/[ff3d]/ff3d/solver/SurfaceMeshGenerator.cpp
ViewVC logotype

Diff of /ff3d/solver/SurfaceMeshGenerator.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.28 by delpinux, Fri Aug 8 08:39:56 2003 UTC revision 1.29 by delpinux, Wed Aug 27 16:49:52 2003 UTC
# Line 60  private: Line 60  private:
60    class IntersectionPoints    class IntersectionPoints
61    //class SurfaceMeshGenerator::Internals::IntersectionPoints    //class SurfaceMeshGenerator::Internals::IntersectionPoints
62    {    {
63    public :    private:
64      TinyVector<2,size_t> numpair;      std::vector<const Vertex* > __listOfPoints;
65      std::vector<const  Vertex* > listofPoints;      TinyVector<3, size_t> __cutTriangle1;
66      TinyVector< 3,size_t > cutTriangle1;      TinyVector<3, size_t> __cutTriangle2;
     TinyVector< 3,size_t > cutTriangle2;  
67    
68      inline const size_t size() const {    public:
69        return listofPoints.size();  
70        TinyVector<3, size_t>& cutTriangle1()
71        {
72          return __cutTriangle1;
73      }      }
74    
75      void addPoint(const Vertex*&  V) {      std::vector<const Vertex*>& listOfPoints()
76        listofPoints.push_back( V );      {
77          return __listOfPoints;
78        }
79    
80        TinyVector<3, size_t>& cutTriangle2()
81        {
82          return __cutTriangle2;
83        }  
84    
85        const TinyVector<3, size_t>& cutTriangle1() const
86        {
87          return __cutTriangle1;
88        }
89    
90        const TinyVector<3, size_t>& cutTriangle2() const
91        {
92          return __cutTriangle2;
93        }
94    
95        inline size_t size() const
96        {
97          return __listOfPoints.size();
98        }
99    
100        const Vertex*& operator[](const size_t& i)
101        {
102          assert(i<__listOfPoints.size());
103          return __listOfPoints[i];
104        }
105    
106        void addPoint(const Vertex*&  V)
107        {
108          __listOfPoints.push_back( V );
109        }
110        const IntersectionPoints& operator=(const IntersectionPoints& P)
111        {
112          __listOfPoints=P.__listOfPoints;
113          __cutTriangle1=P.__cutTriangle1;
114          __cutTriangle2=P.__cutTriangle2;
115          return *this;
116      }      }
     
117    
118      //constructor      //constructor
119      IntersectionPoints() {      IntersectionPoints()
120        numpair=0;        : __cutTriangle1(0),
121        cutTriangle1=0;          __cutTriangle2(0)
122        cutTriangle2=0;      {
123          ;
124      }      }
125    
126      IntersectionPoints(const IntersectionPoints& P) {  
127        numpair=P.numpair;      IntersectionPoints(const IntersectionPoints& P)
128        listofPoints=P.listofPoints;        : __listOfPoints(P.__listOfPoints),
129        cutTriangle1=P.cutTriangle1;          __cutTriangle1(P.__cutTriangle1),
130        cutTriangle2=P.cutTriangle2;          __cutTriangle2(P.__cutTriangle2)
131    
132        {
133          ;
134      }      }
135    
136      //! Destructor.      //! Destructor.
# Line 121  private: Line 165  private:
165      }      }
166    };    };
167    
   
168    class TriangleCut;    class TriangleCut;
169    class ObjectToTreat;    class ObjectToTreat;
170    friend class SurfaceMeshGenerator;    friend class SurfaceMeshGenerator;
171    
172    struct __ltv    typedef std::map<TinyVector<3,real_t>, std::vector<SurfaceMesh*> > VRRefMesh;
   {  
     bool operator()(const TinyVector<3,real_t>& s1, const TinyVector<3,real_t>& s2) const  
     {  
       if (s1[0]<s2[0])  
         return true;  
       if (s1[1]<s2[1])  
         return true;  
       if (s1[2]<s2[2])  
         return true;  
       return false;  
     }  
   };  
   
   typedef std::map<TinyVector<3,real_t>,  
                    std::vector<SurfaceMesh*>, __ltv> VRRefMesh;  
173    
174    typedef std::list<TinyVector<3,Edge::Pair> > EdgePairList;          typedef std::list<TinyVector<3,Edge::Pair> > EdgePairList;      
175    typedef std::list<const Cell*> MotherCellList;    typedef std::list<const Cell*> MotherCellList;
# Line 150  private: Line 178  private:
178                     std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,                     std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,
179                               ReferenceCounting<std::vector<Edge> > > > HexaTetraAssociation;                               ReferenceCounting<std::vector<Edge> > > > HexaTetraAssociation;
180    
181    typedef std::map<TinyVector<3>, Object*, SurfaceMeshGenerator::Internals::__ltv> ObjectReferences;    typedef std::map<TinyVector<3>, Object*> ObjectReferences;
182    
183    //! list of cut edges and intersection vertex.    //! list of cut edges and intersection vertex.
184    typedef std::map<Edge::Pair, size_t> EdgesRef;    typedef std::map<Edge::Pair, size_t> EdgesRef;
185    
   //! This list is used to provide an orientation to the Surfacic Elements.  
   std::list<TinyVector<3> > InsideVertexList;  
   
186    //! list of cut edges and intersection vertex.    //! list of cut edges and intersection vertex.
187    typedef std::map<Edge::Pair, Vertex*> VerticesList;    typedef std::map<Edge::Pair, Vertex*> VerticesList;
188    
   HexaTetraAssociation hexaTetraAssociation;  
   std::set<const Cell*> toTreatHexahedra;  
   
189    typedef std::vector<ReferenceCounting<ObjectToTreat > > ListOfObjectToTreat;    typedef std::vector<ReferenceCounting<ObjectToTreat > > ListOfObjectToTreat;
   
   ListOfObjectToTreat listOfObjectToTreat;  
   ListOfObjectToTreat listOfObjectToTreat2;  
   std::map<const Vertex*,const Vertex*>   listOfVertexMesh;  
190        
191    typedef std::map<const Cell*, std::list<Triangle*> > MapCellTriangle;    typedef std::map<const Cell*, std::list<Triangle*> > MapCellTriangle;
192    
193    typedef std::pair<const Triangle*,const Triangle*> PairTriangleNew;    typedef std::pair<const Triangle*,const Triangle*> PairTriangleNew;
194    //std::vector<Vertex*> listVertexIntersectionNew;    //std::vector<Vertex*> listVertexIntersectionNew;
195    
196    typedef std::map<std::pair<const Vertex*,const Vertex*> ,    typedef std::map<std::pair<const Vertex*,const Vertex*> ,
197                     TinyVector<2,const Vertex* > > PIntersection;                     TinyVector<2,const Vertex* > > PIntersection;
   PIntersection listVertexIntersectionNew;  
198    
199    typedef std::map< PairTriangleNew , IntersectionPoints> PairTriangleToIntersectionPoints;    typedef std::map< PairTriangleNew , IntersectionPoints> PairTriangleToIntersectionPoints;
200    
201      std::map<const Vertex*, const Vertex*>   listOfVertexMesh;
202    
203    PairTriangleToIntersectionPoints pairTriangleToIntersectionPoints;    PairTriangleToIntersectionPoints pairTriangleToIntersectionPoints;
204    
205    std::vector<ReferenceCounting< std::map<const Vertex*,size_t> > > edgesRefVertex;    std::vector<ReferenceCounting< std::map<const Vertex*,size_t> > > edgesRefVertex;
206    
   std::vector<Triangle> triangleListIntersectionNew;  
207    
208    std::map<const Vertex*,TinyVector<2,const Triangle*> > vertexInTriangles;    std::map<const Vertex*,TinyVector<2,const Triangle*> > vertexInTriangles;
209    std::map<const Triangle*,TinyVector<2,const Vertex*> > triangleWithVertex;    std::map<const Triangle*,TinyVector<2,const Vertex*> > triangleWithVertex;
210    
   MotherCellList cellList;  
   
211    TinyVector<3> __splitEdge(const Edge& e,    TinyVector<3> __splitEdge(const Edge& e,
212                              const Shape& S);                              const Shape& S);
213    
# Line 209  private: Line 227  private:
227                          const std::vector<TetrahedronByEdges>::iterator& currentT);                          const std::vector<TetrahedronByEdges>::iterator& currentT);
228    
229    void __calculatePointsIntersection(const ObjectToTreat& O1,    void __calculatePointsIntersection(const ObjectToTreat& O1,
230                                       const ObjectToTreat& O2);                                       const ObjectToTreat& O2,
231                                         const std::set<const Cell*>& toTreatHexahedra,
232                                         PIntersection& listVertexIntersectionNew);
233    
234    void __setBoundaryType(const Boundary& B,    void __setBoundaryType(const Boundary& B,
235                           POVRayReferences& vrRefList,                           POVRayReferences& vrRefList,
# Line 232  private: Line 252  private:
252                            const Vertex* & Vo0,                            const Vertex* & Vo0,
253                            const Vertex* & Vo1,                            const Vertex* & Vo1,
254                            TinyVector<3, const Vertex* > & V,                            TinyVector<3, const Vertex* > & V,
255                            IntersectionPoints& Pinter);                            IntersectionPoints& Pinter,
256    void __CalculPoint(const size_t& objectNumber,                            PIntersection& listVertexIntersectionNew);
257                       const size_t& numobject,  
258                       std::list<Triangle*>::iterator & currentTriangle,    void __CalculPoint(std::list<Triangle*>::iterator & currentTriangle,
259                       std::list<Triangle*>::iterator & currentTriangle1,                       std::list<Triangle*>::iterator & currentTriangle1,
260                       IntersectionPoints & Pinter);                       IntersectionPoints & Pinter,
261                         PIntersection& listVertexIntersectionNew);
262    
263    void __createPointIntersection(const size_t& numobject,    void __createPointIntersection(std::list<Triangle*>::iterator currentTriangle1,
264                                   const size_t& objectNumber,                                   std::list<Triangle*>::iterator currentTriangle2,
265                                   std::list<Triangle*>::iterator currentTriangle1,                                   PIntersection& listVertexIntersectionNew);
                                  std::list<Triangle*>::iterator currentTriangle);  
