// Reconstruction - A tool for building cad models from meshes
// Copyright (C) 2010-2026 Eric Bechet, Borhen Louhichi
//
// See the LICENSE file for license information and contributions.
// Please report all bugs and problems to <bechet@cadxfem.org>.

#include <string.h>
#include "GmshConfig.h"
#include "BrepXfem.h"

#include "GModel.h"
#include "gmshLevelset.h"
#include "math.h"

#include "RecPrimitive.h"
#include "RecBSplineLM.h"
#include "BSplineApproxSurface.h"

#include "outil.h"
#include "outil_spline.h"


std::vector<gp_Pnt > recuperer_Pnts_ls(GModel *_pModelUncut, gLevelset *_ls)
{

  std::vector<MElement *> Elements_ls;

  typedef std::set<GRegion*, GEntityLessThan>::iterator riter;

  for (riter region =_pModelUncut->firstRegion(); region != _pModelUncut->lastRegion(); region++)
  {

    unsigned int numElements[5];
    numElements[0]=0;
    numElements[1]=0;
    numElements[2]=0;
    numElements[3]=0;
    numElements[4]=0;

    (*region)->getNumMeshElements(numElements);

    for (int k=0; k<numElements[0]; k++)
    {

      MElement* tetra = (*region)->getMeshElement(k);

      int testmoin=0;
      int testplus=0;

      for (int v=0; v<tetra->getNumVertices(); v++)
      {

        MVertex* Vertex = tetra->getVertex(v);
        double valeur_level = (*_ls)(Vertex->x(),Vertex->y(),Vertex->z());

        if (valeur_level<=0)
        {
          testmoin++;

        }
        if (valeur_level>=0)
        {
          testplus++;

        }
      }

      if (testmoin>0 && testplus>0) Elements_ls.push_back(tetra);
    }
  }

  std::vector<gp_Pnt > Pnts;

  for (int i=0; i<Elements_ls.size(); i++)
  {

    MElement* tetra = Elements_ls[i];
    std::vector<MEdge > Edges;
    //TColgp_Array1OfPnt thePnts(1,3);

    for (int e=0; e<tetra->getNumEdges(); e++)
    {

      MEdge Edge = tetra->getEdge(e);

      MVertex* Vertex0 = Edge.getVertex(0);
      MVertex* Vertex1 = Edge.getVertex(1);

      double valeur_level0 = (*_ls)(Vertex0->x(),Vertex0->y(),Vertex0->z());
      double valeur_level1 = (*_ls)(Vertex1->x(),Vertex1->y(),Vertex1->z());


      if (valeur_level0*valeur_level1 <=1e-6)
      {
        Edges.push_back(Edge);

        double xx = ((std::abs(valeur_level0)) * Vertex1->x() + (std::abs(valeur_level1)) * Vertex0->x()) / ((std::abs(valeur_level0)) + (std::abs(valeur_level1)));
        double yy = ((std::abs(valeur_level0)) * Vertex1->y() + (std::abs(valeur_level1)) * Vertex0->y()) / ((std::abs(valeur_level0)) + (std::abs(valeur_level1)));
        double zz = ((std::abs(valeur_level0)) * Vertex1->z() + (std::abs(valeur_level1)) * Vertex0->z()) / ((std::abs(valeur_level0)) + (std::abs(valeur_level1)));

        gp_Pnt point(xx, yy, zz);

        int tt=0;
        for (int jj=0; jj<Pnts.size(); jj++)
        {
          if (point.IsEqual(Pnts[jj], 1e-6)) tt=1;

        }
        if (tt==0) Pnts.push_back(gp_Pnt(xx,yy,zz));

      }

    }

  }


  return Pnts;

}

std::vector<gp_Pnt > filtre_pts_arete(std::vector<gp_Pnt > ptsArete, int nbptsmax)
{

  std::vector<double > Parameters;

  double distance=0.0;
  Parameters.push_back(0.0);

  for (int i=1; i<ptsArete.size(); i++)
  {
    distance = distance + ptsArete[i].Distance(ptsArete[i-1]);
    Parameters.push_back(distance);
  }

  std::vector<gp_Pnt > ptsArete_filtree;

  double Uf = Parameters[0];
  double Ul = Parameters[Parameters.size()-1];

  double pas= (Ul-Uf) / (nbptsmax-1);
  double param = 0.0;
  int ini =0;
  ptsArete_filtree.push_back(ptsArete[0]);

  for (int j=0; j< (nbptsmax-1)-1; j++)
  {
    param = param + pas;

    do
    {
      ini = ini+1;


    }
    while (param>Parameters[ini]);

    if (ptsArete[ini-1].IsEqual(ptsArete_filtree[ptsArete_filtree.size()-1], 1e-6) ==0) ptsArete_filtree.push_back(ptsArete[ini-1]);

  }

  ptsArete_filtree.push_back(ptsArete[ptsArete.size()-1]);

  return ptsArete_filtree;

}



std::vector<gp_Pnt > filtre_pts_face(std::vector<gp_Pnt > ptsFace, int nbptsU, int nbptsV)
{



}

