//  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.68 2004/05/15 18:05:40 delpinux Exp $

#include <SurfaceMeshGenerator.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 <HexahedronByEdges.hpp>
#include <TetrahedronByEdges.hpp>

#include <WorkingMesh.hpp>

#include <TinyMatrix.hpp>

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

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)

    {
      ;
    }

    //! Destructor.
    ~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 ObjectToTreat;
  friend class SurfaceMeshGenerator;

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

  typedef std::map<const Cell*, 
    std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,
    ReferenceCounting<std::vector<Edge> > > > HexaTetraAssociation;

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

  //! list of cut edges and intersection vertex.
  typedef std::map<Edge::Pair, size_t> EdgesRef;

  //! list of cut edges and intersection vertex.
  typedef std::map<Edge::Pair, Vertex*> VerticesList;

  typedef std::vector<ReferenceCounting<ObjectToTreat > > ListOfObjectToTreat;
  
  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(ObjectToTreat& O,
				 EdgePairList& triangleListes,
				 VerticesList& verticesListes,
				 MotherCellList& cellListObject);

  void __constructionVerticesList(const ObjectToTreat& O0);

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

  void __createCutEdges(const Scene& scene,
			std::vector<ReferenceCounting< std::vector<const Edge*> > > & cutEdges,
			ObjectReferences& otherReferences,
			const std::vector<TetrahedronByEdges>::iterator& currentT);

  void __calculatePointsIntersection(const ObjectToTreat& O1,
				     const ObjectToTreat& 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 ObjectToTreat& O2,
				     const Triangle*& currentTriangle,
				     TriangleCut & K,
				     const Cell* cell);
  
  TinyVector<3,TinyVector<2,int> > __determinateEdgeCut(const size_t& numobject,
							const ObjectToTreat& O1,
							const ObjectToTreat& O2,
							const Triangle*& currentTriangle,
							const Cell* cell);

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

  void __putRefByFront(const size_t numobject,
		       const ObjectToTreat& O1,
		       const ObjectToTreat& O2,
		       const SurfaceMeshOfTriangles& surfmesh,
		       const Connectivity<Triangle>& connect,
		       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 ObjectToTreat& O1,
				     const ObjectToTreat& 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,
			       MotherCellList& cellList);

  void __createSurface(Structured3DMesh& SMesh,
		       const Shape& S,
		       EdgePairList& triangleListes,
		       EdgePairList& localTriangleListesEdges,
		       VerticesList& verticesListes,
		       MotherCellList& cellListObject);

  bool __isDegenerate(const Triangle& T);

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

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

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

  ReferenceCounting<Vector<Triangle> >
  __marchingTetrahedra(const Structured3DMesh& SMesh,
		       ObjectToTreat& O);
    
  /// @todo FUNCTION TO REMOVE
  static void __getIntersectionReferences(const Intersection& I,
					  std::map<TinyVector<3>, Object*>& otherReferences);

  void plotini(size_t nobjectToTreat);

public:
  Internals()
  {
    ;
  }
};

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

  WorkingMesh __workingMesh;

