// elastic_genTerm - A linear solver for elastic problems using FEM
// 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 _ELASTIC_DOMAIN_H_
#define _ELASTIC_DOMAIN_H_

#include <string>
#include <vector>

#include "genSupports.h"
#include "genTerm.h"
#include "genTensors.h"
#include "genDataIO.h"
#include "genMatrix.h"
#include "savedGenTerm.h"
// #include "genTermCompose.h"
// #include "genTensorFullContrProd.h"

class ElasticDomain : public genSupport
{
public :
  genTerm<genTensor4<double>,0>::ConstHandle C; //elastic tensor

public :
  ElasticDomain(const genDomain &s)
  {
    dim=s.dim;
    tag=s.tag;
    const genMatrix<double> & mC=s.Tensors.find("C")->second;
    genTensor4<double>  elast;
    FromVoigtNotation(mC,elast);
    C = genTerm<genTensor4<double>,0>::Handle(new ConstantField<genTensor4<double> >(elast));
  }
  template <class GT> typename genTerm<typename GT::ValType, GT::Nsp>::Handle Stress(const std::shared_ptr< GT> &Strain);
  template <int n> genTerm<genTensor4<double>,0>::ConstHandle Tangent(const typename genTerm<genTensor2<double>,n>::ConstHandle &Strain);
  template <class GT1, class GT2> typename genTerm<typename FullContrProdOp< typename GT1::ValType, typename GT2::ValType>::ValType, GT1::Nsp+GT2::Nsp>::Handle
  Energy(const std::shared_ptr< GT1> &Strain1, const std::shared_ptr< GT2> &Strain2);
  ElasticDomain() : genSupport() {}
  virtual ~ElasticDomain() {};
};

template <class GT> typename genTerm<typename GT::ValType, GT::Nsp>::Handle ElasticDomain::Stress(const std::shared_ptr<GT> &Strain)
{
  return typename genTerm<typename GT::ValType, GT::Nsp>::Handle(Compose<FullContrProdOp>(C,Strain));
  //return typename genTerm<genTensor2<double>,n>::Handle(FullContractedProductCompose(C,Strain));
}

template <int n> genTerm<genTensor4<double>,0>::ConstHandle ElasticDomain::Tangent(const typename genTerm<genTensor2<double>,n>::ConstHandle &Strain)
{
  return genTerm<genTensor4<double>,0>::ConstHandle(C);
}

template <class GT1, class GT2> typename genTerm<typename FullContrProdOp< typename GT1::ValType, typename GT2::ValType>::ValType, GT1::Nsp + GT2::Nsp>::Handle
ElasticDomain::Energy( const std::shared_ptr<GT1> &Strain1, const std::shared_ptr<GT2> &Strain2)
{
  return typename genTerm<typename FullContrProdOp< typename GT1::ValType, typename GT2::ValType>::ValType ,GT1::Nsp+GT2::Nsp>::Handle (Compose<FullContrProdOp>(Transpose(Strain1),Stress<>(Strain2)));
//  return typename genTerm<double,n1+n2>::Handle (FullContractedProductCompose(Transpose(Strain1),Stress<n2>(Strain2)));
}

class IsotropicElasticDomain : public ElasticDomain
{
public :
  double E, nu; // specific elastic datas (isotropic)

public :
  IsotropicElasticDomain(const genDomain &s)
  {
    dim=s.dim;
    tag=s.tag;
    if(s.Scalars.find("E") != s.Scalars.end() && s.Scalars.find("nu") != s.Scalars.end() ){
	E=s.Scalars.find("E")->second;
	nu=s.Scalars.find("nu")->second;
    }
    else if(s.Scalars.find("lambda") != s.Scalars.end() && s.Scalars.find("mu") != s.Scalars.end() ){
      double lambda = s.Scalars.find("lambda")->second;
      double mu = s.Scalars.find("mu")->second;
      E = mu*(3*lambda+2*mu)/(lambda +mu);
      nu = lambda/(2*(lambda+mu));
    }
    update();
  }
  IsotropicElasticDomain () : ElasticDomain() {}
  ~IsotropicElasticDomain () {}
  void update ()
  {
    double FACT = E / (1 + nu);
    double C11 = FACT * (1 - nu) / (1 - 2 * nu);
    double C12 = FACT * nu / (1 - 2 * nu);
    double C44 = (C11 - C12) / 2;
//       FACT = E / (1 - nu * nu); // plane stress (plates)
//       C11  = FACT;
//       C12  = nu * FACT;
//       C44 = (1. - nu) * .5 * FACT;
    genMatrix<double> mC(6,6);
    for(int i = 0; i < 3; ++i)
    {
      mC(i,i) = C11;
      mC(i+3,i+3) = C44;
    }
    mC(1,0) = mC(0,1) = mC(2,0) = mC(0,2) = mC(1,2) = mC(2,1) = C12;

    genTensor4<double>  elast;
    FromVoigtNotation(mC,elast);

    C = genTerm<genTensor4<double>,0>::Handle(new ConstantField<genTensor4<double> >(elast));
  }
};

class OrthotropicElasticDomain : public ElasticDomain
{
public :
  std::vector<double> E, nu, G; // E1,E2,E3, nu12,nu13,nu23, G12,G13,G23 (orthotropic_3D)

public :
  OrthotropicElasticDomain(const genDomain &s)
  {
    dim=s.dim;
    tag=s.tag;
    const std::vector<double> & E_=s.Vectors.find("E")->second;
    const std::vector<double> & nu_=s.Vectors.find("nu")->second;
    const std::vector<double> & G_=s.Vectors.find("G")->second;
    E = E_; nu = nu_; G = G_;
    update();
  }
  OrthotropicElasticDomain () :  ElasticDomain() {}
  ~OrthotropicElasticDomain () {}
  void update ()
  {
    double S11 = 1. / E[0];
    double S22 = 1. / E[1];
    double S33 = 1. / E[2];
    double S12 = - nu[0] / E[0];
    double S13 = - nu[1] / E[1];
    double S23 = - nu[2] / E[2];

    double dS = S11*S22*S33 - S11*S23*S23 - S22*S13*S13 - S33*S12*S12 + 2*S12*S23*S13;

    genMatrix<double> mC(6,6);
    mC(0,0) = (S22*S33 - S23*S23)/dS;
    mC(1,1) = (S33*S11 - S13*S13)/dS;
    mC(2,2) = (S11*S22 - S12*S12)/dS;
    mC(0,1) = mC(1,0) = (S13*S23 - S12*S33)/dS;
    mC(0,2) = mC(2,0) = (S12*S23 - S13*S22)/dS;
    mC(1,2) = mC(2,1) = (S12*S13 - S23*S11)/dS;
    mC(3,3) = G[2];
    mC(4,4) = G[1];
    mC(5,5) = G[0];

    genTensor4<double>  elast;
    FromVoigtNotation(mC,elast);

    C = genTerm<genTensor4<double>,0>::Handle(new ConstantField<genTensor4<double> >(elast));
  }
};


#endif// _ELASTIC_DOMAIN_H_