// elastic_genTerm - A linear solver for elastic problems using FEM
// 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>.

#include "elastodynamicSolverDisk.h"
#include "genAlgorithms.h"
#include "genSField.h"
#include <stdio.h>
#include "linearSystemCSR.h"
#include "linearSystemPETSc.h"
#include "linearSystemId.h"
#include "genFunction.h"
#include "polarVector.h"
#include "genSimpleFunctionUinc.h"
#include "genSimpleFunctionUtr.h"
#include "genSimpleFunctionUtot.h"
#include "genSimpleFunctionSigmaInc.h"
#include "genSimpleFunctionCase1.h"


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


void ElastodynamicSolverDisk::readInputFile(const std::string &fileName){
    ElastodynamicSolver::readInputFile(fileName);
    fichier = new lectureData("coef_An_disk.txt");
}

void ElastodynamicSolverDisk::BuildFunctionSpaces()
{
  CreateFunctionSpaces(); 
  for ( unsigned int i = 0; i < Data->BCs.size(); i++ )
  {
    if ((Data->BCs[i]->kind=="Dirichlet")||(Data->BCs[i]->kind==""))
    {
      if (Data->BCs[i]->What=="Displacement_Utr")
      {
        std::cout <<  "Dirichlet BC Disp" << std::endl;
        genFilterDofComponent filter ( Data->BCs[i]->comp1 );
	genSimpleFunction<genTensor0<std::complex<double> > >* u_tr = 
	     new genSimpleFunctionUtr<genTensor0<std::complex<double> >, genTensor1<double> >(fichier->vector, fichier->kp_int, fichier->ks_int, Data->BCs[i]->comp1);
        FixNodalDofs ( FSpaceDisp, Data->BCs[i]->begin(), Data->BCs[i]->end(), *pAssembler, *u_tr, filter );
	delete u_tr;
      }
      if (Data->BCs[i]->What=="Displacement_Utot")
      {
        std::cout <<  "Dirichlet BC Disp" << std::endl;
        genFilterDofComponent filter ( Data->BCs[i]->comp1 );
	genSimpleFunction<genTensor1<std::complex<double> > >* u_ref = new genSimpleFunctionCase1<genTensor1<std::complex<double> >, genTensor1<double> >(fichier->vector, fichier->kp,
	fichier->ks);
	size_t inc_size = 100;
	genSimpleFunction<genTensor1<std::complex<double> > >* u_inc = new genSimpleFunctionUinc<genTensor1<std::complex<double> >, genTensor1<double> >(inc_size,fichier->kp);
	genSimpleFunction<genTensor0<std::complex<double> > >* u_tot_ext = new genSimpleFunctionUtot<genTensor0<std::complex<double> > >(u_inc, u_ref, Data->BCs[i]->comp1  );
	FixNodalDofs ( FSpaceDisp, Data->BCs[i]->begin(), Data->BCs[i]->end(), *pAssembler, *u_tot_ext, filter );
	delete u_tot_ext;
	delete u_ref;
	delete u_inc;
      }
      if (Data->BCs[i]->What=="Displacement_Uref")
      {
        std::cout <<  "Dirichlet BC Disp" << std::endl;
        genFilterDofComponent filter ( Data->BCs[i]->comp1 );
        genSimpleFunction<genTensor0<std::complex<double> > >* u_ref = new genSimpleFunctionCase1<genTensor0<std::complex<double> >, genTensor1<double> >(fichier->vector, fichier->kp,
	fichier->ks, Data->BCs[i]->comp1);
	FixNodalDofs ( FSpaceDisp, Data->BCs[i]->begin(), Data->BCs[i]->end(), *pAssembler, *u_ref, filter );
	delete u_ref;
      }
    }
  }
  
  for ( unsigned int i = 0; i < ElastodynamicDomains.size(); ++i )
  {
    ElastodynamicDomain *E = ElastodynamicDomains[i];
    std::cout << "Domain " << i << " : dim=" << E->dim << ", tag=" << E->tag << ", " << E->group()->size() << " elements" << std::endl;
    NumberDofs ( FSpaceDisp, E->begin(), E->end(), *pAssembler );
  }
}

