/[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.15 by dobj, Thu Jul 10 15:37:31 2003 UTC revision 1.16 by delpinux, Sun Jul 13 17:01:27 2003 UTC
# Line 58  class SurfaceMeshGenerator::Internals Line 58  class SurfaceMeshGenerator::Internals
58  private:  private:
59    class IntersectionPoints    class IntersectionPoints
60    //class SurfaceMeshGenerator::Internals::IntersectionPoints    //class SurfaceMeshGenerator::Internals::IntersectionPoints
61  {    {
62  public :    public :
63      TinyVector<2,size_t> numpair;      TinyVector<2,size_t> numpair;
64      std::vector<const  Vertex* > listofPoints;      std::vector<const  Vertex* > listofPoints;
65      TinyVector< 3,size_t > cutTriangle1;      TinyVector< 3,size_t > cutTriangle1;
66      TinyVector< 3,size_t > cutTriangle2;      TinyVector< 3,size_t > cutTriangle2;
67    
68      inline const size_t size() const {      inline const size_t size() const {
69          return listofPoints.size();        return listofPoints.size();
70      }      }
71    
72      void addPoint(const Vertex*&  V) {      void addPoint(const Vertex*&  V) {
73          listofPoints.push_back( V );        listofPoints.push_back( V );
74      }      }
75        
76    
77    //constructor      //constructor
78      IntersectionPoints() {      IntersectionPoints() {
79          numpair=0;        numpair=0;
80          cutTriangle1=0;        cutTriangle1=0;
81          cutTriangle2=0;        cutTriangle2=0;
82      }      }
83    
84    IntersectionPoints(const IntersectionPoints& P) {      IntersectionPoints(const IntersectionPoints& P) {
85        numpair=P.numpair;        numpair=P.numpair;
86        listofPoints=P.listofPoints;        listofPoints=P.listofPoints;
87        cutTriangle1=P.cutTriangle1;        cutTriangle1=P.cutTriangle1;
88        cutTriangle2=P.cutTriangle2;        cutTriangle2=P.cutTriangle2;
89    }      }
   
   //! Destructor.  
   ~IntersectionPoints()  
   {  
     ;  
   }  
       
 };  
90    
91        //! Destructor.
92        ~IntersectionPoints()
93        {
94          ;
95        }
96      };
97    
98    class TriangleCut;    class TriangleCut;
99    friend class SurfaceMeshGenerator;    friend class SurfaceMeshGenerator;
# Line 137  public : Line 135  public :
135    std::vector<ReferenceCounting<EdgePairList> > triangleListes;    std::vector<ReferenceCounting<EdgePairList> > triangleListes;
136    std::vector<ReferenceCounting<EdgePairList> > localTriangleListesEdges;    std::vector<ReferenceCounting<EdgePairList> > localTriangleListesEdges;
137    EdgePairList localTriangleList;    EdgePairList localTriangleList;
   
138    
139    //! This list is used to provide an orientation to the Surfacic Elements.    //! This list is used to provide an orientation to the Surfacic Elements.
140    std::list<TinyVector<3> > InsideVertexList;    std::list<TinyVector<3> > InsideVertexList;
# Line 164  public : Line 161  public :
161    
162    typedef std::pair<const Triangle*,const Triangle*> PairTriangleNew;    typedef std::pair<const Triangle*,const Triangle*> PairTriangleNew;
163    //std::vector<Vertex*> listVertexIntersectionNew;    //std::vector<Vertex*> listVertexIntersectionNew;
164    typedef std::map< std::pair<const Vertex*,const Vertex*> ,    typedef std::map<std::pair<const Vertex*,const Vertex*> ,
165          TinyVector<2,const Vertex* > > PIntersection;                     TinyVector<2,const Vertex* > > PIntersection;
166    PIntersection listVertexIntersectionNew;    PIntersection listVertexIntersectionNew;
167    
168    typedef std::map< PairTriangleNew , IntersectionPoints> PairTriangleToIntersectionPoints;    typedef std::map< PairTriangleNew , IntersectionPoints> PairTriangleToIntersectionPoints;
# Line 237  public : Line 234  public :
234        
235    const Vertex& __determinateNumber(const int& p,    const Vertex& __determinateNumber(const int& p,
236                                      const TinyVector<3,TinyVector<2,const Vertex*> >& pointsint);                                      const TinyVector<3,TinyVector<2,const Vertex*> >& pointsint);
237   void __determinateCase(const int& ncase,    void __determinateCase(const int& ncase,
238                             const TriangleCut & K,                           const TriangleCut & K,
239                             TinyVector<3,int>& tt,                           TinyVector<3,int>& tt,
240                             TinyVector<3,int>& pp);                           TinyVector<3,int>& pp);
241    
242    void __addTriangle(const Vertex& V0,    void __addTriangle(const Vertex& V0,
243                       const Vertex& V1,                       const Vertex& V1,
244                       const Vertex& V2,                       const Vertex& V2,
245                       const Cell* cell);                       const Cell* cell);
246    
247    void __createCase1(int ncase,    void __createCase1(int ncase,
248                       const TriangleCut& K,                       const TriangleCut& K,
249                       const Cell* cell);                       const Cell* cell);
# Line 259  public : Line 258  public :
258                                       const Cell* cell);                                       const Cell* cell);
259        
260    TinyVector<3,TinyVector<2,int> > __determinateEdgeCut(const size_t& numobject,    TinyVector<3,TinyVector<2,int> > __determinateEdgeCut(const size_t& numobject,
261                                           const size_t & objectNumber,                                                          const size_t & objectNumber,
262                                           const Triangle*& currentTriangle,                                                          const Triangle*& currentTriangle,
263                                           const Cell* cell);                                                          const Cell* cell);
264    
265    TinyVector<3,size_t> __putEdgeRef(const size_t& ncase,    TinyVector<3,size_t> __putEdgeRef(const size_t& ncase,
266                                      const Triangle& currentTriangle,                                      const Triangle& currentTriangle,
# Line 271  public : Line 270  public :
270                              const Cell* cell);                              const Cell* cell);
271    
272    void __createByEdges(const TriangleCut& K,const Cell* cell) ;    void __createByEdges(const TriangleCut& K,const Cell* cell) ;
273    
274    void  __createNotIn(const TriangleCut& K,const Cell* cell);    void  __createNotIn(const TriangleCut& K,const Cell* cell);
275    
276    void  __createNotIn2(const TriangleCut& K,const Cell* cell);    void  __createNotIn2(const TriangleCut& K,const Cell* cell);
277    
278    void __createTrianglesIntersection(const size_t& numobject,    void __createTrianglesIntersection(const size_t& numobject,
279                                       const size_t& objectNumber);                                       const size_t& objectNumber);
280    
# Line 282  public : Line 284  public :
284                                 std::vector<const Edge*>& cutEdges,                                 std::vector<const Edge*>& cutEdges,
285                                 const Reference::Representation& ref,                                 const Reference::Representation& ref,
286                                 const Cell& currentCell);                                 const Cell& currentCell);
287    
288    void __createSurface(Structured3DMesh& SMesh,    void __createSurface(Structured3DMesh& SMesh,
289                         const Domain& Omega,                         const Domain& Omega,
290                         const Shape& S ,                         const Shape& S ,
# Line 291  public : Line 294  public :
294                         std::set<const Cell*>& meshedHexahedra);                         std::set<const Cell*>& meshedHexahedra);
295    
296    bool __isDegenerate(const Triangle& T);    bool __isDegenerate(const Triangle& T);
297    
298    void __casePinterNotInList(const TriangleCut& K,const Cell* cell);    void __casePinterNotInList(const TriangleCut& K,const Cell* cell);
299    
300    TinyVector<3,const Vertex*> __determinatePinterNotInList(const TriangleCut& K, TinyVector<3,size_t>);    TinyVector<3,const Vertex*> __determinatePinterNotInList(const TriangleCut& K, TinyVector<3,size_t>);
301  /// @todo FUNCTION TO REMOVE    /// @todo FUNCTION TO REMOVE
302    static void __getIntersectionReferences(const Intersection& I,    static void __getIntersectionReferences(const Intersection& I,
303                                            std::map<TinyVector<3>, Object*, Internals::__ltv>& otherReferences);                                            std::map<TinyVector<3>, Object*, Internals::__ltv>& otherReferences);
304    void plotliste();    void plotliste();
305    
306    void plotini(size_t nobjectToTreat);    void plotini(size_t nobjectToTreat);
307    
308  public:  public:
# Line 309  public: Line 315  public:
315    
316    
317  class SurfaceMeshGenerator::Internals::TriangleCut  class SurfaceMeshGenerator::Internals::TriangleCut
318  : public Triangle    : public Triangle
319  {  {
320  public :  public :
321      TinyVector<3,TinyVector<2,int> > edgecut;    TinyVector<3,TinyVector<2,int> > edgecut;
322      TinyVector<3,size_t> isInTriangle;    TinyVector<3,size_t> isInTriangle;
323      TinyVector<3,TinyVector<2,const Vertex*> > pointsIntersection;    TinyVector<3,TinyVector<2,const Vertex*> > pointsIntersection;
324      std::map< const Vertex* ,const Vertex*> pointsin;    std::map< const Vertex* ,const Vertex*> pointsin;
325      size_t numobject;    size_t numobject;
326      size_t otherobject;    size_t otherobject;
327    
328      //constructor    //constructor
329      TriangleCut() {    TriangleCut()
330      }    {
331      TriangleCut(const TriangleCut& K) {      ;
332          edgecut=K.edgecut;    }
         isInTriangle=K.isInTriangle;  
         pointsIntersection=K.pointsIntersection;  
         pointsin=K.pointsin;  
         numobject=K.numobject;  
         otherobject=K.otherobject;  
     }  
     TriangleCut(const Vertex& v0,  
                 const Vertex& v1,  
                 const Vertex& v2,  
                 size_t num,  
                 size_t num2)  
     {  
         __vertices[0] = (Vertex*)&v0;  
         __vertices[1] = (Vertex*)&v1;  
         __vertices[2] = (Vertex*)&v2;  
         numobject=num;  
         otherobject=num2;  
         edgecut=0;  
         isInTriangle=0;  
     }  
333    
334      //! Destructor.    TriangleCut(const TriangleCut& K) {
335      ~TriangleCut()      edgecut=K.edgecut;
336      {      isInTriangle=K.isInTriangle;
337          ;      pointsIntersection=K.pointsIntersection;
338      }      pointsin=K.pointsin;
339        numobject=K.numobject;
340        otherobject=K.otherobject;
341      }
342    
343      TriangleCut(const Vertex& v0,
344                  const Vertex& v1,
345                  const Vertex& v2,
346                  size_t num,
347                  size_t num2)
348      {
349        __vertices[0] = (Vertex*)&v0;
350        __vertices[1] = (Vertex*)&v1;
351        __vertices[2] = (Vertex*)&v2;
352        numobject=num;
353        otherobject=num2;
354        edgecut=0;
355        isInTriangle=0;
356      }
357    
358      //! Destructor.
359      ~TriangleCut()
360      {
361        ;
362      }
363  };  };
364    
365  void  void
# Line 361  __setBoundaryType(const Boundary& B, Line 371  __setBoundaryType(const Boundary& B,
371    switch (B.type()) {    switch (B.type()) {
372    case (Boundary::povRay): {    case (Boundary::povRay): {
373      const BoundaryPOVRay& bp = dynamic_cast<const BoundaryPOVRay&>(B);      const BoundaryPOVRay& bp = dynamic_cast<const BoundaryPOVRay&>(B);
374      for (size_t l=0; l<bp.nbPOVRef(); ++l)      for (size_t l=0; l<bp.nbPOVRef(); ++l) {
375        vrRefList[bp.POVRef(l)]++;        vrRefList[bp.POVRef(l)]++;
376        }
377      break;      break;
378    }    }
379    case (Boundary::structured3DMesh): {    case (Boundary::structured3DMesh): {
380      const BoundaryStructured3DMesh& bs      const BoundaryStructured3DMesh& bs
381        = dynamic_cast<const BoundaryStructured3DMesh&>(B);        = dynamic_cast<const BoundaryStructured3DMesh&>(B);
382      for (size_t l=0; l<6; ++l)      for (size_t l=0; l<6; ++l) {
383        faces[l] = faces[l] || bs.structMesh(l);        faces[l] = faces[l] || bs.structMesh(l);
384        }
385      break;      break;
386    }    }
387    case (Boundary::surfaceMesh): {    case (Boundary::surfaceMesh): {
# Line 415  __splitEdge(const Edge& e, Line 427  __splitEdge(const Edge& e,
427        vs = (vo+vi)/2;        vs = (vo+vi)/2;
428      }      }
429    }    }
430    
431    if(!S.inside(vs))    if(!S.inside(vs))
432        return vs;      return vs;
433    else    else
434        return vo;      return vo;
   
  // return vs;  
435  }  }
436    
437  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
# Line 428  __createListIntersection(const Shape& S, Line 439  __createListIntersection(const Shape& S,
439                           std::vector<const Edge*>& cutEdges,                           std::vector<const Edge*>& cutEdges,
440                           const Reference::Representation& ref,                           const Reference::Representation& ref,
441                           const size_t& objectNumber,                           const size_t& objectNumber,
442                           const CartesianHexahedron& currentCell) {                           const CartesianHexahedron& currentCell)
443    {
444      switch (cutEdges.size()) {    switch (cutEdges.size()) {
445          case 3: {    case 3: {
446              std::vector<TinyVector<3> > V(cutEdges.size());      std::vector<TinyVector<3> > V(cutEdges.size());
447              for(size_t i=0; i<cutEdges.size(); ++i) {      for(size_t i=0; i<cutEdges.size(); ++i) {
448                  const Edge& e = *cutEdges[i];        const Edge& e = *cutEdges[i];
449                  if(e(0).Ref().isSet(ref))        if(e(0).Ref().isSet(ref)) {
450                      V[i] = (e(0).Vector3() - e(1).Vector3());          V[i] = (e(0).Vector3() - e(1).Vector3());
451                  else                    } else {
452                      V[i] = (e(1).Vector3() - e(0).Vector3());          V[i] = (e(1).Vector3() - e(0).Vector3());
453              }        }
454        }
             if ((V[0]^V[1])*V[2] < 0)  
             std::swap(cutEdges[1], cutEdges[2]);  
   
             //! Addes a triangle to the list.  
             (*localTriangleListesEdges[objectNumber]).push_front( TinyVector<3,Edge::Pair>());  
           //  cellList.push_front(&currentCell);  
   
             //! Sets the triangle vertices (actually just store the edge where it lives).  
             EdgePairList::iterator currentTriangleEdges  
                 = (*localTriangleListesEdges[objectNumber]).begin();  
   
             for(size_t i=0; i<cutEdges.size(); ++i) {  
                 Edge::Pair e=(*cutEdges[i]);  
                 (*currentTriangleEdges)[i] = e;  
             }  
455    
456              //! Stores the Edges on which the vertex will be computed.      if ((V[0]^V[1])*V[2] < 0) {
457              for (size_t i=0; i<cutEdges.size(); ++i)        std::swap(cutEdges[1], cutEdges[2]);
458                  (*verticesListesIntersection[objectNumber])[*cutEdges[i]] =      }
                     new Vertex(Internals::__splitEdge(*cutEdges[i], ref, S));  
459    
460              break;      //! Addes a triangle to the list.
461          }      (*localTriangleListesEdges[objectNumber]).push_front( TinyVector<3,Edge::Pair>());
462          case 4: {      //  cellList.push_front(&currentCell);
             //! Orders the Edges using a map !  
             std::map<Edge::Pair, const Edge*> cutEdgesMap;  
             for(size_t ie=0; ie<cutEdges.size(); ++ie)  
                 cutEdgesMap[(Edge::Pair)*(cutEdges[ie])] = (cutEdges[ie]);  
   
             //! Addes a triangle to the list.  
             (*localTriangleListesEdges[objectNumber]).push_front(TinyVector<3,Edge::Pair>());  
   
          //   cellList.push_front(&currentCell);  
   
             std::vector<TinyVector<3> > V(cutEdges.size());  
             std::map<Edge::Pair, const Edge*>::iterator currentCutEdge=cutEdgesMap.begin();  
             for (size_t i=0; i<V.size(); ++i) {  
                 const Edge& e = *(*currentCutEdge).second;  
                 if(e(0).Ref().isSet(ref))  
                     V[i] = (e(0).Vector3() - e(1).Vector3());  
                 else  
                     V[i] = (e(1).Vector3() - e(0).Vector3());  
                 ++currentCutEdge;  
             }  
463    
464              bool swapEdge = ((V[0]^V[1])*V[2] > 0);      //! Sets the triangle vertices (actually just store the edge where it lives).
465        EdgePairList::iterator currentTriangleEdges
466          = (*localTriangleListesEdges[objectNumber]).begin();
467    
468              //! Sets the triangle vertices (actually just store the edge where it lives).      for(size_t i=0; i<cutEdges.size(); ++i) {
469              EdgePairList::iterator currentTriangleEdges        Edge::Pair e=(*cutEdges[i]);
470                  = (*localTriangleListesEdges[objectNumber]).begin();        (*currentTriangleEdges)[i] = e;
471                            }
             currentCutEdge = cutEdgesMap.begin();  
             if (!swapEdge) {  
                 for(size_t i=0; i<3; ++i) {  
                     Edge::Pair e=(currentCutEdge->first);  
                     (*currentTriangleEdges)[i] = e;  
                     ++currentCutEdge;  
                 }  
             }  
             else {  
                 for(int i=2; i>=0; --i) {  
                     Edge::Pair e=(currentCutEdge->first);  
                     (*currentTriangleEdges)[i] = e;  
                     ++currentCutEdge;  
                 }  
             }  
             for (size_t i=0; i<4; ++i) {  
                 (*verticesListesIntersection[objectNumber])[*cutEdges[i]] =  
                 new Vertex(__splitEdge(*cutEdges[i], ref, S));  
472    
473              }      //! Stores the Edges on which the vertex will be computed.
474        for (size_t i=0; i<cutEdges.size(); ++i) {
475          (*verticesListesIntersection[objectNumber])[*cutEdges[i]] =
476            new Vertex(Internals::__splitEdge(*cutEdges[i], ref, S));
477        }
478        break;
479      }
480      case 4: {
481        //! Orders the Edges using a map !
482        std::map<Edge::Pair, const Edge*> cutEdgesMap;
483        for(size_t ie=0; ie<cutEdges.size(); ++ie) {
484          cutEdgesMap[(Edge::Pair)*(cutEdges[ie])] = (cutEdges[ie]);
485        }
486    
487        //! Addes a triangle to the list.
488        (*localTriangleListesEdges[objectNumber]).push_front(TinyVector<3,Edge::Pair>());
489    
490        //   cellList.push_front(&currentCell);
491    
492        std::vector<TinyVector<3> > V(cutEdges.size());
493        std::map<Edge::Pair, const Edge*>::iterator currentCutEdge=cutEdgesMap.begin();
494        for (size_t i=0; i<V.size(); ++i) {
495          const Edge& e = *(*currentCutEdge).second;
496          if(e(0).Ref().isSet(ref)) {
497            V[i] = (e(0).Vector3() - e(1).Vector3());
498          } else {
499            V[i] = (e(1).Vector3() - e(0).Vector3());
500          }
501          ++currentCutEdge;
502        }
503    
504        bool swapEdge = ((V[0]^V[1])*V[2] > 0);
505    
506        //! Sets the triangle vertices (actually just store the edge where it lives).
507        EdgePairList::iterator currentTriangleEdges
508          = (*localTriangleListesEdges[objectNumber]).begin();
509    
510        currentCutEdge = cutEdgesMap.begin();
511        if (!swapEdge) {
512          for(size_t i=0; i<3; ++i) {
513            Edge::Pair e=(currentCutEdge->first);
514            (*currentTriangleEdges)[i] = e;
515            ++currentCutEdge;
516          }
517        } else {
518          for(int i=2; i>=0; --i) {
519            Edge::Pair e=(currentCutEdge->first);
520            (*currentTriangleEdges)[i] = e;
521            ++currentCutEdge;
522          }
523        }
524        for (size_t i=0; i<4; ++i) {
525          (*verticesListesIntersection[objectNumber])[*cutEdges[i]] =
526            new Vertex(__splitEdge(*cutEdges[i], ref, S));
527        }
528    
529              //! Addes a triangle to the list.      //! Addes a triangle to the list.
530              (*localTriangleListesEdges[objectNumber]).push_front( TinyVector<3,Edge::Pair>());      (*localTriangleListesEdges[objectNumber]).push_front( TinyVector<3,Edge::Pair>());
531    
532              //cellList.push_front(&currentCell);      //cellList.push_front(&currentCell);
533    
534              currentTriangleEdges      currentTriangleEdges
535                  = (*localTriangleListesEdges[objectNumber]).begin();        = (*localTriangleListesEdges[objectNumber]).begin();
536                            
537              //! Sets the triangle vertices (actually just store the edge where it lives).      //! Sets the triangle vertices (actually just store the edge where it lives).
538              currentCutEdge = cutEdgesMap.begin();      currentCutEdge = cutEdgesMap.begin();
539              if (!swapEdge) {      if (!swapEdge) {
540                  for(int i=2; i>=0; --i) {        for(int i=2; i>=0; --i) {
541                      ++currentCutEdge;          ++currentCutEdge;
542                    //  (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);          //  (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
543                      Edge::Pair e=(currentCutEdge->first);          Edge::Pair e=(currentCutEdge->first);
544                      (*currentTriangleEdges)[i] = e;          (*currentTriangleEdges)[i] = e;
545                  }        }
546              }      } else {
547              else {        for(size_t i=0; i<3; ++i) {
548                  for(size_t i=0; i<3; ++i) {          ++currentCutEdge;
549                      ++currentCutEdge;          // (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
550                     // (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);          Edge::Pair e=(currentCutEdge->first);
551                      Edge::Pair e=(currentCutEdge->first);          (*currentTriangleEdges)[i] = e;
552                      (*currentTriangleEdges)[i] = e;        }
                 }  
             }  
             break;  
         }  
         case 1:  
         case 2: {  
             ffout(3) << "?!\n";  
             break;  
         }  
553      }      }
554        break;
555      }
556      case 1:
557      case 2: {
558        ffout(3) << "?!\n";
559        break;
560      }
561      }
562  }  }
563    
564  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
565  __createCutEdges(const Scene& scene,  __createCutEdges(const Scene& scene,
566                   std::vector<ReferenceCounting< std::vector<const Edge*> > > & cutEdges,                   std::vector<ReferenceCounting< std::vector<const Edge*> > > & cutEdges,
567                   ObjectReferences& otherReferences,                   ObjectReferences& otherReferences,
568                   const std::vector<TetrahedronByEdges>::iterator& currentT) {                   const std::vector<TetrahedronByEdges>::iterator& currentT)
569    {
570    
571      //size_t nobjectToTreat = otherReferences.size();      //size_t nobjectToTreat = otherReferences.size();
572      for(size_t ie=0; ie<6; ++ie) {
573      for(size_t ie=0; ie<6; ++ie) {      const Edge& e = currentT->edge(ie);
574          const Edge& e = currentT->edge(ie);      int num=0;
575          int num=0;      for(ObjectReferences::iterator i = otherReferences.begin();
576          for(ObjectReferences::iterator i = otherReferences.begin();          i != otherReferences.end(); ++i) {
577                       i != otherReferences.end(); ++i) {        for(size_t n = 0; n<scene.nbObjects(); ++n) {
578              for(size_t n = 0; n<scene.nbObjects(); ++n) {          Reference::Representation otherRef;
579                  Reference::Representation otherRef;          otherRef.resize(scene.nbObjects());
580                  otherRef.resize(scene.nbObjects());          if (scene.object(n).reference() == ((*i).first)) {
581                  if (scene.object(n).reference() == ((*i).first)) {            otherRef[n] = true;
582                      otherRef[n] = true;            if ((e(0).Ref().isSet(otherRef)) xor (e(1).Ref().isSet(otherRef))) {
583                      if ((e(0).Ref().isSet(otherRef)) xor (e(1).Ref().isSet(otherRef))) {              (*cutEdges[num]).push_back(&e);
584                          (*cutEdges[num]).push_back(&e);            }
                     }  
                 }  
             }  
             num++;  
585          }          }
586          }
587          num++;
588      }      }
589      }
590  }  }
591    
   