public :
  /** 
   * The working mesh
   * 
   * @return __workingMesh
   */
  WorkingMesh& mesh()
  {
    return __workingMesh;
  }

  /** 
   * Read only access to the working mesh
   * 
   * @return __workingMesh
   */
  const WorkingMesh& mesh() const
  {
    return __workingMesh;
  }

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

  const ReferenceCounting<Shape>& shape() const
  {
    assert(__shape != 0);
    return __shape;
  }

  ReferenceCounting<Shape>& shape()
  {
    assert(__shape != 0);
    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 ObjectToTreat& operator=(const ObjectToTreat& O); 
  ObjectToTreat(const ObjectToTreat& O);
public:
  //constructor
  ObjectToTreat()
  {
    ;
  }

  /** 
   * Destructor
   *
   */
  ~ObjectToTreat()
  {
    ;
  }
};

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

  //! Destructor.
  ~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(ObjectToTreat& O0,
		   const ObjectToTreat& O1,
		   const ObjectToTreat& 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;

  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(ObjectToTreat& O0,
			const ObjectToTreat& O1,
			const ObjectToTreat& 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"<<std::flush;
  __createTrianglesIntersection<booleanTest>(1, O2, O1, toTreatHexahedra,
					     triangleListIntersectionNew, listVertexIntersectionNew);
  ffout(4)<<"nb triangles in intersection 1 "<<triangleListIntersectionNew.size()<<"\n"<<std::flush;

  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::
__createCutEdges(const Scene& scene,
		 std::vector<ReferenceCounting< std::vector<const Edge*> > > & cutEdges,
		 ObjectReferences& otherReferences,
		 const std::vector<TetrahedronByEdges>::iterator& currentT)
{

  //size_t nobjectToTreat = otherReferences.size();
  for(size_t ie=0; ie<6; ++ie) {
    const Edge& e = currentT->edge(ie);
    int num=0;
    for(ObjectReferences::iterator i = otherReferences.begin();
	i != otherReferences.end(); ++i) {
      for(size_t n = 0; n<scene.nbObjects(); ++n) {
	Reference::Representation otherRef;
	otherRef.resize(scene.nbObjects());
	if (scene.object(n).reference() == ((*i).first)) {
	  otherRef[n] = true;
	  if (e(0).reference() != e(1).reference()) {
	    (*cutEdges[num]).push_back(&e);
	  }
	}
      }
      num++;
    }
  }
}

void SurfaceMeshGenerator::Internals::__constructionVerticesList(const ObjectToTreat& 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 ObjectToTreat& 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"<<std::flush;
  //! copies vertices into the s_mesh.
  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"<<std::flush;
  
  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());
    listOfTriangleMeshFront[&s_mesh.cell(i)]=&T;
    i++;
    
  }
  // }
  //listOfVertexMesh.clear();
}

void SurfaceMeshGenerator::Internals::__dataStructureConvertion(ObjectToTreat& 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);
      triangleListbis[n] = T;
      hexaToTrianglebis[*itcell].push_back(&triangleListbis[n]);
      ++itcell;
    }
  }
  
  verticesListes.clear();
  triangleListes.clear();
  //}
}

