// 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>.
#include <string>
#include <cstdio>

#include "GmshConfig.h"
#include "elasticitySolverXfem.h"
#include "linearSystemCSR.h"
#include "linearSystemPETSc.h"
#include "linearSystemGmm.h"
#include "Numeric.h"
#include "GModel.h"
#include "functionSpace.h"
#include "terms.h"
#include "solverAlgorithms.h"
#include "quadratureRules.h"
#include "solverField.h"
#include "MPoint.h"
#include "gmshLevelset.h"
#include "FuncGradDisc.h"
#include "FuncHeaviside.h"
#include "filters.h"

#include "SpaceReducerForStableMultipliers.h"

#include "MTriangle.h"
#include "MHexahedron.h"
#include "MTetrahedron.h"
#include "MQuadrangle.h"
#include "MPrism.h"
#include "MPyramid.h"


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

void readTextileAxes ( const std::string &stringAxes, std::vector<double> &axes)
{
  assert(stringAxes.find('[')==0);
  assert(stringAxes.find(']')==stringAxes.size()-1);

  std::string s =stringAxes.substr(1,stringAxes.size()-2);
  std::vector<std::string> vectAxes;
  
  boost::split( vectAxes, s, boost::is_any_of(",") );
  
  assert(vectAxes.size()==2);
  
  axes[0] = atof(vectAxes[0].c_str());
  axes[1] = atof(vectAxes[1].c_str());
}

void readTextileStructure ( const std::string &stringStructure, fullMatrix<int> &structure)
{
  assert(stringStructure.find('[')==0);
  assert(stringStructure.find(']')==stringStructure.size()-1);
  
  std::string s =stringStructure.substr(1,stringStructure.size()-2);
  std::vector<std::string> chaines;
  std::vector< std::vector<std::string> > trames;
  
  boost::split( chaines, s, boost::is_any_of(";") );
  trames.resize(chaines.size());
  
  int sizeTrame = -1;
  for (int i=0; i<chaines.size(); i++)
  {
    std::vector<std::string> temp;
    boost::split( temp, chaines[i], boost::is_any_of(",") );
    trames[i].resize(temp.size());
    for (int j=0; j<temp.size(); j++)
      trames[i][j]=temp[j];
    
    if (i == 0) sizeTrame=trames[i].size();
    else assert(sizeTrame == trames[i].size());
  }
  structure.resize(chaines.size(), sizeTrame);
  
  for (int i=0; i<chaines.size(); i++)
  {
    for (int j=0; j<sizeTrame; j++)
    {
      structure.set( i, j, atoi(trames[i][j].c_str()) );
    }
  }
}

void readTextileMultipliciteChaine ( const std::string &stringMultipliciteChaine, std::vector<int> &MultipliciteChaine)
{
  assert(stringMultipliciteChaine.find('[')==0);
  assert(stringMultipliciteChaine.find(']')==stringMultipliciteChaine.size()-1);
  
  std::string s =stringMultipliciteChaine.substr(1,stringMultipliciteChaine.size()-2);
  std::vector<std::string> chaines;
  
  boost::split( chaines, s, boost::is_any_of(",") );

  MultipliciteChaine.resize(chaines.size());
  
  for (int i=0; i<chaines.size(); i++)
  {
    MultipliciteChaine[i] = atof(chaines[i].c_str()) ;
  }
}


void elasticitySolverXfem::setMesh ( const std::string &meshFileName )
{
  std::cout << "ls list size : " << _lslist.size() << std::endl;

  _modelname = meshFileName;
  // cut model
  _pModel = _pModelUncut;
  for ( unsigned int i = 0 ; i < _lslist.size() ; i++ )
    _pModel =_pModel->buildCutGModel ( _lslist[i] );

  // write in .msh file the cut model
  std::string ModelName = _modelname ;
  ModelName.erase ( ModelName.find ( ".",0 ),4 );
  ModelName.insert ( ModelName.length(),"_cut.msh" );
  _pModel->writeMSH ( ModelName,2.1,false,false );

  GModel::setCurrent ( _pModel );

  // fill the map between cut-uncut physicals
  setOriginalPhysicals ( _pModel,_pModelUncut );

  // fill the map between level sets entries and physicals
  setLevelSetsPhysicals();

}

void elasticitySolverXfem::readLSInInputFile ( const std::string &fn )
{
  FILE *f = fopen ( fn.c_str(), "r" );
  char what[256];
  while ( !feof ( f ) )
  {
    if ( fscanf ( f, "%s", what ) != 1 ) return;
    if ( !strcmp ( what, "CutMeshPlane" ) )
    {
      double a, b, c, d;
      if ( fscanf ( f, "%lf %lf %lf %lf", &a, &b, &c, &d ) != 4 ) return;
      int tag = _lslist.size() + 1;
      gLevelsetPlane *ls = new gLevelsetPlane ( a,b,c,d,tag );
      _lslist.push_back ( ls );
    }
    else if ( !strcmp ( what, "CutMeshSphere" ) )
    {
      double x, y, z, r;
      if ( fscanf ( f, "%lf %lf %lf %lf", &x, &y, &z, &r ) != 4 ) return;
      int tag = _lslist.size() + 1;
      gLevelsetSphere *ls = new gLevelsetSphere ( x,y,z,r,tag );
      _lslist.push_back ( ls );
    }
    else if ( !strcmp ( what, "CutMeshCylinder" ) )
    {
      double pt[3];//={0.,0.1,0.};
      double dir[3];//={0.,0.0,1.0};
      double r;//=0.015;
      int tag = _lslist.size() + 1;
      if ( fscanf ( f, "%lf %lf %lf %lf %lf %lf %lf", &pt[0], &pt[1], &pt[2], &dir[0], &dir[1], &dir[2], &r ) != 7 ) return;
      gLevelsetGenCylinder *ls = new gLevelsetGenCylinder ( pt,dir,r,tag );
      _lslist.push_back ( ls );
    }
    else if ( !strcmp ( what, "CutMeshTextile" ) )
    {
      double epaisseur, dimChaine, dimTrame;
      char charAxesChaine[256], charAxesTrame[256];
      
      int layer;
      char charStructure[256]; // [0,1,2;2,0,2;0,2,1] ou ";" separarteur des lignes et "," separateur des colonnes (valeurs croissantes avec la profondeur)
      char charMultipliciteChaine[256]; // [0,1,2] ou "," separateur des lignes
      
      int nbdx, nbdyz;

      if ( fscanf ( f, "%lf %lf %lf %s %s %d %s %s %d %d", &epaisseur, &dimChaine, &dimTrame, charAxesChaine, charAxesTrame, &layer, charStructure, charMultipliciteChaine, &nbdx, &nbdyz) != 10 ) return;
      
      // Lecture des axes et de la structure
      std::vector<double> axesChaine(2), axesTrame(2);
      fullMatrix<int> structure;
      std::vector<int> multipliciteChaine;
      readTextileAxes(std::string(charAxesChaine), axesChaine);
      readTextileAxes(std::string(charAxesTrame), axesTrame);
      readTextileStructure(std::string(charStructure), structure);
      readTextileMultipliciteChaine(std::string(charMultipliciteChaine), multipliciteChaine);
      
      int tag = _lslist.size() + 1;
      
      gLevelsetTextile *ls = new gLevelsetTextile ( epaisseur, dimChaine, dimTrame, axesChaine, axesTrame, layer, structure, multipliciteChaine, nbdx, nbdyz, tag );

      for (int i=0; i<ls->getNbLsFils(); i++)
        _lslist.push_back( ls->getLsFil(i) );
    }
  }
  fclose ( f );
}

void elasticitySolverXfem::readMeshInInputFile ( const std::string &fn )
{

  FILE *f = fopen ( fn.c_str(), "r" );
  char what[256];
  while ( !feof ( f ) )
  {
    if ( fscanf ( f, "%s", what ) != 1 ) return;
    if ( !strcmp ( what, "MeshFile" ) )
    {
      char name[245];
      if ( fscanf ( f, "%s", name ) != 1 ) return;

#if 1
      // read the uncut model in a .msh file
      _pModelUncut = new GModel();
      _pModelUncut->readMSH ( name );
      _dim = _pModelUncut->getNumRegions() ? 3 : 2;

      setMesh ( name );
#else 
      
            _pModel = new GModel();
      _pModel->readMSH ( name );
      _dim = _pModel->getNumRegions() ? 3 : 2;
#endif

    }
  }
  fclose ( f );

}