592  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
593  __calculatePointsIntersection(const Scene& scene,  __calculatePointsIntersection(const Scene& scene,
594                                size_t no,                                size_t no,
595                                size_t no2) {                                size_t no2)
596      size_t nobject=scene.nbObjects();  {
597      size_t nobject=scene.nbObjects();
598    
599      if(no< nobject and no2<nobject) {    if(no< nobject and no2<nobject) {
600          ListOfHexaToTriangle::iterator hexait=listOfHexaToTriangle.begin();      ListOfHexaToTriangle::iterator hexait=listOfHexaToTriangle.begin();
601          ListOfHexaToTriangle::iterator hexait0=hexait;      ListOfHexaToTriangle::iterator hexait0=hexait;
602          ++hexait;      ++hexait;
603          ListOfHexaToTriangle::iterator hexait1=hexait;      ListOfHexaToTriangle::iterator hexait1=hexait;
604          MapCellTriangle hexaToTriangle0=*(*hexait0);      MapCellTriangle hexaToTriangle0=*(*hexait0);
605          MapCellTriangle hexaToTriangle1=*(*hexait1);      MapCellTriangle hexaToTriangle1=*(*hexait1);
606    
607          for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)      for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin();
608          {          i!=toTreatHexahedra.end(); ++i) {
609              MapCellTriangle::iterator it0=hexaToTriangle0.find(*i);        MapCellTriangle::iterator it0=hexaToTriangle0.find(*i);
610              MapCellTriangle::iterator it1=hexaToTriangle1.find(*i);        MapCellTriangle::iterator it1=hexaToTriangle1.find(*i);
611              if(it0!=hexaToTriangle0.end() and it1!=hexaToTriangle1.end()) {        if(it0!=hexaToTriangle0.end() and it1!=hexaToTriangle1.end()) {
612                     std::list<Triangle*> listTriangle0=(*it0).second;          std::list<Triangle*> listTriangle0=(*it0).second;
613                     std::list<Triangle*> listTriangle1=(*it1).second;          std::list<Triangle*> listTriangle1=(*it1).second;
614              for(std::list<Triangle*>::iterator jt0=listTriangle0.begin() ;          for(std::list<Triangle*>::iterator jt0=listTriangle0.begin() ;
615           jt0!=listTriangle0.end() ; ++jt0) {              jt0!=listTriangle0.end() ; ++jt0) {
616                      for(std::list<Triangle*>::iterator jt1=listTriangle1.begin() ;            for(std::list<Triangle*>::iterator jt1=listTriangle1.begin() ;
617             jt1!=listTriangle1.end() ; ++jt1) {                jt1!=listTriangle1.end() ; ++jt1) {
618                          __createPointIntersection(no,no2,jt0,jt1);              __createPointIntersection(no,no2,jt0,jt1);
619                      }            }
                 }  
                   
             }  
             else {  
                 //ffout(0)<<"pas d'intersection dans l'hexa\n";  
             }  
620          }          }
621          } else {
622            //ffout(0)<<"pas d'intersection dans l'hexa\n";
623          }
624      }      }
625      }
626  }  }
627    
628  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
629  __calculateListesIntersection(const Scene& scene,  __calculateListesIntersection(const Scene& scene,
630                                ObjectReferences& otherReferences) {                                ObjectReferences& otherReferences) {
631      size_t nobject = scene.nbObjects();    size_t nobject = scene.nbObjects();
632      size_t nobjectToTreat = otherReferences.size();    size_t nobjectToTreat = otherReferences.size();
     if(nobjectToTreat==0)  
         nobjectToTreat=1;  
     for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)  
     {  
         SurfaceMeshGenerator::Internals::HexaTetraAssociation::iterator  
                  j=hexaTetraAssociation.find(*i);  
         std::vector<TetrahedronByEdges>& T = *(((*j).second).first);  
         const CartesianHexahedron& currentCell  
             = static_cast<const CartesianHexahedron&>(*(*i));  
   
         triangleListesIntersectionEdgesPerHexa[*i].resize(nobject);  
         for(size_t no=0 ; no<nobjectToTreat ; ++no) {  
             triangleListesIntersectionEdgesPerHexa[*i][no]=  
             new EdgePairList;  
         }  
633    
634          for (std::vector<TetrahedronByEdges>::iterator currentT = T.begin();    if(nobjectToTreat==0) {
635        currentT < T.end(); ++currentT) {      nobjectToTreat=1;
636      }
637    
638              std::vector<ReferenceCounting< std::vector<const Edge*> > > cutEdges;    for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin();
639              cutEdges.resize(nobject);        i!=toTreatHexahedra.end(); ++i) {
640              for(size_t no=0 ; no<nobject ; ++no) {      SurfaceMeshGenerator::Internals::HexaTetraAssociation::iterator
641                  cutEdges[no]= new std::vector<const Edge*>;        j=hexaTetraAssociation.find(*i);
642              }      std::vector<TetrahedronByEdges>& T = *(((*j).second).first);
643        const CartesianHexahedron& currentCell
644          = static_cast<const CartesianHexahedron&>(*(*i));
645    
646              __createCutEdges(scene,cutEdges,otherReferences,currentT);      triangleListesIntersectionEdgesPerHexa[*i].resize(nobject);
647        for(size_t no=0 ; no<nobjectToTreat ; ++no) {
648          triangleListesIntersectionEdgesPerHexa[*i][no]= new EdgePairList;
649        }
650    
651              int num=0;      for (std::vector<TetrahedronByEdges>::iterator currentT = T.begin();
652              for(ObjectReferences::iterator i = otherReferences.begin();           currentT < T.end(); ++currentT) {
          i != otherReferences.end(); ++i) {  
                 for(size_t n = 0; n<scene.nbObjects(); ++n) {  
                     Reference::Representation otherRef;  
                     otherRef.resize(scene.nbObjects());  
                     if (scene.object(n).reference() == ((*i).first)) {  
                         otherRef[n] = true;  
                             __createListIntersection(*(scene.object(n).shape()),  
                                 (*cutEdges[num]),otherRef,num,currentCell);  
                     }  
                 }  
                 num++;  
             }  
         }  
653    
654          for(size_t no=0 ; no<nobjectToTreat ; ++no) {        std::vector<ReferenceCounting< std::vector<const Edge*> > > cutEdges;
655              for (EdgePairList::iterator t = (*localTriangleListesEdges[no]).begin();        cutEdges.resize(nobject);
656            t != (*localTriangleListesEdges[no]).end(); ++t) {        for(size_t no=0 ; no<nobject ; ++no) {
657                  (*triangleListesIntersectionEdgesPerHexa[*i][no]).push_front(*t);          cutEdges[no]= new std::vector<const Edge*>;
658          }
659    
660              }        __createCutEdges(scene,cutEdges,otherReferences,currentT);
         }  
661    
662          for(size_t no=0 ; no<nobject ; ++no) {        int num=0;
663              (*localTriangleListesEdges[no]).clear();        for(ObjectReferences::iterator i = otherReferences.begin();
664              i != otherReferences.end(); ++i) {
665            for(size_t n = 0; n<scene.nbObjects(); ++n) {
666              Reference::Representation otherRef;
667              otherRef.resize(scene.nbObjects());
668              if (scene.object(n).reference() == ((*i).first)) {
669                otherRef[n] = true;
670                __createListIntersection(*(scene.object(n).shape()),
671                                         (*cutEdges[num]),otherRef,num,currentCell);
672              }
673          }          }
674            num++;
675          }
676      }      }
677    
678  }      for(size_t no=0 ; no<nobjectToTreat ; ++no) {
679          for (EdgePairList::iterator t = (*localTriangleListesEdges[no]).begin();
680               t != (*localTriangleListesEdges[no]).end(); ++t) {
681            (*triangleListesIntersectionEdgesPerHexa[*i][no]).push_front(*t);
682          }
683        }
684    
685  bool SurfaceMeshGenerator::Internals::__isDegenerate(const Triangle& T) {      for(size_t no=0 ; no<nobject ; ++no) {
686          (*localTriangleListesEdges[no]).clear();
687        }
688      }
689    }
690    
691      if(&(T(0))==&(T(1)) or &(T(0))==&(T(2)) or &(T(1))==&(T(2)))  bool SurfaceMeshGenerator::Internals::__isDegenerate(const Triangle& T)
692          return true;  {
693      else    return (&(T(0))==&(T(1)) or &(T(0))==&(T(2)) or &(T(1))==&(T(2)));
         return false;  
694  }  }
695    
696  TinyVector<3,const Vertex*> SurfaceMeshGenerator::Internals::__determinatePinterNotInList(const TriangleCut& K,  TinyVector<3,const Vertex*>
697                                                                                            TinyVector<3,size_t> ncase) {  SurfaceMeshGenerator::Internals::__determinatePinterNotInList(const TriangleCut& K,
698      TinyVector<3,const Vertex*> ABC;                                                                TinyVector<3,size_t> ncase)
699      ABC[0]=(K.pointsIntersection[ncase[0]][0]);  {
700      if(K.pointsIntersection[ncase[1]][0]!=ABC[0]) {    TinyVector<3,const Vertex*> ABC;
701          ABC[1]=(K.pointsIntersection[ncase[1]][0]);    ABC[0]=(K.pointsIntersection[ncase[0]][0]);
702      }    if(K.pointsIntersection[ncase[1]][0]!=ABC[0]) {
703      else {      ABC[1]=(K.pointsIntersection[ncase[1]][0]);
704          ABC[1]=(K.pointsIntersection[ncase[1]][1]);    } else {
705      }      ABC[1]=(K.pointsIntersection[ncase[1]][1]);
706      if(K.pointsIntersection[ncase[2]][0]!=ABC[0] and K.pointsIntersection[ncase[2]][0]!=ABC[1]) {    }
707          ABC[2]=(K.pointsIntersection[ncase[2]][0]);  
708      }    if(K.pointsIntersection[ncase[2]][0]!=ABC[0] and K.pointsIntersection[ncase[2]][0]!=ABC[1]) {
709      else {      ABC[2]=(K.pointsIntersection[ncase[2]][0]);
710          ABC[2]=(K.pointsIntersection[ncase[2]][1]);    } else {
711      }      ABC[2]=(K.pointsIntersection[ncase[2]][1]);
712      return ABC;    }
713    
714      return ABC;
715  }  }
716    
717  void SurfaceMeshGenerator::Internals::__casePinterNotInList(const TriangleCut& K,  void SurfaceMeshGenerator::Internals::__casePinterNotInList(const TriangleCut& K,
718                                                              const Cell* cell){                                                              const Cell* cell){
719      TinyVector<3,const Vertex*> ABC;    TinyVector<3,const Vertex*> ABC;
720      TinyVector<3,size_t> ncase;    TinyVector<3,size_t> ncase;
721      if(K.pointsIntersection[0][0]==K.pointsIntersection[0][1]) {    if(K.pointsIntersection[0][0]==K.pointsIntersection[0][1]) {
722          ncase[0]=0;      ncase[0]=0;
723          ncase[1]=1;      ncase[1]=1;
724          ncase[2]=2;      ncase[2]=2;
725          ABC=__determinatePinterNotInList(K,ncase);      ABC=__determinatePinterNotInList(K,ncase);
726      }    } else {
727      else if(K.pointsIntersection[1][0]==K.pointsIntersection[1][1]) {      if(K.pointsIntersection[1][0]==K.pointsIntersection[1][1]) {
728          ncase[0]=1;        ncase[0]=1;
729          ncase[1]=0;        ncase[1]=0;
730          ncase[2]=2;        ncase[2]=2;
731          ABC=__determinatePinterNotInList(K,ncase);        ABC=__determinatePinterNotInList(K,ncase);
732        } else {
733      }        if(K.pointsIntersection[2][0]==K.pointsIntersection[2][1]) {
     else if(K.pointsIntersection[2][0]==K.pointsIntersection[2][1]) {  
734          ncase[0]=2;          ncase[0]=2;
735          ncase[1]=1;          ncase[1]=1;
736          ncase[2]=0;          ncase[2]=0;
737          ABC=__determinatePinterNotInList(K,ncase);          ABC=__determinatePinterNotInList(K,ncase);
738          } else {
     }  
     else {  
739          ABC[0]=K.pointsIntersection[0][0];          ABC[0]=K.pointsIntersection[0][0];
740          ABC[1]=K.pointsIntersection[1][0];          ABC[1]=K.pointsIntersection[1][0];
741          ABC[2]=K.pointsIntersection[2][0];          ABC[2]=K.pointsIntersection[2][0];
742          }
743      }      }
744      __addTriangle(*ABC[0],*ABC[1],*ABC[2],cell);    }
745      __addTriangle(*ABC[0],*ABC[1],*ABC[2],cell);
746  }  }
747    
748  void SurfaceMeshGenerator::Internals::plotini(size_t nobjectToTreat) {  void SurfaceMeshGenerator::Internals::plotini(size_t nobjectToTreat) {
749      std::ofstream sortie2("ini.mesh");    std::ofstream sortie2("ini.mesh");
750      sortie2<<"MeshVersionFormatted 1"<<" \n";    sortie2<<"MeshVersionFormatted 1"<<" \n";
751      sortie2<<" \n";    sortie2<<" \n";
752      sortie2<<"Dimension"<<" \n";    sortie2<<"Dimension"<<" \n";
753      sortie2<<"3"<<" \n";    sortie2<<"3"<<" \n";
754      sortie2<<" \n";    sortie2<<" \n";
755      sortie2<<"Vertices"<<" \n";    sortie2<<"Vertices"<<" \n";
756      size_t nn=0,compt=0;    size_t nn=0,compt=0;
757      TinyVector<10,size_t> num;    TinyVector<10,size_t> num;
758      for(size_t numObj=0 ; numObj<nobjectToTreat ; ++numObj) {    for(size_t numObj=0 ; numObj<nobjectToTreat ; ++numObj) {
759          nn+=(*triangleListes [numObj]).size();      nn+=(*triangleListes [numObj]).size();
760          if(compt<10)      if(compt<10) {
761              num[compt]=nn;        num[compt]=nn;
762          ++compt;      }
763      }      ++compt;
764      sortie2<<3*nn<<"\n";    }
765      sortie2<<3*nn<<"\n";
766      for(size_t numObj=0 ; numObj<nobjectToTreat ; ++numObj) {  
767          for (SurfaceMeshGenerator::Internals::EdgePairList::iterator    for(size_t numObj=0 ; numObj<nobjectToTreat ; ++numObj) {
768        i = (*triangleListes[numObj]).begin();      for (SurfaceMeshGenerator::Internals::EdgePairList::iterator
769        i != (*triangleListes [numObj]).end(); ++i) {             i = (*triangleListes[numObj]).begin();
770              SurfaceMeshGenerator::Internals::VerticesList& V = (*verticesListes[numObj]);           i != (*triangleListes [numObj]).end(); ++i) {
771              const Vertex& V0=* V[(*i)[0]];        SurfaceMeshGenerator::Internals::VerticesList& V = (*verticesListes[numObj]);
772              const Vertex& V1=* V[(*i)[1]];        const Vertex& V0=* V[(*i)[0]];
773              const Vertex& V2=* V[(*i)[2]];        const Vertex& V1=* V[(*i)[1]];
774              sortie2<<V0.Vector3()[0]<<" "<<V0.Vector3()[1]<<" "<<V0.Vector3()[2]<<" 0 \n";        const Vertex& V2=* V[(*i)[2]];
775              sortie2<<V1.Vector3()[0]<<" "<<V1.Vector3()[1]<<" "<<V1.Vector3()[2]<<" 0 \n";        sortie2<<V0.Vector3()[0]<<" "<<V0.Vector3()[1]<<" "<<V0.Vector3()[2]<<" 0 \n";
776              sortie2<<V2.Vector3()[0]<<" "<<V2.Vector3()[1]<<" "<<V2.Vector3()[2]<<" 0 \n";        sortie2<<V1.Vector3()[0]<<" "<<V1.Vector3()[1]<<" "<<V1.Vector3()[2]<<" 0 \n";
777          }        sortie2<<V2.Vector3()[0]<<" "<<V2.Vector3()[1]<<" "<<V2.Vector3()[2]<<" 0 \n";
778      }      }
779      sortie2<<" \n";    }
780      sortie2<<"Triangles"<<"\n";    sortie2<<" \n";
781      sortie2<<nn<<"\n";    sortie2<<"Triangles"<<"\n";
782      sortie2<<nn<<"\n";
783      for(size_t nj=0;nj<nn;++nj) {  
784          sortie2<<3*nj+1<<" "<<3*nj+2<<" "<<3*nj+3;    for(size_t nj=0;nj<nn;++nj) {
785          if(nj<num[0])      sortie2<<3*nj+1<<" "<<3*nj+2<<" "<<3*nj+3;
786              sortie2<<" 2 \n";      if(nj<num[0])
787          else if(nj<num[1])        sortie2<<" 2 \n";
788              sortie2<<" 3 \n";      else if(nj<num[1])
789          else if(nj<num[2])        sortie2<<" 3 \n";
790              sortie2<<" 4 \n";      else if(nj<num[2])
791          else if(nj<num[3])        sortie2<<" 4 \n";
792              sortie2<<" 5 \n";      else if(nj<num[3])
793          else if(nj<num[4])        sortie2<<" 5 \n";
794              sortie2<<" 6 \n";      else if(nj<num[4])
795          else        sortie2<<" 6 \n";
796              sortie2<<" 7 \n";      else
797      }        sortie2<<" 7 \n";
798      sortie2<<"End"<<" \n";    }
799      sortie2<<"End"<<" \n";
800  }  }
801    
802  void SurfaceMeshGenerator::Internals::plotliste() {  void SurfaceMeshGenerator::Internals::plotliste() {
803      std::ofstream sortie2("liste2.mesh");    std::ofstream sortie2("liste2.mesh");
804      sortie2<<"MeshVersionFormatted 1"<<" \n";    sortie2<<"MeshVersionFormatted 1"<<" \n";
805      sortie2<<" \n";    sortie2<<" \n";
806      sortie2<<"Dimension"<<" \n";    sortie2<<"Dimension"<<" \n";
807      sortie2<<"3"<<" \n";    sortie2<<"3"<<" \n";
808      sortie2<<" \n";    sortie2<<" \n";
809      sortie2<<"Vertices"<<" \n";    sortie2<<"Vertices"<<" \n";
810      size_t nn=0,compt=0;    size_t nn=0,compt=0;
811      TinyVector<10,size_t> num;    TinyVector<10,size_t> num;
812      for(std::list<ReferenceCounting<Vector<Triangle> > >::iterator      i=(listOfTriangleListIntersection).begin();    for(std::list<ReferenceCounting<Vector<Triangle> > >::iterator i=(listOfTriangleListIntersection).begin();
813          i!=(listOfTriangleListIntersection).end();++i) {        i!=(listOfTriangleListIntersection).end();++i) {
814          nn+=(*(*i)).size();      nn+=(*(*i)).size();
815          if(compt<10)      if(compt<10)
816              num[compt]=nn;        num[compt]=nn;
817          ++compt;      ++compt;
818      }    }
819      sortie2<<3*nn<<"\n";    sortie2<<3*nn<<"\n";
820      for(std::list<ReferenceCounting<Vector<Triangle> > >::iterator      i=(listOfTriangleListIntersection).begin();    for(std::list<ReferenceCounting<Vector<Triangle> > >::iterator        i=(listOfTriangleListIntersection).begin();
821          i!=(listOfTriangleListIntersection).end();++i) {        i!=(listOfTriangleListIntersection).end();++i) {
822          Vector<Triangle> triangleListIntersection=*(*i);      Vector<Triangle> triangleListIntersection=*(*i);
823          for(size_t size=0 ; size<triangleListIntersection.size() ; ++size) {      for(size_t size=0 ; size<triangleListIntersection.size() ; ++size) {
824              const Triangle& T=triangleListIntersection[size];        const Triangle& T=triangleListIntersection[size];
825              const Vertex& V0=T(0);        const Vertex& V0=T(0);
826              const Vertex& V1=T(1);        const Vertex& V1=T(1);
827              const Vertex& V2=T(2);        const Vertex& V2=T(2);
828              sortie2<<V0.Vector3()[0]<<" "<<V0.Vector3()[1]<<" "<<V0.Vector3()[2]<<" 0 \n";        sortie2<<V0.Vector3()[0]<<" "<<V0.Vector3()[1]<<" "<<V0.Vector3()[2]<<" 0 \n";
829              sortie2<<V1.Vector3()[0]<<" "<<V1.Vector3()[1]<<" "<<V1.Vector3()[2]<<" 0 \n";        sortie2<<V1.Vector3()[0]<<" "<<V1.Vector3()[1]<<" "<<V1.Vector3()[2]<<" 0 \n";
830              sortie2<<V2.Vector3()[0]<<" "<<V2.Vector3()[1]<<" "<<V2.Vector3()[2]<<" 0 \n";        sortie2<<V2.Vector3()[0]<<" "<<V2.Vector3()[1]<<" "<<V2.Vector3()[2]<<" 0 \n";
831          }      }
832      }    }
833      sortie2<<" \n";    sortie2<<" \n";
834      sortie2<<"Triangles"<<"\n";    sortie2<<"Triangles"<<"\n";
835      sortie2<<nn<<"\n";    sortie2<<nn<<"\n";
836    
837      for(size_t nj=0;nj<nn;++nj) {    for(size_t nj=0; nj<nn; ++nj) {
838          sortie2<<3*nj+1<<" "<<3*nj+2<<" "<<3*nj+3;      sortie2 << 3*nj+1 << " " << 3*nj+2 << " " << 3*nj+3;
839          if(nj<num[0])      if(nj<num[0]) {
840              sortie2<<"  "<<2<<" \n";        sortie2 << "  " << 2 << " \n";
841          else if(nj<num[1])      } else {
842              sortie2<<" 3 \n";        if(nj<num[1]) {
843          else if(nj<num[2])          sortie2 << " 3 \n";
844              sortie2<<" 4 \n";        } else {
845          else if(nj<num[3])          if(nj<num[2]) {
846              sortie2<<" 4 \n";
847            } else {
848              if(nj<num[3]) {
849              sortie2<<" 5 \n";              sortie2<<" 5 \n";
850          else if(nj<num[4])            } else {
851              sortie2<<" 6 \n";              if(nj<num[4]) {
852          else                sortie2<<" 6 \n";
853              sortie2<<" 7 \n";              } else {
854                  sortie2<<" 7 \n";
855                }
856              }
857            }
858          }
859      }      }
860      sortie2<<"End"<<" \n";    }
861      sortie2 << "End\n";
862  }  }
863    
864  size_t SurfaceMeshGenerator::Internals::  size_t SurfaceMeshGenerator::Internals::
865  __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 )
866      real_t det,eps=-1e-10;  {
867      const TinyVector<3>& V02=triangleVector[0];    real_t det,eps=-1e-10;
868      const TinyVector<3>& V12=triangleVector[1];    const TinyVector<3>& V02=triangleVector[0];
869      const TinyVector<3>& P0=triangleVector[2];    const TinyVector<3>& V12=triangleVector[1];
870      det=V02[num[0]]*V12[num[1]]-V02[num[1]]*V12[num[0]];    const TinyVector<3>& P0=triangleVector[2];
871      det=V02[num[0]]*V12[num[1]]-V02[num[1]]*V12[num[0]];
872      if(det<eps or det>-eps) {  
873          TinyMatrix<2,2> mat;    if(det<eps or det>-eps) {
874          mat(0,0)=V02[num[0]];      TinyMatrix<2,2> mat;
875          mat(1,0)=V02[num[1]];      mat(0,0)=V02[num[0]];
876          mat(0,1)=V12[num[0]];      mat(1,0)=V02[num[1]];
877          mat(1,1)=V12[num[1]];      mat(0,1)=V12[num[0]];
878          TinyVector<2> bary,secondmb;      mat(1,1)=V12[num[1]];
879          secondmb[0]=P0[num[0]];      TinyVector<2> bary,secondmb;
880          secondmb[1]=P0[num[1]];      secondmb[0]=P0[num[0]];
881          bary=mat.invert()*secondmb;      secondmb[1]=P0[num[1]];
882          if(bary[0]>=eps && bary[1]>=eps && (1.-bary[0]-bary[1])>=eps) {      bary=mat.invert()*secondmb;
883              return 1;      if(bary[0]>=eps && bary[1]>=eps && (1.-bary[0]-bary[1])>=eps) {
884          }        return 1;
885          else {      } else {
886              return 0;        return 0;
         }  
887      }      }
888      else    } else {
889          return 1;      return 1;
890      }
891  }  }
892    
893  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
894  __InTriangle(const Vertex& P,  __InTriangle(const Vertex& P,
895               TinyVector<3,const Vertex* > & V,               TinyVector<3,const Vertex* > & V,
896               size_t& compt) {               size_t& compt)
897    {
898      compt=0;    compt=0;
899      TinyVector<3,TinyVector<3,real_t> >  triangleVector;    TinyVector<3,TinyVector<3,real_t> >  triangleVector;
900      triangleVector[0]=-(*V[2]).Vector3() + (*V[0]).Vector3();    triangleVector[0]=-(*V[2]).Vector3() + (*V[0]).Vector3();
901      triangleVector[1]= -(*V[2]).Vector3() + (*V[1]).Vector3();    triangleVector[1]= -(*V[2]).Vector3() + (*V[1]).Vector3();
902      triangleVector[2]=P.Vector3()-(*V[2]).Vector3();    triangleVector[2]=P.Vector3()-(*V[2]).Vector3();
903      TinyVector<2,size_t>  num;    TinyVector<2,size_t>  num;
904                                                                    
905      num[0]=0;    num[0]=0;
906      num[1]=1;    num[1]=1;
907      size_t ajout=__IsInTriangle(num,triangleVector);  
908      compt+=ajout;    size_t ajout=__IsInTriangle(num,triangleVector);
909      compt+=ajout;
910      num[0]=0;  
911      num[1]=2;    num[0]=0;
912      ajout=__IsInTriangle(num,triangleVector);    num[1]=2;
913      compt+=ajout;    ajout=__IsInTriangle(num,triangleVector);
914      compt+=ajout;
915      num[0]=1;  
916      num[1]=2;    num[0]=1;
917      ajout=__IsInTriangle(num,triangleVector);    num[1]=2;
918      compt+=ajout;    ajout=__IsInTriangle(num,triangleVector);
919          compt+=ajout;  
920  }  }
921    
922  TinyVector<4, real_t>  SurfaceMeshGenerator::Internals::  TinyVector<4, real_t>  SurfaceMeshGenerator::Internals::
923  __EquationPlan(const TinyVector<3,const Vertex *> & V){  __EquationPlan(const TinyVector<3,const Vertex *> & V){
924    
925      //equation of plan      //equation of plan
926      TinyVector<3> V01=(*V[1]).Vector3()-(*V[0]).Vector3();    TinyVector<3> V01=(*V[1]).Vector3()-(*V[0]).Vector3();
927      TinyVector<3> V12=(*V[2]).Vector3()-(*V[1]).Vector3();    TinyVector<3> V12=(*V[2]).Vector3()-(*V[1]).Vector3();
928      TinyVector<3> V02=(*V[2]).Vector3()-(*V[0]).Vector3();    TinyVector<3> V02=(*V[2]).Vector3()-(*V[0]).Vector3();
929    
930      TinyVector<3> coeffo;    TinyVector<3> coeffo;
931      coeffo=V01^V02;    coeffo=V01^V02;
932      real_t d2 = -(coeffo*(*V[0]).Vector3());    real_t d2 = -(coeffo*(*V[0]).Vector3());
933    
934      TinyVector<4> coeff;    TinyVector<4> coeff;
935      coeff[0]=coeffo[0];    coeff[0]=coeffo[0];
936      coeff[1]=coeffo[1];    coeff[1]=coeffo[1];
937      coeff[2]=coeffo[2];    coeff[2]=coeffo[2];
938      coeff[3]=d2;    coeff[3]=d2;
939    
940      return coeff;    return coeff;
       
941  }  }
942    
943  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
# Line 928  __DeterminatePoint(const size_t& ncase, Line 946  __DeterminatePoint(const size_t& ncase,
946                     const Vertex* & Vo1,                     const Vertex* & Vo1,
947                     TinyVector<3,const Vertex* > & V,                     TinyVector<3,const Vertex* > & V,
948                     IntersectionPoints& Pinter) {                     IntersectionPoints& Pinter) {
949      std::pair<const Vertex*,const Vertex*> pV(Vo0,Vo1);    std::pair<const Vertex*,const Vertex*> pV(Vo0,Vo1);
950      std::pair<const Vertex*,const Vertex*> pV2(Vo1,Vo0);    std::pair<const Vertex*,const Vertex*> pV2(Vo1,Vo0);
951    
952      real_t epsilon=1e-10;    real_t epsilon=1e-10;
953      real_t pscalar,k;    real_t pscalar,k;
954      size_t num;    size_t num;
955      Vertex P;    Vertex P;
956      TinyVector<4> coeffo=__EquationPlan(V);    TinyVector<4> coeffo=__EquationPlan(V);
957      TinyVector<3> coeff;    TinyVector<3> coeff;
958      coeff[0]=coeffo[0];    coeff[0]=coeffo[0];
959      coeff[1]=coeffo[1];    coeff[1]=coeffo[1];
960      coeff[2]=coeffo[2];    coeff[2]=coeffo[2];
961      real_t d=coeffo[3];    real_t d=coeffo[3];
962    
963      TinyVector<3> Vo01 = (*Vo1).Vector3() - (*Vo0).Vector3();    TinyVector<3> Vo01 = (*Vo1).Vector3() - (*Vo0).Vector3();
964      TinyVector<3> V01 = (*V[1]).Vector3() - (*V[0]).Vector3();    TinyVector<3> V01 = (*V[1]).Vector3() - (*V[0]).Vector3();
965      TinyVector<3> V12 = (*V[2]).Vector3() - (*V[1]).Vector3();    TinyVector<3> V12 = (*V[2]).Vector3() - (*V[1]).Vector3();
966      TinyVector<3> V02 = (*V[2]).Vector3() - (*V[0]).Vector3();    TinyVector<3> V02 = (*V[2]).Vector3() - (*V[0]).Vector3();
967            
968      real_t norm1=Norm(Vo01);    real_t norm1=Norm(Vo01);
969      pscalar=coeff*Vo01;    pscalar=coeff*Vo01;
970      num=0;    num=0;
971      if(std::abs(pscalar)>1e-10) {    if(std::abs(pscalar)>1e-10) {
972          k=-(coeff*(*Vo0).Vector3()+d)/pscalar;      k=-(coeff*(*Vo0).Vector3()+d)/pscalar;
973          P=(*Vo0)+k*Vo01;      P=(*Vo0)+k*Vo01;
974                    
975          TinyVector<3> V0P=P.Vector3()-(*V[0]).Vector3();      TinyVector<3> V0P=P.Vector3()-(*V[0]).Vector3();
976          TinyVector<3> V1P=P.Vector3()-(*V[1]).Vector3();      TinyVector<3> V1P=P.Vector3()-(*V[1]).Vector3();
977          TinyVector<3> V2P=P.Vector3()-(*V[2]).Vector3();      TinyVector<3> V2P=P.Vector3()-(*V[2]).Vector3();
978    
979          real_t norm2=Norm(P.Vector3()-(*Vo0).Vector3());      real_t norm2=Norm(P.Vector3()-(*Vo0).Vector3());
980          real_t norm3=Norm(P.Vector3()-(*Vo1).Vector3());      real_t norm3=Norm(P.Vector3()-(*Vo1).Vector3());
981          real_t inV01=Norm(V01^V0P);      real_t inV01=Norm(V01^V0P);
982          real_t inV12=Norm(V12^V1P);      real_t inV12=Norm(V12^V1P);
983          real_t inV02=Norm(V02^V2P);      real_t inV02=Norm(V02^V2P);
984    
985          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 ) {
986              if(ncase<3) {        if(ncase<3) {
987                  __InTriangle(P,V,num);          __InTriangle(P,V,num);
988                  if(num==3) {          if(num==3) {
989                      if(!(norm2>epsilon)) {            if(!(norm2>epsilon)) {
990                          Pinter.addPoint(Vo0);              Pinter.addPoint(Vo0);
991                      }            } else
992                      else if(!(norm3>epsilon)) {              if(!(norm3>epsilon)) {
993                          Pinter.addPoint(Vo1);                Pinter.addPoint(Vo1);
994                } else
995                  if(listVertexIntersectionNew.find(pV)==listVertexIntersectionNew.end()
996                     and listVertexIntersectionNew.find(pV2)==listVertexIntersectionNew.end()) {
997                    listVertexIntersectionNew[pV][0]=new Vertex (P.Vector3());
998                    listVertexIntersectionNew[pV][1]=new Vertex (P.Vector3());
999                    Pinter.addPoint(listVertexIntersectionNew[pV][0]);
1000                  } else {
1001                    if(listVertexIntersectionNew.find(pV)!=listVertexIntersectionNew.end()) {
1002                      if(Norm((*(*listVertexIntersectionNew.find(pV)).second[0]).Vector3()
1003                              -P.Vector3())<1e-6) {
1004                        Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[0]);
1005                      } else
1006                        if(Norm((*(*listVertexIntersectionNew.find(pV)).second[1]).Vector3()
1007                                -P.Vector3())<1e-6) {
1008                          Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[1]);
1009                        } else {
1010                          listVertexIntersectionNew[pV][1]=new Vertex (P.Vector3());
1011                          Pinter.addPoint(listVertexIntersectionNew[pV][1]);
1012                        }
1013                    } else {
1014                      if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[0]).Vector3()
1015                              -P.Vector3())<1e-6) {
1016                        Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[0]);
1017                      } else
1018                        if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[1]).Vector3()
1019                                -P.Vector3())<1e-6) {
1020                          Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[1]);
1021                        } else {
1022                          listVertexIntersectionNew[pV2][1]=new Vertex (P.Vector3());
1023                          Pinter.addPoint(listVertexIntersectionNew[pV2][1]);
1024                      }                      }
                     else if(listVertexIntersectionNew.find(pV)==listVertexIntersectionNew.end()  
          and listVertexIntersectionNew.find(pV2)==listVertexIntersectionNew.end()) {  
                         listVertexIntersectionNew[pV][0]=new Vertex (P.Vector3());  
                         listVertexIntersectionNew[pV][1]=new Vertex (P.Vector3());  
                         Pinter.addPoint(listVertexIntersectionNew[pV][0]);  
                     }  
                     else {  
                         if(listVertexIntersectionNew.find(pV)!=listVertexIntersectionNew.end()) {  
                             if(Norm((*(*listVertexIntersectionNew.find(pV)).second[0]).Vector3()  
            -P.Vector3())<1e-6) {  
                                 Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[0]);  
                             }  
                             else if(Norm((*(*listVertexIntersectionNew.find(pV)).second[1]).Vector3()  
                 -P.Vector3())<1e-6) {  
                                 Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[1]);  
                             }  
                             else {  
                                 listVertexIntersectionNew[pV][1]=new Vertex (P.Vector3());  
                                 Pinter.addPoint(listVertexIntersectionNew[pV][1]);  
                             }  
                         }  
                         else {  
                             if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[0]).Vector3()  
              -P.Vector3())<1e-6) {  
                             Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[0]);  
                         }  
                             else if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[1]).Vector3()  
                 -P.Vector3())<1e-6) {  
                                 Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[1]);  
                             }  
                             else {  
                                 listVertexIntersectionNew[pV2][1]=new Vertex (P.Vector3());  
                                 Pinter.addPoint(listVertexIntersectionNew[pV2][1]);  
                             }  
                         }  
                     }  
                     Pinter.cutTriangle1[ncase]=Pinter.size();  
                       
