// 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 "RecXfem.h"

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


std::vector<TopoDS_Edge> Organiser_Edges(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;

}





void RecXfem::Reconstruire()
{

  BRep_Builder B;
  TopoDS_Shell Sh;

  B.MakeShell(Sh);

  //FILE * pFile;
  //pFile = fopen ("myfile1.txt","rt");

  int nb_face;

  nb_face = nb_faces;

  //fscanf(pFile, "%d\n", &nb_face);

  std::cout<<"num Faces : "<<  nb_face <<std::endl;

  Handle(Geom_Surface) Surfs[nb_face];

  int jj=0;

  for (int i=0; i<nb_face; ++i)
  {

    int nb_pts;
    //fscanf(pFile, "%d\n", &nb_pts);
    nb_pts = nb_pts_face[i];


    std::cout<<"num pts : "<<  nb_pts <<std::endl;



    GeomPlate_BuildPlateSurface BPSurf(3,10,3);

    for (int j=0; j<nb_pts; ++j)
    {

      double x, y, z;
      //fscanf(pFile, "%lf %lf %lf\n", &x, &y, &z);
      npoint3 Pt = pts[jj];
      jj=jj+1;
      x = Pt.x();
      y = Pt.y();
      z = Pt.z();
      gp_Pnt P(x,y,z);

      Handle(GeomPlate_PointConstraint) PCont= new GeomPlate_PointConstraint(P,0);
      BPSurf.Add(PCont);

      //std::cout<<"x : "<<  x <<" y : "<<  y <<" z : "<<  z <<std::endl;

    }


    BPSurf.Perform();
    Standard_Integer MaxSeg = 1;
    Standard_Integer MaxDegree=3;
    //Standard_Integer CritOrder=4;
    Standard_Real dmax,Tol;
    Handle(GeomPlate_Surface) PSurf1 = BPSurf.Surface();
    dmax = Max(0.0001,10*BPSurf.G0Error());
    Tol=0.0;
    GeomPlate_MakeApprox Mapp1(PSurf1,Tol,MaxSeg,MaxDegree,dmax,0);
    Handle(Geom_Surface) Surf1 = Mapp1.Surface();

    Surfs[i] = Surf1;

    //TopoDS_Face F = BRepBuilderAPI_MakeFace(Surf1);
    //B.Add(Sh,F);
  }

  std::vector<Handle(Geom_Curve) > courbes_limit[nb_face];

  for (int i=0; i<nb_face; i++)
  {
    Handle(Geom_Surface) Surf = Surfs[i];

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

    for (int j=0; j<nb_face; j++)
    {
      if (i!=j)
      {
        Handle(Geom_Surface) Surf1 = Surfs[j];
        GeomAPI_IntSS SSI(Surf1,Surf, 1e-6);

        //std::cout<<"NbLines()  "<<  SSI.NbLines() <<std::endl;

        for (int l=1; l<=SSI.NbLines(); l++)
        {
          courbes.push_back(SSI.Line(l));

        }

      }

    }

    courbes_limit[i] = courbes;

  }


  for (int i=0; i<nb_face; i++)
  {


    std::vector<Handle(Geom_Curve) > courbes = courbes_limit[i];

    std::cout<<"Nb_courbes_limit  "<<  courbes.size() <<std::endl;

    std::vector<Handle(Geom2d_Curve) > courbes2d;

    for (int j=0; j<courbes.size();j++)
    {
      Handle(Geom2d_Curve) Curve2d = GeomProjLib::Curve2d(courbes[j],Surfs[i]);
      courbes2d.push_back(Curve2d);

    }



    BRepBuilderAPI_MakeWire W;
    TopoDS_Edge E1, E2, E3, E4;

    std::vector<TopoDS_Edge> Edges;
    std::vector<gp_Pnt* > Pnt_lim;

    for (int j=0; j<courbes2d.size(); j++)
    {
      gp_Pnt2d Point_lim[2];
      int kp=0;
      for (int k=0; k<courbes2d.size(); k++)
      {
        if (j!=k)
        {

          Geom2dAPI_InterCurveCurve CCI(courbes2d[k],courbes2d[j]);
          std::cout<<"NbPoints()  "<<  CCI.NbPoints() <<std::endl;
          if (CCI.NbPoints() !=0)
          {

            //for(int np=1; np<=CCI.NbPoints(); np++){

            std::cout<<" kp  "<<  kp  <<std::endl;
            Point_lim[kp] = CCI.Point(1);
            kp++;
            //}

          }
        }
      }


      Geom2dAPI_ProjectPointOnCurve P0(Point_lim[0], courbes2d[j]);
      Geom2dAPI_ProjectPointOnCurve P1(Point_lim[1], courbes2d[j]);

      gp_Pnt P3d_0 = Surfs[i]->Value(Point_lim[0].X(), Point_lim[0].Y());
      gp_Pnt P3d_1 = Surfs[i]->Value(Point_lim[1].X(), Point_lim[1].Y());

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


      double t0 = P0.Parameter(1);
      double t1 = P1.Parameter(1);

      std::cout<<" t0 "<<  t0 <<" t1 "<<  t1 <<"  "<< std::endl;


      //std::cout<<" t0  "<<  t0 <<" t1  "<<  t1 <<std::endl;



      //double Ufirst = std::min(Point_lim[0].X(), Point_lim[1].X());
      //double Ulast = std::max(Point_lim[0].X(), Point_lim[1].X());

      //double Vfirst = std::min(Point_lim[0].Y(), Point_lim[1].Y());
      //double Vlast = std::max(Point_lim[0].Y(), Point_lim[1].Y());


      //std::cout<<"Point_lim[0].X()  "<<  Point_lim[0].X()<<" Point_lim[1].X()  "<<  Point_lim[1].X()<<" Point_lim[0].Y()  "<<  Point_lim[0].Y()<<" Point_lim[1].Y()  "<<  Point_lim[1].Y() <<std::endl;


      //Handle(Geom2d_Curve) Curve2d = GeomProjLib::Curve2d(courbes[j],Surfs[i],Ufirst,Ulast,Vfirst,Vlast);

      TopoDS_Edge E = BRepBuilderAPI_MakeEdge(courbes2d[j],Surfs[i],t0,t1);
      Edges.push_back(E);



      /*
                          if(j==0) E1 = BRepBuilderAPI_MakeEdge(courbes2d[j],Surfs[i],t0,t1);
                          if(j==1) E2 = BRepBuilderAPI_MakeEdge(courbes2d[j],Surfs[i],t0,t1);
                          if(j==2) E3 = BRepBuilderAPI_MakeEdge(courbes2d[j],Surfs[i],t0,t1);
                          if(j==3) E4 = BRepBuilderAPI_MakeEdge(courbes2d[j],Surfs[i],t0,t1);
      */
    }

    std::vector<TopoDS_Edge> EEdges;
    EEdges = Organiser_Edges(Edges);

    std::cout<<" EEdges Sizzzz  "<<  EEdges.size() <<" "<< std::endl;

    for (int j=0; j<EEdges.size(); j++)
    {
      W.Add(EEdges[j]);
    }


    /*W.Add(E1);
    W.Add(E4);
    W.Add(E2);
    W.Add(E3);*/

    TopoDS_Wire W1 = W;

    //W1.Reverse();
    //TopoDS_Face F = BRepBuilderAPI_MakeFace(Surfs[i]);
    TopoDS_Face F = BRepBuilderAPI_MakeFace(Surfs[i],W1);
    B.Add(Sh,F);


  }






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



}

