// 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 "genTerm.h"
#include "genFSpace.h"
#include "genAlgorithms.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<genTensor0<double>,1>::Handle FSpace(new ScalarLagrangeFSpace<genTensor0<double> >(126));
  genTerm<genTensor1<double>,1>::Handle Grad = Gradient(FSpace);
  genTerm<genTensor2<double>,0>::Handle Mat(new ConstantField<genTensor2<double> >(mat));
  genTerm<genTensor0<double>,2>::Handle Term = Compose<FullContrProdOp>(Compose<FullContrProdOp>( Transpose(Grad),Mat), Grad);
  genTerm<genTensor0<double>,2>::Handle Term2 = Compose< FullContrProdOp >(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);
  linearSystemGmm<double>* lsys= new linearSystemGmm<double>;
  dofManager<double>* pAssembler = new dofManager<double> ( lsys );
  NumberDofs ( FSpace, g->begin(), g->end(),*pAssembler );
  std::cout << " numberdofs" << std::endl;

  genMatrix<genTensor0<double > > M;
  genMatrix<genTensor0<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; 
}