void elasticitySolverXfem::readFieldsInInputFile ( const std::string &fn )
{

  FILE *f = fopen ( fn.c_str(), "r" );
  char what[256];
  while ( !feof ( f ) )
  {
    if ( fscanf ( f, "%s", what ) != 1 ) return;
    if ( !strcmp ( what, "UnCutElasticDomain" ) )
    {
      int physical;
      double E,nu;
      if ( fscanf ( f, "%d %lf %lf", &physical, &E, &nu ) != 3 ) return;

      // link with cut physicals
      for ( unsigned int i = 0 ; i < _OriginalPhysicals[_dim][physical].size() ; i++ )
      {
        elasticField field;
        field._tag = _tag;
        field._E = E;
        field._nu = nu;
        field.g = new groupOfElements ( _dim, _OriginalPhysicals[_dim][physical][i] );
        _elasticFields.push_back ( field );

      }

    }
    else if ( ( !strcmp ( what, "CutElasticDomain" ) ) || ( !strcmp ( what, "ElasticDomain" ) ) )
    {
      elasticField field;
      int physical;
      if ( fscanf ( f, "%d %lf %lf", &physical, &field._E, &field._nu ) != 3 ) return;
      field._tag = _tag;
      field.g = new groupOfElements ( _dim, physical );
      _elasticFields.push_back ( field );
    }
    else if ( !strcmp ( what, "LevelSetElasticDomain" ) )
    {

      int numls;
      double E,nu;
      int sign;
      if ( fscanf ( f, "%lf %lf %d %d", &E, &nu,&numls,&sign ) != 4 ) return;
      std::vector<int> listlevelset;
      std::vector<int> listsign;
      listlevelset.clear();
      listsign.clear();
      listlevelset.push_back ( numls );
      listsign.push_back ( sign );


      if ( _dim == 2 )
      {
        MElement *e;
        typedef std::set<GFace*, GEntityLessThan>::iterator fiter;
        for ( fiter it = _pModel->firstFace(); it != _pModel->lastFace(); ++it )
        {

          // take one element per region
          if ( ( *it )->triangles.size() ) e = ( *it )->triangles[0];
          else if ( ( *it )->quadrangles.size() ) e = ( *it )->quadrangles[0];
          else if ( ( *it )->polygons.size() ) e = ( *it )->polygons[0];
          bool isIn = true;
          for ( unsigned int j = 0 ; j < listlevelset.size() ; j++ )
          {
            double val;
            val = ( * ( _lslist[listlevelset[j]] ) ) ( ( e->barycenter() ).x(), ( e->barycenter() ).y(), ( e->barycenter() ).z() );
            if ( ( val*listsign[j] < 0 ) ) isIn = false ;
          }
          if ( isIn && ( *it )->physicals.size() )
          {
            int physical = ( *it )->physicals[0];
            elasticField field;
            field._tag = _tag;
            field._E = E;
            field._nu = nu;
            field.g = new groupOfElements ( _dim, physical );
            _elasticFields.push_back ( field );
          }
        }
      }

      if ( _dim == 3 )
      {
        MElement *e;
        typedef std::set<GRegion*, GEntityLessThan>::iterator riter;
        for ( riter it = _pModel->firstRegion(); it != _pModel->lastRegion(); ++it )
        {

          // take one element per region
          if ( ( *it )->tetrahedra.size() ) e = ( *it )->tetrahedra[0];
          else if ( ( *it )->hexahedra.size() ) e = ( *it )->hexahedra[0];
          else if ( ( *it )->prisms.size() ) e = ( *it )->prisms[0];
          else if ( ( *it )->pyramids.size() ) e = ( *it )->pyramids[0];
          else if ( ( *it )->polyhedra.size() ) e = ( *it )->polyhedra[0];

          // verify if it s in the domain delimited by level sets sign
          bool isIn = true;
          for ( unsigned int j = 0 ; j < listlevelset.size() ; j++ )
          {
            double val;
            val = ( * ( _lslist[listlevelset[j]] ) ) ( ( e->barycenter() ).x(), ( e->barycenter() ).y(), ( e->barycenter() ).z() );
            if ( ( val*listsign[j] < 0 ) ) isIn = false ;
          }
          if ( isIn && ( *it )->physicals.size() )
          {
            int physical = ( *it )->physicals[0];
            elasticField field;
            field._tag = _tag;
            field._E = E;
            field._nu = nu;
            field.g = new groupOfElements ( _dim, physical );
            _elasticFields.push_back ( field );
          }
        }
      }

    }
    else if ( !strcmp ( what, "Solver" ) )
    {
      if ( fscanf ( f, "%d", &_solvertype ) ) return;
    }
  }
  fclose ( f );

}

