// Scalar_Helmholtz_xfem - A linear solver for the scalar helmholtz equation using X-FEM
// Copyright (C) 2012-2026 Eric Bechet
//
// See the LICENSE file for license information.
// Please report all bugs and problems to <bechet@cadxfem.org>.

#include <cstring>
#include <cstdio>
#include <fstream>
#include <complex>

#include "helmholtz_xfem.h"
#include "ConvertEnrichedFSpace.h"
#include "RestrictedFilteredFSpace.h"
#include "ConvertEnrichedFSpace.h"
#include "genAlgorithms.h"
#include "genSField.h"
#include "genFunction.h"
#include "linearSystemCSR.h"
#include "linearSystemPETSc.h"
#include "linearSystemPETSc.hpp"
#include "linearSystemId.h"

#if defined(HAVE_POST)
#include "PView.h"
#include "PViewData.h"
#endif


HelmholtzLevelSetSolver::~HelmholtzLevelSetSolver()
{
  for (int i=0;i<HLSdomains.size();i++) delete HLSdomains[i];
  //if(DataHLS) delete DataHLS;
}


void HelmholtzLevelSetSolver::readInputFile ( const std::string &fileName )
{
  std::cout<<"model name : " << fileName << std::endl ;
  DataHLS->readInputFile(fileName);
  DataHLS->dumpContent(std::cout);
  //*Data=*DataHLS;
  CheckProperties();
  printf("--> %d Supports\n", (int)IDomains.size() );
  printf("--> %d BCs \n", (int)DataHLS->BCs.size());
}


void HelmholtzLevelSetSolver::CheckProperties()
{
  ScalarHelmholtzSolver::CheckProperties();

  for (int i=0; i<DataHLS->LSs.size(); ++i)
  {
    LevelSetDomain* field = new LevelSetDomain(*(DataHLS->LSs[i]));
    LSDomains.push_back(field);
  }

  for (int i=0; i<DataHLS->Domains.size(); ++i)
  {
    if (DataHLS->Domains[i]->type=="Bimaterial")
    {
      int dim = DataHLS->Domains[i]->dim;
      int tag = DataHLS->Domains[i]->tag;
      LevelSetDomain* LS=NULL;

      for (int j=0; j<LSDomains.size(); ++j)
      {
        if ((LSDomains[j]->dim==dim)&&(LSDomains[j]->tag==tag))
          LS=LSDomains[j];
      }

      if (LS)
      {
        BimaterialDomain* field = new BimaterialDomain(*(DataHLS->Domains[i]), LS->Distance());
        BiDomains.push_back(field);
      }
      else
        printf("LevelSetDomain with dim=%d and tag=%d is not already existing.\n", dim, tag);
    }
  }

  for (int i=0; i<DataHLS->HLSsolutions.size(); ++i)
  {
    for (int j=0; j<DataHLS->Domains.size(); j++)
    {
      if (DataHLS->Domains[j]->dim==DataHLS->HLSsolutions[i]->dim
          &&
          DataHLS->Domains[j]->tag==DataHLS->HLSsolutions[i]->tag)
      {
        HelmholtzLevelsetDomain* field = new HelmholtzLevelsetDomain(*(DataHLS->HLSsolutions[i]),*(DataHLS->Domains[j]));
        HLSdomains.push_back(field);
      }
    }
  }
}


void HelmholtzLevelSetSolver::CreateFunctionSpace()
{
  ScalarHelmholtzSolver::CreateFunctionSpace();
  std::cout << "Using_xfem " << DataHLS->Using_xfem << std::endl;
  if(DataHLS->Using_xfem)
  {
    for (int i=0; i<BiDomains.size(); ++i)
    {
      BimaterialDomain* B = BiDomains[i];
      std::cout << "Domain of bi-material interface " << i << " : dim=" << B->dim << ", tag=" << B->tag << ", " << B->group()->psize() << " parent elements" << std::endl;

      std::cout << B->group()->vsize() << " parent vertices" << std::endl;
      std::cout << B->group()->vsize() - B->EnrichedEntities()->size() << " parent vertices along the interface" << std::endl;
      std::cout << B->EnrichedEntities()->size() << " enriched vertices" << std::endl;

      diffTerm<genTensor0<std::complex<double> >,0>::ConstHandle iEF(new ConvertEnrichedFSpace<genTensor0<std::complex<double> > , genTensor0<double> >(B->EnrichmentFunction()));

      genFilterDofMultiEntities* filterDof = new genFilterDofMultiEntities(B->EnrichedEntities());
      genFSpace<genTensor0<std::complex<double> > >::ConstHandle LagrangeFiltered (new FilteredFSpace<genTensor0<std::complex<double> > >(FSpace,filterDof));
      genFSpace<genTensor0<std::complex<double> > >::ConstHandle xFemLagrange (new EnrichedFSpace<genTensor0<std::complex<double> > >(LagrangeFiltered, iEF));

//       genFilterElementGroup* FilterInt = new genFilterElementGroup(B->subDomainInt);
//       genFilterElementGroup* FilterExt = new genFilterElementGroup(B->subDomainExt);
//       genFilterElementOperatorOr* FilterRestricted = new genFilterElementOperatorOr(FilterInt, FilterExt);
//       genFSpace<genTensor0<std::complex<double> > >::ConstHandle restrictedLagrangeFiltered (new RestrictedFilteredFSpace<genTensor0<std::complex<double> > >(FSpace, LagrangeFiltered, FilterRestricted));

      FSpace = genFSpace<genTensor0<std::complex<double> > >::Handle(new CompositeFSpace<genTensor0<std::complex<double> > >(FSpace, xFemLagrange));
    }
  }
}


