// 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 "elastodynamicSolver.h"
#include "genSimpleFunctionCase1.h"
#include "genAlgorithms.h"
#include "genSField.h"
#include <stdio.h>
#include "linearSystemCSR.h"
#include "linearSystemPETSc.h"
#include "linearSystemId.h"
#include "lectureData.h"
#include "genFunction.h"
#include "polarVector.h"
#include "genSimpleFunctionUinc.h"


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

ElastodynamicSolver::~ElastodynamicSolver()
{
  if (pAssembler) delete pAssembler;
  for(int i=0; i< ElastodynamicDomains.size(); i++)
    delete ElastodynamicDomains[i];
    
}

void ElastodynamicSolver::CheckProperties()
{
  for (int i=0; i<Data->Domains.size(); ++i)
  {
    if (Data->Domains[i]->type=="Elastodynamic")
    {
      ElasticDomain *field = new ElasticDomain(*(Data->Domains[i]));
      ElastodynamicDomain *field2 = new ElastodynamicDomain(*field,*(Data->Domains[i]));
      ElastodynamicDomains.push_back (field2);
    }
    else if (Data->Domains[i]->type=="IsotropicElastodynamic")
    {
      IsotropicElasticDomain *field = new IsotropicElasticDomain(*(Data->Domains[i]));
      ElastodynamicDomain *field2 = new ElastodynamicDomain(*field,*(Data->Domains[i]));
      ElastodynamicDomains.push_back (field2);
    }
    else if (Data->Domains[i]->type=="OrthotropicElastodynamic")
    {
      OrthotropicElasticDomain *field = new OrthotropicElasticDomain(*(Data->Domains[i]));
      ElastodynamicDomain *field2 = new ElastodynamicDomain(*field,*(Data->Domains[i]));
      ElastodynamicDomains.push_back (field2);
    }
  }
}


void ElastodynamicSolver::readInputFile ( const std::string &fileName )
{
  std::cout<<"model name : " << fileName << std::endl ;
  Data->readInputFile(fileName);
  Data->dumpContent(std::cout);
  CheckProperties();
  printf("--> %d Domains\n", (int)ElastodynamicDomains.size());
  printf("--> %d BCs \n", (int)Data->BCs.size());
}


void ElastodynamicSolver::CreateFunctionSpaces()
{
  if ( Data->dim==3 )
    FSpaceDisp = genFSpace<genTensor1<std::complex<double> > >::Handle(new VectorLagrangeFSpace<std::complex<double> >());
  if ( Data->dim==2 )
    FSpaceDisp = genFSpace<genTensor1<std::complex<double> > >::Handle(new VectorLagrangeFSpace<std::complex<double> >(VectorLagrangeFSpace<std::complex<double> >::VECTOR_X, 
	VectorLagrangeFSpace<std::complex<double> >::VECTOR_Y));

  std::cout << "Displacement Field : iField=" << FSpaceDisp->getIncidentSpaceTag() << std::endl;

  genTerm<genTensor2<std::complex<double> > ,1>::Handle G(Gradient(FSpaceDisp)*0.5);
  genTerm<genTensor2<std::complex<double> >,1>::Handle CG(new genCache<genTensor2<std::complex<double> > ,1>(G));
  EpsilonDisp = (CG+Transpose(CG));
}


void ElastodynamicSolver::BuildFunctionSpaces()
{
  CreateFunctionSpaces();
  lectureData fichier = lectureData("coef_An.txt");
  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")
      {
        std::cout <<  "Dirichlet BC Disp" << std::endl;
        genFilterDofComponent filter ( Data->BCs[i]->comp1 );
	genSimpleFunction<genTensor0<std::complex<double> > >* u_inc = 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_inc, filter );
      }
    }
  }
  
  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 ElastodynamicSolver::AssembleRHS()
{
  GaussQuadrature Integ_Boundary ( GaussQuadrature::ValVal );

  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")
      {
        std::cout << "Sommerfeld BC" << std::endl;
	double vp = 4000;
	double vs = 2000;
	ElastodynamicDomain* E = 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);
        genTerm<genTensor1<std::complex<double> >,1>::Handle Func=Compose<FullContrProdOp>(Compose<FullContrProdOp>(FSpaceDisp,FuncComplex_n),FuncComplex_n)*vp
        + Compose<FullContrProdOp>(Compose<FullContrProdOp>(FSpaceDisp,FuncComplex_t),FuncComplex_t)*vs; 
	genTerm<genTensor0<std::complex<double> >,2>::Handle f = Compose<FullContrProdOp>(Func, Conj(FSpaceDisp))*std::complex<double>(0.0, -E->omega*E->rho);
        Assemble (f, Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_Boundary, *pAssembler );
	delete n; delete t;
      }
    }
  }
}