void elasticitySolverXfem::readBCInInputFile ( const std::string &fn )
{
  FILE *f = fopen ( fn.c_str(), "r" );
  char what[256];
  std::map<int, std::vector<std::string> > mapExpr;
  std::map<int, std::vector<std::string> > mapVar;
  while ( !feof ( f ) )
  {
    if ( fscanf ( f, "%s", what ) != 1 ) return;
    if ( !strcmp ( what, "NodeDisplacement" ) )
    {
      double val;
      int node, comp;
      if ( fscanf ( f, "%d %d %lf", &node, &comp, &val ) != 3 ) return;
      dirichletBC diri;
      diri.g = new groupOfElements ( 0, node );
      diri._f=  new simpleFunction<double> ( val );
      diri._comp=comp;
      diri._tag=node;
      diri.onWhat=BoundaryCondition::ON_VERTEX;
      _allDirichlet.push_back (diri);
    }
    else if ( !strcmp ( what, "EdgeDisplacement" ) )
    {
      double val;
      int edge, comp;
      if ( fscanf ( f, "%d %d %lf", &edge, &comp, &val ) != 3 ) return;
      
      if (_OriginalPhysicals[1].find(edge) == _OriginalPhysicals[1].end())
      { // link with cut physicals
        dirichletBC diri;
        diri.g = new groupOfElements (1, edge);
        diri._f = new simpleFunction<double> (val);
        diri._comp=comp;
        diri._tag=edge;
        diri.onWhat=BoundaryCondition::ON_EDGE;
        _allDirichlet.push_back (diri);
      }
      else
      { // link with uncut physicals
        for (unsigned int i = 0 ; i < _OriginalPhysicals[1][edge].size() ; i++)
        {
          dirichletBC diri;
          diri.g = new groupOfElements (1, _OriginalPhysicals[1][edge][i]);
          diri._f = new simpleFunction<double> (val);
          diri._comp = comp;
          diri._tag = edge;
          diri.onWhat = BoundaryCondition::ON_EDGE;
          _allDirichlet.push_back (diri);
        }
      }
    }
    else if ( !strcmp ( what, "FaceDisplacement" ) )
    {
      double val;
      int face, comp;
      if ( fscanf ( f, "%d %d %lf", &face, &comp, &val ) != 3 ) return;

      if (_OriginalPhysicals[2].find (face) == _OriginalPhysicals[2].end())
      { // link with cut physicals
        dirichletBC diri;
        diri.g = new groupOfElements (2, face);
        diri._f = new simpleFunction<double> (val);
        diri._comp = comp;
        diri._tag = face;
        diri.onWhat = BoundaryCondition::ON_FACE;
        _allDirichlet.push_back (diri);
      }
      else
      { // link with uncut physicals
        for (unsigned int i = 0 ; i < _OriginalPhysicals[2][face].size() ; i++)
        {
          dirichletBC diri;
          diri.g = new groupOfElements (2, _OriginalPhysicals[2][face][i]);
          diri._f = new simpleFunction<double> (val);
          diri._comp = comp;
          diri._tag = face;
          diri.onWhat = BoundaryCondition::ON_FACE;
          _allDirichlet.push_back (diri);
        }
      }
    }

    else if ( !strcmp ( what, "NodeForce" ) )
    {
      double val1, val2, val3;
      int node;
      if ( fscanf ( f, "%d %lf %lf %lf", &node, &val1, &val2, &val3 ) != 4 ) return;
      neumannBC neu;
      neu.g = new groupOfElements ( 0, node );
      neu._f= new simpleFunction<SVector3> ( SVector3 ( val1, val2, val3 ) );
      neu._tag=node;
      neu.onWhat=BoundaryCondition::ON_VERTEX;
      _allNeumann.push_back ( neu );
    }
    else if ( !strcmp ( what, "EdgeForce" ) )
    {
      double val1, val2, val3;
      int edge;
      if ( fscanf ( f, "%d %lf %lf %lf", &edge, &val1, &val2, &val3 ) != 4 ) return;
      
      if (_OriginalPhysicals[1].find(edge) == _OriginalPhysicals[1].end())
      { // link with cut physicals
        neumannBC neu;
        neu.g = new groupOfElements (1, edge);
        neu._f = new simpleFunction<SVector3> ( SVector3 ( val1, val2, val3 ) );
        neu._tag=edge;
        neu.onWhat=BoundaryCondition::ON_EDGE;
        _allNeumann.push_back ( neu );
      }
      else
      { // link with uncut physicals
        for (unsigned int i = 0 ; i < _OriginalPhysicals[1][edge].size() ; i++)
        {
          neumannBC neu;
          neu.g = new groupOfElements (1, _OriginalPhysicals[1][edge][i]);
          neu._f = new simpleFunction<SVector3> ( SVector3 ( val1, val2, val3 ) );
          neu._tag = edge;
          neu.onWhat = BoundaryCondition::ON_EDGE;
          _allNeumann.push_back (neu);
        }
      }
    }
    else if ( !strcmp ( what, "FaceForce" ) )
    {
      double val1, val2, val3;
      int face;
      if ( fscanf ( f, "%d %lf %lf %lf", &face, &val1, &val2, &val3 ) != 4 ) return;
      
      if (_OriginalPhysicals[2].find (face) == _OriginalPhysicals[2].end())
      { // link with cut physicals
        neumannBC neu;
        neu.g = new groupOfElements (2, face);
        neu._f = new simpleFunction<SVector3> ( SVector3 ( val1, val2, val3 ) );
        neu._tag = face;
        neu.onWhat = BoundaryCondition::ON_FACE;
        _allNeumann.push_back (neu);
      }
      else
      { // link with uncut physicals
      for (unsigned int i = 0 ; i < _OriginalPhysicals[2][face].size() ; i++)
        {
          neumannBC neu;
          neu.g = new groupOfElements (2, _OriginalPhysicals[2][face][i]);
          neu._f = new simpleFunction<SVector3> (SVector3 (val1, val2, val3));
          neu._tag = face;
          neu.onWhat = BoundaryCondition::ON_FACE;
          _allNeumann.push_back (neu);
        }
      }
    }
    else if ( !strcmp ( what, "VolumeForce" ) )
    {
      double val1, val2, val3;
      int volume;
      if ( fscanf ( f, "%d %lf %lf %lf", &volume, &val1, &val2, &val3 ) != 4 ) return;
      neumannBC neu;
      neu.g = new groupOfElements ( 3, volume );
      neu._f= new simpleFunction<SVector3> ( SVector3 ( val1, val2, val3 ) );
      neu._tag=volume;
      neu.onWhat=BoundaryCondition::ON_VOLUME;
      _allNeumann.push_back ( neu );
    }

    else if ( !strcmp ( what, "LagrangeMultipliers" ) )
    {
      LagrangeMultiplierField field;
      int physical;
      double d1, d2, d3, val;
      if ( fscanf ( f, "%d %lf %lf %lf %lf %lf %d", &physical, &field._tau,
                      &d1, &d2, &d3, &val, &field._tag ) != 7 ) return;
      field._tag = _tag;
      field._d = SVector3 ( d1, d2, d3 );
      field._f = new simpleFunction<double> ( val );
      field.g = new groupOfElements ( _dim - 1, physical );
      _LagrangeMultiplierFields.push_back ( field );
    }
    else if ( !strcmp ( what, "StableLagMultField" ) )
    {
      int physical, comp;
      double val;
      if ( fscanf ( f, "%d %d %lf", &physical, &comp, &val ) != 3 ) return;
      LagMultField * oldField = containLagMultField(physical);
      if (oldField==NULL)
      {
        std::vector < groupOfElements * > groupOfElasticFields;
        for ( unsigned int i = 0 ; i < _elasticFields.size() ; i ++ ) groupOfElasticFields.push_back ( _elasticFields[i].g );
        groupOfLagMultElements  * g = new groupOfLagMultElements ( _dim - 1, physical, groupOfElasticFields );
        LagMultField newField ( physical, comp, new simpleFunction<SVector3> (val), g );
        _StableLagrangeMultiplierFields.push_back ( newField );
      }
      else
      {
        oldField->_comp[comp] = true;

        if (mapVar.find(physical) == mapVar.end())
        {
          SVector3 v = (*oldField->_f)(0,0,0);
          double coord[3];
          coord[0]=v.x(); coord[1]=v.y(); coord[2]=v.z();
          coord[comp] = val;
          v = SVector3(coord[0], coord[1], coord[2]) ;
          oldField->changeLagMultField( new simpleFunction<SVector3> ( v ) );
        }
        else
        { // Convert StableLagMultField to StableLagMultFieldFunction
          std::ostringstream oss;
          oss << val;
          mapExpr[physical][comp] = oss.str();
          oldField->changeLagMultField(new MathEvalFunctionSVector3 (mapExpr[physical], mapVar[physical]));
        }

      }
    }
//
//  Complete the readBCInInputFile methode with functions
//
    else if (!strcmp(what, "NodeDisplacementFunction"))
    {
      char exprXYZ[256];
      char x[256], y[256], z[256];
      int node, comp;
      if (fscanf(f, "%d %d %s %s %s %s", &node, &comp, exprXYZ, x, y, z) != 6) return;
      
      std::vector<std::string> expr(1),var(3); 
      expr[0]=exprXYZ; var[0]=x; var[1]=y; var[2]=z;
      
      dirichletBC diri;
      diri.g = new groupOfElements(0, node);
      diri._f= new MathEvalFunctionDouble(expr,var);
      diri._comp=comp;
      diri._tag=node;
      diri.onWhat=BoundaryCondition::ON_VERTEX;
      _allDirichlet.push_back(diri);
    }
    else if (!strcmp(what, "EdgeDisplacementFunction"))
    {
      char exprXYZ[256];
      char x[256], y[256], z[256];
      int edge, comp;
      if (fscanf(f, "%d %d %s %s %s %s", &edge, &comp, exprXYZ, x, y, z) != 6) return;
      
      std::vector<std::string> expr(1),var(3); 
      expr[0]=exprXYZ; var[0]=x; var[1]=y; var[2]=z;
      
      if (_OriginalPhysicals[1].find(edge) == _OriginalPhysicals[1].end())
      { // link with cut physicals
        std::vector<MVertex*> vertices;
        _pModel->getMeshVerticesForPhysicalGroup (1, edge, vertices);

        dirichletBC diri;
        diri.g = new groupOfElements (1, edge);
        diri._f = new MathEvalFunctionDouble (expr, var);
        diri._comp = comp;
        diri._tag = edge;
        diri.onWhat = BoundaryCondition::ON_EDGE;
        _allDirichlet.push_back (diri);
      }
      else
      { // link with cut physicals
      for (unsigned int i = 0 ; i < _OriginalPhysicals[1][edge].size() ; i++)
        {
          std::vector<MVertex*> vertices;
          _pModel->getMeshVerticesForPhysicalGroup (1, edge, vertices);

          dirichletBC diri;
          diri.g = new groupOfElements (1, _OriginalPhysicals[1][edge][i]);
          diri._f = new MathEvalFunctionDouble (expr, var);
          diri._comp = comp;
          diri._tag = edge;
          diri.onWhat = BoundaryCondition::ON_EDGE;
          _allDirichlet.push_back (diri);
        }
      }
    }
    else if (!strcmp(what, "FaceDisplacementFunction"))
    {
      char exprXYZ[256];
      char x[256], y[256], z[256];
      int face, comp;
      if (fscanf(f, "%d %d %s %s %s %s", &face, &comp, exprXYZ, x, y, z) != 6) return;
      
      std::vector<std::string> expr(1),var(3); 
      expr[0]=exprXYZ; var[0]=x; var[1]=y; var[2]=z;
      
      if (_OriginalPhysicals[2].find(face) == _OriginalPhysicals[2].end())
      { // link with cut physicals
        std::vector<MVertex*> vertices;
        _pModel->getMeshVerticesForPhysicalGroup (2, face, vertices);

        dirichletBC diri;
        diri.g = new groupOfElements (2, face);
        diri._f = new MathEvalFunctionDouble (expr, var);
        diri._comp = comp;
        diri._tag = face;
        diri.onWhat = BoundaryCondition::ON_FACE;
        _allDirichlet.push_back (diri);
      }
      else
      { // link with uncut physicals
        for (unsigned int i = 0 ; i < _OriginalPhysicals[2][face].size() ; i++)
        {
          std::vector<MVertex*> vertices;
          _pModel->getMeshVerticesForPhysicalGroup (2, face, vertices);

          dirichletBC diri;
          diri.g = new groupOfElements (2, _OriginalPhysicals[2][face][i]);
          diri._f = new MathEvalFunctionDouble (expr, var);
          diri._comp = comp;
          diri._tag = face;
          diri.onWhat = BoundaryCondition::ON_FACE;
          _allDirichlet.push_back (diri);
        }
      }
    }
    
    else if (!strcmp(what, "NodeForceFunction"))
    {
      char exprX[256],exprY[256],exprZ[256];
      char x[256],y[256],z[256];
      int node;
      if (fscanf(f, "%d %s %s %s %s %s %s", &node, exprX, exprY, exprZ, x, y, z) != 7) return;
      
      std::vector<std::string> expr(3),var(3); 
      expr[0]=exprX; var[0]=x;
      expr[1]=exprY; var[1]=y;
      expr[2]=exprZ; var[2]=z;
      
      neumannBC neu;
      neu.g = new groupOfElements (0, node);
      neu._f= new MathEvalFunctionSVector3(expr,var);
      neu._tag=node;
      neu.onWhat=BoundaryCondition::ON_VERTEX;
      _allNeumann.push_back(neu);
    }
    else if (!strcmp(what, "EdgeForceFunction"))
    {
      char exprX[256],exprY[256],exprZ[256];
      char x[256],y[256],z[256];
      int edge;
      if (fscanf(f, "%d %s %s %s %s %s %s", &edge, exprX, exprY, exprZ, x, y, z) != 7) return;
      
      std::vector<std::string> expr(3),var(3); 
      expr[0]=exprX; var[0]=x;
      expr[1]=exprY; var[1]=y;
      expr[2]=exprZ; var[2]=z;
      
      if (_OriginalPhysicals[1].find (edge) == _OriginalPhysicals[1].end())
      { // link with cut physicals
        std::vector<MVertex*> vertices;
        _pModel->getMeshVerticesForPhysicalGroup(1, edge, vertices);
        
        neumannBC neu;
        neu.g = new groupOfElements (1, edge);
        neu._f= new MathEvalFunctionSVector3(expr,var);
        neu._tag=edge;
        neu.onWhat=BoundaryCondition::ON_EDGE;
        _allNeumann.push_back ( neu );
      }
      else
      { // link with uncut physicals
        for (unsigned int i = 0 ; i < _OriginalPhysicals[1][edge].size() ; i++)
        {
          std::vector<MVertex*> vertices;
          _pModel->getMeshVerticesForPhysicalGroup(1, edge, vertices);
          
          neumannBC neu;
          neu.g = new groupOfElements (1, _OriginalPhysicals[1][edge][i]);
          neu._f= new MathEvalFunctionSVector3(expr,var);
          neu._tag = edge;
          neu.onWhat = BoundaryCondition::ON_EDGE;
          _allNeumann.push_back (neu);
        }
      }
    }
    else if (!strcmp(what, "FaceForceFunction"))
    {
      char exprX[256],exprY[256],exprZ[256];
      char x[256],y[256],z[256];
      int face;
      if (fscanf(f, "%d %s %s %s %s %s %s", &face, exprX, exprY, exprZ, x, y, z) != 7) return;
      
      std::vector<std::string> expr(3),var(3); 
      expr[0]=exprX; var[0]=x;
      expr[1]=exprY; var[1]=y;
      expr[2]=exprZ; var[2]=z;
      
      if (_OriginalPhysicals[2].find (face) == _OriginalPhysicals[2].end())
      { // link with cut physicals
        std::vector<MVertex*> vertices;
        _pModel->getMeshVerticesForPhysicalGroup(2, face, vertices);
        
        neumannBC neu;
        neu.g = new groupOfElements (2, face);
        neu._f= new MathEvalFunctionSVector3(expr,var);
        neu._tag = face;
        neu.onWhat = BoundaryCondition::ON_FACE;
        _allNeumann.push_back (neu);
      }
      else
      { // link with uncut physicals
      for (unsigned int i = 0 ; i < _OriginalPhysicals[2][face].size() ; i++)
        {
          std::vector<MVertex*> vertices;
          _pModel->getMeshVerticesForPhysicalGroup(2, face, vertices);
      
          neumannBC neu;
          neu.g = new groupOfElements (2, _OriginalPhysicals[2][face][i]);
          neu._f= new MathEvalFunctionSVector3(expr,var);
          neu._tag = face;
          neu.onWhat = BoundaryCondition::ON_FACE;
          _allNeumann.push_back (neu);
        }
      }
    }
    else if (!strcmp(what, "VolumeForceFunction"))
    {
      char exprX[256],exprY[256],exprZ[256];
      char x[256],y[256],z[256];
      int volume;
      if (fscanf(f, "%d %s %s %s %s %s %s", &volume, exprX, exprY, exprZ, x, y, z) != 7) return;
      
      std::vector<std::string> expr(3),var(3); 
      expr[0]=exprX; var[0]=x;
      expr[1]=exprY; var[1]=y;
      expr[2]=exprZ; var[2]=z;
      
      neumannBC neu;
      neu.g = new groupOfElements (3, volume);
      neu._f= new MathEvalFunctionSVector3(expr,var);
      neu._tag=volume;
      neu.onWhat=BoundaryCondition::ON_VOLUME;
      _allNeumann.push_back(neu);
    }

    else if ( !strcmp ( what, "LagrangeMultipliersFunction" ) )
    {
      LagrangeMultiplierField field;
      char exprXYZ[256];
      char x[256], y[256], z[256];
      double tau, d1, d2, d3;
      int physical, tag;
      if ( fscanf ( f, "%d %lf %d %lf %lf %lf %s %s %s %s", &physical, &tau, &tag, &d1, &d2, &d3, exprXYZ, x, y, z) != 10 ) return;
      
      
      std::vector<std::string> expr(1),var(3); 
      expr[0]=exprXYZ; var[0]=x; var[1]=y; var[2]=z;
      
      field._tau = tau;
      field._tag = _tag;
      field._d = SVector3 ( d1, d2, d3 );
      field._f = new MathEvalFunctionDouble(expr,var);
      field.g = new groupOfElements ( _dim - 1, physical );
      _LagrangeMultiplierFields.push_back ( field );
    }
    else if ( !strcmp ( what, "StableLagMultFieldFunction" ) )
    {
      int physical;
      int comp;
      char exprXYZ[256];
      char x[256], y[256], z[256];
      if ( fscanf ( f, "%d %d %s %s %s %s", &physical, &comp, exprXYZ, x, y, z) != 6 ) return;

      LagMultField * oldField = containLagMultField(physical);
      if (oldField==NULL)
      {
        std::vector<std::string> expr(3), var(3);
        expr[0]=expr[1]=expr[2]="0";
        expr[comp]=exprXYZ;
        var[0]=x; var[1]=y; var[2]=z;
        
        std::vector < groupOfElements * > groupOfElasticFields;
        for ( unsigned int i = 0 ; i < _elasticFields.size() ; i ++ ) groupOfElasticFields.push_back ( _elasticFields[i].g );
        groupOfLagMultElements  * g = new groupOfLagMultElements ( _dim - 1, physical, groupOfElasticFields );
        LagMultField newField ( physical, comp, new MathEvalFunctionSVector3 (expr,var), g );
        _StableLagrangeMultiplierFields.push_back ( newField );
        mapExpr[physical] = expr;
        mapVar[physical] = var;
      }
      else
      {
        oldField->_comp[comp] = true;
        
        if (mapVar.find(physical) == mapVar.end())
        { // Convert StableLagMultField to StableLagMultFieldFunction
          SVector3 v = (*oldField->_f)(0,0,0);
          std::vector<std::string> expr(3), var(3);
          std::ostringstream oss;
          oss << v.x(); expr[0] = oss.str(); oss.str("");
          oss << v.y(); expr[1] = oss.str(); oss.str("");
          oss << v.z(); expr[2] = oss.str();
          mapExpr[physical] = expr ;

          var[0]=x; var[1]=y; var[2]=z;
          mapVar[physical] = var;
        }

        mapExpr[physical][comp] = exprXYZ; 
        
        oldField->changeLagMultField(new MathEvalFunctionSVector3 (mapExpr[physical], mapVar[physical])); // use the same variables
      }
    }
  }
  fclose(f);
}