// int &nb_faces, std::vector<int> &nb_pts_face, std::vector<npoint3> &pts

void RecXfem::get_pts_faces()
{


//FILE * pFile1;
//pFile1 = fopen ("myfile1.txt","wt");


//std::cout << "level set size : " << gLevelsets.size() << std::endl ;
//fprintf(pFile1, "%d\n", gLevelsets.size());

  nb_faces = gLevelsets.size();

  std::cout << "nb_faces : " << nb_faces << std::endl ;


  for (int ils=0; ils<nb_faces; ils++)

  {

    //_pModel=_pModel->buildCutGModel(gLevelsets[ils]);


    gLevelset *_ls = gLevelsets[ils];


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


    std::set<double* > Verts;

    std::set<MElement *> Elements;

    std::vector<MEdge > Edges;

    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 i=0; i<numElements[0]; i++)
      {

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



        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 ((std::abs(valeur_level0) <=1e-6 || std::abs(valeur_level1) <=1e-6) || (valeur_level0*valeur_level1 <=-1e-6))
          {

            Edges.push_back(Edge);

          }

        }


      }

    }

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

    std::vector<MEdge > Edges1;

    Edges1.push_back(Edges[0]);

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

      for (int j=0; j<Edges1.size(); j++)
      {

        if (Edges[i]==Edges1[j]) eq =1;

      }
      if (eq ==0) Edges1.push_back(Edges[i]);

    }


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



    std::vector<double> vect_x;
    std::vector<double> vect_y;
    std::vector<double> vect_z;

    for (int i=0; i<Edges1.size(); i++)
    {
      MEdge Edge = Edges1[i];

      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 (std::abs(valeur_level0) <=1e-6)
      {
        vect_x.push_back(Vertex0->x());
        vect_y.push_back(Vertex0->y());
        vect_z.push_back(Vertex0->z());
      }

      if (std::abs(valeur_level1) <=1e-6)
      {
        vect_x.push_back(Vertex1->x());
        vect_y.push_back(Vertex1->y());
        vect_z.push_back(Vertex1->z());
      }

      if (valeur_level0*valeur_level1 < -1e-6)
      {

        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)));

        vect_x.push_back(xx);
        vect_y.push_back(yy);
        vect_z.push_back(zz);
      }


    }








    //fprintf(pFile1, "%d\n", vect_x.size());

    nb_pts_face.push_back(vect_x.size());

    std::cout << "nb_pts_face[ils]  : " << nb_pts_face[ils] << std::endl ;


    for (int i=0; i<vect_x.size(); i++)
    {
      //fprintf(pFile1, "%lf %lf %lf\n", vect_x[i], vect_y[i], vect_z[i]);

      pts.push_back(npoint3(vect_x[i], vect_y[i], vect_z[i]));

    }






  }


//fclose(pFile1);


  ReconstruireCylmf();



}




void RecXfem::setMesh(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);

  //setOriginalPhysicals(_pModel,_pModelUncut);

  /*int nbf;
  std::vector<int> nbpts_f;
  std::vector<npoint3> pts;*/


  get_pts_faces();

}


void RecXfem::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;


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

    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
    {
      Msg::Error("Invalid input : %s", what);
//      return;
    }
  }
  fclose(f);
}