void BrepXfem::setMesh0(const std::string &meshFileName)
{
  FILE * pFile;
  pFile = fopen(meshFileName.c_str(),"rt");

  BRep_Builder B;
  TopoDS_Shell Sh;

  B.MakeShell(Sh);

  int nb_pts;
  fscanf(pFile,"%d\n", &nb_pts);
  std::cout<<" nb_pts "<< nb_pts <<std::endl;

  std::vector<gp_Pnt> points;

  double x,y,z;

  for (int p=0; p<nb_pts; p++)
  {
    fscanf(pFile, "%lf %lf %lf\n", &x, &y, &z);
    points.push_back(gp_Pnt(x, y, z));

  }

  int nb_aretes;
  fscanf(pFile,"%d\n", &nb_aretes);
  std::cout<<" nb_aretes "<< nb_aretes <<std::endl;

  std::vector<std::vector<gp_Pnt> > pointsAreteS;

  for (int a=0; a<nb_aretes; a++)
  {

    int ida, nbptsA;
    fscanf(pFile, " %d %d\n", &ida, &nbptsA);

    //std::cout<<" nbpts_aretes "<< nbptsA <<std::endl;

    std::vector<gp_Pnt> pointsArete;

    for (int pa=0; pa<nbptsA; pa++)
    {

      int num;
      fscanf(pFile, " %d", &num);
      gp_Pnt P = points[num];
      pointsArete.push_back(P);

    }
    fscanf(pFile, " \n");

    ///**************** Reconstruction des aretes NURBS BOOK
    int deg =2;
    std::vector<double > Knots;
    std::vector<gp_Pnt > CPnts;

    CurveSpline CS;

    CS.GlobalCurveInterp(pointsArete, deg, Knots, CPnts, 1);

    std::vector<double > KnotsOCC;
    std::vector<int > Multiples;

    Spline S;
    S.OCCKnot_multiplicity(Knots, KnotsOCC, Multiples);

    TColgp_Array1OfPnt Poles(1,CPnts.size());

    for (int i=0; i<=CPnts.size()-1; i++)
    {
      Poles.SetValue(i+1, CPnts[i]);

    }
    TColStd_Array1OfReal Knots0(1,KnotsOCC.size());
    TColStd_Array1OfInteger Multiplicities(1,Multiples.size());


    for (int i=0; i<=KnotsOCC.size()-1; i++)
    {
      Knots0.SetValue(i+1, KnotsOCC[i]);
      Multiplicities.SetValue(i+1, Multiples[i]);
    }

    Handle(Geom_BSplineCurve) SPL =  new Geom_BSplineCurve(Poles, Knots0, Multiplicities, deg);
    ///****************************************************************************************

    std::cout<<" arete "<< a <<std::endl;

    double t0 = SPL->FirstParameter();
    double t1 = SPL->LastParameter();

    std::vector<gp_Pnt > ptsArete_filtree1;

    int decoup = 30;
    for (int j=0; j<=decoup; j++)
    {
      double t=t0+j* ((t1-t0) /decoup);
      gp_Pnt P= SPL->Value(t);

      ptsArete_filtree1.push_back(P);

    }

    pointsAreteS.push_back(ptsArete_filtree1);


  }



  int nb_faces;
  fscanf(pFile,"%d\n", &nb_faces);
  std::cout<<" nb_faces "<< nb_faces <<std::endl;

  for (int f=0; f<nb_faces; f++)
  {
    int idf, nbtriF;
    fscanf(pFile, " %d %d\n", &idf, &nbtriF);

    std::cout<<" idf "<< idf <<std::endl;

    std::vector<gp_Pnt> pointsFace;

    for (int tr=0; tr<nbtriF; tr++)
    {
      double x, y, z;
      fscanf(pFile, "%lf %lf %lf\n", &x, &y, &z);
      pointsFace.push_back(gp_Pnt(x,y,z));
    }

    int nbAF;
    fscanf(pFile, " %d\n", &nbAF);
    std::cout<<" nbAF "<< nbAF <<std::endl;

    std::vector<std::vector<gp_Pnt> > pointsAreteFace;

    for (int af=0; af<nbAF; af++)
    {
      int idaf;
      fscanf(pFile, "%d ", &idaf);
      pointsAreteFace.push_back(pointsAreteS[idaf]);

    }
    fscanf(pFile, " \n");



    std::vector<std::vector<std::vector<gp_Pnt > > > Contour_Aretes_pts;

    int ii = 0;
    do
    {
      int CC=0;
      std::vector<std::vector<gp_Pnt > > Contour_i_Aretes_pts;

      gp_Pnt Pi_contour, Pf_;
      if (((pointsAreteFace[ii])[pointsAreteFace[ii].size()-1].IsEqual((pointsAreteFace[ii+1])[0],1e-6) ==true) ||
          ((pointsAreteFace[ii])[pointsAreteFace[ii].size()-1].IsEqual((pointsAreteFace[ii+1])[pointsAreteFace[ii+1].size()-1],1e-6) ==true))
      {
        Pi_contour= (pointsAreteFace[ii])[0];
        Pf_ = (pointsAreteFace[ii])[pointsAreteFace[ii].size()-1];

      }
      if (((pointsAreteFace[ii])[0].IsEqual((pointsAreteFace[ii+1])[0],1e-6) ==true) ||
          ((pointsAreteFace[ii])[0].IsEqual((pointsAreteFace[ii+1])[pointsAreteFace[ii+1].size()-1],1e-6) ==true))
      {
        Pi_contour= (pointsAreteFace[ii])[pointsAreteFace[ii].size()-1];
        Pf_ = (pointsAreteFace[ii])[0];

      }

      do
      {
        Contour_i_Aretes_pts.push_back(pointsAreteFace[ii]);

        if (CC>0)
        {
          if (Pf_.IsEqual((pointsAreteFace[ii])[0], 1e-6) ==true)
          {
            Pf_= (pointsAreteFace[ii])[pointsAreteFace[ii].size()-1];
          }

          else
          {
            Pf_= (pointsAreteFace[ii])[0];
          }


        }

        CC=CC+1;
        ii = ii+1;


      }
      while (Pf_.IsEqual(Pi_contour, 1e-6) != true);

      Contour_Aretes_pts.push_back(Contour_i_Aretes_pts);


    }
    while (ii<pointsAreteFace.size());


    if (f == 2)
    {
      std::vector<gp_Pnt > Pntss;
      for (int i=0; i<Contour_Aretes_pts.size(); i++)
      {
        for (int j=0; j<Contour_Aretes_pts[i].size(); j++)
        {
          for (int k=1; k< (Contour_Aretes_pts[i])[j].size()-1; k++)
          {
            Pntss.push_back(((Contour_Aretes_pts[i])[j])[k]);

          }
        }
      }
      Pntss.push_back(gp_Pnt(0.0,0.0,0.0));
      Pntss.push_back(gp_Pnt(20.0,0.0,0.0));
      Pntss.push_back(gp_Pnt(0.0,0.0,20.0));
      Pntss.push_back(gp_Pnt(20.0,0.0,20.0));

      std::cout<<" Pntss.size() : "<< Pntss.size() <<std::endl;

      AveragePlane AP(Pntss, Pntss.size(),0.001);

      Handle(Geom_Plane) Plan = AP.Plane();
      gp_Pln Pln1 = Plan->Pln();

      std::cout<<" Position de Plan : "<< Pln1.Location().X() <<"  "<<Pln1.Location().Y() <<"  "<< Pln1.Location().Z() <<std::endl;

      gp_Ax3 Axe = Pln1.Position();

      std::cout<<" Direction : "<< Axe.Direction().X() <<"  "<<Axe.Direction().Y() <<"  "<< Axe.Direction().Z() <<std::endl;

      std::cout<<" Direction X: "<< Axe.XDirection().X() <<"  "<<Axe.XDirection().Y() <<"  "<< Axe.XDirection().Z() << std::endl;
      std::cout<<" Direction Y: "<< Axe.YDirection().X() <<"  "<<Axe.YDirection().Y() <<"  "<< Axe.YDirection().Z() << std::endl;


      ///********************************************

      TColgp_HArray1OfPnt myPts(1, Pntss.size());

      for (int i=0; i<Pntss.size(); i++)
      {
        myPts.SetValue(i+1, Pntss[i]);
      }

      gp_Vec myOX;
      gp_Vec myOY;
      gp_Ax2 Axe1;
      Standard_Boolean IsSingular;
      GeomLib::AxeOfInertia(myPts.Array1(), Axe1, IsSingular, 0.001);

      myOX = Axe1.XDirection();
      myOY = Axe1.YDirection();

      gp_Vec OZ = Axe1.Direction();
      std::cout<<" Direction OZ : "<< Axe1.Direction().X() <<"  "<<Axe1.Direction().Y() <<"  "<< Axe1.Direction().Z() << std::endl;

      std::cout<<" Direction OZ X: "<< Axe1.XDirection().X() <<"  "<<Axe1.XDirection().Y() <<"  "<< Axe1.XDirection().Z() << std::endl;
      std::cout<<" Direction OZ Y: "<< Axe1.YDirection().X() <<"  "<<Axe1.YDirection().Y() <<"  "<< Axe1.YDirection().Z() << std::endl;

      std::cout<<" isSingular: "<< IsSingular << std::endl;

      gp_Pnt GB, myG;

      Standard_Integer i,nb=myPts.Length();
      GB.SetCoord(0.,0.,0.);
      for (i=1; i<=nb; i++)
      {
        GB.SetCoord(1, (GB.Coord(1) +myPts.Value(i).Coord(1)));
        GB.SetCoord(2, (GB.Coord(2) +myPts.Value(i).Coord(2)));
        GB.SetCoord(3, (GB.Coord(3) +myPts.Value(i).Coord(3)));
      }
      myG.SetCoord(1, (GB.Coord(1) /nb));
      myG.SetCoord(2, (GB.Coord(2) /nb));
      myG.SetCoord(3, (GB.Coord(3) /nb));

      std::cout<<" myG :      "<< myG.X() <<"         "<< myG.Y() <<"          "<< myG.Z() <<std::endl;


      ///********************************************

      double Umin, Umax, Vmin, Vmax;
      //AP.MinMaxBox(Umin, Umax, Vmin, Vmax);

      ElSLib::Parameters(Pln1,Pntss[0],Umax,Vmax);
      Umin=Umax;
      Vmin=Vmax;
      double U=0,V=0;
      for (int i=1; i<Pntss.size(); i++)
      {
        ElSLib::Parameters(Pln1,Pntss[i],U,V);
        if (Umax < U) Umax=U;
        if (Umin > U) Umin=U;
        if (Vmax < V) Vmax=V;
        if (Vmin > V) Vmin=V;
      }
      std::cout<<" Umin : "<< Umin <<" Umax : "<< Umax << " Vmin : "<< Vmin << " Vmax : "<< Vmax <<std::endl;




      ///**********
      gp_Vec Normale = gp_Vec(Axe.Direction().X(), Axe.Direction().Y(), Axe.Direction().Z());

      std::vector<gp_Pnt> Pnt_Interpol;
      for (int i=0; i<Pntss.size(); i++)
      {
        GeomAPI_ProjectPointOnSurf Proj(Pntss[i],Plan);
        double distance = Proj.LowerDistance();
        double U,V;
        Proj.LowerDistanceParameters(U,V);

        gp_Pnt P0 = Proj.NearestPoint();
        gp_Vec Vec(P0, Pntss[i]);
        double produit = Vec*Normale;

        double h;
        if (produit>=0) h=distance;
        else
        {
          h=-distance;
        }

        //std::cout<<" h : "<< h <<std::endl;

        Pnt_Interpol.push_back(gp_Pnt(U,h,V));
      }

      std::cout<<" pointsFace.size() "<< pointsFace.size() <<std::endl;

      for (int i=0; i<pointsFace.size(); i++)
      {
        GeomAPI_ProjectPointOnSurf Proj(pointsFace[i],Plan);
        double distance = Proj.LowerDistance();
        double U,V;
        Proj.LowerDistanceParameters(U,V);

        gp_Pnt P0 = Proj.NearestPoint();
        gp_Vec Vec(P0, pointsFace[i]);
        double produit = Vec*Normale;

        double h;
        if (produit>=0) h=distance;
        else
        {
          h=-distance;
        }

        //std::cout<<" h : "<< h <<std::endl;

        Pnt_Interpol.push_back(gp_Pnt(U,h,V));
      }



      std::vector<std::vector<gp_Pnt> > Pnt_C0;
      int nbptsU = 5;
      int nbptsV = 5;
      Tps tps;
      Pnt_C0 = tps.calc_tps(Pnt_Interpol, nbptsU-1, nbptsV-1, Umin-0.1, Umax+0.1, Vmin-0.2, Vmax+0.2);

      FILE *ff0= fopen("pts0.geo", "wt");
      int ii0=1;

      int GRID_W = nbptsU-1;
      int GRID_H = nbptsV-1;
      double pasU= (Umax-Umin) /GRID_W, pasV= (Vmax-Vmin) /GRID_H;

      for (int i=0; i<=GRID_W; ++i)
      {
        for (int j=0; j<=GRID_H; ++j)
        {
          double x=Umin + i*pasU;
          double z=Vmin + j*pasV;
          gp_Pnt PP = Plan->Value(x,z);
          fprintf(ff0,"%s %d %s %lf %s %lf %s %lf %s %lf %s \n","Point(",ii0,") = {",PP.X(),",",PP.Y(),",",PP.Z(),",",10.0,"};");
          ii0++;
        }
      }

      fclose(ff0);


      FILE *ff= fopen("pts1.geo", "wt");
      int ii=1;

      std::vector<std::vector<gp_Pnt> > Pnt_C;
      for (int i=0; i<Pnt_C0.size(); i++)
      {
        std::vector<gp_Pnt> Pts0 = Pnt_C0[i];
        std::vector<gp_Pnt> Pts;
        for (int j=0; j<Pts0.size(); j++)
        {
          gp_Pnt P0 = Plan->Value(Pts0[j].X(), Pts0[j].Z());
          double h = Pts0[j].Y();
          Pts.push_back(gp_Pnt((P0.X() +h*Normale.X()), (P0.Y() +h*Normale.Y()), (P0.Z() +h*Normale.Z())));

          gp_Pnt PP((P0.X() +h*Normale.X()), (P0.Y() +h*Normale.Y()), (P0.Z() +h*Normale.Z()));

          fprintf(ff,"%s %d %s %lf %s %lf %s %lf %s %lf %s \n","Point(",ii,") = {",PP.X(),",",PP.Y(),",",PP.Z(),",",10.0,"};");
          ii++;

          //std::cout<<" P0.X() : "<< P0.X() <<" P0.Y() : "<< P0.Y() <<" P0.Z() : "<< P0.Z() <<" h : "<< h <<std::endl;
          //std::cout<<" Pts[i].X() : "<< Pts[j].X() <<" Pts[i].Y() : "<< Pts[j].Y() <<" Pts[i].Z() : "<< Pts[j].Z() <<std::endl;
        }
        Pnt_C.push_back(Pts);
      }

      fclose(ff);

      std::vector<double > KnotsU;
      std::vector<double > KnotsV;
      std::vector<std::vector<gp_Pnt > > Pij;
      int degU = 2, degV = 2;
      SurfaceSpline SS;
      SS.GlobalSurfInterp(Pnt_C, degU, degV, Pij, KnotsU, KnotsV);

      FILE *ff2= fopen("pts2.geo", "wt");
      int ii2=1;

      TColgp_Array2OfPnt Points_BS0(1,nbptsU,1,nbptsV);
      for (int i=0; i<Pnt_C.size(); i++)
      {
        std::vector<gp_Pnt> Pts = Pij[i];
        for (int j=0; j<Pts.size(); j++)
        {
          Points_BS0.SetValue(i+1,j+1,Pts[j]);

          fprintf(ff2,"%s %d %s %lf %s %lf %s %lf %s %lf %s \n","Point(",ii2,") = {",Pts[j].X(),",",Pts[j].Y(),",",Pts[j].Z(),",",10.0,"};");
          ii2++;
        }
      }

      fclose(ff2);

      std::vector<double > KnotsOCC_U;
      std::vector<int > Multiples_U;

      std::vector<double > KnotsOCC_V;
      std::vector<int > Multiples_V;

      Spline S;
      S.OCCKnot_multiplicity(KnotsU, KnotsOCC_U, Multiples_U);
      S.OCCKnot_multiplicity(KnotsV, KnotsOCC_V, Multiples_V);

      TColStd_Array1OfReal Knots0U(1,KnotsOCC_U.size());
      TColStd_Array1OfInteger MultiplicitiesU(1,Multiples_U.size());
      for (int i=0; i<=KnotsOCC_U.size()-1; i++)
      {
        Knots0U.SetValue(i+1, KnotsOCC_U[i]);
        MultiplicitiesU.SetValue(i+1, Multiples_U[i]);
      }

      TColStd_Array1OfReal Knots0V(1,KnotsOCC_V.size());
      TColStd_Array1OfInteger MultiplicitiesV(1,Multiples_V.size());
      for (int i=0; i<=KnotsOCC_V.size()-1; i++)
      {
        Knots0V.SetValue(i+1, KnotsOCC_V[i]);
        MultiplicitiesV.SetValue(i+1, Multiples_V[i]);
      }

      BRep_Builder B;
      TopoDS_Shell Sh;
      B.MakeShell(Sh);

      Handle(Geom_BSplineSurface) Surf = new Geom_BSplineSurface(Points_BS0,Knots0U,Knots0V,MultiplicitiesU,MultiplicitiesV,degU,degV);

      TopoDS_Face F0 = BRepBuilderAPI_MakeFace(Pln1, Umin-0.1, Umax+0.1, Vmin-0.2, Vmax+0.2);
      B.Add(Sh,F0);

      TopoDS_Face F = BRepBuilderAPI_MakeFace(Surf);
      B.Add(Sh,F);



      TopoDS_Solid Sol;
      B.MakeSolid(Sol);
      B.Add(Sol,Sh);
      STEPControl_Writer writer;
      writer.Transfer(Sol, STEPControl_ManifoldSolidBrep);
      writer.Write("TPS05012011.step");




    }
  }



  /*
      TopoDS_Solid Sol;
      B.MakeSolid(Sol);
      B.Add(Sol,Sh);
      STEPControl_Writer writer;
      writer.Transfer(Sol, STEPControl_ManifoldSolidBrep);

      std::string ModelName = meshFileName.c_str() ;
      ModelName.erase(ModelName.find(".",0),5);
      ModelName.insert(ModelName.length(),".step");

      writer.Write(ModelName.c_str());
  */



}



