/[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.12 by dobj, Fri Jun 20 15:58:50 2003 UTC revision 1.13 by dobj, Wed Jul 2 15:17:02 2003 UTC
# Line 52  Line 52 
52  #include <set>  #include <set>
53  #include <iostream>  #include <iostream>
54    
55    
56  class SurfaceMeshGenerator::Internals  class SurfaceMeshGenerator::Internals
57  {  {
58  private:  private:
59      class IntersectionPoints;    class IntersectionPoints
60      class TriangleCut;    //class SurfaceMeshGenerator::Internals::IntersectionPoints
61      friend class SurfaceMeshGenerator;  {
62    public :
63        TinyVector<2,size_t> numpair;
64        std::vector<const  Vertex* > listofPoints;
65        TinyVector< 3,size_t > cutTriangle1;
66        TinyVector< 3,size_t > cutTriangle2;
67    
68        inline const size_t size() const {
69            return listofPoints.size();
70        }
71    
72        void addPoint(const Vertex*&  V) {
73            listofPoints.push_back( V );
74        }
75      
76    
77      //constructor
78        IntersectionPoints() {
79            numpair=0;
80            cutTriangle1=0;
81            cutTriangle2=0;
82        }
83    
84      IntersectionPoints(const IntersectionPoints& P) {
85          numpair=P.numpair;
86          listofPoints=P.listofPoints;
87          cutTriangle1=P.cutTriangle1;
88          cutTriangle2=P.cutTriangle2;
89      }
90    
91      //! Destructor.
92      ~IntersectionPoints()
93      {
94        ;
95      }
96        
97    };
98    
99    
100      class TriangleCut;
101      friend class SurfaceMeshGenerator;
102    
103    struct __ltv    struct __ltv
104    {    {
# Line 73  private: Line 114  private:
114      }      }
115    };    };
116    
   struct __ltv2  
   {  
       bool operator()(const std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> >& s1,  
                       const std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> >& s2) const  
   {  
       if(s1.first[0][0].first==s2.first[0][0].first) {  
           if(s1.first[0][0].second==s2.first[0][0].second) {  
               if(s1.first[1][0].first==s2.first[1][0].first) {  
                   return s1.first[1][0].second<s2.first[1][0].second;  
               }  
               else {  
                   return s1.first[1][0].first<s2.first[1][0].first;  
               }  
           }  
           else {  
               return s1.first[0][0].second<s2.first[0][0].second;  
           }  
       }  
       else {  
           return s1.first[0][0].first<s2.first[0][0].first;  
       }  
   }  
   };  
   
117    typedef std::map<TinyVector<3,real_t>,    typedef std::map<TinyVector<3,real_t>,
118                     std::vector<SurfaceMesh*>, __ltv> VRRefMesh;                     std::vector<SurfaceMesh*>, __ltv> VRRefMesh;
119    
120    typedef std::list<TinyVector<3,Edge::Pair> > EdgePairList;    typedef std::list<TinyVector<3,Edge::Pair> > EdgePairList;      
   typedef std::list<TinyVector<3,Edge::Pair> * > EdgePairList2;  
   typedef std::pair< const TinyVector<3,Edge::Pair>*,const TinyVector<3,Edge::Pair> *> PairTriangle;  
   typedef std::pair< const TinyVector<3,TinyVector<2,TinyVector<2,Edge::Pair> > >*,  
       const TinyVector<3,TinyVector<2,TinyVector<2,Edge::Pair> > > *> PairTriangleEdges;  
   
   struct __ltvpair  
   {  
     bool operator()(const  PairTriangle & s1, const PairTriangle & s2) const  
     {  
       if(s1.first==s2.first) {  
         return s1.second<s2.second;  
       }  
       else {  
         return s1.first<s2.first;  
       }  
     }  
   };  
   
   struct __ltvpairedges  
   {  
       bool operator()(const  PairTriangleEdges & s1, const PairTriangleEdges & s2) const  
   {  
       if(s1.first==s2.first) {  
           return s1.second<s2.second;  
       }  
       else {  
           return s1.first<s2.first;  
       }  
   }  
   };  
       
121    typedef std::list<const Cell*> MotherCellList;    typedef std::list<const Cell*> MotherCellList;
122    
123    typedef std::map<const Cell*,    typedef std::map<const Cell*,
# Line 139  private: Line 125  private:
125                               ReferenceCounting<std::vector<Edge> > > > HexaTetraAssociation;                               ReferenceCounting<std::vector<Edge> > > > HexaTetraAssociation;
126    
127    typedef std::map<Edge::Pair, Vertex*> VerticesList;    typedef std::map<Edge::Pair, Vertex*> VerticesList;
   typedef std::vector< TinyVector<3,size_t> > TriangleVerticesNumber;  
128    typedef std::map<TinyVector<3>, Object*, SurfaceMeshGenerator::Internals::__ltv> ObjectReferences;    typedef std::map<TinyVector<3>, Object*, SurfaceMeshGenerator::Internals::__ltv> ObjectReferences;
   typedef std::pair<size_t,size_t> ObjectPair;  
129    
130    //! list of cut edges and intersection vertex.    //! list of cut edges and intersection vertex.
131    std::vector<ReferenceCounting<std::map<Edge::Pair, size_t> > > edgesRef;    std::vector<ReferenceCounting<std::map<Edge::Pair, size_t> > > edgesRef;
   std::vector<ReferenceCounting<std::map<Edge::Pair, size_t> > > edgesRefCreate;  
132    
133    typedef  TinyVector<3,TinyVector<2,TinyVector<2,Edge::Pair> > >  TriangleVertices;    typedef std::map<const Cell*, std::vector<ReferenceCounting<EdgePairList> > >
134    typedef std::list< TriangleVertices > TriangleVerticesEdgesList;    TriangleListesIntersectionEdgesPerHexa;
135    typedef std::list< TriangleVertices *> TriangleVerticesEdgesList2;    TriangleListesIntersectionEdgesPerHexa triangleListesIntersectionEdgesPerHexa;
   std::vector<ReferenceCounting<TriangleVerticesEdgesList> > triangleListesIntersectionEdges;  
   typedef std::map<TriangleVertices*,const TriangleVertices*> TriangleVerticesMap;  
   
   std::map<const Cell*, std::vector<ReferenceCounting<TriangleVerticesEdgesList> > >  
       triangleListesIntersectionEdgesPerHexa;  
   std::map<const Cell*, std::vector<ReferenceCounting< TriangleVerticesMap> > >  
       triangleListesIntersectionEdgesTempPerHexa;  
   typedef std::vector< TriangleVertices > TriangleVerticesEdges;  
   typedef std::map<const TriangleVertices*,TriangleVerticesEdges > TriangleTreatMap;  
   TinyVector<2, TriangleTreatMap > triangleTreat;  
136    
137    std::vector<ReferenceCounting<EdgePairList> > triangleListes;    std::vector<ReferenceCounting<EdgePairList> > triangleListes;
138    TinyVector<2,std::vector< TriangleVerticesMap::iterator >  > triangleEraseEdgesTemp;    std::vector<ReferenceCounting<EdgePairList> > localTriangleListesEdges;
   std::vector<ReferenceCounting<TriangleVerticesEdgesList> > localTriangleListesEdges;  
   std::vector<ReferenceCounting<TriangleVerticesEdgesList2> > localTriangleList2Edges;  
139    EdgePairList localTriangleList;    EdgePairList localTriangleList;
140    
141    
142    //! This list is used to provide an orientation to the Surfacic Elements.    //! This list is used to provide an orientation to the Surfacic Elements.
143    std::list<TinyVector<3> > InsideVertexList;    std::list<TinyVector<3> > InsideVertexList;
144    
   //For the intersection  
   std::map< ObjectPair , std::map< PairTriangleEdges , IntersectionPoints, __ltvpairedges > >  
       listVertexIntersectionPerPairEdges;  
   
   typedef std::map<std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> >,Vertex*,__ltv2>             ListVerticesIntersectionPerEdges;  
   ListVerticesIntersectionPerEdges   listVerticesIntersectionPerEdges;  
     
   
   std::map<const Cell*, std::vector<ReferenceCounting< TriangleVerticesEdgesList> > >  
       triangleListesIntersectionPerHexa;  
   std::vector<ReferenceCounting< TriangleVerticesEdges > > listTriangleintersection;  
   TriangleVerticesEdges  localListTriangleintersection;  
     
   
145    //! list of cut edges and intersection vertex.    //! list of cut edges and intersection vertex.
146    std::vector<ReferenceCounting<VerticesList> > verticesListes;    std::vector<ReferenceCounting<VerticesList> > verticesListes;
147    std::vector<ReferenceCounting<VerticesList> > verticesListesIntersection;    std::vector<ReferenceCounting<VerticesList> > verticesListesIntersection;
# Line 192  private: Line 149  private:
149    HexaTetraAssociation hexaTetraAssociation;    HexaTetraAssociation hexaTetraAssociation;
150    std::set<const Cell*> toTreatHexahedra;    std::set<const Cell*> toTreatHexahedra;
151    
152      typedef std::list<ReferenceCounting<Vector<Triangle> > > ListOfTriangleList;
153      typedef std::list<ReferenceCounting<Vector<Vertex*> > > ListOfVertexList;
154    
155    MotherCellList cellList;    ListOfVertexList  listOfVertexList;
156      ListOfVertexList  listOfVertexListIntersection;
157      std::map<const Vertex*,const Vertex*>   listOfVertexMesh;
158      ListOfTriangleList listOfTriangleList;
159      ListOfTriangleList listOfTriangleListIntersection;
160      
161      typedef std::map<const Cell*, std::list<Triangle*> > MapCellTriangle;
162      typedef std::list<ReferenceCounting< MapCellTriangle > > ListOfHexaToTriangle;
163      ListOfHexaToTriangle listOfHexaToTriangle;
164    
165      typedef std::pair<const Triangle*,const Triangle*> PairTriangleNew;
166      //std::vector<Vertex*> listVertexIntersectionNew;
167      typedef std::map< std::pair<const Vertex*,const Vertex*> ,const Vertex* > PIntersection;
168      PIntersection listVertexIntersectionNew;
169    
170      typedef std::map< PairTriangleNew , IntersectionPoints> PairTriangleToIntersectionPoints;
171      PairTriangleToIntersectionPoints pairTriangleToIntersectionPoints;
172    
173      std::vector<ReferenceCounting< std::map<const Vertex*,size_t> > > edgesRefVertex;
174    
175      std::vector<Triangle> triangleListIntersectionNew;
176    
177      std::map<const Vertex*,TinyVector<2,const Triangle*> > vertexInTriangles;
178      std::map<const Triangle*,TinyVector<2,const Vertex*> > triangleWithVertex;
179    
180      MotherCellList cellList;
181      std::vector<ReferenceCounting<MotherCellList> > cellListObject;
182    TinyVector<3> __splitEdge(const Edge& e,    TinyVector<3> __splitEdge(const Edge& e,
183                              const Reference::Representation& ref,                              const Reference::Representation& ref,
184                              const Shape& S);                              const Shape& S);
# Line 214  private: Line 198  private:
198                                       ObjectReferences& otherReferences);                                       ObjectReferences& otherReferences);
199    
200    void __calculatePointsIntersection(const Scene& scene,    void __calculatePointsIntersection(const Scene& scene,
                                      size_t nobjectToTreat,  
201                                       size_t no,                                       size_t no,
202                                       size_t no2);                                       size_t no2);
203    
# Line 226  private: Line 209  private:
209                          const TinyVector<3,TinyVector<3> > & triangleVector );                          const TinyVector<3,TinyVector<3> > & triangleVector );
210        
211    void __InTriangle(const Vertex& P,    void __InTriangle(const Vertex& P,
212                      const TinyVector<3, Vertex* > & V,                      TinyVector<3,const Vertex *> & V,
213                      size_t& ncompt);                      size_t& ncompt);
214    
215    void __inout(const size_t& numobject,    void __inout(const size_t& ncase,
216                   const size_t& numobject,
217                 const int& ttn0,                 const int& ttn0,
218                 const TinyVector<2,Edge::Pair>& e0,                 const Vertex*& V0,
219                 const TinyVector<2,Edge::Pair>& e1,                 const Vertex*& V1,
220                 size_t& ii0,                 TriangleCut& K);
                size_t& ii1,TriangleCut& K);  
221    
222    TinyVector<4, real_t>  __EquationPlan(const TinyVector<3, Vertex* > & V);    TinyVector<4, real_t>  __EquationPlan(const TinyVector<3,const Vertex *> & V);
223    void __DeterminatePoint(const size_t& ncase,    void __DeterminatePoint(const size_t& ncase,
224                            const TinyVector<2,Edge::Pair> & eo0,                            const Vertex* & Vo0,
225                            const TinyVector<2,Edge::Pair> & eo1,                            const Vertex* & Vo1,
226                            const Vertex & Vo0,                            TinyVector<3, const Vertex* > & V,
                           const Vertex & Vo1,  
                           const TinyVector<3, Vertex* > & V,  
227                            IntersectionPoints& Pinter);                            IntersectionPoints& Pinter);
228    void __CalculPoint(const size_t& objectNumber,    void __CalculPoint(const size_t& objectNumber,
229                       const size_t& numobject,                       const size_t& numobject,
230                       SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::const_iterator & currentTriangle,                       std::list<Triangle*>::iterator & currentTriangle,
231                       SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::const_iterator & currentTriangle1,                       std::list<Triangle*>::iterator & currentTriangle1,
232                       IntersectionPoints & Pinter);                       IntersectionPoints & Pinter);
233    
234    void __createPointIntersection(const size_t& numobject,    void __createPointIntersection(const size_t& numobject,
235                                   const size_t& objectNumber,                                   const size_t& objectNumber,
236                                   const Cell*,                                   std::list<Triangle*>::iterator currentTriangle1,
237                                   SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::const_iterator currentTriangle1,                                   std::list<Triangle*>::iterator currentTriangle);
                                  SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::const_iterator currentTriangle);  
238        
239    std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > __determinateNumber(const int& p,    const Vertex& __determinateNumber(const int& p,
240          const TinyVector<3,std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > >& pointsint);                                      const TinyVector<3,TinyVector<2,const Vertex*> >& pointsint);
241    void __determinateCase(const int& ncase,   void __determinateCase(const int& ncase,
242                           const TriangleCut & K,                             const TriangleCut & K,
243                           Vertex & A,Vertex & B,Vertex & C,                             TinyVector<3,int>& tt,
244                           TinyVector<2,TinyVector<2,Edge::Pair> >& e0,                             TinyVector<3,int>& pp);
245                           TinyVector<2,TinyVector<2,Edge::Pair> >& e1,    void __addTriangle(const Vertex& V0,
246                           TinyVector<2,TinyVector<2,Edge::Pair> >& e2,                       const Vertex& V1,
247                           TinyVector<2,int>& tt,                       const Vertex& V2,
248                           TinyVector<2,int>& pp);                       const Cell* cell);
   void __addTriangle(const size_t& numobject,  
                      const TinyVector<2,TinyVector<2,Edge::Pair> >& numA,  
                      const TinyVector<2,TinyVector<2,Edge::Pair> >& numB,  
                      const TinyVector<2,TinyVector<2,Edge::Pair> >& numC,  
                      const Cell* cell,  
                      size_t& kill);  
249    void __createCase1(int ncase,    void __createCase1(int ncase,
250                       const TriangleCut& K,                       const TriangleCut& K,
251                       const Cell* cell,                       const Cell* cell);
                      size_t & kill);  
252    
253    void __createCase2(int ncase,    void __createCase2(int ncase,
254                       const TriangleCut& K,                       const TriangleCut& K,
255                       const Cell* cell,                       const Cell* cell);
                      size_t & kill);  
256    
257    void __createLocalListIntersection(const size_t& objectNumber,    void __createLocalListIntersection(const size_t& objectNumber,
258                                       const TinyVector<3,TinyVector<2,TinyVector<2,Edge::Pair> > >& currentTriangle,                                       const Triangle*& currentTriangle,
259                                       TriangleCut & K,                                       TriangleCut & K,
260                                       const Cell* cell);                                       const Cell* cell);
261        
262    TinyVector<3,int> __determinateEdgeCut(const size_t& numobject,    TinyVector<3,TinyVector<2,int> > __determinateEdgeCut(const size_t& numobject,
263                                           const size_t & objectNumber,                                           const size_t & objectNumber,
264                                  const TinyVector<3,TinyVector<2,TinyVector<2,Edge::Pair> > >& currentTriangle,                                           const Triangle*& currentTriangle,
265                                           const Cell* cell);                                           const Cell* cell);
266    
267    TinyVector<3,size_t> __putEdgeRef(const size_t& ncase,    TinyVector<3,size_t> __putEdgeRef(const size_t& ncase,
268                                      const TinyVector<3,TinyVector<2,TinyVector<2,Edge::Pair> > >& currentTriangle,                                      const Triangle& currentTriangle,
269                                      TriangleCut& K);                                      TriangleCut& K);
270    
271    void __createListTriangle(const TriangleCut& K,    void __createListTriangle(const TriangleCut& K,
272                              const Cell* cell,                              const Cell* cell);
                             size_t & kill);  
273    
274      void __createByEdges(const TriangleCut& K,const Cell* cell) ;
275      void  __createNotIn(const TriangleCut& K,const Cell* cell);
276    void __createTrianglesIntersection(const size_t& numobject,    void __createTrianglesIntersection(const size_t& numobject,
277                                       const size_t& objectNumber,                                       const size_t& objectNumber);
                                      const size_t& nobject);  
278    
279    bool __createTriangleSurface(const size_t& nobject,    bool __createTriangleSurface(const size_t& nobject,
280                                 const Shape& S ,                                 const Shape& S ,
# Line 317  private: Line 289  private:
289                         size_t nobject,                         size_t nobject,
290                         std::set<const Cell*>& cuttedHexahedra,                         std::set<const Cell*>& cuttedHexahedra,
291                         std::set<const Cell*>& meshedHexahedra);                         std::set<const Cell*>& meshedHexahedra);
292    Vertex* __determinateVertex(size_t numobject,  /// @todo FUNCTION TO REMOVE
                              const TinyVector<2,TinyVector<2,Edge::Pair> >& e0);  
   void __createMeshIntersection(size_t nobjectToTreat,  
                                 const Scene & scene);  
   
   /// @todo FUNCTION TO REMOVE  
293    static void __getIntersectionReferences(const Intersection& I,    static void __getIntersectionReferences(const Intersection& I,
294                                            std::map<TinyVector<3>, Object*, Internals::__ltv>& otherReferences);                                            std::map<TinyVector<3>, Object*, Internals::__ltv>& otherReferences);
295    void plotliste(size_t nobject);    void plotliste();
296    void plotlistenew(size_t nobject,size_t nobjectToTreat);    bool triangleplat(const Vertex& V0,
297                        const Vertex& V1,
298    void printpoint(size_t numobject,                      const Vertex& V2);
                   const TinyVector<2,TinyVector<2,Edge::Pair> >& e0 , TriangleCut K);  
   bool triangleplat(size_t numobject,const TinyVector<2,Edge::Pair>& e0 ,  
                     const TinyVector<2,Edge::Pair>& e1,  
                     const TinyVector<2,Edge::Pair>& e2 );  
299    
300  public:  public:
301    Internals()    Internals()
# Line 341  public: Line 304  public:
304    }    }
305  };  };
306    
 class SurfaceMeshGenerator::Internals::IntersectionPoints  
 {  
 public :  
     TinyVector<2,size_t> numpair;  
     std::vector< std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > > listofPoints;  
     TinyVector< 3,size_t > cutTriangle1;  
     TinyVector< 3,size_t > cutTriangle2;  
   
     inline const size_t size() const {  
         return listofPoints.size();  
     }  
   
     void addPoint(std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> >&  V) {  
         listofPoints.push_back( V );  
     }  
     
   
   //constructor  
     IntersectionPoints() {  
         numpair=0;  
         cutTriangle1=0;  
         cutTriangle2=0;  
     }  
307    
   IntersectionPoints(const IntersectionPoints& P) {  
       numpair=P.numpair;  
       listofPoints=P.listofPoints;  
       cutTriangle1=P.cutTriangle1;  
       cutTriangle2=P.cutTriangle2;  
   }  
   
   //! Destructor.  
   ~IntersectionPoints()  
   {  
     ;  
   }  
       
 };  