void SurfaceMeshGenerator::Internals::
__calculatePointsIntersection(const ObjectToTreat& O1,
			      const ObjectToTreat& 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);
  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 ObjectToTreat& 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"<<std::flush;
	  }

	  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 ObjectToTreat& O1,
		     const ObjectToTreat& 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 ObjectToTreat& O1,
		const ObjectToTreat& O2,
		const SurfaceMeshOfTriangles& surfmesh,
		const Connectivity<Triangle>& connect,
		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;
  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(surfmesh.cellNumber(*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 ObjectToTreat& O1,
	     const ObjectToTreat& 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 ObjectToTreat& O1,
			      const ObjectToTreat& 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);
  
  Connectivity<Triangle> connect(nbcell);
  ConnectivityBuilder<SurfaceMeshOfTriangles> C(surfmesh,connect);
  //ffout(4)<<"numero "<<numobject<<" "<<objectNumber<<"\n"<<std::flush;
  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,connect,toTreatHexahedra);
  
  size_t nbinters=listVertexIntersectionNew.size();
  ffout(4)<<"size listVertexIntersectionNew "<<nbinters<<"\n"<<std::flush;
  
  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;
	    
	}
	///////////////////////creation des triangles//////////////////////////////
	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,
			MotherCellList& cellList)
{
  bool jobDone=true;
  switch (cutEdges.size()) {
  case 3: {
    //! Stores the Edges on which the vertex will be computed.
	    
    for (size_t i=0; i<cutEdges.size(); ++i) {
      TinyVector<3> V=__splitEdge(*cutEdges[i], S);
      verticesListes[*cutEdges[i]]
	=  new Vertex(V);
    }
    
    //! Addes a triangle to the list.
    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]);
    }

    //! Sets the triangle vertices (actually just store the edge where it lives).
    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);
    }
    
    //! Orders the Edges using a map !
    std::map<Edge::Pair, const Edge*> cutEdgesMap;
    for(size_t ie=0; ie<cutEdges.size(); ++ie) {
      cutEdgesMap[(Edge::Pair)*(cutEdges[ie])] = (cutEdges[ie]);
    }

    //! Addes a triangle to the list.
    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);

    //! Sets the triangle vertices (actually just store the edge where it lives).
    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;
      }
    }

    //! Addes a triangle to the list.
    localTriangleList.push_front(TinyVector<3,Edge::Pair>());
    cellList.push_front(&currentCell);

    currentTriangle = localTriangleList.begin();

    //! Sets the triangle vertices (actually just store the edge where it lives).
    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(Structured3DMesh& SMesh,
		const Shape& S,
		EdgePairList& triangleListes,
		EdgePairList& localTriangleListesEdges,
		VerticesList& verticesListes,
		MotherCellList& cellListObject)
{
  Vector<unsigned> references(SMesh.nbEdges());

  for (unsigned i=0; i<references.size(); ++i) {
    Edge& e = SMesh.edge(i);
    references[i] = e.reference();
    e.reference() = 0;
  }

  EdgePairList localTriangleList;


  for (size_t i = 0; i<SMesh.nbEdges(); i++) {
    Edge& e = SMesh.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;
    }
  }
    
  if (verticesListes.size()==0) {
    fferr(0) << "\nwarning : "
	     << "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";
  }

  for (size_t icell=0; icell<SMesh.shape().nx()-1; ++icell) {
    for (size_t jcell=0; jcell<SMesh.shape().ny()-1; ++jcell) {
      for (size_t kcell=0; kcell<SMesh.shape().nz()-1; ++kcell) {
	bool intersected = false;

	MotherCellList cellList;

	const Cell& currentCell = SMesh.cell(icell, jcell, kcell);
	    
	const size_t nx
	  = dynamic_cast<const Structured3DMesh&>(SMesh).shape().nx();

	const size_t nx_1 = nx-1;

	const size_t ny
	  = dynamic_cast<const Structured3DMesh&>(SMesh).shape().ny();

	const size_t ny_1 = ny-1;

	const size_t nz
	  = dynamic_cast<const Structured3DMesh&>(SMesh).shape().nz();

	const size_t nz_1 = nz-1;

	HexahedronByEdges H (SMesh.edge(icell*ny*nz + jcell*nz + kcell),
			     SMesh.edge(icell*ny*nz + (jcell+1)*nz + kcell),
			     SMesh.edge(icell*ny*nz + jcell*nz + kcell + 1),
			     SMesh.edge(icell*ny*nz + (jcell+1)*nz + kcell+1),
			     SMesh.edge(nx_1*ny*nz + icell*ny_1*nz + jcell*nz + kcell),
			     SMesh.edge(nx_1*ny*nz + (icell+1)*ny_1*nz + jcell*nz + kcell),
			     SMesh.edge(nx_1*ny*nz + icell*ny_1*nz + jcell*nz + kcell + 1),
			     SMesh.edge(nx_1*ny*nz + (icell+1)*ny_1*nz + jcell*nz + kcell+1),
			     SMesh.edge(nx_1*ny*nz + nx*ny_1*nz + icell*ny*nz_1 + jcell*nz_1 + kcell ),
			     SMesh.edge(nx_1*ny*nz + nx*ny_1*nz + (icell+1)*ny*nz_1 + jcell*nz_1 + kcell ),
			     SMesh.edge(nx_1*ny*nz + nx*ny_1*nz + icell*ny*nz_1 + (jcell+1)*nz_1 + kcell ),
			     SMesh.edge(nx_1*ny*nz + nx*ny_1*nz + (icell+1)*ny*nz_1 + (jcell+1)*nz_1 + kcell ));


	for (size_t l=0; l<H.numberOfEdges(); l++) {
	  Edge& e = H.edge(l);
	  if (verticesListes.find(e) != verticesListes.end()) {
	    intersected = true;
	    break;
	  }
	}

	if (intersected) { // Means that the hexahedron intersects the object.
	  int sens = (icell+jcell+kcell)%2;
	  ReferenceCounting<std::vector<Edge> > pAddedEdges;
	  pAddedEdges = new std::vector<Edge>(6);
	  std::vector<Edge>& addedEdges = *pAddedEdges;

	  std::vector<int> TverticesList(4);
	  switch (sens) {
	  case 0: {
	    TverticesList[0] = 1;
	    TverticesList[1] = 3;
	    TverticesList[2] = 4;
	    TverticesList[3] = 6;
	    break;
	  }
	  case 1: {
	    TverticesList[0] = 0;
	    TverticesList[1] = 2;
	    TverticesList[2] = 5;
	    TverticesList[3] = 7;
	  }
	  }

	  ReferenceCounting<std::vector<TetrahedronByEdges> > pT;
	  pT = new std::vector<TetrahedronByEdges>(5);

	  std::vector<TetrahedronByEdges>& T = *pT;

	  { // Create tetrahedra edge
	    std::vector<Edge>::iterator currentAddedEdge
	      = addedEdges.begin();
	    for (size_t i1=0; i1<3; ++i1) {
	      for (size_t i2=i1; i2<3; ++i2) {
		*(currentAddedEdge) = Edge(currentCell(TverticesList[i1]),
					   currentCell(TverticesList[i2+1]));
		Edge& e = (*currentAddedEdge);
		if (e(0).reference() != e(1).reference()) {
		  e.reference() = 1;
		}
		currentAddedEdge++;
	      }
	    }
	  }

	  switch(sens) {
	  case 0: {
	    T[0] = TetrahedronByEdges(addedEdges[0],
				      addedEdges[1],
				      addedEdges[2],
				      addedEdges[3],
				      addedEdges[4],
				      addedEdges[5]);

	    T[1] = TetrahedronByEdges(H.edge(0),
				      H.edge(4),
				      H.edge(8),
				      addedEdges[0],
				      addedEdges[1],
				      addedEdges[3]);

	    T[2] = TetrahedronByEdges(H.edge(1),
				      H.edge(5),
				      H.edge(11),
				      addedEdges[0],
				      addedEdges[2],
				      addedEdges[4]);

	    T[3] = TetrahedronByEdges(H.edge(2),
				      H.edge(7),
				      H.edge(9),
				      addedEdges[1],
				      addedEdges[2],
				      addedEdges[5]);

	    T[4] = TetrahedronByEdges(H.edge(3),
				      H.edge(6),
				      H.edge(10),
				      addedEdges[3],
				      addedEdges[4],
				      addedEdges[5]);
	    break;
	  }
	  case 1: {
	    T[0] = TetrahedronByEdges(addedEdges[0],
				      addedEdges[1],
				      addedEdges[2],
				      addedEdges[3],
				      addedEdges[4],
				      addedEdges[5]);

	    T[1] = TetrahedronByEdges(H.edge(0),
				      H.edge(5),
				      H.edge(9),
				      addedEdges[0],
				      addedEdges[1],
				      addedEdges[3]);

	    T[2] = TetrahedronByEdges(H.edge(1),
				      H.edge(4),
				      H.edge(10),
				      addedEdges[0],
				      addedEdges[2],
				      addedEdges[4]);

	    T[3] = TetrahedronByEdges(H.edge(2),
				      H.edge(6),
				      H.edge(8),
				      addedEdges[1],
				      addedEdges[2],
				      addedEdges[5]);

	    T[4] = TetrahedronByEdges(H.edge(3),
				      H.edge(7),
				      H.edge(11),
				      addedEdges[3],
				      addedEdges[4],
				      addedEdges[5]);
	  }
	  }

	  bool jobDone = true;
	  for (std::vector<TetrahedronByEdges>::iterator currentT = T.begin();
	       currentT < T.end(); ++currentT) {
	    std::vector<const Edge*> cutEdges;
	    
	    for(size_t ie=0; ie<6; ++ie) {
	      const Edge& e = currentT->edge(ie);
	      if (e.reference() == 1) {
		cutEdges.push_back(&e);
	      }
	    }

	    bool jobDone2=__createTriangleSurface(verticesListes,
						  localTriangleList,
						  S, cutEdges,
						  currentCell,
						  cellList);
	    if(!(jobDone2)) {
	      jobDone=jobDone2;
	    }
	  }
	    
	  if(jobDone) {
	    SurfaceMeshGenerator::Internals::MotherCellList::iterator itcell
	      =cellList.begin();
	    for (SurfaceMeshGenerator::Internals::EdgePairList::iterator t
		   = localTriangleList.begin();
		 t != localTriangleList.end(); ++t) {
	      triangleListes.push_front(*t);
	      cellListObject.push_front(*itcell);
	      ++itcell;
	    }
	  } else {
	    ffout(4)<<"------jobDone false\n";
	    SurfaceMeshGenerator::Internals::MotherCellList::iterator itcell
	      =cellList.begin();
	    for (EdgePairList::iterator t = localTriangleList.begin();
		 t != localTriangleList.end(); ++t) {
	      triangleListes.push_front(*t);
	      cellListObject.push_front(*itcell);
	      ++itcell;
	    }
	  }
	  localTriangleList.clear();
	}
      }
    }
  }

  for (unsigned i=0; i<references.size(); ++i) {
    SMesh.edge(i).reference() = references[i];
  }
}