void BrepXfem::setMesh1(const std::string &meshFileName)
{

  _pModel = new GModel();
  _pModel->readMSH(meshFileName.c_str());
  _dim = _pModel->getNumRegions() ? 3 : 2;

  _pModelUncut = new GModel();
  _pModelUncut=_pModel;


  double pt[3]={50.0,0.0,0.0};
  double dir[3]={0.0,0.0,1.0};

  int tag = 1;
  _ls_cut = new gLevelsetGenCylinder(pt,dir,40.0,tag);
  gLevelsets.push_back(_ls_cut);

  std::vector<gp_Pnt > Pnts = recuperer_Pnts_ls(_pModelUncut, _ls_cut);
  std::vector<gp_Pnt > Pntsbound;

  std::cout<<" Pnts.size() : "<< Pnts.size() <<std::endl;

  for (int i=0; i<Pnts.size(); i++)
  {
    if ((abs(Pnts[i].X()-90.0) <0.0001) || (abs(Pnts[i].X()-10.0) <=0.0001) || (abs(Pnts[i].Z()-100.0) <0.0001) || (abs(Pnts[i].Z()-0.0) <=0.0001))
    {
      Pntsbound.push_back(Pnts[i]);
      std::cout<<" Point : "<< i <<"    :  "<< (Pnts[i].X()-90.0) <<"  "<< Pnts[i].Y() <<"  "<< Pnts[i].Z() <<std::endl;
    }

  }

  for (int i=0; i<Pnts.size(); i++)
  {
    int test =0;
    for (int j=0; j<Pntsbound.size(); j++)
    {
      if (Pnts[i].Distance(Pntsbound[j]) <1e-6) test =1;

    }
    if (test ==0) Pntsbound.push_back(Pnts[i]);
  }

  std::cout<<" Pntsbound.size() : "<< Pntsbound.size() <<std::endl;

  AveragePlane AP(Pntsbound, Pntsbound.size(),0.001);

  Handle(Geom_Plane) Plan = AP.Plane();
  gp_Pln Pln1 = Plan->Pln();

  std::cout<<" Position de Plan : "<< Pln1.Location().X() <<"  "<<Pln1.Location().Y() <<"  "<< Pln1.Location().Z() <<std::endl;

  gp_Ax3 Axe = Pln1.Position();

  std::cout<<" Direction : "<< Axe.Direction().X() <<"  "<<Axe.Direction().Y() <<"  "<< Axe.Direction().Z() <<std::endl;

  double Umin, Umax, Vmin, Vmax;
  //AP.MinMaxBox(Umin, Umax, Vmin, Vmax);

  ElSLib::Parameters(Pln1,Pnts[0],Umax,Vmax);
  Umin=Umax;
  Vmin=Vmax;
  double U=0,V=0;
  for (int i=1; i<Pnts.size(); i++)
  {
    ElSLib::Parameters(Pln1,Pnts[i],U,V);
    if (Umax < U) Umax=U;
    if (Umin > U) Umin=U;
    if (Vmax < V) Vmax=V;
    if (Vmin > V) Vmin=V;
  }

  std::cout<<" Umin : "<< Umin <<" Umax : "<< Umax << " Vmin : "<< Vmin << " Vmax : "<< Vmax <<std::endl;

  gp_Vec Normale = gp_Vec(Axe.Direction().X(), Axe.Direction().Y(), Axe.Direction().Z());

  std::vector<gp_Pnt> Pnt_Interpol;

  for (int i=0; i<Pnts.size(); i++)
  {
    GeomAPI_ProjectPointOnSurf Proj(Pnts[i],Plan);
    double distance = Proj.LowerDistance();
    double U,V;
    Proj.LowerDistanceParameters(U,V);

    gp_Pnt P0 = Proj.NearestPoint();

    gp_Vec Vec(P0, Pnts[i]);

    double produit = Vec*Normale;

    double h;
    if (produit>=0) h=distance;
    else
    {
      h=-distance;
    }

    Pnt_Interpol.push_back(gp_Pnt(U,h,V));
  }

  std::vector<std::vector<gp_Pnt> > Pnt_C0;
  int nbptsU = 10;
  int nbptsV = 10;
  Tps tps;
  Pnt_C0 = tps.calc_tps(Pnt_Interpol, nbptsU-1, nbptsV-1, Umin, Umax, Vmin, Vmax);

  std::vector<std::vector<gp_Pnt> > Pnt_C;
  for (int i=0; i<Pnt_C0.size(); i++)
  {
    std::vector<gp_Pnt> Pts0 = Pnt_C0[i];
    std::vector<gp_Pnt> Pts;
    for (int j=0; j<Pts0.size(); j++)
    {
      gp_Pnt P0 = Plan->Value(Pts0[j].X(), Pts0[j].Z());
      double h = Pts0[j].Y();
      Pts.push_back(gp_Pnt((P0.X() +h*Normale.X()) /100, (P0.Y() +h*Normale.Y()) /100, (P0.Z() +h*Normale.Z()) /100));

      std::cout<<" P0.X() : "<< P0.X() <<" P0.Y() : "<< P0.Y() <<" P0.Z() : "<< P0.Z() <<std::endl;
      std::cout<<" Pts[i].X() : "<< Pts[j].X() <<" Pts[i].Y() : "<< Pts[j].Y() <<" Pts[i].Z() : "<< Pts[j].Z() <<std::endl;
    }
    Pnt_C.push_back(Pts);
  }


  std::vector<double > KnotsU;
  std::vector<double > KnotsV;
  std::vector<std::vector<gp_Pnt > > Pij;
  int degU = 2, degV = 2;
  SurfaceSpline SS;
  SS.GlobalSurfInterp(Pnt_C, degU, degV, Pij, KnotsU, KnotsV);

  TColgp_Array2OfPnt Points_BS0(1,nbptsU,1,nbptsV);
  for (int i=0; i<Pnt_C.size(); i++)
  {
    std::vector<gp_Pnt> Pts = Pij[i];
    for (int j=0; j<Pts.size(); j++)
    {
      Points_BS0.SetValue(i+1,j+1,Pts[j]);

    }

  }

  std::vector<double > KnotsOCC_U;
  std::vector<int > Multiples_U;

  std::vector<double > KnotsOCC_V;
  std::vector<int > Multiples_V;

  Spline S;
  S.OCCKnot_multiplicity(KnotsU, KnotsOCC_U, Multiples_U);
  S.OCCKnot_multiplicity(KnotsV, KnotsOCC_V, Multiples_V);

  TColStd_Array1OfReal Knots0U(1,KnotsOCC_U.size());
  TColStd_Array1OfInteger MultiplicitiesU(1,Multiples_U.size());
  for (int i=0; i<=KnotsOCC_U.size()-1; i++)
  {
    Knots0U.SetValue(i+1, KnotsOCC_U[i]);
    MultiplicitiesU.SetValue(i+1, Multiples_U[i]);
  }

  TColStd_Array1OfReal Knots0V(1,KnotsOCC_V.size());
  TColStd_Array1OfInteger MultiplicitiesV(1,Multiples_V.size());
  for (int i=0; i<=KnotsOCC_V.size()-1; i++)
  {
    Knots0V.SetValue(i+1, KnotsOCC_V[i]);
    MultiplicitiesV.SetValue(i+1, Multiples_V[i]);
  }

  BRep_Builder B;
  TopoDS_Shell Sh;
  B.MakeShell(Sh);

  Handle(Geom_BSplineSurface) Surf = new Geom_BSplineSurface(Points_BS0,Knots0U,Knots0V,MultiplicitiesU,MultiplicitiesV,degU,degV);

  TopoDS_Face F = BRepBuilderAPI_MakeFace(Surf);
  B.Add(Sh,F);

  TopoDS_Solid Sol;
  B.MakeSolid(Sol);
  B.Add(Sol,Sh);
  STEPControl_Writer writer;
  writer.Transfer(Sol, STEPControl_ManifoldSolidBrep);
  writer.Write("TPS15112010.step");



  for (int ils=0; ils<gLevelsets.size(); ils++)
  {
    _pModel=_pModel->buildCutGModel(gLevelsets[ils]);
  }

  std::string ModelName = meshFileName.c_str() ;
  ModelName.erase(ModelName.find(".",0),4);
  ModelName.insert(ModelName.length(),"_cutLS.msh");
  _pModel->writeMSH(ModelName,2.1,false,false);
  GModel::setCurrent(_pModel);



  /*

              gp_Ax2 Axe2(gp_Pnt(50.0, 50.0, -10), gp_Dir(gp_Vec(0.0, 0.0, 100.0)));


      double Rayon =40.0;
      TopoDS_Solid shape = BRepPrimAPI_MakeCylinder(Axe2,Rayon,120.0);

      TopExp_Explorer exp0;

      for(exp0.Init(shape, TopAbs_FACE); exp0.More(); exp0.Next())
      {

          TopoDS_Face Face = TopoDS::Face(exp0.Current());

          std::vector<double> ParamF;

          TopLoc_Location location;
          Handle (Geom_Surface) Surface = BRep_Tool::Surface(Face,location);

          double U1, U2, V1, V2;
          Surface->Bounds (U1, U2, V1, V2);

          if (Surface->DynamicType() == STANDARD_TYPE(Geom_CylindricalSurface))
          {

              std::cout<<" Cylindre : "<<std::endl;
              Handle(Geom_CylindricalSurface) aCylindre = Handle(Geom_CylindricalSurface)::DownCast(Surface);
              gp_Cylinder agpCylindre = aCylindre->Cylinder();

              Convert_CylinderToBSplineSurface Convert_Cyl_BSpline(agpCylindre, -10.0,110.0);

              std::cout<<" U1 : "<< U1 <<" V1 : "<< V1 <<" U2 "<< U2 <<" V2 "<< V2 <<std::endl;

              std::vector<gp_Pnt> Pnt_Interpol;

              for(int i=0; i<Pnts.size(); i++)
              {
                  GeomAPI_ProjectPointOnSurf Proj(Pnts[i],Surface);
                  double distance = Proj.LowerDistance();

                  double U,V;
                  Proj.LowerDistanceParameters(U,V);

                  gp_Pnt P0; gp_Vec D1U, D1V;
                  Surface->D1(U, V, P0, D1U, D1V);
                  gp_Vec Normale = D1U^D1V;

                  gp_Vec Vec(P0, Pnts[i]);

                  double produit = Vec*Normale;
                  //std::cout<<" U : "<< U  <<" Produit "<< produit <<std::endl;
                  double h;
                  if(produit>=0) h=distance;
                  else{h=-distance;}

                  Pnt_Interpol.push_back(gp_Pnt(Rayon*U,h,V));

              }

              std::vector<std::vector<gp_Pnt> > Pnt_C;

              int nbptsU = 10; int nbptsV = 10;

              Tps tps;

              Pnt_C = tps.calc_tps_Revol(Pnt_Interpol, nbptsU-1, nbptsV-1, 0.0, 2*PI, 0.0, 120, Rayon);

              TColgp_Array2OfPnt Points_BS(1,nbptsU,1,nbptsV);

              for(int i=0; i<=nbptsU-1; i++)
              {
                  std::vector<gp_Pnt> PntU=Pnt_C[i];
                  for(int j=0; j<=nbptsV-1; j++)
                  {
                      gp_Pnt P=PntU[j];

                      //std::cout<<" U : "<< P.X()  <<" V : "<< P.Z()  <<" Hauteur : "<< P.Y()  <<std::endl;
                      double h=P.Y();
                      gp_Pnt P0; gp_Vec D1U, D1V;
                      Surface->D1(P.X()/Rayon, P.Z(), P0, D1U, D1V);
                      gp_Vec Normale = D1U^D1V;
                      Normale.Normalize();

                      gp_Pnt Pcalcul(P0.X()+h*Normale.X(), P0.Y()+h*Normale.Y(), P0.Z()+h*Normale.Z());

                      Points_BS.SetValue(i+1,j+1,Pcalcul);

                      //std::cout<<" P0.X() : "<< P0.X()  <<" P0.Y() : "<< P0.Y()  <<" P0.Z() : "<< P0.Z()  <<std::endl;
                      std::cout<<" Pcalcul.X() : "<< Pcalcul.X()  <<" Pcalcul.Y() : "<< Pcalcul.Y()  <<" Pcalcul.Z() : "<< Pcalcul.Z()  <<std::endl;

                  }
              }



              for(int j=0; j<=nbptsV-1; j++)
              {
                  gp_Pnt P((Points_BS(1,j+1).X()+Points_BS(nbptsU,j+1).X())/2,(Points_BS(1,j+1).Y()+Points_BS(nbptsU,j+1).Y())/2,(Points_BS(1,j+1).Z()+Points_BS(nbptsU,j+1).Z())/2);

                  //gp_Pnt P(Points_BS(1,j+1).X(),Points_BS(1,j+1).Y(),Points_BS(1,j+1).Z());

                  //gp_Pnt P1()

                  Points_BS.SetValue(1,j+1,P);
                  Points_BS.SetValue(nbptsU,j+1,P);



              }

              Handle(Geom_BSplineSurface) Surf = GeomAPI_PointsToBSplineSurface(Points_BS).Surface();


              BRep_Builder B;
              TopoDS_Shell Sh;
              B.MakeShell(Sh);

              TopoDS_Face F = BRepBuilderAPI_MakeFace(Surf);
              B.Add(Sh,F);

              TopoDS_Solid Sol;
              B.MakeSolid(Sol);
              B.Add(Sol,Sh);
              STEPControl_Writer writer;
              writer.Transfer(Sol, STEPControl_ManifoldSolidBrep);
              writer.Write("TPS.step");

          }

      }
  */



}

