SurfaceMeshGenerator.cpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stéphane Del Pino

//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2, or (at your option)
//  any later version.

//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

//  $Id: SurfaceMeshGenerator.cpp,v 1.82 2007/06/09 10:37:06 delpinux Exp $

#include <Information.hpp>
#include <SurfaceMeshGenerator.hpp>

#include <MeshOfTetrahedra.hpp>
#include <MeshTetrahedrizor.hpp>

#include <SurfaceMeshOfTriangles.hpp>
#include <Structured3DMesh.hpp>

#include <Scene.hpp>

#include <Union.hpp>
#include <Intersection.hpp>
#include <Difference.hpp>
#include <Not.hpp>

#include <Cube.hpp>
#include <Plane.hpp>
#include <Cylinder.hpp>
#include <InfiniteCylinder.hpp>
#include <Cone.hpp>
#include <InfiniteCone.hpp>
#include <ObjectTransformer.hpp>

#include <Mesh.hpp>
#include <ConnectivityBuilder.hpp>
#include <Connectivity.hpp>
#include <Cell.hpp>

#include <WorkingMesh.hpp>

#include <TinyMatrix.hpp>

#include <triangulation.hpp>
#include <set>
#include <iostream>
#include <fstream>
#include <stack>
#include <algorithm>

inline double SDet(const Vertex & V1, const Vertex & V2, const Vertex & V3) {
    return (V2[0]-V1[0]) * (V3[1]-V1[1]) -  (V2[1]-V1[1]) * (V3[0]-V1[0]);
}

inline bool checkInterbis(const Vertex &P1, const Vertex &P2, const Vertex &Q1, const Vertex &Q2) {
    // retourne true si P1P2 inter Q1Q2 != vide
    ASSERT(!(P1 == P2 || P1 == Q1 || P1 == Q2 || P2 == Q1 || P2 == Q2 || Q1 == Q2));
    //const double eps = 1e-15;
    /*if (std::abs(SDet(Q1,Q2,P1)) <eps||
        std::abs(SDet(Q1,Q2,P2)) <eps) {
        //ffout(0) << " !!!!! NO !!!!!!!" << std::endl;
        return true;
    }*/

    return ((SDet(P1,P2,Q1)*SDet(P1,P2,Q2) < 0) && (SDet(Q1,Q2,P1)*SDet(Q1,Q2,P2) < 0)); // version optimisée
}
class SurfaceMeshGenerator::Internals
{
private:
  class IntersectionPoints
  {
  private:
    std::vector<const Vertex* > __listOfPoints;
    TinyVector<3, size_t> __cutTriangle1;
    TinyVector<3, size_t> __cutTriangle2;

  public:

    TinyVector<3, size_t>& cutTriangle1()
    {
      return __cutTriangle1;
    }

    std::vector<const Vertex*>& listOfPoints()
    {
      return __listOfPoints;
    }

    TinyVector<3, size_t>& cutTriangle2()
    {
      return __cutTriangle2;
    }  

    const TinyVector<3, size_t>& cutTriangle1() const
    {
      return __cutTriangle1;
    }

    const TinyVector<3, size_t>& cutTriangle2() const
    {
      return __cutTriangle2;
    }

    inline size_t size() const
    {
      return __listOfPoints.size();
    }

    const Vertex*& operator[](const size_t& i)
    {
      ASSERT(i<__listOfPoints.size());
      return __listOfPoints[i];
    }

    void addPoint(const Vertex*&  V)
    {
      __listOfPoints.push_back( V );
    }
    const IntersectionPoints& operator=(const IntersectionPoints& P)
    {
      __listOfPoints=P.__listOfPoints;
      __cutTriangle1=P.__cutTriangle1;
      __cutTriangle2=P.__cutTriangle2;
      return *this;
    }

    //constructor
    IntersectionPoints()
      : __cutTriangle1(0),
        __cutTriangle2(0)
    {
      ;
    }


    IntersectionPoints(const IntersectionPoints& P)
      : __listOfPoints(P.__listOfPoints),
        __cutTriangle1(P.__cutTriangle1),
        __cutTriangle2(P.__cutTriangle2)

    {
      ;
    }

    ~IntersectionPoints()
    {
      ;
    }
  };

  struct IntersectionTest
  {
    static bool compute(const Cell& C)
    {
      for (size_t i=0; i<C.numberOfVertices(); ++i) {
        if(C(i).reference()==0)
          return false;
      }
      return true;
    }
  };

  struct UnionTest
  {
    static bool compute(const Cell& C)
    {
      for (size_t i=0; i<C.numberOfVertices(); ++i) {
        if(C(i).reference()==2) {
          return false;
        }
      }
      return true;
    }
  };

  class TriangleCut;
  class MeshedObject;
  friend class SurfaceMeshGenerator;

  typedef std::list<TinyVector<3,Edge::Pair> > EdgePairList;      
  typedef std::list<const Cell*> MotherCellList;

  typedef std::map<TinyVector<3>, Object*> ObjectReferences;

  typedef std::map<Edge::Pair, size_t> EdgesRef;

  typedef std::map<Edge::Pair, Vertex*> VerticesList;

  typedef std::vector<ReferenceCounting<MeshedObject > > ListOfMeshedObject;
  
  typedef std::map<const Cell*, std::list<Triangle*> > MapCellTriangle;
 
  typedef std::pair<const Triangle*,const Triangle*> PairTriangleNew;
  //std::vector<Vertex*> listVertexIntersectionNew;

  typedef std::map<std::pair<const Vertex*,const Vertex*> ,
    TinyVector<2,const Vertex* > > PIntersection;

  typedef std::map< PairTriangleNew , IntersectionPoints> PairTriangleToIntersectionPoints;

  std::map<const Vertex*, const Vertex*>   listOfVertexMesh;
  std::map<const Triangle*, const Triangle*>   listOfTriangleMeshFront;

  PairTriangleToIntersectionPoints pairTriangleToIntersectionPoints;

  std::vector<ReferenceCounting< std::map<const Vertex*,size_t> > > edgesRefVertex;

  std::map<const Vertex*,TinyVector<2,const Triangle*> > vertexInTriangles;
  std::map<const Vertex*,std::vector<const Triangle*> > vertexVectTriangles; 
  std::map<const Triangle*,TinyVector<2,const Vertex*> > triangleWithVertex;

  TinyVector<3> __splitEdge(const Edge& e,
                            const Shape& S);

  void __dataStructureConvertion(MeshedObject& O,
                                 EdgePairList& triangleListes,
                                 VerticesList& verticesListes,
                                 MotherCellList& cellListObject);

  void __constructionVerticesList(const MeshedObject& O0);

  void __constructionFinalMesh(const MeshedObject& O0,
                               SurfaceMeshOfTriangles& s_mesh);

  void __calculatePointsIntersection(const MeshedObject& O1,
                                     const MeshedObject& O2,
                                     const std::set<const Cell*>& toTreatHexahedra,
                                     PIntersection& listVertexIntersectionNew);

  void __InTriangle(const Vertex& P,
                    TinyVector<3,const Vertex *> & V,
                    size_t& ncompt);

  void __inout(const size_t& numobject,
               const int& ttn0,
               const Vertex*& V0,
               const Vertex*& V1);

  TinyVector<4, real_t>  __EquationPlan(const TinyVector<3,const Vertex *> & V);

  void __DeterminatePoint(const size_t& ncase,
                          const Vertex* & Vo0,
                          const Vertex* & Vo1,
                          TinyVector<3, const Vertex* > & V,
                          IntersectionPoints& Pinter,
                          PIntersection& listVertexIntersectionNew);

  void __CalculPoint(std::list<Triangle*>::iterator & currentTriangle,
                     std::list<Triangle*>::iterator & currentTriangle1,
                     IntersectionPoints & Pinter,
                     PIntersection& listVertexIntersectionNew);

  void __createPointIntersection(std::list<Triangle*>::iterator currentTriangle1,
                                 std::list<Triangle*>::iterator currentTriangle2,
                                 PIntersection& listVertexIntersectionNew);
  
  const Vertex& __determinateNumber(const int& p,
                                    const TinyVector<3,TinyVector<2,const Vertex*> >& pointsint);

  void __determinateCase(const int& ncase,
                         const TriangleCut & K,
                         TinyVector<3,int>& tt,
                         TinyVector<3,int>& pp);

  void __addTriangle(const Triangle& T,
                     const Cell* cell,
                     std::vector<Triangle>& triangleListIntersectionNew);

  void __createLocalListIntersection(const size_t& objectNumber,
                                     const MeshedObject& O2,
                                     const Triangle*& currentTriangle,
                                     TriangleCut & K,
                                     const Cell* cell);
  
  TinyVector<3,TinyVector<2,int> > __determinateEdgeCut(const size_t& numobject,
                                                        const MeshedObject& O1,
                                                        const MeshedObject& O2,
                                                        const Triangle*& currentTriangle,
                                                        const Cell* cell);

  void __putEdgeRef(const Triangle*& Tmesh,
                    const Triangle*& T,
                    const size_t numobject,
                    const MeshedObject& O1,
                    const MeshedObject& O2);

  void __putRefByFront(const size_t numobject,
                       const MeshedObject& O1,
                       const MeshedObject& O2,
                       const SurfaceMeshOfTriangles& surfmesh,
                       const std::set<const Cell*>& toTreatHexahedra);

  void __createCase(const TriangleCut& K ,
                    size_t& in, size_t& in1,
                    TinyVector<3,size_t>& place);

  void __transformVertex(const TriangleCut& K,
                         std::vector<const Vertex*>& points,
                         std::vector<Vertex> &pointstrian);
  void __addPoints(const Vertex*& Pcut,
                   size_t& stopcut,
                   std::vector<const Vertex*>& pointslist1,
                   std::vector<const Vertex*>& pointslist2);
  
  void __addList(const TriangleCut& K);
  void __createList(const TriangleCut& K,
                    std::vector<const Vertex*>& points,
                    size_t& in,
                    size_t& in1,
                    TinyVector<3,size_t> place);

  bool __verifExist(const Vertex* Ptemp,
                    std::vector<const Vertex*>& points);
  bool  __findOneEdge(const TriangleCut& K,
                      const size_t & number,
                      TinyVector<2,const Vertex*>& P);
  
  void __findInEdges(const TriangleCut& K,
                     TinyVector<2,const Vertex*>& P,
                     TinyVector<2,const Triangle*>& T,
                     size_t& iInter,
                     size_t& stop,
                     size_t& in,
                     size_t& in1,
                     TinyVector<3,size_t> place,
                     std::vector<const Vertex*>& points);
      
  void __findVertex(const TriangleCut& K,
                    TinyVector<2,const Vertex*>& P,
                    TinyVector<2,const Triangle*>& T,
                    size_t& iCut,
                    size_t& iInter,
                    size_t& stop,
                    size_t& in,
                    size_t& in1,
                    TinyVector<3,size_t> place,
                    std::vector<const Vertex*>& points);

  void __findTriangle(TinyVector<2,const Vertex*>& P,
                      TinyVector<2,const Triangle*>& T,
                      size_t& stop);
    
  void __createGeneral(const TriangleCut& K,const Cell* cell,
                       std::vector<Triangle>& triangleListIntersectionNew);
  
  void __create2SD(const TriangleCut& K,const Cell* cell,
                         std::vector<const Vertex*>& points,
                         std::vector<Vertex> &pointstrian,
                         const TinyVector<3,size_t>& place,
                         std::vector<Triangle>& triangleListIntersectionNew);
  
  void __create3in(const TriangleCut& K,const Cell* cell,
                         std::vector<const Vertex*>& points,
                         const size_t in,const size_t in1,
                         const TinyVector<3,size_t>& place,
                         std::vector<Triangle>& triangleListIntersectionNew);

  void __createTriangle(const TriangleCut& K,const Cell* cell,
                        std::vector<Triangle>& triangleListIntersectionNew,
                        std::vector<const Vertex*>& points,
                        std::vector<Vertex>& pointtrian);

  template <typename booleanTest>
  void __createTrianglesIntersection(const size_t numobject,
                                     const MeshedObject& O1,
                                     const MeshedObject& O2,
                                     const std::set<const Cell*>& toTreatHexahedra,
                                     std::vector<Triangle>& triangleListIntersectionNew,
                                     PIntersection& listVertexIntersectionNew);

  bool __createTriangleSurface(VerticesList& verticesListes,
                               EdgePairList& localTriangleList,
                               const Shape& S ,
                               std::vector<const Edge*>& cutEdges,
                               const Cell& currentCell);

  void __createSurface(MeshedObject& O,
                       EdgePairList& triangleListes,
                       EdgePairList& localTriangleListesEdges,
                       VerticesList& verticesListes,
                       MotherCellList& cellListObject);

  bool __isDegenerate(const Triangle& T);

  template <typename booleanTest>
  void __calculateIntersection(MeshedObject& O0,
                               const MeshedObject& O1,
                               const MeshedObject& O2,
                               const std::set<const Cell*>& toTreatHexahedra);

  template <typename booleanTest>
  void __operationBoolean(MeshedObject& O0,
                          const MeshedObject& O1,
                          const MeshedObject& O2);

  void __generateMesh(const Domain& omega, const Object& object, const size_t& level,
                      std::stack<ReferenceCounting<MeshedObject> >& objectStack);

  void __generateCoarseMesh(const Domain& omega, const Object& object);

  void __setMotherCells(SurfaceMeshOfTriangles& surfmesh);

  ReferenceCounting<Vector<Triangle> >
  __marchingTetrahedra(MeshedObject& O);
    
  static void __getIntersectionReferences(const Intersection& I,
                                          std::map<TinyVector<3>, Object*>& otherReferences);

  void plotini(size_t nobjectToTreat);

  ReferenceCounting<MeshTetrahedrizor::CellMapping>
  __motherCells;        
  ReferenceCounting<MeshOfTetrahedra>
  __backgroundMesh;     
  const Mesh*
  __lastBackgoundMesh;  
public:
  Internals()
    : __motherCells(0),
      __backgroundMesh(0),
      __lastBackgoundMesh(0)
  {
    ;
  }
};

class SurfaceMeshGenerator::Internals::MeshedObject
{
private:
  ReferenceCounting<Shape> __shape;

public :

  void setShape(ReferenceCounting<Shape> s)
  {
    __shape = s;
  }

  ReferenceCounting<Shape> shapeReference() const
  {
    ASSERT(__shape != 0);
    return __shape;
  }

  const Shape& shape() const
  {
    return *__shape;
  }

  bool inside(const Vertex& V) const
  {
    return __shape->inside(V);
  }

  ReferenceCounting<Vector<Triangle> > trianglelist;
  ReferenceCounting<Vector<Vertex*> > verticeslist;
  ReferenceCounting< MapCellTriangle > hexalist;

private:  
  // Copying objects is forbidden
  const MeshedObject& operator=(const MeshedObject& O); 
  MeshedObject(const MeshedObject& O);
public:
  //constructor
  MeshedObject()
    : __shape(0)
  {
    ;
  }

  ~MeshedObject()
  {
    ;
  }
};

