// Piezo - A linear solver for piezo problem 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 _PIEZODOMAIN_H_
#define _PIEZODOMAIN_H_

#include <string>
#include <vector>

#include "genSupports.h"
#include "genTerm.h"
#include "genTensors.h"
#include "genDataIO.h"
#include "genMatrix.h"

class PiezoDomain : public genSupport
{
  public :
    genTerm<genTensor4<double>,0>::Handle C; //elastic tensor
    genTerm<genTensor3<double>,0>::Handle E; //piezoelectric tensor
    genTerm<genTensor2<double>,0>::Handle K; //dielectric tensor

  public :
    PiezoDomain(const genDomain &s) : genSupport(s)
    {
      const genMatrix<double> & mC=s.Tensors.find("C")->second;
      const genMatrix<double> & mE=s.Tensors.find("E")->second;
      const genMatrix<double> & mK=s.Tensors.find("K")->second;
      genTensor4<double> elast;
      genTensor3<double> piezo;
      genTensor2<double> dielec;
      FromVoigtNotation(mC,elast);
      FromVoigtNotation(mE,piezo);
      FromVoigtNotation(mK,dielec);
      C=genTerm<genTensor4<double>,0>::Handle(new ConstantField<genTensor4<double> >(elast));
      E=genTerm<genTensor3<double>,0>::Handle(new ConstantField<genTensor3<double> >(piezo));
      K=genTerm<genTensor2<double>,0>::Handle(new ConstantField<genTensor2<double> >(dielec));
    }
    template <int n> typename genTerm<genTensor2<double>,n>::Handle Stress(const typename genTerm<genTensor2<double>,n>::ConstHandle &Strain,const typename genTerm<genTensor1<double>,n>::ConstHandle &EField);
    template <int n> typename genTerm<genTensor1<double>,n>::Handle Edisp(const typename genTerm<genTensor2<double>,n>::ConstHandle &Strain,const typename genTerm<genTensor1<double>,n>::ConstHandle &EField);
    template <int  n1, int n2> typename genTerm<genTensor0<double>, n1 + n2  >::Handle 
    Enthalpy( const typename genTerm<genTensor2<double>,n1>::ConstHandle &Strain1,const typename genTerm<genTensor1<double>,n1>::ConstHandle &EField1,
              const typename genTerm<genTensor2<double>,n2>::ConstHandle &Strain2,const typename genTerm<genTensor1<double>,n2>::ConstHandle &EField2
            );
    template <int  n1, int n2> typename genTerm<genTensor0<double>, n1 + n2  >::Handle             
    Energy( const typename genTerm<genTensor2<double>,n1>::ConstHandle &Strain1,const typename genTerm<genTensor1<double>,n1>::ConstHandle &EField1,
            const typename genTerm<genTensor2<double>,n2>::ConstHandle &Strain2,const typename genTerm<genTensor1<double>,n2>::ConstHandle &EField2
          );

    PiezoDomain() : genSupport() {}
    virtual ~PiezoDomain() {};
};


template <int n> typename genTerm<genTensor2<double>,n>::Handle
PiezoDomain::Stress(const typename genTerm<genTensor2<double>,n>::ConstHandle &Strain,const typename genTerm<genTensor1<double>,n>::ConstHandle &EField)
{
  return typename genTerm<genTensor2<double>,n>::Handle(Compose<FullContrProdOp>(C,Strain)+Compose<FullContrProdOp>(Transpose(E),EField)*-1.0);
}

template <int n> typename genTerm<genTensor1<double>,n>::Handle
PiezoDomain::Edisp(const typename genTerm<genTensor2<double>,n>::ConstHandle &Strain,const typename genTerm<genTensor1<double>,n>::ConstHandle &EField)
{
  return typename genTerm<genTensor1<double>,n>::Handle(Compose<FullContrProdOp>(E,Strain)+Compose<FullContrProdOp>(K,EField));
}

template <int  n1, int n2> typename genTerm< genTensor0<double>, n1 + n2  >::Handle 
PiezoDomain::Enthalpy( const typename genTerm<genTensor2<double>,n1>::ConstHandle &Strain1,const typename genTerm<genTensor1<double>,n1>::ConstHandle &EField1,
              const typename genTerm<genTensor2<double>,n2>::ConstHandle &Strain2,const typename genTerm<genTensor1<double>,n2>::ConstHandle &EField2)
{
  typename genTerm<genTensor0<double>,n1+n2>::Handle Mech(Compose<FullContrProdOp>(Transpose(Strain1),Stress<n2>(Strain2,EField2)));
  typename genTerm<genTensor0<double>,n1+n2>::Handle Elec(Compose<FullContrProdOp>(Transpose(EField1),Edisp<n2>(Strain2,EField2)));
  return Mech+Elec*-1.0;
}


template <int  n1, int n2> typename genTerm< genTensor0<double>, n1 + n2  >::Handle 
PiezoDomain::Energy( const typename genTerm<genTensor2<double>,n1>::ConstHandle &Strain1,const typename genTerm<genTensor1<double>,n1>::ConstHandle &EField1,
            const typename genTerm<genTensor2<double>,n2>::ConstHandle &Strain2,const typename genTerm<genTensor1<double>,n2>::ConstHandle &EField2)

{
  typename genTerm<genTensor0<double>,n1+n2>::Handle Mech(Compose<FullContrProdOp>(Transpose(Strain1),Stress<n2>(Strain2,EField2)));
  typename genTerm<genTensor0<double>,n1+n2>::Handle Elec(Compose<FullContrProdOp>(Transpose(EField1),Edisp<n2>(Strain2,EField2)));
  return Mech+Elec;
}

#endif// _PIEZODOMAIN_H_