// 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 "elastodynamicMultLagSolver.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"
#include "polarVector.h"
#include "spaceReducer.h"

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

ElastodynamicMultLagSolver::~ElastodynamicMultLagSolver()
{
     if(fichier) delete fichier;
     for(int i= 0; i<ElastodynamicMultLagDomains.size(); i++){
       delete ElastodynamicMultLagDomains[i];
     }
     for(int i=0; i< InterfaceElastoDomains.size(); i++){
       delete InterfaceElastoDomains[i];
     }
     if(DataLS) delete DataLS;
}



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

void ElastodynamicMultLagSolver::CheckProperties()
{
   ElastodynamicSolver::CheckProperties();
   
   for(int i=0; i< DataLS->LSs.size(); i++){
     InterfaceElastoDomain* field = new InterfaceElastoDomain(*DataLS->LSs[i]); // for levelSet
     InterfaceElastoDomains.push_back(field);
   }
   
   for(int i= 0; i<Data->BCs.size(); i++){
     if (Data->BCs[i]->What=="Interface"){
       InterfaceElastoDomain* field = new InterfaceElastoDomain(*Data->BCs[i]); //No levelSet
       InterfaceElastoDomains.push_back(field);
     }
   }
   
   for(int i = 0; i< ElastodynamicDomains.size() ; i++){
     ElastodynamicMultLagDomain* field =  new ElastodynamicMultLagDomain(*ElastodynamicDomains[i]);
     ElastodynamicMultLagDomains.push_back(field);
   }
}	


void ElastodynamicMultLagSolver::CreateFunctionSpaces()
{
  FSpacesDisp.resize(ElastodynamicDomains.size());
  EpsilonSpacesDisp.resize(ElastodynamicDomains.size());
  for(int i=0; i< FSpacesDisp.size(); i++){
    if(Data->dim == 3){
        FSpacesDisp[i] = genFSpace<genTensor1<std::complex<double> > >::Handle(new VectorLagrangeFSpace<std::complex<double> >());
        genTerm<genTensor2<std::complex<double> > ,1>::Handle G(Gradient(FSpacesDisp[i])*0.5);
	genTerm<genTensor2<std::complex<double> >,1>::Handle CG(new genCache<genTensor2<std::complex<double> > ,1>(G));
	EpsilonSpacesDisp[i] = (CG+Transpose(CG));
  }
    if ( Data->dim==2 ){
	FSpacesDisp[i] = genFSpace<genTensor1<std::complex<double> > >::Handle(new VectorLagrangeFSpace<std::complex<double> >(VectorLagrangeFSpace<std::complex<double> >::VECTOR_X, 
	VectorLagrangeFSpace<std::complex<double> >::VECTOR_Y));
        genTerm<genTensor2<std::complex<double> > ,1>::Handle G(Gradient(FSpacesDisp[i])*0.5);
	genTerm<genTensor2<std::complex<double> >,1>::Handle CG(new genCache<genTensor2<std::complex<double> > ,1>(G));
	EpsilonSpacesDisp[i] = (CG+Transpose(CG));
    }
  }
  
  FSpacesMultLag.resize(InterfaceElastoDomains.size());
  for(int i=0; i< FSpacesMultLag.size(); i++){
  if(Data->dim == 3)
      FSpacesMultLag[i] = genFSpace<genTensor1<std::complex<double> > >::Handle(new VectorLagrangeFSpace<std::complex<double> >());
  if ( Data->dim==2 )
      FSpacesMultLag[i] = genFSpace<genTensor1<std::complex<double> > >::Handle(new VectorLagrangeFSpace<std::complex<double> >(VectorLagrangeFSpace<std::complex<double> >::VECTOR_X, 
      VectorLagrangeFSpace<std::complex<double> >::VECTOR_Y));
  }
  
}