class SurfaceMeshGenerator::Internals::TriangleCut
  : public Triangle
{
public :
  TinyVector<3,TinyVector<2,int> > edgecut;
  TinyVector<3,size_t> isInTriangle;
  TinyVector<3,TinyVector<2,const Vertex*> > pointsIntersection;
  std::map< const Vertex* ,const Vertex*> pointsin;
  size_t numobject;
  size_t otherobject;

  //constructor
  TriangleCut()
  {
    ;
  }

  TriangleCut(const TriangleCut& K) {
    edgecut=K.edgecut;
    isInTriangle=K.isInTriangle;
    pointsIntersection=K.pointsIntersection;
    pointsin=K.pointsin;
    numobject=K.numobject;
    otherobject=K.otherobject;
  }

  TriangleCut(const Vertex& v0,
              const Vertex& v1,
              const Vertex& v2,
              size_t num,
              size_t num2)
    : Triangle(v0,v1,v2)
  {
    numobject=num;
    otherobject=num2;
    edgecut=0;
    isInTriangle=0;
  }

  ~TriangleCut()
  {
    ;
  }
};


TinyVector<3, real_t>
SurfaceMeshGenerator::Internals::
__splitEdge(const Edge& e,
            const Shape& S)
{
  int in = -1;          // number of the vertex inside
  int out = -1;         // same thing for outside

  for (int j=0; j<2; j++) {
    if (e(j).reference() != 1) {
      out = j;
    } else {
      in = j;
    }
  }

  ASSERT((in != -1)&&(out != -1));

  TinyVector<3> vi = e(in);
  TinyVector<3> vo = e(out);
  TinyVector<3> vs = 0.5*(vi+vo);

  for (int l = 0; l<20/*20*/; l++) { // 3% error after 5 iterations.
    if(S.inside(vs)) {
      vi = vs;
      vs = 0.5*(vo+vi);
    } else {
      vo = vs;
      vs = 0.5*(vo+vi);
    }
  }

  return vs;
}

template <typename booleanTest>
void SurfaceMeshGenerator::Internals::
__operationBoolean(MeshedObject& O0,
                   const MeshedObject& O1,
                   const MeshedObject& O2)
{
  std::set<const Cell*> hexaList1;
  std::set<const Cell*> hexaList2;

  for (MapCellTriangle::const_iterator i=(*O1.hexalist).begin(); i!= (*O1.hexalist).end(); ++i) {
    hexaList1.insert((*i).first);
  }
  for (MapCellTriangle::const_iterator i=(*O2.hexalist).begin(); i!= (*O2.hexalist).end(); ++i) {
    hexaList2.insert((*i).first);
  }

  std::set<const Cell*> toTreatHexahedra;

  if (not(Information::instance().coarseMesh())) {
    std::set_intersection (hexaList1.begin(), hexaList1.end(),
                           hexaList2.begin(), hexaList2.end(),
                           inserter(toTreatHexahedra, toTreatHexahedra.begin()));
    ffout(4)<<"nb hexa to treat "<<toTreatHexahedra.size()<<"\n";
  }
  
  __calculateIntersection<booleanTest>(O0,O1,O2, toTreatHexahedra);
}

template <typename booleanTest>
void SurfaceMeshGenerator::Internals::
__calculateIntersection(MeshedObject& O0,
                        const MeshedObject& O1,
                        const MeshedObject& O2,
                        const std::set<const Cell*>& toTreatHexahedra)
{
  PIntersection listVertexIntersectionNew;
  std::vector<Triangle> triangleListIntersectionNew;

  // CALCUL DES INTERSECTIONS
  __calculatePointsIntersection(O1, O2, toTreatHexahedra, listVertexIntersectionNew);

  //create mesh intersection for the object 0 and 1
  // GENERATION DES MAILLAGES
  __createTrianglesIntersection<booleanTest>(0, O1, O2, toTreatHexahedra,
                                             triangleListIntersectionNew, listVertexIntersectionNew);
  ffout(4)<<"nb triangles in intersection 0 "<<triangleListIntersectionNew.size()<<'\n';
  __createTrianglesIntersection<booleanTest>(1, O2, O1, toTreatHexahedra,
                                             triangleListIntersectionNew, listVertexIntersectionNew);
  ffout(4)<<"nb triangles in intersection 1 "<<triangleListIntersectionNew.size()<<'\n';

  pairTriangleToIntersectionPoints.clear();
  
  //triangle from intersection of objects 0-1 in O0
  size_t sizeinter=triangleListIntersectionNew.size();
  
  O0.trianglelist = new Vector<Triangle> (sizeinter);
  O0.hexalist = new std::map<const Cell*, std::list<Triangle*> >;
  for(size_t size=0 ; size<sizeinter ; ++size) {
    (*O0.trianglelist)[size]=triangleListIntersectionNew[size];
    (*O0.hexalist)[&(*O0.trianglelist)[size].mother()].push_back(&(*O0.trianglelist)[size]);
  }
  triangleListIntersectionNew.clear();
  listVertexIntersectionNew.clear();
}


void SurfaceMeshGenerator::Internals::__constructionVerticesList(const MeshedObject& O0) {

  const Vector<Triangle>& triangleList=(*O0.trianglelist);

  for(size_t size=0 ; size< triangleList.size() ; ++size) {
    const Triangle& T=triangleList[size];
    if (listOfVertexMesh.find(&T(0))==listOfVertexMesh.end()) {
      listOfVertexMesh[&T(0)] = new Vertex(T(0));;
    }
    if (listOfVertexMesh.find(&T(1))==listOfVertexMesh.end()) {
      listOfVertexMesh[&T(1)] = new Vertex(T(1));;
    }
    if (listOfVertexMesh.find(&T(2))==listOfVertexMesh.end()) {
      listOfVertexMesh[&T(2)] = new Vertex(T(2));;
    }
  }
}

void SurfaceMeshGenerator::Internals::__constructionFinalMesh(const MeshedObject& O0,
                                                              SurfaceMeshOfTriangles& s_mesh) {

  __constructionVerticesList(O0);

  size_t ntr=0;
  
  ntr+=listOfVertexMesh.size();
  
  s_mesh.setNumberOfVertices(ntr);
  ffout(4)<<"nb tot vertices "<<s_mesh.numberOfVertices()<<'\n';
  int i = 0;
  /* for(size_t v=0 ; v<ntr ; ++v) {
     s_mesh.vertex(i) = *listVertexFinal[v];
     listVertexFinal[v] = &s_mesh.vertex(i);
     i++;
     }*/

  std::set<const Vertex*> newVertices;
  std::set<const Vertex*> oldVertices;

  for(std::map<const Vertex*, const Vertex*>::iterator it=listOfVertexMesh.begin() ;
      it!=listOfVertexMesh.end() ; it++) {
    s_mesh.vertex(i) = *((*it).second);

    delete ((*it).second);
    ((*it).second) = &s_mesh.vertex(i);
    i++;
  }

  const Vector<Triangle>& triangleList=(*O0.trianglelist);
  i = 0;
  size_t ncell=0;
  //if(nobjectToTreat==0){
  ncell+=triangleList.size();
  //}
  //  ffout(4)<<"ncell "<<ncell<<" n0 "<<n0<<" n1 "<<n1<<"\n";
  
  s_mesh.setNumberOfCells(ncell);
  ffout(4)<<"nb tot cells "<<s_mesh.numberOfCells()<<'\n';
  
  for(size_t size=0 ; size< triangleList.size() ; ++size) {
    const Triangle& T=triangleList[size];
    std::map<const Vertex*, const Vertex*>::iterator
      it0=listOfVertexMesh.find(&T(0));
    const Vertex& V0=*((*it0).second);
    std::map<const Vertex*, const Vertex*>::iterator
      it1=listOfVertexMesh.find(&T(1));
    const Vertex& V1=*((*it1).second);
    std::map<const Vertex*, const Vertex*>::iterator
      it2=listOfVertexMesh.find(&T(2));
    const Vertex& V2=*((*it2).second);
    s_mesh.cell(i) = Triangle(V0,V1,V2, T.reference());
    s_mesh.cell(i).setMother(&T.mother(),
                             std::numeric_limits<size_t>::max()); // Le triangle n'est pas une face de la maille
    listOfTriangleMeshFront[&s_mesh.cell(i)]=&T;
    i++;
    
  }
  // }
  //listOfVertexMesh.clear();
}

void SurfaceMeshGenerator::Internals::
__dataStructureConvertion(MeshedObject& O,
                          EdgePairList& triangleListes,
                          VerticesList& verticesListes,
                          MotherCellList& cellListObject)
{
  //for (size_t numObj = 0; numObj<nb; ++numObj) {
  // Data structure convertion
  {
#warning keep this ?
    size_t n=0;
    O.verticeslist
      = new Vector<Vertex*>(verticesListes.size());
    Vector<Vertex*>& vertexListbis = (*O.verticeslist);
    for (SurfaceMeshGenerator::Internals::VerticesList::iterator
           i = verticesListes.begin();
         i != verticesListes.end(); ++i, ++n) {
      vertexListbis[n] = (*i).second;
    }
  }
  
  {
    size_t n=0;
    O.hexalist
      = new std::map<const Cell*, std::list<Triangle*> >;
    std::map<const Cell*, std::list<Triangle*> >& hexaToTrianglebis=(*O.hexalist);
    
    O.trianglelist
      = new Vector<Triangle> (triangleListes.size());
    Vector<Triangle>& triangleListbis = (*O.trianglelist);
    SurfaceMeshGenerator::Internals::MotherCellList::iterator
      itcell=cellListObject.begin();
    
    for (SurfaceMeshGenerator::Internals::EdgePairList::iterator
           i = triangleListes.begin();
         i != triangleListes.end(); ++i, ++n) {
      SurfaceMeshGenerator::Internals::VerticesList& V = verticesListes;
      Triangle T(* V[(*i)[0]],
                 * V[(*i)[1]],
                 * V[(*i)[2]], 1);
      
      T.setMother(*itcell,
                  std::numeric_limits<size_t>::max()); // Le triangle n'est pas une face de la maille
      triangleListbis[n] = T;
      hexaToTrianglebis[*itcell].push_back(&triangleListbis[n]);
      ++itcell;
    }
  }
  
  verticesListes.clear();
  triangleListes.clear();
  //}
}

void SurfaceMeshGenerator::Internals::
__calculatePointsIntersection(const MeshedObject& O1,
                              const MeshedObject& O2,
                              const std::set<const Cell*>& toTreatHexahedra,
                              PIntersection& listVertexIntersectionNew)
{
  // calculate points intersection between object O1 and O2
  const MapCellTriangle& hexaToTriangle0=*(O1.hexalist);
  const MapCellTriangle& hexaToTriangle1=*(O2.hexalist);
    
  for(std::set<const Cell*>::const_iterator i=toTreatHexahedra.begin();
      i!=toTreatHexahedra.end(); ++i) {
    MapCellTriangle::const_iterator it0=hexaToTriangle0.find(*i);
    MapCellTriangle::const_iterator it1=hexaToTriangle1.find(*i);

    ASSERT(it0!=hexaToTriangle0.end() and it1!=hexaToTriangle1.end());

    std::list<Triangle*> listTriangle0=(*it0).second;
    std::list<Triangle*> listTriangle1=(*it1).second;

    for(std::list<Triangle*>::iterator jt0=listTriangle0.begin() ;
        jt0!=listTriangle0.end() ; ++jt0) {
      for(std::list<Triangle*>::iterator jt1=listTriangle1.begin() ;
          jt1!=listTriangle1.end() ; ++jt1) {
        __createPointIntersection(jt0,jt1,listVertexIntersectionNew);
      }
    }
  }
}

bool SurfaceMeshGenerator::Internals::__isDegenerate(const Triangle& T)
{
  return (&(T(0))==&(T(1)) or &(T(0))==&(T(2)) or &(T(1))==&(T(2)));
}

void SurfaceMeshGenerator::Internals::
__InTriangle(const Vertex& P,
             TinyVector<3,const Vertex* > & V,
             size_t& compt)
{
  compt=0;
  TinyVector<3,TinyVector<3,real_t> >  triangleVector;
  triangleVector[0]= -(*V[2]) + (*V[0]);
  triangleVector[1]= -(*V[2]) + (*V[1]);
  triangleVector[2]= P - (*V[2]);

  TinyMatrix<3,3, real_t> A;
  for (unsigned i=0; i<3; ++i) {
    for (unsigned j=0; j<3; ++j) {
      A(i,j)=(*V[j])[i];
    }
  }

  TinyVector<3, real_t> lambda;

  gaussPivot(A,P,lambda);

  bool ok=true;
  for (size_t i=0; i<3; ++i) {
    ok = (lambda[i]>-1e-7)?ok:false;
  }

  compt = (ok) ? 3 : 0;
}

TinyVector<4, real_t>  SurfaceMeshGenerator::Internals::
__EquationPlan(const TinyVector<3,const Vertex *> & V){

  //equation of plan
  TinyVector<3> V01=(*V[1])-(*V[0]);
  TinyVector<3> V12=(*V[2])-(*V[1]);
  TinyVector<3> V02=(*V[2])-(*V[0]);

  TinyVector<3> coeffo;
  coeffo=V01^V02;
  real_t d2 = -(coeffo*(*V[0]));

  TinyVector<4> coeff;
  coeff[0]=coeffo[0];
  coeff[1]=coeffo[1];
  coeff[2]=coeffo[2];
  coeff[3]=d2;

  return coeff;
}

