// 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 "GmshGlobal.h"
#include "GModel.h"
#include "BrepXfem.h"
#include <iterator>
#include "RecIncludes.h"
#include "RecXfem.h"
#include "RecBSplineLM.h"
#include "RecPrimitive.h"
#include "BSplineApproxSurface.h"

#include "outil_spline.h"




int main(int argc, char* argv[])
{


  if (argc != 2)
  {
    printf("usage : Reconstruction_xfem input_file_name\n");
    return -1;
  }

  GmshInitialize(argc, argv);
  // globals are still present in Gmsh
  // instanciate a solver
  BrepXfem recons(1000);
  // read some input file
  recons.readInputFile(argv[1]);


  /*
  GmshInitialize(argc, argv);

   GModel *m = new GModel;

   m->readGEO("suspension.geo");
   m->readMSH("suspension.msh");

   std::cout<<"num Faces : "<<  m->getNumFaces() <<std::endl;
   std::cout<<"num Aretes : "<<  m->getNumEdges() <<std::endl;
   std::cout<<"num Sommets : "<<  m->getNumVertices() <<std::endl;


   BRep_Builder B;
   TopoDS_Shell Sh;

   B.MakeShell(Sh);

   int cont = 0;


   std::vector<gp_Pnt > Pnts;
   std::map<int, gp_Pnt> mapPnts;

   int nbface=0;///****
   std::vector<int > nb_areteparface; ///****
   std::vector<int > nb_ptspararete; ///****



   for( GModel::eiter ite = m->firstEdge(); ite != m->lastEdge(); ++ite){

          std::cout<<"Tag    : "<<  (*ite)->tag() <<std::endl;

          int nbe = (*ite)->getNumMeshVertices();
          std::cout<<"nbe : "<<  nbe <<std::endl;


          GVertex *V1 = (*ite)->getBeginVertex();
          int nbv1 = V1->getNumMeshVertices();
          std::cout<<"nbv1 : "<<  nbv1 <<std::endl;
          MVertex *Nodv1 = V1->mesh_vertices[0];
          std::cout<<" Num : "<<  Nodv1->getNum() <<std::endl;
          std::cout<<" xNv1 : "<<  Nodv1->x()<<" yNv1 : "<<  Nodv1->y()<<" zNv1 : "<<  Nodv1->z() <<std::endl;

          //



          Pnts.push_back(gp_Pnt(Nodv1->x(),Nodv1->y(),Nodv1->z())); ///*****************
          mapPnts.insert(std::pair<int, gp_Pnt>(Nodv1->getNum(),gp_Pnt(Nodv1->x(),Nodv1->y(),Nodv1->z())));

          for (int ii = 0; ii < nbe; ii++){

          MVertex *Node = (*ite)->mesh_vertices[ii];
          std::cout<<" Num : "<<  Node->getNum() <<std::endl;
          std::cout<<" xNe : "<<  Node->x()<<" yNe : "<<  Node->y()<<" zNe : "<<  Node->z() <<std::endl;



          mapPnts.insert(std::pair<int, gp_Pnt>(Node->getNum(),gp_Pnt(Node->x(),Node->y(),Node->z())));
          }

          GVertex *V2 = (*ite)->getEndVertex();
          int nbv2 = V2->getNumMeshVertices();
          std::cout<<"nbv2 : "<<  nbv2 <<std::endl;
          MVertex *Nodv2 = V2->mesh_vertices[0];
          std::cout<<" Num : "<<  Nodv2->getNum() <<std::endl;
          std::cout<<" xNv2 : "<<  Nodv2->x()<<" yNv2 : "<<  Nodv2->y()<<" zNv2 : "<<  Nodv2->z() <<std::endl;

          mapPnts.insert(std::pair<int, gp_Pnt>(Nodv2->getNum(),gp_Pnt(Nodv2->x(),Nodv2->y(),Nodv2->z())));

      }

      //getNumEdges()

        FILE * pFile0;
        pFile0 = fopen("suspension.surf","wt");

        fprintf(pFile0, "%d \n",mapPnts.size());

      std::map<int,gp_Pnt>::iterator it0;
      for ( it0=mapPnts.begin() ; it0 != mapPnts.end(); it0++ )
      {
          fprintf(pFile0, "%lf %lf %lf\n", (*it0).second.X(), (*it0).second.Y(), (*it0).second.Z());
      }

      fprintf(pFile0, "%d \n",m->getNumEdges());

      //std::map<int,gp_Pnt>::iterator itmap;

      for( GModel::eiter ite = m->firstEdge(); ite != m->lastEdge(); ++ite){


          int nbe = (*ite)->getNumMeshVertices();
          fprintf(pFile0, "%d %d\n",(*ite)->tag()-1,nbe+2);


          GVertex *V1 = (*ite)->getBeginVertex();
          MVertex *Nodv1 = V1->mesh_vertices[0];

          fprintf(pFile0, "%d ",std::distance(mapPnts.begin(), mapPnts.find(Nodv1->getNum())));

          for (int ii = 0; ii < nbe; ii++){

          MVertex *Node = (*ite)->mesh_vertices[ii];

          fprintf(pFile0, "%d ",std::distance(mapPnts.begin(), mapPnts.find(Node->getNum())));
          }

          GVertex *V2 = (*ite)->getEndVertex();

          MVertex *Nodv2 = V2->mesh_vertices[0];

          fprintf(pFile0, "%d \n",std::distance(mapPnts.begin(), mapPnts.find(Nodv2->getNum())));
      }

      int ccc=0;

      fprintf(pFile0, "%d \n",m->getNumFaces());

      for (GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
      {
          std::list<GEdge*> edgs = (*it)->edges();

          if(edgs.size()>0)
          {
              ccc++;

              fprintf(pFile0, "%d %d \n",(*it)->tag(), 0);

              fprintf(pFile0, "%d \n",edgs.size());
              for( std::list<GEdge*>::const_iterator ite = edgs.begin(); ite != edgs.end(); ite++)
              {
                  fprintf(pFile0, "%d ",(*ite)->tag()-1);
              }
              fprintf(pFile0, "\n");
          }

      }

      fclose(pFile0);

      std::cout<<" number Faces : "<<  m->getNumFaces() <<std::endl;
      std::cout<<" ccc :          "<<  ccc <<std::endl;


  /*
    std::map<int, gp_Pnt> myset;
    std::map<int, gp_Pnt>::iterator it;

  myset.insert(std::pair<int, gp_Pnt>(60,gp_Pnt(0.0,0.0,0.0)));
    for (int i=1; i<=5; i++) myset.insert(std::pair<int, gp_Pnt>(i*10,gp_Pnt(0.0,0.0,0.0)));    // set: 10 20 30 40 50

  myset.insert(std::pair<int,gp_Pnt>(20,gp_Pnt(0.0,0.0,0.0)));
  //myset.insert(60);
  myset.insert(std::pair<int,gp_Pnt>(30,gp_Pnt(0.0,0.0,0.0)));


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






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

    control_points.push_back(gp_Pnt(-5.0,-5.0,-5.0));
    control_points.push_back(gp_Pnt(-5.0,-5.0,5.0));
    control_points.push_back(gp_Pnt(5.0,-5.0,5.0));
    control_points.push_back(gp_Pnt(5.0,-5.0,-5.0));
    control_points.push_back(gp_Pnt(0.0,5.0,0.0));

    control_points.push_back(gp_Pnt(-5.0,-5.0,0.0));
    control_points.push_back(gp_Pnt(5.0,-5.0,0.0));
    control_points.push_back(gp_Pnt(0.0,-5.0,-5.0));
    control_points.push_back(gp_Pnt(0.0,-5.0,5.0));

    control_points.push_back(gp_Pnt(2.0,-20.0,2.0));

    calc_tps();


    Handle(Geom_BSplineSurface) Surf = GeomAPI_PointsToBSplineSurface(Points_BS, 2, 8, GeomAbs_C1,1.0e-6).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");


      for(int i=0; i<control_points.size(); i++)
      {
          GeomAPI_ProjectPointOnSurf Proj(control_points[i],Surf);
          double distance = Proj.LowerDistance();
          std::cout<<" erreur : "<<  distance <<std::endl;


      }








  /*
      std::vector<npoint3> points;

      double y =0.0;
   for (int j0 = 0; j0 < 10; j0++)
      {

    gp_Vec Vecteur(1,0,0);

    for (int i0 = 0; i0 < 10; i0++)
    {
     gp_Vec Vec;

     if(i0%2==0) Vec = gp_Vec(100*Vecteur.X()+y, 100*Vecteur.Y(), 100*Vecteur.Z());
     if(i0%2==1) Vec = gp_Vec(100*Vecteur.X()+y, 100*Vecteur.Y(), 100*Vecteur.Z());

     points.push_back(npoint3(Vec.X(),Vec.Y(),Vec.Z()));

     Vecteur.Rotate(gp_Ax1(gp_Pnt(0,0,0),gp_Dir(0,1,0)),PI/10);

    }
   y=y+10;
   }


   RecPrimitive Recp(points);
   std::vector<double> Params = Recp.Plan();

   std::cout<<"centre[0] : "<<  Params[0] <<std::endl;
      std::cout<<"centre[1] : "<<  Params[1] <<std::endl;
      std::cout<<"centre[2] : "<<  Params[2] <<std::endl;

      std::cout<<"axe[0] : "<<  Params[3] <<std::endl;
      std::cout<<"axe[1] : "<<  Params[4] <<std::endl;
      std::cout<<"axe[2] : "<<  Params[5] <<std::endl;
  */












  /*



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

      TopoDS_Solid shape = BRepPrimAPI_MakeCylinder(Axe2,40.0,100.0);

      TopExp_Explorer exp0;

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

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

          //Faces.push_back(Face);

          std::vector<double> ParamF;

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


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

              std::cout<<" Plane : "<<std::endl;
              Handle(Geom_Plane) aPlane = Handle(Geom_Plane)::DownCast(Surface);
              gp_Pln agpPlane = aPlane->Pln();
              gp_Ax1 normale = agpPlane.Axis();

              ParamF.push_back(normale.Location().X());
              ParamF.push_back(normale.Location().Y());
              ParamF.push_back(normale.Location().Z());
              ParamF.push_back(normale.Direction().X());
              ParamF.push_back(normale.Direction().Y());
              ParamF.push_back(normale.Direction().Z());

          }


          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();
              gp_Ax1 axe = agpCylindre.Axis();
              double R = agpCylindre.Radius();
              ParamF.push_back(axe.Location().X());
              ParamF.push_back(axe.Location().Y());
              ParamF.push_back(axe.Location().Z());
              ParamF.push_back(axe.Direction().X());
              ParamF.push_back(axe.Direction().Y());
              ParamF.push_back(axe.Direction().Z());
              ParamF.push_back(R);

              Convert_CylinderToBSplineSurface Convert_Cyl_BSpline(agpCylindre, 0.0,100.0);

              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,1) = Pnts1(1,1);
              W1(nbUPoles+1,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;

              }

              //Mu(4) = 3;

              BRep_Builder B;
              TopoDS_Shell Sh;

              B.MakeShell(Sh);

              Handle(Geom_BSplineSurface) Surf = new Geom_BSplineSurface(Pnts, W,UKnots,VKnots, Mu, Mv, udeg, vdeg, Standard_True,Standard_False);

              //Surf->InsertVKnot (25.0, 2, 0.0);
              //Surf->InsertVKnot (50.0, 2, 0.0);
              //Surf->InsertVKnot (75.0, 2, 0.0);

              //Surf->InsertUKnot ((2*PI)/6, 2, PI/6);
              //Surf->InsertUKnot (PI, 2, PI/6);
              //Surf->InsertUKnot ((10*PI)/6, 2, PI/6);

              Surf->IncreaseDegree (2,2);

              //npoint3 PointsContraint[3];
              std::vector<npoint3> points;

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

              int nb;
              double x, y, z;

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

              for(int i=0; i<nb; i++)
              {
                  fscanf(pFile, "%lf %lf %lf\n", &x, &y, &z);
                  npoint3 p0(x,y,z);
                  points.push_back(p0);

                  //std::cout<<"p0 : "<<  p0.x() <<"  "<<  p0.y() <<"  "<<  p0.z() <<std::endl;

              }
              fclose(pFile);

  /*            npoint3 p0(80.0,0.0,10.0);
              points.push_back(p0);
              p0 = npoint3(0.0,80.0,10.0);
              points.push_back(p0);
              p0 = npoint3(-80.0,0.0,10.0);
              points.push_back(p0);
              p0 = npoint3(0.0,-80.0,10.0);
              points.push_back(p0);*/


  /*
              RecBSplineLM RecS(Surf, points);
              Handle(Geom_BSplineSurface) Surf1 = RecS.Surface();

              double emax = 1e-6;
              double somme=0.0;

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

                  GeomAPI_ProjectPointOnSurf Proj(gp_Pnt(points[i].x(),points[i].y(),points[i].z()),Surf1);
                  double distance = Proj.LowerDistance();

                  somme = somme + distance;
                  if(distance>emax) emax=distance;
              }

              std::cout<<"emax : "<<  emax <<std::endl;

              std::cout<<"emoy : "<<  somme/nb <<std::endl;

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

              //TopoDS_Face F1 = BRepBuilderAPI_MakeFace(agpCylindre, 0.0, 2*PI, 0.0,100.0);
              //B.Add(Sh,F1);






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




             // TColgp_Array2OfPnt Pnts_Controle(1,nbUPoles,1,nbVPoles);
              //Convert_Cyl_BSpline.Poles(Pnts_Controle);


          }

      }
  */


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

  FILE * pFile2;
  pFile2 = fopen ("myfile2.txt","wt");
  */
  /*
  FILE * pFile3;
  pFile3 = fopen ("myfile3.txt","wt");

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

              int nb;
              double x, y, z;

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

              std::cout<<"nb : "<< nb <<std::endl;

              int ii=0;

              for(int i=0; i<nb; i++)
              {
                  fscanf(pFile, "%lf %lf %lf\n", &x, &y, &z);
                  npoint3 p0(x,y,z);
                  //points.push_back(p0);

                  std::cout<<"p0 : "<<  p0.x() <<"  "<<  p0.y() <<"  "<<  p0.z() <<std::endl;

                  //if(((z+10.0))<=1e-6 && (y-50.0)>=1e-6) fprintf (pFile1, "%lf %lf %lf\n", x, y, z);
                  //if(abs((z-110.0))<=1e-6 && (y-50.0)>=1e-6) fprintf (pFile2, "%lf %lf %lf\n", x, y, z);
                  if(x>=1e-6 && (x)<=100.0 && z>=1e-6 && (z)<=100.0 && (y)>=50.0)
                  {
                      fprintf (pFile3, "%lf %lf %lf\n", x, y, z);
                      ii = ii+1;
                  }

              }
              std::cout<<"ii : "<< ii <<std::endl;

              fclose(pFile);

  fclose(pFile3);
  /*
  fclose(pFile2);
  fclose(pFile1);
  */
  /*

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

  int nb1, nb2;
  double x, y, z;

  fscanf (pFile1, "%d\n",&nb1);
  TColgp_Array1OfPnt array1 (1,nb1);

  for(int i=1; i<=nb1; i++)
  {
      fscanf(pFile1, "%lf %lf %lf\n", &x, &y, &z);
      gp_Pnt P(x,y,z);
      array1.SetValue(i,P);

  }
  fclose(pFile1);


  FILE * pFile2;
  pFile2 = fopen ("myfile2.txt","rt");

  fscanf (pFile2, "%d\n",&nb2);
  TColgp_Array1OfPnt array2 (1,nb2);

  for(int i=1; i<=nb2; i++)
  {
      fscanf(pFile2, "%lf %lf %lf\n", &x, &y, &z);
      gp_Pnt P(x,y,z);
      array2.SetValue(i,P);

  }
  fclose(pFile2);


  Handle(Geom_BSplineCurve) SPL1 = GeomAPI_PointsToBSpline(array1).Curve();
  Handle(Geom_BSplineCurve) SPL2 = GeomAPI_PointsToBSpline(array2).Curve();

  TColgp_Array1OfPnt array3 (1,2);
  gp_Pnt P31(0.0, 0.0, 0.0);
  gp_Pnt P32(0.0, 100.0, 0.0);
  array3.SetValue(1,P31);
  array3.SetValue(2,P32);
  Handle(Geom_BSplineCurve) SPL3 = GeomAPI_PointsToBSpline(array3).Curve();

  TColgp_Array1OfPnt array4 (1,2);
  gp_Pnt P41(100.0, 0.0, 0.0);
  gp_Pnt P42(100.0, 100.0, 0.0);
  array4.SetValue(1,P41);
  array4.SetValue(2,P42);
  Handle(Geom_BSplineCurve) SPL4 = GeomAPI_PointsToBSpline(array4).Curve();

  GeomPlate_BuildPlateSurface BPSurf(3,10,3);

  Handle(GeomAdaptor_HCurve) SPL1Adaptor =  new GeomAdaptor_HCurve(SPL1);
   Handle(BRepFill_CurveConstraint) Cont1= new BRepFill_CurveConstraint(SPL1Adaptor,0);

  Handle(GeomAdaptor_HCurve) SPL2Adaptor =  new GeomAdaptor_HCurve(SPL2);
   Handle(BRepFill_CurveConstraint) Cont2= new BRepFill_CurveConstraint(SPL2Adaptor,0);

  Handle(GeomAdaptor_HCurve) SPL3Adaptor =  new GeomAdaptor_HCurve(SPL3);
   Handle(BRepFill_CurveConstraint) Cont3= new BRepFill_CurveConstraint(SPL3Adaptor,0);

  Handle(GeomAdaptor_HCurve) SPL4Adaptor =  new GeomAdaptor_HCurve(SPL4);
   Handle(BRepFill_CurveConstraint) Cont4= new BRepFill_CurveConstraint(SPL4Adaptor,0);

  BPSurf.Add(Cont1);
  BPSurf.Add(Cont3);
  BPSurf.Add(Cont2);
  BPSurf.Add(Cont4);
  /*
  BPSurf.Perform();



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

              int nb;
              double x1, y1, z1;
              fscanf (pFile, "%d\n",&nb);
              int ii=0;

              for(int i=0; i<nb; i++)
              {
                  fscanf(pFile, "%lf %lf %lf\n", &x1, &y1, &z1);
                  gp_Pnt P(x1,y1,z1);
                  Handle(GeomPlate_PointConstraint) PCont= new GeomPlate_PointConstraint(P,0);
                  BPSurf.Add(PCont);

              }
              fclose(pFile);
  */
  /*

  gp_Pnt P0(50,50,50);
  Handle(GeomPlate_PointConstraint) PCont0= new GeomPlate_PointConstraint(P0,0);
  BPSurf.Add(PCont0);

  gp_Pnt P1(10,10,-10);
  Handle(GeomPlate_PointConstraint) PCont1= new GeomPlate_PointConstraint(P1,0);
  BPSurf.Add(PCont1);

  gp_Pnt P2(10,90,-10);
  Handle(GeomPlate_PointConstraint) PCont2= new GeomPlate_PointConstraint(P2,0);
  BPSurf.Add(PCont2);

  gp_Pnt P3(90,90,-10);
  Handle(GeomPlate_PointConstraint) PCont3= new GeomPlate_PointConstraint(P3,0);
  BPSurf.Add(PCont3);

  gp_Pnt P4(90,10,-10);
  Handle(GeomPlate_PointConstraint) PCont4= new GeomPlate_PointConstraint(P4,0);
  BPSurf.Add(PCont4);

  BPSurf.Perform();

  /*
      Handle (Geom_Surface) Surf = new Geom_CylindricalSurface (gp_Ax3(gp_Pnt(50.0, 50.0, -10), gp_Dir(gp_Vec(0.0, 0.0, 120.0))), 30.0);

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

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

              int nb;
              double x, y, z;

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

              for(int i=0; i<nb; i++)
              {
                  fscanf(pFile, "%lf %lf %lf\n", &x, &y, &z);
                  gp_Pnt P(x,y,z);
                  Handle(GeomPlate_PointConstraint) PCont= new GeomPlate_PointConstraint(P,0);
                  BPSurf.Add(PCont);


              }
              fclose(pFile);
  */


  /*



           Standard_Integer MaxSeg = 4;
           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;

           BSplineApproxSurface Mapp1(PSurf1,Tol,MaxSeg,MaxDegree,dmax,0);
           //GeomPlate_MakeApprox Mapp1(PSurf1,Tol,MaxSeg,MaxDegree,dmax,0);
           Handle (Geom_BSplineSurface) Surf1 = Mapp1.Surface();


           BRep_Builder B;
           TopoDS_Shell Sh;

           B.MakeShell(Sh);

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

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



    std::cout<<"ApproxError : "<< Mapp1.ApproxError() <<std::endl;
    std::cout<<"CriterionError : "<< Mapp1.CriterionError() <<std::endl;


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

              int nb;
              double x1, y1, z1;
              fscanf (pFile, "%d\n",&nb);
              int ii=0;

              double somme=0.0;

              for(int i=0; i<nb; i++)
              {
                  fscanf(pFile, "%lf %lf %lf\n", &x1, &y1, &z1);
                  gp_Pnt P(x1,y1,z1);
                  GeomAPI_ProjectPointOnSurf PPS(P,Surf1);
                  double distance = PPS.LowerDistance();

                  std::cout<<" distance : "<< distance <<std::endl;

                  somme = somme+distance;

              }
              fclose(pFile);


              std::cout<<" moyenne : "<< somme/nb <<std::endl;





  /*

      std::vector<npoint3> points;

   gp_Vec Vecteur(1,0,0);


   for (int j0 = 0; j0 <= 20; j0++)
      {
    for (int i0 = 0; i0 <= 20; i0++)
    {
     points.push_back(npoint3(100*Vecteur.X(), 100*Vecteur.Y(), 100*Vecteur.Z()));
     Vecteur.Rotate(gp_Ax1(gp_Pnt(0,0,0),gp_Dir(0,1,0)),PI/10);
    }
    Vecteur.Rotate(gp_Ax1(gp_Pnt(0,0,0),gp_Dir(0,0,1)),(PI/10));
   }



      RecPrimitive RecPrimi(points);

      std::vector<double> Param = RecPrimi.Sphere();

  std::cout<<"centre[0] : "<<  Param[0] <<std::endl;
  std::cout<<"centre[1] : "<<  Param[1] <<std::endl;
  std::cout<<"centre[2] : "<<  Param[2] <<std::endl;
  std::cout<<"Rayon : "<<  Param[3] <<std::endl;









  /*

  GmshInitialize(argc, argv);

   GModel *m = new GModel;

   m->readMSH("PP.msh");

   //std::cout<<"num Faces : "<<  m->getNumFaces() <<std::endl;
   //std::cout<<"num Aretes : "<<  m->getNumEdges() <<std::endl;
   //std::cout<<"num Sommets : "<<  m->getNumVertices() <<std::endl;

                  //////////level Set
                  int tag = 1;
                  const double pt[3]={0.5,0.5,0.0};
                  const double dir[3]={0.0,0.0,1.0};
                  double r=0.3;

                  gLevelset *_ls_cut = new gLevelsetGenCylinder(pt,dir,r,tag);
                  m = m->buildCutGModel(_ls_cut);


                  int tag1 = 1;
                  const double pt1[3]={0.5,0.5,0.5};
                  const double dir1[3]={0.0,0.0,1.0};

                  gLevelset *_ls = new gLevelsetPlane(pt1,dir1,tag1);
                  m = m->buildCutGModel(_ls);


      m->writeMSH("PP_cut.msh",2.1,false,false);
      GModel::setCurrent(m);

      GmshFinalize();

  /*


   GmshInitialize(argc, argv);

   GModel *m = new GModel;

   //m->readGEO("PP.geo");

   //m->readOCCSTEP("P1.STEP");

   //m->writeGEO("cube.geo", true);

   m->readGEO("cube.geo");


   std::cout<<"num Faces : "<<  m->getNumFaces() <<std::endl;
   std::cout<<"num Aretes : "<<  m->getNumEdges() <<std::endl;
   std::cout<<"num Sommets : "<<  m->getNumVertices() <<std::endl;

   BRep_Builder B;
   TopoDS_Shell Sh;

   B.MakeShell(Sh);

   int cont = 0;

   for (GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){



       GeomPlate_BuildPlateSurface BPSurf(3,10,3);



          std::list<GVertex*> verts = (*it)->vertices();

          int i=0;

          for( std::list<GVertex*>::const_iterator itv = verts.begin(); itv != verts.end(); itv++){

              i=i+1;


          }

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


      std::list<GEdge*> edgs = (*it)->edges();

      for( std::list<GEdge*>::const_iterator ite = edgs.begin(); ite != edgs.end(); ite++){

          /*Range<double> bounds = (*ite)->parBounds(0);
          double umin = bounds.low();
          double umax = bounds.high();

          for(int i = 1; i < (*ite)->minimumDrawSegments(); i++){
          double u = umin + (double)i / (*ite)->minimumDrawSegments() * (umax - umin);
          GPoint p = (*ite)->point(u);

          std::cout<<"  coordx : "<<  p.x() <<"  coordy : "<< p.y() << " coordz : "<< p.z() <<std::endl;

          }

          //std::cout<<"  nbpts : "<<  (*ite)->nbPointsTransfinite() <<std::endl;

  /*
          std::list<GVertex*> verts = (*ite)->vertices();

          int i=0;

          for( std::list<GVertex*>::const_iterator itv = verts.begin(); itv != verts.end(); itv++){

              i=i+1;
              std::cout<<"nbre sommets : "<<  i <<std::endl;

          }
  */


  /*
   int nbe = (*ite)->getNumMeshVertices();
   std::cout<<"nbe : "<<  nbe <<std::endl;

   TColgp_Array1OfPnt Array(1,nbe+2);
   int cc = 1;


   GVertex *V1 = (*ite)->getBeginVertex();
   int nbv1 = V1->getNumMeshVertices();
   std::cout<<"nbv1 : "<<  nbv1 <<std::endl;
   MVertex *Nodv1 = V1->mesh_vertices[0];
   std::cout<<" xNv1 : "<<  Nodv1->x()<<" yNv1 : "<<  Nodv1->y()<<" zNv1 : "<<  Nodv1->z() <<std::endl;

   Array.SetValue(cc,gp_Pnt(Nodv1->x(),Nodv1->y(),Nodv1->z()));
   cc++;




   for (int ii = 0; ii < nbe; ii++){

   MVertex *Node = (*ite)->mesh_vertices[ii];
   std::cout<<" xNe : "<<  Node->x()<<" yNe : "<<  Node->y()<<" zNe : "<<  Node->z() <<std::endl;

   Array.SetValue(cc,gp_Pnt(Node->x(),Node->y(),Node->z()));
   cc++;

   }


   GVertex *V2 = (*ite)->getEndVertex();
   int nbv2 = V2->getNumMeshVertices();
   std::cout<<"nbv2 : "<<  nbv2 <<std::endl;
   MVertex *Nodv2 = V2->mesh_vertices[0];
   std::cout<<" xNv2 : "<<  Nodv2->x()<<" yNv2 : "<<  Nodv2->y()<<" zNv2 : "<<  Nodv2->z() <<std::endl;

   Array.SetValue(cc,gp_Pnt(Nodv2->x(),Nodv2->y(),Nodv2->z()));


   Handle(Geom_BSplineCurve) SPL = GeomAPI_PointsToBSpline (Array,3,8,GeomAbs_C1,0).Curve();
   Handle(GeomAdaptor_HCurve) SPL1Adaptor =  new GeomAdaptor_HCurve(SPL);
   Handle(BRepFill_CurveConstraint) Cont= new BRepFill_CurveConstraint(SPL1Adaptor,0);
   BPSurf.Add(Cont);




  }
  */
  /*

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

           MVertex *Nod = (*it)->mesh_vertices[i];

           //std::cout<<" xN : "<<  Nod->x()<<" yN : "<<  Nod->y()<<" zN : "<<  Nod->z() <<std::endl;

           gp_Pnt P(Nod->x(),Nod->y(),Nod->z());
           Handle(GeomPlate_PointConstraint) PCont= new GeomPlate_PointConstraint(P,0);
           BPSurf.Add(PCont);

           }

           BPSurf.Perform();

           Standard_Integer MaxSeg = 50;
           Standard_Integer MaxDegree=8;
           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();


          // Construction de Contour
          BRepBuilderAPI_MakeWire W1;

          for( std::list<GEdge*>::const_iterator ite = edgs.begin(); ite != edgs.end(); ite++){

          int nbe = (*ite)->getNumMeshVertices();

          TColgp_Array1OfPnt2d Array01(1,nbe+2);
          int cc = 1;


          GVertex *V1 = (*ite)->getBeginVertex();
          int nbv1 = V1->getNumMeshVertices();

          MVertex *Nodv1 = V1->mesh_vertices[0];


          Standard_Real U,V;
          GeomAPI_ProjectPointOnSurf Proj(gp_Pnt(Nodv1->x(),Nodv1->y(),Nodv1->z()),Surf1);
          Proj.LowerDistanceParameters(U, V);
          Array01.SetValue(cc,gp_Pnt2d(U,V));

          cc++;




          for (int ii = 0; ii < nbe; ii++){

          MVertex *Node = (*ite)->mesh_vertices[ii];
          GeomAPI_ProjectPointOnSurf Proj(gp_Pnt(Node->x(),Node->y(),Node->z()),Surf1);
          Proj.LowerDistanceParameters(U, V);
          Array01.SetValue(cc,gp_Pnt2d(U,V));

          cc++;

          }


          GVertex *V2 = (*ite)->getEndVertex();
          int nbv2 = V2->getNumMeshVertices();

          MVertex *Nodv2 = V2->mesh_vertices[0];

          GeomAPI_ProjectPointOnSurf Proj1(gp_Pnt(Nodv2->x(),Nodv2->y(),Nodv2->z()),Surf1);
          Proj1.LowerDistanceParameters(U, V);
          Array01.SetValue(cc,gp_Pnt2d(U,V));


          Handle(Geom2d_BSplineCurve) SSPL1 = Geom2dAPI_PointsToBSpline(Array01,3,8,GeomAbs_C1,0).Curve();
          TopoDS_Edge E = BRepBuilderAPI_MakeEdge(SSPL1,Surf1);
          W1.Add(E);



      }

           TopoDS_Face F1 = BRepBuilderAPI_MakeFace(Surf1,W1,Standard_True);

           B.Add(Sh,F1);
  */

  /*
       }

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

  */


// stop gmsh
//  GmshFinalize();






  /*



      GmshInitialize(argc, argv);

    //****** lire le fichier

      GModel *m = new GModel();

      m->readMSH("PP.msh");

      std::cout<<"num Faces : "<<  m->getNumFaces() <<std::endl;
      std::cout<<"num Aretes : "<<  m->getNumEdges() <<std::endl;
      std::cout<<"num Sommets : "<<  m->getNumVertices() <<std::endl;

      //////////level Set

                  int tag2 =1;
                  const double pt2[3]={1.0,0.0,1.0};
                  const double dir2[3]={0.0,1.0,0.0};
                  double r2 = 0.3;

                  gLevelset *_ls_cut1 = new gLevelsetGenCylinder(pt2,dir2,r2,tag2);
                  GModel *m1 = new GModel();
                  m1 = m->buildCutGModel(_ls_cut1);



                  int tag =2;
                  const double pt1[3]={0.0,0.0,0.0};
                  const double dir1[3]={0.0,1.0,0.0};
                  double r = 0.3;

                  gLevelset *_ls_cut = new gLevelsetGenCylinder(pt1,dir1,r,tag);
                  GModel *m2 = new GModel();
                  m2 = m1->buildCutGModel(_ls_cut);


  /*
                  double a(0.), b(0), c(1), d(-0.5);
                  int n(2); // pour level set
                  gLevelset *_ls = new gLevelsetPlane(a, b, c, d, n);
                  m = m->buildCutGModel(_ls);
  */

  //m2->writeMSH("PP_cut.msh",2.1,false,false);

  // stop gmsh
  //GmshFinalize();


  /*
      GmshInitialize(argc, argv);

      GModel *m = new GModel;


      STEPControl_Reader readerStep;
      readerStep.ReadFile("P.STEP");
      readerStep.TransferRoots ();
      TopoDS_Shape shape = readerStep.OneShape();

      FILE * pFile;
      pFile = fopen ("myfile0.txt","wt");


      std::vector<TopoDS_Vertex > Vertices;

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

          TopoDS_Vertex Vertex = TopoDS::Vertex(exp0.Current());

          int test =0;
          for (int ii = 0; ii <  Vertices.size(); ii++){

              if(Vertex.IsSame(Vertices[ii])){

              test =1;

              }
          }

          if(test ==0){

          Vertices.push_back(Vertex);
          }

      }

      fprintf (pFile, "Points\n");
      fprintf (pFile, "%d\n", Vertices.size());

      int tagV = 1;

      for (int ii = 0; ii <  Vertices.size(); ii++){

          gp_Pnt Pt = BRep_Tool::Pnt(Vertices[ii]);
          fprintf (pFile, "%d %lf %lf %lf %d\n",ii,Pt.X(),Pt.Y(),Pt.Z(), tagV);
      }


      std::vector<TopoDS_Edge > Edges;

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

          TopoDS_Edge Edge = TopoDS::Edge(exp0.Current());

          int test =0;
          for (int ii = 0; ii <  Edges.size(); ii++){

              if(Edge.IsSame(Edges[ii])){
              test =1;
              }
          }

          if(test ==0){
          Edges.push_back(Edge);
          }

      }

      fprintf (pFile, "Edges\n");
      fprintf (pFile, "%d\n", Edges.size());

      int tagE =1;
      for (int ii = 0; ii <  Edges.size(); ii++){

          int idv1, idv2;

          TopExp_Explorer exp1;
          int id12 =1;
          for(exp1.Init(Edges[ii], TopAbs_VERTEX); exp1.More(); exp1.Next()){

              TopoDS_Vertex Vertex = TopoDS::Vertex(exp1.Current());

              for (int jj = 0; jj <  Vertices.size(); jj++){

                  if(Vertex.IsSame(Vertices[jj])&&id12 == 1){

                      idv1 =jj;

                  }
                  if(Vertex.IsSame(Vertices[jj])&&id12 == 2){

                      idv2 =jj;

                  }
              }
              id12++;
          }

          fprintf (pFile, "%d %d %d %d\n",ii,idv1,idv2,tagE);
      }

      std::vector<TopoDS_Face > Faces;

      int idF=0;

      std::vector<std::vector<double> > Parametres;

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

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

          Faces.push_back(Face);

          std::vector<double> ParamF;

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

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

              std::cout<<" Plane Plane Plane : "<<std::endl;
              Handle(Geom_Plane) aPlane = Handle(Geom_Plane)::DownCast(Surface);
              gp_Pln agpPlane = aPlane->Pln();
              gp_Ax1 normale = agpPlane.Axis();

              ParamF.push_back(normale.Location().X());
              ParamF.push_back(normale.Location().Y());
              ParamF.push_back(normale.Location().Z());
              ParamF.push_back(normale.Direction().X());
              ParamF.push_back(normale.Direction().Y());
              ParamF.push_back(normale.Direction().Z());

          }

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

              std::cout<<" Cylindre Cylindre Cylindre : "<<std::endl;
              Handle(Geom_CylindricalSurface) aCylindre = Handle(Geom_CylindricalSurface)::DownCast(Surface);
              gp_Cylinder agpCylindre = aCylindre->Cylinder();
              gp_Ax1 axe = agpCylindre.Axis();
              double R = agpCylindre.Radius();
              ParamF.push_back(axe.Location().X());
              ParamF.push_back(axe.Location().Y());
              ParamF.push_back(axe.Location().Z());
              ParamF.push_back(axe.Direction().X());
              ParamF.push_back(axe.Direction().Y());
              ParamF.push_back(axe.Direction().Z());
              ParamF.push_back(R);
          }

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

              std::cout<<" Cone Cone Cone : "<<std::endl;
              Handle(Geom_ConicalSurface) aCone = Handle(Geom_ConicalSurface)::DownCast(Surface);
              gp_Cone agpCone = aCone->Cone();
              gp_Ax1 axe = agpCone.Axis();
              double Rf = agpCone.RefRadius();
              double Angle = agpCone.SemiAngle();
              ParamF.push_back(agpCone.Location().X());
              ParamF.push_back(agpCone.Location().Y());
              ParamF.push_back(agpCone.Location().Z());
              ParamF.push_back(axe.Direction().X());
              ParamF.push_back(axe.Direction().Y());
              ParamF.push_back(axe.Direction().Z());
              ParamF.push_back(Rf);
              ParamF.push_back(Angle);

          }

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

              std::cout<<" Tore Tore Tore : "<<std::endl;
              Handle(Geom_ToroidalSurface) aTore = Handle(Geom_ToroidalSurface)::DownCast(Surface);
              gp_Torus agpTore = aTore->Torus();
              gp_Ax1 axe = agpTore.Axis();

              double gR = agpTore.MajorRadius();
              double pR = agpTore.MinorRadius();

              ParamF.push_back(agpTore.Location().X());
              ParamF.push_back(agpTore.Location().Y());
              ParamF.push_back(agpTore.Location().Z());
              ParamF.push_back(axe.Direction().X());
              ParamF.push_back(axe.Direction().Y());
              ParamF.push_back(axe.Direction().Z());
              ParamF.push_back(gR);
              ParamF.push_back(pR);


          }

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

              std::cout<<" Sphere Sphere Sphere : "<<std::endl;
              Handle(Geom_SphericalSurface) aSphere = Handle(Geom_SphericalSurface)::DownCast(Surface);
              gp_Sphere agpSphere = aSphere->Sphere();
              //gp_Ax3 axe = agpSphere.Position();
              double R = agpSphere.Radius();

              ParamF.push_back(agpSphere.Location().X());
              ParamF.push_back(agpSphere.Location().Y());
              ParamF.push_back(agpSphere.Location().Z());
              //ParamF.push_back(axe.Direction().X());
              //ParamF.push_back(axe.Direction().Y());
              //ParamF.push_back(axe.Direction().Z());
              ParamF.push_back(R);
  dindin.step
          }

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

              std::cout<<" Revolution Revolution Revolution  : "<<std::endl;
              Handle(Geom_SurfaceOfRevolution) aRevolSurf = Handle(Geom_SurfaceOfRevolution)::DownCast(Surface);

              Handle(Geom_Curve) Courbe = aRevolSurf->BasisCurve();
              gp_Ax1 Axe = aRevolSurf->Axis();
              gp_Pnt ptPosition = aRevolSurf->Location();

          }

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

              Handle(Geom_SurfaceOfLinearExtrusion) aExtrudSurf = Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(Surface);

              Handle(Geom_Curve) Courbe = aExtrudSurf->BasisCurve();
              gp_Dir dir = aExtrudSurf->Direction();

          }

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

              Handle(Geom_BezierSurface) aBezierSurf = Handle(Geom_BezierSurface)::DownCast(Surface);

              int nbCPntU =  aBezierSurf->NbUPoles();
              int nbCPntV =  aBezierSurf->NbVPoles();

              TColgp_Array2OfPnt PtsControl(0,nbCPntU,0,nbCPntV);
              aBezierSurf->Poles(PtsControl);
              int Udeg =  aBezierSurf->UDegree();
              int Vdeg =  aBezierSurf->VDegree();

          }

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

              Handle(Geom_BSplineSurface) aBSplineSurf = Handle(Geom_BSplineSurface)::DownCast(Surface);

              int nbCPntU =  aBSplineSurf->NbUPoles();
              int nbCPntV =  aBSplineSurf->NbVPoles();

              int nbUknots = aBSplineSurf->NbUKnots();
              int nbVknots = aBSplineSurf->NbVKnots();

              TColStd_Array1OfReal Ku(0,nbUknots-1);
              TColStd_Array1OfReal Kv(0,nbVknots-1);


              TColgp_Array2OfPnt PtsControl(0,nbCPntU-1,0,nbCPntV-1);
              aBSplineSurf->Poles(PtsControl);
              int Udeg =  aBSplineSurf->UDegree();
              int Vdeg =  aBSplineSurf->VDegree();

              }
          Parametres.push_back(ParamF);
          idF++;

      }

      fprintf (pFile, "Faces\n");
      fprintf (pFile, "%d", Faces.size());

      int tagF =1;
      int typeGeomSuface;
      for (int ii = 0; ii <  Faces.size(); ii++){

          fprintf (pFile, "\n");

          TopoDS_Face Face = Faces[ii];

          std::vector<double> ParamF = Parametres[ii];

          fprintf (pFile,"%d ", ii);

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

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

              typeGeomSuface =100;
              fprintf (pFile,"Plan ");
          }

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

              typeGeomSuface =101;
              fprintf (pFile,"Cylindre ");
          }

          for (int jj = 0; jj <  ParamF.size(); jj++){

              fprintf (pFile,"%lf ", ParamF[jj]);

              }

          fprintf (pFile,"%d\n", tagF);

          TopExp_Explorer exp1;

          int nbAreteFace = 0;
          for(exp1.Init(Face, TopAbs_EDGE); exp1.More(); exp1.Next()){

             nbAreteFace = nbAreteFace + 1;
          }

          fprintf (pFile,"%d ", nbAreteFace);

          for(exp1.Init(Face, TopAbs_EDGE); exp1.More(); exp1.Next()){

              TopoDS_Edge Edge = TopoDS::Edge(exp1.Current());

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

                  if(Edge.IsSame(Edges[jj])){

                      fprintf (pFile,"%d ", jj);

                  }

              }

          }

      }




      fclose (pFile);


  */



  /*
  TopExp_Explorer exp0, exp1, exp2, exp3, exp4, exp5, exp01;

      TopTools_IndexedMapOfShape fmap, emap, vmap, somap, shmap, wmap;

    for(exp0.Init(shape, TopAbs_SOLID); exp0.More(); exp0.Next()){
      TopoDS_Solid solid = TopoDS::Solid(exp0.Current());




            for(exp2.Init(solid, TopAbs_FACE); exp2.More(); exp2.Next()){
              TopoDS_Face face = TopoDS::Face(exp2.Current().Composed(solid.Orientation()));

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

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

                  std::cout<<" Plane Plane Plane : "<<std::endl;
                  Handle(Geom_Plane) aPlane = Handle(Geom_Plane)::DownCast(Surface);
                  gp_Pln agpPlane = aPlane->Pln();
                  gp_Ax1 norm = agpPlane.Axis();
                  gp_Dir dir = norm.Direction();
                  //std::cout<<" X_Normal : "<<  dir.X() <<" Y_Normal : "<<  dir.Y() <<" Z_Normal : "<<  dir.Z()<<std::endl;

              }

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

                  std::cout<<" Cylindre Cylindre Cylindre : "<<std::endl;
                  Handle(Geom_CylindricalSurface) aCylindre = Handle(Geom_CylindricalSurface)::DownCast(Surface);
                  gp_Cylinder agpCylindre = aCylindre->Cylinder();
                  gp_Ax1 norm = agpCylindre.Axis();
                  gp_Dir dir = norm.Direction();
                  double R = agpCylindre.Radius();
                  //std::cout<<" Rayon Cylindre : "<<  R <<std::endl;
              }

              if (Surface->DynamicType() == STANDARD_TYPE(Geom_ConicalSurface))
              {
                  std::cout<<" Cone Cone Cone : "<<std::endl;
                  Handle(Geom_ConicalSurface) aCone = Handle(Geom_ConicalSurface)::DownCast(Surface);
                  gp_Cone agpCone = aCone->Cone();
                  gp_Ax1 norm = agpCone.Axis();
                  gp_Dir dir = norm.Direction();
                  double Rf = agpCone.RefRadius();
                  gp_Ax3 Position = agpCone.Position();
                  double Angle = agpCone.SemiAngle();

                  //std::cout<<" Rayon Cone : "<<  Rf <<std::endl;
                  //std::cout<<" Angle : "<<  Angle <<std::endl;
              }

              if (Surface->DynamicType() == STANDARD_TYPE(Geom_ToroidalSurface))
              {
                  std::cout<<" Tore Tore Tore : "<<std::endl;
                  Handle(Geom_ToroidalSurface) aTore = Handle(Geom_ToroidalSurface)::DownCast(Surface);
                  gp_Torus agpTore = aTore->Torus();
                  gp_Ax1 norm = agpTore.Axis();

                  double gR = agpTore.MajorRadius();
                  double pR = agpTore.MinorRadius();
                  gp_Ax3 Position = agpTore.Position();

                  //std::cout<<" Major Radius : "<<  gR <<" Minor Radius : "<<  pR <<std::endl;

              }

              if (Surface->DynamicType() == STANDARD_TYPE(Geom_SphericalSurface))
              {
                  std::cout<<" Sphere Sphere Sphere : "<<std::endl;
                  Handle(Geom_SphericalSurface) aSphere = Handle(Geom_SphericalSurface)::DownCast(Surface);
                  gp_Sphere agpSphere = aSphere->Sphere();

                  gp_Pnt ptInsertion = agpSphere.Location();
                  gp_Ax3 Axe = agpSphere.Position();
                  double R = agpSphere.Radius();

                  //std::cout<<" Rayon : "<<  R <<" X_Insertion : "<<  ptInsertion.X() <<" Y_Insertion : "<<  ptInsertion.Y() <<" Z_Insertion : "<<  ptInsertion.Z()<<std::endl;

              }

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

                  std::cout<<" Revolution Revolution Revolution  : "<<std::endl;
                  Handle(Geom_SurfaceOfRevolution) aRevolSurf = Handle(Geom_SurfaceOfRevolution)::DownCast(Surface);

                  Handle(Geom_Curve) Courbe = aRevolSurf->BasisCurve();
                  gp_Ax1 Axe = aRevolSurf->Axis();
                  gp_Pnt ptPosition = aRevolSurf->Location();

                  //std::cout<<" X_Position : "<<  ptPosition.X() <<" Y_Position : "<<  ptPosition.Y() <<" Z_Position : "<<  ptPosition.Z()<<std::endl;

              }

              if (Surface->DynamicType() == STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion))
              {
                  Handle(Geom_SurfaceOfLinearExtrusion) aExtrudSurf = Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(Surface);

                  Handle(Geom_Curve) Courbe = aExtrudSurf->BasisCurve();
                  gp_Dir dir = aExtrudSurf->Direction();

              }

              if (Surface->DynamicType() == STANDARD_TYPE(Geom_BezierSurface))
              {
                  Handle(Geom_BezierSurface) aBezierSurf = Handle(Geom_BezierSurface)::DownCast(Surface);

                  int nbCPntU =  aBezierSurf->NbUPoles();
                  int nbCPntV =  aBezierSurf->NbVPoles();


                  TColgp_Array2OfPnt PtsControl(0,nbCPntU,0,nbCPntV);
                  aBezierSurf->Poles(PtsControl);
                  int Udeg =  aBezierSurf->UDegree();
                  int Vdeg =  aBezierSurf->VDegree();

              }

              if (Surface->DynamicType() == STANDARD_TYPE(Geom_BSplineSurface))
              {
                  Handle(Geom_BSplineSurface) aBSplineSurf = Handle(Geom_BSplineSurface)::DownCast(Surface);

                  int nbCPntU =  aBSplineSurf->NbUPoles();
                  int nbCPntV =  aBSplineSurf->NbVPoles();

                  int nbUknots = aBSplineSurf->NbUKnots();
                  int nbVknots = aBSplineSurf->NbVKnots();

                  TColStd_Array1OfReal Ku(0,nbUknots-1);
                  TColStd_Array1OfReal Kv(0,nbVknots-1);


                  TColgp_Array2OfPnt PtsControl(0,nbCPntU-1,0,nbCPntV-1);
                  aBSplineSurf->Poles(PtsControl);
                  int Udeg =  aBSplineSurf->UDegree();
                  int Vdeg =  aBSplineSurf->VDegree();

              }





              std::vector<TopoDS_Face > Faces;


                //for(exp3.Init(exp2.Current(), TopAbs_WIRE); exp3.More(); exp3.Next()){
                  //TopoDS_Wire wire = TopoDS::Wire(exp3.Current().Composed(face.Orientation()));


                    for(exp4.Init(exp2.Current(), TopAbs_EDGE); exp4.More(); exp4.Next()){
                      TopoDS_Edge edge = TopoDS::Edge(exp4.Current().Composed(face.Orientation()));

                          int numface = 0;
                          for(exp01.Init(solid, TopAbs_FACE); exp01.More(); exp01.Next()){
                          TopoDS_Face face01 = TopoDS::Face(exp01.Current());


                              if(face01.IsSame(face)==false){
                                  for (TopExp_Explorer expEdge(face01, TopAbs_EDGE); expEdge.More(); expEdge.Next())
                                  {
                                      TopoDS_Edge edge01 = TopoDS::Edge(expEdge.Current()); // this is the edge on the face
                                      if(edge01.IsSame(edge))// same edge, but maybe different orientation
                                      {

                                          int test =0;
                                          for (int ii = 0; ii <  Faces.size(); ii++){

                                              TopoDS_Face facetest = Faces[ii];

                                              if(facetest == face01){

                                                  test =1;

                                                  }


                                              }

                                              if(test ==0){

                                                  Faces.push_back(face01);
                                                  }


                                          //Faces.push_back(&face01);

                                          //if(E1.IsEqual(E2))sense = true; // same orientation too
                                          //break; // same edge ( ignore the face's other edges )
                                      }
                                  }
                              }

                          }







                    }

                    std::cout<<" Nombre de faces qui limitent cette face :  "<< Faces.size() <<std::endl;

                    for(exp4.Init(exp2.Current(), TopAbs_EDGE); exp4.More(); exp4.Next()){
                      TopoDS_Edge edge = TopoDS::Edge(exp4.Current());

                      TColgp_Array1OfPnt Array(1,2);
                      int cc = 1;

                      double pfirst, plast;

                      Handle (Geom_Curve) Courbe = BRep_Tool::Curve(edge,pfirst,plast);
                      std::cout<<" pfist :  "<< pfirst <<" plast :  "<< plast <<std::endl;

  /*
                      if (Courbe->DynamicType() == STANDARD_TYPE(Geom_Line)){

                          for(exp5.Init(exp4.Current(), TopAbs_VERTEX); exp5.More(); exp5.Next()){
                          TopoDS_Vertex Vertex = TopoDS::Vertex(exp5.Current());

                          gp_Pnt Point = BRep_Tool::Pnt(Vertex);

                          Array.SetValue(cc,Point);
                          cc++;

                          }

                          }


                      if (Courbe->DynamicType() != STANDARD_TYPE(Geom_Line)){

                          for(exp5.Init(exp4.Current(), TopAbs_VERTEX); exp5.More(); exp5.Next()){
                          TopoDS_Vertex Vertex = TopoDS::Vertex(exp5.Current());

                          gp_Pnt Point = BRep_Tool::Pnt(Vertex);

                          Array.SetValue(cc,Point);
                          cc++;

                          }

                          }

  */
  /*
                    }



                //}


            }



    }






      // stop gmsh
      GmshFinalize();

  */

  /*
    GmshInitialize(argc, argv);

    //****** lire le fichier

      GModel *m = new GModel();

      //m->readGEO("p.geo");
      m->readMSH("p.msh");

      std::cout<<"num Faces : "<<  m->getNumFaces() <<std::endl;
      std::cout<<"num Aretes : "<<  m->getNumEdges() <<std::endl;
      std::cout<<"num Sommets : "<<  m->getNumVertices() <<std::endl;



      int tag =2;



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

     double R = 1;
     double r = 0.025;
     double ir = 0.02;
     double a = 0.20;
     double b = 0.40;
     double c = 0.05;
     double d = 0.10;
     double e = 0.50;

      gp_Circ cc = gp_Circ(gp_Ax2(gp_Pnt(0.,0.,0.),gp_Dir(0.,0.,1.)), R);
   TopoDS_Edge Ec = BRepBuilderAPI_MakeEdge(cc);
   TopoDS_Wire Wc = BRepBuilderAPI_MakeWire(Ec);
   TopoDS_Face F = BRepBuilderAPI_MakeFace(gp_Pln(gp::XOY()),Wc);


      double start =0;
      double end = 0;
      int nb =0;
      double dir = 0;

      dir =1;
      start = b/2 + c;
      end = start + e;
      int nb_y = floor((end-start)/((2*r+ir)*(sqrt(3)/2)));
      start = b/2 + c + e/2 -(nb_y*(2*r+ir)*(sqrt(3)/2)/2);
      double y_0 = start;

      start = -sqrt( (R-d)*(R-d) - (y_0)*(y_0) );
      end = sqrt( (R-d)*(R-d) - (y_0)*(y_0) );
      int nb_x = Floor((end-start)/(2*r + ir));
      start = -(nb_x*(2*r + ir)/2);
      double x_0 = start;

      //////////////////
      // Bloc 1
      //////////////////

      int i = 0;
      double x = x_0;
      double y = y_0;

   //createHole
   for (int ii = 1; ii <=nb_x+1; ii++){

       if(Sqrt (x*x + y*y) < (R - r - d)){
       gp_Circ c1 = gp_Circ(gp_Ax2(gp_Pnt(x,y,0.),gp_Dir(0.,0.,1.)), r);
          TopoDS_Edge Ec1 = BRepBuilderAPI_MakeEdge(c1);
          TopoDS_Wire Wc1 = BRepBuilderAPI_MakeWire(Ec1);
          F = BRepBuilderAPI_MakeFace(F,Wc1);


                  //////////level Set
                  const double pt1[3]={x,y,0.0};
                  const double dir1[3]={0.0,0.0,1.0};

                  gLevelset *_ls_cut = new gLevelsetGenCylinder(pt1,dir1,r,tag);
                  m = m->buildCutGModel(_ls_cut);


       }

       if (dir == 0){
              x -= 2*r + ir;
       }

          if (dir == 1){
              x += 2*r + ir;
          }
   }

   for (int ii = 2; ii <=nb_y; ii++){

       i++;
       if (i%2 == 0){
              x = x_0;
       }

          if (i%2 == 1){
              x = x_0 - ((2*r + ir)/2);
          }
          y += (2*r + ir)*(Sqrt(3)/2);

          //createHole
          for (int jj = 1; jj <=nb_x+1; jj++){

              if(Sqrt (x*x + y*y) < (R - r - d)){
              gp_Circ c1 = gp_Circ(gp_Ax2(gp_Pnt(x,y,0.),gp_Dir(0.,0.,1.)), r);
              TopoDS_Edge Ec1 = BRepBuilderAPI_MakeEdge(c1);
              TopoDS_Wire Wc1 = BRepBuilderAPI_MakeWire(Ec1);
              F = BRepBuilderAPI_MakeFace(F,Wc1);

                              //////////level Set
                  const double pt1[3]={x,y,0.0};
                  const double dir1[3]={0.0,0.0,1.0};

                  gLevelset *_ls_cut = new gLevelsetGenCylinder(pt1,dir1,r,tag);
                  m = m->buildCutGModel(_ls_cut);
              }

              if (dir == 0){
                  x -= 2*r + ir;
              }

              if (dir == 1){
                  x += 2*r + ir;
              }
          }
          /////////////


   }

      //////////////////
      //* Bloc 4
      //////////////////

       i = 0;
       x = x_0;
       y = -y_0;

   //createHole
   for (int ii = 1; ii <=nb_x+1; ii++){

       if(Sqrt (x*x + y*y) < (R - r - d)){
       gp_Circ c1 = gp_Circ(gp_Ax2(gp_Pnt(x,y,0.),gp_Dir(0.,0.,1.)), r);
          TopoDS_Edge Ec1 = BRepBuilderAPI_MakeEdge(c1);
          TopoDS_Wire Wc1 = BRepBuilderAPI_MakeWire(Ec1);
          F = BRepBuilderAPI_MakeFace(F,Wc1);

                          //////////level Set
                  const double pt1[3]={x,y,0.0};
                  const double dir1[3]={0.0,0.0,1.0};
                  double r1=r;

                  gLevelset *_ls_cut = new gLevelsetGenCylinder(pt1,dir1,r1,tag);
                  m = m->buildCutGModel(_ls_cut);

       }

       if (dir == 0){
              x -= 2*r + ir;
       }

          if (dir == 1){
              x += 2*r + ir;
          }
   }

   for (int ii = 2; ii <=nb_y; ii++){

       i++;
       if (i%2 == 0){
              x = x_0;
       }

          if (i%2 == 1){
              x = x_0 - ((2*r + ir)/2);
          }
          y -= (2*r + ir)*(Sqrt(3)/2);

          //createHole
          for (int jj = 1; jj <=nb_x+1; jj++){

              if(Sqrt (x*x + y*y) < (R - r - d)){
              gp_Circ c1 = gp_Circ(gp_Ax2(gp_Pnt(x,y,0.),gp_Dir(0.,0.,1.)), r);
              TopoDS_Edge Ec1 = BRepBuilderAPI_MakeEdge(c1);
              TopoDS_Wire Wc1 = BRepBuilderAPI_MakeWire(Ec1);
              F = BRepBuilderAPI_MakeFace(F,Wc1);

                              //////////level Set
                  const double pt1[3]={x,y,0.0};
                  const double dir1[3]={0.0,0.0,1.0};
                  double r1=r;

                  gLevelset *_ls_cut = new gLevelsetGenCylinder(pt1,dir1,r1,tag);
                  m = m->buildCutGModel(_ls_cut);

              }

              if (dir == 0){
                  x -= 2*r + ir;
              }

              if (dir == 1){
                  x += 2*r + ir;
              }
          }
          /////////////

   }



  ////////////////////////
  // bloc 2 et 3
  ////////////////////////
  start = -b/2;
  end = b/2;
  nb_y = Floor((end-start)/((2*r + ir)*(Sqrt(3)/2)));
  start = -(nb_y*(2*r + ir)*(Sqrt(3)/2)/2);
  y_0 = start;

  start = a/2;
  end = R-d;
  nb_x = Floor((end-start)/(2*r + ir));
  x_0 = start + r;
  y_0 += r;
  ////////////////////////
  // bloc 2
  ////////////////////////
  i = 0;
  x = x_0;
  y = y_0;

  //createHole
          for (int jj = 1; jj <=nb_x+1; jj++){

              if(Sqrt (x*x + y*y) < (R - r - d)){
              gp_Circ c1 = gp_Circ(gp_Ax2(gp_Pnt(x,y,0.),gp_Dir(0.,0.,1.)), r);
              TopoDS_Edge Ec1 = BRepBuilderAPI_MakeEdge(c1);
              TopoDS_Wire Wc1 = BRepBuilderAPI_MakeWire(Ec1);
              F = BRepBuilderAPI_MakeFace(F,Wc1);

                              //////////level Set
                  const double pt1[3]={x,y,0.0};
                  const double dir1[3]={0.0,0.0,1.0};
                  double r1=r;

                  gLevelset *_ls_cut = new gLevelsetGenCylinder(pt1,dir1,r1,tag);
                  m = m->buildCutGModel(_ls_cut);

              }

              if (dir == 0){
                  x -= 2*r + ir;
              }

              if (dir == 1){
                  x += 2*r + ir;
              }
          }
          /////////////

  for(int ii = 2; ii <=nb_y; ii++){
    i++;
    if(i%2 == 0){
    x = x_0;
    }
    if (i%2 == 1){
    x = x_0 + ((2*r + ir)/2);
    }
    y += (2*r + ir)*(Sqrt(3)/2);
    //createHole
          for (int jj = 1; jj <=nb_x+1; jj++){

              if(Sqrt (x*x + y*y) < (R - r - d)){
              gp_Circ c1 = gp_Circ(gp_Ax2(gp_Pnt(x,y,0.),gp_Dir(0.,0.,1.)), r);
              TopoDS_Edge Ec1 = BRepBuilderAPI_MakeEdge(c1);
              TopoDS_Wire Wc1 = BRepBuilderAPI_MakeWire(Ec1);
              F = BRepBuilderAPI_MakeFace(F,Wc1);



              }

              if (dir == 0){
                  x -= 2*r + ir;
              }

              if (dir == 1){
                  x += 2*r + ir;
              }
          }
          /////////////
  }


  ////////////////////////
  // bloc 3
  ////////////////////////
  dir = 0;
  i = 0;
  x = -x_0;
  y = y_0;

  //createHole
          for (int jj = 1; jj <=nb_x+1; jj++){

              if(Sqrt (x*x + y*y) < (R - r - d)){
              gp_Circ c1 = gp_Circ(gp_Ax2(gp_Pnt(x,y,0.),gp_Dir(0.,0.,1.)), r);
              TopoDS_Edge Ec1 = BRepBuilderAPI_MakeEdge(c1);
              TopoDS_Wire Wc1 = BRepBuilderAPI_MakeWire(Ec1);
              F = BRepBuilderAPI_MakeFace(F,Wc1);

                              //////////level Set
                  const double pt1[3]={x,y,0.0};
                  const double dir1[3]={0.0,0.0,1.0};
                  double r1=r;

                  gLevelset *_ls_cut = new gLevelsetGenCylinder(pt1,dir1,r1,tag);
                  m = m->buildCutGModel(_ls_cut);

              }

              if (dir == 0){
                  x -= 2*r + ir;
              }

              if (dir == 1){
                  x += 2*r + ir;
              }
          }
          /////////////

  for(int ii = 2; ii <=nb_y; ii++){
    i++;
    if(i%2 == 0){
    x = -x_0;
    }
    if (i%2 == 1){
    x = -x_0 - ((2*r + ir)/2);
    }
    y += (2*r + ir)*(Sqrt(3)/2);
    //createHole
          for (int jj = 1; jj <=nb_x+1; jj++){

              if(Sqrt (x*x + y*y) < (R - r - d)){
              gp_Circ c1 = gp_Circ(gp_Ax2(gp_Pnt(x,y,0.),gp_Dir(0.,0.,1.)), r);
              TopoDS_Edge Ec1 = BRepBuilderAPI_MakeEdge(c1);
              TopoDS_Wire Wc1 = BRepBuilderAPI_MakeWire(Ec1);
              F = BRepBuilderAPI_MakeFace(F,Wc1);



              }

              if (dir == 0){
                  x -= 2*r + ir;
              }

              if (dir == 1){
                  x += 2*r + ir;
              }
          }
          /////////////
  }



      m->writeMSH("p_cut.msh",2.1,false,false);
      //GModel::setCurrent(m);


   double H = 0.1;
   TopoDS_Shape S = BRepPrimAPI_MakePrism(F,gp_Vec(0.,0.,H));






      STEPControl_Writer writer;
   writer.Transfer(S, STEPControl_ManifoldSolidBrep);
   writer.Write("plaque.step");

  */

  /*

   GmshInitialize(argc, argv);

   GModel *m = new GModel;

   m->readGEO("Chemise.geo");
   m->readMSH("Chemise.msh");

   std::cout<<"num Faces : "<<  m->getNumFaces() <<std::endl;
   std::cout<<"num Aretes : "<<  m->getNumEdges() <<std::endl;
   std::cout<<"num Sommets : "<<  m->getNumVertices() <<std::endl;


   BRep_Builder B;
   TopoDS_Shell Sh;

   B.MakeShell(Sh);

   int cont = 0;


   std::vector<gp_Pnt > Pnts;
   std::map<int, gp_Pnt> mapPnts;

   int nbface=0;///****
   std::vector<int > nb_areteparface; ///****
   std::vector<int > nb_ptspararete; ///****



   for( GModel::eiter ite = m->firstEdge(); ite != m->lastEdge(); ++ite){

          std::cout<<"Tag    : "<<  (*ite)->tag() <<std::endl;

          int nbe = (*ite)->getNumMeshVertices();
          std::cout<<"nbe : "<<  nbe <<std::endl;


          GVertex *V1 = (*ite)->getBeginVertex();
          int nbv1 = V1->getNumMeshVertices();
          std::cout<<"nbv1 : "<<  nbv1 <<std::endl;
          MVertex *Nodv1 = V1->mesh_vertices[0];
          std::cout<<" Num : "<<  Nodv1->getNum() <<std::endl;
          std::cout<<" xNv1 : "<<  Nodv1->x()<<" yNv1 : "<<  Nodv1->y()<<" zNv1 : "<<  Nodv1->z() <<std::endl;

          //



          Pnts.push_back(gp_Pnt(Nodv1->x(),Nodv1->y(),Nodv1->z())); ///*****************
          mapPnts.insert(std::pair<int, gp_Pnt>(Nodv1->getNum(),gp_Pnt(Nodv1->x(),Nodv1->y(),Nodv1->z())));

          for (int ii = 0; ii < nbe; ii++){

          MVertex *Node = (*ite)->mesh_vertices[ii];
          std::cout<<" Num : "<<  Node->getNum() <<std::endl;
          std::cout<<" xNe : "<<  Node->x()<<" yNe : "<<  Node->y()<<" zNe : "<<  Node->z() <<std::endl;



          mapPnts.insert(std::pair<int, gp_Pnt>(Node->getNum(),gp_Pnt(Node->x(),Node->y(),Node->z())));
          }

          GVertex *V2 = (*ite)->getEndVertex();
          int nbv2 = V2->getNumMeshVertices();
          std::cout<<"nbv2 : "<<  nbv2 <<std::endl;
          MVertex *Nodv2 = V2->mesh_vertices[0];
          std::cout<<" Num : "<<  Nodv2->getNum() <<std::endl;
          std::cout<<" xNv2 : "<<  Nodv2->x()<<" yNv2 : "<<  Nodv2->y()<<" zNv2 : "<<  Nodv2->z() <<std::endl;

          mapPnts.insert(std::pair<int, gp_Pnt>(Nodv2->getNum(),gp_Pnt(Nodv2->x(),Nodv2->y(),Nodv2->z())));

      }

      //getNumEdges()

        FILE * pFile0;
        pFile0 = fopen("Chemise0.surf","wt");

        fprintf(pFile0, "%d \n",mapPnts.size());

      std::map<int,gp_Pnt>::iterator it0;
      for ( it0=mapPnts.begin() ; it0 != mapPnts.end(); it0++ )
      {
          fprintf(pFile0, "%lf %lf %lf\n", (*it0).second.X(), (*it0).second.Y(), (*it0).second.Z());
      }

      fprintf(pFile0, "%d \n",m->getNumEdges());

      //std::map<int,gp_Pnt>::iterator itmap;

      for( GModel::eiter ite = m->firstEdge(); ite != m->lastEdge(); ++ite){


          int nbe = (*ite)->getNumMeshVertices();
          fprintf(pFile0, "%d %d\n",(*ite)->tag()-1,nbe+2);


          GVertex *V1 = (*ite)->getBeginVertex();
          MVertex *Nodv1 = V1->mesh_vertices[0];

          fprintf(pFile0, "%d ",std::distance(mapPnts.begin(), mapPnts.find(Nodv1->getNum())));

          for (int ii = 0; ii < nbe; ii++){

          MVertex *Node = (*ite)->mesh_vertices[ii];

          fprintf(pFile0, "%d ",std::distance(mapPnts.begin(), mapPnts.find(Node->getNum())));
          }

          GVertex *V2 = (*ite)->getEndVertex();

          MVertex *Nodv2 = V2->mesh_vertices[0];

          fprintf(pFile0, "%d \n",std::distance(mapPnts.begin(), mapPnts.find(Nodv2->getNum())));
      }

      fprintf(pFile0, "%d \n",m->getNumFaces());

      for (GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
      {
          fprintf(pFile0, "%d %d \n",(*it)->tag(), 0);

          std::list<GEdge*> edgs = (*it)->edges();
          fprintf(pFile0, "%d \n",edgs.size());
          for( std::list<GEdge*>::const_iterator ite = edgs.begin(); ite != edgs.end(); ite++)
          {
              fprintf(pFile0, "%d ",(*ite)->tag()-1);
          }
          fprintf(pFile0, "\n");

      }

      fclose(pFile0);

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


    std::vector<std::vector<gp_Pnt > > ptsfaces;

   for (GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){

       //if(nbface<=4)
       {

      nbface=nbface+1;///*************

      GeomPlate_BuildPlateSurface BPSurf(3,10,3);

       int nb = (*it)->getNumMeshVertices();
       std::cout<<"nb : "<<  nb <<std::endl;

       int nbareteface=0;///*************

       std::list<GEdge*> edgs = (*it)->edges();

      for( std::list<GEdge*>::const_iterator ite = edgs.begin(); ite != edgs.end(); ite++){

          nbareteface=nbareteface+1;///*******************

          std::cout<<"Tag    : "<<  (*ite)->tag() <<std::endl;

          int nbe = (*ite)->getNumMeshVertices();
          std::cout<<"nbe : "<<  nbe <<std::endl;

          TColgp_Array1OfPnt Array(1,nbe+2);
          int cc = 1;

          GVertex *V1 = (*ite)->getBeginVertex();
          int nbv1 = V1->getNumMeshVertices();
          std::cout<<"nbv1 : "<<  nbv1 <<std::endl;
          MVertex *Nodv1 = V1->mesh_vertices[0];
          std::cout<<" Num : "<<  Nodv1->getNum() <<std::endl;
          std::cout<<" xNv1 : "<<  Nodv1->x()<<" yNv1 : "<<  Nodv1->y()<<" zNv1 : "<<  Nodv1->z() <<std::endl;

          //



          Array.SetValue(cc,gp_Pnt(Nodv1->x(),Nodv1->y(),Nodv1->z()));
          cc++;

          Pnts.push_back(gp_Pnt(Nodv1->x(),Nodv1->y(),Nodv1->z())); ///*****************
          //mapPnts.insert(std::pair<int, gp_Pnt>(Nodv1->getNum(),gp_Pnt(Nodv1->x(),Nodv1->y(),Nodv1->z())));

          for (int ii = 0; ii < nbe; ii++){

          MVertex *Node = (*ite)->mesh_vertices[ii];
          std::cout<<" Num : "<<  Node->getNum() <<std::endl;
          std::cout<<" xNe : "<<  Node->x()<<" yNe : "<<  Node->y()<<" zNe : "<<  Node->z() <<std::endl;

          Array.SetValue(cc,gp_Pnt(Node->x(),Node->y(),Node->z()));
          cc++;
          Pnts.push_back(gp_Pnt(Node->x(),Node->y(),Node->z())); ///*****************
          //mapPnts.insert(std::pair<int, gp_Pnt>(Node->getNum(),gp_Pnt(Node->x(),Node->y(),Node->z())));
          }

          GVertex *V2 = (*ite)->getEndVertex();
          int nbv2 = V2->getNumMeshVertices();
          std::cout<<"nbv2 : "<<  nbv2 <<std::endl;
          MVertex *Nodv2 = V2->mesh_vertices[0];
          std::cout<<" Num : "<<  Nodv2->getNum() <<std::endl;
          std::cout<<" xNv2 : "<<  Nodv2->x()<<" yNv2 : "<<  Nodv2->y()<<" zNv2 : "<<  Nodv2->z() <<std::endl;

          Array.SetValue(cc,gp_Pnt(Nodv2->x(),Nodv2->y(),Nodv2->z()));

          Pnts.push_back(gp_Pnt(Nodv2->x(),Nodv2->y(),Nodv2->z())); ///*****************
          //mapPnts.insert(std::pair<int, gp_Pnt>(Nodv2->getNum(),gp_Pnt(Nodv2->x(),Nodv2->y(),Nodv2->z())));

          nb_ptspararete.push_back(cc);///*****************

          Handle(Geom_BSplineCurve) SPL = GeomAPI_PointsToBSpline (Array,2,3,GeomAbs_C1,0).Curve();
          Handle(GeomAdaptor_HCurve) SPL1Adaptor =  new GeomAdaptor_HCurve(SPL);
          Handle(BRepFill_CurveConstraint) Cont= new BRepFill_CurveConstraint(SPL1Adaptor,0);
          BPSurf.Add(Cont);


      }

      nb_areteparface.push_back(edgs.size());///****************

  /*    std::vector<gp_Pnt > ptsface;

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

           MVertex *Nod = (*it)->mesh_vertices[i];
           gp_Pnt P(Nod->x(),Nod->y(),Nod->z());

           ptsface.push_back(P);

           }

       ptsfaces.push_back(ptsface);
  /*
           BPSurf.Perform();

           Standard_Integer MaxSeg = 8;
           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();

           std::cout<<" **************  num face : "<<  nbface <<std::endl;


          // Construction de Contour
          BRepBuilderAPI_MakeWire W1;

          for( std::list<GEdge*>::const_iterator ite = edgs.begin(); ite != edgs.end(); ite++){

          int nbe = (*ite)->getNumMeshVertices();

          TColgp_Array1OfPnt2d Array01(1,nbe+2);
          int cc = 1;


          GVertex *V1 = (*ite)->getBeginVertex();
          int nbv1 = V1->getNumMeshVertices();

          MVertex *Nodv1 = V1->mesh_vertices[0];


          Standard_Real U,V;
          GeomAPI_ProjectPointOnSurf Proj(gp_Pnt(Nodv1->x(),Nodv1->y(),Nodv1->z()),Surf1);
          Proj.LowerDistanceParameters(U, V);
          Array01.SetValue(cc,gp_Pnt2d(U,V));

          cc++;




          for (int ii = 0; ii < nbe; ii++){

          MVertex *Node = (*ite)->mesh_vertices[ii];
          GeomAPI_ProjectPointOnSurf Proj(gp_Pnt(Node->x(),Node->y(),Node->z()),Surf1);
          Proj.LowerDistanceParameters(U, V);
          Array01.SetValue(cc,gp_Pnt2d(U,V));

          cc++;

          }


          GVertex *V2 = (*ite)->getEndVertex();
          int nbv2 = V2->getNumMeshVertices();

          MVertex *Nodv2 = V2->mesh_vertices[0];

          GeomAPI_ProjectPointOnSurf Proj1(gp_Pnt(Nodv2->x(),Nodv2->y(),Nodv2->z()),Surf1);
          Proj1.LowerDistanceParameters(U, V);
          Array01.SetValue(cc,gp_Pnt2d(U,V));


          Handle(Geom2d_BSplineCurve) SSPL1 = Geom2dAPI_PointsToBSpline(Array01,2,3,GeomAbs_C1,0).Curve();
          TopoDS_Edge E = BRepBuilderAPI_MakeEdge(SSPL1,Surf1);
          W1.Add(E);



      }

  */

  //TopoDS_Face F1 = BRepBuilderAPI_MakeFace(Surf1  /*,W1,Standard_True*/);

  //B.Add(Sh,F1);
  /*

       }



       }


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

        FILE * pFile;
        pFile = fopen("cube.surf","wt");

        fprintf(pFile, "%d \n",Pnts.size());

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

        fprintf(pFile, "%d \n",nb_ptspararete.size());

        int cc_pt=0;

        for(int i=0; i<nb_ptspararete.size(); i++)
        {
            fprintf(pFile, "%d %d \n",i, nb_ptspararete[i]);
            for(int j=0; j<nb_ptspararete[i]; j++)
            {
                fprintf(pFile, "%d ",cc_pt);
                cc_pt++;
            }
            fprintf(pFile, "\n");

        }

        fprintf(pFile, "%d \n",nbface);
        int cc_arete=0;

        int nbf = nbface;

        for(int i=0; i<nbface; i++)
        {
            //std::vector<gp_Pnt > ptsface = ptsfaces[i];
            if(nb_areteparface[i] > 0)
            {
              fprintf(pFile, "%d %d \n",i, 0);


              fprintf(pFile, "%d \n",nb_areteparface[i]);
              for(int j=0; j<nb_areteparface[i]; j++)
              {
                  fprintf(pFile, "%d ",cc_arete);
                  cc_arete++;
              }
              fprintf(pFile, "\n");

            }
            else
            {
                nbf--;

            }

        }

        std::cout<<"nbf : "<<  nbf <<std::endl;


        fclose(pFile);


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



// stop gmsh
  GmshFinalize();





  return 0;
}