void BrepXfem::setMesh2(const std::string &meshFileName)
{

  std::vector<gp_Pnt> Pnt_Interpol;

  Pnt_Interpol.push_back(gp_Pnt(0.0,-5.0,0.0));
  Pnt_Interpol.push_back(gp_Pnt(10.0,-5.0,0.0));
  Pnt_Interpol.push_back(gp_Pnt(10.0,-5.0,10.0));
  Pnt_Interpol.push_back(gp_Pnt(0.0,-5.0,10.0));

  Pnt_Interpol.push_back(gp_Pnt(2.0,3.0,2.0));
  Pnt_Interpol.push_back(gp_Pnt(8.0,3.0,2.0));
  Pnt_Interpol.push_back(gp_Pnt(8.0,3.0,8.0));
  Pnt_Interpol.push_back(gp_Pnt(2.0,3.0,8.0));

  Pnt_Interpol.push_back(gp_Pnt(5.0,-5.0,5.0));



  std::vector<std::vector<gp_Pnt> > Pnt_C;
  int nbptsU = 5;
  int nbptsV = 5;
  Tps tps;
  Pnt_C = tps.calc_tps(Pnt_Interpol, nbptsU-1, nbptsV-1, 0.0, 10.0, 0.0, 10.0);

  TColgp_Array2OfPnt Points_BS(1,nbptsU,1,nbptsV);

  for (int i=0; i<Pnt_C.size(); i++)
  {
    std::vector<gp_Pnt> Pts = Pnt_C[i];
    for (int j=0; j<Pts.size(); j++)
    {
      Points_BS.SetValue(i+1,j+1,Pts[j]);
      std::cout<<" Point "<< i <<"  "<< j <<" : "<< Pts[j].X() <<"  "<< Pts[j].Y() <<"  "<< Pts[j].Z() <<std::endl;
    }

  }

  std::vector<double > KnotsU;
  std::vector<double > KnotsV;
  std::vector<std::vector<gp_Pnt > > Pij;
  int degU = 2, degV = 2;
  SurfaceSpline SS;
  SS.GlobalSurfInterp(Pnt_C, degU, degV, Pij, KnotsU, KnotsV);

  TColgp_Array2OfPnt Points_BS0(1,nbptsU,1,nbptsV);
  for (int i=0; i<Pnt_C.size(); i++)
  {
    std::vector<gp_Pnt> Pts = Pij[i];
    for (int j=0; j<Pts.size(); j++)
    {
      Points_BS0.SetValue(i+1,j+1,Pts[j]);
      std::cout<<" Point "<< i <<"  "<< j <<" : "<< Pts[j].X() <<"  "<< Pts[j].Y() <<"  "<< Pts[j].Z() <<std::endl;
    }

  }

  std::vector<double > KnotsOCC_U;
  std::vector<int > Multiples_U;

  std::vector<double > KnotsOCC_V;
  std::vector<int > Multiples_V;

  Spline S;
  S.OCCKnot_multiplicity(KnotsU, KnotsOCC_U, Multiples_U);
  S.OCCKnot_multiplicity(KnotsV, KnotsOCC_V, Multiples_V);

  TColStd_Array1OfReal Knots0U(1,KnotsOCC_U.size());
  TColStd_Array1OfInteger MultiplicitiesU(1,Multiples_U.size());
  for (int i=0; i<=KnotsOCC_U.size()-1; i++)
  {
    Knots0U.SetValue(i+1, KnotsOCC_U[i]);
    MultiplicitiesU.SetValue(i+1, Multiples_U[i]);
  }

  TColStd_Array1OfReal Knots0V(1,KnotsOCC_V.size());
  TColStd_Array1OfInteger MultiplicitiesV(1,Multiples_V.size());
  for (int i=0; i<=KnotsOCC_V.size()-1; i++)
  {
    Knots0V.SetValue(i+1, KnotsOCC_V[i]);
    MultiplicitiesV.SetValue(i+1, Multiples_V[i]);
  }

  BRep_Builder B;
  TopoDS_Shell Sh;
  B.MakeShell(Sh);

  Handle(Geom_BSplineSurface) Surf = new Geom_BSplineSurface(Points_BS0,Knots0U,Knots0V,MultiplicitiesU,MultiplicitiesV,degU,degV);

  TopoDS_Face F = BRepBuilderAPI_MakeFace(Surf);
  B.Add(Sh,F);

  TopoDS_Solid Sol;
  B.MakeSolid(Sol);
  B.Add(Sol,Sh);
  STEPControl_Writer writer;
  writer.Transfer(Sol, STEPControl_ManifoldSolidBrep);
  writer.Write("TPS.step");


}