1025                  }                  }
1026              }                }
1027              else if(inV01> epsilon and inV02> epsilon and inV12> epsilon ) {            Pinter.cutTriangle1[ncase]=Pinter.size();
1028                  __InTriangle(P,V,num);          }
1029                  if(num==3) {        } else {
1030                      //ffout(0)<<"ajout du point cas 2\n";          if(inV01 > epsilon and inV02 > epsilon and inV12 > epsilon ) {
1031                      //ffout(0)<<P.Vector3()[0]<<" "<<P.Vector3()[1]<<" "<<P.Vector3()[2]<<" \n";            __InTriangle(P,V,num);
1032                      //listVertexIntersectionNew.push_back(new Vertex (P.Vector3()));            if(num==3) {
1033                      PIntersection PIn;              //ffout(0)<<"ajout du point cas 2\n";
1034                      if(!(norm2>epsilon)) {              //ffout(0)<<P.Vector3()[0]<<" "<<P.Vector3()[1]<<" "<<P.Vector3()[2]<<" \n";
1035                          Pinter.addPoint(Vo0);              //listVertexIntersectionNew.push_back(new Vertex (P.Vector3()));
1036                      }              PIntersection PIn;
1037                      else if(!(norm3>epsilon)) {              if(!(norm2>epsilon)) {
1038                          Pinter.addPoint(Vo1);                Pinter.addPoint(Vo0);
1039                      }              } else {
1040                      else if(listVertexIntersectionNew.find(pV)==listVertexIntersectionNew.end()                if(!(norm3>epsilon)) {
1041           and listVertexIntersectionNew.find(pV2)==listVertexIntersectionNew.end()) {                  Pinter.addPoint(Vo1);
1042                          listVertexIntersectionNew[pV][0]=new Vertex (P.Vector3());                } else {
1043                    if(listVertexIntersectionNew.find(pV)==listVertexIntersectionNew.end()
1044                       and listVertexIntersectionNew.find(pV2)==listVertexIntersectionNew.end()) {
1045                      listVertexIntersectionNew[pV][0]=new Vertex (P.Vector3());
1046                      listVertexIntersectionNew[pV][1]=new Vertex (P.Vector3());
1047                      Pinter.addPoint(listVertexIntersectionNew[pV][0]);
1048                    } else {
1049                      if(listVertexIntersectionNew.find(pV)!=listVertexIntersectionNew.end()) {
1050                        if(Norm((*(*listVertexIntersectionNew.find(pV)).second[0]).Vector3()
1051                                -P.Vector3())<1e-6) {
1052                          Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[0]);
1053                        } else {
1054                          if(Norm((*(*listVertexIntersectionNew.find(pV)).second[1]).Vector3()
1055                                  -P.Vector3())<1e-6) {
1056                            Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[1]);
1057                          } else {
1058                          listVertexIntersectionNew[pV][1]=new Vertex (P.Vector3());                          listVertexIntersectionNew[pV][1]=new Vertex (P.Vector3());
1059                          Pinter.addPoint(listVertexIntersectionNew[pV][0]);                          Pinter.addPoint(listVertexIntersectionNew[pV][1]);
1060                          }
1061                      }                      }
1062                      else {                    } else {
1063                          if(listVertexIntersectionNew.find(pV)!=listVertexIntersectionNew.end()) {                      if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[0]).Vector3()
1064                              if(Norm((*(*listVertexIntersectionNew.find(pV)).second[0]).Vector3()                              -P.Vector3())<1e-6) {
1065                 -P.Vector3())<1e-6) {                        Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[0]);
1066                                  Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[0]);                      } else {
1067                              }                        if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[1]).Vector3()
1068                              else if(Norm((*(*listVertexIntersectionNew.find(pV)).second[1]).Vector3()                                -P.Vector3())<1e-6) {
1069                      -P.Vector3())<1e-6) {                          Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[1]);
1070                                  Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[1]);                        } else {
1071                              }                          listVertexIntersectionNew[pV2][1]=new Vertex (P.Vector3());
1072                              else {                          Pinter.addPoint(listVertexIntersectionNew[pV2][1]);
1073                                  listVertexIntersectionNew[pV][1]=new Vertex (P.Vector3());                        }
                                 Pinter.addPoint(listVertexIntersectionNew[pV][1]);  
                             }  
                         }  
                         else {  
                             if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[0]).Vector3()  
                -P.Vector3())<1e-6) {  
                                 Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[0]);  
                             }  
                             else if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[1]).Vector3()  
                     -P.Vector3())<1e-6) {  
                                 Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[1]);  
                             }  
                             else {  
                                 listVertexIntersectionNew[pV2][1]=new Vertex (P.Vector3());  
                                 Pinter.addPoint(listVertexIntersectionNew[pV2][1]);  
                             }  
                         }  
                           
1074                      }                      }
1075                     // Pinter.addPoint(listVertexIntersectionNew[listVertexIntersectionNew.size()-1]);                    }
1076                  }                  }
1077                  }
1078              }              }
1079                // Pinter.addPoint(listVertexIntersectionNew[listVertexIntersectionNew.size()-1]);
1080              }
1081          }          }
1082          }
1083      }      }
1084      }
1085  }  }
1086    
1087  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1088  __CalculPoint(const size_t& objectNumber,const size_t& numobject,  __CalculPoint(const size_t& objectNumber,const size_t& numobject,
1089                std::list<Triangle*>::iterator & currentTriangle,                std::list<Triangle*>::iterator & currentTriangle,
1090                std::list<Triangle*>::iterator & currentTriangle1,                std::list<Triangle*>::iterator & currentTriangle1,
1091                IntersectionPoints& Pinter) {                IntersectionPoints& Pinter)
1092    {
1093      size_t num1=objectNumber,num0=numobject;    size_t num1=objectNumber,num0=numobject;
       
     TinyVector<3,const Vertex * > V012;  
     const Triangle& T0=*(*currentTriangle1);  
     V012[0] = &T0(0);  
     V012[1] = &T0(1);  
     V012[2] = &T0(2);  
     TinyVector<3,const Vertex* > Vo012;  
     Triangle& To0=*(*currentTriangle);  
     Vo012[0] = &To0(0);  
     Vo012[1] = &To0(1);  
     Vo012[2] = &To0(2);  
     //for each edge of the triangle V0V1V2  
     Pinter.numpair[0]=num0;  
     Pinter.numpair[1]=num1;  
   
     __DeterminatePoint(0,(V012[0]),(V012[1]),Vo012,Pinter);  
     __DeterminatePoint(1,(V012[0]),(V012[2]),Vo012,Pinter);  
     __DeterminatePoint(2,(V012[1]),(V012[2]),Vo012,Pinter);  
   
     __DeterminatePoint(3,(Vo012[0]),(Vo012[1]),V012,Pinter);  
     __DeterminatePoint(3,(Vo012[0]),(Vo012[2]),V012,Pinter);  
     __DeterminatePoint(3,(Vo012[1]),(Vo012[2]),V012,Pinter);  
1094    
1095      TinyVector<3,const Vertex * > V012;
1096      const Triangle& T0=*(*currentTriangle1);
1097      V012[0] = &T0(0);
1098      V012[1] = &T0(1);
1099      V012[2] = &T0(2);
1100      TinyVector<3,const Vertex* > Vo012;
1101      Triangle& To0=*(*currentTriangle);
1102      Vo012[0] = &To0(0);
1103      Vo012[1] = &To0(1);
1104      Vo012[2] = &To0(2);
1105      //for each edge of the triangle V0V1V2
1106      Pinter.numpair[0]=num0;
1107      Pinter.numpair[1]=num1;
1108    
1109      __DeterminatePoint(0,(V012[0]),(V012[1]),Vo012,Pinter);
1110      __DeterminatePoint(1,(V012[0]),(V012[2]),Vo012,Pinter);
1111      __DeterminatePoint(2,(V012[1]),(V012[2]),Vo012,Pinter);
1112    
1113      __DeterminatePoint(3,(Vo012[0]),(Vo012[1]),V012,Pinter);
1114      __DeterminatePoint(3,(Vo012[0]),(Vo012[2]),V012,Pinter);
1115      __DeterminatePoint(3,(Vo012[1]),(Vo012[2]),V012,Pinter);
1116  }  }
1117    
1118  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1119  __inout(const size_t& numobject,  __inout(const size_t& numobject,
1120          const int & ttn0,          const int & ttn0,
1121          const Vertex*& V0,          const Vertex*& V0,
1122          const Vertex*& V1) {          const Vertex*& V1)
1123    {
1124      size_t& re0=(*edgesRefVertex[numobject])[V0];    size_t& re0=(*edgesRefVertex[numobject])[V0];
1125      size_t& re1=(*edgesRefVertex[numobject])[V1];    size_t& re1=(*edgesRefVertex[numobject])[V1];
1126    
1127      if(re0+re1==4 and ttn0==1) {    if(re0+re1==4 and ttn0==1) {
1128          ffout(0)<<"warning edge ref 2 cut \n";      ffout(0)<<"warning edge ref 2 cut \n";
1129          ffout(0)<<(*V0).Vector3()[0]<<" "<<(*V0).Vector3()[1]<<" "<<(*V0).Vector3()[2]<<" \n";      ffout(0)<<(*V0).Vector3()[0]<<" "<<(*V0).Vector3()[1]<<" "<<(*V0).Vector3()[2]<<" \n";
1130          ffout(0)<<(*V1).Vector3()[0]<<" "<<(*V1).Vector3()[1]<<" "<<(*V1).Vector3()[2]<<" \n";      ffout(0)<<(*V1).Vector3()[0]<<" "<<(*V1).Vector3()[1]<<" "<<(*V1).Vector3()[2]<<" \n";
1131          re0=0;      re0=0;
1132          re1=0;      re1=0;
1133      }    } else {
1134      else if(re0+re1==2 and ttn0==0) {      // CECILE: ne serait-ce pas plutôt une assertion ?
1135          re0=2;      if(re0+re1==2 and ttn0==0) {
1136          re1=2;        re0=2;
1137          re1=2;
1138      }      }
1139      }
1140  }  }
1141    
1142  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1143  __createPointIntersection(const size_t& numobject,const size_t& objectNumber,  __createPointIntersection(const size_t& numobject,const size_t& objectNumber,
1144                            std::list<Triangle*>::iterator currentTriangle1,                            std::list<Triangle*>::iterator currentTriangle1,
1145                            std::list<Triangle*>::iterator currentTriangle) {                            std::list<Triangle*>::iterator currentTriangle)
1146    {    
1147          //ffout(0)<< "createpointintersection "<<numobject<<" "<<objectNumber<<"\n";
1148      //ffout(0)<< "createpointintersection "<<numobject<<" "<<objectNumber<<"\n";    //ffout(0)<<"*************************\n";
     //ffout(0)<<"*************************\n";  
1149    
1150      PairTriangleNew pTriangle;    PairTriangleNew pTriangle;
1151      pTriangle.first=(*currentTriangle1);    pTriangle.first=(*currentTriangle1);
1152      pTriangle.second=(*currentTriangle);    pTriangle.second=(*currentTriangle);
1153    
1154      IntersectionPoints Pinter;    IntersectionPoints Pinter;
1155            
1156      __CalculPoint(objectNumber,numobject,currentTriangle,currentTriangle1,Pinter);    __CalculPoint(objectNumber,numobject,currentTriangle,currentTriangle1,Pinter);
1157    
1158      const TinyVector<3,size_t >& pintercut=Pinter.cutTriangle1;    const TinyVector<3,size_t >& pintercut=Pinter.cutTriangle1;
1159    
1160      if(pintercut[0]+pintercut[1]+pintercut[2]!=0) {    if(pintercut[0]+pintercut[1]+pintercut[2]!=0) {
1161          //ffout(0)<<"le tr est cut\n";      //ffout(0)<<"le tr est cut\n";
1162          if(pairTriangleToIntersectionPoints.find(pTriangle)==pairTriangleToIntersectionPoints.end()) {      if(pairTriangleToIntersectionPoints.find(pTriangle)==pairTriangleToIntersectionPoints.end()) {
1163              //ffout(0) << "rajout de pTriangle\n";        //ffout(0) << "rajout de pTriangle\n";
1164              pairTriangleToIntersectionPoints[pTriangle] = Pinter;        pairTriangleToIntersectionPoints[pTriangle] = Pinter;
         }  
1165      }      }
1166      else if(Pinter.size()!=0 and numobject<objectNumber) {    } else {
1167          //ffout(0)<<"rajout2 \n";      if(Pinter.size()!=0 and numobject<objectNumber) {
1168          pairTriangleToIntersectionPoints[pTriangle] = Pinter;        //ffout(0)<<"rajout2 \n";
1169          pairTriangleToIntersectionPoints[pTriangle] = Pinter;
1170      }      }
1171      }
1172    
1173      if(pairTriangleToIntersectionPoints.find(pTriangle)!= pairTriangleToIntersectionPoints.end()    if(pairTriangleToIntersectionPoints.find(pTriangle)!= pairTriangleToIntersectionPoints.end()
1174         and numobject<objectNumber) {       and numobject<objectNumber) {
1175          TinyVector<3,size_t >& cutTr2= pairTriangleToIntersectionPoints[pTriangle].cutTriangle2;      TinyVector<3,size_t >& cutTr2= pairTriangleToIntersectionPoints[pTriangle].cutTriangle2;
1176    
1177          for(size_t i=0 ; i<Pinter.size() ; ++i) {      for(size_t i=0 ; i<Pinter.size() ; ++i) {
1178              real_t epsilon=1e-10;        real_t epsilon=1e-10;
1179              const Triangle& To=(*(*currentTriangle));        const Triangle& To=(*(*currentTriangle));
1180              const Vertex& Vo0 = To(0);        const Vertex& Vo0 = To(0);
1181              const Vertex& Vo1 = To(1);        const Vertex& Vo1 = To(1);
1182              const Vertex& Vo2 = To(2);        const Vertex& Vo2 = To(2);
1183              const Vertex& P=*(pairTriangleToIntersectionPoints[pTriangle].listofPoints[i]);        const Vertex& P=*(pairTriangleToIntersectionPoints[pTriangle].listofPoints[i]);
1184                            
1185              real_t Vo01=Norm(Vo1.Vector3()-Vo0.Vector3());        real_t Vo01=Norm(Vo1.Vector3()-Vo0.Vector3());
1186              real_t Vo12=Norm(Vo2.Vector3()-Vo1.Vector3());        real_t Vo12=Norm(Vo2.Vector3()-Vo1.Vector3());
1187              real_t Vo02=Norm(Vo2.Vector3()-Vo0.Vector3());        real_t Vo02=Norm(Vo2.Vector3()-Vo0.Vector3());
1188              real_t norm0=Norm(P.Vector3()-Vo0.Vector3());        real_t norm0=Norm(P.Vector3()-Vo0.Vector3());
1189              real_t norm1=Norm(P.Vector3()-Vo1.Vector3());        real_t norm1=Norm(P.Vector3()-Vo1.Vector3());
1190              real_t norm2=Norm(P.Vector3()-Vo2.Vector3());        real_t norm2=Norm(P.Vector3()-Vo2.Vector3());
1191              //for filter points        //for filter points
1192              real_t eps=1e-10;        real_t eps=1e-10;
1193              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) {
1194                  cutTr2[0]=i+1;          cutTr2[0]=i+1;
1195              }        }
1196              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) {
1197                  cutTr2[1]=i+1;          cutTr2[1]=i+1;
1198              }        }
1199              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) {
1200                  cutTr2[2]=i+1;          cutTr2[2]=i+1;
1201              }        }
         }  
1202      }      }
1203      }
1204  }        }      
1205    
1206    
1207  const Vertex& SurfaceMeshGenerator::Internals::  const Vertex& SurfaceMeshGenerator::Internals::
1208  __determinateNumber(const int& p,  __determinateNumber(const int& p,
1209                      const TinyVector<3,TinyVector<2,const Vertex*> >& pointsint) {                      const TinyVector<3,TinyVector<2,const Vertex*> >& pointsint)
1210      if(p<=2)  {
1211          return *pointsint[p][0];    if(p<=2) {
1212      else      return *pointsint[p][0];
1213          return *pointsint[2][0];    } else {
1214        return *pointsint[2][0];
1215      }
1216  }  }
1217    
1218  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1219  __addTriangle(const Vertex& V0,  __addTriangle(const Vertex& V0,
1220                const Vertex& V1,                const Vertex& V1,
1221                const Vertex& V2,                const Vertex& V2,
1222                const Cell* cell) {                const Cell* cell)
1223      Triangle T(V0,V1,V2);  {
1224      T.setMother(cell);    Triangle T(V0,V1,V2);
1225      if(!__isDegenerate(T)) {    T.setMother(cell);
1226          triangleListIntersectionNew.push_back(T);    if(!__isDegenerate(T)) {
1227      }      triangleListIntersectionNew.push_back(T);
1228      }
1229  }  }
1230    
1231  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1232  __determinateCase(const int& ncase,  __determinateCase(const int& ncase,
1233                    const TriangleCut & K,                    const TriangleCut & K,
1234                    TinyVector<3,int>& tt,                    TinyVector<3,int>& tt,
1235                    TinyVector<3,int>& pp) {                    TinyVector<3,int>& pp)
1236    {
1237      const TinyVector<3,TinyVector<2,int> >& tt2=K.edgecut;    const TinyVector<3,TinyVector<2,int> >& tt2=K.edgecut;
1238      if(ncase==0) {    if(ncase==0) {
1239          tt[0]=tt2[0][0];      tt[0]=tt2[0][0];
1240          tt[1]=tt2[1][0];      tt[1]=tt2[1][0];
1241          tt[2]=tt2[2][0];      tt[2]=tt2[2][0];
1242          pp[0]=0;      pp[0]=0;
1243          pp[1]=1;      pp[1]=1;
1244          pp[2]=2;      pp[2]=2;
1245      }    } else {
1246      else if(ncase==1) {      if(ncase==1) {
1247          tt[0]=tt2[0][0];        tt[0]=tt2[0][0];
1248          tt[1]=tt2[2][0];        tt[1]=tt2[2][0];
1249          tt[2]=tt2[1][0];        tt[2]=tt2[1][0];
1250          pp[0]=0;        pp[0]=0;
1251          pp[1]=2;        pp[1]=2;
1252          pp[2]=1;        pp[2]=1;
1253      }      } else {
1254      else if(ncase==2) {        if(ncase==2) {
1255          tt[0]=tt2[1][0];          tt[0]=tt2[1][0];
1256          tt[1]=tt2[2][0];          tt[1]=tt2[2][0];
1257          tt[2]=tt2[0][0];          tt[2]=tt2[0][0];
1258          pp[0]=1;          pp[0]=1;
1259          pp[1]=2;          pp[1]=2;
1260            // CECILE: est-il normal que pp[2] ne soit pas défini ?
1261          }
1262      }      }
1263      }
1264  }  }
1265    
1266  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1267  __createCase1(int ncase,  __createCase1(int ncase,
1268                const TriangleCut& K,                const TriangleCut& K,
1269                const Cell* cell) {                const Cell* cell)
1270      const Vertex& V0=K(0);  {
1271      const Vertex& V1=K(1);    const Vertex& V0=K(0);
1272      const Vertex& V2=K(2);    const Vertex& V1=K(1);
1273      const Vertex& V2=K(2);
1274            
1275      const TinyVector<3,TinyVector<2,const Vertex*> >& pointsint=K.pointsIntersection;    const TinyVector<3,TinyVector<2,const Vertex*> >& pointsint=K.pointsIntersection;
1276      const std::map< const Vertex* ,const Vertex* >& pointsin=K.pointsin;    const std::map< const Vertex* ,const Vertex* >& pointsin=K.pointsin;
1277      const size_t& numberin=pointsin.size();    const size_t& numberin=pointsin.size();
1278      const Vertex& A=K(ncase);    const Vertex& A=K(ncase);
1279    
1280      TinyVector<3,int> tt;    TinyVector<3,int> tt;
1281      TinyVector<3,int> pp;    TinyVector<3,int> pp;
1282    
1283      __determinateCase(ncase,K,tt,pp);    __determinateCase(ncase,K,tt,pp);
1284    
1285      if(numberin==0) {    if(numberin==0) {
1286          if(tt[0]==1) {      if(tt[0]==1) {
1287              const Vertex& B=__determinateNumber(pp[0],pointsint);        const Vertex& B=__determinateNumber(pp[0],pointsint);
1288              if(tt[1]==1) {        if(tt[1]==1) {
1289                  const Vertex& C=__determinateNumber(pp[1],pointsint);          const Vertex& C=__determinateNumber(pp[1],pointsint);
1290                  __addTriangle(A,B,C,cell);          __addTriangle(A,B,C,cell);
1291              }        } else {
1292          else {          ffout(0) <<"pbs 10\n"<<std::flush;
1293              ffout(0) <<"pbs 10\n"<<std::flush;          ffout(0)<<V0.Vector3()<<"\n";
1294              ffout(0)<<V0.Vector3()<<"\n";          ffout(0)<<V1.Vector3()<<"\n";
1295              ffout(0)<<V1.Vector3()<<"\n";          ffout(0)<<V2.Vector3()<<"\n";
1296              ffout(0)<<V2.Vector3()<<"\n";        }
1297          }      } else {
1298      }        if(tt[1]==1) {
1299      else  {          ffout(0)<<"pbs 11\n"<<std::flush;
1300          if(tt[1]==1) {        } else {
1301              ffout(0)<<"pbs 11\n"<<std::flush;          ffout(0) <<"pbs 111\n";
1302          }        }
         else {  
             ffout(0) <<"pbs 111\n";  
         }  
1303      }      }
1304    }    } else {
1305    else if(numberin==1) {      if(numberin==1) {
1306        std::map<const Vertex*,const Vertex*>::const_iterator it=pointsin.begin();        std::map<const Vertex*,const Vertex*>::const_iterator it=pointsin.begin();
1307        const Vertex&  D=*((*it).second);        const Vertex&  D=*((*it).second);
1308        if(tt[0]==1) {        if(tt[0]==1) {
1309            const Vertex& B=__determinateNumber(pp[0],pointsint);          const Vertex& B=__determinateNumber(pp[0],pointsint);
1310            if(tt[1]==1) {          if(tt[1]==1) {
1311                const Vertex& C=__determinateNumber(pp[1],pointsint);            const Vertex& C=__determinateNumber(pp[1],pointsint);
1312                __addTriangle(A,B,D,cell);            __addTriangle(A,B,D,cell);
1313                __addTriangle(A,D,C,cell);            __addTriangle(A,D,C,cell);
1314            }          } else {
1315            else {            ffout(0)<<"-------- warning trou in 0\n"<<std::flush;
1316                ffout(0)<<"-------- warning trou in 0\n"<<std::flush;            ffout(0)<<V0.Vector3()<<"\n";
1317                ffout(0)<<V0.Vector3()<<"\n";            ffout(0)<<V1.Vector3()<<"\n";
1318                ffout(0)<<V1.Vector3()<<"\n";            ffout(0)<<V2.Vector3()<<"\n";
1319                ffout(0)<<V2.Vector3()<<"\n";            __addTriangle(A,B,D,cell);
1320                __addTriangle(A,B,D,cell);          }
1321            }        } else {
1322        }          if(tt[1]==1) {
       else if(tt[1]==1) {  
1323            ffout(0)<<"-------- warning trou in\n"<<std::flush;            ffout(0)<<"-------- warning trou in\n"<<std::flush;
1324            const Vertex& B=__determinateNumber(pp[1],pointsint);            const Vertex& B=__determinateNumber(pp[1],pointsint);
1325            __addTriangle(A,B,D,cell);            __addTriangle(A,B,D,cell);
1326            } else {
       }  
       else {  
1327            ffout(0) <<"pbs 1\n"<<std::flush;            ffout(0) <<"pbs 1\n"<<std::flush;
1328            }
1329        }        }
1330    }      } else {
   else {  
1331  #warning is it true?  #warning is it true?
1332      //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
1333      ffout(0) << numberin<<" be careful 1\n"<<std::flush;        ffout(0) << numberin<<" be careful 1\n"<<std::flush;
1334        //POUR STEPH :  peut generer des tr plats        //POUR STEPH :  peut generer des tr plats
1335      //  ffout(0)<<V0.Vector3()<<"\n";        //  ffout(0)<<V0.Vector3()<<"\n";
1336        //ffout(0)<<V1.Vector3()<<"\n";        //ffout(0)<<V1.Vector3()<<"\n";
1337        //ffout(0)<<V2.Vector3()<<"\n";        //ffout(0)<<V2.Vector3()<<"\n";
1338        TinyVector<10,const Vertex*> P,Ptemp;        TinyVector<10,const Vertex*> P,Ptemp;
1339        size_t i=0;        size_t i=0;
1340        for(std::map<const Vertex*,const Vertex*>::const_iterator it=pointsin.begin();        for(std::map<const Vertex*,const Vertex*>::const_iterator it=pointsin.begin();
1341            it!=pointsin.end() ; ++it) {            it!=pointsin.end() ; ++it) {
1342            P[i+1]=(*it).second;          P[i+1]=(*it).second;
1343            ++i;          ++i;
1344        }        }
1345    
1346        if(tt[0]==1) {        if(tt[0]==1) {
1347            const Vertex& B=__determinateNumber(pp[0],pointsint);          const Vertex& B=__determinateNumber(pp[0],pointsint);
1348            P[0]=&B;          P[0]=&B;
1349            if(tt[1]==1) {          if(tt[1]==1) {
1350                const Vertex& C=__determinateNumber(pp[1],pointsint);            const Vertex& C=__determinateNumber(pp[1],pointsint);
1351                size_t size=numberin+1;            size_t size=numberin+1;
1352              for(size_t boucle=0 ; boucle<numberin-1 ; ++boucle) {
1353                size_t num=1;
1354                real_t norm=Norm((*P[0]).Vector3()-(*P[1]).Vector3());
1355                for(size_t i=2 ; i<size ; ++i) {
1356                  real_t normtemp=Norm((*P[0]).Vector3()-(*P[i]).Vector3());
1357                  if(normtemp<=norm) {
1358                    norm=normtemp;
1359                    Ptemp[i-1]=P[num];
1360                    num=i;
1361                  } else {
1362                    Ptemp[i-1]=P[i];
1363                  }
1364                }
1365                __addTriangle(A,*P[0],*P[num],cell);
1366                Ptemp[0]=P[num];
1367                P=Ptemp;
1368                --size;
1369                // ffout(0)<<"size = "<<size<<"\n";
1370              }
1371              __addTriangle(A,*P[0],*P[1],cell);
1372              __addTriangle(A,*P[1],C,cell);
1373              /*  if(tt[2]==1) {
1374                  const Vertex& D=__determinateNumber(pp[2],pointsint);
1375                  for(size_t i=0 ; i<numberin ; ++i) {
1376                  P[i+1]=(pointsin[i]);
1377                  }
1378                  P[0]=&B;
1379                  size=numberin+1;
1380                for(size_t boucle=0 ; boucle<numberin-1 ; ++boucle) {                for(size_t boucle=0 ; boucle<numberin-1 ; ++boucle) {
1381                    size_t num=1;                size_t num=1;
1382                    real_t norm=Norm((*P[0]).Vector3()-(*P[1]).Vector3());                real_t norm=Norm((*P[0]).Vector3()-(*P[1]).Vector3());
1383                    for(size_t i=2 ; i<size ; ++i) {                for(size_t i=2 ; i<size ; ++i) {
1384                        real_t normtemp=Norm((*P[0]).Vector3()-(*P[i]).Vector3());                real_t normtemp=Norm((*P[0]).Vector3()-(*P[i]).Vector3());
1385                        if(normtemp<=norm) {                if(normtemp<=norm) {
1386                            norm=normtemp;                norm=normtemp;
1387                            Ptemp[i-1]=P[num];                Ptemp[i-1]=P[num];
1388                            num=i;                num=i;
                       }  
                       else {  
                           Ptemp[i-1]=P[i];  
                       }  
                   }  
                   __addTriangle(A,*P[0],*P[num],cell);  
                   Ptemp[0]=P[num];  
                   P=Ptemp;  
                   --size;  
                  // ffout(0)<<"size = "<<size<<"\n";  
1389                }                }
1390                __addTriangle(A,*P[0],*P[1],cell);                else {
1391                __addTriangle(A,*P[1],C,cell);                Ptemp[i-1]=P[i];
1392              /*  if(tt[2]==1) {                }
1393                    const Vertex& D=__determinateNumber(pp[2],pointsint);                }
1394                    for(size_t i=0 ; i<numberin ; ++i) {                __addTriangle(D,*P[0],*P[num],cell);
1395                        P[i+1]=(pointsin[i]);                Ptemp[0]=P[num];
1396                    }                P=Ptemp;
1397                    P[0]=&B;                --size;
1398                    size=numberin+1;                // ffout(0)<<"size = "<<size<<"\n";
1399                    for(size_t boucle=0 ; boucle<numberin-1 ; ++boucle) {                }
1400                        size_t num=1;                __addTriangle(D,*P[0],*P[1],cell);
1401                        real_t norm=Norm((*P[0]).Vector3()-(*P[1]).Vector3());                //  __addTriangle(D,*P[1],C,cell);
                       for(size_t i=2 ; i<size ; ++i) {  
                           real_t normtemp=Norm((*P[0]).Vector3()-(*P[i]).Vector3());  
                           if(normtemp<=norm) {  
                               norm=normtemp;  
                               Ptemp[i-1]=P[num];  
                               num=i;  
                           }  
                           else {  
                               Ptemp[i-1]=P[i];  
                           }  
                       }  
                       __addTriangle(D,*P[0],*P[num],cell);  
                       Ptemp[0]=P[num];  
                       P=Ptemp;  
                       --size;  
                       // ffout(0)<<"size = "<<size<<"\n";  
                   }  
                   __addTriangle(D,*P[0],*P[1],cell);  
                 //  __addTriangle(D,*P[1],C,cell);  
1402                }*/                }*/
1403            }          } else {
1404            else {            ffout(0)<<"une seule arete cut pbs\n";
1405                ffout(0)<<"une seule arete cut pbs\n";          }
1406            }        } else {
1407        }          ffout(0)<<"une seule arete cut pbs2\n";
       else {  
           ffout(0)<<"une seule arete cut pbs2\n";  
1408        }        }
1409        }
1410    }    }
1411  }  }
1412    
1413  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1414  __createCase2(int ncase,  __createCase2(int ncase,
1415                const TriangleCut& K,                const TriangleCut& K,
1416                const Cell* cell) {                const Cell* cell)
1417      {  
1418    int ttn0,ttn1;    int ttn0,ttn1;
1419    int p0,p1;    int p0,p1;
1420    const Vertex& V0=K(0);    const Vertex& V0=K(0);
# Line 1420  __createCase2(int ncase, Line 1433  __createCase2(int ncase,
1433        ttn1=tt2[2][0];        ttn1=tt2[2][0];
1434        p0=1;        p0=1;
1435        p1=2;        p1=2;
1436    }    } else {
1437    else if(ncase==1) {      if(ncase==1) {
1438      ttn0=tt2[0][0];        ttn0=tt2[0][0];
1439      ttn1=tt2[2][0];        ttn1=tt2[2][0];
1440      p0=0;        p0=0;
1441      p1=2;        p1=2;
1442      ncase0=0;        ncase0=0;
1443      ncase1=2;        ncase1=2;
1444    }      } else {
1445    else if(ncase==2) {        if(ncase==2) {
1446      ttn0=tt2[0][0];          ttn0=tt2[0][0];
1447      ttn1=tt2[1][0];          ttn1=tt2[1][0];
1448      p0=0;          p0=0;
1449      p1=1;          p1=1;
1450      ncase0=1;          ncase0=1;
1451      ncase1=2;          ncase1=2;
1452          }
1453        }
1454    }    }
1455    
1456    const Vertex& A=K(ncase0);    const Vertex& A=K(ncase0);
1457    const Vertex& B=K(ncase1);    const Vertex& B=K(ncase1);
1458    
1459    if(numberin==0) {    if(numberin==0) {
1460        if(ttn0==1) {      if(ttn0==1) {
1461            const Vertex& C=__determinateNumber(p0,pointsint);        const Vertex& C=__determinateNumber(p0,pointsint);
1462            if(ttn1==1) {        if(ttn1==1) {
1463                const Vertex& D=__determinateNumber(p1,pointsint);          const Vertex& D=__determinateNumber(p1,pointsint);
1464                __addTriangle(A,B,C,cell);          __addTriangle(A,B,C,cell);
1465                __addTriangle(C,B,D,cell);          __addTriangle(C,B,D,cell);
1466            }        } else {
1467            else {          // ffout(0) <<"pbs 20\n";
1468               // ffout(0) <<"pbs 20\n";          __addTriangle(A,B,C,cell);
               __addTriangle(A,B,C,cell);  
           }  
1469        }        }
1470        else {      } else {
1471            if(ttn1==1) {        if(ttn1==1) {
1472               // ffout(0) <<"pbs 210\n";          // ffout(0) <<"pbs 210\n";
1473                const Vertex& C=__determinateNumber(p1,pointsint);          const Vertex& C=__determinateNumber(p1,pointsint);
1474                __addTriangle(A,B,C,cell);          __addTriangle(A,B,C,cell);
1475          } else {
1476            }          ffout(0) <<"pbs 21\n"<<std::flush;
1477            else {          ffout(0)<<V0.Vector3()<<"\n";
1478                ffout(0) <<"pbs 21\n"<<std::flush;          ffout(0)<<V1.Vector3()<<"\n";
1479                ffout(0)<<V0.Vector3()<<"\n";          ffout(0)<<V2.Vector3()<<"\n";
               ffout(0)<<V1.Vector3()<<"\n";  
               ffout(0)<<V2.Vector3()<<"\n";  
           }  
1480        }        }
1481    }      }
1482    else if(numberin==1 /*or numberin==2*/) {    } else
1483        if(numberin==1 /*or numberin==2*/) {
1484        std::map<const Vertex*,const Vertex*>::const_iterator it=pointsin.begin();        std::map<const Vertex*,const Vertex*>::const_iterator it=pointsin.begin();
1485        const Vertex& E=*((*it).second);        const Vertex& E=*((*it).second);
1486    
1487        if(ttn0==1) {        if(ttn0==1) {
1488            const Vertex& C=__determinateNumber(p0,pointsint);          const Vertex& C=__determinateNumber(p0,pointsint);
1489            if(ttn1==1) {          if(ttn1==1) {
1490                const Vertex& D=__determinateNumber(p1,pointsint);            const Vertex& D=__determinateNumber(p1,pointsint);
1491                __addTriangle(A,C,E,cell);            __addTriangle(A,C,E,cell);
1492                __addTriangle(A,E,B,cell);            __addTriangle(A,E,B,cell);
1493                __addTriangle(B,E,D,cell);            __addTriangle(B,E,D,cell);
1494            }          } else {
1495            else {            __addTriangle(A,B,E,cell);
1496                __addTriangle(A,B,E,cell);            __addTriangle(C,E,B,cell);
1497                __addTriangle(C,E,B,cell);          }
1498            }        } else {
1499        }          if(ttn1==1) {
1500        else {            const Vertex& D=__determinateNumber(p1,pointsint);
1501            if(ttn1==1) {            __addTriangle(D,E,B,cell);
1502                const Vertex& D=__determinateNumber(p1,pointsint);            __addTriangle(A,B,E,cell);
1503                __addTriangle(D,E,B,cell);          } else {
1504                __addTriangle(A,B,E,cell);            __addTriangle(A,B,E,cell);
1505            }          }
1506            else {        }
1507                __addTriangle(A,B,E,cell);      } else
1508            }        if(numberin>=2) {
1509        }          //POUR STEPH :  peut generer des tr plats
1510    }          ffout(0) << numberin <<"be careful 2\n"<<std::flush;
1511    else if(numberin>=2) {          //     ffout(0)<<V0.Vector3()<<"\n";
1512        //POUR STEPH :  peut generer des tr plats          //   ffout(0)<<V1.Vector3()<<"\n";
1513        ffout(0) << numberin <<"be careful 2\n"<<std::flush;          // ffout(0)<<V2.Vector3()<<"\n";
1514   //     ffout(0)<<V0.Vector3()<<"\n";          TinyVector<10,const Vertex*> P,Ptemp;
1515     //   ffout(0)<<V1.Vector3()<<"\n";          size_t i=0;
1516       // ffout(0)<<V2.Vector3()<<"\n";          for(std::map<const Vertex*,const Vertex*>::const_iterator it=pointsin.begin();
1517        TinyVector<10,const Vertex*> P,Ptemp;              it!=pointsin.end() ; ++it) {
       size_t i=0;  
       for(std::map<const Vertex*,const Vertex*>::const_iterator it=pointsin.begin();  
           it!=pointsin.end() ; ++it) {  
1518            P[i+1]=(*it).second;            P[i+1]=(*it).second;
1519            ++i;            ++i;
1520        }          }
1521        if(ttn0==1) {  
1522            // CECILE: les deux 'if' qui suivent en sont-ils vraiment ou
1523            // veux-tu vérifier simplement que l'assertion est vérifiée ?
1524            // Dans ce cas, il serait bien préférable d'utiliser assert()
1525    
1526            if(ttn0==1) {
1527            const Vertex& C=__determinateNumber(p0,pointsint);            const Vertex& C=__determinateNumber(p0,pointsint);
1528            __addTriangle(A,B,C,cell);            __addTriangle(A,B,C,cell);
1529            P[0]=&C;            P[0]=&C;
1530            if(ttn1==1) {            if(ttn1==1) {
1531                const Vertex& D=__determinateNumber(p1,pointsint);              const Vertex& D=__determinateNumber(p1,pointsint);
1532                size_t size=numberin+1;              size_t size=numberin+1;
1533                for(size_t boucle=0 ; boucle<numberin-1 ; ++boucle) {              for(size_t boucle=0 ; boucle<numberin-1 ; ++boucle) {
1534                    size_t num=1;                size_t num=1;
1535                    real_t norm=Norm((*P[0]).Vector3()-(*P[1]).Vector3());                real_t norm=Norm((*P[0]).Vector3()-(*P[1]).Vector3());
1536                    for(size_t i=2 ; i<size ; ++i) {                for(size_t i=2 ; i<size ; ++i) {
1537                        real_t normtemp=Norm((*P[0]).Vector3()-(*P[i]).Vector3());                  real_t normtemp=Norm((*P[0]).Vector3()-(*P[i]).Vector3());
1538                        if(normtemp<=norm) {                  if(normtemp<=norm) {
1539                            norm=normtemp;                    norm=normtemp;
1540                            Ptemp[i-1]=P[num];                    Ptemp[i-1]=P[num];
1541                            num=i;                    num=i;
1542                        }                  } else {
1543                        else {                    Ptemp[i-1]=P[i];
1544                            Ptemp[i-1]=P[i];                  }
                       }  
                   }  
                   __addTriangle(B,*P[0],*P[num],cell);  
                   Ptemp[0]=P[num];  
                   P=Ptemp;  
                   --size;  
                   // ffout(0)<<"size = "<<size<<"\n";  
1545                }                }
1546                __addTriangle(B,*P[0],*P[1],cell);                __addTriangle(B,*P[0],*P[num],cell);
1547                __addTriangle(B,*P[1],D,cell);                Ptemp[0]=P[num];
1548                //__addTriangle(A,B,D,cell);                P=Ptemp;
1549            }                --size;
1550            else {                // ffout(0)<<"size = "<<size<<"\n";
1551                ffout(0)<<"une seule arete cut pbs\n";              }
1552                __addTriangle(B,*P[0],*P[1],cell);
1553                __addTriangle(B,*P[1],D,cell);
1554                //__addTriangle(A,B,D,cell);
1555              } else {
1556                ffout(0)<<"une seule arete cut pbs\n";
1557            }            }
1558        }          } else {
       else {  
1559            ffout(0)<<"une seule arete cut pbs2\n";            ffout(0)<<"une seule arete cut pbs2\n";
1560            }
1561        }        }
   }  
