// Autodiff - automatic differentiation tools for GenFem
// Copyright (C) 2012-2026 Eric Bechet
//
// See the LICENSE file for license information and contributions.
// Please report all bugs and problems to <bechet@cadxfem.org>.

#include "dfloat.h"
#include "genTensors.h"
#include "genTraits.h"
#include "genFSpace.h"
#include "genTerm.h"
#include "genDataIO.h"
#include "genDeriveOp.h"
#include "genSolver.h"
#include "GmshGlobal.h"
#include "GModel.h"
#include "genGroupOfElements.h"
#include "linearSystemCSR.h"
#include "genAlgorithms.h"
#include "genTermCompose.h"

#ifdef USE_FTENSOR
#include "FTensor.hpp"
using namespace FTensor;
#endif // USE_FTENSOR


#include <iostream>
using namespace std;

dfloat<double> f(dfloat<double> in)
{
  return in*in+in-100.;
}

dfloat<double> functionA(dfloat<double> x,dfloat<double> y,dfloat<double> z)
{
  return sqrt(x*x+y*y+z*z)/(x+y+z);
}

dfloat<double> functionX(dfloat<double> x,dfloat<double> y,dfloat<double> z)
{
  return sqrt(x*x+y*y+z*z)/(x+y+z);
}

dfloat<double> functionY(dfloat<double> x,dfloat<double> y,dfloat<double> z)
{
  return sqrt(x*x-y*y+z*z)/(x+y+z);
}

dfloat<double> functionZ(dfloat<double> x,dfloat<double> y,dfloat<double> z)
{
  return sqrt(x*x+y*y-z*z)/(x+y+z);
}


#ifdef USE_FTENSOR



template <class TA,class TB,class TC,char T3,char T4,int N>
class Op
{
public:
  Op()  {}
  void operator()(const TA &tens1, const TB &tens2, TC &res)
  {
    Index<'i',N> t1;
    Index<'j',N> t2;
    Index<T3,N> t3;
    Index<T4,N> t4;
    res(t1,t2)=tens1(t3)*tens2(t4);
  }
  void operator()(const TA &tens1, TC &res)
  {
    Index<'i',N> t1;
    Index<'j',N> t2;
    Index<T3,N> t3;
    Index<T4,N> t4;
    res(t1,t2)=tens1(t3)*tens1(t4);
  }
};
#endif // USE_FTENSOR


#include "genSimpleFunctionExtended.h"