266        
267    const Vertex& __determinateNumber(const int& p,    const Vertex& __determinateNumber(const int& p,
268                                      const TinyVector<3,TinyVector<2,const Vertex*> >& pointsint);                                      const TinyVector<3,TinyVector<2,const Vertex*> >& pointsint);
# Line 254  private: Line 274  private:
274    void __addTriangle(const Vertex& V0,    void __addTriangle(const Vertex& V0,
275                       const Vertex& V1,                       const Vertex& V1,
276                       const Vertex& V2,                       const Vertex& V2,
277                       const Cell* cell);                       const Cell* cell,
278                         std::vector<Triangle>& triangleListIntersectionNew);
279    
280    void __createCase1(int ncase,    void __createCase1(int ncase,
281                       const TriangleCut& K,                       const TriangleCut& K,
282                       const Cell* cell);                       const Cell* cell,
283                         std::vector<Triangle>& triangleListIntersectionNew);
284    
285    void __createCase2(int ncase,    void __createCase2(int ncase,
286                       const TriangleCut& K,                       const TriangleCut& K,
287                       const Cell* cell);                       const Cell* cell,
288                         std::vector<Triangle>& triangleListIntersectionNew);
289    
290    void __createLocalListIntersection(const size_t& objectNumber,    void __createLocalListIntersection(const size_t& objectNumber,
291                                       const ObjectToTreat& O2,                                       const ObjectToTreat& O2,
# Line 281  private: Line 304  private:
304                                       TriangleCut& K);                                       TriangleCut& K);
305    
306    void __createListTriangle(const TriangleCut& K,    void __createListTriangle(const TriangleCut& K,
307                              const Cell* cell);                              const Cell* cell,
308                                std::vector<Triangle>& triangleListIntersectionNew);
309    
310    void __createByEdges(const TriangleCut& K,const Cell* cell) ;    void __createByEdges(const TriangleCut& K,const Cell* cell,
311                           std::vector<Triangle>& triangleListIntersectionNew);
312    
313    void  __createNotIn(const TriangleCut& K,const Cell* cell);    void  __createNotIn(const TriangleCut& K,const Cell* cell,
314                          std::vector<Triangle>& triangleListIntersectionNew);
315    
316    void  __createNotIn2(const TriangleCut& K,const Cell* cell);    void  __createNotIn2(const TriangleCut& K,const Cell* cell,
317                           std::vector<Triangle>& triangleListIntersectionNew);
318    
319    template <typename booleanTest>    template <typename booleanTest>
320    void __createTrianglesIntersection(const size_t numobject,    void __createTrianglesIntersection(const size_t numobject,
321                                       const ObjectToTreat& O1,                                       const ObjectToTreat& O1,
322                                       const ObjectToTreat& O2);                                       const ObjectToTreat& O2,
323                                         const std::set<const Cell*>& toTreatHexahedra,
324                                         std::vector<Triangle>& triangleListIntersectionNew,
325                                         PIntersection& listVertexIntersectionNew);
326    
327    bool __createTriangleSurface(VerticesList& verticesListes,    bool __createTriangleSurface(VerticesList& verticesListes,
328                                 EdgePairList& localTriangleList,                                 EdgePairList& localTriangleList,
329                                 EdgesRef& edgesRef,                                 EdgesRef& edgesRef,
330                                 const Shape& S ,                                 const Shape& S ,
331                                 std::vector<const Edge*>& cutEdges,                                 std::vector<const Edge*>& cutEdges,
332                                 const Cell& currentCell);                                 const Cell& currentCell,
333                                   MotherCellList& cellList);
334    
335    void __createSurface(Structured3DMesh& SMesh,    void __createSurface(Structured3DMesh& SMesh,
336                         const Shape& S,                         const Shape& S,
# Line 311  private: Line 342  private:
342    
343    bool __isDegenerate(const Triangle& T);    bool __isDegenerate(const Triangle& T);
344    
345    void __casePinterNotInList(const TriangleCut& K,const Cell* cell);    void __casePinterNotInList(const TriangleCut& K,const Cell* cell, std::vector<Triangle>& triangleListIntersectionNew);
346    
347    TinyVector<3,const Vertex*> __determinatePinterNotInList(const TriangleCut& K, TinyVector<3,size_t>);    TinyVector<3,const Vertex*> __determinatePinterNotInList(const TriangleCut& K, TinyVector<3,size_t>);
348    
349    template <typename booleanTest>    template <typename booleanTest>
350    void __calculateIntersection(ObjectToTreat& O0,    void __calculateIntersection(ObjectToTreat& O0,
351                                 const ObjectToTreat& O1,                                 const ObjectToTreat& O1,
352                                 const ObjectToTreat& O2);                                 const ObjectToTreat& O2,
353                                   const std::set<const Cell*>& toTreatHexahedra);
354    
355    template <typename booleanTest>    template <typename booleanTest>
356    void __operationBoolean(ObjectToTreat& O0,    void __operationBoolean(ObjectToTreat& O0,
357                            const ObjectToTreat& O1,                            const ObjectToTreat& O1,
358                            ObjectToTreat& O2);                            const ObjectToTreat& O2);
359    
360    void __generateMesh(Structured3DMesh& SMesh,    void __generateMesh(Structured3DMesh& SMesh,
361                        const Object& object, const size_t& level,                        const Object& object, const size_t& level,
# Line 339  private: Line 371  private:
371            
372    /// @todo FUNCTION TO REMOVE    /// @todo FUNCTION TO REMOVE
373    static void __getIntersectionReferences(const Intersection& I,    static void __getIntersectionReferences(const Intersection& I,
374                                            std::map<TinyVector<3>, Object*, Internals::__ltv>& otherReferences);                                            std::map<TinyVector<3>, Object*>& otherReferences);
375    void plotini(size_t nobjectToTreat);    void plotini(size_t nobjectToTreat);
376    
377  public:  public:
# Line 350  public: Line 382  public:
382  };  };
383    
384  class SurfaceMeshGenerator::Internals::ObjectToTreat  class SurfaceMeshGenerator::Internals::ObjectToTreat
 //   : public Object  
385  {  {
386  private:  private:
387    ReferenceCounting<Shape> __shape;    ReferenceCounting<Shape> __shape;
# Line 544  template <typename booleanTest> Line 575  template <typename booleanTest>
575  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
576  __operationBoolean(ObjectToTreat& O0,  __operationBoolean(ObjectToTreat& O0,
577                     const ObjectToTreat& O1,                     const ObjectToTreat& O1,
578                     ObjectToTreat& O2)                     const ObjectToTreat& O2)
579  {  {
580    std::set<const Cell*> hexaList1;    std::set<const Cell*> hexaList1;
581    std::set<const Cell*> hexaList2;    std::set<const Cell*> hexaList2;
# Line 556  __operationBoolean(ObjectToTreat& O0, Line 587  __operationBoolean(ObjectToTreat& O0,
587      hexaList2.insert((*i).first);      hexaList2.insert((*i).first);
588    }    }
589    
590    toTreatHexahedra.clear();    std::set<const Cell*> toTreatHexahedra;
591    
592    std::set_intersection (hexaList1.begin(), hexaList1.end(),    std::set_intersection (hexaList1.begin(), hexaList1.end(),
593                           hexaList2.begin(), hexaList2.end(),                           hexaList2.begin(), hexaList2.end(),
# Line 564  __operationBoolean(ObjectToTreat& O0, Line 595  __operationBoolean(ObjectToTreat& O0,
595    
596    ffout(0)<<"nb hexa to treat "<<toTreatHexahedra.size()<<"\n";    ffout(0)<<"nb hexa to treat "<<toTreatHexahedra.size()<<"\n";
597        
598    __calculateIntersection<booleanTest>(O0,O1,O2);    __calculateIntersection<booleanTest>(O0,O1,O2, toTreatHexahedra);
599  }  }
600    
601  template <typename booleanTest>  template <typename booleanTest>
602  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
603  __calculateIntersection(ObjectToTreat& O0,  __calculateIntersection(ObjectToTreat& O0,
604                          const ObjectToTreat& O1,                          const ObjectToTreat& O1,
605                          const ObjectToTreat& O2)                          const ObjectToTreat& O2,
606                            const std::set<const Cell*>& toTreatHexahedra)
607  {  {
608    __calculatePointsIntersection(O1,O2);    PIntersection listVertexIntersectionNew;
609      std::vector<Triangle> triangleListIntersectionNew;
610    
611      __calculatePointsIntersection(O1, O2, toTreatHexahedra, listVertexIntersectionNew);
612    
613    //create mesh intersection for the object 0 and 1    //create mesh intersection for the object 0 and 1
614    __createTrianglesIntersection<booleanTest>(0,O1,O2);    __createTrianglesIntersection<booleanTest>(0, O1, O2, toTreatHexahedra,
615                                                 triangleListIntersectionNew, listVertexIntersectionNew);
616    ffout(0)<<"nb triangles in intersection 0 "<<triangleListIntersectionNew.size()<<"\n"<<std::flush;    ffout(0)<<"nb triangles in intersection 0 "<<triangleListIntersectionNew.size()<<"\n"<<std::flush;
617    __createTrianglesIntersection<booleanTest>(1,O2,O1);    __createTrianglesIntersection<booleanTest>(1, O2, O1, toTreatHexahedra,
618                                                 triangleListIntersectionNew, listVertexIntersectionNew);
619    ffout(0)<<"nb triangles in intersection 1 "<<triangleListIntersectionNew.size()<<"\n"<<std::flush;    ffout(0)<<"nb triangles in intersection 1 "<<triangleListIntersectionNew.size()<<"\n"<<std::flush;
620    
621      pairTriangleToIntersectionPoints.clear();
622        
623    //triangle from intersection of objects 0-1 in O0    //triangle from intersection of objects 0-1 in O0
624    size_t sizeinter=triangleListIntersectionNew.size();    size_t sizeinter=triangleListIntersectionNew.size();
# Line 591  __calculateIntersection(ObjectToTreat& O Line 630  __calculateIntersection(ObjectToTreat& O
630      (*O0.hexalist)[&(*O0.trianglelist)[size].mother()].push_back(&(*O0.trianglelist)[size]);      (*O0.hexalist)[&(*O0.trianglelist)[size].mother()].push_back(&(*O0.trianglelist)[size]);
631    }    }
632    triangleListIntersectionNew.clear();    triangleListIntersectionNew.clear();
633      listVertexIntersectionNew.clear();
634  }  }
635    
636  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
# Line 624  __createCutEdges(const Scene& scene, Line 664  __createCutEdges(const Scene& scene,
664  void SurfaceMeshGenerator::Internals::__constructionVerticesList(const ObjectToTreat& O0) {  void SurfaceMeshGenerator::Internals::__constructionVerticesList(const ObjectToTreat& O0) {
665    
666    const Vector<Triangle>& triangleList=(*O0.trianglelist);    const Vector<Triangle>& triangleList=(*O0.trianglelist);
667    //if (nobjectToTreat==0){  
668    for(size_t size=0 ; size< triangleList.size() ; ++size) {    for(size_t size=0 ; size< triangleList.size() ; ++size) {
669      const Triangle& T=triangleList[size];      const Triangle& T=triangleList[size];
670      if (listOfVertexMesh.find(&T(0))==listOfVertexMesh.end()) {      if (listOfVertexMesh.find(&T(0))==listOfVertexMesh.end()) {
671        listOfVertexMesh[&T(0)]=(&T(0));        listOfVertexMesh[&T(0)] = &T(0);
672      }      }
673      if (listOfVertexMesh.find(&T(1))==listOfVertexMesh.end()) {      if (listOfVertexMesh.find(&T(1))==listOfVertexMesh.end()) {
674        listOfVertexMesh[&T(1)]=(&T(1));        listOfVertexMesh[&T(1)] = &T(1);
675      }      }
676      if (listOfVertexMesh.find(&T(2))==listOfVertexMesh.end()) {      if (listOfVertexMesh.find(&T(2))==listOfVertexMesh.end()) {
677        listOfVertexMesh[&T(2)]=(&T(2));        listOfVertexMesh[&T(2)] = &T(2);
678      }      }
679    }    }
   //}  
   /* Vector<const Vertex*> listVertexFinal(listOfVertexMesh.size());  
      Vector<Vertex> listVertexFinal2(listOfVertexMesh.size());  
      size_t number=0;  
      for(std::map<const Vertex*,const Vertex*>::iterator it=listOfVertexMesh.begin() ;  
      it!=listOfVertexMesh.end() ; it++) {  
      listVertexFinal[number]=(*it).second;  
      listVertexFinal2[number]=*((*it).second);  
      ++number;  
      }*/  
680  }  }
681    
682  void SurfaceMeshGenerator::Internals::__constructionFinalMesh(const ObjectToTreat& O0,  void SurfaceMeshGenerator::Internals::__constructionFinalMesh(const ObjectToTreat& O0,
# Line 656  void SurfaceMeshGenerator::Internals::__ Line 686  void SurfaceMeshGenerator::Internals::__
686    
687    size_t ntr=0;    size_t ntr=0;
688        
   //ntr+=listVertexFinal.size();  
689    ntr+=listOfVertexMesh.size();    ntr+=listOfVertexMesh.size();
690        
691    s_mesh.setNumberOfVertices(ntr);    s_mesh.setNumberOfVertices(ntr);
# Line 669  void SurfaceMeshGenerator::Internals::__ Line 698  void SurfaceMeshGenerator::Internals::__
698       listVertexFinal[v] = &s_mesh.vertex(i);       listVertexFinal[v] = &s_mesh.vertex(i);
699       i++;       i++;
700       }*/       }*/
701    for(std::map<const Vertex*,const Vertex*>::iterator it=listOfVertexMesh.begin() ;  
702      std::set<const Vertex*> newVertices;
703      std::set<const Vertex*> oldVertices;
704    
705      for(std::map<const Vertex*, const Vertex*>::iterator it=listOfVertexMesh.begin() ;
706        it!=listOfVertexMesh.end() ; it++) {        it!=listOfVertexMesh.end() ; it++) {
707      s_mesh.vertex(i) = *((*it).second);      s_mesh.vertex(i) = *((*it).second);
708      delete (*it).second;  
709      (*it).second = &s_mesh.vertex(i);      delete ((*it).second);
710        ((*it).second) = &s_mesh.vertex(i);
711      i++;      i++;
712    }    }
713      
714    const Vector<Triangle>& triangleList=(*O0.trianglelist);    const Vector<Triangle>& triangleList=(*O0.trianglelist);
715    i = 0;    i = 0;
716    size_t ncell=0;    size_t ncell=0;
# Line 688  void SurfaceMeshGenerator::Internals::__ Line 722  void SurfaceMeshGenerator::Internals::__
722    s_mesh.setNumberOfCells(ncell);    s_mesh.setNumberOfCells(ncell);
723    ffout(0)<<"nb tot cells "<<s_mesh.numberOfCells()<<" \n"<<std::flush;    ffout(0)<<"nb tot cells "<<s_mesh.numberOfCells()<<" \n"<<std::flush;
724        
   /*   SurfaceMeshGenerator::Internals::MotherCellList::iterator motherCell  
        = cellList.begin();  
         
        for(size_t size=0 ; size< (*(*itcell)).size() ; ++size) {  
        const Triangle& T=(*(*itcell))[size];  
        Vertex P=T(0);  
        size_t n0=listVertexFinal2.number(P);  
        size_t n1=listVertexFinal2.number(T(1));  
        size_t n2=listVertexFinal2.number(T(2));  
        const Vertex& V0=*listVertexFinal[n0];  
        const Vertex& V1=*listVertexFinal[n1];  
        const Vertex& V2=*listVertexFinal[n2];  
        s_mesh.cell(i) = Triangle(V0,V1,V2);  
        #warning SHOULD USE BC NUMBER (ie OBJECT NUMBER IN SCENE INSTEAD)  
        s_mesh.cell(i).reference() = 1;  
        s_mesh.cell(i).setMother(&T.mother());  
        i++;  
        motherCell++;  
        }        */  
   //  if(nobjectToTreat==0){  
     
   //SurfaceMeshGenerator::Internals::MotherCellList::iterator motherCell  
     //= cellList.begin();  
725    for(size_t size=0 ; size< triangleList.size() ; ++size) {    for(size_t size=0 ; size< triangleList.size() ; ++size) {
726      const Triangle& T=triangleList[size];      const Triangle& T=triangleList[size];
727      std::map<const Vertex*,const Vertex*>::iterator      std::map<const Vertex*, const Vertex*>::iterator
728        it0=listOfVertexMesh.find(&T(0));        it0=listOfVertexMesh.find(&T(0));
729      const Vertex& V0=*((*it0).second);      const Vertex& V0=*((*it0).second);
730      std::map<const Vertex*,const Vertex*>::iterator      std::map<const Vertex*, const Vertex*>::iterator
731        it1=listOfVertexMesh.find(&T(1));        it1=listOfVertexMesh.find(&T(1));
732      const Vertex& V1=*((*it1).second);      const Vertex& V1=*((*it1).second);
733      std::map<const Vertex*,const Vertex*>::iterator      std::map<const Vertex*, const Vertex*>::iterator
734        it2=listOfVertexMesh.find(&T(2));        it2=listOfVertexMesh.find(&T(2));
735      const Vertex& V2=*((*it2).second);      const Vertex& V2=*((*it2).second);
736      s_mesh.cell(i) = Triangle(V0,V1,V2);      s_mesh.cell(i) = Triangle(V0,V1,V2, T.reference());
737  #warning SHOULD USE BC NUMBER (ie OBJECT NUMBER IN SCENE INSTEAD)  #warning SHOULD USE BC NUMBER (ie OBJECT NUMBER IN SCENE INSTEAD)
738      s_mesh.cell(i).reference() = 1;      //    s_mesh.cell(i).reference() = 1;
739      s_mesh.cell(i).setMother(&T.mother());      s_mesh.cell(i).setMother(&T.mother());
740      i++;      i++;
741    }    }
# Line 785  void SurfaceMeshGenerator::Internals::__ Line 796  void SurfaceMeshGenerator::Internals::__
796    
797  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
798  __calculatePointsIntersection(const ObjectToTreat& O1,  __calculatePointsIntersection(const ObjectToTreat& O1,
799                                const ObjectToTreat& O2)                                const ObjectToTreat& O2,
800                                  const std::set<const Cell*>& toTreatHexahedra,
801                                  PIntersection& listVertexIntersectionNew)
802  {  {
803        // calculate points intersection between object O1 and O2
804    MapCellTriangle hexaToTriangle0=*(O1.hexalist);    MapCellTriangle hexaToTriangle0=*(O1.hexalist);
805    MapCellTriangle hexaToTriangle1=*(O2.hexalist);    MapCellTriangle hexaToTriangle1=*(O2.hexalist);
806            
# Line 801  __calculatePointsIntersection(const Obje Line 815  __calculatePointsIntersection(const Obje
815            jt0!=listTriangle0.end() ; ++jt0) {            jt0!=listTriangle0.end() ; ++jt0) {
816          for(std::list<Triangle*>::iterator jt1=listTriangle1.begin() ;          for(std::list<Triangle*>::iterator jt1=listTriangle1.begin() ;
817              jt1!=listTriangle1.end() ; ++jt1) {              jt1!=listTriangle1.end() ; ++jt1) {
818            __createPointIntersection(0,1,jt0,jt1);            __createPointIntersection(jt0,jt1,listVertexIntersectionNew);
819          }          }
820        }        }
821      } else {      } else {
# Line 837  SurfaceMeshGenerator::Internals::__deter Line 851  SurfaceMeshGenerator::Internals::__deter
851  }  }
852    
853  void SurfaceMeshGenerator::Internals::__casePinterNotInList(const TriangleCut& K,  void SurfaceMeshGenerator::Internals::__casePinterNotInList(const TriangleCut& K,
854                                                              const Cell* cell){                                                              const Cell* cell,
855                                                                std::vector<Triangle>& triangleListIntersectionNew)
856    {
857        //build triangles if K is outside but the 3 egdes are cut
858    TinyVector<3,const Vertex*> triangle;    TinyVector<3,const Vertex*> triangle;
859    TinyVector<3,size_t> ncase;    TinyVector<3,size_t> ncase;
860    if(K.pointsIntersection[0][0]==K.pointsIntersection[0][1]) {    if(K.pointsIntersection[0][0]==K.pointsIntersection[0][1]) {
# Line 864  void SurfaceMeshGenerator::Internals::__ Line 881  void SurfaceMeshGenerator::Internals::__
881        }        }
882      }      }
883    }    }
884    __addTriangle(*triangle[0],*triangle[1],*triangle[2],cell);    __addTriangle(*triangle[0],*triangle[1],*triangle[2],cell, triangleListIntersectionNew);
885  }  }
886    
 // void SurfaceMeshGenerator::Internals::plotini(size_t nobjectToTreat) {  
 //   std::ofstream sortie2("ini.mesh");  
 //   sortie2<<"MeshVersionFormatted 1"<<" \n";  
 //   sortie2<<" \n";  
 //   sortie2<<"Dimension"<<" \n";  
 //   sortie2<<"3"<<" \n";  
 //   sortie2<<" \n";  
 //   sortie2<<"Vertices"<<" \n";  
 //   size_t nn=0,compt=0;  
 //   TinyVector<10,size_t> num;  
 //   for(size_t numObj=0 ; numObj<nobjectToTreat ; ++numObj) {  
 //     nn+=triangleListes.size();  
 //     if(compt<10) {  
 //       num[compt]=nn;  
 //     }  
 //     ++compt;  
 //   }  
 //   sortie2<<3*nn<<"\n";  
   
 //   for(size_t numObj=0 ; numObj<nobjectToTreat ; ++numObj) {  
 //     for (SurfaceMeshGenerator::Internals::EdgePairList::iterator  
 //         i = triangleListes.begin();  
 //       i != triangleListes.end(); ++i) {  
 //       SurfaceMeshGenerator::Internals::VerticesList& V = verticesListes;  
 //       const Vertex& V0=* V[(*i)[0]];  
 //       const Vertex& V1=* V[(*i)[1]];  
 //       const Vertex& V2=* V[(*i)[2]];  
 //       sortie2<<V0.Vector3()[0]<<" "<<V0.Vector3()[1]<<" "<<V0.Vector3()[2]<<" 0 \n";  
 //       sortie2<<V1.Vector3()[0]<<" "<<V1.Vector3()[1]<<" "<<V1.Vector3()[2]<<" 0 \n";  
 //       sortie2<<V2.Vector3()[0]<<" "<<V2.Vector3()[1]<<" "<<V2.Vector3()[2]<<" 0 \n";  
 //     }  
 //   }  
 //   sortie2<<" \n";  
 //   sortie2<<"Triangles"<<"\n";  
 //   sortie2<<nn<<"\n";  
   
 //   for(size_t nj=0;nj<nn;++nj) {  
 //     sortie2<<3*nj+1<<" "<<3*nj+2<<" "<<3*nj+3;  
 //     if(nj<num[0])  
 //       sortie2<<" 2 \n";  
 //     else if(nj<num[1])  
 //       sortie2<<" 3 \n";  
 //     else if(nj<num[2])  
 //       sortie2<<" 4 \n";  
 //     else if(nj<num[3])  
 //       sortie2<<" 5 \n";  
 //     else if(nj<num[4])  
 //       sortie2<<" 6 \n";  
 //     else  
 //       sortie2<<" 7 \n";  
 //   }  
 //   sortie2<<"End"<<" \n";  
 // }  
   