void ElastodynamicSolverDisk::AssembleRHS()
{
  ElastodynamicSolver::AssembleRHS();
  for ( unsigned int i = 0; i < Data->BCs.size(); i++ )
  {
    if ((Data->BCs[i]->kind=="Neumann")||(Data->BCs[i]->kind==""))
    {
      if (Data->BCs[i]->What=="Sommerfeld_tot")
      {
        std::cout << "Sommerfeld BC" << std::endl;
	double vp = 4000;
	double vs = 2000;
	int inc_size = 100;
	ElastodynamicDomain* P = ElastodynamicDomains[0];
	genSimpleFunction<genTensor1<double> >* n = new radialVector<genTensor1<double> >();
	genSimpleFunction<genTensor1<double> >* t = new tangentVector<genTensor1<double> >();
	genTerm<genTensor1<double>,0>::Handle Func_n( new genFunction<genTensor1<double> >( *n ) );
	genTerm<genTensor1<double>,0>::Handle Func_t( new genFunction<genTensor1<double> >( *t ) );
        genTerm<genTensor1<std::complex<double> >,0>::Handle FuncComplex_n=Build<genTensor1<std::complex<double> > >(Func_n)*std::complex<double>(1.0,0.0);
	genTerm<genTensor1<std::complex<double> >,0>::Handle FuncComplex_t=Build<genTensor1<std::complex<double> > >(Func_t)*std::complex<double>(1.0,0.0);
	
	GaussQuadrature Integ_Boundary_2 ( GaussQuadrature::Val );
	genSimpleFunction<genTensor1<std::complex<double> > >* u_inc = 
	    new genSimpleFunctionUinc<genTensor1<std::complex<double> >, genTensor1<double> >(inc_size,fichier->kp);
	genTerm<genTensor1<std::complex<double> >,0>::Handle FuncInc( new genFunction<genTensor1<std::complex<double> > >( *u_inc) );
	genTerm<genTensor1<std::complex<double> >,0>::Handle FuncInc_Neum=(Compose<FullContrProdOp>(Compose<FullContrProdOp>(FuncInc,FuncComplex_n),FuncComplex_n)*vp
	    + Compose<FullContrProdOp>(Compose<FullContrProdOp>(FuncInc,FuncComplex_t),FuncComplex_t)*vs)*std::complex<double>(0.0, -P->omega*P->rho); 
	const double E = P->E;
	const double nu = P->nu;
	double lambda = E*nu/((1+nu)*(1-2*nu));
	double mu= E/(2*(1+nu));
	genSimpleFunction<genTensor1<std::complex<double> > >* sigma_inc =
	      new genSimpleFunctionSigmaInc<genTensor1<std::complex<double> >, genTensor1<double> >(inc_size,fichier->kp, lambda, mu);
	genTerm<genTensor1<std::complex<double> >,0>::Handle Sigma_Inc( new genFunction<genTensor1<std::complex<double> > >( *sigma_inc) );  
	genTerm<genTensor0<std::complex<double> >,1>::Handle FuncTot_Neum = Compose<FullContrProdOp>(Sigma_Inc+FuncInc_Neum, Conj(FSpaceDisp));
        Assemble (FuncTot_Neum, Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_Boundary_2, *pAssembler );
	
	delete n;
	delete t;
	delete u_inc;
	delete sigma_inc;
      }
    }
  }
}