308    
309  class SurfaceMeshGenerator::Internals::TriangleCut  class SurfaceMeshGenerator::Internals::TriangleCut
310    : public Triangle
311  {  {
312  public :  public :
313      TinyVector<3,TinyVector<2,TinyVector<2,Edge::Pair> > > listpairedges;      TinyVector<3,TinyVector<2,int> > edgecut;
     TinyVector<3,int> edgecut;  
314      TinyVector<3,size_t> isInTriangle;      TinyVector<3,size_t> isInTriangle;
315      TinyVector<3,std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > > pointsIntersection;      TinyVector<3,TinyVector<2,const Vertex*> > pointsIntersection;
316      std::vector< std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > > pointsin;      std::map< const Vertex* ,const Vertex*> pointsin;
317      size_t numobject;      size_t numobject;
318      size_t otherobject;      size_t otherobject;
319    
# Line 395  public : Line 321  public :
321      TriangleCut() {      TriangleCut() {
322      }      }
323      TriangleCut(const TriangleCut& K) {      TriangleCut(const TriangleCut& K) {
         listpairedges=K.listpairedges;  
324          edgecut=K.edgecut;          edgecut=K.edgecut;
325          isInTriangle=K.isInTriangle;          isInTriangle=K.isInTriangle;
326          pointsIntersection=K.pointsIntersection;          pointsIntersection=K.pointsIntersection;
# Line 403  public : Line 328  public :
328          numobject=K.numobject;          numobject=K.numobject;
329          otherobject=K.otherobject;          otherobject=K.otherobject;
330      }      }
331      TriangleCut(const TinyVector<2,TinyVector<2,Edge::Pair> > & e0,      TriangleCut(const Vertex& v0,
332                  const TinyVector<2,TinyVector<2,Edge::Pair> > & e1,                  const Vertex& v1,
333                  const TinyVector<2,TinyVector<2,Edge::Pair> > & e2,                  const Vertex& v2,
334                  size_t num,size_t num2) {                  size_t num,
335          listpairedges[0]=e0;                  size_t num2)
336          listpairedges[1]=e1;      {
337          listpairedges[2]=e2;          __vertices[0] = (Vertex*)&v0;
338            __vertices[1] = (Vertex*)&v1;
339            __vertices[2] = (Vertex*)&v2;
340          numobject=num;          numobject=num;
341          otherobject=num2;          otherobject=num2;
342          edgecut=0;          edgecut=0;
343          isInTriangle=0;          isInTriangle=0;
344      }      }
345        
346      //! Destructor.      //! Destructor.
347      ~TriangleCut()      ~TriangleCut()
348      {      {
# Line 514  __createListIntersection(const Shape& S, Line 441  __createListIntersection(const Shape& S,
441              std::swap(cutEdges[1], cutEdges[2]);              std::swap(cutEdges[1], cutEdges[2]);
442    
443              //! Addes a triangle to the list.              //! Addes a triangle to the list.
444              (*localTriangleListesEdges[objectNumber]).push_front(              (*localTriangleListesEdges[objectNumber]).push_front( TinyVector<3,Edge::Pair>());
445                                                            TinyVector<3,TinyVector<2,TinyVector<2,Edge::Pair> > >());            //  cellList.push_front(&currentCell);
             cellList.push_front(&currentCell);  
446    
447              //! Sets the triangle vertices (actually just store the edge where it lives).              //! Sets the triangle vertices (actually just store the edge where it lives).
448              TriangleVerticesEdgesList::iterator currentTriangleEdges              EdgePairList::iterator currentTriangleEdges
449                  = (*localTriangleListesEdges[objectNumber]).begin();                  = (*localTriangleListesEdges[objectNumber]).begin();
450    
451              for(size_t i=0; i<cutEdges.size(); ++i) {              for(size_t i=0; i<cutEdges.size(); ++i) {
452                  Edge::Pair e=(*cutEdges[i]);                  Edge::Pair e=(*cutEdges[i]);
453                  TinyVector<2,TinyVector<2,Edge::Pair> > pe;                  (*currentTriangleEdges)[i] = e;
                 pe[0][0]=e;  
                 pe[1][0]=e;  
                 pe[0][1]=e;  
                 pe[1][1]=e;  
                 (*currentTriangleEdges)[i] = pe;  
454              }              }
455    
456              //! Stores the Edges on which the vertex will be computed.              //! Stores the Edges on which the vertex will be computed.
# Line 546  __createListIntersection(const Shape& S, Line 467  __createListIntersection(const Shape& S,
467                  cutEdgesMap[(Edge::Pair)*(cutEdges[ie])] = (cutEdges[ie]);                  cutEdgesMap[(Edge::Pair)*(cutEdges[ie])] = (cutEdges[ie]);
468    
469              //! Addes a triangle to the list.              //! Addes a triangle to the list.
470              (*localTriangleListesEdges[objectNumber]).push_front(              (*localTriangleListesEdges[objectNumber]).push_front(TinyVector<3,Edge::Pair>());
                                                           TinyVector<3,TinyVector<2,TinyVector<2,Edge::Pair> > >());  
471    
472              cellList.push_front(&currentCell);           //   cellList.push_front(&currentCell);
473    
474              std::vector<TinyVector<3> > V(cutEdges.size());              std::vector<TinyVector<3> > V(cutEdges.size());
475              std::map<Edge::Pair, const Edge*>::iterator currentCutEdge=cutEdgesMap.begin();              std::map<Edge::Pair, const Edge*>::iterator currentCutEdge=cutEdgesMap.begin();
# Line 565  __createListIntersection(const Shape& S, Line 485  __createListIntersection(const Shape& S,
485              bool swapEdge = ((V[0]^V[1])*V[2] > 0);              bool swapEdge = ((V[0]^V[1])*V[2] > 0);
486    
487              //! Sets the triangle vertices (actually just store the edge where it lives).              //! Sets the triangle vertices (actually just store the edge where it lives).
488              TriangleVerticesEdgesList::iterator currentTriangleEdges              EdgePairList::iterator currentTriangleEdges
489                  = (*localTriangleListesEdges[objectNumber]).begin();                  = (*localTriangleListesEdges[objectNumber]).begin();
490                                            
491              currentCutEdge = cutEdgesMap.begin();              currentCutEdge = cutEdgesMap.begin();
492              if (!swapEdge) {              if (!swapEdge) {
493                  for(size_t i=0; i<3; ++i) {                  for(size_t i=0; i<3; ++i) {
494                      Edge::Pair e=(currentCutEdge->first);                      Edge::Pair e=(currentCutEdge->first);
495                      TinyVector<2,TinyVector<2,Edge::Pair> > pe;                      (*currentTriangleEdges)[i] = e;
                     pe[0][0]=e;  
                     pe[1][0]=e;  
                     pe[0][1]=e;  
                     pe[1][1]=e;  
                     (*currentTriangleEdges)[i] = pe;  
496                      ++currentCutEdge;                      ++currentCutEdge;
497                  }                  }
498              }              }
499              else {              else {
500                  for(int i=2; i>=0; --i) {                  for(int i=2; i>=0; --i) {
501                      Edge::Pair e=(currentCutEdge->first);                      Edge::Pair e=(currentCutEdge->first);
502                      TinyVector<2,Edge::Pair> pe;                      (*currentTriangleEdges)[i] = e;
                     pe[0]=e;  
                     pe[1]=e;  
                     (*currentTriangleEdges)[i] = pe;  
503                      ++currentCutEdge;                      ++currentCutEdge;
504                  }                  }
505              }              }
# Line 598  __createListIntersection(const Shape& S, Line 510  __createListIntersection(const Shape& S,
510              }              }
511    
512              //! Addes a triangle to the list.              //! Addes a triangle to the list.
513              (*localTriangleListesEdges[objectNumber]).push_front(              (*localTriangleListesEdges[objectNumber]).push_front( TinyVector<3,Edge::Pair>());
                                                           TinyVector<3,TinyVector<2,TinyVector<2,Edge::Pair> > >());  
514    
515              cellList.push_front(&currentCell);              //cellList.push_front(&currentCell);
516    
517              currentTriangleEdges              currentTriangleEdges
518                  = (*localTriangleListesEdges[objectNumber]).begin();                  = (*localTriangleListesEdges[objectNumber]).begin();
# Line 613  __createListIntersection(const Shape& S, Line 524  __createListIntersection(const Shape& S,
524                      ++currentCutEdge;                      ++currentCutEdge;
525                    //  (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);                    //  (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
526                      Edge::Pair e=(currentCutEdge->first);                      Edge::Pair e=(currentCutEdge->first);
527                      TinyVector<2,TinyVector<2,Edge::Pair> > pe;                      (*currentTriangleEdges)[i] = e;
                     pe[0][0]=e;  
                     pe[1][0]=e;  
                     pe[0][1]=e;  
                     pe[1][1]=e;  
                     (*currentTriangleEdges)[i] = pe;  
528                  }                  }
529              }              }
530              else {              else {
# Line 626  __createListIntersection(const Shape& S, Line 532  __createListIntersection(const Shape& S,
532                      ++currentCutEdge;                      ++currentCutEdge;
533                     // (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);                     // (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
534                      Edge::Pair e=(currentCutEdge->first);                      Edge::Pair e=(currentCutEdge->first);
535                      TinyVector<2,TinyVector<2,Edge::Pair> > pe;                      (*currentTriangleEdges)[i] = e;
                     pe[0][0]=e;  
                     pe[1][0]=e;  
                     pe[0][1]=e;  
                     pe[1][1]=e;  
                     (*currentTriangleEdges)[i] = pe;  
536                  }                  }
537              }              }
538              break;              break;
# Line 674  __createCutEdges(const Scene& scene, Line 575  __createCutEdges(const Scene& scene,
575    
576  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
577  __calculatePointsIntersection(const Scene& scene,  __calculatePointsIntersection(const Scene& scene,
                               size_t nobjectTotreat,  
578                                size_t no,                                size_t no,
579                                size_t no2) {                                size_t no2) {
580      size_t nobject=scene.nbObjects();      size_t nobject=scene.nbObjects();
581    
582      if(no< nobject and no2<nobject) {      if(no< nobject and no2<nobject) {
583            ListOfHexaToTriangle::iterator hexait=listOfHexaToTriangle.begin();
584            ListOfHexaToTriangle::iterator hexait0=hexait;
585            ++hexait;
586            ListOfHexaToTriangle::iterator hexait1=hexait;
587            MapCellTriangle hexaToTriangle0=*(*hexait0);
588            MapCellTriangle hexaToTriangle1=*(*hexait1);
589    
590          for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)          for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)
591          {          {
592              for(SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::const_iterator currentTriangle1              MapCellTriangle::iterator it0=hexaToTriangle0.find(*i);
593           =  (*triangleListesIntersectionEdgesPerHexa[*i][no]).begin();              MapCellTriangle::iterator it1=hexaToTriangle1.find(*i);
594           currentTriangle1 !=  (*triangleListesIntersectionEdgesPerHexa[*i][no]).end();              if(it0!=hexaToTriangle0.end() and it1!=hexaToTriangle1.end()) {
595           ++currentTriangle1) {                     std::list<Triangle*> listTriangle0=(*it0).second;
596                  for(SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::const_iterator currentTriangle                     std::list<Triangle*> listTriangle1=(*it1).second;
597        =  (*triangleListesIntersectionEdgesPerHexa[*i][no2]).begin();              for(std::list<Triangle*>::iterator jt0=listTriangle0.begin() ;
598        currentTriangle !=  (*triangleListesIntersectionEdgesPerHexa[*i][no2]).end();           jt0!=listTriangle0.end() ; ++jt0) {
599        ++currentTriangle) {                      for(std::list<Triangle*>::iterator jt1=listTriangle1.begin() ;
600                      __createPointIntersection(no,no2,*i,currentTriangle1,currentTriangle);             jt1!=listTriangle1.end() ; ++jt1) {
601                  }                          __createPointIntersection(no,no2,jt0,jt1);
             }  
         }  
     }  
     else {  
         size_t n0,n1;  
         if(no< nobject) {  
             n0=no;  
             n1=no2-nobject;  
             for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)  
             {  
                 if((*triangleListesIntersectionEdgesPerHexa[*i][n0]).size()!=0 and  
      (*triangleListesIntersectionPerHexa[*i][n1]).size()!=0) {  
                     for(SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::const_iterator currentTriangle1  
           =  (*triangleListesIntersectionEdgesPerHexa[*i][n0]).begin();  
           currentTriangle1 !=  (*triangleListesIntersectionEdgesPerHexa[*i][n0]).end();  
           ++currentTriangle1) {  
                         for(SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::const_iterator currentTriangle  
        =  (*triangleListesIntersectionPerHexa[*i][n1]).begin();  
        currentTriangle !=  (*triangleListesIntersectionEdgesPerHexa[*i][n1]).end();  
        ++currentTriangle) {  
                             __createPointIntersection(n0,n1,*i,currentTriangle1,currentTriangle);  
                         }  
                     }  
                 }  
             }  
         }  
         else if(no2<nobject) {  
             n1=no2;  
             n0=no-nobject;  
             //ffout(0)<<"******************* "<<n0<<" "<<n1<<"\n";  
             for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)  
             {  
                 if((*triangleListesIntersectionEdgesPerHexa[*i][n1]).size()!=0 and  
      (*triangleListesIntersectionPerHexa[*i][n0]).size()!=0) {  
                     for(SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::const_iterator currentTriangle1  
           =  (*triangleListesIntersectionPerHexa[*i][n0]).begin();  
           currentTriangle1 !=  (*triangleListesIntersectionPerHexa[*i][n0]).end();  
           ++currentTriangle1) {  
                         for(SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::const_iterator currentTriangle  
        =  (*triangleListesIntersectionEdgesPerHexa[*i][n1]).begin();  
        currentTriangle !=  (*triangleListesIntersectionEdgesPerHexa[*i][n1]).end();  
        ++currentTriangle) {  
                             __createPointIntersection(n0,n1,*i,currentTriangle1,currentTriangle);  
                         }  
602                      }                      }
603                  }                  }
604                  else {                  
                   //  ffout(0)<<"pas inter dans hexa\n";  
                 }  
605              }              }
606          }              else {
607          else {                  //ffout(0)<<"pas d'intersection dans l'hexa\n";
             n1=no2-nobject;  
             n0=no-nobject;  
             for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)  
             {  
                 if((*triangleListesIntersectionPerHexa[*i][n1]).size()!=0 and  
      (*triangleListesIntersectionPerHexa[*i][n0]).size()!=0) {  
                     for(SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::const_iterator  
           currentTriangle1 =  (*triangleListesIntersectionPerHexa[*i][n0]).begin();  
           currentTriangle1 !=  (*triangleListesIntersectionPerHexa[*i][n0]).end();  
           ++currentTriangle1) {  
                         for(SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::const_iterator  
        currentTriangle =  (*triangleListesIntersectionPerHexa[*i][n1]).begin();  
        currentTriangle !=  (*triangleListesIntersectionPerHexa[*i][n1]).end();  
        ++currentTriangle) {  
                             __createPointIntersection(n0,n1,*i,currentTriangle1,currentTriangle);  
                         }  
                     }  
                 }  
                 else {  
                     //  ffout(0)<<"pas inter dans hexa\n";  
                 }  
608              }              }
609          }          }
610      }      }
# Line 774  __calculateListesIntersection(const Scen Line 615  __calculateListesIntersection(const Scen
615                                ObjectReferences& otherReferences) {                                ObjectReferences& otherReferences) {
616      size_t nobject = scene.nbObjects();      size_t nobject = scene.nbObjects();
617      size_t nobjectToTreat = otherReferences.size();      size_t nobjectToTreat = otherReferences.size();
618        TinyVector<2,size_t> objpair;
619        objpair[0]=nobject;
620        objpair[1]=nobject;
621        
622    
623      for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)      for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)
624      {      {
# Line 786  __calculateListesIntersection(const Scen Line 631  __calculateListesIntersection(const Scen
631          triangleListesIntersectionEdgesPerHexa[*i].resize(nobject);          triangleListesIntersectionEdgesPerHexa[*i].resize(nobject);
632          for(size_t no=0 ; no<nobjectToTreat ; ++no) {          for(size_t no=0 ; no<nobjectToTreat ; ++no) {
633              triangleListesIntersectionEdgesPerHexa[*i][no]=              triangleListesIntersectionEdgesPerHexa[*i][no]=
634              new SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList;              new EdgePairList;
635          }          }
636    
637          for (std::vector<TetrahedronByEdges>::iterator currentT = T.begin();          for (std::vector<TetrahedronByEdges>::iterator currentT = T.begin();
# Line 817  __calculateListesIntersection(const Scen Line 662  __calculateListesIntersection(const Scen
662          }          }
663    
664          for(size_t no=0 ; no<nobjectToTreat ; ++no) {          for(size_t no=0 ; no<nobjectToTreat ; ++no) {
665              for (TriangleVerticesEdgesList::iterator t = (*localTriangleListesEdges[no]).begin();              for (EdgePairList::iterator t = (*localTriangleListesEdges[no]).begin();
666            t != (*localTriangleListesEdges[no]).end(); ++t) {            t != (*localTriangleListesEdges[no]).end(); ++t) {
                 (*triangleListesIntersectionEdges[no]).push_front(*t);  
                 (*localTriangleList2Edges[no]).push_front(&(*(*triangleListesIntersectionEdges[no]).begin()));  
667                  (*triangleListesIntersectionEdgesPerHexa[*i][no]).push_front(*t);                  (*triangleListesIntersectionEdgesPerHexa[*i][no]).push_front(*t);
668    
669              }              }
# Line 828  __calculateListesIntersection(const Scen Line 671  __calculateListesIntersection(const Scen
671    
672          for(size_t no=0 ; no<nobject ; ++no) {          for(size_t no=0 ; no<nobject ; ++no) {
673              (*localTriangleListesEdges[no]).clear();              (*localTriangleListesEdges[no]).clear();
             (*localTriangleList2Edges[no]).clear();  
674          }          }
675      }      }
676    
677  }  }
678    
679  void SurfaceMeshGenerator::Internals::printpoint(size_t numobject,  bool SurfaceMeshGenerator::Internals::triangleplat(const Vertex& V0 ,
680                                                   const TinyVector<2,TinyVector<2,Edge::Pair> >& e0 , TriangleCut K) {                                                     const Vertex& V1,
681      Vertex V0;                                                     const Vertex& V2 ) {
682      if(e0[0][0]==e0[1][0])      real_t eps=1e-6;
683          V0=(*(*verticesListesIntersection[numobject])[e0[0][0]]);      real_t normV0V1=Norm(V0.Vector3()-V1.Vector3());
684      else {      real_t normV0V2=Norm(V0.Vector3()-V2.Vector3());
685          std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > p0;      real_t normV1V2=Norm(V2.Vector3()-V1.Vector3());
686          p0.first=e0;      if(normV0V1>eps and normV0V2>eps and normV1V2>eps) {
687          TinyVector<2,size_t> num;          return false;
         num[0]=K.numobject;  
         num[1]=K.otherobject;  
         p0.second=num;  
         SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator  
             ip0=listVerticesIntersectionPerEdges.find(p0);  
         if(ip0!=listVerticesIntersectionPerEdges.end())  
             V0 =*(*ip0).second;  
         else  
             ffout(0)<<"point exi pas\n";  
     }  
     ffout(0)<<V0.Vector3()[0]<<" "<<V0.Vector3()[1]<<" "<<V0.Vector3()[2]<<"\n";  
 }  
   
 bool SurfaceMeshGenerator::Internals::triangleplat(size_t numobject,const TinyVector<2,Edge::Pair>& e0 ,  
                                                    const TinyVector<2,Edge::Pair>& e1,  
                                                    const TinyVector<2,Edge::Pair>& e2 ) {  
   
   
     Vertex V0,V1,V2;  
     real_t eps=1e-13;  
     if(e0[0]==e0[1] and e1[0]==e1[1] and e2[0]==e2[1]) {  
         V0=(*(*verticesListesIntersection[numobject])[e0[0]]);  
         V1=(*(*verticesListesIntersection[numobject])[e1[0]]);  
         V2=(*(*verticesListesIntersection[numobject])[e2[0]]);  
         real_t normV0V1=Norm(V0.Vector3()-V1.Vector3());  
         real_t normV0V2=Norm(V0.Vector3()-V2.Vector3());  
         real_t normV1V2=Norm(V2.Vector3()-V1.Vector3());  
         if(normV0V1>eps and normV0V2>eps and normV1V2>eps) {  
             return false;  
         }  
         else  
             return true;  
688      }      }
689      else      else
690          return false;          return true;
691  }  }
692    
693  void SurfaceMeshGenerator::Internals::plotliste(size_t nobject) {  void SurfaceMeshGenerator::Internals::plotliste() {
694    
695  std::ofstream sortie2("liste2.mesh");      std::ofstream sortie2("liste2.mesh");
696    
697      sortie2<<"MeshVersionFormatted 1"<<" \n";      sortie2<<"MeshVersionFormatted 1"<<" \n";
698      sortie2<<" \n";      sortie2<<" \n";
# Line 890  std::ofstream sortie2("liste2.mesh"); Line 700  std::ofstream sortie2("liste2.mesh");
700      sortie2<<"3"<<" \n";      sortie2<<"3"<<" \n";
701      sortie2<<" \n";      sortie2<<" \n";
702      sortie2<<"Vertices"<<" \n";      sortie2<<"Vertices"<<" \n";
703      size_t nn=0;      size_t nn=0,n0=0,n1=0,compt=0;
704      for(size_t no=0;no<nobject;++no)      for(std::list<ReferenceCounting<Vector<Triangle> > >::iterator i=(listOfTriangleListIntersection).begin();
705          nn+=(*triangleListesIntersectionEdges[no]).size();          i!=(listOfTriangleListIntersection).end();++i) {
706            nn+=(*(*i)).size();
707            if(compt==0)
708                n0=nn;
709            else if(compt==1)
710                n1=nn;
711            ++compt;
712        }
713      sortie2<<3*nn<<"\n";      sortie2<<3*nn<<"\n";
714      for(size_t nz=0;nz<nobject;++nz) {      for(std::list<ReferenceCounting<Vector<Triangle> > >::iterator i=(listOfTriangleListIntersection).begin();
715          for(std::list<TinyVector<3,TinyVector<2,TinyVector<2,Edge::Pair> > > >::iterator currentTriangle          i!=(listOfTriangleListIntersection).end();++i) {
716       = (*triangleListesIntersectionEdges[nz]).begin();          Vector<Triangle> triangleListIntersection=*(*i);
717       currentTriangle != (*triangleListesIntersectionEdges[nz]).end();          for(size_t size=0 ; size<triangleListIntersection.size() ; ++size) {
718       ++currentTriangle) {              const Triangle& T=triangleListIntersection[size];
719              const Edge& e0=(*currentTriangle)[0][0][0];              const Vertex& V0=T(0);
720              const Edge& e1=(*currentTriangle)[1][0][0];              const Vertex& V1=T(1);
721              const Edge& e2=(*currentTriangle)[2][0][0];              const Vertex& V2=T(2);
             const Vertex& V0=(*(*verticesListesIntersection[nz])[e0]);  
             const Vertex& V1=(*(*verticesListesIntersection[nz])[e1]);  
             const Vertex& V2=(*(*verticesListesIntersection[nz])[e2]);  
722    
723              sortie2<<V0.Vector3()[0]<<" "<<V0.Vector3()[1]<<" "<<V0.Vector3()[2]<<" 0 \n";              sortie2<<V0.Vector3()[0]<<" "<<V0.Vector3()[1]<<" "<<V0.Vector3()[2]<<" 0 \n";
724              sortie2<<V1.Vector3()[0]<<" "<<V1.Vector3()[1]<<" "<<V1.Vector3()[2]<<" 0 \n";              sortie2<<V1.Vector3()[0]<<" "<<V1.Vector3()[1]<<" "<<V1.Vector3()[2]<<" 0 \n";
# Line 917  std::ofstream sortie2("liste2.mesh"); Line 731  std::ofstream sortie2("liste2.mesh");
731    
732      for(size_t nj=0;nj<nn;++nj) {      for(size_t nj=0;nj<nn;++nj) {
733          sortie2<<3*nj+1<<" "<<3*nj+2<<" "<<3*nj+3;          sortie2<<3*nj+1<<" "<<3*nj+2<<" "<<3*nj+3;
734          if(nj<(*triangleListesIntersectionEdges[0]).size())          if(nj<n0)
735              sortie2<<" 2 \n";              sortie2<<" 2 \n";
736          else if(nj<(*triangleListesIntersectionEdges[0]).size()+(*triangleListesIntersectionEdges[1]).size())          else if(nj<n1)
737              sortie2<<" 3 \n";              sortie2<<" 3 \n";
738          else          else
739              sortie2<<" 4 \n";              sortie2<<" 4 \n";
# Line 927  std::ofstream sortie2("liste2.mesh"); Line 741  std::ofstream sortie2("liste2.mesh");
741      sortie2<<"End"<<" \n";      sortie2<<"End"<<" \n";
742  }  }
743    
 void SurfaceMeshGenerator::Internals::plotlistenew(size_t nobject,size_t nobjectToTreat) {  
   
     std::ofstream sortie2("newinter.mesh");  
   
     sortie2<<"MeshVersionFormatted 1"<<" \n";  
     sortie2<<" \n";  
     sortie2<<"Dimension"<<" \n";  
     sortie2<<"3"<<" \n";  
     sortie2<<" \n";  
     sortie2<<"Vertices"<<" \n";  
     size_t nn=0;  
     for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)  
     {  
         for(size_t no=0;no<nobjectToTreat;++no) {  
             nn+=(*triangleListesIntersectionPerHexa[*i][no]).size();  
         }  
     }  
     sortie2<<3*nn<<"\n";  
   
     for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)  
     {  
         for(size_t nz=0;nz<nobjectToTreat;++nz) {  
             for(TriangleVerticesEdgesList::iterator currentTriangle  
          = (*triangleListesIntersectionPerHexa[*i][nz]).begin();  
          currentTriangle != (*triangleListesIntersectionPerHexa[*i][nz]).end();  
          ++currentTriangle) {  
                 const TinyVector<2,TinyVector<2,Edge::Pair> >& e0=(*currentTriangle)[0];  
                 const TinyVector<2,TinyVector<2,Edge::Pair> >& e1=(*currentTriangle)[1];  
                 const TinyVector<2,TinyVector<2,Edge::Pair> >& e2=(*currentTriangle)[2];  
                 Vertex V0=(*__determinateVertex(nz,e0));  
                 Vertex V1=(*__determinateVertex(nz,e1));  
                 Vertex V2=(*__determinateVertex(nz,e2));  
                 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;  
         sortie2<<" 4 \n";  
     }  
     sortie2<<"End"<<" \n";  
 }  
   
   
 Vertex* SurfaceMeshGenerator::Internals::__determinateVertex(size_t numobject,  
                                                             const TinyVector<2,TinyVector<2,Edge::Pair> >& e0) {  
     Vertex* V0;  
     if(e0[0][0]==e0[0][1] and e0[1][0]==e0[1][1] and e0[1][0]==e0[0][1]) {  
         //ffout(0)<<"I search in verticesListesIntersection[numobject]\n";  
         V0=((*verticesListesIntersection[numobject])[e0[0][0]]);  
     }  
     else {  
         std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > p0;  
         p0.first=e0;  
         TinyVector<2,size_t> num;  
         num[0]=0;  
         num[1]=0;  
         p0.second=num;  
         SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator  
             ip0=listVerticesIntersectionPerEdges.find(p0);  
         if(ip0!=listVerticesIntersectionPerEdges.end())  
             V0 =(*ip0).second;  
         else  
             ffout(0)<<"point exi pas\n";  
     }  
     return V0;  
 }  
   
744  size_t SurfaceMeshGenerator::Internals::  size_t SurfaceMeshGenerator::Internals::
745  __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 ) {
746    
# Line 1034  __IsInTriangle(const TinyVector<2,size_t Line 774  __IsInTriangle(const TinyVector<2,size_t
774    
775  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
776  __InTriangle(const Vertex& P,  __InTriangle(const Vertex& P,
777               const TinyVector<3, Vertex* > & V,               TinyVector<3,const Vertex* > & V,
778               size_t& compt) {               size_t& compt) {
779    
780      compt=0;      compt=0;
781        TinyVector<3,TinyVector<3,real_t> >  triangleVector;
     TinyVector<3,TinyVector<3> >  triangleVector;  
782      triangleVector[0]=-(*V[2]).Vector3() + (*V[0]).Vector3();      triangleVector[0]=-(*V[2]).Vector3() + (*V[0]).Vector3();
783      triangleVector[1]= -(*V[2]).Vector3() + (*V[1]).Vector3();      triangleVector[1]= -(*V[2]).Vector3() + (*V[1]).Vector3();
784      triangleVector[2]=P.Vector3()-(*V[2]).Vector3();      triangleVector[2]=P.Vector3()-(*V[2]).Vector3();
# Line 1063  __InTriangle(const Vertex& P, Line 802  __InTriangle(const Vertex& P,
802  }  }
803    
804  TinyVector<4, real_t>  SurfaceMeshGenerator::Internals::  TinyVector<4, real_t>  SurfaceMeshGenerator::Internals::
805  __EquationPlan(const TinyVector<3, Vertex* > & V){  __EquationPlan(const TinyVector<3,const Vertex *> & V){
806    
807      //equation of plan      //equation of plan
808    
# Line 1087  __EquationPlan(const TinyVector<3, Verte Line 826  __EquationPlan(const TinyVector<3, Verte
826    
827  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
828  __DeterminatePoint(const size_t& ncase,  __DeterminatePoint(const size_t& ncase,
829                     const TinyVector<2,Edge::Pair> & eo0,                     const Vertex* & Vo0,
830                     const TinyVector<2,Edge::Pair> & eo1,                     const Vertex* & Vo1,
831                     const Vertex & Vo0,                     TinyVector<3,const Vertex* > & V,
                    const Vertex & Vo1,  
                    const TinyVector<3, Vertex* > & V,  
832                     IntersectionPoints& Pinter) {                     IntersectionPoints& Pinter) {
833      std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > p1,p0;      std::pair<const Vertex*,const Vertex*> pV(Vo0,Vo1);
834      TinyVector<2,TinyVector<2,Edge::Pair> > pt1,pt0;      std::pair<const Vertex*,const Vertex*> pV2(Vo1,Vo0);
     pt0[0]=eo0;  
     pt0[1]=eo1;  
     pt1[0]=eo1;  
     pt1[1]=eo0;  
     p1.first=pt1;  
     p0.first=pt0;  
     TinyVector<2,size_t> nuum;  
     nuum[0]=Pinter.numpair[0];  
     nuum[1]=Pinter.numpair[1];  
     p1.second=nuum;  
     p0.second=nuum;  
835    
836      real_t epsilon=1e-10;      real_t epsilon=1e-10;
837      real_t pscalar,k;      real_t pscalar,k;
# Line 1118  __DeterminatePoint(const size_t& ncase, Line 844  __DeterminatePoint(const size_t& ncase,
844      coeff[2]=coeffo[2];      coeff[2]=coeffo[2];
845      real_t d=coeffo[3];      real_t d=coeffo[3];
846    
847      TinyVector<3> Vo01 = Vo1.Vector3() - Vo0.Vector3();      TinyVector<3> Vo01 = (*Vo1).Vector3() - (*Vo0).Vector3();
848      TinyVector<3> V01 = (*V[1]).Vector3() - (*V[0]).Vector3();      TinyVector<3> V01 = (*V[1]).Vector3() - (*V[0]).Vector3();
849      TinyVector<3> V12 = (*V[2]).Vector3() - (*V[1]).Vector3();      TinyVector<3> V12 = (*V[2]).Vector3() - (*V[1]).Vector3();
850      TinyVector<3> V02 = (*V[2]).Vector3() - (*V[0]).Vector3();      TinyVector<3> V02 = (*V[2]).Vector3() - (*V[0]).Vector3();
# Line 1127  __DeterminatePoint(const size_t& ncase, Line 853  __DeterminatePoint(const size_t& ncase,
853      pscalar=coeff*Vo01;      pscalar=coeff*Vo01;
854      num=0;      num=0;
855      if(std::abs(pscalar)>1e-10) {      if(std::abs(pscalar)>1e-10) {
856          k=-(coeff*Vo0.Vector3()+d)/pscalar;          k=-(coeff*(*Vo0).Vector3()+d)/pscalar;
857          P=Vo0+k*Vo01;          P=(*Vo0)+k*Vo01;
858                    
859          TinyVector<3> V0P=P.Vector3()-(*V[0]).Vector3();          TinyVector<3> V0P=P.Vector3()-(*V[0]).Vector3();
860          TinyVector<3> V1P=P.Vector3()-(*V[1]).Vector3();          TinyVector<3> V1P=P.Vector3()-(*V[1]).Vector3();
861          TinyVector<3> V2P=P.Vector3()-(*V[2]).Vector3();          TinyVector<3> V2P=P.Vector3()-(*V[2]).Vector3();
862    
863          real_t norm2=Norm(P.Vector3()-Vo0.Vector3());          real_t norm2=Norm(P.Vector3()-(*Vo0).Vector3());
864          real_t norm3=Norm(P.Vector3()-Vo1.Vector3());          real_t norm3=Norm(P.Vector3()-(*Vo1).Vector3());
865          real_t inV01=Norm(V01^V0P);          real_t inV01=Norm(V01^V0P);
866          real_t inV12=Norm(V12^V1P);          real_t inV12=Norm(V12^V1P);
867          real_t inV02=Norm(V02^V2P);          real_t inV02=Norm(V02^V2P);
868    
869          if(norm2>epsilon and norm3>epsilon and std::abs((norm2+norm3)-norm1) < epsilon ) {          if(/*norm2>epsilon and norm3>epsilon and*/ std::abs((norm2+norm3)-norm1) < epsilon ) {
870              if(ncase<3) {              if(ncase<3) {
871                  __InTriangle(P,V,num);                  __InTriangle(P,V,num);
872                  if(num==3) {                  if(num==3) {
873                      if(p1.first[0][0]!=p1.first[0][1] or p1.first[1][0]!=p1.first[1][1])                      if(!(norm2>epsilon)) {
874                          ffout(0)<<"un new intersect============\n";                          Pinter.addPoint(Vo0);
875                      if (listVerticesIntersectionPerEdges.find(p1) == listVerticesIntersectionPerEdges.end()                      }
876            and listVerticesIntersectionPerEdges.find(p0) == listVerticesIntersectionPerEdges.end()) {                      else if(!(norm3>epsilon)) {
877                          listVerticesIntersectionPerEdges[p0]=new Vertex (P.Vector3());                          Pinter.addPoint(Vo1);
878                          Pinter.addPoint(p0);                      }
879                      }                      else if(listVertexIntersectionNew.find(pV)==listVertexIntersectionNew.end()
880                      else if(listVerticesIntersectionPerEdges.find(p1) != listVerticesIntersectionPerEdges.end()){           and listVertexIntersectionNew.find(pV2)==listVertexIntersectionNew.end()) {
881                          Vertex Ptemp=(*listVerticesIntersectionPerEdges[p1]);                          listVertexIntersectionNew[pV]=new Vertex (P.Vector3());
882                          Pinter.addPoint(p1);                          Pinter.addPoint(listVertexIntersectionNew[pV]);
                     }  
                     else if(listVerticesIntersectionPerEdges.find(p0) != listVerticesIntersectionPerEdges.end()){  
                         Vertex Ptemp=(*listVerticesIntersectionPerEdges[p0]);  
                         Pinter.addPoint(p0);  
883                      }                      }
884                      else {                      else {
885                          ffout(0)<<"pbbbbb\n";                          if(listVertexIntersectionNew.find(pV)!=listVertexIntersectionNew.end())
886                                Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second);
887                            else
888                                 Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second);
889                      }                      }
                     //if(norm2>epsilon and norm3>epsilon)  
890                      Pinter.cutTriangle1[ncase]=Pinter.size();                      Pinter.cutTriangle1[ncase]=Pinter.size();
891                        
892                  }                  }
893              }              }
894              else if(inV01> epsilon and inV02> epsilon and inV12> epsilon ) {              else if(inV01> epsilon and inV02> epsilon and inV12> epsilon ) {
895                  __InTriangle(P,V,num);                  __InTriangle(P,V,num);
896                  if(num==3) {                  if(num==3) {
897                      if(p1.first[0][0]!=p1.first[0][1] or p1.first[1][0]!=p1.first[1][1])                      //ffout(0)<<"ajout du point cas 2\n";
898                          ffout(0)<<"un new intersect============\n";                      //ffout(0)<<P.Vector3()[0]<<" "<<P.Vector3()[1]<<" "<<P.Vector3()[2]<<" \n";
899                      if (listVerticesIntersectionPerEdges.find(p1) == listVerticesIntersectionPerEdges.end()                      //listVertexIntersectionNew.push_back(new Vertex (P.Vector3()));
900            and listVerticesIntersectionPerEdges.find(p0) == listVerticesIntersectionPerEdges.end()) {                      PIntersection PIn;
901                          listVerticesIntersectionPerEdges[p0]=new Vertex (P.Vector3());                      if(!(norm2>epsilon)) {
902                          Pinter.addPoint(p0);                          Pinter.addPoint(Vo0);
                     }  
                     else if(listVerticesIntersectionPerEdges.find(p1) != listVerticesIntersectionPerEdges.end()){  
                         Vertex Ptemp=(*listVerticesIntersectionPerEdges[p1]);  
                         if(Ptemp.Vector3()==P.Vector3())  
                             Pinter.addPoint(p1);  
                         else {  
                             TinyVector<2,size_t> temp;  
                             temp[0]=Pinter.numpair[1];  
                             temp[1]=Pinter.numpair[0];  
                             p1.second=temp;  
                             listVerticesIntersectionPerEdges[p1]=new Vertex (P.Vector3());  
                             Pinter.addPoint(p1);  
                         }  
903                      }                      }
904                      else if(listVerticesIntersectionPerEdges.find(p0) != listVerticesIntersectionPerEdges.end()){                      else if(!(norm3>epsilon)) {
905                          Vertex Ptemp=(*listVerticesIntersectionPerEdges[p0]);                          Pinter.addPoint(Vo1);
906                          if(Ptemp.Vector3()==P.Vector3())                      }
907                              Pinter.addPoint(p0);                      else if(listVertexIntersectionNew.find(pV)==listVertexIntersectionNew.end()
908                          else {           and listVertexIntersectionNew.find(pV2)==listVertexIntersectionNew.end()) {
909                              TinyVector<2,size_t> temp;                          listVertexIntersectionNew[pV]=new Vertex (P.Vector3());
910                              temp[0]=Pinter.numpair[1];                          Pinter.addPoint(listVertexIntersectionNew[pV]);
                             temp[1]=Pinter.numpair[0];  
                             p0.second=temp;  
                             listVerticesIntersectionPerEdges[p0]=new Vertex (P.Vector3());  
                             Pinter.addPoint(p0);  
                         }  
   
911                      }                      }
912                      else {                      else {
913                          ffout(0)<<"pbbbbb\n";                          if(listVertexIntersectionNew.find(pV)!=listVertexIntersectionNew.end())
914                                Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second);
915                            else
916                                Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second);    
917                      }                      }
918                       // Pinter.addPoint(listVertexIntersectionNew[listVertexIntersectionNew.size()-1]);
919                  }                  }
920              }              }
921          }          }
# Line 1214  __DeterminatePoint(const size_t& ncase, Line 924  __DeterminatePoint(const size_t& ncase,
924    
925  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
926  __CalculPoint(const size_t& objectNumber,const size_t& numobject,  __CalculPoint(const size_t& objectNumber,const size_t& numobject,
927                SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::const_iterator & currentTriangle,                std::list<Triangle*>::iterator & currentTriangle,
928                SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::const_iterator & currentTriangle1,                std::list<Triangle*>::iterator & currentTriangle1,
929                IntersectionPoints& Pinter) {                IntersectionPoints& Pinter) {
930    
931      //ffout(0) <<"here ok\n";      //ffout(0) <<"here ok\n";
932      size_t num1=objectNumber,num0=numobject;      size_t num1=objectNumber,num0=numobject;
933            
934      TinyVector<3, Vertex* > V012;      TinyVector<3,const Vertex * > V012;
935      const TinyVector<2,TinyVector<2,Edge::Pair> >& e0=((*currentTriangle1))[0];      const Triangle& T0=*(*currentTriangle1);
936      const TinyVector<2,TinyVector<2,Edge::Pair> >& e1=((*currentTriangle1))[1];      V012[0] = &T0(0);
937      const TinyVector<2,TinyVector<2,Edge::Pair> >& e2=((*currentTriangle1))[2];      V012[1] = &T0(1);
938      V012[0] = __determinateVertex(num0,e0);      V012[2] = &T0(2);
939      V012[1] = __determinateVertex(num0,e1);      TinyVector<3,const Vertex* > Vo012;
940      V012[2] = __determinateVertex(num0,e2);      Triangle& To0=*(*currentTriangle);
941        Vo012[0] = &To0(0);
942      TinyVector<3, Vertex* > Vo012;      Vo012[1] = &To0(1);
943      const TinyVector<2,TinyVector<2,Edge::Pair> >& eo0=((*currentTriangle))[0];      Vo012[2] = &To0(2);
     const TinyVector<2,TinyVector<2,Edge::Pair> >& eo1=((*currentTriangle))[1];  
     const TinyVector<2,TinyVector<2,Edge::Pair> >& eo2=((*currentTriangle))[2];  
     Vo012[0] = __determinateVertex(num1,eo0);  
     Vo012[1] = __determinateVertex(num1,eo1);  
     Vo012[2] = __determinateVertex(num1,eo2);  
   
944      //for each edge of the triangle V0V1V2      //for each edge of the triangle V0V1V2
945      Pinter.numpair[0]=num0;      Pinter.numpair[0]=num0;
946      Pinter.numpair[1]=num1;      Pinter.numpair[1]=num1;
947      TinyVector<2,Edge::Pair> ee0,ee1,ee2,eeo0,eeo1,eeo2;  
948      ee0[0]=e0[0][0];      __DeterminatePoint(0,(V012[0]),(V012[1]),Vo012,Pinter);
949      ee0[1]=e0[1][0];      __DeterminatePoint(1,(V012[0]),(V012[2]),Vo012,Pinter);
950      ee1[0]=e1[0][0];      __DeterminatePoint(2,(V012[1]),(V012[2]),Vo012,Pinter);
951      ee1[1]=e1[1][0];  
952      ee2[0]=e2[0][0];      __DeterminatePoint(3,(Vo012[0]),(Vo012[1]),V012,Pinter);
953      ee2[1]=e2[1][0];      __DeterminatePoint(3,(Vo012[0]),(Vo012[2]),V012,Pinter);
954      eeo0[0]=eo0[0][0];      __DeterminatePoint(3,(Vo012[1]),(Vo012[2]),V012,Pinter);
955      eeo0[1]=eo0[1][0];  
     eeo1[0]=eo1[0][0];  
     eeo1[1]=eo1[1][0];  
     eeo2[0]=eo2[0][0];  
     eeo2[1]=eo2[1][0];  
       
     __DeterminatePoint(0,ee0,ee1,(*V012[0]),(*V012[1]),Vo012,Pinter);  
     __DeterminatePoint(1,ee0,ee2,(*V012[0]),(*V012[2]),Vo012,Pinter);  
     __DeterminatePoint(2,ee1,ee2,(*V012[1]),(*V012[2]),Vo012,Pinter);  
   
   
     __DeterminatePoint(3,eeo0,eeo1,(*Vo012[0]),(*Vo012[1]),V012,Pinter);  
     __DeterminatePoint(3,eeo0,eeo2,(*Vo012[0]),(*Vo012[2]),V012,Pinter);  
     __DeterminatePoint(3,eeo1,eeo2,(*Vo012[1]),(*Vo012[2]),V012,Pinter);  
       
956  }  }
957    
958  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
959  __inout(const size_t& numobject,const int & ttn0,const TinyVector<2,Edge::Pair>& e0,  __inout(const size_t& ncase,
960          const TinyVector<2,Edge::Pair>& e1,size_t& ii0,size_t& ii1,TriangleCut& K) {          const size_t& numobject,
961            const int & ttn0,
962      size_t& re0=(*edgesRefCreate[numobject])[e0[0]];          const Vertex*& V0,
963      size_t& re1=(*edgesRefCreate[numobject])[e1[0]];          const Vertex*& V1,
964                TriangleCut& K) {
965    ii0=0;  
966    ii1=0;      size_t& re0=(*edgesRefVertex[numobject])[V0];
967    //if(e0[0]!=e0[1])      size_t& re1=(*edgesRefVertex[numobject])[V1];
968      //  re0=2;  
969    //if(e1[0]!=e1[1])      if(re0==1 and re1==1) {
970      //  re1=2;          ffout(0)<<"--------------- re ==1 1\n";
971              if(ttn0!=0)
972    if(re0==1 and re1==1) {              ffout(0)<<"-- cas indetermine\n";
973        ffout(0)<<"--------------- re ==1 1\n";      }
974        if(ttn0!=0)      else if((re0+re1==1)) {
975            ffout(0)<<"-- cas indetermine\n";          ffout(0)<<"-------------- re ==1 2\n";
976    }          if(ttn0!=0) {
977    else if((re0+re1==1)) {              //ffout(0)<<"-- 1 in 0 out on sait traite\n";
978        ffout(0)<<"-------------- re ==1 2\n";          }
979      if(ttn0!=0) {      }
980        //ffout(0)<<"-- 1 in 0 out on sait traite\n";      else if((re0+re1==3)) {
981        ii0=re0;          ffout(0)<<"------------- re ==1 3\n";
982        ii1=re1;  
983      }          if(ttn0==0) {
984    }              //ffout(0)<<"-- 1 2 pas coupe\n";
985    else if((re0+re1==3)) {              re0=2;
986        ffout(0)<<"------------- re ==1 3\n";              re1=2;
987            }
988      if(ttn0==0) {          else {
989        //ffout(0)<<"-- 1 2 pas coupe\n";              //ffout(0)<<"-- 1 2 on sait traite\n";
990        re0=2;              if(re0==1) {
991        re1=2;                  re0=0;
992                }
993                else {
994                    re1=0;
995                }
996            }
997      }      }
998      else {      else if(re0+re1==4 and ttn0==1) {
999        //ffout(0)<<"-- 1 2 on sait traite\n";          ffout(0)<<"warning edge ref 2 cut \n";
1000        if(re0==1) {          ffout(0)<<(*V0).Vector3()[0]<<" "<<(*V0).Vector3()[1]<<" "<<(*V0).Vector3()[2]<<" \n";
1001          ii0=0;          ffout(0)<<(*V1).Vector3()[0]<<" "<<(*V1).Vector3()[1]<<" "<<(*V1).Vector3()[2]<<" \n";
         ii1=1;  
1002          re0=0;          re0=0;
       }  
       else {  
         ii0=1;  
         ii1=0;  
1003          re1=0;          re1=0;
       }  
1004      }      }
1005    }      else if(re0+re1==2 and ttn0==0) {
1006    else if( ttn0==0 and ( (re0==0 and re1==2) or (re0==2 and re1==0) ) ) {          re0=2;
1007        //ffout(0)<<"new rajout------------\n";          re1=2;
1008        re0=2;      }
1009        re1=2;      //     else
   }  
   else if(re0+re1==4 and ttn0==1) {  
       ffout(0)<<"warning edge ref 2 cut \n";  
       printpoint(numobject,e0,K);  
       printpoint(numobject,e1,K);  
       re0=0;  
       re1=0;  
   }  
   //     else  
1010    //ffout(0)<<"--0 out 0 out rien a faire\n";    //ffout(0)<<"--0 out 0 out rien a faire\n";
1011    
1012  }  }
1013    
1014  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1015  __createPointIntersection(const size_t& numobject,const size_t& objectNumber,  __createPointIntersection(const size_t& numobject,const size_t& objectNumber,
1016                            const Cell* cell,                            std::list<Triangle*>::iterator currentTriangle1,
1017                            SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::const_iterator currentTriangle1,                            std::list<Triangle*>::iterator currentTriangle) {
                           SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::const_iterator currentTriangle) {  
     ObjectPair objectp(numobject,objectNumber);  
1018    
1019        
1020      //ffout(0)<< "createpointintersection "<<numobject<<" "<<objectNumber<<"\n";      //ffout(0)<< "createpointintersection "<<numobject<<" "<<objectNumber<<"\n";
1021      //ffout(0)<<"*************************\n";      //ffout(0)<<"*************************\n";
1022    
1023      PairTriangleEdges pTriangle;      PairTriangleNew pTriangle;
1024      pTriangle.first=&(*currentTriangle1);      pTriangle.first=(*currentTriangle1);
1025      pTriangle.second=&(*currentTriangle);      pTriangle.second=(*currentTriangle);
1026    
1027      IntersectionPoints Pinter;      IntersectionPoints Pinter;
1028            
1029      __CalculPoint(objectNumber,numobject,currentTriangle,currentTriangle1,Pinter);      __CalculPoint(objectNumber,numobject,currentTriangle,currentTriangle1,Pinter);
1030    
1031      TinyVector<3,size_t> pintercut=Pinter.cutTriangle1;      const TinyVector<3,size_t >& pintercut=Pinter.cutTriangle1;
1032    
1033      if(pintercut[0]+pintercut[1]+pintercut[2]!=0) {      if(pintercut[0]+pintercut[1]+pintercut[2]!=0) {
1034          //ffout(0)<<"le tr est cut\n";          //ffout(0)<<"le tr est cut\n";
1035          if((listVertexIntersectionPerPairEdges[objectp]).find(pTriangle)==          if(pairTriangleToIntersectionPoints.find(pTriangle)==pairTriangleToIntersectionPoints.end()) {
     (listVertexIntersectionPerPairEdges[objectp]).end()) {  
1036              //ffout(0) << "rajout de pTriangle\n";              //ffout(0) << "rajout de pTriangle\n";
1037              (listVertexIntersectionPerPairEdges[objectp])[pTriangle] = Pinter;              pairTriangleToIntersectionPoints[pTriangle] = Pinter;
1038          }          }
1039      }      }
1040      else if(Pinter.size()!=0 and numobject<objectNumber) {      else if(Pinter.size()!=0 and numobject<objectNumber) {
1041          (listVertexIntersectionPerPairEdges[objectp])[pTriangle] = Pinter;          //ffout(0)<<"rajout2 \n";
1042            pairTriangleToIntersectionPoints[pTriangle] = Pinter;
1043      }      }
1044    
1045      if((listVertexIntersectionPerPairEdges[objectp]).find(pTriangle)!=      if(pairTriangleToIntersectionPoints.find(pTriangle)!= pairTriangleToIntersectionPoints.end() and
        (listVertexIntersectionPerPairEdges[objectp]).end() and  
1046         numobject<objectNumber) {         numobject<objectNumber) {
1047          TinyVector<3,size_t>& cutTr2=          TinyVector<3,size_t >& cutTr2= pairTriangleToIntersectionPoints[pTriangle].cutTriangle2;
         (listVertexIntersectionPerPairEdges[objectp])[pTriangle].cutTriangle2;  
1048    
1049          for(size_t i=0 ; i<Pinter.size() ; ++i) {          for(size_t i=0 ; i<Pinter.size() ; ++i) {
1050              real_t epsilon=1e-10;              real_t epsilon=1e-10;
1051              const TinyVector<2,TinyVector<2,Edge::Pair> >& eo0=((*currentTriangle))[0];              const Triangle& To=(*(*currentTriangle));
1052              const TinyVector<2,TinyVector<2,Edge::Pair> >& eo1=((*currentTriangle))[1];              const Vertex& Vo0 = To(0);
1053              const TinyVector<2,TinyVector<2,Edge::Pair> >& eo2=((*currentTriangle))[2];              const Vertex& Vo1 = To(1);
1054              const Vertex& Vo0 = *__determinateVertex(objectNumber,eo0);              const Vertex& Vo2 = To(2);
1055              const Vertex& Vo1 = *__determinateVertex(objectNumber,eo1);              const Vertex& P=*(pairTriangleToIntersectionPoints[pTriangle].listofPoints[i]);
1056              const Vertex& Vo2 = *__determinateVertex(objectNumber,eo2);              
             std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> >  
                 num=((listVertexIntersectionPerPairEdges[objectp])[pTriangle].listofPoints[i]);  
             Vertex P;  
             if(listVerticesIntersectionPerEdges.find(num)!=listVerticesIntersectionPerEdges.end()) {  
                 P=(*listVerticesIntersectionPerEdges[num]);  
             }  
   
1057              real_t Vo01=Norm(Vo1.Vector3()-Vo0.Vector3());              real_t Vo01=Norm(Vo1.Vector3()-Vo0.Vector3());
1058              real_t Vo12=Norm(Vo2.Vector3()-Vo1.Vector3());              real_t Vo12=Norm(Vo2.Vector3()-Vo1.Vector3());
1059              real_t Vo02=Norm(Vo2.Vector3()-Vo0.Vector3());              real_t Vo02=Norm(Vo2.Vector3()-Vo0.Vector3());
# Line 1392  __createPointIntersection(const size_t& Line 1062  __createPointIntersection(const size_t&
1062              real_t norm2=Norm(P.Vector3()-Vo2.Vector3());              real_t norm2=Norm(P.Vector3()-Vo2.Vector3());
1063              //for filter points              //for filter points
1064              real_t eps=1e-10;              real_t eps=1e-10;
1065              if(norm0>eps and norm1>eps and std::abs((norm0+norm1)-Vo01)<epsilon) {              if(/*norm0>eps and norm1>eps and*/ std::abs((norm0+norm1)-Vo01)<epsilon) {
1066                  cutTr2[0]=i+1;                  cutTr2[0]=i+1;
1067              }              }
1068              if(norm0>eps and norm2>eps and std::abs((norm0+norm2)-Vo02)<epsilon) {              if(/*norm0>eps and norm2>eps and */std::abs((norm0+norm2)-Vo02)<epsilon) {
1069                  cutTr2[1]=i+1;                  cutTr2[1]=i+1;
1070              }              }
1071              if(norm2>eps and norm1>eps and std::abs((norm2+norm1)-Vo12)<epsilon) {              if(/*norm2>eps and norm1>eps and */std::abs((norm2+norm1)-Vo12)<epsilon) {
1072                  cutTr2[2]=i+1;                  cutTr2[2]=i+1;
1073              }              }
1074          }          }
# Line 1406  __createPointIntersection(const size_t& Line 1076  __createPointIntersection(const size_t&
1076  }        }      
1077    
1078    
1079  std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > SurfaceMeshGenerator::Internals::  const Vertex& SurfaceMeshGenerator::Internals::
1080  __determinateNumber(const int& p,  __determinateNumber(const int& p,
1081      const TinyVector<3,std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > >& pointsint) {      const TinyVector<3,TinyVector<2,const Vertex*> >& pointsint) {
1082    
1083      if(p==0) {      if(p==0) {
1084          return pointsint[0];          return *pointsint[0][0];
1085      }      }
1086      else if(p==1) {      else if(p==1) {
1087          return pointsint[1];          return *pointsint[1][0];
1088      }      }
1089      else {      else {
1090          return pointsint[2];          return *pointsint[2][0];
1091      }      }
1092  }  }
1093    
1094  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1095  __addTriangle(const size_t& numobject,  __addTriangle(const Vertex& V0,
1096                const TinyVector<2,TinyVector<2,Edge::Pair> >& numA,                const Vertex& V1,
1097                const TinyVector<2,TinyVector<2,Edge::Pair> >& numB,                const Vertex& V2,
1098                const TinyVector<2,TinyVector<2,Edge::Pair> >& numC,                const Cell* cell) {
1099                const Cell* cell,      Triangle T(V0,V1,V2);
1100                size_t& kill) {      T.setMother(cell);
1101        triangleListIntersectionNew.push_back(T);
     TinyVector<3,TinyVector<2,TinyVector<2,Edge::Pair> > > trianglebis;  
     trianglebis[0]=numA;  
     trianglebis[1]=numB;  
     trianglebis[2]=numC;  
     (*listTriangleintersection[numobject]).push_back(trianglebis);  
     localListTriangleintersection.push_back(trianglebis);  
     (*triangleListesIntersectionPerHexa[cell][numobject]).push_front(trianglebis);  
     kill=1;  
1102  }  }
1103    
1104  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1105  __determinateCase(const int& ncase,  __determinateCase(const int& ncase,
1106                    const TriangleCut & K,                    const TriangleCut & K,
1107                    Vertex & A,Vertex & B,Vertex & C,                    TinyVector<3,int>& tt,
1108                    TinyVector<2,TinyVector<2,Edge::Pair> >& e0,                    TinyVector<3,int>& pp) {
                   TinyVector<2,TinyVector<2,Edge::Pair> >& e1,  
                   TinyVector<2,TinyVector<2,Edge::Pair> >& e2,  
                   TinyVector<2,int>& tt,  
                   TinyVector<2,int>& pp) {  
   
     const TinyVector<3,int>& tt2=K.edgecut;  
     const size_t& numobject=K.numobject;  
     Vertex V0,V1,V2;  
     if(K.listpairedges[0][0][0]==K.listpairedges[0][1][0])  
         V0=(*(*verticesListesIntersection[numobject])[K.listpairedges[0][0][0]]);  
     else {  
         std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > p0;  
         p0.first=K.listpairedges[0];  
         TinyVector<2,size_t> num;  
         num[0]=K.numobject;  
         num[1]=K.otherobject;  
         p0.second=num;  
         SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator  
         ip0=listVerticesIntersectionPerEdges.find(p0);  
         if(ip0!=listVerticesIntersectionPerEdges.end())  
             V0 =*(*ip0).second;  
     }  
     if(K.listpairedges[1][0][0]==K.listpairedges[1][1][0])  
         V1=(*(*verticesListesIntersection[numobject])[K.listpairedges[1][0][0]]);  
     else {  
         std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > p0;  
         p0.first=K.listpairedges[1];  
         TinyVector<2,size_t> num;  
         num[0]=0;  
         num[1]=0;  
         p0.second=num;  
         SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator  
         ip1=listVerticesIntersectionPerEdges.find(p0);  
         if(ip1!=listVerticesIntersectionPerEdges.end())  
             V1 =*(*ip1).second;  
     }  
       
     if(K.listpairedges[2][0][0]==K.listpairedges[2][1][0])  
         V2=(*(*verticesListesIntersection[numobject])[K.listpairedges[2][0][0]]);  
     else {  
         std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > p0;  
         p0.first=K.listpairedges[2];  
         TinyVector<2,size_t> num;  
         num[0]=0;  
         num[1]=0;  
         p0.second=num;  
         SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator  
         ip2=listVerticesIntersectionPerEdges.find(p0);  
         if(ip2!=listVerticesIntersectionPerEdges.end())  
             V2 =*(*ip2).second;  
     }  
       
1109    
1110        const TinyVector<3,TinyVector<2,int> >& tt2=K.edgecut;
1111      if(ncase==0) {      if(ncase==0) {
1112          tt[0]=tt2[0];          tt[0]=tt2[0][0];
1113          tt[1]=tt2[1];          tt[1]=tt2[1][0];
1114            tt[2]=tt2[2][0];
1115          pp[0]=0;          pp[0]=0;
1116          pp[1]=1;          pp[1]=1;
1117          A=V0;          pp[2]=2;
1118          e0=K.listpairedges[0];          /*if(tt[0]==0 and tt[1]==1) {
         if(tt[0]==0 and tt[1]==1) {  
1119              B=V2;              B=V2;
             e1=K.listpairedges[2];  
1120          }          }
1121          if(tt[0]==1 and tt[1]==0) {          if(tt[0]==1 and tt[1]==0) {
1122              C=V1;              C=V1;
1123              e2=K.listpairedges[1];          }*/
         }  
1124      }      }
1125      else if(ncase==1) {      else if(ncase==1) {
1126          tt[0]=tt2[0];          tt[0]=tt2[0][0];
1127          tt[1]=tt2[2];          tt[1]=tt2[2][0];
1128            tt[2]=tt2[1][0];
1129          pp[0]=0;          pp[0]=0;
1130          pp[1]=2;          pp[1]=2;
1131          A=V1;          pp[2]=1;
1132          e0=K.listpairedges[1];          /*if(tt[0]==0 and tt[1]==1) {
         if(tt[0]==0 and tt[1]==1) {  
1133              B=V2;              B=V2;
             e1=K.listpairedges[2];  
1134          }          }
1135          if(tt[0]==1 and tt[1]==0) {          if(tt[0]==1 and tt[1]==0) {
1136              C=V0;              C=V0;
1137              e2=K.listpairedges[0];          }*/
         }  
1138      }      }
1139      else if(ncase==2) {      else if(ncase==2) {
1140          tt[0]=tt2[1];          tt[0]=tt2[1][0];
1141          tt[1]=tt2[2];          tt[1]=tt2[2][0];
1142            tt[2]=tt2[0][0];
1143          pp[0]=1;          pp[0]=1;
1144          pp[1]=2;          pp[1]=2;
1145          A=V2;          pp[2]=0;
1146          e0=K.listpairedges[2];          /*if(tt[0]==0 and tt[1]==1) {
         if(tt[0]==0 and tt[1]==1) {  
1147              B=V0;              B=V0;
             e1=K.listpairedges[0];  
1148          }          }
1149          if(tt[0]==1 and tt[1]==0) {          if(tt[0]==1 and tt[1]==0) {
1150              C=V1;              C=V1;
1151              e2=K.listpairedges[1];          }*/
         }  
1152      }      }
1153    
1154    
# Line 1551  __determinateCase(const int& ncase, Line 1156  __determinateCase(const int& ncase,
1156  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1157  __createCase1(int ncase,  __createCase1(int ncase,
1158                const TriangleCut& K,                const TriangleCut& K,
1159                const Cell* cell,                const Cell* cell) {
1160                size_t & kill) {      const Vertex& V0=K(0);
1161      const size_t& numobject=K.numobject;      const Vertex& V1=K(1);
1162      const TinyVector<2,TinyVector<2,Edge::Pair> >& e0=K.listpairedges[0];      const Vertex& V2=K(2);
1163      const TinyVector<2,TinyVector<2,Edge::Pair> >& e1=K.listpairedges[1];      
1164      const TinyVector<2,TinyVector<2,Edge::Pair> >& e2=K.listpairedges[2];      const TinyVector<3,TinyVector<2,const Vertex*> >& pointsint=K.pointsIntersection;
1165      Vertex V0,V1,V2;      const std::map< const Vertex* ,const Vertex* >& pointsin=K.pointsin;
1166      std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > p0,p1,p2;      const size_t& numberin=pointsin.size();
1167        const Vertex& A=K(ncase);
1168    
1169        TinyVector<3,int> tt;
1170        TinyVector<3,int> pp;
1171    
1172        __determinateCase(ncase,K,tt,pp);
1173    
     if(e0[0][0]==e0[1][0])  
         V0 = *((*verticesListesIntersection[numobject])[e0[0][0]]);  
     else {  
         p0.first=e0;  
         TinyVector<2,size_t> num;  
         num[0]=K.numobject;  
         num[1]=K.otherobject;  
         p0.second=num;  
         SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator  
         ip0=listVerticesIntersectionPerEdges.find(p0);  
         if(ip0!=listVerticesIntersectionPerEdges.end())  
             V0 =*(*ip0).second;  
     }  
     if(e1[0][0]==e1[1][0])  
         V1 = *((*verticesListesIntersection[numobject])[e1[0][0]]);  
     else {  
         p1.first=e1;  
         TinyVector<2,size_t> num;  
         num[0]=K.numobject;  
         num[1]=K.otherobject;  
         p0.second=num;  
         SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator  
         ip1=listVerticesIntersectionPerEdges.find(p1);  
         if(ip1!=listVerticesIntersectionPerEdges.end())  
             V1 =*(*ip1).second;  
     }  
     if(e2[0][0]==e2[1][0])  
         V2 = *((*verticesListesIntersection[numobject])[e2[0][0]]);  
     else {  
         p2.first=e2;  
         TinyVector<2,size_t> num;  
         num[0]=K.numobject;  
         num[1]=K.otherobject;  
         p0.second=num;  
         SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator  
         ip2=listVerticesIntersectionPerEdges.find(p2);  
         if(ip2!=listVerticesIntersectionPerEdges.end())  
             V2 =*(*ip2).second;  
     }  
     const TinyVector<3,std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > >&  
         pointsint=K.pointsIntersection;  
     const std::vector< std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> >  >&  
         pointsin=K.pointsin;  
     const int& numberin=pointsin.size();  
     Vertex A=V0,B=V1,C=V2;  
     TinyVector<2,TinyVector<2,Edge::Pair> > eA=e0,eB=e1,eC=e2;  
     std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > numA,numB,numC,numD;  
   
     TinyVector<2,int> tt;  
     TinyVector<2,int> pp;  
   
     __determinateCase(ncase,K,A,B,C,eA,eB,eC,tt,pp);  
   
     std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > pA,pB;  
     TinyVector<2,size_t> nuum;  
     nuum[0]=0;  
     nuum[1]=0;  
     pA.second=nuum;  
     pB.second=nuum;  
     pA.first=eA;  
     pB.first=eB;  
     numA=pA;  
1174      if(numberin==0) {      if(numberin==0) {
1175          if(tt[0]==1) {          if(tt[0]==1) {
1176              numB=__determinateNumber(pp[0],pointsint);              const Vertex& B=__determinateNumber(pp[0],pointsint);
   
             if(listVerticesIntersectionPerEdges.find(numB)!=listVerticesIntersectionPerEdges.end()) {  
                 B=(*listVerticesIntersectionPerEdges[numB]);  
             }  
             else {  
                 ffout(0)<<"1 "<<pp[0]<<" pbs point inexistant\n";  
                 std::exit(1);  
             }  
1177              if(tt[1]==1) {              if(tt[1]==1) {
1178                    const Vertex& C=__determinateNumber(pp[1],pointsint);
1179                  numC=__determinateNumber(pp[1],pointsint);                  __addTriangle(A,B,C,cell);
                 if(listVerticesIntersectionPerEdges.find(numC)!=listVerticesIntersectionPerEdges.end()) {  
                     C=(*listVerticesIntersectionPerEdges[numC]);  
                 }  
                 else {  
                     ffout(0)<<"2 "<<pp[1]<<"  pbs point inexistant\n";  
                     std::exit(1);  
                 }  
                 __addTriangle(numobject,numA.first,numB.first,numC.first,cell,kill);  
   
1180              }              }
1181          else {          else {
1182              ffout(0) <<"pbs 10\n";              ffout(0) <<"pbs 10\n";
1183              printpoint(numobject,K.listpairedges[0],K);              ffout(0)<<V0.Vector3()<<"\n";
1184              printpoint(numobject,K.listpairedges[1],K);              ffout(0)<<V1.Vector3()<<"\n";
1185              printpoint(numobject,K.listpairedges[2],K);              ffout(0)<<V2.Vector3()<<"\n";
1186              /*listpoints.push_back(new Vertex(C.Vector3()));              /*listpoints.push_back(new Vertex(C.Vector3()));
1187              numC=(listpoints.size()-1);              numC=(listpoints.size()-1);
1188              __addTriangle(numobject,numA,numB,numC,kill);              __addTriangle(numobject,numA,numB,numC);*/
1189              */              
1190          }          }
1191      }      }
1192      else  {      else  {
1193          if(tt[1]==1) {          if(tt[1]==1) {
1194              ffout(0)<<"pbs 11\n";              ffout(0)<<"pbs 11\n";
             printpoint(numobject,K.listpairedges[0],K);  
             printpoint(numobject,K.listpairedges[1],K);  
             printpoint(numobject,K.listpairedges[2],K);  
1195              /* numC=__determinateNumber(pp[1],pointsint);              /* numC=__determinateNumber(pp[1],pointsint);
1196            C=(*listpoints[numC]);            C=(*listpoints[numC]);
1197            listpoints.push_back(new Vertex(A.Vector3()));            listpoints.push_back(new Vertex(A.Vector3()));
1198            numA=(listpoints.size()-1);            numA=(listpoints.size()-1);
1199            listpoints.push_back(new Vertex(B.Vector3()));            listpoints.push_back(new Vertex(B.Vector3()));
1200            numC=(listpoints.size()-1);            numC=(listpoints.size()-1);
1201            __addTriangle(numobject,numA,numB,numC,kill);*/            __addTriangle(numobject,numA,numB,numC);*/
1202          }          }
1203          else {          else {
1204              ffout(0) <<"pbs 111\n";              ffout(0) <<"pbs 111\n";
1205              printpoint(numobject,K.listpairedges[0],K);              //__addTriangle(numobject,K.otherobject,e0,e1,e2,cell);
             printpoint(numobject,K.listpairedges[1],K);  
             printpoint(numobject,K.listpairedges[2],K);  
               
1206          }          }
1207      }      }
1208    }    }
1209    else if(numberin==1 or numberin==2) {    else if(numberin==1 /*or numberin==2 or numberin==3*/) {
1210                std::map<const Vertex*,const Vertex*>::const_iterator it=pointsin.begin();
1211        numD=(pointsin[0]);        const Vertex&  D=*((*it).second);
1212        if(listVerticesIntersectionPerEdges.find(numD)!=listVerticesIntersectionPerEdges.end()) {            
       Vertex D=(*listVerticesIntersectionPerEdges[numD]);  
       }  
       else {  
           ffout(0)<<"3 "<<pp[0]<<"  pbs point inexistant\n";  
           std::exit(1);  
       }  
         
1213        if(tt[0]==1) {        if(tt[0]==1) {
1214            numB=__determinateNumber(pp[0],pointsint);            const Vertex& B=__determinateNumber(pp[0],pointsint);
           if(listVerticesIntersectionPerEdges.find(numB)!=listVerticesIntersectionPerEdges.end()) {  
               B=(*listVerticesIntersectionPerEdges[numB]);  
           }  
           else {  
               ffout(0)<<"4 "<<pp[0]<<"  pbs point inexistant\n";  
               std::exit(1);  
           }  
   
1215            if(tt[1]==1) {            if(tt[1]==1) {
1216                numC=__determinateNumber(pp[1],pointsint);                const Vertex& C=__determinateNumber(pp[1],pointsint);
1217                if(listVerticesIntersectionPerEdges.find(numC)!=listVerticesIntersectionPerEdges.end()) {                
1218                C=(*listVerticesIntersectionPerEdges[numC]);                __addTriangle(A,B,D,cell);
1219                }                __addTriangle(A,D,C,cell);
               else {  
                   ffout(0)<<"5 "<<pp[1]<<"  pbs point inexistant\n";  
                   std::exit(1);  
               }  
               __addTriangle(numobject,numA.first,numB.first,numD.first,cell,kill);  
               __addTriangle(numobject,numA.first,numD.first,numC.first,cell,kill);  
1220            }            }
1221            else {            else {
1222                ffout(0)<<"-------- warning trou in 0\n";                ffout(0)<<"-------- warning trou in 0\n";
1223                printpoint(numobject,K.listpairedges[0],K);                ffout(0)<<V0.Vector3()<<"\n";
1224                printpoint(numobject,K.listpairedges[1],K);                ffout(0)<<V1.Vector3()<<"\n";
1225                printpoint(numobject,K.listpairedges[2],K);                      ffout(0)<<V2.Vector3()<<"\n";
1226                __addTriangle(numobject,numA.first,numB.first,numD.first,cell,kill);                __addTriangle(A,B,D,cell);
1227            }            }
1228        }        }
1229        else if(tt[1]==1) {        else if(tt[1]==1) {
1230            ffout(0)<<"-------- warning trou in\n";            ffout(0)<<"-------- warning trou in\n";
1231            printpoint(numobject,K.listpairedges[0],K);            const Vertex& B=__determinateNumber(pp[1],pointsint);
1232            printpoint(numobject,K.listpairedges[1],K);            __addTriangle(A,B,D,cell);
           printpoint(numobject,K.listpairedges[2],K);    
           numB=__determinateNumber(pp[1],pointsint);  
           if(listVerticesIntersectionPerEdges.find(numB)!=listVerticesIntersectionPerEdges.end()) {  
               B=(*listVerticesIntersectionPerEdges[numB]);  
           }  
           else {  
               ffout(0)<<"6 "<<pp[0]<<"  pbs point inexistant\n";  
               std::exit(1);  
           }  
           __addTriangle(numobject,numA.first,numB.first,numD.first,cell,kill);  
1233    
1234        }        }
1235        else {        else {
1236            ffout(0) <<"pbs 1\n";            ffout(0) <<"pbs 1\n";
1237            printpoint(numobject,K.listpairedges[0],K);           //__addTriangle(numobject,K.otherobject,e0,e1,e2,cell);
           printpoint(numobject,K.listpairedges[1],K);  
           printpoint(numobject,K.listpairedges[2],K);      
1238        }        }
1239    }    }
1240    else {    else {
1241    #warning is it true?
1242      //a faire si plus d' 1 point d'intersection pas dans le triangle      //a faire si plus d' 1 point d'intersection pas dans le triangle
1243      ffout(0) << numberin<<" be careful 1\n";      ffout(0) << numberin<<" be careful 1\n";
1244        //  ffout(0)<<V0.Vector3()<<"\n";
1245        //  ffout(0)<<V1.Vector3()<<"\n";
1246        //  ffout(0)<<V2.Vector3()<<"\n";
1247          TinyVector<10,const Vertex*> P,Ptemp;
1248          size_t i=0;
1249          for(std::map<const Vertex*,const Vertex*>::const_iterator it=pointsin.begin();
1250              it!=pointsin.end() ; ++it) {
1251              P[i+1]=(*it).second;
1252              ++i;
1253          }
1254    
1255          if(tt[0]==1) {
1256              const Vertex& B=__determinateNumber(pp[0],pointsint);
1257              P[0]=&B;
1258              if(tt[1]==1) {
1259                  const Vertex& C=__determinateNumber(pp[1],pointsint);
1260                  size_t size=numberin+1;
1261                  for(size_t boucle=0 ; boucle<numberin-1 ; ++boucle) {
1262                      size_t num=1;
1263                      real_t norm=Norm((*P[0]).Vector3()-(*P[1]).Vector3());
1264                      for(size_t i=2 ; i<size ; ++i) {
1265                          real_t normtemp=Norm((*P[0]).Vector3()-(*P[i]).Vector3());
1266                          if(normtemp<=norm) {
1267                              norm=normtemp;
1268                              Ptemp[i-1]=P[num];
1269                              num=i;
1270                          }
1271                          else {
1272                              Ptemp[i-1]=P[i];
1273                          }
1274                      }
1275                      __addTriangle(A,*P[0],*P[num],cell);
1276                      Ptemp[0]=P[num];
1277                      P=Ptemp;
1278                      --size;
1279                     // ffout(0)<<"size = "<<size<<"\n";
1280                  }
1281                  __addTriangle(A,*P[0],*P[1],cell);
1282                  __addTriangle(A,*P[1],C,cell);
1283                /*  if(tt[2]==1) {
1284                      const Vertex& D=__determinateNumber(pp[2],pointsint);
1285                      for(size_t i=0 ; i<numberin ; ++i) {
1286                          P[i+1]=(pointsin[i]);
1287                      }
1288                      P[0]=&B;
1289                      size=numberin+1;
1290                      for(size_t boucle=0 ; boucle<numberin-1 ; ++boucle) {
1291                          size_t num=1;
1292                          real_t norm=Norm((*P[0]).Vector3()-(*P[1]).Vector3());
1293                          for(size_t i=2 ; i<size ; ++i) {
1294                              real_t normtemp=Norm((*P[0]).Vector3()-(*P[i]).Vector3());
1295                              if(normtemp<=norm) {
1296                                  norm=normtemp;
1297                                  Ptemp[i-1]=P[num];
1298                                  num=i;
1299                              }
1300                              else {
1301                                  Ptemp[i-1]=P[i];
1302                              }
1303                          }
1304                          __addTriangle(D,*P[0],*P[num],cell);
1305                          Ptemp[0]=P[num];
1306                          P=Ptemp;
1307                          --size;
1308                          // ffout(0)<<"size = "<<size<<"\n";
1309                      }
1310                      __addTriangle(D,*P[0],*P[1],cell);
1311                    //  __addTriangle(D,*P[1],C,cell);
1312                  }*/
1313              }
1314              else {
1315                  ffout(0)<<"une seule arete cut pbs\n";
1316              }
1317          }
1318          else {
1319              ffout(0)<<"une seule arete cut pbs2\n";
1320          }
1321    
1322    }    }
1323  }  }
1324    
1325  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1326  __createCase2(int ncase,  __createCase2(int ncase,
1327                const TriangleCut& K,                const TriangleCut& K,
1328                const Cell* cell,                const Cell* cell) {
1329                size_t & kill) {    
   
1330    int ttn0,ttn1;    int ttn0,ttn1;
1331    int p0,p1;    int p0,p1;
1332    const size_t& numobject=K.numobject;    const Vertex& V0=K(0);
1333    const TinyVector<2,TinyVector<2,Edge::Pair> >& e0=K.listpairedges[0];    const Vertex& V1=K(1);
1334    const TinyVector<2,TinyVector<2,Edge::Pair> >& e1=K.listpairedges[1];    const Vertex& V2=K(2);
1335    const TinyVector<2,TinyVector<2,Edge::Pair> >& e2=K.listpairedges[2];  
1336    std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > pp0,pp1,pp2;    const TinyVector<3,TinyVector<2,const Vertex*> >& pointsint=K.pointsIntersection;
1337      const std::map<const Vertex* ,const Vertex*>& pointsin=K.pointsin;
1338    Vertex V0,V1,V2;    const TinyVector<3,TinyVector<2,int> >&  tt2=K.edgecut;
   if(e0[0][0]==e0[1][0]) {  
       V0 = *((*verticesListesIntersection[numobject])[e0[0][0]]);  
   }  
   else {  
       pp0.first=e0;  
       TinyVector<2,size_t> num;  
       num[0]=K.numobject;  
       num[1]=K.otherobject;  
       pp0.second=num;  
         SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator  
         ip0=listVerticesIntersectionPerEdges.find(pp0);  
         if(ip0!=listVerticesIntersectionPerEdges.end())  
             V0 =*(*ip0).second;  
     }  
   if(e1[0][0]==e1[1][0]) {  
       V1 = *((*verticesListesIntersection[numobject])[e1[0][0]]);  
   }  
   else {  
       pp1.first=e1;  
       TinyVector<2,size_t> num;  
       num[0]=K.numobject;  
       num[1]=K.otherobject;  
       pp1.second=num;  
       SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator  
       ip1=listVerticesIntersectionPerEdges.find(pp1);  
       if(ip1!=listVerticesIntersectionPerEdges.end())  
           V1 =*(*ip1).second;  
   }  
   if(e2[0][0]==e2[1][0]) {  
       V2 = *((*verticesListesIntersection[numobject])[e2[0][0]]);  
   }  
   else {  
       pp2.first=e2;  
       TinyVector<2,size_t> num;  
       num[0]=K.numobject;  
       num[1]=K.otherobject;  
       pp2.second=num;  
       SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator  
       ip2=listVerticesIntersectionPerEdges.find(pp2);  
       if(ip2!=listVerticesIntersectionPerEdges.end())  
           V2 =*(*ip2).second;  
   }  
   const TinyVector<3,std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > >&  
       pointsint=K.pointsIntersection;  
   const std::vector< std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> >  >&  
       pointsin=K.pointsin;  
   const TinyVector<3,int>&  tt2=K.edgecut;  
1339    const int& numberin=pointsin.size();    const int& numberin=pointsin.size();
1340    
1341    Vertex A=V0,B=V1;    size_t ncase0=0,ncase1=1;
1342    TinyVector<2,TinyVector<2,Edge::Pair> > eA=e0,eB=e1;    
   std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > pA=pp0,pB=pp1;  
   std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > numA,numB,numC,numD,numE;  
   
1343    if(ncase==0) {    if(ncase==0) {
1344        ttn0=tt2[1];        ttn0=tt2[1][0];
1345        ttn1=tt2[2];        ttn1=tt2[2][0];
1346        p0=1;        p0=1;
1347        p1=2;        p1=2;
1348    }    }
1349    else if(ncase==1) {    else if(ncase==1) {
1350      ttn0=tt2[0];      ttn0=tt2[0][0];
1351      ttn1=tt2[2];      ttn1=tt2[2][0];
1352      p0=0;      p0=0;
1353      p1=2;      p1=2;
1354      A=V0;      ncase0=0;
1355      eA=K.listpairedges[0];      ncase1=2;
     B=V2;  
     eB=K.listpairedges[2];  
1356    }    }
1357    else if(ncase==2) {    else if(ncase==2) {
1358      ttn0=tt2[0];      ttn0=tt2[0][0];
1359      ttn1=tt2[1];      ttn1=tt2[1][0];
1360      p0=0;      p0=0;
1361      p1=1;      p1=1;
1362      A=V1;      ncase0=1;
1363      eA=K.listpairedges[1];      ncase1=2;
     B=V2;  
     eB=K.listpairedges[2];  
1364    }    }
1365    
1366    pA.first=eA;    const Vertex& A=K(ncase0);
1367    pB.first=eB;    const Vertex& B=K(ncase1);
   numA=pA;  
   numB=pB;  
     
1368    if(numberin==0) {    if(numberin==0) {
1369        if(ttn0==1) {        if(ttn0==1) {
1370            Vertex C;            const Vertex& C=__determinateNumber(p0,pointsint);
           numC=__determinateNumber(p0,pointsint);  
           if(listVerticesIntersectionPerEdges.find(numC)!=listVerticesIntersectionPerEdges.end()) {  
               C=(*listVerticesIntersectionPerEdges[numC]);  
           }  
1371            if(ttn1==1) {            if(ttn1==1) {
1372                Vertex D;                const Vertex& D=__determinateNumber(p1,pointsint);
1373                numD=__determinateNumber(p1,pointsint);                __addTriangle(A,B,C,cell);
1374                if(listVerticesIntersectionPerEdges.find(numD)!=listVerticesIntersectionPerEdges.end()) {                __addTriangle(C,B,D,cell);
                   D=(*listVerticesIntersectionPerEdges[numD]);  
               }  
               __addTriangle(numobject,numA.first,numB.first,numC.first,cell,kill);  
               __addTriangle(numobject,numC.first,numB.first,numD.first,cell,kill);  
1375            }            }
1376            else {            else {
1377               // ffout(0) <<"pbs 20\n";               // ffout(0) <<"pbs 20\n";
1378                __addTriangle(numobject,numA.first,numB.first,numC.first,cell,kill);                __addTriangle(A,B,C,cell);
1379            }            }
1380        }        }
1381        else {        else {
1382            if(ttn1==1) {            if(ttn1==1) {
1383               // ffout(0) <<"pbs 210\n";               // ffout(0) <<"pbs 210\n";
1384                Vertex C;                const Vertex& C=__determinateNumber(p1,pointsint);
1385                numC=__determinateNumber(p1,pointsint);                __addTriangle(A,B,C,cell);
               if(listVerticesIntersectionPerEdges.find(numC)!=listVerticesIntersectionPerEdges.end()) {  
                   C=(*listVerticesIntersectionPerEdges[numC]);  
               }  
               __addTriangle(numobject,numA.first,numB.first,numC.first,cell,kill);  
1386    
1387            }            }
1388            else {            else {
1389                ffout(0) <<"pbs 21\n";                ffout(0) <<"pbs 21\n";
1390                printpoint(numobject,K.listpairedges[0],K);                ffout(0)<<V0.Vector3()<<"\n";
1391                printpoint(numobject,K.listpairedges[1],K);                ffout(0)<<V1.Vector3()<<"\n";
1392                printpoint(numobject,K.listpairedges[2],K);                ffout(0)<<V2.Vector3()<<"\n";
1393                               // __addTriangle(V0,V1,V2,cell);
1394            }            }
1395        }        }
1396    }    }
1397    else if(numberin==1 or numberin==2) {    else if(numberin==1 or numberin==2) {
1398        numE=(pointsin[0]);        std::map<const Vertex*,const Vertex*>::const_iterator it=pointsin.begin();
1399          const Vertex& E=*((*it).second);
1400    
1401        if(ttn0==1) {        if(ttn0==1) {
1402            Vertex C;            const Vertex& C=__determinateNumber(p0,pointsint);
           numC=__determinateNumber(p0,pointsint);  
           if(listVerticesIntersectionPerEdges.find(numC)!=listVerticesIntersectionPerEdges.end()) {  
               C=(*listVerticesIntersectionPerEdges[numC]);  
           }  
1403            if(ttn1==1) {            if(ttn1==1) {
1404                Vertex D;                const Vertex& D=__determinateNumber(p1,pointsint);
1405                numD=__determinateNumber(p1,pointsint);                
1406                if(listVerticesIntersectionPerEdges.find(numD)!=listVerticesIntersectionPerEdges.end()) {                __addTriangle(A,C,E,cell);
1407                    D=(*listVerticesIntersectionPerEdges[numD]);                __addTriangle(A,E,B,cell);
1408                }                __addTriangle(B,E,D,cell);
               __addTriangle(numobject,numA.first,numC.first,numE.first,cell,kill);  
               __addTriangle(numobject,numA.first,numE.first,numB.first,cell,kill);  
               __addTriangle(numobject,numB.first,numE.first,numD.first,cell,kill);  
1409            }            }
1410            else {            else {
1411                __addTriangle(numobject,numA.first,numB.first,numE.first,cell,kill);                __addTriangle(A,B,E,cell);
1412                __addTriangle(numobject,numC.first,numE.first,numB.first,cell,kill);                __addTriangle(C,E,B,cell);
1413            }            }
1414        }        }
1415        else {        else {
1416            if(ttn1==1) {            if(ttn1==1) {
1417                Vertex D;                const Vertex& D=__determinateNumber(p1,pointsint);
1418                numD=__determinateNumber(p1,pointsint);                __addTriangle(D,E,B,cell);
1419                if(listVerticesIntersectionPerEdges.find(numD)!=listVerticesIntersectionPerEdges.end()) {                __addTriangle(A,B,E,cell);
                   D=(*listVerticesIntersectionPerEdges[numD]);  
               }  
               __addTriangle(numobject,numD.first,numE.first,numB.first,cell,kill);  
               __addTriangle(numobject,numA.first,numB.first,numE.first,cell,kill);  
1420            }            }
1421            else {            else {
1422                __addTriangle(numobject,numA.first,numB.first,numE.first,cell,kill);                __addTriangle(A,B,E,cell);
1423            }            }
1424        }        }
1425    }    }
1426    else if(numberin>2) {    else if(numberin>2) {
1427        ffout(0) << numberin <<"be careful 2\n";        ffout(0) << numberin <<"be careful 2\n";
1428    }        TinyVector<10,const Vertex*> P,Ptemp;
1429          size_t i=0;
1430          for(std::map<const Vertex*,const Vertex*>::const_iterator it=pointsin.begin();
1431              it!=pointsin.end() ; ++it) {
1432              P[i+1]=(*it).second;
1433              ++i;
1434          }
1435          if(ttn0==1) {
1436              const Vertex& C=__determinateNumber(p0,pointsint);
1437              __addTriangle(A,B,C,cell);
1438              P[0]=&C;
1439              if(ttn1==1) {
1440                  const Vertex& D=__determinateNumber(p1,pointsint);
1441                  size_t size=numberin+1;
1442                  for(size_t boucle=0 ; boucle<numberin-1 ; ++boucle) {
1443                      size_t num=1;
1444                      real_t norm=Norm((*P[0]).Vector3()-(*P[1]).Vector3());
1445                      for(size_t i=2 ; i<size ; ++i) {
1446                          real_t normtemp=Norm((*P[0]).Vector3()-(*P[i]).Vector3());
1447                          if(normtemp<=norm) {
1448                              norm=normtemp;
1449                              Ptemp[i-1]=P[num];
1450                              num=i;
1451                          }
1452                          else {
1453                              Ptemp[i-1]=P[i];
1454                          }
1455                      }
1456                      __addTriangle(B,*P[0],*P[num],cell);
1457                      Ptemp[0]=P[num];
1458                      P=Ptemp;
1459                      --size;
1460                      // ffout(0)<<"size = "<<size<<"\n";
1461                  }
1462                  __addTriangle(B,*P[0],*P[1],cell);
1463                  __addTriangle(B,*P[1],D,cell);
1464                  //__addTriangle(A,B,D,cell);
1465    
1466              }
1467              else {
1468                  ffout(0)<<"une seule arete cut pbs\n";
1469              }
1470          }
1471          else {
1472              ffout(0)<<"une seule arete cut pbs2\n";
1473          }
1474    
1475      }
1476  }  }
1477    
1478  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1479  __createLocalListIntersection(const size_t& objectNumber,  __createLocalListIntersection(const size_t& objectNumber,
1480                                const TinyVector<3,TinyVector<2,TinyVector<2,Edge::Pair> > >& currentTriangle,                                const Triangle*& currentTriangle,
1481                                TriangleCut & K,                                TriangleCut & K,
1482                                const Cell* cell) {                                const Cell* cell) {
1483    
1484      TinyVector<3,int>& tt=K.edgecut;      TinyVector<3,TinyVector<2,int> >& tt=K.edgecut;
1485      TinyVector<3,std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > >&      TinyVector<3,TinyVector<2,const Vertex*> > & pointsIntersection=K.pointsIntersection;
1486          pointsIntersection=K.pointsIntersection;      std::map<const Vertex* ,const Vertex*>& pointsin=K.pointsin;
     std::vector< std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > >&  
         pointsin=K.pointsin;  
1487    
1488      int numberin=pointsin.size();      int numberin=pointsin.size();
1489      const size_t& numobject=K.numobject;      const size_t& numobject=K.numobject;
1490    
1491      ObjectPair objectp(numobject,objectNumber);      ListOfHexaToTriangle::iterator hexait=listOfHexaToTriangle.begin();
1492      if(numobject>objectNumber) {      MapCellTriangle hexaToTriangle1;
1493          objectp.first=objectNumber;      if(numobject!=0)
1494          objectp.second=numobject;          hexaToTriangle1=*(*hexait);
1495        else {
1496            ++hexait;
1497            hexaToTriangle1=*(*hexait);
1498      }      }
1499    
1500      for(SurfaceMeshGenerator::Internals::TriangleVerticesMap::iterator currentTriangle1      MapCellTriangle::iterator it1=hexaToTriangle1.find(cell);
1501          = (*triangleListesIntersectionEdgesTempPerHexa[cell][objectNumber]).begin();      if(it1!=hexaToTriangle1.end()) {
1502          currentTriangle1 != (*triangleListesIntersectionEdgesTempPerHexa[cell][objectNumber]).end();      std::list<Triangle*> listTriangle1=(*it1).second;
1503          ++currentTriangle1) {          for(std::list<Triangle*>::iterator jt0=listTriangle1.begin() ; jt0!=listTriangle1.end() ; ++jt0) {
1504    
1505         const SurfaceMeshGenerator::Internals::TriangleVertices& tv=(*(*currentTriangle1).second);              
1506                        PairTriangleNew pT;
         PairTriangleEdges pT;  
         if(numobject<objectNumber) {  
             pT.first=&(currentTriangle);  
             pT.second=&(tv);  
         }  
         else {  
             pT.first=&(tv);  
             pT.second=&(currentTriangle);  
         }  
         size_t temp=0;  
         if((listVertexIntersectionPerPairEdges[objectp]).find(pT)!=  
     (listVertexIntersectionPerPairEdges[objectp]).end()) {  
             TinyVector<3,size_t> cc;  
1507              if(numobject<objectNumber) {              if(numobject<objectNumber) {
1508                  cc=((listVertexIntersectionPerPairEdges[objectp])[pT]).cutTriangle1;                  pT.first=(currentTriangle);
1509                    pT.second=(*jt0);
1510              }              }
1511              else {              else {
1512                  cc=((listVertexIntersectionPerPairEdges[objectp])[pT]).cutTriangle2;                  pT.second=(currentTriangle);
1513                    pT.first=(*jt0);
1514              }              }
1515              if(cc[0]!=0) {  
1516                  temp+=1;              size_t temp=0;
1517                  tt[0]=1;              
1518                  pointsIntersection[0]=              if(pairTriangleToIntersectionPoints.find(pT)!=pairTriangleToIntersectionPoints.end()) {
1519                      ((listVertexIntersectionPerPairEdges[objectp])[pT].listofPoints[cc[0]-1]);                  TinyVector<3,size_t > cc;
1520              }                  if(numobject<objectNumber) {
1521              if(cc[1]!=0) {                      cc=(pairTriangleToIntersectionPoints[pT]).cutTriangle1;
1522                  temp+=1;                  }
                 tt[1]=1;  
                 pointsIntersection[1]=  
                     ((listVertexIntersectionPerPairEdges[objectp])[pT].listofPoints[cc[1]-1]);  
             }  
             if(cc[2]!=0) {  
                 temp+=1;  
                 tt[2]=1;  
                 pointsIntersection[2]=  
                     ((listVertexIntersectionPerPairEdges[objectp])[pT].listofPoints[cc[2]-1]);  
             }  
             if(temp<(listVertexIntersectionPerPairEdges[objectp])[pT].size()) {  
                 numberin+=1;  
                 size_t r;  
                 size_t sum=cc[0]+cc[1]+cc[2];  
                 //ffout(0)<<(listVertexIntersectionPerPairEdges[objectp])[pT].size()<<" "<<sum<<"\n";  
                 if(sum==0 or sum==2 or sum==5)  
                     r=0;  
                 else if(sum==4 or sum==1)  
                     r=1;  
                 else if(sum==3)  
                     r=2;  
1523                  else {                  else {
1524                      ffout(0)<<"pbs sum = "<<sum<<" size "<<(listVertexIntersectionPerPairEdges[objectp])[pT].size()<<" \n";                      cc=(pairTriangleToIntersectionPoints[pT]).cutTriangle2;
1525                      ffout(0)<<numobject<<" "<<objectNumber<<"\n";                  }
1526                      ffout(0)<<cc[0]<<" "<<cc[1]<<" "<<cc[2]<<"\n";                  if(pairTriangleToIntersectionPoints[pT].size()==2) {
1527                    for(size_t b=0; b<pairTriangleToIntersectionPoints[pT].size() ; ++b) {
1528                        const Vertex*& V=pairTriangleToIntersectionPoints[pT].listofPoints[b];
1529                        if(vertexInTriangles.find(V)==vertexInTriangles.end()) {
1530                            TinyVector<2,const Triangle*> TT;
1531                            TT[0]=*jt0;
1532                            TT[1]=*jt0;
1533                            vertexInTriangles[V]=TT;
1534                        }
1535                        else {
1536                            (*vertexInTriangles.find(V)).second[1]=*jt0;
1537                        }
1538                        if(triangleWithVertex.find(*jt0)==triangleWithVertex.end()) {
1539                            TinyVector<2,const Vertex*> VV;
1540                            VV[0]=V;
1541                            VV[1]=V;
1542                            triangleWithVertex[*jt0]=VV;
1543                        }
1544                        else {
1545                            (*triangleWithVertex.find(*jt0)).second[1]=V;
1546                        }
1547                    }
1548                    }
1549                    if(cc[0]!=0) {
1550                        temp+=1;
1551                        if(tt[0][0]==0) {
1552                            tt[0][0]=1;
1553                            pointsIntersection[0][0]=
1554                            (pairTriangleToIntersectionPoints[pT].listofPoints[cc[0]-1]);
1555                        }
1556                        else {
1557                            tt[0][1]=1;
1558                            pointsIntersection[0][1]=
1559                                (pairTriangleToIntersectionPoints[pT].listofPoints[cc[0]-1]);
1560                        }
1561                        
1562                    }
1563                    if(cc[1]!=0) {
1564                        temp+=1;
1565                        if(tt[1][0]==0) {
1566                            tt[1][0]=1;
1567                            pointsIntersection[1][0]=
1568                                (pairTriangleToIntersectionPoints[pT].listofPoints[cc[1]-1]);
1569                        }
1570                        else{
1571                            tt[1][1]=1;
1572                            pointsIntersection[1][1]=
1573                                (pairTriangleToIntersectionPoints[pT].listofPoints[cc[1]-1]);
1574                        }
1575                        
1576                    }
1577                    if(cc[2]!=0) {
1578                        temp+=1;
1579                        if(tt[2][0]==0) {
1580                            tt[2][0]=1;
1581                            pointsIntersection[2][0]=
1582                                (pairTriangleToIntersectionPoints[pT].listofPoints[cc[2]-1]);
1583                        }
1584                        else{
1585                            tt[2][1]=1;
1586                            pointsIntersection[2][1]=
1587                                (pairTriangleToIntersectionPoints[pT].listofPoints[cc[2]-1]);
1588                        }
1589                        
1590                    }
1591                    if(temp<pairTriangleToIntersectionPoints[pT].size()) {
1592                        numberin+=1;
1593                        size_t r=0;
1594                        size_t sum=cc[0]+cc[1]+cc[2];
1595                        if(sum>3)
1596                            ffout(0)<<"sum "<<sum<<"\n";
1597                        if(sum==0) {
1598                            r=0;
1599                            for(size_t i=1 ; i<pairTriangleToIntersectionPoints[pT].size() ;++i)
1600                                pointsin[(pairTriangleToIntersectionPoints[pT].listofPoints[i])]=
1601                                    (pairTriangleToIntersectionPoints[pT].listofPoints[i]);
1602                        }
1603                        else if(sum==2 or sum==5)
1604                            r=0;
1605                        else if(sum==4 or sum==1)
1606                            r=1;
1607                        else if(sum==3)
1608                            r=2;
1609                        else {
1610                            ffout(0)<<"pbs sum = "<<sum<<" size "<<
1611                            pairTriangleToIntersectionPoints[pT].size()<<" \n";
1612                            ffout(0)<<K(0)<<"\n";
1613                            ffout(0)<<K(1)<<"\n";
1614                            ffout(0)<<K(2)<<"\n";
1615                            ffout(0)<<cc[0]<<" "<<cc[1]<<" "<<cc[2]<<"\n";
1616                        }
1617                        if(pointsin.find(pairTriangleToIntersectionPoints[pT].listofPoints[r])
1618             ==pointsin.end())
1619                            pointsin[(pairTriangleToIntersectionPoints[pT].listofPoints[r])]
1620                                =(pairTriangleToIntersectionPoints[pT].listofPoints[r]);
1621                  }                  }
                 pointsin.push_back((listVertexIntersectionPerPairEdges[objectp])[pT].listofPoints[r]);  
1622              }              }
1623              //ffout(0)<<cc[0]<<" "<<cc[1]<<" "<<cc[2]<<"\n";              //ffout(0)<<cc[0]<<" "<<cc[1]<<" "<<cc[2]<<"\n";
1624       //ffout(0)<<((listVertexIntersectionPerPair[objectp])[pT]).size()<<"\n";;       //ffout(0)<<((listVertexIntersectionPerPair[objectp])[pT]).size()<<"\n";;
1625          }              else {
         else {  
1626             // ffout(0) <<"pasintersecte\n";             // ffout(0) <<"pasintersecte\n";
1627                }
1628          }          }
1629      }      }
1630  }  }
1631    
1632  TinyVector<3,int> SurfaceMeshGenerator::Internals::  TinyVector<3,TinyVector<2,int> > SurfaceMeshGenerator::Internals::
1633  __determinateEdgeCut(const size_t& numobject,  __determinateEdgeCut(const size_t& numobject,
1634                       const size_t & objectNumber,                       const size_t & objectNumber,
1635                       const TinyVector<3,TinyVector<2,TinyVector<2,Edge::Pair> > >& currentTriangle,                       const Triangle*& currentTriangle,
1636                       const Cell* cell){                       const Cell* cell){
1637        
1638      TinyVector<3,int> tt;      TinyVector<3,TinyVector<2,int> > tt;
1639      tt=0;      tt=0;
     ObjectPair objectp(numobject,objectNumber);  
     if(numobject>objectNumber) {  
         objectp.first=objectNumber;  
         objectp.second=numobject;  
     }  
   
     for(SurfaceMeshGenerator::Internals::TriangleVerticesMap::iterator currentTriangle1  
         = (*triangleListesIntersectionEdgesTempPerHexa[cell][objectNumber]).begin();  
         currentTriangle1 != (*triangleListesIntersectionEdgesTempPerHexa[cell][objectNumber]).end();  
         ++currentTriangle1) {  
   
         const SurfaceMeshGenerator::Internals::TriangleVertices& tv=(*(*currentTriangle1).second);  
   
         PairTriangleEdges pT;  
         if(numobject<objectNumber) {  
             pT.first=&(currentTriangle);  
             pT.second=&(tv);  
         }  
         else {  
             pT.first=&(tv);  
             pT.second=&(currentTriangle);  
         }  
1640    
1641          if((listVertexIntersectionPerPairEdges[objectp]).find(pT)!=      ListOfHexaToTriangle::iterator hexait=listOfHexaToTriangle.begin();
1642      (listVertexIntersectionPerPairEdges[objectp]).end()) {      MapCellTriangle hexaToTriangle1;
1643            //  if((pT.first[0])[0][0]!=(pT.first[0])[1][0] or (pT.first[0])[0][1]!=(pT.first[0])[1][1])      if(numobject!=0)
1644            //    ffout(0)<<"un determinate new intersecte ============\n";          hexaToTriangle1=*(*hexait);
1645                    else {
1646              TinyVector<3,size_t> cc;          ++hexait;
1647            hexaToTriangle1=*(*hexait);
1648        }
1649        
1650        MapCellTriangle::iterator it1=hexaToTriangle1.find(cell);
1651        
1652        if(it1!=hexaToTriangle1.end()) {
1653            std::list<Triangle*> listTriangle1=(*it1).second;
1654            for(std::list<Triangle*>::iterator jt0=listTriangle1.begin() ; jt0!=listTriangle1.end() ; ++jt0) {
1655    
1656                PairTriangleNew pT;
1657              if(numobject<objectNumber) {              if(numobject<objectNumber) {
1658                  cc=((listVertexIntersectionPerPairEdges[objectp])[pT]).cutTriangle1;                  pT.first=(currentTriangle);
1659                    pT.second=(*jt0);
1660              }              }
1661              else {              else {
1662                  cc=((listVertexIntersectionPerPairEdges[objectp])[pT]).cutTriangle2;                  pT.second=(currentTriangle);
1663              }                  pT.first=(*jt0);
             if(cc[0]!=0) {  
                 tt[0]=1;  
1664              }              }
1665              if(cc[1]!=0) {  
1666                  tt[1]=1;              if(pairTriangleToIntersectionPoints.find(pT)!= pairTriangleToIntersectionPoints.end()) {
1667    
1668                    TinyVector<3,size_t >cc;
1669                    if(numobject<objectNumber) {
1670                        cc=(pairTriangleToIntersectionPoints[pT]).cutTriangle1;
1671                    }
1672                    else {
1673                        cc=(pairTriangleToIntersectionPoints[pT]).cutTriangle2;
1674                    }
1675                    if(cc[0]!=0) {
1676                        tt[0]=1;
1677                    }
1678                    if(cc[1]!=0) {
1679                        tt[1]=1;
1680                    }
1681                    if(cc[2]!=0) {
1682                        tt[2]=1;
1683                    }
1684              }              }
1685              if(cc[2]!=0) {              else {
1686                  tt[2]=1;                  // ffout(0)<<"la paire de triangle n'existe pas!\n";
1687              }              }
1688          }          }
         else {  
            // ffout(0)<<"la paire de triangle n'existe pas!\n";  
         }  
1689      }      }
1690      return tt;      return tt;
1691  }  }
1692    
1693  TinyVector<3,size_t> SurfaceMeshGenerator::Internals::  TinyVector<3,size_t> SurfaceMeshGenerator::Internals::
1694  __putEdgeRef(const size_t& ncase,  __putEdgeRef(const size_t& ncase,
1695               const TinyVector<3,TinyVector<2,TinyVector<2,Edge::Pair> > >& currentTriangle,               const Triangle& currentTriangle,
1696               TriangleCut& K) {               TriangleCut& K) {
1697    
1698      const TinyVector<3,int>  tt=K.edgecut;      const TinyVector<3,TinyVector<2,int> >&  tt=K.edgecut;
1699      TinyVector<3,size_t>& ii=K.isInTriangle;      TinyVector<3,size_t>& ii=K.isInTriangle;
1700      size_t numobject=K.numobject;      size_t numobject=K.numobject;
1701        const Vertex* V0=&currentTriangle(0);
1702        const Vertex* V1=&currentTriangle(1);
1703        const Vertex* V2=&currentTriangle(2);
1704    
1705      const TinyVector<2,Edge::Pair>& e0=(currentTriangle)[0][0];      
1706      const TinyVector<2,Edge::Pair>& e1=(currentTriangle)[1][0];      if(tt[0][0]+tt[1][0]+tt[2][0]!=0 and ((*edgesRefVertex[numobject])[V0]+
1707      const TinyVector<2,Edge::Pair>& e2=(currentTriangle)[2][0];                                   (*edgesRefVertex[numobject])[V1]+(*edgesRefVertex[numobject])[V2])==6) {
   
     if(tt[0]+tt[1]+tt[2]!=0 and ((*edgesRefCreate[numobject])[e0[0]]+  
                                  (*edgesRefCreate[numobject])[e1[0]]+(*edgesRefCreate[numobject])[e2[0]])==6) {  
1708          ffout(0)<<"************************ ---------tr a pbs\n";          ffout(0)<<"************************ ---------tr a pbs\n";
         printpoint(numobject,e0,K);  
         printpoint(numobject,e1,K);  
         printpoint(numobject,e2,K);  
1709          ffout(0)<<tt[0]<<" "<<tt[1]<<" "<<tt[2]<<"\n";          ffout(0)<<tt[0]<<" "<<tt[1]<<" "<<tt[2]<<"\n";
1710      }      }
1711    
1712        __inout(ncase,numobject,tt[0][0],V0,V1,K);
1713        __inout(ncase,numobject,tt[1][0],V0,V2,K);
1714        __inout(ncase,numobject,tt[2][0],V1,V2,K);
1715            
     __inout(numobject,tt[0],e0,e1,ii[0],ii[1],K);  
     __inout(numobject,tt[1],e0,e2,ii[0],ii[2],K);  
     __inout(numobject,tt[2],e1,e2,ii[1],ii[2],K);  
1716      if(ncase==0) {      if(ncase==0) {
1717          if(tt[0]+tt[1]+tt[2]==0 and ((*edgesRefCreate[numobject])[e0[0]]+          if(tt[0][0]+tt[1][0]+tt[2][0]==0 and ((*edgesRefVertex[numobject])[V0]+
1718                                (*edgesRefCreate[numobject])[e1[0]]+(*edgesRefCreate[numobject])[e2[0]])>=4) {                                (*edgesRefVertex[numobject])[V1]+(*edgesRefVertex[numobject])[V2])>=4) {
1719              //if triangle pas plat.........              (*edgesRefVertex[numobject])[V0]=2;
1720             // if(!(triangleplat(numobject,e0,e1,e2)) ) {              (*edgesRefVertex[numobject])[V1]=2;
1721                  (*edgesRefCreate[numobject])[e0[0]]=2;              (*edgesRefVertex[numobject])[V2]=2;
                 (*edgesRefCreate[numobject])[e1[0]]=2;  
                 (*edgesRefCreate[numobject])[e2[0]]=2;  
            // }  
             /*else {  
                 real_t eps=1e-13;  
                 Vertex V0,V1,V2;  
                 V0=(*(*verticesListesIntersection[numobject])[e0[0]]);  
                 V1=(*(*verticesListesIntersection[numobject])[e1[0]]);  
                 V2=(*(*verticesListesIntersection[numobject])[e2[0]]);  
                 real_t normV0V1=Norm(V0.Vector3()-V1.Vector3());  
                 real_t normV0V2=Norm(V0.Vector3()-V2.Vector3());  
                 real_t normV1V2=Norm(V2.Vector3()-V1.Vector3());  
                 if(normV0V1>eps and (*edgesRefCreate[numobject])[e0[0]]+  
      (*edgesRefCreate[numobject])[e1[0]]==4)  
                     (*edgesRefCreate[numobject])[e2[0]]=2;  
                 if(normV0V2>eps and (*edgesRefCreate[numobject])[e0[0]]+  
      (*edgesRefCreate[numobject])[e2[0]]==4)  
                     (*edgesRefCreate[numobject])[e1[0]]=2;  
                 if(normV1V2>eps and (*edgesRefCreate[numobject])[e1[0]]+  
      (*edgesRefCreate[numobject])[e2[0]]==4)  
                     (*edgesRefCreate[numobject])[e0[0]]=2;  
                 ffout(0)<<"triangle plat\n";  
             }*/  
1722          }                }      
1723      }      }
1724      else if(ncase==1) {      else if(ncase==1 or ncase==0) {
1725            
1726          if((*edgesRefCreate[numobject])[e0[0]]==2)          if((*edgesRefVertex[numobject])[V0]==2 )
1727              ii[0]=1;              ii[0]=1;
1728          if((*edgesRefCreate[numobject])[e1[0]]==2)          else
1729                ii[0]=0;
1730            if((*edgesRefVertex[numobject])[V1]==2 )
1731              ii[1]=1;              ii[1]=1;
1732          if((*edgesRefCreate[numobject])[e2[0]]==2)          else
1733                ii[1]=0;
1734            if((*edgesRefVertex[numobject])[V2]==2 )
1735              ii[2]=1;              ii[2]=1;
1736            else
1737                ii[2]=0;
1738      }      }
1739        
1740    
# Line 2161  __putEdgeRef(const size_t& ncase, Line 1743  __putEdgeRef(const size_t& ncase,
1743    
1744  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1745  __createListTriangle(const TriangleCut& K,  __createListTriangle(const TriangleCut& K,
1746                       const Cell* cell,                       const Cell* cell){
1747                       size_t & kill){      const TinyVector<3,TinyVector<2,int> >&  tt=K.edgecut;
     const TinyVector<3,int>&  tt=K.edgecut;  
1748      const TinyVector<3,size_t>& ii=K.isInTriangle;      const TinyVector<3,size_t>& ii=K.isInTriangle;
1749      const size_t& numobject=K.numobject;      const Vertex& V0=K(0);
1750      const TinyVector<2,TinyVector<2,Edge::Pair> >& e0=K.listpairedges[0];      const Vertex& V1=K(1);
1751      const TinyVector<2,TinyVector<2,Edge::Pair> >& e1=K.listpairedges[1];      const Vertex& V2=K(2);
1752      const TinyVector<2,TinyVector<2,Edge::Pair> >& e2=K.listpairedges[2];    
1753      TinyVector<2,TinyVector<2,Edge::Pair> > pe0(e0),pe1(e1),pe2(e2);     if(((tt[0][0]+tt[1][0]+tt[2][0]==0 and ii[0]+ii[1]+ii[2]==3) /*or ii[0]+ii[1]+ii[2]==3*/)) {
1754           __addTriangle(V0,V1,V2,cell);
    if(/*kill==3 and*/ ((tt[0]+tt[1]+tt[2]==0 and ii[0]+ii[1]+ii[2]==3) /*or ii[0]+ii[1]+ii[2]==3*/)) {  
        TinyVector<2,TinyVector<2,Edge::Pair> > num0=pe0,num1=pe1,num2=pe2;  
         __addTriangle(numobject,num0,num1,num2,cell,kill);  
1755      }      }
1756      else if(ii[0]+ii[1]+ii[2]==1) {      else if(ii[0]+ii[1]+ii[2]==1) {
1757          if(ii[0]==1) {          if(ii[0]==1) {
1758              __createCase1(0,K,cell,kill);              __createCase1(0,K,cell);
1759          }          }
1760          if(ii[1]==1) {          if(ii[1]==1) {
1761              __createCase1(1,K,cell,kill);              __createCase1(1,K,cell);
1762          }          }
1763          if(ii[2]==1) {          if(ii[2]==1) {
1764              __createCase1(2,K,cell,kill);              __createCase1(2,K,cell);
1765          }          }
1766      }      }
1767      else if(ii[0]+ii[1]+ii[2]==2) {      else if(ii[0]+ii[1]+ii[2]==2) {
1768          if(ii[0]+ii[1]==2) {          if(ii[0]+ii[1]==2) {
1769              __createCase2(0,K,cell,kill);              __createCase2(0,K,cell);
1770          }          }
1771          if(ii[0]+ii[2]==2) {          if(ii[0]+ii[2]==2) {
1772              __createCase2(1,K,cell,kill);              __createCase2(1,K,cell);
1773          }          }
1774          if(ii[1]+ii[2]==2) {          if(ii[1]+ii[2]==2) {
1775              __createCase2(2,K,cell,kill);              __createCase2(2,K,cell);
1776            }
1777        }
1778    }
1779    
1780    void SurfaceMeshGenerator::Internals::
1781    __createNotIn(const TriangleCut& K,const Cell* cell) {
1782        size_t place0=0;
1783        if(K.edgecut[1][0]==1 and K.edgecut[1][1]==1) {
1784            place0=1;
1785        }
1786        else if(K.edgecut[2][0]==1 and K.edgecut[2][1]==1) {
1787            place0=2;
1788        }
1789        const Vertex* Pend=K.pointsIntersection[place0][1];
1790        TinyVector<10,const Vertex*> P;
1791        TinyVector<10,const Triangle*> T;
1792        P[0]=K.pointsIntersection[place0][0];
1793        if(vertexInTriangles.find(P[0])!=vertexInTriangles.end()) {
1794            T[0]=(*vertexInTriangles.find(P[0])).second[0];
1795        }
1796        else {
1797            ffout(0)<<"pbs 0 existe pas\n";
1798        }
1799        size_t b=0;
1800        while(P[b]!=Pend and b<9) {
1801            if(triangleWithVertex.find(T[b])!=triangleWithVertex.end()) {
1802                TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[b])).second;
1803                if(VV[0]==P[b]) {
1804                    if(Pend!=VV[1])
1805                        __addTriangle(*Pend,*P[b],*VV[1],cell);
1806                    P[b+1]=VV[1];
1807                }
1808                else {
1809                    if(Pend!=VV[0])
1810                        __addTriangle(*Pend,*P[b],*VV[0],cell);
1811                    P[b+1]=VV[0];
1812                }
1813            }
1814            else
1815                ffout(0)<<"pbs 1 existe pas\n";
1816            TinyVector<2,const Triangle*> T2;
1817            if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {
1818                T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];
1819                T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];
1820            }
1821            if(T2[0]==T[b]) {
1822                if(T2[1]==T[b]) {
1823                    //point edge
1824                }
1825                else
1826                    T[b+1]=T2[1];
1827            }
1828            else {
1829                T[b+1]=T2[0];
1830            }
1831            ++b;
1832        }
1833    }
1834    
1835    
1836    void SurfaceMeshGenerator::Internals::
1837    __createByEdges(const TriangleCut& K,const Cell* cell) {
1838    
1839        size_t in0=0,place0=0,place1=1,place2=2;
1840        if(K.isInTriangle[1]==1) {
1841            in0=1;
1842            place1=2;
1843            place2=1;
1844        }
1845        else if(K.isInTriangle[2]==1) {
1846            in0=2;
1847            place0=1;
1848            place1=2;
1849            place2=0;
1850        }
1851        const Vertex* Pend=K.pointsIntersection[place1][0];
1852        TinyVector<10,const Vertex*> P;
1853        TinyVector<10,const Triangle*> T;
1854        P[0]=K.pointsIntersection[place0][0];
1855        if(vertexInTriangles.find(P[0])!=vertexInTriangles.end()) {
1856            T[0]=(*vertexInTriangles.find(P[0])).second[0];
1857        }
1858        else {
1859            ffout(0)<<"pbs 0 existe pas\n";
1860        }
1861        size_t b=0;
1862        while(P[b]!=Pend and b<9) {
1863            if(triangleWithVertex.find(T[b])!=triangleWithVertex.end()) {
1864                TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[b])).second;
1865                if(VV[0]==P[b]) {
1866                    if(VV[1]==P[b]) {
1867                        ffout(0)<<"egal\n";
1868                        if(K.pointsIntersection[place2][0]==P[b]) {
1869                            P[b+1]=K.pointsIntersection[place2][1];
1870                            __addTriangle(K(in0),*P[b],*P[b+1],cell);
1871                        }
1872                        else if(K.pointsIntersection[place2][1]==P[b]){
1873                            P[b+1]=K.pointsIntersection[place2][0];
1874                            __addTriangle(K(in0),*P[b],*P[b+1],cell);
1875                        }
1876                        else {
1877                            ffout(0)<<"pbs pbs\n";
1878                            ffout(0)<<*P[b]<<"\n";
1879                            std::exit(0);
1880                        }
1881                    }
1882                    else {
1883                    __addTriangle(K(in0),*P[b],*VV[1],cell);
1884                    P[b+1]=VV[1];
1885                    }
1886                }
1887                else {
1888                    __addTriangle(K(in0),*P[b],*VV[0],cell);
1889                    P[b+1]=VV[0];
1890                }
1891            }
1892            else
1893                ffout(0)<<"pbs 1 existe pas\n";
1894            TinyVector<2,const Triangle*> T2;
1895            if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {
1896                T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];
1897                T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];
1898            }
1899            else {
1900                ffout(0)<<"pbs tr existe pas\n";
1901          }          }
1902            if(T2[0]==T[b] and T2[1]==T[b]) {
1903               // ffout(0)<<"here0\n";
1904                //point edge
1905                if(P[b+1]!=Pend) {
1906                    if(K.pointsIntersection[place2][0]==P[b+1]) {
1907                        const Vertex* A=K.pointsIntersection[place2][1];
1908                        __addTriangle(K(in0),*P[b+1],*A,cell);
1909                        if(vertexInTriangles.find(A)!=vertexInTriangles.end()) {
1910                            T[b+1]=(*vertexInTriangles.find(A)).second[0];
1911                            P[b+1]=A;
1912                        }
1913                    }
1914                    else {
1915                        const Vertex* A=K.pointsIntersection[place2][0];
1916                        __addTriangle(K(in0),*P[b+1],*A,cell);
1917                        if(vertexInTriangles.find(A)!=vertexInTriangles.end()) {
1918                            T[b+1]=(*vertexInTriangles.find(A)).second[0];
1919                            P[b+1]=A;
1920                        }
1921                    }
1922                }
1923            }
1924            else if(T2[0]==T[b]) {
1925                T[b+1]=T2[1];
1926            }
1927            else {
1928                T[b+1]=T2[0];
1929            }
1930            ++b;
1931      }      }
1932        
1933  }  }
1934    
1935    
1936  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1937  __createTrianglesIntersection(const size_t& numobject,  __createTrianglesIntersection(const size_t& numobject,
1938                                const size_t& objectNumber,                                const size_t& objectNumber) {
                               const size_t& nobject) {  
1939    
1940      ffout(0)<<"numero "<<numobject<<" "<<objectNumber<<"\n";      ffout(0)<<"numero "<<numobject<<" "<<objectNumber<<"\n";
1941    
1942      for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)      ListOfHexaToTriangle::iterator hexait=listOfHexaToTriangle.begin();
1943      {      MapCellTriangle hexaToTriangle0;
1944          for(SurfaceMeshGenerator::Internals::TriangleVerticesMap::iterator currentTriangle      if(numobject==0)
1945       = (*triangleListesIntersectionEdgesTempPerHexa[*i][numobject]).begin();          hexaToTriangle0=*(*hexait);
1946       currentTriangle != (*triangleListesIntersectionEdgesTempPerHexa[*i][numobject]).end();      else {
1947       ++currentTriangle) {          ++hexait;
1948              const SurfaceMeshGenerator::Internals::TriangleVertices& tv=(*(*currentTriangle).second);          hexaToTriangle0=*(*hexait);
1949              const SurfaceMeshGenerator::Internals::TriangleVertices& tv2=(*(*currentTriangle).first);      }
1950              const TinyVector<2,Edge::Pair>& e0=tv2[0][0];  //Initialize edgesRefVertex
1951              const TinyVector<2,Edge::Pair>& e1=tv2[1][0];      size_t taille=0;
1952              const TinyVector<2,Edge::Pair>& e2=tv2[2][0];      for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i) {
1953            MapCellTriangle::iterator it0=hexaToTriangle0.find(*i);
1954              (*edgesRefCreate[numobject])[e0[0]]=0;          if(it0!=hexaToTriangle0.end()) {
1955              (*edgesRefCreate[numobject])[e1[0]]=0;              std::list<Triangle*> listTriangle0=(*it0).second;
1956              (*edgesRefCreate[numobject])[e2[0]]=0;              taille +=listTriangle0.size();
1957                for(std::list<Triangle*>::iterator jt0=listTriangle0.begin() ;
1958             jt0!=listTriangle0.end() ; ++jt0) {
1959                    const Triangle& T = *(*jt0);
1960                    
1961                    (*edgesRefVertex[numobject])[&(T(0))]=0;
1962                    (*edgesRefVertex[numobject])[&(T(1))]=0;
1963                    (*edgesRefVertex[numobject])[&(T(2))]=0;
1964                }
1965          }          }
1966      }      }
1967    
1968      for(SurfaceMeshGenerator::Internals::EdgePairList::const_iterator currentTriangle      ListOfTriangleList::iterator triangleit=listOfTriangleList.begin();
1969          = (*triangleListes[numobject]).begin();      if(numobject!=0) {
1970          currentTriangle != (*triangleListes[numobject]).end();          ++triangleit;
1971          ++currentTriangle) {      }
1972          const Edge& e0=(*currentTriangle)[0];      for(size_t size=0 ; size<(*(*triangleit)).size() ; ++size)  {
1973          const Edge& e1=(*currentTriangle)[1];          const Triangle& T=(*(*triangleit))[size];
1974          const Edge& e2=(*currentTriangle)[2];          (*edgesRefVertex[numobject])[&T(0)]=2;
1975          (*edgesRefCreate[numobject])[e0]=2;          (*edgesRefVertex[numobject])[&T(1)]=2;
1976          (*edgesRefCreate[numobject])[e1]=2;          (*edgesRefVertex[numobject])[&T(2)]=2;
         (*edgesRefCreate[numobject])[e2]=2;  
     }  
   
     ObjectPair objectp(numobject,objectNumber);  
     if(numobject>objectNumber) {  
         objectp.first=objectNumber;  
         objectp.second=numobject;  
1977      }      }
1978      size_t nbinters=listVertexIntersectionPerPairEdges[objectp].size();      
1979      ffout(0)<<"size listVertexIntersectionPerPair "<<nbinters<<"\n";     size_t nbinters=listVertexIntersectionNew.size();
1980        ffout(0)<<"size listVertexIntersectionNew "<<nbinters<<"\n";
1981    
1982      size_t dd_max=3;      size_t dd_max=4;//4
1983        ffout(0)<<"size liste intersection "<<taille<<"\n";
1984      for(size_t dd=0 ; dd<dd_max ; ++dd) {      for(size_t dd=0 ; dd<dd_max ; ++dd) {
1985          //ffout(0)<<"num passage "<<dd<<"\n";          ffout(0)<<"num passage "<<dd<<"\n";
           
1986          for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)          for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)
1987          {          {
1988              for(SurfaceMeshGenerator::Internals::TriangleVerticesMap::iterator currentTriangle              MapCellTriangle::iterator it0=hexaToTriangle0.find(*i);
1989           = (*triangleListesIntersectionEdgesTempPerHexa[*i][numobject]).begin();              if(it0!=hexaToTriangle0.end()) {
1990           currentTriangle != (*triangleListesIntersectionEdgesTempPerHexa[*i][numobject]).end();                  std::list<Triangle*> listTriangle0=(*it0).second;
1991           ++currentTriangle) {                  for(std::list<Triangle*>::iterator jt0=listTriangle0.begin()
1992          ; jt0!=listTriangle0.end() ; ++jt0) {
1993                  const SurfaceMeshGenerator::Internals::TriangleVertices& tv=(*(*currentTriangle).second);                      const Triangle* T=(*jt0);
1994                  const SurfaceMeshGenerator::Internals::TriangleVertices& tv2=(*(*currentTriangle).first);                      const Vertex V0=(*T)(0);
1995                  const TinyVector<2,Edge::Pair>& e0=tv2[0][0];                      const Vertex V1=(*T)(1);
1996                  const TinyVector<2,Edge::Pair>& e1=tv2[1][0];                      const Vertex V2=(*T)(2);
1997                  const TinyVector<2,Edge::Pair>& e2=tv2[2][0];                      TriangleCut K(V0,V1,V2,numobject,objectNumber);
   
                 TriangleCut K(e0,e1,e2,numobject,objectNumber);  
1998    
1999                  K.edgecut =__determinateEdgeCut(numobject,objectNumber,(*(*currentTriangle).second),*i);                      K.edgecut =__determinateEdgeCut(numobject,objectNumber,T,*i);
2000                  K.isInTriangle =__putEdgeRef(0,(*(*currentTriangle).second),K);                      K.isInTriangle =__putEdgeRef(0,(*T),K);
2001                        }
2002              }              }
2003          }          }
2004      }      }  
   
     size_t kill=0;  