void ElastodynamicMultLagSolver::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;
	for(int j =0; j< ElastodynamicMultLagDomains.size(); j++){
	  ElastodynamicMultLagDomain* P = ElastodynamicMultLagDomains[j];
	  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 ( FSpacesDisp[j], Data->BCs[i]->begin(), Data->BCs[i]->end(), *pAssembler, *u_tr, *P->FilterDof(Data->BCs[i]->comp1) );
	  delete u_tr;
	}
      }
      if (Data->BCs[i]->What=="Displacement_Utot")
      {
        std::cout <<  "Dirichlet BC Disp" << std::endl;
	for(int j =0; j< ElastodynamicDomains.size(); j++){
	  ElastodynamicMultLagDomain* P = ElastodynamicMultLagDomains[j];
	  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 ( FSpacesDisp[j], Data->BCs[i]->begin(), Data->BCs[i]->end(), *pAssembler, *u_tot_ext, *P->FilterDof(Data->BCs[i]->comp1) );
	  delete u_tot_ext;
	  delete u_ref;
	  delete u_inc;
	}
      }
      if (Data->BCs[i]->What=="Displacement_Uref")
      {
        std::cout <<  "Dirichlet BC Disp" << std::endl;
	for(int j =0; j< ElastodynamicDomains.size(); j++){
	  ElastodynamicMultLagDomain* P = ElastodynamicMultLagDomains[j];
	  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 ( FSpacesDisp[j], Data->BCs[i]->begin(), Data->BCs[i]->end(), *pAssembler, *u_ref, *P->FilterDof(Data->BCs[i]->comp1) );
	  delete u_ref;
	}
      }
    }
  }
  
  for ( unsigned int i = 0; i < ElastodynamicMultLagDomains.size(); ++i )
  {
    ElastodynamicMultLagDomain *E = ElastodynamicMultLagDomains[i];
    NumberDofs ( FSpacesDisp[i], E->begin(), E->end(), *pAssembler );
    std::cout << "Domain " << i << " : dim=" << E->dim << ", tag=" << E->tag << ", " << E->group()->size() << " elements" << std::endl;
  }
  // /* Only if using levelSet
  for ( unsigned int i = 0; i < InterfaceElastoDomains.size(); ++i){
    SpaceReducer spaceReducer(InterfaceElastoDomains[i]->group());
    spaceReducer.BuildLinearConstraints(FSpacesMultLag[i], pAssembler);
  }
  // */
  for ( unsigned int i=0; i < InterfaceElastoDomains.size() ; i++){
    InterfaceElastoDomain* E = InterfaceElastoDomains[i];
    NumberDofs ( FSpacesMultLag[i], E->begin(), E->end(), *pAssembler );
    std::cout << "InterfaceDomain " << i << " : dim=" << E->dim << ", tag=" << E->tag << ", " << E->group()->size() << " elements" << std::endl;
  }
  
}

void ElastodynamicMultLagSolver::AssembleRHS()
{
  for (unsigned int i=0; i<InterfaceElastoDomains.size(); ++i){
      GaussQuadrature Integ_Boundary ( GaussQuadrature::Val );
      size_t v = 100;
      genSimpleFunction<genTensor1<std::complex<double> > >* u_inc = 
	  new genSimpleFunctionUinc<genTensor1<std::complex<double> >, genTensor1<double> >(v,fichier->kp);
      genTerm<genTensor1<std::complex<double> >,0>::Handle Func( new genFunction<genTensor1<std::complex<double> > >( *u_inc ) );   
      genTerm<genTensor0<std::complex<double> >,1>::Handle loadterm = Compose<FullContrProdOp>(FSpacesMultLag[0],Func);
      Assemble ( loadterm , InterfaceElastoDomains[i]->begin(), InterfaceElastoDomains[i]->end(), Integ_Boundary, *pAssembler );
  }
}