void BrepXfem::setMesh(const std::string &meshFileName)
{
  _pModel = new GModel();
  _pModel->readMSH(meshFileName.c_str());
  _dim = _pModel->getNumRegions() ? 3 : 2;

  _pModelUncut = new GModel();
  _pModelUncut=_pModel;


  FILE * pFile;
  pFile = fopen(meshFileName.c_str(),"rt");

  BRep_Builder B;
  TopoDS_Shell Sh;

  B.MakeShell(Sh);

  int nb_pts;
  fscanf(pFile,"%d\n", &nb_pts);
  std::cout<<" nb_pts "<< nb_pts <<std::endl;

  std::vector<gp_Pnt> points;

  double x,y,z;

  for (int p=0; p<nb_pts; p++)
  {
    fscanf(pFile, "%lf %lf %lf\n", &x, &y, &z);
    points.push_back(gp_Pnt(x, y, z));

  }



  ///********************** paraview
  /*char filename[100];
  sprintf(filename, "Triangulation.vtk");
  //tprintf("Writing %s ....\n", filename);
  FILE *fp = fopen (filename, "w");
  fprintf (fp, "# vtk DataFile Version 2.0\n");
  fprintf (fp, "Really cool data\n");
  fprintf (fp, "ASCII\n");
  fprintf (fp, "DATASET POLYDATA\n");
  fprintf (fp, "POINTS %d double\n", (int)nb_pts);

      for(int j=0; j<nb_pts; j++)
      {
          fprintf (fp, "%lf %lf %lf\n", points[j].X(), points[j].Y(), points[j].Z());
      }

  int nblines=0;*/





  int nb_aretes;
  fscanf(pFile,"%d\n", &nb_aretes);
  std::cout<<" nb_aretes "<< nb_aretes <<std::endl;

  std::vector<std::vector<gp_Pnt> > pointsAreteS;

  std::map<int, gp_Pnt> mapVertex;

  std::vector<std::vector<int > > num_Vertex_edge;

  for (int a=0; a<nb_aretes; a++)
  {

    int ida, nbptsA;
    fscanf(pFile, " %d %d\n", &ida, &nbptsA);

    //std::cout<<" nbpts_aretes "<< nbptsA <<std::endl;
    std::vector<int > num_Vertex_e;
    std::vector<gp_Pnt> pointsArete;

    for (int pa=0; pa<nbptsA; pa++)
    {

      int num;
      fscanf(pFile, " %d", &num);
      gp_Pnt P = points[num];
      pointsArete.push_back(P);
      if (pa==0 || pa==nbptsA-1)
      {
        mapVertex.insert(std::pair<int, gp_Pnt> (num,P));
        num_Vertex_e.push_back(num);
        std::cout<<"  num      "<< num <<std::endl;
      }


    }
    fscanf(pFile, " \n");

    num_Vertex_edge.push_back(num_Vertex_e);

    /*

    std::vector<gp_Pnt > ptsArete_filtree;

    //if(pointsArete.size()>10)
    {
        ptsArete_filtree = filtre_pts_arete(pointsArete, 3);
    }
    //else
    {
       // ptsArete_filtree = pointsArete;
    }
    //pointsAreteS.push_back(ptsArete_filtree);

    //std::cout<<" ptsArete_filtree.size() "<< ptsArete_filtree.size() <<std::endl;

        TColgp_Array1OfPnt array(1,pointsArete.size());
        for(int j=0; j<pointsArete.size(); j++)
        {
            array.SetValue(j+1,pointsArete[j]);

        }
        Handle(Geom_BSplineCurve) SPL = GeomAPI_PointsToBSpline(array,2,3,GeomAbs_C0,1.0e-3).Curve();
    */

    ///**************** Reconstruction des aretes NURBS BOOK
    int deg =2;
    std::vector<double > Knots;
    std::vector<gp_Pnt > CPnts;

    CurveSpline CS;

    CS.GlobalCurveInterp(pointsArete, deg, Knots, CPnts, 1);

    std::vector<double > KnotsOCC;
    std::vector<int > Multiples;

    Spline S;
    S.OCCKnot_multiplicity(Knots, KnotsOCC, Multiples);

    TColgp_Array1OfPnt Poles(1,CPnts.size());

    for (int i=0; i<=CPnts.size()-1; i++)
    {
      Poles.SetValue(i+1, CPnts[i]);

    }
    TColStd_Array1OfReal Knots0(1,KnotsOCC.size());
    TColStd_Array1OfInteger Multiplicities(1,Multiples.size());


    for (int i=0; i<=KnotsOCC.size()-1; i++)
    {
      Knots0.SetValue(i+1, KnotsOCC[i]);
      Multiplicities.SetValue(i+1, Multiples[i]);
    }

    Handle(Geom_BSplineCurve) SPL =  new Geom_BSplineCurve(Poles, Knots0, Multiplicities, deg);
    ///****************************************************************************************

    std::cout<<" arete "<< a <<std::endl;

    double t0 = SPL->FirstParameter();
    double t1 = SPL->LastParameter();

    std::vector<gp_Pnt > ptsArete_filtree1;

    int decoup = 10;
    for (int j=0; j<=decoup; j++)
    {
      double t=t0+j* ((t1-t0) /decoup);
      gp_Pnt P= SPL->Value(t);

      ptsArete_filtree1.push_back(P);

    }

    pointsAreteS.push_back(ptsArete_filtree1);


  }

  double precision = Precision::Confusion();
  BRep_Builder B1;

  std::vector<TopoDS_Vertex > Vertices(mapVertex.size());

  std::map<int, int> indVertex;

  int iV=0;
  std::map<int,gp_Pnt>::iterator it;

  for (it=mapVertex.begin() ; it != mapVertex.end(); it++)
  {
    std::cout<<"  (*it).first "<< (*it).first <<std::endl;
    std::cout<<"  (*it).second "<< ((*it).second).X() <<"   "<< ((*it).second).Y() <<"  "<< ((*it).second).Z() <<std::endl;
  }

  for (it=mapVertex.begin() ; it != mapVertex.end(); it++)
  {
    B1.MakeVertex(Vertices[iV], (*it).second,precision);
    indVertex.insert(std::pair<int, int> ((*it).first,iV));
    iV++;

    std::cout<<"  (*it).first "<< (*it).first <<"  iV "<< iV <<std::endl;
  }

  std::vector<TopoDS_Edge > Edges(pointsAreteS.size());

  Handle(Geom_BSplineCurve) SPL;
  for (int i=0; i<pointsAreteS.size(); i++)
  {
    TColgp_Array1OfPnt array(1,pointsAreteS[i].size());
    for (int j=0; j<pointsAreteS[i].size(); j++)
    {
      array.SetValue(j+1, (pointsAreteS[i])[j]);
    }
    SPL = GeomAPI_PointsToBSpline(array,2,3,GeomAbs_C1,0).Curve();
    B1.MakeEdge(Edges[i],SPL,precision);

    Vertices[indVertex[(num_Vertex_edge[i])[0]]].Orientation(TopAbs_FORWARD);
    Vertices[indVertex[(num_Vertex_edge[i])[1]]].Orientation(TopAbs_REVERSED);

    B1.Add(Edges[i],Vertices[indVertex[(num_Vertex_edge[i])[0]]]);
    B1.Add(Edges[i],Vertices[indVertex[(num_Vertex_edge[i])[1]]]);

    std::vector<int > num_Vertex_e = num_Vertex_edge[i];

    std::cout<<num_Vertex_e.size() <<"  (num_Vertex_edge[i])[0] "<< num_Vertex_e[0] <<" (num_Vertex_edge[i])[1] "<< num_Vertex_e[1] <<std::endl;
  }

  int nb_faces;
  fscanf(pFile,"%d\n", &nb_faces);
  std::cout<<" nb_faces "<< nb_faces <<std::endl;

  std::vector<TopoDS_Face > Faces(nb_faces);
  TopoDS_Wire W10;

  int indicefaceReconstruite[nb_faces];

  for (int f=0; f<nb_faces; f++)
  {
    indicefaceReconstruite[f]=0;
    int idf, nbtriF;
    fscanf(pFile, " %d %d\n", &idf, &nbtriF);

    std::cout<<" idf "<< idf <<std::endl;

    std::vector<gp_Pnt> pointsFace;

    for (int tr=0; tr<nbtriF; tr++)
    {
      int p1,p2,p3;
      fscanf(pFile, "%d %d %d\n", &p1, &p2, &p3);

      pointsFace.push_back(points[p1]);
      pointsFace.push_back(points[p2]);
      pointsFace.push_back(points[p3]);

      /// paraview
      /*if(f!=19 && f!=20 && f!=27 && f!=32)
      {
      fprintf (fp, "%d %d %d\n", 2, p1, p2);
      fprintf (fp, "%d %d %d\n", 2, p2, p3);
      fprintf (fp, "%d %d %d\n", 2, p3, p1);
      nblines=nblines+3;
      }*/
    }

    int nbAF;
    fscanf(pFile, " %d\n", &nbAF);
    std::cout<<" nbAF "<< nbAF <<std::endl;

    std::vector<std::vector<gp_Pnt> > pointsAreteFace;
    std::vector<int > id_AreteS_face;

    for (int af=0; af<nbAF; af++)
    {
      int idaf;
      fscanf(pFile, "%d ", &idaf);
      pointsAreteFace.push_back(pointsAreteS[idaf]);

      id_AreteS_face.push_back(idaf);

    }
    fscanf(pFile, " \n");




    ///********** vtk (paraview) ************///
    /*
            /// paraview
        char filename1[100];
        sprintf(filename1, "Surface_%d.vtk", f);
        //tprintf("Writing %s ....\n", filename);
        FILE *fp1 = fopen (filename1, "w");
        fprintf (fp1, "# vtk DataFile Version 2.0\n");
        fprintf (fp1, "Really cool data\n");
        fprintf (fp1, "ASCII\n");
        fprintf (fp1, "DATASET POLYDATA\n");
        fprintf (fp1, "POINTS %d double\n", (int)pointsFace.size());

        for(int i=0; i<pointsFace.size(); i++)
        {
            fprintf (fp1, "%lf %lf %lf\n", pointsFace[i].X(), pointsFace[i].Y(), pointsFace[i].Z());

        }

        fprintf (fp1, "LINES %d %d\n", (int)pointsFace.size(), (int)pointsFace.size()*3);

        int iii=0;
        for(int i=0; i<pointsFace.size()/3; i++)
        {
            fprintf (fp1, "%d %d %d\n", 2, iii, iii+1);
            fprintf (fp1, "%d %d %d\n", 2, iii+1, iii+2);
            fprintf (fp1, "%d %d %d\n", 2, iii+2, iii);
            iii=iii+3;
        }
    */

    /*
            int nbptParaview=0;
            for(int i=0; i<pointsAreteFace.size(); i++) nbptParaview = nbptParaview + pointsAreteFace[i].size();

        /// paraview
            char filename1[100];
        sprintf(filename1, "Aretes_%d.vtk", f);
        //tprintf("Writing %s ....\n", filename);
        FILE *fp1 = fopen (filename1, "w");
        fprintf (fp1, "# vtk DataFile Version 2.0\n");
        fprintf (fp1, "Really cool data\n");
        fprintf (fp1, "ASCII\n");
        fprintf (fp1, "DATASET POLYDATA\n");
        fprintf (fp1, "POINTS %d double\n", (int)nbptParaview);


        for(int i=0; i<pointsAreteFace.size(); i++)
        {

            for(int j=0; j<pointsAreteFace[i].size(); j++)
            {
                fprintf (fp1, "%lf %lf %lf\n", (pointsAreteFace[i])[j].X(), (pointsAreteFace[i])[j].Y(), (pointsAreteFace[i])[j].Z());

            }

        }
        int nb_lines = nbptParaview-pointsAreteFace.size();

        fprintf (fp1, "LINES %d %d\n", (int)nb_lines, (int)nb_lines*3);

        int cc=0;
        for(int i=0; i<pointsAreteFace.size(); i++)
        {
            for(int j=0; j<pointsAreteFace[i].size()-1; j++)
            {
                fprintf (fp1, "%d %d %d\n", 2, cc, cc+1);
                cc=cc+1;

            }
            cc=cc+1;

        }
    */
///*************************************************************************************************************

    std::vector<std::vector<std::vector<gp_Pnt > > > Contour_Aretes_pts;
    std::vector<std::vector<int> > id_arete_contour;



    int ii = 0;
    do
    {
      int CC=0;
      std::vector<std::vector<gp_Pnt > > Contour_i_Aretes_pts;

      std::vector<int > id_arete_i_contour;

      gp_Pnt Pi_contour, Pf_;
      if (((pointsAreteFace[ii])[pointsAreteFace[ii].size()-1].IsEqual((pointsAreteFace[ii+1])[0],1e-6) ==true) ||
          ((pointsAreteFace[ii])[pointsAreteFace[ii].size()-1].IsEqual((pointsAreteFace[ii+1])[pointsAreteFace[ii+1].size()-1],1e-6) ==true))
      {
        Pi_contour= (pointsAreteFace[ii])[0];
        Pf_ = (pointsAreteFace[ii])[pointsAreteFace[ii].size()-1];

      }
      if (((pointsAreteFace[ii])[0].IsEqual((pointsAreteFace[ii+1])[0],1e-6) ==true) ||
          ((pointsAreteFace[ii])[0].IsEqual((pointsAreteFace[ii+1])[pointsAreteFace[ii+1].size()-1],1e-6) ==true))
      {
        Pi_contour= (pointsAreteFace[ii])[pointsAreteFace[ii].size()-1];
        Pf_ = (pointsAreteFace[ii])[0];
      }

      do
      {
        Contour_i_Aretes_pts.push_back(pointsAreteFace[ii]);
        id_arete_i_contour.push_back(id_AreteS_face[ii]);


        if (CC>0)
        {
          if (Pf_.IsEqual((pointsAreteFace[ii])[0], 1e-6) ==true)
          {
            Pf_= (pointsAreteFace[ii])[pointsAreteFace[ii].size()-1];
          }

          else
          {
            Pf_= (pointsAreteFace[ii])[0];
          }
        }

        CC=CC+1;
        ii = ii+1;

      }
      while (Pf_.IsEqual(Pi_contour, 1e-6) != true);


      Contour_Aretes_pts.push_back(Contour_i_Aretes_pts);

      id_arete_contour.push_back(id_arete_i_contour);

    }
    while (ii<pointsAreteFace.size());





///********************************************************************************************************************************
    /// Reconstruction des Surfaces en utilisant le contour exterieur (3 ou 4 aretes)
///********************************************************************************************************************************

    if (Contour_Aretes_pts.size() ==1 && (nbAF==4 || nbAF==3))
    {
      indicefaceReconstruite[f]=1;

      Handle(Geom_Surface) Surface;

      TColGeom_Array1OfBSplineCurve ArrayCurve(1,nbAF);

      if (nbAF==4)
      {
        //TColGeom_Array1OfBSplineCurve ArrayCurve(1,nbAF);
        for (int i=0; i<nbAF; i++)
        {
          ///**************** Reconstruction des aretes NURBS BOOK
          /*int deg =2;
          std::vector<double > Knots;
          std::vector<gp_Pnt > CPnts;

          CurveSpline CS;

          CS.GlobalCurveInterp(pointsAreteFace[i], deg, Knots, CPnts, 1);

          std::vector<double > KnotsOCC;
          std::vector<int > Multiples;

          Spline S;
          S.OCCKnot_multiplicity(Knots, KnotsOCC, Multiples);

          TColgp_Array1OfPnt Poles(1,CPnts.size());

          for(int i=0; i<=CPnts.size()-1; i++)
          {
              Poles.SetValue(i+1, CPnts[i]);

          }
          TColStd_Array1OfReal Knots0(1,KnotsOCC.size());
          TColStd_Array1OfInteger Multiplicities(1,Multiples.size());


          for(int i=0; i<=KnotsOCC.size()-1; i++)
          {
              Knots0.SetValue(i+1, KnotsOCC[i]);
              Multiplicities.SetValue(i+1, Multiples[i]);
          }

          Handle(Geom_BSplineCurve) SPL =  new Geom_BSplineCurve(Poles, Knots0, Multiplicities, deg);*/
          ///****************************************************************************************

          TColgp_Array1OfPnt array(1,pointsAreteFace[i].size());
          for (int j=0; j<pointsAreteFace[i].size(); j++)
          {
            array.SetValue(j+1, (pointsAreteFace[i])[j]);
          }
          Handle(Geom_BSplineCurve) SPL = GeomAPI_PointsToBSpline(array,2,3,GeomAbs_C1,0).Curve();

          ArrayCurve.SetValue(i+1, SPL);
        }

        Handle(GeomAdaptor_HCurve) SPL1Adaptor = new GeomAdaptor_HCurve(ArrayCurve(1));
        //Handle(GeomFill_SimpleBound) B1 = new GeomFill_SimpleBound(SPL1Adaptor,0.0,0.0);
        Handle(GeomFill_SimpleBound) B1 = new GeomFill_SimpleBound(SPL1Adaptor,Precision::Approximation(),Precision::Angular());
        Handle(GeomAdaptor_HCurve) SPL2Adaptor = new GeomAdaptor_HCurve(ArrayCurve(2));
        //Handle(GeomFill_SimpleBound) B2 = new GeomFill_SimpleBound(SPL2Adaptor,0.0,0.0);
        Handle(GeomFill_SimpleBound) B2 = new GeomFill_SimpleBound(SPL2Adaptor,Precision::Approximation(),Precision::Angular());
        Handle(GeomAdaptor_HCurve) SPL3Adaptor = new GeomAdaptor_HCurve(ArrayCurve(3));
        //Handle(GeomFill_SimpleBound) B3 = new GeomFill_SimpleBound(SPL3Adaptor,0.0,0.0);
        Handle(GeomFill_SimpleBound) B3 = new GeomFill_SimpleBound(SPL3Adaptor,Precision::Approximation(),Precision::Angular());
        Handle(GeomAdaptor_HCurve) SPL4Adaptor = new GeomAdaptor_HCurve(ArrayCurve(4));
        //Handle(GeomFill_SimpleBound) B4 = new GeomFill_SimpleBound(SPL4Adaptor,0.0,0.0);
        Handle(GeomFill_SimpleBound) B4 = new GeomFill_SimpleBound(SPL4Adaptor,Precision::Approximation(),Precision::Angular());

        Standard_Boolean NoCheck= Standard_False;
        Standard_Integer MaxDeg = 8;
        Standard_Integer MaxSeg = 8;

        GeomFill_ConstrainedFilling aConstrained(MaxDeg, MaxSeg);
        aConstrained.Init(B1,B2,B3,B4, Standard_False);
        Surface = aConstrained.Surface();
      }


      if (nbAF==3)
      {
        //TColGeom_Array1OfBSplineCurve ArrayCurve(1,nbAF);
        for (int i=0; i<nbAF; i++)
        {
          ///**************** Reconstruction des aretes NURBS BOOK
          /*int deg =2;
          std::vector<double > Knots;
          std::vector<gp_Pnt > CPnts;

          CurveSpline CS;

          CS.GlobalCurveInterp(pointsAreteFace[i], deg, Knots, CPnts, 1);

          std::vector<double > KnotsOCC;
          std::vector<int > Multiples;

          Spline S;
          S.OCCKnot_multiplicity(Knots, KnotsOCC, Multiples);

          TColgp_Array1OfPnt Poles(1,CPnts.size());

          for(int i=0; i<=CPnts.size()-1; i++)
          {
              Poles.SetValue(i+1, CPnts[i]);

          }
          TColStd_Array1OfReal Knots0(1,KnotsOCC.size());
          TColStd_Array1OfInteger Multiplicities(1,Multiples.size());

          for(int i=0; i<=KnotsOCC.size()-1; i++)
          {
              Knots0.SetValue(i+1, KnotsOCC[i]);
              Multiplicities.SetValue(i+1, Multiples[i]);
          }

          Handle(Geom_BSplineCurve) SPL =  new Geom_BSplineCurve(Poles, Knots0, Multiplicities, deg);*/
          ///****************************************************************************************

          TColgp_Array1OfPnt array(1,pointsAreteFace[i].size());
          for (int j=0; j<pointsAreteFace[i].size(); j++)
          {
            array.SetValue(j+1, (pointsAreteFace[i])[j]);
          }
          Handle(Geom_BSplineCurve) SPL = GeomAPI_PointsToBSpline(array,2,3,GeomAbs_C1,0).Curve();

          ArrayCurve.SetValue(i+1, SPL);
        }

        Handle(GeomAdaptor_HCurve) SPL1Adaptor = new GeomAdaptor_HCurve(ArrayCurve(1));
        Handle(GeomFill_SimpleBound) B1 = new GeomFill_SimpleBound(SPL1Adaptor,Precision::Approximation(),Precision::Angular());
        Handle(GeomAdaptor_HCurve) SPL2Adaptor = new GeomAdaptor_HCurve(ArrayCurve(2));
        Handle(GeomFill_SimpleBound) B2 = new GeomFill_SimpleBound(SPL2Adaptor,Precision::Approximation(),Precision::Angular());
        Handle(GeomAdaptor_HCurve) SPL3Adaptor = new GeomAdaptor_HCurve(ArrayCurve(3));
        Handle(GeomFill_SimpleBound) B3 = new GeomFill_SimpleBound(SPL3Adaptor,Precision::Approximation(),Precision::Angular());

        Standard_Boolean NoCheck= Standard_False;
        Standard_Integer MaxDeg = 8;
        Standard_Integer MaxSeg = 8;

        GeomFill_ConstrainedFilling aConstrained(MaxDeg, MaxSeg);
        aConstrained.Init(B1,B2,B3, Standard_False);
        Surface = aConstrained.Surface();
      }



      TopoDS_Face F1 = BRepBuilderAPI_MakeFace(Surface);
      B.Add(Sh,F1);


      B1.MakeFace(Faces[f],Surface,precision);
      B1.MakeWire(W10);


      //std::cout<<" ************************************** Contour_Aretes_pts.size() :   "<< Contour_Aretes_pts.size() <<std::endl;

      for (int i=0; i<Contour_Aretes_pts.size(); i++)
      {

        std::vector<std::vector<gp_Pnt > > Aretes_pts = Contour_Aretes_pts[i];

        B1.MakeWire(W10);

        for (int k=0; k<Aretes_pts.size(); k++)
        {

          std::vector<gp_Pnt > pts = Aretes_pts[k];
          TColgp_Array1OfPnt2d Array01(1,pts.size());

          for (int j=0; j<pts.size(); j++)
          {
            Standard_Real U,V;
            GeomAPI_ProjectPointOnSurf Proj(pts[j],Surface);
            Proj.LowerDistanceParameters(U, V);
            Array01.SetValue(j+1,gp_Pnt2d(U,V));

            //std::cout<<" U "<< U <<" V "<< V <<std::endl;
          }

          Handle(Geom2d_BSplineCurve) SSPL1 = Geom2dAPI_PointsToBSpline(Array01,2,3,GeomAbs_C1,0).Curve();

          TopoDS_Edge E = BRepBuilderAPI_MakeEdge(SSPL1,Surface);

          B1.Add(W10,Edges[(id_arete_contour[i])[k]]);
          B1.UpdateEdge(Edges[(id_arete_contour[i])[k]],SSPL1,Faces[f],precision);

        }
      }

      B1.Add(Faces[f],W10);

      std::cout<<" face *** "<< f <<" reconstruite "<<std::endl;

    }


///********************************************************************************************************************************
    /// FIN Reconstruction des Surfaces en utilisant le contour exterieur
///********************************************************************************************************************************

    else
    {
//if(nbAF>2)
//if(Contour_Aretes_pts.size()==1)
      {
        indicefaceReconstruite[f]=1;

        GeomPlate_BuildPlateSurface BPSurf(3,10,3);

        for (int i=0; i<nbAF; i++)
        {
          ///**************** Reconstruction des aretes NURBS BOOK
          /*int deg =2;
          std::vector<double > Knots;
          std::vector<gp_Pnt > CPnts;

          CurveSpline CS;

          CS.GlobalCurveInterp(pointsAreteFace[i], deg, Knots, CPnts, 1);

          std::vector<double > KnotsOCC;
          std::vector<int > Multiples;

          Spline S;
          S.OCCKnot_multiplicity(Knots, KnotsOCC, Multiples);

          TColgp_Array1OfPnt Poles(1,CPnts.size());

          for(int i=0; i<=CPnts.size()-1; i++)
          {
              Poles.SetValue(i+1, CPnts[i]);

          }
          TColStd_Array1OfReal Knots0(1,KnotsOCC.size());
          TColStd_Array1OfInteger Multiplicities(1,Multiples.size());


          for(int i=0; i<=KnotsOCC.size()-1; i++)
          {
              Knots0.SetValue(i+1, KnotsOCC[i]);
              Multiplicities.SetValue(i+1, Multiples[i]);
          }

          Handle(Geom_BSplineCurve) SPL =  new Geom_BSplineCurve(Poles, Knots0, Multiplicities, deg);*/
          ///****************************************************************************************

          TColgp_Array1OfPnt array(1,pointsAreteFace[i].size());
          for (int j=0; j<pointsAreteFace[i].size(); j++)
          {
            array.SetValue(j+1, (pointsAreteFace[i])[j]);

          }
          Handle(Geom_BSplineCurve) SPL = GeomAPI_PointsToBSpline(array,2,3,GeomAbs_C1).Curve();


          Handle(GeomAdaptor_HCurve) SPLAdaptor =  new GeomAdaptor_HCurve(SPL);
          Handle(BRepFill_CurveConstraint) Cont= new BRepFill_CurveConstraint(SPLAdaptor,1e-6);
          BPSurf.Add(Cont);


        }


        BPSurf.Perform();

        Standard_Integer MaxSeg = 8;
        Standard_Integer MaxDegree=8;
        Standard_Real dmax,Tol;
        Handle(GeomPlate_Surface) PSurf = BPSurf.Surface();
        dmax = Max(0.0001,10*BPSurf.G0Error());
        Tol=1e-6;

        BSplineApproxSurface Mapp(PSurf,Tol,MaxSeg,MaxDegree,dmax,0);
        Handle(Geom_BSplineSurface) Surf = Mapp.Surface();


        ///**********************************************************************************************************
        /// AJOUT DES CONTOURS ///
        ///**********************************************************************************************************



        //EY10.Orientation(TopAbs_FORWARD);

        //pcurve
        //L2d = new Geom2d_Line(gp_Pnt2d(0,0),gp_Dir2d(1,0));
        //

        B1.MakeFace(Faces[f],Surf,precision);
        B1.MakeWire(W10);


        TopoDS_Face F1; // = BRepBuilderAPI_MakeFace(Surf);

        for (int i=0; i<Contour_Aretes_pts.size(); i++)
        {

          std::vector<std::vector<gp_Pnt > > Aretes_pts = Contour_Aretes_pts[i];
          BRepBuilderAPI_MakeWire W1;

          B1.MakeWire(W10);

          for (int k=0; k<Aretes_pts.size(); k++)
          {

            std::vector<gp_Pnt > pts = Aretes_pts[k];
            TColgp_Array1OfPnt2d Array01(1,pts.size());

            for (int j=0; j<pts.size(); j++)
            {
              Standard_Real U,V;
              GeomAPI_ProjectPointOnSurf Proj(pts[j],Surf);
              Proj.LowerDistanceParameters(U, V);
              Array01.SetValue(j+1,gp_Pnt2d(U,V));
            }

            Handle(Geom2d_BSplineCurve) SSPL1 = Geom2dAPI_PointsToBSpline(Array01,2,3,GeomAbs_C1,0).Curve();

            TopoDS_Edge E = BRepBuilderAPI_MakeEdge(SSPL1,Surf);
            W1.Add(E);

            B1.Add(W10,Edges[(id_arete_contour[i])[k]]);
            B1.UpdateEdge(Edges[(id_arete_contour[i])[k]],SSPL1,Faces[f],precision);

          }

          B1.Add(Faces[f],W10);

          TopoDS_Wire W11 = W1;

          if (i==0)
          {
            F1 = BRepBuilderAPI_MakeFace(Surf,W11,Standard_True);
          }
          if (i > 0)
          {
            BRepBuilderAPI_MakeFace F2(Surf);
            F2.Init(Surf,Standard_False);
            F2.Add(W11);

            if (F2.IsDone())
            {
              W11.Reverse();
            }
            F1 = BRepBuilderAPI_MakeFace(F1,W11);
          }

        }

        B.Add(Sh,F1);

        //TopoDS_Face F01 = BRepBuilderAPI_MakeFace(Surf);
        //B.Add(Sh,F01);

        //std::cout<<" face "<< f <<" reconstruite "<<std::endl;


        ///****************************************************************************************************************
        /// FIN AJOUT DES CONTOURS ///
        ///****************************************************************************************************************
      }
    }




  }

  ///paraview
  //fprintf (fp, "%d\n", nblines);


  TopoDS_Shell Sh1;
  B1.MakeShell(Sh1);
  for (int i=0; i<nb_faces; i++)
  {
    if (indicefaceReconstruite[i]==1)
      B1.Add(Sh1,Faces[i]);
  }

  // Solid
  TopoDS_Solid Sol1;
  B1.MakeSolid(Sol1);
  B1.Add(Sol1,Sh1);

  STEPControl_Writer writer1;
  writer1.Transfer(Sol1, STEPControl_ManifoldSolidBrep);

  std::string ModelName1 = meshFileName.c_str() ;
  ModelName1.erase(ModelName1.find(".",0),5);
  ModelName1.insert(ModelName1.length(),"T.step");

  writer1.Write(ModelName1.c_str());






  TopoDS_Solid Sol;
  B.MakeSolid(Sol);
  B.Add(Sol,Sh);
  STEPControl_Writer writer;
  writer.Transfer(Sol, STEPControl_ManifoldSolidBrep);

  std::string ModelName = meshFileName.c_str() ;
  ModelName.erase(ModelName.find(".",0),5);
  ModelName.insert(ModelName.length(),".step");

  writer.Write(ModelName.c_str());



}