887  size_t SurfaceMeshGenerator::Internals::  size_t SurfaceMeshGenerator::Internals::
888  __IsInTriangle(const TinyVector<2,size_t> & num,const TinyVector<3,TinyVector<3> > & triangleVector )  __IsInTriangle(const TinyVector<2,size_t> & num,const TinyVector<3,TinyVector<3> > & triangleVector )
889  {  {
# Line 1004  __DeterminatePoint(const size_t& ncase, Line 967  __DeterminatePoint(const size_t& ncase,
967                     const Vertex* & Vo0,                     const Vertex* & Vo0,
968                     const Vertex* & Vo1,                     const Vertex* & Vo1,
969                     TinyVector<3,const Vertex* > & V,                     TinyVector<3,const Vertex* > & V,
970                     IntersectionPoints& Pinter) {                     IntersectionPoints& Pinter,
971                       PIntersection& listVertexIntersectionNew)
972    {
973        //calculate intersection between edge Vo0Vo1 and triangle V
974        // if ncase<3 the edge is in the first object else the edge is in the second object
975        // P = point intersection
976    std::pair<const Vertex*,const Vertex*> pV(Vo0,Vo1);    std::pair<const Vertex*,const Vertex*> pV(Vo0,Vo1);
977    std::pair<const Vertex*,const Vertex*> pV2(Vo1,Vo0);    std::pair<const Vertex*,const Vertex*> pV2(Vo1,Vo0);
978    
# Line 1045  __DeterminatePoint(const size_t& ncase, Line 1013  __DeterminatePoint(const size_t& ncase,
1013        if(ncase<3) {        if(ncase<3) {
1014          __InTriangle(P,V,num);          __InTriangle(P,V,num);
1015          if(num==3) {          if(num==3) {
1016            if(!(norm2>epsilon)) {  //        if(!(norm2>epsilon)) {
1017              Pinter.addPoint(Vo0);  //          Pinter.addPoint(Vo0);
1018            } else  //        } else
1019              if(!(norm3>epsilon)) {  //          if(!(norm3>epsilon)) {
1020                Pinter.addPoint(Vo1);  //            Pinter.addPoint(Vo1);
1021              } else  //          } else
1022                if(listVertexIntersectionNew.find(pV)==listVertexIntersectionNew.end()                if(listVertexIntersectionNew.find(pV)==listVertexIntersectionNew.end()
1023                   and listVertexIntersectionNew.find(pV2)==listVertexIntersectionNew.end()) {                   and listVertexIntersectionNew.find(pV2)==listVertexIntersectionNew.end()) {
1024                  listVertexIntersectionNew[pV][0]=new Vertex (P);                  listVertexIntersectionNew[pV][0]=new Vertex (P);
# Line 1062  __DeterminatePoint(const size_t& ncase, Line 1030  __DeterminatePoint(const size_t& ncase,
1030                            -P)<1e-6) {                            -P)<1e-6) {
1031                      Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[0]);                      Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[0]);
1032                    } else                    } else
1033                      if(Norm((*(*listVertexIntersectionNew.find(pV)).second[1])                      if(Norm((*(*listVertexIntersectionNew.find(pV)).second[1]) - P) < 1e-6) {
                             -P)<1e-6) {  
1034                        Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[1]);                        Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[1]);
1035                      } else {                      } else {
1036                        listVertexIntersectionNew[pV][1]=new Vertex (P);                        listVertexIntersectionNew[pV][1]=new Vertex (P);
1037                        Pinter.addPoint(listVertexIntersectionNew[pV][1]);                        Pinter.addPoint(listVertexIntersectionNew[pV][1]);
1038                      }                      }
1039                  } else {                  } else {
1040                    if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[0])                    if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[0]) - P) < 1e-6) {
                           -P)<1e-6) {  
1041                      Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[0]);                      Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[0]);
1042                    } else                    } else
1043                      if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[1])                      if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[1]) - P) < 1e-6) {
                             -P)<1e-6) {  
1044                        Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[1]);                        Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[1]);
1045                      } else {                      } else {
1046                        listVertexIntersectionNew[pV2][1]=new Vertex (P);                        listVertexIntersectionNew[pV2][1]=new Vertex (P);
# Line 1083  __DeterminatePoint(const size_t& ncase, Line 1048  __DeterminatePoint(const size_t& ncase,
1048                      }                      }
1049                  }                  }
1050                }                }
1051            Pinter.cutTriangle1[ncase]=Pinter.size();            Pinter.cutTriangle1()[ncase]=Pinter.size();
1052          }          }
1053        } else {        } else {
1054          if(inV01 > epsilon and inV02 > epsilon and inV12 > epsilon ) {          if(inV01 > epsilon and inV02 > epsilon and inV12 > epsilon ) {
# Line 1093  __DeterminatePoint(const size_t& ncase, Line 1058  __DeterminatePoint(const size_t& ncase,
1058              //ffout(0)<<P[0]<<" "<<P[1]<<" "<<P[2]<<" \n";              //ffout(0)<<P[0]<<" "<<P[1]<<" "<<P[2]<<" \n";
1059              //listVertexIntersectionNew.push_back(new Vertex (P));              //listVertexIntersectionNew.push_back(new Vertex (P));
1060              PIntersection PIn;              PIntersection PIn;
1061              if(!(norm2>epsilon)) {  //          if(!(norm2>epsilon)) {
1062                Pinter.addPoint(Vo0);  //            Pinter.addPoint(Vo0);
1063              } else {  //          } else {
1064                if(!(norm3>epsilon)) {  //            if(!(norm3>epsilon)) {
1065                  Pinter.addPoint(Vo1);  //              Pinter.addPoint(Vo1);
1066                } else {  //            } else
1067                    {
1068                  if(listVertexIntersectionNew.find(pV)==listVertexIntersectionNew.end()                  if(listVertexIntersectionNew.find(pV)==listVertexIntersectionNew.end()
1069                     and listVertexIntersectionNew.find(pV2)==listVertexIntersectionNew.end()) {                     and listVertexIntersectionNew.find(pV2)==listVertexIntersectionNew.end()) {
1070                    listVertexIntersectionNew[pV][0]=new Vertex (P);                    listVertexIntersectionNew[pV][0]=new Vertex (P);
# Line 1134  __DeterminatePoint(const size_t& ncase, Line 1100  __DeterminatePoint(const size_t& ncase,
1100                    }                    }
1101                  }                  }
1102                }                }
1103              }                //            }
1104              // Pinter.addPoint(listVertexIntersectionNew[listVertexIntersectionNew.size()-1]);              // Pinter.addPoint(listVertexIntersectionNew[listVertexIntersectionNew.size()-1]);
1105            }            }
1106          }          }
# Line 1144  __DeterminatePoint(const size_t& ncase, Line 1110  __DeterminatePoint(const size_t& ncase,
1110  }  }
1111    
1112  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1113  __CalculPoint(const size_t& objectNumber,const size_t& numobject,  __CalculPoint(std::list<Triangle*>::iterator & currentTriangle,
               std::list<Triangle*>::iterator & currentTriangle,  
1114                std::list<Triangle*>::iterator & currentTriangle1,                std::list<Triangle*>::iterator & currentTriangle1,
1115                IntersectionPoints& Pinter)                IntersectionPoints& Pinter,
1116                  PIntersection& listVertexIntersectionNew)
1117  {  {
1118    size_t num1=objectNumber,num0=numobject;      // calculate points intersection between triangle "currentTriangle" and "currentTriangle1"
1119        // this points are in "Pinter"
1120    TinyVector<3,const Vertex * > V012;    TinyVector<3,const Vertex * > V012;
1121    const Triangle& T0=*(*currentTriangle1);    const Triangle& T0=*(*currentTriangle1);
1122    V012[0] = &T0(0);    V012[0] = &T0(0);
# Line 1161  __CalculPoint(const size_t& objectNumber Line 1127  __CalculPoint(const size_t& objectNumber
1127    Vo012[0] = &To0(0);    Vo012[0] = &To0(0);
1128    Vo012[1] = &To0(1);    Vo012[1] = &To0(1);
1129    Vo012[2] = &To0(2);    Vo012[2] = &To0(2);
1130    //for each edge of the triangle V0V1V2  
1131    Pinter.numpair[0]=num0;    __DeterminatePoint(0,(V012[0]),(V012[1]),Vo012,Pinter, listVertexIntersectionNew);
1132    Pinter.numpair[1]=num1;    __DeterminatePoint(1,(V012[0]),(V012[2]),Vo012,Pinter, listVertexIntersectionNew);
1133      __DeterminatePoint(2,(V012[1]),(V012[2]),Vo012,Pinter, listVertexIntersectionNew);
1134    __DeterminatePoint(0,(V012[0]),(V012[1]),Vo012,Pinter);    
1135    __DeterminatePoint(1,(V012[0]),(V012[2]),Vo012,Pinter);    __DeterminatePoint(3,(Vo012[0]),(Vo012[1]),V012,Pinter, listVertexIntersectionNew);
1136    __DeterminatePoint(2,(V012[1]),(V012[2]),Vo012,Pinter);    __DeterminatePoint(3,(Vo012[0]),(Vo012[2]),V012,Pinter, listVertexIntersectionNew);
1137        __DeterminatePoint(3,(Vo012[1]),(Vo012[2]),V012,Pinter, listVertexIntersectionNew);
   __DeterminatePoint(3,(Vo012[0]),(Vo012[1]),V012,Pinter);  
   __DeterminatePoint(3,(Vo012[0]),(Vo012[2]),V012,Pinter);  
   __DeterminatePoint(3,(Vo012[1]),(Vo012[2]),V012,Pinter);  
1138  }  }
1139    
1140  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
# Line 1179  __inout(const size_t& numobject, Line 1142  __inout(const size_t& numobject,
1142          const int & edgeCut,          const int & edgeCut,
1143          const Vertex*& V0,          const Vertex*& V0,
1144          const Vertex*& V1)          const Vertex*& V1)
1145  {  {// if the edge is not cut , put ref=2 for the 2 vertex
1146    size_t& refVertex0=(*edgesRefVertex[numobject])[V0];    size_t& refVertex0=(*edgesRefVertex[numobject])[V0];
1147    size_t& refVertex1=(*edgesRefVertex[numobject])[V1];    size_t& refVertex1=(*edgesRefVertex[numobject])[V1];
1148    
# Line 1199  __inout(const size_t& numobject, Line 1162  __inout(const size_t& numobject,
1162  }  }
1163    
1164  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1165  __createPointIntersection(const size_t& numobject,const size_t& objectNumber,  __createPointIntersection(std::list<Triangle*>::iterator currentTriangle1,
1166                            std::list<Triangle*>::iterator currentTriangle1,                            std::list<Triangle*>::iterator currentTriangle,
1167                            std::list<Triangle*>::iterator currentTriangle)                            PIntersection& listVertexIntersectionNew)
1168  {      {    
1169    //ffout(0)<< "createpointintersection "<<numobject<<" "<<objectNumber<<"\n";    //ffout(0)<< "createpointintersection "<<numobject<<" "<<objectNumber<<"\n";
1170    //ffout(0)<<"*************************\n";    //ffout(0)<<"*************************\n";
# Line 1212  __createPointIntersection(const size_t& Line 1175  __createPointIntersection(const size_t&
1175    
1176    IntersectionPoints Pinter;    IntersectionPoints Pinter;
1177            
1178    __CalculPoint(objectNumber,numobject,currentTriangle,currentTriangle1,Pinter);    __CalculPoint(currentTriangle, currentTriangle1, Pinter, listVertexIntersectionNew);
1179    
1180    const TinyVector<3,size_t >& pintercut=Pinter.cutTriangle1;    const TinyVector<3,size_t >& pintercut=Pinter.cutTriangle1();
1181    
1182    if(pintercut[0]+pintercut[1]+pintercut[2]!=0) {    if(pintercut[0]+pintercut[1]+pintercut[2]!=0) {
1183      //ffout(0)<<"le tr est cut\n";      //ffout(0)<<"le tr est cut\n";
1184      if(pairTriangleToIntersectionPoints.find(pTriangle)==pairTriangleToIntersectionPoints.end()) {      assert(pairTriangleToIntersectionPoints.find(pTriangle)==pairTriangleToIntersectionPoints.end());
1185        //ffout(0) << "rajout de pTriangle\n";      pairTriangleToIntersectionPoints[pTriangle] = Pinter;
       pairTriangleToIntersectionPoints[pTriangle] = Pinter;  
     }  
1186    } else {    } else {
1187      if(Pinter.size()!=0 and numobject<objectNumber) {      if(Pinter.size()!=0) {
1188        //ffout(0)<<"rajout2 \n";        //ffout(0)<<"rajout2 \n";
1189        pairTriangleToIntersectionPoints[pTriangle] = Pinter;        pairTriangleToIntersectionPoints[pTriangle] = Pinter;
1190      }      }
1191    }    }
1192    
1193    if(pairTriangleToIntersectionPoints.find(pTriangle)!= pairTriangleToIntersectionPoints.end()    if(pairTriangleToIntersectionPoints.find(pTriangle)!= pairTriangleToIntersectionPoints.end()) {
1194       and numobject<objectNumber) {      TinyVector<3,size_t >& cutTr2= pairTriangleToIntersectionPoints[pTriangle].cutTriangle2();
     TinyVector<3,size_t >& cutTr2= pairTriangleToIntersectionPoints[pTriangle].cutTriangle2;  
1195    
1196      for(size_t i=0 ; i<Pinter.size() ; ++i) {      for(size_t i=0 ; i<Pinter.size() ; ++i) {
1197        real_t epsilon=1e-10;        real_t epsilon=1e-10;
# Line 1239  __createPointIntersection(const size_t& Line 1199  __createPointIntersection(const size_t&
1199        const Vertex& Vo0 = To(0);        const Vertex& Vo0 = To(0);
1200        const Vertex& Vo1 = To(1);        const Vertex& Vo1 = To(1);
1201        const Vertex& Vo2 = To(2);        const Vertex& Vo2 = To(2);
1202        const Vertex& P=*(pairTriangleToIntersectionPoints[pTriangle].listofPoints[i]);        const Vertex& P=*(pairTriangleToIntersectionPoints[pTriangle][i]);
1203                            
1204        real_t Vo01=Norm(Vo1-Vo0);        real_t Vo01=Norm(Vo1-Vo0);
1205        real_t Vo12=Norm(Vo2-Vo1);        real_t Vo12=Norm(Vo2-Vo1);
# Line 1266  __createPointIntersection(const size_t& Line 1226  __createPointIntersection(const size_t&
1226  const Vertex& SurfaceMeshGenerator::Internals::  const Vertex& SurfaceMeshGenerator::Internals::
1227  __determinateNumber(const int& place,  __determinateNumber(const int& place,
1228                      const TinyVector<3,TinyVector<2,const Vertex*> >& pointsint)                      const TinyVector<3,TinyVector<2,const Vertex*> >& pointsint)
1229  {  {//return the point intersection
1230    if(place<=2) {    if(place<=2) {
1231      return *pointsint[place][0];      return *pointsint[place][0];
1232    } else {    } else {
# Line 1278  void SurfaceMeshGenerator::Internals:: Line 1238  void SurfaceMeshGenerator::Internals::
1238  __addTriangle(const Vertex& V0,  __addTriangle(const Vertex& V0,
1239                const Vertex& V1,                const Vertex& V1,
1240                const Vertex& V2,                const Vertex& V2,
1241                const Cell* cell)                const Cell* cell,
1242                  std::vector<Triangle>& triangleListIntersectionNew)
1243  {  {
1244    Triangle T(V0,V1,V2);    Triangle T(V0,V1,V2);
1245    T.setMother(cell);    T.setMother(cell);
1246    if(!__isDegenerate(T)) {    if(!__isDegenerate(T)) {
1247      triangleListIntersectionNew.push_back(T);      triangleListIntersectionNew.push_back(T);
1248        T.reference() = 40;
1249      } else {
1250        fferr(0) << "\t\tTriangle dégénéré :"
1251                 << &V0 << ':' << &V1 << ':' << &V2 << '\n';
1252        T.reference() = 2;
1253        triangleListIntersectionNew.push_back(T);
1254    }    }
1255  }  }
1256    
# Line 1293  __determinateCase(const int& ncase, Line 1260  __determinateCase(const int& ncase,
1260                    TinyVector<3,int>& edgesCut,                    TinyVector<3,int>& edgesCut,
1261                    TinyVector<3,int>& vertexNum)                    TinyVector<3,int>& vertexNum)
1262  {  {
1263        //determinate which point is inside
1264    const TinyVector<3,TinyVector<2,int> >& edgesCut2=K.edgecut;    const TinyVector<3,TinyVector<2,int> >& edgesCut2=K.edgecut;
1265    switch(ncase) {        switch(ncase) {    
1266    case 0: {    case 0: {
# Line 1328  __determinateCase(const int& ncase, Line 1296  __determinateCase(const int& ncase,
1296  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1297  __createCase1(int ncase,  __createCase1(int ncase,
1298                const TriangleCut& K,                const TriangleCut& K,
1299                const Cell* cell)                const Cell* cell,
1300                  std::vector<Triangle>& triangleListIntersectionNew)
1301  {  {
1302        //case 1 points is inside
1303        //ncase=number of the point which is inside
1304    const Vertex& V0=K(0);    const Vertex& V0=K(0);
1305    const Vertex& V1=K(1);    const Vertex& V1=K(1);
1306    const Vertex& V2=K(2);    const Vertex& V2=K(2);
# Line 1350  __createCase1(int ncase, Line 1321  __createCase1(int ncase,
1321        const Vertex& B=__determinateNumber(vertexNum[0],pointsIntersection);        const Vertex& B=__determinateNumber(vertexNum[0],pointsIntersection);
1322        if(edgesCut[1]==1) {        if(edgesCut[1]==1) {
1323          const Vertex& C=__determinateNumber(vertexNum[1],pointsIntersection);          const Vertex& C=__determinateNumber(vertexNum[1],pointsIntersection);
1324          __addTriangle(A,B,C,cell);          __addTriangle(A,B,C,cell,triangleListIntersectionNew);
1325        } else {        } else {
1326          ffout(0) <<"pbs 10\n"<<std::flush;          ffout(0) <<"pbs 10\n"<<std::flush;
1327          ffout(0)<<V0<<"\n";          ffout(0)<<V0<<"\n";
# Line 1378  __createCase1(int ncase, Line 1349  __createCase1(int ncase,
1349        const Vertex& B=__determinateNumber(vertexNum[0],pointsIntersection);        const Vertex& B=__determinateNumber(vertexNum[0],pointsIntersection);
1350        if(edgesCut[1]==1) {        if(edgesCut[1]==1) {
1351          const Vertex& C=__determinateNumber(vertexNum[1],pointsIntersection);          const Vertex& C=__determinateNumber(vertexNum[1],pointsIntersection);
1352          __addTriangle(A,B,D,cell);          __addTriangle(A,B,D,cell,triangleListIntersectionNew);
1353          __addTriangle(A,D,C,cell);          __addTriangle(A,D,C,cell,triangleListIntersectionNew);
1354        } else {        } else {
1355          ffout(0)<<"-------- warning trou in 0\n"<<std::flush;          ffout(0)<<"-------- warning trou in 0\n"<<std::flush;
1356          ffout(0)<<V0<<"\n";          ffout(0)<<V0<<"\n";
1357          ffout(0)<<V1<<"\n";          ffout(0)<<V1<<"\n";
1358          ffout(0)<<V2<<"\n";          ffout(0)<<V2<<"\n";
1359          __addTriangle(A,B,D,cell);          __addTriangle(A,B,D,cell,triangleListIntersectionNew);
1360        }        }
1361      } else {      } else {
1362        if(edgesCut[1]==1) {        if(edgesCut[1]==1) {
1363          ffout(0)<<"-------- warning trou in\n"<<std::flush;          ffout(0)<<"-------- warning trou in\n"<<std::flush;
1364          const Vertex& B=__determinateNumber(vertexNum[1],pointsIntersection);          const Vertex& B=__determinateNumber(vertexNum[1],pointsIntersection);
1365          __addTriangle(A,B,D,cell);          __addTriangle(A,B,D,cell,triangleListIntersectionNew);
1366        } else {        } else {
1367          ffout(0) <<"pbs 1\n"<<std::flush;          ffout(0) <<"pbs 1\n"<<std::flush;
1368        }        }
# Line 1399  __createCase1(int ncase, Line 1370  __createCase1(int ncase,
1370      break;      break;
1371    }    }
1372    default: {    default: {
1373  #warning is it true?  #warning is it true? case where only one edge is cut?
     //a faire si plus d' 1 point d'intersection pas dans le triangle  
1374      ffout(0) << numberPointsInterior<<" be careful 1\n"<<std::flush;      ffout(0) << numberPointsInterior<<" be careful 1\n"<<std::flush;
1375      //POUR STEPH :  peut generer des tr plats        //case 1 point inside and more than 1 point in the triangle
1376       // ffout(0)<<V0<<"\n";      __createByEdges(K, cell, triangleListIntersectionNew);
      // ffout(0)<<V1<<"\n";  
      // ffout(0)<<V2<<"\n";  
     __createByEdges(K,cell);  
1377    }    }
1378    }    }
1379  }  }
# Line 1414  __createCase1(int ncase, Line 1381  __createCase1(int ncase,
1381  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1382  __createCase2(int ncase,  __createCase2(int ncase,
1383                const TriangleCut& K,                const TriangleCut& K,
1384                const Cell* cell)                const Cell* cell,
1385                  std::vector<Triangle>& triangleListIntersectionNew)
1386  {  {
1387        // case 2 points inside
1388    TinyVector<2,int> edgesCut;    TinyVector<2,int> edgesCut;
1389    TinyVector<2,int> vertexNum;    TinyVector<2,int> vertexNum;
1390    const Vertex& V0=K(0);    const Vertex& V0=K(0);
# Line 1464  __createCase2(int ncase, Line 1433  __createCase2(int ncase,
1433        const Vertex& C=__determinateNumber(vertexNum[0],pointsIntersection);        const Vertex& C=__determinateNumber(vertexNum[0],pointsIntersection);
1434        if(edgesCut[1]==1) {        if(edgesCut[1]==1) {
1435          const Vertex& D=__determinateNumber(vertexNum[1],pointsIntersection);          const Vertex& D=__determinateNumber(vertexNum[1],pointsIntersection);
1436          __addTriangle(A,B,C,cell);          __addTriangle(A,B,C,cell,triangleListIntersectionNew);
1437          __addTriangle(C,B,D,cell);          __addTriangle(C,B,D,cell,triangleListIntersectionNew);
1438        } else {        } else {
1439          // ffout(0) <<"pbs 20\n";          // ffout(0) <<"pbs 20\n";
1440          __addTriangle(A,B,C,cell);          __addTriangle(A,B,C,cell,triangleListIntersectionNew);
1441        }        }
1442      } else {      } else {
1443        if(edgesCut[1]==1) {        if(edgesCut[1]==1) {
1444          // ffout(0) <<"pbs 210\n";          // ffout(0) <<"pbs 210\n";
1445          const Vertex& C=__determinateNumber(vertexNum[1],pointsIntersection);          const Vertex& C=__determinateNumber(vertexNum[1],pointsIntersection);
1446          __addTriangle(A,B,C,cell);          __addTriangle(A,B,C,cell,triangleListIntersectionNew);
1447        } else {        } else {
1448          ffout(0) <<"pbs 21\n"<<std::flush;          ffout(0) <<"pbs 21\n"<<std::flush;
1449          ffout(0)<<V0<<"\n";          ffout(0)<<V0<<"\n";
# Line 1492  __createCase2(int ncase, Line 1461  __createCase2(int ncase,
1461        const Vertex& C=__determinateNumber(vertexNum[0],pointsIntersection);        const Vertex& C=__determinateNumber(vertexNum[0],pointsIntersection);
1462        if(edgesCut[1]==1) {        if(edgesCut[1]==1) {
1463          const Vertex& D=__determinateNumber(vertexNum[1],pointsIntersection);          const Vertex& D=__determinateNumber(vertexNum[1],pointsIntersection);
1464          __addTriangle(A,C,E,cell);          __addTriangle(A,C,E,cell,triangleListIntersectionNew);
1465          __addTriangle(A,E,B,cell);          __addTriangle(A,E,B,cell,triangleListIntersectionNew);
1466          __addTriangle(B,E,D,cell);          __addTriangle(B,E,D,cell,triangleListIntersectionNew);
1467        } else {        } else {
1468          __addTriangle(A,B,E,cell);          __addTriangle(A,B,E,cell,triangleListIntersectionNew);
1469          __addTriangle(C,E,B,cell);          __addTriangle(C,E,B,cell,triangleListIntersectionNew);
1470        }        }
1471      } else {      } else {
1472        if(edgesCut[1]==1) {        if(edgesCut[1]==1) {
1473          const Vertex& D=__determinateNumber(vertexNum[1],pointsIntersection);          const Vertex& D=__determinateNumber(vertexNum[1],pointsIntersection);
1474          __addTriangle(D,E,B,cell);          __addTriangle(D,E,B,cell,triangleListIntersectionNew);
1475          __addTriangle(A,B,E,cell);          __addTriangle(A,B,E,cell,triangleListIntersectionNew);
1476        } else {        } else {
1477          __addTriangle(A,B,E,cell);          __addTriangle(A,B,E,cell,triangleListIntersectionNew);
1478        }        }
1479      }      }
1480      break;      break;
# Line 1530  __createCase2(int ncase, Line 1499  __createCase2(int ncase,
1499    
1500      if(edgesCut[0]==1) {      if(edgesCut[0]==1) {
1501        const Vertex& C=__determinateNumber(vertexNum[0],pointsIntersection);        const Vertex& C=__determinateNumber(vertexNum[0],pointsIntersection);
1502        __addTriangle(A,B,C,cell);        __addTriangle(A,B,C,cell,triangleListIntersectionNew);
1503        P[0]=&C;        P[0]=&C;
1504        if(edgesCut[1]==1) {        if(edgesCut[1]==1) {
1505          const Vertex& D=__determinateNumber(vertexNum[1],pointsIntersection);          const Vertex& D=__determinateNumber(vertexNum[1],pointsIntersection);
# Line 1548  __createCase2(int ncase, Line 1517  __createCase2(int ncase,
1517                Ptemp[i-1]=P[i];                Ptemp[i-1]=P[i];
1518              }              }
1519            }            }
1520            __addTriangle(B,*P[0],*P[num],cell);            __addTriangle(B,*P[0],*P[num],cell,triangleListIntersectionNew);
1521            Ptemp[0]=P[num];            Ptemp[0]=P[num];
1522            P=Ptemp;            P=Ptemp;
1523            --size;            --size;
1524            // ffout(0)<<"size = "<<size<<"\n";            // ffout(0)<<"size = "<<size<<"\n";
1525          }          }
1526          __addTriangle(B,*P[0],*P[1],cell);          __addTriangle(B,*P[0],*P[1],cell,triangleListIntersectionNew);
1527          __addTriangle(B,*P[1],D,cell);          __addTriangle(B,*P[1],D,cell,triangleListIntersectionNew);
1528          //__addTriangle(A,B,D,cell);          //__addTriangle(A,B,D,cell);
1529        } else {        } else {
1530          ffout(0)<<"une seule arete cut pbs\n";          ffout(0)<<"une seule arete cut pbs\n";
# Line 1574  __createLocalListIntersection(const size Line 1543  __createLocalListIntersection(const size
1543                                TriangleCut & K,                                TriangleCut & K,
1544                                const Cell* cell)                                const Cell* cell)
1545  {  {
1546        //build the vectors K.pointsIntersection, K.pointsin with points intersection include in K
1547        //put K.edgecut[i]=1 if the edge i is cut
1548        //build vertexInTriangles, triangleWithVertex (for build mesh by front)
1549    TinyVector<3,TinyVector<2,int> >& edgesCut=K.edgecut;    TinyVector<3,TinyVector<2,int> >& edgesCut=K.edgecut;
1550    TinyVector<3,TinyVector<2,const Vertex*> > & pointsIntersection=K.pointsIntersection;    TinyVector<3,TinyVector<2,const Vertex*> > & pointsIntersection=K.pointsIntersection;
1551    std::map<const Vertex* ,const Vertex*>& pointsInterior=K.pointsin;    std::map<const Vertex* ,const Vertex*>& pointsInterior=K.pointsin;
# Line 1601  __createLocalListIntersection(const size Line 1573  __createLocalListIntersection(const size
1573        if(pairTriangleToIntersectionPoints.find(pairTriangle)!=pairTriangleToIntersectionPoints.end()) {        if(pairTriangleToIntersectionPoints.find(pairTriangle)!=pairTriangleToIntersectionPoints.end()) {
1574          TinyVector<3,size_t > cutTriangle;          TinyVector<3,size_t > cutTriangle;
1575          if(numobject<objectNumber) {          if(numobject<objectNumber) {
1576            cutTriangle=(pairTriangleToIntersectionPoints[pairTriangle]).cutTriangle1;            cutTriangle=(pairTriangleToIntersectionPoints[pairTriangle]).cutTriangle1();
1577          } else {          } else {
1578            cutTriangle=(pairTriangleToIntersectionPoints[pairTriangle]).cutTriangle2;            cutTriangle=(pairTriangleToIntersectionPoints[pairTriangle]).cutTriangle2();
1579          }          }
1580          //detection tr plat          //detection tr plat
1581          size_t isDegenerate=0;          size_t isDegenerate=0;
# Line 1612  __createLocalListIntersection(const size Line 1584  __createLocalListIntersection(const size
1584            isDegenerate=1;            isDegenerate=1;
1585          }          }
1586          if(pairTriangleToIntersectionPoints[pairTriangle].size()==2 and isDegenerate==0) {          if(pairTriangleToIntersectionPoints[pairTriangle].size()==2 and isDegenerate==0) {
1587            const Vertex*& V0=pairTriangleToIntersectionPoints[pairTriangle].listofPoints[0];            const Vertex*& V0=pairTriangleToIntersectionPoints[pairTriangle][0];
1588            const Vertex*& V1=pairTriangleToIntersectionPoints[pairTriangle].listofPoints[1];            const Vertex*& V1=pairTriangleToIntersectionPoints[pairTriangle][1];
1589            if((V0==&(*(*jt0))(0) or V0==&(*(*jt0))(1) or V0==&(*(*jt0))(2))            if((V0==&(*(*jt0))(0) or V0==&(*(*jt0))(1) or V0==&(*(*jt0))(2))
1590               and (V1==&(*(*jt0))(0) or V1==&(*(*jt0))(1) or V1==&(*(*jt0))(2))) {               and (V1==&(*(*jt0))(0) or V1==&(*(*jt0))(1) or V1==&(*(*jt0))(2))) {
1591              //ffout(0)<<"inter confondu avec arete\n";              //ffout(0)<<"inter confondu avec arete\n";
# Line 1622  __createLocalListIntersection(const size Line 1594  __createLocalListIntersection(const size
1594          }          }
1595          if(pairTriangleToIntersectionPoints[pairTriangle].size()>=2 and isDegenerate==0) {          if(pairTriangleToIntersectionPoints[pairTriangle].size()>=2 and isDegenerate==0) {
1596            for(size_t b=0; b<pairTriangleToIntersectionPoints[pairTriangle].size() ; ++b) {            for(size_t b=0; b<pairTriangleToIntersectionPoints[pairTriangle].size() ; ++b) {
1597              const Vertex*& V=pairTriangleToIntersectionPoints[pairTriangle].listofPoints[b];              const Vertex*& V=pairTriangleToIntersectionPoints[pairTriangle][b];
1598              if(!(V==&(*(*jt0))(0) or V==&(*(*jt0))(1) or V==&(*(*jt0))(2)) or b==0) {              if(!(V==&(*(*jt0))(0) or V==&(*(*jt0))(1) or V==&(*(*jt0))(2)) or b==0) {
1599                if(vertexInTriangles.find(V)==vertexInTriangles.end()) {                if(vertexInTriangles.find(V)==vertexInTriangles.end()) {
1600                  TinyVector<2,const Triangle*> triangles;                  TinyVector<2,const Triangle*> triangles;
# Line 1650  __createLocalListIntersection(const size Line 1622  __createLocalListIntersection(const size
1622              if(edgesCut[n][0]==0) {              if(edgesCut[n][0]==0) {
1623                edgesCut[n][0]=1;                edgesCut[n][0]=1;
1624                pointsIntersection[n][0]=                pointsIntersection[n][0]=
1625                  (pairTriangleToIntersectionPoints[pairTriangle].listofPoints[cutTriangle[n]-1]);                  (pairTriangleToIntersectionPoints[pairTriangle][cutTriangle[n]-1]);
1626              } else {              } else {
1627                edgesCut[n][1]=1;                edgesCut[n][1]=1;
1628                pointsIntersection[n][1]=                pointsIntersection[n][1]=
1629                  (pairTriangleToIntersectionPoints[pairTriangle].listofPoints[cutTriangle[n]-1]);                  (pairTriangleToIntersectionPoints[pairTriangle][cutTriangle[n]-1]);
1630              }              }
1631            }            }
1632          }          }
# Line 1671  __createLocalListIntersection(const size Line 1643  __createLocalListIntersection(const size
1643            case 0: {            case 0: {
1644              r=0;              r=0;
1645              for(size_t i=1 ; i<pairTriangleToIntersectionPoints[pairTriangle].size() ;++i)              for(size_t i=1 ; i<pairTriangleToIntersectionPoints[pairTriangle].size() ;++i)
1646                pointsInterior[(pairTriangleToIntersectionPoints[pairTriangle].listofPoints[i])]=                pointsInterior[(pairTriangleToIntersectionPoints[pairTriangle][i])]=
1647                  (pairTriangleToIntersectionPoints[pairTriangle].listofPoints[i]);                  (pairTriangleToIntersectionPoints[pairTriangle][i]);
1648              break;              break;
1649            }            }
1650            case 2:            case 2:
# Line 1697  __createLocalListIntersection(const size Line 1669  __createLocalListIntersection(const size
1669              ffout(0) << cutTriangle[0] << " " << cutTriangle[1] << " " << cutTriangle[2] << "\n";              ffout(0) << cutTriangle[0] << " " << cutTriangle[1] << " " << cutTriangle[2] << "\n";
1670            }            }
1671            }            }
1672            if(pointsInterior.find(pairTriangleToIntersectionPoints[pairTriangle].listofPoints[r])            if(pointsInterior.find(pairTriangleToIntersectionPoints[pairTriangle][r])
1673               ==pointsInterior.end()) {               ==pointsInterior.end()) {
1674              pointsInterior[(pairTriangleToIntersectionPoints[pairTriangle].listofPoints[r])]              pointsInterior[(pairTriangleToIntersectionPoints[pairTriangle][r])]
1675                =(pairTriangleToIntersectionPoints[pairTriangle].listofPoints[r]);                =(pairTriangleToIntersectionPoints[pairTriangle][r]);
1676            }            }
1677          }          }
1678        } else {        } else {
# Line 1715  __determinateEdgeCut(const size_t& numob Line 1687  __determinateEdgeCut(const size_t& numob
1687                       const ObjectToTreat& O1,                       const ObjectToTreat& O1,
1688                       const ObjectToTreat& O2,                       const ObjectToTreat& O2,
1689                       const Triangle*& currentTriangle,                       const Triangle*& currentTriangle,
1690                       const Cell* cell){                       const Cell* cell)
1691        {
1692      // put edgecut=1 if the edge is cut  
1693    TinyVector<3,TinyVector<2,int> > edgesCut;    TinyVector<3,TinyVector<2,int> > edgesCut;
1694    edgesCut=0;    edgesCut=0;
1695    size_t objectNumber=0;    size_t objectNumber=0;
# Line 1746  __determinateEdgeCut(const size_t& numob Line 1719  __determinateEdgeCut(const size_t& numob
1719    
1720          TinyVector<3,size_t> cutTriangle;          TinyVector<3,size_t> cutTriangle;
1721          if(numobject<objectNumber) {          if(numobject<objectNumber) {
1722            cutTriangle=(pairTriangleToIntersectionPoints[pairTriangle]).cutTriangle1;            cutTriangle=(pairTriangleToIntersectionPoints[pairTriangle]).cutTriangle1();
1723          } else {          } else {
1724            cutTriangle=(pairTriangleToIntersectionPoints[pairTriangle]).cutTriangle2;            cutTriangle=(pairTriangleToIntersectionPoints[pairTriangle]).cutTriangle2();
1725          }          }
1726    
1727          for (size_t n=0; n<3; ++n) {          for (size_t n=0; n<3; ++n) {
# Line 1770  __putEdgeRef(const size_t& ncase, Line 1743  __putEdgeRef(const size_t& ncase,
1743               const Triangle& currentTriangle,               const Triangle& currentTriangle,
1744               TriangleCut& K)               TriangleCut& K)
1745  {  {
1746        //put edgesRefVertex=2 if the vertex is in the domain
1747    const TinyVector<3,TinyVector<2,int> >&  cutTriangle=K.edgecut;    const TinyVector<3,TinyVector<2,int> >&  cutTriangle=K.edgecut;
1748    TinyVector<3,size_t>& isIn=K.isInTriangle;    TinyVector<3,size_t>& isIn=K.isInTriangle;
1749    size_t numobject=K.numobject;    size_t numobject=K.numobject;
# Line 1815  __putEdgeRef(const size_t& ncase, Line 1789  __putEdgeRef(const size_t& ncase,
1789    
1790  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1791  __createListTriangle(const TriangleCut& K,  __createListTriangle(const TriangleCut& K,
1792                       const Cell* cell)                       const Cell* cell,
1793                         std::vector<Triangle>& triangleListIntersectionNew)
1794  {  {
1795        //create triangles "include" in the triangle K
1796    const TinyVector<3,TinyVector<2,int> >&  edgesCut=K.edgecut;    const TinyVector<3,TinyVector<2,int> >&  edgesCut=K.edgecut;
1797    const TinyVector<3,size_t>& isIn=K.isInTriangle;    const TinyVector<3,size_t>& isIn=K.isInTriangle;
1798    const Vertex& V0=K(0);    const Vertex& V0=K(0);
# Line 1824  __createListTriangle(const TriangleCut& Line 1800  __createListTriangle(const TriangleCut&
1800    const Vertex& V2=K(2);    const Vertex& V2=K(2);
1801        
1802    if(((edgesCut[0][0]+edgesCut[1][0]+edgesCut[2][0]==0 and isIn[0]+isIn[1]+isIn[2]==3) /*or isIn[0]+isIn[1]+isIn[2]==3*/)) {    if(((edgesCut[0][0]+edgesCut[1][0]+edgesCut[2][0]==0 and isIn[0]+isIn[1]+isIn[2]==3) /*or isIn[0]+isIn[1]+isIn[2]==3*/)) {
1803      __addTriangle(V0,V1,V2,cell);      __addTriangle(V0,V1,V2,cell,triangleListIntersectionNew);
1804    } else {    } else {
1805      if(isIn[0]+isIn[1]+isIn[2]==1) {      if(isIn[0]+isIn[1]+isIn[2]==1) {
1806        if(isIn[0]==1) {        if(isIn[0]==1) {
1807          __createCase1(0,K,cell);          __createCase1(0,K,cell,triangleListIntersectionNew);
1808        }        }
1809        if(isIn[1]==1) {        if(isIn[1]==1) {
1810          __createCase1(1,K,cell);          __createCase1(1,K,cell,triangleListIntersectionNew);
1811        }        }
1812        if(isIn[2]==1) {        if(isIn[2]==1) {
1813          __createCase1(2,K,cell);          __createCase1(2,K,cell,triangleListIntersectionNew);
1814        }        }
1815      } else {      } else {
1816        if(isIn[0]+isIn[1]+isIn[2]==2) {        if(isIn[0]+isIn[1]+isIn[2]==2) {
1817          if(isIn[0]+isIn[1]==2) {          if(isIn[0]+isIn[1]==2) {
1818            __createCase2(0,K,cell);            __createCase2(0,K,cell,triangleListIntersectionNew);
1819          }          }
1820          if(isIn[0]+isIn[2]==2) {          if(isIn[0]+isIn[2]==2) {
1821            __createCase2(1,K,cell);            __createCase2(1,K,cell,triangleListIntersectionNew);
1822          }          }
1823          if(isIn[1]+isIn[2]==2) {          if(isIn[1]+isIn[2]==2) {
1824            __createCase2(2,K,cell);            __createCase2(2,K,cell,triangleListIntersectionNew);
1825          }          }
1826        }        }
1827      }      }
# Line 1853  __createListTriangle(const TriangleCut& Line 1829  __createListTriangle(const TriangleCut&
1829  }  }
1830    
1831  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1832  __createNotIn(const TriangleCut& K,const Cell* cell)  __createNotIn(const TriangleCut& K,const Cell* cell,
1833                  std::vector<Triangle>& triangleListIntersectionNew)
1834  {  {
1835        //create triangles if only one edge is cut
1836        // more than 3 points intersection?
1837    size_t place0=0;    size_t place0=0;
1838    if(K.edgecut[1][0]==1 and K.edgecut[1][1]==1) {    if(K.edgecut[1][0]==1 and K.edgecut[1][1]==1) {
1839      place0=1;      place0=1;
# Line 1863  __createNotIn(const TriangleCut& K,const Line 1842  __createNotIn(const TriangleCut& K,const
1842    }    }
1843    
1844    const Vertex* Pend=K.pointsIntersection[place0][1];    const Vertex* Pend=K.pointsIntersection[place0][1];
1845      #warning size TinyVector
1846    TinyVector<10,const Vertex*> P;    TinyVector<10,const Vertex*> P;
1847    TinyVector<10,const Triangle*> T;    TinyVector<10,const Triangle*> T;
1848    P[0]=K.pointsIntersection[place0][0];    P[0]=K.pointsIntersection[place0][0];
# Line 1880  __createNotIn(const TriangleCut& K,const Line 1860  __createNotIn(const TriangleCut& K,const
1860        TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[b])).second;        TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[b])).second;
1861        if(VV[0]==P[b]) {        if(VV[0]==P[b]) {
1862          if(Pend!=VV[1])          if(Pend!=VV[1])
1863            __addTriangle(*Pend,*P[b],*VV[1],cell);            __addTriangle(*Pend,*P[b],*VV[1],cell,triangleListIntersectionNew);
1864          P[b+1]=VV[1];          P[b+1]=VV[1];
1865        } else {        } else {
1866          if(Pend!=VV[0])          if(Pend!=VV[0])
1867            __addTriangle(*Pend,*P[b],*VV[0],cell);            __addTriangle(*Pend,*P[b],*VV[0],cell,triangleListIntersectionNew);
1868          P[b+1]=VV[0];          P[b+1]=VV[0];
1869        }        }
1870      } else {      } else {
# Line 1909  __createNotIn(const TriangleCut& K,const Line 1889  __createNotIn(const TriangleCut& K,const
1889  }  }
1890    
1891  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1892  __createNotIn2(const TriangleCut& K,const Cell* cell)  __createNotIn2(const TriangleCut& K,const Cell* cell,
1893  {                 std::vector<Triangle>& triangleListIntersectionNew)
1894    {// create triangles if 2 edges is cut
1895        // and more than 3 points intersection?
1896    size_t place0=0;    size_t place0=0;
1897    if(K.edgecut[1][0]==1 and K.edgecut[1][1]==1) {    if(K.edgecut[1][0]==1 and K.edgecut[1][1]==1) {
1898      place0=1;      place0=1;
# Line 1937  __createNotIn2(const TriangleCut& K,cons Line 1919  __createNotIn2(const TriangleCut& K,cons
1919        TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[b])).second;        TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[b])).second;
1920        if(VV[0]==P[b]) {        if(VV[0]==P[b]) {
1921          if(Pend!=VV[1]) {          if(Pend!=VV[1]) {
1922            __addTriangle(*Pend,*P[b],*VV[1],cell);            __addTriangle(*Pend,*P[b],*VV[1],cell,triangleListIntersectionNew);
1923            //  ffout(0)<<"1 "<<*Pend<<"\n";            //  ffout(0)<<"1 "<<*Pend<<"\n";
1924            //  ffout(0)<<*P[b]<<"\n";            //  ffout(0)<<*P[b]<<"\n";
1925            //  ffout(0)<<*VV[1]<<"\n";            //  ffout(0)<<*VV[1]<<"\n";
# Line 1945  __createNotIn2(const TriangleCut& K,cons Line 1927  __createNotIn2(const TriangleCut& K,cons
1927          P[b+1]=VV[1];          P[b+1]=VV[1];
1928        } else {        } else {
1929          if(Pend!=VV[0]) {          if(Pend!=VV[0]) {
1930            __addTriangle(*Pend,*P[b],*VV[0],cell);            __addTriangle(*Pend,*P[b],*VV[0],cell,triangleListIntersectionNew);
1931            //ffout(0)<<"2 "<<*Pend<<"\n";            //ffout(0)<<"2 "<<*Pend<<"\n";
1932            //ffout(0)<<*P[b]<<"\n";            //ffout(0)<<*P[b]<<"\n";
1933            //ffout(0)<<*VV[0]<<"\n";            //ffout(0)<<*VV[0]<<"\n";
# Line 1982  __createNotIn2(const TriangleCut& K,cons Line 1964  __createNotIn2(const TriangleCut& K,cons
1964          if(K.edgecut[0][0]==1 and K.edgecut[0][1]==1) {          if(K.edgecut[0][0]==1 and K.edgecut[0][1]==1) {
1965            //ffout(0)<<"here2\n";            //ffout(0)<<"here2\n";
1966            if(P[b+1]==K.pointsIntersection[0][1]) {            if(P[b+1]==K.pointsIntersection[0][1]) {
1967              __addTriangle(*Pend,*P[b+1],*K.pointsIntersection[0][0],cell);              __addTriangle(*Pend,*P[b+1],*K.pointsIntersection[0][0],cell,triangleListIntersectionNew);
1968              P[b+1]=K.pointsIntersection[0][0];                                P[b+1]=K.pointsIntersection[0][0];                  
1969              if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {              if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {
1970                T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];                T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];
# Line 1991  __createNotIn2(const TriangleCut& K,cons Line 1973  __createNotIn2(const TriangleCut& K,cons
1973              T[b+1]=T2[0];              T[b+1]=T2[0];
1974            } else            } else
1975              if(P[b+1]==K.pointsIntersection[0][0]){              if(P[b+1]==K.pointsIntersection[0][0]){
1976                __addTriangle(*Pend,*P[b+1],*K.pointsIntersection[0][1],cell);                __addTriangle(*Pend,*P[b+1],*K.pointsIntersection[0][1],cell,triangleListIntersectionNew);
1977                P[b+1]=K.pointsIntersection[0][1];                P[b+1]=K.pointsIntersection[0][1];
1978                if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {                if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {
1979                  T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];                  T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];
# Line 2011  __createNotIn2(const TriangleCut& K,cons Line 1993  __createNotIn2(const TriangleCut& K,cons
1993            if(K.edgecut[2][0]==1 and K.edgecut[2][1]==1) {            if(K.edgecut[2][0]==1 and K.edgecut[2][1]==1) {
1994              //ffout(0)<<"here3\n";              //ffout(0)<<"here3\n";
1995              if(P[b+1]==K.pointsIntersection[2][1]) {              if(P[b+1]==K.pointsIntersection[2][1]) {
1996                __addTriangle(*Pend,*P[b+1],*K.pointsIntersection[2][0],cell);                __addTriangle(*Pend,*P[b+1],*K.pointsIntersection[2][0],cell,triangleListIntersectionNew);
1997                P[b+1]=K.pointsIntersection[2][0];                P[b+1]=K.pointsIntersection[2][0];
1998                if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {                if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {
1999                  T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];                  T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];
# Line 2020  __createNotIn2(const TriangleCut& K,cons Line 2002  __createNotIn2(const TriangleCut& K,cons
2002                T[b+1]=T2[0];                T[b+1]=T2[0];
2003              } else {              } else {
2004                if(P[b+1]==K.pointsIntersection[2][0]){                if(P[b+1]==K.pointsIntersection[2][0]){
2005                __addTriangle(*Pend,*P[b+1],*K.pointsIntersection[2][1],cell);                __addTriangle(*Pend,*P[b+1],*K.pointsIntersection[2][1],cell,triangleListIntersectionNew);
2006                P[b+1]=K.pointsIntersection[2][1];                P[b+1]=K.pointsIntersection[2][1];
2007                if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {                if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {
2008                  T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];                  T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];
# Line 2044  __createNotIn2(const TriangleCut& K,cons Line 2026  __createNotIn2(const TriangleCut& K,cons
2026  }  }
2027    
2028  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
2029  __createByEdges(const TriangleCut& K,const Cell* cell)  __createByEdges(const TriangleCut& K,const Cell* cell,
2030                    std::vector<Triangle>& triangleListIntersectionNew)
2031  {  {
2032        //createtriangles if one point is in Triangle
2033        // and two edges are cuts
2034    size_t in0=0,place0=0,place1=1,place2=2;    size_t in0=0,place0=0,place1=1,place2=2;
2035    
2036    if(K.isInTriangle[1]==1) {    if(K.isInTriangle[1]==1) {
# Line 2062  __createByEdges(const TriangleCut& K,con Line 2047  __createByEdges(const TriangleCut& K,con
2047    }    }
2048    
2049    const Vertex* Pend=K.pointsIntersection[place1][0];    const Vertex* Pend=K.pointsIntersection[place1][0];
2050      #warning size of Tiny
2051    TinyVector<10,const Vertex*> P;    TinyVector<10,const Vertex*> P;
2052    TinyVector<10,const Triangle*> T;    TinyVector<10,const Triangle*> T;
2053    P[0]=K.pointsIntersection[place0][0];    P[0]=K.pointsIntersection[place0][0];
# Line 2080  __createByEdges(const TriangleCut& K,con Line 2066  __createByEdges(const TriangleCut& K,con
2066            ffout(0)<<"egal\n";            ffout(0)<<"egal\n";
2067            if(K.pointsIntersection[place2][0]==P[b]) {            if(K.pointsIntersection[place2][0]==P[b]) {
2068              P[b+1]=K.pointsIntersection[place2][1];              P[b+1]=K.pointsIntersection[place2][1];
2069              __addTriangle(K(in0),*P[b],*P[b+1],cell);              __addTriangle(K(in0),*P[b],*P[b+1],cell,triangleListIntersectionNew);
2070            } else {            } else {
2071              if(K.pointsIntersection[place2][1]==P[b]){              if(K.pointsIntersection[place2][1]==P[b]){
2072                P[b+1]=K.pointsIntersection[place2][0];                P[b+1]=K.pointsIntersection[place2][0];
2073                __addTriangle(K(in0),*P[b],*P[b+1],cell);                __addTriangle(K(in0),*P[b],*P[b+1],cell,triangleListIntersectionNew);
2074              } else {              } else {
2075                ffout(0)<<"pbs pbs\n";                ffout(0)<<"pbs pbs\n";
2076                ffout(0)<<*P[b]<<"\n";                ffout(0)<<*P[b]<<"\n";
# Line 2092  __createByEdges(const TriangleCut& K,con Line 2078  __createByEdges(const TriangleCut& K,con
2078              }              }
2079            }            }
2080          } else {          } else {
2081            __addTriangle(K(in0),*P[b],*VV[1],cell);            __addTriangle(K(in0),*P[b],*VV[1],cell,triangleListIntersectionNew);
2082            P[b+1]=VV[1];            P[b+1]=VV[1];
2083          }          }
2084        } else {        } else {
2085          __addTriangle(K(in0),*P[b],*VV[0],cell);          __addTriangle(K(in0),*P[b],*VV[0],cell,triangleListIntersectionNew);
2086          P[b+1]=VV[0];          P[b+1]=VV[0];
2087        }        }
2088      } else {      } else {
# Line 2117  __createByEdges(const TriangleCut& K,con Line 2103  __createByEdges(const TriangleCut& K,con
2103        if(P[b+1]!=Pend) {        if(P[b+1]!=Pend) {
2104          if(K.pointsIntersection[place2][0]==P[b+1]) {          if(K.pointsIntersection[place2][0]==P[b+1]) {
2105            const Vertex* A=K.pointsIntersection[place2][1];            const Vertex* A=K.pointsIntersection[place2][1];
2106            __addTriangle(K(in0),*P[b+1],*A,cell);            __addTriangle(K(in0),*P[b+1],*A,cell,triangleListIntersectionNew);
2107            if(vertexInTriangles.find(A)!=vertexInTriangles.end()) {            if(vertexInTriangles.find(A)!=vertexInTriangles.end()) {
2108              T[b+1]=(*vertexInTriangles.find(A)).second[0];              T[b+1]=(*vertexInTriangles.find(A)).second[0];
2109              P[b+1]=A;              P[b+1]=A;
2110            }            }
2111          } else {          } else {
2112            const Vertex* A=K.pointsIntersection[place2][0];            const Vertex* A=K.pointsIntersection[place2][0];
2113            __addTriangle(K(in0),*P[b+1],*A,cell);            __addTriangle(K(in0),*P[b+1],*A,cell,triangleListIntersectionNew);
2114            if(vertexInTriangles.find(A)!=vertexInTriangles.end()) {            if(vertexInTriangles.find(A)!=vertexInTriangles.end()) {
2115              T[b+1]=(*vertexInTriangles.find(A)).second[0];              T[b+1]=(*vertexInTriangles.find(A)).second[0];
2116              P[b+1]=A;              P[b+1]=A;
# Line 2147  template <typename booleanTest> Line 2133  template <typename booleanTest>
2133  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
2134  __createTrianglesIntersection(const size_t numobject,  __createTrianglesIntersection(const size_t numobject,
2135                                const ObjectToTreat& O1,                                const ObjectToTreat& O1,
2136                                const ObjectToTreat& O2)                                const ObjectToTreat& O2,
2137                                  const std::set<const Cell*>& toTreatHexahedra,
2138                                  std::vector<Triangle>& triangleListIntersectionNew,
2139                                  PIntersection& listVertexIntersectionNew)
2140  {  {
2141        //build mesh of object O1 (which is intersected with O2)
2142        
2143    //ffout(0)<<"numero "<<numobject<<" "<<objectNumber<<"\n"<<std::flush;    //ffout(0)<<"numero "<<numobject<<" "<<objectNumber<<"\n"<<std::flush;
2144    size_t objectNumber=0;    size_t objectNumber=0;
2145    if(numobject==0)    if(numobject==0)
# Line 2196  __createTrianglesIntersection(const size Line 2187  __createTrianglesIntersection(const size
2187          if((*edgesRefVertex[numobject])[&T(0)]==2          if((*edgesRefVertex[numobject])[&T(0)]==2
2188             and (*edgesRefVertex[numobject])[&T(1)]==2             and (*edgesRefVertex[numobject])[&T(1)]==2
2189             and (*edgesRefVertex[numobject])[&T(2)]==2) {             and (*edgesRefVertex[numobject])[&T(2)]==2) {
2190            __addTriangle(T(0),T(1),T(2),&T.mother());            __addTriangle(T(0),T(1),T(2),&T.mother(),triangleListIntersectionNew);
2191          }          }
2192        }        }
2193      }      }
# Line 2224  __createTrianglesIntersection(const size Line 2215  __createTrianglesIntersection(const size
2215    
2216              K.edgecut =__determinateEdgeCut(numobject,O1,O2,T,*i);              K.edgecut =__determinateEdgeCut(numobject,O1,O2,T,*i);
2217    
             // CECILE: Je commente cette affectation qui ne sert à rien (selon moi)  
             //  K.isInTriangle =  
2218              __putEdgeRef(0,(*T),K);              __putEdgeRef(0,(*T),K);
2219            }            }
2220          }          }
2221        }        }
2222    }    }
2223    
   size_t err=0,err2=0;  
2224    for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin();    for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin();
2225        i!=toTreatHexahedra.end() ; ++i) {        i!=toTreatHexahedra.end() ; ++i) {
2226      MapCellTriangle::iterator it0=hexaToTriangle0.find(*i);      MapCellTriangle::iterator it0=hexaToTriangle0.find(*i);
# Line 2251  __createTrianglesIntersection(const size Line 2239  __createTrianglesIntersection(const size
2239    
2240          __createLocalListIntersection(objectNumber,O2,T,K,*i);          __createLocalListIntersection(objectNumber,O2,T,K,*i);
2241    
         // CECILE: Je commente cette affectation qui ne sert à rien (selon moi)  
         //  K.isInTriangle =  
2242          __putEdgeRef(1,*T,K);          __putEdgeRef(1,*T,K);
2243    
2244          if(K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==0 ) {          // K.isInTriangle[i] vaut 1 si le point est a garder pour construire le maillage 0 sinon.
2245            ++err;          // K.edgecut[i] vaut 1 si l'arete est coupee
         }  
         if(K.edgecut[0]+K.edgecut[1]+K.edgecut[2]==0 ) {  
           ++err2;  
         }  
         /*size_t test=0;  
         TinyVector<3> A(0.0483333, 0.232544, -0.0263158);  
           TinyVector<3> B(-0.0263158, 0.319152, -0.187573);  
           TinyVector<3> C(0.0315867, 0.179403, -0.0627097);  
           real_t ee=1e-5;  
           if((Norm(V0-A)<=ee or Norm(V1-A)<=ee or Norm(V2-A)<=ee)  
           and (Norm(V0-B)<=ee or Norm(V1-B)<=ee or Norm(V2-B)<=ee)  
           and (Norm(V0-C)<=ee or Norm(V1-C)<=ee or Norm(V2-C)<=ee)) {  
           ffout(0)<<"*****************\n";  
           ffout(0)<<V0<<"\n";  
           ffout(0)<<V1<<"\n";  
           ffout(0)<<V2<<"\n";  
           ffout(0)<<K.edgecut<<"\n";  
           ffout(0)<<K.isInTriangle<<"\n";  
           ffout(0)<<vertexInTriangles.size()<<"\n";  
           for(std::map<const Vertex*,TinyVector<2,const Triangle*> >::iterator  
           it=vertexInTriangles.begin() ; it!=vertexInTriangles.end() ; ++it) {  
               ffout(0)<<*(*it).first<<"\n";  
           }  
           test=1;  
           }*/  
2246          if(K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==0          if(K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==0
2247             and K.edgecut[0]+K.edgecut[1]+K.edgecut[2]==3 and vertexInTriangles.size()<3) {             and K.edgecut[0]+K.edgecut[1]+K.edgecut[2]==3 and vertexInTriangles.size()<3) {
2248                // K is outside but the 3 edges are cut "and there are less than 3 points in K"
2249    
2250              // On ne peut normalement pas se trouver ICI.
2251            ffout(0)<<"tous les points d'inter ne sont pas dans la liste\n";            ffout(0)<<"tous les points d'inter ne sont pas dans la liste\n";
2252            //ffout(0)<<V0<<"\n";            __casePinterNotInList(K,*i, triangleListIntersectionNew);
2253            //ffout(0)<<V1<<"\n";            std::exit(1);
           //ffout(0)<<V2<<"\n";  
           //ffout(0)<<K.edgecut<<"\n";  
           //ffout(0)<<K.isInTriangle<<"\n";  
           __casePinterNotInList(K,*i);  
2254          }          }
2255    
2256          if(K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==0          if(K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==0
2257             and K.edgecut[0]+K.edgecut[1]+K.edgecut[2]!=0) {             and K.edgecut[0]+K.edgecut[1]+K.edgecut[2]!=0) {
2258                //K is outside
2259            if(vertexInTriangles.size()==3) {            if(vertexInTriangles.size()==3) {
2260                  //there are 3 points intersection
2261              TinyVector<3,const Vertex*> V;              TinyVector<3,const Vertex*> V;
2262              size_t num=0;              size_t num=0;
2263              for(std::map<const Vertex*,TinyVector<2,const Triangle*> >::iterator              for(std::map<const Vertex*,TinyVector<2,const Triangle*> >::iterator
# Line 2303  __createTrianglesIntersection(const size Line 2265  __createTrianglesIntersection(const size
2265                V[num]=(*it).first;                V[num]=(*it).first;
2266              }              }
2267              if(V[0]!=V[1] and V[1]!=V[2] and V[2]!=V[0]) {              if(V[0]!=V[1] and V[1]!=V[2] and V[2]!=V[0]) {
2268                __addTriangle(*V[0],*V[1],*V[2],*i);                __addTriangle(*V[0],*V[1],*V[2],*i,triangleListIntersectionNew);
2269                ffout(0)<<"on comble le trou\n"<<std::flush;                ffout(0)<<"on comble le trou\n"<<std::flush;
               //ffout(0)<<V0<<"\n";  
               //ffout(0)<<V1<<"\n";  
               //ffout(0)<<V2<<"\n";  
2270              } else {              } else {
2271                ffout(0)<<"----- deux inter =\n";                ffout(0)<<"----- deux inter =\n";
2272                ffout(0)<<V[0]<<"\n";                ffout(0)<<V[0]<<"\n";
# Line 2320  __createTrianglesIntersection(const size Line 2279  __createTrianglesIntersection(const size
2279                //one edge cut but twice                //one edge cut but twice
2280                ffout(0)<<"createNotIn\n";                ffout(0)<<"createNotIn\n";
2281                //POUR STEPH : createNotIn peut generer des tr plats                //POUR STEPH : createNotIn peut generer des tr plats
2282                __createNotIn(K,*i);                __createNotIn(K,*i, triangleListIntersectionNew);
2283              } else {              } else {
2284                if(vertexInTriangles.size()>3 and                if(vertexInTriangles.size()>3 and
2285                   K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==2) {                   K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==2) {
2286                //two edge cut but twice                //two edge cut but twice
2287                ffout(0)<<"createNotIn2\n";                ffout(0)<<"createNotIn2\n";
2288                __createNotIn2(K,*i);                __createNotIn2(K,*i, triangleListIntersectionNew);
2289                } else {                } else {
2290                  if(K.edgecut[0]!=0) {                  if(K.edgecut[0]!=0) {
2291                    if(K.pointsIntersection[0][0]==&V0                    if(K.pointsIntersection[0][0]==&V0
# Line 2351  __createTrianglesIntersection(const size Line 2310  __createTrianglesIntersection(const size
2310                        if(vertexInTriangles.size()>3 and                        if(vertexInTriangles.size()>3 and
2311                           K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==2){                           K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==2){
2312                          ffout(0)<<"hhhhhhhh\n";                          ffout(0)<<"hhhhhhhh\n";
2313                          __createNotIn2(K,*i);                          __createNotIn2(K,*i,triangleListIntersectionNew);
2314                        } else {                        } else {
2315                          if(vertexInTriangles.size()>3 and                          if(vertexInTriangles.size()>3 and
2316                             K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==1) {                             K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==1) {
2317                            ffout(0)<<"case\n";                            ffout(0)<<"case\n";
2318                            __createNotIn(K,*i);                            __createNotIn(K,*i,triangleListIntersectionNew);
2319                          }                          }
2320                        }                        }
2321                      } else {                      } else {
# Line 2408  __createTrianglesIntersection(const size Line 2367  __createTrianglesIntersection(const size
2367               and K.pointsIntersection[1][0]!=K.pointsIntersection[2][0]               and K.pointsIntersection[1][0]!=K.pointsIntersection[2][0]
2368               and K.pointsIntersection[2][0]!=K.pointsIntersection[0][0]) {               and K.pointsIntersection[2][0]!=K.pointsIntersection[0][0]) {
2369              ffout(0)<<"---------------- 3 edge cut \n";              ffout(0)<<"---------------- 3 edge cut \n";
2370              //ffout(0)<<V0<<"\n";              __createByEdges(K,*i,triangleListIntersectionNew);
             //ffout(0)<<V1<<"\n";  
             //ffout(0)<<V2<<"\n";  
             //ffout(0)<<"createByEdges\n";  
             __createByEdges(K,*i);  
2371            } else {            } else {
2372              __createListTriangle(K,*i);              __createListTriangle(K,*i,triangleListIntersectionNew);
2373            }            }
2374          } else {          } else {
2375            if(nbinters!=0) {            if(nbinters!=0) {
2376              __createListTriangle(K,*i);              __createListTriangle(K,*i,triangleListIntersectionNew);
2377            } else {            } else {
2378              //ffout(0)<<"j ajoute le tr car pas d'intersection dans les 2 listes\n";              //ffout(0)<<"j ajoute le tr car pas d'intersection dans les 2 listes\n";
2379              __addTriangle(V0,V1,V2,*i);              __addTriangle(V0,V1,V2,*i,triangleListIntersectionNew);
2380            }            }
2381          }          }
2382        }        }
# Line 2437  __createTriangleSurface(VerticesList& ve Line 2392  __createTriangleSurface(VerticesList& ve
2392                          EdgesRef& edgesRef,                          EdgesRef& edgesRef,
2393                          const Shape& S ,                          const Shape& S ,
2394                          std::vector<const Edge*>& cutEdges,                          std::vector<const Edge*>& cutEdges,
2395                          const Cell& currentCell)                          const Cell& currentCell,
2396                            MotherCellList& cellList)
2397  {  {
2398    bool jobDone=true;    bool jobDone=true;
2399    switch (cutEdges.size()) {    switch (cutEdges.size()) {
# Line 2590  __createSurface(Structured3DMesh& SMesh, Line 2546  __createSurface(Structured3DMesh& SMesh,
2546      std::exit(1);      std::exit(1);
2547    }    }
2548    
2549    for (size_t icell=0; icell<SMesh.shape().nx()-1; ++icell)    for (size_t icell=0; icell<SMesh.shape().nx()-1; ++icell) {
2550      for (size_t jcell=0; jcell<SMesh.shape().ny()-1; ++jcell)      for (size_t jcell=0; jcell<SMesh.shape().ny()-1; ++jcell) {
2551        for (size_t kcell=0; kcell<SMesh.shape().nz()-1; ++kcell) {        for (size_t kcell=0; kcell<SMesh.shape().nz()-1; ++kcell) {
2552          bool intersected = false;          bool intersected = false;
2553    
2554            MotherCellList cellList;
2555    
2556          const Cell& currentCell = SMesh.cell(icell, jcell, kcell);          const Cell& currentCell = SMesh.cell(icell, jcell, kcell);
2557                            
2558          const size_t nx          const size_t nx
# Line 2771  __createSurface(Structured3DMesh& SMesh, Line 2729  __createSurface(Structured3DMesh& SMesh,
2729                                                    localTriangleList,                                                    localTriangleList,
2730                                                    edgesRef,                                                    edgesRef,
2731                                                    S, cutEdges,                                                    S, cutEdges,
2732                                                    currentCell);                                                    currentCell,
2733                                                      cellList);
2734              if(!(jobDone2)) {              if(!(jobDone2)) {
2735                  jobDone=jobDone2;                  jobDone=jobDone2;
2736              }              }
# Line 2787  __createSurface(Structured3DMesh& SMesh, Line 2746  __createSurface(Structured3DMesh& SMesh,
2746                cellListObject.push_front(*itcell);                cellListObject.push_front(*itcell);
2747                ++itcell;                ++itcell;
2748              }              }
             hexaTetraAssociation[&currentCell]  
               = std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,  
               ReferenceCounting<std::vector<Edge> > >  
               (pT,pAddedEdges);  
2749            } else {            } else {
2750                ffout(0)<<"------jobDone false\n";                ffout(0)<<"------jobDone false\n";
             hexaTetraAssociation[&currentCell]  
               = std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,  
               ReferenceCounting<std::vector<Edge> > >  
               (pT,pAddedEdges);  
2751              SurfaceMeshGenerator::Internals::MotherCellList::iterator itcell              SurfaceMeshGenerator::Internals::MotherCellList::iterator itcell
2752                =cellList.begin();                =cellList.begin();
2753              for (EdgePairList::iterator t = localTriangleList.begin();              for (EdgePairList::iterator t = localTriangleList.begin();
# Line 2807  __createSurface(Structured3DMesh& SMesh, Line 2758  __createSurface(Structured3DMesh& SMesh,
2758              }              }
2759            }            }
2760            localTriangleList.clear();            localTriangleList.clear();
           cellList.clear();  
2761          }          }
2762        }        }
2763        }
2764      }
2765  }  }
2766    
2767  void  void
# Line 3116  __generateMesh(Structured3DMesh& SMesh, Line 3068  __generateMesh(Structured3DMesh& SMesh,
3068    default: {    default: {
3069      ReferenceCounting<ObjectToTreat> o = new ObjectToTreat();      ReferenceCounting<ObjectToTreat> o = new ObjectToTreat();
3070      (*o).setShape(object.shape());      (*o).setShape(object.shape());
3071      __marchingTetrahedra(SMesh, *o,      __marchingTetrahedra(SMesh, *o, 0);
                          0);  
3072    
3073      objectStack.push(o);      objectStack.push(o);
3074      std::cout << "Marching cube sur " << objectStack.top() << '\n';      std::cout << "Marching cube sur " << objectStack.top() << '\n';
# Line 3129  ReferenceCounting<Vector<Triangle> > Line 3080  ReferenceCounting<Vector<Triangle> >
3080  SurfaceMeshGenerator::Internals::  SurfaceMeshGenerator::Internals::
3081  __marchingTetrahedra(Structured3DMesh& SMesh,  __marchingTetrahedra(Structured3DMesh& SMesh,
3082                       ObjectToTreat& O,                       ObjectToTreat& O,
                      //                        const Shape& S,  
                      //                        const Reference::Representation& ref,  
3083                       size_t nobject)                       size_t nobject)
3084  {  {
3085    EdgePairList triangleListes;    EdgePairList triangleListes;
# Line 3158  generateSurfacicMesh(const Domain& Omega Line 3107  generateSurfacicMesh(const Domain& Omega
3107                       SurfaceMeshOfTriangles& s_mesh,                       SurfaceMeshOfTriangles& s_mesh,
3108                       Structured3DMesh& SMesh)                       Structured3DMesh& SMesh)
3109  {  {
3110    typedef std::map<TinyVector<3>, Object*,    typedef std::map<TinyVector<3>, Object*> ObjectReferences;
     SurfaceMeshGenerator::Internals::__ltv> ObjectReferences;  
3111    
3112    const Object& object = Omega.object();    const Object& object = Omega.object();
3113    const Scene& scene = Omega.scene();    const Scene& scene = Omega.scene();
# Line 3178  generateSurfacicMesh(const Domain& Omega Line 3126  generateSurfacicMesh(const Domain& Omega
3126                                           s_mesh);                                           s_mesh);
3127    
3128    int n0=0, n1=0;    int n0=0, n1=0;
   //  plot2(s_mesh, (*__internals).toTreatHexahedra, SMesh);  
   plotmedit(scene.nbObjects(),s_mesh, (*__internals).toTreatHexahedra,n0,n1);  
3129    
3130    #warning REMOVE THIS CALL
3131      std::set<const Cell*> toTreatHexahedra;
3132      plotmedit(scene.nbObjects(), s_mesh, toTreatHexahedra, n0, n1);
3133  }  }
3134    
3135  SurfaceMeshGenerator::SurfaceMeshGenerator()  SurfaceMeshGenerator::SurfaceMeshGenerator()
# Line 3260  void plotmedit(const size_t& nobject, Line 3209  void plotmedit(const size_t& nobject,
3209      sortie<<s_mesh.vertex(i)[0]<<" "<<s_mesh.vertex(i)[1]<<" "<<s_mesh.vertex(i)[2]<<" 0"<<" \n";      sortie<<s_mesh.vertex(i)[0]<<" "<<s_mesh.vertex(i)[1]<<" "<<s_mesh.vertex(i)[2]<<" 0"<<" \n";
3210    }    }
3211    
   /*  for (size_t i=0; i<s_mesh.numberOfCells(); ++i)  
       {  
       const Triangle& T = s_mesh.cell(i);  
       sortie<<T(0)[0]<<" "<<T(0)[1]<<" "<<T(0)[2]<<" 0"<<" \n";  
       sortie<<T(1)[0]<<" "<<T(1)[1]<<" "<<T(1)[2]<<" 0"<<" \n";  
       sortie<<T(2)[0]<<" "<<T(2)[1]<<" "<<T(2)[2]<<" 0"<<" \n";  
       }*/  
3212    sortie<<" \n";    sortie<<" \n";
3213    sortie<<"Triangles"<<"\n";    sortie<<"Triangles"<<"\n";
3214    sortie<<s_mesh.numberOfCells()<<"\n";    sortie<<s_mesh.numberOfCells()<<"\n";
# Line 3277  void plotmedit(const size_t& nobject, Line 3219  void plotmedit(const size_t& nobject,
3219      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<<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;
3220            
3221      // sortie<<3*i+1+8*toTreatHexahedra.size()<<" "<<3*i+2+8*toTreatHexahedra.size()<<" //"<<3*i+3+8*toTreatHexahedra.size();//<<" "<<1<<"\n";      // sortie<<3*i+1+8*toTreatHexahedra.size()<<" "<<3*i+2+8*toTreatHexahedra.size()<<" //"<<3*i+3+8*toTreatHexahedra.size();//<<" "<<1<<"\n";
3222      if(i<n0)  //     if(i<n0)
3223        sortie<<" 1"<<"\n";  //       sortie<<" 1"<<"\n";
3224      else if(i<n1)  //     else if(i<n1)
3225        sortie<<" 2"<<"\n";  //       sortie<<" 2"<<"\n";
3226      else  //     else
3227        sortie<<" 4"<<"\n";  //       sortie<<" 4"<<"\n";
3228        sortie << " " << T.reference() << '\n';
3229    }    }
3230    ffout(0)<<"fin ecriture triangles\n";    ffout(0)<<"fin ecriture triangles\n";
3231    sortie<<" \n";    sortie<<" \n";
# Line 3294  void plotmedit(const size_t& nobject, Line 3237  void plotmedit(const size_t& nobject,
3237    sortie<<"End"<<" \n";    sortie<<"End"<<" \n";
3238    ffout(0)<<"fin ecriture intersection.mesh\n";    ffout(0)<<"fin ecriture intersection.mesh\n";
3239  }  }
   
   
 // #include <vtkPolyDataMapper.h>  
 // #include <vtkRenderWindow.h>  
 // #include <vtkRenderWindowInteractor.h>  
 // #include <vtkTriangle.h>  
 // #include <vtkHexahedron.h>  
 // #include <vtkTetra.h>  
 // #include <vtkUnstructuredGrid.h>  
 // #include <vtkDataSetMapper.h>  
 // #include <vtkDataSetToPolyDataFilter.h>  
 // #include <vtkDoubleArray.h>  
   
 // void plot2(const SurfaceMeshOfTriangles& s,  
 //         std::set<const Cell*>& hlist,  
 //         Structured3DMesh& m)  
 // {  
 //   vtkRenderer *ren1= vtkRenderer::New();  
   
 //   vtkActor *anActor1 = vtkActor::New();  
 //   vtkActor *anActor2 = vtkActor::New();  
   
 //   {  
 //     vtkPoints* vertices = vtkPoints::New();  
   
 //     vertices->SetNumberOfPoints(m.numberOfVertices());  
   
 //     for (size_t i=0; i<m.numberOfVertices(); ++i) {  
 //       const Vertex& X = m.vertex(i);  
 //       vertices->InsertPoint(i,X[0],X[1],X[2]);  
 //     }  
   
 //     vtkHexahedron* hexahedron = vtkHexahedron::New();  
   
 //     vtkUnstructuredGrid* anHexahedraGrid = vtkUnstructuredGrid::New();  
 //     anHexahedraGrid->Allocate(m.numberOfCells(),1);  
   
 //     for (std::set<const Cell*>::iterator i = hlist.begin();  
 //       i != hlist.end(); ++i) {  
 //       const CartesianHexahedron& H = static_cast<const CartesianHexahedron&>(*(*i));  
 //       for (size_t j=0; j<8; ++j)  
 //      hexahedron->GetPointIds()->SetId(j,m.vertexNumber(H(j)));  
   
 //       anHexahedraGrid->InsertNextCell(hexahedron->GetCellType(),  
 //                                    hexahedron->GetPointIds());  
 //     }  
   
 //     anHexahedraGrid->SetPoints(vertices);  
   
 //     vtkDoubleArray* values = vtkDoubleArray::New();  
 //     values->Resize(m.numberOfVertices());  
 //     for (size_t i=0; i<m.numberOfVertices(); ++i) {  
 //       values->InsertNextValue(1);  
 //     }  
   
 //     anHexahedraGrid->GetPointData()->SetScalars(values);  
   
 //     // We create an instance of vtkPolyDataMapper to map the polygonal data  
 //     // into graphics primitives. We connect the output of the cone souece  
 //     // to the input of this mapper  
   
 //     vtkDataSetMapper* mapper = vtkDataSetMapper::New();  
 //     mapper->SetInput(anHexahedraGrid);  
   
 //     // create an actor to represent the cone. The actor coordinates rendering of  
 //     // the graphics primitives for a mapper. We set this actor's mapper to be  
 //     // coneMapper which we created above.  
 //     //  
   
 //     vtkProperty* prop = vtkProperty::New();  
 //     prop->SetRepresentationToWireframe();  
   
 //     anActor1->SetMapper( mapper );  
 //     anActor1->SetProperty(prop);  
   
 //   }  
   
   
 //   {  
 //     vtkPoints* vertices = vtkPoints::New();  
   
 //     vertices->SetNumberOfPoints(s.numberOfVertices());  
   
 //     for (size_t i=0; i<s.numberOfVertices(); ++i) {  
 //       const Vertex& X = s.vertex(i);  
 //       vertices->InsertPoint(i,X[0],X[1],X[2]);  
 //     }  
   
 //     vtkTriangle* triangle = vtkTriangle::New();  
   
 //     vtkUnstructuredGrid* aTriangleGrid = vtkUnstructuredGrid::New();  
 //     aTriangleGrid->Allocate(s.numberOfCells(),1);  
   
 //     for (size_t i=0; i<s.numberOfCells(); ++i) {  
 //       const Triangle& T = s.cell(i);  
 //       triangle->GetPointIds()->SetId(0,s.vertexNumber(T(0)));  
 //       triangle->GetPointIds()->SetId(1,s.vertexNumber(T(1)));  
 //       triangle->GetPointIds()->SetId(2,s.vertexNumber(T(2)));  
   
 //       aTriangleGrid->InsertNextCell(triangle->GetCellType(),  
 //                                  triangle->GetPointIds());  
 //     }  
   
 //     aTriangleGrid->SetPoints(vertices);  
   
 //     vtkDoubleArray* values = vtkDoubleArray::New();  
 //     values->Resize(s.numberOfVertices());  
 //     for (size_t i=0; i<s.numberOfVertices(); ++i) {  
 //       values->InsertNextValue(0);  
 //     }  
   
 //     aTriangleGrid->GetPointData()->SetScalars(values);  
   
 //     // We create an instance of vtkPolyDataMapper to map the polygonal data  
 //     // into graphics primitives. We connect the output of the cone souece  
 //     // to the input of this mapper  
   
 //     vtkDataSetMapper* mapper = vtkDataSetMapper::New();  
 //     mapper->SetInput(aTriangleGrid);  
   
   
 //     // create an actor to represent the cone. The actor coordinates rendering of  
 //     // the graphics primitives for a mapper. We set this actor's mapper to be  
 //     // coneMapper which we created above.  
 //     //  
 //     anActor2->SetMapper( mapper );  
   
 //   }  
   
 //   // Create the Renderer and assign actors to it. A renderer is like a  
 //   // viewport. It is part or all of a window on the screen and it is  
 //   // responsible for drawing the actors it has.  We also set the background  
 //   // color here  
 //   //  
 //   //  ren1->AddActor( anActor1 );  
 //   ren1->AddActor( anActor2 );  
 //   ren1->SetBackground( 0.1, 0.2, 0.4 );  
   
 //   //  
 //   // Finally we create the render window which will show up on the screen  
 //   // We put our renderer into the render window using AddRenderer. We also  
 //   // set the size to be 300 pixels by 300  
 //   //  
 //   vtkRenderWindow *renWin = vtkRenderWindow::New();  
 //   renWin->AddRenderer( ren1 );  
 //   renWin->SetSize( 300, 300 );  
   
 //   vtkRenderWindowInteractor* iren = vtkRenderWindowInteractor::New();  
 //   iren->SetRenderWindow(renWin);  
   
 //   // now we loop over 360 degreeees and render the cone each time  
 //   //  
 //   // render the image  
 //   renWin->Render();  
 //   iren->Start();  
 //   //  
 //   // Free up any objects we created  
 //   //  
 //   //   cone->Delete();  
 //   //   coneMapper->Delete();  
 //   //   coneActor->Delete();  
 //   //   ren1->Delete();  
 // }  

Legend:
Removed from v.1.28  
changed lines
  Added in v.1.29

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26