1562  }  }
1563    
1564  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1565  __createLocalListIntersection(const size_t& objectNumber,  __createLocalListIntersection(const size_t& objectNumber,
1566                                const Triangle*& currentTriangle,                                const Triangle*& currentTriangle,
1567                                TriangleCut & K,                                TriangleCut & K,
1568                                const Cell* cell) {                                const Cell* cell)
1569    {
1570      TinyVector<3,TinyVector<2,int> >& tt=K.edgecut;    TinyVector<3,TinyVector<2,int> >& tt=K.edgecut;
1571      TinyVector<3,TinyVector<2,const Vertex*> > & pointsIntersection=K.pointsIntersection;    TinyVector<3,TinyVector<2,const Vertex*> > & pointsIntersection=K.pointsIntersection;
1572      std::map<const Vertex* ,const Vertex*>& pointsin=K.pointsin;    std::map<const Vertex* ,const Vertex*>& pointsin=K.pointsin;
1573    
1574      int numberin=pointsin.size();    int numberin=pointsin.size();
1575      const size_t& numobject=K.numobject;    const size_t& numobject=K.numobject;
1576    
1577      ListOfHexaToTriangle::iterator hexait=listOfHexaToTriangle.begin();    ListOfHexaToTriangle::iterator hexait=listOfHexaToTriangle.begin();
1578      MapCellTriangle hexaToTriangle1;    MapCellTriangle hexaToTriangle1;
1579      if(numobject!=0)  
1580          hexaToTriangle1=*(*hexait);    if(numobject!=0) {
1581      else {      hexaToTriangle1=*(*hexait);
1582          ++hexait;    } else {
1583          hexaToTriangle1=*(*hexait);      ++hexait;
1584      }      hexaToTriangle1=*(*hexait);
1585      }
1586                    
1587      MapCellTriangle::iterator it1=hexaToTriangle1.find(cell);    MapCellTriangle::iterator it1=hexaToTriangle1.find(cell);
1588      if(it1!=hexaToTriangle1.end()) {    if(it1!=hexaToTriangle1.end()) {
1589      std::list<Triangle*> listTriangle1=(*it1).second;      std::list<Triangle*> listTriangle1=(*it1).second;
1590          for(std::list<Triangle*>::iterator jt0=listTriangle1.begin() ;      for(std::list<Triangle*>::iterator jt0=listTriangle1.begin() ;
1591       jt0!=listTriangle1.end() ; ++jt0) {          jt0!=listTriangle1.end() ; ++jt0) {
1592    
1593              PairTriangleNew pT;        PairTriangleNew pT;
1594              if(numobject<objectNumber) {        if(numobject<objectNumber) {
1595                  pT.first=(currentTriangle);          pT.first=(currentTriangle);
1596                  pT.second=(*jt0);          pT.second=(*jt0);
1597              }        } else {
1598              else {          pT.second=(currentTriangle);
1599                  pT.second=(currentTriangle);          pT.first=(*jt0);
1600                  pT.first=(*jt0);        }
1601              }        size_t temp=0;
             size_t temp=0;  
1602                            
1603              if(pairTriangleToIntersectionPoints.find(pT)!=pairTriangleToIntersectionPoints.end()) {        if(pairTriangleToIntersectionPoints.find(pT)!=pairTriangleToIntersectionPoints.end()) {
1604                  TinyVector<3,size_t > cc;          TinyVector<3,size_t > cc;
1605                  if(numobject<objectNumber) {          if(numobject<objectNumber) {
1606                      cc=(pairTriangleToIntersectionPoints[pT]).cutTriangle1;            cc=(pairTriangleToIntersectionPoints[pT]).cutTriangle1;
1607                  }          } else {
1608                  else {            cc=(pairTriangleToIntersectionPoints[pT]).cutTriangle2;
1609                      cc=(pairTriangleToIntersectionPoints[pT]).cutTriangle2;          }
1610                  }          //detection tr plat
1611                  //detection tr plat          size_t isDegenerate=0;
1612                  size_t isDegenerate=0;          if(__isDegenerate(*(*jt0))) {
1613                  if(__isDegenerate(*(*jt0))) {            // ffout(0)<<"degenerate triangle\n";
1614                     // ffout(0)<<"degenerate triangle\n";            isDegenerate=1;
1615                      isDegenerate=1;          }
1616                  }          if(pairTriangleToIntersectionPoints[pT].size()==2 and isDegenerate==0) {
1617                  if(pairTriangleToIntersectionPoints[pT].size()==2 and isDegenerate==0) {            const Vertex*& V0=pairTriangleToIntersectionPoints[pT].listofPoints[0];
1618                      const Vertex*& V0=pairTriangleToIntersectionPoints[pT].listofPoints[0];            const Vertex*& V1=pairTriangleToIntersectionPoints[pT].listofPoints[1];
1619                      const Vertex*& V1=pairTriangleToIntersectionPoints[pT].listofPoints[1];            if((V0==&(*(*jt0))(0) or V0==&(*(*jt0))(1) or V0==&(*(*jt0))(2))
1620                      if((V0==&(*(*jt0))(0) or V0==&(*(*jt0))(1) or V0==&(*(*jt0))(2))               and (V1==&(*(*jt0))(0) or V1==&(*(*jt0))(1) or V1==&(*(*jt0))(2))) {
1621           and (V1==&(*(*jt0))(0) or V1==&(*(*jt0))(1) or V1==&(*(*jt0))(2))) {              //ffout(0)<<"inter confondu avec arete\n";
1622                          //ffout(0)<<"inter confondu avec arete\n";              isDegenerate=1;
1623                          isDegenerate=1;            }
1624                      }          }
1625                  }          if(pairTriangleToIntersectionPoints[pT].size()>=2 and isDegenerate==0) {
1626                  if(pairTriangleToIntersectionPoints[pT].size()>=2 and isDegenerate==0) {            for(size_t b=0; b<pairTriangleToIntersectionPoints[pT].size() ; ++b) {
1627                      for(size_t b=0; b<pairTriangleToIntersectionPoints[pT].size() ; ++b) {              const Vertex*& V=pairTriangleToIntersectionPoints[pT].listofPoints[b];
1628                          const Vertex*& V=pairTriangleToIntersectionPoints[pT].listofPoints[b];              if(!(V==&(*(*jt0))(0) or V==&(*(*jt0))(1) or V==&(*(*jt0))(2)) or b==0) {
1629                          if(!(V==&(*(*jt0))(0) or V==&(*(*jt0))(1) or V==&(*(*jt0))(2)) or b==0) {                if(vertexInTriangles.find(V)==vertexInTriangles.end()) {
1630                          if(vertexInTriangles.find(V)==vertexInTriangles.end()) {                  TinyVector<2,const Triangle*> TT;
1631                              TinyVector<2,const Triangle*> TT;                  TT[0]=*jt0;
1632                              TT[0]=*jt0;                  TT[1]=*jt0;
1633                              TT[1]=*jt0;                  vertexInTriangles[V]=TT;
1634                              vertexInTriangles[V]=TT;                } else {
1635                          }                  (*vertexInTriangles.find(V)).second[1]=*jt0;
1636                          else {                }
1637                              (*vertexInTriangles.find(V)).second[1]=*jt0;                if(triangleWithVertex.find(*jt0)==triangleWithVertex.end()) {
1638                          }                  TinyVector<2,const Vertex*> VV;
1639                          if(triangleWithVertex.find(*jt0)==triangleWithVertex.end()) {                  VV[0]=V;
1640                              TinyVector<2,const Vertex*> VV;                  VV[1]=V;
1641                              VV[0]=V;                  triangleWithVertex[*jt0]=VV;
1642                              VV[1]=V;                } else {
1643                              triangleWithVertex[*jt0]=VV;                  (*triangleWithVertex.find(*jt0)).second[1]=V;
1644                          }                }
                         else {  
                             (*triangleWithVertex.find(*jt0)).second[1]=V;  
                         }  
                         }  
                     }  
                 }  
                 if(cc[0]!=0) {  
                     temp+=1;  
                     if(tt[0][0]==0) {  
                         tt[0][0]=1;  
                         pointsIntersection[0][0]=  
                         (pairTriangleToIntersectionPoints[pT].listofPoints[cc[0]-1]);  
                     }  
                     else {  
                         tt[0][1]=1;  
                         pointsIntersection[0][1]=  
                             (pairTriangleToIntersectionPoints[pT].listofPoints[cc[0]-1]);  
                     }  
                 }  
                 if(cc[1]!=0) {  
                     temp+=1;  
                     if(tt[1][0]==0) {  
                         tt[1][0]=1;  
                         pointsIntersection[1][0]=  
                             (pairTriangleToIntersectionPoints[pT].listofPoints[cc[1]-1]);  
                     }  
                     else{  
                         tt[1][1]=1;  
                         pointsIntersection[1][1]=  
                             (pairTriangleToIntersectionPoints[pT].listofPoints[cc[1]-1]);  
                     }  
                       
                 }  
                 if(cc[2]!=0) {  
                     temp+=1;  
                     if(tt[2][0]==0) {  
                         tt[2][0]=1;  
                         pointsIntersection[2][0]=  
                             (pairTriangleToIntersectionPoints[pT].listofPoints[cc[2]-1]);  
                     }  
                     else{  
                         tt[2][1]=1;  
                         pointsIntersection[2][1]=  
                             (pairTriangleToIntersectionPoints[pT].listofPoints[cc[2]-1]);  
                     }  
                       
                 }  
                 if(temp<pairTriangleToIntersectionPoints[pT].size()) {  
                     numberin+=1;  
                     size_t r=0;  
                     size_t sum=cc[0]+cc[1]+cc[2];  
                     if(sum>3)  
                         ffout(0)<<"sum "<<sum<<"\n"<<std::flush;  
                     if(sum==0) {  
                         r=0;  
                         for(size_t i=1 ; i<pairTriangleToIntersectionPoints[pT].size() ;++i)  
                             pointsin[(pairTriangleToIntersectionPoints[pT].listofPoints[i])]=  
                                 (pairTriangleToIntersectionPoints[pT].listofPoints[i]);  
                     }  
                     else if(sum==2 or sum==5)  
                         r=0;  
                     else if(sum==4 or sum==1)  
                         r=1;  
                     else if(sum==3)  
                         r=2;  
                     else {  
                         ffout(0)<<"pbs sum = "<<sum<<" size "<<  
                         pairTriangleToIntersectionPoints[pT].size()<<" \n";  
                         ffout(0)<<K(0)<<"\n";  
                         ffout(0)<<K(1)<<"\n";  
                         ffout(0)<<K(2)<<"\n";  
                         ffout(0)<<cc[0]<<" "<<cc[1]<<" "<<cc[2]<<"\n";  
                     }  
                     if(pointsin.find(pairTriangleToIntersectionPoints[pT].listofPoints[r])  
          ==pointsin.end())  
                         pointsin[(pairTriangleToIntersectionPoints[pT].listofPoints[r])]  
                             =(pairTriangleToIntersectionPoints[pT].listofPoints[r]);  
                 }  
1645              }              }
1646              //ffout(0)<<cc[0]<<" "<<cc[1]<<" "<<cc[2]<<"\n";            }
1647       //ffout(0)<<((listVertexIntersectionPerPair[objectp])[pT]).size()<<"\n";;          }
1648              else {          if(cc[0]!=0) {
1649             // ffout(0) <<"pasintersecte\n";            temp+=1;
1650              if(tt[0][0]==0) {
1651                tt[0][0]=1;
1652                pointsIntersection[0][0]=
1653                  (pairTriangleToIntersectionPoints[pT].listofPoints[cc[0]-1]);
1654              } else {
1655                tt[0][1]=1;
1656                pointsIntersection[0][1]=
1657                  (pairTriangleToIntersectionPoints[pT].listofPoints[cc[0]-1]);
1658              }
1659            }
1660            if(cc[1]!=0) {
1661              temp+=1;
1662              if(tt[1][0]==0) {
1663                tt[1][0]=1;
1664                pointsIntersection[1][0]=
1665                  (pairTriangleToIntersectionPoints[pT].listofPoints[cc[1]-1]);
1666              } else{
1667                tt[1][1]=1;
1668                pointsIntersection[1][1]=
1669                  (pairTriangleToIntersectionPoints[pT].listofPoints[cc[1]-1]);
1670              }
1671            }
1672            if(cc[2]!=0) {
1673              temp+=1;
1674              if(tt[2][0]==0) {
1675                tt[2][0]=1;
1676                pointsIntersection[2][0]=
1677                  (pairTriangleToIntersectionPoints[pT].listofPoints[cc[2]-1]);
1678              } else{
1679                tt[2][1]=1;
1680                pointsIntersection[2][1]=
1681                  (pairTriangleToIntersectionPoints[pT].listofPoints[cc[2]-1]);
1682              }
1683            }
1684    
1685            if(temp<pairTriangleToIntersectionPoints[pT].size()) {
1686              numberin+=1;
1687              size_t r=0;
1688              size_t sum=cc[0]+cc[1]+cc[2];
1689              if(sum>3) {
1690                ffout(0)<<"sum "<<sum<<"\n"<<std::flush;
1691              }
1692              if(sum==0) {
1693                r=0;
1694                for(size_t i=1 ; i<pairTriangleToIntersectionPoints[pT].size() ;++i)
1695                  pointsin[(pairTriangleToIntersectionPoints[pT].listofPoints[i])]=
1696                    (pairTriangleToIntersectionPoints[pT].listofPoints[i]);
1697              } else {
1698                if(sum==2 or sum==5) {
1699                  r=0;
1700                } else {
1701                  if(sum==4 or sum==1) {
1702                    r=1;
1703                  } else {
1704                    if(sum==3) {
1705                      r=2;
1706                    } else {
1707                      ffout(0) << "pbs sum = "<<sum<<" size "
1708                               << pairTriangleToIntersectionPoints[pT].size()<<" \n";
1709                      ffout(0) << K(0) << "\n";
1710                      ffout(0) << K(1) << "\n";
1711                      ffout(0) << K(2) << "\n";
1712                      ffout(0) << cc[0] << " " << cc[1] << " " << cc[2] << "\n";
1713                    }
1714                  }
1715              }              }
1716              }
1717              if(pointsin.find(pairTriangleToIntersectionPoints[pT].listofPoints[r])
1718                 ==pointsin.end()) {
1719                pointsin[(pairTriangleToIntersectionPoints[pT].listofPoints[r])]
1720                  =(pairTriangleToIntersectionPoints[pT].listofPoints[r]);
1721              }
1722          }          }
1723          } else {
1724            // ffout(0) <<"pasintersecte\n";
1725          }
1726      }      }
1727      }
1728  }  }
1729    
1730  TinyVector<3,TinyVector<2,int> > SurfaceMeshGenerator::Internals::  TinyVector<3,TinyVector<2,int> > SurfaceMeshGenerator::Internals::
# Line 1725  __determinateEdgeCut(const size_t& numob Line 1733  __determinateEdgeCut(const size_t& numob
1733                       const Triangle*& currentTriangle,                       const Triangle*& currentTriangle,
1734                       const Cell* cell){                       const Cell* cell){
1735            
1736      TinyVector<3,TinyVector<2,int> > tt;    TinyVector<3,TinyVector<2,int> > tt;
1737      tt=0;    tt=0;
   
     ListOfHexaToTriangle::iterator hexait=listOfHexaToTriangle.begin();  
     MapCellTriangle hexaToTriangle1;  
     if(numobject!=0)  
         hexaToTriangle1=*(*hexait);  
     else {  
         ++hexait;  
         hexaToTriangle1=*(*hexait);  
     }  
       
     MapCellTriangle::iterator it1=hexaToTriangle1.find(cell);  
       
     if(it1!=hexaToTriangle1.end()) {  
         std::list<Triangle*> listTriangle1=(*it1).second;  
         for(std::list<Triangle*>::iterator jt0=listTriangle1.begin() ;  
      jt0!=listTriangle1.end() ; ++jt0) {  
1738    
1739              PairTriangleNew pT;    ListOfHexaToTriangle::iterator hexait=listOfHexaToTriangle.begin();
1740              if(numobject<objectNumber) {    MapCellTriangle hexaToTriangle1;
1741                  pT.first=(currentTriangle);    if(numobject!=0) {
1742                  pT.second=(*jt0);      hexaToTriangle1=*(*hexait);
1743              }    } else {
1744              else {      ++hexait;
1745                  pT.second=(currentTriangle);      hexaToTriangle1=*(*hexait);
1746                  pT.first=(*jt0);    }
             }  
1747    
1748              if(pairTriangleToIntersectionPoints.find(pT)!=    MapCellTriangle::iterator it1=hexaToTriangle1.find(cell);
1749          pairTriangleToIntersectionPoints.end()) {      
1750      if(it1!=hexaToTriangle1.end()) {
1751        std::list<Triangle*> listTriangle1=(*it1).second;
1752        for(std::list<Triangle*>::iterator jt0=listTriangle1.begin() ;
1753            jt0!=listTriangle1.end() ; ++jt0) {
1754    
1755                  TinyVector<3,size_t >cc;        PairTriangleNew pT;
1756                  if(numobject<objectNumber) {        if(numobject<objectNumber) {
1757                      cc=(pairTriangleToIntersectionPoints[pT]).cutTriangle1;          pT.first=(currentTriangle);
1758                  }          pT.second=(*jt0);
1759                  else {        } else {
1760                      cc=(pairTriangleToIntersectionPoints[pT]).cutTriangle2;          pT.second=(currentTriangle);
1761                  }          pT.first=(*jt0);
1762                  if(cc[0]!=0) {        }
1763                      tt[0]=1;  
1764                  }        if(pairTriangleToIntersectionPoints.find(pT)!=
1765                  if(cc[1]!=0) {           pairTriangleToIntersectionPoints.end()) {
1766                      tt[1]=1;  
1767                  }          TinyVector<3,size_t >cc;
1768                  if(cc[2]!=0) {          if(numobject<objectNumber) {
1769                      tt[2]=1;            cc=(pairTriangleToIntersectionPoints[pT]).cutTriangle1;
1770                  }          } else {
1771              }            cc=(pairTriangleToIntersectionPoints[pT]).cutTriangle2;
1772              else {          }
1773                  // ffout(0)<<"la paire de triangle n'existe pas!\n";          if(cc[0]!=0) {
1774              }            tt[0]=1;
1775            }
1776            if(cc[1]!=0) {
1777              tt[1]=1;
1778            }
1779            if(cc[2]!=0) {
1780              tt[2]=1;
1781          }          }
1782          }
1783          else {
1784            // ffout(0)<<"la paire de triangle n'existe pas!\n";
1785          }
1786      }      }
1787      return tt;    }
1788      return tt;
1789  }  }
1790    
1791  TinyVector<3,size_t> SurfaceMeshGenerator::Internals::  TinyVector<3,size_t> SurfaceMeshGenerator::Internals::
1792  __putEdgeRef(const size_t& ncase,  __putEdgeRef(const size_t& ncase,
1793               const Triangle& currentTriangle,               const Triangle& currentTriangle,
1794               TriangleCut& K) {               TriangleCut& K)
1795    {
1796      const TinyVector<3,TinyVector<2,int> >&  tt=K.edgecut;    const TinyVector<3,TinyVector<2,int> >&  tt=K.edgecut;
1797      TinyVector<3,size_t>& ii=K.isInTriangle;    TinyVector<3,size_t>& ii=K.isInTriangle;
1798      size_t numobject=K.numobject;    size_t numobject=K.numobject;
1799      const Vertex* V0=&currentTriangle(0);    const Vertex* V0=&currentTriangle(0);
1800      const Vertex* V1=&currentTriangle(1);    const Vertex* V1=&currentTriangle(1);
1801      const Vertex* V2=&currentTriangle(2);    const Vertex* V2=&currentTriangle(2);
1802    
1803        
1804      if(tt[0][0]+tt[1][0]+tt[2][0]!=0 and ((*edgesRefVertex[numobject])[V0]+
1805                                            (*edgesRefVertex[numobject])[V1]+(*edgesRefVertex[numobject])[V2])==6) {
1806        ffout(0)<<"************************ ---------tr a pbs\n"<<std::flush;
1807        ffout(0)<<tt[0]<<" "<<tt[1]<<" "<<tt[2]<<"\n";
1808        ffout(0)<<*V0<<"\n";
1809        ffout(0)<<*V1<<"\n";
1810        ffout(0)<<*V2<<"\n";
1811      }
1812    
1813      __inout(numobject,tt[0][0],V0,V1);
1814      __inout(numobject,tt[1][0],V0,V2);
1815      __inout(numobject,tt[2][0],V1,V2);
1816            
1817      if(tt[0][0]+tt[1][0]+tt[2][0]!=0 and ((*edgesRefVertex[numobject])[V0]+    if(ncase==0) {
1818                                   (*edgesRefVertex[numobject])[V1]+(*edgesRefVertex[numobject])[V2])==6) {      if(tt[0][0]+tt[1][0]+tt[2][0]==0 and ((*edgesRefVertex[numobject])[V0]+
1819          ffout(0)<<"************************ ---------tr a pbs\n"<<std::flush;                                            (*edgesRefVertex[numobject])[V1]+(*edgesRefVertex[numobject])[V2])>=4) {
1820          ffout(0)<<tt[0]<<" "<<tt[1]<<" "<<tt[2]<<"\n";        (*edgesRefVertex[numobject])[V0]=2;
1821          ffout(0)<<*V0<<"\n";        (*edgesRefVertex[numobject])[V1]=2;
1822          ffout(0)<<*V1<<"\n";        (*edgesRefVertex[numobject])[V2]=2;
1823          ffout(0)<<*V2<<"\n";      }  
1824      }    } else {
1825        if(ncase==1 or ncase==0) {
1826          if((*edgesRefVertex[numobject])[V0]==2 ) {
1827            ii[0]=1;
1828          } else {
1829            ii[0]=0;
1830          }
1831    
1832      __inout(numobject,tt[0][0],V0,V1);        if((*edgesRefVertex[numobject])[V1]==2 ) {
1833      __inout(numobject,tt[1][0],V0,V2);          ii[1]=1;
1834      __inout(numobject,tt[2][0],V1,V2);        } else {
1835                ii[1]=0;
1836      if(ncase==0) {        }
1837          if(tt[0][0]+tt[1][0]+tt[2][0]==0 and ((*edgesRefVertex[numobject])[V0]+  
1838                                (*edgesRefVertex[numobject])[V1]+(*edgesRefVertex[numobject])[V2])>=4) {        if((*edgesRefVertex[numobject])[V2]==2 ) {
1839              (*edgesRefVertex[numobject])[V0]=2;          ii[2]=1;
1840              (*edgesRefVertex[numobject])[V1]=2;        } else {
1841              (*edgesRefVertex[numobject])[V2]=2;          ii[2]=0;
1842          }              }
     }  
     else if(ncase==1 or ncase==0) {  
         if((*edgesRefVertex[numobject])[V0]==2 )  
             ii[0]=1;  
         else  
             ii[0]=0;  
         if((*edgesRefVertex[numobject])[V1]==2 )  
             ii[1]=1;  
         else  
             ii[1]=0;  
         if((*edgesRefVertex[numobject])[V2]==2 )  
             ii[2]=1;  
         else  
             ii[2]=0;  
1843      }      }
1844        }
1845    
1846      return ii;    return ii;
1847  }  }
1848    
1849  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1850  __createListTriangle(const TriangleCut& K,  __createListTriangle(const TriangleCut& K,
1851                       const Cell* cell){                       const Cell* cell)
1852      const TinyVector<3,TinyVector<2,int> >&  tt=K.edgecut;  {
1853      const TinyVector<3,size_t>& ii=K.isInTriangle;    const TinyVector<3,TinyVector<2,int> >&  tt=K.edgecut;
1854      const Vertex& V0=K(0);    const TinyVector<3,size_t>& ii=K.isInTriangle;
1855      const Vertex& V1=K(1);    const Vertex& V0=K(0);
1856      const Vertex& V2=K(2);    const Vertex& V1=K(1);
1857      const Vertex& V2=K(2);
1858        
1859     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*/)) {    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*/)) {
1860         __addTriangle(V0,V1,V2,cell);      __addTriangle(V0,V1,V2,cell);
1861      }    } else {
1862      else if(ii[0]+ii[1]+ii[2]==1) {      if(ii[0]+ii[1]+ii[2]==1) {
1863          if(ii[0]==1) {        if(ii[0]==1) {
1864              __createCase1(0,K,cell);          __createCase1(0,K,cell);
1865          }        }
1866          if(ii[1]==1) {        if(ii[1]==1) {
1867              __createCase1(1,K,cell);          __createCase1(1,K,cell);
1868          }        }
1869          if(ii[2]==1) {        if(ii[2]==1) {
1870              __createCase1(2,K,cell);          __createCase1(2,K,cell);
1871          }        }
1872      }      } else {
1873      else if(ii[0]+ii[1]+ii[2]==2) {        if(ii[0]+ii[1]+ii[2]==2) {
1874          if(ii[0]+ii[1]==2) {          if(ii[0]+ii[1]==2) {
1875              __createCase2(0,K,cell);            __createCase2(0,K,cell);
1876          }          }
1877          if(ii[0]+ii[2]==2) {          if(ii[0]+ii[2]==2) {
1878              __createCase2(1,K,cell);            __createCase2(1,K,cell);
1879          }          }
1880          if(ii[1]+ii[2]==2) {          if(ii[1]+ii[2]==2) {
1881              __createCase2(2,K,cell);            __createCase2(2,K,cell);
1882          }          }
1883          }
1884      }      }
1885      }
1886  }  }
1887    
1888  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1889  __createNotIn(const TriangleCut& K,const Cell* cell) {  __createNotIn(const TriangleCut& K,const Cell* cell)
1890      size_t place0=0;  {
1891      if(K.edgecut[1][0]==1 and K.edgecut[1][1]==1) {    size_t place0=0;
1892          place0=1;    if(K.edgecut[1][0]==1 and K.edgecut[1][1]==1) {
1893      }      place0=1;
1894      else if(K.edgecut[2][0]==1 and K.edgecut[2][1]==1) {    } else if(K.edgecut[2][0]==1 and K.edgecut[2][1]==1) {
1895          place0=2;      place0=2;
1896      }    }
1897      const Vertex* Pend=K.pointsIntersection[place0][1];  
1898      TinyVector<10,const Vertex*> P;    const Vertex* Pend=K.pointsIntersection[place0][1];
1899      TinyVector<10,const Triangle*> T;    TinyVector<10,const Vertex*> P;
1900      P[0]=K.pointsIntersection[place0][0];    TinyVector<10,const Triangle*> T;
1901      if(vertexInTriangles.find(P[0])!=vertexInTriangles.end()) {    P[0]=K.pointsIntersection[place0][0];
1902          T[0]=(*vertexInTriangles.find(P[0])).second[0];  
1903      }    if(vertexInTriangles.find(P[0])!=vertexInTriangles.end()) {
1904      else {      T[0]=(*vertexInTriangles.find(P[0])).second[0];
1905          ffout(0)<<"pbs 0 existe pas\n";    } else {
1906      }      ffout(0)<<"pbs 0 existe pas\n";
1907      size_t b=0;    }
1908      while(P[b]!=Pend and b<9) {  
1909          if(triangleWithVertex.find(T[b])!=triangleWithVertex.end()) {    size_t b=0;
1910              TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[b])).second;  
1911              if(VV[0]==P[b]) {    while(P[b]!=Pend and b<9) {
1912                  if(Pend!=VV[1])      if(triangleWithVertex.find(T[b])!=triangleWithVertex.end()) {
1913                      __addTriangle(*Pend,*P[b],*VV[1],cell);        TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[b])).second;
1914                  P[b+1]=VV[1];        if(VV[0]==P[b]) {
1915              }          if(Pend!=VV[1])
1916              else {            __addTriangle(*Pend,*P[b],*VV[1],cell);
1917                  if(Pend!=VV[0])          P[b+1]=VV[1];
1918                      __addTriangle(*Pend,*P[b],*VV[0],cell);        } else {
1919                  P[b+1]=VV[0];          if(Pend!=VV[0])
1920              }            __addTriangle(*Pend,*P[b],*VV[0],cell);
1921          }          P[b+1]=VV[0];
1922          else        }
1923              ffout(0)<<"pbs 1 existe pas\n";      } else {
1924          TinyVector<2,const Triangle*> T2;        ffout(0)<<"pbs 1 existe pas\n";
1925          if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {      }
1926              T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];      TinyVector<2,const Triangle*> T2;
1927              T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];      if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {
1928          }        T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];
1929          if(T2[0]==T[b]) {        T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];
1930              if(T2[1]==T[b]) {      }
1931                  //point edge      if(T2[0]==T[b]) {
1932              }        if(T2[1]==T[b]) {
1933              else          //point edge
1934                  T[b+1]=T2[1];        } else {
1935          }          T[b+1]=T2[1];
1936          else {        }
1937              T[b+1]=T2[0];      } else {
1938          }        T[b+1]=T2[0];
         ++b;  