void BrepXfem::setMesh00(const std::string &meshFileName)
{

  _pModel = new GModel();
  _pModel->readMSH(meshFileName.c_str());
  _dim = _pModel->getNumRegions() ? 3 : 2;

  _pModelUncut = new GModel();
  _pModelUncut=_pModel;

  for (int ils=0; ils<gLevelsets.size(); ils++)
  {
    _pModel=_pModel->buildCutGModel(gLevelsets[ils]);
  }


  std::string ModelName = meshFileName.c_str() ;
  ModelName.erase(ModelName.find(".",0),4);
  ModelName.insert(ModelName.length(),"_cutLS.msh");
  _pModel->writeMSH(ModelName,2.1,false,false);
  //GModel::setCurrent(_pModel);


}





void BrepXfem::readInputFile(const std::string &fn)
{
  FILE *f = fopen(fn.c_str(), "r");
  char what[256];
  while (!feof(f))
  {
    if (fscanf(f, "%s", what) != 1) return;

    else if (!strcmp(what, "levelSetCyl"))
    {

      double x, y, z, dx, dy, dz;
      double r;
      int tag;
      if (fscanf(f, "%d %lf %lf %lf %lf %lf %lf %lf", &tag, &r, &x, &y, &z, &dx, &dy, &dz) != 8) return;
      double pt[3]={x,y,z};
      double dir[3]={dx,dy,dz};

      _ls_cut = new gLevelsetGenCylinder(pt,dir,r,tag);
      gLevelsets.push_back(_ls_cut);

    }

    else if (!strcmp(what, "levelSetPlan"))
    {

      double x, y, z, nx, ny, nz;
      int tag;
      if (fscanf(f, "%d %lf %lf %lf %lf %lf %lf", &tag, &x, &y, &z, &nx, &ny, &nz) != 7) return;
      double pt[3]={x,y,z};
      double dir[3]={nx,ny,nz};

      _ls = new gLevelsetPlane(pt,dir,tag);
      gLevelsets.push_back(_ls);

    }

    else if (!strcmp(what, "RecFile"))
    {
      char name[245];
      if (fscanf(f, "%s", name) != 1) return;
      setMesh(name);
    }

    else if (!strcmp(what, "RecFile0"))
    {
      char name[245];
      if (fscanf(f, "%s", name) != 1) return;
      setMesh0(name);
    }

    else if (!strcmp(what, "CutMesh0"))
    {
      char name[245];
      if (fscanf(f, "%s", name) != 1) return;
      setMesh00(name);
    }

    else
    {
      Msg::Error("Invalid input : %s", what);
//      return;
    }
  }
  fclose(f);
}