int main(int argc, char** argv)
{
/*  std::vector<std::string> expr(1), var(3);
  expr[0]="1.*x+2.*y+3.*z";
  var[0]="x"; var[1]="y"; var[2]="z";
  genSimpleFunctionAnalytical<genTensor0<double> > f(expr,expr,expr,expr,var);

  double x=-1., y=1., z=1.;
  std::cout << "Test Mathematical function : f(x,y,z) = " << f(x,y,z) << std::endl;
  std::cout << "Test Derivative function : df(x,y,z) = " << f.grad(x,y,z) << std::endl;
  std::cout << "Test output : std::cout << f = " << f << std::endl;

  genSimpleFunctionOperatorProduct<genTensor0<double> > prodf(2., &f);
  std::cout << "Test Product : 2f(x,y,z) = " << prodf(x,y,z) << std::endl;

  genSimpleFunctionOperatorTensorize<genTensor1<double> > Tf(&f,&prodf,&f);
  std::cout << "Test Tensorize : [f,2f,f](x,y,z) = " << Tf(x,y,z) << std::endl;

  genSimpleFunctionOperatorComposition<genTensor0<double>,genTensor1<double> > cf(&f,&Tf);
  std::cout << "Test Composition : (foTf)(x,y,z) = " << cf(x,y,z) << std::endl;

  //Cylindrical Coordinates System
  genSimpleFunctionCylindricalCoordSyst<genTensor1<double> >fCylindricalCS;
  std::cout << "Test CylindricalCoordSyst : fCylindricalCS(x,y,z) = " << fCylindricalCS(x,y,z) << std::endl;

  std::vector<std::string> expr1(3), var1(3);
  expr1[0]="rho*cos(theta)"; expr1[1]="rho*sin(theta)"; expr1[2]="z";
  var1[0]="rho"; var1[1]="theta"; var1[2]="z";
  genSimpleFunctionAnalytical<genTensor1<double> > f1(expr1,var1);

  genSimpleFunctionOperatorComposition<genTensor1<double>,genTensor1<double> > f1_o_fCylindricalCS(&f1,&fCylindricalCS);
  std::cout << "Test Function in Cylindrical Coord Syst : f1_o_fCylindricalCS(x,y,z) = " << f1_o_fCylindricalCS(x,y,z) << std::endl;

  //Spherical Coordinates System
  genSimpleFunctionSphericalCoordSyst<genTensor1<double> >fSphericalCS;
  std::cout << "Test SphericalCoordSyst : fSphericalCS(x,y,z) = " << fSphericalCS(x,+y,z) << std::endl;

  std::vector<std::string> expr3(3), var3(3);
  expr3[0]="r*cos(theta)*sin(phi)"; expr3[1]="r*sin(theta)*sin(phi)"; expr3[2]="r*cos(phi)";
  var3[0]="r"; var3[1]="theta"; var3[2]="phi";
  genSimpleFunctionAnalytical<genTensor1<double> > f2(expr3,var3);

  genSimpleFunctionOperatorComposition<genTensor1<double>,genTensor1<double> > f2_o_fSphericalCS(&f2,&fSphericalCS);
  std::cout << "Test Function in Spherical Coord Syst : f2_o_fSphericalCS(x,y,z) = " << f2_o_fSphericalCS(x,y,z) << std::endl;
*/

/*  std::vector<double> v(1,9.);
  double d;
  tensorize(v,d);
  std::cout << d << std::endl;

  genTensor0<double> d0;
  tensorize(v,d0);
  std::cout << d0 << std::endl;

  std::vector<double> v1(3,9.);
  genTensor1<double> d1;
  tensorize(v1,d1);
  std::cout << d1 << std::endl;

  std::vector<double> v2(9,9.);
  genTensor2<double> d2;
  tensorize(v2,d2);
  std::cout << d2 << std::endl;

  std::vector<double> v3(27,9.);
  genTensor3<double> d3;
  tensorize(v3,d3);
  std::cout << d3 << std::endl;

  std::vector<double> v4(81,9.);
  genTensor4<double> d4;
  tensorize(v4,d4);
  std::cout << d4 << std::endl;*/
#if 1
  
  for (unsigned long i=0;i<524;++i)
  {
    std::cout << i << " " << nbits(i) << std::endl;
  }
  DofIdManager* ddd = DofIdManager::getInstance();
  int iField=0, iComp=0, iIndex=1, iLocalIndex=0;
//   int iField=32767, iComp=15, iIndex=15, iLocalIndex=255;
//   int iField=32767, iComp=255, iIndex=255;
//   int iType;
//   iType = ddd->createTypeWithTwoInts(iField,iComp);
//   std::cout << "Create Type With Two Ints : " << " (" << iField << "," << iComp << ") => " << iType << std::endl;
//   ddd->getTwoIntsFromType(iType, iField, iComp);
//   std::cout << "Get Two Ints From Type : " << iType << " => (" << iField << "," << iComp << ")" << std::endl;
//   std::cout << std::endl;

//   iType = ddd->createTypeWithThreeInts(iField, iComp, iIndex);
//   std::cout << "Create Type With Three Ints : " << " (" << iField << "," << iComp << "," << iIndex << ") => " << iType << std::endl;
//   ddd->getThreeIntsFromType(iType, iField, iComp, iIndex);
//   std::cout << "Get Three Ints From Type : " << iType << " => (" << iField << "," << iComp << "," << iIndex << ")" << std::endl;
//   std::cout << std::endl;
// 
//   iType = ddd->createTypeWithFourInts(iField, iComp, iIndex, iLocalIndex);
//   std::cout << "Create Type With Four Ints : " << " (" << iField << "," << iComp << "," << iIndex << "," << iLocalIndex << ") => " << iType << std::endl;
//   ddd->getFourIntsFromType(iType, iField, iComp, iIndex, iLocalIndex);
//   std::cout << "Get Four Ints From Type : " << iType << " => (" << iField << "," << iComp << "," << iIndex << "," << iLocalIndex << ")" << std::endl;
//   std::cout << std::endl;
// 
//   ddd->setBinaryShift(511,1023,15);
//   int Shift_Comp, Shift_Index, Shift_LocalIndex;
//   ddd->getBinaryShift(Shift_Comp, Shift_Index, Shift_LocalIndex);
//   std::cout << "Binary Shift : " << "Shift_Comp=" << Shift_Comp << ", Shift_Index=" << Shift_Index << ", Shift_LocalIndex=" << Shift_LocalIndex << std::endl;
//   int Max_Field, Max_Comp, Max_Index, Max_LocalIndex;
//   ddd->getMaxValues(Max_Field, Max_Comp, Max_Index, Max_LocalIndex);
//   std::cout << "Max Values : " << "Max_Field=" << Max_Field << ", Max_Comp=" << Max_Comp << ", Max_Index=" << Max_Index << ", Max_LocalIndex=" << Max_LocalIndex << std::endl;
//   std::cout << std::endl;
// 
//   iField=255, iComp=511, iIndex=1023, iLocalIndex=15;
//   iType = ddd->createTypeWithFourInts(iField, iComp, iIndex, iLocalIndex);
//   std::cout << "Create Type With Four Ints : " << " (" << iField << "," << iComp << "," << iIndex << "," << iLocalIndex << ") => " << iType << std::endl;
//   ddd->getFourIntsFromType(iType, iField, iComp, iIndex, iLocalIndex);
//   std::cout << "Get Four Ints From Type : " << iType << " => (" << iField << "," << iComp << "," << iIndex << "," << iLocalIndex << ")" << std::endl;
//   std::cout << std::endl;

  std::cout << std::endl;
  std::cout << "Create Field n°" << ddd->createField()<< std::endl;
  std::cout << "Create Field n°" << ddd->createField()<< std::endl;
  std::cout << "Create Field n°" << ddd->createField()<< std::endl;
  std::cout << "Create Field n°" << ddd->createField()<< std::endl;
  std::cout << "Create Field n°" << ddd->createField()<< std::endl;
  std::cout << "Create Field n°" << ddd->createField()<< std::endl;
  std::cout << std::endl;
  std::cout << "Create Index n°" << ddd->createIndex(1) << " in Field n°" << 1 << std::endl;
  std::cout << "Create Index n°" << ddd->createIndex(2) << " in Field n°" << 2 << std::endl;
  std::cout << "Create Index n°" << ddd->createIndex(2) << " in Field n°" << 2 << std::endl;
  std::cout << "Create Index n°" << ddd->createIndex(3) << " in Field n°" << 3 << std::endl;
  std::cout << "Create Index n°" << ddd->createIndex(3) << " in Field n°" << 3 << std::endl;
  std::cout << "Create Index n°" << ddd->createIndex(3) << " in Field n°" << 3 << std::endl;
  std::cout << std::endl;
  std::cout << *ddd;
  std::cout << std::endl;

// #if 1
  dfloat<double> df1(10.0,1.0);
  dfloat<double> df2=f(df1);
//  cout << df2 << endl;
  cout << "f=" << df2.v << " df/dx=" << df2.dv/df1.dv << endl;
  double x(1.0),y(1.0),z(1.0);
//  double diff=diff(x,y,z);
  dfloat<double> res=functionA(x,y,z);
  double grad[3];
  double cur[3];
  double lf=laplacian(functionA,x,y,z);
  double div=divergence(functionX,functionY,functionZ,x,y,z);
  curl(functionX,functionY,functionZ,x,y,z,cur);
  gradient(functionA,x,y,z,grad);
  
//  double lf=diff(complex_function,x,y,z,x);
  
  cout << "f=" << res.v << " df/dx=" << res.dv << endl;
  cout << "gradient D(f)=[" << grad[0] << "," <<grad[1] << "," << grad[2] << "]" << endl;
  cout << "laplacian L(f)=" << lf <<  endl;
  cout << "divergence D(ff)=" << div <<  endl;
  cout << "curl C(ff)=[" << cur[0] << "," <<cur[1] << "," << cur[2] << "]" << endl;
  
  genTensor2<dfloat<double> ,2> mat(1.0);
  
  genTensor1<dfloat<double> ,2> m(1.0);
  genTensor1<dfloat<double> ,2> m2(dfloat<double>(1.0,0.0));
  m2=mat*m;
  m2[0]=dfloat<double>(1.0,1.0);
  m2*=10;
  cout << norm(m) << endl;
  cout << m[0] << endl;
  cout << (norm(m2))<< endl;
  cout << (norm(m2)).dv << endl;

  //  test 

  genFSpace<genTensor0<double> >::Handle FSpaceDisp2(new ScalarLagrangeFSpace<genTensor0<double> >(126));
  genTerm<genTensor0<double>,1>::Handle FSpaceDisp3(FSpaceDisp2);
  genTerm<genTensor1<double>,1>::Handle G2 = Gradient(FSpaceDisp2);
  genFSpace<genTensor1<dfloat<double> > >::Handle FSpaceDisp(new VectorLagrangeFSpace<dfloat<double> >(125));
  genTerm<genTensor1<double>,1>::Handle FSd= Build<genTensor1<double> >(FSpaceDisp); // type cast over the scalar type 
  genTerm<genTensor2<dfloat<double> >,1>::Handle GGGG = Gradient(FSpaceDisp);
  genTerm<genTensor2<dfloat<double> ,2>,1>::Handle GG = Build<genTensor2<dfloat<double>,2 > > (GGGG); // type cast over the tensorial type
  genTerm<genTensor0<double,2>,1>::Handle GGG = Build<genTensor0<double,2 > > (FSpaceDisp3); // type cast over the tensorial type
  genTerm<genTensor1<double,2>,1>::Handle FSd2= Build<genTensor1<double,2> >(FSpaceDisp); // does compile now !!!

  genTerm<genTensor2<dfloat<double> >,1>::Handle G = Gradient(FSpaceDisp);
  dfloat<double> val(465,1);
  genTerm<genTensor0<dfloat<double> >,0>::Handle Mat(new ConstantField<genTensor0<dfloat<double> > >(val));
  genTerm<genTensor0<dfloat<double> >,2>::Handle Term = Compose<FullContrProdOp>(Compose<FullContrProdOp>( Transpose(G),Mat), G );
  genTerm<genTensor0<dfloat<double> >,2>::Handle Term2 = Compose< FullContrProdOp >(G, G);
  genTerm<genTensor0<dfloat<double> >,2>::Handle Termadd = Term+Term2;
  genTerm<genTensor0<double>,1>::Handle Termadd2 = FSpaceDisp2+FSpaceDisp3;
  genTerm<genTensor0<double>,2>::Handle Term3 = Compose<FullContrProdOp>(G2,G2);
                            
  genTerm<genTensor0<double>,2>::Handle NominalTerm = Nominal(Term);
  genTerm<genTensor0<double>,2>::Handle TangentTerm = Derive(Term);
  genTerm<genTensor0<dfloat<double> >,2>::Handle RebuildTerm=Build<genTensor0<dfloat<double> > >(NominalTerm,TangentTerm); // should be automatic
 

  

  GmshInitialize(argc, argv);
  GModel* Model=new GModel();
  Model->readMSH("mesh.msh");
  int dim = Model->getNumRegions() ? 3 : 2;
  std::map<int, std::vector<GEntity*> > groups[4];
  GModel::current()->getPhysicalGroups(groups);
 
  std::cout << "There are " << groups[0].size() << " 0D groups, "  << groups[1].size() << " 1D groups, " << groups[2].size() << " 2D groups and " << groups[3].size() << " 3D groups" << std::endl;
  for (int i=0;i<4;++i)
  {
    if (groups[i].size())
    {
      std::cout << i << "D groups : ";
      for (std::map<int, std::vector<GEntity*> >::iterator j=groups[i].begin();j!=groups[i].end();j++)
        std::cout << j->first << " " ;
      std::cout << std::endl;
    }
  }
  genGroupOfElements* g = new genGroupOfElements(dim, 7);
  linearSystemCSRGmm<double>* lsys= new linearSystemCSRGmm<double>;
  dofManager<double>* pAssembler = new dofManager<double> ( lsys );
  NumberDofs ( FSpaceDisp, g->begin(), g->end(),*pAssembler );
  NumberDofs ( FSpaceDisp2, g->begin(), g->end(),*pAssembler );
  cout << " numberdofs" << endl;
  

  genMatrix<genTensor0<dfloat<double> > > M;
  genMatrix<genTensor0<dfloat<double> > > MM;
  genMatrix<genTensor0<double> > M2;
  genMatrix<genTensor0<double> > M3;

  genMatrix<genTensor0<double> > MN;
  genMatrix<genTensor0<double> > MT;
  genMatrix<genTensor0<dfloat<double> > > MR;
  
  IntPt pt;pt.pt[0]=0.333;pt.pt[1]=0.333;pt.pt[2]=0;pt.weight=0.5;

  Term->get(*(g->begin()),1,&pt,M);
  NominalTerm->get(*(g->begin()),1,&pt,MN);
  TangentTerm->get(*(g->begin()),1,&pt,MT);
//  RebuildTerm->get(*(g->begin()),1,&pt,MR);
  for (int i=0;i<M.size1();++i)
  {
    for (int j=0;j<M.size2();++j)
      cout << M(i,j) << " " ;
    cout << endl;
  }

  MN.print();
  MT.print();
  for (int i=0;i<MR.size1();++i)
  {
    for (int j=0;j<MR.size2();++j)
      cout << MR(i,j) << " " ;
    cout << endl;
  }

      for (int i=0;i<M.size1();++i)
  {
    for (int j=0;j<M.size2();++j)
      cout << M(i,j) << " " ;
    cout << endl;
  }

  
  Term2->get(*(g->begin()),1,&pt,MM);
  TangentTerm->get(*(g->begin()),1,&pt,M2);
  Term->get(*(g->begin()),1,&pt,M);
  Term3->get(*(g->begin()),1,&pt,M2);
  for (int i=0;i<M.size1();++i)
  {
    for (int j=0;j<M.size2();++j)
      cout << M(i,j) << " " ;
    cout << endl;
  }
    cout << endl;
  for (int i=0;i<MM.size1();++i)
  {
    for (int j=0;j<MM.size2();++j)
      cout << MM(i,j) << " " ;
    cout << endl;
  }
    cout << endl;
      for (int i=0;i<M2.size1();++i)
  {
    for (int j=0;j<M2.size2();++j)
      cout << M2(i,j) << " " ;
    cout << endl;
  }  cout << endl;
  NominalTerm->get(*(g->begin()),1,&pt,M2);
  for (int i=0;i<M2.size1();++i)
  {
    for (int j=0;j<M2.size2();++j)
      cout << M2(i,j) << " " ;
    cout << endl;
  }
  genTensor1<dfloat<> > xx(dfloat<>(1.0,0.0));
  xx(0)=dfloat<>(1.0,1.0);
  dfloat<> ttt=norm(xx);
  std::cout << xx << ttt <<  " " << ttt.v*ttt.dv << std::endl;
  return 0;
#ifdef USE_FTENSOR
  Tensor1<double,3> T1(1,2,3);
  Tensor2<double,3,3> T2;
  Tensor3_dg<double ,3,3> T3;
  Tensor4_ddg<double,3,3> T4;
  Index<'i',3> i,ii;
  Index<'i',3> iii(i);
  ii=i;
  Index<'j',3> j;

  Op<Tensor1<double,3>,Tensor1<double,3>,Tensor2<double,3,3>,'j', 'i', 3> Op;

  Op(T1,T2);
  //T2(i,j)=T1(i)*T1(j);
  //T4(i,j,k,l)=T2(i,j)*T2(k,l);

  for(int i=0;i<3;++i)
  {
  for(int j=0;j<3;++j)
          std::cout << T2(i,j) << " ";
          std::cout << std::endl;
  }
  genFSpace<Scalar0<double,3> >::Handle FSpaceDispFT(new ScalarLagrangeFSpace<Scalar0<double,3> >(126));
//  genTerm<Tensor1<double,3>,1>::Handle graddisp=Gradient(FSpaceDispFT);
  Tensor1<dfloat<double>,3> ntval(dfloat<double>(465,0.5),dfloat<double>(1,1),dfloat<double>(465,0));
  genTerm<Tensor1<dfloat<double>,3>,0>::Handle ntMat(new ConstantField<Tensor1<dfloat<double>,3> >(ntval));
  genTerm<Tensor1<dfloat<double>,3>,0>::Handle kkk= ntMat + ntMat;
  genScalar<Tensor1<dfloat<double>,3> > Res;
  kkk->get(*(g->begin()),1,&pt,Res);
  std::cout << Res()(0) << std::endl;
#endif // USE_FTENSOR
#endif
}