void ElastodynamicSolver::AssembleLHS()
{
  genTerm<genTensor0<std::complex<double> >,2>::Handle M=Compose<FullContrProdOp>(FSpaceDisp,Conj(FSpaceDisp));
  for (int i=0;i<ElastodynamicDomains.size();++i)
  {
    ElastodynamicDomain *E = ElastodynamicDomains[i];
    genTerm<genTensor0<std::complex<double> >,2>::Handle Total = E->Energy<>(EpsilonDisp,Conj(EpsilonDisp))+ M*std::complex<double>(-pow(E->omega,2)*E->rho,0.0);
    GaussQuadrature Integ_Bulk ( GaussQuadrature::GradGrad );
    Assemble ( Total, E->begin(), E->end(), Integ_Bulk, *pAssembler );
  }
}


void ElastodynamicSolver::BuildLinearSystem()
{
  AssembleRHS();
  AssembleLHS();

  printf ( "nDofs=%d\n",pAssembler->sizeOfR() );
  printf ( "-- done assembling!\n" );

  if ( Data->solvertype == 1 )
  {
    exportKbCompressed();
//     exportKbToMatrixMarketFormat();
  }
}

void ElastodynamicSolver::solve()
{
  linearSystem<std::complex<double> > * lsys=NULL;
  if ( Data->solvertype == 2 )
  {
#if defined(HAVE_TAUCS)
    lsys = new linearSystemCSRTaucs<std::complex<double> >;
#else
    printf ( "Taucs is not installed : Gmm is chosen to solve but not available with complex\n" );
//    linearSystemCSRGmm<std::complex<double> > *buf= new linearSystemCSRGmm<std::complex<double> >;
//    buf->setNoisy ( 2 );
//    lsys=buf;
//    Data->solvertype = 1;
#endif
  }
  else if ( Data->solvertype == 3 )
  {
//#if defined(HAVE_PETSC)
//    lsys = new linearSystemPETSc<std::complex<double> >;
//#else
    printf ( "Petsc is not installed : Gmm is chosen to solve but not available with complex\n" );
//    linearSystemCSRGmm<std::complex<double> > *buf= new linearSystemCSRGmm<std::complex<double> >;
//    buf->setNoisy ( 2 );
//    lsys=buf;
//    Data->solvertype = 1;
//#endif

  }
  else if ( Data->solvertype == 1 )
  {
    printf ( "Gmm is chosen to solve but not available with complex\n" );
//    linearSystemCSRGmm<std::complex<double> > *buf= new linearSystemCSRGmm<std::complex<double> >;
//    buf->setNoisy ( 2 );
//    lsys=buf;
  }

  if ( pAssembler ) delete pAssembler;
  if (lsys)
  {
  pAssembler = new dofManager<std::complex<double> > ( lsys );

  printf ( "-- start solving\n" );
  BuildFunctionSpaces();
  BuildLinearSystem();
  printf("ElastodynamicDomains size %d \n", (int)ElastodynamicDomains.size());
  pAssembler->systemSolve();
  printf ( "-- done solving! \n" );
  }
  else
  {
   printf ( "no solver compatible with complex\n" );
   }
}