#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::
__generateMesh(const Domain& omega,
	       const Structured3DMesh& SMesh,
	       const Object& object, const size_t& level,
	       std::stack<ReferenceCounting<ObjectToTreat> >& objectStack)
{
  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, SMesh, *s, level+1, objectStack);
      if (j>0) {
	// get first object of the intersection
	ReferenceCounting<ObjectToTreat> O1 = objectStack.top();
	objectStack.pop();

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

	ReferenceCounting<ObjectToTreat> O0 = new ObjectToTreat();
	__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).shape()));
	  u->push_back(new Object((*O2).shape()));
	  (*O0).setShape(u);
	} else { //use the union shape contained in O2
	  Union& u = dynamic_cast<Union&>(*(*O2).shape());
	  u.push_back(new Object((*O1).shape()));
	  (*O0).setShape((*O2).shape());
	}

	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, SMesh, *s, level+1, objectStack);
      if (j>0) {
	// get first object of the intersection
	ReferenceCounting<ObjectToTreat> O1 = objectStack.top();
	objectStack.pop();

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

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

	if (j == 1) { // create an intersection shape
	  Intersection* newI = new Intersection();
	  newI->push_back(new Object((*O1).shape()));
	  newI->push_back(new Object((*O2).shape()));
	  (*O0).setShape(newI);
	} else { //use the intersection shape contained in O2
	  Intersection& newI = dynamic_cast<Intersection&>(*(*O2).shape());
	  newI.push_back(new Object((*O1).shape()));
	  (*O0).setShape((*O2).shape());
	}
	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, SMesh, o,
		   level+1, objectStack);
    break;
  }
  case Shape::not_: {
    const Object& notObject = *(static_cast<const Not&>(*object.shape()).object());

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

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

    break;
  }
  case Shape::cube: {
    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, SMesh, o,
		   level+1, objectStack);
    break;
  }
  case Shape::cylinder: {
    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, SMesh, o,
		   level+1, objectStack);
    break;
  }
  case Shape::cone: {
    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, SMesh, o,
		   level+1, objectStack);
    break;
  }
  default: {
    ReferenceCounting<ObjectToTreat> o = new ObjectToTreat();
    (*o).setShape(object.shape());
    __marchingTetrahedra(SMesh, *o);

    objectStack.push(o);
    ffout(4) << "Marching cube sur " << objectStack.top() << '\n';
  }
  }

  ObjectToTreat& 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(const Structured3DMesh& mesh,
		     ObjectToTreat& O)
{
  ffout(3) << "Marching cube ..." << std::flush;
  EdgePairList triangleListes;
  EdgePairList localTriangleListesEdges;
  VerticesList verticesListes;
  MotherCellList cellListObject;

#warning VERY BAD HACK!
  Structured3DMesh& SMesh = const_cast<Structured3DMesh&>(mesh);

  SMesh.createReferences(*(O.shape()));

  __createSurface(SMesh, *(O.shape()), triangleListes, localTriangleListesEdges,
		  verticesListes, cellListObject);

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

/*! Generates a surfacic mesh using only characteristic function of an object
  and the Structured 3D Mesh.

  \warning Implementation could be improved and function size reduced.
*/
void SurfaceMeshGenerator::
generateSurfacicMesh(const Domain& omega,
		     const Mesh& M,
		     SurfaceMeshOfTriangles& s_mesh)
{
  switch (M.type()) {
  case Mesh::cartesianHexahedraMesh: {
    break;
  }
  default: {
    fferr(0) << __FILE__ << ':' << __LINE__ << ": Only cartesian structured 3d mesh are allowed\n";
    fferr(0) << " for surface mesh generation. Others are not already supported\n";
    std::exit(1);
  }
  }

  const Structured3DMesh& SMesh = dynamic_cast<const Structured3DMesh&>(M);

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

  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::ObjectToTreat> > objectStack;
  (*__internals).__generateMesh(omega, SMesh,object, 0, objectStack);

  SurfaceMeshGenerator::Internals::ObjectToTreat& 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);
  s_mesh.setBackgroundMesh(&M);

/*#warning _________________________________A ENLEVER_________
  Connectivity<Triangle> connect(s_mesh.numberOfCells());
  ConnectivityBuilder<SurfaceMeshOfTriangles> C(s_mesh,connect);
  size_t compt=0;
  for(size_t i=0 ; i<s_mesh.numberOfCells() ; ++i) {
    const Triangle& T=s_mesh.cell(i);
    TinyVector<3,const Triangle*> V=connect(s_mesh.cellNumber(T));
    for(size_t j=0 ; j<3 ; ++j) {
      if(V[j]==0) {
	//ffout(4)<<"attention un des voisins n'existe pas\n";
	++compt;
      }
    }
  }
  ffout(4)<<"****************************il manque "<<compt<<" sommets\n";*/
  /////////////////////////////////////////////////////////////////////////////

  int n0=0, n1=0;

#warning REMOVE THIS CALL
  std::set<const Cell*> toTreatHexahedra;
  plotmedit(scene.nbObjects(), s_mesh, toTreatHexahedra, n0, n1);
}

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