void HelmholtzLevelSetSolver::error()
{
  genTensor0<std::complex<double> > sum_diff_L2;
  genTensor0<std::complex<double> > sum_ana_L2;
  genTensor0<std::complex<double> > sum_diff_H1;
  genTensor0<std::complex<double> > sum_ana_H1;
  for (int i=0; i< HLSdomains.size(); ++i)
  {
    HelmholtzLevelsetDomain* P = HLSdomains[i];
    diffTerm<genTensor0<std::complex<double> >,0>::Handle Field_sol(new genSField<genTensor0<std::complex<double> > >( pAssembler, FSpace));
    diffTerm<genTensor0<double>,0>::Handle Func_Real( new genFunction<genTensor0<double> >(*(P->fscalar_real)) );
    diffTerm<genTensor0<double>,0>::Handle Func_Im( new genFunction<genTensor0<double> >(*(P->fscalar_im)) );
    genTerm<genTensor1<double>,0>::Handle Grad_Func_Real = Gradient(Func_Real);
    genTerm<genTensor1<double>,0>::Handle Grad_Func_Im = Gradient(Func_Im);
    genTerm<genTensor1<std::complex<double> >, 0>::Handle Field_Grad_ana = Build<genTensor1<std::complex<double> > >(Grad_Func_Real)*std::complex<double>(1.0,0)
                                                                         + Build<genTensor1<std::complex<double> > >(Grad_Func_Im)*std::complex<double>(0.0,1.0);
    genTerm<genTensor0<std::complex<double> >,0>::Handle Field_ana = Build<genTensor0<std::complex<double> > >(Func_Real)*std::complex<double>(1.0,0)
                                                                    + Build<genTensor0<std::complex<double> > >(Func_Im)*std::complex<double>(0.0,1.0);

    GaussQuadrature Integ_Bulk ( GaussQuadrature::ValVal );

    // norme L2
    genTerm<genTensor0<std::complex<double > >,0>::Handle Field_diff = Field_ana - Field_sol;
    genTerm<genTensor0<std::complex<double > >,0>::Handle Err_L2=Compose<FullContrProdOp>(Field_diff,Conj(Field_diff));
    genTerm<genTensor0<std::complex<double > >,0>::Handle Ref_L2=Compose<FullContrProdOp>(Field_ana,Conj(Field_ana));
    Assemble(Err_L2,P->begin(),P->end(),Integ_Bulk,sum_diff_L2);
    Assemble(Ref_L2,P->begin(),P->end(),Integ_Bulk,sum_ana_L2);

    // norme H1
    genTerm<genTensor1<std::complex<double> >, 0>::Handle Field_Grad_diff = Field_Grad_ana - Gradient(Field_sol);
    genTerm<genTensor0<std::complex<double> >, 0>::Handle Err_H1 = Compose<FullContrProdOp>(Field_Grad_diff,Conj(Field_Grad_diff));
    //- (P->k)*(P->k)*Err_L2;
    genTerm<genTensor0<std::complex<double> >, 0>::Handle Ref_H1 = Compose<FullContrProdOp>(Field_Grad_ana,Conj(Field_Grad_ana)); 
    //- (P->k)*(P->k)*Ref_L2;
    Assemble(Err_H1,P->begin(),P->end(),Integ_Bulk,sum_diff_H1);
    Assemble(Ref_H1,P->begin(),P->end(),Integ_Bulk,sum_ana_H1);

    Field_ANA.push_back(Re(Field_ana));
    Field_DIFF.push_back(Re(Field_diff));
  }
  std::cout << "sum_ana_L2 " << sum_ana_L2() << "\n";
  std::cout << "sum_diff_L2 " << sum_diff_L2() << "\n";
  std::complex<double> error_L2=sum_diff_L2()/sum_ana_L2();
  std::cout << "displacement error " << pow(std::real(error_L2),0.5) << std::endl;

  std::cout << "sum_ana_H1 " << sum_ana_H1() << "\n";
  std::cout << "sum_diff_H1 " << sum_diff_H1() << "\n";
  std::complex<double> error_H1=sum_diff_H1()/sum_ana_H1();
  std::cout << "energy error " << pow(std::real(error_H1), 0.5) << std::endl;
}


PView* HelmholtzLevelSetSolver::buildViewDiff( const std::string &postFileName )
{
  return buildViewNodal(Field_DIFF,HLSdomains,postFileName);
}

PView* HelmholtzLevelSetSolver::buildViewAna( const std::string &postFileName )
{
  return buildViewNodal(Field_ANA,HLSdomains,postFileName);
}

PView* HelmholtzLevelSetSolver::buildViewSol( const std::string &postFileName )
{
  diffTerm<genTensor0<std::complex<double> >,0>::Handle Field_sol(new genSField<genTensor0<std::complex<double> > >( pAssembler, FSpace));
  genTerm<genTensor0<double>,0>::Handle RealField(Re(Field_sol));
  return buildViewNodal(RealField,HLSdomains,postFileName);
}