// Incompressible - A linear solver for incompressible 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 _INCOMPDOMAIN_H_
#define _INCOMPDOMAIN_H_

#include <string>
#include <vector>

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


inline double kron(int i,int j)
{
  if (i==j) return 1.0;
  return 0.0 ;
}

class IncompDomain : public genSupport
{
  public :
    genTerm<genTensor4<double>,0>::Handle C; //elastic tensor
    genTerm<genTensor4<double>,0>::Handle CmCtilde; // obtained from C or E and nu
    genTerm<genTensor2<double>,0>::Handle Itilde; // obtained from C or E and nu
    genTerm<genTensor0<double>,0>::Handle invk; //inverse of the compressibility obtained from C or E and nu
    double factor; // <>1 to avoid an ill-conditionned linear system
  public :
    IncompDomain(const genDomain &s) : genSupport(s)
    {
      std::map<std::string,genMatrix<double> >::const_iterator it;
      it=s.Tensors.find("C");

      if (it!=s.Tensors.end())
      {
        factor=1.;
        const genMatrix<double> & mC=it->second;
        genTensor4<double> elast;
        genTensor4<double> elast_til;
        genTensor4<double> elast_m_til;
        genTensor2<double> elast_itil;
        genTensor2<double> elast_itil2;
        double kk;
        double Caabb=0;
        FromVoigtNotation(mC,elast);

        for (int a=0;a<3;++a)
          for (int b=0;b<3;++b)
            Caabb+=elast(a,a,b,b);
        std::cout << "compressibility=" << 9./Caabb << std::endl;
        double norm=0.;
        for (int i=0;i<3;++i)
          for (int j=0;j<3;++j)
            for (int k=0;k<3;++k)
              for (int l=0;l<3;++l)
              {
                elast_til(i,j,k,l)=0.;
                for (int a=0;a<3;++a)
                  for (int b=0;b<3;++b)
                    elast_til(i,j,k,l)+=elast(a,a,i,j)*elast(b,b,k,l)/Caabb;
              }

        for (int i=0;i<3;++i)
          for (int j=0;j<3;++j)
          {
            elast_itil(i,j)=0.;
            for (int a=0;a<3;++a)
              elast_itil(i,j)+=3*elast(a,a,i,j)/Caabb;
          }

        elast_m_til=elast+elast_til*(-1.);
        factor=0;
        for (int i=0;i<3;++i)
          for (int j=0;j<3;++j)
            for (int k=0;k<3;++k)
              for (int l=0;l<3;++l)
              {
                if ((i!=j)&&(k!=l)) factor+=elast_m_til(i,j,k,l);
              }
        
        for (int i=0;i<3;++i)
          for (int j=0;j<3;++j)
             elast_itil2(i,j)+=elast_itil(i,j)*factor;
   
        
        C=genTerm<genTensor4<double>,0>::Handle(new ConstantField<genTensor4<double> >(elast));
        CmCtilde= genTerm<genTensor4<double>,0>::Handle(new ConstantField<genTensor4<double> >(elast_m_til));
        Itilde= genTerm<genTensor2<double>,0>::Handle(new ConstantField<genTensor2<double> >(elast_itil2));
        invk=genTerm<genTensor0<double>,0>::Handle(new ConstantField<genTensor0<double> >(factor*factor*9/Caabb));
      }
      else
      {
        double E = s.Scalars.find("E")->second;
        double nu = s.Scalars.find("nu")->second;
        genTensor2<double> elast_itil;
        genTensor2<double> elast_itil2;
        genTensor4<double> elast_m_til;
        for (int i=0;i<3;++i)
          for (int j=0;j<3;++j)
          {
            elast_itil(i,j)=kron(i,j);
          }

          
        double mu=E/(2*(1+nu));
        factor=mu;
//        factor=1.0;
        for (int i=0;i<3;++i)
          for (int j=0;j<3;++j)
            for (int k=0;k<3;++k)
              for (int l=0;l<3;++l)
              {
                elast_m_til(i,j,k,l)=mu*(kron(i,k)*kron(j,l)+kron(i,l)*kron(j,k)-kron(i,j)*kron(k,l)*(2./3.));
              }
        for (int i=0;i<3;++i)
          for (int j=0;j<3;++j)
          {
            elast_itil2(i,j)=elast_itil(i,j)*factor;
          }
                
                
          std::cout << "compressibility=" << 3*(1-2*nu)/E << std::endl;
          
          Itilde= genTerm<genTensor2<double>,0>::Handle(new ConstantField<genTensor2<double> >(elast_itil2));
          invk=genTerm<genTensor0<double>,0>::Handle(new ConstantField<genTensor0<double> >(factor*factor*3*(1-2*nu)/E));
          CmCtilde= genTerm<genTensor4<double>,0>::Handle(new ConstantField<genTensor4<double> >(elast_m_til));
      }
    }
    /*
    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< 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< 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
          );*/
    IncompDomain() : genSupport() {}
    virtual ~IncompDomain() {}
};

/*
template <int n> typename genTerm<genTensor2<double>,n>::Handle
PiezoSupport::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>(EField,E)*-1.0);
  // now working !!
}

template <int n> typename genTerm<genTensor1<double>,n>::Handle
PiezoSupport::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>(EField,K));
// now working !!
}

template <int  n1, int n2> typename genTerm< double, n1 + n2  >::Handle 
PiezoSupport::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<double,n1+n2>::Handle Mech(Compose<FullContrProdOp>(Strain1,Stress<n2>(Strain2,EField2)));
    typename genTerm<double,n1+n2>::Handle Elec(Compose<FullContrProdOp>(EField1,Edisp<n2>(Strain2,EField2)));
    return Mech+Elec*-1.0;
}


template <int  n1, int n2> typename genTerm< double, n1 + n2  >::Handle 
PiezoSupport::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<double,n1+n2>::Handle Mech(Compose<FullContrProdOp>(Strain1,Stress<n2>(Strain2,EField2)));
    typename genTerm<double,n1+n2>::Handle Elec(Compose<FullContrProdOp>(EField1,Edisp<n2>(Strain2,EField2)));
    return Mech+Elec;
}
*/
#endif// _INCOMPDOMAIN_H_