void elasticitySolverXfem::readInputFile ( const std::string &fn )
{

  _modelname = fn;
  std::cout<<"model name : " << _modelname << std::endl ;
  readLSInInputFile ( fn );
  readMeshInInputFile ( fn );
  readFieldsInInputFile ( fn );
  readBCInInputFile ( fn );
}

void elasticitySolverXfem::BuildFunctionSpaces()
{
  mechanicsSolver::BuildFunctionSpaces();

//  if ( _dim==2 )
//     _StableLagrangeMultipliersSpace = new VectorLagrangeFunctionSpace ( _tag+2,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y );
//   if ( _dim==3 )
//     _StableLagrangeMultipliersSpace = new VectorLagrangeFunctionSpace ( _tag+2 );
  _StableLagrangeMultipliersSpace.resize(_StableLagrangeMultiplierFields.size());
  for ( unsigned int i = 0 ; i < _StableLagrangeMultiplierFields.size(); i++ )
  {
    VectorLagrangeFunctionSpace::Along vectLag[3];
    vectLag[0] = VectorLagrangeFunctionSpace::VECTOR_X;
    vectLag[1] = VectorLagrangeFunctionSpace::VECTOR_Y;
    vectLag[2] = VectorLagrangeFunctionSpace::VECTOR_Z;

    std::vector<int> vect;
    for ( unsigned int j = 0; j < _dim; j++)
      if (_StableLagrangeMultiplierFields[i]._comp[j]) vect.push_back(j);

    switch (vect.size())
    {
      case 1 :
        _StableLagrangeMultipliersSpace[i] = new VectorLagrangeFunctionSpace ( _tag+2, vectLag[vect[0]]);
        break;
      case 2 :
        _StableLagrangeMultipliersSpace[i] = new VectorLagrangeFunctionSpace ( _tag+2, vectLag[vect[0]], vectLag[vect[1]]);
        break;
      case 3 :
        _StableLagrangeMultipliersSpace[i] = new VectorLagrangeFunctionSpace ( _tag+2, vectLag[vect[0]], vectLag[vect[1]], vectLag[vect[2]]);
        break;
      default :
        Msg::Error("StableLagrangeMultipliersSpace size error");
    }
  }

  // build linear constraints, (to do fix dof when no bulk element associate with)
  for ( unsigned int i = 0 ; i < _StableLagrangeMultiplierFields.size(); i++ )
  {
    SpaceReducerForStableMultipliers *spacereducer = new SpaceReducerForStableMultipliers ( _StableLagrangeMultiplierFields[i]._g ,/* _lslist[i],*/_tag+2 );
    spacereducer->BuildLinearConstraints ( _StableLagrangeMultiplierFields[i]._comp, _pAssembler );
  }


//    SpaceReducerForStableMultipliers *spacereducer;
//
//    spacereducer = new SpaceReducerForStableMultipliers(std::pair < int, int >(2,2), _lslist[1],_tag+2);
//    spacereducer->BuildLinearConstraints(_pAssembler);
//    delete spacereducer;
//    spacereducer = new SpaceReducerForStableMultipliers(std::pair < int, int >(2,3), _lslist[2],_tag+2);
//    spacereducer->BuildLinearConstraints(_pAssembler);
//    delete spacereducer;



  // Gaetan method
  if ( _LagrangeMultiplierSpace ) delete _LagrangeMultiplierSpace;
  //LagrangeMultiplierSpace = new ScalarLagrangeFunctionSpace(_tag+1);  // here we take parent vertex for dof...
  _LagrangeMultiplierSpace = new ScalarLagrangeFunctionSpaceOfElement ( _tag+1 ); // here level set vertices ?

  /*

      // fill mean hanging nodes

  if (_ListHangingNodes)
  {
      std::map<int,std::vector <int> >::iterator ith;
      ith = _ListHangingNodes->begin();
      int compt = 1;
      while (ith!=_ListHangingNodes->end())
      {
          float fac;
          fac = 1.0/(ith->second).size();  // mean hanging nodes
          for (int j=0;j<_dim;j++)
          {
              DofAffineConstraint<double> constraint;
              int type = Dof::createTypeWithTwoInts(j, _tag);
              Dof hgnd(ith->first,type);
              for (int i=0;i<(ith->second).size();i++)
              {
                  Dof key((ith->second)[i],type);
                  std::pair<Dof, double >  linDof(key,fac);
                  constraint.linear.push_back(linDof);
              }
              constraint.shift = 0;
              _pAssembler->setLinearConstraint (hgnd, constraint);
          }
          ith++;
          compt++;
      }
  }
  */


  // we number the LagrangeMultipliers dofs
  for ( unsigned int i = 0; i < _StableLagrangeMultiplierFields.size(); ++i ) // E.Bechet Method
  {
    NumberDofs ( *_StableLagrangeMultipliersSpace[i], _StableLagrangeMultiplierFields[i]._g->begin(), _StableLagrangeMultiplierFields[i]._g->end(), *_pAssembler );
  }
  for ( unsigned int i = 0; i < _LagrangeMultiplierFields.size(); ++i ) // Gaetan method
  {
    NumberDofs ( *_LagrangeMultiplierSpace, _LagrangeMultiplierFields[i].g->begin(), _LagrangeMultiplierFields[i].g->end(), *_pAssembler );
  }

}