SurfaceMeshGenerator::~SurfaceMeshGenerator()
{
  ;
}

void plothexa(size_t ncase,std::set<const Cell*>& cuttedHexahedra) {
  char* cccc="hhhhhhhhhhhhh";
  if(ncase==0) {
    cccc="meshed1.mesh";
  }
  else if(ncase==1) {
    cccc="meshed2.mesh";
  }
  else {
    cccc="meshed3.mesh";
  }
  std::ofstream sortie3(cccc);
  sortie3<<"MeshVersionFormatted 1"<<" \n";
  sortie3<<" \n";
  sortie3<<"Dimension"<<" \n";
  sortie3<<"3"<<" \n";
  sortie3<<" \n";
  sortie3<<"Vertices"<<" \n";
  sortie3<<8*cuttedHexahedra.size()<<" \n";
  for (std::set<const Cell*>::iterator i = cuttedHexahedra.begin();
       i != cuttedHexahedra.end(); ++i) {
    const CartesianHexahedron& H = static_cast<const CartesianHexahedron&>(*(*i));
    for (size_t j=0; j<8; ++j)
      {
	sortie3<<H(j)[0]<<" "<<H(j)[1]<<" "<<H(j)[2]<<" 0"<<" \n";
      }
  }
  sortie3<<" \n";

  sortie3<<"End"<<" \n"; 
}

    


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

  std::ofstream sortie("intersection.mesh");
     
  sortie<<"MeshVersionFormatted 1"<<" \n";
  sortie<<" \n";
  sortie<<"Dimension"<<" \n";
  sortie<<"3"<<" \n";
  sortie<<" \n";
  sortie<<"Vertices"<<" \n";
  sortie<<8*toTreatHexahedra.size()/*+3*s_mesh.numberOfCells()*/+s_mesh.numberOfVertices()<<" \n";
  for (std::set<const Cell*>::iterator i = toTreatHexahedra.begin();
       i != toTreatHexahedra.end(); ++i) {
    const CartesianHexahedron& H = static_cast<const CartesianHexahedron&>(*(*i));
    for (size_t j=0; j<8; ++j) {
      sortie<<H(j)[0]<<" "<<H(j)[1]<<" "<<H(j)[2]<<" 0"<<" \n";
    }
  }
  for(size_t i=0;i<s_mesh.numberOfVertices();++i) {
    sortie<<s_mesh.vertex(i)[0]<<" "<<s_mesh.vertex(i)[1]<<" "<<s_mesh.vertex(i)[2]<<" 0"<<" \n";
  }

  sortie<<" \n";
  sortie<<"Triangles"<<"\n";
  sortie<<s_mesh.numberOfCells()<<"\n";
  ffout(4)<<"fin ecriture vertex\n";
  for (size_t i=0; i<s_mesh.numberOfCells(); ++i) {
    const Triangle& T = s_mesh.cell(i);
    // ffout(4)<<s_mesh.vertexNumber(T(0))<<" "<<s_mesh.vertexNumber(T(1))<<" "<<s_mesh.//vertexNumber(T(2))<<"\n";
    sortie<<s_mesh.vertexNumber(T(0))+8*toTreatHexahedra.size()+1<<" "<<s_mesh.vertexNumber(T(1))+8*toTreatHexahedra.size()+1<<" "<<s_mesh.vertexNumber(T(2))+8*toTreatHexahedra.size()+1;
    
    // sortie<<3*i+1+8*toTreatHexahedra.size()<<" "<<3*i+2+8*toTreatHexahedra.size()<<" //"<<3*i+3+8*toTreatHexahedra.size();//<<" "<<1<<"\n";
    //     if(i<n0)
    //       sortie<<" 1"<<"\n";
    //     else if(i<n1)
    //       sortie<<" 2"<<"\n";
    //     else
    //       sortie<<" 4"<<"\n";
    sortie << " " << T.reference() << '\n';
  }
  ffout(4)<<"fin ecriture triangles\n";
  sortie<<" \n";

  sortie<<"End"<<" \n";
  ffout(4)<<"fin ecriture intersection.mesh\n";
}