1939      }      }
1940        ++b;
1941      }
1942  }  }
1943    
1944  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1945  __createNotIn2(const TriangleCut& K,const Cell* cell) {  __createNotIn2(const TriangleCut& K,const Cell* cell)
1946      size_t place0=0;  {
1947      if(K.edgecut[1][0]==1 and K.edgecut[1][1]==1) {    size_t place0=0;
1948          place0=1;    if(K.edgecut[1][0]==1 and K.edgecut[1][1]==1) {
1949      }      place0=1;
1950      else if(K.edgecut[2][0]==1 and K.edgecut[2][1]==1) {    } else
1951          place0=2;      if(K.edgecut[2][0]==1 and K.edgecut[2][1]==1) {
1952      }        place0=2;
1953      //ffout(0)<<"place = "<<place0<<"\n";      }
1954      const Vertex* Pend=K.pointsIntersection[place0][1];  
1955      TinyVector<10,const Vertex*> P;    //ffout(0)<<"place = "<<place0<<"\n";
1956      TinyVector<10,const Triangle*> T;    const Vertex* Pend=K.pointsIntersection[place0][1];
1957      P[0]=K.pointsIntersection[place0][0];    TinyVector<10,const Vertex*> P;
1958      if(vertexInTriangles.find(P[0])!=vertexInTriangles.end()) {    TinyVector<10,const Triangle*> T;
1959          T[0]=(*vertexInTriangles.find(P[0])).second[0];    P[0]=K.pointsIntersection[place0][0];
1960      }    if(vertexInTriangles.find(P[0])!=vertexInTriangles.end()) {
1961      else {      T[0]=(*vertexInTriangles.find(P[0])).second[0];
1962          ffout(0)<<"pbs 0 existe pas\n";    } else {
1963      }      ffout(0)<<"pbs 0 existe pas\n";
1964      size_t b=0;    }
1965      while(P[b]!=Pend and b<9) {  
1966          //ffout(0)<<"------------------- b "<<b<<"\n";    size_t b=0;
1967          if(triangleWithVertex.find(T[b])!=triangleWithVertex.end()) {    while(P[b]!=Pend and b<9) {
1968              TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[b])).second;      //ffout(0)<<"------------------- b "<<b<<"\n";
1969              if(VV[0]==P[b]) {      if(triangleWithVertex.find(T[b])!=triangleWithVertex.end()) {
1970                  if(Pend!=VV[1]) {        TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[b])).second;
1971                      __addTriangle(*Pend,*P[b],*VV[1],cell);        if(VV[0]==P[b]) {
1972                    //  ffout(0)<<"1 "<<*Pend<<"\n";          if(Pend!=VV[1]) {
1973                    //  ffout(0)<<*P[b]<<"\n";            __addTriangle(*Pend,*P[b],*VV[1],cell);
1974                    //  ffout(0)<<*VV[1]<<"\n";            //  ffout(0)<<"1 "<<*Pend<<"\n";
1975                  }            //  ffout(0)<<*P[b]<<"\n";
1976                  P[b+1]=VV[1];            //  ffout(0)<<*VV[1]<<"\n";
1977              }          }
1978              else {          P[b+1]=VV[1];
1979                  if(Pend!=VV[0]) {        } else {
1980                      __addTriangle(*Pend,*P[b],*VV[0],cell);          if(Pend!=VV[0]) {
1981                      //ffout(0)<<"2 "<<*Pend<<"\n";            __addTriangle(*Pend,*P[b],*VV[0],cell);
1982                      //ffout(0)<<*P[b]<<"\n";            //ffout(0)<<"2 "<<*Pend<<"\n";
1983                      //ffout(0)<<*VV[0]<<"\n";            //ffout(0)<<*P[b]<<"\n";
1984                  }            //ffout(0)<<*VV[0]<<"\n";
                 P[b+1]=VV[0];  
             }  
         }  
         else {  
             ffout(0)<<"pbs 1 existe pas\n";  
         }  
         TinyVector<2,const Triangle*> T2;  
         if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {  
             T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];  
             T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];  
         }  
         else {  
             ffout(0)<<"points exisste pas \n";  
1985          }          }
1986            P[b+1]=VV[0];
1987          }
1988        } else {
1989          ffout(0)<<"pbs 1 existe pas\n";
1990        }
1991    
1992        TinyVector<2,const Triangle*> T2;
1993        if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {
1994          T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];
1995          T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];
1996        } else {
1997          ffout(0)<<"points exisste pas \n";
1998        }
1999        if(K.edgecut[0][0]==1 and K.edgecut[0][1]==1) {
2000          if(P[b+1]==K.pointsIntersection[0][1] or P[b+1]==K.pointsIntersection[0][0]) {
2001            T2[0]=T[b];
2002            T2[1]=T[b];
2003          }
2004        } else
2005          if(K.edgecut[2][0]==1 and K.edgecut[2][1]==1) {
2006            if(P[b+1]==K.pointsIntersection[2][1] or P[b+1]==K.pointsIntersection[2][0]) {
2007              T2[0]=T[b];
2008              T2[1]=T[b];
2009            }
2010          }
2011        if(T2[0]==T[b]) {
2012          if(T2[1]==T[b]) {
2013            //point edge
2014            //ffout(0)<<"here\n";
2015          if(K.edgecut[0][0]==1 and K.edgecut[0][1]==1) {          if(K.edgecut[0][0]==1 and K.edgecut[0][1]==1) {
2016              if(P[b+1]==K.pointsIntersection[0][1] or P[b+1]==K.pointsIntersection[0][0]) {            //ffout(0)<<"here2\n";
2017                  T2[0]=T[b];            if(P[b+1]==K.pointsIntersection[0][1]) {
2018                  T2[1]=T[b];              __addTriangle(*Pend,*P[b+1],*K.pointsIntersection[0][0],cell);
2019              }              P[b+1]=K.pointsIntersection[0][0];                  
2020          }              if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {
2021          else if(K.edgecut[2][0]==1 and K.edgecut[2][1]==1) {                T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];
2022              if(P[b+1]==K.pointsIntersection[2][1] or P[b+1]==K.pointsIntersection[2][0]) {                T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];
                 T2[0]=T[b];  
                 T2[1]=T[b];  
             }  
         }  
         if(T2[0]==T[b]) {  
             if(T2[1]==T[b]) {  
                 //point edge  
                 //ffout(0)<<"here\n";  
                 if(K.edgecut[0][0]==1 and K.edgecut[0][1]==1) {  
                     //ffout(0)<<"here2\n";  
                     if(P[b+1]==K.pointsIntersection[0][1]) {  
                         __addTriangle(*Pend,*P[b+1],*K.pointsIntersection[0][0],cell);  
                         P[b+1]=K.pointsIntersection[0][0];                        
                         if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {  
                             T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];  
                             T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];  
                         }  
                         T[b+1]=T2[0];  
                     }  
                     else if(P[b+1]==K.pointsIntersection[0][0]){  
                         __addTriangle(*Pend,*P[b+1],*K.pointsIntersection[0][1],cell);  
                         P[b+1]=K.pointsIntersection[0][1];  
                         if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {  
                             T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];  
                             T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];  
                         }  
                         else {  
                             ffout(0)<<"le tr e pas\n";  
                         }  
                         T[b+1]=T2[0];  
                     }  
                     else if(P[b+1]!=Pend) {  
                         ffout(0)<<"point pas sur edge trou???? \n";  
                         ffout(0)<<*P[b+1]<<"\n";  
                         ffout(0)<<*Pend<<"\n";  
                         P[b+1]=Pend;  
                     }  
                 }  
                 else if(K.edgecut[2][0]==1 and K.edgecut[2][1]==1) {  
                     ffout(0)<<"here3\n";  
                     if(P[b+1]==K.pointsIntersection[2][1]) {  
                         __addTriangle(*Pend,*P[b+1],*K.pointsIntersection[2][0],cell);  
                         P[b+1]=K.pointsIntersection[2][0];  
                         if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {  
                             T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];  
                             T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];  
                         }  
                         T[b+1]=T2[0];  
                           
                     }  
                     else if(P[b+1]==K.pointsIntersection[2][0]){  
                         __addTriangle(*Pend,*P[b+1],*K.pointsIntersection[2][1],cell);  
                         P[b+1]=K.pointsIntersection[2][1];  
                         if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {  
                             T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];  
                             T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];  
                         }  
                         T[b+1]=T2[0];  
                     }  
                 }  
                 else  
                     ffout(0)<<"pas trouve l'egde et le point\n";  
2023              }              }
             else  
                 T[b+1]=T2[1];  
         }  
         else {  
2024              T[b+1]=T2[0];              T[b+1]=T2[0];
2025              } else
2026                if(P[b+1]==K.pointsIntersection[0][0]){
2027                  __addTriangle(*Pend,*P[b+1],*K.pointsIntersection[0][1],cell);
2028                  P[b+1]=K.pointsIntersection[0][1];
2029                  if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {
2030                    T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];
2031                    T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];
2032                  } else {
2033                    ffout(0)<<"le tr e pas\n";
2034                  }
2035                  T[b+1]=T2[0];
2036                } else
2037                  if(P[b+1]!=Pend) {
2038                    ffout(0)<<"point pas sur edge trou???? \n";
2039                    ffout(0)<<*P[b+1]<<"\n";
2040                    ffout(0)<<*Pend<<"\n";
2041                    P[b+1]=Pend;
2042                  }
2043            } else {
2044              if(K.edgecut[2][0]==1 and K.edgecut[2][1]==1) {
2045                ffout(0)<<"here3\n";
2046                if(P[b+1]==K.pointsIntersection[2][1]) {
2047                  __addTriangle(*Pend,*P[b+1],*K.pointsIntersection[2][0],cell);
2048                  P[b+1]=K.pointsIntersection[2][0];
2049                  if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {
2050                    T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];
2051                    T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];
2052                  }
2053                  T[b+1]=T2[0];
2054                } else {
2055                  if(P[b+1]==K.pointsIntersection[2][0]){
2056                  __addTriangle(*Pend,*P[b+1],*K.pointsIntersection[2][1],cell);
2057                  P[b+1]=K.pointsIntersection[2][1];
2058                  if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {
2059                    T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];
2060                    T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];
2061                  }
2062                  T[b+1]=T2[0];
2063                  }
2064                }
2065              } else {
2066                ffout(0)<<"pas trouve l'egde et le point\n";
2067              }
2068          }          }
2069          ++b;        } else {
2070            T[b+1]=T2[1];
2071          }
2072        } else {
2073          T[b+1]=T2[0];
2074      }      }
2075        ++b;
2076      }
2077  }  }
2078    
2079  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
2080  __createByEdges(const TriangleCut& K,const Cell* cell) {  __createByEdges(const TriangleCut& K,const Cell* cell)
2081    {
2082      size_t in0=0,place0=0,place1=1,place2=2;
2083    
2084      size_t in0=0,place0=0,place1=1,place2=2;    if(K.isInTriangle[1]==1) {
2085      if(K.isInTriangle[1]==1) {      in0=1;
2086          in0=1;      place1=2;
2087          place1=2;      place2=1;
2088          place2=1;    } else {
2089      }      if(K.isInTriangle[2]==1) {
2090      else if(K.isInTriangle[2]==1) {        in0=2;
2091          in0=2;        place0=1;
2092          place0=1;        place1=2;
2093          place1=2;        place2=0;
2094          place2=0;      }
2095      }    }
2096      const Vertex* Pend=K.pointsIntersection[place1][0];  
2097      TinyVector<10,const Vertex*> P;    const Vertex* Pend=K.pointsIntersection[place1][0];
2098      TinyVector<10,const Triangle*> T;    TinyVector<10,const Vertex*> P;
2099      P[0]=K.pointsIntersection[place0][0];    TinyVector<10,const Triangle*> T;
2100      if(vertexInTriangles.find(P[0])!=vertexInTriangles.end()) {    P[0]=K.pointsIntersection[place0][0];
2101          T[0]=(*vertexInTriangles.find(P[0])).second[0];    if(vertexInTriangles.find(P[0])!=vertexInTriangles.end()) {
2102      }      T[0]=(*vertexInTriangles.find(P[0])).second[0];
2103      else {    } else {
2104          ffout(0)<<"pbs 0 existe pas\n";      ffout(0)<<"pbs 0 existe pas\n";
2105      }    }
2106      size_t b=0;  
2107      while(P[b]!=Pend and b<9) {    size_t b=0;
2108          if(triangleWithVertex.find(T[b])!=triangleWithVertex.end()) {    while(P[b]!=Pend and b<9) {
2109              TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[b])).second;      if(triangleWithVertex.find(T[b])!=triangleWithVertex.end()) {
2110              if(VV[0]==P[b]) {        TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[b])).second;
2111                  if(VV[1]==P[b]) {        if(VV[0]==P[b]) {
2112                      ffout(0)<<"egal\n";          if(VV[1]==P[b]) {
2113                      if(K.pointsIntersection[place2][0]==P[b]) {            ffout(0)<<"egal\n";
2114                          P[b+1]=K.pointsIntersection[place2][1];            if(K.pointsIntersection[place2][0]==P[b]) {
2115                          __addTriangle(K(in0),*P[b],*P[b+1],cell);              P[b+1]=K.pointsIntersection[place2][1];
2116                      }              __addTriangle(K(in0),*P[b],*P[b+1],cell);
2117                      else if(K.pointsIntersection[place2][1]==P[b]){            } else {
2118                          P[b+1]=K.pointsIntersection[place2][0];              if(K.pointsIntersection[place2][1]==P[b]){
2119                          __addTriangle(K(in0),*P[b],*P[b+1],cell);                P[b+1]=K.pointsIntersection[place2][0];
2120                      }                __addTriangle(K(in0),*P[b],*P[b+1],cell);
2121                      else {              } else {
2122                          ffout(0)<<"pbs pbs\n";                ffout(0)<<"pbs pbs\n";
2123                          ffout(0)<<*P[b]<<"\n";                ffout(0)<<*P[b]<<"\n";
2124                          std::exit(0);                std::exit(0);
2125                      }              }
2126                  }            }
2127                  else {          } else {
2128                  __addTriangle(K(in0),*P[b],*VV[1],cell);            __addTriangle(K(in0),*P[b],*VV[1],cell);
2129                  P[b+1]=VV[1];            P[b+1]=VV[1];
2130                  }          }
2131              }        } else {
2132              else {          __addTriangle(K(in0),*P[b],*VV[0],cell);
2133                  __addTriangle(K(in0),*P[b],*VV[0],cell);          P[b+1]=VV[0];
2134                  P[b+1]=VV[0];        }
2135              }      } else {
2136          }        ffout(0)<<"pbs 1 existe pas\n";
2137          else      }
2138              ffout(0)<<"pbs 1 existe pas\n";      TinyVector<2,const Triangle*> T2;
2139          TinyVector<2,const Triangle*> T2;  
2140          if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {      if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {
2141              T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];        T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];
2142              T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];        T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];
2143          }      } else {
2144          else {        ffout(0)<<"pbs tr existe pas\n";
2145              ffout(0)<<"pbs tr existe pas\n";      }
2146          }  
2147          if(T2[0]==T[b] and T2[1]==T[b]) {      if(T2[0]==T[b] and T2[1]==T[b]) {
2148             // ffout(0)<<"here0\n";        // ffout(0)<<"here0\n";
2149              //point edge        //point edge
2150              if(P[b+1]!=Pend) {        if(P[b+1]!=Pend) {
2151                  if(K.pointsIntersection[place2][0]==P[b+1]) {          if(K.pointsIntersection[place2][0]==P[b+1]) {
2152                      const Vertex* A=K.pointsIntersection[place2][1];            const Vertex* A=K.pointsIntersection[place2][1];
2153                      __addTriangle(K(in0),*P[b+1],*A,cell);            __addTriangle(K(in0),*P[b+1],*A,cell);
2154                      if(vertexInTriangles.find(A)!=vertexInTriangles.end()) {            if(vertexInTriangles.find(A)!=vertexInTriangles.end()) {
2155                          T[b+1]=(*vertexInTriangles.find(A)).second[0];              T[b+1]=(*vertexInTriangles.find(A)).second[0];
2156                          P[b+1]=A;              P[b+1]=A;
2157                      }            }
2158                  }          } else {
2159                  else {            const Vertex* A=K.pointsIntersection[place2][0];
2160                      const Vertex* A=K.pointsIntersection[place2][0];            __addTriangle(K(in0),*P[b+1],*A,cell);
2161                      __addTriangle(K(in0),*P[b+1],*A,cell);            if(vertexInTriangles.find(A)!=vertexInTriangles.end()) {
2162                      if(vertexInTriangles.find(A)!=vertexInTriangles.end()) {              T[b+1]=(*vertexInTriangles.find(A)).second[0];
2163                          T[b+1]=(*vertexInTriangles.find(A)).second[0];              P[b+1]=A;
2164                          P[b+1]=A;            }
                     }  
                 }  
             }  
         }  
         else if(T2[0]==T[b]) {  
             T[b+1]=T2[1];  
         }  
         else {  
             T[b+1]=T2[0];  
2165          }          }
2166          ++b;        }
2167        } else {
2168          if(T2[0]==T[b]) {
2169            T[b+1]=T2[1];
2170          } else {
2171            T[b+1]=T2[0];
2172          }
2173      }      }
2174            ++b;
2175      }
2176  }  }
2177    
2178    
2179  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
2180  __createTrianglesIntersection(const size_t& numobject,  __createTrianglesIntersection(const size_t& numobject,
2181                                const size_t& objectNumber) {                                const size_t& objectNumber)
2182    {
2183      ffout(0)<<"numero "<<numobject<<" "<<objectNumber<<"\n"<<std::flush;    ffout(0)<<"numero "<<numobject<<" "<<objectNumber<<"\n"<<std::flush;
2184    
2185      ListOfHexaToTriangle::iterator hexait=listOfHexaToTriangle.begin();    ListOfHexaToTriangle::iterator hexait=listOfHexaToTriangle.begin();
2186      MapCellTriangle hexaToTriangle0;    MapCellTriangle hexaToTriangle0;
2187      if(numobject==0)    if(numobject==0) {
2188          hexaToTriangle0=*(*hexait);      hexaToTriangle0=*(*hexait);
2189      else {    } else {
2190          ++hexait;      ++hexait;
2191          hexaToTriangle0=*(*hexait);      hexaToTriangle0=*(*hexait);
2192      }    }
2193  //Initialize edgesRefVertex  //Initialize edgesRefVertex
2194      size_t taille=0;    size_t taille=0;
2195      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) {
2196          MapCellTriangle::iterator it0=hexaToTriangle0.find(*i);      MapCellTriangle::iterator it0=hexaToTriangle0.find(*i);
2197          if(it0!=hexaToTriangle0.end()) {      if(it0!=hexaToTriangle0.end()) {
2198              std::list<Triangle*> listTriangle0=(*it0).second;        std::list<Triangle*> listTriangle0=(*it0).second;
2199              taille +=listTriangle0.size();        taille +=listTriangle0.size();
2200              for(std::list<Triangle*>::iterator jt0=listTriangle0.begin() ;        for(std::list<Triangle*>::iterator jt0=listTriangle0.begin() ;
2201           jt0!=listTriangle0.end() ; ++jt0) {            jt0!=listTriangle0.end() ; ++jt0) {
2202                  const Triangle& T = *(*jt0);          const Triangle& T = *(*jt0);
2203                            
2204                  (*edgesRefVertex[numobject])[&(T(0))]=0;          (*edgesRefVertex[numobject])[&(T(0))]=0;
2205                  (*edgesRefVertex[numobject])[&(T(1))]=0;          (*edgesRefVertex[numobject])[&(T(1))]=0;
2206                  (*edgesRefVertex[numobject])[&(T(2))]=0;          (*edgesRefVertex[numobject])[&(T(2))]=0;
2207              }        }
         }  
2208      }      }
2209      }
2210    
2211      ListOfTriangleList::iterator triangleit=listOfTriangleList.begin();    ListOfTriangleList::iterator triangleit=listOfTriangleList.begin();
2212      if(numobject!=0) {    if(numobject!=0) {
2213          ++triangleit;      ++triangleit;
2214      }    }
2215      for(size_t size=0 ; size<(*(*triangleit)).size() ; ++size)  {    for(size_t size=0 ; size<(*(*triangleit)).size() ; ++size)  {
2216          const Triangle& T=(*(*triangleit))[size];      const Triangle& T=(*(*triangleit))[size];
2217          (*edgesRefVertex[numobject])[&T(0)]=2;      (*edgesRefVertex[numobject])[&T(0)]=2;
2218          (*edgesRefVertex[numobject])[&T(1)]=2;      (*edgesRefVertex[numobject])[&T(1)]=2;
2219          (*edgesRefVertex[numobject])[&T(2)]=2;      (*edgesRefVertex[numobject])[&T(2)]=2;
2220      }    }
2221            
2222     size_t nbinters=listVertexIntersectionNew.size();    size_t nbinters=listVertexIntersectionNew.size();
2223      ffout(0)<<"size listVertexIntersectionNew "<<nbinters<<"\n"<<std::flush;    ffout(0)<<"size listVertexIntersectionNew "<<nbinters<<"\n"<<std::flush;
   
     size_t dd_max=4;//4  
     ffout(0)<<"size liste intersection "<<taille<<"\n"<<std::flush;  
     for(size_t dd=0 ; dd<dd_max ; ++dd) {  
         ffout(0)<<"num passage "<<dd<<"\n"<<std::flush;  
         for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)  
         {  
             MapCellTriangle::iterator it0=hexaToTriangle0.find(*i);  
             if(it0!=hexaToTriangle0.end()) {  
                 std::list<Triangle*> listTriangle0=(*it0).second;  
                 for(std::list<Triangle*>::iterator jt0=listTriangle0.begin()  
       ; jt0!=listTriangle0.end() ; ++jt0) {  
                     const Triangle* T=(*jt0);  
                     const Vertex V0=(*T)(0);  
                     const Vertex V1=(*T)(1);  
                     const Vertex V2=(*T)(2);  
                     TriangleCut K(V0,V1,V2,numobject,objectNumber);  
   
                     K.edgecut =__determinateEdgeCut(numobject,objectNumber,T,*i);  
                     K.isInTriangle =__putEdgeRef(0,(*T),K);  
                 }  
             }  
         }  
     }    