void elasticitySolverXfem::BuildLinearSystem()
{
  mechanicsSolver::BuildLinearSystem();
  GaussQuadrature Integ_Boundary ( GaussQuadrature::Val );

//   // assemble stable lag mult terms  E. Bechet Method
  std::cout <<  "Cross Term, and neumann stable lag mult"<< std::endl;
  GaussQuadrature Integ_LagrangeMult ( GaussQuadrature::ValVal );
  for ( unsigned int i = 0; i < _StableLagrangeMultiplierFields.size(); i++ )
  {
    double eqfac = _elasticFields[0]._E ;  // ? equations scaling ?
    // eqfac = 1;
    double k = 0; // stiffness inverse for the condition (0 = infinite stiffness)
    std::cout << "stiffness inverse k = " << k << std::endl ;
    LagMultTerm CrossLagTerm ( *_FSpace, *_StableLagrangeMultipliersSpace[i], eqfac );
    LagMultTerm SelfLagTerm ( *_StableLagrangeMultipliersSpace[i], *_StableLagrangeMultipliersSpace[i], k );
    LoadTermOnBorder<SVector3> Lterm ( *_StableLagrangeMultipliersSpace[i], *(_StableLagrangeMultiplierFields[i]._f),eqfac );
    Assemble ( CrossLagTerm, *_FSpace, *_StableLagrangeMultipliersSpace[i],
               _StableLagrangeMultiplierFields[i]._g->begin(), _StableLagrangeMultiplierFields[i]._g->end(), Integ_LagrangeMult, *_pAssembler );
    Assemble ( SelfLagTerm,*_StableLagrangeMultipliersSpace[i],
               _StableLagrangeMultiplierFields[i]._g->begin(), _StableLagrangeMultiplierFields[i]._g->end(), Integ_LagrangeMult, *_pAssembler );
    Assemble ( Lterm, *_StableLagrangeMultipliersSpace[i],
               _StableLagrangeMultiplierFields[i]._g->begin(), _StableLagrangeMultiplierFields[i]._g->end(), Integ_Boundary, *_pAssembler );
  }

  // assemble cross term, laplace term and rhs term for LM. Gaetan Method
  GaussQuadrature Integ_Laplace ( GaussQuadrature::GradGrad );
  for ( unsigned int i = 0; i < _LagrangeMultiplierFields.size(); i++ )
  {
    LagrangeMultiplierTerm LagTerm ( *_FSpace, *_LagrangeMultiplierSpace, _LagrangeMultiplierFields[i]._d );
    Assemble ( LagTerm, *_FSpace, *_LagrangeMultiplierSpace,
               _LagrangeMultiplierFields[i].g->begin(), _LagrangeMultiplierFields[i].g->end(), Integ_LagrangeMult, *_pAssembler );

    LaplaceTerm<double,double> LapTerm ( *_LagrangeMultiplierSpace, _LagrangeMultiplierFields[i]._tau );
    Assemble ( LapTerm, *_LagrangeMultiplierSpace,
               _LagrangeMultiplierFields[i].g->begin(), _LagrangeMultiplierFields[i].g->end(), Integ_Laplace, *_pAssembler );

    LoadTermOnBorder<double> Lterm ( *_LagrangeMultiplierSpace, *(_LagrangeMultiplierFields[i]._f) );
    Assemble ( Lterm, *_LagrangeMultiplierSpace,
               _LagrangeMultiplierFields[i].g->begin(),_LagrangeMultiplierFields[i].g->end(), Integ_Boundary, *_pAssembler );
  }

  printf ( "nDofs=%d\n",_pAssembler->sizeOfR() );

  printf ( "-- done assembling!\n" );


  if ( _solvertype == 1 )
    exportKb();
}