std::vector<std::vector<gp_Pnt> > Surf2Surf(std::vector<npoint3> pts)
{

  TColgp_Array1OfPnt Array_Point(1,pts.size());

  for (int i = 1; i<=pts.size();i++)
  {

    Array_Point.SetValue(i,gp_Pnt(pts[i-1].x(), pts[i-1].y(), pts[i-1].z()));
  }

  GProp_PEquation Boite = GProp_PEquation(Array_Point,0);

  gp_Pnt P_def;
  gp_Vec V11_def;
  gp_Vec V22_def;
  gp_Vec V33_def;

  Boite.Box(P_def, V11_def, V22_def, V33_def);

  //V11_def.Normalize(); V22_def.Normalize(); V33_def.Normalize();

  std::cout<<" P_def : "<<  P_def.X() <<"  "<<  P_def.Y() <<"  "<<  P_def.Z() <<std::endl;
  std::cout<<" V11_def : "<<  V11_def.X() <<"  "<<  V11_def.Y() <<"  "<<  V11_def.Z() <<std::endl;
  std::cout<<" V22_def : "<<  V22_def.X() <<"  "<<  V22_def.Y() <<"  "<<  V22_def.Z() <<std::endl;
  std::cout<<" V33_def : "<<  V33_def.X() <<"  "<<  V33_def.Y() <<"  "<<  V33_def.Z() <<std::endl;

  gp_Vec V = (V11_def+V22_def) +V33_def;

  gp_Pnt Pt(P_def.X() +V.X(),P_def.Y() +V.Y(),P_def.Z() +V.Z());

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

  std::cout<<" V11_def.Magnitude() : "<<  V11_def.Magnitude() <<"  "<<std::endl;
  std::cout<<" V22_def.Magnitude() : "<<  V22_def.Magnitude() <<"  "<<std::endl;
  std::cout<<" V33_def.Magnitude() : "<<  V33_def.Magnitude() <<"  "<<std::endl;

  std::vector<gp_Pnt> demiPnt1;
  std::vector<gp_Pnt> demiPnt2;

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

  for (int i = 1; i<=pts.size();i++)
  {
    gp_Vec V(Array_Point(i).X()-P_def.X(),Array_Point(i).Y()-P_def.Y(),Array_Point(i).Z()-P_def.Z());

    double Produit = V*V22_def.Normalized();
    std::cout<<" Produit : "<<  Produit <<"  "<<std::endl;

    if (Produit <= V22_def.Magnitude() /2) demiPnt1.push_back(Array_Point(i));
    if (Produit >= V22_def.Magnitude() /2) demiPnt2.push_back(Array_Point(i));

  }

  std::vector<std::vector<gp_Pnt> > deux_demiPnt;

  deux_demiPnt.push_back(demiPnt1);
  deux_demiPnt.push_back(demiPnt2);

  std::cout<<" demiPnt1.size() : "<<  demiPnt1.size() <<"  "<<std::endl;
  std::cout<<" demiPnt2.size() : "<<  demiPnt2.size() <<"  "<<std::endl;

  return deux_demiPnt;



  /*
   Array_Boite[1] = P_def.X();
   Array_Boite[2] = P_def.Y();
   Array_Boite[3] = P_def.Z();

   Array_Boite[4] = V11_def.X();
   Array_Boite[5] = V11_def.Y();
   Array_Boite[6] = V11_def.Z();

   Array_Boite[7] = V22_def.X();
   Array_Boite[8] = V22_def.Y();
   Array_Boite[9] = V22_def.Z();

   Array_Boite[10] = V33_def.X();
   Array_Boite[11] = V33_def.Y();
   Array_Boite[12] = V33_def.Z();
  */




}


Handle(Geom_BSplineSurface) convertCylBSpline(gp_Cylinder agpCylindre, double u1, double u2, double v1, double v2)
{

  Convert_CylinderToBSplineSurface Convert_Cyl_BSpline(agpCylindre, u1, u2, v1, v2);

  int udeg = Convert_Cyl_BSpline.UDegree();
  int vdeg = Convert_Cyl_BSpline.VDegree();

  int nbknotsU = Convert_Cyl_BSpline.NbUKnots();
  int nbknotsV = Convert_Cyl_BSpline.NbVKnots();

  int nbUPoles = Convert_Cyl_BSpline.NbUPoles();
  int nbVPoles = Convert_Cyl_BSpline.NbVPoles();



  std::cout<<" Degre U : "<< udeg <<" Degre V : "<< vdeg <<std::endl;
  std::cout<<" nbknots U : "<< nbknotsU <<" nbknots V : "<< nbknotsV <<std::endl;
  std::cout<<" nb Poles U : "<< nbUPoles <<" nb Poles V : "<< nbVPoles <<std::endl;

  TColgp_Array2OfPnt Pnts(1,nbUPoles,1,nbVPoles);
  TColStd_Array2OfReal W(1,nbUPoles,1,nbVPoles);

  //TColgp_Array2OfPnt Pnts1(1,nbUPoles+1,1,nbVPoles);
  //TColStd_Array2OfReal W1(1,nbUPoles+1,1,nbVPoles);

  for (int i=1; i<=nbUPoles; i++)
  {

    for (int j=1; j<=nbVPoles; j++)
    {

      gp_Pnt Pt_Controle = Convert_Cyl_BSpline.Pole(i,j);

      double W2 = Convert_Cyl_BSpline.Weight(i,j);

      std::cout<<" X : "<< Pt_Controle.X() <<" Y : "<< Pt_Controle.Y() <<" Z : "<< Pt_Controle.Z() <<std::endl;
      std::cout<<" Wij : "<< W2 <<std::endl;

      Pnts(i,j) = Pt_Controle;
      W(i,j) = Convert_Cyl_BSpline.Weight(i,j);

      //Pnts1(i,j) = Pt_Controle;
      //W1(i,j) = Convert_Cyl_BSpline.Weight(i,j);

    }
  }
  //Pnts1(nbUPoles,1) = Pnts1(1,1);
  //W1(nbUPoles,1) = W1(1,1);
  //Pnts1(nbUPoles+1,2) = Pnts1(1,2);
  //W1(nbUPoles+1,2) = W1(1,2);

  TColStd_Array1OfReal UKnots(1,nbknotsU);
  TColStd_Array1OfReal VKnots(1,nbknotsV);
  TColStd_Array1OfInteger Mu(1,nbknotsU);
  TColStd_Array1OfInteger Mv(1,nbknotsV);

  for (int i=1; i<=nbknotsU; i++)
  {

    UKnots(i) = Convert_Cyl_BSpline.UKnot(i);
    std::cout<<" UKnot : "<< UKnots(i) <<std::endl;

  }
  for (int i=1; i<=nbknotsV; i++)
  {

    VKnots(i) = Convert_Cyl_BSpline.VKnot(i);
    std::cout<<" UKnots(i) : "<< VKnots(i) <<std::endl;

  }

  for (int i=1; i<=nbknotsU; i++)
  {

    Mu(i) = Convert_Cyl_BSpline.UMultiplicity(i);
    std::cout<<" Mu(i) : "<< Mu(i) <<std::endl;

  }
  for (int i=1; i<=nbknotsV; i++)
  {

    Mv(i) = Convert_Cyl_BSpline.VMultiplicity(i);
    std::cout<<" Mv(i) : "<< Mv(i) <<std::endl;

  }



  Handle(Geom_BSplineSurface) Surf = new Geom_BSplineSurface(Pnts, W,UKnots,VKnots, Mu, Mv, udeg, vdeg, Convert_Cyl_BSpline.IsUPeriodic(),Convert_Cyl_BSpline.IsVPeriodic());

  return Surf;
}