2224    
2225      size_t err=0,err2=0;    size_t dd_max=4;//4
2226      ffout(0)<<"size liste intersection "<<taille<<"\n"<<std::flush;
2227      for(size_t dd=0 ; dd<dd_max ; ++dd) {
2228        ffout(0)<<"num passage "<<dd<<"\n"<<std::flush;
2229      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)
2230      {        {
2231          MapCellTriangle::iterator it0=hexaToTriangle0.find(*i);          MapCellTriangle::iterator it0=hexaToTriangle0.find(*i);
                   
2232          if(it0!=hexaToTriangle0.end()) {          if(it0!=hexaToTriangle0.end()) {
2233              std::list<Triangle*> listTriangle0=(*it0).second;            std::list<Triangle*> listTriangle0=(*it0).second;
2234              for(std::list<Triangle*>::iterator jt0=listTriangle0.begin() ;            for(std::list<Triangle*>::iterator jt0=listTriangle0.begin();
2235           jt0!=listTriangle0.end() ; ++jt0) {                jt0!=listTriangle0.end() ; ++jt0) {
2236                  const Triangle* T=(*jt0);              const Triangle* T=(*jt0);
2237                  const Vertex& V0=(*T)(0);              const Vertex V0=(*T)(0);
2238                  const Vertex& V1=(*T)(1);              const Vertex V1=(*T)(1);
2239                  const Vertex& V2=(*T)(2);              const Vertex V2=(*T)(2);
2240                  TriangleCut K(V0,V1,V2,numobject,objectNumber);              TriangleCut K(V0,V1,V2,numobject,objectNumber);
2241                  triangleWithVertex.clear();  
2242                  vertexInTriangles.clear();              K.edgecut =__determinateEdgeCut(numobject,objectNumber,T,*i);
2243                K.isInTriangle =__putEdgeRef(0,(*T),K);
2244              }
2245            }
2246          }
2247      }
2248    
2249                  __createLocalListIntersection(objectNumber,T,K,*i);    size_t err=0,err2=0;
2250                  K.isInTriangle =__putEdgeRef(1,*T,K);    for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin();
2251          i!=toTreatHexahedra.end() ; ++i) {
2252        MapCellTriangle::iterator it0=hexaToTriangle0.find(*i);
2253                                    
2254                  if(K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==0 ) {      if(it0!=hexaToTriangle0.end()) {
2255                      ++err;        std::list<Triangle*> listTriangle0=(*it0).second;
2256                  }        for(std::list<Triangle*>::iterator jt0=listTriangle0.begin() ;
2257                  if(K.edgecut[0]+K.edgecut[1]+K.edgecut[2]==0 ) {            jt0!=listTriangle0.end() ; ++jt0) {
2258                      ++err2;          const Triangle* T=(*jt0);
2259                  }          const Vertex& V0=(*T)(0);
2260                  /*TinyVector<3> A(1.000000, -1.800000, -0.800000);          const Vertex& V1=(*T)(1);
2261                  TinyVector<3> B(1.000000 ,-1.800000 ,-1.200000);          const Vertex& V2=(*T)(2);
2262                  TinyVector<3> C(1.000000 ,-1.600000 ,-1.000000 );          TriangleCut K(V0,V1,V2,numobject,objectNumber);
2263                  real_t ee=1e-5;          triangleWithVertex.clear();
2264                  if((Norm(V0.Vector3()-A)<=ee or Norm(V1.Vector3()-A)<=ee or Norm(V2.Vector3()-A)<=ee)          vertexInTriangles.clear();
2265       and (Norm(V0.Vector3()-B)<=ee or Norm(V1.Vector3()-B)<=ee or Norm(V2.Vector3()-B)<=ee)  
2266       and (Norm(V0.Vector3()-C)<=ee or Norm(V1.Vector3()-C)<=ee or Norm(V2.Vector3()-C)<=ee)) {          __createLocalListIntersection(objectNumber,T,K,*i);
2267                      ffout(0)<<"*****************\n";          K.isInTriangle =__putEdgeRef(1,*T,K);
2268                      ffout(0)<<V0<<"\n";  
2269                      ffout(0)<<V1<<"\n";          if(K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==0 ) {
2270                      ffout(0)<<V2<<"\n";            ++err;
2271                      ffout(0)<<K.edgecut<<"\n";          }
2272                      ffout(0)<<K.isInTriangle<<"\n";          if(K.edgecut[0]+K.edgecut[1]+K.edgecut[2]==0 ) {
2273                      ffout(0)<<vertexInTriangles.size()<<"\n";            ++err2;
2274                      for(std::map<const Vertex*,TinyVector<2,const Triangle*> >::iterator          }
2275            /*TinyVector<3> A(1.000000, -1.800000, -0.800000);
2276              TinyVector<3> B(1.000000 ,-1.800000 ,-1.200000);
2277              TinyVector<3> C(1.000000 ,-1.600000 ,-1.000000 );
2278              real_t ee=1e-5;
2279              if((Norm(V0.Vector3()-A)<=ee or Norm(V1.Vector3()-A)<=ee or Norm(V2.Vector3()-A)<=ee)
2280              and (Norm(V0.Vector3()-B)<=ee or Norm(V1.Vector3()-B)<=ee or Norm(V2.Vector3()-B)<=ee)
2281              and (Norm(V0.Vector3()-C)<=ee or Norm(V1.Vector3()-C)<=ee or Norm(V2.Vector3()-C)<=ee)) {
2282              ffout(0)<<"*****************\n";
2283              ffout(0)<<V0<<"\n";
2284              ffout(0)<<V1<<"\n";
2285              ffout(0)<<V2<<"\n";
2286              ffout(0)<<K.edgecut<<"\n";
2287              ffout(0)<<K.isInTriangle<<"\n";
2288              ffout(0)<<vertexInTriangles.size()<<"\n";
2289              for(std::map<const Vertex*,TinyVector<2,const Triangle*> >::iterator
2290            it=vertexInTriangles.begin() ; it!=vertexInTriangles.end() ; ++it) {            it=vertexInTriangles.begin() ; it!=vertexInTriangles.end() ; ++it) {
2291                          ffout(0)<<*(*it).first<<"\n";            ffout(0)<<*(*it).first<<"\n";
2292                      }            }
2293                  }*/            }*/
2294                  if(K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==0          if(K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==0
2295       and K.edgecut[0]+K.edgecut[1]+K.edgecut[2]==3 and vertexInTriangles.size()<3){             and K.edgecut[0]+K.edgecut[1]+K.edgecut[2]==3 and vertexInTriangles.size()<3) {
2296                      ffout(0)<<"tous les points d'inter ne sont pas dans la liste\n";            ffout(0)<<"tous les points d'inter ne sont pas dans la liste\n";
2297                      //ffout(0)<<V0<<"\n";            //ffout(0)<<V0<<"\n";
2298                      //ffout(0)<<V1<<"\n";            //ffout(0)<<V1<<"\n";
2299                      //ffout(0)<<V2<<"\n";            //ffout(0)<<V2<<"\n";
2300                      //ffout(0)<<K.edgecut<<"\n";            //ffout(0)<<K.edgecut<<"\n";
2301                      //ffout(0)<<K.isInTriangle<<"\n";            //ffout(0)<<K.isInTriangle<<"\n";
2302                      __casePinterNotInList(K,*i);            __casePinterNotInList(K,*i);
2303                                }
2304                  }  
2305                  if(K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==0          if(K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==0
2306       and K.edgecut[0]+K.edgecut[1]+K.edgecut[2]!=0) {             and K.edgecut[0]+K.edgecut[1]+K.edgecut[2]!=0) {
2307                      if(vertexInTriangles.size()==3) {            if(vertexInTriangles.size()==3) {
2308                          TinyVector<3,const Vertex*> V;              TinyVector<3,const Vertex*> V;
2309                          size_t num=0;              size_t num=0;
2310                          for(std::map<const Vertex*,TinyVector<2,const Triangle*> >::iterator              for(std::map<const Vertex*,TinyVector<2,const Triangle*> >::iterator
2311         it=vertexInTriangles.begin() ; it!=vertexInTriangles.end() ; ++it , ++num) {                    it=vertexInTriangles.begin() ; it!=vertexInTriangles.end() ; ++it , ++num) {
2312                              V[num]=(*it).first;                V[num]=(*it).first;
2313                          }              }
2314                          if(V[0]!=V[1] and V[1]!=V[2] and V[2]!=V[0]) {              if(V[0]!=V[1] and V[1]!=V[2] and V[2]!=V[0]) {
2315                              __addTriangle(*V[0],*V[1],*V[2],*i);                __addTriangle(*V[0],*V[1],*V[2],*i);
2316                              ffout(0)<<"on comble le trou\n"<<std::flush;                ffout(0)<<"on comble le trou\n"<<std::flush;
2317                              //ffout(0)<<V0<<"\n";                //ffout(0)<<V0<<"\n";
2318                              //ffout(0)<<V1<<"\n";                //ffout(0)<<V1<<"\n";
2319                              //ffout(0)<<V2<<"\n";                //ffout(0)<<V2<<"\n";
2320                          }              } else {
2321                          else {                ffout(0)<<"----- deux inter =\n";
2322                              ffout(0)<<"----- deux inter =\n";                ffout(0)<<V[0]<<"\n";
2323                              ffout(0)<<V[0]<<"\n";                ffout(0)<<V[1]<<"\n";
2324                              ffout(0)<<V[1]<<"\n";                ffout(0)<<V[2]<<"\n";
2325                              ffout(0)<<V[2]<<"\n";              }
2326                          }            } else {
2327                      }              if(vertexInTriangles.size()>3 and
2328                      else if(vertexInTriangles.size()>3 and                 K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==1){
2329                K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==1){                //one edge cut but twice
2330                          //one edge cut but twice                //ffout(0)<<"createNotIn\n";
2331                          //ffout(0)<<"createNotIn\n";                //POUR STEPH : createNotIn peut generer des tr plats
2332                          //POUR STEPH : createNotIn peut generer des tr plats                __createNotIn(K,*i);
2333                          __createNotIn(K,*i);              } else {
2334                      }                if(vertexInTriangles.size()>3 and
2335                      else if(vertexInTriangles.size()>3 and                   K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==2) {
2336                K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==2){                //two edge cut but twice
2337                          //two edge cut but twice                //ffout(0)<<"createNotIn2\n";
2338                          //ffout(0)<<"createNotIn2\n";                  __createNotIn2(K,*i);
2339                  } else {
2340                    if(K.edgecut[0]!=0) {
2341                      if(K.pointsIntersection[0][0]==&V0
2342                         or K.pointsIntersection[0][0]==&V1
2343                         or K.pointsIntersection[0][0]==&V2 or vertexInTriangles.size()==1) {
2344                        // ffout(0)<<"**********normal\n";
2345                      } else {
2346                        if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==3
2347                           and K.edgecut[0][1]+K.edgecut[1][1]+K.edgecut[2][1]==3) {
2348                          if(K.pointsIntersection[0][1]==K.pointsIntersection[0][0]) {
2349                            K.edgecut[0][0]=0;
2350                            K.edgecut[0][1]=0;
2351                          }
2352                          if(K.pointsIntersection[1][1]==K.pointsIntersection[1][0]) {
2353                            K.edgecut[1][0]=0;
2354                            K.edgecut[1][1]=0;
2355                          }
2356                          if(K.pointsIntersection[2][1]==K.pointsIntersection[2][0]) {
2357                            K.edgecut[2][0]=0;
2358                            K.edgecut[2][1]=0;
2359                          }
2360                          if(vertexInTriangles.size()>3 and
2361                             K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==2){
2362                            ffout(0)<<"hhhhhhhh\n";
2363                          __createNotIn2(K,*i);                          __createNotIn2(K,*i);
2364                      }                        } else {
2365                      else if(K.edgecut[0]!=0) {                          if(vertexInTriangles.size()>3 and
2366                          if(K.pointsIntersection[0][0]==&V0                             K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==1) {
2367        or K.pointsIntersection[0][0]==&V1                            ffout(0)<<"case\n";
2368        or K.pointsIntersection[0][0]==&V2 or vertexInTriangles.size()==1) {                            __createNotIn(K,*i);
                            // ffout(0)<<"**********normal\n";  
                         }  
                         else {  
                             if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==3  
           and K.edgecut[0][1]+K.edgecut[1][1]+K.edgecut[2][1]==3) {  
                                 if(K.pointsIntersection[0][1]==K.pointsIntersection[0][0]) {  
                                     K.edgecut[0][0]=0;  
                                     K.edgecut[0][1]=0;  
                                 }  
                                 if(K.pointsIntersection[1][1]==K.pointsIntersection[1][0]) {  
                                     K.edgecut[1][0]=0;  
                                     K.edgecut[1][1]=0;  
                                 }  
                                 if(K.pointsIntersection[2][1]==K.pointsIntersection[2][0]) {  
                                     K.edgecut[2][0]=0;  
                                     K.edgecut[2][1]=0;  
                                 }  
                                 if(vertexInTriangles.size()>3 and  
               K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==2){  
                                     ffout(0)<<"hhhhhhhh\n";  
                                     __createNotIn2(K,*i);  
                                 }  
                                 else if(vertexInTriangles.size()>3 and  
             K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==1) {  
                                     ffout(0)<<"case\n";  
                                     __createNotIn(K,*i);  
                                 }  
                             }  
                             else if(vertexInTriangles.size()>=3) {  
                             ffout(0)<<"triangle cut but outside1\n";  
                             ffout(0)<<V0<<"\n";  
                             ffout(0)<<V1<<"\n";  
                             ffout(0)<<V2<<"\n";  
                             ffout(0)<<"nb points inter "<<vertexInTriangles.size()<<"\n";  
                             }  
2369                          }                          }
2370                          }
2371                        } else {
2372                          if(vertexInTriangles.size()>=3) {
2373                            ffout(0)<<"triangle cut but outside1\n";
2374                            ffout(0)<<V0<<"\n";
2375                            ffout(0)<<V1<<"\n";
2376                            ffout(0)<<V2<<"\n";
2377                            ffout(0)<<"nb points inter "<<vertexInTriangles.size()<<"\n";
2378                          }
2379                      }                      }
2380                      else if(K.edgecut[1]!=0) {                    }
2381                          if(K.pointsIntersection[1][0]==&V0                  } else {
2382        or K.pointsIntersection[1][0]==&V1                    if(K.edgecut[1]!=0) {
2383        or K.pointsIntersection[1][0]==&V2 or vertexInTriangles.size()<3) {                      if(K.pointsIntersection[1][0]==&V0
2384                            //  ffout(0)<<"**********normal\n";                         or K.pointsIntersection[1][0]==&V1
2385                          }                         or K.pointsIntersection[1][0]==&V2 or vertexInTriangles.size()<3) {
2386                          else {                        //  ffout(0)<<"**********normal\n";
2387                              ffout(0)<<"triangle cut but outside2\n";                      } else {
2388                              ffout(0)<<V0<<"\n";                        ffout(0)<<"triangle cut but outside2\n";
2389                              ffout(0)<<V1<<"\n";                        ffout(0)<<V0<<"\n";
2390                              ffout(0)<<V2<<"\n";                        ffout(0)<<V1<<"\n";
2391                              ffout(0)<<"nb points inter "<<vertexInTriangles.size()<<"\n";                        ffout(0)<<V2<<"\n";
2392                          }                        ffout(0)<<"nb points inter "<<vertexInTriangles.size()<<"\n";
2393                      }                      }
2394                      else if(K.edgecut[2]!=0) {                    } else {
2395                          if(K.pointsIntersection[2][0]==&V0                      if(K.edgecut[2]!=0) {
2396        or K.pointsIntersection[2][0]==&V1                        if(K.pointsIntersection[2][0]==&V0
2397        or K.pointsIntersection[2][0]==&V2 or vertexInTriangles.size()<3) {                           or K.pointsIntersection[2][0]==&V1
2398                              //ffout(0)<<"**********normal\n";                           or K.pointsIntersection[2][0]==&V2 or vertexInTriangles.size()<3) {
2399                          }                          //ffout(0)<<"**********normal\n";
2400                          else {                        } else {
2401                              ffout(0)<<"triangle cut but outside3\n";                          ffout(0)<<"triangle cut but outside3\n";
2402                              ffout(0)<<V0<<"\n";                          ffout(0)<<V0<<"\n";
2403                              ffout(0)<<V1<<"\n";                          ffout(0)<<V1<<"\n";
2404                              ffout(0)<<V2<<"\n";                          ffout(0)<<V2<<"\n";
2405                              ffout(0)<<"nb points inter "<<vertexInTriangles.size()<<"\n";                          ffout(0)<<"nb points inter "<<vertexInTriangles.size()<<"\n";
2406                          }                        }
                     }  
                 }  
                 if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==3  
      and K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==1) {  
                     if(K.pointsIntersection[0][0]!=K.pointsIntersection[1][0]  
          and K.pointsIntersection[1][0]!=K.pointsIntersection[2][0]  
          and K.pointsIntersection[2][0]!=K.pointsIntersection[0][0]) {  
                         //ffout(0)<<"---------------- 3 edge cut \n";  
                         //ffout(0)<<V0<<"\n";  
                         //ffout(0)<<V1<<"\n";  
                         //ffout(0)<<V2<<"\n";  
                         //ffout(0)<<"createByEdges\n";  
                         __createByEdges(K,*i);  
                     }  
                     else {  
                         __createListTriangle(K,*i);  
2407                      }                      }
2408                      }
2409                  }                  }
2410                  else if(nbinters!=0) {                }
                     __createListTriangle(K,*i);  
                 }  
                 else {  
                         //ffout(0)<<"j ajoute le tr car pas d'intersection dans les 2 listes\n";  
                         __addTriangle(V0,V1,V2,*i);  
                 }  
2411              }              }
2412              }
2413            }
2414            if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==3
2415               and K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==1) {
2416              if(K.pointsIntersection[0][0]!=K.pointsIntersection[1][0]
2417                 and K.pointsIntersection[1][0]!=K.pointsIntersection[2][0]
2418                 and K.pointsIntersection[2][0]!=K.pointsIntersection[0][0]) {
2419                //ffout(0)<<"---------------- 3 edge cut \n";
2420                //ffout(0)<<V0<<"\n";
2421                //ffout(0)<<V1<<"\n";
2422                //ffout(0)<<V2<<"\n";
2423                //ffout(0)<<"createByEdges\n";
2424                __createByEdges(K,*i);
2425              } else {
2426                __createListTriangle(K,*i);
2427              }
2428            } else {
2429              if(nbinters!=0) {
2430                __createListTriangle(K,*i);
2431              } else {
2432                //ffout(0)<<"j ajoute le tr car pas d'intersection dans les 2 listes\n";
2433                __addTriangle(V0,V1,V2,*i);
2434              }
2435          }          }
2436          }
2437      }      }
2438      if(err==taille)    }
2439          ffout(0)<<"tous les triangles de la liste d'intersection sont dehors\n";  
2440      if(err2==taille)    if(err==taille) {
2441          ffout(0)<<"aucun triangle de la liste d'intersection n'est coupe\n";      ffout(0)<<"tous les triangles de la liste d'intersection sont dehors\n";
2442      }
2443    
2444      if(err2==taille) {
2445        ffout(0)<<"aucun triangle de la liste d'intersection n'est coupe\n";
2446      }
2447    
2448      (*edgesRefVertex[numobject]).clear();    (*edgesRefVertex[numobject]).clear();
2449  }  }
2450    
2451  bool SurfaceMeshGenerator::Internals::  bool SurfaceMeshGenerator::Internals::
# Line 2436  __createTriangleSurface(const size_t& no Line 2455  __createTriangleSurface(const size_t& no
2455                          std::vector<const Edge*>& cutEdges,                          std::vector<const Edge*>& cutEdges,
2456                          const Reference::Representation& ref,                          const Reference::Representation& ref,
2457                          const Cell& currentCell){                          const Cell& currentCell){
2458      bool jobDone=true;    bool jobDone=true;
2459      switch (cutEdges.size()) {    switch (cutEdges.size()) {
2460          case 3: {    case 3: {
2461              //! Stores the Edges on which the vertex will be computed.      //! Stores the Edges on which the vertex will be computed.
2462                            
2463              for (size_t i=0; i<cutEdges.size(); ++i) {      for (size_t i=0; i<cutEdges.size(); ++i) {
2464                  TinyVector<3> V=__splitEdge(*cutEdges[i], ref, S);        TinyVector<3> V=__splitEdge(*cutEdges[i], ref, S);
2465                  if(Omega.inside(V)) {        if(Omega.inside(V)) {
2466                      (*verticesListes[nobject])[*cutEdges[i]]          (*verticesListes[nobject])[*cutEdges[i]]
2467                          =  new Vertex(V);            =  new Vertex(V);
2468                      (*edgesRef[nobject])[*cutEdges[i]] = 1;          (*edgesRef[nobject])[*cutEdges[i]] = 1;
2469                  }        } else {
2470                  else {          // ffout(0)<<"outside\n";
2471                     // ffout(0)<<"outside\n";          jobDone=false;
2472                      jobDone=false;          break;
2473                      break;        }
2474                  }      }
             }  
             //! Addes a triangle to the list.  
             localTriangleList.push_front(TinyVector<3,Edge::Pair>());  
             cellList.push_front(&currentCell);  
   
             std::vector<TinyVector<3> > V(cutEdges.size());  
             for(size_t i=0; i<cutEdges.size(); ++i) {  
                 const Edge& e = *cutEdges[i];  
                 if(e(0).Ref().isSet(ref))  
                     V[i] = (e(0).Vector3() - e(1).Vector3());  
                 else  
                     V[i] = (e(1).Vector3() - e(0).Vector3());  
             }  
2475    
2476              if ((V[0]^V[1])*V[2] < 0)      //! Addes a triangle to the list.
2477                   std::swap(cutEdges[1], cutEdges[2]);      localTriangleList.push_front(TinyVector<3,Edge::Pair>());
2478        cellList.push_front(&currentCell);
2479    
2480              //! Sets the triangle vertices (actually just store the edge where it lives).      std::vector<TinyVector<3> > V(cutEdges.size());
2481              SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle      for(size_t i=0; i<cutEdges.size(); ++i) {
2482                  = localTriangleList.begin();        const Edge& e = *cutEdges[i];
2483          if(e(0).Ref().isSet(ref)) {
2484            V[i] = (e(0).Vector3() - e(1).Vector3());
2485          } else {
2486            V[i] = (e(1).Vector3() - e(0).Vector3());
2487          }
2488        }
2489    
2490              for(size_t i=0; i<cutEdges.size(); ++i)      if ((V[0]^V[1])*V[2] < 0) {
2491                  (*currentTriangle)[i] = (Edge::Pair)(*cutEdges[i]);        std::swap(cutEdges[1], cutEdges[2]);
2492                    }
             break;  
         }  
         case 4: {  
2493    
2494              for (size_t i=0; i<4; ++i) {      //! Sets the triangle vertices (actually just store the edge where it lives).
2495                  TinyVector<3> V=__splitEdge(*cutEdges[i], ref, S);      SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle
2496                  if(Omega.inside(V)) {        = localTriangleList.begin();
                     (*verticesListes[nobject])[*cutEdges[i]]  
                     =  new Vertex(V);  
                     (*edgesRef[nobject])[*cutEdges[i]] = 1;  
                 }  
                 else {  
                     //ffout(0)<<"outside4\n";  
                     jobDone=false;  
                     break;  
                 }  
             }  
             //! Orders the Edges using a map !  
             std::map<Edge::Pair, const Edge*> cutEdgesMap;  
             for(size_t ie=0; ie<cutEdges.size(); ++ie)  
                 cutEdgesMap[(Edge::Pair)*(cutEdges[ie])] = (cutEdges[ie]);  
   
             //! Addes a triangle to the list.  
             localTriangleList.push_front(TinyVector<3,Edge::Pair>());  
             cellList.push_front(&currentCell);  
   
             std::vector<TinyVector<3> > V(cutEdges.size());  
             std::map<Edge::Pair, const Edge*>::iterator currentCutEdge=cutEdgesMap.begin();  
             for (size_t i=0; i<V.size(); ++i) {  
                 const Edge& e = *(*currentCutEdge).second;  
                 if(e(0).Ref().isSet(ref))  
                     V[i] = (e(0).Vector3() - e(1).Vector3());  
                 else  
                     V[i] = (e(1).Vector3() - e(0).Vector3());  
                 ++currentCutEdge;  
             }  
2497    
2498              bool swapEdge = ((V[0]^V[1])*V[2] > 0);      for(size_t i=0; i<cutEdges.size(); ++i) {
2499          (*currentTriangle)[i] = (Edge::Pair)(*cutEdges[i]);
2500        }
2501    
2502              //! Sets the triangle vertices (actually just store the edge where it lives).      break;
2503              SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle    }
                 = localTriangleList.begin();  
   
             currentCutEdge = cutEdgesMap.begin();  
             if (!swapEdge)  
                 for(size_t i=0; i<3; ++i) {  
                     (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);  
                     ++currentCutEdge;  
                 }  
             else  
                 for(int i=2; i>=0; --i) {  
                     (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);  
                     ++currentCutEdge;  
                 }  
2504    
2505                  case 4: {
2506        for (size_t i=0; i<4; ++i) {
2507          TinyVector<3> V=__splitEdge(*cutEdges[i], ref, S);
2508          if(Omega.inside(V)) {
2509            (*verticesListes[nobject])[*cutEdges[i]]
2510              =  new Vertex(V);
2511            (*edgesRef[nobject])[*cutEdges[i]] = 1;
2512          } else {
2513            //ffout(0)<<"outside4\n";
2514            jobDone=false;
2515            break;
2516          }
2517        }
2518    
2519              //! Addes a triangle to the list.      //! Orders the Edges using a map !
2520              localTriangleList.push_front(TinyVector<3,Edge::Pair>());      std::map<Edge::Pair, const Edge*> cutEdgesMap;
2521              cellList.push_front(&currentCell);      for(size_t ie=0; ie<cutEdges.size(); ++ie) {
2522          cutEdgesMap[(Edge::Pair)*(cutEdges[ie])] = (cutEdges[ie]);
2523              currentTriangle = localTriangleList.begin();      }
2524    
2525              //! Sets the triangle vertices (actually just store the edge where it lives).      //! Addes a triangle to the list.
2526              currentCutEdge = cutEdgesMap.begin();      localTriangleList.push_front(TinyVector<3,Edge::Pair>());
2527              if (!swapEdge)      cellList.push_front(&currentCell);
2528                  for(int i=2; i>=0; --i) {  
2529                      ++currentCutEdge;      std::vector<TinyVector<3> > V(cutEdges.size());
2530                      (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);      std::map<Edge::Pair, const Edge*>::iterator currentCutEdge=cutEdgesMap.begin();
2531                  }      for (size_t i=0; i<V.size(); ++i) {
2532              else        const Edge& e = *(*currentCutEdge).second;
2533                  for(size_t i=0; i<3; ++i) {        if(e(0).Ref().isSet(ref)) {
2534                      ++currentCutEdge;          V[i] = (e(0).Vector3() - e(1).Vector3());
2535                      (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);        } else {
2536                  }          V[i] = (e(1).Vector3() - e(0).Vector3());
2537              break;        }
2538          }        ++currentCutEdge;
2539          case 1:      }
2540          case 2: {  
2541              jobDone = false;      bool swapEdge = ((V[0]^V[1])*V[2] > 0);
2542              break;  
2543          }      //! Sets the triangle vertices (actually just store the edge where it lives).
2544        SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle
2545          = localTriangleList.begin();
2546    
2547        currentCutEdge = cutEdgesMap.begin();
2548        if (!swapEdge) {
2549          for(size_t i=0; i<3; ++i) {
2550            (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
2551            ++currentCutEdge;
2552          }
2553        } else {
2554          for(int i=2; i>=0; --i) {
2555            (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
2556            ++currentCutEdge;
2557          }
2558      }      }
2559    
2560      return jobDone;      //! Addes a triangle to the list.
2561            localTriangleList.push_front(TinyVector<3,Edge::Pair>());
2562        cellList.push_front(&currentCell);
2563    
2564        currentTriangle = localTriangleList.begin();
2565    
2566        //! Sets the triangle vertices (actually just store the edge where it lives).
2567        currentCutEdge = cutEdgesMap.begin();
2568        if (!swapEdge) {
2569          for(int i=2; i>=0; --i) {
2570            ++currentCutEdge;
2571            (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
2572          }
2573        } else {
2574          for(size_t i=0; i<3; ++i) {
2575            ++currentCutEdge;
2576            (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
2577          }
2578        }
2579        break;
2580      }
2581      case 1:
2582      case 2: {
2583        jobDone = false;
2584        break;
2585      }
2586      }
2587    
2588      return jobDone;
2589  }  }
2590    
2591  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
# Line 2577  __createSurface(Structured3DMesh& SMesh, Line 2600  __createSurface(Structured3DMesh& SMesh,
2600      Edge& e = SMesh.edge(i);      Edge& e = SMesh.edge(i);
2601  #warning THIS TEST IS NOT ENOUGH! THINK TO THE BORDER!  #warning THIS TEST IS NOT ENOUGH! THINK TO THE BORDER!
2602      if(/*(e(0).Ref().insideDomain() xor e(1).Ref().insideDomain())      if(/*(e(0).Ref().insideDomain() xor e(1).Ref().insideDomain())
2603         and*/((e(0).Ref().isSet(ref)) xor (e(1).Ref().isSet(ref)))) {           and*/((e(0).Ref().isSet(ref)) xor (e(1).Ref().isSet(ref)))) {
2604        (*verticesListes[nobject])[e] = new Vertex(__splitEdge(e, ref, S));        (*verticesListes[nobject])[e] = new Vertex(__splitEdge(e, ref, S));
2605        (*edgesRef[nobject])[e] = 1;        (*edgesRef[nobject])[e] = 1;
               
2606      }      }
2607    }    }
2608            
# Line 2600  __createSurface(Structured3DMesh& SMesh, Line 2622  __createSurface(Structured3DMesh& SMesh,
2622          bool intersected = false;          bool intersected = false;
2623    
2624          const Cell& currentCell = SMesh.cell(icell, jcell, kcell);          const Cell& currentCell = SMesh.cell(icell, jcell, kcell);
2625    
2626          const size_t nx          const size_t nx
2627            = dynamic_cast<const Structured3DMesh&>(SMesh).shape().nx();            = dynamic_cast<const Structured3DMesh&>(SMesh).shape().nx();
2628    
2629          const size_t nx_1 = nx-1;          const size_t nx_1 = nx-1;
2630    
2631          const size_t ny          const size_t ny
2632            = dynamic_cast<const Structured3DMesh&>(SMesh).shape().ny();            = dynamic_cast<const Structured3DMesh&>(SMesh).shape().ny();
2633    
2634          const size_t ny_1 = ny-1;          const size_t ny_1 = ny-1;
2635    
2636          const size_t nz          const size_t nz
2637            = dynamic_cast<const Structured3DMesh&>(SMesh).shape().nz();            = dynamic_cast<const Structured3DMesh&>(SMesh).shape().nz();
2638    
2639          const size_t nz_1 = nz-1;          const size_t nz_1 = nz-1;
2640    
2641          HexahedronByEdges H (SMesh.edge(icell*ny*nz + jcell*nz + kcell),          HexahedronByEdges H (SMesh.edge(icell*ny*nz + jcell*nz + kcell),
# Line 2629  __createSurface(Structured3DMesh& SMesh, Line 2655  __createSurface(Structured3DMesh& SMesh,
2655          for (size_t l=0; l<H.numberOfEdges(); l++) {          for (size_t l=0; l<H.numberOfEdges(); l++) {
2656            Edge& e = H.edge(l);            Edge& e = H.edge(l);
2657            if ((*verticesListes[nobject]).count(e) > 0) {            if ((*verticesListes[nobject]).count(e) > 0) {
2658                intersected = true;              intersected = true;
2659                cuttedHexahedra.insert(&currentCell);              cuttedHexahedra.insert(&currentCell);
2660              break;              break;
2661            }            }
2662          }          }
# Line 2754  __createSurface(Structured3DMesh& SMesh, Line 2780  __createSurface(Structured3DMesh& SMesh,
2780                                        addedEdges[3],                                        addedEdges[3],
2781                                        addedEdges[4],                                        addedEdges[4],
2782                                        addedEdges[5]);                                        addedEdges[5]);
   