void elasticitySolverXfem::setLevelSetsPhysicals()
{
  std::map <int,int> lsphysicals;
  double fac = 1e-2; // warning : too mesh size dependent
  if ( _dim == 3 )
  {
    MElement *e;
    typedef std::set<GFace*, GEntityLessThan>::iterator fiter;
    for ( fiter it = _pModel->firstFace(); it != _pModel->lastFace(); ++it )
    {
      // take one element per face
      if ( ( *it )->triangles.size() ) e = ( *it )->triangles[0];
      else if ( ( *it )->quadrangles.size() ) e = ( *it )->quadrangles[0];
      else if ( ( *it )->polygons.size() ) e = ( *it )->polygons[0];
      for ( unsigned int j = 0 ; j < _lslist.size() ; j++ )
      {
        double val;
        val = ( * ( _lslist[j] ) ) ( ( e->barycenter() ).x(), ( e->barycenter() ).y(), ( e->barycenter() ).z() );
        if ( ( std::abs ( val ) < e->getInnerRadius() *fac ) && ( *it )->physicals.size() )   // too arbitrary condition
        {
          lsphysicals[j] = ( *it )->physicals[0];
          break; // one level set per physical ; out of the for (lslist) loop
        }
      }
    }
  }
  if ( _dim == 2 )
  {
    MElement *e;
    typedef std::set<GEdge*, GEntityLessThan>::iterator eiter;
    for ( eiter it = _pModel->firstEdge(); it != _pModel->lastEdge(); ++it )
    {
      // take one element per face
      if ( ( *it )->lines.size() ) e = ( *it )->lines[0];
      for ( unsigned int j = 0 ; j < _lslist.size() ; j++ )
      {
        double val;
        val = ( * ( _lslist[j] ) ) ( ( e->barycenter() ).x(), ( e->barycenter() ).y(), ( e->barycenter() ).z() );
        if ( ( std::abs ( val ) < e->getInnerRadius() *fac ) && ( *it )->physicals.size() )  // too arbitrary condition
        {
          lsphysicals[j] = ( *it )->physicals[0];
          break; // one level set per physical ; out of the for (lslist) loop
        }
      }
    }
  }


  for ( unsigned int i = 0 ; i < _lslist.size() ; i++ )
  {
    std::pair < int,int > LevelSetEntity;
    LevelSetEntity.first = _dim - 1; // level set dim
    // LevelSetEntity.second = lsphysicals[i] + 1 ;  // tag warning !!
    LevelSetEntity.second = lsphysicals[i] ;  // tag
    _lsEntitylist.push_back ( LevelSetEntity );
  }

  std::cout<<"-----------------------------------------"<<std::endl;
  std::cout<<"Warning : setLevelSetsPhysicals results..."<<std::endl;

  std::map < int, int >::iterator it;
  it = lsphysicals.begin();
  for ( ;it!=lsphysicals.end();it++ )
  {
    std::cout << ( *it ).first << " -> " << ( *it ).second << std::endl ;
  }
  std::cout<<"-----------------------------------------"<<std::endl;

}

