// elastic_genTerm - A linear solver for elastic problems using FEM & new genTensors
// 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_NGT_H_
#define _ELASTIC_DOMAIN_NGT_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_ngt : public genSupport
{
public :
  genTerm<genTensor<4,double>,0>::ConstHandle C; //elastic tensor

public :
  ElasticDomain_ngt(const genDomain &s)
  {
    dim=s.dim;
    tag=s.tag;
    const genMatrix<double> & mC=s.Tensors.find("C")->second;
    genTensor<4,double>  elast;
    FromVoigtNotation(mC,elast);
    C = genTerm<genTensor<4,double>,0>::Handle(new ConstantField<genTensor<4,double> >(elast));
  }
  template <int n> typename genTerm<genTensor<2,double>,n>::Handle Stress(const typename genTerm<genTensor<2,double>,n>::ConstHandle &Strain);
  template <int n> genTerm<genTensor<4,double>,0>::ConstHandle Tangent(const typename genTerm<genTensor<2,double>,n>::ConstHandle &Strain);
  template <int n1, int n2> typename genTerm<genTensor<0,double>, n1+n2>::Handle
  Energy( const typename genTerm<genTensor<2,double>,n1>::ConstHandle &Strain1, const typename genTerm<genTensor<2,double>,n2>::ConstHandle &Strain2 );
  ElasticDomain_ngt() : genSupport() {}
  virtual ~ElasticDomain_ngt() {};
};

template <int n> typename genTerm<genTensor<2,double>,n>::Handle ElasticDomain_ngt::Stress(const typename genTerm<genTensor<2,double>,n>::ConstHandle &Strain)
{
  return typename genTerm<genTensor<2,double>,n>::Handle(Compose<FullContrProdOp>(C,Strain));
  //return typename genTerm<genTensor<2,double>,n>::Handle(FullContractedProductCompose(C,Strain));
}

template <int n> genTerm<genTensor<4,double>,0>::ConstHandle ElasticDomain_ngt::Tangent(const typename genTerm<genTensor<2,double>,n>::ConstHandle &Strain)
{
  return genTerm<genTensor<4,double>,0>::ConstHandle(C);
}

template <int n1, int n2> typename genTerm<genTensor<0,double>,n1+n2>::Handle
ElasticDomain_ngt::Energy( const typename genTerm<genTensor<2,double>,n1>::ConstHandle &Strain1, const typename genTerm<genTensor<2,double>,n2>::ConstHandle &Strain2)
{
  return typename genTerm<genTensor<0,double>,n1+n2>::Handle (Compose<FullContrProdOp>(Transpose(Strain1),Stress<n2>(Strain2)));
//  return typename genTerm<double,n1+n2>::Handle (FullContractedProductCompose(Transpose(Strain1),Stress<n2>(Strain2)));
}

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

public :
  IsotropicElasticDomain_ngt(const genDomain &s)
  {
    dim=s.dim;
    tag=s.tag;
    const double & E_=s.Scalars.find("E")->second;
    const double & nu_=s.Scalars.find("nu")->second;
    E = E_; nu = nu_;
    update();
  }
  IsotropicElasticDomain_ngt () : ElasticDomain_ngt() {}
  ~IsotropicElasticDomain_ngt () {}
  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;

    genTensor<4,double>  elast;
    FromVoigtNotation(mC,elast);

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

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

public :
  OrthotropicElasticDomain_ngt(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_ngt () :  ElasticDomain_ngt() {}
  ~OrthotropicElasticDomain_ngt () {}
  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];

    genTensor<4,double>  elast;
    FromVoigtNotation(mC,elast);

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


#endif// _ELASTIC_DOMAIN_NGT_H_