// GenFem - A high-level finite element library
// 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>.

// check terms
#include "GmshGlobal.h"
#include "GModel.h"
#include "linearSystemCSR.h"
#include "genTerm.h"
#include "genFSpace.h"
#include "genAlgorithms.h"
#include "genTermCompose.h"


int main(int argc,char**argv)
{
  const double mat[]={1.0,0.0,0.0,
                     0.0,2.0,0.0,
                     0.0,0.0,3.0};
//  mat(0,0)=2.0;mat(0,1)=0.0;mat(0,2)=0.0
  diffTerm<genTensor<0,double>,1>::Handle FSpace(new ScalarLagrangeFSpace<genTensor<0,double> >(126));
  genTerm<genTensor<1,double>,1>::Handle Grad = Gradient(FSpace);
  genTerm<genTensor<2,double>,0>::Handle Mat(new ConstantField<genTensor<2,double> >(mat));
  genTerm<genTensor<0,double>,2>::Handle Term = Compose<FullContrProdOp>(Compose<FullContrProdOp>( Transpose(Grad),Mat), Grad);
//  genTerm<genTensor<0,double>,2>::Handle Term_new = Compose<NewContrProdOp<genTerm<genTensor<1,double>,1>,genTerm<genTensor<1,double>,1>,0> >(Compose<FullContrProdOp>( Transpose(Grad),Mat), Grad);

  genTerm<genTensor<0,double>,2>::Handle Term2 = Compose< FullContrProdOp >(Grad, Grad);
  /*
  NewContrProdOp<genTensor<1,double>,genTensor<1,double>,genTensor<0,double> > Op;
  FullContrProdOp<genTensor<1,double>,genTensor<1,double> > Op2;
  genTerm<genTensor<0,double>,2>::Handle Term2new;
  Term2new = Compose3<NewContrProdOp<genTensor<1,double>,genTensor<1,double>,genTensor<0,double> >, genTerm<genTensor<1,double>,1>,genTerm<genTensor<1,double>,1> >(Grad, Grad,Op);
*/
  /*  genTerm<genTensor<0,double>,2>::Handle Term2new = 
  Compose< 
  NewContrProdOp<genTerm<genTensor<1,double>,1>,
                 genTerm<genTensor<1,double>,1>,
                 genTerm<genTensor<0,double>,1> > >(Grad, Grad);
*/
  GmshInitialize(argc, argv);
  GModel* Model=new GModel();
  std::cout << 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 ( FSpace, g->begin(), g->end(),*pAssembler );
  std::cout << " numberdofs" << std::endl;

  genMatrix<genTensor<0,double > > M;
  genMatrix<genTensor<0,double> > M2;

  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);

  std::cout.precision(15);
  std::cout << std::showpos;
  std::cout << std::scientific;

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


  Term2->get(*(g->begin()),1,&pt,M2);
  std::cout << std::endl;
  for (int i=0;i<M2.size1();++i)
  {
    for (int j=0;j<M2.size2();++j)
      std::cout << M2(i,j)() << " " ;
    std::cout << std::endl;
  }
  std::cout << std::endl;

  std::cout << M.norm()() << std::endl;
  std::cout << M2.norm()() << std::endl;
  if (fabs( M.norm()() - 2.121320343559643e+00) > 1e-12) return 1;
  if (fabs( M2.norm()() - 1.581138830084189e+00) > 1e-12) return 1;
  
  return 0; 
}