void SurfaceMeshGenerator::Internals::
__DeterminatePoint(const size_t& ncase,
                   const Vertex* & Vo0,
                   const Vertex* & Vo1,
                   TinyVector<3,const Vertex* > & V,
                   IntersectionPoints& Pinter,
                   PIntersection& listVertexIntersectionNew)
{
  //calculate intersection between edge Vo0Vo1 and triangle V
  // if ncase<3 the edge is in the first object else the edge is in the second object
  // P = point intersection
  std::pair<const Vertex*,const Vertex*> pV(Vo0,Vo1);
  std::pair<const Vertex*,const Vertex*> pV2(Vo1,Vo0);

    // rq: epsilon a 10-6 on loupe des points!!!!! (testStep2.txt)
  real_t epsilon=1e-10;
  real_t eps=1e-6;
  real_t pscalar,k;
  size_t num;
  Vertex P;
  TinyVector<4> coeffo=__EquationPlan(V);
  TinyVector<3> coeff;
  coeff[0]=coeffo[0];
  coeff[1]=coeffo[1];
  coeff[2]=coeffo[2];
  real_t d=coeffo[3];

  TinyVector<3> Vo01 = (*Vo1) - (*Vo0);
  TinyVector<3> V01 = (*V[1]) - (*V[0]);
  TinyVector<3> V12 = (*V[2]) - (*V[1]);
  TinyVector<3> V02 = (*V[2]) - (*V[0]);
    
  real_t norm1=Norm(Vo01);
  pscalar=coeff*Vo01;
  num=0;
  if(std::abs(pscalar)>1e-10) {
    k=-(coeff*(*Vo0)+d)/pscalar;
    P=(*Vo0)+k*Vo01;
        
    TinyVector<3> V0P=P-(*V[0]);
    TinyVector<3> V1P=P-(*V[1]);
    TinyVector<3> V2P=P-(*V[2]);

    real_t norm2=Norm(P-(*Vo0));
    real_t norm3=Norm(P-(*Vo1));
    real_t inV01=Norm(V01^V0P);
    real_t inV12=Norm(V12^V1P);
    real_t inV02=Norm(V02^V2P);

    if(/*norm2>epsilon and norm3>epsilon and */std::abs((norm2+norm3)-norm1) < epsilon ) {
      if(ncase<3) {
        __InTriangle(P,V,num);
        if(num==3) {
          if(!(norm2>epsilon)) {
            //point confondus avec un des sommets
            Pinter.addPoint(Vo0);
          } else
            if(!(norm3>epsilon)) {
              //point confondus avec l'autre sommet
              Pinter.addPoint(Vo1);
            } else {
              // CONSIDÈRE QU'IL PEUT Y AVOIR DEUX POINTS D'INTERSECTION SUR L'ARÊTE, MAIS PAS PLUS
              if(listVertexIntersectionNew.find(pV)==listVertexIntersectionNew.end()
                 and listVertexIntersectionNew.find(pV2)==listVertexIntersectionNew.end()) {
                listVertexIntersectionNew[pV][0]=new Vertex (P);
                listVertexIntersectionNew[pV][1]=new Vertex (P);
                Pinter.addPoint(listVertexIntersectionNew[pV][0]);
              } else {
                if(listVertexIntersectionNew.find(pV)!=listVertexIntersectionNew.end()) {
                  if(Norm((*(*listVertexIntersectionNew.find(pV)).second[0])
                          -P)<eps) {
                    Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[0]);
                  } else
                    if(Norm((*(*listVertexIntersectionNew.find(pV)).second[1]) - P) < eps) {
                      Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[1]);
                    } else {
                      listVertexIntersectionNew[pV][1]=new Vertex (P);
                      Pinter.addPoint(listVertexIntersectionNew[pV][1]);
                    }
                } else {
                  if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[0]) - P) < eps) {
                    Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[0]);
                  } else
                    if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[1]) - P) < eps) {
                      Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[1]);
                    } else {
                      listVertexIntersectionNew[pV2][1]=new Vertex (P);
                      Pinter.addPoint(listVertexIntersectionNew[pV2][1]);
                    }
                }
              }
            }
          Pinter.cutTriangle1()[ncase]=Pinter.size();
        }
      } else {
        //il ne faut pas que le point soit sur une arete du tr 1
        if(inV01 > epsilon and inV02 > epsilon and inV12 > epsilon ) {
          __InTriangle(P,V,num);
          if(num==3) {
            //ffout(4)<<"ajout du point cas 2\n";
            //ffout(4)<<P[0]<<" "<<P[1]<<" "<<P[2]<<" \n";
            //listVertexIntersectionNew.push_back(new Vertex (P));
            PIntersection PIn;
            if(!(norm2>epsilon)) {
              Pinter.addPoint(Vo0);
            } else {
              if(!(norm3>epsilon)) {
                Pinter.addPoint(Vo1);
              } else 
                {
                  if(listVertexIntersectionNew.find(pV)==listVertexIntersectionNew.end()
                     and listVertexIntersectionNew.find(pV2)==listVertexIntersectionNew.end()) {
                    listVertexIntersectionNew[pV][0]=new Vertex (P);
                    listVertexIntersectionNew[pV][1]=new Vertex (P);
                    Pinter.addPoint(listVertexIntersectionNew[pV][0]);
                  } else {
                    if(listVertexIntersectionNew.find(pV)!=listVertexIntersectionNew.end()) {
                      if(Norm((*(*listVertexIntersectionNew.find(pV)).second[0])
                              -P)<eps) {
                        Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[0]);
                      } else {
                        if(Norm((*(*listVertexIntersectionNew.find(pV)).second[1])
                                -P)<eps) {
                          Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[1]);
                        } else {
                          listVertexIntersectionNew[pV][1]=new Vertex (P);
                          Pinter.addPoint(listVertexIntersectionNew[pV][1]);
                        }
                      }
                    } else {
                      if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[0])
                              -P)<eps) {
                        Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[0]);
                      } else {
                        if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[1])
                                -P)<eps) {
                          Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[1]);
                        } else {
                          listVertexIntersectionNew[pV2][1]=new Vertex (P);
                          Pinter.addPoint(listVertexIntersectionNew[pV2][1]);
                        }
                      }
                    }
                  }
                }
            }
            // Pinter.addPoint(listVertexIntersectionNew[listVertexIntersectionNew.size()-1]);
          }
        }
      }
    }
  }
}

void SurfaceMeshGenerator::Internals::
__CalculPoint(std::list<Triangle*>::iterator & currentTriangle,
              std::list<Triangle*>::iterator & currentTriangle1,
              IntersectionPoints& Pinter,
              PIntersection& listVertexIntersectionNew)
{
  // calculate points intersection between triangle "currentTriangle" and "currentTriangle1"
  // this points are in "Pinter"

  // NECESSITE DE V012 et Vo012 ?

  TinyVector<3,const Vertex * > V012;
  const Triangle& T0=*(*currentTriangle1);
  V012[0] = &T0(0);
  V012[1] = &T0(1);
  V012[2] = &T0(2);
  TinyVector<3,const Vertex* > Vo012;
  Triangle& To0=*(*currentTriangle);
  Vo012[0] = &To0(0);
  Vo012[1] = &To0(1);
  Vo012[2] = &To0(2);

  __DeterminatePoint(0,(V012[0]),(V012[1]),Vo012,Pinter, listVertexIntersectionNew);
  __DeterminatePoint(1,(V012[0]),(V012[2]),Vo012,Pinter, listVertexIntersectionNew);
  __DeterminatePoint(2,(V012[1]),(V012[2]),Vo012,Pinter, listVertexIntersectionNew);
  
  __DeterminatePoint(3,(Vo012[0]),(Vo012[1]),V012,Pinter, listVertexIntersectionNew);
  __DeterminatePoint(3,(Vo012[0]),(Vo012[2]),V012,Pinter, listVertexIntersectionNew);
  __DeterminatePoint(3,(Vo012[1]),(Vo012[2]),V012,Pinter, listVertexIntersectionNew);
}

void SurfaceMeshGenerator::Internals::
__inout(const size_t& numobject,
        const int & edgeCut,
        const Vertex*& V0,
        const Vertex*& V1)
{// if the edge is not cut , put ref=2 for the 2 vertex
  size_t& refVertex0=(*edgesRefVertex[numobject])[V0];
  size_t& refVertex1=(*edgesRefVertex[numobject])[V1];

  if(refVertex0+refVertex1==4 and edgeCut==1) {
    ffout(4)<<"warning edge ref 2 cut \n";
    ffout(4)<<(*V0)[0]<<" "<<(*V0)[1]<<" "<<(*V0)[2]<<" \n";
    ffout(4)<<(*V1)[0]<<" "<<(*V1)[1]<<" "<<(*V1)[2]<<" \n";
    refVertex0=0;
    refVertex1=0;
  } else {
    // CECILE: ne serait-ce pas plutôt une ASSERTION ?
    if(refVertex0+refVertex1==2 and edgeCut==0) {
      refVertex0=2;
      refVertex1=2;
    }
  }
}

void SurfaceMeshGenerator::Internals::
__createPointIntersection(std::list<Triangle*>::iterator currentTriangle1,
                          std::list<Triangle*>::iterator currentTriangle,
                          PIntersection& listVertexIntersectionNew)
{    
  //ffout(4)<< "createpointintersection "<<numobject<<" "<<objectNumber<<"\n";
  //ffout(4)<<"*************************\n";

  PairTriangleNew pTriangle;
  pTriangle.first=(*currentTriangle1);
  pTriangle.second=(*currentTriangle);

  IntersectionPoints Pinter;
  // CALCUL DU POINT D'INTERSECTION
  __CalculPoint(currentTriangle, currentTriangle1, Pinter, listVertexIntersectionNew);

  //const TinyVector<3,size_t >& pintercut=Pinter.cutTriangle1();

  if(Pinter.size()!=0) {
    //ffout(4)<<"rajout2 \n";
    pairTriangleToIntersectionPoints[pTriangle] = Pinter;
  }

  if(pairTriangleToIntersectionPoints.find(pTriangle)!= pairTriangleToIntersectionPoints.end()) {
    TinyVector<3,size_t >& cutTr2= pairTriangleToIntersectionPoints[pTriangle].cutTriangle2();

    for(size_t i=0 ; i<Pinter.size() ; ++i) {
      real_t epsilon=1e-7;
      const Triangle& To=(*(*currentTriangle));
      const Vertex& Vo0 = To(0);
      const Vertex& Vo1 = To(1);
      const Vertex& Vo2 = To(2);
      const Vertex& P=*(pairTriangleToIntersectionPoints[pTriangle][i]);

      // TESTE SI LE POINT EST SUR UNE DES ARETES EN TESTANT TROIS "EGALITÉ" TRIANGULAIRE
            
      real_t Vo01=Norm(Vo1-Vo0);
      real_t Vo12=Norm(Vo2-Vo1);
      real_t Vo02=Norm(Vo2-Vo0);
      real_t norm0=Norm(P-Vo0);
      real_t norm1=Norm(P-Vo1);
      real_t norm2=Norm(P-Vo2);
      //for filter points
      //real_t eps=1e-6;
      if(/*norm0>eps and norm1>eps and */std::abs((norm0+norm1)-Vo01)<epsilon) {
        cutTr2[0]=i+1;
      }
      if(/*norm0>eps and norm2>eps and */std::abs((norm0+norm2)-Vo02)<epsilon) {
        cutTr2[1]=i+1;
      }
      if(/*norm2>eps and norm1>eps and */std::abs((norm2+norm1)-Vo12)<epsilon) {
        cutTr2[2]=i+1;
      }
    }
  }
}       


const Vertex& SurfaceMeshGenerator::Internals::
__determinateNumber(const int& place,
                    const TinyVector<3,TinyVector<2,const Vertex*> >& pointsint)
{//return the point intersection 
  if(place<=2) {
    return *pointsint[place][0];
  } else {
    return *pointsint[2][0];
  }
}

void SurfaceMeshGenerator::Internals::
__addTriangle(const Triangle& t,
              const Cell* cell,
              std::vector<Triangle>& triangleListIntersectionNew)
{
  Triangle T(t);
  T.setMother(cell,
              std::numeric_limits<size_t>::max()); // Le triangle n'est pas une face de la cellule
  if(!__isDegenerate(T)) {
    triangleListIntersectionNew.push_back(T);
  } else {
    fferr(3) << "\t\tTriangle dégénéré :"
             << &T(0) << ':' << &T(1) << ':' << &T(2) << '\n';
    fferr(3)
      << T(0) << ':' << T(1) << ':' << T(2) << '\n';
    // triangleListIntersectionNew.push_back(T);
  }
}

void SurfaceMeshGenerator::Internals::
__determinateCase(const int& ncase,
                  const TriangleCut & K,
                  TinyVector<3,int>& edgesCut,
                  TinyVector<3,int>& vertexNum)
{
  //determinate which point is inside
  const TinyVector<3,TinyVector<2,int> >& edgesCut2=K.edgecut;
  switch(ncase) {    
  case 0: {
    edgesCut[0]=edgesCut2[0][0];
    edgesCut[1]=edgesCut2[1][0];
    edgesCut[2]=edgesCut2[2][0];
    vertexNum[0]=0;
    vertexNum[1]=1;
    vertexNum[2]=2;
    break;
  }
  case 1: {
    edgesCut[0]=edgesCut2[0][0];
    edgesCut[1]=edgesCut2[2][0];
    edgesCut[2]=edgesCut2[1][0];
    vertexNum[0]=0;
    vertexNum[1]=2;
    vertexNum[2]=1;
    break;
  }
  case 2: {
    edgesCut[0]=edgesCut2[1][0];
    edgesCut[1]=edgesCut2[2][0];
    edgesCut[2]=edgesCut2[0][0];
    vertexNum[0]=1;
    vertexNum[1]=2;
    vertexNum[2]=0;
    break;
  }
  }
}