2005    
2006        size_t err=0,err2=0;
2007      for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)      for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)
2008      {      {
2009          for(SurfaceMeshGenerator::Internals::TriangleVerticesMap::iterator currentTriangle          MapCellTriangle::iterator it0=hexaToTriangle0.find(*i);
2010       = (*triangleListesIntersectionEdgesTempPerHexa[*i][numobject]).begin();                  
2011       currentTriangle != (*triangleListesIntersectionEdgesTempPerHexa[*i][numobject]).end();          if(it0!=hexaToTriangle0.end()) {
2012       ++currentTriangle) {              std::list<Triangle*> listTriangle0=(*it0).second;
2013              const SurfaceMeshGenerator::Internals::TriangleVertices& tv=(*(*currentTriangle).second);              for(std::list<Triangle*>::iterator jt0=listTriangle0.begin() ;
2014              const SurfaceMeshGenerator::Internals::TriangleVertices& tv2=(*(*currentTriangle).first);           jt0!=listTriangle0.end() ; ++jt0) {
2015              TinyVector<2,TinyVector<2,Edge::Pair> > e0,e1,e2;                  const Triangle* T=(*jt0);
2016              if(tv2[0][0]==tv2[0][1]) {                  const Vertex& V0=(*T)(0);
2017                  e0=tv[0];}                  const Vertex& V1=(*T)(1);
2018              else {                  const Vertex& V2=(*T)(2);
2019                  e0=tv2[0];                  TriangleCut K(V0,V1,V2,numobject,objectNumber);
2020              }                  triangleWithVertex.clear();
2021              if(tv2[1][0]==tv2[1][1]) {                  vertexInTriangles.clear();
2022                  e1=tv[1];}  
2023              else {                  __createLocalListIntersection(objectNumber,T,K,*i);
2024                  e1=tv2[1];                  K.isInTriangle =__putEdgeRef(1,*T,K);
2025              }                  if(K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==0 ) {
2026              if(tv2[2][0]==tv2[2][1]) {                      ++err;
2027                  e2=tv[2];}                  }
2028              else {                  if(K.edgecut[0]+K.edgecut[1]+K.edgecut[2]==0 ) {
2029                  e2=tv2[2];                      ++err2;
2030              }                  }
2031                    if(K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==0
2032              TriangleCut K(e0,e1,e2,numobject,objectNumber);       and K.edgecut[0]+K.edgecut[1]+K.edgecut[2]!=0) {
2033                        if(vertexInTriangles.size()==3) {
2034              __createLocalListIntersection(objectNumber,(*(*currentTriangle).second),K,*i);                          TinyVector<3,const Vertex*> V;
2035              K.isInTriangle =__putEdgeRef(1,(*(*currentTriangle).second),K);                          size_t num=0;
2036                                      for(std::map<const Vertex*,TinyVector<2,const Triangle*> >::iterator
2037              kill=0;         it=vertexInTriangles.begin() ; it!=vertexInTriangles.end() ; ++it , ++num) {
2038                                V[num]=(*it).first;
2039              if(nbinters!=0) {                          }
2040                  //if(!(triangleplat(numobject,e0,e1,e2)))                          if(V[0]!=V[1] and V[1]!=V[2] and V[2]!=V[0]) {
2041                  __createListTriangle(K,*i,kill);                              __addTriangle(*V[0],*V[1],*V[2],*i);
2042              }                              ffout(0)<<"on comble le trou\n";
2043              else {                          }
2044                  //ffout(0)<<"j ajoute le tr car pas d'intersection dans les 2 listes\n";                      }
2045                  __addTriangle(numobject,e0,e1,e2,*i,kill);                      else if(vertexInTriangles.size()>3 and
2046                  K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==1){
2047                            __createNotIn(K,*i);
2048                        }
2049                        else if(K.edgecut[0]!=0) {
2050                            if(K.pointsIntersection[0][0]==&V0
2051          or K.pointsIntersection[0][0]==&V1
2052          or K.pointsIntersection[0][0]==&V2) {
2053                                ffout(0)<<"**********normal\n";
2054                            }
2055                            else {
2056                                ffout(0)<<"triangle cut but outside\n";
2057                                ffout(0)<<V0<<"\n";
2058                                ffout(0)<<V1<<"\n";
2059                                ffout(0)<<V2<<"\n";
2060                                ffout(0)<<"nb points inter "<<K.pointsIntersection.size()<<"\n";
2061                            }
2062                        }
2063                        else if(K.edgecut[1]!=0) {
2064                            if(K.pointsIntersection[1][0]==&V0
2065          or K.pointsIntersection[1][0]==&V1
2066          or K.pointsIntersection[1][0]==&V2) {
2067                                ffout(0)<<"**********normal\n";
2068                            }
2069                            else {
2070                                ffout(0)<<"triangle cut but outside\n";
2071                                ffout(0)<<V0<<"\n";
2072                                ffout(0)<<V1<<"\n";
2073                                ffout(0)<<V2<<"\n";
2074                                ffout(0)<<"nb points inter "<<K.pointsIntersection.size()<<"\n";
2075                            }
2076                        }
2077                        else if(K.edgecut[2]!=0) {
2078                            if(K.pointsIntersection[2][0]==&V0
2079          or K.pointsIntersection[2][0]==&V1
2080          or K.pointsIntersection[2][0]==&V2) {
2081                                ffout(0)<<"**********normal\n";
2082                            }
2083                            else {
2084                                ffout(0)<<"triangle cut but outside\n";
2085                                ffout(0)<<V0<<"\n";
2086                                ffout(0)<<V1<<"\n";
2087                                ffout(0)<<V2<<"\n";
2088                                ffout(0)<<"nb points inter "<<K.pointsIntersection.size()<<"\n";
2089                            }
2090                        }
2091                    }
2092                    if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==3
2093         and K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==1) {
2094                        if(K.pointsIntersection[0][0]!=K.pointsIntersection[1][0]
2095             and K.pointsIntersection[1][0]!=K.pointsIntersection[2][0]
2096             and K.pointsIntersection[2][0]!=K.pointsIntersection[0][0]) {
2097                            ffout(0)<<"---------------- 3 edge cut \n";
2098                            ffout(0)<<V0<<"\n";
2099                            ffout(0)<<V1<<"\n";
2100                            ffout(0)<<V2<<"\n";
2101                        __createByEdges(K,*i);
2102                        }
2103                        else {
2104                            __createListTriangle(K,*i);
2105                        }
2106                    }
2107                    else if(nbinters!=0) {
2108                            //if(!(triangleplat(numobject,e0,e1,e2)))
2109                        __createListTriangle(K,*i);
2110                    }
2111                    else {
2112                            //ffout(0)<<"j ajoute le tr car pas d'intersection dans les 2 listes\n";
2113                            __addTriangle(V0,V1,V2,*i);
2114                    }
2115              }              }
2116          }          }
2117      }      }
2118        if(err==taille)
2119            ffout(0)<<"tous les triangles de la liste d'intersection sont dehors\n";
2120        if(err2==taille)
2121            ffout(0)<<"aucun triangle de la liste d'intersection n'est coupe\n";
2122    
2123    
2124        (*edgesRefVertex[numobject]).clear();
2125  }  }
2126    
2127  bool SurfaceMeshGenerator::Internals::  bool SurfaceMeshGenerator::Internals::
# Line 2664  __createSurface(Structured3DMesh& SMesh, Line 2472  __createSurface(Structured3DMesh& SMesh,
2472            }            }
2473                            
2474            if(jobDone) {            if(jobDone) {
2475                meshedHexahedra.insert(&currentCell);               // meshedHexahedra.insert(&currentCell);
2476                  SurfaceMeshGenerator::Internals::MotherCellList::iterator itcell
2477                  =cellList.begin();
2478                for (SurfaceMeshGenerator::Internals::EdgePairList::iterator t                for (SurfaceMeshGenerator::Internals::EdgePairList::iterator t
2479              = localTriangleList.begin();              = localTriangleList.begin();
2480              t != localTriangleList.end(); ++t) {              t != localTriangleList.end(); ++t) {
2481                    (*triangleListes[nobject]).push_front(*t);                    (*triangleListes[nobject]).push_front(*t);
2482                      (*cellListObject[nobject]).push_front(*itcell);
2483                      ++itcell;
2484                }                }
2485                  hexaTetraAssociation[&currentCell]
2486                  = std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,
2487                  ReferenceCounting<std::vector<Edge> > >
2488                  (pT,pAddedEdges);
2489            } else {            } else {
2490                hexaTetraAssociation[&currentCell]                hexaTetraAssociation[&currentCell]
2491                = std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,                = std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,
# Line 2677  __createSurface(Structured3DMesh& SMesh, Line 2493  __createSurface(Structured3DMesh& SMesh,
2493                (pT,pAddedEdges);                (pT,pAddedEdges);
2494            }            }
2495            localTriangleList.clear();            localTriangleList.clear();
2496              cellList.clear();
2497          }          }
2498        }        }
2499    
2500  }  }
2501    
 void SurfaceMeshGenerator::Internals::  
 __createMeshIntersection(size_t nobjectToTreat,  
                          const Scene & scene) {  
   
     size_t nobject=scene.nbObjects();  
     for(size_t no=0 ; no<nobjectToTreat ; ++no) {  
         size_t begin=2;  
         if(no>=2)  
             begin=no+1;  
         for(size_t no2=begin ; no2<nobjectToTreat ; ++no2) {  
             if(no==0)  
                 __calculatePointsIntersection(scene,nobjectToTreat,nobject+no,no2);  
             else  
                 __calculatePointsIntersection(scene,nobjectToTreat,nobject+no,nobject+no2);  
   
             for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ;  
          i!=toTreatHexahedra.end() ; ++i)  
             {  
                 (*triangleListesIntersectionEdgesTempPerHexa[*i][no]).clear();  
                 (*listTriangleintersection[no]).clear();  
                 (*triangleListesIntersectionEdgesTempPerHexa[*i][no2]).clear();  
                 (*listTriangleintersection[no2]).clear();  
   
                 for(SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::iterator currentTriangle1  
       = (*triangleListesIntersectionPerHexa[*i][no]).begin();  
       currentTriangle1 != (*triangleListesIntersectionPerHexa[*i][no]).end();  
       ++currentTriangle1) {  
                     (*  
        triangleListesIntersectionEdgesTempPerHexa[*i][no])[&(*currentTriangle1)]=&(*currentTriangle1);  
   
                 }  
                 if(no==0) {  
                     for(SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::iterator currentTriangle1  
           = (*triangleListesIntersectionEdgesPerHexa[*i][no2]).begin();  
           currentTriangle1 != (*triangleListesIntersectionEdgesPerHexa[*i][no2]).end();  
           ++currentTriangle1) {  
                         (*  
     triangleListesIntersectionEdgesTempPerHexa[*i][no2])[&(*currentTriangle1)]=&(*currentTriangle1);  
                     }  
                 }  
                 else {  
                     for(SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::iterator currentTriangle1  
           = (*triangleListesIntersectionPerHexa[*i][no2]).begin();  
           currentTriangle1 != (*triangleListesIntersectionPerHexa[*i][no2]).end();  
           ++currentTriangle1) {  
                         (*  
     triangleListesIntersectionEdgesTempPerHexa[*i][no2])[&(*currentTriangle1)]=&(*currentTriangle1);  
   
                     }  
                 }  
             }  
             __createTrianglesIntersection(no,no2,nobjectToTreat);  
             __createTrianglesIntersection(no2,no,nobjectToTreat);  
         }  
     }      
 }  
   