int test_intersect_2ls_tetra(MElement* tetra, gLevelset *_ls1, gLevelset *_ls2)
{
  std::vector<gp_Pnt> Pnts_Plan;

  for (int i=0; i<tetra->getNumFaces(); ++i)
  {
    MFace Face = tetra->getFace(i);


    std::vector<MVertex *> vmoin;
    std::vector<MVertex *> vplus;

    for (int j=0; j<Face.getNumVertices(); j++)
    {
      if ((*_ls1)(Face.getVertex(j)->x(),Face.getVertex(j)->y(),Face.getVertex(j)->z()) <-1e-3)
      {
        vmoin.push_back(Face.getVertex(j));

      }
      if ((*_ls1)(Face.getVertex(j)->x(),Face.getVertex(j)->y(),Face.getVertex(j)->z()) >1e-3)
      {
        vplus.push_back(Face.getVertex(j));

      }
    }

    if (vmoin.size() >=1 && vplus.size() >=1)
    {
      for (int ii=0; ii<vmoin.size(); ii++)
      {
        for (int jj=0; jj<vplus.size(); jj++)
        {

          MVertex* Vertex0 = vmoin[ii];
          MVertex* Vertex1 = vplus[jj];
          double valeur_level0 = (*_ls1)(vmoin[ii]->x(),vmoin[ii]->y(),vmoin[ii]->z());
          double valeur_level1 = (*_ls1)(vplus[jj]->x(),vplus[jj]->y(),vplus[jj]->z());

          double xx = ((std::abs(valeur_level0)) * Vertex1->x() + (std::abs(valeur_level1)) * Vertex0->x()) / ((std::abs(valeur_level0)) + (std::abs(valeur_level1)));
          double yy = ((std::abs(valeur_level0)) * Vertex1->y() + (std::abs(valeur_level1)) * Vertex0->y()) / ((std::abs(valeur_level0)) + (std::abs(valeur_level1)));
          double zz = ((std::abs(valeur_level0)) * Vertex1->z() + (std::abs(valeur_level1)) * Vertex0->z()) / ((std::abs(valeur_level0)) + (std::abs(valeur_level1)));

          int tt=0;
          for (int kk=0; kk<Pnts_Plan.size(); kk++)
          {
            if (gp_Vec(Pnts_Plan[kk],gp_Pnt(xx,yy,zz)).Magnitude() <=1e-9) tt=1;
          }

          if (tt ==0) Pnts_Plan.push_back(gp_Pnt(xx,yy,zz));

        }
      }

    }

  }

  //TColgp_Array1OfPnt Pnts(1,Pnts_Plan.size());

  //std::cout<<" --------------------Pnts_Plan : "<<  Pnts_Plan.size() <<"  "<<std::endl;


  gp_Vec V0P(Pnts_Plan[0], Pnts_Plan[1]);
  V0P^=gp_Vec(Pnts_Plan[0], Pnts_Plan[2]);

  gp_Pln Pln(Pnts_Plan[0],V0P);

  //GProp_PEquation GP(Pnts, 1e-7);
  //gp_Pln Pln = GP.Plane();

  std::vector<gp_Pnt> Pnts_Proj;

  for (int i=0; i<tetra->getNumFaces(); ++i)
  {
    MFace Face = tetra->getFace(i);

    std::vector<MVertex *> vmoin;
    std::vector<MVertex *> vplus;

    for (int j=0; j<Face.getNumVertices(); j++)
    {
      if ((*_ls2)(Face.getVertex(j)->x(),Face.getVertex(j)->y(),Face.getVertex(j)->z()) <-1e-3)
      {
        vmoin.push_back(Face.getVertex(j));

      }
      if ((*_ls2)(Face.getVertex(j)->x(),Face.getVertex(j)->y(),Face.getVertex(j)->z()) >1e-3)
      {
        vplus.push_back(Face.getVertex(j));

      }
    }

    if (vmoin.size() >=1 && vplus.size() >=1)
    {

      TColgp_Array1OfPnt thePnts(1,2);
      int cc = 0;
      for (int ii=0; ii<vmoin.size(); ii++)
      {
        for (int jj=0; jj<vplus.size(); jj++)
        {

          MVertex* Vertex0 = vmoin[ii];
          MVertex* Vertex1 = vplus[jj];
          double valeur_level0 = (*_ls2)(vmoin[ii]->x(),vmoin[ii]->y(),vmoin[ii]->z());
          double valeur_level1 = (*_ls2)(vplus[jj]->x(),vplus[jj]->y(),vplus[jj]->z());

          double xx = ((std::abs(valeur_level0)) * Vertex1->x() + (std::abs(valeur_level1)) * Vertex0->x()) / ((std::abs(valeur_level0)) + (std::abs(valeur_level1)));
          double yy = ((std::abs(valeur_level0)) * Vertex1->y() + (std::abs(valeur_level1)) * Vertex0->y()) / ((std::abs(valeur_level0)) + (std::abs(valeur_level1)));
          double zz = ((std::abs(valeur_level0)) * Vertex1->z() + (std::abs(valeur_level1)) * Vertex0->z()) / ((std::abs(valeur_level0)) + (std::abs(valeur_level1)));

          int tt=0;
          for (int kk=0; kk<Pnts_Proj.size(); kk++)
          {
            if (gp_Vec(Pnts_Proj[kk],gp_Pnt(xx,yy,zz)).Magnitude() <=1e-9) tt=1;
          }

          if (tt == 0) Pnts_Proj.push_back(gp_Pnt(xx,yy,zz));

        }
      }

    }

  }

  //std::cout<<" --------------------Pnts_Proj : "<<  Pnts_Proj.size() <<"  "<<std::endl;

  /// A améliorer *************

  Handle(Geom_Surface) Surf = new Geom_Plane(Pln);

  GeomAPI_ProjectPointOnSurf Proj0(Pnts_Proj[0], Surf);
  gp_Vec Vec0(Pnts_Proj[0],Proj0.NearestPoint());

  Vec0.Normalize();

  //std::cout<<" --------------------Vec0.Magnitude() : "<<  Vec0.Magnitude() <<"  "<<std::endl;

  int test = 0;

  for (int i=1; i<Pnts_Proj.size(); i++)
  {
    GeomAPI_ProjectPointOnSurf Proj(Pnts_Proj[i], Surf);
    gp_Vec Vec(Pnts_Proj[i],Proj.NearestPoint());

    //std::cout<<" --------------------Vec.Magnitude() : "<<  Vec.Magnitude() <<"  "<<std::endl;

    //std::cout<<" --------------------Vec0*Vec : "<<  Vec0*Vec <<"  "<<std::endl;

    //*************   A Vérifier  et à Améliorer **************//
    if (Vec0*Vec<-0.05) test = 1;
    //*************   A Vérifier  et à Améliorer **************//

  }



  return test;

}



Handle(Geom_BSplineSurface) face_intersect_ls_tetra(MElement* tetra, gLevelset *_ls, int nbIntersct)
{
  std::vector<Handle(GeomFill_SimpleBound) > curveBound;

  for (int i=0; i<tetra->getNumFaces(); ++i)
  {
    MFace Face = tetra->getFace(i);

    std::vector<MVertex *> vmoin;
    std::vector<MVertex *> vplus;

    for (int j=0; j<Face.getNumVertices(); j++)
    {
      if ((*_ls)(Face.getVertex(j)->x(),Face.getVertex(j)->y(),Face.getVertex(j)->z()) <-0.001)
      {
        vmoin.push_back(Face.getVertex(j));

      }
      if ((*_ls)(Face.getVertex(j)->x(),Face.getVertex(j)->y(),Face.getVertex(j)->z()) >0.001)
      {
        vplus.push_back(Face.getVertex(j));

      }
    }

    if (vmoin.size() >=1 && vplus.size() >=1)
    {

      TColgp_Array1OfPnt thePnts(1,2);
      int cc = 0;
      for (int ii=0; ii<vmoin.size(); ii++)
      {
        for (int jj=0; jj<vplus.size(); jj++)
        {

          MVertex* Vertex0 = vmoin[ii];
          MVertex* Vertex1 = vplus[jj];
          double valeur_level0 = (*_ls)(vmoin[ii]->x(),vmoin[ii]->y(),vmoin[ii]->z());
          double valeur_level1 = (*_ls)(vplus[jj]->x(),vplus[jj]->y(),vplus[jj]->z());

          double xx = ((std::abs(valeur_level0)) * Vertex1->x() + (std::abs(valeur_level1)) * Vertex0->x()) / ((std::abs(valeur_level0)) + (std::abs(valeur_level1)));
          double yy = ((std::abs(valeur_level0)) * Vertex1->y() + (std::abs(valeur_level1)) * Vertex0->y()) / ((std::abs(valeur_level0)) + (std::abs(valeur_level1)));
          double zz = ((std::abs(valeur_level0)) * Vertex1->z() + (std::abs(valeur_level1)) * Vertex0->z()) / ((std::abs(valeur_level0)) + (std::abs(valeur_level1)));

          thePnts(cc+1) = gp_Pnt(xx,yy,zz);
          cc++;

        }
      }


      Handle(Geom_BSplineCurve) SPL11 = GeomAPI_PointsToBSpline(thePnts).Curve();
      Handle(GeomAdaptor_HCurve) SPL11Adaptor = new GeomAdaptor_HCurve(SPL11);
      Handle(GeomFill_SimpleBound) B11 = new GeomFill_SimpleBound(SPL11Adaptor,Precision::Approximation(),Precision::Angular());

      curveBound.push_back(B11);

    }

  }

  GeomFill_ConstrainedFilling aConstrainedFilling1(8, 2);

  if (nbIntersct == 3)
  {
    aConstrainedFilling1.Init(curveBound[0],curveBound[1],curveBound[2],Standard_False);
  }

  if (nbIntersct == 4)
  {
    aConstrainedFilling1.Init(curveBound[0],curveBound[1],curveBound[2],curveBound[3],Standard_False);
  }

  Handle(Geom_BSplineSurface) aBSplineSurface1 = aConstrainedFilling1.Surface();

  return aBSplineSurface1;

}