void SurfaceMeshGenerator::Internals::
__createLocalListIntersection(const size_t& objectNumber,
                              const MeshedObject& O2,
                              const Triangle*& currentTriangle,
                              TriangleCut & K,
                              const Cell* cell)
{
  //build the vectors K.pointsIntersection, K.pointsin with points intersection include in K
  //put K.edgecut[i]=1 if the edge i is cut
  //build vertexInTriangles, triangleWithVertex (for build mesh by front)
  //build vertexVectTriangles 
  TinyVector<3,TinyVector<2,int> >& edgesCut=K.edgecut;
  TinyVector<3,TinyVector<2,const Vertex*> > & pointsIntersection=K.pointsIntersection;
  std::map<const Vertex* ,const Vertex*>& pointsInterior=K.pointsin;
  
  int numberPointsInterior=pointsInterior.size();
  const size_t& numobject=K.numobject;
  const MapCellTriangle& hexaToTriangle1=*(O2.hexalist);
        
  MapCellTriangle::const_iterator it1=hexaToTriangle1.find(cell);
  if(it1!=hexaToTriangle1.end()) {
    std::list<Triangle*> listTriangle1=(*it1).second;
    for(std::list<Triangle*>::iterator jt0=listTriangle1.begin() ;
        jt0!=listTriangle1.end() ; ++jt0) {
 
      PairTriangleNew pairTriangle;
      if(numobject<objectNumber) {
        pairTriangle.first=(currentTriangle);
        pairTriangle.second=(*jt0);
      } else {
        pairTriangle.second=(currentTriangle);
        pairTriangle.first=(*jt0);
      }
      size_t temp=0;
            
      if(pairTriangleToIntersectionPoints.find(pairTriangle)!=pairTriangleToIntersectionPoints.end()) {
        TinyVector<3,size_t > cutTriangle;
        if(numobject<objectNumber) {
          cutTriangle=(pairTriangleToIntersectionPoints[pairTriangle]).cutTriangle1();
        } else {
          cutTriangle=(pairTriangleToIntersectionPoints[pairTriangle]).cutTriangle2();
        }
        //detection tr plat
        size_t isDegenerate=0;
        if(__isDegenerate(*(*jt0))) {
          // ffout(4)<<"degenerate triangle\n";
          isDegenerate=1;
        }
        if(pairTriangleToIntersectionPoints[pairTriangle].size()==2 and isDegenerate==0) {
          const Vertex*& V0=pairTriangleToIntersectionPoints[pairTriangle][0];
          const Vertex*& V1=pairTriangleToIntersectionPoints[pairTriangle][1];
          const Triangle T=(*(*jt0));
          if((V0==&(*(*jt0))(0) or V0==&(*(*jt0))(1) or V0==&(*(*jt0))(2))
             and (V1==&(*(*jt0))(0) or V1==&(*(*jt0))(1) or V1==&(*(*jt0))(2))) {
            //ffout(4)<<"inter confondu avec arete\n";
            isDegenerate=0;
          }
        }
        if(pairTriangleToIntersectionPoints[pairTriangle].size()>=2 and isDegenerate==0) {
          for(size_t b=0; b<pairTriangleToIntersectionPoints[pairTriangle].size() ; ++b) {
            const Vertex*& V=pairTriangleToIntersectionPoints[pairTriangle][b];
            const Triangle T=(*(*jt0));
            // si on a le point est confondu avec un sommet on n'a pas besoin du triangle pour la connectivite
            if(!(V==&(*(*jt0))(0) or V==&(*(*jt0))(1) or V==&(*(*jt0))(2)) and
        Norm(*V-(*(*jt0))(0))>1e-6 and Norm(*V-(*(*jt0))(1))>1e-6 and Norm(*V-(*(*jt0))(2))>1e-6 
               /*or b==0*/) {
              vertexVectTriangles[V].push_back(*jt0);
              if(vertexInTriangles.find(V)==vertexInTriangles.end()) {
                TinyVector<2,const Triangle*> triangles;
                triangles[0]=*jt0;
                triangles[1]=*jt0;
                vertexInTriangles[V]=triangles;
              } else {
                (*vertexInTriangles.find(V)).second[1]=*jt0;
              }
              if(triangleWithVertex.find(*jt0)==triangleWithVertex.end()) {
                TinyVector<2,const Vertex*> vertexTriangle;
                vertexTriangle[0]=V;
                vertexTriangle[1]=V;
                triangleWithVertex[*jt0]=vertexTriangle;
              } else {
                (*triangleWithVertex.find(*jt0)).second[1]=V;
              }
            } else {
              //le point est confondu avec un sommet du tr mais dans certains cas on en a besoin qd meme!!
              vertexVectTriangles[V].push_back(*jt0);
            }
          }
        } else {
            if(pairTriangleToIntersectionPoints[pairTriangle].size()!=0) {
                //ffout(4)<<"on a un point non pris... triangle plat??\n";
                //ffout(4)<<(*(*jt0))(0)<<"\n";
                //ffout(4)<<(*(*jt0))(1)<<"\n";
                //ffout(4)<<(*(*jt0))(2)<<"\n";
                //ffout(4)<<*pairTriangleToIntersectionPoints[pairTriangle][0]<<"\n";
                vertexVectTriangles[pairTriangleToIntersectionPoints[pairTriangle][0]].push_back(*jt0);
                
            } 
        }

        for (size_t n=0; n < cutTriangle.size(); ++n) {
          if(cutTriangle[n]!=0) {
            temp+=1;
            if(edgesCut[n][0]==0) {
              edgesCut[n][0]=1;
              pointsIntersection[n][0]=
                (pairTriangleToIntersectionPoints[pairTriangle][cutTriangle[n]-1]);
            } else {
              edgesCut[n][1]=1;
              pointsIntersection[n][1]=
                (pairTriangleToIntersectionPoints[pairTriangle][cutTriangle[n]-1]);
            }
          }
        }

        if(temp<pairTriangleToIntersectionPoints[pairTriangle].size()) {
          numberPointsInterior+=1;
          size_t r=0;
          size_t sum=cutTriangle[0]+cutTriangle[1]+cutTriangle[2];
          if(sum>5) {
            ffout(4)<<"sum "<<sum<<'\n';
          }

          switch (sum) {
          case 0: {
            r=0;
            for(size_t i=1 ; i<pairTriangleToIntersectionPoints[pairTriangle].size() ;++i)
              pointsInterior[(pairTriangleToIntersectionPoints[pairTriangle][i])]=
                (pairTriangleToIntersectionPoints[pairTriangle][i]);
            break;
          }
          case 2:
          case 5: {
            r=0;
            break;
          case 4:
          case 1:
            r=1;
            break;
          }
          case 3: {
            r=2;
            break;
          }
          default: {
            ffout(4) << "pbs sum = "<<sum<<" size "
                     << pairTriangleToIntersectionPoints[pairTriangle].size()<<" \n";
            ffout(4) << K(0) << "\n";
            ffout(4) << K(1) << "\n";
            ffout(4) << K(2) << "\n";
            ffout(4) << "cutTriangle "<<cutTriangle[0] << " " << cutTriangle[1] << " " << cutTriangle[2] << "\n";
            ffout(4) << "edgesCut "<<edgesCut[0] << " " << edgesCut[1] << " " << edgesCut[2] << "\n";
            //edgesCut[0]=cutTriangle[0];
            //edgesCut[1]=cutTriangle[1];
            //edgesCut[2]=cutTriangle[2];
            
          }
          }
          if(pointsInterior.find(pairTriangleToIntersectionPoints[pairTriangle][r])
             ==pointsInterior.end() and sum<=5) {
            pointsInterior[(pairTriangleToIntersectionPoints[pairTriangle][r])]
              =(pairTriangleToIntersectionPoints[pairTriangle][r]);
          }
        }
      } else {
        // ffout(4) <<"pasintersecte\n";
      }
    }
  }
}

TinyVector<3,TinyVector<2,int> > SurfaceMeshGenerator::Internals::
__determinateEdgeCut(const size_t& numobject,
                     const MeshedObject& O1,
                     const MeshedObject& O2,
                     const Triangle*& currentTriangle,
                     const Cell* cell)
{
  // put edgecut=1 if the edge is cut  
  TinyVector<3,TinyVector<2,int> > edgesCut;
  edgesCut=0;
  size_t objectNumber=0;
  if(numobject==0)
    objectNumber=1;
  const MapCellTriangle& hexaToTriangle1=*(O2.hexalist);

  MapCellTriangle::const_iterator it1=hexaToTriangle1.find(cell);
    
  if(it1!=hexaToTriangle1.end()) {
    std::list<Triangle*> listTriangle1=(*it1).second;
    for(std::list<Triangle*>::iterator jt0=listTriangle1.begin() ;
        jt0!=listTriangle1.end() ; ++jt0) {

      PairTriangleNew pairTriangle;
      if(numobject<objectNumber) {
        pairTriangle.first=(currentTriangle);
        pairTriangle.second=(*jt0);
      } else {
        pairTriangle.second=(currentTriangle);
        pairTriangle.first=(*jt0);
      }
      
      std::map< PairTriangleNew , IntersectionPoints>::iterator i = pairTriangleToIntersectionPoints.find(pairTriangle);
      if(i != pairTriangleToIntersectionPoints.end()) {
        TinyVector<3,size_t> cutTriangle;
        if(numobject<objectNumber) {
          cutTriangle=(*i).second.cutTriangle1();
        } else {
          cutTriangle=(*i).second.cutTriangle2();
        }

        for (size_t n=0; n<3; ++n) {
          if(cutTriangle[n]!=0) {
            edgesCut[n]=1;
          }
        }
      }
      else {
        // ffout(4)<<"la paire de triangle n'existe pas!\n";
      }
    }
  }
  return edgesCut;
}

void SurfaceMeshGenerator::Internals::
__putRefByFront(const size_t numobject,
                const MeshedObject& O1,
                const MeshedObject& O2,
                const SurfaceMeshOfTriangles& surfmesh,
                const std::set<const Cell*>& toTreatHexahedra) {
  //pour etre sur qu'on traite bien tous les tr qu'il faut
  std::map<const Triangle*,size_t> listTreat;
  //to put the reference by front
  std::list<const Triangle*> front;
  int frontSize = 0;

  size_t objectNumber=0;
  if(numobject==0)
    objectNumber=1;
  size_t nbcell=surfmesh.numberOfCells();
  //initialize the front
  ffout(4)<<"taille de totreatHexahedra = "<<toTreatHexahedra.size()<<"\n";
  size_t i=0;
  while(i<nbcell) {
    const Triangle* T=&(surfmesh.cell(i));
    const Cell* cell=&((*T).mother());
    if(toTreatHexahedra.find(cell)!=toTreatHexahedra.end()
       and ((*edgesRefVertex[numobject])[&(*T)(0)]==2
            or (*edgesRefVertex[numobject])[&(*T)(1)]==2
            or (*edgesRefVertex[numobject])[&(*T)(2)]==2)) {
      front.push_front(T);
      listTreat[T]=1;
      frontSize++;
      const Triangle*& TT=(*listOfTriangleMeshFront.find(T)).second;
      // on traite le triangle du front      
      __putEdgeRef(T,TT,numobject,O1,O2);
    }
    ++i;
  }
  
  size_t iter=0;
  ffout(4)<<"taille initial du front ="<< frontSize <<"\n";
  //if(front.size()==0)
  if(nbcell != 0) {
    ffout(4)<<"ref de surfmesh = "<<surfmesh.cell(0).reference()<<"\n";
  }
  //put reference by front
  size_t erase=0;

  const Connectivity<SurfaceMeshOfTriangles>& connect
    = surfmesh.connectivity();

  while(frontSize !=0 and iter<nbcell) {
    std::list<const Triangle*>::iterator er=front.begin() ;
    for(std::list<const Triangle*>::iterator jt0=front.begin() ;
        jt0!=front.end() ; ++jt0) {
      if(jt0!=front.begin() and erase==1) {
        er=front.erase(er);
        frontSize--;
      }
      er=jt0;
      erase=0;
      const Triangle* T=*jt0;
      const Triangle*& TTo=(*listOfTriangleMeshFront.find(T)).second;
      const Vertex& Vo0 = (*TTo)(0);
      const Vertex& Vo1 = (*TTo)(1);
      const Vertex& Vo2 = (*TTo)(2);
      TriangleCut Kneigho(Vo0,Vo1,Vo2,numobject,objectNumber);
      const Cell* cello=&((*TTo).mother());
      Kneigho.edgecut =__determinateEdgeCut(numobject,O1,O2,TTo,cello);
      TinyVector<2,int> edgeo0=Kneigho.edgecut[0];
      TinyVector<2,int> edgeo1=Kneigho.edgecut[1];
      TinyVector<2,int> edgeo2=Kneigho.edgecut[2];
      //__putEdgeRef(T,TTo,numobject,O1,O2);
      size_t& refo0=(*edgesRefVertex[numobject])[&(*T)(0)];
      size_t& refo1=(*edgesRefVertex[numobject])[&(*T)(1)];
      size_t& refo2=(*edgesRefVertex[numobject])[&(*T)(2)];
      if(((refo0==2 or refo1==2 or refo2==2) and  refo0+refo1+refo2<6 and edgeo0+edgeo1+edgeo2==0)
          /*or (refo0+refo1+refo2!=6 and edgeo0+edgeo1+edgeo2==0 and refo0+refo1+refo2>0)*/
         /*or (edgeo0[0]==1 and refo0+refo1!=2 and refo0+refo1+refo2!=0)
           or (edgeo1[0]==1 and refo0+refo2!=2 and refo0+refo1+refo2!=0)
           or (edgeo2[0]==1 and refo2+refo1!=2 and refo0+refo1+refo2!=0)*/) {
        //ffout(4)<<"chouette on garde "<<refo0+refo1+refo2<<"\n";
        refo0=2;
        refo1=2;
        refo2=2;
        //erase=1;
      } else {
        erase=1;
      }

      const TinyVector<3,const Triangle*>& V=connect.cells(*T);
      
      for(size_t j=0 ; j<3 ; ++j) {
        //check the neighbours exist
        if(V[j]!=0){
          const Triangle* triangle=V[j];
          const Triangle*& TT=(*listOfTriangleMeshFront.find(triangle)).second;
          const Vertex& V0 = (*TT)(0);
          const Vertex& V1 = (*TT)(1);
          const Vertex& V2 = (*TT)(2);
          TriangleCut Kneigh(V0,V1,V2,numobject,objectNumber);
          const Cell* cell=&((*TT).mother());
          Kneigh.edgecut =__determinateEdgeCut(numobject,O1,O2,TT,cell);
          size_t& ref0=(*edgesRefVertex[numobject])[&(*triangle)(0)];
          size_t& ref1=(*edgesRefVertex[numobject])[&(*triangle)(1)];
          size_t& ref2=(*edgesRefVertex[numobject])[&(*triangle)(2)];
          if(ref0+ref1+ref2!=6) {
            //ffout(0)<<"le voisin n'est pas de ref 2\n";
            //ffout(0)<<"avt "<<ref0<<" "<<ref1<<" "<<ref2<<"\n";
            size_t test=0;
            TinyVector<2,int> edge0=Kneigh.edgecut[0];
            TinyVector<2,int> edge1=Kneigh.edgecut[1];
            TinyVector<2,int> edge2=Kneigh.edgecut[2];
            size_t sum=edge0[0]+edge1[0]+edge2[0];
            if((ref0==2 or ref1==2 or ref2==2)  and
               (ref0!=1 and ref1!=1 and ref2!=1) and sum==0) {
              //ffout(4)<<"ref 2 partout\n";
              test=1;
              ref0=2;
              ref1=2;
              ref2=2;
            }  else {
              if(((ref1==2) xor (ref0==2)) and (edge0==0) and (ref1!=1) and (ref0!=1)) {
                  ref1=2;
                  ref0=2;
                  test=1;
              }
              if(((ref2==2) xor (ref0==2)) and (edge1==0) and (ref2!=1) and (ref0!=1)) {
                    ref2=2;
                    ref0=2;
                    test=1;
              }
              if(((ref1==2) xor (ref2==2)) and (edge2==0) and (ref1!=1) and (ref2!=1)) {
                    ref1=2;
                    ref2=2;
                    test=1;
              }
              if(((ref1==2) xor (ref0==2)) and (edge0==1) and (sum>1)){
                //ffout(4)<<"ref 1 2 0 partout\n";
                //ffout(4)<<ref0<<" "<<ref1<<"\n";
                if(ref1+ref0!=3) {
                  test=1;
                }
                if(ref0==2) {
                  ref1=1;
                } else {
                  ref0=1;
                }
              }
              if(((ref2==2) xor (ref0==2)) and (edge1==1) and (sum>1)) {
                //ffout(0)<<"ref 1 2 1 partout\n";
                //ffout(4)<<ref0<<" "<<ref2<<"\n";
                if(ref2+ref0!=3) {
                  test=1;
                }
                if(ref0==2) {
                  ref2=1;
                } else {
                  ref0=1;
                }
              }
              if(((ref2==2) xor (ref1==2)) and (edge2==1) and (sum>1)) {
                //ffout(0)<<"ref 1 2 2 partout\n";
                //ffout(4)<<ref1<<" "<<ref2<<"\n";
                if(ref1+ref2!=3) {
                  test=1;
                }
                if(ref2==2) {
                  ref1=1;
                } else {
                  ref2=1;
                }
              }
            }       
            if(test==1/*ref0+ref1+ref2==6 *//*or edge0[0]+edge1[0]+edge2[0]==1*/) {
              front.push_front(triangle);
              listTreat[triangle]=1;
              frontSize++;
            } else if(listTreat.find(triangle)==listTreat.end()) {
                //ffout(4)<<"pas encore traite\n";
                front.push_front(triangle);
                listTreat[triangle]=1;
                frontSize++;
            }
            //ffout(0)<<ref0<<" "<<ref1<<" "<<ref2<<"\n";
          }
        }
      }
    }
    if(erase==1 and frontSize==1) {
      er=front.erase(er);
      frontSize--;
    }
    ++iter;
    ffout(4)<<"taille du front = "<<front.size()<<" "<<iter<<" "<<nbcell<<"\n";
  }
  
}