2502  void  void
2503  SurfaceMeshGenerator::  SurfaceMeshGenerator::
2504  generatesAllSufaceMeshes(const Problem& problem,  generatesAllSufaceMeshes(const Problem& problem,
# Line 2962  void plotmedit(const size_t& nobject, Line 2722  void plotmedit(const size_t& nobject,
2722                 const SurfaceMeshOfTriangles& s_mesh,                 const SurfaceMeshOfTriangles& s_mesh,
2723                 std::set<const Cell*>& toTreatHexahedra,                 std::set<const Cell*>& toTreatHexahedra,
2724                 std::set<const Cell*>& cuttedHexahedra,                 std::set<const Cell*>& cuttedHexahedra,
2725                 std::vector<ReferenceCounting<std::map<Edge::Pair, Vertex*> > > verticesListesIntersection,                 size_t nobjother,size_t n0,size_t n1);
                std::vector<ReferenceCounting<std::list<TinyVector<3,TinyVector<2,TinyVector<2,Edge::Pair> > > > > >  
                triangleListesIntersectionEdges,  
                size_t nobjother,size_t ninter,size_t ninter2,size_t ntr1);  
2726    
2727  void plothexa(size_t ncase,std::set<const Cell*>& cuttedHexahedra);  void plothexa(size_t ncase,std::set<const Cell*>& cuttedHexahedra);
2728    
# Line 3021  generateSurfacicMesh(const Object& O, Line 2778  generateSurfacicMesh(const Object& O,
2778    
2779    size_t nobject = scene.nbObjects();    size_t nobject = scene.nbObjects();
2780    size_t nobjectToTreat = otherReferences.size();    size_t nobjectToTreat = otherReferences.size();
2781    //ffout(0)<<"size other ref "<<nobjectToTreat<<"\n";    ffout(0)<<"size other ref "<<nobjectToTreat<<"\n";
2782    
2783    (*__internals).verticesListes.resize(nobject);    (*__internals).verticesListes.resize(nobject);
2784    (*__internals).triangleListes.resize(nobject);    (*__internals).triangleListes.resize(nobject);
   (*__internals).triangleListesIntersectionEdges.resize(nobject);  
