/[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.10 by delpinux, Tue Mar 25 20:09:57 2003 UTC revision 1.11 by dobj, Wed Jun 18 10:26:41 2003 UTC
# Line 55  Line 55 
55  class SurfaceMeshGenerator::Internals  class SurfaceMeshGenerator::Internals
56  {  {
57  private:  private:
58    friend class SurfaceMeshGenerator;      class IntersectionPoints;
59        class TriangleCut;
60        friend class SurfaceMeshGenerator;
61    
62    struct __ltv    struct __ltv
63    {    {
# Line 71  private: Line 73  private:
73      }      }
74    };    };
75    
76      struct __ltv2
77      {
78          bool operator()(const std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> >& s1,
79                          const std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> >& s2) const
80      {
81          if(s1.first[0].first==s2.first[0].first) {
82              if(s1.first[0].second==s2.first[0].second) {
83                  if(s1.first[1].first==s2.first[1].first) {
84                      return s1.first[1].second<s2.first[1].second;
85                  }
86                  else {
87                      return s1.first[1].first<s2.first[1].first;
88                  }
89              }
90              else {
91                  return s1.first[0].second<s2.first[0].second;
92              }
93          }
94          else {
95              return s1.first[0].first<s2.first[0].first;
96          }
97      }
98      };
99    
100    typedef std::map<TinyVector<3,real_t>,    typedef std::map<TinyVector<3,real_t>,
101                     std::vector<SurfaceMesh*>, __ltv> VRRefMesh;                     std::vector<SurfaceMesh*>, __ltv> VRRefMesh;
102    
103    typedef std::list<TinyVector<3,Edge::Pair> > EdgePairList;    typedef std::list<TinyVector<3,Edge::Pair> > EdgePairList;
104      typedef std::list<TinyVector<3,Edge::Pair> * > EdgePairList2;
105      typedef std::pair< const TinyVector<3,Edge::Pair>*,const TinyVector<3,Edge::Pair> *> PairTriangle;
106      typedef std::pair< const TinyVector<3,TinyVector<2,Edge::Pair> >*,
107          const TinyVector<3,TinyVector<2,Edge::Pair> > *> PairTriangleEdges;
108    
109      struct __ltvpair
110      {
111        bool operator()(const  PairTriangle & s1, const PairTriangle & s2) const
112        {
113          if(s1.first==s2.first) {
114            return s1.second<s2.second;
115          }
116          else {
117            return s1.first<s2.first;
118          }
119        }
120      };
121    
122      struct __ltvpairedges
123      {
124          bool operator()(const  PairTriangleEdges & s1, const PairTriangleEdges & s2) const
125      {
126          if(s1.first==s2.first) {
127              return s1.second<s2.second;
128          }
129          else {
130              return s1.first<s2.first;
131          }
132      }
133      };
134        
135    typedef std::list<const Cell*> MotherCellList;    typedef std::list<const Cell*> MotherCellList;
136    
137    typedef std::map<const Cell*,    typedef std::map<const Cell*,
# Line 83  private: Line 140  private:
140    HexaTetraAssociation;    HexaTetraAssociation;
141    
142    typedef std::map<Edge::Pair, Vertex*> VerticesList;    typedef std::map<Edge::Pair, Vertex*> VerticesList;
143      typedef std::vector< TinyVector<3,size_t> > TriangleVerticesNumber;
144      typedef std::map<TinyVector<3>, Object*, SurfaceMeshGenerator::Internals::__ltv> ObjectReferences;
145      typedef std::pair<size_t,size_t> ObjectPair;
146    
   
     
147    //! list of cut edges and intersection vertex.    //! list of cut edges and intersection vertex.
148    std::map<Edge::Pair, size_t> edgesRef;    std::vector<ReferenceCounting<std::map<Edge::Pair, size_t> > > edgesRef;
149      std::vector<ReferenceCounting<std::map<Edge::Pair, size_t> > > edgesRefCreate;
150    
151    std::vector<ReferenceCounting<EdgePairList> > triangleListes;    typedef  TinyVector<3,TinyVector<2,Edge::Pair> >  TriangleVertices;
152      typedef std::list< TriangleVertices > TriangleVerticesEdgesList;
153      typedef std::list< TriangleVertices *> TriangleVerticesEdgesList2;
154      std::vector<ReferenceCounting<TriangleVerticesEdgesList> > triangleListesIntersectionEdges;
155      typedef std::map<TriangleVertices*,const TriangleVertices*> TriangleVerticesMap;
156      std::vector<ReferenceCounting< TriangleVerticesMap> > triangleListesIntersectionEdgesTemp;
157    
158      typedef std::vector< TriangleVertices > TriangleVerticesEdges;
159      typedef std::map<const TriangleVertices*,TriangleVerticesEdges > TriangleTreatMap;
160      TinyVector<2, TriangleTreatMap > triangleTreat;
161    
162      std::vector<ReferenceCounting<EdgePairList> > triangleListes;
163      TinyVector<2,std::vector< TriangleVerticesMap::iterator >  > triangleEraseEdgesTemp;
164      std::vector<ReferenceCounting<TriangleVerticesEdgesList> > localTriangleListesEdges;
165      std::vector<ReferenceCounting<TriangleVerticesEdgesList2> > localTriangleList2Edges;
166    EdgePairList localTriangleList;    EdgePairList localTriangleList;
167    EdgePairList localTriangleListother;  
168    
169    //! This list is used to provide an orientation to the Surfacic Elements.    //! This list is used to provide an orientation to the Surfacic Elements.
170    std::list<TinyVector<3> > InsideVertexList;    std::list<TinyVector<3> > InsideVertexList;
171    
172      //For the intersection
173      std::map< ObjectPair , std::map< PairTriangleEdges , IntersectionPoints, __ltvpairedges > >
174          listVertexIntersectionPerPairEdges;
175    
176      typedef std::map<std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> >,Vertex*,__ltv2>            ListVerticesIntersectionPerEdges;
177      ListVerticesIntersectionPerEdges   listVerticesIntersectionPerEdges;
178      
179      
180      std::vector<ReferenceCounting< TriangleVerticesEdges > > listTriangleintersection;
181      TriangleVerticesEdges  localListTriangleintersection;
182      
183    
184    //! list of cut edges and intersection vertex.    //! list of cut edges and intersection vertex.
  //VerticesList verticesList;  
185    std::vector<ReferenceCounting<VerticesList> > verticesListes;    std::vector<ReferenceCounting<VerticesList> > verticesListes;
186      std::vector<ReferenceCounting<VerticesList> > verticesListesIntersection;
187    
188    HexaTetraAssociation hexaTetraAssociation;    HexaTetraAssociation hexaTetraAssociation;
189      std::set<const Cell*> toTreatHexahedra;
190    
191    
192    MotherCellList cellList;    MotherCellList cellList;
193    
# Line 108  private: Line 195  private:
195                              const Reference::Representation& ref,                              const Reference::Representation& ref,
196                              const Shape& S);                              const Shape& S);
197    
198    void __createTriangles(const Shape& S,    void __createListIntersection(const Shape& S,
199                           const Reference::Representation& ref,const int objectNumber);                                  std::vector<const Edge*>& cutEdges,
200                                    const Reference::Representation& ref,
201                                    const size_t& objectNumber,
202                                    const CartesianHexahedron& currentCell);
203    
204      void __createCutEdges(const Scene& scene,
205                            std::vector<ReferenceCounting< std::vector<const Edge*> > > & cutEdges,
206                            ObjectReferences& otherReferences,
207                            const std::vector<TetrahedronByEdges>::iterator& currentT);
208    
209      void __calculatePointsIntersection(const Scene& scene,
210                                         ObjectReferences& otherReferences);
211    
212    void __setBoundaryType(const Boundary& B,    void __setBoundaryType(const Boundary& B,
213                           POVRayReferences& vrRefList,                           POVRayReferences& vrRefList,
214                           TinyVector<6,bool>& faces);                           TinyVector<6,bool>& faces);
215    
216      size_t __IsInTriangle(const TinyVector<2,size_t> & num,
217                            const TinyVector<3,TinyVector<3> > & triangleVector );
218      
219    void __InTriangle(const Vertex& P,    void __InTriangle(const Vertex& P,
220              TinyVector<3>& P0,                      const TinyVector<3, Vertex* > & V,
221              TinyVector<3>& V12,                      size_t& ncompt);
222              TinyVector<3>& V02,  
223              size_t& ncompt);    void __inout(const size_t& numobject,
224    void __inout(size_t iter,                 const int& ttn0,
225               const size_t nn,                 const TinyVector<2,Edge::Pair>& e0,
226               int& ttn0,                 const TinyVector<2,Edge::Pair>& e1,
227                 int& ttn1,                 size_t& ii0,
228                 int& ttn2,                 size_t& ii1,TriangleCut& K);
229               size_t& re0,  
230               size_t& re1,    TinyVector<4, real_t>  __EquationPlan(const TinyVector<3, Vertex* > & V);
231               size_t& re2,    void __DeterminatePoint(const size_t& ncase,
232               size_t& ii0,                            const Edge::Pair & e0,
233               size_t& ii1,                            const Edge::Pair & e1,
234               size_t& ii2);                            const Vertex & V0,
235                              const Vertex & V1,
236                              const TinyVector<3, Vertex* > & V,
237                              IntersectionPoints& Pinter);
238      void __CalculPoint(const size_t& objectNumber,
239                         const size_t& numobject,
240                         SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList2::const_iterator & currentTriangle,
241                         SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList2::const_iterator & currentTriangle1,
242                         IntersectionPoints & Pinter);
243    
244      void __createPointIntersection(const size_t& numobject,
245                                     const size_t& objectNumber);
246        
247    void __CalculPoint(const Vertex& V0, const Vertex& V1 , const Vertex& V2,    std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > __determinateNumber(const int& p,
248                       const Vertex& Vo0, const Vertex& Vo1 , const Vertex& Vo2,                                                    const TinyVector<3,std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > >& pointsint);
249                       int nn,int& ttn0,int& ttn1 ,int& ttn2,    void __determinateCase(const int& ncase,
250                       Vertex *& vvn0 , Vertex*& vvn1 , Vertex*& vvn2);                           const TriangleCut & K,
251                             Vertex & A,Vertex & B,Vertex & C,
252    void __IntersectionPoint(const Scene& scene ,const Reference::Representation& ref,                           TinyVector<2,Edge::Pair> & e0,TinyVector<2,Edge::Pair> & e1,
253                             const Reference::Representation& otherRef,                           TinyVector<2,Edge::Pair> & e2,
254                             std::map<Edge::Pair,int>& triangleedge,                           TinyVector<2,int>& tt,
255                             std::map<Edge::Pair, Vertex*>&  listofVertex);                           TinyVector<2,int>& pp);
256      void __addTriangle(const size_t& numobject,
257                            const TinyVector<2,Edge::Pair>& numA,
258                            const TinyVector<2,Edge::Pair>& numB,
259                            const TinyVector<2,Edge::Pair>& numC,
260                            size_t& kill);
261      void __createCase1(int ncase,
262                         const TriangleCut& K,
263                         size_t & kill);
264    
265      void __createCase2(int ncase,
266                         const TriangleCut& K,
267                         size_t & kill);
268    
269      void __createLocalListIntersection(const size_t& objectNumber,
270                                         const TinyVector<3,TinyVector<2,Edge::Pair> >& currentTriangle,
271                                         TriangleCut & K);
272        
273      TinyVector<3,int> __determinateEdgeCut(const size_t& numobject,
274                                             const size_t & objectNumber,
275                                             const TinyVector<3,TinyVector<2,Edge::Pair> >& currentTriangle);
276    
277      TinyVector<3,size_t> __putEdgeRef(const size_t& ncase,
278                                        const TinyVector<3,TinyVector<2,Edge::Pair> >& currentTriangle,
279                                        TriangleCut& K);
280    
281      void __createListTriangle(const TriangleCut& K,
282                                size_t & kill);
283    
284      void __createTrianglesIntersection(const size_t& numobject,
285                                         const size_t& objectNumber,
286                                         const size_t& nobject);
287    
288      bool __createTriangleSurface(const size_t& nobject,
289                                   const Shape& S ,
290                                   const Domain& Omega,
291                                   std::vector<const Edge*>& cutEdges,
292                                   const Reference::Representation& ref,
293                                   const Cell& currentCell);
294      void __createSurface(Structured3DMesh& SMesh,
295                           const Domain& Omega,
296                           const Shape& S ,
297                           const Reference::Representation& ref,
298                           size_t nobject,
299                           std::set<const Cell*>& cuttedHexahedra,
300                           std::set<const Cell*>& meshedHexahedra);
301    
302      /// @todo FUNCTION TO REMOVE
303      static void __getIntersectionReferences(const Intersection& I,
304                                              std::map<TinyVector<3>, Object*, Internals::__ltv>& otherReferences);
305      void plotliste(size_t nobject);
306    
307      void printpoint(size_t numobject,const TinyVector<2,Edge::Pair>& e0 , TriangleCut K);
308      bool triangleplat(size_t numobject,const TinyVector<2,Edge::Pair>& e0 ,
309                        const TinyVector<2,Edge::Pair>& e1,
310                        const TinyVector<2,Edge::Pair>& e2 );
311    
312  public:  public:
313    Internals()    Internals()
314    {    {
# Line 148  public: Line 316  public:
316    }    }
317  };  };
318    
319    class SurfaceMeshGenerator::Internals::IntersectionPoints
320    {
321    public :
322        TinyVector<2,size_t> numpair;
323        std::vector< std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > > listofPoints;
324        TinyVector< 3,size_t > cutTriangle1;
325        TinyVector< 3,size_t > cutTriangle2;
326    
327        inline const size_t size() const {
328            return listofPoints.size();
329        }
330    
331        void addPoint(std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> >&  V) {
332            listofPoints.push_back( V );
333        }
334      
335    
336      //constructor
337        IntersectionPoints() {
338            numpair=0;
339            cutTriangle1=0;
340            cutTriangle2=0;
341        }
342    
343      IntersectionPoints(const IntersectionPoints& P) {
344          numpair=P.numpair;
345          listofPoints=P.listofPoints;
346          cutTriangle1=P.cutTriangle1;
347          cutTriangle2=P.cutTriangle2;
348      }
349    
350      //! Destructor.
351      ~IntersectionPoints()
352      {
353        ;
354      }
355        
356    };
357    
358    class SurfaceMeshGenerator::Internals::TriangleCut
359    {
360    public :
361        TinyVector<3,TinyVector<2,Edge::Pair> > listpairedges;
362        TinyVector<3,int> edgecut;
363        TinyVector<3,size_t> isInTriangle;
364        TinyVector<3,std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > > pointsIntersection;
365        std::vector< std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > > pointsin;
366        size_t numobject;
367        size_t otherobject;
368    
369        //constructor
370        TriangleCut() {
371        }
372        TriangleCut(const TriangleCut& K) {
373            listpairedges=K.listpairedges;
374            edgecut=K.edgecut;
375            isInTriangle=K.isInTriangle;
376            pointsIntersection=K.pointsIntersection;
377            pointsin=K.pointsin;
378            numobject=K.numobject;
379            otherobject=K.otherobject;
380        }
381        TriangleCut(const TinyVector<2,Edge::Pair> & e0, const TinyVector<2,Edge::Pair> & e1,
382                    const TinyVector<2,Edge::Pair> & e2,
383                    size_t num,size_t num2) {
384            listpairedges[0]=e0;
385            listpairedges[1]=e1;
386            listpairedges[2]=e2;
387            numobject=num;
388            otherobject=num2;
389            edgecut=0;
390            isInTriangle=0;
391        }
392        
393        //! Destructor.
394        ~TriangleCut()
395        {
396            ;
397        }
398    };
399    
400  void  void
401  SurfaceMeshGenerator::Internals::  SurfaceMeshGenerator::Internals::
402  __setBoundaryType(const Boundary& B,  __setBoundaryType(const Boundary& B,
# Line 202  __splitEdge(const Edge& e, Line 451  __splitEdge(const Edge& e,
451    TinyVector<3> vo = e(out);    TinyVector<3> vo = e(out);
452    TinyVector<3> vs = (vi+vo)/2;    TinyVector<3> vs = (vi+vo)/2;
453    
454    for (int l = 0; l<5; l++) { // 3% error after 5 iterations.    for (int l = 0; l<10; l++) { // 3% error after 5 iterations.
455      if(S.inside(vs)) {      if(S.inside(vs)) {
456        vi = vs;        vi = vs;
457        vs = (vo+vi)/2;        vs = (vo+vi)/2;
# Line 211  __splitEdge(const Edge& e, Line 460  __splitEdge(const Edge& e,
460        vs = (vo+vi)/2;        vs = (vo+vi)/2;
461      }      }
462    }    }
463      if(!S.inside(vs))
464    return vs;        return vs;
465      else
466          return vo;
467  }  }
468    
469  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
470  __createTriangles(const Shape& S,  __createListIntersection(const Shape& S,
471                    const Reference::Representation& ref,const int objectNumber)                           std::vector<const Edge*>& cutEdges,
472  {                           const Reference::Representation& ref,
473                                 const size_t& objectNumber,
474    for(SurfaceMeshGenerator::Internals::HexaTetraAssociation::iterator i                           const CartesianHexahedron& currentCell) {
475          = hexaTetraAssociation.begin();  
476        i != hexaTetraAssociation.end(); ++i) {      switch (cutEdges.size()) {
477                case 3: {
478      const CartesianHexahedron& currentCell              std::vector<TinyVector<3> > V(cutEdges.size());
479        = static_cast<const CartesianHexahedron&>(*(*i).first);              for(size_t i=0; i<cutEdges.size(); ++i) {
480      std::vector<TetrahedronByEdges>& T = *(((*i).second).first);                  const Edge& e = *cutEdges[i];
481                          if(e(0).Ref().isSet(ref))
482      for (std::vector<TetrahedronByEdges>::iterator currentT = T.begin();                      V[i] = (e(0).Vector3() - e(1).Vector3());
483           currentT < T.end(); ++currentT) {                  else            
484        std::vector<const Edge*> cutEdges;                      V[i] = (e(1).Vector3() - e(0).Vector3());
485                }
486    
487                if ((V[0]^V[1])*V[2] < 0)
488                std::swap(cutEdges[1], cutEdges[2]);
489    
490                //! Addes a triangle to the list.
491                (*localTriangleListesEdges[objectNumber]).push_front(TinyVector<3,TinyVector<2,Edge::Pair> >());
492                cellList.push_front(&currentCell);
493    
494                //! Sets the triangle vertices (actually just store the edge where it lives).
495                TriangleVerticesEdgesList::iterator currentTriangleEdges
496                    = (*localTriangleListesEdges[objectNumber]).begin();
497    
498                for(size_t i=0; i<cutEdges.size(); ++i) {
499                    Edge::Pair e=(*cutEdges[i]);
500                    TinyVector<2,Edge::Pair> pe;
501                    pe[0]=e;
502                    pe[1]=e;
503                    (*currentTriangleEdges)[i] = pe;
504                }
505    
506                //! Stores the Edges on which the vertex will be computed.
507                for (size_t i=0; i<cutEdges.size(); ++i)
508                    (*verticesListesIntersection[objectNumber])[*cutEdges[i]] =
509                        new Vertex(Internals::__splitEdge(*cutEdges[i], ref, S));
510    
511        for(size_t ie=0; ie<6; ++ie) {              break;
         const Edge& e = currentT->edge(ie);  
         if ((e(0).Ref().isSet(ref)) xor (e(1).Ref().isSet(ref))) {  
           cutEdges.push_back(&e);  
512          }          }
513        }          case 4: {
514                //! Orders the Edges using a map !
515                std::map<Edge::Pair, const Edge*> cutEdgesMap;
516                for(size_t ie=0; ie<cutEdges.size(); ++ie)
517                    cutEdgesMap[(Edge::Pair)*(cutEdges[ie])] = (cutEdges[ie]);
518    
519        switch (cutEdges.size()) {              //! Addes a triangle to the list.
520        case 3: {              (*localTriangleListesEdges[objectNumber]).push_front(TinyVector<3,TinyVector<2,Edge::Pair> >());
         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());  
         }  
521    
522          if ((V[0]^V[1])*V[2] < 0)              cellList.push_front(&currentCell);
           std::swap(cutEdges[1], cutEdges[2]);  
523    
524          //! Addes a triangle to the list.              std::vector<TinyVector<3> > V(cutEdges.size());
525          localTriangleList.push_front(TinyVector<3,Edge::Pair>());              std::map<Edge::Pair, const Edge*>::iterator currentCutEdge=cutEdgesMap.begin();
526          cellList.push_front(&currentCell);              for (size_t i=0; i<V.size(); ++i) {
527                    const Edge& e = *(*currentCutEdge).second;
528                    if(e(0).Ref().isSet(ref))
529                        V[i] = (e(0).Vector3() - e(1).Vector3());
530                    else
531                        V[i] = (e(1).Vector3() - e(0).Vector3());
532                    ++currentCutEdge;
533                }
534    
535          //! Sets the triangle vertices (actually just store the edge where it lives).              bool swapEdge = ((V[0]^V[1])*V[2] > 0);
         EdgePairList::iterator currentTriangle = localTriangleList.begin();  
         for(size_t i=0; i<cutEdges.size(); ++i)  
           (*currentTriangle)[i] = (Edge::Pair)(*cutEdges[i]);  
536    
537          //! Stores the Edges on which the vertex will be computed.              //! Sets the triangle vertices (actually just store the edge where it lives).
538          for (size_t i=0; i<cutEdges.size(); ++i)              TriangleVerticesEdgesList::iterator currentTriangleEdges
539          //verticesList[*cutEdges[i]] =  new Vertex(Internals::__splitEdge(*cutEdges[i], ref, S));                  = (*localTriangleListesEdges[objectNumber]).begin();
540              (*verticesListes[objectNumber])[*cutEdges[i]] =  new Vertex(Internals::__splitEdge(*cutEdges[i], ref, S));                      
541                currentCutEdge = cutEdgesMap.begin();
542                if (!swapEdge) {
543                    for(size_t i=0; i<3; ++i) {
544                        Edge::Pair e=(currentCutEdge->first);
545                        TinyVector<2,Edge::Pair> pe;
546                        pe[0]=e;
547                        pe[1]=e;
548                        (*currentTriangleEdges)[i] = pe;
549                        ++currentCutEdge;
550                    }
551                }
552                else {
553                    for(int i=2; i>=0; --i) {
554                        Edge::Pair e=(currentCutEdge->first);
555                        TinyVector<2,Edge::Pair> pe;
556                        pe[0]=e;
557                        pe[1]=e;
558                        (*currentTriangleEdges)[i] = pe;
559                        ++currentCutEdge;
560                    }
561                }
562                for (size_t i=0; i<4; ++i) {
563                    (*verticesListesIntersection[objectNumber])[*cutEdges[i]] =
564                    new Vertex(__splitEdge(*cutEdges[i], ref, S));
565    
566          break;              }
       }  
       case 4: {  
         //! 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;  
         }  
567    
568          bool swapEdge = ((V[0]^V[1])*V[2] > 0);              //! Addes a triangle to the list.
569                (*localTriangleListesEdges[objectNumber]).push_front(TinyVector<3,TinyVector<2,Edge::Pair> >());
570    
571          //! Sets the triangle vertices (actually just store the edge where it lives).              cellList.push_front(&currentCell);
         EdgePairList::iterator currentTriangle = localTriangleList.begin();  
572    
573          currentCutEdge = cutEdgesMap.begin();              currentTriangleEdges
574          if (!swapEdge)                  = (*localTriangleListesEdges[objectNumber]).begin();
575            for(size_t i=0; i<3; ++i) {              
576              (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);              //! Sets the triangle vertices (actually just store the edge where it lives).
577              ++currentCutEdge;              currentCutEdge = cutEdgesMap.begin();
578            }              if (!swapEdge) {
579          else                  for(int i=2; i>=0; --i) {
580            for(int i=2; i>=0; --i) {                      ++currentCutEdge;
581              (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);                    //  (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
582              ++currentCutEdge;                      Edge::Pair e=(currentCutEdge->first);
583            }                      TinyVector<2,Edge::Pair> pe;
584                        pe[0]=e;
585                        pe[1]=e;
586                        (*currentTriangleEdges)[i] = pe;
587                    }
588                }
589                else {
590                    for(size_t i=0; i<3; ++i) {
591                        ++currentCutEdge;
592                       // (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
593                        Edge::Pair e=(currentCutEdge->first);
594                        TinyVector<2,Edge::Pair> pe;
595                        pe[0]=e;
596                        pe[1]=e;
597                        (*currentTriangleEdges)[i] = pe;
598                    }
599                }
600                break;
601            }
602            case 1:
603            case 2: {
604                ffout(3) << "?!\n";
605                break;
606            }
607        }
608    }
609    void SurfaceMeshGenerator::Internals::
610    __createCutEdges(const Scene& scene,
611                     std::vector<ReferenceCounting< std::vector<const Edge*> > > & cutEdges,
612                     ObjectReferences& otherReferences,
613                     const std::vector<TetrahedronByEdges>::iterator& currentT) {
614    
615          for (size_t i=0; i<4; ++i) {      //size_t nobjectToTreat = otherReferences.size();
           //verticesList[*cutEdges[i]] =  new Vertex(__splitEdge(*cutEdges[i], ref, S));  
             (*verticesListes[objectNumber])[*cutEdges[i]] =  new Vertex(__splitEdge(*cutEdges[i], ref, S));  
616    
617        for(size_t ie=0; ie<6; ++ie) {
618            const Edge& e = currentT->edge(ie);
619            int num=0;
620            for(ObjectReferences::iterator i = otherReferences.begin();
621                         i != otherReferences.end(); ++i) {
622                for(size_t n = 0; n<scene.nbObjects(); ++n) {
623                    Reference::Representation otherRef;
624                    otherRef.resize(scene.nbObjects());
625                    if (scene.object(n).reference() == ((*i).first)) {
626                        otherRef[n] = true;
627                        if ((e(0).Ref().isSet(otherRef)) xor (e(1).Ref().isSet(otherRef))) {
628                            (*cutEdges[num]).push_back(&e);
629                        }
630                    }
631                }
632                num++;
633          }          }
634        }
635    }
636    
637          //! Addes a triangle to the list.  void SurfaceMeshGenerator::Internals::
638          localTriangleList.push_front(TinyVector<3,Edge::Pair>());  __calculatePointsIntersection(const Scene& scene,
639          cellList.push_front(&currentCell);                                ObjectReferences& otherReferences) {
640        size_t nobject = scene.nbObjects();
641        size_t nobjectToTreat = otherReferences.size();
642    
643        for(std::set<const Cell*>::iterator i=toTreatHexahedra.begin() ; i!=toTreatHexahedra.end() ; ++i)
644        {
645            SurfaceMeshGenerator::Internals::HexaTetraAssociation::iterator
646                     j=hexaTetraAssociation.find(*i);
647            std::vector<TetrahedronByEdges>& T = *(((*j).second).first);
648            const CartesianHexahedron& currentCell
649                = static_cast<const CartesianHexahedron&>(*(*i));
650    
651            for (std::vector<TetrahedronByEdges>::iterator currentT = T.begin();
652          currentT < T.end(); ++currentT) {
653    
654                std::vector<ReferenceCounting< std::vector<const Edge*> > > cutEdges;
655                cutEdges.resize(nobject);
656                for(size_t no=0 ; no<nobject ; ++no) {
657                    cutEdges[no]= new std::vector<const Edge*>;
658                }
659    
660                __createCutEdges(scene,cutEdges,otherReferences,currentT);
661    
662                int num=0;
663                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 (TriangleVerticesEdgesList::iterator t = (*localTriangleListesEdges[no]).begin();
680              t != (*localTriangleListesEdges[no]).end(); ++t) {
681                    (*triangleListesIntersectionEdges[no]).push_front(*t);
682                    (*localTriangleList2Edges[no]).push_front(&(*(*triangleListesIntersectionEdges[no]).begin()));
683                }
684            }
685    
686            size_t intersec=0;
687            for(size_t no=0 ; no<nobjectToTreat ; ++no) {
688                for(size_t no2=no+1 ; no2<nobject ; ++no2) {
689                    ObjectPair objectp(no,no2);
690                    size_t sizetemp=(listVertexIntersectionPerPairEdges[objectp]).size();
691                    __createPointIntersection(no,no2);
692                    if(sizetemp<(listVertexIntersectionPerPairEdges[objectp]).size())
693                        intersec=1;
694                }
695            }
696            if(intersec==0) {
697               // ffout(0)<<"pas d'intersection dans l hexa \n";
698            }
699            for(size_t no=0 ; no<nobject ; ++no) {
700                (*localTriangleListesEdges[no]).clear();
701                (*localTriangleList2Edges[no]).clear();
702            }
703        }
704    
705          currentTriangle = localTriangleList.begin();  }
706    
707          //! Sets the triangle vertices (actually just store the edge where it lives).  void SurfaceMeshGenerator::Internals::printpoint(size_t numobject,
708          currentCutEdge = cutEdgesMap.begin();                                                   const TinyVector<2,Edge::Pair>& e0 , TriangleCut K) {
709          if (!swapEdge)      Vertex V0;
710            for(int i=2; i>=0; --i) {      if(e0[0]==e0[1])
711              ++currentCutEdge;          V0=(*(*verticesListesIntersection[numobject])[e0[0]]);
712              (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);      else {
713            }          std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > p0;
714            p0.first=e0;
715            TinyVector<2,size_t> num;
716            num[0]=K.numobject;
717            num[1]=K.otherobject;
718            p0.second=num;
719            SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator
720                ip0=listVerticesIntersectionPerEdges.find(p0);
721            if(ip0!=listVerticesIntersectionPerEdges.end())
722                V0 =*(*ip0).second;
723          else          else
724            for(size_t i=0; i<3; ++i) {              ffout(0)<<"point exi pas\n";
             ++currentCutEdge;  
             (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);  
           }  
         break;  
       }  
       case 1:  
       case 2: {  
         ffout(3) << "?!\n";  
         break;  
       }  
       }  
725      }      }
726        ffout(0)<<V0.Vector3()[0]<<" "<<V0.Vector3()[1]<<" "<<V0.Vector3()[2]<<"\n";
727    }
728    
729      for (EdgePairList::iterator t = localTriangleList.begin();  bool SurfaceMeshGenerator::Internals::triangleplat(size_t numobject,const TinyVector<2,Edge::Pair>& e0 ,
730           t != localTriangleList.end(); ++t) {                                                     const TinyVector<2,Edge::Pair>& e1,
731          (*triangleListes[objectNumber]).push_front(*t);                                                     const TinyVector<2,Edge::Pair>& e2 ) {
732    
733    
734        Vertex V0,V1,V2;
735        real_t eps=1e-13;
736        if(e0[0]==e0[1] and e1[0]==e1[1] and e2[0]==e2[1]) {
737            V0=(*(*verticesListesIntersection[numobject])[e0[0]]);
738            V1=(*(*verticesListesIntersection[numobject])[e1[0]]);
739            V2=(*(*verticesListesIntersection[numobject])[e2[0]]);
740            real_t normV0V1=Norm(V0.Vector3()-V1.Vector3());
741            real_t normV0V2=Norm(V0.Vector3()-V2.Vector3());
742            real_t normV1V2=Norm(V2.Vector3()-V1.Vector3());
743            if(normV0V1>eps and normV0V2>eps and normV1V2>eps) {
744                return false;
745            }
746            else
747                return true;
748      }      }
749      localTriangleList.clear();      else
750    }          return false;
751  }  }
752    
753    void SurfaceMeshGenerator::Internals::plotliste(size_t nobject) {
754    
755  void SurfaceMeshGenerator::Internals::__InTriangle(const Vertex& P,TinyVector<3>& P0,TinyVector<3>& V12,TinyVector<3>& V02,size_t& compt)  std::ofstream sortie2("liste2.mesh");
756  {  
757        sortie2<<"MeshVersionFormatted 1"<<" \n";
758        sortie2<<" \n";
759        sortie2<<"Dimension"<<" \n";
760        sortie2<<"3"<<" \n";
761        sortie2<<" \n";
762        sortie2<<"Vertices"<<" \n";
763        size_t nn=0;
764        for(size_t no=0;no<nobject;++no)
765            nn+=(*triangleListesIntersectionEdges[no]).size();
766        sortie2<<3*nn<<"\n";
767        for(size_t nz=0;nz<nobject;++nz) {
768            for(std::list<TinyVector<3,TinyVector<2,Edge::Pair> > >::iterator currentTriangle
769         = (*triangleListesIntersectionEdges[nz]).begin();
770         currentTriangle != (*triangleListesIntersectionEdges[nz]).end();
771         ++currentTriangle) {
772                const Edge& e0=(*currentTriangle)[0][0];
773                const Edge& e1=(*currentTriangle)[1][0];
774                const Edge& e2=(*currentTriangle)[2][0];
775                const Vertex& V0=(*(*verticesListesIntersection[nz])[e0]);
776                const Vertex& V1=(*(*verticesListesIntersection[nz])[e1]);
777                const Vertex& V2=(*(*verticesListesIntersection[nz])[e2]);
778    
779                sortie2<<V0.Vector3()[0]<<" "<<V0.Vector3()[1]<<" "<<V0.Vector3()[2]<<" 0 \n";
780                sortie2<<V1.Vector3()[0]<<" "<<V1.Vector3()[1]<<" "<<V1.Vector3()[2]<<" 0 \n";
781                sortie2<<V2.Vector3()[0]<<" "<<V2.Vector3()[1]<<" "<<V2.Vector3()[2]<<" 0 \n";
782            }
783        }
784        sortie2<<" \n";
785        sortie2<<"Triangles"<<"\n";
786        sortie2<<nn<<"\n";
787    
788        for(size_t nj=0;nj<nn;++nj) {
789            sortie2<<3*nj+1<<" "<<3*nj+2<<" "<<3*nj+3;
790            if(nj<(*triangleListesIntersectionEdges[0]).size())
791                sortie2<<" 2 \n";
792            else if(nj<(*triangleListesIntersectionEdges[0]).size()+(*triangleListesIntersectionEdges[1]).size())
793                sortie2<<" 3 \n";
794            else
795                sortie2<<" 4 \n";
796        }
797        sortie2<<"End"<<" \n";
798    }    
799    size_t SurfaceMeshGenerator::Internals::
800    __IsInTriangle(const TinyVector<2,size_t> & num,const TinyVector<3,TinyVector<3> > & triangleVector ) {
801    
     compt=0;  
802      real_t det,eps=-1e-10;      real_t det,eps=-1e-10;
803      det=V02[0]*V12[1]-V02[1]*V12[0];      const TinyVector<3>& V02=triangleVector[0];
804      if(det<eps or det>-eps)      const TinyVector<3>& V12=triangleVector[1];
805      {      const TinyVector<3>& P0=triangleVector[2];
806          //std::cout<<"det "<<det<<"\n";      det=V02[num[0]]*V12[num[1]]-V02[num[1]]*V12[num[0]];
807    
808        if(det<eps or det>-eps) {
809          TinyMatrix<2,2> mat;          TinyMatrix<2,2> mat;
810          mat(0,0)=V02[0];          mat(0,0)=V02[num[0]];
811          mat(1,0)=V02[1];          mat(1,0)=V02[num[1]];
812          mat(0,1)=V12[0];          mat(0,1)=V12[num[0]];
813          mat(1,1)=V12[1];          mat(1,1)=V12[num[1]];
814          TinyVector<2> bary,secondmb;          TinyVector<2> bary,secondmb;
815          secondmb[0]=P0[0];          secondmb[0]=P0[num[0]];
816          secondmb[1]=P0[1];          secondmb[1]=P0[num[1]];
817          bary=mat.invert()*secondmb;          bary=mat.invert()*secondmb;
818          real_t temp=(1.-bary[0]-bary[1]);          real_t temp=(1.-bary[0]-bary[1]);
819          if(bary[0]>=eps && bary[1]>=eps && temp>=eps)          if(bary[0]>=eps && bary[1]>=eps && temp>=eps) {
820          {              return 1;
821              compt+=1;          }
822              //  std::cout<<bary[0]<<" "<<bary[1]<<" "<<1.-bary[0]-bary[1]<<" "<<compt<<"\n";          else {
823       //std::cout<<verif[0]<<" "<<verif[1]<<" \n";              return 0;
      //std::cout<<P0[0]<<" "<<P0[2]<<" \n";  
824          }          }
825      }      }
826      else      else
827          compt+=1;          return 1;
828      {  }
         det=V02[0]*V12[2]-V02[2]*V12[0];  
         if(det<eps or det>-eps)  
         {  
             TinyMatrix<2,2> mat;  
             mat(0,0)=V02[0];  
             mat(1,0)=V02[2];  
             mat(1,1)=V12[2];  
             mat(0,1)=V12[0];  
             TinyVector<2> bary,secondmb;  
             secondmb[0]=P0[0];  
             secondmb[1]=P0[2];  
             bary=mat.invert()*secondmb;  
             TinyVector<2> verif=mat*bary;  
             if(bary[0]>=eps && bary[1]>=eps && (1-bary[0]-bary[1])>=eps)  
             {  
                 compt+=1;  
                 //std::cout<<bary[0]<<" "<<bary[1]<<" "<<1.-bary[0]-bary[1]<<"\n";  
   //std::cout<<verif[0]<<" "<<verif[1]<<" \n";  
   //std::cout<<P0[0]<<" "<<P0[2]<<" \n";  
829    
830    void SurfaceMeshGenerator::Internals::
831    __InTriangle(const Vertex& P,
832                 const TinyVector<3, Vertex* > & V,
833                 size_t& compt) {
834    
835              }      compt=0;
         }  
         else  
             compt+=1;  
     }  
     {  
         det=V02[1]*V12[2]-V02[2]*V12[1];  
         if(det<eps or det>-eps)  
         {  
             TinyMatrix<2,2> mat;  
             mat(0,0)=V02[1];  
             mat(1,0)=V02[2];  
             mat(1,1)=V12[2];  
             mat(0,1)=V12[1];  
             TinyVector<2> bary,secondmb;  
             secondmb[0]=P0[1];  
             secondmb[1]=P0[2];  
             bary=mat.invert()*secondmb;  
             TinyVector<2> verif=mat*bary;  
             if(bary[0]>=eps && bary[1]>=eps && (1-bary[0]-bary[1])>=eps)  
             {  
                 compt+=1;  
                 //std::cout<<bary[0]<<" "<<bary[1]<<" "<<1.-bary[0]-bary[1]<<" "<<compt<<"\n";  
   //std::cout<<verif[0]<<" "<<verif[1]<<" \n";  
   //std::cout<<P0[0]<<" "<<P0[1]<<" \n";  
836    
837              }      TinyVector<3,TinyVector<3> >  triangleVector;
838          }      triangleVector[0]=-(*V[2]).Vector3() + (*V[0]).Vector3();
839          else      triangleVector[1]= -(*V[2]).Vector3() + (*V[1]).Vector3();
840              compt+=1;      triangleVector[2]=P.Vector3()-(*V[2]).Vector3();
841      }      TinyVector<2,size_t>  num;
842                                    
843        num[0]=0;
844        num[1]=1;
845        size_t ajout=__IsInTriangle(num,triangleVector);
846        compt+=ajout;
847    
848        num[0]=0;
849        num[1]=2;
850        ajout=__IsInTriangle(num,triangleVector);
851        compt+=ajout;
852    
853        num[0]=1;
854        num[1]=2;
855        ajout=__IsInTriangle(num,triangleVector);
856        compt+=ajout;
857        
858  }  }
859    
860  void SurfaceMeshGenerator::Internals::__CalculPoint(const Vertex& V0, const Vertex& V1 , const Vertex& V2,  TinyVector<4, real_t>  SurfaceMeshGenerator::Internals::
861                                                      const Vertex& Vo0, const Vertex& Vo1 , const Vertex& Vo2,  __EquationPlan(const TinyVector<3, Vertex* > & V){
862                                                      int nn, int& ttn0, int& ttn1 ,int& ttn2,  
                                                     Vertex *& vvn0 , Vertex *& vvn1 , Vertex *& vvn2)  
 {  
     TinyVector<3>  V01,V12,V02;  
     V01 = V1.Vector3() - V0.Vector3();  
     V12 = V2.Vector3() - V1.Vector3();  
     V02 = V2.Vector3() - V0.Vector3();  
       
863      //equation of plan      //equation of plan
864      TinyVector<3>  Vo01,Vo12,Vo02;  
865      Vo01=Vo1.Vector3()-Vo0.Vector3();      TinyVector<3> V01=(*V[1]).Vector3()-(*V[0]).Vector3();
866      Vo12=Vo2.Vector3()-Vo1.Vector3();      TinyVector<3> V12=(*V[2]).Vector3()-(*V[1]).Vector3();
867      Vo02=Vo2.Vector3()-Vo0.Vector3();      TinyVector<3> V02=(*V[2]).Vector3()-(*V[0]).Vector3();
868    
869      TinyVector<3> coeffo;      TinyVector<3> coeffo;
870      coeffo=Vo01^Vo02;      coeffo=V01^V02;
871      real_t d2 = -(coeffo*Vo0.Vector3());      real_t d2 = -(coeffo*(*V[0]).Vector3());
872    
873        TinyVector<4> coeff;
874        coeff[0]=coeffo[0];
875        coeff[1]=coeffo[1];
876        coeff[2]=coeffo[2];
877        coeff[3]=d2;
878    
879        return coeff;
880        
881    }
882    
883    void SurfaceMeshGenerator::Internals::
884    __DeterminatePoint(const size_t& ncase,
885                       const Edge::Pair & eo0,
886                       const Edge::Pair & eo1,
887                       const Vertex & Vo0,
888                       const Vertex & Vo1,
889                       const TinyVector<3, Vertex* > & V,
890                       IntersectionPoints& Pinter) {
891        std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > p1,p0;
892        TinyVector<2,Edge::Pair> pt1,pt0;
893        pt0[0]=eo0;
894        pt0[1]=eo1;
895        pt1[0]=eo1;
896        pt1[1]=eo0;
897        p1.first=pt1;
898        p0.first=pt0;
899        TinyVector<2,size_t> nuum;
900        nuum[0]=Pinter.numpair[0];
901        nuum[1]=Pinter.numpair[1];
902        p1.second=nuum;
903        p0.second=nuum;
904    
905        real_t epsilon=1e-10;
906      real_t pscalar,k;      real_t pscalar,k;
907      size_t num;      size_t num;
908      {      Vertex P;
909          Vertex P;      TinyVector<4> coeffo=__EquationPlan(V);
910          real_t norm1=Norm(V01);      TinyVector<3> coeff;
911          pscalar=coeffo*V01;      coeff[0]=coeffo[0];
912          num=0;      coeff[1]=coeffo[1];
913          if(std::abs(pscalar)>1e-13) {      coeff[2]=coeffo[2];
914              k=-(coeffo[0]*V0[0]+coeffo[1]*V0[1]+coeffo[2]*V0[2]+d2)/pscalar;      real_t d=coeffo[3];
915              P=V0+k*V01;  
916              TinyVector<3> P0=P.Vector3()-Vo2.Vector3(),P1=Vo12*(-1.),P2=Vo02*(-1.);      TinyVector<3> Vo01 = Vo1.Vector3() - Vo0.Vector3();
917              real_t norm2=Norm(P.Vector3()-V0.Vector3());      TinyVector<3> V01 = (*V[1]).Vector3() - (*V[0]).Vector3();
918              real_t norm3=Norm(P.Vector3()-V1.Vector3());      TinyVector<3> V12 = (*V[2]).Vector3() - (*V[1]).Vector3();
919              if(norm2<=norm1 and norm3<=norm1)   {      TinyVector<3> V02 = (*V[2]).Vector3() - (*V[0]).Vector3();
920                  __InTriangle(P,P0,P1,P2,num);      
921                  if(num==3)      real_t norm1=Norm(Vo01);
922                  {      pscalar=coeff*Vo01;
923                      //std::cout<<"droite1 "<<P.Vector3()[0]<<" "<<P.Vector3()[1]<<" "<<P.Vector3()[2]<<"  \n";      num=0;
924                      ttn0=1;      if(std::abs(pscalar)>1e-10) {
925                      vvn0=new Vertex(P.Vector3());          k=-(coeff*Vo0.Vector3()+d)/pscalar;
926                      //std::cout<<"P "<<nn<<" "<< vv[0][0].Vector3()[0]<<" "<<(vv[0][0]).Vector3()[1]<<" "<<(vv[0][0]).Vector3()[2]<<"\n";          P=Vo0+k*Vo01;
927        //std::cout<<"verif "<<coeffo[0]*P.Vector3()[0]+coeffo[1]*P.Vector3()[1]+coeffo[2]*P.Vector3()[2]+d2<<"\n";          
928        //std::cout<<"verif "<<coeff[0]*P.Vector3()[0]+coeff[1]*P.Vector3()[1]+coeff[2]*P.Vector3()[2]+d<<"\n";          TinyVector<3> V0P=P.Vector3()-(*V[0]).Vector3();
929            TinyVector<3> V1P=P.Vector3()-(*V[1]).Vector3();
930            TinyVector<3> V2P=P.Vector3()-(*V[2]).Vector3();
931    
932            real_t norm2=Norm(P.Vector3()-Vo0.Vector3());
933            real_t norm3=Norm(P.Vector3()-Vo1.Vector3());
934            real_t inV01=Norm(V01^V0P);
935            real_t inV12=Norm(V12^V1P);
936            real_t inV02=Norm(V02^V2P);
937    
938            if(norm2>epsilon and norm3>epsilon and std::abs((norm2+norm3)-norm1) < epsilon ) {
939                if(ncase<3) {
940                    __InTriangle(P,V,num);
941                    if(num==3) {
942                        if (listVerticesIntersectionPerEdges.find(p1) == listVerticesIntersectionPerEdges.end()
943              and listVerticesIntersectionPerEdges.find(p0) == listVerticesIntersectionPerEdges.end()) {
944                            listVerticesIntersectionPerEdges[p0]=new Vertex (P.Vector3());
945                            Pinter.addPoint(p0);
946                        }
947                        else if(listVerticesIntersectionPerEdges.find(p1) != listVerticesIntersectionPerEdges.end()){
948                            Vertex Ptemp=(*listVerticesIntersectionPerEdges[p1]);
949                            Pinter.addPoint(p1);
950                        }
951                        else if(listVerticesIntersectionPerEdges.find(p0) != listVerticesIntersectionPerEdges.end()){
952                            Vertex Ptemp=(*listVerticesIntersectionPerEdges[p0]);
953                            Pinter.addPoint(p0);
954                        }
955                        else {
956                            ffout(0)<<"pbbbbb\n";
957                        }
958                        //if(norm2>epsilon and norm3>epsilon)
959                        Pinter.cutTriangle1[ncase]=Pinter.size();
960                  }                  }
961              }              }
962          }              else if(inV01> epsilon and inV02> epsilon and inV12> epsilon ) {
963      }                  __InTriangle(P,V,num);
     {  
         Vertex P;  
         pscalar=coeffo*V02;  
         real_t norm1=Norm(V02);  
         if(std::abs(pscalar)>1e-13) {  
             k=-(coeffo[0]*V0[0]+coeffo[1]*V0[1]+coeffo[2]*V0[2]+d2)/pscalar;  
             P=V0+k*V02;  
             TinyVector<3> P0=P.Vector3()-Vo2.Vector3(),P1=Vo12*(-1.),P2=Vo02*(-1.);  
             real_t norm2=Norm(P.Vector3()-V0.Vector3());  
             real_t norm3=Norm(P.Vector3()-V2.Vector3());  
             //std::cout<<"norm "<<norm3<<" "<<norm2<<" "<<norm1<<"\n";  
             if(norm2<=norm1 and norm3<=norm1) {  
                 __InTriangle(P,P0,P1,P2,num);  
964                  if(num==3) {                  if(num==3) {
965                     // std::cout<<"droite2 "<<P.Vector3()[0]<<" "<<P.Vector3()[1]<<" "<<P.Vector3()[2]<<"  \n";                      if (listVerticesIntersectionPerEdges.find(p1) == listVerticesIntersectionPerEdges.end()
966                      ttn1=1;            and listVerticesIntersectionPerEdges.find(p0) == listVerticesIntersectionPerEdges.end()) {
967                      vvn1=new Vertex(P.Vector3());                          listVerticesIntersectionPerEdges[p0]=new Vertex (P.Vector3());
968                            Pinter.addPoint(p0);
969                        }
970                        else if(listVerticesIntersectionPerEdges.find(p1) != listVerticesIntersectionPerEdges.end()){
971                            Vertex Ptemp=(*listVerticesIntersectionPerEdges[p1]);
972                            if(Ptemp.Vector3()==P.Vector3())
973                                Pinter.addPoint(p1);
974                            else {
975                                TinyVector<2,size_t> temp;
976                                temp[0]=Pinter.numpair[1];
977                                temp[1]=Pinter.numpair[0];
978                                p1.second=temp;
979                                listVerticesIntersectionPerEdges[p1]=new Vertex (P.Vector3());
980                                Pinter.addPoint(p1);
981                            }
982                        }
983                        else if(listVerticesIntersectionPerEdges.find(p0) != listVerticesIntersectionPerEdges.end()){
984                            Vertex Ptemp=(*listVerticesIntersectionPerEdges[p0]);
985                            if(Ptemp.Vector3()==P.Vector3())
986                                Pinter.addPoint(p0);
987                            else {
988                                TinyVector<2,size_t> temp;
989                                temp[0]=Pinter.numpair[1];
990                                temp[1]=Pinter.numpair[0];
991                                p0.second=temp;
992                                listVerticesIntersectionPerEdges[p0]=new Vertex (P.Vector3());
993                                Pinter.addPoint(p0);
994                            }
995    
996                        }
997                        else {
998                            ffout(0)<<"pbbbbb\n";
999                        }
1000                  }                  }
1001              }              }
1002          }          }
1003      }      }
1004      {  }
         Vertex P;  
         pscalar=coeffo*V12;  
         real_t norm1=Norm(V12);  
         if(std::abs(pscalar)>1e-13) {  
             k=-(coeffo*V1.Vector3() + d2)/pscalar;  
             P=V1+k*V12;  
             TinyVector<3> P0=P.Vector3()-Vo2.Vector3(),P1=Vo12*(-1.),P2=Vo02*(-1.);  
             real_t norm2=Norm(P.Vector3()-V1.Vector3());  
             real_t norm3=Norm(P.Vector3()-V2.Vector3());  
1005    
1006              if(norm2<=norm1 and norm3<=norm1) {  void SurfaceMeshGenerator::Internals::
1007                  __InTriangle(P,P0,P1,P2,num);  __CalculPoint(const size_t& objectNumber,const size_t& numobject,
1008                  if(num==3) {                SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList2::const_iterator & currentTriangle,
1009  //std::cout<<"droite3 "<<P.Vector3()[0]<<" "<<P.Vector3()[1]<<" "<<P.Vector3()[2]<<"  \n";                SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList2::const_iterator & currentTriangle1,
1010                  IntersectionPoints& Pinter) {
1011    
1012        //ffout(0) <<"here ok\n";
1013        size_t num1=objectNumber,num0=numobject;
1014        
1015        TinyVector<3, Vertex* > V012;
1016        const Edge::Pair& e0=(*(*currentTriangle1))[0][0];
1017        const Edge::Pair& e1=(*(*currentTriangle1))[1][0];
1018        const Edge::Pair& e2=(*(*currentTriangle1))[2][0];
1019        V012[0] = ((*verticesListesIntersection[num0])[e0]);
1020        V012[1] = ((*verticesListesIntersection[num0])[e1]);
1021        V012[2] = ((*verticesListesIntersection[num0])[e2]);
1022    
1023        TinyVector<3, Vertex* > Vo012;
1024        const Edge::Pair& eo0=(*(*currentTriangle))[0][0];
1025        const Edge::Pair& eo1=(*(*currentTriangle))[1][0];
1026        const Edge::Pair& eo2=(*(*currentTriangle))[2][0];
1027        Vo012[0] = ((*verticesListesIntersection[num1])[eo0]);
1028        Vo012[1] = ((*verticesListesIntersection[num1])[eo1]);
1029        Vo012[2] = ((*verticesListesIntersection[num1])[eo2]);
1030    
1031        //for each edge of the triangle V0V1V2
1032        Pinter.numpair[0]=num0;
1033        Pinter.numpair[1]=num1;
1034        __DeterminatePoint(0,e0,e1,(*V012[0]),(*V012[1]),Vo012,Pinter);
1035        __DeterminatePoint(1,e0,e2,(*V012[0]),(*V012[2]),Vo012,Pinter);
1036        __DeterminatePoint(2,e1,e2,(*V012[1]),(*V012[2]),Vo012,Pinter);
1037    
1038    
1039        __DeterminatePoint(3,eo0,eo1,(*Vo012[0]),(*Vo012[1]),V012,Pinter);
1040        __DeterminatePoint(3,eo0,eo2,(*Vo012[0]),(*Vo012[2]),V012,Pinter);
1041        __DeterminatePoint(3,eo1,eo2,(*Vo012[1]),(*Vo012[2]),V012,Pinter);
1042        
1043    }
1044    
1045    void SurfaceMeshGenerator::Internals::
1046    __inout(const size_t& numobject,const int & ttn0,const TinyVector<2,Edge::Pair>& e0,
1047            const TinyVector<2,Edge::Pair>& e1,size_t& ii0,size_t& ii1,TriangleCut& K) {
1048    
1049        size_t& re0=(*edgesRefCreate[numobject])[e0[0]];
1050        size_t& re1=(*edgesRefCreate[numobject])[e1[0]];
1051        
1052      ii0=0;
1053      ii1=0;
1054      //if(e0[0]!=e0[1])
1055        //  re0=2;
1056      //if(e1[0]!=e1[1])
1057        //  re1=2;
1058      
1059      if(re0==1 and re1==1) {
1060          ffout(0)<<"--------------- re ==1 1\n";
1061          if(ttn0!=0)
1062              ffout(0)<<"-- cas indetermine\n";
1063      }
1064      else if((re0+re1==1)) {
1065          ffout(0)<<"-------------- re ==1 2\n";
1066        if(ttn0!=0) {
1067          //ffout(0)<<"-- 1 in 0 out on sait traite\n";
1068          ii0=re0;
1069          ii1=re1;
1070        }
1071      }
1072      else if((re0+re1==3)) {
1073          ffout(0)<<"------------- re ==1 3\n";
1074    
1075        if(ttn0==0) {
1076          //ffout(0)<<"-- 1 2 pas coupe\n";
1077          re0=2;
1078          re1=2;
1079        }
1080        else {
1081          //ffout(0)<<"-- 1 2 on sait traite\n";
1082          if(re0==1) {
1083            ii0=0;
1084            ii1=1;
1085            re0=0;
1086          }
1087          else {
1088            ii0=1;
1089            ii1=0;
1090            re1=0;
1091          }
1092        }
1093      }
1094      else if(re0+re1==4 and ttn0==1) {
1095          ffout(0)<<"warning edge ref 2 cut \n";
1096          printpoint(numobject,e0,K);
1097          printpoint(numobject,e1,K);
1098          re0=0;
1099          re1=0;
1100      }
1101      //     else
1102      //ffout(0)<<"--0 out 0 out rien a faire\n";
1103    
1104    }
1105    
1106    void SurfaceMeshGenerator::Internals::
1107    __createPointIntersection(const size_t& numobject,const size_t& objectNumber) {
1108        ObjectPair objectp(numobject,objectNumber);
1109    
1110        //ffout(0)<< "createpointintersection "<<numobject<<" "<<objectNumber<<"\n";
1111        for(SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList2::const_iterator currentTriangle1
1112            =  (*localTriangleList2Edges[numobject]).begin();
1113            currentTriangle1 !=  (*localTriangleList2Edges[numobject]).end();
1114            ++currentTriangle1) {
1115    
1116            PairTriangleEdges pTriangle;
1117            pTriangle.first=(*currentTriangle1);
1118    
1119            for(SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList2::const_iterator currentTriangle
1120         =  (*localTriangleList2Edges[objectNumber]).begin();
1121         currentTriangle !=  (*localTriangleList2Edges[objectNumber]).end();
1122         ++currentTriangle) {
1123                //ffout(0)<<"*************************\n";
1124    
1125                      ttn2=1;              IntersectionPoints Pinter;
1126                      vvn2=new Vertex(P.Vector3());              pTriangle.second=(*currentTriangle);
1127                
1128                __CalculPoint(objectNumber,numobject,currentTriangle,currentTriangle1,Pinter);
1129                
1130                TinyVector<3,size_t> pintercut=Pinter.cutTriangle1;
1131    
1132                if(pintercut[0]+pintercut[1]+pintercut[2]!=0) {            
1133                    if((listVertexIntersectionPerPairEdges[objectp]).find(pTriangle)==
1134         (listVertexIntersectionPerPairEdges[objectp]).end()) {
1135                        //ffout(0) << "rajout de pTriangle\n";
1136                        (listVertexIntersectionPerPairEdges[objectp])[pTriangle] = Pinter;
1137                    }
1138                }
1139                else if(Pinter.size()!=0 and numobject<objectNumber) {
1140                    (listVertexIntersectionPerPairEdges[objectp])[pTriangle] = Pinter;
1141                }
1142    
1143                if((listVertexIntersectionPerPairEdges[objectp]).find(pTriangle)!=
1144            (listVertexIntersectionPerPairEdges[objectp]).end() and
1145            numobject<objectNumber) {
1146                    TinyVector<3,size_t>& cutTr2=
1147                        (listVertexIntersectionPerPairEdges[objectp])[pTriangle].cutTriangle2;
1148    
1149                    for(size_t i=0 ; i<Pinter.size() ; ++i) {
1150                        real_t epsilon=1e-10;
1151                        const Edge::Pair& eo0=(*(*currentTriangle))[0][0];
1152                        const Edge::Pair& eo1=(*(*currentTriangle))[1][0];
1153                        const Edge::Pair& eo2=(*(*currentTriangle))[2][0];
1154                        const Vertex& Vo0 = *((*verticesListesIntersection[objectNumber])[eo0]);
1155                        const Vertex& Vo1 = *((*verticesListesIntersection[objectNumber])[eo1]);
1156                        const Vertex& Vo2 = *((*verticesListesIntersection[objectNumber])[eo2]);
1157                        std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> >
1158                            num=((listVertexIntersectionPerPairEdges[objectp])[pTriangle].listofPoints[i]);
1159                        Vertex P;
1160                        if(listVerticesIntersectionPerEdges.find(num)!=listVerticesIntersectionPerEdges.end()) {
1161                            P=(*listVerticesIntersectionPerEdges[num]);
1162                        }
1163    
1164                        real_t Vo01=Norm(Vo1.Vector3()-Vo0.Vector3());
1165                        real_t Vo12=Norm(Vo2.Vector3()-Vo1.Vector3());
1166                        real_t Vo02=Norm(Vo2.Vector3()-Vo0.Vector3());
1167                        real_t norm0=Norm(P.Vector3()-Vo0.Vector3());
1168                        real_t norm1=Norm(P.Vector3()-Vo1.Vector3());
1169                        real_t norm2=Norm(P.Vector3()-Vo2.Vector3());
1170                        //for filter points
1171                        real_t eps=1e-10;
1172                        if(norm0>eps and norm1>eps and std::abs((norm0+norm1)-Vo01)<epsilon) {
1173                            cutTr2[0]=i+1;
1174                        }
1175                        if(norm0>eps and norm2>eps and std::abs((norm0+norm2)-Vo02)<epsilon) {
1176                            cutTr2[1]=i+1;
1177                        }
1178                        if(norm2>eps and norm1>eps and std::abs((norm2+norm1)-Vo12)<epsilon) {
1179                            cutTr2[2]=i+1;
1180                        }
1181                  }                  }
1182              }              }
1183    
1184          }          }
1185      }      }
1186  }  }
1187    
1188    
1189    std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > SurfaceMeshGenerator::Internals::
1190    __determinateNumber(const int& p,const TinyVector<3,std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > >& pointsint) {
1191    
1192        if(p==0) {
1193            return pointsint[0];
1194        }
1195        else if(p==1) {
1196            return pointsint[1];
1197        }
1198        else {
1199            return pointsint[2];
1200        }
1201    }
1202    
1203  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1204  __inout(size_t iter,const size_t nn,int& ttn0,  __addTriangle(const size_t& numobject,
1205                 int& ttn1,                   const TinyVector<2,Edge::Pair>& numA,
1206                 int& ttn2,                   const TinyVector<2,Edge::Pair>& numB, const TinyVector<2,Edge::Pair>& numC,
1207          size_t& re0,size_t& re1,size_t& re2,size_t& ii0,size_t& ii1,size_t& ii2)                   size_t& kill) {
1208  {  
1209        TinyVector<3,TinyVector<2,Edge::Pair> > trianglebis;
1210        trianglebis[0]=numA;
1211        trianglebis[1]=numB;
1212        trianglebis[2]=numC;
1213        (*listTriangleintersection[numobject]).push_back(trianglebis);
1214        localListTriangleintersection.push_back(trianglebis);
1215        kill=1;
1216    }
1217    
1218      if(re0==1 and re1==1)  void SurfaceMeshGenerator::Internals::
1219      {  __determinateCase(const int& ncase,
1220          if(ttn0!=0 and iter>1)                    const TriangleCut & K,
1221  std::cout<<"-- cas indetermine\n";                    Vertex & A,Vertex & B,Vertex & C,
1222                      TinyVector<2,Edge::Pair> & e0,TinyVector<2,Edge::Pair> & e1,
1223                      TinyVector<2,Edge::Pair> & e2,
1224                      TinyVector<2,int>& tt,
1225                      TinyVector<2,int>& pp) {
1226    
1227        const TinyVector<3,int>& tt2=K.edgecut;
1228        const size_t& numobject=K.numobject;
1229        Vertex V0,V1,V2;
1230        if(K.listpairedges[0][0]==K.listpairedges[0][1])
1231            V0=(*(*verticesListesIntersection[numobject])[K.listpairedges[0][0]]);
1232        else {
1233            std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > p0;
1234            p0.first=K.listpairedges[0];
1235            TinyVector<2,size_t> num;
1236            num[0]=K.numobject;
1237            num[1]=K.otherobject;
1238            p0.second=num;
1239            SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator
1240            ip0=listVerticesIntersectionPerEdges.find(p0);
1241            if(ip0!=listVerticesIntersectionPerEdges.end())
1242                V0 =*(*ip0).second;
1243        }
1244        if(K.listpairedges[1][0]==K.listpairedges[1][1])
1245            V1=(*(*verticesListesIntersection[numobject])[K.listpairedges[1][0]]);
1246        else {
1247            std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > p0;
1248            p0.first=K.listpairedges[1];
1249            TinyVector<2,size_t> num;
1250            num[0]=0;
1251            num[1]=0;
1252            p0.second=num;
1253            SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator
1254            ip1=listVerticesIntersectionPerEdges.find(p0);
1255            if(ip1!=listVerticesIntersectionPerEdges.end())
1256                V1 =*(*ip1).second;
1257      }      }
1258      else if((re0+re1==1))      
1259      {      if(K.listpairedges[2][0]==K.listpairedges[2][1])
1260          if(ttn0==0)          V2=(*(*verticesListesIntersection[numobject])[K.listpairedges[2][0]]);
1261          {      else {
1262              //std::cout<<"-- 1 0\n";          std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > p0;
1263              re0=0;          p0.first=K.listpairedges[2];
1264              re1=0;          TinyVector<2,size_t> num;
1265          }          num[0]=0;
1266          else          num[1]=0;
1267          {          p0.second=num;
1268              //std::cout<<"-- 1 in 0 out on sait traite\n";          SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator
1269              ii0=re0;          ip2=listVerticesIntersectionPerEdges.find(p0);
1270              ii1=re1;          if(ip2!=listVerticesIntersectionPerEdges.end())
1271                V2 =*(*ip2).second;
1272        }
1273        
1274    
1275        if(ncase==0) {
1276            tt[0]=tt2[0];
1277            tt[1]=tt2[1];
1278            pp[0]=0;
1279            pp[1]=1;
1280            A=V0;
1281            e0=K.listpairedges[0];
1282            if(tt[0]==0 and tt[1]==1) {
1283                B=V2;
1284                e1=K.listpairedges[2];
1285            }
1286            if(tt[0]==1 and tt[1]==0) {
1287                C=V1;
1288                e2=K.listpairedges[1];
1289            }
1290        }
1291        else if(ncase==1) {
1292            tt[0]=tt2[0];
1293            tt[1]=tt2[2];
1294            pp[0]=0;
1295            pp[1]=2;
1296            A=V1;
1297            e0=K.listpairedges[1];
1298            if(tt[0]==0 and tt[1]==1) {
1299                B=V2;
1300                e1=K.listpairedges[2];
1301            }
1302            if(tt[0]==1 and tt[1]==0) {
1303                C=V0;
1304                e2=K.listpairedges[0];
1305            }
1306        }
1307        else if(ncase==2) {
1308            tt[0]=tt2[1];
1309            tt[1]=tt2[2];
1310            pp[0]=1;
1311            pp[1]=2;
1312            A=V2;
1313            e0=K.listpairedges[2];
1314            if(tt[0]==0 and tt[1]==1) {
1315                B=V0;
1316                e1=K.listpairedges[0];
1317            }
1318            if(tt[0]==1 and tt[1]==0) {
1319                C=V1;
1320                e2=K.listpairedges[1];
1321          }          }
1322      }      }
1323      else if((re0==2 and re1==0) or  
1324              (re0==0 and re1==2))  
1325      {  }
1326          if(ttn0==0)  void SurfaceMeshGenerator::Internals::
1327          {  __createCase1(int ncase,
1328              if(ttn0+ttn1+ttn2==0)                const TriangleCut& K,
1329              {                size_t & kill) {
1330                  re0=2;      const size_t& numobject=K.numobject;
1331                  re1=2;      const TinyVector<2,Edge::Pair>& e0=K.listpairedges[0];
1332        const TinyVector<2,Edge::Pair>& e1=K.listpairedges[1];
1333        const TinyVector<2,Edge::Pair>& e2=K.listpairedges[2];
1334        Vertex V0,V1,V2;
1335        std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > p0,p1,p2;
1336    
1337        if(e0[0]==e0[1])
1338            V0 = *((*verticesListesIntersection[numobject])[e0[0]]);
1339        else {
1340            p0.first=e0;
1341            TinyVector<2,size_t> num;
1342            num[0]=K.numobject;
1343            num[1]=K.otherobject;
1344            p0.second=num;
1345            SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator
1346            ip0=listVerticesIntersectionPerEdges.find(p0);
1347            if(ip0!=listVerticesIntersectionPerEdges.end())
1348                V0 =*(*ip0).second;
1349        }
1350        if(e1[0]==e1[1])
1351            V1 = *((*verticesListesIntersection[numobject])[e1[0]]);
1352        else {
1353            p1.first=e1;
1354            TinyVector<2,size_t> num;
1355            num[0]=K.numobject;
1356            num[1]=K.otherobject;
1357            p0.second=num;
1358            SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator
1359            ip1=listVerticesIntersectionPerEdges.find(p1);
1360            if(ip1!=listVerticesIntersectionPerEdges.end())
1361                V1 =*(*ip1).second;
1362        }
1363        if(e2[0]==e2[1])
1364            V2 = *((*verticesListesIntersection[numobject])[e2[0]]);
1365        else {
1366            p2.first=e2;
1367            TinyVector<2,size_t> num;
1368            num[0]=K.numobject;
1369            num[1]=K.otherobject;
1370            p0.second=num;
1371            SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator
1372            ip2=listVerticesIntersectionPerEdges.find(p2);
1373            if(ip2!=listVerticesIntersectionPerEdges.end())
1374                V2 =*(*ip2).second;
1375        }
1376        const TinyVector<3,std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > >& pointsint=K.pointsIntersection;
1377        const std::vector< std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> >  >& pointsin=K.pointsin;
1378        const int& numberin=pointsin.size();
1379        Vertex A=V0,B=V1,C=V2;
1380        TinyVector<2,Edge::Pair> eA=e0,eB=e1,eC=e2;
1381        std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > numA,numB,numC,numD;
1382    
1383        TinyVector<2,int> tt;
1384        TinyVector<2,int> pp;
1385    
1386        __determinateCase(ncase,K,A,B,C,eA,eB,eC,tt,pp);
1387    
1388        std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > pA,pB;
1389        TinyVector<2,size_t> nuum;
1390        nuum[0]=0;
1391        nuum[1]=0;
1392        pA.second=nuum;
1393        pB.second=nuum;
1394        pA.first=eA;
1395        pB.first=eB;
1396        numA=pA;
1397        if(numberin==0) {
1398            if(tt[0]==1) {
1399                numB=__determinateNumber(pp[0],pointsint);
1400    
1401                if(listVerticesIntersectionPerEdges.find(numB)!=listVerticesIntersectionPerEdges.end()) {
1402                    B=(*listVerticesIntersectionPerEdges[numB]);
1403                }
1404                else {
1405                    ffout(0)<<"1 "<<pp[0]<<" pbs point inexistant\n";
1406                    std::exit(1);
1407                }
1408                if(tt[1]==1) {
1409    
1410                    numC=__determinateNumber(pp[1],pointsint);
1411                    if(listVerticesIntersectionPerEdges.find(numC)!=listVerticesIntersectionPerEdges.end()) {
1412                        C=(*listVerticesIntersectionPerEdges[numC]);
1413                    }
1414                    else {
1415                        ffout(0)<<"2 "<<pp[1]<<"  pbs point inexistant\n";
1416                        std::exit(1);
1417                    }
1418                    __addTriangle(numobject,numA.first,numB.first,numC.first,kill);
1419    
1420              }              }
1421              else          else {
1422              {              ffout(0) <<"pbs 10\n";
1423  std::cout<<"-- 0 2 pas coupe\n";              printpoint(numobject,K.listpairedges[0],K);
1424  std::cout<<ttn0<<" "<<ttn1<<" "<<ttn2<<"\n";              printpoint(numobject,K.listpairedges[1],K);
1425  std::cout<<re0<<" "<<re1<<"\n";              printpoint(numobject,K.listpairedges[2],K);
1426                  re0=2;              /*listpoints.push_back(new Vertex(C.Vector3()));
1427                  re1=2;              numC=(listpoints.size()-1);
1428                __addTriangle(numobject,numA,numB,numC,kill);
1429                */
1430            }
1431        }
1432        else  {
1433            if(tt[1]==1) {
1434                ffout(0)<<"pbs 11\n";
1435                printpoint(numobject,K.listpairedges[0],K);
1436                printpoint(numobject,K.listpairedges[1],K);
1437                printpoint(numobject,K.listpairedges[2],K);
1438                /* numC=__determinateNumber(pp[1],pointsint);
1439              C=(*listpoints[numC]);
1440              listpoints.push_back(new Vertex(A.Vector3()));
1441              numA=(listpoints.size()-1);
1442              listpoints.push_back(new Vertex(B.Vector3()));
1443              numC=(listpoints.size()-1);
1444              __addTriangle(numobject,numA,numB,numC,kill);*/
1445            }
1446            else {
1447                ffout(0) <<"pbs 111\n";
1448                printpoint(numobject,K.listpairedges[0],K);
1449                printpoint(numobject,K.listpairedges[1],K);
1450                printpoint(numobject,K.listpairedges[2],K);
1451                
1452            }
1453        }
1454      }
1455      else if(numberin==1 or numberin==2) {
1456          
1457          numD=(pointsin[0]);
1458          if(listVerticesIntersectionPerEdges.find(numD)!=listVerticesIntersectionPerEdges.end()) {
1459          Vertex D=(*listVerticesIntersectionPerEdges[numD]);
1460          }
1461          else {
1462              ffout(0)<<"3 "<<pp[0]<<"  pbs point inexistant\n";
1463              std::exit(1);
1464          }
1465          
1466          if(tt[0]==1) {
1467              numB=__determinateNumber(pp[0],pointsint);
1468              if(listVerticesIntersectionPerEdges.find(numB)!=listVerticesIntersectionPerEdges.end()) {
1469                  B=(*listVerticesIntersectionPerEdges[numB]);
1470              }
1471              else {
1472                  ffout(0)<<"4 "<<pp[0]<<"  pbs point inexistant\n";
1473                  std::exit(1);
1474              }
1475    
1476              if(tt[1]==1) {
1477                  numC=__determinateNumber(pp[1],pointsint);
1478                  if(listVerticesIntersectionPerEdges.find(numC)!=listVerticesIntersectionPerEdges.end()) {
1479                  C=(*listVerticesIntersectionPerEdges[numC]);
1480                  }
1481                  else {
1482                      ffout(0)<<"5 "<<pp[1]<<"  pbs point inexistant\n";
1483                      std::exit(1);
1484                  }
1485                  __addTriangle(numobject,numA.first,numB.first,numD.first,kill);
1486                  __addTriangle(numobject,numA.first,numD.first,numC.first,kill);
1487              }
1488              else {
1489                  ffout(0)<<"-------- warning trou in 0\n";
1490                  printpoint(numobject,K.listpairedges[0],K);
1491                  printpoint(numobject,K.listpairedges[1],K);
1492                  printpoint(numobject,K.listpairedges[2],K);      
1493                  __addTriangle(numobject,numA.first,numB.first,numD.first,kill);
1494              }
1495          }
1496          else if(tt[1]==1) {
1497              ffout(0)<<"-------- warning trou in\n";
1498              printpoint(numobject,K.listpairedges[0],K);
1499              printpoint(numobject,K.listpairedges[1],K);
1500              printpoint(numobject,K.listpairedges[2],K);  
1501              numB=__determinateNumber(pp[1],pointsint);
1502              if(listVerticesIntersectionPerEdges.find(numB)!=listVerticesIntersectionPerEdges.end()) {
1503                  B=(*listVerticesIntersectionPerEdges[numB]);
1504              }
1505              else {
1506                  ffout(0)<<"6 "<<pp[0]<<"  pbs point inexistant\n";
1507                  std::exit(1);
1508              }
1509              __addTriangle(numobject,numA.first,numB.first,numD.first,kill);
1510    
1511          }
1512          else {
1513              ffout(0) <<"pbs 1\n";
1514              printpoint(numobject,K.listpairedges[0],K);
1515              printpoint(numobject,K.listpairedges[1],K);
1516              printpoint(numobject,K.listpairedges[2],K);    
1517          }
1518      }
1519      else {
1520        //a faire si plus d' 1 point d'intersection pas dans le triangle
1521        ffout(0) << numberin<<" be careful 1\n";
1522      }
1523    }
1524    
1525    void SurfaceMeshGenerator::Internals::
1526    __createCase2(int ncase,
1527                  const TriangleCut& K,
1528                  size_t & kill) {
1529    
1530      int ttn0,ttn1;
1531      int p0,p1;
1532      const size_t& numobject=K.numobject;
1533      const TinyVector<2,Edge::Pair>& e0=K.listpairedges[0];
1534      const TinyVector<2,Edge::Pair>& e1=K.listpairedges[1];
1535      const TinyVector<2,Edge::Pair>& e2=K.listpairedges[2];
1536      std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > pp0,pp1,pp2;
1537    
1538      Vertex V0,V1,V2;
1539      if(e0[0]==e0[1]) {
1540          V0 = *((*verticesListesIntersection[numobject])[e0[0]]);
1541      }
1542      else {
1543          pp0.first=e0;
1544          TinyVector<2,size_t> num;
1545          num[0]=K.numobject;
1546          num[1]=K.otherobject;
1547          pp0.second=num;
1548            SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator
1549            ip0=listVerticesIntersectionPerEdges.find(pp0);
1550            if(ip0!=listVerticesIntersectionPerEdges.end())
1551                V0 =*(*ip0).second;
1552        }
1553      if(e1[0]==e1[1]) {
1554          V1 = *((*verticesListesIntersection[numobject])[e1[0]]);
1555      }
1556      else {
1557          pp1.first=e1;
1558          TinyVector<2,size_t> num;
1559          num[0]=K.numobject;
1560          num[1]=K.otherobject;
1561          pp1.second=num;
1562          SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator
1563          ip1=listVerticesIntersectionPerEdges.find(pp1);
1564          if(ip1!=listVerticesIntersectionPerEdges.end())
1565              V1 =*(*ip1).second;
1566      }
1567      if(e2[0]==e2[1]) {
1568          V2 = *((*verticesListesIntersection[numobject])[e2[0]]);
1569      }
1570      else {
1571          pp2.first=e2;
1572          TinyVector<2,size_t> num;
1573          num[0]=K.numobject;
1574          num[1]=K.otherobject;
1575          pp2.second=num;
1576          SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator
1577          ip2=listVerticesIntersectionPerEdges.find(pp2);
1578          if(ip2!=listVerticesIntersectionPerEdges.end())
1579              V2 =*(*ip2).second;
1580      }
1581      const TinyVector<3,std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > >& pointsint=K.pointsIntersection;
1582      const std::vector< std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> >  >& pointsin=K.pointsin;
1583      const TinyVector<3,int>&  tt2=K.edgecut;
1584      const int& numberin=pointsin.size();
1585    
1586      Vertex A=V0,B=V1;
1587      TinyVector<2,Edge::Pair> eA=e0,eB=e1;
1588      std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > pA=pp0,pB=pp1;
1589      std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > numA,numB,numC,numD,numE;
1590    
1591      if(ncase==0) {
1592          ttn0=tt2[1];
1593          ttn1=tt2[2];
1594          p0=1;
1595          p1=2;
1596      }
1597      else if(ncase==1) {
1598        ttn0=tt2[0];
1599        ttn1=tt2[2];
1600        p0=0;
1601        p1=2;
1602        A=V0;
1603        eA=K.listpairedges[0];
1604        B=V2;
1605        eB=K.listpairedges[2];
1606      }
1607      else if(ncase==2) {
1608        ttn0=tt2[0];
1609        ttn1=tt2[1];
1610        p0=0;
1611        p1=1;
1612        A=V1;
1613        eA=K.listpairedges[1];
1614        B=V2;
1615        eB=K.listpairedges[2];
1616      }
1617    
1618      pA.first=eA;
1619      pB.first=eB;
1620      numA=pA;
1621      numB=pB;
1622      
1623      if(numberin==0) {
1624          if(ttn0==1) {
1625              Vertex C;
1626              numC=__determinateNumber(p0,pointsint);
1627              if(listVerticesIntersectionPerEdges.find(numC)!=listVerticesIntersectionPerEdges.end()) {
1628                  C=(*listVerticesIntersectionPerEdges[numC]);
1629              }
1630              if(ttn1==1) {
1631                  Vertex D;
1632                  numD=__determinateNumber(p1,pointsint);
1633                  if(listVerticesIntersectionPerEdges.find(numD)!=listVerticesIntersectionPerEdges.end()) {
1634                      D=(*listVerticesIntersectionPerEdges[numD]);
1635                  }
1636                  __addTriangle(numobject,numA.first,numB.first,numC.first,kill);
1637                  __addTriangle(numobject,numC.first,numB.first,numD.first,kill);
1638              }
1639              else {
1640                 // ffout(0) <<"pbs 20\n";
1641                  __addTriangle(numobject,numA.first,numB.first,numC.first,kill);
1642              }
1643          }
1644          else {
1645              if(ttn1==1) {
1646                 // ffout(0) <<"pbs 210\n";
1647                  Vertex C;
1648                  numC=__determinateNumber(p1,pointsint);
1649                  if(listVerticesIntersectionPerEdges.find(numC)!=listVerticesIntersectionPerEdges.end()) {
1650                      C=(*listVerticesIntersectionPerEdges[numC]);
1651                  }
1652                  __addTriangle(numobject,numA.first,numB.first,numC.first,kill);
1653    
1654              }
1655              else {
1656                  ffout(0) <<"pbs 21\n";
1657                  printpoint(numobject,K.listpairedges[0],K);
1658                  printpoint(numobject,K.listpairedges[1],K);
1659                  printpoint(numobject,K.listpairedges[2],K);
1660                  
1661              }
1662          }
1663      }
1664      else if(numberin==1 or numberin==2) {
1665          numE=(pointsin[0]);
1666    
1667          if(ttn0==1) {
1668              Vertex C;
1669              numC=__determinateNumber(p0,pointsint);
1670              if(listVerticesIntersectionPerEdges.find(numC)!=listVerticesIntersectionPerEdges.end()) {
1671                  C=(*listVerticesIntersectionPerEdges[numC]);
1672              }
1673              if(ttn1==1) {
1674                  Vertex D;
1675                  numD=__determinateNumber(p1,pointsint);
1676                  if(listVerticesIntersectionPerEdges.find(numD)!=listVerticesIntersectionPerEdges.end()) {
1677                      D=(*listVerticesIntersectionPerEdges[numD]);
1678                  }
1679                  __addTriangle(numobject,numA.first,numC.first,numE.first,kill);
1680                  __addTriangle(numobject,numA.first,numE.first,numB.first,kill);
1681                  __addTriangle(numobject,numB.first,numE.first,numD.first,kill);
1682              }
1683              else {
1684                  __addTriangle(numobject,numA.first,numB.first,numE.first,kill);
1685                  __addTriangle(numobject,numC.first,numE.first,numB.first,kill);
1686              }
1687          }
1688          else {
1689              if(ttn1==1) {
1690                  Vertex D;
1691                  numD=__determinateNumber(p1,pointsint);
1692                  if(listVerticesIntersectionPerEdges.find(numD)!=listVerticesIntersectionPerEdges.end()) {
1693                      D=(*listVerticesIntersectionPerEdges[numD]);
1694                  }
1695                  __addTriangle(numobject,numD.first,numE.first,numB.first,kill);
1696                  __addTriangle(numobject,numA.first,numB.first,numE.first,kill);
1697              }
1698              else {
1699                  __addTriangle(numobject,numA.first,numB.first,numE.first,kill);
1700              }
1701          }
1702      }
1703      else if(numberin>2) {
1704          ffout(0) << numberin <<"be careful 2\n";
1705      }
1706    
1707    }
1708    
1709    void SurfaceMeshGenerator::Internals::
1710    __createLocalListIntersection(const size_t& objectNumber,
1711                                  const TinyVector<3,TinyVector<2,Edge::Pair> >& currentTriangle,
1712                                  TriangleCut & K) {
1713    
1714        TinyVector<3,int>& tt=K.edgecut;
1715        TinyVector<3,std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > >& pointsIntersection=K.pointsIntersection;
1716        std::vector< std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > >& pointsin=K.pointsin;
1717    
1718        int numberin=pointsin.size();
1719        const size_t& numobject=K.numobject;
1720    
1721        ObjectPair objectp(numobject,objectNumber);
1722        if(numobject>objectNumber) {
1723            objectp.first=objectNumber;
1724            objectp.second=numobject;
1725        }
1726    
1727        for(SurfaceMeshGenerator::Internals::TriangleVerticesMap::iterator currentTriangle1
1728            = (*triangleListesIntersectionEdgesTemp[objectNumber]).begin();
1729            currentTriangle1 != (*triangleListesIntersectionEdgesTemp[objectNumber]).end();
1730            ++currentTriangle1) {
1731    
1732           const SurfaceMeshGenerator::Internals::TriangleVertices& tv=(*(*currentTriangle1).second);
1733            
1734            PairTriangleEdges pT;
1735            if(numobject<objectNumber) {
1736                pT.first=&(currentTriangle);
1737                pT.second=&(tv);
1738            }
1739            else {
1740                pT.first=&(tv);
1741                pT.second=&(currentTriangle);
1742            }
1743            size_t temp=0;
1744            if((listVertexIntersectionPerPairEdges[objectp]).find(pT)!=
1745        (listVertexIntersectionPerPairEdges[objectp]).end()) {
1746                TinyVector<3,size_t> cc;
1747                if(numobject<objectNumber) {
1748                    cc=((listVertexIntersectionPerPairEdges[objectp])[pT]).cutTriangle1;
1749                }
1750                else {
1751                    cc=((listVertexIntersectionPerPairEdges[objectp])[pT]).cutTriangle2;
1752                }
1753                if(cc[0]!=0) {
1754                    temp+=1;
1755                    tt[0]=1;
1756                    pointsIntersection[0]=
1757                        ((listVertexIntersectionPerPairEdges[objectp])[pT].listofPoints[cc[0]-1]);
1758                }
1759                if(cc[1]!=0) {
1760                    temp+=1;
1761                    tt[1]=1;
1762                    pointsIntersection[1]=
1763                        ((listVertexIntersectionPerPairEdges[objectp])[pT].listofPoints[cc[1]-1]);
1764                }
1765                if(cc[2]!=0) {
1766                    temp+=1;
1767                    tt[2]=1;
1768                    pointsIntersection[2]=
1769                        ((listVertexIntersectionPerPairEdges[objectp])[pT].listofPoints[cc[2]-1]);
1770                }
1771                if(temp<(listVertexIntersectionPerPairEdges[objectp])[pT].size()) {
1772                    numberin+=1;
1773                    size_t r;
1774                    size_t sum=cc[0]+cc[1]+cc[2];
1775                    //ffout(0)<<(listVertexIntersectionPerPairEdges[objectp])[pT].size()<<" "<<sum<<"\n";
1776                    if(sum==0 or sum==2 or sum==5)
1777                        r=0;
1778                    else if(sum==4 or sum==1)
1779                        r=1;
1780                    else if(sum==3)
1781                        r=2;
1782                    else {
1783                        ffout(0)<<"pbs sum = "<<sum<<" size "<<(listVertexIntersectionPerPairEdges[objectp])[pT].size()<<" \n";
1784                        ffout(0)<<numobject<<" "<<objectNumber<<"\n";
1785                        ffout(0)<<cc[0]<<" "<<cc[1]<<" "<<cc[2]<<"\n";
1786                    }
1787                    pointsin.push_back((listVertexIntersectionPerPairEdges[objectp])[pT].listofPoints[r]);
1788              }              }
1789            /*      if(re0==2 and (ttn1==1 or ttn2==1))              //ffout(0)<<cc[0]<<" "<<cc[1]<<" "<<cc[2]<<"\n";
1790                    {       //ffout(0)<<((listVertexIntersectionPerPair[objectp])[pT]).size()<<"\n";;
                   re0=-nn;  
                   ii0=-nn;  
                   }  
                   else  
                   {  
                   re1=-nn;  
                   ii1=-nn;  
                   }*/  
           //  re0=0;  
           //re1=0;  
1791          }          }
1792          else          else {
1793          {             // ffout(0) <<"pasintersecte\n";
           //std::cout<<"-- 0 2 on sait traite\n";  
             ii0=re0;  
             ii1=re1;  
1794          }          }
1795      }      }
1796      else if((re0+re1==3))  }
1797      {  
1798          if(ttn0==0)  TinyVector<3,int> SurfaceMeshGenerator::Internals::
1799          {  __determinateEdgeCut(const size_t& numobject,
1800            //std::cout<<"-- 1 2 pas coupe\n";                       const size_t & objectNumber,
1801              re0=2;                       const TinyVector<3,TinyVector<2,Edge::Pair> >& currentTriangle){
1802              re1=2;  
1803        TinyVector<3,int> tt;
1804        tt=0;
1805        ObjectPair objectp(numobject,objectNumber);
1806        if(numobject>objectNumber) {
1807            objectp.first=objectNumber;
1808            objectp.second=numobject;
1809        }
1810    
1811        for(SurfaceMeshGenerator::Internals::TriangleVerticesMap::iterator currentTriangle1
1812            = (*triangleListesIntersectionEdgesTemp[objectNumber]).begin();
1813            currentTriangle1 != (*triangleListesIntersectionEdgesTemp[objectNumber]).end();
1814            ++currentTriangle1) {
1815    
1816            const SurfaceMeshGenerator::Internals::TriangleVertices& tv=(*(*currentTriangle1).second);
1817    
1818            PairTriangleEdges pT;
1819            if(numobject<objectNumber) {
1820                pT.first=&(currentTriangle);
1821                pT.second=&(tv);
1822          }          }
1823          else          else {
1824          {              pT.first=&(tv);
1825            //std::cout<<"-- 1 2 on sait traite\n";              pT.second=&(currentTriangle);
1826              if(re0==1)          }
1827              {          if((listVertexIntersectionPerPairEdges[objectp]).find(pT)!=
1828                  ii0=0;      (listVertexIntersectionPerPairEdges[objectp]).end()) {
1829                  ii1=1;              TinyVector<3,size_t> cc;
1830                  re0=0;              if(numobject<objectNumber) {
1831                    cc=((listVertexIntersectionPerPairEdges[objectp])[pT]).cutTriangle1;
1832              }              }
1833              else              else {
1834              {                  cc=((listVertexIntersectionPerPairEdges[objectp])[pT]).cutTriangle2;
1835                  ii0=1;              }
1836                  ii1=0;              if(cc[0]!=0) {
1837                  re1=0;                  tt[0]=1;
1838                }
1839                if(cc[1]!=0) {
1840                    tt[1]=1;
1841                }
1842                if(cc[2]!=0) {
1843                    tt[2]=1;
1844              }              }
1845          }          }
1846      }      }
1847    //     else      return tt;
1848    //std::cout<<"--0 out 0 out rien a faire\n";  }
1849    
1850      if(re0==1 and re2==1)  TinyVector<3,size_t> SurfaceMeshGenerator::Internals::
1851      {  __putEdgeRef(const size_t& ncase,
1852          if(ttn1!=0 and iter>1)               const TinyVector<3,TinyVector<2,Edge::Pair> >& currentTriangle,
1853  std::cout<<"--* cas indetermine\n";               TriangleCut& K) {
1854    
1855        const TinyVector<3,int>  tt=K.edgecut;
1856        TinyVector<3,size_t>& ii=K.isInTriangle;
1857        size_t numobject=K.numobject;
1858    
1859        const TinyVector<2,Edge::Pair>& e0=(currentTriangle)[0];
1860        const TinyVector<2,Edge::Pair>& e1=(currentTriangle)[1];
1861        const TinyVector<2,Edge::Pair>& e2=(currentTriangle)[2];
1862    
1863        if(tt[0]+tt[1]+tt[2]!=0 and ((*edgesRefCreate[numobject])[e0[0]]+
1864                                     (*edgesRefCreate[numobject])[e1[0]]+(*edgesRefCreate[numobject])[e2[0]])==6) {
1865            ffout(0)<<"************************ ---------tr a pbs\n";
1866            printpoint(numobject,e0,K);
1867            printpoint(numobject,e1,K);
1868            printpoint(numobject,e2,K);
1869            ffout(0)<<tt[0]<<" "<<tt[1]<<" "<<tt[2]<<"\n";
1870      }      }
1871      else if((re0==1 and re2==0) or      
1872              (re0==0 and re2==1))      __inout(numobject,tt[0],e0,e1,ii[0],ii[1],K);
1873      {      __inout(numobject,tt[1],e0,e2,ii[0],ii[2],K);
1874          if(ttn1==0)      __inout(numobject,tt[2],e1,e2,ii[1],ii[2],K);
1875          {      if(ncase==0) {
1876              //std::cout<<"--* 1 O non inter\n";          if(tt[0]+tt[1]+tt[2]==0 and ((*edgesRefCreate[numobject])[e0[0]]+
1877              re0=0;                                (*edgesRefCreate[numobject])[e1[0]]+(*edgesRefCreate[numobject])[e2[0]])>=4) {
1878              re2=0;              //if triangle pas plat.........
1879          }              if(!(triangleplat(numobject,e0,e1,e2)) ) {
1880          else                  (*edgesRefCreate[numobject])[e0[0]]=2;
1881          {                  (*edgesRefCreate[numobject])[e1[0]]=2;
1882            //std::cout<<"--* 1 in 0 out on sait traite\n";                  (*edgesRefCreate[numobject])[e2[0]]=2;
1883              ii0=re0;              }
1884              ii2=re2;              else {
1885          }                  real_t eps=1e-13;
1886      }                  Vertex V0,V1,V2;
1887      else if((re0==2 and re2==0) or                  V0=(*(*verticesListesIntersection[numobject])[e0[0]]);
1888              (re0==0 and re2==2))                  V1=(*(*verticesListesIntersection[numobject])[e1[0]]);
1889      {                  V2=(*(*verticesListesIntersection[numobject])[e2[0]]);
1890          if(ttn1==0)                  real_t normV0V1=Norm(V0.Vector3()-V1.Vector3());
1891          {                  real_t normV0V2=Norm(V0.Vector3()-V2.Vector3());
1892              if(ttn0+ttn1+ttn2==0)                  real_t normV1V2=Norm(V2.Vector3()-V1.Vector3());
1893              {                  if(normV0V1>eps and (*edgesRefCreate[numobject])[e0[0]]+
1894                  re0=2;       (*edgesRefCreate[numobject])[e1[0]]==4)
1895                  re2=2;                      (*edgesRefCreate[numobject])[e2[0]]=2;
1896                    if(normV0V2>eps and (*edgesRefCreate[numobject])[e0[0]]+
1897         (*edgesRefCreate[numobject])[e2[0]]==4)
1898                        (*edgesRefCreate[numobject])[e1[0]]=2;
1899                    if(normV1V2>eps and (*edgesRefCreate[numobject])[e1[0]]+
1900         (*edgesRefCreate[numobject])[e2[0]]==4)
1901                        (*edgesRefCreate[numobject])[e0[0]]=2;
1902                    ffout(0)<<"triangle plat\n";
1903              }              }
1904              else          }      
1905              {      }
1906        else if(ncase==1) {
1907        
1908            if((*edgesRefCreate[numobject])[e0[0]]==2)
1909                ii[0]=1;
1910            if((*edgesRefCreate[numobject])[e1[0]]==2)
1911                ii[1]=1;
1912            if((*edgesRefCreate[numobject])[e2[0]]==2)
1913                ii[2]=1;
1914        }
1915      
1916    
1917  std::cout<<"--* 0 2 pas coupe\n";      return ii;
1918  std::cout<<ttn0<<" "<<ttn1<<" "<<ttn2<<"\n";  }
 std::cout<<re0<<" "<<re2<<"\n";  
                 re0=2;  
                 re2=2;  
             }  
           /*   if(re0==2 and (ttn0==1 or ttn2==1))  
                {  
                re0=-nn;  
                ii0=-nn;  
                }  
                else  
                {  
                re2=-nn;  
                ii2=-nn;  
                }*/  
1919    
1920            //      re0=0;  void SurfaceMeshGenerator::Internals::
1921            //    re2=0;  __createListTriangle(const TriangleCut& K,
1922                         size_t & kill){
1923        const TinyVector<3,int>&  tt=K.edgecut;
1924        const TinyVector<3,size_t>& ii=K.isInTriangle;
1925        const size_t& numobject=K.numobject;
1926        const TinyVector<2,Edge::Pair>& e0=K.listpairedges[0];
1927        const TinyVector<2,Edge::Pair>& e1=K.listpairedges[1];
1928        const TinyVector<2,Edge::Pair>& e2=K.listpairedges[2];
1929        TinyVector<2,Edge::Pair> pe0(e0),pe1(e1),pe2(e2);
1930    
1931       if(kill==3 and ((tt[0]+tt[1]+tt[2]==0 and ii[0]+ii[1]+ii[2]==3) or ii[0]+ii[1]+ii[2]==3)) {
1932           TinyVector<2,Edge::Pair> num0=pe0,num1=pe1,num2=pe2;
1933            __addTriangle(numobject,num0,num1,num2,kill);
1934        }
1935        else if(ii[0]+ii[1]+ii[2]==1) {
1936            if(ii[0]==1) {
1937                __createCase1(0,K,kill);
1938            }
1939            if(ii[1]==1) {
1940                __createCase1(1,K,kill);
1941            }
1942            if(ii[2]==1) {
1943                __createCase1(2,K,kill);
1944            }
1945        }
1946        else if(ii[0]+ii[1]+ii[2]==2) {
1947            if(ii[0]+ii[1]==2) {
1948                __createCase2(0,K,kill);
1949          }          }
1950          else          if(ii[0]+ii[2]==2) {
1951          {              __createCase2(1,K,kill);
           //std::cout<<"--* 0 2 on sait traite\n";  
             ii0=re0;  
             ii2=re2;  
1952          }          }
1953      }          if(ii[1]+ii[2]==2) {
1954      else if((re0+re2==3))              __createCase2(2,K,kill);
     {  
         if(ttn1==0)  
         {  
           //std::cout<<"--* 1 2 pas coupe\n";  
             re0=2;  
             re2=2;  
1955          }          }
1956          else      }
1957          {  }
           //std::cout<<"--* 1 2 on sait traite\n";  
             if(re0==1)  
             {  
                 ii0=0;  
                 ii2=1;  
                 re0=0;  
             }  
             else  
             {  
                 ii0=1;  
                 ii2=0;  
                 re2=0;  
             }  
         }        
     }        
   //       else  
   //std::cout<<"--* 0 out 0 out rien a faire\n";  