void SurfaceMeshGenerator::Internals::
__putEdgeRef(const Triangle*& Tmesh,
             const Triangle*& T,
             const size_t numobject,
             const MeshedObject& O1,
             const MeshedObject& O2 )
{
  //put edgesRefVertex=2 if the vertex is in the domain
  //ffout(0)<<"j'initialise le front\n";

  TinyVector<3,TinyVector<2,int> > edgecut;
  const Cell* cell=&((*T).mother());
  edgecut =__determinateEdgeCut(numobject,O1,O2,T,cell);
  size_t& ref0=(*edgesRefVertex[numobject])[&(*Tmesh)(0)];
  size_t& ref1=(*edgesRefVertex[numobject])[&(*Tmesh)(1)];
  size_t& ref2=(*edgesRefVertex[numobject])[&(*Tmesh)(2)];
    
  if(ref0+ref1+ref2!=6) {
    TinyVector<2,int> edge0=edgecut[0];
    TinyVector<2,int> edge1=edgecut[1];
    TinyVector<2,int> edge2=edgecut[2];
    size_t sum=edge0[0]+edge1[0]+edge2[0];
    if((ref0==2 or ref1==2 or ref2==2)  and sum==0) {
      ref0=2;
      ref1=2;
      ref2=2;
    } else {
      if((ref1==2 xor ref0==2) and edge0==1 and sum>1){
        if(ref0==2) {
          ref1=1;
        } else {
          ref0=1;
        }
      }
      if((ref2==2 xor ref0==2) and edge1==1 and sum>1 ){
        if(ref0==2) {
          ref2=1;
        } else {
          ref0=1;
        }
      }
      if((ref2==2 xor ref1==2) and edge2==1 and sum>1 ){
        if(ref2==2) {
          ref1=1;
        } else {
          ref2=1;
        }
      }
    }
  }
}

void SurfaceMeshGenerator::Internals::
__createCase(const TriangleCut& K ,
             size_t& in, size_t& in1,
             TinyVector<3,size_t>& place)
{
  size_t numberIn=K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2];
    
  switch (numberIn) {
  case 0: {
    if(K.edgecut[0][0]==1) {
      place[0]=0;
      if(K.edgecut[1][0]==1) {
        place[1]=1;
        place[2]=2;
      } else {
        place[1]=2;
        place[2]=1;
      }
    } else {
      if(K.edgecut[1][0]==1) {
        place[0]=1;
        place[1]=2;
        place[2]=0;
      } else {
        place[0]=2;
        place[1]=0;
        place[2]=1;
      }
    }
    break;
  }
  case 1: {
    if(K.isInTriangle[1]==1) {
      in=1;
      in1=in;
      place[0]=2;
      place[1]=0;
      place[2]=1;
    } else if(K.isInTriangle[2]==1) {
      in=2;
      in1=in;
      place[0]=1;
      place[1]=2;
      place[2]=0;
    } else {
      in1=in;
      place[0]=0;
      place[1]=1;
      place[2]=2;
    }
    break;
  }
  case 2: {
    if(K.isInTriangle[1]==1 and K.isInTriangle[0]==1) {
      in=1;
      in1=0;
      place[0]=2;
      place[1]=1;
      place[2]=0;
    } else if(K.isInTriangle[2]==1 and K.isInTriangle[1]==1) {
      in=2;
      in1=1;
      place[0]=1;
      place[1]=0;
      place[2]=2;
    } else {
      in=0;
      in1=2;
      place[0]=0;
      place[1]=2;
      place[2]=1;
    }
    break;
  }
  case 3: {
      if(K.edgecut[0]==1) {
          place[0]=0;
          in1=2;
          in=0;
          if(K.edgecut[1]==1) {
              place[1]=1;
              place[2]=2;
          } else {
              place[1]=2;
              place[2]=1;
          }
      } else if(K.edgecut[1]==1) {
          in1=1;
          in=2;
          place[0]=1;
          place[1]=2;
          place[2]=0;
      } else {
          in1=0;
          in=1;
          place[0]=2;
          place[1]=1;
          place[2]=0;
      }
      break;
  }
  }
    
}

void SurfaceMeshGenerator::Internals::
__transformVertex(const TriangleCut& K,
                  std::vector<const Vertex*>& points,
                  std::vector<Vertex> &pointstrian)
{
    size_t taille=points.size();
    for(unsigned ie=0;ie<taille;++ie) {
        const Vertex& P=*points[ie];
        TinyMatrix<3,3, real_t> A;
        for (size_t i=0; i<3; ++i) {
            for (size_t j=0; j<3; ++j) {
                A(i,j)=(K(j))[i];
            }
        }
        TinyVector<3, real_t> lambda;
        gaussPivot(A,P,lambda);
        Vertex Pi;
        Pi[0]=lambda[0];
        Pi[1]=lambda[1];
        Pi[2]=0;
        pointstrian[ie]=Pi;
    }
    
}
void SurfaceMeshGenerator::Internals::
__addPoints(const Vertex*& Pcut,
            size_t& stopcut,
            std::vector<const Vertex*>& pointslist1,
            std::vector<const Vertex*>& pointslist2)
{
    if(stopcut==0) {
        //ffout(4)<<"00\n";
        pointslist1.push_back(Pcut);
        stopcut=3;
    } else {
        //ffout(4)<<"10\n";
        pointslist2.push_back(Pcut);
    }
}

void SurfaceMeshGenerator::Internals::
__addList(const TriangleCut& K)
{

// sdp: not used?    
//     size_t numberInEdges=K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]
//     +K.edgecut[0][1]+K.edgecut[1][1]+K.edgecut[2][1];
    size_t nbInter=vertexInTriangles.size();
    //ffout(4)<<"nb Inter avt = "<<nbInter<<"\n";
    //if(vertexVectTriangles.size()!=0) {
    //  ffout(4)<<"nb points total "<<vertexVectTriangles.size()<<"\n";
    //}

    //premiere etage : verifier que vertexInTriangles contient bien tous les points
    if(nbInter<vertexVectTriangles.size()) {
        //il manque des points
        for(std::map<const Vertex*,std::vector<const Triangle*> >::iterator it=vertexVectTriangles.begin()
            ; it!=vertexVectTriangles.end() ; ++it) {
            const Vertex* PP=(*it).first;
            //ffout(4)<<"-----\n";
            //ffout(4)<<*PP<<"\n";

            std::vector<const Triangle*> VT=(*it).second;
            if(vertexInTriangles.find(PP)==vertexInTriangles.end()) {
                //ffout(4)<<"je l'ajoute\n";
                //le point n'existe pas il faut le rajouter
                //on cherche a quel triangle l'associer
                if(nbInter==0) {
                    TinyVector<2,const Triangle*> triangles;
                    triangles[0]=VT[0];
                    triangles[1]=VT[0];
                    vertexInTriangles[PP]=triangles;
                } else {
                    size_t nn=0,stop2=0;
                    size_t taille=VT.size();
                    //test si le point est confondu avec un point in
                    if(K.isInTriangle[0]==1 and Norm(K(0)-*PP)<1e-6) {
                        ffout(4)<<"trop proche 1\n";
                        stop2=1;
                    }
                    if(K.isInTriangle[1]==1 and Norm(K(1)-*PP)<1e-6) {
                        ffout(4)<<"trop proche 2\n";
                        stop2=1;
                    }
                    if(K.isInTriangle[2]==1 and Norm(K(2)-*PP)<1e-6) {
                        ffout(4)<< Norm(K(2)-*PP)<<" trop proche 3\n";
                        ffout(4)<<*PP<<"\n";
                        ffout(4)<<K(2)<<"\n";
                        stop2=1;
                    }
                    while(nn<taille and stop2==0) {
                        const Triangle* TT=VT[nn];
                        if(triangleWithVertex.find(TT)!=triangleWithVertex.end()) {
                            (*triangleWithVertex.find(TT)).second[1]=PP;
                            TinyVector<2,const Triangle*> triangles;
                            triangles[0]=TT;
                            triangles[1]=TT;
                            vertexInTriangles[PP]=triangles;
                            stop2=1;
                        }
                        ++nn;
                    }
                    if(nn==taille and stop2==0) {
                        (*triangleWithVertex.find(&K)).second[0]=PP;
                        TinyVector<2,const Triangle*> triangles;
                        triangles[0]=&K;
                        triangles[1]=&K;
                        vertexInTriangles[PP]=triangles;
                        ffout(4)<<"bizarreeeeee on a pas trouve le tr\n";
                        //ffout(4)<<nn<<" =? "<<taille<<" "<<stop2<<" =? 0\n";
                        //ffout(4)<<*PP<<"\n";
                        //ffout(4)<<K(0)<<"\n";
                        //ffout(4)<<K(1)<<"\n";
                        //ffout(4)<<K(2)<<"\n";
                    }
                }
            } else {

            }
        }
    }
}

void SurfaceMeshGenerator::Internals::
__createList(const TriangleCut& K,
             std::vector<const Vertex*>& points,
             size_t& in,
             size_t& in1,
             TinyVector<3,size_t> place)
{
  TinyVector<2,const Vertex*> P;
  TinyVector<2,const Triangle*> T;
  size_t nbInter=vertexInTriangles.size();
  size_t nbCut=K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0];
  size_t numberIn=K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2];
  size_t stop=0,iInter=0,iCut=0;
  //initialisation
  //if(numberIn==3) {
  //    points.push_back(&K(1));
  //    points.push_back(&K(2));
  //}
  if((numberIn==3 or numberIn==2) and stop==0) {
    points.push_back(&K(in1));
  }
  const Vertex* A=&K(in);
  if(numberIn>0) {
    points.push_back(A);
    //ffout(4)<<*A<<"\n";
  }
  if(nbCut>0) {
    P[0]=K.pointsIntersection[place[0]][0];
    //attention si on a 2 points d'inter et un point in choisir le point le plus proche de A
    if(K.edgecut[place[0]][1]!=0) {
      P[1]=K.pointsIntersection[place[0]][1];
      if(numberIn>0) {
        real_t norm0=Norm(*A-*P[0]);
        real_t norm1=Norm(*A-*P[1]);
        if(norm1<norm0) {
          ffout(4)<<"--------------- cas inverse\n";
          P[0]=P[1];
        }
      }
    }
    if(P[0]!=A  and nbInter!=0) {
      points.push_back(P[0]);
    }
  } else if(nbCut!=0) {
    if(P[0]==A) {
      ffout(4)<<"pbs point confondu A\n";
    }
    std::map<const Vertex*,TinyVector<2,const Triangle*> >::iterator itvert=vertexInTriangles.begin();
    P[0]=(*itvert).first;
  }
  if(nbCut!=0) {
    if(vertexInTriangles.find(P[0])!=vertexInTriangles.end()) {
      T[0]=(*vertexInTriangles.find(P[0])).second[0];
      iInter=1;
    } else {
      ffout(4)<<"******\n";
      ffout(4)<<"nb Inter "<<nbInter<<"\n";
      ffout(4)<<"P0 = "<<P[0]<<*P[0]<<"\n";
      ffout(4)<<"K(2) = "<<K(2)<<&K(2)<<"\n";
      ffout(4)<<"pbs 0 existe pas\n";
      ffout(4)<<place<<"\n";
      //ffout(4)<<*K.pointsIntersection[0][0]<<"\n";
      //ffout(4)<<*K.pointsIntersection[1][0]<<"\n";
      ffout(4)<<K(0)<<"\n";
      ffout(4)<<K(1)<<"\n";
      ffout(4)<<K(2)<<"\n";
      points.clear();
      stop=1;
    }
  }
  
  //ffout(4)<<"nb de points d'inter "<<nbInter<<"\n";
  if(stop==0 and nbCut==2 ) {
    if(K.edgecut[place[1]][0]==1 and K.edgecut[place[0]][0]==1) {
      if(Norm(*K.pointsIntersection[place[1]][0]-*K.pointsIntersection[place[0]][0])<1e-12) {
        ffout(4)<<"les points sont confondus\n";
        ffout(4)<<*K.pointsIntersection[place[0]][0]<<"\n";
        stop=1;
      }
    }
  }

  if(iInter==nbInter and nbInter==1 and numberIn==1) {
    //ATTENTION PEUT ETRE PLUS DE POINTS A PRENDRE
    if(nbCut>1) {
      const Vertex* PP=K.pointsIntersection[place[1]][0];
      if(PP!=P[0]) {
        points.push_back(PP);
      } else {
        if(K.edgecut[place[2]]!=0) {
          const Vertex* PP1=K.pointsIntersection[place[2]][0];
          if(PP1!=P[0]) {
            points.push_back(PP1);
          }
        }
      }
    }
  }
  
  while(iInter<nbInter and stop==0) {
    //ffout(4)<<"createGeneral "<<iInter<<" "<<nbInter<<"\n";
    //on cherche le voisin de Pi
    __findVertex(K,P,T,iCut,iInter,stop,in,in1,place,points);

    //on cherche le triangle associe au nouveau Pi pour relancer
    size_t k=0;
    while(stop==3 or k==0) {
      //ffout(4)<<k<<" findTriangle\n";
      __findTriangle(P,T,stop);
      if(stop==3) {
        //ffout(4)<<"stop3\n";
        __findInEdges(K,P,T,iInter,stop,in,in1,place,points);
      }
      ++k;
    }
    
  }
    
}

bool SurfaceMeshGenerator::Internals::
__verifExist(const Vertex* Pt,
             std::vector<const Vertex*>& points)
{
  bool exist=false;
  for(size_t k=0 ; k<points.size() ; ++k) {
    if(points[k]==Pt) {
      //ffout(4)<<"je l'ai deja!!\n";
      exist=true;
    }
  }
  return exist;
}

bool SurfaceMeshGenerator::Internals::
__findOneEdge(const TriangleCut& K,
              const size_t & number,
              TinyVector<2,const Vertex*>& P)
{   bool find=false;
 if(K.edgecut[number][0]!=0) {
   if(P[0]==K.pointsIntersection[number][0]) {
     //ffout(4)<<"0"<<number<<"\n";
     find=true;
   }
 }
 if(K.edgecut[number][1]!=0 and find==false) {
   if(P[0]==K.pointsIntersection[number][1]) {
     //ffout(4)<<"1"<<number<<"\n";
     find=true;
   }
 }
 return find;
}


void SurfaceMeshGenerator::Internals::
__findInEdges(const TriangleCut& K,
              TinyVector<2,const Vertex*>& P,
              TinyVector<2,const Triangle*>& T,
              size_t& iInter,
              size_t& stop,
              size_t& in,
              size_t& in1,
              TinyVector<3,size_t> place,
              std::vector<const Vertex*>& points)
{
  //size_t nbCut=K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0];
  size_t nbInter=vertexInTriangles.size();
  size_t iCut=0;
  //ffout(4)<<*P[0]<<"\n";
  //find on which edge we are
  bool find=false;
  while(find==false and iCut<3) {
    find=__findOneEdge(K,iCut,P);
    ++iCut;
  }
  if(iCut==3 and find==false) {
    if(iInter!=iInter) {
      ffout(4)<<iInter<<" "<<iInter<<"\n";
      ffout(4)<<"pbs on a plus 3\n";
    }
    stop=1;
    iCut=0;
  } else {
    --iCut;
  }
    
  //ffout(4)<<"iCut "<<iCut<<"\n";
  if(K.edgecut[iCut][1]!=0 and stop!=1) {
    //verifier que le point n'est pas deja dans la liste
    const Vertex* Ptemp=K.pointsIntersection[iCut][1];
    bool exist=__verifExist(Ptemp,points);
    if(!exist) {
      P[0]=Ptemp;
      points.push_back(P[0]);
      //ffout(4)<<P[0]<<"\n";
      //ffout(4)<<*P[0]<<"\n";
      ++iInter;
    } else {
      exist=__verifExist(K.pointsIntersection[iCut][0],points);
      if(!exist) {
        P[0]=K.pointsIntersection[iCut][0];
        points.push_back(P[0]);
        //ffout(4)<<P[0]<<"\n";
        //ffout(4)<<*P[0]<<"\n";
        ++iInter;
      } else {
        if(nbInter!=iInter) {
          ffout(4)<<"pbs on a plus de points\n";
          ffout(4)<<nbInter<<" "<<iInter<<"\n";
          ffout(4)<<K(0)<<"\n";
          ffout(4)<<K(1)<<"\n";
          ffout(4)<<K(2)<<"\n";
          ffout(4)<<"dernier "<<iCut<<" "<<*K.pointsIntersection[iCut][0]<<"\n";
        }
        stop=1;
      }
    }
  } else {
    if(nbInter!=iInter) {
      ffout(4)<<stop<<" pbs on a plus de points2\n";
      ffout(4)<<nbInter<<" "<<iInter<<"\n";
      ffout(4)<<K(0)<<"\n";
      ffout(4)<<K(1)<<"\n";
      ffout(4)<<K(2)<<"\n";
    }
    stop=1;
  }
}