2783            }            }
2784            }            }
2785    
# Line 2765  __createSurface(Structured3DMesh& SMesh, Line 2790  __createSurface(Structured3DMesh& SMesh,
2790    
2791              for(size_t ie=0; ie<6; ++ie) {              for(size_t ie=0; ie<6; ++ie) {
2792                const Edge& e = currentT->edge(ie);                const Edge& e = currentT->edge(ie);
2793                if ((*edgesRef[nobject]).find(e) != (*edgesRef[nobject]).end())                if ((*edgesRef[nobject]).find(e) != (*edgesRef[nobject]).end()) {
2794                  cutEdges.push_back(&e);                  cutEdges.push_back(&e);
2795                  }
2796              }              }
2797    
2798              bool jobDone2=__createTriangleSurface(nobject,S,Omega,cutEdges,              bool jobDone2=__createTriangleSurface(nobject,S,Omega,cutEdges,
2799                                       ref,currentCell);                                       ref,currentCell);
2800              if(!(jobDone2))              if(!(jobDone2)) {
2801                  jobDone=jobDone2;                  jobDone=jobDone2;
2802                }
2803            }            }
2804                            
2805            if(jobDone) {            if(jobDone) {
2806               // meshedHexahedra.insert(&currentCell);              // meshedHexahedra.insert(&currentCell);
2807                SurfaceMeshGenerator::Internals::MotherCellList::iterator itcell              SurfaceMeshGenerator::Internals::MotherCellList::iterator itcell
2808                =cellList.begin();                =cellList.begin();
2809                for (SurfaceMeshGenerator::Internals::EdgePairList::iterator t              for (SurfaceMeshGenerator::Internals::EdgePairList::iterator t
2810              = localTriangleList.begin();                     = localTriangleList.begin();
2811              t != localTriangleList.end(); ++t) {                   t != localTriangleList.end(); ++t) {
2812                    (*triangleListes[nobject]).push_front(*t);                (*triangleListes[nobject]).push_front(*t);
2813                    (*cellListObject[nobject]).push_front(*itcell);                (*cellListObject[nobject]).push_front(*itcell);
2814                    ++itcell;                ++itcell;
2815                }              }
2816                hexaTetraAssociation[&currentCell]              hexaTetraAssociation[&currentCell]
2817                = std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,                = std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,
2818                ReferenceCounting<std::vector<Edge> > >                ReferenceCounting<std::vector<Edge> > >
2819                (pT,pAddedEdges);                (pT,pAddedEdges);
2820            } else {            } else {
2821                hexaTetraAssociation[&currentCell]              hexaTetraAssociation[&currentCell]
2822                = std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,                = std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,
2823                ReferenceCounting<std::vector<Edge> > >                ReferenceCounting<std::vector<Edge> > >
2824                (pT,pAddedEdges);                (pT,pAddedEdges);
# Line 2801  __createSurface(Structured3DMesh& SMesh, Line 2827  __createSurface(Structured3DMesh& SMesh,
2827            cellList.clear();            cellList.clear();
2828          }          }
2829        }        }
   
2830  }  }
2831    
2832  void  void
# Line 2907  generatesAllSufaceMeshes(const Problem& Line 2932  generatesAllSufaceMeshes(const Problem&
2932            = dynamic_cast<const BoundaryStructured3DMesh&>(B);            = dynamic_cast<const BoundaryStructured3DMesh&>(B);
2933          for (size_t k=0; k<6; ++k) {          for (size_t k=0; k<6; ++k) {
2934            Structured3DMesh& SMesh = dynamic_cast<Structured3DMesh&>(M);            Structured3DMesh& SMesh = dynamic_cast<Structured3DMesh&>(M);
2935            if(bs.structMesh(k))            if(bs.structMesh(k)) {
2936              bcMeshAssociation.addPDEBoundaryConditionAndMesh(i,&BC,&SMesh.surfaceMeshOfQuadrangles(k));              bcMeshAssociation.addPDEBoundaryConditionAndMesh(i,&BC,&SMesh.surfaceMeshOfQuadrangles(k));
2937              }
2938          }          }
2939          break;          break;
2940        }        }
# Line 3049  generateSurfacicMesh(const Object& O, Line 3075  generateSurfacicMesh(const Object& O,
3075    const Shape& S = *(O.shape());    const Shape& S = *(O.shape());
3076            
3077    typedef std::map<TinyVector<3>, Object*,    typedef std::map<TinyVector<3>, Object*,
3078    SurfaceMeshGenerator::Internals::__ltv> ObjectReferences;      SurfaceMeshGenerator::Internals::__ltv> ObjectReferences;
3079    
3080    ObjectReferences otherReferences;    ObjectReferences otherReferences;
3081    
# Line 3069  generateSurfacicMesh(const Object& O, Line 3095  generateSurfacicMesh(const Object& O,
3095      Internals::__getIntersectionReferences(I, otherReferences);      Internals::__getIntersectionReferences(I, otherReferences);
3096      break;      break;
3097    } case Shape::difference: {    } case Shape::difference: {
3098        fferr(0) << __FILE__ << ':' << __LINE__ << ": Not implemented\n";      fferr(0) << __FILE__ << ':' << __LINE__ << ": Not implemented\n";
3099        break;      break;
3100    } case Shape::not_: {    } case Shape::not_: {
3101      //  fferr(0) << __FILE__ << ':' << __LINE__ << ": Not implemented\n";      //  fferr(0) << __FILE__ << ':' << __LINE__ << ": Not implemented\n";
3102        break;      break;
3103    } default: {    } default: {
3104    //std::exit(1);      //std::exit(1);
3105    }    }
3106    }        }    
3107    
# Line 3095  generateSurfacicMesh(const Object& O, Line 3121  generateSurfacicMesh(const Object& O,
3121    (*__internals).edgesRefVertex.resize(nobject);    (*__internals).edgesRefVertex.resize(nobject);
3122    
3123    for(size_t no=0 ; no<nobject ; ++no) {    for(size_t no=0 ; no<nobject ; ++no) {
3124        (*__internals).triangleListes[no]= new SurfaceMeshGenerator::Internals::EdgePairList;      (*__internals).triangleListes[no]= new SurfaceMeshGenerator::Internals::EdgePairList;
3125        (*__internals).localTriangleListesEdges[no]= new SurfaceMeshGenerator::Internals::EdgePairList;      (*__internals).localTriangleListesEdges[no]= new SurfaceMeshGenerator::Internals::EdgePairList;
3126        (*__internals).verticesListes[no]= new SurfaceMeshGenerator::Internals::VerticesList;      (*__internals).verticesListes[no]= new SurfaceMeshGenerator::Internals::VerticesList;
3127        (*__internals).verticesListesIntersection[no]= new SurfaceMeshGenerator::Internals::VerticesList;      (*__internals).verticesListesIntersection[no]= new SurfaceMeshGenerator::Internals::VerticesList;
3128        (*__internals).cellListObject[no]= new SurfaceMeshGenerator::Internals::MotherCellList;      (*__internals).cellListObject[no]= new SurfaceMeshGenerator::Internals::MotherCellList;
3129    
3130        (*__internals).edgesRef[no]= new std::map<Edge::Pair, size_t> ;      (*__internals).edgesRef[no]= new std::map<Edge::Pair, size_t> ;
3131        (*__internals).edgesRefVertex[no] = new std::map<const Vertex*,size_t> ;      (*__internals).edgesRefVertex[no] = new std::map<const Vertex*,size_t> ;
3132    }    }
3133    
3134        
# Line 3114  generateSurfacicMesh(const Object& O, Line 3140  generateSurfacicMesh(const Object& O,
3140    
3141    int num=0;    int num=0;
3142    if(nobjectToTreat<2) {    if(nobjectToTreat<2) {
3143        (*__internals).__createSurface(SMesh,Omega,S,ref,num,      (*__internals).__createSurface(SMesh,Omega,S,ref,num,
                                      cuttedHexahedra,meshedHexahedra);  
   }  
   else {  
       for(ObjectReferences::iterator i = otherReferences.begin();  
           i != otherReferences.end(); ++i) {  
           for(size_t n = 0; n<scene.nbObjects(); ++n) {  
               Reference::Representation otherRef;  
               otherRef.resize(scene.nbObjects());  
               //ffout(0) << scene.object(n).reference() << ":"<< ((*i).first) << '\n';  
               if (scene.object(n).reference() == ((*i).first)) {  
                   otherRef[n] = true;  
                   //create surface objects  
                   (*__internals).__createSurface(SMesh,Omega,  
                                    *(scene.object(n).shape()),otherRef,num,  
3144                                     cuttedHexahedra,meshedHexahedra);                                     cuttedHexahedra,meshedHexahedra);
3145                }    } else {
3146            }      for(ObjectReferences::iterator i = otherReferences.begin();
3147            plothexa(num,cuttedHexahedra);          i != otherReferences.end(); ++i) {
3148            if(num==0)        for(size_t n = 0; n<scene.nbObjects(); ++n) {
3149                plothexa(num-1,cuttedHexahedra);          Reference::Representation otherRef;
3150            num++;          otherRef.resize(scene.nbObjects());
3151        }          //ffout(0) << scene.object(n).reference() << ":"<< ((*i).first) << '\n';
3152            if (scene.object(n).reference() == ((*i).first)) {
3153              otherRef[n] = true;
3154              //create surface objects
3155              (*__internals).__createSurface(SMesh,Omega,
3156                                             *(scene.object(n).shape()),otherRef,num,
3157                                             cuttedHexahedra,meshedHexahedra);
3158            }
3159          }
3160          plothexa(num,cuttedHexahedra);
3161          if(num==0) {
3162            plothexa(num-1,cuttedHexahedra);
3163          }
3164          num++;
3165        }
3166    }    }
3167    
3168  std::set_difference (cuttedHexahedra.begin(), cuttedHexahedra.end(),    std::set_difference (cuttedHexahedra.begin(), cuttedHexahedra.end(),
3169                       meshedHexahedra.begin(), meshedHexahedra.end(),                         meshedHexahedra.begin(), meshedHexahedra.end(),
3170                       inserter((*__internals).toTreatHexahedra, (*__internals).toTreatHexahedra.begin()));                         inserter((*__internals).toTreatHexahedra, (*__internals).toTreatHexahedra.begin()));
3171    
3172  //  ffout(0)<<(*__internals).toTreatHexahedra.size()<<"\n";    //  ffout(0)<<(*__internals).toTreatHexahedra.size()<<"\n";
3173  //  (*__internals).toTreatHexahedra=meshedHexahedra;    //  (*__internals).toTreatHexahedra=meshedHexahedra;
3174  /*  std::set<const Cell*>::iterator iterat=(*__internals).toTreatHexahedra.begin();    /*  std::set<const Cell*>::iterator iterat=(*__internals).toTreatHexahedra.begin();
3175    size_t treatsize=(*__internals).toTreatHexahedra.size();        size_t treatsize=(*__internals).toTreatHexahedra.size();
3176    for(size_t yy=0 ; yy<treatsize ; ++yy)        for(size_t yy=0 ; yy<treatsize ; ++yy)
3177    {        {
3178        SurfaceMeshGenerator::Internals::HexaTetraAssociation::iterator        SurfaceMeshGenerator::Internals::HexaTetraAssociation::iterator
3179        j=(*__internals).hexaTetraAssociation.find(*iterat);        j=(*__internals).hexaTetraAssociation.find(*iterat);
3180        std::vector<TetrahedronByEdges>& T = *(((*j).second).first);        std::vector<TetrahedronByEdges>& T = *(((*j).second).first);
3181        const CartesianHexahedron& currentCell        const CartesianHexahedron& currentCell
3182            = static_cast<const CartesianHexahedron&>(*(*iterat));        = static_cast<const CartesianHexahedron&>(*(*iterat));
3183       // ffout(0)<<"*********************************\n";        // ffout(0)<<"*********************************\n";
3184        size_t test=0;        size_t test=0;
3185        std::set<const Cell*>::iterator temp=iterat++;        std::set<const Cell*>::iterator temp=iterat++;
3186        for (std::vector<TetrahedronByEdges>::iterator currentT = T.begin();        for (std::vector<TetrahedronByEdges>::iterator currentT = T.begin();
3187             currentT < T.end(); ++currentT) {        currentT < T.end(); ++currentT) {
3188    
3189            std::vector<ReferenceCounting< std::vector<const Edge*> > > cutEdges;        std::vector<ReferenceCounting< std::vector<const Edge*> > > cutEdges;
3190            cutEdges.resize(nobject);        cutEdges.resize(nobject);
3191            for(size_t no=0 ; no<nobject ; ++no) {        for(size_t no=0 ; no<nobject ; ++no) {
3192                cutEdges[no]= new std::vector<const Edge*>;        cutEdges[no]= new std::vector<const Edge*>;
3193            }        }
3194            (*__internals).__createCutEdges(scene,cutEdges,otherReferences,currentT);        (*__internals).__createCutEdges(scene,cutEdges,otherReferences,currentT);
3195            for(size_t no=0 ; no<nobject ; ++no) {        for(size_t no=0 ; no<nobject ; ++no) {
3196                if((*cutEdges[no]).size()>0) {        if((*cutEdges[no]).size()>0) {
3197                    test++;        test++;
3198                }        }
3199            }        }
3200            if(test<2)        if(test<2)
3201                test=0;        test=0;
3202            else        else
3203                test=2;        test=2;
3204        }        }
3205        if(test!=2) {        if(test!=2) {
3206            (*__internals).toTreatHexahedra.erase(temp);        (*__internals).toTreatHexahedra.erase(temp);
3207        }        }
3208    }*/        }*/
3209    ffout(0)<<"nb tr objet 0 complet "<<(*(*__internals).triangleListes[0]).size()<<" \n";    ffout(0)<<"nb tr objet 0 complet "<<(*(*__internals).triangleListes[0]).size()<<" \n";
3210            
3211    size_t nobjother=0;    size_t nobjother=0;
3212    for(size_t no=0;no<nobjectToTreat;++no) {    for(size_t no=0;no<nobjectToTreat;++no) {
3213        nobjother+=(*(*__internals).triangleListes[no]).size();      nobjother+=(*(*__internals).triangleListes[no]).size();
3214    }    }
3215   // ffout(0)<<"nb tr objet1 complet "<<nobjother<<" \n";   // ffout(0)<<"nb tr objet1 complet "<<nobjother<<" \n";
3216    
# Line 3225  std::set_difference (cuttedHexahedra.beg Line 3251  std::set_difference (cuttedHexahedra.beg
3251        Vector<Triangle>& triangleList = *pTriangleList;        Vector<Triangle>& triangleList = *pTriangleList;
3252    
3253        SurfaceMeshGenerator::Internals::MotherCellList::iterator        SurfaceMeshGenerator::Internals::MotherCellList::iterator
3254            itcell=(*(*__internals).cellListObject[numObj]).begin();          itcell=(*(*__internals).cellListObject[numObj]).begin();
3255        for (SurfaceMeshGenerator::Internals::EdgePairList::iterator        for (SurfaceMeshGenerator::Internals::EdgePairList::iterator
3256               i = (*(*__internals).triangleListes[numObj]).begin();               i = (*(*__internals).triangleListes[numObj]).begin();
3257             i != (*(*__internals).triangleListes [numObj]).end(); ++i, ++n) {             i != (*(*__internals).triangleListes [numObj]).end(); ++i, ++n) {
# Line 3245  std::set_difference (cuttedHexahedra.beg Line 3271  std::set_difference (cuttedHexahedra.beg
3271      (*(*__internals).triangleListes[numObj]).clear();      (*(*__internals).triangleListes[numObj]).clear();
3272    
3273      // Intersection lists      // Intersection lists
3274          
     
3275      {      {
3276  #warning keep this ?  #warning keep this ?
3277        size_t n=0;        size_t n=0;
# Line 3256  std::set_difference (cuttedHexahedra.beg Line 3281  std::set_difference (cuttedHexahedra.beg
3281        Vector<Vertex*>& vertexListIntersection = *pVertexListIntersection;        Vector<Vertex*>& vertexListIntersection = *pVertexListIntersection;
3282    
3283        //for not have double points        //for not have double points
3284          for (SurfaceMeshGenerator::Internals::VerticesList::iterator        for (SurfaceMeshGenerator::Internals::VerticesList::iterator
3285        i = (*(*__internals).verticesListesIntersection[numObj]).begin();               i = (*(*__internals).verticesListesIntersection[numObj]).begin();
3286        i != (*(*__internals).verticesListesIntersection[numObj]).end(); ++i, ++n) {             i != (*(*__internals).verticesListesIntersection[numObj]).end(); ++i, ++n) {
3287              if((*(*__internals).verticesListes[numObj]).find((*i).first)          if ((*(*__internals).verticesListes[numObj]).find((*i).first)
3288          !=(*(*__internals).verticesListes[numObj]).end()) {              !=(*(*__internals).verticesListes[numObj]).end()) {
3289                  //ffout(0)<<"point existe deja\n";            //ffout(0)<<"point existe deja\n";
3290                  SurfaceMeshGenerator::Internals::VerticesList::iterator            SurfaceMeshGenerator::Internals::VerticesList::iterator
3291                  j=(*(*__internals).verticesListes[numObj]).find((*i).first);              j=(*(*__internals).verticesListes[numObj]).find((*i).first);
3292                  vertexListIntersection[n] = (*j).second;            vertexListIntersection[n] = (*j).second;
3293              }          } else {
3294              else {            vertexListIntersection[n] = (*i).second;
3295                  vertexListIntersection[n] = (*i).second;            (*(*__internals).verticesListes[numObj])[(*i).first]
3296                  (*(*__internals).verticesListes[numObj])[(*i).first]              =(*i).second;              
                                       =(*i).second;              
             }  
3297          }          }
3298          }
3299      }      }
3300      (*(*__internals).verticesListesIntersection[numObj]).clear();      (*(*__internals).verticesListesIntersection[numObj]).clear();
3301      {      {
# Line 3302  std::set_difference (cuttedHexahedra.beg Line 3326  std::set_difference (cuttedHexahedra.beg
3326                 j = (*(*i).second[numObj]).begin();                 j = (*(*i).second[numObj]).begin();
3327               j !=  (*(*i).second[numObj]).end(); ++j) {               j !=  (*(*i).second[numObj]).end(); ++j) {
3328            SurfaceMeshGenerator::Internals::VerticesList& V = (*(*__internals).verticesListes[numObj]);            SurfaceMeshGenerator::Internals::VerticesList& V = (*(*__internals).verticesListes[numObj]);
3329              Triangle T(* V[(*j)[0]],            Triangle T(* V[(*j)[0]],
3330                  * V[(*j)[1]],                       * V[(*j)[1]],
3331                  * V[(*j)[2]], 1);                       * V[(*j)[2]], 1);
3332                            
3333              T.setMother((*i).first);            T.setMother((*i).first);
3334    
3335              triangleListIntersection[n] = T;            triangleListIntersection[n] = T;
3336    
3337              hexaToTriangle[(*i).first].push_back(&triangleListIntersection[n]);            hexaToTriangle[(*i).first].push_back(&triangleListIntersection[n]);
3338              ++n;            ++n;
3339          }          }
3340        }        }
3341      }      }
3342    
3343      (*(*__internals).verticesListes[numObj]).clear();      (*(*__internals).verticesListes[numObj]).clear();
3344    }    }
   (*__internals).triangleListesIntersectionEdgesPerHexa.clear();  