void elasticitySolverXfem::setOriginalPhysicals ( GModel *pModel, GModel *pModelCut )
{

  typedef std::set<GFace*, GEntityLessThan>::iterator fiter;
  typedef std::set<GEdge*, GEntityLessThan>::iterator eiter;
  typedef std::set<GRegion*, GEntityLessThan>::iterator riter;

  double eps = 1e-12;

  for ( eiter it = _pModel->firstEdge(); it != _pModel->lastEdge(); ++it )
  {

    //int ElementNum;
    SPoint3 bary1;
    std::pair<int,int> physicalpair;

    bool found;
    found = false;
    // take the first MLine in the Edge entity
    MElement *ec = ( *it )->lines[0];
    while ( ec->getParent() ) ec = ec->getParent();
    //ElementNum = ec->getNum();
    bary1 = ec->barycenter();

    // one physical per entity assumption
    if ( ( *it )->physicals[0] )
    {
      physicalpair.first = ( *it )->physicals[0];  // cut element physical tag
      // search in uncut entities
      for ( eiter itc = _pModelUncut->firstEdge(); itc != _pModelUncut->lastEdge() && !found; ++itc )
      {
        if ( ( *itc )->physicals[0] )
          for ( unsigned int i = 0; i < ( *itc )->lines.size(); i++ )
          {
            MElement *ep;
            ep = ( *itc )->lines[i];
            //if (ep->getNum() == ElementNum)
            if ( bary1.distance ( ep->barycenter() ) < eps )
            {
              physicalpair.second = ( *itc )->physicals[0];  // original element physical tag
              std::vector<int>::iterator itv;
              itv = _OriginalPhysicals[1][physicalpair.second].begin();
              // research in vector...
              while ( itv != _OriginalPhysicals[1][physicalpair.second].end() )
              {
                if ( ( *itv ) == physicalpair.first ) break;
                else itv++;
              }
              if ( itv == _OriginalPhysicals[1][physicalpair.second].end() )
              {
                _OriginalPhysicals[1][physicalpair.second].push_back ( physicalpair.first );
                found = true;
                break;
              }
            }
          }
      }
    }
  }


  for ( fiter it = _pModel->firstFace(); it != _pModel->lastFace(); ++it )
  {

    //int ElementNum;
    SPoint3 bary1;
    std::pair<int,int> physicalpair;
    unsigned int numElements[3];
    numElements[0]=0;
    numElements[1]=0;
    numElements[2]=0;

    bool found;
    found = false;

    // get num mesh element in entity by type
    ( *it )->getNumMeshElements ( numElements );

    MElement *ec;
    // take the first MElement in the face entity
    if ( numElements[2] )
    {
      ec = ( *it )->polygons[0];
      while ( ec->getParent() ) ec = ec->getParent();
      //ElementNum = ep->getNum();
      bary1 = ec->barycenter();
    }
    else if ( numElements[0] )
    {
      ec = ( *it )->triangles[0] ;
      //ElementNum = (*it)->triangles[0]->getNum();
      bary1 = ec->barycenter();
    }
    else if ( numElements[1] )
    {
      ec = ( *it )->quadrangles[0] ;
      //ElementNum = (*it)->quadrangles[0]->getNum();
      bary1 = ec->barycenter();
    }

    // one physical per entity assumption
    if ( ( *it )->physicals[0] )
    {
      physicalpair.first = ( *it )->physicals[0];  // cut tag
      // search in uncut entities
      for ( fiter itc = _pModelUncut->firstFace(); itc != _pModelUncut->lastFace() && !found ; ++itc )
      {
        if ( ( *itc )->physicals[0] )
        {
          if ( numElements[0] )
          {
            for ( unsigned int i = 0; i < ( *itc )->triangles.size(); i++ )
            {
              MElement *ep = ( *itc )->triangles[i];
              if ( bary1.distance ( ep->barycenter() ) < eps )
              {
                physicalpair.second = ( *itc )->physicals[0];  // original physical
                std::vector<int>::iterator itv;
                itv = _OriginalPhysicals[2][physicalpair.second].begin();
                // research in vector...
                while ( itv != _OriginalPhysicals[2][physicalpair.second].end() )
                {
                  if ( ( *itv ) == physicalpair.first ) break;
                  else itv++;
                }
                if ( itv == _OriginalPhysicals[2][physicalpair.second].end() )
                {
                  _OriginalPhysicals[2][physicalpair.second].push_back ( physicalpair.first );
                  found = true;
                  break ;
                }
              }
            }
          }
          else
          {
            if ( numElements[1] )
            {
              for ( unsigned int i = 0; i < ( *itc )->quadrangles.size(); i++ )
              {
                MElement *ep = ( *itc )->quadrangles[i];
                if ( bary1.distance ( ep->barycenter() ) < eps )
                {
                  physicalpair.second = ( *itc )->physicals[0];  // original physical
                  std::vector<int>::iterator itv;
                  itv = _OriginalPhysicals[2][physicalpair.second].begin();
                  // research in vector...
                  while ( itv != _OriginalPhysicals[2][physicalpair.second].end() )
                  {
                    if ( ( *itv ) == physicalpair.first ) break;
                    else itv++;
                  }
                  if ( itv == _OriginalPhysicals[2][physicalpair.second].end() )
                  {
                    _OriginalPhysicals[2][physicalpair.second].push_back ( physicalpair.first );
                    found = true;
                    break ;
                  }
                }
              }
            }
          }
        }
      }
    }
  }


  for ( riter it = _pModel->firstRegion(); it != _pModel->lastRegion(); ++it )
  {

    //int ElementNum;
    SPoint3 bary1;
    std::pair<int,int> physicalpair;
    unsigned int numElements[5];
    numElements[0]=0;
    numElements[1]=0;
    numElements[2]=0;
    numElements[3]=0;
    numElements[4]=0;

    bool found;
    found = false;

    // get num mesh element in entity by type
    ( *it )->getNumMeshElements ( numElements );

    MElement *ec;
    // take the first MElement in the face entity
    if ( numElements[4] )
    {
      ec = ( *it )->polyhedra[0];
      while ( ec->getParent() ) ec = ec->getParent();
      //ElementNum = ep->getNum();
      bary1 = ec->barycenter();
    }
    else if ( numElements[0] )
    {
      ec = ( *it )->tetrahedra[0] ;
      //ElementNum = (*it)->triangles[0]->getNum();
      bary1 = ec->barycenter();
    }
    else if ( numElements[1] )
    {
      ec = ( *it )->hexahedra[0] ;
      //ElementNum = (*it)->quadrangles[0]->getNum();
      bary1 = ec->barycenter();
    } // to be continued...

    // one physical per entity assumption
    if ( ( *it )->physicals[0] )
    {
      physicalpair.first = ( *it )->physicals[0];  // cut tag
      // search in uncut entities
      for ( riter itc = _pModelUncut->firstRegion(); itc != _pModelUncut->lastRegion() && !found ; ++itc )
      {
        if ( ( *itc )->physicals[0] )
        {
          if ( numElements[0] )
          {
            for ( unsigned int i = 0; i < ( *itc )->tetrahedra.size(); i++ )
            {
              MElement *ep = ( *itc )->tetrahedra[i];
              if ( bary1.distance ( ep->barycenter() ) < eps )
              {
                physicalpair.second = ( *itc )->physicals[0];  // original physical
                std::vector<int>::iterator itv;
                itv = _OriginalPhysicals[3][physicalpair.second].begin();
                // research in vector...
                while ( itv != _OriginalPhysicals[3][physicalpair.second].end() )
                {
                  if ( ( *itv ) == physicalpair.first ) break;
                  else itv++;
                }
                if ( itv == _OriginalPhysicals[3][physicalpair.second].end() )
                {
                  _OriginalPhysicals[3][physicalpair.second].push_back ( physicalpair.first );
                  found = true;
                  break ;
                }
              }
            }
          }
          else
          {
            if ( numElements[1] )
            {
              for ( unsigned int i = 0; i < ( *itc )->hexahedra.size(); i++ )
              {
                MElement *ep = ( *itc )->hexahedra[i];
                if ( bary1.distance ( ep->barycenter() ) < eps )
                {
                  physicalpair.second = ( *itc )->physicals[0];  // original physical
                  std::vector<int>::iterator itv;
                  itv = _OriginalPhysicals[3][physicalpair.second].begin();
                  // research in vector...
                  while ( itv != _OriginalPhysicals[3][physicalpair.second].end() )
                  {
                    if ( ( *itv ) == physicalpair.first ) break;
                    else itv++;
                  }
                  if ( itv == _OriginalPhysicals[3][physicalpair.second].end() )
                  {
                    _OriginalPhysicals[3][physicalpair.second].push_back ( physicalpair.first );
                    found = true;
                    break ;
                  }
                }
              }
            }
          }
        }
      }
    }
  }

  std::cout<<"-----------------------------------------"<<std::endl;
  std::cout<<"Warning : setOriginalphysicals results..."<<std::endl;
  std::cout<<"lines : "<<std::endl;
  std::map < int, std::vector < int > >::iterator it;
  it = _OriginalPhysicals[1].begin();
  for ( ;it!=_OriginalPhysicals[1].end();it++ )
  {
    for ( unsigned int i = 0 ;  i < ( *it ).second.size() ; i++ )
      std::cout << ( *it ).first << " -> " << ( *it ).second[i] << std::endl ;
  }

  std::cout<<"faces :"<<std::endl;

  it = _OriginalPhysicals[2].begin();
  for ( ;it!=_OriginalPhysicals[2].end();it++ )
  {
    for ( unsigned int i = 0 ;  i < ( *it ).second.size() ; i++ )
      std::cout<< ( *it ).first << " -> " << ( *it ).second[i] << std::endl ;
  }

  std::cout<<"regions :"<<std::endl;

  it = _OriginalPhysicals[3].begin();
  for ( ;it!=_OriginalPhysicals[3].end();it++ )
  {
    for ( unsigned int i = 0 ;  i < ( *it ).second.size() ; i++ )
      std::cout<< ( *it ).first << " -> " << ( *it ).second[i] << std::endl ;
  }

  std::cout<<"-----------------------------------------"<<std::endl;

}



#if defined(HAVE_POST)

static double vonMises ( dofManager<double> *a, MElement *e,
                         double u, double v, double w,
                         double E, double nu, int _tag )
{
  double valx[256];
  double valy[256];
  double valz[256];
  for ( int k = 0; k < e->getNumVertices(); k++ )
  {
    a->getDofValue ( e->getVertex ( k ), 0, _tag, valx[k] );
    a->getDofValue ( e->getVertex ( k ), 1, _tag, valy[k] );
    a->getDofValue ( e->getVertex ( k ), 2, _tag, valz[k] );
  }
  double gradux[3];
  double graduy[3];
  double graduz[3];
  e->interpolateGrad ( valx, u, v, w, gradux );
  e->interpolateGrad ( valy, u, v, w, graduy );
  e->interpolateGrad ( valz, u, v, w, graduz );

  double eps[6] = {gradux[0], graduy[1], graduz[2],
                   0.5 * ( gradux[1] + graduy[0] ),
                   0.5 * ( gradux[2] + graduz[0] ),
                   0.5 * ( graduy[2] + graduz[1] )
                  };
  double A = E / ( 1. + nu );
  double B = A * ( nu / ( 1. - 2 * nu ) );
  double trace = eps[0] + eps[1] + eps[2] ;
  double sxx = A * eps[0] + B * trace;
  double syy = A * eps[1] + B * trace;
  double szz = A * eps[2] + B * trace;
  double sxy = A * eps[3];
  double sxz = A * eps[4];
  double syz = A * eps[5];

  double s[9] = {sxx, sxy, sxz,
                 sxy, syy, syz,
                 sxz, syz, szz
                };

  return ComputeVonMises ( s );
}