void SurfaceMeshGenerator::Internals::
__findVertex(const TriangleCut& K,
             TinyVector<2,const Vertex*>& P,
             TinyVector<2,const Triangle*>& T,
             size_t& iCut,
             size_t& iInter,
             size_t& stop,
             size_t& in,
             size_t& in1,
             TinyVector<3,size_t> place,
             std::vector<const Vertex*>& points)
{
  //size_t nbCut=K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0];
  //size_t nbInter=vertexInTriangles.size();
  
  if(triangleWithVertex.find(T[0])!=triangleWithVertex.end()) {
    TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[0])).second;
    if(VV[0]==P[0]) {
      if(VV[1]==P[0]) {
        //on est a priori sur une arete
          ffout(4)<<"cas findinedges\n";
        __findInEdges(K,P,T,iInter,stop,in,in1,place,points);
      } else {
        //ffout(4)<<"cas 2\n";
        P[1]=VV[1];
        P[0]=P[1];
        //ffout(4)<<P[0]<<"\n";
        //ffout(4)<<*P[0]<<"\n";
        points.push_back(P[0]);
        ++iInter;
      }
    } else {
      //ffout(4)<<"cas 3\n";
      P[1]=VV[0];
      P[0]=P[1];
      //ffout(4)<<P[0]<<"\n";
      //ffout(4)<<*P[0]<<"\n";
      bool exist=__verifExist(P[0],points);
      if(!exist) {
        points.push_back(P[0]);
        ++iInter;
      } else {
        ffout(4)<<"on a un pbs le point existe deja!!\n";
        stop=1;
      }
      

    }
  }
}

void SurfaceMeshGenerator::Internals::
__findTriangle(TinyVector<2,const Vertex*>& P,
               TinyVector<2,const Triangle*>& T,
               size_t& stop)
{
  if(stop==3) {
    stop=0;
  }
  TinyVector<2,const Triangle*> T2;
  if(vertexInTriangles.find(P[0])!=vertexInTriangles.end()) {
    T2[0]=(*vertexInTriangles.find(P[0])).second[0];
    T2[1]=(*vertexInTriangles.find(P[0])).second[1];
  } else {
    ffout(4)<<"pbs tr existe pas\n";
    stop=1;
  }

  if(T2[0]==T[0]) {
    if(T2[1]==T[0]) {
      stop=3;
    } else {
      T[0]=T2[1];
    }
  } else {
    T[0]=T2[0];
  }
    
}

void SurfaceMeshGenerator::Internals::
__createTriangle(const TriangleCut& K,const Cell* cell,
                 std::vector<Triangle>& triangleListIntersectionNew,
                 std::vector<const Vertex*>& points,
                 std::vector<Vertex>& pointstrian) {
//     static int toto=0;
    size_t taille=points.size();
    Triangulation trgen;
    Triangulation::CurveVertex envVertex(taille+1);
    std::ofstream sortie("testy2mb");
    sortie.precision(30);
    sortie<<taille<<" 0 0 \n";
    for(unsigned ie=0;ie<taille;++ie) {
        sortie<<pointstrian[ie][0]<<" "<<pointstrian[ie][1]<<"\n";
        //ffout(4)<<points[ie]<<"\n";
 //ffout(4)<<*points[ie]<<"\n";
 //ffout(4)<<pointstrian[ie]<<"\n";
 //ffout(4)<<"----\n";
        envVertex[ie] = &pointstrian[ie];
    }
    envVertex[taille] = envVertex[0];
    sortie<<taille<<" ";
    for(unsigned ie=0;ie<taille;++ie) {
        sortie<<envVertex[ie]-envVertex[0]<<" ";
    }

    sortie.close();

    //dans ce cas pas de trou
std::vector<Triangulation::CurveVertex> holes(0);
std::vector<Triangulation::CurveVertex> curves(0);
    
    if(taille>3) {
        //ffout(4)<<"je triangule\n";
       // ffout(4)<<"num = "<<++toto<<"\n";
        //verifions que les points ne s'intersectent pas avant de lancer la triangulation
        bool isOK=false;
        //ffout(4)<<"nb de points "<<taille<<"\n";
        size_t ie=1;
        while(ie<taille-2 and !isOK) {
            isOK=checkInterbis(pointstrian[0],pointstrian[taille-1],pointstrian[ie],pointstrian[ie+1]);
            ++ie;
        }
        if(!isOK) {
            //verification que les points ne sont pas trop proches!!
            ie=0;
            while(ie<taille-1 and !isOK) {
                const Vertex& A=pointstrian[ie];
                const Vertex& B=pointstrian[ie+1];
                real_t norm=Norm(A-B);
                isOK=(norm<1e-6);
                ++ie;
            }
        }
        if(isOK) {
            ffout(4)<<"isok est vraie-----------\n";
        }
        if(!isOK) {
          //  ffout(4)<<"je triangule\n";
            trgen.triangulize(pointstrian,
                       envVertex,
                       holes,
                       curves);
           // ffout(4)<<"fin triangule-----\n";
            std::ofstream o("full.mesh");
            trgen.export_mesh(o);
            o.close();
            std::list<ConnectedTriangle> L=trgen.getTriangles();
            for(std::list<ConnectedTriangle>::iterator itL=L.begin() ; itL!= L.end() ; itL++) {
                Triangle T=(*itL);
                Vertex* P0=&T(0);
                Vertex* P1=&T(1);
                Vertex* P2=&T(2);
                __addTriangle(Triangle(*points[P0-&pointstrian[0]],
                         *points[P1-&pointstrian[0]],
                         *points[P2-&pointstrian[0]],
                         K.reference()),cell,triangleListIntersectionNew);
            }
        } else {
            ffout(4)<<"isOK est vraie on a un pbs.....\n";
            ffout(4)<<K(0)<<"\n";
            ffout(4)<<K(1)<<"\n";
            ffout(4)<<K(2)<<"\n";
        }
        //ffout(4)<<L.size()<<" ecrit\n";
    } else if(taille==3) {
        //ffout(4)<<"je cree un triangle\n";
        //ffout(4)<<*points[0]<<"\n";
        //ffout(4)<<*points[1]<<"\n";
        //ffout(4)<<*points[2]<<"\n";
        __addTriangle(Triangle(*points[0],
                        *points[1],
                        *points[2],
                        K.reference()),cell,triangleListIntersectionNew);
    } else {
        ffout(4)<<"la taille est inferieur a 3 = "<<taille<<"\n";
        ffout(4)<<K(0)<<"\n";
        ffout(4)<<K(1)<<"\n";
        ffout(4)<<K(2)<<"\n";
        ffout(4)<<"¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡\n";
    }
}

void SurfaceMeshGenerator::Internals::
__create3in(const TriangleCut& K,const Cell* cell,
                std::vector<const Vertex*>& points,
                const size_t in,const size_t in1,
                const  TinyVector<3,size_t>& place,
                std::vector<Triangle>& triangleListIntersectionNew)
{
  size_t taille=points.size();
  size_t number=K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0];
  TinyVector<2,const Vertex*> P;
  size_t numberP=K.edgecut[place[0]][0]+K.edgecut[place[0]][1];
  switch (numberP) {
    case 1: {
        P[0]=K.pointsIntersection[place[0]][0];
        break;
    }
    case 2: {
        P[0]=K.pointsIntersection[place[0]][0];
        P[1]=K.pointsIntersection[place[0]][1];
        break;
    }
  }
  std::vector<const Vertex*> pointslist1;
  std::vector<const Vertex*> pointslist2; 
  size_t stopcut=0;
  for(unsigned ie=1;ie<taille;++ie) {
    const Vertex* Pcut=points[ie];
    //ffout(4)<<*Pcut<<"\n";
    if(Pcut==P[0]) {
        __addPoints(Pcut,stopcut,pointslist1,pointslist2);
        if(stopcut==3) {
          stopcut=1;
        }
    } else {
        if(number==2) {
          if(Pcut==P[1]) {
            __addPoints(Pcut,stopcut,pointslist1,pointslist2);
            if(stopcut==3) {
                stopcut=1;
            }
          } else {
            __addPoints(Pcut,stopcut,pointslist1,pointslist2);
            if(stopcut==3) {
                stopcut=0;
            }
          }
        } else {
          __addPoints(Pcut,stopcut,pointslist1,pointslist2);
          if(stopcut==3) {
            stopcut=0;
          }
        }
    }
  }
  
  //on  rajoute le 3e sommets
  size_t in2=0;
  if(in+in1==1) {
    in2=2;
  } else if(in+in1==2) {
    in2=1;
  }
  pointslist2.push_back(&K(in2));
  //on rajoute le point oppose a l'arete cut
  pointslist2.push_back(points[0]);
  pointslist1.push_back(points[0]);
  //on rajoute le premier point d'inter a la 2e liste
  if(!__verifExist(P[0],pointslist2)) {
    pointslist2.push_back(P[0]);
  } else {
    pointslist2.push_back(P[1]);
  }
  
  //verif des listes
  ffout(4)<<"taille de list1 = "<<pointslist1.size()<<"\n";
  for(unsigned ie=0;ie<pointslist1.size();++ie) {
    ffout(4)<<*pointslist1[ie]<<"\n";
  }
  ffout(4)<<"taille de list2 = "<<pointslist2.size()<<"\n";
  for(unsigned ie=0;ie<pointslist2.size();++ie) {
    ffout(4)<<*pointslist2[ie]<<"\n";
  }
  
  //on cree les 2 maillages
  if(pointslist1.size()>=3) {
    std::vector<Vertex> pointslist1plan(pointslist1.size());
    __transformVertex(K,pointslist1,pointslist1plan);
    __createTriangle(K,cell,triangleListIntersectionNew,pointslist1,pointslist1plan);
  }
  if(pointslist2.size()>=3) {
    std::vector<Vertex> pointslist2plan(pointslist2.size());
    __transformVertex(K,pointslist2,pointslist2plan);
    __createTriangle(K,cell,triangleListIntersectionNew,pointslist2,pointslist2plan);
  }
         
}

void SurfaceMeshGenerator::Internals::
__create2SD(const TriangleCut& K,const Cell* cell,
                std::vector<const Vertex*>& points,
                std::vector<Vertex> &pointstrian,
                const  TinyVector<3,size_t>& place,
                std::vector<Triangle>& triangleListIntersectionNew)
{
  // on veut decomposer la liste en 2 sous listes pour creer 2 sous domaines
  size_t number=K.edgecut[place[2]][0]+K.edgecut[place[2]][1];
  size_t taille=points.size();

  // pour determiner le ou les points ou on separe les listes
  TinyVector<2,const Vertex*> P;
  switch (number) {
  case 1: {
    P[0]=K.pointsIntersection[place[2]][0];
    break;
  }
  case 2: {
    P[0]=K.pointsIntersection[place[2]][0];
    P[1]=K.pointsIntersection[place[2]][1];
    break;
  }
  }
  std::vector<const Vertex*> pointslist1;
  std::vector<const Vertex*> pointslist2;
  
  size_t stopcut=0;
  for(unsigned ie=1;ie<taille;++ie) {
    const Vertex* Pcut=points[ie];
    //ffout(4)<<*Pcut<<"\n";
//     const Vertex& Pcutplan=pointstrian[ie];
    if(Pcut==P[0]) {
      __addPoints(Pcut,stopcut,pointslist1,pointslist2);
      if(stopcut==3) {
        stopcut=1;
      }
    } else {
      if(number==2) {
        if(Pcut==P[1]) {
          __addPoints(Pcut,stopcut,pointslist1,pointslist2);
          if(stopcut==3) {
            stopcut=1;
          }
        } else {
          __addPoints(Pcut,stopcut,pointslist1,pointslist2);
          if(stopcut==3) {
            stopcut=0;
          }
        }
      } else {
        __addPoints(Pcut,stopcut,pointslist1,pointslist2);
        if(stopcut==3) {
          stopcut=0;
        }
      }
    }
  }
  pointslist2.push_back(points[0]);
  //verif des listes
  //ffout(4)<<"taille de list1 = "<<pointslist1.size()<<"\n";
  //for(unsigned ie=0;ie<pointslist1.size();++ie) {
  //    ffout(4)<<*pointslist1[ie]<<"\n";
  //}
  //ffout(4)<<"taille de list2 = "<<pointslist2.size()<<"\n";
  //for(unsigned ie=0;ie<pointslist2.size();++ie) {
  //    ffout(4)<<*pointslist2[ie]<<"\n";
  //}
  //on cree les 2 maillages
  if(pointslist1.size()>=3) {
    std::vector<Vertex> pointslist1plan(pointslist1.size());
    __transformVertex(K,pointslist1,pointslist1plan);
    __createTriangle(K,cell,triangleListIntersectionNew,pointslist1,pointslist1plan);
  }
  if(pointslist2.size()>=3) {
    std::vector<Vertex> pointslist2plan(pointslist2.size());
    __transformVertex(K,pointslist2,pointslist2plan);
    __createTriangle(K,cell,triangleListIntersectionNew,pointslist2,pointslist2plan);
  }
  
}