Handle(Geom_Curve) Organiser_Segments(std::vector<Handle(Geom_Curve) > Segments_Intersection)
{

  for (int i=0; i<Segments_Intersection.size(); i++)
  {
    Handle(Geom_Curve) curve1 = Segments_Intersection[i];

    double t10 = curve1->FirstParameter();
    double t11 = curve1->LastParameter();
    gp_Pnt P10 = curve1->Value(t10);
    gp_Pnt P11 = curve1->Value(t11);

    std::cout<<" iii  :      "<<   i  <<"  "<<std::endl;

    std::cout<<" ------------ Segments_Intersection[i]  :      "<<   P10.X()  <<"  "<<   P10.Y()  <<"  "<<   P10.Z()  <<"  "<<std::endl;
    std::cout<<" ------------ Segments_Intersection[i]  :      "<<   P11.X()  <<"  "<<   P11.Y()  <<"  "<<   P11.Z()  <<"  "<<std::endl;

  }

  std::vector<Handle(Geom_Curve) > temp;

  // Chercher le premier Segment
  Handle(Geom_Curve) curve0;
  int test0, test1;
  int cc=0;
  do
  {
    test0 = 0; ///test0 et test1 sont utilisés pour voir si le segment est lié à un autre par son premier sommet (test0) ou son deuxieme sommet (test1)
    test1 = 0;
    curve0 = Segments_Intersection[cc];

    double t00 = curve0->FirstParameter();
    double t01 = curve0->LastParameter();
    gp_Pnt P00 = curve0->Value(t00);
    gp_Pnt P01 = curve0->Value(t01);

    for (int j=0; j<Segments_Intersection.size(); j++)
    {
      if (j!=cc)
      {
        Handle(Geom_Curve) curve1 = Segments_Intersection[j];

        double t10 = curve1->FirstParameter();
        double t11 = curve1->LastParameter();
        gp_Pnt P10 = curve1->Value(t10);
        gp_Pnt P11 = curve1->Value(t11);

        if ((gp_Vec(P00, P10).Magnitude() <1e-1 || gp_Vec(P00, P11).Magnitude() <1e-1)) test0 = 1;          /// Vérifier et Améliorer
        if ((gp_Vec(P01, P10).Magnitude() <1e-1 || gp_Vec(P01, P11).Magnitude() <1e-1)) test1 = 1;          /// Vérifier et Améliorer
      }

    }

    cc++;
  }
  while (cc<Segments_Intersection.size() && (test0==1 && test1==1));

  std::cout<<" ------------ ccc  :      "<<   cc  <<"  test0    "<<   test0  <<"  test1    "<<   test1  <<"  "<<std::endl;

  temp.push_back(curve0);

  double t00 = curve0->FirstParameter();
  double t01 = curve0->LastParameter();
  gp_Pnt P00 = curve0->Value(t00);
  gp_Pnt P01 = curve0->Value(t01);

  std::cout<<" ------------ Curve0.first  :      "<<   P00.X()  <<"  "<<   P00.Y()  <<"  "<<   P00.Z()  <<"  "<<std::endl;
  std::cout<<" ------------ Curve0.last  :      "<<   P01.X()  <<"  "<<   P01.Y()  <<"  "<<   P01.Z()  <<"  "<<std::endl;


  //std::cout<<" --------------------              iiiii  :           "<<  i <<"  "<<std::endl;

  int test[Segments_Intersection.size()];

  for (int i=0; i<Segments_Intersection.size(); i++)
  {
    test[i] =0;
  }

  test[cc-1] = 1;

  //int cc0=0;
  gp_Pnt Pnt_Candidat;

  std::vector<gp_Pnt > liste_Pnt;

  do
  {
    //curve0 = temp[temp.size()-1];

    if (temp.size() == 1)
    {
      if (test0 ==1)
      {
        Pnt_Candidat = temp[0]->Value(temp[0]->FirstParameter());
        liste_Pnt.push_back(temp[0]->Value(temp[0]->LastParameter()));
        liste_Pnt.push_back(temp[0]->Value(temp[0]->FirstParameter()));
      }
      if (test1 ==1&&test0==0)
      {
        Pnt_Candidat = temp[0]->Value(temp[0]->LastParameter());
        liste_Pnt.push_back(temp[0]->Value(temp[0]->FirstParameter()));
        liste_Pnt.push_back(temp[0]->Value(temp[0]->LastParameter()));
      }

    }

    double dmin = 1e6;

    int indice;
    for (int i=0; i<Segments_Intersection.size(); i++)
    {
      Handle(Geom_Curve) curve1 = Segments_Intersection[i];

      double t10 = curve1->FirstParameter();
      double t11 = curve1->LastParameter();
      gp_Pnt P10 = curve1->Value(t10);
      gp_Pnt P11 = curve1->Value(t11);

      if (gp_Vec(Pnt_Candidat, P10).Magnitude() <dmin)
      {
        if (test[i] == 0)
        {
          //Pnt_Candidat = P11;
          indice = i;
          dmin = gp_Vec(Pnt_Candidat, P10).Magnitude();

        }
      }
      if (gp_Vec(Pnt_Candidat, P11).Magnitude() <dmin)
      {
        if (test[i] == 0)
        {
          //Pnt_Candidat = P10;
          indice = i;
          dmin = gp_Vec(Pnt_Candidat, P11).Magnitude();

        }
      }
    }

    temp.push_back(Segments_Intersection[indice]);
    test[indice]=1;

    if (gp_Vec(Pnt_Candidat, Segments_Intersection[indice]->Value(Segments_Intersection[indice]->FirstParameter())).Magnitude()- dmin <= 1e-32)
    {
      Pnt_Candidat = Segments_Intersection[indice]->Value(Segments_Intersection[indice]->LastParameter());
    }
    if (gp_Vec(Pnt_Candidat, Segments_Intersection[indice]->Value(Segments_Intersection[indice]->LastParameter())).Magnitude() - dmin <=1e-32)
    {
      Pnt_Candidat = Segments_Intersection[indice]->Value(Segments_Intersection[indice]->FirstParameter());
    }

    liste_Pnt.push_back(Pnt_Candidat);



    //temp.push_back();
    //cc0++;

  }
  while (temp.size() <Segments_Intersection.size());



  for (int i=0; i<liste_Pnt.size(); i++)
  {
    std::cout<<" ------------ liste_Pnt[i]  :      "<<   liste_Pnt[i].X()  <<"  "<<   liste_Pnt[i].Y()  <<"  "<<   liste_Pnt[i].Z()  <<"  "<<std::endl;

  }


  std::cout<<" ------------ Segments_Intersection.size()  :      "<<   Segments_Intersection.size()  <<"  "<<std::endl;
  std::cout<<" ------------ liste_Pnt.size  :      "<<   liste_Pnt.size()  <<"  "<<std::endl;


  std::vector<gp_Pnt > liste_Pnt_A;

  liste_Pnt_A.push_back(liste_Pnt[0]);

  gp_Pnt Pnt_test = liste_Pnt[0];
  for (int i=1; i<liste_Pnt.size(); i++)
  {
    if (Pnt_test.Distance(liste_Pnt[i]) >1e-3)      /// a verifier
    {
      liste_Pnt_A.push_back(liste_Pnt[i]);
      Pnt_test = liste_Pnt[i];

    }

  }

  std::cout<<" ------------ liste_Pnt_A.size  :      "<<   liste_Pnt_A.size()  <<"  "<<std::endl;

  TColgp_Array1OfPnt Array(1,liste_Pnt_A.size());

  for (int i=1; i<=liste_Pnt_A.size(); i++)
  {
    Array.SetValue(i,liste_Pnt_A[i-1]);
    std::cout<<" ------------ liste_Pnt_A[i]  :      "<<   liste_Pnt_A[i-1].X()  <<"  "<<   liste_Pnt_A[i-1].Y()  <<"  "<<   liste_Pnt_A[i-1].Z()  <<"  "<<std::endl;

  }

  Handle(Geom_BSplineCurve) Spline = GeomAPI_PointsToBSpline(Array,3,8,GeomAbs_C1,0).Curve();


  return Spline;


}



Handle(Geom_Curve) Courbe_porteuse_Edge(Handle(Geom_Curve) Curve, gp_Pnt Vertex0, gp_Pnt Vertex1)
{

  double tmin = Curve->FirstParameter();
  double tmax = Curve->LastParameter();


  GeomAPI_ProjectPointOnCurve Proj0(Vertex0, Curve);
  GeomAPI_ProjectPointOnCurve Proj1(Vertex1, Curve);

  double t0 = Proj0.LowerDistanceParameter();
  double t1 = Proj1.LowerDistanceParameter();

  std::cout<<" ------------ Proj0.Parameter(1)  :      "<<   t0  <<"  "<<std::endl;
  std::cout<<" ------------ Proj1.Parameter(1)  :      "<<   t1  <<"  "<<std::endl;

  std::cout<<" ------------ Vertex0  :      "<<   Vertex0.X()   <<"  "<<   Vertex0.Y()   <<"  "<<   Vertex0.Z()   <<"  "<<std::endl;
  std::cout<<" ------------ Vertex1  :      "<<   Vertex1.X()   <<"  "<<   Vertex1.Y()   <<"  "<<   Vertex1.Z()   <<"  "<<std::endl;

  gp_Pnt V0, V1;

  if (t0<t1)
  {
    V0 = gp_Pnt(Vertex0.X(), Vertex0.Y(), Vertex0.Z());
    V1 = gp_Pnt(Vertex1.X(), Vertex1.Y(), Vertex1.Z());
  }
  else
  {
    V0 = gp_Pnt(Vertex1.X(), Vertex1.Y(), Vertex1.Z());
    V1 = gp_Pnt(Vertex0.X(), Vertex0.Y(), Vertex0.Z());
  }

  GeomAPI_ProjectPointOnCurve ProjV0(V0, Curve);
  GeomAPI_ProjectPointOnCurve ProjV1(V1, Curve);

  int i=0;
  double pas = (tmax-tmin) /10;
  double tin, tax;
  do
  {
    tin = tmin + i*pas;
    std::cout<<" ------------  tin  :      "<<   tin  <<"  "<<std::endl;
    i++;
  }
  while (tin<=ProjV0.LowerDistanceParameter());


  i=0;
  do
  {
    tax = tmax - i*pas;
    std::cout<<" ------------  tax  :      "<<   tax  <<"  "<<std::endl;
    i++;
  }
  while (tax>=ProjV1.LowerDistanceParameter());

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

  std::vector<gp_Pnt > liste_Pnt;

  liste_Pnt.push_back(V0);
  for (int i=0; i<10; i++)
  {
    double t = tin + i* ((tax-tin) /10);
    liste_Pnt.push_back(Curve->Value(t));

  }
  liste_Pnt.push_back(V1);

  for (int i=0; i<liste_Pnt.size(); i++)
  {
    std::cout<<" ------------  liste_Pnt  :      "<<   liste_Pnt[i].X()  <<"  "<<   liste_Pnt[i].Y()  <<"  "<<   liste_Pnt[i].Z()  <<"  "<<std::endl;

  }

  TColgp_Array1OfPnt Array(1,liste_Pnt.size());

  for (int i=1; i<=liste_Pnt.size(); i++)
  {
    Array.SetValue(i,liste_Pnt[i-1]);

  }

  Handle(Geom_BSplineCurve) Spline = GeomAPI_PointsToBSpline(Array,3,8,GeomAbs_C1,0).Curve();


  return Spline;


}


std::vector<TopoDS_Edge> Organiser_Aretes(std::vector<TopoDS_Edge> Edges)
{

  int test[Edges.size()];
  for (int i=0; i<Edges.size(); i++)
  {
    test[i]=0;
  }
  std::vector<TopoDS_Edge> Edges_temp;
  Edges_temp.push_back(Edges[0]);
  test[0] = 1;
  int compteur =1;

  TopoDS_Edge EdgeC = Edges[0];
  TopoDS_Vertex V0c = TopExp::FirstVertex(EdgeC);
  TopoDS_Vertex V1c = TopExp::LastVertex(EdgeC);

  do
  {

    for (int i=0; i<Edges.size(); i++)
    {
      TopoDS_Vertex V0 = TopExp::FirstVertex(Edges[i]);
      TopoDS_Vertex V1 = TopExp::LastVertex(Edges[i]);
      if ((test[i]==0) && ((BRepTools::Compare(V0,V0c)) || (BRepTools::Compare(V0,V1c)) || (BRepTools::Compare(V1,V0c)) || (BRepTools::Compare(V1,V1c))))
      {
        Edges_temp.push_back(Edges[i]);
        test[i]=1;
        V0c=V0;
        V1c=V1;
        compteur++;

      }
    }

  }
  while (compteur!=Edges.size());



  return Edges_temp;

}