PView* elasticitySolverXfem::buildDisplacementView ( const std::string &postFileName )
{

  SolverField<SVector3> Field ( _pAssembler, _FSpace );

  // store de displacement values <vertex,value>
  std::map<int, std::vector<double> > data;

  // for new elements
  std::vector<MVertex*> vecVertices;
  std::vector <MElement *> vecElements;
  std::map<int, std::vector<MElement*> > Cut_Elements;

  // copy of _pModel
  GModel *pModelPost = new GModel ( *_pModel );
  GModel::setCurrent ( pModelPost );

  // tolerance for discontinuity point calcul
  double eps = 1e-6;

  // entity associate to cut elements
  std::map<MElement*, GEntity*> ElementCutEntity;
  typedef std::set<GFace*, GEntityLessThan>::iterator fiter;
  for ( fiter it = pModelPost->firstFace(); it != pModelPost->lastFace(); ++it )
    for ( unsigned int i = 0; i < ( *it )->polygons.size(); i++ )
      ElementCutEntity[ ( *it )->polygons[i]] = ( *it );

  // REGION
  typedef std::set<GRegion*, GEntityLessThan>::iterator riter;
  for ( riter it = pModelPost->firstRegion(); it != pModelPost->lastRegion(); ++it )
    for ( unsigned int i = 0; i < ( *it )->polyhedra.size(); i++ )
      ElementCutEntity[ ( *it )->polyhedra[i]] = ( *it );


  for ( unsigned int i = 0; i < _elasticFields.size(); ++i )
  {
    for ( groupOfElements::elementContainer::const_iterator it = _elasticFields[i].g->begin(); it != _elasticFields[i].g->end(); ++it )
    {
      MElement *e=*it;
      MElement *ep=e->getParent(); // cut
      if ( !ep )
      {
        for ( int j = 0; j < e->getNumVertices(); ++j )
        {
          SVector3 val;
          SPoint3 p ( e->getVertex ( j )->x(),e->getVertex ( j )->y(),e->getVertex ( j )->z() );
          double pnt[3];
          pnt[0]=p.x();
          pnt[1]=p.y();
          pnt[2]=p.z();
          double uvw[3];
          e->xyz2uvw ( pnt,uvw );
          Field.f ( e,uvw[0],uvw[1],uvw[2],val );
          std::vector<double> vec ( 3 );
          vec[0]=val ( 0 );
          vec[1]=val ( 1 );
          vec[2]=val ( 2 );
          data[e->getVertex ( j )->getNum() ]=vec;
        }
      }
      else
      {
        // determine tag
        GEntity *g = ElementCutEntity[e];
        // take all the parts
        for ( int i = 0; i < e->getNumChildren() ; i ++ )
        {
          MElement *ec;
          ec = e->getChild ( i );  // get part
          std::vector<MVertex *> vertices;
          vertices.clear();
          for ( int j = 0; j < ec->getNumVertices(); ++j )
          {

            // duplicate vertices
            MVertex *v = new MVertex ( ec->getVertex ( j )->x(),ec->getVertex ( j )->y(),ec->getVertex ( j )->z(),g );

            vertices.push_back ( v );
            vecVertices.push_back ( v );

            SVector3 val;
            SPoint3 p ( ec->getVertex ( j )->x() , ec->getVertex ( j )->y(), ec->getVertex ( j )->z() );

            // if discontinuity point, 'limit calcul'
            if ( false )
            {
              if ( ( *_lslist[0] ) ( p.x() , p.y(), p.z() ) < eps ) // fix for all lslist[0]
              {
                SPoint3 b = ec->barycenter();
                b = b - p;
                b = b*eps;
                p = p + b;
              }
            }
            double pnt[3];
            pnt[0]=p.x();
            pnt[1]=p.y();
            pnt[2]=p.z();
            double uvw[3];
            // the parent owned the field
            ep->xyz2uvw ( pnt,uvw );
            Field.f ( ep,uvw[0],uvw[1],uvw[2],val );
            std::vector<double> vec ( 3 );
            vec[0]=val ( 0 );
            vec[1]=val ( 1 );
            vec[2]=val ( 2 );
            data[ec->getVertex ( j )->getNum() ]=vec;
            data[v->getNum() ]=vec;
          }

          // create elements associate with entities
          MElement *ecopy;
          if ( ec->getType() != 3 && ec->getType() != 4 && ec->getType() != 5 )
            std::cout<<"warning.. element type : " << ec->getType()  << std::endl;
          if ( ec->getType() == 3 ) ecopy = new MTriangle ( vertices );
          if ( ec->getType() == 4 ) ecopy = new MQuadrangle ( vertices );
          if ( ec->getType() == 5 ) ecopy = new MTetrahedron ( vertices );
          // fill the new elements to be created map
          Cut_Elements[g->tag() ].push_back ( ecopy );
        }
      }
    }
  }

  std::map < int, std::map<int, std::string > > physical; // empty

  // store new elements and vertices
  pModelPost->storeChain ( vecVertices,_dim, Cut_Elements, physical );

//  pModelPost->writeMSH("ModelCut.msh",2.1,false,false);

  PView *pv = new PView ( postFileName, "NodeData", pModelPost, data, 0.0 );
  return pv;
}

PView *elasticitySolverXfem::buildElasticEnergyView ( const std::string &postFileName )
{
  std::map<int, std::vector<double> > data;
  GaussQuadrature Integ_Bulk ( GaussQuadrature::GradGrad );
  for ( unsigned int i = 0; i < _elasticFields.size(); ++i )
  {
    SolverField<SVector3> Field ( _pAssembler, _FSpace );
    IsotropicElasticTerm Eterm ( Field,_elasticFields[i]._E,_elasticFields[i]._nu );
    BilinearTermToScalarTerm Elastic_Energy_Term ( Eterm );
    ScalarTermConstant<double> One ( 1.0 );
    for ( groupOfElements::elementContainer::const_iterator it = _elasticFields[i].g->begin(); it != _elasticFields[i].g->end(); ++it )
    {
      MElement *e=*it;
//      if ( e->getParent() ) e=e->getParent();
      double energ;
      double vol;
      IntPt *GP;
      int npts=Integ_Bulk.getIntPoints ( e,&GP );
      Elastic_Energy_Term.get ( e,npts,GP,energ );
      One.get ( e,npts,GP,vol );
      std::vector<double> vec;
      vec.push_back ( energ/vol );
      data[ ( e )->getNum() ]=vec;
    }
  }
  PView *pv = new PView ( postFileName, "ElementData", _pModel, data, 0.0 );
  return pv;
}



PView *elasticitySolverXfem::buildVonMisesView ( const std::string &postFileName )
{
  std::map<int, std::vector<double> > data;
  GaussQuadrature Integ_Bulk ( GaussQuadrature::GradGrad );
  SolverField<SVector3> Field ( _pAssembler, _FSpace );
  GradTerm<SVector3> GTerm(Field);
  GradTerm<SVector3> GTerm2(Field);
  const LinearTermBase<STensor3> &somme=GTerm+GTerm2+GTerm;
  const BilinearTermBase &bbb=GTerm|GTerm2;
  for ( unsigned int i = 0; i < _elasticFields.size(); ++i )
  {
    IsotropicElasticTerm Eterm ( Field,_elasticFields[i]._E,_elasticFields[i]._nu );
//    BilinearTermToScalarTerm Elastic_Energy_Term ( Eterm );
    const double lam = _elasticFields[i]._nu *_elasticFields[i]._E/((1.+_elasticFields[i]._nu )*(1.-2.*_elasticFields[i]._nu ));
    const double mu  = _elasticFields[i]._E/(2.*(1.+_elasticFields[i]._nu ));
    int ii, jj, kk, ll;
    STensor43 elast;
    for (ii = 0; ii < 3; ++ii)
    {
      for (jj = 0; jj < 3; ++jj)
      {
        for (kk = 0; kk < 3; ++kk)
        {
          for (ll = 0; ll < 3; ++ll)
          {
            elast(ii,jj,kk,ll) = lam * delta(ii,jj) * delta(kk,ll) + mu * (delta(ii,kk) * delta(jj,ll) + delta(ii,ll) * delta(jj,kk));
          }
        }
      }
    }
  
    ScalarTermConstant<STensor43> t(elast);
    BilinearTermContractWithLaw<STensor3> bbb2(GTerm,t,GTerm);

    for ( groupOfElements::elementContainer::const_iterator it = _elasticFields[i].g->begin(); it != _elasticFields[i].g->end(); ++it )
    {
      MElement *e=*it;
      fullMatrix< double> s;
      fullMatrix< double> ss;
      fullVector<STensor3> s2;
      fullMatrix< double> sss;
      IntPt *GP;
      int npts=Integ_Bulk.getIntPoints ( e,&GP );
  //    Elastic_Energy_Term.get ( e,npts,GP,energ );
      Eterm.get(e,npts,GP,s);
      GTerm.get(e,npts,GP,s2);
      bbb.get(e,npts,GP,ss);
      bbb2.get(e,npts,GP,sss);
      s.print();
      s2(0).print("");
      ss.print();
      sss.print();
//      std::vector<double> vec;
//      vec.push_back ( energ );
//      data[ ( e )->getNum() ]=vec;
    }
  }
  PView *pv = new PView ( postFileName, "ElementData", _pModel, data, 0.0 );
  return pv;
}


#else
PView* elasticitySolverXfem::buildDisplacementView ( const std::string &postFileName )
{
  Msg::Error ( "Post-pro module not available" );
  return 0;
}

PView* elasticitySolverXfem::buildElasticEnergyView ( const std::string &postFileName )
{
  Msg::Error ( "Post-pro module not available" );
  return 0;
}


PView* elasticitySolverXfem::buildVonMisesView ( const std::string &postFileName )
{
  Msg::Error ( "Post-pro module not available" );
  return 0;
}

#endif