void SurfaceMeshGenerator::Internals::
__createGeneral(const TriangleCut& K,const Cell* cell,
                std::vector<Triangle>& triangleListIntersectionNew)
{
  
  // pour l'instant on suppose edge cut une seule fois
  size_t in=0,in1=1;
  TinyVector<3,size_t> place;
  std::vector<const Vertex*> points;
  
  real_t epsilon=1e-6;
  real_t norm01=Norm(K(0)-K(1));
  real_t norm02=Norm(K(0)-K(2));
  real_t norm21=Norm(K(2)-K(1));
  //test si le triangle est plat!! cad si 2 points sont "confondus"
  if(norm01<epsilon or norm02<epsilon or norm21<epsilon) {
    ffout(4)<<"le triangle est tres tres tres plat\n";
  } else {
    //pour determiner le numero des edges cut et ceux des points in
    __createCase(K,in,in1,place);
    
    //pour 3 edges cut et un in il faut changer place
    if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==3 and
       K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==1) {
      size_t temp=place[1];
      place[1]=place[2];
      place[2]=temp;
    }
    
    //Vertex A(-0.000141,0.129744,0.110641);
    //if(Norm(A-K(0))<1e-6 or Norm(A-K(1))<1e-6 or Norm(A-K(2))<1e-6 ){
    //  ffout(4)<<"=============\n";
    //  ffout(4)<<K(0)<<"\n";
    //  ffout(4)<<K(1)<<"\n";
    //  ffout(4)<<K(2)<<"\n";
    //  ffout(4)<<K.isInTriangle<<"\n";
    //  ffout(4)<<K.edgecut<<"\n";
    //  ffout(4)<<place<<"\n";
    //    }
    
    
    //add vertex to vertexInTriangles if need
    __addList(K);
    
    //create the list "points"
    size_t numberOfPoints=vertexInTriangles.size()+K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2];
    if(numberOfPoints>3) {
      __createList(K,points,in,in1,place);
    } else if(numberOfPoints==3) {
      //dans ce cas on a rien besoin de faire on cree un triangle
      //ffout(4)<<"on a juste 3 points\n";
      std::vector<const Vertex*> PP;
      for(std::map<const Vertex*,TinyVector<2,const Triangle*> >::iterator jt=vertexInTriangles.begin()
            ; jt!=vertexInTriangles.end() ; jt++) {
        PP.push_back((*jt).first);
        //ffout(4)<<*(*jt).first<<"\n";
      }
      for(size_t i=0 ; i<3 ; ++i) {
        if(K.isInTriangle[i]==1) {
          PP.push_back(&K(i));
        }
      }
      ASSERT(PP.size()==3);
      real_t norm0=Norm(*PP[0]-*PP[1]);
      real_t norm1=Norm(*PP[0]-*PP[2]);
      real_t norm2=Norm(*PP[2]-*PP[1]);
      if(norm0>1e-6 and norm1>1e-6 and norm2>1e-6) {
        __addTriangle(Triangle(*PP[0],*PP[1],*PP[2],K.reference()),cell,triangleListIntersectionNew);
      }
      
    } else {
      ffout(4)<<"moins de 3 points\n";
    }
    
    size_t taille=points.size();
    // met les points dans le plan du triangle
    std::vector<Vertex> pointstrian(taille);
    //if(taille!=0) {
    //  ffout(4)<<"points size "<<taille<<"\n";
    //}
    //pour mettre les points dans un plan
    __transformVertex(K,points,pointstrian);
    
    if(K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==2 and
       K.edgecut[place[2]][0]!=0 and taille!=0) {
      ffout(4)<<"on est dans un cas a traiter pour le mailleur\n";
      ffout(4)<<K(0)<<"\n";
      ffout(4)<<K(1)<<"\n";
      ffout(4)<<K(2)<<"\n";
      // on veut decomposer la liste en 2 sous listes pour creer 2 sous domaines
      __create2SD(K,cell,points,pointstrian,place,triangleListIntersectionNew);
      
    } else if(K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==3) {
      ffout(4)<<"----------------- cas 3 in\n";
      ffout(4)<<"place "<<place<<"\n";
      ffout(4)<<K(0)<<"\n";
      ffout(4)<<K(1)<<"\n";
      ffout(4)<<K(2)<<"\n";
      
      size_t number=K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0];
      if(number==1) {
          __create3in(K,cell,points,in,in1,place,triangleListIntersectionNew);
        
      } else {
        ffout(4)<<"il faut le traiter\n";
      }
      //for(unsigned ie=0;ie<points.size();++ie) {
      //    ffout(4)<<*points[ie]<<"\n";
      //}
    } else if(taille!=0) {
      __createTriangle(K,cell,triangleListIntersectionNew,points,pointstrian);
    }
  }
}

template <typename booleanTest>
void SurfaceMeshGenerator::Internals::
__createTrianglesIntersection(const size_t numobject,
                              const MeshedObject& O1,
                              const MeshedObject& O2,
                              const std::set<const Cell*>& toTreatHexahedra,
                              std::vector<Triangle>& triangleListIntersectionNew,
                              PIntersection& listVertexIntersectionNew)
{
  //build mesh of object O1 (which is intersected with O2)

  size_t nbcell=(*O1.trianglelist).size();
  SurfaceMeshOfTriangles surfmesh(nbcell);
  // construction du maillage surfmesh a partir des listes connues
  // construction de listOfTriangleMeshFront
  __constructionFinalMesh(O1,surfmesh);
  
  ConnectivityBuilder<SurfaceMeshOfTriangles> builder(surfmesh);
  builder.generates(Connectivity<SurfaceMeshOfTriangles>::CellToCells);
  //ffout(4)<<"numero "<<numobject<<" "<<objectNumber<<'\n';
  size_t objectNumber=0;
  if(numobject==0)
    objectNumber=1;
  // (*edgesRefVertex[numobject]).clear();
  const MapCellTriangle& hexaToTriangle0=*(O1.hexalist);

  
  std::set<const Cell*> hexaList1;
  std::set<const Cell*> hexaListIn;

  for (MapCellTriangle::const_iterator i=(*O1.hexalist).begin(); i!= (*O1.hexalist).end(); ++i) {
    hexaList1.insert((*i).first);
  }
  //hexaListIn : hexa ou se passent les intersections
  std::set_difference (hexaList1.begin(), hexaList1.end(),
                       toTreatHexahedra.begin(), toTreatHexahedra.end(),
                       inserter(hexaListIn, hexaListIn.begin()));

  ffout(4)<<"taille hexalistIn "<<hexaListIn.size()<<"\n";
  for (std::set<const Cell*>::iterator icell = hexaListIn.begin(); icell != hexaListIn.end(); ++icell) {
    
    MapCellTriangle::const_iterator it0=hexaToTriangle0.find(*icell);
    
    if(it0!=hexaToTriangle0.end()) {
      //ffout(4)<<"ici\n";
      Cell& C = const_cast<Cell&>(*(*icell));
      for (size_t i=0; i<C.numberOfVertices(); ++i) {
        C(i).reference() = (O1.shape().inside(C(i))?1:0)
          + (O2.shape().inside(C(i))?1:0);
      }

      size_t ref;
      if (booleanTest::compute(C)) {
        ref=2;
      } else {
        ref=0;
      }
        
      std::list<Triangle*> listTriangle0=(*it0).second;
      for(std::list<Triangle*>::iterator jt0=listTriangle0.begin() ;
          jt0!=listTriangle0.end() ; ++jt0) {

        Triangle& T = (*(*jt0));
        for (size_t n=0; n<3; ++n) {
          // if(ref==2)
          // ffout(4)<<"ref ="<<ref<<"\n";
          (*edgesRefVertex[numobject])[&T(n)]=ref;
          //put ref 2 for points of surfmesh
          const Vertex* P0=(*listOfVertexMesh.find(&T(n))).second;
          (*edgesRefVertex[numobject])[P0]=ref;
        }

        if((*edgesRefVertex[numobject])[&T(0)]==2
           and (*edgesRefVertex[numobject])[&T(1)]==2
           and (*edgesRefVertex[numobject])[&T(2)]==2) {
          //si le triangle est dans le domaine on le construit
          __addTriangle(T, &T.mother(),triangleListIntersectionNew);
        }
      }
    }
  }

  //put the reference of points
  __putRefByFront(numobject,O1,O2,surfmesh,toTreatHexahedra);
  
  size_t nbinters=listVertexIntersectionNew.size();
  ffout(4)<<"size listVertexIntersectionNew "<<nbinters<<'\n';
  
  for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin();
      i!=toTreatHexahedra.end() ; ++i) {
    MapCellTriangle::const_iterator it0=hexaToTriangle0.find(*i);
    
    if(it0!=hexaToTriangle0.end()) {
      std::list<Triangle*> listTriangle0=(*it0).second;
      for(std::list<Triangle*>::iterator jt0=listTriangle0.begin() ;
          jt0!=listTriangle0.end() ; ++jt0) {
        const Triangle* T=(*jt0);
        const Vertex& V0 = (*T)(0);
        const Vertex& V1 = (*T)(1);
        const Vertex& V2 = (*T)(2);
        TriangleCut K(V0,V1,V2,numobject,objectNumber);
        K.reference()=(*T).reference();
        triangleWithVertex.clear();
        vertexInTriangles.clear();
        vertexVectTriangles.clear();
        
        //build the vectors K.pointsIntersection, K.pointsin with points intersection include in K
        //put K.edgecut[i]=1 if the edge i is cut
        //build vertexInTriangles, triangleWithVertex (for build mesh by front)
        //build vertexVectTriangles
        __createLocalListIntersection(objectNumber,O2,T,K,*i);
        std::map< const Vertex* ,const Vertex*> pointsin=K.pointsin;

        const Vertex* P0=(*listOfVertexMesh.find(&K(0))).second;
        const Vertex* P1=(*listOfVertexMesh.find(&K(1))).second;
        const Vertex* P2=(*listOfVertexMesh.find(&K(2))).second;
        (*edgesRefVertex[numobject])[&K(0)]=(*edgesRefVertex[numobject])[P0];
        (*edgesRefVertex[numobject])[&K(1)]=(*edgesRefVertex[numobject])[P1];
        (*edgesRefVertex[numobject])[&K(2)]=(*edgesRefVertex[numobject])[P2];
        for(size_t k=0 ; k<3 ; ++k) {
          if((*edgesRefVertex[numobject])[&K(k)]==2) {
            K.isInTriangle[k]=1;
          } else {
            K.isInTriangle[k]=0;
          }
        }
        if( K.isInTriangle[0]+  K.isInTriangle[1]+  K.isInTriangle[2]!=0 and
            K.isInTriangle[0]+  K.isInTriangle[1]+  K.isInTriangle[2]!=3
            and  K.edgecut[0]+ K.edgecut[1]+ K.edgecut[2]==0) {
          //cad qu'on a au moins un point in et un point out mais qu'aucune arete n'est coupee
          ffout(4)<<"cas pbssssss\n";
          ffout(4)<<(*edgesRefVertex[numobject])[&K(0)]<<" "<<(*edgesRefVertex[numobject])[&K(1)]<<
            " "<<(*edgesRefVertex[numobject])[&K(2)]<<" "<<K.edgecut<<"\n";
          ffout(4)<<K(0)<<"\n";
          ffout(4)<<K(1)<<"\n";
          ffout(4)<<K(2)<<"\n";
          
          // __addTriangle(Triangle(K(0),K(1),K(2),K.reference()),*i,triangleListIntersectionNew);

        }
        if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]>=6) {
          ffout(4)<<"on veut rien faire car triangle plat\n";
          ffout(4)<<"K.isInTriangle "<< K.isInTriangle[0]<<" "<< K.isInTriangle[1]<<" "<< K.isInTriangle[2]<<"\n";
          ffout(4)<<K(0)<<"\n";
          ffout(4)<<K(1)<<"\n";
          ffout(4)<<K(2)<<"\n";
          ffout(4)<<"edge cut "<<K.edgecut<<"\n";
          K.isInTriangle[0]=1;
          K.isInTriangle[1]=1;
          K.isInTriangle[2]=1;
            
        }
        size_t numberVertexInDomain=K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2];
        switch (numberVertexInDomain) {
        case 0: {
          if(K.edgecut[0]+K.edgecut[1]+K.edgecut[2]==3 and vertexVectTriangles.size()<3) {
            ffout(4)<<"tous les points d'inter ne sont pas dans la liste\n";
          } else {
            //ffout(4)<<"------ case 0 in\n";
            __createGeneral(K,*i,triangleListIntersectionNew);
          }   
          break;  
        } 
        case 1: {
          if(K.edgecut[0]+K.edgecut[1]+K.edgecut[2]==0) {
            ffout(4)<<"pbbbbbbbbbbbbbbb\n";
          }
          if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==3) {
            if(K.pointsIntersection[0][0]!=K.pointsIntersection[1][0]
               and K.pointsIntersection[1][0]!=K.pointsIntersection[2][0]
               and K.pointsIntersection[2][0]!=K.pointsIntersection[0][0]) {
              //ffout(4)<<"---------------- 3 edge cut \n";
              __createGeneral(K,*i,triangleListIntersectionNew);
            } else {
              //ffout(4)<<"1------------\n";
              //ffout(4)<<K.edgecut[0]<<" "<<K.edgecut[1]<<" "<<K.edgecut[2]<<"\n";
              __createGeneral(K, *i, triangleListIntersectionNew);
            }
          } else {
            //ffout(4)<<"2------------\n";
            //ffout(4)<<K.edgecut[0]<<" "<<K.edgecut[1]<<" "<<K.edgecut[2]<<"\n";
            __createGeneral(K, *i, triangleListIntersectionNew);
          }
          break;
        }
        case 2: {
          if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]<2) {
            ffout(4)<<"pbs pas assez d'edge cut \n";
          } else {
            //ffout(4)<<"---------- case 2\n";
            //ffout(4)<<K.edgecut[0]<<" "<<K.edgecut[1]<<" "<<K.edgecut[2]<<"\n";
            __createGeneral(K, *i, triangleListIntersectionNew);
          }
          break;
        }
        case 3: {
          if(K.edgecut[0]+K.edgecut[1]+K.edgecut[2]==0) {
            __addTriangle(*T,*i,triangleListIntersectionNew);
          } else {
              //ffout(4)<<"------------ 3 in cas traite\n";
              //ffout(4)<<K.edgecut[0]<<" "<<K.edgecut[1]<<" "<<K.edgecut[2]<<"\n";
              __createGeneral(K, *i, triangleListIntersectionNew);
          }
          break;
        }
        }
        //      // K.isInTriangle[i] vaut 1 si le point est a garder pour construire le maillage 0 sinon.
        //      // K.edgecut[i] vaut 1 si l'arete est coupee

      }
    }
  }
  (*edgesRefVertex[numobject]).clear();
  listOfVertexMesh.clear();
  listOfTriangleMeshFront.clear();
}