void RecXfem::ReconstruireCyl()
{

  BRep_Builder B;
  TopoDS_Shell Sh;

  B.MakeShell(Sh);

  //FILE * pFile;
  //pFile = fopen ("myfile1.txt","rt");

  int nb_face;
  nb_face = nb_faces;

  //fscanf(pFile, "%d\n", &nb_face);
  std::cout<<"num Faces : "<<  nb_face <<std::endl;
  //Handle(Geom_Surface) Surfs[nb_face];

  std::vector<Handle(Geom_Surface) > Surfs;
  int ii=0;

  int jj=0;
  Handle(Geom_CylindricalSurface) SS;

  for (int i=0; i<nb_face; ++i)
  {

    int nb_pts;
    nb_pts = nb_pts_face[i];



    if (i==0)
    {

      std::cout<<"num pts : "<<  nb_pts <<std::endl;
      SS =  new Geom_CylindricalSurface(gp_Ax3(gp_Pnt(50.0,50.0,0.0),gp_Vec(gp_Pnt(50.0,50.0,0.0),gp_Pnt(50.0,50.0,100.0))),40.0);
      //GeomPlate_BuildPlateSurface BPSurf(SS,3,10,3);

      std::vector<npoint3> pts1;
      for (int j=0; j<nb_pts; ++j)
      {

        double x, y, z;
        //fscanf(pFile, "%lf %lf %lf\n", &x, &y, &z);
        npoint3 Pt = pts[jj];
        jj=jj+1;
        x = Pt.x();
        y = Pt.y();
        z = Pt.z();
        gp_Pnt P(x,y,z);
        pts1.push_back(Pt);

        //GeomAPI_ProjectPointOnSurf Proj(P,SS);
        //double U,V;
        //Proj.LowerDistanceParameters(U,V);
        //***************************umax; umin********************//

        //std::cout<<" U : "<<  U <<" V : "<<  V <<std::endl;
        //Handle(GeomPlate_PointConstraint) PCont= new GeomPlate_PointConstraint(P,0);
        //BPSurf.Add(PCont);
        //std::cout<<"x : "<<  x <<" y : "<<  y <<" z : "<<  z <<std::endl;

      }

      std::vector<gp_Pnt> demiPnt1;
      std::vector<gp_Pnt> demiPnt2;

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

      for (int ip = 0; ip<pts1.size();ip++)
      {
        gp_Pnt P(pts1[ip].x(),pts1[ip].y(),pts1[ip].z());
        GeomAPI_ProjectPointOnSurf Proj(P,SS);
        double U,V;
        Proj.LowerDistanceParameters(U,V);

        if ((U <= (PI))) demiPnt1.push_back(P);
        if ((U >= (PI))) demiPnt2.push_back(P);

      }

      std::vector<std::vector<gp_Pnt> > deux_demiPnt;

      deux_demiPnt.push_back(demiPnt1);
      deux_demiPnt.push_back(demiPnt2);

      std::cout<<" demiPnt1.size() : "<<  demiPnt1.size() <<std::endl;
      std::cout<<" demiPnt2.size() : "<<  demiPnt2.size() <<std::endl;


      //****************************************************************************
      /*gp_Cylinder agpCylindre = SS->Cylinder();
      Handle(Geom_BSplineSurface) SS1 = convertCylBSpline(agpCylindre, 0.0, 100.0);
      RecBSplineLM Rec(SS1, pts1);
      Handle(Geom_BSplineSurface) Surf2 = Rec.Surface();


      TopoDS_Face F = BRepBuilderAPI_MakeFace(Surf2);
      B.Add(Sh,F);*/
      //****************************************************************************

      //std::vector<std::vector<gp_Pnt> > ptsss = Surf2Surf(pts1);
      //std::cout<<" ptsss.size() : "<<  ptsss.size() <<std::endl;

      for (int j=0; j<deux_demiPnt.size(); j++)
      {

        Handle(Geom_BSplineSurface) SS1;
        if (j==0)
        {

          gp_Cylinder agpCylindre = SS->Cylinder();
          SS1 = convertCylBSpline(agpCylindre, 0.0, PI, 0.0, 100.0);


        }
        if (j==1)
        {

          gp_Cylinder agpCylindre = SS->Cylinder();
          SS1 = convertCylBSpline(agpCylindre, PI, 2.0*PI, 0.0, 100.0);
        }


        //Handle(Geom_BoundedSurface) SSS = Handle(Geom_BoundedSurface)::DownCast(Surface);

        //GeomLib::ExtendSurfByLength(SS1, 20*PI,  2 ,Standard_True,  Standard_False);
        //GeomLib::ExtendSurfByLength(SS1, 20*PI,  2 ,Standard_True,  Standard_True);

        //Handle (Geom_Surface) Surf1 = Handle(Geom_Surface)::DownCast(SS1);

        GeomPlate_BuildPlateSurface BPSurf(3,10,3);

        for (int k=0; k<deux_demiPnt[j].size(); k++)
        {

          gp_Pnt P = (deux_demiPnt[j])[k];
          Handle(GeomPlate_PointConstraint) PCont= new GeomPlate_PointConstraint(P,0);
          BPSurf.Add(PCont);

        }

        BPSurf.Perform();
        Standard_Integer MaxSeg = 10;
        Standard_Integer MaxDegree=2;
        //Standard_Integer CritOrder=4;
        Standard_Real dmax,Tol;
        Handle(GeomPlate_Surface) PSurf1 = BPSurf.Surface();
        dmax = Max(0.0001,10*BPSurf.G0Error());
        Tol=0.0;
        GeomPlate_MakeApprox Mapp1(PSurf1,Tol,MaxSeg,MaxDegree,dmax,0);
        Handle(Geom_BSplineSurface) Surf1 = Mapp1.Surface();

        //GeomLib::ExtendSurfByLength(Surf1, 30*PI,  2 ,Standard_True,  Standard_False);
        //GeomLib::ExtendSurfByLength(Surf1, 30*PI,  2 ,Standard_True,  Standard_True);

        Surfs.push_back(Surf1);
        ii=ii+1;

        //TopoDS_Face F = BRepBuilderAPI_MakeFace(Surf1);
        //B.Add(Sh,F);
      }


    }

    else
    {

      GeomPlate_BuildPlateSurface BPSurf(3,10,3);

      std::cout<<"nb_pts : "<<  nb_pts <<std::endl;

      for (int j=0; j<nb_pts; ++j)
      {

        double x, y, z;
        //fscanf(pFile, "%lf %lf %lf\n", &x, &y, &z);
        npoint3 Pt = pts[jj];
        jj=jj+1;
        x = Pt.x();
        y = Pt.y();
        z = Pt.z();
        gp_Pnt P(x,y,z);

        Handle(GeomPlate_PointConstraint) PCont= new GeomPlate_PointConstraint(P,0);
        BPSurf.Add(PCont);

        //std::cout<<"x : "<<  x <<" y : "<<  y <<" z : "<<  z <<std::endl;

      }


      BPSurf.Perform();
      Standard_Integer MaxSeg = 1;
      Standard_Integer MaxDegree=2;
      //Standard_Integer CritOrder=4;
      Standard_Real dmax,Tol;
      Handle(GeomPlate_Surface) PSurf1 = BPSurf.Surface();
      dmax = Max(0.0001,10*BPSurf.G0Error());
      Tol=0.0;
      GeomPlate_MakeApprox Mapp1(PSurf1,Tol,MaxSeg,MaxDegree,dmax,0);
      Handle(Geom_Surface) Surf1 = Mapp1.Surface();

      Surfs.push_back(Surf1);
      ii=ii+1;

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

    }

  }

  std::vector<Handle(Geom_Curve) > courbes_limit[Surfs.size()];

  for (int i=0; i<Surfs.size(); i++)
  {
    Handle(Geom_Surface) Surf = Surfs[i];

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

    for (int j=0; j<Surfs.size(); j++)
    {
      if (i!=j)
      {
        Handle(Geom_Surface) Surf1 = Surfs[j];
        GeomAPI_IntSS SSI(Surf1,Surf, 1e-6);

        //std::cout<<"NbLines()  "<<  SSI.NbLines() <<std::endl;

        for (int l=1; l<=SSI.NbLines(); l++)
        {
          courbes.push_back(SSI.Line(l));

        }

      }

    }

    courbes_limit[i] = courbes;

  }


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


    std::vector<Handle(Geom_Curve) > courbes = courbes_limit[i];

    std::cout<<"Nb_courbes_limit  "<<  courbes.size() <<std::endl;

    std::vector<Handle(Geom2d_Curve) > courbes2d;

    for (int j=0; j<courbes.size();j++)
    {
      Handle(Geom2d_Curve) Curve2d = GeomProjLib::Curve2d(courbes[j],Surfs[i]);

      //Handle(GeomAdaptor_HCurve) SPLAdaptor =  new GeomAdaptor_HCurve(courbes[j]);
      //Adaptor3d_CurveOnSurface (const Handle(Adaptor3d_HCurve2d)& C,  const Handle(Adaptor3d_HSurface)& S);

      TopoDS_Edge E = BRepBuilderAPI_MakeEdge(courbes[j]);

      TopoDS_Vertex V0 = TopExp::FirstVertex(E);
      TopoDS_Vertex V1 = TopExp::LastVertex(E);

      gp_Pnt Pt0 = BRep_Tool::Pnt(V0);
      gp_Pnt Pt1 = BRep_Tool::Pnt(V1);

      std::cout<<" X0 pt Edge : "<<  Pt0.X() <<" Y0 pt Edge : "<<  Pt0.Y() <<" Z0 pt Edge : "<<  Pt0.Z() <<std::endl;
      std::cout<<" X1 pt Edge : "<<  Pt1.X() <<" Y1 pt Edge : "<<  Pt1.Y() <<" Z1 pt Edge : "<<  Pt1.Z() <<std::endl;

      GeomAPI_ProjectPointOnSurf Proj0(Pt0, Surfs[i], 1e-6);
      GeomAPI_ProjectPointOnSurf Proj1(Pt1, Surfs[i], 1e-6);

      double t0 = courbes[j]->FirstParameter();
      double t1 = courbes[j]->LastParameter();

      std::cout<<" FirstParameter t0 : "<<  t0 <<" LastParameter t1 : "<<  t1 <<std::endl;

      TColgp_Array1OfPnt2d Array0(1,41);

      for (int p=0; p<=40; p++)
      {
        double tt = t0 + p*0.025;

        gp_Pnt Ptcourbe = courbes[j]->Value(tt);

        std::cout<<" x(t) : "<<  Ptcourbe.X() <<" y(t) : "<<  Ptcourbe.Y() <<" z(t) : "<<  Ptcourbe.Z() <<std::endl;

        GeomAPI_ProjectPointOnSurf Proj(Ptcourbe,Surfs[i]);
        double U,V;
        Proj.LowerDistanceParameters(U, V);
        Array0.SetValue(p+1,gp_Pnt2d(U,V));

      }

      Handle(Geom2d_BSplineCurve) SSPL = Geom2dAPI_PointsToBSpline(Array0,3,8,GeomAbs_C1,0).Curve();

      //double u0, v0, u1, v1;

      //Proj0.Parameters(1, u0, v0);
      //Proj1.Parameters(1, u1, v1);

      //std::cout<<" U0 : "<<  u0 <<" V0 : "<<  v0 <<" U1 : "<<  u1 <<" V1 : "<<  v1 <<std::endl;

      //TopLoc_Location L;
      //Bnd_Box2d Boite2d = BRep_Tool::UVBox(E, Surfs[i],L);

      //courbes2d.push_back(Curve2d);

      courbes2d.push_back(SSPL);

    }



    BRepBuilderAPI_MakeWire W;
    TopoDS_Edge E1, E2, E3, E4;

    std::vector<TopoDS_Edge> Edges;

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

    for (int j=0; j<courbes2d.size(); j++)
    {

      gp_Pnt2d Point_lim[2];
      int kp=0;
      for (int k=0; k<courbes2d.size(); k++)
      {

        if (j!=k)
        {

          Geom2dAPI_InterCurveCurve CCI(courbes2d[k],courbes2d[j]);
          std::cout<<"NbPoints()  "<<  CCI.NbPoints() <<std::endl;
          if (CCI.NbPoints() !=0)
          {

            for (int np=1; np<=CCI.NbPoints(); np++)
            {

              std::cout<<" kp  "<<  kp  <<std::endl;
              Point_lim[kp] = CCI.Point(np);
              kp++;
            }

          }
        }
      }


      Geom2dAPI_ProjectPointOnCurve P0(Point_lim[0], courbes2d[j]);
      Geom2dAPI_ProjectPointOnCurve P1(Point_lim[1], courbes2d[j]);

      gp_Pnt P3d_0 = Surfs[i]->Value(Point_lim[0].X(), Point_lim[0].Y());
      gp_Pnt P3d_1 = Surfs[i]->Value(Point_lim[1].X(), Point_lim[1].Y());

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


      double t0 = P0.Parameter(1);
      double t1 = P1.Parameter(1);

      std::cout<<" t0 "<<  t0 <<" t1 "<<  t1 <<"  "<< std::endl;


      //std::cout<<" t0  "<<  t0 <<" t1  "<<  t1 <<std::endl;



      //double Ufirst = std::min(Point_lim[0].X(), Point_lim[1].X());
      //double Ulast = std::max(Point_lim[0].X(), Point_lim[1].X());

      //double Vfirst = std::min(Point_lim[0].Y(), Point_lim[1].Y());
      //double Vlast = std::max(Point_lim[0].Y(), Point_lim[1].Y());


      //std::cout<<"Point_lim[0].X()  "<<  Point_lim[0].X()<<" Point_lim[1].X()  "<<  Point_lim[1].X()<<" Point_lim[0].Y()  "<<  Point_lim[0].Y()<<" Point_lim[1].Y()  "<<  Point_lim[1].Y() <<std::endl;


      //Handle(Geom2d_Curve) Curve2d = GeomProjLib::Curve2d(courbes[j],Surfs[i],Ufirst,Ulast,Vfirst,Vlast);

      TopoDS_Edge E = BRepBuilderAPI_MakeEdge(courbes2d[j],Surfs[i],t0,t1);
      Edges.push_back(E);



      /*
                          if(j==0) E1 = BRepBuilderAPI_MakeEdge(courbes2d[j],Surfs[i],t0,t1);
                          if(j==1) E2 = BRepBuilderAPI_MakeEdge(courbes2d[j],Surfs[i],t0,t1);
                          if(j==2) E3 = BRepBuilderAPI_MakeEdge(courbes2d[j],Surfs[i],t0,t1);
                          if(j==3) E4 = BRepBuilderAPI_MakeEdge(courbes2d[j],Surfs[i],t0,t1);
      */
    }

    std::vector<TopoDS_Edge> EEdges;
    //if(Edges.size()>2)
    //{
    EEdges = Organiser_Edges(Edges);


    std::cout<<" EEdges Sizzzz  "<<  EEdges.size() <<" "<< std::endl;

    for (int j=0; j<EEdges.size(); j++)
    {
      W.Add(EEdges[j]);
    }
    //}

    //else
    //{
    //for(int j=0; j<Edges.size(); j++)
    //{
    //W.Add(Edges[0]);
    //}

    //}





    TopoDS_Wire W1 = W;

    //W1.Reverse();
    //TopoDS_Face F = BRepBuilderAPI_MakeFace(Surfs[i]);
    TopoDS_Face F = BRepBuilderAPI_MakeFace(Surfs[i],W1);
    B.Add(Sh,F);


  }






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



}