2785    (*__internals).localTriangleListesEdges.resize(nobject);    (*__internals).localTriangleListesEdges.resize(nobject);
   (*__internals).localTriangleList2Edges.resize(nobject);  
2786    (*__internals).verticesListesIntersection.resize(nobject);    (*__internals).verticesListesIntersection.resize(nobject);
2787    (*__internals).edgesRef.resize(nobject);    (*__internals).edgesRef.resize(nobject);
2788    (*__internals).edgesRefCreate.resize(nobject);    (*__internals).cellListObject.resize(nobject);
2789    (*__internals).listTriangleintersection.resize(nobject);    (*__internals).edgesRefVertex.resize(nobject);
2790    
       
2791    for(size_t no=0 ; no<nobject ; ++no) {    for(size_t no=0 ; no<nobject ; ++no) {
2792        (*__internals).triangleListes[no]= new SurfaceMeshGenerator::Internals::EdgePairList;        (*__internals).triangleListes[no]= new SurfaceMeshGenerator::Internals::EdgePairList;
2793        (*__internals).triangleListesIntersectionEdges[no]=        (*__internals).localTriangleListesEdges[no]= new SurfaceMeshGenerator::Internals::EdgePairList;
           new SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList;  
       (*__internals).localTriangleListesEdges[no]= new SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList;  
       (*__internals).localTriangleList2Edges[no]=  
           new SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList2;  
2794        (*__internals).verticesListes[no]= new SurfaceMeshGenerator::Internals::VerticesList;        (*__internals).verticesListes[no]= new SurfaceMeshGenerator::Internals::VerticesList;
2795        (*__internals).verticesListesIntersection[no]= new SurfaceMeshGenerator::Internals::VerticesList;        (*__internals).verticesListesIntersection[no]= new SurfaceMeshGenerator::Internals::VerticesList;
2796          (*__internals).cellListObject[no]= new SurfaceMeshGenerator::Internals::MotherCellList;
2797    
2798        (*__internals).edgesRef[no]= new std::map<Edge::Pair, size_t> ;        (*__internals).edgesRef[no]= new std::map<Edge::Pair, size_t> ;
2799        (*__internals).edgesRefCreate[no]= new std::map<Edge::Pair, size_t> ;        (*__internals).edgesRefVertex[no] = new std::map<const Vertex*,size_t> ;
       (*__internals).listTriangleintersection[no]=  
           new SurfaceMeshGenerator::Internals::TriangleVerticesEdges;  
2800    }    }
2801    
2802        
# Line 3076  generateSurfacicMesh(const Object& O, Line 2825  generateSurfacicMesh(const Object& O,
2825                                     cuttedHexahedra,meshedHexahedra);                                     cuttedHexahedra,meshedHexahedra);
2826                }                }
2827            }            }
2828            plothexa(num,meshedHexahedra);            plothexa(num,cuttedHexahedra);
2829            if(num==0)            if(num==0)
2830                plothexa(num-1,cuttedHexahedra);                plothexa(num-1,cuttedHexahedra);
2831            num++;            num++;
2832        }        }
2833    }    }
2834    
2835    std::set_difference (cuttedHexahedra.begin(), cuttedHexahedra.end(),  std::set_difference (cuttedHexahedra.begin(), cuttedHexahedra.end(),
2836                       meshedHexahedra.begin(), meshedHexahedra.end(),                       meshedHexahedra.begin(), meshedHexahedra.end(),
2837                       inserter((*__internals).toTreatHexahedra, (*__internals).toTreatHexahedra.begin()));                       inserter((*__internals).toTreatHexahedra, (*__internals).toTreatHexahedra.begin()));
2838    
2839    std::set<const Cell*>::iterator iterat=(*__internals).toTreatHexahedra.begin();  //  ffout(0)<<(*__internals).toTreatHexahedra.size()<<"\n";
2840    //  (*__internals).toTreatHexahedra=meshedHexahedra;
2841    /*  std::set<const Cell*>::iterator iterat=(*__internals).toTreatHexahedra.begin();
2842    size_t treatsize=(*__internals).toTreatHexahedra.size();    size_t treatsize=(*__internals).toTreatHexahedra.size();
2843    for(size_t yy=0 ; yy<treatsize ; ++yy)    for(size_t yy=0 ; yy<treatsize ; ++yy)
2844    {    {
# Line 3121  generateSurfacicMesh(const Object& O, Line 2872  generateSurfacicMesh(const Object& O,
2872        if(test!=2) {        if(test!=2) {
2873            (*__internals).toTreatHexahedra.erase(temp);            (*__internals).toTreatHexahedra.erase(temp);
2874        }        }
2875    }    }*/
2876    //ffout(0)<<"nb tr objet 0 complet "<<(*(*__internals).triangleListes[0]).size()<<" \n";    ffout(0)<<"nb tr objet 0 complet "<<(*(*__internals).triangleListes[0]).size()<<" \n";
2877          
   //  ffout(0)<<"liste points "<<(*__internals).listVerticesIntersectionPerEdges.size()<<"\n";  
   //ffout(0)<<"nb tr intersection objet 0 "<<(*(*__internals).triangleListesIntersectionEdges[0]).size()<<" \n";  
   