void ElastodynamicSolverDisk::error_norm_L2()
{
  genTensor0<std::complex<double> > sum_diff_L2;
  genTensor0<std::complex<double> > sum_ana_L2;
  for (int i=0;i< ElastodynamicDomains.size();++i)
  {
      ElastodynamicDomain *P = ElastodynamicDomains[i];
      genTerm<genTensor1<std::complex<double> >,0>::Handle Field_num(new genSField<genTensor1<std::complex<double> > >( pAssembler, FSpaceDisp));
      genSimpleFunction<genTensor1<std::complex<double> > >* f = new genSimpleFunctionCase1<genTensor1<std::complex<double> >, genTensor1<double> >(fichier->vector, 
	fichier->kp, fichier->ks);
      genTerm<genTensor1<std::complex<double> >,0>::Handle Func_ana( new genFunction<genTensor1<std::complex<double> > >( *f ) );
      GaussQuadrature Integ_Bulk ( GaussQuadrature::ValVal );
      
      
      // norme L2
      genTerm<genTensor1<std::complex<double > >,0>::Handle Field_diff = Func_ana - Field_num;
      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>(Func_ana,Conj(Func_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);
      
      delete f;
  }
  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 << "displacement error imag part :" << pow(std::imag(error_L2),0.5) << std::endl;
  FILE *fi;
  fi = fopen ( "error_L2.txt", "at" );
  if(!fi) fi=fopen ( "error_L2.txt", "w" );
  fprintf(fi, " # displacement error \n");
  fprintf(fi, " partie reelle de l'erreur : %f \n", pow(std::real(error_L2),0.5) );
  fprintf(fi, " partie imaginaire de l'erreur : %f \n", pow(std::imag(error_L2),0.5) );
  fclose(fi);
}

PView* ElastodynamicSolverDisk::buildViewDispInc_x( const std::string &postFileName ){
    size_t v = 100;
    int comp =  0;
    genSimpleFunction<genTensor0<std::complex<double> > > *f_inc = new
    genSimpleFunctionUinc<genTensor0<std::complex<double> >, genTensor1<double> >(v,fichier->kp, comp);
    genTerm<genTensor0<std::complex<double> >,0>::Handle Func( new genFunction<genTensor0<std::complex<double> > >( *f_inc ) );
    genTerm<genTensor0<double>,0>::Handle Func_Real(Re(Func));
    return buildViewNodal(Func_Real,ElastodynamicDomains,postFileName);
}


PView* ElastodynamicSolverDisk::buildViewDispAna_x( const std::string &postFileName ){
    genSimpleFunction<genTensor1<std::complex<double> > >* u_ref = new genSimpleFunctionCase1<genTensor1<std::complex<double> >, genTensor1<double> >(fichier->vector, fichier->kp, fichier->ks);
    size_t v = 100;
    genSimpleFunction<genTensor1<std::complex<double> > >* u_inc = 
	    new genSimpleFunctionUinc<genTensor1<std::complex<double> >, genTensor1<double> >(v,fichier->kp);
    genSimpleFunction<genTensor1<std::complex<double> > >* u_tot_ext = new genSimpleFunctionUtot<genTensor1<std::complex<double> > >(u_inc, u_ref );
    genTerm<genTensor1<std::complex<double> >,0>::Handle FuncExt( new genFunction<genTensor1<std::complex<double> > >( *u_tot_ext ) );
    genTerm<genTensor1<double>,0>::Handle FuncExt_Real(Re(FuncExt));
    FuncANA.push_back(FuncExt_Real);
    genSimpleFunction<genTensor1<std::complex<double> > >* f_int = new genSimpleFunctionUtr<genTensor1<std::complex<double> >, genTensor1<double> >(fichier->vector, fichier->kp_int, fichier->ks_int);
    genTerm<genTensor1<std::complex<double> >,0>::Handle FuncInt( new genFunction<genTensor1<std::complex<double> > >( *f_int ) );
    genTerm<genTensor1<double>,0>::Handle FuncInt_Real(Re(FuncInt));
    FuncANA.push_back(FuncInt_Real);
    std::map<int, std::vector<double> > data;
    for(int i=0; i<ElastodynamicDomains.size();  i++){
      ElastodynamicDomain* E = ElastodynamicDomains[i]; 
      genTerm<genTensor1<double>,0>::Handle Func = FuncANA[i];
      for ( genGroupOfElements::elementContainer::const_iterator it = E->begin(); it != E->end(); ++it )
      {
	MElement* e = *it;

	for (int j=0; j<e->getNumVertices(); ++j)
	{
	  double pnt[3];
	  pnt[0]=e->getVertex(j)->x();
	  pnt[1]=e->getVertex(j)->y();
	  pnt[2]=e->getVertex(j)->z();
	  

	  double uvw[3];
	  e->xyz2uvw(pnt, uvw);

	  
	  IntPt pt;
	  pt.pt[0]=uvw[0];
	  pt.pt[1]=uvw[1];
	  pt.pt[2]=uvw[2];
	  pt.weight=1.;

	  std::vector<genScalar<genTensor1<double> > > val(1);
	  Func->get(e, 1, &pt, val);

	  int vNum = e->getVertex(j)->getNum();
	  genTensor0<double> val_bis = val[0]()(0);
	  data[vNum] = vectorize(val_bis);
	}
      }
    }
    int numComp = vectorize(genTensor0<double>()).size();
    PView* pv = new PView (postFileName, "NodeData", Data->Model, data, 0., numComp);
    delete u_ref;
    delete u_inc;
    delete u_tot_ext;
    delete f_int;
  return pv;
}



PView* ElastodynamicSolverDisk::buildDispNum_x(const std::string &postFileName){
    genTerm<genTensor1<std::complex<double> >,0>::ConstHandle Field(new genSField<genTensor1<std::complex<double> > >( pAssembler, FSpaceDisp ));
    genTerm<genTensor1<double>,0>::Handle Field_Real(Re(Field));
    std::map<int, std::vector<double> > data;
    for(int i=0; i<ElastodynamicDomains.size();  i++){
      ElastodynamicDomain* E = ElastodynamicDomains[i]; 
      for ( genGroupOfElements::elementContainer::const_iterator it = E->begin(); it != E->end(); ++it )
      {
	MElement* e = *it;

	for (int j=0; j<e->getNumVertices(); ++j)
	{
	  double pnt[3];
	  pnt[0]=e->getVertex(j)->x();
	  pnt[1]=e->getVertex(j)->y();
	  pnt[2]=e->getVertex(j)->z();
	  

	  double uvw[3];
	  e->xyz2uvw(pnt, uvw);

	  
	  IntPt pt;
	  pt.pt[0]=uvw[0];
	  pt.pt[1]=uvw[1];
	  pt.pt[2]=uvw[2];
	  pt.weight=1.;

	  std::vector<genScalar<genTensor1<double> > > val(1);
	  Field_Real->get(e, 1, &pt, val);

	  int vNum = e->getVertex(j)->getNum();
	  genTensor0<double> val_bis = val[0]()(0);
	  data[vNum] = vectorize(val_bis);
	}
      }
    }
    int numComp = vectorize(genTensor0<double>()).size();
    PView* pv = new PView (postFileName, "NodeData", Data->Model, data, 0., numComp);

  return pv;
  
}

PView* ElastodynamicSolverDisk::buildViewDispDiff_x(const std::string &postFileName){
    //Ana
    genSimpleFunction<genTensor0<std::complex<double> > >* f_x = new genSimpleFunctionCase1<genTensor0<std::complex<double> >, genTensor1<double> >(fichier->vector, 
	fichier->kp, fichier->ks);
    genTerm<genTensor0<std::complex<double> >,0>::Handle Func( new genFunction<genTensor0<std::complex<double> > >( *f_x ) );
    genTerm<genTensor0<double>,0>::Handle Func_Real(Im(Func));
    //Num
    genTerm<genTensor1<std::complex<double> >,0>::ConstHandle Field(new genSField<genTensor1<std::complex<double> > >( pAssembler, FSpaceDisp ));
    genTerm<genTensor1<double>,0>::Handle Field_Real(Im(Field));
    
    std::map<int, std::vector<double> > data;
    ElastodynamicDomain* E = ElastodynamicDomains[0]; 
    for ( genGroupOfElements::elementContainer::const_iterator it = E->begin(); it != E->end(); ++it )
    {
      MElement* e = *it;

      for (int j=0; j<e->getNumVertices(); ++j)
      {
        double pnt[3];
        pnt[0]=e->getVertex(j)->x();
        pnt[1]=e->getVertex(j)->y();
        pnt[2]=e->getVertex(j)->z();
	

        double uvw[3];
        e->xyz2uvw(pnt, uvw);

	
        IntPt pt;
        pt.pt[0]=uvw[0];
        pt.pt[1]=uvw[1];
        pt.pt[2]=uvw[2];
        pt.weight=1.;

        std::vector<genScalar<genTensor1<double> > > val(1);
        Field_Real->get(e, 1, &pt, val);
	
	std::vector<genScalar<genTensor0<double> > > val2(1);
	Func_Real->get(e, 1, &pt, val2);
	

        int vNum = e->getVertex(j)->getNum();
	genTensor0<double> val_bis = val[0]()(0);
	genTensor0<double> val2_bis  = val2[0]();
	
	genTensor0<double> val_diff = val_bis - val2_bis;
	data[vNum] = vectorize(val_diff);
      }
    }
    int numComp = vectorize(genTensor0<double>()).size();
    PView* pv = new PView (postFileName, "NodeData", Data->Model, data, 0., numComp);
    delete f_x;
  return pv;
  
}