// GenFem - A high-level finite element library
// Copyright (C) 2010-2026 Eric Bechet
//
// See the LICENSE file for license information and contributions.
// Please report all bugs and problems to <bechet@cadxfem.org>.

#ifndef _GEN_SIMPLE_FUNCTION_SIGMA_INC_H_
#define _GEN_SIMPLE_FUNCTION_SIGMA_INC_H_


#include "genSimpleFunction.h"
#include "genSimpleFunctionExtended.h"
#include <boost/math/special_functions/bessel.hpp>
#include <cmath>

template<class T1,class T2> class genSimpleFunctionSigmaInc : public genSimpleFunction<T1>
{
public :
  typedef typename TensorialTraits<T1>::Scalar Scalar;
  typedef typename TensorialTraits<T1>::ScalarType ScalarType;
  typedef typename TensorialTraits<T1>::ValType ValType;
  typedef typename TensorialTraits<T1>::GradType GradType;
  typedef typename TensorialTraits<T1>::HessType HessType;
    
  genSimpleFunction<T2> *polar;
  double kp, lambda, mu;
  std::vector<double>::size_type size;
public :
  genSimpleFunctionSigmaInc( std::vector<double>::size_type size_, double &kp_, double &lambda_, double & mu_ ) : kp(kp_), lambda(lambda_), mu(mu_), size(size_) {
        polar = new genSimpleFunctionCylindricalCoordSyst<T2>();
  };
  virtual ValType operator()(double x, double y, double z) const {
   
    std::vector<double> polar_coor = vectorize((*polar)(x,y,z));
    double r = polar_coor[0];
    double theta = polar_coor[1];
    Scalar sigma_rr, sigma_rtheta, sigma_x, sigma_y;  
    std::complex<double> i(0,1);
    for(int n=0; n< (int)size; n++){
    if(n==0){
    	sigma_rr += pow(i,n)*kp*( (-lambda*kp*kp-2*mu/(r*r)*(n-n*n+(kp*r)*(kp*r)) )*boost::math::cyl_bessel_j(n,kp*r) + 2*mu*kp/r*boost::math::cyl_bessel_j(n+1,kp*r) )*cos(n*theta);

	sigma_rtheta += (-2*n/r*mu*pow(i,n)*(-(1-n)/r*boost::math::cyl_bessel_j(n,kp*r) -kp*boost::math::cyl_bessel_j(n+1,kp*r) ) )*sin(n*theta);
    }
    else{
	sigma_rr += 2.0*pow(i,n)*kp*( (-lambda*kp*kp-2*mu/(r*r)*(n-n*n+(kp*r)*(kp*r)) )*boost::math::cyl_bessel_j(n,kp*r) + 2*mu*kp/r*boost::math::cyl_bessel_j(n+1,kp*r) )*cos(n*theta);

	sigma_rtheta += (-4.0*n/r*mu*pow(i,n)*(-(1-n)/r*boost::math::cyl_bessel_j(n,kp*r) -kp*boost::math::cyl_bessel_j(n+1,kp*r) ) )*sin(n*theta);
	} 
    }
    sigma_x = sigma_rr*cos(theta) - sigma_rtheta*sin(theta);
    sigma_y = sigma_rr*sin(theta) + sigma_rtheta*cos(theta);
    std::vector<Scalar> res;
    res.push_back(sigma_x);
    res.push_back(sigma_y);
    return ValType(res);
  };
  
  virtual GradType grad(double x, double y, double z) const { return GradType(); };
  virtual HessType hess(double x, double y, double z) const { return HessType(); };
  ~genSimpleFunctionSigmaInc(){ delete polar;}
};



#endif // _GEN_SIMPLE_FUNCTION_UINC_H_