1958    
1959      if(re2==1 and re1==1)  void SurfaceMeshGenerator::Internals::
1960      {  __createTrianglesIntersection(const size_t& numobject,
1961          if(ttn2!=0 and iter>1)                                const size_t& objectNumber,
1962  std::cout<<"--+ cas indetermine\n";                                const size_t& nobject) {
1963    
1964        ffout(0)<<"numero "<<numobject<<" "<<objectNumber<<"\n";
1965    
1966        for(SurfaceMeshGenerator::Internals::TriangleVerticesMap::iterator currentTriangle
1967            = (*triangleListesIntersectionEdgesTemp[numobject]).begin();
1968            currentTriangle != (*triangleListesIntersectionEdgesTemp[numobject]).end();
1969            ++currentTriangle) {
1970    
1971            const SurfaceMeshGenerator::Internals::TriangleVertices& tv=(*(*currentTriangle).second);
1972            const SurfaceMeshGenerator::Internals::TriangleVertices& tv2=(*(*currentTriangle).first);
1973            const TinyVector<2,Edge::Pair>& e0=tv2[0];
1974            const TinyVector<2,Edge::Pair>& e1=tv2[1];
1975            const TinyVector<2,Edge::Pair>& e2=tv2[2];
1976            
1977            (*edgesRefCreate[numobject])[e0[0]]=0;
1978            (*edgesRefCreate[numobject])[e1[0]]=0;
1979            (*edgesRefCreate[numobject])[e2[0]]=0;
1980      }      }
1981      else if((re2+re1==1))  
1982      {      for(SurfaceMeshGenerator::Internals::EdgePairList::const_iterator currentTriangle
1983          if(ttn2==0)          = (*triangleListes[numobject]).begin();
1984          {          currentTriangle != (*triangleListes[numobject]).end();
1985            //std::cout<<"--+ 1 0 non inter\n";          ++currentTriangle) {
1986              re2=0;          const Edge& e0=(*currentTriangle)[0];
1987              re1=0;          const Edge& e1=(*currentTriangle)[1];
1988          }          const Edge& e2=(*currentTriangle)[2];
1989          else          (*edgesRefCreate[numobject])[e0]=2;
1990          {          (*edgesRefCreate[numobject])[e1]=2;
1991            //std::cout<<"--+ 1 in 0 out on sait traite\n";          (*edgesRefCreate[numobject])[e2]=2;
1992              ii1=re1;      }
1993              ii2=re2;  
1994        ObjectPair objectp(numobject,objectNumber);
1995        if(numobject>objectNumber) {
1996            objectp.first=objectNumber;
1997            objectp.second=numobject;
1998        }
1999        size_t nbinters=listVertexIntersectionPerPairEdges[objectp].size();
2000        ffout(0)<<"size listVertexIntersectionPerPair "<<nbinters<<"\n";
2001    
2002        size_t dd_max=2;
2003        for(size_t dd=0 ; dd<dd_max ; ++dd) {
2004            ffout(0)<<"num passage "<<dd<<"\n";
2005            for(SurfaceMeshGenerator::Internals::TriangleVerticesMap::iterator currentTriangle
2006         = (*triangleListesIntersectionEdgesTemp[numobject]).begin();
2007         currentTriangle != (*triangleListesIntersectionEdgesTemp[numobject]).end();
2008         ++currentTriangle) {
2009    
2010                const SurfaceMeshGenerator::Internals::TriangleVertices& tv=(*(*currentTriangle).second);
2011                const SurfaceMeshGenerator::Internals::TriangleVertices& tv2=(*(*currentTriangle).first);
2012                const TinyVector<2,Edge::Pair>& e0=tv2[0];
2013                const TinyVector<2,Edge::Pair>& e1=tv2[1];
2014                const TinyVector<2,Edge::Pair>& e2=tv2[2];
2015                
2016                TriangleCut K(e0,e1,e2,numobject,objectNumber);
2017                //si c un tr pas encore traite
2018                if(tv2[0][0]==tv2[0][1]
2019            and tv2[1][0]==tv2[1][1]
2020            and tv2[2][0]==tv2[2][1]) {
2021                    K.edgecut =__determinateEdgeCut(numobject,objectNumber,(*(*currentTriangle).second));
2022                    K.isInTriangle =__putEdgeRef(0,(*(*currentTriangle).second),K);
2023                }
2024                else {
2025                    ffout(0)<<"----------------------------- point intersection1\n";
2026                    K.edgecut =__determinateEdgeCut(numobject,objectNumber,(*(*currentTriangle).second));
2027                    K.isInTriangle =__putEdgeRef(0,(*(*currentTriangle).first),K);
2028                }
2029          }          }
2030      }      }
     else if((re2==2 and re1==0) or  
             (re2==0 and re1==2))  
     {  
         if(ttn2==0)  
         {  
             if(ttn0+ttn1+ttn2==0)  
             {  
                 re1=2;  
                 re2=2;  
             }  
             else  
             {  
2031    
2032  std::cout<<"--+ 0 2 pas coupe\n";      size_t kill=0;
2033  std::cout<<ttn0<<" "<<ttn1<<" "<<ttn2<<"\n";  
2034  std::cout<<re1<<" "<<re2<<"\n";      for(SurfaceMeshGenerator::Internals::TriangleVerticesMap::iterator currentTriangle
2035                  re1=2;          = (*triangleListesIntersectionEdgesTemp[numobject]).begin();
2036                  re2=2;          currentTriangle != (*triangleListesIntersectionEdgesTemp[numobject]).end();
2037              }          ++currentTriangle) {
2038            /*   if(re1==2 and (ttn1==1 or ttn0==1))  
2039                 {          const SurfaceMeshGenerator::Internals::TriangleVertices& tv=(*(*currentTriangle).second);
2040                 re1=-nn;          const SurfaceMeshGenerator::Internals::TriangleVertices& tv2=(*(*currentTriangle).first);
2041                 ii1=-nn;          TinyVector<2,Edge::Pair> e0,e1,e2;
2042                 }          if(tv2[0][0]==tv2[0][1]) {
2043                 else              e0=tv[0];}
2044                 {          else {
2045                 re2=-nn;              e0=tv2[0];
2046                 ii2=-nn;          }
2047                 }*/          if(tv2[1][0]==tv2[1][1]) {
2048            //re1=0;              e1=tv[1];}
2049            //re2=0;          else {
2050          }              e1=tv2[1];
2051          else          }
2052          {          if(tv2[2][0]==tv2[2][1]) {
2053            //std::cout<<"--+ 0 2 on sait traite\n";              e2=tv[2];}
2054              ii1=re1;          else {
2055              ii2=re2;              e2=tv2[2];
2056          }          }
2057            /*const TinyVector<2,Edge::Pair>& e0=tv2[0];
2058            const TinyVector<2,Edge::Pair>& e1=tv2[1];
2059            const TinyVector<2,Edge::Pair>& e2=tv2[2];*/
2060    
2061            TriangleCut K(e0,e1,e2,numobject,objectNumber);
2062            
2063            if(tv2[0][0]==tv2[0][1]
2064        and tv2[1][0]==tv2[1][1]
2065        and tv2[2][0]==tv2[2][1]) {
2066    
2067                __createLocalListIntersection(objectNumber,(*(*currentTriangle).second),K);
2068                K.isInTriangle =__putEdgeRef(1,(*(*currentTriangle).second),K);
2069            }
2070            else {
2071                ffout(0)<<"----------------------------- point intersection2\n";
2072                __createLocalListIntersection(objectNumber,(*(*currentTriangle).second),K);
2073                K.isInTriangle =__putEdgeRef(1,(*(*currentTriangle).first),K);
2074            }
2075                kill=0;
2076                if(objectNumber==nobject-1 or (objectNumber==nobject-2 and numobject==nobject-1))
2077                    kill=3;
2078                else if(nbinters==0)
2079                    kill=4;
2080    
2081                if(kill!=4) {
2082                    //if(!(triangleplat(numobject,e0,e1,e2)))
2083                        __createListTriangle(K,kill);
2084                }
2085                if(kill==1) {
2086                    //ffout(0)<<"construction d'un triangle\n";
2087                    if(numobject<objectNumber) {
2088                        triangleEraseEdgesTemp[0].push_back(currentTriangle);
2089                        for(size_t kk=0 ; kk<localListTriangleintersection.size() ; ++kk) {
2090                            triangleTreat[0][&tv].push_back(localListTriangleintersection[kk]);
2091                        }
2092                    }
2093                    else {
2094                        triangleEraseEdgesTemp[1].push_back(currentTriangle);
2095                        for(size_t kk=0 ; kk<localListTriangleintersection.size() ; ++kk) {
2096                            triangleTreat[1][&tv].push_back(localListTriangleintersection[kk]);
2097                        }
2098                    }
2099                }
2100                localListTriangleintersection.clear();
2101            
2102      }      }
2103      else if((re2+re1==3))      ffout(0)<<" triangleListesIntersection size "<<(*triangleListesIntersectionEdges[numobject]).size()<<"\n";
2104      {  }
2105          if(ttn2==0)  
2106          {  bool SurfaceMeshGenerator::Internals::
2107              //std::cout<<"--+ 1 2 pas coupe\n";  __createTriangleSurface(const size_t& nobject,
2108              re2=2;                          const Shape& S ,
2109              re1=2;                          const Domain& Omega,
2110          }                          std::vector<const Edge*>& cutEdges,
2111          else                          const Reference::Representation& ref,
2112          {                          const Cell& currentCell){
2113            //std::cout<<"--+ 1 2 on sait traite\n";      bool jobDone=true;
2114              if(re1==1)      switch (cutEdges.size()) {
2115              {          case 3: {
2116                  ii1=0;              //! Stores the Edges on which the vertex will be computed.
2117                  ii2=1;              
2118                  re1=0;              for (size_t i=0; i<cutEdges.size(); ++i) {
2119                    TinyVector<3> V=__splitEdge(*cutEdges[i], ref, S);
2120                    if(Omega.inside(V)) {
2121                        (*verticesListes[nobject])[*cutEdges[i]]
2122                            =  new Vertex(V);
2123                        (*edgesRef[nobject])[*cutEdges[i]] = 1;
2124                    }
2125                    else {
2126                       // ffout(0)<<"outside\n";
2127                        jobDone=false;
2128                        break;
2129                    }
2130              }              }
2131              else              //! Addes a triangle to the list.
2132              {              localTriangleList.push_front(TinyVector<3,Edge::Pair>());
2133                  ii1=1;              cellList.push_front(&currentCell);
2134                  ii2=0;  
2135                  re2=0;              std::vector<TinyVector<3> > V(cutEdges.size());
2136                for(size_t i=0; i<cutEdges.size(); ++i) {
2137                    const Edge& e = *cutEdges[i];
2138                    if(e(0).Ref().isSet(ref))
2139                        V[i] = (e(0).Vector3() - e(1).Vector3());
2140                    else
2141                        V[i] = (e(1).Vector3() - e(0).Vector3());
2142              }              }
2143    
2144                if ((V[0]^V[1])*V[2] < 0)
2145                     std::swap(cutEdges[1], cutEdges[2]);
2146    
2147                //! Sets the triangle vertices (actually just store the edge where it lives).
2148                SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle
2149                    = localTriangleList.begin();
2150    
2151                for(size_t i=0; i<cutEdges.size(); ++i)
2152                    (*currentTriangle)[i] = (Edge::Pair)(*cutEdges[i]);
2153                
2154                break;
2155          }          }
2156      }          case 4: {
   //       else  
   //std::cout<<"--+ 0 out 0 out rien a faire\n";  
2157    
2158      if(/*iter>1 and*/ (((re2==1 and re1==1 and re0==1) or re0+re1+re2==6)              for (size_t i=0; i<4; ++i) {
2159                         and ttn0+ttn1+ttn2==0))                  TinyVector<3> V=__splitEdge(*cutEdges[i], ref, S);
2160      {                  if(Omega.inside(V)) {
2161          re0=2;                      (*verticesListes[nobject])[*cutEdges[i]]
2162          re1=2;                      =  new Vertex(V);
2163          re2=2;                      (*edgesRef[nobject])[*cutEdges[i]] = 1;
2164          ii0=1;                  }
2165          ii1=1;                  else {
2166          ii2=1;                      //ffout(0)<<"outside4\n";
2167      }                      jobDone=false;
2168      if(iter>2 and re0+re1+re2==0 and ttn0+ttn1+ttn2!=0)                      break;
2169      {                  }
2170  std::cout<<"pbs de precision?? tr out cut\n";              }
2171          ii0=0;              //! Orders the Edges using a map !
2172          ii1=0;              std::map<Edge::Pair, const Edge*> cutEdgesMap;
2173          ii2=0;              for(size_t ie=0; ie<cutEdges.size(); ++ie)
2174                    cutEdgesMap[(Edge::Pair)*(cutEdges[ie])] = (cutEdges[ie]);
2175    
2176                //! Addes a triangle to the list.
2177                localTriangleList.push_front(TinyVector<3,Edge::Pair>());
2178                cellList.push_front(&currentCell);
2179    
2180                std::vector<TinyVector<3> > V(cutEdges.size());
2181                std::map<Edge::Pair, const Edge*>::iterator currentCutEdge=cutEdgesMap.begin();
2182                for (size_t i=0; i<V.size(); ++i) {
2183                    const Edge& e = *(*currentCutEdge).second;
2184                    if(e(0).Ref().isSet(ref))
2185                        V[i] = (e(0).Vector3() - e(1).Vector3());
2186                    else
2187                        V[i] = (e(1).Vector3() - e(0).Vector3());
2188                    ++currentCutEdge;
2189                }
2190    
2191                bool swapEdge = ((V[0]^V[1])*V[2] > 0);
2192    
2193                //! Sets the triangle vertices (actually just store the edge where it lives).
2194                SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle
2195                    = localTriangleList.begin();
2196    
2197                currentCutEdge = cutEdgesMap.begin();
2198                if (!swapEdge)
2199                    for(size_t i=0; i<3; ++i) {
2200                        (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
2201                        ++currentCutEdge;
2202                    }
2203                else
2204                    for(int i=2; i>=0; --i) {
2205                        (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
2206                        ++currentCutEdge;
2207                    }
2208    
2209                
2210    
2211                //! Addes a triangle to the list.
2212                localTriangleList.push_front(TinyVector<3,Edge::Pair>());
2213                cellList.push_front(&currentCell);
2214    
2215                currentTriangle = localTriangleList.begin();
2216    
2217                //! Sets the triangle vertices (actually just store the edge where it lives).
2218                currentCutEdge = cutEdgesMap.begin();
2219                if (!swapEdge)
2220                    for(int i=2; i>=0; --i) {
2221                        ++currentCutEdge;
2222                        (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
2223                    }
2224                else
2225                    for(size_t i=0; i<3; ++i) {
2226                        ++currentCutEdge;
2227                        (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);
2228                    }
2229                break;
2230            }
2231            case 1:
2232            case 2: {
2233                jobDone = false;
2234                break;
2235            }
2236      }      }
   /*if(iter==1)  
     {  
       if(ttn0+ttn1+ttn2==0 and ((re0==0 and re1==1 and re2==1) or (re0==1 and re1==0 and re2==1)  
                                       or (re0==1 and re1==1 and re2==0) or (re0==0 and re1==0 and re2==1)  
                                       or (re0==0 and re1==1 and re2==0) or (re0==1 and re1==0 and re2==0)))  
         {  
           re0=0;  
           re1=0;  
           re2=0;  
         }  
     }*/  
   
   if(ii0==2)  
     ii0=1;  
   if(ii1==2)  
     ii1=1;  
   if(ii2==2)  
     ii2=1;  
2237    
2238        return jobDone;
2239        
2240  }  }
2241    
2242  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
2243  __IntersectionPoint(const Scene& scene ,const Reference::Representation& ref,  __createSurface(Structured3DMesh& SMesh,const Domain& Omega,
2244                      const Reference::Representation& otherRef,                  const Shape& S ,
2245                      std::map<Edge::Pair,int>& triangleedge,                  const Reference::Representation& ref,
2246                      std::map<Edge::Pair, Vertex*>&  listofVertex)                  size_t nobject,
2247  {                  std::set<const Cell*>& cuttedHexahedra,
2248                    std::set<const Cell*>& meshedHexahedra) {
2249    
2250    int nn=0;    for (size_t i = 0; i<SMesh.nbEdges(); i++) {
2251    EdgePairList::iterator currentTriangle=(*triangleListes[0]).end();      Edge& e = SMesh.edge(i);
2252    EdgePairList::iterator currentTriangleother=(*triangleListes[1]).end();  #warning THIS TEST IS NOT ENOUGH! THINK TO THE BORDER!
2253    --currentTriangle;      if(/*(e(0).Ref().insideDomain() xor e(1).Ref().insideDomain())
2254    --currentTriangleother;         and*/((e(0).Ref().isSet(ref)) xor (e(1).Ref().isSet(ref)))) {
2255          (*verticesListes[nobject])[e] = new Vertex(__splitEdge(e, ref, S));
2256          (*edgesRef[nobject])[e] = 1;
2257                
2258        }
2259      }
2260            
2261    for(SurfaceMeshGenerator::Internals::HexaTetraAssociation::iterator i    if ((*verticesListes[nobject]).size()==0) {
         = hexaTetraAssociation.begin();  
       i != hexaTetraAssociation.end(); ++i) {  
2262    
2263            fferr(2) << __FILE__ << ':' << __LINE__ << ": "
2264      std::vector<TetrahedronByEdges>& T = *(((*i).second).first);               << "meshing object " << S << " a problem occured ...\n"
2265  //      ffout(0) << " hexa*****************************************\n";               << "Check that:\n"
2266                 << "- the object is not to small compare to the background mesh\n"
2267                 << "- the object intersectes the background mesh\n\n";
2268        std::exit(1);
2269      }
2270    
2271      for (std::vector<TetrahedronByEdges>::iterator currentT = T.begin();    for (size_t icell=0; icell<SMesh.shape().nx()-1; ++icell)
2272           currentT < T.end(); ++currentT) {      for (size_t jcell=0; jcell<SMesh.shape().ny()-1; ++jcell)
2273          for (size_t kcell=0; kcell<SMesh.shape().nz()-1; ++kcell) {
2274            bool intersected = false;
2275    
2276        std::vector<const Edge*> cutEdges;          const Cell& currentCell = SMesh.cell(icell, jcell, kcell);
2277        std::vector<const Edge*> cutEdgesother;          const size_t nx
2278              = dynamic_cast<const Structured3DMesh&>(SMesh).shape().nx();
2279            const size_t nx_1 = nx-1;
2280    
2281        for(size_t ie=0; ie<6; ++ie) {          const size_t ny
2282          const Edge& e = currentT->edge(ie);            = dynamic_cast<const Structured3DMesh&>(SMesh).shape().ny();
2283          if ((e(0).Ref().isSet(ref)) xor (e(1).Ref().isSet(ref))) {          const size_t ny_1 = ny-1;
2284            cutEdges.push_back(&e);  
2285          }          const size_t nz
2286          if ((e(0).Ref().isSet(otherRef)) xor (e(1).Ref().isSet(otherRef))) {            = dynamic_cast<const Structured3DMesh&>(SMesh).shape().nz();
2287            cutEdgesother.push_back(&e);          const size_t nz_1 = nz-1;
2288    
2289            HexahedronByEdges H (SMesh.edge(icell*ny*nz + jcell*nz + kcell),
2290                                 SMesh.edge(icell*ny*nz + (jcell+1)*nz + kcell),
2291                                 SMesh.edge(icell*ny*nz + jcell*nz + kcell + 1),
2292                                 SMesh.edge(icell*ny*nz + (jcell+1)*nz + kcell+1),
2293                                 SMesh.edge(nx_1*ny*nz + icell*ny_1*nz + jcell*nz + kcell),
2294                                 SMesh.edge(nx_1*ny*nz + (icell+1)*ny_1*nz + jcell*nz + kcell),
2295                                 SMesh.edge(nx_1*ny*nz + icell*ny_1*nz + jcell*nz + kcell + 1),
2296                                 SMesh.edge(nx_1*ny*nz + (icell+1)*ny_1*nz + jcell*nz + kcell+1),
2297                                 SMesh.edge(nx_1*ny*nz + nx*ny_1*nz + icell*ny*nz_1 + jcell*nz_1 + kcell ),
2298                                 SMesh.edge(nx_1*ny*nz + nx*ny_1*nz + (icell+1)*ny*nz_1 + jcell*nz_1 + kcell ),
2299                                 SMesh.edge(nx_1*ny*nz + nx*ny_1*nz + icell*ny*nz_1 + (jcell+1)*nz_1 + kcell ),
2300                                 SMesh.edge(nx_1*ny*nz + nx*ny_1*nz + (icell+1)*ny*nz_1 + (jcell+1)*nz_1 + kcell ));
2301    
2302    
2303            for (size_t l=0; l<H.numberOfEdges(); l++) {
2304              Edge& e = H.edge(l);
2305              if ((*verticesListes[nobject]).count(e) > 0) {
2306                  intersected = true;
2307                  cuttedHexahedra.insert(&currentCell);
2308                break;
2309              }
2310          }          }
       }  
 //          ffout(0) <<cutEdges.size()<<" "<<cutEdgesother.size()<< " tetra__________*****************************************\n";  
2311    
2312        switch (cutEdges.size()) {          if (intersected) { // Means that the hexahedron intersectes the object.
2313        case 3: {            int sens = (icell+jcell+kcell)%2;
2314          localTriangleList.push_front(TinyVector<3,Edge::Pair>());            ReferenceCounting<std::vector<Edge> > pAddedEdges;
2315          EdgePairList::iterator currentT = localTriangleList.begin();            pAddedEdges = new std::vector<Edge>(6);
2316          for(size_t i=0; i<3; ++i)            std::vector<Edge>& addedEdges = *pAddedEdges;
           (*currentT)[i] = (*currentTriangle)[i];  
         --currentTriangle;  
         break;  
       }  
       case 4: {  
         {  
           localTriangleList.push_front(TinyVector<3,Edge::Pair>());  
           EdgePairList::iterator currentT = localTriangleList.begin();  
           for(size_t i=0; i<3; ++i)  
             (*currentT)[i] = (*currentTriangle)[i];  
         }  
         --currentTriangle;  
         localTriangleList.push_front(TinyVector<3,Edge::Pair>());  
         EdgePairList::iterator currentT = localTriangleList.begin();  
   
         for(size_t i=0; i<3; ++i)  
           (*currentT)[i] = (*currentTriangle)[i];  
         --currentTriangle;  
         break;  
       }  
       case 1:  
       case 2: {  
         ffout(0) << "?!\n";  
         break;  
       }  
       }  
2317    
2318        switch (cutEdgesother.size()) {            std::vector<int> TverticesList(4);
2319        case 3: {            switch (sens) {
2320          localTriangleListother.push_front(TinyVector<3,Edge::Pair>());            case 0: {
2321          EdgePairList::iterator currentT = localTriangleListother.begin();              TverticesList[0] = 1;
2322          for(size_t i=0; i<3; ++i)              TverticesList[1] = 3;
2323            (*currentT)[i] = (*currentTriangleother)[i];              TverticesList[2] = 4;
2324          --currentTriangleother;              TverticesList[3] = 6;
2325          break;              break;
2326        }            }
2327        case 4: {            case 1: {
2328          {              TverticesList[0] = 0;
2329            localTriangleListother.push_front(TinyVector<3,Edge::Pair>());              TverticesList[1] = 2;
2330            EdgePairList::iterator currentT = localTriangleListother.begin();              TverticesList[2] = 5;
2331            for(size_t i=0; i<3; ++i)              TverticesList[3] = 7;
2332              (*currentT)[i] = (*currentTriangleother)[i];            }
2333          }            }
         --currentTriangleother;  
   
         localTriangleListother.push_front(TinyVector<3,Edge::Pair>());  
         EdgePairList::iterator currentT = localTriangleListother.begin();  
   
         for(size_t i=0; i<3; ++i)  
           (*currentT)[i] = (*currentTriangleother)[i];  
         --currentTriangleother;  
         break;  
       }  
       case 1:  
       case 2: {  
         ffout(0) << "?!\n";  
         break;  
       }  
       }  
     }  
2334    
2335              ReferenceCounting<std::vector<TetrahedronByEdges> > pT;
2336              pT = new std::vector<TetrahedronByEdges>(5);
2337    
2338      for(SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle1            std::vector<TetrahedronByEdges>& T = *pT;
2339            = localTriangleList.begin();  
2340          currentTriangle1 != localTriangleList.end();            { // Create tetrahedra edge
2341          ++currentTriangle1) {              std::vector<Edge>::iterator currentAddedEdge
2342        const Edge::Pair& e0=(*currentTriangle1)[0];                = addedEdges.begin();
2343        const Edge::Pair& e1=(*currentTriangle1)[1];              for (size_t i1=0; i1<3; ++i1) {
2344        const Edge::Pair& e2=(*currentTriangle1)[2];                for (size_t i2=i1; i2<3; ++i2) {
2345        //ffout(0) << (*e0.first).Vector3()[0]<<"\n";                  *(currentAddedEdge) = Edge(currentCell(TverticesList[i1]),
2346        const Vertex& V0 = *((*verticesListes[0])[e0]);                                             currentCell(TverticesList[i2+1]));
2347        const Vertex& V1 = *((*verticesListes[0])[e1]);                  Edge& e = (*currentAddedEdge);
2348        const Vertex& V2 = *((*verticesListes[0])[e2]);                  if (((e(0).Ref().isSet(ref)) xor (e(1).Ref().isSet(ref)))
2349        //tt[nn] = 0;                      /*and (e(0).Ref().insideDomain() xor e(1).Ref().insideDomain())*/) {
2350        (triangleedge)[e0]=0;                    e.Ref() = 1;
2351        (triangleedge)[e1]=0;                    (*edgesRef[nobject])[e] = 1;
2352        (triangleedge)[e2]=0;                  }
2353        for(SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle                  currentAddedEdge++;
2354              = localTriangleListother.begin();                }
2355            currentTriangle != localTriangleListother.end();              }
2356            ++currentTriangle) {            }
2357          //std::cout<<"*************************\n";  
2358          const Edge::Pair& eo0=(*currentTriangle)[0];            switch(sens) {
2359          const Edge::Pair& eo1=(*currentTriangle)[1];            case 0: {
2360          const Edge::Pair& eo2=(*currentTriangle)[2];              T[0] = TetrahedronByEdges(addedEdges[0],
2361          const Vertex& Vo0 = *((*verticesListes[1])[eo0]);                                        addedEdges[1],
2362          const Vertex& Vo1 = *((*verticesListes[1])[eo1]);                                        addedEdges[2],
2363          const Vertex& Vo2 = *((*verticesListes[1])[eo2]);                                        addedEdges[3],
2364                                          addedEdges[4],
2365          __CalculPoint(V0, V1, V2, Vo0, Vo1, Vo2, nn,                                        addedEdges[5]);
2366                        triangleedge[e0], triangleedge[e1], triangleedge[e2],  
2367                        listofVertex[e0], listofVertex[e1], listofVertex[e2]);              T[1] = TetrahedronByEdges(H.edge(0),
2368        }                                        H.edge(4),
2369        ++nn;                                        H.edge(8),
2370        //ffout(0) << "nn "<<nn<<"\n";                                        addedEdges[0],
2371      }                                        addedEdges[1],
2372      localTriangleList.clear();                                        addedEdges[3]);
2373      localTriangleListother.clear();  
2374    }              T[2] = TetrahedronByEdges(H.edge(1),
2375  /*    nn=0;                                        H.edge(5),
2376                                          H.edge(11),
2377      for(SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle1                                        addedEdges[0],
2378          = (*triangleListes[0]).begin();                                        addedEdges[2],
2379          currentTriangle1 != (*triangleListes[0]).end();                                        addedEdges[4]);
2380          ++currentTriangle1)  
2381      {              T[3] = TetrahedronByEdges(H.edge(2),
2382          const Edge::Pair& e0=(*currentTriangle1)[0];                                        H.edge(7),
2383          const Edge::Pair& e1=(*currentTriangle1)[1];                                        H.edge(9),
2384          const Edge::Pair& e2=(*currentTriangle1)[2];                                        addedEdges[1],
2385          const Vertex& V0 = *((*verticesListes[0])[e0]);                                        addedEdges[2],
2386          const Vertex& V1 = *((*verticesListes[0])[e1]);                                        addedEdges[5]);
2387          const Vertex& V2 = *((*verticesListes[0])[e2]);  
2388                T[4] = TetrahedronByEdges(H.edge(3),
2389          //TinyVector<3>  V01,V12,V02;                                        H.edge(6),
2390          //V01 = V1.Vector3() - V0.Vector3();                                        H.edge(10),
2391          //V12 = V2.Vector3() - V1.Vector3();                                        addedEdges[3],
2392          //V02 = V2.Vector3() - V0.Vector3();                                        addedEdges[4],
2393          //  const TinyVector<3> coeff=V01^V02;                                        addedEdges[5]);
2394   //  const real_t d=-(coeff[0]*V0[0]+coeff[1]*V0[1]+coeff[2]*V0[2]);              break;
2395   //      std::cout<<"--------------------------------------------------------------\n";            }
2396              case 1: {
2397          tt[nn] = 0;              T[0] = TetrahedronByEdges(addedEdges[0],
2398                                          addedEdges[1],
2399          // Pour Cecile :                                        addedEdges[2],
2400   // Je crois que tu ne pouvais pas espérer que ça marche pour plus de 2 objets                                        addedEdges[3],
2401   // puisque tu écris explicitement d'utiliser la liste [1]                                        addedEdges[4],
2402                                          addedEdges[5]);
2403          //       for(SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle  
2404   //         = (*(*__internals).triangleListes[1]).begin();              T[1] = TetrahedronByEdges(H.edge(0),
2405   //       currentTriangle != (*(*__internals).triangleListes[1]).end();                                        H.edge(5),
2406   //       ++currentTriangle)                                        H.edge(9),
2407                                          addedEdges[0],
2408          for (size_t currentObjectNumber = 1; currentObjectNumber < scene.nbObjects(); ++currentObjectNumber) {                                        addedEdges[1],
2409                                          addedEdges[3]);
 SurfaceMeshGenerator::Internals::EdgePairList& otherTrianglesList  
             = *triangleListes[currentObjectNumber];  
   
 SurfaceMeshGenerator::Internals::VerticesList&  
             otherVerticesList = (*verticesListes[currentObjectNumber]);  
   
             for(SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle  
          = otherTrianglesList.begin();  
          currentTriangle != otherTrianglesList.end();  
          ++currentTriangle)  
             {  
                 //std::cout<<"*************************\n";  
                 const Vertex& Vo0 = *(otherVerticesList[(*currentTriangle)[0]]);  
                 const Vertex& Vo1 = *(otherVerticesList[(*currentTriangle)[1]]);  
                 const Vertex& Vo2 = *(otherVerticesList[(*currentTriangle)[2]]);  
2410    
2411                  __CalculPoint(V0,V1,V2,Vo0,Vo1,Vo2,nn,tt[nn],vv[nn]);              T[2] = TetrahedronByEdges(H.edge(1),
2412                                          H.edge(4),
2413                                          H.edge(10),
2414                                          addedEdges[0],
2415                                          addedEdges[2],
2416                                          addedEdges[4]);
2417    
2418                T[3] = TetrahedronByEdges(H.edge(2),
2419                                          H.edge(6),
2420                                          H.edge(8),
2421                                          addedEdges[1],
2422                                          addedEdges[2],
2423                                          addedEdges[5]);
2424    
2425                T[4] = TetrahedronByEdges(H.edge(3),
2426                                          H.edge(7),
2427                                          H.edge(11),
2428                                          addedEdges[3],
2429                                          addedEdges[4],
2430                                          addedEdges[5]);
2431    
2432              }
2433              }
2434    
2435              bool jobDone = true;
2436              for (std::vector<TetrahedronByEdges>::iterator currentT = T.begin();
2437                   currentT < T.end(); ++currentT) {
2438                std::vector<const Edge*> cutEdges;
2439    
2440                for(size_t ie=0; ie<6; ++ie) {
2441                  const Edge& e = currentT->edge(ie);
2442                  if ((*edgesRef[nobject]).find(e) != (*edgesRef[nobject]).end())
2443                    cutEdges.push_back(&e);
2444              }              }
2445    
2446                bool jobDone2=__createTriangleSurface(nobject,S,Omega,cutEdges,
2447                                         ref,currentCell);
2448                if(!(jobDone2))
2449                    jobDone=jobDone2;
2450    
2451              }
2452                
2453              if(jobDone) {
2454                //  if(nobject==1)
2455                      meshedHexahedra.insert(&currentCell);
2456                  for (SurfaceMeshGenerator::Internals::EdgePairList::iterator t
2457                = localTriangleList.begin();
2458                t != localTriangleList.end(); ++t) {
2459                      (*triangleListes[nobject]).push_front(*t);
2460                  }
2461              } else {
2462                 // if(nobject==1)
2463                      hexaTetraAssociation[&currentCell]
2464                          = std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,
2465                          ReferenceCounting<std::vector<Edge> > >
2466                          (pT,pAddedEdges);
2467              }
2468              localTriangleList.clear();
2469          }          }
2470          nn++;        }
     }*/  
 }  
2471    
2472    }
2473    
2474  void  void
2475  SurfaceMeshGenerator::  SurfaceMeshGenerator::
# Line 1048  generatesAllSufaceMeshes(const Problem& Line 2486  generatesAllSufaceMeshes(const Problem&
2486    TinyVector<6,bool> faces = false;    TinyVector<6,bool> faces = false;
2487    
2488    for (size_t i=0; i<problem.numberOfUnknown(); ++i) {    for (size_t i=0; i<problem.numberOfUnknown(); ++i) {
2489        const BoundaryConditionSet& bcSet      const BoundaryConditionSet& bcSet
2490          = problem.boundaryConditionSet(i);        = problem.boundaryConditionSet(i);
2491    
2492      for (size_t j=0; j<bcSet.nbBoundaryCondition(); ++j) {      for (size_t j=0; j<bcSet.nbBoundaryCondition(); ++j) {
2493        const Boundary& B = bcSet[j].boundary();        const Boundary& B = bcSet[j].boundary();
# Line 1248  void plot(const SurfaceMeshOfTriangles& Line 2686  void plot(const SurfaceMeshOfTriangles&
2686            std::set<const Cell*>& hlist,            std::set<const Cell*>& hlist,
2687            Structured3DMesh& SMesh);            Structured3DMesh& SMesh);
2688  void plot2(const SurfaceMeshOfTriangles& s,  void plot2(const SurfaceMeshOfTriangles& s,
2689            std::set<const Cell*>& hlist,             std::set<const Cell*>& hlist,
2690            Structured3DMesh& SMesh);             Structured3DMesh& SMesh);
2691    
2692    
2693  void plotmedit(const SurfaceMeshOfTriangles& s_mesh,  void plotmedit(const size_t& nobject,
2694                   const SurfaceMeshOfTriangles& s_mesh,
2695                 std::set<const Cell*>& toTreatHexahedra,                 std::set<const Cell*>& toTreatHexahedra,
2696                 size_t nobj);                 std::set<const Cell*>& cuttedHexahedra,
2697                   std::vector<ReferenceCounting<std::map<Edge::Pair, Vertex*> > > verticesListesIntersection,
2698                   std::vector<ReferenceCounting<std::list<TinyVector<3,TinyVector<2,Edge::Pair> > > > >
2699                   triangleListesIntersectionEdges,
2700                   size_t nobjother,size_t ninter,size_t ninter2,size_t ntr1);
2701    
2702    void plothexa(size_t ncase,std::set<const Cell*>& cuttedHexahedra);
2703    
   /*! Generates a surfacic mesh using only characteristic function of an object  
     and the Structured 3D Mesh.  
2704    
2705      \warning Implementation could be improved and function size reduced.  /*! Generates a surfacic mesh using only characteristic function of an object
2706    */    and the Structured 3D Mesh.
2707  void  
2708  SurfaceMeshGenerator::    \warning Implementation could be improved and function size reduced.
2709    */
2710    void SurfaceMeshGenerator::
2711  generateSurfacicMesh(const Object& O,  generateSurfacicMesh(const Object& O,
2712                       const Domain& Omega,                       const Domain& Omega,
2713                       SurfaceMeshOfTriangles& s_mesh,                       SurfaceMeshOfTriangles& s_mesh,
2714                       Structured3DMesh& SMesh,                       Structured3DMesh& SMesh,
2715                       const Reference::Representation& ref,                       const Reference::Representation& ref,
2716                       const std::set<Reference::Representation>& referencesList)                       const std::set<Reference::Representation>& referencesList) {
 {  
   const Shape& S = *(O.shape());  
2717    
2718      const Shape& S = *(O.shape());
2719        
2720    typedef std::map<TinyVector<3>, Object*,    typedef std::map<TinyVector<3>, Object*,
2721                     SurfaceMeshGenerator::Internals::__ltv> ObjectReferences;    SurfaceMeshGenerator::Internals::__ltv> ObjectReferences;
2722    
2723    ObjectReferences otherReferences;    ObjectReferences otherReferences;
2724    
# Line 1290  generateSurfacicMesh(const Object& O, Line 2735  generateSurfacicMesh(const Object& O,
2735      break;      break;
2736    }  case Shape::intersection: {    }  case Shape::intersection: {
2737      const Intersection& I = static_cast<const Intersection&>(*(object.shape()));      const Intersection& I = static_cast<const Intersection&>(*(object.shape()));
2738      for(Intersection::const_iterator i = I.begin();      Internals::__getIntersectionReferences(I, otherReferences);
         i != I.end(); ++i) {  
         otherReferences[(*(*i)).reference()] = *i;  
         //ffout(0) << (*(*i)).reference() << "\n";  
     }  
2739      break;      break;
2740      } case Shape::difference: {
2741          fferr(0) << __FILE__ << ':' << __LINE__ << ": Not implemented\n";
2742          break;
2743      } case Shape::not_: {
2744        //  fferr(0) << __FILE__ << ':' << __LINE__ << ": Not implemented\n";
2745          break;
2746      } default: {  
2747      //std::exit(1);
2748    }    }
2749    default: {    fferr(0) << __FILE__ << ':' << __LINE__ << ": Not implemented\n";    }    
     std::exit(1);  
   }  
   }  
   
   otherReferences.erase(O.reference());  
2750    
2751  //  for(ObjectReferences::iterator i=otherReferences.begin();i!=otherReferences.end();++i)    //  for(ObjectReferences::iterator i=otherReferences.begin();i!=otherReferences.end();++i)
2752    //    ffout(0) << (*i).first <<"\n";    //    ffout(0) << (*i).first <<"\n";
2753      
2754    size_t nobject = scene.nbObjects();    size_t nobject = scene.nbObjects();
2755      size_t nobjectToTreat = otherReferences.size();
2756      ffout(9)<<"size other ref "<<nobjectToTreat<<"\n";
2757    
2758    (*__internals).verticesListes.resize(nobject);    (*__internals).verticesListes.resize(nobject);
2759    (*__internals).triangleListes.resize(nobject);    (*__internals).triangleListes.resize(nobject);
2760      (*__internals).triangleListesIntersectionEdges.resize(nobject);
2761      (*__internals).triangleListesIntersectionEdgesTemp.resize(nobject);
2762      (*__internals).localTriangleListesEdges.resize(nobject);
2763      (*__internals).localTriangleList2Edges.resize(nobject);
2764      (*__internals).verticesListesIntersection.resize(nobject);
2765      (*__internals).edgesRef.resize(nobject);
2766      (*__internals).edgesRefCreate.resize(nobject);
2767      (*__internals).listTriangleintersection.resize(nobject);
2768    
2769    std::cout<<"nb object "<<nobject<<"\n";      
2770      for(size_t no=0 ; no<nobject ; ++no) {
   for(size_t no=0 ; no<nobject ; ++no)  
     {  
2771        (*__internals).triangleListes[no]= new SurfaceMeshGenerator::Internals::EdgePairList;        (*__internals).triangleListes[no]= new SurfaceMeshGenerator::Internals::EdgePairList;
2772          (*__internals).triangleListesIntersectionEdges[no]=
2773              new SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList;
2774          (*__internals).triangleListesIntersectionEdgesTemp[no]=
2775              new std::map<SurfaceMeshGenerator::Internals::TriangleVertices*,
2776                      const SurfaceMeshGenerator::Internals::TriangleVertices*>;
2777    
2778          (*__internals).localTriangleListesEdges[no]= new SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList;
2779          (*__internals).localTriangleList2Edges[no]=
2780              new SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList2;
2781        (*__internals).verticesListes[no]= new SurfaceMeshGenerator::Internals::VerticesList;        (*__internals).verticesListes[no]= new SurfaceMeshGenerator::Internals::VerticesList;
2782      }        (*__internals).verticesListesIntersection[no]= new SurfaceMeshGenerator::Internals::VerticesList;
2783    
2784    for (size_t i = 0; i<SMesh.nbEdges(); i++) {        (*__internals).edgesRef[no]= new std::map<Edge::Pair, size_t> ;
2785      Edge& e = SMesh.edge(i);        (*__internals).edgesRefCreate[no]= new std::map<Edge::Pair, size_t> ;
2786  #warning THIS TEST IS NOT ENOUGH! THINK TO THE BORDER!        (*__internals).listTriangleintersection[no]=
2787      if((e(0).Ref().insideDomain() xor e(1).Ref().insideDomain())            new SurfaceMeshGenerator::Internals::TriangleVerticesEdges;
        and((e(0).Ref().isSet(ref)) xor (e(1).Ref().isSet(ref))))  
       {  
         (*(*__internals).verticesListes[0])[e] = new Vertex((*__internals).__splitEdge(e, ref, S));  
         //(*__internals).verticesList[e] = new Vertex((*__internals).__splitEdge(e, ref, S));  
         (*__internals).edgesRef[e] = 1;  
       }  
2788    }    }
2789    
2790    /// at least one edge of the hexa has reference 1    
2791      /// at least one edge of the hexa has reference 1
2792    std::set<const Cell*> cuttedHexahedra;    std::set<const Cell*> cuttedHexahedra;
2793    
2794    /// a triangulation has been generated inside the hexahedron    /// a triangulation has been generated inside the hexahedron
2795    std::set<const Cell*> meshedHexahedra;    std::set<const Cell*> meshedHexahedra;
2796    
2797    if ((*(*__internals).verticesListes[0]).size()==0) {    int num=0;
2798      if(nobjectToTreat==0) {
2799      fferr(2) << __FILE__ << ':' << __LINE__ << ": "        (*__internals).__createSurface(SMesh,Omega,S,ref,num,
2800               << "meshing object " << S << " a problem occured ...\n"                                       cuttedHexahedra,meshedHexahedra);
2801               << "Check that:\n"    }
2802               << "- the object is not to small compare to the background mesh\n"    else {
2803               << "- the object intersectes the background mesh\n\n";        for(ObjectReferences::iterator i = otherReferences.begin();
2804      std::exit(1);            i != otherReferences.end(); ++i) {
2805    }            for(size_t n = 0; n<scene.nbObjects(); ++n) {
2806                  Reference::Representation otherRef;
2807    for (size_t icell=0; icell<SMesh.shape().nx()-1; ++icell)                otherRef.resize(scene.nbObjects());
2808      for (size_t jcell=0; jcell<SMesh.shape().ny()-1; ++jcell)                //ffout(0) << scene.object(n).reference() << ":"<< ((*i).first) << '\n';
2809        for (size_t kcell=0; kcell<SMesh.shape().nz()-1; ++kcell) {                if (scene.object(n).reference() == ((*i).first)) {
2810          bool intersected = false;                    otherRef[n] = true;
2811                      //create surface objects
2812          const Cell& currentCell = SMesh.cell(icell, jcell, kcell);                    (*__internals).__createSurface(SMesh,Omega,*(scene.object(n).shape()),otherRef,num,
2813          const size_t nx                                     cuttedHexahedra,meshedHexahedra);
           = dynamic_cast<const Structured3DMesh&>(SMesh).shape().nx();  
         const size_t nx_1 = nx-1;  
   
         const size_t ny  
           = dynamic_cast<const Structured3DMesh&>(SMesh).shape().ny();  
         const size_t ny_1 = ny-1;  
   
         const size_t nz  
           = dynamic_cast<const Structured3DMesh&>(SMesh).shape().nz();  
         const size_t nz_1 = nz-1;  
   
         HexahedronByEdges H (SMesh.edge(icell*ny*nz + jcell*nz + kcell),  
                              SMesh.edge(icell*ny*nz + (jcell+1)*nz + kcell),  
                              SMesh.edge(icell*ny*nz + jcell*nz + kcell + 1),  
                              SMesh.edge(icell*ny*nz + (jcell+1)*nz + kcell+1),  
                              SMesh.edge(nx_1*ny*nz + icell*ny_1*nz + jcell*nz + kcell),  
                              SMesh.edge(nx_1*ny*nz + (icell+1)*ny_1*nz + jcell*nz + kcell),  
                              SMesh.edge(nx_1*ny*nz + icell*ny_1*nz + jcell*nz + kcell + 1),  
                              SMesh.edge(nx_1*ny*nz + (icell+1)*ny_1*nz + jcell*nz + kcell+1),  
                              SMesh.edge(nx_1*ny*nz + nx*ny_1*nz + icell*ny*nz_1 + jcell*nz_1 + kcell ),  
                              SMesh.edge(nx_1*ny*nz + nx*ny_1*nz + (icell+1)*ny*nz_1 + jcell*nz_1 + kcell ),  
                              SMesh.edge(nx_1*ny*nz + nx*ny_1*nz + icell*ny*nz_1 + (jcell+1)*nz_1 + kcell ),  
                              SMesh.edge(nx_1*ny*nz + nx*ny_1*nz + (icell+1)*ny*nz_1 + (jcell+1)*nz_1 + kcell ));  
   
   
         for (size_t l=0; l<H.numberOfEdges(); l++) {  
           Edge& e = H.edge(l);  
           if ((*(*__internals).verticesListes[0]).count(e) > 0) {  
             intersected = true;  
             cuttedHexahedra.insert(&currentCell);  
             break;  
           }  
         }  
   
         if (intersected) { // Means that the hexahedron intersectes the object.  
           int sens = (icell+jcell+kcell)%2;  
           ReferenceCounting<std::vector<Edge> > pAddedEdges;  
           pAddedEdges = new std::vector<Edge>(6);  
           std::vector<Edge>& addedEdges = *pAddedEdges;  
   
           std::vector<int> TverticesList(4);  
           switch (sens) {  
           case 0: {  
             TverticesList[0] = 1;  
             TverticesList[1] = 3;  
             TverticesList[2] = 4;  
             TverticesList[3] = 6;  
             break;  
           }  
           case 1: {  
             TverticesList[0] = 0;  
             TverticesList[1] = 2;  
             TverticesList[2] = 5;  
             TverticesList[3] = 7;  
           }  
           }  
   
           ReferenceCounting<std::vector<TetrahedronByEdges> > pT;  
           pT = new std::vector<TetrahedronByEdges>(5);  
   
           std::vector<TetrahedronByEdges>& T = *pT;  
   
           { // Create tetrahedra edge  
             std::vector<Edge>::iterator currentAddedEdge  
               = addedEdges.begin();  
             for (size_t i1=0; i1<3; ++i1) {  
               for (size_t i2=i1; i2<3; ++i2) {  
                 *(currentAddedEdge) = Edge(currentCell(TverticesList[i1]),  
                                            currentCell(TverticesList[i2+1]));  
                 Edge& e = (*currentAddedEdge);  
                 if (((e(0).Ref().isSet(ref)) xor (e(1).Ref().isSet(ref)))  
                     and (e(0).Ref().insideDomain() xor e(1).Ref().insideDomain())) {  
                   e.Ref() = 1;  
                   (*__internals).edgesRef[e] = 1;  
                 }  
                 currentAddedEdge++;  
2814                }                }
             }  
2815            }            }
2816              plothexa(num,meshedHexahedra);
2817              if(num==0)
2818                  plothexa(num-1,cuttedHexahedra);
2819              num++;
2820          }
2821      }
2822    
2823            switch(sens) {    std::set_difference (cuttedHexahedra.begin(), cuttedHexahedra.end(),
2824            case 0: {                       meshedHexahedra.begin(), meshedHexahedra.end(),
2825              T[0] = TetrahedronByEdges(addedEdges[0],                       inserter((*__internals).toTreatHexahedra, (*__internals).toTreatHexahedra.begin()));
                                       addedEdges[1],  
                                       addedEdges[2],  
                                       addedEdges[3],  
                                       addedEdges[4],  
                                       addedEdges[5]);  
   
             T[1] = TetrahedronByEdges(H.edge(0),  
                                       H.edge(4),  
                                       H.edge(8),  
                                       addedEdges[0],  
                                       addedEdges[1],  
                                       addedEdges[3]);  
   
             T[2] = TetrahedronByEdges(H.edge(1),  
                                       H.edge(5),  
                                       H.edge(11),  
                                       addedEdges[0],  
                                       addedEdges[2],  
                                       addedEdges[4]);  
   
             T[3] = TetrahedronByEdges(H.edge(2),  
                                       H.edge(7),  
                                       H.edge(9),  
                                       addedEdges[1],  
                                       addedEdges[2],  
                                       addedEdges[5]);  
   
             T[4] = TetrahedronByEdges(H.edge(3),  
                                       H.edge(6),  
                                       H.edge(10),  
                                       addedEdges[3],  
                                       addedEdges[4],  
                                       addedEdges[5]);  
             break;  
           }  
           case 1: {  
             T[0] = TetrahedronByEdges(addedEdges[0],  
                                       addedEdges[1],  
                                       addedEdges[2],  
                                       addedEdges[3],  
                                       addedEdges[4],  
                                       addedEdges[5]);  
   
             T[1] = TetrahedronByEdges(H.edge(0),  
                                       H.edge(5),  
                                       H.edge(9),  
                                       addedEdges[0],  
                                       addedEdges[1],  
                                       addedEdges[3]);  
   
             T[2] = TetrahedronByEdges(H.edge(1),  
                                       H.edge(4),  
                                       H.edge(10),  
                                       addedEdges[0],  
                                       addedEdges[2],  
                                       addedEdges[4]);  
   
             T[3] = TetrahedronByEdges(H.edge(2),  
                                       H.edge(6),  
                                       H.edge(8),  
                                       addedEdges[1],  
                                       addedEdges[2],  
                                       addedEdges[5]);  
   
             T[4] = TetrahedronByEdges(H.edge(3),  
                                       H.edge(7),  
                                       H.edge(11),  
                                       addedEdges[3],  
                                       addedEdges[4],  
                                       addedEdges[5]);  
2826    
2827      std::set<const Cell*>::iterator iterat=(*__internals).toTreatHexahedra.begin();
2828      size_t treatsize=(*__internals).toTreatHexahedra.size();
2829      for(size_t yy=0 ; yy<treatsize ; ++yy)
2830      {
2831          SurfaceMeshGenerator::Internals::HexaTetraAssociation::iterator
2832          j=(*__internals).hexaTetraAssociation.find(*iterat);
2833          std::vector<TetrahedronByEdges>& T = *(((*j).second).first);
2834          const CartesianHexahedron& currentCell
2835              = static_cast<const CartesianHexahedron&>(*(*iterat));
2836         // ffout(0)<<"*********************************\n";
2837          size_t test=0;
2838          std::set<const Cell*>::iterator temp=iterat++;
2839          for (std::vector<TetrahedronByEdges>::iterator currentT = T.begin();
2840               currentT < T.end(); ++currentT) {
2841    
2842              std::vector<ReferenceCounting< std::vector<const Edge*> > > cutEdges;
2843              cutEdges.resize(nobject);
2844              for(size_t no=0 ; no<nobject ; ++no) {
2845                  cutEdges[no]= new std::vector<const Edge*>;
2846            }            }
2847            }            (*__internals).__createCutEdges(scene,cutEdges,otherReferences,currentT);
2848              for(size_t no=0 ; no<nobject ; ++no) {
2849            bool jobDone = true;                if((*cutEdges[no]).size()>0) {
2850                      test++;
           for (std::vector<TetrahedronByEdges>::iterator currentT = T.begin();  
                currentT < T.end(); ++currentT) {  
             std::vector<const Edge*> cutEdges;  
   
             for(size_t ie=0; ie<6; ++ie) {  
               const Edge& e = currentT->edge(ie);  
               if ((*__internals).edgesRef.find(e) != (*__internals).edgesRef.end())  
                 cutEdges.push_back(&e);  
             }  
   
             switch (cutEdges.size()) {  
             case 3: {  
               //! Addes a triangle to the list.  
               (*__internals).localTriangleList.push_front(TinyVector<3,Edge::Pair>());  
               (*__internals).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());  
               }  
   
               if ((V[0]^V[1])*V[2] < 0)  
                 std::swap(cutEdges[1], cutEdges[2]);  
   
               //! Sets the triangle vertices (actually just store the edge where it lives).  
               SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle  
                 = (*__internals).localTriangleList.begin();  
   
               for(size_t i=0; i<cutEdges.size(); ++i)  
                 (*currentTriangle)[i] = (Edge::Pair)(*cutEdges[i]);  
   
               //! Stores the Edges on which the vertex will be computed.  
               for (size_t i=0; i<cutEdges.size(); ++i) {  
                 //(*__internals).verticesList[*cutEdges[i]]  
                   (*(*__internals).verticesListes[0])[*cutEdges[i]]  
                   =  new Vertex((*__internals).__splitEdge(*cutEdges[i], ref, S));  
                   (*__internals).edgesRef[*cutEdges[i]] = 1;  
               }  
               break;  
             }  
             case 4: {  
               //! 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.  
               (*__internals).localTriangleList.push_front(TinyVector<3,Edge::Pair>());  
               (*__internals).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;  
2851                }                }
   
               bool swapEdge = ((V[0]^V[1])*V[2] > 0);  
   
               //! Sets the triangle vertices (actually just store the edge where it lives).  
               SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle  
                 = (*__internals).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;  
                 }  
   
               for (size_t i=0; i<4; ++i) {  
                 //(*__internals).verticesList[*cutEdges[i]]  
                   (*(*__internals).verticesListes[0])[*cutEdges[i]]  
                   =  new Vertex((*__internals).__splitEdge(*cutEdges[i], ref, S));  
                         (*__internals).edgesRef[*cutEdges[i]] = 1;  
               }  
   
               //! Addes a triangle to the list.  
               (*__internals).localTriangleList.push_front(TinyVector<3,Edge::Pair>());  
               (*__internals).cellList.push_front(&currentCell);  
   
               currentTriangle = (*__internals).localTriangleList.begin();  
   
               //! Sets the triangle vertices (actually just store the edge where it lives).  
               currentCutEdge = cutEdgesMap.begin();  
               if (!swapEdge)  
                 for(int i=2; i>=0; --i) {  
                   ++currentCutEdge;  
                   (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);  
                 }  
               else  
                 for(size_t i=0; i<3; ++i) {  
                   ++currentCutEdge;  
                   (*currentTriangle)[i] = (Edge::Pair)(currentCutEdge->first);  
                 }  
               break;  
             }  
             case 1:  
             case 2: {  
                 jobDone = false;  
               break;  
             }  
             }  
           }  
   
           if(jobDone) {  
             meshedHexahedra.insert(&currentCell);  
             for (SurfaceMeshGenerator::Internals::EdgePairList::iterator t  
                    = (*__internals).localTriangleList.begin();  
                  t != (*__internals).localTriangleList.end(); ++t) {  
                 (*(*__internals).triangleListes[0]).push_front(*t);  
   
             }  
 //          (*__internals).hexaTetraAssociation[&currentCell]  
 //              = std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,  
 //              ReferenceCounting<std::vector<Edge> > >  
 //              (pT,pAddedEdges);  
               
           } else {  
             (*__internals).hexaTetraAssociation[&currentCell]  
               = std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,  
                           ReferenceCounting<std::vector<Edge> > >  
               (pT,pAddedEdges);  
               
2852            }            }
2853              if(test<2)
2854            (*__internals).localTriangleList.clear();                test=0;
2855          }            else
2856                  test=2;
2857        }        }
2858          if(test!=2) {
2859  /*  (*(*__internals).triangleListes[0]).clear();            (*__internals).toTreatHexahedra.erase(temp);
   (*(*__internals).verticesListes[0]).clear();  
   
   (*__internals).__createTriangles(S,ref,0);  
   int num=0;  
   
   size_t nobj= (*(*__internals).triangleListes[0]).size();  
   std::cout<<"nb tr objet1 "<<nobj<<" \n";  
   
   for(ObjectReferences::iterator i = otherReferences.begin();  
       i != otherReferences.end(); ++i) {  
     num++;  
     std::cout<<"num "<<num<<"\n";  
     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;  
         //std::cout<<"here\n";  
         (*__internals).__createTriangles(*(scene.object(n).shape()),  
                                          otherRef,num);  
       }  
     }  
   }*/  
   
   size_t nobjother=(*(*__internals).triangleListes[1]).size();  
   std::cout<<"nb tr objet2 "<<nobjother<<" \n";  
   if(scene.nbObjects()>2)  
     std::cout<<"nb tr objet3 "<<(*(*__internals).triangleListes[2]).size()<<" \n";  
 /*  
   // *******************************************************************************  
   
 SurfaceMeshGenerator::Internals::VerticesList listofVertex;  
 //  std::map<Edge::Pair, Vertex*>  listofVertex;  
   std::map<Edge::Pair,int> triangleedge;  
   
 std::vector< TinyVector<3,int> > tt;  
   std::vector< std::vector<Vertex> > vv;  
   vv.resize((*(*__internals).triangleListes[0]).size());  
   tt.resize((*(*__internals).triangleListes[0]).size());  
   for(size_t kk=0 ;kk<nobj;kk++)  
     vv[kk].resize(3);  
   
   num=0;  
   for(ObjectReferences::iterator i = otherReferences.begin();  
       i != otherReferences.end(); ++i) {  
       num++;  
       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;  
               (*__internals).__IntersectionPoint(scene,ref,otherRef,triangleedge,listofVertex);  
           }  
2860        }        }
2861    }    }
2862      ffout(0)<<"nb tr objet 0 complet "<<(*(*__internals).triangleListes[0]).size()<<" \n";
2863        
2864    // create triangles    ffout(9)<<"calculate points intersection\n";
2865    std::vector<Vertex> vertexplot;    (*__internals).__calculatePointsIntersection(scene,otherReferences);
2866    {    
2867      int nn=0;    ffout(9)<<"liste points "<<(*__internals).listVerticesIntersectionPerEdges.size()<<"\n";
2868      size_t ii0,ii1,ii2;    ffout(9)<<"nb tr intersection objet 0 "<<(*(*__internals).triangleListesIntersectionEdges[0]).size()<<" \n";
     for(SurfaceMeshGenerator::Internals::EdgePairList::const_iterator currentTriangle  
           = (*(*__internals).triangleListes[0]).begin();  
         currentTriangle != (*(*__internals).triangleListes[0]).end();  
         ++currentTriangle) {  
       const Edge& e0=(*currentTriangle)[0];  
       const Edge& e1=(*currentTriangle)[1];  
       const Edge& e2=(*currentTriangle)[2];  
         
      // if(tn0+tn1+tn2!=0) {  
 //std::cout<<"-------------------------------------------------"<<"\n";  
 //std::cout<<"tt "<<nn<<" "<<tn0<<" "<<tn1<<" "<<tn2<<"\n";  
 //std::cout<<"avt "<<nn<<" "<<(*__internals).edgesRef[e0]<<" "<<(*__internals).edgesRef[e1]<<" "<<(*__internals).edgesRef[e2]<<"\n";  
       //}  
       (*__internals).__inout(1,nn,triangleedge[e0],triangleedge[e1],triangleedge[e2],  
                              (*__internals).edgesRef[e0],(*__internals).edgesRef[e1]  
                              ,(*__internals).edgesRef[e2],ii0,ii1,ii2);  
       //if(tn0+tn1+tn2!=0) {  
 //std::cout<<"apres "<<nn<<" "<<(*__internals).edgesRef[e0]<<" "<<(*__internals).edgesRef[e1]<<" //"<<(*__internals).edgesRef[e2]<<"\n";  
       //}  
       nn++;  
   
       
     }  
 ffout(0) << "fin premiere boucle\n";  
     nn=0;  
     for(SurfaceMeshGenerator::Internals::EdgePairList::const_iterator currentTriangle  
           = (*(*__internals).triangleListes[0]).begin();  
         currentTriangle != (*(*__internals).triangleListes[0]).end();  
         ++currentTriangle) {  
       const Edge& e0=(*currentTriangle)[0];  
       const Edge& e1=(*currentTriangle)[1];  
       const Edge& e2=(*currentTriangle)[2];  
       if(tt[nn][0]+tt[nn][1]+tt[nn][2]!=0) {  
 //std::cout<<"-------------------------------------------------"<<"\n";  
 //std::cout<<"tt "<<nn<<" "<<tt[nn][0]<<" "<<tt[nn][1]<<" "<<tt[nn][2]<<"\n";  
 //std::cout<<"avt "<<nn<<" "<<(*__internals).edgesRef[e0]<<" "<<(*__internals).edgesRef[e1]<<" "<<(*__internals).edgesRef[e2]<<"\n";  
       }  
       (*__internals).__inout(2,nn,triangleedge[e0],triangleedge[e1],triangleedge[e2],  
                              (*__internals).edgesRef[e0],(*__internals).edgesRef[e1]  
                              ,(*__internals).edgesRef[e2],ii0,ii1,ii2);  
       if(tt[nn][0]+tt[nn][1]+tt[nn][2]!=0) {  
 //std::cout<<"apres "<<nn<<" "<<(*__internals).edgesRef[e0]<<" "<<(*__internals).edgesRef[e1]<<" //"<<(*__internals).edgesRef[e2]<<"\n";  
       }  
         
       nn++;  
     }  
 ffout(0) << "fin deuxieme boucle\n";  
2869    
2870      nn=0;    size_t nobjother=0;
2871      for(SurfaceMeshGenerator::Internals::EdgePairList::const_iterator currentTriangle    for(size_t no=0;no<nobjectToTreat;++no) {
2872            = (*(*__internals).triangleListes[0]).begin();        nobjother+=(*(*__internals).triangleListes[no]).size();
2873          currentTriangle != (*(*__internals).triangleListes[0]).end();    }
2874          ++currentTriangle) {  //  ffout(9)<<"nb tr objet1 complet "<<nobjother<<" \n";
2875        const Edge& e0=(*currentTriangle)[0];  //  ffout(9)<<"nb tr objet1 intersection "<<(*(*__internals).triangleListesIntersectionEdges[1]).size()<<" \n";
2876        const Edge& e1=(*currentTriangle)[1];  
2877        const Edge& e2=(*currentTriangle)[2];    if(scene.nbObjects()>2) {
2878        if(tt[nn][0]+tt[nn][1]+tt[nn][2]!=0) {        ffout(9)<<"nb tr objet2 intersection "<<(*(*__internals).triangleListesIntersectionEdges[2]).size()<<" \n";
2879  //std::cout<<"-------------------------------------------------"<<"\n";    }
2880  //std::cout<<"tt "<<nn<<" "<<tt[nn][0]<<" "<<tt[nn][1]<<" "<<tt[nn][2]<<"\n";  
2881  //std::cout<<"avt "<<nn<<" "<<(*__internals).edgesRef[e0]<<" "<<(*__internals).edgesRef[e1]<<" "<<(*__internals).edgesRef[e2]<<"\n";     for(size_t no=0;no<nobjectToTreat;++no) {
2882        }        for(SurfaceMeshGenerator::Internals::TriangleVerticesEdgesList::iterator currentTriangle1
2883        (*__internals).__inout(3,nn,triangleedge[e0],triangleedge[e1],triangleedge[e2],            = (*(*__internals).triangleListesIntersectionEdges[no]).begin();
2884                               (*__internals).edgesRef[e0],(*__internals).edgesRef[e1]            currentTriangle1 != (*(*__internals).triangleListesIntersectionEdges[no]).end();
2885                               ,(*__internals).edgesRef[e2],ii0,ii1,ii2);            ++currentTriangle1) {
2886        if(tt[nn][0]+tt[nn][1]+tt[nn][2]!=0) {            (*(*__internals).triangleListesIntersectionEdgesTemp[no])[&(*currentTriangle1)]=&(*currentTriangle1);
 //std::cout<<"apres "<<nn<<" "<<(*__internals).edgesRef[e0]<<" "<<(*__internals).edgesRef[e1]<<" //"<<(*__internals).edgesRef[e2]<<"\n";  
2887        }        }
2888            }
       nn++;  
 }  
 ffout(0) << "fin troisieme boucle\n";  
2889    
2890      nn=0;    //create triangle with intersection points
2891      for(SurfaceMeshGenerator::Internals::EdgePairList::const_iterator currentTriangle   // nobjectToTreat=2;
2892            = (*(*__internals).triangleListes[0]).begin();   for(size_t no=0 ; no<nobjectToTreat ; ++no) {
2893          currentTriangle != (*(*__internals).triangleListes[0]).end();        for(size_t no2=no+1 ; no2<nobjectToTreat ; ++no2) {
2894          ++currentTriangle) {            (*__internals).__createTrianglesIntersection(no,no2,nobjectToTreat);
2895        //if the triangle is intersected            (*__internals).__createTrianglesIntersection(no2,no,nobjectToTreat);
2896        //std::cout<<"**************************************************"<<"\n";            for(size_t r=0 ; r<2 ; ++r) {
2897        const Edge& e0=(*currentTriangle)[0];                if(r==0)
2898        const Edge& e1=(*currentTriangle)[1];                    num=no;
2899        const Edge& e2=(*currentTriangle)[2];                else
2900        int tn0=triangleedge[e0];                    num=no2;
2901        int tn1=triangleedge[e1];                ffout(9)<<"triangleErase size "<<(*__internals).triangleEraseEdgesTemp[r].size()<<"\n";
2902        int tn2=triangleedge[e2];                for(size_t k=0 ; k<(*__internals).triangleEraseEdgesTemp[r].size() ; ++k) {
2903        //std::cout<<"  avt i123 "<<nn<<" "<<ii0<<" "<<ii1<<" "<<ii2<<" \n";                    //save adresse
2904        if(tn0+tn1+tn2!=0) {                    const SurfaceMeshGenerator::Internals::TriangleVertices* constantAdresse=
2905          std::cout<<"-------------------------------------------------"<<"\n";                    (*(*__internals).triangleEraseEdgesTemp[r][k]).second;
2906          std::cout<<"tt "<<nn<<" "<<tn0<<" "<<tn1<<" "<<tn2<<"\n";                    //erase triangle
2907            std::cout<<"avt "<<nn<<" "<<(*__internals).edgesRef[e0]<<" "<<(*__internals).edgesRef[e1]<<" "<<(*__internals).edgesRef[e2]<<"\n";                    (*(*__internals).triangleListesIntersectionEdgesTemp[num]).
2908                          erase((*__internals).triangleEraseEdgesTemp[r][k]);
2909                      //add new triangle
2910                       SurfaceMeshGenerator::Internals::TriangleTreatMap::iterator Ttmiterator=
2911                          ((*__internals).triangleTreat[r]).find(constantAdresse);
2912                      for(size_t hh=0 ; hh<((*Ttmiterator).second).size();++hh){
2913                          SurfaceMeshGenerator::Internals::TriangleVertices tv=((*Ttmiterator).second)[hh];
2914                          (*(*__internals).triangleListesIntersectionEdgesTemp[num])[&tv]
2915                          =constantAdresse;
2916                      }
2917                  }
2918                  (*__internals).triangleEraseEdgesTemp[r].clear();
2919                  (*__internals).triangleTreat[r].clear();
2920              }
2921        }        }
2922            }
2923    
2924        //std::cout<<"ref "<<(*__internals).edgesRef[e0]<<" "<<(*__internals).edgesRef[e1]<<" "<<(*__internals).edgesRef[e2]<<"\n";    //ajout des points de verticesListesIntersection qui sont dans le maillage final
2925        //std::cout<<"inter "<<nn<<" "<<(vv[0][0]).Vector3()[0]<<" "<<(vv[0][0]).Vector3()[1]<<" "<<(vv[0][0]).Vector3()[2]<<" \n";    ffout(9)<<"nb vertices object 0 "<<(*(*__internals).verticesListes[0]).size()<<"\n";
2926        // if(triangleedge[(*currentTriangle)[0]]+triangleedge[(*currentTriangle)[1]]+ triangleedge[(*currentTriangle)[2]]==2)    ffout(9)<<"nb vertices object 1 "<<(*(*__internals).verticesListes[1]).size()<<"\n";
2927        //std::cout<<"avt "<<(*__internals).edgesRef[e0]<<" "<<(*__internals).edgesRef[e1]<<" "<<(*__internals).edgesRef[e2]<<"\n";  
2928      for(size_t kk=0 ; kk<nobject ; ++kk) {
2929        //std::cout<<"-------------------------------------------------"<<"\n";  
2930          for(size_t n=0 ; n<(*(*__internals).listTriangleintersection[kk]).size(); ++n) {
2931        (*__internals).__inout(4,nn,tn0,tn1,tn2,            TinyVector<2,Edge::Pair> p0=((*(*__internals).listTriangleintersection[kk]))[n][0];
2932                               (*__internals).edgesRef[e0],(*__internals).edgesRef[e1]            TinyVector<2,Edge::Pair> p1=((*(*__internals).listTriangleintersection[kk]))[n][1];
2933                               ,(*__internals).edgesRef[e2],ii0,ii1,ii2);            TinyVector<2,Edge::Pair> p2=((*(*__internals).listTriangleintersection[kk]))[n][2];
2934        if(tn0+tn1+tn2!=0) {            if(p0[0]==p0[1]) {
2935                  if((*(*__internals).verticesListesIntersection[kk]).find(p0[1])!=
2936          std::cout<<"apres "<<nn<<" "<<(*__internals).edgesRef[e0]<<" "<<(*__internals).edgesRef[e1]<<" "<<(*__internals).edgesRef[e2]<<"\n";            (*(*__internals).verticesListesIntersection[kk]).end()) {
2937          std::cout<<"ii "<<ii0<<" "<<ii1<<" "<<ii2<<"\n";                    if((*(*__internals).verticesListes[kk]).find(p0[1])==(*(*__internals).verticesListes[kk]).end())
2938        }                    {
2939        if(tn0+tn1+tn2==2) {                        (*(*__internals).verticesListes[kk])[p0[1]]
2940          std::cout<<"-------------------------------------------------"<<"\n";                        =((*(*__internals).verticesListesIntersection[kk])[p0[1]]);
2941          Vertex& V0 = *((*(*__internals).verticesListes[0])[(*currentTriangle)[0]]);                    }
2942          Vertex& V1 = *((*(*__internals).verticesListes[0])[(*currentTriangle)[1]]);                }
2943          Vertex& V2 = *((*(*__internals).verticesListes[0])[(*currentTriangle)[2]]);                else
2944          std::cout<<"VO "<<V0.Vector3()[0]<<"  "<<V0.Vector3()[1]<<" "<<V0.Vector3()[2]<<"\n";                    ffout(0)<<"pbs point 0 pas dans vertices\n";
         std::cout<<"V1 "<<V1.Vector3()[0]<<"  "<<V1.Vector3()[1]<<" "<<V1.Vector3()[2]<<"\n";  
         std::cout<<"V2 "<<V2.Vector3()[0]<<"  "<<V2.Vector3()[1]<<" "<<V2.Vector3()[2]<<"\n";  
         std::cout<<" "<<"\n";  
   
         // Pour Cécile :  
         //  Cette partie m'a l'air assez dangereuse : tu es en train  
         // de modifier un caractère essentiel de certains points du  
         // maillage : s'ils sont ou non à l'intérieur du domaine !  
         // (à cause de la référence bool&)  
   
         bool i0=(ii0==1);  
         bool i1 = (ii1==1);  
         bool i2 = (ii2==1);  
   
         std::cout<<"inside? "<<i0<<" "<<i1<<" "<<i2<<"\n";  
   
         if(tn0+tn1==2) {  
           if(i0) {  
             Vertex A=V0;  
 //std::cout<<"A "<<A.Vector3()[0]<<"  "<<A.Vector3()[1]<<" "<<A.Vector3()[2]<<"\n";  
             Vertex& B=(*(listofVertex[e0]));  
             Vertex& C=(*(listofVertex[e1]));  
 //std::cout<<"B "<<B.Vector3()[0]<<"  "<<B.Vector3()[1]<<" "<<B.Vector3()[2]<<"\n";  
 //std::cout<<"C "<<C.Vector3()[0]<<"  "<<C.Vector3()[1]<<" "<<C.Vector3()[2]<<"\n";  
             ffout(0) << "----------\n";  
             vertexplot.push_back(A);  
             vertexplot.push_back(B);  
             vertexplot.push_back(C);  
               
           } else if(i1 and i2) {  
             Vertex A=V1;  
             Vertex B=V2;  
             Vertex& C=(*(listofVertex[e0]));  
             Vertex& D=(*(listofVertex[e1]));  
             vertexplot.push_back(A);  
             vertexplot.push_back(C);  
             vertexplot.push_back(D);  
             vertexplot.push_back(A);  
             vertexplot.push_back(B);  
             vertexplot.push_back(D);  
           } else {  
             std::cout<<"erreur1\n";  
2945            }            }
2946          }            if(p1[0]==p1[1]) {
2947          // else if(triangleedge[(*currentTriangle)[0]]+triangleedge[(*currentTriangle)[2]]==2)                if((*(*__internals).verticesListesIntersection[kk]).find(p1[1])!=
2948          else if(tn0+tn2==2) {            (*(*__internals).verticesListesIntersection[kk]).end()) {
2949            if(i1) {                    if((*(*__internals).verticesListes[kk]).find(p1[1])==(*(*__internals).verticesListes[kk]).end())
2950              Vertex A=V1;                    {
2951              Vertex& B=(*(listofVertex[e0]));                        (*(*__internals).verticesListes[kk])[p1[1]]=
2952              Vertex& C=(*(listofVertex[e2]));                        ((*(*__internals).verticesListesIntersection[kk])[p1[1]]);
2953              vertexplot.push_back(A);                    }
2954              vertexplot.push_back(B);                }
2955              vertexplot.push_back(C);                else
2956              //std::cout<<"A "<<A.Vector3()[0]<<"  "<<A.Vector3()[1]<<" "<<A.Vector3()[2]<<"\n";                    ffout(0)<<"pbs point 1 pas dans vertices\n";
             //std::cout<<"B "<<B.Vector3()[0]<<"  "<<B.Vector3()[1]<<" "<<B.Vector3()[2]<<"\n";  
             //std::cout<<"C "<<C.Vector3()[0]<<"  "<<C.Vector3()[1]<<" "<<C.Vector3()[2]<<"\n";  
           } else if(i0 and i2) {  
             Vertex A=V0;  
             Vertex B=V2;  
             Vertex& C=(*(listofVertex[e0]));  
             Vertex& D=(*(listofVertex[e2]));  
             vertexplot.push_back(A);  
             vertexplot.push_back(C);  
             vertexplot.push_back(D);  
             vertexplot.push_back(A);  
             vertexplot.push_back(B);  
             vertexplot.push_back(D);  
           } else {  
             std::cout<<"erreur2\n";  
2957            }            }
2958          } else {            if(p2[0]==p2[1]) {
2959            if(i2) {                if((*(*__internals).verticesListesIntersection[kk]).find(p2[1])!=
2960              Vertex A=V2;            (*(*__internals).verticesListesIntersection[kk]).end()) {
2961              Vertex& B=(*(listofVertex[e1]));                    if((*(*__internals).verticesListes[kk]).find(p2[1])==(*(*__internals).verticesListes[kk]).end())
2962              Vertex& C=(*(listofVertex[e2]));                    {
2963              vertexplot.push_back(A);                        (*(*__internals).verticesListes[kk])[p2[1]]=
2964              vertexplot.push_back(B);                        ((*(*__internals).verticesListesIntersection[kk])[p2[1]]);
2965              vertexplot.push_back(C);                    }
2966           //   std::cout<<"A "<<A.Vector3()[0]<<"  "<<A.Vector3()[1]<<" "<<A.Vector3()[2]<<"\n";                }
2967              //std::cout<<"B "<<B.Vector3()[0]<<"  "<<B.Vector3()[1]<<" "<<B.Vector3()[2]<<"\n";                else
2968              //std::cout<<"C "<<C.Vector3()[0]<<"  "<<C.Vector3()[1]<<" "<<C.Vector3()[2]<<"\n";                    ffout(0)<<"pbs point 0 pas dans vertices\n";
           } else if(i1 and i0) {  
             Vertex A=V0;  
             Vertex B=V1;  
             std::cout<<"A "<<A.Vector3()[0]<<"  "<<A.Vector3()[1]<<" "<<A.Vector3()[2]<<"\n";  
             std::cout<<"B "<<B.Vector3()[0]<<"  "<<B.Vector3()[1]<<" "<<B.Vector3()[2]<<"\n";  
             Vertex& C=(*(listofVertex[e1]));  
             Vertex& D=(*(listofVertex[e2]));  
             std::cout<<"C "<<C.Vector3()[0]<<"  "<<C.Vector3()[1]<<" "<<C.Vector3()[2]<<"\n";  
             vertexplot.push_back(A);  
             vertexplot.push_back(C);  
             vertexplot.push_back(D);  
             vertexplot.push_back(A);  
             vertexplot.push_back(B);  
             vertexplot.push_back(D);  
           } else {  
             std::cout<<"erreur3\n";  
2969            }            }
         }  
2970        }        }
2971        if(tn0+tn1+tn2==3) {    }
         Vertex& V0 = *((*(*__internals).verticesListes[0])[(*currentTriangle)[0]]);  
         Vertex& V1 = *((*(*__internals).verticesListes[0])[(*currentTriangle)[1]]);  
         Vertex& V2 = *((*(*__internals).verticesListes[0])[(*currentTriangle)[2]]);  
2972    
2973          // Pour Cécile : comme plus haut...    ffout(9)<<"apres nb vertices object 0 "<<(*(*__internals).verticesListes[0]).size()<<"\n";
2974      ffout(9)<<"apres nb vertices object 1 "<<(*(*__internals).verticesListes[1]).size()<<"\n";
2975    
2976          bool i0=(ii0==1);  #warning remove this
2977          bool i1 = (ii1==1);    (*__internals).plotliste(nobject);
         bool i2 = (ii2==1);  
   
         std::cout<<"be careful triangle \n";  
         //std::cout<<"inside? "<<V0.Ref().insideDomain()<<" "<<V1.Ref().insideDomain()<<" "<<V2.Ref().insideDomain()<<"\n";  
         if(i0+i1+i2==1)  
         {  
         if(i0) {  
           std::cout<<"be careful triangle 0\n";  
           Vertex A=V0;  
           Vertex& B=(*(listofVertex[e0]));  
           Vertex& C=(*(listofVertex[e1]));  
           Vertex& D=(*(listofVertex[e2]));  
           vertexplot.push_back(A);  
           vertexplot.push_back(B);  
           vertexplot.push_back(D);  
           vertexplot.push_back(A);  
           vertexplot.push_back(C);  
           vertexplot.push_back(D);  
           //std::cout<<"A "<<A.Vector3()[0]<<"  "<<A.Vector3()[1]<<" "<<A.Vector3()[2]<<"\n";  
           //std::cout<<"B "<<B.Vector3()[0]<<"  "<<B.Vector3()[1]<<" "<<B.Vector3()[2]<<"\n";  
           //std::cout<<"C "<<C.Vector3()[0]<<"  "<<C.Vector3()[1]<<" "<<C.Vector3()[2]<<"\n";  
           //std::cout<<"D "<<D.Vector3()[0]<<"  "<<D.Vector3()[1]<<" "<<D.Vector3()[2]<<"\n";  
         } else if(i1) {  
           std::cout<<"be careful triangle 1 \n";  
           Vertex A=V1;  
           Vertex& B=(*(listofVertex[e0]));  
           Vertex& C=(*(listofVertex[e1]));  
           Vertex& D=(*(listofVertex[e2]));  
           vertexplot.push_back(A);  
           vertexplot.push_back(C);  
           vertexplot.push_back(D);  
           vertexplot.push_back(A);  
           vertexplot.push_back(C);  
           vertexplot.push_back(B);  
         } else if(i2) {  
           std::cout<<"be careful triangle 2 \n";  
           Vertex A=V2;  
           Vertex& B=(*(listofVertex[e0]));  
           Vertex& C=(*(listofVertex[e1]));  
           Vertex& D=(*(listofVertex[e2]));  
           vertexplot.push_back(A);  
           vertexplot.push_back(C);  
           vertexplot.push_back(D);  
           vertexplot.push_back(A);  
           vertexplot.push_back(B);  
           vertexplot.push_back(D);  
           //std::cout<<"A "<<A.Vector3()[0]<<"  "<<A.Vector3()[1]<<" "<<A.Vector3()[2]<<"\n";  
           //std::cout<<"B "<<B.Vector3()[0]<<"  "<<B.Vector3()[1]<<" "<<B.Vector3()[2]<<"\n";  
           //std::cout<<"C "<<C.Vector3()[0]<<"  "<<C.Vector3()[1]<<" "<<C.Vector3()[2]<<"\n";  
           //std::cout<<"D "<<D.Vector3()[0]<<"  "<<D.Vector3()[1]<<" "<<D.Vector3()[2]<<"\n";  
         }  
         }  
         else  
         {  
             if(i0+i1==2)  
             {  
                 Vertex A=V0;  
                 Vertex B=V1;  
                 Vertex& C=(*(listofVertex[e1]));  
                 Vertex& D=(*(listofVertex[e2]));  
                 vertexplot.push_back(A);  
                 vertexplot.push_back(C);  
                 vertexplot.push_back(D);  
                 vertexplot.push_back(A);  
                 vertexplot.push_back(B);  
                 vertexplot.push_back(D);  
             }  
             else if(i0+i2==2)  
             {  
                 Vertex A=V0;  
                 Vertex B=V2;  
                 Vertex& C=(*(listofVertex[e0]));  
                 Vertex& D=(*(listofVertex[e2]));  
                 vertexplot.push_back(A);  
                 vertexplot.push_back(C);  
                 vertexplot.push_back(D);  
                 vertexplot.push_back(A);  
                 vertexplot.push_back(B);  
                 vertexplot.push_back(D);  
             }  
             else if(i1+i2==2)  
             {  
                 Vertex A=V1;  
                 Vertex B=V2;  
                 Vertex& C=(*(listofVertex[e0]));  
                 Vertex& D=(*(listofVertex[e1]));  
                 vertexplot.push_back(A);  
                 vertexplot.push_back(C);  
                 vertexplot.push_back(D);  
                 vertexplot.push_back(A);  
                 vertexplot.push_back(B);  
                 vertexplot.push_back(D);  
             }  
             else  
 std::cout<<"3 points in\n";  
         }  
       }  
       nn++;  
     }  
   }*/  
2978    
2979    size_t ntr=0;    size_t ntr=0;
2980    for(size_t n = 0; n<scene.nbObjects(); ++n)  
2981      for(size_t n = 0; n<nobject; ++n)
2982      ntr+=(*(*__internals).verticesListes[n]).size();      ntr+=(*(*__internals).verticesListes[n]).size();
2983    
2984      ntr+=(*__internals).listVerticesIntersectionPerEdges.size();
2985      
2986    s_mesh.setNumberOfVertices(ntr);    s_mesh.setNumberOfVertices(ntr);
2987      ffout(9)<<"nb vertices intersection "<<(*__internals).listVerticesIntersectionPerEdges.size()<<"\n";
2988    // //  s_mesh.setNumberOfVertices((*__internals).verticesList.size());    ffout(9)<<"nb tot vertices "<<s_mesh.numberOfVertices()<<"\n";
   
   std::cout<<"nb tot "<<s_mesh.numberOfVertices()<<"\n";  
2989    
2990    //! copies vertices into the s_mesh.    //! copies vertices into the s_mesh.
2991    int i = 0;    int i = 0;
2992    ////   for(SurfaceMeshGenerator::Internals::VerticesList::iterator splitVertex    for(size_t n = 0; n<scene.nbObjects(); ++n) {
   ////  =  (*__internals).verticesList.begin();  
   ////  splitVertex != (*__internals).verticesList.end(); ++splitVertex) {  
   
   // //s_mesh.vertex(i) = (*(*splitVertex).second);  
   ////delete (*splitVertex).second;  
   ////(*splitVertex).second = &s_mesh.vertex(i);  
   ////i++;  
   ////}  
   for(size_t n = 0; n<scene.nbObjects(); ++n)  
     {  
2993        for(SurfaceMeshGenerator::Internals::VerticesList::iterator splitVertex        for(SurfaceMeshGenerator::Internals::VerticesList::iterator splitVertex
2994              =  (*(*__internals).verticesListes[n]).begin();              =  (*(*__internals).verticesListes[n]).begin();
2995            splitVertex != (*(*__internals).verticesListes[n]).end(); ++splitVertex)            splitVertex != (*(*__internals).verticesListes[n]).end(); ++splitVertex) {
         {  
2996            s_mesh.vertex(i) = (*(*splitVertex).second);            s_mesh.vertex(i) = (*(*splitVertex).second);
2997            delete (*splitVertex).second;            delete (*splitVertex).second;
2998            (*splitVertex).second = &s_mesh.vertex(i);            (*splitVertex).second = &s_mesh.vertex(i);
2999            i++;            i++;
3000          }          }
3001      }      }
   
   
   
   /*std::set<const Cell*> toTreatHexahedra;  
   std::set_difference (cuttedHexahedra.begin(), cuttedHexahedra.end(),  
                        meshedHexahedra.begin(), meshedHexahedra.end(),  
                        inserter(toTreatHexahedra, toTreatHexahedra.begin()));*/  
3002        
3003      for(SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator
3004          n=(*__internals).listVerticesIntersectionPerEdges.begin() ;
3005          n!=(*__internals).listVerticesIntersectionPerEdges.end() ; ++n) {
3006          s_mesh.vertex(i)=(*((*n).second));
3007          delete (*n).second;
3008          (*n).second = &s_mesh.vertex(i);
3009          i++;
3010      }
3011    
3012    i = 0;    i = 0;
3013    size_t ncell=0;    size_t ncell=0;
3014    
3015    for(size_t n = 0; n<scene.nbObjects(); ++n)    for(size_t n = 0; n<scene.nbObjects(); ++n)
3016      ncell+=(*(*__internals).triangleListes[n]).size();      ncell+=(*(*__internals).triangleListes[n]).size();
3017    
3018    s_mesh.setNumberOfCells(ncell);    size_t ntr1=(*(*__internals).triangleListes[0]).size();
3019    std::cout<<"nb tot "<<s_mesh.numberOfCells()<<" \n";    for(size_t no=0 ; no<nobject ; ++no) {
3020          ncell+=(*(*__internals).listTriangleintersection[no]).size();
3021      }
3022      size_t ninter=(*(*__internals).listTriangleintersection[0]).size();
3023      size_t ninter2=(*(*__internals).listTriangleintersection[1]).size();
3024    
3025      s_mesh.setNumberOfCells(ncell);
3026      ffout(9)<<"nb tot cells "<<s_mesh.numberOfCells()<<" \n";
3027    
3028    SurfaceMeshGenerator::Internals::MotherCellList::iterator motherCell    SurfaceMeshGenerator::Internals::MotherCellList::iterator motherCell
3029      = (*__internals).cellList.begin();      = (*__internals).cellList.begin();
3030    
3031    
3032    for(size_t n = 0; n<scene.nbObjects(); ++n)    for(size_t n = 0; n<scene.nbObjects(); ++n) {
3033      {      for(SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle
3034              = (*(*__internals).triangleListes[n]).begin();
3035        for(SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle          currentTriangle != (*(*__internals).triangleListes[n]).end();
3036              = (*(*__internals).triangleListes[n]).begin();          ++currentTriangle) {
3037            currentTriangle != (*(*__internals).triangleListes[n]).end();        Vertex& V0 = *((*(*__internals).verticesListes[n])[(*currentTriangle)[0]]);
3038            ++currentTriangle)        Vertex& V1 = *((*(*__internals).verticesListes[n])[(*currentTriangle)[1]]);
3039          {        Vertex& V2 = *((*(*__internals).verticesListes[n])[(*currentTriangle)[2]]);
           //Vertex& V0 = *((*__internals).verticesList[(*currentTriangle)[0]]);  
           //Vertex& V1 = *((*__internals).verticesList[(*currentTriangle)[1]]);  
           //Vertex& V2 = *((*__internals).verticesList[(*currentTriangle)[2]]);  
           Vertex& V0 = *((*(*__internals).verticesListes[n])[(*currentTriangle)[0]]);  
           Vertex& V1 = *((*(*__internals).verticesListes[n])[(*currentTriangle)[1]]);  
           Vertex& V2 = *((*(*__internals).verticesListes[n])[(*currentTriangle)[2]]);  
3040    
3041            s_mesh.cell(i) = Triangle(V0,V1,V2);          
3042          s_mesh.cell(i) = Triangle(V0,V1,V2);
3043  #warning SHOULD USE BC NUMBER (ie OBJECT NUMBER IN SCENE INSTEAD)  #warning SHOULD USE BC NUMBER (ie OBJECT NUMBER IN SCENE INSTEAD)
3044          s_mesh.cell(i).reference() = 1;
3045          s_mesh.cell(i).setMother(*motherCell);
3046          i++;
3047          motherCell++;
3048        }
3049      }
3050    
3051      for(size_t kk=0 ; kk<nobject ; ++kk) {
3052          for(size_t n=0 ; n<(*(*__internals).listTriangleintersection[kk]).size(); ++n) {
3053    
3054              TinyVector<2,Edge::Pair> p0=((*(*__internals).listTriangleintersection[kk]))[n][0];
3055              TinyVector<2,Edge::Pair> p1=((*(*__internals).listTriangleintersection[kk]))[n][1];
3056              TinyVector<2,Edge::Pair> p2=((*(*__internals).listTriangleintersection[kk]))[n][2];
3057              Vertex *V0,*V1,*V2;
3058              std::pair<TinyVector<2,Edge::Pair>,TinyVector<2,size_t> > pp0,pp1,pp2;
3059              TinyVector<2,size_t> nuum;
3060              nuum[0]=0;
3061              nuum[1]=0;
3062              pp0.second=nuum;
3063              pp1.second=nuum;
3064              pp2.second=nuum;
3065              pp0.first=p0;
3066              pp1.first=p1;
3067              pp2.first=p2;
3068    
3069              if(p0[0]==p0[1]) {
3070                  V0=((*(*__internals).verticesListes[kk])[p0[0]]);
3071              }
3072              else {
3073                  SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator
3074                  ip0=(*__internals).listVerticesIntersectionPerEdges.find(pp0);
3075                  if(ip0!=(*__internals).listVerticesIntersectionPerEdges.end())
3076                      V0 =(*ip0).second;
3077              }
3078              if(p1[0]==p1[1]) {
3079                  V1=((*(*__internals).verticesListes[kk])[p1[0]]);
3080              }
3081              else {
3082                  SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator
3083                  ip0=(*__internals).listVerticesIntersectionPerEdges.find(pp1);
3084                  if(ip0!=(*__internals).listVerticesIntersectionPerEdges.end())
3085                      V1 =(*ip0).second;
3086              }
3087              if(p2[0]==p2[1]) {
3088                  V2=((*(*__internals).verticesListes[kk])[p2[0]]);
3089              }
3090              else {
3091                  SurfaceMeshGenerator::Internals::ListVerticesIntersectionPerEdges::iterator
3092                  ip0=(*__internals).listVerticesIntersectionPerEdges.find(pp2);
3093                  if(ip0!=(*__internals).listVerticesIntersectionPerEdges.end())
3094                      V2 =(*ip0).second;
3095              }
3096              s_mesh.cell(i) = Triangle(*V0,*V1,*V2);
3097    
3098    #warning SHOULD USE BC NUMBER (ie OBJECT NUMBER IN SCENE INSTEAD)
3099            s_mesh.cell(i).reference() = 1;            s_mesh.cell(i).reference() = 1;
3100            s_mesh.cell(i).setMother(*motherCell);            s_mesh.cell(i).setMother(*motherCell);
   
3101            i++;            i++;
3102            motherCell++;            motherCell++;
3103          }        }
3104      }    }
   
   //  plot2(s_mesh, toTreatHexahedra, SMesh);  
 /*ffout(0) << "plot medit \n";  
   plotmedit(s_mesh, toTreatHexahedra,nobj);  
 ffout(0) << "plot sortie \n";  
   {  
     std::ofstream sortie("plot.mesh");  
3105    
3106      sortie<<"MeshVersionFormatted 1"<<" \n";    //  plot2(s_mesh, (*__internals).toTreatHexahedra, SMesh);
3107      sortie<<" \n";    ffout(9) << "plot medit \n";
3108      sortie<<"Dimension"<<" \n";    plotmedit(nobject,s_mesh, (*__internals).toTreatHexahedra,cuttedHexahedra,
3109      sortie<<"3"<<" \n";              (*__internals).verticesListesIntersection,(*__internals).triangleListesIntersectionEdges,
3110      sortie<<" \n";              nobjother,ninter,ninter2,ntr1);
     sortie<<"Vertices"<<" \n";  
     sortie<<vertexplot.size()<<"\n";  
     for(size_t i=0;i<vertexplot.size();++i)  
       {  
         sortie<<(*(&vertexplot[i])).Vector3()[0]<<" "<<(*(&vertexplot[i])).Vector3()[1]<<" "<<(*(&vertexplot[i])).Vector3()[2]<<" 0"<<" \n";  
       }  
   
     sortie<<" \n";  
     sortie<<"Triangles"<<"\n";  
     sortie<<vertexplot.size()/3<<"\n";  
   
     for (size_t i=0; i<vertexplot.size()/3; ++i)  
       {  
         sortie<<3*i+1<<" "<<3*i+2<<" "<<3*i+3<<" "<<1<<"\n";  
       }  
     sortie<<" \n";  
     sortie<<"End"<<"\n";  
   }*/  
3111  }  }
3112    
   
   
3113  SurfaceMeshGenerator::SurfaceMeshGenerator()  SurfaceMeshGenerator::SurfaceMeshGenerator()
3114    : __internals(new SurfaceMeshGenerator::Internals())    : __internals(new SurfaceMeshGenerator::Internals())
3115  {  {
# Line 2171  SurfaceMeshGenerator::~SurfaceMeshGenera Line 3121  SurfaceMeshGenerator::~SurfaceMeshGenera
3121    ;    ;
3122  }  }
3123    
3124    void plothexa(size_t ncase,std::set<const Cell*>& cuttedHexahedra) {
3125        char* cccc="hhhhhhhhhhhhh";
3126        if(ncase==0) {
3127            cccc="meshed1.mesh";
3128        }
3129        else if(ncase==1) {
3130            cccc="meshed2.mesh";
3131        }
3132        else {
3133            cccc="meshed3.mesh";
3134        }
3135        std::ofstream sortie3(cccc);
3136        sortie3<<"MeshVersionFormatted 1"<<" \n";
3137        sortie3<<" \n";
3138        sortie3<<"Dimension"<<" \n";
3139        sortie3<<"3"<<" \n";
3140        sortie3<<" \n";
3141        sortie3<<"Vertices"<<" \n";
3142        sortie3<<8*cuttedHexahedra.size()<<" \n";
3143        for (std::set<const Cell*>::iterator i = cuttedHexahedra.begin();
3144             i != cuttedHexahedra.end(); ++i)
3145        {
3146            const CartesianHexahedron& H = static_cast<const CartesianHexahedron&>(*(*i));
3147            for (size_t j=0; j<8; ++j)
3148            {
3149                sortie3<<H(j)[0]<<" "<<H(j)[1]<<" "<<H(j)[2]<<" 0"<<" \n";
3150            }
3151        }
3152        sortie3<<" \n";
3153        sortie3<<"Hexahedra"<<"\n";
3154        sortie3<<cuttedHexahedra.size()<<" \n";
3155    
3156        for(size_t i=0;i<cuttedHexahedra.size();i++)
3157            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";
3158        sortie3<<"End"<<" \n";
3159        
3160    }
3161    
3162        
3163    
3164    
3165  void plotmedit(const SurfaceMeshOfTriangles& s_mesh,std::set<const Cell*>& toTreatHexahedra,size_t nobj)  void plotmedit(const size_t& nobject,
3166  {                 const SurfaceMeshOfTriangles& s_mesh,
3167                   std::set<const Cell*>& toTreatHexahedra,
3168                   std::set<const Cell*>& cuttedHexahedra,
3169                   std::vector<ReferenceCounting<std::map<Edge::Pair, Vertex*> > > verticesListesIntersection,
3170                   std::vector<ReferenceCounting<std::list<TinyVector<3,TinyVector<2,Edge::Pair> > > > >
3171                   triangleListesIntersectionEdges,
3172                   size_t nobjother,size_t ninter,size_t ninter2,size_t ntr1) {
3173    
3174    // **********************************************************************    // **********************************************************************
3175        std::ofstream sortie3("cutted.mesh");
3176    
3177        sortie3<<"MeshVersionFormatted 1"<<" \n";
3178        sortie3<<" \n";
3179        sortie3<<"Dimension"<<" \n";
3180        sortie3<<"3"<<" \n";
3181        sortie3<<" \n";
3182        sortie3<<"Vertices"<<" \n";
3183        sortie3<<8*cuttedHexahedra.size()<<" \n";
3184        for (std::set<const Cell*>::iterator i = cuttedHexahedra.begin();
3185             i != cuttedHexahedra.end(); ++i)
3186        {
3187            const CartesianHexahedron& H = static_cast<const CartesianHexahedron&>(*(*i));
3188            for (size_t j=0; j<8; ++j)
3189            {
3190                sortie3<<H(j)[0]<<" "<<H(j)[1]<<" "<<H(j)[2]<<" 0"<<" \n";
3191            }
3192        }
3193        sortie3<<" \n";
3194        sortie3<<"Hexahedra"<<"\n";
3195        sortie3<<cuttedHexahedra.size()<<" \n";
3196    
3197        for(size_t i=0;i<cuttedHexahedra.size();i++)
3198            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";
3199        sortie3<<"End"<<" \n";
3200        
3201    
3202        
3203    std::ofstream sortie("intersection.mesh");    std::ofstream sortie("intersection.mesh");
3204                
3205    sortie<<"MeshVersionFormatted 1"<<" \n";    sortie<<"MeshVersionFormatted 1"<<" \n";
# Line 2202  void plotmedit(const SurfaceMeshOfTriang Line 3222  void plotmedit(const SurfaceMeshOfTriang
3222      {      {
3223        sortie<<s_mesh.vertex(i)[0]<<" "<<s_mesh.vertex(i)[1]<<" "<<s_mesh.vertex(i)[2]<<" 0"<<" \n";        sortie<<s_mesh.vertex(i)[0]<<" "<<s_mesh.vertex(i)[1]<<" "<<s_mesh.vertex(i)[2]<<" 0"<<" \n";
3224      }      }
3225    //       for (size_t i=0; i<s_mesh.numberOfCells(); ++i)  
3226    //       {    /*        for (size_t i=0; i<s_mesh.numberOfCells(); ++i)
3227    //       const Triangle& T = s_mesh.cell(i);              {
3228    //       sortie<<T(0)[0]<<" "<<T(0)[1]<<" "<<T(0)[2]<<" 0"<<" \n";              const Triangle& T = s_mesh.cell(i);
3229    //       sortie<<T(1)[0]<<" "<<T(1)[1]<<" "<<T(1)[2]<<" 0"<<" \n";              sortie<<T(0)[0]<<" "<<T(0)[1]<<" "<<T(0)[2]<<" 0"<<" \n";
3230    //       sortie<<T(2)[0]<<" "<<T(2)[1]<<" "<<T(2)[2]<<" 0"<<" \n";              sortie<<T(1)[0]<<" "<<T(1)[1]<<" "<<T(1)[2]<<" 0"<<" \n";
3231    //              sortie<<T(2)[0]<<" "<<T(2)[1]<<" "<<T(2)[2]<<" 0"<<" \n";
3232    //       }    
3233                }*/
3234    sortie<<" \n";    sortie<<" \n";
3235    sortie<<"Triangles"<<"\n";    sortie<<"Triangles"<<"\n";
3236    sortie<<s_mesh.numberOfCells()<<"\n";    sortie<<s_mesh.numberOfCells()<<"\n";
3237      ffout(9)<<"fin ecriture vertex\n";
3238    for (size_t i=0; i<s_mesh.numberOfCells(); ++i)    for (size_t i=0; i<s_mesh.numberOfCells(); ++i)
3239      {      {
3240        const Triangle& T = s_mesh.cell(i);        const Triangle& T = s_mesh.cell(i);
3241         // ffout(0)<<s_mesh.vertexNumber(T(0))<<" "<<s_mesh.vertexNumber(T(1))<<" "<<s_mesh.vertexNumber(T(2))<<"\n";
3242        sortie<<s_mesh.vertexNumber(T(0))+8*toTreatHexahedra.size()+1<<" "<<s_mesh.vertexNumber(T(1))+8*toTreatHexahedra.size()+1<<" "<<s_mesh.vertexNumber(T(2))+8*toTreatHexahedra.size()+1;        sortie<<s_mesh.vertexNumber(T(0))+8*toTreatHexahedra.size()+1<<" "<<s_mesh.vertexNumber(T(1))+8*toTreatHexahedra.size()+1<<" "<<s_mesh.vertexNumber(T(2))+8*toTreatHexahedra.size()+1;
3243        //          sortie<<3*i+1+8*toTreatHexahedra.size()<<" "<<3*i+2+8*toTreatHexahedra.size()<<" //"<<3*i+3+8*toTreatHexahedra.size()<<" "<<1<<"\n";        
3244        if(i<(nobj))  //                sortie<<3*i+1+8*toTreatHexahedra.size()<<" "<<3*i+2+8*toTreatHexahedra.size()<<" //"<<3*i+3+8*toTreatHexahedra.size()<<" "<<1<<"\n";
3245          sortie<<" 1"<<"\n";        if(i<ntr1)
3246              sortie<<" 1"<<"\n";
3247          else if(i<(nobjother))
3248              sortie<<" 2"<<"\n";
3249          else if(i<(nobjother+ninter))
3250              sortie<<" 3"<<"\n";
3251          else if(i<(nobjother+ninter+ninter2))
3252              sortie<<" 4"<<"\n";
3253        else        else
3254          sortie<<" 2"<<"\n";            sortie<<" 5"<<"\n";
3255      }      }
3256      ffout(9)<<"fin ecriture triangles\n";
3257    sortie<<" \n";    sortie<<" \n";
3258    sortie<<"Hexahedra"<<"\n";    sortie<<"Hexahedra"<<"\n";
3259    sortie<<toTreatHexahedra.size()<<" \n";    sortie<<toTreatHexahedra.size()<<" \n";
3260    
   /*  for (std::set<const Cell*>::iterator i = toTreatHexahedra.begin();  
       i != toTreatHexahedra.end(); ++i)  
       {  
       const CartesianHexahedron& H = static_cast<const CartesianHexahedron&>(*(*i));  
       for(size_t j=0; j<8; ++j)  
       sortie<<SMesh.vertexNumber(H(j))<<" \n";  
       }*/  
3261    for(size_t i=0;i<toTreatHexahedra.size();i++)    for(size_t i=0;i<toTreatHexahedra.size();i++)
3262      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";
3263    sortie<<"End"<<" \n";    sortie<<"End"<<" \n";
3264      ffout(9)<<"fin ecriture intersection.mesh\n";
3265    
3266  }  }
3267    
3268    
3269    #warning REMOVE THIS
3270    void SurfaceMeshGenerator::Internals::__getIntersectionReferences(const Intersection& I,
3271                                                                      std::map<TinyVector<3>, Object*, SurfaceMeshGenerator::Internals::__ltv>& otherReferences)
3272    {
3273      for(Intersection::const_iterator i = I.begin();
3274          i != I.end(); ++i) {
3275        if ((*(*i)).shapeType() == Shape::intersection) {
3276          const Intersection& I2 = static_cast<const Intersection&>(*(*(*i)).shape());
3277          SurfaceMeshGenerator::Internals::__getIntersectionReferences(I2, otherReferences);
3278        } else {
3279          otherReferences[(*(*i)).reference()] = *i;
3280          ffout(0) << (*(*i)).reference() << "\n";
3281        }
3282      }
3283    }
3284    
3285  // #include <vtkPolyDataMapper.h>  // #include <vtkPolyDataMapper.h>
3286  // #include <vtkRenderWindow.h>  // #include <vtkRenderWindow.h>
3287  // #include <vtkRenderWindowInteractor.h>  // #include <vtkRenderWindowInteractor.h>

Legend:
Removed from v.1.10  
changed lines
  Added in v.1.11

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