void RecXfem::ReconstruireCylmf()
{

  BRep_Builder B;
  TopoDS_Shell Sh;

  B.MakeShell(Sh);

  //FILE * pFile;
  //pFile = fopen ("myfile1.txt","rt");

  int nb_face;
  nb_face = nb_faces;

  //fscanf(pFile, "%d\n", &nb_face);
  std::cout<<"num Faces : "<<  nb_face <<std::endl;
  //Handle(Geom_Surface) Surfs[nb_face];

  std::vector<Handle(Geom_Surface) > Surfs;
  int ii=0;

  int jj=0;
  Handle(Geom_CylindricalSurface) SS;

  for (int i=0; i<1; ++i)
  {

    int nb_pts;
    nb_pts = nb_pts_face[i];



    if (i==0)
    {

      std::cout<<"num pts : "<<  nb_pts <<std::endl;
      SS =  new Geom_CylindricalSurface(gp_Ax3(gp_Pnt(50.0,50.0,0.0),gp_Vec(gp_Pnt(50.0,50.0,0.0),gp_Pnt(50.0,50.0,100.0))),30.0);
      //GeomPlate_BuildPlateSurface BPSurf(SS,3,10,3);

      std::vector<npoint3> pts1;
      for (int j=0; j<nb_pts; ++j)
      {

        double x, y, z;
        //fscanf(pFile, "%lf %lf %lf\n", &x, &y, &z);
        npoint3 Pt = pts[jj];
        jj=jj+1;
        x = Pt.x();
        y = Pt.y();
        z = Pt.z();
        gp_Pnt P(x,y,z);
        pts1.push_back(Pt);

        //GeomAPI_ProjectPointOnSurf Proj(P,SS);
        //double U,V;
        //Proj.LowerDistanceParameters(U,V);
        //***************************umax; umin********************//

        //std::cout<<" U : "<<  U <<" V : "<<  V <<std::endl;
        //Handle(GeomPlate_PointConstraint) PCont= new GeomPlate_PointConstraint(P,0);
        //BPSurf.Add(PCont);
        //std::cout<<"x : "<<  x <<" y : "<<  y <<" z : "<<  z <<std::endl;

      }

      std::vector<gp_Pnt> demiPnt1;
      std::vector<gp_Pnt> demiPnt2;

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

      for (int ip = 0; ip<pts1.size();ip++)
      {
        gp_Pnt P(pts1[ip].x(),pts1[ip].y(),pts1[ip].z());
        GeomAPI_ProjectPointOnSurf Proj(P,SS);
        double U,V;
        Proj.LowerDistanceParameters(U,V);

        if ((U <= (PI))) demiPnt1.push_back(P);
        if ((U >= (PI))) demiPnt2.push_back(P);

      }

      std::vector<std::vector<gp_Pnt> > deux_demiPnt;

      deux_demiPnt.push_back(demiPnt1);
      deux_demiPnt.push_back(demiPnt2);

      std::cout<<" demiPnt1.size() : "<<  demiPnt1.size() <<std::endl;
      std::cout<<" demiPnt2.size() : "<<  demiPnt2.size() <<std::endl;


      //****************************************************************************
      /*gp_Cylinder agpCylindre = SS->Cylinder();
      Handle(Geom_BSplineSurface) SS1 = convertCylBSpline(agpCylindre, 0.0, 100.0);
      RecBSplineLM Rec(SS1, pts1);
      Handle(Geom_BSplineSurface) Surf2 = Rec.Surface();


      TopoDS_Face F = BRepBuilderAPI_MakeFace(Surf2);
      B.Add(Sh,F);*/
      //****************************************************************************

      //std::vector<std::vector<gp_Pnt> > ptsss = Surf2Surf(pts1);
      //std::cout<<" ptsss.size() : "<<  ptsss.size() <<std::endl;

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

      for (int j=0; j<deux_demiPnt.size(); j++)
      {

        Handle(Geom_BSplineSurface) SS1;
        if (j==0)
        {

          gp_Cylinder agpCylindre = SS->Cylinder();
          SS1 = convertCylBSpline(agpCylindre, 0.0, 2*PI, 0.0, 100.0);
        }

        if (j==1)
        {

          gp_Cylinder agpCylindre = SS->Cylinder();
          SS1 = convertCylBSpline(agpCylindre, 0.0, 2*PI, 0.0, 100.0);
        }


        //Handle(Geom_BoundedSurface) SSS = Handle(Geom_BoundedSurface)::DownCast(Surface);

        //GeomLib::ExtendSurfByLength(SS1, 20*PI,  2 ,Standard_True,  Standard_False);
        //GeomLib::ExtendSurfByLength(SS1, 20*PI,  2 ,Standard_True,  Standard_True);

        //Handle (Geom_Surface) Surf1 = Handle(Geom_Surface)::DownCast(SS1);

        GeomPlate_BuildPlateSurface BPSurf(SS1,3,10,3);

        for (int k=0; k<deux_demiPnt[j].size(); k++)
        {

          gp_Pnt P = (deux_demiPnt[j])[k];
          Handle(GeomPlate_PointConstraint) PCont= new GeomPlate_PointConstraint(P,0);
          BPSurf.Add(PCont);

        }

        BPSurf.Perform();
        Standard_Integer MaxSeg = 10;
        Standard_Integer MaxDegree=2;
        //Standard_Integer CritOrder=4;
        Standard_Real dmax,Tol;
        Handle(GeomPlate_Surface) PSurf1 = BPSurf.Surface();
        dmax = Max(0.0001,10*BPSurf.G0Error());
        Tol=0.0;
        GeomPlate_MakeApprox Mapp1(PSurf1,Tol,MaxSeg,MaxDegree,dmax,0);
        Handle(Geom_BSplineSurface) Surf1 = Mapp1.Surface();

        //GeomLib::ExtendSurfByLength(Surf1, 30*PI,  2 ,Standard_True,  Standard_False);
        //GeomLib::ExtendSurfByLength(Surf1, 30*PI,  2 ,Standard_True,  Standard_True);

        Surfs.push_back(Surf1);
        ii=ii+1;

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


    }

    else
    {

      GeomPlate_BuildPlateSurface BPSurf(3,10,3);

      std::cout<<"nb_pts : "<<  nb_pts <<std::endl;

      for (int j=0; j<nb_pts; ++j)
      {

        double x, y, z;
        //fscanf(pFile, "%lf %lf %lf\n", &x, &y, &z);
        npoint3 Pt = pts[jj];
        jj=jj+1;
        x = Pt.x();
        y = Pt.y();
        z = Pt.z();
        gp_Pnt P(x,y,z);

        Handle(GeomPlate_PointConstraint) PCont= new GeomPlate_PointConstraint(P,0);
        BPSurf.Add(PCont);

        //std::cout<<"x : "<<  x <<" y : "<<  y <<" z : "<<  z <<std::endl;

      }


      BPSurf.Perform();
      Standard_Integer MaxSeg = 1;
      Standard_Integer MaxDegree=2;
      //Standard_Integer CritOrder=4;
      Standard_Real dmax,Tol;
      Handle(GeomPlate_Surface) PSurf1 = BPSurf.Surface();
      dmax = Max(0.0001,10*BPSurf.G0Error());
      Tol=0.0;
      GeomPlate_MakeApprox Mapp1(PSurf1,Tol,MaxSeg,MaxDegree,dmax,0);
      Handle(Geom_Surface) Surf1 = Mapp1.Surface();

      Surfs.push_back(Surf1);
      ii=ii+1;

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

    }

  }


  /*
           std::vector<Handle(Geom_Curve)> courbes_limit[Surfs.size()];

           for(int i=0; i<Surfs.size(); i++)
           {
               Handle (Geom_Surface) Surf = Surfs[i];

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

               for(int j=0; j<Surfs.size(); j++)
               {
                   if(i!=j)
                   {
                       Handle (Geom_Surface) Surf1 = Surfs[j];
                       GeomAPI_IntSS SSI(Surf1,Surf, 1e-6);

                       //std::cout<<"NbLines()  "<<  SSI.NbLines() <<std::endl;

                       for(int l=1; l<=SSI.NbLines(); l++)
                       {
                           courbes.push_back(SSI.Line(l));

                           }

                       }

                   }

                   courbes_limit[i] = courbes;

               }


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


                   std::vector<Handle(Geom_Curve)> courbes = courbes_limit[i];

                   std::cout<<"Nb_courbes_limit  "<<  courbes.size() <<std::endl;

                   std::vector<Handle(Geom2d_Curve)> courbes2d;

                   for(int j=0; j<courbes.size();j++)
                   {
                       Handle(Geom2d_Curve) Curve2d = GeomProjLib::Curve2d (courbes[j],Surfs[i]);

                       //Handle(GeomAdaptor_HCurve) SPLAdaptor =  new GeomAdaptor_HCurve(courbes[j]);
                       //Adaptor3d_CurveOnSurface (const Handle(Adaptor3d_HCurve2d)& C,  const Handle(Adaptor3d_HSurface)& S);

                      TopoDS_Edge E = BRepBuilderAPI_MakeEdge(courbes[j]);

                      TopoDS_Vertex V0 = TopExp::FirstVertex(E);
                      TopoDS_Vertex V1 = TopExp::LastVertex(E);

                      gp_Pnt Pt0 = BRep_Tool::Pnt(V0);
                      gp_Pnt Pt1 = BRep_Tool::Pnt(V1);

                      std::cout<<" X0 pt Edge : "<<  Pt0.X() <<" Y0 pt Edge : "<<  Pt0.Y() <<" Z0 pt Edge : "<<  Pt0.Z() <<std::endl;
                      std::cout<<" X1 pt Edge : "<<  Pt1.X() <<" Y1 pt Edge : "<<  Pt1.Y() <<" Z1 pt Edge : "<<  Pt1.Z() <<std::endl;

                      GeomAPI_ProjectPointOnSurf Proj0(Pt0, Surfs[i], 1e-6);
                      GeomAPI_ProjectPointOnSurf Proj1(Pt1, Surfs[i], 1e-6);

                      double t0 = courbes[j]->FirstParameter();
                      double t1 = courbes[j]->LastParameter();

                      std::cout<<" FirstParameter t0 : "<<  t0 <<" LastParameter t1 : "<<  t1 <<std::endl;

                      TColgp_Array1OfPnt2d Array0(1,41);

                      for(int p=0; p<=40; p++)
                      {
                          double tt = t0 + p*0.025;

                          gp_Pnt Ptcourbe = courbes[j]->Value(tt);

                          std::cout<<" x(t) : "<<  Ptcourbe.X() <<" y(t) : "<<  Ptcourbe.Y() <<" z(t) : "<<  Ptcourbe.Z() <<std::endl;

                          GeomAPI_ProjectPointOnSurf Proj(Ptcourbe,Surfs[i]);
                          double U,V;
                          Proj.LowerDistanceParameters(U, V);
                          Array0.SetValue(p+1,gp_Pnt2d(U,V));

                          }

                      Handle(Geom2d_BSplineCurve) SSPL = Geom2dAPI_PointsToBSpline(Array0,3,8,GeomAbs_C1,0).Curve();

                      //double u0, v0, u1, v1;

                      //Proj0.Parameters(1, u0, v0);
                      //Proj1.Parameters(1, u1, v1);

                      //std::cout<<" U0 : "<<  u0 <<" V0 : "<<  v0 <<" U1 : "<<  u1 <<" V1 : "<<  v1 <<std::endl;

                      //TopLoc_Location L;
                      //Bnd_Box2d Boite2d = BRep_Tool::UVBox(E, Surfs[i],L);

                      //courbes2d.push_back(Curve2d);

                      courbes2d.push_back(SSPL);

                      }



                  BRepBuilderAPI_MakeWire W;
                  TopoDS_Edge E1, E2, E3, E4;

                  std::vector<TopoDS_Edge> Edges;

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

                  for(int j=0; j<courbes2d.size(); j++)
                  {

                      gp_Pnt2d Point_lim[2];
                      int kp=0;
                      for(int k=0; k<courbes2d.size(); k++)
                      {

                          if(j!=k){

                              Geom2dAPI_InterCurveCurve CCI(courbes2d[k],courbes2d[j]);
                              std::cout<<"NbPoints()  "<<  CCI.NbPoints() <<std::endl;
                              if(CCI.NbPoints() !=0){

                                  for(int np=1; np<=CCI.NbPoints(); np++){

                                      std::cout<<" kp  "<<  kp  <<std::endl;
                                      Point_lim[kp] = CCI.Point(np);
                                      kp++;
                                      }

                                  }
                              }
                          }


                          Geom2dAPI_ProjectPointOnCurve P0(Point_lim[0], courbes2d[j]);
                          Geom2dAPI_ProjectPointOnCurve P1(Point_lim[1], courbes2d[j]);

                          gp_Pnt P3d_0 = Surfs[i]->Value (Point_lim[0].X(), Point_lim[0].Y());
                          gp_Pnt P3d_1 = Surfs[i]->Value (Point_lim[1].X(), Point_lim[1].Y());

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


                          double t0 = P0.Parameter(1);
                          double t1 = P1.Parameter(1);

                          std::cout<<" t0 "<<  t0 <<" t1 "<<  t1 <<"  "<< std::endl;


                          //std::cout<<" t0  "<<  t0 <<" t1  "<<  t1 <<std::endl;



                          //double Ufirst = std::min(Point_lim[0].X(), Point_lim[1].X());
                          //double Ulast = std::max(Point_lim[0].X(), Point_lim[1].X());

                          //double Vfirst = std::min(Point_lim[0].Y(), Point_lim[1].Y());
                          //double Vlast = std::max(Point_lim[0].Y(), Point_lim[1].Y());


                      //std::cout<<"Point_lim[0].X()  "<<  Point_lim[0].X()<<" Point_lim[1].X()  "<<  Point_lim[1].X()<<" Point_lim[0].Y()  "<<  Point_lim[0].Y()<<" Point_lim[1].Y()  "<<  Point_lim[1].Y() <<std::endl;


                      //Handle(Geom2d_Curve) Curve2d = GeomProjLib::Curve2d(courbes[j],Surfs[i],Ufirst,Ulast,Vfirst,Vlast);

                      TopoDS_Edge E = BRepBuilderAPI_MakeEdge(courbes2d[j],Surfs[i],t0,t1);
                      Edges.push_back(E);




                      }

                      std::vector<TopoDS_Edge> EEdges;
                      //if(Edges.size()>2)
                      //{
                          EEdges = Organiser_Edges(Edges);


                          std::cout<<" EEdges Sizzzz  "<<  EEdges.size() <<" "<< std::endl;

                          for(int j=0; j<EEdges.size(); j++)
                          {
                              W.Add(EEdges[j]);
                          }
                      //}

                      //else
                      //{
                          //for(int j=0; j<Edges.size(); j++)
                          //{
                              //W.Add(Edges[0]);
                          //}

                      //}





                      TopoDS_Wire W1 = W;

                      //W1.Reverse();
                      //TopoDS_Face F = BRepBuilderAPI_MakeFace(Surfs[i]);
                      TopoDS_Face F = BRepBuilderAPI_MakeFace(Surfs[i],W1);
                      B.Add(Sh,F);


                   }

  */




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



}