void ElastodynamicSolver::exportKb()
{
  FILE *f = fopen("K.txt", "w");
  std::complex<double> valeur;
  std::string sysname = "A";
  for (int i=0; i < pAssembler->sizeOfR(); ++i)
  {
    for (int j=0; j < pAssembler->sizeOfR(); ++j)
    {
      pAssembler->getLinearSystem(sysname)->getFromMatrix(i,j,valeur);
      if (valeur.real() || valeur.imag() )
      fprintf ( f,"(%+.16e, %+.16e) ",valeur.real(),valeur.imag() );
    }
    fprintf(f ,"\n");
  }
  fclose(f);

  f = fopen("b.txt", "w");
  for ( int i=0 ; i < pAssembler->sizeOfR(); ++i)
  {
    pAssembler->getLinearSystem(sysname)->getFromRightHandSide(i,valeur);
    fprintf ( f,"(%+.16e, %+.16e) ",valeur.real(),valeur.imag() );
  }
  fclose(f);
}

void ElastodynamicSolver::exportKbCompressed()
{
  FILE *f = fopen("K.txt", "w");
  std::complex<double> valeur;
  std::string sysname = "A";
  for (int i=0; i < pAssembler->sizeOfR(); ++i)
  {
    bool first=true;
    int nz=0;
    for (int j=0; j < pAssembler->sizeOfR(); ++j)
    {
      pAssembler->getLinearSystem(sysname)->getFromMatrix(i,j,valeur);
      if (valeur!=0.)
      {
        if (nz>3) fprintf(f, "Z%d ", nz);
        else
          for (int k=0;k<nz;++k) fprintf(f, "%d ", 0);
        nz=0;
        fprintf ( f,"(%+.16e, %+.16e) ",valeur.real(),valeur.imag() );
      }
      else
        nz++;
    }
    if (nz>3) fprintf(f, "Z%d ", nz);
     else
      for (int k=0;k<nz;++k) fprintf(f, "%d ", 0);
    fprintf(f ,"\n");
  }
  fclose(f);

  f = fopen("b.txt", "w");
  for ( int i=0 ; i < pAssembler->sizeOfR(); ++i)
  {
    pAssembler->getLinearSystem(sysname)->getFromRightHandSide(i,valeur);
    fprintf ( f,"(%+.16e, %+.16e) ",valeur.real(),valeur.imag() );
  }
  fclose(f);
}

void ElastodynamicSolver::exportb()
{
  FILE *f;
  std::complex<double> valeur;
  std::string sysname = "A";
  f = fopen ( "b.txt", "w" );
  fprintf(f,"# name: b\n");
  fprintf(f,"# type: complex matrix\n");
  fprintf(f,"# rows: %d\n",pAssembler->sizeOfR());
  fprintf(f,"# columns: 1\n");
  for ( int i = 0 ; i < pAssembler->sizeOfR() ; i ++ )
  {
    pAssembler->getLinearSystem ( sysname )->getFromRightHandSide ( i,valeur );
    fprintf ( f,"(%+.16e, %+.16e)\n",valeur.real(),valeur.imag() ) ;
  }
  fclose ( f );
}

void ElastodynamicSolver::exportx()
{
  FILE *f;
  std::complex<double> valeur;
  std::string sysname = "A";
  f = fopen ( "x.txt", "w" );
  fprintf(f,"# name: x\n");
  fprintf(f,"# type: complex matrix\n");
  fprintf(f,"# rows: %d\n",pAssembler->sizeOfR());
  fprintf(f,"# columns: 1\n");
  for ( int i = 0 ; i < pAssembler->sizeOfR() ; i ++ )
  {
    pAssembler->getLinearSystem ( sysname )->getFromSolution ( i,valeur );
    fprintf ( f,"(%+.16e, %+.16e)\n",valeur.real(),valeur.imag() ) ;
  }
  fclose ( f );
}