3345    
3346      (*__internals).triangleListesIntersectionEdgesPerHexa.clear();
3347    (*__internals).plotliste();    (*__internals).plotliste();
3348    
3349    size_t n0=0,n1=0;    size_t n0=0,n1=0;
3350    if(nobjectToTreat>1) {    if(nobjectToTreat>1) {
3351        //for(size_t ncompt=0 ; ncompt<=nobjectToTreat-5; ++ncompt) {      //for(size_t ncompt=0 ; ncompt<=nobjectToTreat-5; ++ncompt) {
3352        for(size_t ncompt=0 ; ncompt<=nobjectToTreat-2; ++ncompt) {      for(size_t ncompt=0 ; ncompt<=nobjectToTreat-2; ++ncompt) {
3353            ffout(0)<<"treat object "<<ncompt+1<<"\n"<<std::flush;        ffout(0)<<"treat object "<<ncompt+1<<"\n"<<std::flush;
3354            //calculate points intersection between object 0 and object 1        //calculate points intersection between object 0 and object 1
3355            (*__internals).__calculatePointsIntersection(scene,0,1);        (*__internals).__calculatePointsIntersection(scene,0,1);
3356            //  ffout(0)<<"end calculatePointsIntersection\n";        //  ffout(0)<<"end calculatePointsIntersection\n";
3357     //create mesh intersection for the object 0 and 1        //create mesh intersection for the object 0 and 1
3358            (*__internals).__createTrianglesIntersection(0,1);        (*__internals).__createTrianglesIntersection(0,1);
3359            ffout(0)<<"nb triangles in intersection 0 "<<(*__internals).triangleListIntersectionNew.size()<<"\n"<<std::flush;        ffout(0)<<"nb triangles in intersection 0 "<<(*__internals).triangleListIntersectionNew.size()<<"\n"<<std::flush;
3360            (*__internals).__createTrianglesIntersection(1,0);        (*__internals).__createTrianglesIntersection(1,0);
3361            ffout(0)<<"nb triangles in intersection 1 "<<(*__internals).triangleListIntersectionNew.size()<<"\n"<<std::flush;        ffout(0)<<"nb triangles in intersection 1 "<<(*__internals).triangleListIntersectionNew.size()<<"\n"<<std::flush;
3362    
3363    //erase triangles from object 0 and 1 and add triangle from intersection of objects 0-1        //erase triangles from object 0 and 1 and add triangle from intersection of objects 0-1
3364            SurfaceMeshGenerator::Internals::ListOfTriangleList::iterator        SurfaceMeshGenerator::Internals::ListOfTriangleList::iterator
3365                it=(*__internals).listOfTriangleListIntersection.begin();          it=(*__internals).listOfTriangleListIntersection.begin();
3366            SurfaceMeshGenerator::Internals::ListOfTriangleList::iterator        SurfaceMeshGenerator::Internals::ListOfTriangleList::iterator
3367                er=(*__internals).listOfTriangleListIntersection.erase(it);          er=(*__internals).listOfTriangleListIntersection.erase(it);
3368            it=(*__internals).listOfTriangleListIntersection.begin();        it=(*__internals).listOfTriangleListIntersection.begin();
3369            er=(*__internals).listOfTriangleListIntersection.erase(it);        er=(*__internals).listOfTriangleListIntersection.erase(it);
3370    
3371            size_t sizeinter=(*__internals).triangleListIntersectionNew.size();        size_t sizeinter=(*__internals).triangleListIntersectionNew.size();
3372            ReferenceCounting<Vector<Triangle> > pTriangleListIntersection = new Vector<Triangle>(sizeinter);        ReferenceCounting<Vector<Triangle> > pTriangleListIntersection = new Vector<Triangle>(sizeinter);
3373            (*__internals).listOfTriangleListIntersection.push_front(pTriangleListIntersection);        (*__internals).listOfTriangleListIntersection.push_front(pTriangleListIntersection);
3374            Vector<Triangle>& triangleListIntersection = *pTriangleListIntersection;        Vector<Triangle>& triangleListIntersection = *pTriangleListIntersection;
3375    
3376            SurfaceMeshGenerator::Internals::ListOfHexaToTriangle::iterator        SurfaceMeshGenerator::Internals::ListOfHexaToTriangle::iterator
3377                ithexa=(*__internals).listOfHexaToTriangle.begin();          ithexa=(*__internals).listOfHexaToTriangle.begin();
3378            SurfaceMeshGenerator::Internals::ListOfHexaToTriangle::iterator        SurfaceMeshGenerator::Internals::ListOfHexaToTriangle::iterator
3379                erhexa=(*__internals).listOfHexaToTriangle.erase(ithexa);          erhexa=(*__internals).listOfHexaToTriangle.erase(ithexa);
3380            ithexa=(*__internals).listOfHexaToTriangle.begin();        ithexa=(*__internals).listOfHexaToTriangle.begin();
3381            erhexa=(*__internals).listOfHexaToTriangle.erase(ithexa);        erhexa=(*__internals).listOfHexaToTriangle.erase(ithexa);
3382            (*__internals).listOfHexaToTriangle.push_front(new std::map<const Cell*, std::list<Triangle*> >);        (*__internals).listOfHexaToTriangle.push_front(new std::map<const Cell*, std::list<Triangle*> >);
3383            ithexa=(*__internals).listOfHexaToTriangle.begin();        ithexa=(*__internals).listOfHexaToTriangle.begin();
3384    
3385            for(size_t size=0 ; size<sizeinter ; ++size) {        for(size_t size=0 ; size<sizeinter ; ++size) {
3386                triangleListIntersection[size]=(*__internals).triangleListIntersectionNew[size];          triangleListIntersection[size]=(*__internals).triangleListIntersectionNew[size];
3387                (*(*ithexa))[&triangleListIntersection[size].mother()].push_back(&triangleListIntersection[size]);          (*(*ithexa))[&triangleListIntersection[size].mother()].push_back(&triangleListIntersection[size]);
3388            }        }
3389            (*__internals).triangleListIntersectionNew.clear();        (*__internals).triangleListIntersectionNew.clear();
3390    
3391            //copy meshes not intersected of object 0 and 1        //copy meshes not intersected of object 0 and 1
3392            ReferenceCounting<Vector<Triangle> > triangleListTemp;        ReferenceCounting<Vector<Triangle> > triangleListTemp;
3393            size_t size=0;        size_t size=0;
3394            SurfaceMeshGenerator::Internals::ListOfTriangleList::iterator        SurfaceMeshGenerator::Internals::ListOfTriangleList::iterator
3395                triangleit=(*__internals).listOfTriangleList.begin();          triangleit=(*__internals).listOfTriangleList.begin();
3396            size+=(*(*triangleit)).size();        size+=(*(*triangleit)).size();
3397            n0=size;        n0=size;
3398            ++triangleit;        ++triangleit;
3399            size+=(*(*triangleit)).size();        size+=(*(*triangleit)).size();
3400            n1=size;        n1=size;
3401            triangleListTemp=new Vector<Triangle>(size);        triangleListTemp=new Vector<Triangle>(size);
3402            triangleit=(*__internals).listOfTriangleList.begin();        triangleit=(*__internals).listOfTriangleList.begin();
3403            for(size_t s=0 ; s<(*(*triangleit)).size() ; ++s) {        for(size_t s=0 ; s<(*(*triangleit)).size() ; ++s) {
3404                (*triangleListTemp)[s]=(*(*triangleit))[s];          (*triangleListTemp)[s]=(*(*triangleit))[s];
3405            }        }
3406            size_t numsize=(*(*triangleit)).size();        size_t numsize=(*(*triangleit)).size();
3407            ++triangleit;        ++triangleit;
3408            for(size_t s=0 ; s<(*(*triangleit)).size() ; ++s) {        for(size_t s=0 ; s<(*(*triangleit)).size() ; ++s) {
3409                (*triangleListTemp)[numsize+s]=(*(*triangleit))[s];          (*triangleListTemp)[numsize+s]=(*(*triangleit))[s];
3410            }        }
3411    
3412            it=(*__internals).listOfTriangleList.begin();        it=(*__internals).listOfTriangleList.begin();
3413            er=(*__internals).listOfTriangleList.erase(it);        er=(*__internals).listOfTriangleList.erase(it);
3414            it=(*__internals).listOfTriangleList.begin();        it=(*__internals).listOfTriangleList.begin();
3415            er=(*__internals).listOfTriangleList.erase(it);        er=(*__internals).listOfTriangleList.erase(it);
3416            (*__internals).listOfTriangleList.push_front(new Vector<Triangle>(size));        (*__internals).listOfTriangleList.push_front(new Vector<Triangle>(size));
3417            it=(*__internals).listOfTriangleList.begin();        it=(*__internals).listOfTriangleList.begin();
3418            for(size_t s=0 ; s<(*triangleListTemp).size() ; ++s) {        for(size_t s=0 ; s<(*triangleListTemp).size() ; ++s) {
3419                (*(*it))[s]=(*triangleListTemp)[s];          (*(*it))[s]=(*triangleListTemp)[s];
           }  
3420        }        }
3421        }
3422    }    }
3423    //end intersection mesh    //end intersection mesh
3424    
3425    //construction of list vertices    //construction of list vertices
3426        
3427    SurfaceMeshGenerator::Internals::ListOfTriangleList::const_iterator    SurfaceMeshGenerator::Internals::ListOfTriangleList::const_iterator
3428        itcell=((*__internals).listOfTriangleList).begin();      itcell=((*__internals).listOfTriangleList).begin();
3429    SurfaceMeshGenerator::Internals::ListOfTriangleList::const_iterator    SurfaceMeshGenerator::Internals::ListOfTriangleList::const_iterator
3430        intertr=((*__internals).listOfTriangleListIntersection).begin();      intertr=((*__internals).listOfTriangleListIntersection).begin();
3431    if(nobjectToTreat==0){    if (nobjectToTreat==0){
3432    for(size_t size=0 ; size< (*(*itcell)).size() ; ++size) {      for(size_t size=0 ; size< (*(*itcell)).size() ; ++size) {
3433        const Triangle& T=(*(*itcell))[size];        const Triangle& T=(*(*itcell))[size];
3434        if((*__internals).listOfVertexMesh.find(&T(0))==(*__internals).listOfVertexMesh.end()) {        if ((*__internals).listOfVertexMesh.find(&T(0))==(*__internals).listOfVertexMesh.end()) {
3435            (*__internals).listOfVertexMesh[&T(0)]=(&T(0));          (*__internals).listOfVertexMesh[&T(0)]=(&T(0));
3436        }        }
3437        if((*__internals).listOfVertexMesh.find(&T(1))==(*__internals).listOfVertexMesh.end()) {        if ((*__internals).listOfVertexMesh.find(&T(1))==(*__internals).listOfVertexMesh.end()) {
3438            (*__internals).listOfVertexMesh[&T(1)]=(&T(1));          (*__internals).listOfVertexMesh[&T(1)]=(&T(1));
3439        }        }
3440        if((*__internals).listOfVertexMesh.find(&T(2))==(*__internals).listOfVertexMesh.end()) {        if ((*__internals).listOfVertexMesh.find(&T(2))==(*__internals).listOfVertexMesh.end()) {
3441            (*__internals).listOfVertexMesh[&T(2)]=(&T(2));          (*__internals).listOfVertexMesh[&T(2)]=(&T(2));
3442        }        }
3443        }
3444    }    }
3445    }  
3446    for(size_t size=0 ; size< (*(*intertr)).size() ; ++size) {    for(size_t size=0 ; size< (*(*intertr)).size() ; ++size) {
3447        const Triangle& T=(*(*intertr))[size];      const Triangle& T=(*(*intertr))[size];
3448        if((*__internals).listOfVertexMesh.find(&T(0))==(*__internals).listOfVertexMesh.end()) {      if ((*__internals).listOfVertexMesh.find(&T(0))==(*__internals).listOfVertexMesh.end()) {
3449            (*__internals).listOfVertexMesh[&T(0)]=(&T(0));        (*__internals).listOfVertexMesh[&T(0)]=(&T(0));
3450        }      }
3451        if((*__internals).listOfVertexMesh.find(&T(1))==(*__internals).listOfVertexMesh.end()) {      if ((*__internals).listOfVertexMesh.find(&T(1))==(*__internals).listOfVertexMesh.end()) {
3452            (*__internals).listOfVertexMesh[&T(1)]=(&T(1));        (*__internals).listOfVertexMesh[&T(1)]=(&T(1));
3453        }      }
3454        if((*__internals).listOfVertexMesh.find(&T(2))==(*__internals).listOfVertexMesh.end()) {      if ((*__internals).listOfVertexMesh.find(&T(2))==(*__internals).listOfVertexMesh.end()) {
3455            (*__internals).listOfVertexMesh[&T(2)]=(&T(2));        (*__internals).listOfVertexMesh[&T(2)]=(&T(2));
3456        }      }
3457    }    }
3458   /* Vector<const Vertex*> listVertexFinal((*__internals).listOfVertexMesh.size());    /* Vector<const Vertex*> listVertexFinal((*__internals).listOfVertexMesh.size());
3459    Vector<Vertex> listVertexFinal2((*__internals).listOfVertexMesh.size());       Vector<Vertex> listVertexFinal2((*__internals).listOfVertexMesh.size());
3460    size_t number=0;       size_t number=0;
3461    for(std::map<const Vertex*,const Vertex*>::iterator it=(*__internals).listOfVertexMesh.begin() ;       for(std::map<const Vertex*,const Vertex*>::iterator it=(*__internals).listOfVertexMesh.begin() ;
3462        it!=(*__internals).listOfVertexMesh.end() ; it++) {       it!=(*__internals).listOfVertexMesh.end() ; it++) {
3463        listVertexFinal[number]=(*it).second;       listVertexFinal[number]=(*it).second;
3464        listVertexFinal2[number]=*((*it).second);       listVertexFinal2[number]=*((*it).second);
3465        ++number;       ++number;
3466    }*/       }*/
3467     size_t ntr=0;    size_t ntr=0;
3468    
3469     //ntr+=listVertexFinal.size();    //ntr+=listVertexFinal.size();
3470    ntr+=(*__internals).listOfVertexMesh.size();    ntr+=(*__internals).listOfVertexMesh.size();
3471        
3472     s_mesh.setNumberOfVertices(ntr);    s_mesh.setNumberOfVertices(ntr);
3473     ffout(0)<<"nb tot vertices "<<s_mesh.numberOfVertices()<<"\n"<<std::flush;    ffout(0)<<"nb tot vertices "<<s_mesh.numberOfVertices()<<"\n"<<std::flush;
3474    
3475     //! copies vertices into the s_mesh.    //! copies vertices into the s_mesh.
3476     int i = 0;    int i = 0;
3477    /* for(size_t v=0 ; v<ntr ; ++v) {    /* for(size_t v=0 ; v<ntr ; ++v) {
3478         s_mesh.vertex(i) = *listVertexFinal[v];       s_mesh.vertex(i) = *listVertexFinal[v];
3479         listVertexFinal[v] = &s_mesh.vertex(i);       listVertexFinal[v] = &s_mesh.vertex(i);
3480         i++;       i++;
3481       }*/       }*/
3482     for(std::map<const Vertex*,const Vertex*>::iterator it=(*__internals).listOfVertexMesh.begin() ;    for(std::map<const Vertex*,const Vertex*>::iterator it=(*__internals).listOfVertexMesh.begin() ;
3483         it!=(*__internals).listOfVertexMesh.end() ; it++) {        it!=(*__internals).listOfVertexMesh.end() ; it++) {
3484         s_mesh.vertex(i) = *((*it).second);      s_mesh.vertex(i) = *((*it).second);
3485         delete (*it).second;      delete (*it).second;
3486         (*it).second = &s_mesh.vertex(i);      (*it).second = &s_mesh.vertex(i);
3487         i++;      i++;
3488     }    }
3489    
3490     i = 0;    i = 0;
3491     size_t ncell=0;    size_t ncell=0;
3492     if(nobjectToTreat==0){    if(nobjectToTreat==0){
3493         ncell+=(*(*itcell)).size();      ncell+=(*(*itcell)).size();
3494     }    }
3495   //  ffout(0)<<"ncell "<<ncell<<" n0 "<<n0<<" n1 "<<n1<<"\n";    //  ffout(0)<<"ncell "<<ncell<<" n0 "<<n0<<" n1 "<<n1<<"\n";
3496      ncell+=(*(*intertr)).size();    ncell+=(*(*intertr)).size();
3497    
3498      s_mesh.setNumberOfCells(ncell);    s_mesh.setNumberOfCells(ncell);
3499      ffout(0)<<"nb tot cells "<<s_mesh.numberOfCells()<<" \n"<<std::flush;    ffout(0)<<"nb tot cells "<<s_mesh.numberOfCells()<<" \n"<<std::flush;
3500    
3501  //   SurfaceMeshGenerator::Internals::MotherCellList::iterator motherCell    //   SurfaceMeshGenerator::Internals::MotherCellList::iterator motherCell
3502  //     = (*__internals).cellList.begin();    //     = (*__internals).cellList.begin();
3503    
3504     /* for(size_t size=0 ; size< (*(*itcell)).size() ; ++size) {    /* for(size_t size=0 ; size< (*(*itcell)).size() ; ++size) {
3505          const Triangle& T=(*(*itcell))[size];       const Triangle& T=(*(*itcell))[size];
3506          Vertex P=T(0);       Vertex P=T(0);
3507          size_t n0=listVertexFinal2.number(P);       size_t n0=listVertexFinal2.number(P);
3508          size_t n1=listVertexFinal2.number(T(1));       size_t n1=listVertexFinal2.number(T(1));
3509          size_t n2=listVertexFinal2.number(T(2));       size_t n2=listVertexFinal2.number(T(2));
3510          const Vertex& V0=*listVertexFinal[n0];       const Vertex& V0=*listVertexFinal[n0];
3511          const Vertex& V1=*listVertexFinal[n1];       const Vertex& V1=*listVertexFinal[n1];
3512          const Vertex& V2=*listVertexFinal[n2];       const Vertex& V2=*listVertexFinal[n2];
3513          s_mesh.cell(i) = Triangle(V0,V1,V2);       s_mesh.cell(i) = Triangle(V0,V1,V2);
3514  #warning SHOULD USE BC NUMBER (ie OBJECT NUMBER IN SCENE INSTEAD)       #warning SHOULD USE BC NUMBER (ie OBJECT NUMBER IN SCENE INSTEAD)
3515          s_mesh.cell(i).reference() = 1;       s_mesh.cell(i).reference() = 1;
3516          s_mesh.cell(i).setMother(&T.mother());       s_mesh.cell(i).setMother(&T.mother());
3517          i++;       i++;
3518          //motherCell++;       //motherCell++;
3519      }   */       }  */
3520       if(nobjectToTreat==0){    if(nobjectToTreat==0){
3521       SurfaceMeshGenerator::Internals::MotherCellList::iterator motherCell      SurfaceMeshGenerator::Internals::MotherCellList::iterator motherCell
3522                = (*__internals).cellList.begin();        = (*__internals).cellList.begin();
3523     for(size_t size=0 ; size< (*(*itcell)).size() ; ++size) {      for(size_t size=0 ; size< (*(*itcell)).size() ; ++size) {
3524          const Triangle& T=(*(*itcell))[size];        const Triangle& T=(*(*itcell))[size];
3525          std::map<const Vertex*,const Vertex*>::iterator        std::map<const Vertex*,const Vertex*>::iterator
         it0=(*__internals).listOfVertexMesh.find(&T(0));  
         const Vertex& V0=*((*it0).second);  
         std::map<const Vertex*,const Vertex*>::iterator  
             it1=(*__internals).listOfVertexMesh.find(&T(1));  
         const Vertex& V1=*((*it1).second);  
         std::map<const Vertex*,const Vertex*>::iterator  
             it2=(*__internals).listOfVertexMesh.find(&T(2));  
         const Vertex& V2=*((*it2).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(&T.mother());  
         i++;  
         //motherCell++;  
     }  
      }  
     for(size_t size=0 ; size< (*(*intertr)).size() ; ++size) {  
         const Triangle& T=(*(*intertr))[size];  
         std::map<const Vertex*,const Vertex*>::iterator  
3526          it0=(*__internals).listOfVertexMesh.find(&T(0));          it0=(*__internals).listOfVertexMesh.find(&T(0));
3527          const Vertex& V0=*((*it0).second);        const Vertex& V0=*((*it0).second);
3528          std::map<const Vertex*,const Vertex*>::iterator        std::map<const Vertex*,const Vertex*>::iterator
3529          it1=(*__internals).listOfVertexMesh.find(&T(1));          it1=(*__internals).listOfVertexMesh.find(&T(1));
3530          const Vertex& V1=*((*it1).second);        const Vertex& V1=*((*it1).second);
3531          std::map<const Vertex*,const Vertex*>::iterator        std::map<const Vertex*,const Vertex*>::iterator
3532          it2=(*__internals).listOfVertexMesh.find(&T(2));          it2=(*__internals).listOfVertexMesh.find(&T(2));
3533          const Vertex& V2=*((*it2).second);        const Vertex& V2=*((*it2).second);
3534          s_mesh.cell(i) = Triangle(V0,V1,V2);        s_mesh.cell(i) = Triangle(V0,V1,V2);
           
3535  #warning SHOULD USE BC NUMBER (ie OBJECT NUMBER IN SCENE INSTEAD)  #warning SHOULD USE BC NUMBER (ie OBJECT NUMBER IN SCENE INSTEAD)
3536          s_mesh.cell(i).reference() = 1;        s_mesh.cell(i).reference() = 1;
3537          s_mesh.cell(i).setMother(&T.mother());        s_mesh.cell(i).setMother(&T.mother());
3538          i++;        i++;
3539          //motherCell++;        //motherCell++;
3540      }      }
3541      }
3542    
3543      for(size_t size=0 ; size< (*(*intertr)).size() ; ++size) {
3544        const Triangle& T=(*(*intertr))[size];
3545        std::map<const Vertex*,const Vertex*>::iterator
3546          it0=(*__internals).listOfVertexMesh.find(&T(0));
3547        const Vertex& V0=*((*it0).second);
3548        std::map<const Vertex*,const Vertex*>::iterator
3549          it1=(*__internals).listOfVertexMesh.find(&T(1));
3550        const Vertex& V1=*((*it1).second);
3551        std::map<const Vertex*,const Vertex*>::iterator
3552          it2=(*__internals).listOfVertexMesh.find(&T(2));
3553        const Vertex& V2=*((*it2).second);
3554        s_mesh.cell(i) = Triangle(V0,V1,V2);
3555    
3556    #warning SHOULD USE BC NUMBER (ie OBJECT NUMBER IN SCENE INSTEAD)
3557        s_mesh.cell(i).reference() = 1;
3558        s_mesh.cell(i).setMother(&T.mother());
3559        i++;
3560        //motherCell++;
3561      }
3562            
3563            
3564  //   //  plot2(s_mesh, (*__internals).toTreatHexahedra, SMesh);  //   //  plot2(s_mesh, (*__internals).toTreatHexahedra, SMesh);
3565  //   //ffout(0) << "plot medit \n";  //   //ffout(0) << "plot medit \n";
3566     plotmedit(nobject,s_mesh, (*__internals).toTreatHexahedra,cuttedHexahedra,    plotmedit(nobject,s_mesh, (*__internals).toTreatHexahedra,cuttedHexahedra,
3567              nobjother,n0,n1);              nobjother,n0,n1);
3568  }  }
3569    
# Line 3553  SurfaceMeshGenerator::~SurfaceMeshGenera Line 3579  SurfaceMeshGenerator::~SurfaceMeshGenera
3579  }  }
3580    
3581  void plothexa(size_t ncase,std::set<const Cell*>& cuttedHexahedra) {  void plothexa(size_t ncase,std::set<const Cell*>& cuttedHexahedra) {
3582      char* cccc="hhhhhhhhhhhhh";    char* cccc="hhhhhhhhhhhhh";
3583      if(ncase==0) {    if(ncase==0) {
3584          cccc="meshed1.mesh";      cccc="meshed1.mesh";
3585      }    }
3586      else if(ncase==1) {    else if(ncase==1) {
3587          cccc="meshed2.mesh";      cccc="meshed2.mesh";
3588      }    }
3589      else {    else {
3590          cccc="meshed3.mesh";      cccc="meshed3.mesh";
3591      }    }
3592      std::ofstream sortie3(cccc);    std::ofstream sortie3(cccc);
3593      sortie3<<"MeshVersionFormatted 1"<<" \n";    sortie3<<"MeshVersionFormatted 1"<<" \n";
3594      sortie3<<" \n";    sortie3<<" \n";
3595      sortie3<<"Dimension"<<" \n";    sortie3<<"Dimension"<<" \n";
3596      sortie3<<"3"<<" \n";    sortie3<<"3"<<" \n";
3597      sortie3<<" \n";    sortie3<<" \n";
3598      sortie3<<"Vertices"<<" \n";    sortie3<<"Vertices"<<" \n";
3599      sortie3<<8*cuttedHexahedra.size()<<" \n";    sortie3<<8*cuttedHexahedra.size()<<" \n";
3600      for (std::set<const Cell*>::iterator i = cuttedHexahedra.begin();    for (std::set<const Cell*>::iterator i = cuttedHexahedra.begin();
3601           i != cuttedHexahedra.end(); ++i)         i != cuttedHexahedra.end(); ++i) {
3602      {      const CartesianHexahedron& H = static_cast<const CartesianHexahedron&>(*(*i));
3603          const CartesianHexahedron& H = static_cast<const CartesianHexahedron&>(*(*i));      for (size_t j=0; j<8; ++j)
3604          for (size_t j=0; j<8; ++j)        {
3605          {          sortie3<<H(j)[0]<<" "<<H(j)[1]<<" "<<H(j)[2]<<" 0"<<" \n";
3606              sortie3<<H(j)[0]<<" "<<H(j)[1]<<" "<<H(j)[2]<<" 0"<<" \n";        }
3607          }    }
3608      }    sortie3<<" \n";
3609      sortie3<<" \n";    sortie3<<"Hexahedra"<<"\n";
3610      sortie3<<"Hexahedra"<<"\n";    sortie3<<cuttedHexahedra.size()<<" \n";
3611      sortie3<<cuttedHexahedra.size()<<" \n";  
3612      for(size_t i=0;i<cuttedHexahedra.size();i++)
3613      for(size_t i=0;i<cuttedHexahedra.size();i++)      sortie3<<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";
3614          sortie3<<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";    sortie3<<"End"<<" \n";
     sortie3<<"End"<<" \n";  
       
3615  }  }
3616    
3617            
# Line 3597  void plotmedit(const size_t& nobject, Line 3621  void plotmedit(const size_t& nobject,
3621                 const SurfaceMeshOfTriangles& s_mesh,                 const SurfaceMeshOfTriangles& s_mesh,
3622                 std::set<const Cell*>& toTreatHexahedra,                 std::set<const Cell*>& toTreatHexahedra,
3623                 std::set<const Cell*>& cuttedHexahedra,                 std::set<const Cell*>& cuttedHexahedra,
3624                 size_t nobjother,size_t n0,size_t n1) {                 size_t nobjother,size_t n0,size_t n1)
3625    {
3626    // **********************************************************************    // **********************************************************************
3627      std::ofstream sortie3("cutted.mesh");    std::ofstream sortie3("cutted.mesh");
3628    
3629      sortie3<<"MeshVersionFormatted 1"<<" \n";    sortie3<<"MeshVersionFormatted 1"<<" \n";
3630      sortie3<<" \n";    sortie3<<" \n";
3631      sortie3<<"Dimension"<<" \n";    sortie3<<"Dimension"<<" \n";
3632      sortie3<<"3"<<" \n";    sortie3<<"3"<<" \n";
3633      sortie3<<" \n";    sortie3<<" \n";
3634      sortie3<<"Vertices"<<" \n";    sortie3<<"Vertices"<<" \n";
3635      sortie3<<8*cuttedHexahedra.size()<<" \n";    sortie3<<8*cuttedHexahedra.size()<<" \n";
3636      for (std::set<const Cell*>::iterator i = cuttedHexahedra.begin();    for (std::set<const Cell*>::iterator i = cuttedHexahedra.begin();
3637           i != cuttedHexahedra.end(); ++i)         i != cuttedHexahedra.end(); ++i) {
3638      {      const CartesianHexahedron& H = static_cast<const CartesianHexahedron&>(*(*i));
3639          const CartesianHexahedron& H = static_cast<const CartesianHexahedron&>(*(*i));      for (size_t j=0; j<8; ++j) {
3640          for (size_t j=0; j<8; ++j)        sortie3<<H(j)[0]<<" "<<H(j)[1]<<" "<<H(j)[2]<<" 0"<<" \n";
3641          {      }
3642              sortie3<<H(j)[0]<<" "<<H(j)[1]<<" "<<H(j)[2]<<" 0"<<" \n";    }
3643          }    sortie3<<" \n";
3644      }    sortie3<<"Hexahedra"<<"\n";
3645      sortie3<<" \n";    sortie3<<cuttedHexahedra.size()<<" \n";
3646      sortie3<<"Hexahedra"<<"\n";  
3647      sortie3<<cuttedHexahedra.size()<<" \n";    for(size_t i=0;i<cuttedHexahedra.size();i++)
3648        sortie3<<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";
3649      for(size_t i=0;i<cuttedHexahedra.size();i++)    sortie3<<"End"<<" \n";
         sortie3<<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";  
     sortie3<<"End"<<" \n";  
       
3650    
       
3651    std::ofstream sortie("intersection.mesh");    std::ofstream sortie("intersection.mesh");
3652              
3653    sortie<<"MeshVersionFormatted 1"<<" \n";    sortie<<"MeshVersionFormatted 1"<<" \n";
3654    sortie<<" \n";    sortie<<" \n";
3655    sortie<<"Dimension"<<" \n";    sortie<<"Dimension"<<" \n";
# Line 3638  void plotmedit(const size_t& nobject, Line 3658  void plotmedit(const size_t& nobject,
3658    sortie<<"Vertices"<<" \n";    sortie<<"Vertices"<<" \n";
3659    sortie<<8*toTreatHexahedra.size()/*+3*s_mesh.numberOfCells()*/+s_mesh.numberOfVertices()<<" \n";    sortie<<8*toTreatHexahedra.size()/*+3*s_mesh.numberOfCells()*/+s_mesh.numberOfVertices()<<" \n";
3660    for (std::set<const Cell*>::iterator i = toTreatHexahedra.begin();    for (std::set<const Cell*>::iterator i = toTreatHexahedra.begin();
3661         i != toTreatHexahedra.end(); ++i)         i != toTreatHexahedra.end(); ++i) {
3662      {      const CartesianHexahedron& H = static_cast<const CartesianHexahedron&>(*(*i));
3663        const CartesianHexahedron& H = static_cast<const CartesianHexahedron&>(*(*i));      for (size_t j=0; j<8; ++j) {
3664        for (size_t j=0; j<8; ++j)        sortie<<H(j)[0]<<" "<<H(j)[1]<<" "<<H(j)[2]<<" 0"<<" \n";
         {  
           sortie<<H(j)[0]<<" "<<H(j)[1]<<" "<<H(j)[2]<<" 0"<<" \n";  
         }  
     }  
   for(size_t i=0;i<s_mesh.numberOfVertices();++i)  
     {  
       sortie<<s_mesh.vertex(i)[0]<<" "<<s_mesh.vertex(i)[1]<<" "<<s_mesh.vertex(i)[2]<<" 0"<<" \n";  
3665      }      }
3666      }
3667      for(size_t i=0;i<s_mesh.numberOfVertices();++i) {
3668        sortie<<s_mesh.vertex(i)[0]<<" "<<s_mesh.vertex(i)[1]<<" "<<s_mesh.vertex(i)[2]<<" 0"<<" \n";
3669      }
3670    
3671       /*  for (size_t i=0; i<s_mesh.numberOfCells(); ++i)    /*  for (size_t i=0; i<s_mesh.numberOfCells(); ++i)
3672              {        {
3673              const Triangle& T = s_mesh.cell(i);        const Triangle& T = s_mesh.cell(i);
3674              sortie<<T(0)[0]<<" "<<T(0)[1]<<" "<<T(0)[2]<<" 0"<<" \n";        sortie<<T(0)[0]<<" "<<T(0)[1]<<" "<<T(0)[2]<<" 0"<<" \n";
3675              sortie<<T(1)[0]<<" "<<T(1)[1]<<" "<<T(1)[2]<<" 0"<<" \n";        sortie<<T(1)[0]<<" "<<T(1)[1]<<" "<<T(1)[2]<<" 0"<<" \n";
3676              sortie<<T(2)[0]<<" "<<T(2)[1]<<" "<<T(2)[2]<<" 0"<<" \n";        sortie<<T(2)[0]<<" "<<T(2)[1]<<" "<<T(2)[2]<<" 0"<<" \n";
3677            }*/
             }*/  
3678    sortie<<" \n";    sortie<<" \n";
3679    sortie<<"Triangles"<<"\n";    sortie<<"Triangles"<<"\n";
3680    sortie<<s_mesh.numberOfCells()<<"\n";    sortie<<s_mesh.numberOfCells()<<"\n";
3681    ffout(0)<<"fin ecriture vertex\n";    ffout(0)<<"fin ecriture vertex\n";
3682    for (size_t i=0; i<s_mesh.numberOfCells(); ++i)    for (size_t i=0; i<s_mesh.numberOfCells(); ++i) {
3683      {      const Triangle& T = s_mesh.cell(i);
3684        const Triangle& T = s_mesh.cell(i);      // ffout(0)<<s_mesh.vertexNumber(T(0))<<" "<<s_mesh.vertexNumber(T(1))<<" "<<s_mesh.//vertexNumber(T(2))<<"\n";
3685       // ffout(0)<<s_mesh.vertexNumber(T(0))<<" "<<s_mesh.vertexNumber(T(1))<<" "<<s_mesh.//vertexNumber(T(2))<<"\n";      sortie<<s_mesh.vertexNumber(T(0))+8*toTreatHexahedra.size()+1<<" "<<s_mesh.vertexNumber(T(1))+8*toTreatHexahedra.size()+1<<" "<<s_mesh.vertexNumber(T(2))+8*toTreatHexahedra.size()+1;
3686        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;      
3687              // sortie<<3*i+1+8*toTreatHexahedra.size()<<" "<<3*i+2+8*toTreatHexahedra.size()<<" //"<<3*i+3+8*toTreatHexahedra.size();//<<" "<<1<<"\n";
3688       // sortie<<3*i+1+8*toTreatHexahedra.size()<<" "<<3*i+2+8*toTreatHexahedra.size()<<" //"<<3*i+3+8*toTreatHexahedra.size();//<<" "<<1<<"\n";      if(i<n0)
3689        if(i<n0)        sortie<<" 1"<<"\n";
3690            sortie<<" 1"<<"\n";      else if(i<n1)
3691        else if(i<n1)        sortie<<" 2"<<"\n";
3692            sortie<<" 2"<<"\n";      /* else if(i<(nobjother+ninter))
3693       /* else if(i<(nobjother+ninter))         sortie<<" 3"<<"\n";
3694            sortie<<" 3"<<"\n";         else if(i<(nobjother+ninter+ninter2))
3695        else if(i<(nobjother+ninter+ninter2))         sortie<<" 4"<<"\n";*/
3696            sortie<<" 4"<<"\n";*/      else
3697        else        sortie<<" 4"<<"\n";
3698            sortie<<" 4"<<"\n";    }
     }  
3699    ffout(0)<<"fin ecriture triangles\n";    ffout(0)<<"fin ecriture triangles\n";
3700    sortie<<" \n";    sortie<<" \n";
3701    sortie<<"Hexahedra"<<"\n";    sortie<<"Hexahedra"<<"\n";
# Line 3690  void plotmedit(const size_t& nobject, Line 3705  void plotmedit(const size_t& nobject,
3705      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";
3706    sortie<<"End"<<" \n";    sortie<<"End"<<" \n";
3707    ffout(0)<<"fin ecriture intersection.mesh\n";    ffout(0)<<"fin ecriture intersection.mesh\n";
   
3708  }  }
3709    
3710    

Legend:
Removed from v.1.15  
changed lines
  Added in v.1.16

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