2878    size_t nobjother=0;    size_t nobjother=0;
2879    for(size_t no=0;no<nobjectToTreat;++no) {    for(size_t no=0;no<nobjectToTreat;++no) {
2880        nobjother+=(*(*__internals).triangleListes[no]).size();        nobjother+=(*(*__internals).triangleListes[no]).size();
2881    }    }
2882   // ffout(0)<<"nb tr objet1 complet "<<nobjother<<" \n";   // ffout(0)<<"nb tr objet1 complet "<<nobjother<<" \n";
  // ffout(0)<<"nb tr objet1 intersection "<<(*(*__internals).triangleListesIntersectionEdges[1]).size()<<" \n";  
2883    
2884   /* if(scene.nbObjects()>2) {    ffout(0)<<"calculate listes intersection\n";
       ffout(0)<<"nb tr objet2 intersection "<<(*(*__internals).triangleListesIntersectionEdges[2]).size()<<" \n";  
   }*/  
   
   ffout(0)<<"calculate points intersection\n";  
2885        
2886    //create triangleListesIntersectionEdgesPerHexa for all objects    //create triangleListesIntersectionEdgesPerHexa for all objects
2887    (*__internals).__calculateListesIntersection(scene,otherReferences);    (*__internals).__calculateListesIntersection(scene,otherReferences);
2888    
2889    //triangleListesIntersectionPerHexa is for the triangle creating by intersection    size_t nb=nobjectToTreat;
2890    for(std::set<const Cell*>::iterator i=(*__internals).toTreatHexahedra.begin() ;    if(nb==0)
2891        i!=(*__internals).toTreatHexahedra.end() ; ++i)        nb=1;
2892    {    for (size_t numObj = 0; numObj<nb; ++numObj) {
2893        (*__internals).triangleListesIntersectionPerHexa[*i].resize(nobject);      // Data structure convertion
2894        for(size_t no=0 ; no<nobjectToTreat; ++no) {      {
2895            (*__internals).triangleListesIntersectionPerHexa[*i][no]=  #warning keep this ?
2896                new SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList;        size_t n=0;
       }  
   }  
   if(nobjectToTreat!=0) {  
       //calculate points intersection between object 0 and object 1  
       (*__internals).__calculatePointsIntersection(scene,nobjectToTreat,0,1);  
2897    
2898        for(std::set<const Cell*>::iterator i=(*__internals).toTreatHexahedra.begin() ;        ReferenceCounting<Vector<Vertex*> > pVertexList
2899            i!=(*__internals).toTreatHexahedra.end() ; ++i)          = new Vector<Vertex*>((*(*__internals).verticesListes[numObj]).size());
2900        {        (*__internals).listOfVertexList.push_back(pVertexList);
2901            (*__internals).triangleListesIntersectionEdgesTempPerHexa[*i].resize(nobject);        Vector<Vertex*>& vertexList = *pVertexList;
2902            for(size_t no=0 ; no<nobjectToTreat ; ++no) {  
2903                (*__internals).triangleListesIntersectionEdgesTempPerHexa[*i][no]=        for (SurfaceMeshGenerator::Internals::VerticesList::iterator
2904                new SurfaceMeshGenerator::Internals::TriangleVerticesMap;               i = (*(*__internals).verticesListes[numObj]).begin();
2905            }             i != (*(*__internals).verticesListes[numObj]).end(); ++i, ++n) {
2906            for(size_t no=0 ; no<nobjectToTreat ; ++no) {          vertexList[n] = (*i).second;
               for(SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::iterator currentTriangle1  
            = (*(*__internals).triangleListesIntersectionEdgesPerHexa[*i][no]).begin();  
            currentTriangle1 != (*(*__internals).triangleListesIntersectionEdgesPerHexa[*i][no]).end();  
            ++currentTriangle1) {  
                   (*(*__internals).  
      triangleListesIntersectionEdgesTempPerHexa[*i][no])[&(*currentTriangle1)]=&(*currentTriangle1);  
               }  
           }  
2907        }        }
2908        }
2909    
2910        //create mesh intersection for the object 0 and 1      {
2911        (*__internals).__createTrianglesIntersection(0,1,nobjectToTreat);        size_t n=0;
       (*__internals).__createTrianglesIntersection(1,0,nobjectToTreat);  
   }  
     
   if(nobjectToTreat>2) {  
       (*__internals).__createMeshIntersection(nobjectToTreat,scene);  
   }  
   
   (*__internals).plotlistenew(nobject,nobjectToTreat);  
2912    
2913    //ajout des points de verticesListesIntersection qui sont dans le maillage final        ReferenceCounting<Vector<Triangle> > pTriangleList
2914    //ffout(0)<<"nb vertices object 0 "<<(*(*__internals).verticesListes[0]).size()<<"\n";          = new Vector<Triangle> ((*(*__internals).triangleListes[numObj]).size());
2915   // ffout(0)<<"nb vertices object 1 "<<(*(*__internals).verticesListes[1]).size()<<"\n";        (*__internals).listOfTriangleList.push_back(pTriangleList);
2916          Vector<Triangle>& triangleList = *pTriangleList;
2917    for(size_t kk=0 ; kk<nobject ; ++kk) {  
2918          SurfaceMeshGenerator::Internals::MotherCellList::iterator
2919        for(size_t n=0 ; n<(*(*__internals).listTriangleintersection[kk]).size(); ++n) {            itcell=(*(*__internals).cellListObject[numObj]).begin();
2920            TinyVector<2,TinyVector<2,Edge::Pair> > p0=((*(*__internals).listTriangleintersection[kk]))[n][0];        for (SurfaceMeshGenerator::Internals::EdgePairList::iterator
2921            TinyVector<2,TinyVector<2,Edge::Pair> > p1=((*(*__internals).listTriangleintersection[kk]))[n][1];               i = (*(*__internals).triangleListes[numObj]).begin();
2922            TinyVector<2,TinyVector<2,Edge::Pair> > p2=((*(*__internals).listTriangleintersection[kk]))[n][2];             i != (*(*__internals).triangleListes [numObj]).end(); ++i, ++n) {
2923            if(p0[0]==p0[1]) {          SurfaceMeshGenerator::Internals::VerticesList& V = (*(*__internals).verticesListes[numObj]);
2924                if((*(*__internals).verticesListesIntersection[kk]).find(p0[1][0])!=          Triangle T(* V[(*i)[0]],
2925            (*(*__internals).verticesListesIntersection[kk]).end()) {                     * V[(*i)[1]],
2926                    if((*(*__internals).verticesListes[kk]).find(p0[1][0])==(*(*__internals).verticesListes[kk]).end())                     * V[(*i)[2]], 1);
2927                    {  
2928                        (*(*__internals).verticesListes[kk])[p0[1][0]]  #warning Store mother cell here
2929                        =((*(*__internals).verticesListesIntersection[kk])[p0[1][0]]);          T.setMother(*itcell);
2930                    }          triangleList[n] = T;
2931                }          ++itcell;
               else  
                   ffout(0)<<"pbs point 0 pas dans vertices\n";  
           }  
           if(p1[0]==p1[1]) {  
               if((*(*__internals).verticesListesIntersection[kk]).find(p1[1][0])!=  
           (*(*__internals).verticesListesIntersection[kk]).end()) {  
                   if((*(*__internals).verticesListes[kk]).find(p1[1][0])==(*(*__internals).verticesListes[kk]).end())  
                   {  
                       (*(*__internals).verticesListes[kk])[p1[1][0]]=  
                       ((*(*__internals).verticesListesIntersection[kk])[p1[1][0]]);  
                   }  
               }  
               else  
                   ffout(0)<<"pbs point 1 pas dans vertices\n";  
           }  
           if(p2[0]==p2[1]) {  
               if((*(*__internals).verticesListesIntersection[kk]).find(p2[1][0])!=  
           (*(*__internals).verticesListesIntersection[kk]).end()) {  
                   if((*(*__internals).verticesListes[kk]).find(p2[1][0])==(*(*__internals).verticesListes[kk]).end())  
                   {  
                       (*(*__internals).verticesListes[kk])[p2[1][0]]=  
                       ((*(*__internals).verticesListesIntersection[kk])[p2[1][0]]);  
                   }  
               }  
               else  
                   ffout(0)<<"pbs point 0 pas dans vertices\n";  
           }  
2932        }        }
2933    }      }
2934    
2935   // ffout(0)<<"apres nb vertices object 0 "<<(*(*__internals).verticesListes[0]).size()<<"\n";      //(*(*__internals).verticesListes[numObj]).clear();
2936   // ffout(0)<<"apres nb vertices object 1 "<<(*(*__internals).verticesListes[1]).size()<<"\n";      (*(*__internals).triangleListes[numObj]).clear();
2937    
2938  #warning remove this      // Intersection lists
2939    (*__internals).plotliste(nobject);        
2940      
2941        {
2942    #warning keep this ?
2943          size_t n=0;
2944          ReferenceCounting<Vector<Vertex*> > pVertexListIntersection
2945            = new Vector<Vertex*> ((*(*__internals).verticesListesIntersection[numObj]).size());
2946          (*__internals).listOfVertexListIntersection.push_back(pVertexListIntersection);
2947          Vector<Vertex*>& vertexListIntersection = *pVertexListIntersection;
2948    
2949          //for not have double points
2950            for (SurfaceMeshGenerator::Internals::VerticesList::iterator
2951          i = (*(*__internals).verticesListesIntersection[numObj]).begin();
2952          i != (*(*__internals).verticesListesIntersection[numObj]).end(); ++i, ++n) {
2953                if((*(*__internals).verticesListes[numObj]).find((*i).first)
2954            !=(*(*__internals).verticesListes[numObj]).end()) {
2955                    //ffout(0)<<"point existe deja\n";
2956                    SurfaceMeshGenerator::Internals::VerticesList::iterator
2957                    j=(*(*__internals).verticesListes[numObj]).find((*i).first);
2958                    vertexListIntersection[n] = (*j).second;
2959                }
2960                else {
2961                    vertexListIntersection[n] = (*i).second;
2962                    (*(*__internals).verticesListes[numObj])[(*i).first]
2963                                          =(*i).second;            
2964                }
2965            }
2966        }
2967        (*(*__internals).verticesListesIntersection[numObj]).clear();
2968        {
2969          size_t n=0;
2970    
2971    size_t ntr=0;  #warning precompute this:
2972          size_t size = 0;
2973          for (SurfaceMeshGenerator::Internals::TriangleListesIntersectionEdgesPerHexa::iterator
2974                 i = (*__internals).triangleListesIntersectionEdgesPerHexa.begin();
2975               i != (*__internals).triangleListesIntersectionEdgesPerHexa.end(); ++i) {
2976            size += (*(*i).second[numObj]).size();
2977          }
2978    
2979          ReferenceCounting<Vector<Triangle> > pTriangleListIntersection
2980            = new Vector<Triangle>(size);
2981          (*__internals).listOfTriangleListIntersection.push_back(pTriangleListIntersection);
2982          Vector<Triangle>& triangleListIntersection = *pTriangleListIntersection;
2983    
2984          ReferenceCounting<std::map<const Cell*, std::list<Triangle*> > > pHexaToTriangle
2985            = new std::map<const Cell*, std::list<Triangle*> >;
2986          (*__internals).listOfHexaToTriangle.push_back(pHexaToTriangle);
2987          std::map<const Cell*, std::list<Triangle*> >& hexaToTriangle = *pHexaToTriangle;
2988    
2989          for (SurfaceMeshGenerator::Internals::TriangleListesIntersectionEdgesPerHexa::iterator
2990                 i = (*__internals).triangleListesIntersectionEdgesPerHexa.begin();
2991               i != (*__internals).triangleListesIntersectionEdgesPerHexa.end(); ++i) {
2992            for (SurfaceMeshGenerator::Internals::EdgePairList::iterator
2993                   j = (*(*i).second[numObj]).begin();
2994                 j !=  (*(*i).second[numObj]).end(); ++j) {
2995              SurfaceMeshGenerator::Internals::VerticesList& V = (*(*__internals).verticesListes[numObj]);
2996                Triangle T(* V[(*j)[0]],
2997                    * V[(*j)[1]],
2998                    * V[(*j)[2]], 1);
2999                
3000                T.setMother((*i).first);
3001    
3002    for(size_t n = 0; n<nobject; ++n)              triangleListIntersection[n] = T;
     ntr+=(*(*__internals).verticesListes[n]).size();  
3003    
3004    ntr+=(*__internals).listVerticesIntersectionPerEdges.size();              hexaToTriangle[(*i).first].push_back(&triangleListIntersection[n]);
3005                  ++n;
   s_mesh.setNumberOfVertices(ntr);  
  // ffout(0)<<"nb vertices intersection "<<(*__internals).listVerticesIntersectionPerEdges.size()<<"\n";  
   ffout(0)<<"nb tot vertices "<<s_mesh.numberOfVertices()<<"\n";  
   
   //! copies vertices into the s_mesh.  
   int i = 0;  
   for(size_t n = 0; n<scene.nbObjects(); ++n) {  
       for(SurfaceMeshGenerator::Internals::VerticesList::iterator splitVertex  
             =  (*(*__internals).verticesListes[n]).begin();  
           splitVertex != (*(*__internals).verticesListes[n]).end(); ++splitVertex) {  
           s_mesh.vertex(i) = (*(*splitVertex).second);  
           delete (*splitVertex).second;  
           (*splitVertex).second = &s_mesh.vertex(i);  
           i++;  
3006          }          }
3007          }
3008      }      }
     
   for(SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator  
       n=(*__internals).listVerticesIntersectionPerEdges.begin() ;  
       n!=(*__internals).listVerticesIntersectionPerEdges.end() ; ++n) {  
       s_mesh.vertex(i)=(*((*n).second));  
       delete (*n).second;  
       (*n).second = &s_mesh.vertex(i);  
       i++;  
   }  