void ElastodynamicSolver::exportSolNum(){
    lectureData fichier = lectureData("coef_An.txt");
    
    genSimpleFunction<genTensor0<std::complex<double> > > *f = 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 ) );
    std::map<std::vector<double>, std::complex<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<genTensor0<std::complex<double> > > > val(1);
        Func->get(e, 1, &pt, val);
	
	std::vector<double> pnt_bis(3);
	pnt_bis[0] = pnt[0] ; pnt_bis[1] = pnt[1] ; pnt_bis[2] = pnt[2] ; 

        data.insert(std::pair<std::vector<double>, std::complex<double> >(pnt_bis, val[0]()() ) );
      }
    }
     
  FILE *fi;
  fi = fopen ( "NumSol.txt", "w" );
  fprintf(fi, " # valeur de ux \n");
      std::cout << "data size " << data.size() << std::endl;
      
  for (std::map<std::vector<double>, std::complex<double> >::iterator it = data.begin() ; it != data.end() ; it++){
    fprintf(fi,"x : %f, y : %f, z : %f, value : (%+.16e, %+.16e)\n", it->first[0] , it->first[1], it->first[2] , it->second.real(), it->second.imag() );  
  }
  
  fclose(fi);
  delete f;
}

void ElastodynamicSolver::error_norm_L2()
{
  lectureData fichier = lectureData("coef_An.txt");
  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* ElastodynamicSolver::buildViewDispInc_x( const std::string &postFileName ){
    lectureData fichier = lectureData("coef_An.txt");
    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* ElastodynamicSolver::buildViewDispAna_x( const std::string &postFileName ){
    lectureData fichier = lectureData("coef_An.txt");
    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( new genFunction<genTensor1<std::complex<double> > >( *f ) );
    genTerm<genTensor1<double>,0>::Handle Func_Real(Im(Func));
        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);
        Func_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);
    delete f;
  return pv;
}



PView* ElastodynamicSolver::buildViewDispNum_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;
    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);

        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* ElastodynamicSolver::buildViewDispDiff_x(const std::string &postFileName){
    //Ana
    lectureData fichier = lectureData("coef_An.txt");
    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;
  
}

PView* ElastodynamicSolver::buildDisplacementView(const std::string &postFileName)
{
  genTerm<genTensor1<std::complex<double> >,0>::ConstHandle Field(new genSField<genTensor1<std::complex<double> > >( pAssembler, FSpaceDisp ));
  return buildViewNodal(Field,ElastodynamicDomains,postFileName);
}


PView *ElastodynamicSolver::buildStressView(const std::string &postFileName)
{
  diffTerm<genTensor1<std::complex<double> >,0>::Handle Field(new genSField<genTensor1<std::complex<double> > >( pAssembler, FSpaceDisp ));
  genTerm<genTensor2<std::complex<double> >,0>::Handle G = Gradient(Field);
  genTerm<genTensor2<std::complex<double> >,0>::Handle CG(new genCache<genTensor2<std::complex<double> >,0>(G));
  genTerm<genTensor2<std::complex<double> >,0>::Handle B = (CG+Transpose(CG))*0.5;

  std::vector<genTerm<genTensor2<std::complex<double> >,0>::ConstHandle> Stress;
  for( unsigned int i = 0; i < ElastodynamicDomains.size(); i++ )
  {
    ElastodynamicDomain *E = ElastodynamicDomains[i];
    Stress.push_back(E->Stress<>(B));
  }

  return buildViewElementNode(Stress, ElastodynamicDomains, postFileName);
//   return buildViewElement(Stress, ElastodynamicDomains, postFileName);
}

PView *ElastodynamicSolver::buildStrainView(const std::string &postFileName)
{
  diffTerm<genTensor1<std::complex<double> >,0>::Handle Field(new genSField<genTensor1<std::complex<double> > >( pAssembler, FSpaceDisp ));
  genTerm<genTensor2<std::complex<double> >,0>::Handle G = Gradient(Field);
  genTerm<genTensor2<std::complex<double> >,0>::Handle CG(new genCache<genTensor2<std::complex<double> >,0>(G));
  genTerm<genTensor2<std::complex<double> >,0>::ConstHandle Strain = (CG+Transpose(CG))*0.5;

  return buildViewElementNode(Strain, ElastodynamicDomains, postFileName);
//   return buildViewElement(*Strain, ElastodynamicDomains, postFileName);
}