void ElastodynamicMultLagSolver::AssembleLHS()
{
  GaussQuadrature Integ_Bulk_val ( GaussQuadrature::ValVal );
  GaussQuadrature Integ_Bulk_grad ( GaussQuadrature::GradGrad );
  for(int i=0; i<ElastodynamicDomains.size(); i++){
    ElastodynamicDomain* E = ElastodynamicDomains[i];
    genTerm<genTensor0<std::complex<double> >,2>::Handle M=Compose<FullContrProdOp>(FSpacesDisp[i],Conj(FSpacesDisp[i]))*std::complex<double>(-pow(E->omega,2)*E->rho,0.0);
    Assemble ( M, E->begin(), E->end(), Integ_Bulk_val, *pAssembler );
    
    genTerm<genTensor0<std::complex<double> >,2>::Handle K = E->Energy<>(EpsilonSpacesDisp[i],Conj(EpsilonSpacesDisp[i]));
    Assemble ( K, E->begin(), E->end(), Integ_Bulk_grad, *pAssembler );
    
    for( int j= 0; j<FSpacesMultLag.size(); j++){
      genTerm<genTensor0<std::complex<double> >,2>::Handle Q= Compose<FullContrProdOp>(FSpacesDisp[i], FSpacesMultLag[j])*pow(-1,i+j+1); // pow argument depend on data in main.dat
      genTerm<genTensor0<std::complex<double> >,2>::Handle TrQ = Compose<FullContrProdOp>(FSpacesMultLag[j], FSpacesDisp[i])*pow(-1,i+j+1); // pow argument depend on data in main.dat
      Assemble ( Q ,InterfaceElastoDomains[j]->begin(),InterfaceElastoDomains[j]->end(),Integ_Bulk_val,*pAssembler );
      Assemble ( TrQ ,InterfaceElastoDomains[j]->begin(),InterfaceElastoDomains[j]->end(),Integ_Bulk_val,*pAssembler );
    }   
  }
  
  for( int i= 0; i<FSpacesMultLag.size(); i++){
      genTerm<genTensor0<std::complex<double> >,2>::Handle NullQ= Compose<FullContrProdOp>(FSpacesMultLag[i], FSpacesMultLag[i])*0.0;
      Assemble ( NullQ ,InterfaceElastoDomains[i]->begin(),InterfaceElastoDomains[i]->end(),Integ_Bulk_val,*pAssembler );
  }
  
  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];  //depend on data in main.dat
	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>(FSpacesDisp[0],FuncComplex_n),FuncComplex_n)*vp
        + Compose<FullContrProdOp>(Compose<FullContrProdOp>(FSpacesDisp[0],FuncComplex_t),FuncComplex_t)*vs; 
	genTerm<genTensor0<std::complex<double> >,2>::Handle f = Compose<FullContrProdOp>(Func, Conj(FSpacesDisp[0]))*std::complex<double>(0.0, -E->omega*E->rho);
        Assemble (f, Data->BCs[i]->begin(), Data->BCs[i]->end(), Integ_Bulk_val, *pAssembler );
	delete n; delete t;
      }
    }
  }

}


void ElastodynamicMultLagSolver::error_norm_L2()
{
  genTensor0<std::complex<double> > sum_diff_L2;
  genTensor0<std::complex<double> > sum_ana_L2;
  for (int i=0;i< ElastodynamicMultLagDomains.size();++i)
  {
      ElastodynamicMultLagDomain *P = ElastodynamicMultLagDomains[i];
      FSpaceDisp = FSpacesDisp[i];
      genTerm<genTensor1<std::complex<double> >,0>::Handle Field_num(new genSField<genTensor1<std::complex<double> > >( pAssembler, FSpaceDisp));
      genSimpleFunction<genTensor1<std::complex<double> > >* f;
      if(i==0){
	f = new genSimpleFunctionCase1<genTensor1<std::complex<double> >, genTensor1<double> >(fichier->vector, 
	  fichier->kp, fichier->ks);
      }
      else if(i==1){
	f = new genSimpleFunctionUtr<genTensor1<std::complex<double> >, 
	  genTensor1<double> >(fichier->vector, fichier->kp_int, fichier->ks_int);
      }
      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* ElastodynamicMultLagSolver::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* ElastodynamicMultLagSolver::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);
    genTerm<genTensor1<std::complex<double> >,0>::Handle FuncExt( new genFunction<genTensor1<std::complex<double> > >( *u_ref ) );
    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 f_int;
  return pv;
}



PView* ElastodynamicMultLagSolver::buildViewDispNum_x(const std::string &postFileName){
    /*
    std::vector<genTerm<genTensor1<double>,0>::Handle> Field_NUM;
    for(int k =0; k< FSpacesDisp.size(); k++){
      genTerm<genTensor1<std::complex<double> >,0>::ConstHandle Field(new genSField<genTensor1<std::complex<double> > >( pAssembler, FSpacesDisp[k] ));
      genTerm<genTensor1<double>,0>::Handle Field_Real(Re(Field));
      Field_NUM.push_back(Re(Field_Real));
    }
    return buildViewNodal(Field_NUM,ElastodynamicDomains,postFileName);
    */
    std::map<int, std::vector<double> > data;
    for(int i=0; i<ElastodynamicDomains.size();  i++){
      ElastodynamicDomain* E = ElastodynamicDomains[i]; 
      genTerm<genTensor1<std::complex<double> >,0>::ConstHandle Field(new genSField<genTensor1<std::complex<double> > >( pAssembler, FSpacesDisp[i] ));
      genTerm<genTensor1<double>,0>::Handle Field_Real(Re(Field));
      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* ElastodynamicMultLagSolver::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;
  
}