3009    
3010    i = 0;      (*(*__internals).verticesListes[numObj]).clear();
3011    size_t ncell=0;    }
3012      (*__internals).triangleListesIntersectionEdgesPerHexa.clear();
3013    
3014    for(size_t n = 0; n<scene.nbObjects(); ++n)    (*__internals).plotliste();
     ncell+=(*(*__internals).triangleListes[n]).size();  
3015    
3016    size_t ntr1=(*(*__internals).triangleListes[0]).size();    size_t n0=0,n1=0;
3017    for(size_t no=0 ; no<nobject ; ++no) {    if(nobjectToTreat!=0) {
3018        ncell+=(*(*__internals).listTriangleintersection[no]).size();        for(size_t ncompt=0 ; ncompt<=nobjectToTreat-2; ++ncompt) {
3019              ffout(0)<<"treat object "<<ncompt+1<<"\n";
3020              //calculate points intersection between object 0 and object 1
3021              (*__internals).__calculatePointsIntersection(scene,0,1);
3022              //  ffout(0)<<"end calculatePointsIntersection\n";
3023       //create mesh intersection for the object 0 and 1
3024              (*__internals).__createTrianglesIntersection(0,1);
3025              ffout(0)<<"nb triangles in intersection 0 "<<(*__internals).triangleListIntersectionNew.size()<<"\n";
3026              (*__internals).__createTrianglesIntersection(1,0);
3027              ffout(0)<<"nb triangles in intersection 1 "<<(*__internals).triangleListIntersectionNew.size()<<"\n";
3028    
3029      //erase triangles from object 0 and 1 and add triangle from intersection of objects 0-1
3030              SurfaceMeshGenerator::Internals::ListOfTriangleList::iterator
3031                  it=(*__internals).listOfTriangleListIntersection.begin();
3032              SurfaceMeshGenerator::Internals::ListOfTriangleList::iterator
3033                  er=(*__internals).listOfTriangleListIntersection.erase(it);
3034              it=(*__internals).listOfTriangleListIntersection.begin();
3035              er=(*__internals).listOfTriangleListIntersection.erase(it);
3036    
3037              size_t sizeinter=(*__internals).triangleListIntersectionNew.size();
3038              ReferenceCounting<Vector<Triangle> > pTriangleListIntersection = new Vector<Triangle>(sizeinter);
3039              (*__internals).listOfTriangleListIntersection.push_front(pTriangleListIntersection);
3040              Vector<Triangle>& triangleListIntersection = *pTriangleListIntersection;
3041    
3042              SurfaceMeshGenerator::Internals::ListOfHexaToTriangle::iterator
3043                  ithexa=(*__internals).listOfHexaToTriangle.begin();
3044              SurfaceMeshGenerator::Internals::ListOfHexaToTriangle::iterator
3045                  erhexa=(*__internals).listOfHexaToTriangle.erase(ithexa);
3046              ithexa=(*__internals).listOfHexaToTriangle.begin();
3047              erhexa=(*__internals).listOfHexaToTriangle.erase(ithexa);
3048              (*__internals).listOfHexaToTriangle.push_front(new std::map<const Cell*, std::list<Triangle*> >);
3049              ithexa=(*__internals).listOfHexaToTriangle.begin();
3050    
3051              for(size_t size=0 ; size<sizeinter ; ++size) {
3052                  triangleListIntersection[size]=(*__internals).triangleListIntersectionNew[size];
3053                  (*(*ithexa))[&triangleListIntersection[size].mother()].push_back(&triangleListIntersection[size]);
3054              }
3055              (*__internals).triangleListIntersectionNew.clear();
3056    
3057              //copy meshes not intersected of object 0 and 1
3058              ReferenceCounting<Vector<Triangle> > triangleListTemp;
3059              size_t size=0;
3060              SurfaceMeshGenerator::Internals::ListOfTriangleList::iterator
3061                  triangleit=(*__internals).listOfTriangleList.begin();
3062              size+=(*(*triangleit)).size();
3063              n0=size;
3064              ++triangleit;
3065              size+=(*(*triangleit)).size();
3066              n1=size;
3067              triangleListTemp=new Vector<Triangle>(size);
3068              triangleit=(*__internals).listOfTriangleList.begin();
3069              for(size_t s=0 ; s<(*(*triangleit)).size() ; ++s) {
3070                  (*triangleListTemp)[s]=(*(*triangleit))[s];
3071              }
3072              size_t numsize=(*(*triangleit)).size();
3073              ++triangleit;
3074              for(size_t s=0 ; s<(*(*triangleit)).size() ; ++s) {
3075                  (*triangleListTemp)[numsize+s]=(*(*triangleit))[s];
3076              }
3077    
3078              it=(*__internals).listOfTriangleList.begin();
3079              er=(*__internals).listOfTriangleList.erase(it);
3080              it=(*__internals).listOfTriangleList.begin();
3081              er=(*__internals).listOfTriangleList.erase(it);
3082              (*__internals).listOfTriangleList.push_front(new Vector<Triangle>(size));
3083              it=(*__internals).listOfTriangleList.begin();
3084              for(size_t s=0 ; s<(*triangleListTemp).size() ; ++s) {
3085                  (*(*it))[s]=(*triangleListTemp)[s];
3086              }
3087          }
3088    }    }
3089    size_t ninter=(*(*__internals).listTriangleintersection[0]).size();    //end intersection mesh
   size_t ninter2=(*(*__internals).listTriangleintersection[1]).size();  