bool SurfaceMeshGenerator::Internals::
__createTriangleSurface(VerticesList& verticesListes,
                        EdgePairList& localTriangleList,
                        const Shape& S ,
                        std::vector<const Edge*>& cutEdges,
                        const Cell& currentCell)
{
  bool jobDone=true;
  switch (cutEdges.size()) {
  case 3: {
            
    for (size_t i=0; i<cutEdges.size(); ++i) {
      TinyVector<3> V=__splitEdge(*cutEdges[i], S);
      verticesListes[*cutEdges[i]]
        =  new Vertex(V);
    }
    
    localTriangleList.push_front(TinyVector<3,Edge::Pair>());
//     cellList.push_front(&currentCell);

    std::vector<TinyVector<3> > V(cutEdges.size());
    for(size_t i=0; i<cutEdges.size(); ++i) {
      const Edge& e = *cutEdges[i];
      if(e(0).reference() == 1) {
        V[i] = (e(0) - e(1));
      } else {
        V[i] = (e(1) - e(0));
      }
    }

    if ((V[0]^V[1])*V[2] < 0) {
      std::swap(cutEdges[1], cutEdges[2]);
    }

    SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle
      = localTriangleList.begin();

    for(size_t i=0; i<cutEdges.size(); ++i) {
      (*currentTriangle)[i] = (Edge::Pair)(*cutEdges[i]);
    }

    break;
  }

  case 4: {
    for (size_t i=0; i<4; ++i) {
      TinyVector<3> V=__splitEdge(*cutEdges[i], S);
      verticesListes[*cutEdges[i]]
        =  new Vertex(V);
    }
    
    std::map<Edge::Pair, const Edge*> cutEdgesMap;
    for(size_t ie=0; ie<cutEdges.size(); ++ie) {
      cutEdgesMap[(Edge::Pair)*(cutEdges[ie])] = (cutEdges[ie]);
    }

    localTriangleList.push_front(TinyVector<3,Edge::Pair>());
//     cellList.push_front(&currentCell);

    std::vector<TinyVector<3> > V(cutEdges.size());
    std::map<Edge::Pair, const Edge*>::iterator currentCutEdge=cutEdgesMap.begin();
    for (size_t i=0; i<V.size(); ++i) {
      const Edge& e = *(*currentCutEdge).second;
      if(e(0).reference() == 1) {
        V[i] = (e(0) - e(1));
      } else {
        V[i] = (e(1) - e(0));
      }
      ++currentCutEdge;
    }

    bool swapEdge = ((V[0]^V[1])*V[2] > 0);

    SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle
      = localTriangleList.begin();

    currentCutEdge = cutEdgesMap.begin();
    if (!swapEdge) {
      for(size_t i=0; i<3; ++i) {
        (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
        ++currentCutEdge;
      }
    } else {
      for(int i=2; i>=0; --i) {
        (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
        ++currentCutEdge;
      }
    }

    localTriangleList.push_front(TinyVector<3,Edge::Pair>());
//     cellList.push_front(&currentCell);

    currentTriangle = localTriangleList.begin();

    currentCutEdge = cutEdgesMap.begin();
    if (!swapEdge) {
      for(int i=2; i>=0; --i) {
        ++currentCutEdge;
        (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
      }
    } else {
      for(size_t i=0; i<3; ++i) {
        ++currentCutEdge;
        (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
      }
    }
    break;
  }
  case 1:
  case 2: {
    //jobDone = false;
    break;
  }
  }

  return jobDone;
}

void SurfaceMeshGenerator::Internals::
__createSurface(MeshedObject& O,
                EdgePairList& triangleListes,
                EdgePairList& localTriangleListesEdges,
                VerticesList& verticesListes,
                MotherCellList& cellListObject)
{
  MeshOfTetrahedra& tetrahedraMesh = (*__backgroundMesh);
  const Shape& S = O.shape();

  // Sets vertices references
  for (size_t i=0; i<tetrahedraMesh.numberOfVertices(); ++i) {
    Vertex& V = tetrahedraMesh.vertex(i);
    V.reference() = S.inside(V);
  }

  for (size_t i = 0; i<tetrahedraMesh.numberOfEdges(); i++) {
    Edge& e = tetrahedraMesh.edge(i);
#warning THIS TEST IS NOT ENOUGH! THINK TO THE BORDER!
    if(e(0).reference() != e(1).reference()) {
      verticesListes[e] = new Vertex(__splitEdge(e, S));
      e.reference() = 1;
    } else {
      e.reference() = 0;
    }
  }
    
  if (verticesListes.size()==0) {
    fferr(1) << "\nwarning :\n"
             << "meshing object " << S << " a problem occured ...\n"
             << "Check that:\n"
             << "- the object is not to small compare to the background mesh\n"
             << "- the object intersects the background mesh\n\n";
  }

#warning A garder ?
  EdgePairList localTriangleList;

  for (MeshOfTetrahedra::iterator i(tetrahedraMesh);
       not(i.end()); ++i) {
    Tetrahedron& T = (*i);
    size_t numberInside
      = T(0).reference()
      + T(1).reference()
      + T(2).reference()
      + T(3).reference();

    bool jobDone = true;
    if ((numberInside)%4 != 0) {
      std::vector<const Edge*> cutEdges;
      const Connectivity<MeshOfTetrahedra>::CellToEdgesType&
        tetrahedronEdge = tetrahedraMesh.connectivity().edges(T);

      for (unsigned l=0; l<Tetrahedron::NumberOfEdges; ++l) {
        const Edge& e = *(tetrahedronEdge[l]);
        if (e.reference() == 1) {
          cutEdges.push_back(&e);
        }
      }

      bool jobDone2=__createTriangleSurface(verticesListes,
                                            localTriangleList,
                                            S, cutEdges,
                                            T);

      if(!(jobDone2)) {
        jobDone=jobDone2;
      }
    }
            
    if(jobDone) {
      for (SurfaceMeshGenerator::Internals::EdgePairList::iterator t
             = localTriangleList.begin();
           t != localTriangleList.end(); ++t) {
        triangleListes.push_front(*t);
        cellListObject.push_front(&T);
      }
    } else {
      throw ErrorHandler(__FILE__,__LINE__,
                         "jobDone false",
                         ErrorHandler::unexpected);
    }
    localTriangleList.clear();
  }
}

#warning REMOVE THIS NOW !
void plot(const SurfaceMeshOfTriangles& s,
          std::set<const Cell*>& hlist,
          Structured3DMesh& SMesh);
void plot2(const SurfaceMeshOfTriangles& s,
           std::set<const Cell*>& hlist,
           Structured3DMesh& SMesh);

void plotmedit(const size_t& nobject,
               const SurfaceMeshOfTriangles& s_mesh,
               std::set<const Cell*>& toTreatHexahedra,
               size_t n0,size_t n1);

void plothexa(size_t ncase,std::set<const Cell*>& cuttedHexahedra);

void SurfaceMeshGenerator::Internals::
__setMotherCells(SurfaceMeshOfTriangles& surfaceMesh)
{
  for (SurfaceMeshOfTriangles::iterator i(surfaceMesh);
       not(i.end()); ++i) {
    Triangle& triangle = *i;
    const Tetrahedron& tetrahedron
      = static_cast<const Tetrahedron&>(triangle.mother());

    const size_t tetrahedronNumber
      = (*__backgroundMesh).cellNumber(tetrahedron);

    triangle.setMother((*__motherCells)[tetrahedronNumber],
                       std::numeric_limits<size_t>::max()); // Le triangle n'est pas une face de la maille
  }
}


void SurfaceMeshGenerator::Internals::
__generateMesh(const Domain& omega,
               const Object& object, const size_t& level,
               std::stack<ReferenceCounting<MeshedObject> >& objectStack)
{
  if(Information::instance().coarseMesh() and object.hasReference()) {
    ReferenceCounting<MeshedObject> o = new MeshedObject();
    (*o).setShape(object.shape()->getCopy());
    __marchingTetrahedra(*o);

    objectStack.push(o);
    ffout(4) << "Marching cube sur " << objectStack.top() << '\n';
  } else {
    switch ((*object.shape()).type()) {
    case Shape::union_: {
      const Union& U = static_cast<const Union&>(*(object.shape()));

      size_t j=0;
      ObjectTransformer objectTransformer(U.transformationsList());
      for (Union::const_iterator i = U.begin();
           i != U.end(); ++i,++j) {
        ReferenceCounting<Object> s = objectTransformer(*(*i));
        __generateMesh(omega, *s, level+1, objectStack);
        if (j>0) {
          // get first object of the intersection
          ReferenceCounting<MeshedObject> O1 = objectStack.top();
          objectStack.pop();

          // get second object of the intersection
          ReferenceCounting<MeshedObject> O2 = objectStack.top();
          objectStack.pop();

          ReferenceCounting<MeshedObject> O0 = new MeshedObject();
          __operationBoolean<UnionTest>(*O0,*O1, *O2);

          // Set shape function
          if (j == 1) { // create an union shape
            Union* u = new Union();
            u->push_back(new Object((*O1).shapeReference()));
            u->push_back(new Object((*O2).shapeReference()));
            (*O0).setShape(u);
          } else { //use the union shape contained in O2
            Union& u = dynamic_cast<Union&>(*(*O2).shapeReference());
            u.push_back(new Object((*O1).shapeReference()));
            (*O0).setShape((*O2).shapeReference());
          }

          objectStack.push(O0);
        }
      }
      break;
    }
    case Shape::intersection: {
      const Intersection& I = static_cast<const Intersection&>(*(object.shape()));

      size_t j=0;
      ObjectTransformer objectTransformer(I.transformationsList());
      for(Intersection::const_iterator i = I.begin();
          i != I.end(); ++i,++j) {
        ReferenceCounting<Object> s = objectTransformer(*(*i));
        __generateMesh(omega, *s, level+1, objectStack);
        if (j>0) {
          // get first object of the intersection
          ReferenceCounting<MeshedObject> O1 = objectStack.top();
          objectStack.pop();

          // get second object of the intersection
          ReferenceCounting<MeshedObject> O2 = objectStack.top();
          objectStack.pop();

          ReferenceCounting<MeshedObject> O0 = new MeshedObject();
          __operationBoolean<IntersectionTest>(*O0,*O1, *O2);

          if (j == 1) { // create an intersection shape
            Intersection* newI = new Intersection();
            newI->push_back(new Object((*O1).shapeReference()));
            newI->push_back(new Object((*O2).shapeReference()));
            (*O0).setShape(newI);
          } else { //use the intersection shape contained in O2
            Intersection& newI = dynamic_cast<Intersection&>(*(*O2).shapeReference());
            newI.push_back(new Object((*O1).shapeReference()));
            (*O0).setShape((*O2).shapeReference());
          }
          objectStack.push(O0);
        }
      }
      break;
    }
    case Shape::difference: {
      const Difference& D = static_cast<const Difference&>(*(object.shape()));

      ObjectTransformer objectTransformer(D.transformationsList());
      Intersection* I = new Intersection;

      Difference::const_iterator i = D.begin();
      ReferenceCounting<Object> s = objectTransformer(*(*i));
      I->push_back(s);
 
      Union* U = new Union;
      ++i;
      for(; i != D.end(); ++i) {
        s = objectTransformer(*(*i));
        U->push_back(s);
      }

      I->push_back(new Object(new Not(new Object(U))));
      Object o(I);
      __generateMesh(omega,  o, level+1, objectStack);
      break;
    }
    case Shape::not_: {
      const Object& notObject = *(static_cast<const Not&>(*object.shape()).object());

      __generateMesh(omega, notObject, level+1, objectStack);

      MeshedObject& o = (*objectStack.top());
      ReferenceCounting<Shape> s = o.shapeReference();
      o.setShape(new Not(new Object(s)));

      break;
    }
    case Shape::cube: {
      if (not(Information::instance().coarseMesh())) {
        const Cube& c = static_cast<const Cube&>(*object.shape());
        Intersection* I = new Intersection;
        for (size_t i=0; i<6; ++i) {
          I->push_back(new Object(new Plane(c, i)));
        }

        Object o(I);
        __generateMesh(omega, o, level+1, objectStack);
        break;
      }
    }
    case Shape::cylinder: {
      if (not(Information::instance().coarseMesh())) {
        const Cylinder& c = static_cast<const Cylinder&>(*object.shape());
        Intersection* I = new Intersection;
        I->push_back(new Object(new InfiniteCylinder(c)));
        I->push_back(new Object(new Plane(c, 0)));
        I->push_back(new Object(new Plane(c, 1)));

        Object o(I);
        __generateMesh(omega, o, level+1, objectStack);
        break;
      }
    }
    case Shape::cone: {
      if (not(Information::instance().coarseMesh())) {
        const Cone& c = static_cast<const Cone&>(*object.shape());
        Intersection* I = new Intersection;
        I->push_back(new Object(new InfiniteCone(c)));
        I->push_back(new Object(new Plane(c, 0)));
        I->push_back(new Object(new Plane(c, 1)));

        Object o(I);
        __generateMesh(omega, o, level+1, objectStack);
        break;
      }
    }
    default: {
      ReferenceCounting<MeshedObject> o = new MeshedObject();
      (*o).setShape(object.shape()->getCopy());
      __marchingTetrahedra(*o);

      objectStack.push(o);
      ffout(4) << "Marching cube sur " << objectStack.top() << '\n';
    }
    }
  }
  MeshedObject& o = (*objectStack.top());
  if(object.hasReference()) {
    const size_t ref = omega.reference(object.reference());
    for (size_t i=0; i<(*o.trianglelist).size(); ++i) {
      Triangle& t = (*o.trianglelist)[i];
      t.reference()=ref;
    }
  }
}

ReferenceCounting<Vector<Triangle> >
SurfaceMeshGenerator::Internals::
__marchingTetrahedra(MeshedObject& O)
{
  ffout(3) << "Marching tetrahedra\n";
  EdgePairList triangleListes;
  EdgePairList localTriangleListesEdges;
  VerticesList verticesListes;
  MotherCellList cellListObject;

  __createSurface(O, triangleListes, localTriangleListesEdges,
                  verticesListes, cellListObject);

  __dataStructureConvertion(O, triangleListes, verticesListes, cellListObject);
  ffout(3) << "Marching tetrahedra: done\n";
  return 0;
}

void SurfaceMeshGenerator::
generateSurfacicMesh(const Domain& omega,
                     const Mesh& M,
                     SurfaceMeshOfTriangles& s_mesh)
{
  switch (M.type()) {
  case Mesh::cartesianHexahedraMesh: {
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
                       "only cartesian meshes are supported for surface mesh generation",
                       ErrorHandler::unexpected);
  }
  }

  Internals& internals = (*__internals);

  // Check if a new Tetrahedral mesh is to generate
  if (&M != internals.__lastBackgoundMesh) {
    // puts a mark
    internals.__lastBackgoundMesh = &M; 
    
    // now computes the correct mesh
    MeshTetrahedrizor tetrahedrizor(&M);

    // the mesh is converted, but we do not need octree and connectivity
    tetrahedrizor.run(false);

    internals.__backgroundMesh
      = & static_cast<MeshOfTetrahedra&>(*tetrahedrizor.mesh());

    // Finally get the cell mapping: needed to impose boundary conditions.
    internals.__motherCells
      = tetrahedrizor.motherCells();

    ConnectivityBuilder<MeshOfTetrahedra> builder(*internals.__backgroundMesh);
    builder.generates(Connectivity<MeshOfTetrahedra>::CellToEdges);
  }

  const Object& object = omega.object();

//   const Scene& scene = omega.scene();
//   (*__internals).edgesRefVertex.resize(scene.nbObjects());
  (*__internals).edgesRefVertex.resize(2);

  for(size_t no=0 ; no<2/*scene.nbObjects()*/ ; ++no) {
    (*__internals).edgesRefVertex[no] = new std::map<const Vertex*,size_t> ;
  }

  std::stack<ReferenceCounting<SurfaceMeshGenerator::Internals::MeshedObject> > objectStack;

  (*__internals).__generateMesh(omega, object, 0, objectStack);

  SurfaceMeshGenerator::Internals::MeshedObject& o = (*objectStack.top());
  std::set<int> refset;
  for (size_t i=0; i<(*o.trianglelist).size(); ++i) {
    const Triangle& t = (*o.trianglelist)[i];
    refset.insert(t.reference());
  }

  ffout(0) << "Sets " << refset.size() << " References\n";

  //construction of s_mesh
  (*__internals).__constructionFinalMesh((*objectStack.top()),
                                         s_mesh);

  (*__internals).__setMotherCells(s_mesh);

  s_mesh.setBackgroundMesh(&M);
}

SurfaceMeshGenerator::SurfaceMeshGenerator()
  : __internals(new SurfaceMeshGenerator::Internals())
{
  ;
}

SurfaceMeshGenerator::~SurfaceMeshGenerator()
{
  ;
}

---