// Periodic_genTerm - A solver for periodic problems using FEM
// Copyright (C) 2010-2026 Eric Bechet, Frederic Duboeuf
//
// See the LICENSE file for license information and contributions.
// Please report all bugs and problems to <bechet@cadxfem.org> or <duboeuf@outlook.com>.

#include "periodicSolver.h"
#include "genAlgorithms.h"
#include "genSField.h"

template<class Iterator, class T> void AssembleInteg(const genTerm<T,0> &term,
                                                     Iterator itbegin, Iterator itend,
                                                     QuadratureBase &integrator, T &val)
{
  for (Iterator it = itbegin; it != itend; ++it)
  {
    MElement* e = *it;
    AssembleInteg(term,e,integrator,val);
  }
}

template<class T> void AssembleInteg(const genTerm<T,0> &term, MElement* e,
                                     QuadratureBase &integrator, T &val)
{
  genScalar<T> localval;
  IntPt* GP;
  int npts = integrator.getIntPoints(e, &GP);
  term.get(e, npts, GP, localval);
  val += localval;
}




PeriodicSolver::~PeriodicSolver()
{
}

std::vector<genBoundaryCondition*> PeriodicSolver::copyBCs()
{
  std::vector<genBoundaryCondition*> fullBCs=Data->BCs;
  return fullBCs;
}

void PeriodicSolver::selectBCs(const std::string &fileName, int cas, std::vector<genBoundaryCondition*> &fullBCs)
{
  // copy boundary conditions
  Data->BCs=fullBCs;
  printf("--> %d BCs copy \n", (int)Data->BCs.size());
  int NbFaces = Data->User.Integers["NbFaces"];
  //cas = which BC, select BC associated
  if(cas==0)
  {
    Data->BCs.resize(NbFaces);
    printf("--> %d BCs calcul \n", (int)Data->BCs.size());
  }
  else if(cas>0)
  {
  for ( unsigned int i=0 ; i < (cas*NbFaces); i++ )
  {
    Data->BCs.erase(Data->BCs.begin());
  }
  Data->BCs.resize(NbFaces);
  printf("--> %d BCs calcul \n", (int)Data->BCs.size());
  }
}

void PeriodicSolver::buildStiffness(int cas, genMatrix<double> &Stiff)
{
  double volume=0;
  GaussQuadrature Integ_Bulk ( GaussQuadrature::GradGrad );
  for ( unsigned int i = 0; i < EDomains.size(); i++ )
  {
    ConstantField<genTensor0<double> > Volume(1.);
    Assemble ( Volume, EDomains[i]->begin(), EDomains[i]->end(), Integ_Bulk, volume );
  }
  printf ( "volume=%f\n",volume );


  genTensor2<double> stressH,strainH;
  genTensor0<double> energyH=0.;
  for ( unsigned int i = 0; i < EDomains.size(); i++ )
  {
    ElasticDomain* E = EDomains[i];
    diffTerm<genTensor1<double>,0>::Handle Field (new genSField<genTensor1<double> >(pAssembler, FSpaceDisp) );

    genTerm<genTensor2<double>,0>::Handle G = Gradient(Field);
    genTerm<genTensor2<double>,0>::Handle Strain = (G+Transpose(G))*0.5;

    genTerm<genTensor2<double>,0>::Handle Stress = E->Stress<>(Strain);
//     genTerm<genTensor4<double>,0>::Handle C = E->Tangent<0>(Strain);
//     genTerm<genTensor2<double>,0>::Handle Stress = FullContractedProductCompose(C,Strain);

    genTerm<genTensor0<double>,0>::Handle Energy = E->Energy<>(Strain,Strain);
//     genTerm<double,0>::Handle Energy = FullContractedProductCompose(Transpose(Strain),Stress);

    GaussQuadrature Integ_Grad ( GaussQuadrature::Grad );
    GaussQuadrature Integ_GradGrad ( GaussQuadrature::GradGrad );
    AssembleInteg ( *Strain, E->begin(), E->end(), Integ_Grad, strainH );
    AssembleInteg ( *Stress, E->begin(), E->end(), Integ_Grad, stressH );
    AssembleInteg ( *Energy, E->begin(), E->end(), Integ_GradGrad, energyH );
  }

  double energyEq=0;
  energyEq = (stressH * strainH)/volume;
  strainH.print("strainH");
  stressH.print("stressH");
  printf ( "energyH =%f\n",energyH() );
  printf ( "energyEq=%f\n",energyEq );

  //int cas=0;
  switch (cas)
  {
    case 0 :
      //cout << "C1111 = "<< stressH(0,0)<< endl;
      //cout << "C1122 = "<< stressH(1,1)<< endl;
      //cout << "C1133 = "<< stressH(2,2)<< endl;
      Stiff(0,0) = stressH(0,0)/volume;
      Stiff(0,1) = stressH(1,1)/volume;
      Stiff(1,0) = stressH(1,1)/volume;
      Stiff(0,2) = stressH(2,2)/volume;
      Stiff(2,0) = stressH(2,2)/volume;
      break;
    case 1 :
      //cout << "C2222 = "<< stressH(1,1)<< endl;
      //cout << "C2233 = "<< stressH(2,2)<< endl;
      Stiff(1,1) = stressH(1,1)/volume;
      Stiff(1,2) = stressH(2,2)/volume;
      Stiff(2,1) = stressH(2,2)/volume;
      break;
    case 2 :
      //cout << "C3333 = "<< stressH(2,2)<< endl;
      Stiff(2,2) = stressH(2,2)/volume;
      break;
    case 3 :
      //cout << "C2323 = "<< stressH(1,2)<< endl;
      Stiff(3,3) = stressH(1,2)/(2*volume);
      break;
    case 4 :
      //cout << "C1313 = "<< stressH(1,3)<< endl;
      Stiff(4,4) = stressH(0,2)/(2*volume);
      break;
    case 5 :
      //cout << "C1212 = "<< stressH(0,1)<< endl;
      Stiff(5,5) = stressH(0,1)/(2*volume);
      break;
  }
}