3090    
3091    s_mesh.setNumberOfCells(ncell);    //construction of list vertices
3092    ffout(0)<<"nb tot cells "<<s_mesh.numberOfCells()<<" \n";    
3093     // SurfaceMeshGenerator::Internals::ListOfTriangleList::const_iterator
3094       //   itcell=((*__internals).listOfTriangleList).begin();
3095      SurfaceMeshGenerator::Internals::ListOfTriangleList::const_iterator
3096          intertr=((*__internals).listOfTriangleListIntersection).begin();
3097    /*  for(size_t size=0 ; size< (*(*itcell)).size() ; ++size) {
3098          const Triangle& T=(*(*itcell))[size];
3099          if((*__internals).listOfVertexMesh.find(&T(0))==(*__internals).listOfVertexMesh.end()) {
3100              (*__internals).listOfVertexMesh[&T(0)]=(&T(0));
3101          }
3102          if((*__internals).listOfVertexMesh.find(&T(1))==(*__internals).listOfVertexMesh.end()) {
3103              (*__internals).listOfVertexMesh[&T(1)]=(&T(1));
3104          }
3105          if((*__internals).listOfVertexMesh.find(&T(2))==(*__internals).listOfVertexMesh.end()) {
3106              (*__internals).listOfVertexMesh[&T(2)]=(&T(2));
3107          }
3108      }*/
3109      for(size_t size=0 ; size< (*(*intertr)).size() ; ++size) {
3110          const Triangle& T=(*(*intertr))[size];
3111          if((*__internals).listOfVertexMesh.find(&T(0))==(*__internals).listOfVertexMesh.end()) {
3112              (*__internals).listOfVertexMesh[&T(0)]=(&T(0));
3113          }
3114          if((*__internals).listOfVertexMesh.find(&T(1))==(*__internals).listOfVertexMesh.end()) {
3115              (*__internals).listOfVertexMesh[&T(1)]=(&T(1));
3116          }
3117          if((*__internals).listOfVertexMesh.find(&T(2))==(*__internals).listOfVertexMesh.end()) {
3118              (*__internals).listOfVertexMesh[&T(2)]=(&T(2));
3119          }
3120      }
3121     /* Vector<const Vertex*> listVertexFinal((*__internals).listOfVertexMesh.size());
3122      Vector<Vertex> listVertexFinal2((*__internals).listOfVertexMesh.size());
3123      size_t number=0;
3124      for(std::map<const Vertex*,const Vertex*>::iterator it=(*__internals).listOfVertexMesh.begin() ;
3125          it!=(*__internals).listOfVertexMesh.end() ; it++) {
3126          listVertexFinal[number]=(*it).second;
3127          listVertexFinal2[number]=*((*it).second);
3128          ++number;
3129      }*/
3130       size_t ntr=0;
3131    
3132    SurfaceMeshGenerator::Internals::MotherCellList::iterator motherCell     //ntr+=listVertexFinal.size();
3133      = (*__internals).cellList.begin();    ntr+=(*__internals).listOfVertexMesh.size();
3134      
3135       s_mesh.setNumberOfVertices(ntr);
3136    for(size_t n = 0; n<scene.nbObjects(); ++n) {     ffout(0)<<"nb tot vertices "<<s_mesh.numberOfVertices()<<"\n";
     for(SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle  
           = (*(*__internals).triangleListes[n]).begin();  
         currentTriangle != (*(*__internals).triangleListes[n]).end();  
         ++currentTriangle) {  
       Vertex& V0 = *((*(*__internals).verticesListes[n])[(*currentTriangle)[0]]);  
       Vertex& V1 = *((*(*__internals).verticesListes[n])[(*currentTriangle)[1]]);  
       Vertex& V2 = *((*(*__internals).verticesListes[n])[(*currentTriangle)[2]]);  
3137    
3138       //! copies vertices into the s_mesh.
3139       int i = 0;
3140      /* for(size_t v=0 ; v<ntr ; ++v) {
3141           s_mesh.vertex(i) = *listVertexFinal[v];
3142           listVertexFinal[v] = &s_mesh.vertex(i);
3143           i++;
3144         }*/
3145       for(std::map<const Vertex*,const Vertex*>::iterator it=(*__internals).listOfVertexMesh.begin() ;
3146           it!=(*__internals).listOfVertexMesh.end() ; it++) {
3147           s_mesh.vertex(i) = *((*it).second);
3148           delete (*it).second;
3149           (*it).second = &s_mesh.vertex(i);
3150           i++;
3151       }
3152    
3153       i = 0;
3154       size_t ncell=0;
3155     //  ncell+=(*(*itcell)).size();
3156     //  ffout(0)<<"ncell "<<ncell<<" n0 "<<n0<<" n1 "<<n1<<"\n";
3157        ncell+=(*(*intertr)).size();
3158    
3159        s_mesh.setNumberOfCells(ncell);
3160        ffout(0)<<"nb tot cells "<<s_mesh.numberOfCells()<<" \n";
3161    
3162    //   SurfaceMeshGenerator::Internals::MotherCellList::iterator motherCell
3163    //     = (*__internals).cellList.begin();
3164    
3165       /* for(size_t size=0 ; size< (*(*itcell)).size() ; ++size) {
3166            const Triangle& T=(*(*itcell))[size];
3167            Vertex P=T(0);
3168            size_t n0=listVertexFinal2.number(P);
3169            size_t n1=listVertexFinal2.number(T(1));
3170            size_t n2=listVertexFinal2.number(T(2));
3171            const Vertex& V0=*listVertexFinal[n0];
3172            const Vertex& V1=*listVertexFinal[n1];
3173            const Vertex& V2=*listVertexFinal[n2];
3174            s_mesh.cell(i) = Triangle(V0,V1,V2);
3175    #warning SHOULD USE BC NUMBER (ie OBJECT NUMBER IN SCENE INSTEAD)
3176            s_mesh.cell(i).reference() = 1;
3177            s_mesh.cell(i).setMother(&T.mother());
3178            i++;
3179            //motherCell++;
3180        }   */
3181    /*    for(size_t size=0 ; size< (*(*itcell)).size() ; ++size) {
3182            const Triangle& T=(*(*itcell))[size];
3183            std::map<const Vertex*,const Vertex*>::iterator
3184            it0=(*__internals).listOfVertexMesh.find(&T(0));
3185            const Vertex& V0=*((*it0).second);
3186            std::map<const Vertex*,const Vertex*>::iterator
3187                it1=(*__internals).listOfVertexMesh.find(&T(1));
3188            const Vertex& V1=*((*it1).second);
3189            std::map<const Vertex*,const Vertex*>::iterator
3190                it2=(*__internals).listOfVertexMesh.find(&T(2));
3191            const Vertex& V2=*((*it2).second);
3192            s_mesh.cell(i) = Triangle(V0,V1,V2);
3193    #warning SHOULD USE BC NUMBER (ie OBJECT NUMBER IN SCENE INSTEAD)
3194            s_mesh.cell(i).reference() = 1;
3195            s_mesh.cell(i).setMother(&T.mother());
3196            i++;
3197            //motherCell++;
3198        }*/
3199        for(size_t size=0 ; size< (*(*intertr)).size() ; ++size) {
3200            const Triangle& T=(*(*intertr))[size];
3201            std::map<const Vertex*,const Vertex*>::iterator
3202            it0=(*__internals).listOfVertexMesh.find(&T(0));
3203            const Vertex& V0=*((*it0).second);
3204            std::map<const Vertex*,const Vertex*>::iterator
3205            it1=(*__internals).listOfVertexMesh.find(&T(1));
3206            const Vertex& V1=*((*it1).second);
3207            std::map<const Vertex*,const Vertex*>::iterator
3208            it2=(*__internals).listOfVertexMesh.find(&T(2));
3209            const Vertex& V2=*((*it2).second);
3210            s_mesh.cell(i) = Triangle(V0,V1,V2);
3211                    
       s_mesh.cell(i) = Triangle(V0,V1,V2);  
3212  #warning SHOULD USE BC NUMBER (ie OBJECT NUMBER IN SCENE INSTEAD)  #warning SHOULD USE BC NUMBER (ie OBJECT NUMBER IN SCENE INSTEAD)
3213        s_mesh.cell(i).reference() = 1;          s_mesh.cell(i).reference() = 1;
3214        s_mesh.cell(i).setMother(*motherCell);          s_mesh.cell(i).setMother(&T.mother());
3215        i++;          i++;
3216        motherCell++;          //motherCell++;
3217      }      }
3218    }      
3219        
3220    for(size_t kk=0 ; kk<nobject ; ++kk) {  //   //  plot2(s_mesh, (*__internals).toTreatHexahedra, SMesh);
3221        for(size_t n=0 ; n<(*(*__internals).listTriangleintersection[kk]).size(); ++n) {  //   //ffout(0) << "plot medit \n";
3222       plotmedit(nobject,s_mesh, (*__internals).toTreatHexahedra,cuttedHexahedra,
3223            TinyVector<2,TinyVector<2,Edge::Pair> > p0=((*(*__internals).listTriangleintersection[kk]))[n][0];              nobjother,n0,n1);
           TinyVector<2,TinyVector<2,Edge::Pair> > p1=((*(*__internals).listTriangleintersection[kk]))[n][1];  
           TinyVector<2,TinyVector<2,Edge::Pair> > p2=((*(*__internals).listTriangleintersection[kk]))[n][2];  
           Vertex *V0,*V1,*V2;  
           std::pair<TinyVector<2,TinyVector<2,Edge::Pair> >,TinyVector<2,size_t> > pp0,pp1,pp2;  
           TinyVector<2,size_t> nuum;  
           nuum[0]=0;  
           nuum[1]=0;  
           pp0.second=nuum;  
           pp1.second=nuum;  
           pp2.second=nuum;  
           pp0.first=p0;  
           pp1.first=p1;  
           pp2.first=p2;  
   
           if(p0[0][0]==p0[1][0]) {  
               V0=((*(*__internals).verticesListes[kk])[p0[0][0]]);  
           }  
           else {  
               SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator  
               ip0=(*__internals).listVerticesIntersectionPerEdges.find(pp0);  
               if(ip0!=(*__internals).listVerticesIntersectionPerEdges.end())  
                   V0 =(*ip0).second;  
           }  
           if(p1[0][0]==p1[1][0]) {  
               V1=((*(*__internals).verticesListes[kk])[p1[0][0]]);  
           }  
           else {  
               SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator  
               ip0=(*__internals).listVerticesIntersectionPerEdges.find(pp1);  
               if(ip0!=(*__internals).listVerticesIntersectionPerEdges.end())  
                   V1 =(*ip0).second;  
           }  
           if(p2[0][0]==p2[1][0]) {  
               V2=((*(*__internals).verticesListes[kk])[p2[0][0]]);  
           }  
           else {  
               SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator  
               ip0=(*__internals).listVerticesIntersectionPerEdges.find(pp2);  
               if(ip0!=(*__internals).listVerticesIntersectionPerEdges.end())  
                   V2 =(*ip0).second;  
           }  
           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(*motherCell);  
           i++;  
           motherCell++;  
       }  
   }  
   
   //  plot2(s_mesh, (*__internals).toTreatHexahedra, SMesh);  
   //ffout(0) << "plot medit \n";  
   plotmedit(nobject,s_mesh, (*__internals).toTreatHexahedra,cuttedHexahedra,  
             (*__internals).verticesListesIntersection,(*__internals).triangleListesIntersectionEdges,  
             nobjother,ninter,ninter2,ntr1);  
3224  }  }
3225    
3226  SurfaceMeshGenerator::SurfaceMeshGenerator()  SurfaceMeshGenerator::SurfaceMeshGenerator()
# Line 3432  void plotmedit(const size_t& nobject, Line 3279  void plotmedit(const size_t& nobject,
3279                 const SurfaceMeshOfTriangles& s_mesh,                 const SurfaceMeshOfTriangles& s_mesh,
3280                 std::set<const Cell*>& toTreatHexahedra,                 std::set<const Cell*>& toTreatHexahedra,
3281                 std::set<const Cell*>& cuttedHexahedra,                 std::set<const Cell*>& cuttedHexahedra,
3282                 std::vector<ReferenceCounting<std::map<Edge::Pair, Vertex*> > > verticesListesIntersection,                 size_t nobjother,size_t n0,size_t n1) {
                std::vector<ReferenceCounting<std::list<TinyVector<3,TinyVector<2,TinyVector<2,Edge::Pair> > > > > >  
                triangleListesIntersectionEdges,  
                size_t nobjother,size_t ninter,size_t ninter2,size_t ntr1) {  
3283    
3284    // **********************************************************************    // **********************************************************************
3285  std::ofstream sortie3("cutted.mesh");      std::ofstream sortie3("cutted.mesh");
3286    
3287      sortie3<<"MeshVersionFormatted 1"<<" \n";      sortie3<<"MeshVersionFormatted 1"<<" \n";
3288      sortie3<<" \n";      sortie3<<" \n";
# Line 3474  std::ofstream sortie3("cutted.mesh"); Line 3318  std::ofstream sortie3("cutted.mesh");
3318    sortie<<"3"<<" \n";    sortie<<"3"<<" \n";
3319    sortie<<" \n";    sortie<<" \n";
3320    sortie<<"Vertices"<<" \n";    sortie<<"Vertices"<<" \n";
3321    sortie<<8*toTreatHexahedra.size()+s_mesh.numberOfVertices()<<" \n";    sortie<<8*toTreatHexahedra.size()/*+3*s_mesh.numberOfCells()*/+s_mesh.numberOfVertices()<<" \n";
3322    for (std::set<const Cell*>::iterator i = toTreatHexahedra.begin();    for (std::set<const Cell*>::iterator i = toTreatHexahedra.begin();
3323         i != toTreatHexahedra.end(); ++i)         i != toTreatHexahedra.end(); ++i)
3324      {      {
# Line 3489  std::ofstream sortie3("cutted.mesh"); Line 3333  std::ofstream sortie3("cutted.mesh");
3333        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";
3334      }      }
3335    
3336    /*        for (size_t i=0; i<s_mesh.numberOfCells(); ++i)       /*  for (size_t i=0; i<s_mesh.numberOfCells(); ++i)
3337              {              {
3338              const Triangle& T = s_mesh.cell(i);              const Triangle& T = s_mesh.cell(i);
3339              sortie<<T(0)[0]<<" "<<T(0)[1]<<" "<<T(0)[2]<<" 0"<<" \n";              sortie<<T(0)[0]<<" "<<T(0)[1]<<" "<<T(0)[2]<<" 0"<<" \n";
# Line 3500  std::ofstream sortie3("cutted.mesh"); Line 3344  std::ofstream sortie3("cutted.mesh");
3344    sortie<<" \n";    sortie<<" \n";
3345    sortie<<"Triangles"<<"\n";    sortie<<"Triangles"<<"\n";
3346    sortie<<s_mesh.numberOfCells()<<"\n";    sortie<<s_mesh.numberOfCells()<<"\n";
3347   // ffout(0)<<"fin ecriture vertex\n";    ffout(0)<<"fin ecriture vertex\n";
3348    for (size_t i=0; i<s_mesh.numberOfCells(); ++i)    for (size_t i=0; i<s_mesh.numberOfCells(); ++i)
3349      {      {
3350        const Triangle& T = s_mesh.cell(i);        const Triangle& T = s_mesh.cell(i);
3351       // ffout(0)<<s_mesh.vertexNumber(T(0))<<" "<<s_mesh.vertexNumber(T(1))<<" "<<s_mesh.vertexNumber(T(2))<<"\n";       // ffout(0)<<s_mesh.vertexNumber(T(0))<<" "<<s_mesh.vertexNumber(T(1))<<" "<<s_mesh.//vertexNumber(T(2))<<"\n";
3352        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;
3353                
3354  //                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";
3355        if(i<ntr1)        if(i<n0)
3356            sortie<<" 1"<<"\n";            sortie<<" 1"<<"\n";
3357        else if(i<(nobjother))        else if(i<n1)
3358            sortie<<" 2"<<"\n";            sortie<<" 2"<<"\n";
3359        else if(i<(nobjother+ninter))       /* else if(i<(nobjother+ninter))
3360            sortie<<" 3"<<"\n";            sortie<<" 3"<<"\n";
3361        else if(i<(nobjother+ninter+ninter2))        else if(i<(nobjother+ninter+ninter2))
3362            sortie<<" 4"<<"\n";            sortie<<" 4"<<"\n";*/
3363        else        else
3364            sortie<<" 5"<<"\n";            sortie<<" 4"<<"\n";
3365      }      }
3366   // ffout(0)<<"fin ecriture triangles\n";    ffout(0)<<"fin ecriture triangles\n";
3367    sortie<<" \n";    sortie<<" \n";
3368    sortie<<"Hexahedra"<<"\n";    sortie<<"Hexahedra"<<"\n";
3369    sortie<<toTreatHexahedra.size()<<" \n";    sortie<<toTreatHexahedra.size()<<" \n";
# Line 3527  std::ofstream sortie3("cutted.mesh"); Line 3371  std::ofstream sortie3("cutted.mesh");
3371    for(size_t i=0;i<toTreatHexahedra.size();i++)    for(size_t i=0;i<toTreatHexahedra.size();i++)
3372      sortie<<8*i+1<<" "<<8*i+2<<" "<<8*i+3<<" "<<8*i+4<<" "<<8*i+5<<" "<<8*i+6<<" "<<8*i+7<<" "<<8*i+8<<" 3"<<"\n";      sortie<<8*i+1<<" "<<8*i+2<<" "<<8*i+3<<" "<<8*i+4<<" "<<8*i+5<<" "<<8*i+6<<" "<<8*i+7<<" "<<8*i+8<<" 3"<<"\n";
3373    sortie<<"End"<<" \n";    sortie<<"End"<<" \n";
3374   // ffout(0)<<"fin ecriture intersection.mesh\n";    ffout(0)<<"fin ecriture intersection.mesh\n";
3375    
3376  }  }
3377    

Legend:
Removed from v.1.12  
changed lines
  Added in v.1.13

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