/[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.51 by delpinux, Wed Oct 22 17:08:27 2003 UTC revision 1.52 by dobj, Tue Oct 28 18:42:10 2003 UTC
# Line 55  Line 55 
55    
56  #include <MeshExtractor.hpp>  #include <MeshExtractor.hpp>
57    
58    #include <triangulation.hpp>
59  #include <set>  #include <set>
60  #include <iostream>  #include <iostream>
61  #include <stack>  #include <stack>
62    
63    
64    inline double SDet(const Vertex & V1, const Vertex & V2, const Vertex & V3) {
65        return (V2[0]-V1[0]) * (V3[1]-V1[1]) -  (V2[1]-V1[1]) * (V3[0]-V1[0]);
66    }
67    
68    inline bool checkInterbis(const Vertex &P1, const Vertex &P2, const Vertex &Q1, const Vertex &Q2) {
69        // retourne true si P1P2 inter Q1Q2 != vide
70        assert(!(P1 == P2 || P1 == Q1 || P1 == Q2 || P2 == Q1 || P2 == Q2 || Q1 == Q2));
71        const double eps = 1e-15;
72        /*if (std::abs(SDet(Q1,Q2,P1)) <eps||
73            std::abs(SDet(Q1,Q2,P2)) <eps) {
74            //ffout(0) << " !!!!! NO !!!!!!!" << std::endl;
75            return true;
76        }*/
77    
78        return ((SDet(P1,P2,Q1)*SDet(P1,P2,Q2) < 0) && (SDet(Q1,Q2,P1)*SDet(Q1,Q2,P2) < 0)); // version optimisée
79    }
80  class SurfaceMeshGenerator::Internals  class SurfaceMeshGenerator::Internals
81  {  {
82  private:  private:
# Line 179  private: Line 197  private:
197    typedef std::list<const Cell*> MotherCellList;    typedef std::list<const Cell*> MotherCellList;
198    
199    typedef std::map<const Cell*,    typedef std::map<const Cell*,
200                     std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,      std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,
201                               ReferenceCounting<std::vector<Edge> > > > HexaTetraAssociation;      ReferenceCounting<std::vector<Edge> > > > HexaTetraAssociation;
202    
203    typedef std::map<TinyVector<3>, Object*> ObjectReferences;    typedef std::map<TinyVector<3>, Object*> ObjectReferences;
204    
# Line 198  private: Line 216  private:
216    //std::vector<Vertex*> listVertexIntersectionNew;    //std::vector<Vertex*> listVertexIntersectionNew;
217    
218    typedef std::map<std::pair<const Vertex*,const Vertex*> ,    typedef std::map<std::pair<const Vertex*,const Vertex*> ,
219                     TinyVector<2,const Vertex* > > PIntersection;      TinyVector<2,const Vertex* > > PIntersection;
220    
221    typedef std::map< PairTriangleNew , IntersectionPoints> PairTriangleToIntersectionPoints;    typedef std::map< PairTriangleNew , IntersectionPoints> PairTriangleToIntersectionPoints;
222    
# Line 210  private: Line 228  private:
228    std::vector<ReferenceCounting< std::map<const Vertex*,size_t> > > edgesRefVertex;    std::vector<ReferenceCounting< std::map<const Vertex*,size_t> > > edgesRefVertex;
229    
230    std::map<const Vertex*,TinyVector<2,const Triangle*> > vertexInTriangles;    std::map<const Vertex*,TinyVector<2,const Triangle*> > vertexInTriangles;
231      std::map<const Vertex*,std::vector<const Triangle*> > vertexVectTriangles;
232    std::map<const Triangle*,TinyVector<2,const Vertex*> > triangleWithVertex;    std::map<const Triangle*,TinyVector<2,const Vertex*> > triangleWithVertex;
233    
234    TinyVector<3> __splitEdge(const Edge& e,    TinyVector<3> __splitEdge(const Edge& e,
# Line 302  private: Line 321  private:
321                         const SurfaceMeshOfTriangles& surfmesh,                         const SurfaceMeshOfTriangles& surfmesh,
322                         const Connectivity<Triangle>& connect,                         const Connectivity<Triangle>& connect,
323                         const std::set<const Cell*>& toTreatHexahedra);                         const std::set<const Cell*>& toTreatHexahedra);
     
   void __createByEdges(const TriangleCut& K,const Cell* cell,  
                        std::vector<Triangle>& triangleListIntersectionNew);  
324    
325      void __createCase(const TriangleCut& K ,
326                        size_t& in, size_t& in1,
327                        TinyVector<3,size_t>& place);
328    
329      void __transformVertex(const TriangleCut& K,
330                             std::vector<const Vertex*>& points,
331                             std::vector<Vertex> &pointstrian);
332      void __addPoints(const Vertex*& Pcut,
333                       size_t& stopcut,
334                       std::vector<const Vertex*>& pointslist1,
335                       std::vector<const Vertex*>& pointslist2);
336      
337      void __addList(const TriangleCut& K);
338      void __createList(const TriangleCut& K,
339                        std::vector<const Vertex*>& points,
340                        size_t& in,
341                        size_t& in1,
342                        TinyVector<3,size_t> place);
343    
344      bool __verifExist(const Vertex* Ptemp,
345                        std::vector<const Vertex*>& points);
346      bool  __findOneEdge(const TriangleCut& K,
347                          const size_t & number,
348                          TinyVector<2,const Vertex*>& P);
349      
350      void __findInEdges(const TriangleCut& K,
351                         TinyVector<2,const Vertex*>& P,
352                         TinyVector<2,const Triangle*>& T,
353                         size_t& iInter,
354                         size_t& stop,
355                         size_t& in,
356                         size_t& in1,
357                         TinyVector<3,size_t> place,
358                         std::vector<const Vertex*>& points);
359          
360      void __findVertex(const TriangleCut& K,
361                        TinyVector<2,const Vertex*>& P,
362                        TinyVector<2,const Triangle*>& T,
363                        size_t& iCut,
364                        size_t& iInter,
365                        size_t& stop,
366                        size_t& in,
367                        size_t& in1,
368                        TinyVector<3,size_t> place,
369                        std::vector<const Vertex*>& points);
370    
371      void __findTriangle(TinyVector<2,const Vertex*>& P,
372                          TinyVector<2,const Triangle*>& T,
373                          size_t& stop);
374        
375    void __createGeneral(const TriangleCut& K,const Cell* cell,    void __createGeneral(const TriangleCut& K,const Cell* cell,
376                         std::vector<Triangle>& triangleListIntersectionNew);                         std::vector<Triangle>& triangleListIntersectionNew);
     
   void  __createNotIn(const TriangleCut& K,const Cell* cell,  
                       std::vector<Triangle>& triangleListIntersectionNew);  
   
   void  __createNotIn2(const TriangleCut& K,const Cell* cell,  
                        std::vector<Triangle>& triangleListIntersectionNew);  
377    
378    void __3pointsintersection(const TriangleCut& K,const Cell* cell,    void __createTriangle(const TriangleCut& K,const Cell* cell,
379                               std::vector<Triangle>& triangleListIntersectionNew);                          std::vector<Triangle>& triangleListIntersectionNew,
380                            std::vector<const Vertex*>& points,
381                            std::vector<Vertex>& pointtrian);
382    
383    template <typename booleanTest>    template <typename booleanTest>
384    void __createTrianglesIntersection(const size_t numobject,    void __createTrianglesIntersection(const size_t numobject,
# Line 342  private: Line 404  private:
404    
405    bool __isDegenerate(const Triangle& T);    bool __isDegenerate(const Triangle& T);
406    
   void __casePinterNotInList(const TriangleCut& K,const Cell* cell, std::vector<Triangle>& triangleListIntersectionNew);  
   
   TinyVector<3,const Vertex*> __determinatePinterNotInList(const TriangleCut& K, TinyVector<3,size_t>);  
   
407    template <typename booleanTest>    template <typename booleanTest>
408    void __calculateIntersection(ObjectToTreat& O0,    void __calculateIntersection(ObjectToTreat& O0,
409                                 const ObjectToTreat& O1,                                 const ObjectToTreat& O1,
# Line 364  private: Line 422  private:
422    ReferenceCounting<Vector<Triangle> >    ReferenceCounting<Vector<Triangle> >
423    __marchingTetrahedra(Structured3DMesh& SMesh,    __marchingTetrahedra(Structured3DMesh& SMesh,
424                         ObjectToTreat& O,                         ObjectToTreat& O,
425  //                     const Shape& S,                         //                      const Shape& S,
426  //                     const size_t& ref,                         //                      const size_t& ref,
427                         size_t nobject);                         size_t nobject);
428            
429    /// @todo FUNCTION TO REMOVE    /// @todo FUNCTION TO REMOVE
# Line 475  public : Line 533  public :
533                const Vertex& v2,                const Vertex& v2,
534                size_t num,                size_t num,
535                size_t num2)                size_t num2)
536        : Triangle(v0,v1,v2)      : Triangle(v0,v1,v2)
537    {    {
538      numobject=num;      numobject=num;
539      otherobject=num2;      otherobject=num2;
# Line 545  __splitEdge(const Edge& e, Line 603  __splitEdge(const Edge& e,
603    TinyVector<3> vo = e(out);    TinyVector<3> vo = e(out);
604    TinyVector<3> vs = 0.5*(vi+vo);    TinyVector<3> vs = 0.5*(vi+vo);
605    
606    for (int l = 0; l<20; l++) { // 3% error after 5 iterations.    for (int l = 0; l<20/*20*/; l++) { // 3% error after 5 iterations.
607      if(S.inside(vs)) {      if(S.inside(vs)) {
608        vi = vs;        vi = vs;
609        vs = 0.5*(vo+vi);        vs = 0.5*(vo+vi);
# Line 629  __createCutEdges(const Scene& scene, Line 687  __createCutEdges(const Scene& scene,
687                   const std::vector<TetrahedronByEdges>::iterator& currentT)                   const std::vector<TetrahedronByEdges>::iterator& currentT)
688  {  {
689    
690      //size_t nobjectToTreat = otherReferences.size();    //size_t nobjectToTreat = otherReferences.size();
691    for(size_t ie=0; ie<6; ++ie) {    for(size_t ie=0; ie<6; ++ie) {
692      const Edge& e = currentT->edge(ie);      const Edge& e = currentT->edge(ie);
693      int num=0;      int num=0;
# Line 731  void SurfaceMeshGenerator::Internals::__ Line 789  void SurfaceMeshGenerator::Internals::__
789            
790    }    }
791    // }    // }
792  //listOfVertexMesh.clear();    //listOfVertexMesh.clear();
793  }  }
794    
795  void SurfaceMeshGenerator::Internals::__dataStructureConvertion(ObjectToTreat& O,  void SurfaceMeshGenerator::Internals::__dataStructureConvertion(ObjectToTreat& O,
# Line 792  __calculatePointsIntersection(const Obje Line 850  __calculatePointsIntersection(const Obje
850                                const std::set<const Cell*>& toTreatHexahedra,                                const std::set<const Cell*>& toTreatHexahedra,
851                                PIntersection& listVertexIntersectionNew)                                PIntersection& listVertexIntersectionNew)
852  {  {
853      // calculate points intersection between object O1 and O2    // calculate points intersection between object O1 and O2
854    const MapCellTriangle& hexaToTriangle0=*(O1.hexalist);    const MapCellTriangle& hexaToTriangle0=*(O1.hexalist);
855    const MapCellTriangle& hexaToTriangle1=*(O2.hexalist);    const MapCellTriangle& hexaToTriangle1=*(O2.hexalist);
856            
# Line 821  bool SurfaceMeshGenerator::Internals::__ Line 879  bool SurfaceMeshGenerator::Internals::__
879    return (&(T(0))==&(T(1)) or &(T(0))==&(T(2)) or &(T(1))==&(T(2)));    return (&(T(0))==&(T(1)) or &(T(0))==&(T(2)) or &(T(1))==&(T(2)));
880  }  }
881    
 TinyVector<3,const Vertex*>  
 SurfaceMeshGenerator::Internals::__determinatePinterNotInList(const TriangleCut& K,  
                                                               TinyVector<3,size_t> ncase)  
 {  
   TinyVector<3,const Vertex*> triangle;  
   triangle[0]=(K.pointsIntersection[ncase[0]][0]);  
   if(K.pointsIntersection[ncase[1]][0]!=triangle[0]) {  
     triangle[1]=(K.pointsIntersection[ncase[1]][0]);  
   } else {  
     triangle[1]=(K.pointsIntersection[ncase[1]][1]);  
   }  
   
   if(K.pointsIntersection[ncase[2]][0]!=triangle[0] and K.pointsIntersection[ncase[2]][0]!=triangle[1]) {  
     triangle[2]=(K.pointsIntersection[ncase[2]][0]);  
   } else {  
     triangle[2]=(K.pointsIntersection[ncase[2]][1]);  
   }  
   
   return triangle;  
 }  
   
 void SurfaceMeshGenerator::Internals::__casePinterNotInList(const TriangleCut& K,  
                                                             const Cell* cell,  
                                                             std::vector<Triangle>& triangleListIntersectionNew)  
 {  
     //build triangles if K is outside but the 3 egdes are cut  
   TinyVector<3,const Vertex*> triangle;  
   TinyVector<3,size_t> ncase;  
   if(K.pointsIntersection[0][0]==K.pointsIntersection[0][1]) {  
     ncase[0]=0;  
     ncase[1]=1;  
     ncase[2]=2;  
     triangle=__determinatePinterNotInList(K,ncase);  
   } else {  
     if(K.pointsIntersection[1][0]==K.pointsIntersection[1][1]) {  
       ncase[0]=1;  
       ncase[1]=0;  
       ncase[2]=2;  
       triangle=__determinatePinterNotInList(K,ncase);  
     } else {  
       if(K.pointsIntersection[2][0]==K.pointsIntersection[2][1]) {  
         ncase[0]=2;  
         ncase[1]=1;  
         ncase[2]=0;  
         triangle=__determinatePinterNotInList(K,ncase);  
       } else {  
         triangle[0]=K.pointsIntersection[0][0];  
         triangle[1]=K.pointsIntersection[1][0];  
         triangle[2]=K.pointsIntersection[2][0];  
       }  
     }  
   }  
 //#warning Reference not set.  
   Triangle T(*triangle[0], *triangle[1], *triangle[2],K.reference());  
   __addTriangle(T, cell, triangleListIntersectionNew);  
 }  
   
882  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
883  __InTriangle(const Vertex& P,  __InTriangle(const Vertex& P,
884               TinyVector<3,const Vertex* > & V,               TinyVector<3,const Vertex* > & V,
# Line 911  __InTriangle(const Vertex& P, Line 912  __InTriangle(const Vertex& P,
912  TinyVector<4, real_t>  SurfaceMeshGenerator::Internals::  TinyVector<4, real_t>  SurfaceMeshGenerator::Internals::
913  __EquationPlan(const TinyVector<3,const Vertex *> & V){  __EquationPlan(const TinyVector<3,const Vertex *> & V){
914    
915      //equation of plan    //equation of plan
916    TinyVector<3> V01=(*V[1])-(*V[0]);    TinyVector<3> V01=(*V[1])-(*V[0]);
917    TinyVector<3> V12=(*V[2])-(*V[1]);    TinyVector<3> V12=(*V[2])-(*V[1]);
918    TinyVector<3> V02=(*V[2])-(*V[0]);    TinyVector<3> V02=(*V[2])-(*V[0]);
# Line 937  __DeterminatePoint(const size_t& ncase, Line 938  __DeterminatePoint(const size_t& ncase,
938                     IntersectionPoints& Pinter,                     IntersectionPoints& Pinter,
939                     PIntersection& listVertexIntersectionNew)                     PIntersection& listVertexIntersectionNew)
940  {  {
941      //calculate intersection between edge Vo0Vo1 and triangle V    //calculate intersection between edge Vo0Vo1 and triangle V
942      // if ncase<3 the edge is in the first object else the edge is in the second object    // if ncase<3 the edge is in the first object else the edge is in the second object
943      // P = point intersection    // P = point intersection
944    std::pair<const Vertex*,const Vertex*> pV(Vo0,Vo1);    std::pair<const Vertex*,const Vertex*> pV(Vo0,Vo1);
945    std::pair<const Vertex*,const Vertex*> pV2(Vo1,Vo0);    std::pair<const Vertex*,const Vertex*> pV2(Vo1,Vo0);
946    
947    real_t epsilon=1e-6;      // rq: a 10-6 on loupe des points!!!!! (testStep2.txt)
948      real_t epsilon=1e-10;
949    real_t pscalar,k;    real_t pscalar,k;
950    size_t num;    size_t num;
951    Vertex P;    Vertex P;
# Line 1022  __DeterminatePoint(const size_t& ncase, Line 1024  __DeterminatePoint(const size_t& ncase,
1024            Pinter.cutTriangle1()[ncase]=Pinter.size();            Pinter.cutTriangle1()[ncase]=Pinter.size();
1025          }          }
1026        } else {        } else {
1027            //il ne faut pas que le point soit sur une arete du tr 1          //il ne faut pas que le point soit sur une arete du tr 1
1028          if(inV01 > epsilon and inV02 > epsilon and inV12 > epsilon ) {          if(inV01 > epsilon and inV02 > epsilon and inV12 > epsilon ) {
1029            __InTriangle(P,V,num);            __InTriangle(P,V,num);
1030            if(num==3) {            if(num==3) {
# Line 1037  __DeterminatePoint(const size_t& ncase, Line 1039  __DeterminatePoint(const size_t& ncase,
1039                  Pinter.addPoint(Vo1);                  Pinter.addPoint(Vo1);
1040                } else                } else
1041                  {                  {
1042                  if(listVertexIntersectionNew.find(pV)==listVertexIntersectionNew.end()                    if(listVertexIntersectionNew.find(pV)==listVertexIntersectionNew.end()
1043                     and listVertexIntersectionNew.find(pV2)==listVertexIntersectionNew.end()) {                       and listVertexIntersectionNew.find(pV2)==listVertexIntersectionNew.end()) {
1044                    listVertexIntersectionNew[pV][0]=new Vertex (P);                      listVertexIntersectionNew[pV][0]=new Vertex (P);
1045                    listVertexIntersectionNew[pV][1]=new Vertex (P);                      listVertexIntersectionNew[pV][1]=new Vertex (P);
1046                    Pinter.addPoint(listVertexIntersectionNew[pV][0]);                      Pinter.addPoint(listVertexIntersectionNew[pV][0]);
1047                  } else {                    } else {
1048                    if(listVertexIntersectionNew.find(pV)!=listVertexIntersectionNew.end()) {                      if(listVertexIntersectionNew.find(pV)!=listVertexIntersectionNew.end()) {
1049                      if(Norm((*(*listVertexIntersectionNew.find(pV)).second[0])                        if(Norm((*(*listVertexIntersectionNew.find(pV)).second[0])
                             -P)<1e-6) {  
                       Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[0]);  
                     } else {  
                       if(Norm((*(*listVertexIntersectionNew.find(pV)).second[1])  
1050                                -P)<1e-6) {                                -P)<1e-6) {
1051                          Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[1]);                          Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[0]);
1052                        } else {                        } else {
1053                          listVertexIntersectionNew[pV][1]=new Vertex (P);                          if(Norm((*(*listVertexIntersectionNew.find(pV)).second[1])
1054                          Pinter.addPoint(listVertexIntersectionNew[pV][1]);                                  -P)<1e-6) {
1055                              Pinter.addPoint((*listVertexIntersectionNew.find(pV)).second[1]);
1056                            } else {
1057                              listVertexIntersectionNew[pV][1]=new Vertex (P);
1058                              Pinter.addPoint(listVertexIntersectionNew[pV][1]);
1059                            }
1060                        }                        }
                     }  
                   } else {  
                     if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[0])  
                             -P)<1e-6) {  
                       Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[0]);  
1061                      } else {                      } else {
1062                        if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[1])                        if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[0])
1063                                -P)<1e-6) {                                -P)<1e-6) {
1064                          Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[1]);                          Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[0]);
1065                        } else {                        } else {
1066                          listVertexIntersectionNew[pV2][1]=new Vertex (P);                          if(Norm((*(*listVertexIntersectionNew.find(pV2)).second[1])
1067                          Pinter.addPoint(listVertexIntersectionNew[pV2][1]);                                  -P)<1e-6) {
1068                              Pinter.addPoint((*listVertexIntersectionNew.find(pV2)).second[1]);
1069                            } else {
1070                              listVertexIntersectionNew[pV2][1]=new Vertex (P);
1071                              Pinter.addPoint(listVertexIntersectionNew[pV2][1]);
1072                            }
1073                        }                        }
1074                      }                      }
1075                    }                    }
1076                  }                  }
1077                }              }
                     }  
1078              // Pinter.addPoint(listVertexIntersectionNew[listVertexIntersectionNew.size()-1]);              // Pinter.addPoint(listVertexIntersectionNew[listVertexIntersectionNew.size()-1]);
1079            }            }
1080          }          }
# Line 1087  __CalculPoint(std::list<Triangle*>::iter Line 1089  __CalculPoint(std::list<Triangle*>::iter
1089                IntersectionPoints& Pinter,                IntersectionPoints& Pinter,
1090                PIntersection& listVertexIntersectionNew)                PIntersection& listVertexIntersectionNew)
1091  {  {
1092      // calculate points intersection between triangle "currentTriangle" and "currentTriangle1"    // calculate points intersection between triangle "currentTriangle" and "currentTriangle1"
1093      // this points are in "Pinter"    // this points are in "Pinter"
1094    
1095    // NECESSITE DE V012 et Vo012 ?    // NECESSITE DE V012 et Vo012 ?
1096    
# Line 1152  __createPointIntersection(std::list<Tria Line 1154  __createPointIntersection(std::list<Tria
1154    // CALCUL DU POINT D'INTERSECTION    // CALCUL DU POINT D'INTERSECTION
1155    __CalculPoint(currentTriangle, currentTriangle1, Pinter, listVertexIntersectionNew);    __CalculPoint(currentTriangle, currentTriangle1, Pinter, listVertexIntersectionNew);
1156    
1157    const TinyVector<3,size_t >& pintercut=Pinter.cutTriangle1();    //const TinyVector<3,size_t >& pintercut=Pinter.cutTriangle1();
1158    
1159    if(Pinter.size()!=0) {    if(Pinter.size()!=0) {
1160      //ffout(4)<<"rajout2 \n";      //ffout(4)<<"rajout2 \n";
# Line 1163  __createPointIntersection(std::list<Tria Line 1165  __createPointIntersection(std::list<Tria
1165      TinyVector<3,size_t >& cutTr2= pairTriangleToIntersectionPoints[pTriangle].cutTriangle2();      TinyVector<3,size_t >& cutTr2= pairTriangleToIntersectionPoints[pTriangle].cutTriangle2();
1166    
1167      for(size_t i=0 ; i<Pinter.size() ; ++i) {      for(size_t i=0 ; i<Pinter.size() ; ++i) {
1168        real_t epsilon=1e-6;        real_t epsilon=1e-7;
1169        const Triangle& To=(*(*currentTriangle));        const Triangle& To=(*(*currentTriangle));
1170        const Vertex& Vo0 = To(0);        const Vertex& Vo0 = To(0);
1171        const Vertex& Vo1 = To(1);        const Vertex& Vo1 = To(1);
# Line 1179  __createPointIntersection(std::list<Tria Line 1181  __createPointIntersection(std::list<Tria
1181        real_t norm1=Norm(P-Vo1);        real_t norm1=Norm(P-Vo1);
1182        real_t norm2=Norm(P-Vo2);        real_t norm2=Norm(P-Vo2);
1183        //for filter points        //for filter points
1184        real_t eps=1e-6;        //real_t eps=1e-6;
1185        if(/*norm0>eps and norm1>eps and */std::abs((norm0+norm1)-Vo01)<epsilon) {        if(/*norm0>eps and norm1>eps and */std::abs((norm0+norm1)-Vo01)<epsilon) {
1186          cutTr2[0]=i+1;          cutTr2[0]=i+1;
1187        }        }
# Line 1217  __addTriangle(const Triangle& t, Line 1219  __addTriangle(const Triangle& t,
1219    } else {    } else {
1220      fferr(3) << "\t\tTriangle dégénéré :"      fferr(3) << "\t\tTriangle dégénéré :"
1221               << &T(0) << ':' << &T(1) << ':' << &T(2) << '\n';               << &T(0) << ':' << &T(1) << ':' << &T(2) << '\n';
1222        fferr(3)      fferr(3)
1223            << T(0) << ':' << T(1) << ':' << T(2) << '\n';        << T(0) << ':' << T(1) << ':' << T(2) << '\n';
1224     // triangleListIntersectionNew.push_back(T);      // triangleListIntersectionNew.push_back(T);
1225    }    }
1226  }  }
1227    
# Line 1229  __determinateCase(const int& ncase, Line 1231  __determinateCase(const int& ncase,
1231                    TinyVector<3,int>& edgesCut,                    TinyVector<3,int>& edgesCut,
1232                    TinyVector<3,int>& vertexNum)                    TinyVector<3,int>& vertexNum)
1233  {  {
1234      //determinate which point is inside    //determinate which point is inside
1235    const TinyVector<3,TinyVector<2,int> >& edgesCut2=K.edgecut;    const TinyVector<3,TinyVector<2,int> >& edgesCut2=K.edgecut;
1236    switch(ncase) {        switch(ncase) {    
1237    case 0: {    case 0: {
# Line 1269  __createLocalListIntersection(const size Line 1271  __createLocalListIntersection(const size
1271                                TriangleCut & K,                                TriangleCut & K,
1272                                const Cell* cell)                                const Cell* cell)
1273  {  {
1274      //build the vectors K.pointsIntersection, K.pointsin with points intersection include in K    //build the vectors K.pointsIntersection, K.pointsin with points intersection include in K
1275      //put K.edgecut[i]=1 if the edge i is cut    //put K.edgecut[i]=1 if the edge i is cut
1276      //build vertexInTriangles, triangleWithVertex (for build mesh by front)    //build vertexInTriangles, triangleWithVertex (for build mesh by front)
1277      //build vertexVectTriangles
1278    TinyVector<3,TinyVector<2,int> >& edgesCut=K.edgecut;    TinyVector<3,TinyVector<2,int> >& edgesCut=K.edgecut;
1279    TinyVector<3,TinyVector<2,const Vertex*> > & pointsIntersection=K.pointsIntersection;    TinyVector<3,TinyVector<2,const Vertex*> > & pointsIntersection=K.pointsIntersection;
1280    std::map<const Vertex* ,const Vertex*>& pointsInterior=K.pointsin;    std::map<const Vertex* ,const Vertex*>& pointsInterior=K.pointsin;
# Line 1324  __createLocalListIntersection(const size Line 1327  __createLocalListIntersection(const size
1327              const Vertex*& V=pairTriangleToIntersectionPoints[pairTriangle][b];              const Vertex*& V=pairTriangleToIntersectionPoints[pairTriangle][b];
1328              const Triangle T=(*(*jt0));              const Triangle T=(*(*jt0));
1329              // si on a le point est confondu avec un sommet on n'a pas besoin du triangle pour la connectivite              // si on a le point est confondu avec un sommet on n'a pas besoin du triangle pour la connectivite
1330              if(!(V==&(*(*jt0))(0) or V==&(*(*jt0))(1) or V==&(*(*jt0))(2))              if(!(V==&(*(*jt0))(0) or V==&(*(*jt0))(1) or V==&(*(*jt0))(2)) and
1331                  /*or b==0*/) {          Norm(*V-(*(*jt0))(0))>1e-6 and Norm(*V-(*(*jt0))(1))>1e-6 and Norm(*V-(*(*jt0))(2))>1e-6
1332                   /*or b==0*/) {
1333                  vertexVectTriangles[V].push_back(*jt0);
1334                if(vertexInTriangles.find(V)==vertexInTriangles.end()) {                if(vertexInTriangles.find(V)==vertexInTriangles.end()) {
1335                  TinyVector<2,const Triangle*> triangles;                  TinyVector<2,const Triangle*> triangles;
1336                  triangles[0]=*jt0;                  triangles[0]=*jt0;
# Line 1342  __createLocalListIntersection(const size Line 1347  __createLocalListIntersection(const size
1347                } else {                } else {
1348                  (*triangleWithVertex.find(*jt0)).second[1]=V;                  (*triangleWithVertex.find(*jt0)).second[1]=V;
1349                }                }
1350                } else {
1351                  //le point est confondu avec un sommet du tr mais dans certains cas on en a besoin qd meme!!
1352                  vertexVectTriangles[V].push_back(*jt0);
1353              }              }
1354            }            }
1355            } else {
1356                if(pairTriangleToIntersectionPoints[pairTriangle].size()!=0) {
1357                    ffout(4)<<"on a un point non pris... triangle plat??\n";
1358                    //ffout(4)<<(*(*jt0))(0)<<"\n";
1359                    //ffout(4)<<(*(*jt0))(1)<<"\n";
1360                    //ffout(4)<<(*(*jt0))(2)<<"\n";
1361                    ffout(4)<<*pairTriangleToIntersectionPoints[pairTriangle][0]<<"\n";
1362                    vertexVectTriangles[pairTriangleToIntersectionPoints[pairTriangle][0]].push_back(*jt0);
1363                }
1364          }          }
1365    
1366          for (size_t n=0; n < cutTriangle.size(); ++n) {          for (size_t n=0; n < cutTriangle.size(); ++n) {
# Line 1396  __createLocalListIntersection(const size Line 1413  __createLocalListIntersection(const size
1413              ffout(4) << K(0) << "\n";              ffout(4) << K(0) << "\n";
1414              ffout(4) << K(1) << "\n";              ffout(4) << K(1) << "\n";
1415              ffout(4) << K(2) << "\n";              ffout(4) << K(2) << "\n";
1416              ffout(4) << cutTriangle[0] << " " << cutTriangle[1] << " " << cutTriangle[2] << "\n";              ffout(4) << "cutTriangle "<<cutTriangle[0] << " " << cutTriangle[1] << " " << cutTriangle[2] << "\n";
1417              ffout(4) << edgesCut[0] << " " << edgesCut[1] << " " << edgesCut[2] << "\n";              ffout(4) << "edgesCut "<<edgesCut[0] << " " << edgesCut[1] << " " << edgesCut[2] << "\n";
1418              edgesCut[0]=cutTriangle[0];              //edgesCut[0]=cutTriangle[0];
1419              edgesCut[1]=cutTriangle[1];              //edgesCut[1]=cutTriangle[1];
1420              edgesCut[2]=cutTriangle[2];              //edgesCut[2]=cutTriangle[2];
1421                            
1422            }            }
1423            }            }
# Line 1429  __determinateEdgeCut(const size_t& numob Line 1446  __determinateEdgeCut(const size_t& numob
1446    edgesCut=0;    edgesCut=0;
1447    size_t objectNumber=0;    size_t objectNumber=0;
1448    if(numobject==0)    if(numobject==0)
1449        objectNumber=1;      objectNumber=1;
1450    const MapCellTriangle& hexaToTriangle1=*(O2.hexalist);    const MapCellTriangle& hexaToTriangle1=*(O2.hexalist);
1451    
1452    MapCellTriangle::const_iterator it1=hexaToTriangle1.find(cell);    MapCellTriangle::const_iterator it1=hexaToTriangle1.find(cell);
# Line 1494  __putRefByFront(const size_t numobject, Line 1511  __putRefByFront(const size_t numobject,
1511      const Cell* cell=&((*T).mother());      const Cell* cell=&((*T).mother());
1512      if(toTreatHexahedra.find(cell)!=toTreatHexahedra.end()      if(toTreatHexahedra.find(cell)!=toTreatHexahedra.end()
1513         and ((*edgesRefVertex[numobject])[&(*T)(0)]==2         and ((*edgesRefVertex[numobject])[&(*T)(0)]==2
1514        or (*edgesRefVertex[numobject])[&(*T)(1)]==2              or (*edgesRefVertex[numobject])[&(*T)(1)]==2
1515        or (*edgesRefVertex[numobject])[&(*T)(2)]==2)) {              or (*edgesRefVertex[numobject])[&(*T)(2)]==2)) {
1516        front.push_front(T);        front.push_front(T);
1517        frontSize++;        frontSize++;
1518        const Triangle*& TT=(*listOfTriangleMeshFront.find(T)).second;        const Triangle*& TT=(*listOfTriangleMeshFront.find(T)).second;
# Line 1508  __putRefByFront(const size_t numobject, Line 1525  __putRefByFront(const size_t numobject,
1525    size_t iter=0;    size_t iter=0;
1526    ffout(4)<<"taille initial du front ="<< frontSize <<"\n";    ffout(4)<<"taille initial du front ="<< frontSize <<"\n";
1527    //if(front.size()==0)    //if(front.size()==0)
1528        ffout(4)<<"ref de surfmesh = "<<surfmesh.cell(0).reference()<<"\n";    ffout(4)<<"ref de surfmesh = "<<surfmesh.cell(0).reference()<<"\n";
1529    //put reference by front    //put reference by front
1530    size_t erase=0;    size_t erase=0;
1531    while(frontSize !=0 and iter<nbcell) {    while(frontSize !=0 and iter<nbcell) {
# Line 1538  __putRefByFront(const size_t numobject, Line 1555  __putRefByFront(const size_t numobject,
1555        size_t& refo2=(*edgesRefVertex[numobject])[&(*T)(2)];        size_t& refo2=(*edgesRefVertex[numobject])[&(*T)(2)];
1556                
1557        if(((refo0==2 or refo1==2 or refo2==2) and  refo0+refo1+refo2<6 and edgeo0+edgeo1+edgeo2==0)        if(((refo0==2 or refo1==2 or refo2==2) and  refo0+refo1+refo2<6 and edgeo0+edgeo1+edgeo2==0)
1558          /* or (refo0+refo1+refo2!=6 and edgeo0+edgeo1+edgeo2==0 and refo0+refo1+refo2>0)*/           /* or (refo0+refo1+refo2!=6 and edgeo0+edgeo1+edgeo2==0 and refo0+refo1+refo2>0)*/
1559           /*or (edgeo0[0]==1 and refo0+refo1!=2 and refo0+refo1+refo2!=0)           /*or (edgeo0[0]==1 and refo0+refo1!=2 and refo0+refo1+refo2!=0)
1560           or (edgeo1[0]==1 and refo0+refo2!=2 and refo0+refo1+refo2!=0)             or (edgeo1[0]==1 and refo0+refo2!=2 and refo0+refo1+refo2!=0)
1561           or (edgeo2[0]==1 and refo2+refo1!=2 and refo0+refo1+refo2!=0)*/) {             or (edgeo2[0]==1 and refo2+refo1!=2 and refo0+refo1+refo2!=0)*/) {
1562            ffout(4)<<"chouette on garde "<<refo0+refo1+refo2<<"\n";          ffout(4)<<"chouette on garde "<<refo0+refo1+refo2<<"\n";
1563            refo0=2;          refo0=2;
1564            refo1=2;          refo1=2;
1565            refo2=2;          refo2=2;
1566            erase=1;          erase=1;
1567        } else {        } else {
1568            erase=1;          erase=1;
1569        }        }
1570        const TinyVector<3,const Triangle*>& V=connect(surfmesh.cellNumber(*T));        const TinyVector<3,const Triangle*>& V=connect(surfmesh.cellNumber(*T));
1571                
# Line 1576  __putRefByFront(const size_t numobject, Line 1593  __putRefByFront(const size_t numobject,
1593              TinyVector<2,int> edge2=Kneigh.edgecut[2];              TinyVector<2,int> edge2=Kneigh.edgecut[2];
1594              size_t sum=edge0[0]+edge1[0]+edge2[0];              size_t sum=edge0[0]+edge1[0]+edge2[0];
1595              if((ref0==2 or ref1==2 or ref2==2)  and              if((ref0==2 or ref1==2 or ref2==2)  and
1596          (ref0!=1 and ref1!=1 and ref2!=1) and sum==0) {                 (ref0!=1 and ref1!=1 and ref2!=1) and sum==0) {
1597                  ffout(4)<<"ref 2 partout\n";                ffout(4)<<"ref 2 partout\n";
1598                  test=1;                test=1;
1599                  ref0=2;                ref0=2;
1600                  ref1=2;                ref1=2;
1601                  ref2=2;                ref2=2;
1602              } else {              } else {
1603                  if((ref1==2 or ref0==2) and edge0==1){                if((ref1==2 or ref0==2) and edge0==1){
1604                      //ffout(4)<<"ref 1 2 0 partout\n";                  //ffout(4)<<"ref 1 2 0 partout\n";
1605                      //ffout(4)<<ref0<<" "<<ref1<<"\n";                  //ffout(4)<<ref0<<" "<<ref1<<"\n";
1606                      if(ref1+ref0!=3) {                  if(ref1+ref0!=3) {
1607                          test=1;                    test=1;
                     }  
                     if(ref0==2) {  
                         ref1=1;  
                     } else {  
                         ref0=1;  
                     }  
1608                  }                  }
1609                  if((ref2==2 or ref0==2) and edge1==1){                  if(ref0==2) {
1610                      //ffout(0)<<"ref 1 2 1 partout\n";                    ref1=1;
1611                      //ffout(4)<<ref0<<" "<<ref2<<"\n";                  } else {
1612                      if(ref2+ref0!=3) {                    ref0=1;
                         test=1;  
                     }  
                     if(ref0==2) {  
                         ref2=1;  
                     } else {  
                         ref0=1;  
                     }  
1613                  }                  }
1614                  if((ref2==2 or ref1==2) and edge2==1){                }
1615                      //ffout(0)<<"ref 1 2 2 partout\n";                if((ref2==2 or ref0==2) and edge1==1){
1616                      //ffout(4)<<ref1<<" "<<ref2<<"\n";                  //ffout(0)<<"ref 1 2 1 partout\n";
1617                      if(ref1+ref2!=3) {                  //ffout(4)<<ref0<<" "<<ref2<<"\n";
1618                          test=1;                  if(ref2+ref0!=3) {
1619                      }                    test=1;
1620                      if(ref2==2) {                  }
1621                          ref1=1;                  if(ref0==2) {
1622                      } else {                    ref2=1;
1623                          ref2=1;                  } else {
1624                      }                    ref0=1;
1625                    }
1626                  }
1627                  if((ref2==2 or ref1==2) and edge2==1){
1628                    //ffout(0)<<"ref 1 2 2 partout\n";
1629                    //ffout(4)<<ref1<<" "<<ref2<<"\n";
1630                    if(ref1+ref2!=3) {
1631                      test=1;
1632                    }
1633                    if(ref2==2) {
1634                      ref1=1;
1635                    } else {
1636                      ref2=1;
1637                  }                  }
1638                  }
1639              }                    }      
1640              if(test==1/*ref0+ref1+ref2==6 *//*or edge0[0]+edge1[0]+edge2[0]==1*/) {              if(test==1/*ref0+ref1+ref2==6 *//*or edge0[0]+edge1[0]+edge2[0]==1*/) {
1641                front.push_front(triangle);                front.push_front(triangle);
# Line 1630  __putRefByFront(const size_t numobject, Line 1647  __putRefByFront(const size_t numobject,
1647        }        }
1648      }      }
1649      if(erase==1 and frontSize==1) {      if(erase==1 and frontSize==1) {
1650          er=front.erase(er);        er=front.erase(er);
1651          frontSize--;        frontSize--;
1652      }      }
1653      ++iter;      ++iter;
1654      ffout(4)<<"taille du front = "<<front.size()<<" "<<iter<<" "<<nbcell<<"\n";      ffout(4)<<"taille du front = "<<front.size()<<" "<<iter<<" "<<nbcell<<"\n";
# Line 1646  __putEdgeRef(const Triangle*& Tmesh, Line 1663  __putEdgeRef(const Triangle*& Tmesh,
1663               const ObjectToTreat& O1,               const ObjectToTreat& O1,
1664               const ObjectToTreat& O2 )               const ObjectToTreat& O2 )
1665  {  {
1666      //put edgesRefVertex=2 if the vertex is in the domain    //put edgesRefVertex=2 if the vertex is in the domain
1667      //ffout(0)<<"j'initialise le front\n";    //ffout(0)<<"j'initialise le front\n";
1668    
1669      TinyVector<3,TinyVector<2,int> > edgecut;    TinyVector<3,TinyVector<2,int> > edgecut;
1670      const Cell* cell=&((*T).mother());    const Cell* cell=&((*T).mother());
1671      edgecut =__determinateEdgeCut(numobject,O1,O2,T,cell);    edgecut =__determinateEdgeCut(numobject,O1,O2,T,cell);
1672      size_t& ref0=(*edgesRefVertex[numobject])[&(*Tmesh)(0)];    size_t& ref0=(*edgesRefVertex[numobject])[&(*Tmesh)(0)];
1673      size_t& ref1=(*edgesRefVertex[numobject])[&(*Tmesh)(1)];    size_t& ref1=(*edgesRefVertex[numobject])[&(*Tmesh)(1)];
1674      size_t& ref2=(*edgesRefVertex[numobject])[&(*Tmesh)(2)];    size_t& ref2=(*edgesRefVertex[numobject])[&(*Tmesh)(2)];
1675            
1676      if(ref0+ref1+ref2!=6) {    if(ref0+ref1+ref2!=6) {
1677          TinyVector<2,int> edge0=edgecut[0];      TinyVector<2,int> edge0=edgecut[0];
1678          TinyVector<2,int> edge1=edgecut[1];      TinyVector<2,int> edge1=edgecut[1];
1679          TinyVector<2,int> edge2=edgecut[2];      TinyVector<2,int> edge2=edgecut[2];
1680          size_t sum=edge0[0]+edge1[0]+edge2[0];      size_t sum=edge0[0]+edge1[0]+edge2[0];
1681          if((ref0==2 or ref1==2 or ref2==2)  and sum==0) {      if((ref0==2 or ref1==2 or ref2==2)  and sum==0) {
1682              ref0=2;        ref0=2;
1683              ref1=2;        ref1=2;
1684              ref2=2;        ref2=2;
1685        } else {
1686          if((ref1==2 or ref0==2) and edge0==1){
1687            if(ref0==2) {
1688              ref1=1;
1689          } else {          } else {
1690              if((ref1==2 or ref0==2) and edge0==1){            ref0=1;
1691                  if(ref0==2) {          }
1692                      ref1=1;        }
1693                  } else {        if((ref2==2 or ref0==2) and edge1==1){
1694                      ref0=1;          if(ref0==2) {
1695                  }            ref2=1;
1696              }          } else {
1697              if((ref2==2 or ref0==2) and edge1==1){            ref0=1;
1698                  if(ref0==2) {          }
1699                      ref2=1;        }
1700                  } else {        if((ref2==2 or ref1==2) and edge2==1){
1701                      ref0=1;          if(ref2==2) {
1702                  }            ref1=1;
1703              }          } else {
1704              if((ref2==2 or ref1==2) and edge2==1){            ref2=1;
                 if(ref2==2) {  
                     ref1=1;  
                 } else {  
                     ref2=1;  
                 }  
             }  
1705          }          }
1706          }
1707      }      }
1708   }    }
1709    }
1710    
1711  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1712  __createGeneral(const TriangleCut& K,const Cell* cell,  __createCase(const TriangleCut& K ,
1713                  std::vector<Triangle>& triangleListIntersectionNew)               size_t& in, size_t& in1,
1714                 TinyVector<3,size_t>& place)
1715  {  {
1716  // pour l'instant on suppose edge cut une seule fois    size_t numberIn=K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2];
1717      size_t iCut=0, iInter=0, stop=0,in=0,in1=0;      
1718      size_t nbInter=triangleWithVertex.size();    switch (numberIn) {
1719      size_t nbCut=K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0];    case 0: {
1720      TinyVector<2,const Vertex*> P;      if(K.edgecut[0][0]==1) {
1721      TinyVector<2,const Triangle*> T;        place[0]=0;
1722      TinyVector<3,size_t> place;        if(K.edgecut[1][0]==1) {
1723      size_t numberIn=K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2];          place[1]=1;
1724      switch (numberIn) {          place[2]=2;
1725        case 1: {        } else {
1726            if(K.isInTriangle[1]==1) {          place[1]=2;
1727                in=1;          place[2]=1;
               place[0]=2;  
               place[1]=0;  
               place[2]=1;  
           } else if(K.isInTriangle[2]==1) {  
               in=2;  
               place[0]=1;  
               place[1]=2;  
               place[2]=0;  
           } else {  
               place[0]=0;  
               place[1]=1;  
               place[2]=2;  
           }  
           break;  
1728        }        }
1729        case 2: {      } else {
1730            if(K.isInTriangle[1]==1 and K.isInTriangle[0]==1) {        if(K.edgecut[1][0]==1) {
1731                in=1;          place[0]=1;
1732                in1=0;          place[1]=2;
1733                place[0]=1;          place[2]=0;
1734                place[1]=2;        } else {
1735                place[2]=0;          place[0]=2;
1736            } else if(K.isInTriangle[2]==1 and K.isInTriangle[1]==1) {          place[1]=0;
1737                in=2;          place[2]=1;
1738                in1=1;        }
1739                place[0]=0;      }
1740        break;
1741      }
1742      case 1: {
1743        if(K.isInTriangle[1]==1) {
1744          in=1;
1745          in1=in;
1746          place[0]=2;
1747          place[1]=0;
1748          place[2]=1;
1749        } else if(K.isInTriangle[2]==1) {
1750          in=2;
1751          in1=in;
1752          place[0]=1;
1753          place[1]=2;
1754          place[2]=0;
1755        } else {
1756          in1=in;
1757          place[0]=0;
1758          place[1]=1;
1759          place[2]=2;
1760        }
1761        break;
1762      }
1763      case 2: {
1764        if(K.isInTriangle[1]==1 and K.isInTriangle[0]==1) {
1765          in=1;
1766          in1=0;
1767          place[0]=2;
1768          place[1]=1;
1769          place[2]=0;
1770        } else if(K.isInTriangle[2]==1 and K.isInTriangle[1]==1) {
1771          in=2;
1772          in1=1;
1773          place[0]=1;
1774          place[1]=0;
1775          place[2]=2;
1776        } else {
1777          in=0;
1778          in1=2;
1779          place[0]=0;
1780          place[1]=2;
1781          place[2]=1;
1782        }
1783        break;
1784      }
1785      case 3: {
1786          if(K.isInTriangle[0]==1) {
1787              place[0]=0;
1788              if(K.isInTriangle[1]==1) {
1789                place[1]=1;                place[1]=1;
1790                place[2]=2;                place[2]=2;
1791            } else {            } else {
1792                in=0;                place[1]=2;
               in1=2;  
               place[0]=2;  
               place[1]=0;  
1793                place[2]=1;                place[2]=1;
1794            }            }
1795            break;        } else if(K.isInTriangle[1]==1) {
1796        }              place[0]=1;
1797              place[1]=2;
1798              place[2]=0;
1799          } else {
1800              place[0]=2;
1801              place[1]=1;
1802              place[2]=0;
1803          }
1804          break;
1805      }
1806      }
1807        
1808    }
1809    
1810    void SurfaceMeshGenerator::Internals::
1811    __transformVertex(const TriangleCut& K,
1812                      std::vector<const Vertex*>& points,
1813                      std::vector<Vertex> &pointstrian)
1814    {
1815        size_t taille=points.size();
1816        for(unsigned ie=0;ie<taille;++ie) {
1817            const Vertex& P=*points[ie];
1818            TinyMatrix<3,3, real_t> A;
1819            for (size_t i=0; i<3; ++i) {
1820                for (size_t j=0; j<3; ++j) {
1821                    A(i,j)=(K(j))[i];
1822                }
1823            }
1824            TinyVector<3, real_t> lambda;
1825            gaussPivot(A,P,lambda);
1826            Vertex Pi;
1827            Pi[0]=lambda[0];
1828            Pi[1]=lambda[1];
1829            Pi[2]=0;
1830            pointstrian[ie]=Pi;
1831      }      }
1832      const Vertex* A=&K(in);      
1833      P[0]=K.pointsIntersection[place[0]][0];  }
1834      ++iCut;  void SurfaceMeshGenerator::Internals::
1835      if(vertexInTriangles.find(P[0])!=vertexInTriangles.end()) {  __addPoints(const Vertex*& Pcut,
1836          T[0]=(*vertexInTriangles.find(P[0])).second[0];              size_t& stopcut,
1837                std::vector<const Vertex*>& pointslist1,
1838                std::vector<const Vertex*>& pointslist2)
1839    {
1840        if(stopcut==0) {
1841            //ffout(4)<<"00\n";
1842            pointslist1.push_back(Pcut);
1843            stopcut=3;
1844      } else {      } else {
1845          stop=1;          //ffout(4)<<"10\n";
1846          ffout(4)<<"pbs 0 existe pas\n";          pointslist2.push_back(Pcut);
1847      }      }
1848      if(numberIn==2 and stop==0) {  }
         __addTriangle(Triangle(K(in1),*A,*P[0],K.reference()),cell,triangleListIntersectionNew);  
     }  
     //iInter=1;  
     while(iInter<nbInter and stop==0) {  
         //ffout(4)<<"createGeneral "<<iInter<<" "<<nbInter<<"\n";  
         //on cherche le voisin de Pi  
         if(triangleWithVertex.find(T[0])!=triangleWithVertex.end()) {  
             TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[0])).second;  
             if(VV[0]==P[0]) {  
                 if(VV[1]==P[0]) {  
                     //on est a priori sur une arete  
                     if(iCut<nbCut) {  
                         #warning CAS edge cut twice  
                         //ffout(4)<<"cas 1\n";  
                         P[1]=K.pointsIntersection[place[iCut]][0];  
                         __addTriangle(Triangle(*A,*P[0],*P[1],K.reference()),cell,triangleListIntersectionNew);  
                         P[0]=P[1];  
                         ++iInter;  
                         ++iCut;  
                     } else {  
                         stop=1;  
                         ffout(4)<<"on s'arrete en plein milieu\n";  
                     }  
                 } else {  
                     //ffout(4)<<"cas 2\n";  
                     P[1]=VV[1];  
                     __addTriangle(Triangle(*A,*P[0],*P[1],K.reference()),cell,triangleListIntersectionNew);  
                     P[0]=P[1];  
                     ++iInter;  
                 }  
             } else {  
                 //ffout(4)<<"cas 3\n";  
                 P[1]=VV[0];  
                 __addTriangle(Triangle(*A,*P[0],*P[1],K.reference()),cell,triangleListIntersectionNew);  
                 P[0]=P[1];  
                 ++iInter;  
                   
             }  
         }  
1849    
1850          //on cherche le triangle associe au nouveau Pi pour relancer  void SurfaceMeshGenerator::Internals::
1851          if(stop==0 and iInter<nbInter) {  __addList(const TriangleCut& K)
1852              TinyVector<2,const Triangle*> T2;  {
1853              if(vertexInTriangles.find(P[0])!=vertexInTriangles.end()) {      
1854                  T2[0]=(*vertexInTriangles.find(P[0])).second[0];      size_t numberInEdges=K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]
1855                  T2[1]=(*vertexInTriangles.find(P[0])).second[1];      +K.edgecut[0][1]+K.edgecut[1][1]+K.edgecut[2][1];
1856              } else {      size_t nbInter=vertexInTriangles.size();
1857                  ffout(4)<<"pbs tr existe pas\n";      ffout(4)<<"nb Inter avt = "<<nbInter<<"\n";
1858                  stop=1;      ffout(4)<<"nb points total "<<vertexVectTriangles.size()<<"\n";
1859              }      //premiere etage : verifier que vertexInTriangles contient bien tous les points
1860              if(T2[0]==T[0]) {      if(nbInter<vertexVectTriangles.size()) {
1861                  if(T2[1]==T[0]) {          //il manque des points
1862                      ++iCut;          for(std::map<const Vertex*,std::vector<const Triangle*> >::iterator it=vertexVectTriangles.begin()
1863                      if(iCut<nbCut) {              ; it!=vertexVectTriangles.end() ; ++it) {
1864                          #warning CAS edge cut twice              const Vertex* PP=(*it).first;
1865                          ffout(4)<<iCut<<" cas new\n";              //ffout(4)<<"-----\n";
1866                          if(P[0]==K.pointsIntersection[place[iCut]][0]) {              //ffout(4)<<*PP<<"\n";
1867                              ffout(4)<<"j'ai le point1\n";  
1868                              --iCut;              std::vector<const Triangle*> VT=(*it).second;
1869                          }              if(vertexInTriangles.find(PP)==vertexInTriangles.end()) {
1870                          if(P[0]==K.pointsIntersection[place[1]][0]) {                  //ffout(4)<<"je l'ajoute\n";
1871                              ffout(4)<<"j'ai le point2\n";                  //le point n'existe pas il faut le rajouter
1872                              --iCut;    //on cherche a quel triangle l'associer
1873                          }                  if(nbInter==0) {
1874                          if(P[0]==K.pointsIntersection[place[0]][0]) {                      TinyVector<2,const Triangle*> triangles;
1875                              ffout(4)<<"j'ai le point3\n";                      triangles[0]=VT[0];
1876                              --iCut;                      triangles[1]=VT[0];
1877                        vertexInTriangles[PP]=triangles;
1878                    } else {
1879                        size_t nn=0,stop2=0;
1880                        size_t taille=VT.size();
1881                        //test si le point est confondu avec un point in
1882                        if(K.isInTriangle[0]==1 and Norm(K(0)-*PP)<1e-6) {
1883                            stop2=1;
1884                        }
1885                        if(K.isInTriangle[1]==1 and Norm(K(1)-*PP)<1e-6) {
1886                            stop2=1;
1887                        }
1888                        if(K.isInTriangle[2]==1 and Norm(K(2)-*PP)<1e-6) {
1889                            stop2=1;
1890                        }
1891                        while(nn<taille and stop2==0) {
1892                            const Triangle* TT=VT[nn];
1893                            if(triangleWithVertex.find(TT)!=triangleWithVertex.end()) {
1894                                (*triangleWithVertex.find(TT)).second[1]=PP;
1895                                TinyVector<2,const Triangle*> triangles;
1896                                triangles[0]=TT;
1897                                triangles[1]=TT;
1898                                vertexInTriangles[PP]=triangles;
1899                                stop2=1;
1900                          }                          }
1901                                                    ++nn;
1902                                                }
1903                          //P[1]=K.pointsIntersection[place[iCut]][0];                      if(nn==taille and stop2==0) {
1904                          //__addTriangle(Triangle(*A,*P[0],*P[1],K.reference()),cell,triangleListIntersectionNew);                          ffout(4)<<"bizarreeeeee on a pas trouve le tr\n";
                         //P[0]=P[1];  
                         //++iInter;  
                         //++iCut;  
                     } else {  
                         stop=1;  
                         ffout(4)<<"nbCut icut "<<iCut<<" "<<nbCut<<"\n";  
                         ffout(4)<<"pbs triangle\n";  
1905                      }                      }
                       
                 } else {  
                     T[0]=T2[1];  
1906                  }                  }
1907              } else {              } else {
1908                  T[0]=T2[0];  
1909              }              }
1910          }          }
1911      }      }
       
1912  }  }
1913    
1914  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
1915  __createNotIn(const TriangleCut& K,const Cell* cell,  __createList(const TriangleCut& K,
1916                std::vector<Triangle>& triangleListIntersectionNew)               std::vector<const Vertex*>& points,
1917  {               size_t& in,
1918      //create triangles if only one edge is cut               size_t& in1,
1919      // more than 3 points intersection?               TinyVector<3,size_t> place)
1920    size_t place0=0;  {
1921    if(K.edgecut[1][0]==1 and K.edgecut[1][1]==1) {    TinyVector<2,const Vertex*> P;
1922      place0=1;    TinyVector<2,const Triangle*> T;
1923    } else if(K.edgecut[2][0]==1 and K.edgecut[2][1]==1) {    size_t nbInter=vertexInTriangles.size();
1924      place0=2;    size_t nbCut=K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0];
1925      size_t numberIn=K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2];
1926      size_t stop=0,iInter=0,iCut=0;
1927      //initialisation
1928      if(numberIn==3) {
1929          points.push_back(&K(1));
1930          points.push_back(&K(2));
1931      }
1932      if(numberIn==2 and stop==0) {
1933        points.push_back(&K(in1));
1934      }
1935      const Vertex* A=&K(in);
1936      if(numberIn>0) {
1937        points.push_back(A);
1938        //ffout(4)<<*A<<"\n";
1939    }    }
1940      if(nbCut>0) {
1941    const Vertex* Pend=K.pointsIntersection[place0][1];      P[0]=K.pointsIntersection[place[0]][0];
1942    #warning size TinyVector      //attention si on a 2 points d'inter et un point in choisir le point le plus proche de A
1943    TinyVector<30,const Vertex*> P;      if(K.edgecut[place[0]][1]!=0) {
1944    TinyVector<30,const Triangle*> T;        P[1]=K.pointsIntersection[place[0]][1];
1945    P[0]=K.pointsIntersection[place0][0];        if(numberIn>0) {
1946            real_t norm0=Norm(*A-*P[0]);
1947    if(vertexInTriangles.find(P[0])!=vertexInTriangles.end()) {          real_t norm1=Norm(*A-*P[1]);
1948      T[0]=(*vertexInTriangles.find(P[0])).second[0];          if(norm1<norm0) {
1949    } else {            ffout(4)<<"--------------- cas inverse\n";
1950      ffout(4)<<"pbs 0 existe pas\n";            //P[0]=P[1];
1951            }
1952          }
1953        }
1954        if(P[0]!=A  and nbInter!=0) {
1955          points.push_back(P[0]);
1956        }
1957      } else if(nbCut!=0) {
1958        if(P[0]==A) {
1959          ffout(4)<<"pbs point confondu A\n";
1960        }
1961        std::map<const Vertex*,TinyVector<2,const Triangle*> >::iterator itvert=vertexInTriangles.begin();
1962        P[0]=(*itvert).first;
1963      }
1964      if(nbCut!=0) {
1965        if(vertexInTriangles.find(P[0])!=vertexInTriangles.end()) {
1966          T[0]=(*vertexInTriangles.find(P[0])).second[0];
1967          iInter=1;
1968        } else {
1969          ffout(4)<<"******\n";
1970          ffout(4)<<"nb Inter "<<nbInter<<"\n";
1971          ffout(4)<<"P0 = "<<P[0]<<*P[0]<<"\n";
1972          ffout(4)<<"K(2) = "<<K(1)<<&K(1)<<"\n";
1973          ffout(4)<<"pbs 0 existe pas\n";
1974          ffout(4)<<*K.pointsIntersection[1][0]<<"\n";
1975          ffout(4)<<*K.pointsIntersection[2][0]<<"\n";
1976          ffout(4)<<*K.pointsIntersection[2][1]<<"\n";
1977          stop=1;
1978        }
1979      }
1980      
1981      //ffout(4)<<"nb de points d'inter "<<nbInter<<"\n";
1982      if(stop==0 and nbCut==2 ) {
1983        if(K.edgecut[place[1]][0]==1 and K.edgecut[place[0]][0]==1) {
1984          if(Norm(*K.pointsIntersection[place[1]][0]-*K.pointsIntersection[place[0]][0])<1e-12) {
1985            ffout(4)<<"les points sont confondus\n";
1986            ffout(4)<<*K.pointsIntersection[place[0]][0]<<"\n";
1987            stop=1;
1988          }
1989        }
1990    }    }
1991    
1992    size_t b=0;    if(iInter==nbInter and nbInter==1 and numberIn==1) {
1993        //ATTENTION PEUT ETRE PLUS DE POINTS A PRENDRE
1994    while(P[b]!=Pend and b<29) {      if(nbCut>1) {
1995      if(triangleWithVertex.find(T[b])!=triangleWithVertex.end()) {        const Vertex* PP=K.pointsIntersection[place[1]][0];
1996        TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[b])).second;        if(PP!=P[0]) {
1997        if(VV[0]==P[b]) {          points.push_back(PP);
         if(Pend!=VV[1]) {  
 //#warning Reference not set.  
           __addTriangle(Triangle(*Pend,*P[b],*VV[1],K.reference()),cell,triangleListIntersectionNew);  
         }  
         P[b+1]=VV[1];  
1998        } else {        } else {
1999          if(Pend!=VV[0]) {          if(K.edgecut[place[2]]!=0) {
2000  //#warning Reference not set.            const Vertex* PP1=K.pointsIntersection[place[2]][0];
2001            __addTriangle(Triangle(*Pend,*P[b],*VV[0],K.reference()),cell,triangleListIntersectionNew);            if(PP1!=P[0]) {
2002                points.push_back(PP1);
2003              }
2004          }          }
         P[b+1]=VV[0];  
2005        }        }
     } else {  
       ffout(4)<<"pbs 1 existe pas\n";  
       b=28;  
2006      }      }
2007      TinyVector<2,const Triangle*> T2;    }
2008      if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {    
2009        T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];    while(iInter<nbInter and stop==0) {
2010        T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];      //ffout(4)<<"createGeneral "<<iInter<<" "<<nbInter<<"\n";
2011      }      //on cherche le voisin de Pi
2012      if(T2[0]==T[b]) {      __findVertex(K,P,T,iCut,iInter,stop,in,in1,place,points);
2013        if(T2[1]==T[b]) {  
2014          //point edge      //on cherche le triangle associe au nouveau Pi pour relancer
2015        } else {      size_t k=0;
2016          T[b+1]=T2[1];      while(stop==3 or k==0) {
2017          //ffout(4)<<k<<" findTriangle\n";
2018          __findTriangle(P,T,stop);
2019          if(stop==3) {
2020            //ffout(4)<<"stop3\n";
2021            __findInEdges(K,P,T,iInter,stop,in,in1,place,points);
2022        }        }
2023      } else {        ++k;
       T[b+1]=T2[0];  
2024      }      }
2025      ++b;      
2026    }    }
2027    if(b==29) ffout(4)<<"createnotin b ="<<b<<"\n";      
2028  }  }
2029    
2030  void SurfaceMeshGenerator::Internals::  bool SurfaceMeshGenerator::Internals::
2031  __createNotIn2(const TriangleCut& K,const Cell* cell,  __verifExist(const Vertex* Pt,
2032                 std::vector<Triangle>& triangleListIntersectionNew)               std::vector<const Vertex*>& points)
2033  {// create triangles if 2 edges is cut  {
2034      // and more than 3 points intersection?    bool exist=false;
2035    size_t place0=0;    for(size_t k=0 ; k<points.size() ; ++k) {
2036    if(K.edgecut[1][0]==1 and K.edgecut[1][1]==1) {      if(points[k]==Pt) {
2037      place0=1;        //ffout(4)<<"je l'ai deja!!\n";
2038    } else        exist=true;
     if(K.edgecut[2][0]==1 and K.edgecut[2][1]==1) {  
       place0=2;  
2039      }      }
2040      }
2041      return exist;
2042    }
2043    
2044    //ffout(4)<<"place = "<<place0<<"\n";  bool SurfaceMeshGenerator::Internals::
2045    const Vertex* Pend=K.pointsIntersection[place0][1];  __findOneEdge(const TriangleCut& K,
2046    TinyVector<30,const Vertex*> P;                const size_t & number,
2047    TinyVector<30,const Triangle*> T;                TinyVector<2,const Vertex*>& P)
2048    P[0]=K.pointsIntersection[place0][0];  {   bool find=false;
2049    if(vertexInTriangles.find(P[0])!=vertexInTriangles.end()) {   if(K.edgecut[number][0]!=0) {
2050      T[0]=(*vertexInTriangles.find(P[0])).second[0];     if(P[0]==K.pointsIntersection[number][0]) {
2051         //ffout(4)<<"0"<<number<<"\n";
2052         find=true;
2053       }
2054     }
2055     if(K.edgecut[number][1]!=0 and find==false) {
2056       if(P[0]==K.pointsIntersection[number][1]) {
2057         //ffout(4)<<"1"<<number<<"\n";
2058         find=true;
2059       }
2060     }
2061     return find;
2062    }
2063    
2064    
2065    void SurfaceMeshGenerator::Internals::
2066    __findInEdges(const TriangleCut& K,
2067                  TinyVector<2,const Vertex*>& P,
2068                  TinyVector<2,const Triangle*>& T,
2069                  size_t& iInter,
2070                  size_t& stop,
2071                  size_t& in,
2072                  size_t& in1,
2073                  TinyVector<3,size_t> place,
2074                  std::vector<const Vertex*>& points)
2075    {
2076      //size_t nbCut=K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0];
2077      size_t nbInter=vertexInTriangles.size();
2078      size_t iCut=0;
2079      //ffout(4)<<*P[0]<<"\n";
2080      //find on which edge we are
2081      bool find=false;
2082      while(find==false and iCut<3) {
2083        find=__findOneEdge(K,iCut,P);
2084        ++iCut;
2085      }
2086      if(iCut==3 and find==false) {
2087        if(iInter!=iInter) {
2088          ffout(4)<<iInter<<" "<<iInter<<"\n";
2089          ffout(4)<<"pbs on a plus 3\n";
2090        }
2091        stop=1;
2092        iCut=0;
2093    } else {    } else {
2094      ffout(4)<<"pbs 0 existe pas\n";      --iCut;
2095    }    }
2096        
2097    size_t b=0;    //ffout(4)<<"iCut "<<iCut<<"\n";
2098    while(P[b]!=Pend and b<29) {    if(K.edgecut[iCut][1]!=0 and stop!=1) {
2099      //ffout(4)<<"------------------- b "<<b<<"\n";      //verifier que le point n'est pas deja dans la liste
2100      if(triangleWithVertex.find(T[b])!=triangleWithVertex.end()) {      const Vertex* Ptemp=K.pointsIntersection[iCut][1];
2101        TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[b])).second;      bool exist=__verifExist(Ptemp,points);
2102        if(VV[0]==P[b]) {      if(!exist) {
2103          if(Pend!=VV[1]) {        P[0]=Ptemp;
2104            __addTriangle(Triangle(*Pend,*P[b],*VV[1],K.reference()),cell,triangleListIntersectionNew);        points.push_back(P[0]);
2105            //  ffout(4)<<"1 "<<*Pend<<"\n";        //ffout(4)<<P[0]<<"\n";
2106            //  ffout(4)<<*P[b]<<"\n";        //ffout(4)<<*P[0]<<"\n";
2107            //  ffout(4)<<*VV[1]<<"\n";        ++iInter;
2108          }      } else {
2109          P[b+1]=VV[1];        exist=__verifExist(K.pointsIntersection[iCut][0],points);
2110          if(!exist) {
2111            P[0]=K.pointsIntersection[iCut][0];
2112            points.push_back(P[0]);
2113            //ffout(4)<<P[0]<<"\n";
2114            //ffout(4)<<*P[0]<<"\n";
2115            ++iInter;
2116        } else {        } else {
2117          if(Pend!=VV[0]) {          if(nbInter!=iInter) {
2118              __addTriangle(Triangle(*Pend,*P[b],*VV[0],K.reference()),cell,triangleListIntersectionNew);            ffout(4)<<nbInter<<" "<<iInter<<"\n";
2119            //ffout(4)<<"2 "<<*Pend<<"\n";            ffout(4)<<"pbs on a plus de points\n";
           //ffout(4)<<*P[b]<<"\n";  
           //ffout(4)<<*VV[0]<<"\n";  
2120          }          }
2121          P[b+1]=VV[0];          stop=1;
2122        }        }
     } else {  
       ffout(4)<<"pbs 1 existe pas\n";  
       b=28;  
     }  
         
     TinyVector<2,const Triangle*> T2;  
     if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {  
       T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];  
       T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];  
     } else {  
       ffout(4)<<"points exisste pas \n";  
       b=28;  
2123      }      }
2124          } else {
2125      if(K.edgecut[0][0]==1 and K.edgecut[0][1]==1) {      if(nbInter!=iInter) {
2126        if(P[b+1]==K.pointsIntersection[0][1] or P[b+1]==K.pointsIntersection[0][0]) {        ffout(4)<<nbInter<<" "<<iInter<<"\n";
2127          T2[0]=T[b];        ffout(4)<<"pbs on a plus de points2\n";
2128          T2[1]=T[b];      }
2129        }      stop=1;
2130      } else    }
2131        if(K.edgecut[2][0]==1 and K.edgecut[2][1]==1) {  }
2132          if(P[b+1]==K.pointsIntersection[2][1] or P[b+1]==K.pointsIntersection[2][0]) {  
2133            T2[0]=T[b];  void SurfaceMeshGenerator::Internals::
2134            T2[1]=T[b];  __findVertex(const TriangleCut& K,
2135          }               TinyVector<2,const Vertex*>& P,
2136        }               TinyVector<2,const Triangle*>& T,
2137                     size_t& iCut,
2138      if(T2[0]==T[b]) {               size_t& iInter,
2139        if(T2[1]==T[b]) {               size_t& stop,
2140          //point edge               size_t& in,
2141          //ffout(4)<<"here\n";               size_t& in1,
2142          if(K.edgecut[0][0]==1 and K.edgecut[0][1]==1 and place0!=0) {               TinyVector<3,size_t> place,
2143            //ffout(4)<<"here2\n";               std::vector<const Vertex*>& points)
2144            if(P[b+1]==K.pointsIntersection[0][1]) {  {
2145              __addTriangle(Triangle(*Pend,*P[b+1],*K.pointsIntersection[0][0],    //size_t nbCut=K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0];
2146                              K.reference()),cell,triangleListIntersectionNew);    //size_t nbInter=vertexInTriangles.size();
2147              P[b+1]=K.pointsIntersection[0][0];                      
2148              if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {    if(triangleWithVertex.find(T[0])!=triangleWithVertex.end()) {
2149                T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];      TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[0])).second;
2150                T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];      if(VV[0]==P[0]) {
2151              }        if(VV[1]==P[0]) {
2152              T[b+1]=T2[0];          //on est a priori sur une arete
2153            } else            ffout(4)<<"cas findinedges\n";
2154              if(P[b+1]==K.pointsIntersection[0][0]){          __findInEdges(K,P,T,iInter,stop,in,in1,place,points);
               __addTriangle(Triangle(*Pend,*P[b+1],*K.pointsIntersection[0][1],  
                               K.reference()),cell,triangleListIntersectionNew);  
               P[b+1]=K.pointsIntersection[0][1];  
               if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {  
                 T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];  
                 T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];  
               } else {  
                 ffout(4)<<"le tr e pas\n";  
               }  
               T[b+1]=T2[0];  
             } else  
               if(P[b+1]!=Pend) {  
                 ffout(4)<<"point pas sur edge trou???? \n";  
                 ffout(4)<<*P[b+1]<<"\n";  
                 ffout(4)<<*Pend<<"\n";  
                 P[b+1]=Pend;  
               }  
         } else {  
           if(K.edgecut[2][0]==1 and K.edgecut[2][1]==1 and place0!=2) {  
             if(P[b+1]==K.pointsIntersection[2][1]) {  
               __addTriangle(Triangle(*Pend,*P[b+1],*K.pointsIntersection[2][0],  
                               K.reference()),cell,triangleListIntersectionNew);  
               P[b+1]=K.pointsIntersection[2][0];  
               if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {  
                 T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];  
                 T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];  
               }  
               T[b+1]=T2[0];  
             } else {  
               if(P[b+1]==K.pointsIntersection[2][0]){  
                 __addTriangle(Triangle(*Pend,*P[b+1],*K.pointsIntersection[2][1],  
                          K.reference()),cell,triangleListIntersectionNew);  
                 P[b+1]=K.pointsIntersection[2][1];  
                 if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {  
                   T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];  
                   T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];  
                 }  
                 T[b+1]=T2[0];  
               }  
             }  
           } else {  
             ffout(4)<<"pas trouve l'egde et le point\n";  
           }  
         }  
2155        } else {        } else {
2156          T[b+1]=T2[1];          //ffout(4)<<"cas 2\n";
2157            P[1]=VV[1];
2158            P[0]=P[1];
2159            ffout(4)<<P[0]<<"\n";
2160            ffout(4)<<*P[0]<<"\n";
2161            points.push_back(P[0]);
2162            ++iInter;
2163        }        }
2164      } else {      } else {
2165        T[b+1]=T2[0];        //ffout(4)<<"cas 3\n";
2166          P[1]=VV[0];
2167          P[0]=P[1];
2168          //ffout(4)<<P[0]<<"\n";
2169          //ffout(4)<<*P[0]<<"\n";
2170          bool exist=__verifExist(P[0],points);
2171          if(!exist) {
2172            points.push_back(P[0]);
2173            ++iInter;
2174          } else {
2175            ffout(4)<<"on a un pbs le point existe deja!!\n";
2176            stop=1;
2177          }
2178          
2179    
2180      }      }
     ++b;  
2181    }    }
   if(b==29) ffout(4)<<"createnotin2 b ="<<b<<"\n";  
2182  }  }
2183    
2184  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
2185  __createByEdges(const TriangleCut& K,const Cell* cell,  __findTriangle(TinyVector<2,const Vertex*>& P,
2186                  std::vector<Triangle>& triangleListIntersectionNew)                 TinyVector<2,const Triangle*>& T,
2187                   size_t& stop)
2188  {  {
2189      //createtriangles if one point is in Triangle    if(stop==3) {
2190      // and two or 3? edges are cuts      stop=0;
2191    size_t in0=0,place0=0,place1=1,place2=2;    }
2192      TinyVector<2,const Triangle*> T2;
2193    if(K.isInTriangle[1]==1) {    if(vertexInTriangles.find(P[0])!=vertexInTriangles.end()) {
2194      in0=1;      T2[0]=(*vertexInTriangles.find(P[0])).second[0];
2195      place1=2;      T2[1]=(*vertexInTriangles.find(P[0])).second[1];
     place2=1;  
2196    } else {    } else {
2197      if(K.isInTriangle[2]==1) {      ffout(4)<<"pbs tr existe pas\n";
2198        in0=2;      stop=1;
       place0=1;  
       place1=2;  
       place2=0;  
     }  
2199    }    }
2200    
2201    const Vertex* Pend=K.pointsIntersection[place1][0];    if(T2[0]==T[0]) {
2202    #warning size of Tiny      if(T2[1]==T[0]) {
2203    size_t b=0;        stop=3;
2204    TinyVector<30,const Vertex*> P;      } else {
2205    TinyVector<30,const Triangle*> T;        T[0]=T2[1];
2206    P[0]=K.pointsIntersection[place0][0];      }
   if(vertexInTriangles.find(P[0])!=vertexInTriangles.end()) {  
     T[0]=(*vertexInTriangles.find(P[0])).second[0];  
2207    } else {    } else {
2208      b=29;        T[0]=T2[0];
     ffout(4)<<"pbs 0 existe pas\n";  
2209    }    }
2210        
2211    }
2212    
2213    while(P[b]!=Pend and b<29) {  void SurfaceMeshGenerator::Internals::
2214      if(triangleWithVertex.find(T[b])!=triangleWithVertex.end()) {  __createTriangle(const TriangleCut& K,const Cell* cell,
2215        TinyVector<2,const Vertex*> VV=(*triangleWithVertex.find(T[b])).second;                   std::vector<Triangle>& triangleListIntersectionNew,
2216        if(VV[0]==P[b]) {                   std::vector<const Vertex*>& points,
2217          if(VV[1]==P[b]) {                   std::vector<Vertex>& pointstrian) {
2218            ffout(4)<<"b ="<<b<<" egal\n";  
2219            if(K.pointsIntersection[place2][0]==P[b]) {      size_t taille=points.size();
2220              P[b+1]=K.pointsIntersection[place2][1];      Triangulation trgen;
2221              __addTriangle(Triangle(K(in0),*P[b],*P[b+1],K.reference()),cell,triangleListIntersectionNew);      Triangulation::CurveVertex envVertex(taille+1);
2222            } else {      std::ofstream sortie("testy2mb");
2223              if(K.pointsIntersection[place2][1]==P[b]){      sortie.precision(30);
2224                P[b+1]=K.pointsIntersection[place2][0];      sortie<<taille<<" 0 0 \n";
2225                __addTriangle(Triangle(K(in0),*P[b],*P[b+1],K.reference()),cell,triangleListIntersectionNew);      for(unsigned ie=0;ie<taille;++ie) {
2226              } else {          sortie<<pointstrian[ie][0]<<" "<<pointstrian[ie][1]<<"\n";
2227                ffout(4)<<"pbs pbs\n";          //ffout(4)<<points[ie]<<"\n";
2228                ffout(4)<<*P[b]<<"\n";   //ffout(4)<<*points[ie]<<"\n";
2229                P[b+1]=K.pointsIntersection[place2][1];   //ffout(4)<<pointstrian[ie]<<"\n";
2230                //std::exit(0);   //ffout(4)<<"----\n";
2231            envVertex[ie] = &pointstrian[ie];
2232        }
2233        envVertex[taille] = envVertex[0];
2234        sortie<<taille<<" ";
2235        for(unsigned ie=0;ie<taille;++ie) {
2236            sortie<<envVertex[ie]-envVertex[0]<<" ";
2237        }
2238    
2239        sortie.close();
2240    
2241        //dans ce cas pas de trou
2242    std::vector<Triangulation::CurveVertex> holes(0);
2243    std::vector<Triangulation::CurveVertex> curves(0);
2244        
2245        if(taille>3) {
2246            //ffout(4)<<"je triangule\n";
2247            //ffout(4)<<++toto<<"\n";
2248            //verifions que les points ne s'intersectent pas avant de lancer la triangulation
2249            bool isOK=false;
2250            //ffout(4)<<"nb de points "<<taille<<"\n";
2251            size_t ie=1;
2252            while(ie<taille-2 and !isOK) {
2253                isOK=checkInterbis(pointstrian[0],pointstrian[taille-1],pointstrian[ie],pointstrian[ie+1]);
2254                ++ie;
2255            }
2256            if(!isOK) {
2257                //verification que les points ne sont pas trop proches!!
2258                ie=0;
2259                while(ie<taille-1 and !isOK) {
2260                    const Vertex& A=pointstrian[ie];
2261                    const Vertex& B=pointstrian[ie+1];
2262                    real_t norm=Norm(A-B);
2263                    isOK=(norm<1e-6);
2264                    ++ie;
2265                }
2266            }
2267            if(!isOK) {
2268                trgen.triangulize(pointstrian,
2269                           envVertex,
2270                           holes,
2271                           curves);
2272                //ffout(4)<<"fin triangule-----\n";
2273                std::ofstream o("full.mesh");
2274                trgen.export_mesh(o);
2275                std::list<ConnectedTriangle> L=trgen.getTriangles();
2276                for(std::list<ConnectedTriangle>::iterator itL=L.begin() ; itL!= L.end() ; itL++) {
2277                    Triangle T=(*itL);
2278                    Vertex* P0=&T(0);
2279                    Vertex* P1=&T(1);
2280                    Vertex* P2=&T(2);
2281                    __addTriangle(Triangle(*points[P0-&pointstrian[0]],
2282                             *points[P1-&pointstrian[0]],
2283                             *points[P2-&pointstrian[0]],
2284                             K.reference()),cell,triangleListIntersectionNew);
2285              }              }
           }  
2286          } else {          } else {
2287            __addTriangle(Triangle(K(in0),*P[b],*VV[1],K.reference()),cell,triangleListIntersectionNew);              ffout(4)<<"isOK est vraie on a un pbs.....\n";
2288            P[b+1]=VV[1];              ffout(4)<<K(0)<<"\n";
2289          }              ffout(4)<<K(1)<<"\n";
2290        } else {              ffout(4)<<K(2)<<"\n";
2291          __addTriangle(Triangle(K(in0),*P[b],*VV[0],K.reference()),cell,triangleListIntersectionNew);          }
2292          P[b+1]=VV[0];          //ffout(4)<<L.size()<<" ecrit\n";
2293        }      } else if(taille==3) {
2294            ffout(4)<<"je cree un triangle\n";
2295            __addTriangle(Triangle(*points[0],
2296                            *points[1],
2297                            *points[2],
2298                            K.reference()),cell,triangleListIntersectionNew);
2299      } else {      } else {
2300        ffout(4)<<"pbs 1 existe pas\n";          ffout(4)<<"la taille est inferieur a 3 = "<<taille<<"\n";
2301        b=28;          ffout(4)<<K(0)<<"\n";
2302            ffout(4)<<K(1)<<"\n";
2303            ffout(4)<<K(2)<<"\n";
2304            ffout(4)<<"¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡\n";
2305      }      }
2306      TinyVector<2,const Triangle*> T2;  }
2307    
2308      if(vertexInTriangles.find(P[b+1])!=vertexInTriangles.end()) {  void SurfaceMeshGenerator::Internals::
2309        T2[0]=(*vertexInTriangles.find(P[b+1])).second[0];  __createGeneral(const TriangleCut& K,const Cell* cell,
2310        T2[1]=(*vertexInTriangles.find(P[b+1])).second[1];                  std::vector<Triangle>& triangleListIntersectionNew)
2311      } else {  {
2312        ffout(4)<<"pbs tr existe pas\n";    static int toto=0;
2313        b=28;    // pour l'instant on suppose edge cut une seule fois
2314      }    size_t in=0,in1=1;
2315      TinyVector<3,size_t> place;
2316      std::vector<const Vertex*> points;
2317      
2318      real_t epsilon=1e-7;
2319      real_t norm01=Norm(K(0)-K(1));
2320      real_t norm02=Norm(K(0)-K(2));
2321      real_t norm21=Norm(K(2)-K(1));
2322      //test si le triangle est plat!! cad si 2 points sont "confondus"
2323      if(norm01<epsilon or norm02<epsilon or norm21<epsilon) {
2324        ffout(4)<<"le triangle est tres tres tres plat\n";
2325      } else {
2326    
2327      if(T2[0]==T[b] and T2[1]==T[b]) {      //pour determiner le numero des edges cut et ceux des points in
2328        // ffout(4)<<"here0\n";      __createCase(K,in,in1,place);
2329        //point edge  
2330        if(P[b+1]!=Pend) {      //pour 3 edges cut et un in il faut changer place
2331          if(K.pointsIntersection[place2][0]==P[b+1]) {      if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==3 and
2332            const Vertex* A=K.pointsIntersection[place2][1];         K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==1) {
2333  //#warning Reference not set.        size_t temp=place[1];
2334            __addTriangle(Triangle(K(in0),*P[b+1],*A,K.reference()),cell,triangleListIntersectionNew);        place[1]=place[2];
2335            if(vertexInTriangles.find(A)!=vertexInTriangles.end()) {        place[2]=temp;
2336              T[b+1]=(*vertexInTriangles.find(A)).second[0];      }
2337              P[b+1]=A;  
2338            }         /////////////////////////////////////////////////////////////////////////////
2339       /* if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==3
2340           and K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==1 ) {
2341          ffout(4)<<norm01<<" "<<norm02<<" "<<norm21<<"\n";
2342          ffout(4)<<in<<"\n";
2343          ffout(4)<<place[0]<<" "<<place[1]<<" "<<place[2]<<"\n";
2344          ffout(4)<<K(0)<<"\n";
2345          ffout(4)<<K(1)<<"\n";
2346          ffout(4)<<K(2)<<"\n";
2347          ffout(4)<<"points in "<<&K(in)<<"\n";
2348          ffout(4)<<K.pointsIntersection[place[0]][0]<<"\n";
2349          ffout(4)<<K.pointsIntersection[place[1]][0]<<"\n";
2350            ffout(4)<<K.pointsIntersection[place[2]][0]<<"\n";
2351          ffout(4)<<"inter 0 "<<*K.pointsIntersection[place[0]][0]<<"\n";
2352          ffout(4)<<"inter 1"<<*K.pointsIntersection[place[1]][0]<<"\n";
2353          ffout(4)<<*K.pointsIntersection[place[2]][0]<<"\n";
2354          ffout(4)<<"inter 0 1 "<<K.pointsIntersection[place[0]][1]<<"\n";
2355          ffout(4)<<"inter 1 1"<<K.pointsIntersection[place[1]][1]<<"\n";
2356          ffout(4)<<"inter 0 1 "<<*K.pointsIntersection[place[0]][1]<<"\n";
2357          ffout(4)<<"inter 1 1"<<*K.pointsIntersection[place[1]][1]<<"\n";
2358    
2359        }*/
2360        /////////////////////////////////////////////////////////////////////////////
2361        
2362        //add vertex to VertexInTriangles if need
2363        __addList(K);
2364    
2365        //create the list "points"
2366        __createList(K,points,in,in1,place);
2367    
2368        size_t taille=points.size();
2369        // met les points dans le plan du triangle
2370        std::vector<Vertex> pointstrian(taille);
2371        __transformVertex(K,points,pointstrian);
2372        
2373        if(K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==2 and
2374           K.edgecut[place[2]][0]!=0) {
2375            ffout(4)<<"on est dans un cas a traiter pour le mailleur\n";
2376            ffout(4)<<K(0)<<"\n";
2377            ffout(4)<<K(1)<<"\n";
2378            ffout(4)<<K(2)<<"\n";
2379            // on veut decomposer la liste en 2 sous listes pour creer 2 sous domaines
2380            size_t number=K.edgecut[place[2]][0]+K.edgecut[place[2]][1];
2381            // pour determiner le ou les points ou on separe les listes
2382            TinyVector<2,const Vertex*> P;
2383            switch (number) {
2384                case 1: {
2385                    P[0]=K.pointsIntersection[place[2]][0];
2386                    break;
2387                }
2388                case 2: {
2389                    P[0]=K.pointsIntersection[place[2]][0];
2390                    P[1]=K.pointsIntersection[place[2]][1];
2391                    break;
2392                }
2393            }
2394            std::vector<const Vertex*> pointslist1;
2395            std::vector<const Vertex*> pointslist2;
2396            
2397            size_t stopcut=0;
2398            for(unsigned ie=1;ie<taille;++ie) {
2399                const Vertex* Pcut=points[ie];
2400                //ffout(4)<<*Pcut<<"\n";
2401                const Vertex& Pcutplan=pointstrian[ie];
2402                if(Pcut==P[0]) {
2403                    __addPoints(Pcut,stopcut,pointslist1,pointslist2);
2404                    if(stopcut==3) {
2405                        stopcut=1;
2406                    }
2407                } else {
2408                    if(number==2) {
2409                        if(Pcut==P[1]) {
2410                            __addPoints(Pcut,stopcut,pointslist1,pointslist2);
2411                            if(stopcut==3) {
2412                                stopcut=1;
2413                            }
2414                        } else {
2415                            __addPoints(Pcut,stopcut,pointslist1,pointslist2);
2416                            if(stopcut==3) {
2417                                stopcut=0;
2418                            }
2419                        }
2420                    } else {
2421                        __addPoints(Pcut,stopcut,pointslist1,pointslist2);
2422                        if(stopcut==3) {
2423                            stopcut=0;
2424                        }
2425                    }
2426                }
2427            }
2428            pointslist2.push_back(points[0]);
2429            /*//verif des listes
2430            ffout(4)<<"taille de list1 = "<<pointslist1.size()<<"\n";
2431            for(unsigned ie=0;ie<pointslist1.size();++ie) {
2432                ffout(4)<<*pointslist1[ie]<<"\n";
2433            }
2434            ffout(4)<<"taille de list2 = "<<pointslist2.size()<<"\n";
2435            for(unsigned ie=0;ie<pointslist2.size();++ie) {
2436                ffout(4)<<*pointslist2[ie]<<"\n";
2437            }*/
2438            //on cree les 2 maillages
2439            if(pointslist1.size()>=3) {
2440                std::vector<Vertex> pointslist1plan(pointslist1.size());
2441                __transformVertex(K,pointslist1,pointslist1plan);
2442                __createTriangle(K,cell,triangleListIntersectionNew,pointslist1,pointslist1plan);
2443            }
2444            if(pointslist2.size()>=3) {
2445                 std::vector<Vertex> pointslist2plan(pointslist2.size());
2446                __transformVertex(K,pointslist2,pointslist2plan);
2447                __createTriangle(K,cell,triangleListIntersectionNew,pointslist2,pointslist2plan);
2448            }
2449        } else if(K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2]==3) {
2450            ffout(4)<<"cas 3 in\n";
2451            size_t number=K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0];
2452            if(number==1) {
2453                TinyVector<2,const Vertex*> P;
2454                size_t numberP=K.edgecut[place[0]][0]+K.edgecut[place[0]][1];
2455                switch (numberP) {
2456                    case 1: {
2457                        P[0]=K.pointsIntersection[place[0]][0];
2458                        break;
2459                    }
2460                    case 2: {
2461                        P[0]=K.pointsIntersection[place[0]][0];
2462                        P[1]=K.pointsIntersection[place[0]][1];
2463                        break;
2464                    }
2465                }
2466                std::vector<const Vertex*> pointslist1;
2467                std::vector<const Vertex*> pointslist2;
2468                size_t stopcut=0;
2469                for(unsigned ie=1;ie<taille;++ie) {
2470                    const Vertex* Pcut=points[ie];
2471                    //ffout(4)<<*Pcut<<"\n";
2472                    if(Pcut==P[0]) {
2473                        __addPoints(Pcut,stopcut,pointslist1,pointslist2);
2474                        if(stopcut==3) {
2475                            stopcut=1;
2476                        }
2477                    } else {
2478                        if(number==2) {
2479                            if(Pcut==P[1]) {
2480                                __addPoints(Pcut,stopcut,pointslist1,pointslist2);
2481                                if(stopcut==3) {
2482                                    stopcut=1;
2483                                }
2484                            } else {
2485                                __addPoints(Pcut,stopcut,pointslist1,pointslist2);
2486                                if(stopcut==3) {
2487                                    stopcut=0;
2488                                }
2489                            }
2490                        } else {
2491                            __addPoints(Pcut,stopcut,pointslist1,pointslist2);
2492                            if(stopcut==3) {
2493                                stopcut=0;
2494                            }
2495                        }
2496                    }
2497                }
2498                pointslist2.push_back(points[0]);
2499                pointslist2.push_back(points[1]);
2500                pointslist2.push_back(P[0]);
2501    
2502                //on cree les 2 maillages
2503                if(pointslist1.size()>=3) {
2504                     std::vector<Vertex> pointslist1plan(pointslist1.size());
2505                    __transformVertex(K,pointslist1,pointslist1plan);
2506                    __createTriangle(K,cell,triangleListIntersectionNew,pointslist1,pointslist1plan);
2507                }
2508                if(pointslist2.size()>=3) {
2509                    std::vector<Vertex> pointslist2plan(pointslist2.size());
2510                    __transformVertex(K,pointslist2,pointslist2plan);
2511                    __createTriangle(K,cell,triangleListIntersectionNew,pointslist2,pointslist2plan);
2512                }
2513                
2514                
2515          } else {          } else {
2516            const Vertex* A=K.pointsIntersection[place2][0];              ffout(4)<<"il faut le traiter\n";
 //#warning Reference not set.  
           __addTriangle(Triangle(K(in0),*P[b+1],*A,K.reference()),cell,triangleListIntersectionNew);  
           if(vertexInTriangles.find(A)!=vertexInTriangles.end()) {  
             T[b+1]=(*vertexInTriangles.find(A)).second[0];  
             P[b+1]=A;  
           }  
2517          }          }
2518        }          //for(unsigned ie=0;ie<points.size();++ie) {
2519      } else {          //    ffout(4)<<*points[ie]<<"\n";
2520        if(T2[0]==T[b]) {          //}
2521          T[b+1]=T2[1];      } else if(taille!=0) {
2522        } else {          __createTriangle(K,cell,triangleListIntersectionNew,points,pointstrian);
         T[b+1]=T2[0];  
       }  
2523      }      }
     ++b;  
2524    }    }
   if(b==29) ffout(4)<<"createbyedges b ="<<b<<"\n";  
2525  }  }
2526    
 void SurfaceMeshGenerator::Internals::__3pointsintersection(const TriangleCut& K,const Cell* cell,  
                                                             std::vector<Triangle>& triangleListIntersectionNew)  
 {  
     TinyVector<3,const Vertex*> V;  
     size_t num=0;  
     for(std::map<const Vertex*,TinyVector<2,const Triangle*> >::iterator  
         it=vertexInTriangles.begin() ; it!=vertexInTriangles.end() ; ++it , ++num) {  
         V[num]=(*it).first;  
     }  
     if(V[0]!=V[1] and V[1]!=V[2] and V[2]!=V[0]) {  
         __addTriangle(Triangle(*V[0],*V[1],*V[2],K.reference()),cell,triangleListIntersectionNew);  
         ffout(4)<<"on comble le trou\n"<<std::flush;  
     } else {  
         ffout(4)<<"----- deux inter =\n";  
         ffout(4)<<V[0]<<"\n";  
         ffout(4)<<V[1]<<"\n";  
         ffout(4)<<V[2]<<"\n";  
     }  
 }  
2527  template <typename booleanTest>  template <typename booleanTest>
2528  void SurfaceMeshGenerator::Internals::  void SurfaceMeshGenerator::Internals::
2529  __createTrianglesIntersection(const size_t numobject,  __createTrianglesIntersection(const size_t numobject,
# Line 2196  __createTrianglesIntersection(const size Line 2533  __createTrianglesIntersection(const size
2533                                std::vector<Triangle>& triangleListIntersectionNew,                                std::vector<Triangle>& triangleListIntersectionNew,
2534                                PIntersection& listVertexIntersectionNew)                                PIntersection& listVertexIntersectionNew)
2535  {  {
2536      //build mesh of object O1 (which is intersected with O2)    //build mesh of object O1 (which is intersected with O2)
2537    
2538    size_t nbcell=(*O1.trianglelist).size();    size_t nbcell=(*O1.trianglelist).size();
2539    SurfaceMeshOfTriangles surfmesh(nbcell);    SurfaceMeshOfTriangles surfmesh(nbcell);
# Line 2207  __createTrianglesIntersection(const size Line 2544  __createTrianglesIntersection(const size
2544    size_t objectNumber=0;    size_t objectNumber=0;
2545    if(numobject==0)    if(numobject==0)
2546      objectNumber=1;      objectNumber=1;
2547   // (*edgesRefVertex[numobject]).clear();    // (*edgesRefVertex[numobject]).clear();
2548    const MapCellTriangle& hexaToTriangle0=*(O1.hexalist);    const MapCellTriangle& hexaToTriangle0=*(O1.hexalist);
2549    
2550        
# Line 2228  __createTrianglesIntersection(const size Line 2565  __createTrianglesIntersection(const size
2565      MapCellTriangle::const_iterator it0=hexaToTriangle0.find(*icell);      MapCellTriangle::const_iterator it0=hexaToTriangle0.find(*icell);
2566            
2567      if(it0!=hexaToTriangle0.end()) {      if(it0!=hexaToTriangle0.end()) {
2568          //ffout(4)<<"ici\n";        //ffout(4)<<"ici\n";
2569          Cell& C = const_cast<Cell&>(*(*icell));        Cell& C = const_cast<Cell&>(*(*icell));
2570          for (size_t i=0; i<C.numberOfVertices(); ++i) {        for (size_t i=0; i<C.numberOfVertices(); ++i) {
2571              C(i).reference() = ((*O1.shape()).inside(C(i))?1:0)          C(i).reference() = ((*O1.shape()).inside(C(i))?1:0)
2572              + ((*O2.shape()).inside(C(i))?1:0);            + ((*O2.shape()).inside(C(i))?1:0);
2573          }        }
2574    
2575          size_t ref;        size_t ref;
2576          if (booleanTest::compute(C)) {        if (booleanTest::compute(C)) {
2577              ref=2;          ref=2;
2578          } else {        } else {
2579              ref=0;          ref=0;
2580          }        }
2581                    
2582        std::list<Triangle*> listTriangle0=(*it0).second;        std::list<Triangle*> listTriangle0=(*it0).second;
2583        for(std::list<Triangle*>::iterator jt0=listTriangle0.begin() ;        for(std::list<Triangle*>::iterator jt0=listTriangle0.begin() ;
# Line 2248  __createTrianglesIntersection(const size Line 2585  __createTrianglesIntersection(const size
2585    
2586          Triangle& T = (*(*jt0));          Triangle& T = (*(*jt0));
2587          for (size_t n=0; n<3; ++n) {          for (size_t n=0; n<3; ++n) {
2588              // if(ref==2)            // if(ref==2)
2589             // ffout(4)<<"ref ="<<ref<<"\n";            // ffout(4)<<"ref ="<<ref<<"\n";
2590            (*edgesRefVertex[numobject])[&T(n)]=ref;            (*edgesRefVertex[numobject])[&T(n)]=ref;
2591              //put ref 2 for points of surfmesh            //put ref 2 for points of surfmesh
2592              const Vertex* P0=(*listOfVertexMesh.find(&T(n))).second;            const Vertex* P0=(*listOfVertexMesh.find(&T(n))).second;
2593              (*edgesRefVertex[numobject])[P0]=ref;            (*edgesRefVertex[numobject])[P0]=ref;
2594          }          }
2595    
2596          if((*edgesRefVertex[numobject])[&T(0)]==2          if((*edgesRefVertex[numobject])[&T(0)]==2
# Line 2287  __createTrianglesIntersection(const size Line 2624  __createTrianglesIntersection(const size
2624          K.reference()=(*T).reference();          K.reference()=(*T).reference();
2625          triangleWithVertex.clear();          triangleWithVertex.clear();
2626          vertexInTriangles.clear();          vertexInTriangles.clear();
2627            vertexVectTriangles.clear();
2628                    
2629          __createLocalListIntersection(objectNumber,O2,T,K,*i);          __createLocalListIntersection(objectNumber,O2,T,K,*i);
2630          std::map< const Vertex* ,const Vertex*> pointsin=K.pointsin;          std::map< const Vertex* ,const Vertex*> pointsin=K.pointsin;
2631          size_t numberP=pointsin.size();  
           
2632          const Vertex* P0=(*listOfVertexMesh.find(&K(0))).second;          const Vertex* P0=(*listOfVertexMesh.find(&K(0))).second;
2633          const Vertex* P1=(*listOfVertexMesh.find(&K(1))).second;          const Vertex* P1=(*listOfVertexMesh.find(&K(1))).second;
2634          const Vertex* P2=(*listOfVertexMesh.find(&K(2))).second;          const Vertex* P2=(*listOfVertexMesh.find(&K(2))).second;
# Line 2299  __createTrianglesIntersection(const size Line 2636  __createTrianglesIntersection(const size
2636          (*edgesRefVertex[numobject])[&K(1)]=(*edgesRefVertex[numobject])[P1];          (*edgesRefVertex[numobject])[&K(1)]=(*edgesRefVertex[numobject])[P1];
2637          (*edgesRefVertex[numobject])[&K(2)]=(*edgesRefVertex[numobject])[P2];          (*edgesRefVertex[numobject])[&K(2)]=(*edgesRefVertex[numobject])[P2];
2638          for(size_t k=0 ; k<3 ; ++k) {          for(size_t k=0 ; k<3 ; ++k) {
2639          if((*edgesRefVertex[numobject])[&K(k)]==2) {            if((*edgesRefVertex[numobject])[&K(k)]==2) {
2640              K.isInTriangle[k]=1;              K.isInTriangle[k]=1;
2641          } else {            } else {
2642              K.isInTriangle[k]=0;              K.isInTriangle[k]=0;
2643          }            }
2644          }          }
2645          if( K.isInTriangle[0]+  K.isInTriangle[1]+  K.isInTriangle[2]!=0 and          if( K.isInTriangle[0]+  K.isInTriangle[1]+  K.isInTriangle[2]!=0 and
2646       K.isInTriangle[0]+  K.isInTriangle[1]+  K.isInTriangle[2]!=3              K.isInTriangle[0]+  K.isInTriangle[1]+  K.isInTriangle[2]!=3
2647       and  K.edgecut[0]+ K.edgecut[1]+ K.edgecut[2]==0) {              and  K.edgecut[0]+ K.edgecut[1]+ K.edgecut[2]==0) {
2648              ffout(4)<<"cas pbssssss\n";            ffout(4)<<"cas pbssssss\n";
2649              ffout(4)<<(*edgesRefVertex[numobject])[&K(0)]<<" "<<(*edgesRefVertex[numobject])[&K(1)]<<            ffout(4)<<(*edgesRefVertex[numobject])[&K(0)]<<" "<<(*edgesRefVertex[numobject])[&K(1)]<<
2650                  " "<<(*edgesRefVertex[numobject])[&K(2)]<<" "<<K.edgecut<<"\n";              " "<<(*edgesRefVertex[numobject])[&K(2)]<<" "<<K.edgecut<<"\n";
2651              // __addTriangle(Triangle(K(0),K(1),K(2),K.reference()),*i,triangleListIntersectionNew);
2652    
2653          }          }
2654          if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]>=6) {          if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]>=6) {
2655              ffout(4)<<"on veut rien faire car triangle plat\n";            ffout(4)<<"on veut rien faire car triangle plat\n";
2656              ffout(4)<< K.isInTriangle[0]<<" "<< K.isInTriangle[1]<<" "<< K.isInTriangle[2]<<"\n";            ffout(4)<<"K.isInTriangle "<< K.isInTriangle[0]<<" "<< K.isInTriangle[1]<<" "<< K.isInTriangle[2]<<"\n";
2657              K.isInTriangle[0]=1;            ffout(4)<<K(0)<<"\n";
2658               K.isInTriangle[1]=1;            ffout(4)<<K(1)<<"\n";
2659                K.isInTriangle[2]=1;            ffout(4)<<K(2)<<"\n";
2660              ffout(4)<<"edge cut "<<K.edgecut<<"\n";
2661              K.isInTriangle[0]=1;
2662              K.isInTriangle[1]=1;
2663              K.isInTriangle[2]=1;
2664                            
2665          }          }
2666  ///////////////////////creation des triangles//////////////////////////////          ///////////////////////creation des triangles//////////////////////////////
2667          size_t numberVertexInDomain=K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2];          size_t numberVertexInDomain=K.isInTriangle[0]+K.isInTriangle[1]+K.isInTriangle[2];
2668          switch (numberVertexInDomain) {          switch (numberVertexInDomain) {
2669          case 0: {          case 0: {
2670            if(K.edgecut[0]+K.edgecut[1]+K.edgecut[2]==3 and vertexInTriangles.size()<3) {            if(K.edgecut[0]+K.edgecut[1]+K.edgecut[2]==3 and vertexVectTriangles.size()<3) {
2671              ffout(4)<<"tous les points d'inter ne sont pas dans la liste\n";              ffout(4)<<"tous les points d'inter ne sont pas dans la liste\n";
2672              //__casePinterNotInList(K,*i, triangleListIntersectionNew);            } else {
2673              //std::exit(1);              ffout(4)<<"------ case 0 in\n";
2674            } else if(K.edgecut[0]+K.edgecut[1]+K.edgecut[2]!=0) {              __createGeneral(K,*i,triangleListIntersectionNew);
             if(vertexInTriangles.size()==3) {  
               //there are 3 points intersection so we create a triangle  
                 //__3pointsintersection(K,*i,triangleListIntersectionNew);  
             } else {  
                 size_t numberEdgesCut=K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0];  
               if(vertexInTriangles.size()>3 and numberEdgesCut==1){  
                 //one edge cut but twice  
                 ffout(4)<<"createNotIn\n";  
                 //POUR STEPH : createNotIn peut generer des tr plats  
                 //__createNotIn(K,*i, triangleListIntersectionNew);  
               } else {  
                 if(vertexInTriangles.size()>3 and numberEdgesCut==2) {  
                   //two edge cut but twice  
                   ffout(4)<<"createNotIn2\n";  
                   //__createNotIn2(K,*i, triangleListIntersectionNew);  
                 } else {  
                   if(K.edgecut[0]!=0) {  
                     if(K.pointsIntersection[0][0]==&V0  
                        or K.pointsIntersection[0][0]==&V1  
                        or K.pointsIntersection[0][0]==&V2 or vertexInTriangles.size()==1) {  
                       // ffout(4)<<"**********normal\n";  
                     } else {  
                       if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==3  
                          and K.edgecut[0][1]+K.edgecut[1][1]+K.edgecut[2][1]==3) {  
                         if(K.pointsIntersection[0][1]==K.pointsIntersection[0][0]) {  
                           K.edgecut[0][0]=0;  
                           K.edgecut[0][1]=0;  
                         }  
                         if(K.pointsIntersection[1][1]==K.pointsIntersection[1][0]) {  
                           K.edgecut[1][0]=0;  
                           K.edgecut[1][1]=0;  
                         }  
                         if(K.pointsIntersection[2][1]==K.pointsIntersection[2][0]) {  
                           K.edgecut[2][0]=0;  
                           K.edgecut[2][1]=0;  
                         }  
                         if(vertexInTriangles.size()>3 and  
                            K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==2){  
                            ffout(4)<<"hhhhhhhh\n";  
                             ffout(4)<<K(0)<<"\n";  
                             ffout(4)<<K(1)<<"\n";  
                             ffout(4)<<K(2)<<"\n";  
                           //__createNotIn2(K,*i,triangleListIntersectionNew);  
                         } else {  
                           if(vertexInTriangles.size()>3 and  
                              K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==1) {  
                             ffout(4)<<"case\n";  
                             //__createNotIn(K,*i,triangleListIntersectionNew);  
                           }  
                         }  
                       } else {  
                         if(vertexInTriangles.size()>=3) {  
                           ffout(4)<<"triangle cut but outside1\n";  
                           ffout(4)<<V0<<"\n";  
                           ffout(4)<<V1<<"\n";  
                           ffout(4)<<V2<<"\n";  
                           ffout(4)<<"nb points inter "<<vertexInTriangles.size()<<"\n";  
                         }  
                       }  
                     }  
                   } else {  
                     if(K.edgecut[1]!=0) {  
                       if(K.pointsIntersection[1][0]==&V0  
                          or K.pointsIntersection[1][0]==&V1  
                          or K.pointsIntersection[1][0]==&V2 or vertexInTriangles.size()<3) {  
                         //  ffout(4)<<"**********normal\n";  
                       } else {  
                         ffout(4)<<"triangle cut but outside2\n";  
                         ffout(4)<<V0<<"\n";  
                         ffout(4)<<V1<<"\n";  
                         ffout(4)<<V2<<"\n";  
                         ffout(4)<<"nb points inter "<<vertexInTriangles.size()<<"\n";  
                       }  
                     } else {  
                       if(K.edgecut[2]!=0) {  
                         if(K.pointsIntersection[2][0]==&V0  
                            or K.pointsIntersection[2][0]==&V1  
                            or K.pointsIntersection[2][0]==&V2 or vertexInTriangles.size()<3) {  
                           //ffout(4)<<"**********normal\n";  
                         } else {  
                           ffout(4)<<"triangle cut but outside3\n";  
                           ffout(4)<<V0<<"\n";  
                           ffout(4)<<V1<<"\n";  
                           ffout(4)<<V2<<"\n";  
                           ffout(4)<<"nb points inter "<<vertexInTriangles.size()<<"\n";  
                         }  
                       }  
                     }  
                   }  
                 }  
               }  
             }  
2675            }              }  
2676            break;              break;  
2677          }          }
2678          case 1: {          case 1: {
2679              if(K.edgecut[0]+K.edgecut[1]+K.edgecut[2]==0) {            if(K.edgecut[0]+K.edgecut[1]+K.edgecut[2]==0) {
2680                  ffout(4)<<"pbbbbbbbbbbbbbbb\n";              ffout(4)<<"pbbbbbbbbbbbbbbb\n";
2681              }            }
2682            if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==3) {            if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==3) {
2683              if(K.pointsIntersection[0][0]!=K.pointsIntersection[1][0]              if(K.pointsIntersection[0][0]!=K.pointsIntersection[1][0]
2684                 and K.pointsIntersection[1][0]!=K.pointsIntersection[2][0]                 and K.pointsIntersection[1][0]!=K.pointsIntersection[2][0]
2685                  and K.pointsIntersection[2][0]!=K.pointsIntersection[0][0]) {                 and K.pointsIntersection[2][0]!=K.pointsIntersection[0][0]) {
2686                   ffout(4)<<"---------------- 3 edge cut \n";                ffout(4)<<"---------------- 3 edge cut \n";
2687                   //__createGeneral(K,*i,triangleListIntersectionNew);                __createGeneral(K,*i,triangleListIntersectionNew);
2688              } else {              } else {
2689                  if(numberP<=2) {                ffout(4)<<"1------------\n";
2690                   __createGeneral(K, *i, triangleListIntersectionNew);                ffout(4)<<K.edgecut[0]<<" "<<K.edgecut[1]<<" "<<K.edgecut[2]<<"\n";
2691                  } else {                __createGeneral(K, *i, triangleListIntersectionNew);
                     ffout(4)<<"numberP "<<numberP<<"\n";  
                 }  
2692              }              }
2693            } else {            } else {
2694                if(numberP<=2) {              ffout(4)<<"2------------\n";
2695                    __createGeneral(K, *i, triangleListIntersectionNew);              ffout(4)<<K.edgecut[0]<<" "<<K.edgecut[1]<<" "<<K.edgecut[2]<<"\n";
2696                } else {              __createGeneral(K, *i, triangleListIntersectionNew);
                   ffout(4)<<"numberP "<<numberP<<"\n";  
               }  
2697            }            }
2698            break;            break;
2699          }          }
2700          case 2: {          case 2: {
2701              if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==3) {            if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]<2) {
2702                  ffout(4)<<"case 2 pas traite \n";              ffout(4)<<"pbs pas assez d'edge cut \n";
2703              } else {            } else {
2704                  if(numberP<=2) {              ffout(4)<<"---------- case 2\n";
2705                      __createGeneral(K, *i, triangleListIntersectionNew);              ffout(4)<<K.edgecut[0]<<" "<<K.edgecut[1]<<" "<<K.edgecut[2]<<"\n";
2706                  } else {              __createGeneral(K, *i, triangleListIntersectionNew);
2707                      ffout(4)<<"numberP "<<numberP<<"\n";            }
                 }  
             }  
2708            break;            break;
2709          }          }
2710          case 3: {          case 3: {
2711              if(K.edgecut[0]+K.edgecut[1]+K.edgecut[2]==0) {            if(K.edgecut[0]+K.edgecut[1]+K.edgecut[2]==0) {
2712                __addTriangle(*T,*i,triangleListIntersectionNew);              __addTriangle(*T,*i,triangleListIntersectionNew);
2713              } else {
2714                if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==3) {
2715                    __createGeneral(K, *i, triangleListIntersectionNew);
2716                  ffout(4)<<"cas a traiter-------------\n";
2717              } else {              } else {
2718                  if(K.edgecut[0][0]+K.edgecut[1][0]+K.edgecut[2][0]==3) {                  ffout(4)<<"cas pas traite\n";
2719                      ffout(4)<<"cas a traiter-------------\n";                  ffout(4)<<K.edgecut[0][0]<<" "<<K.edgecut[1][0]<<" "<<K.edgecut[2][0]<<"\n";
2720                  } else {                  __createGeneral(K, *i, triangleListIntersectionNew);
                     ffout(4)<<K.edgecut[0][0]<<" "<<K.edgecut[1][0]<<" "<<K.edgecut[2][0]<<"\n";  
                     ffout(4)<<"cas pas traite\n";  
                 }  
2721              }              }
2722              }
2723            break;            break;
2724          }          }
2725          }          }
2726  //////////////////////////////////////////////////////////////////////////          //////////////////////////////////////////////////////////////////////////
2727  //      // K.isInTriangle[i] vaut 1 si le point est a garder pour construire le maillage 0 sinon.          //      // K.isInTriangle[i] vaut 1 si le point est a garder pour construire le maillage 0 sinon.
2728  //      // K.edgecut[i] vaut 1 si l'arete est coupee          //      // K.edgecut[i] vaut 1 si l'arete est coupee
2729    
2730          }
2731      }      }
2732    }    }
2733  }    (*edgesRefVertex[numobject]).clear();
2734      (*edgesRefVertex[numobject]).clear();    listOfVertexMesh.clear();
2735      listOfVertexMesh.clear();    listOfTriangleMeshFront.clear();
     listOfTriangleMeshFront.clear();  
2736  }  }
2737    
2738  bool SurfaceMeshGenerator::Internals::  bool SurfaceMeshGenerator::Internals::
# Line 2838  __createSurface(Structured3DMesh& SMesh, Line 3084  __createSurface(Structured3DMesh& SMesh,
3084                                                    currentCell,                                                    currentCell,
3085                                                    cellList);                                                    cellList);
3086              if(!(jobDone2)) {              if(!(jobDone2)) {
3087                  jobDone=jobDone2;                jobDone=jobDone2;
3088              }              }
3089            }            }
3090                            
# Line 2853  __createSurface(Structured3DMesh& SMesh, Line 3099  __createSurface(Structured3DMesh& SMesh,
3099                ++itcell;                ++itcell;
3100              }              }
3101            } else {            } else {
3102                ffout(4)<<"------jobDone false\n";              ffout(4)<<"------jobDone false\n";
3103              SurfaceMeshGenerator::Internals::MotherCellList::iterator itcell              SurfaceMeshGenerator::Internals::MotherCellList::iterator itcell
3104                =cellList.begin();                =cellList.begin();
3105              for (EdgePairList::iterator t = localTriangleList.begin();              for (EdgePairList::iterator t = localTriangleList.begin();
# Line 2961  generatesAllSufaceMeshes(const Problem& Line 3207  generatesAllSufaceMeshes(const Problem&
3207            Structured3DMesh& SMesh = dynamic_cast<Structured3DMesh&>(M);            Structured3DMesh& SMesh = dynamic_cast<Structured3DMesh&>(M);
3208            if(bs.structMesh(k)) {            if(bs.structMesh(k)) {
3209              bcMeshAssociation.addPDEBoundaryConditionAndMesh(i,&BC,&SMesh.surfaceMeshOfQuadrangles(k));              bcMeshAssociation.addPDEBoundaryConditionAndMesh(i,&BC,&SMesh.surfaceMeshOfQuadrangles(k));
3210              }            }
3211          }          }
3212          break;          break;
3213        }        }
# Line 3251  generateSurfacicMesh(const Domain& omega Line 3497  generateSurfacicMesh(const Domain& omega
3497    (*__internals).__constructionFinalMesh((*objectStack.top()),    (*__internals).__constructionFinalMesh((*objectStack.top()),
3498                                           s_mesh);                                           s_mesh);
3499    
3500  /*  #warning _________________________________A ENLEVER_________  #warning _________________________________A ENLEVER_________
3501    Connectivity<Triangle> connect(s_mesh.numberOfCells());    Connectivity<Triangle> connect(s_mesh.numberOfCells());
3502    ConnectivityBuilder<SurfaceMeshOfTriangles> C(s_mesh,connect);    ConnectivityBuilder<SurfaceMeshOfTriangles> C(s_mesh,connect);
3503    size_t compt=0;    size_t compt=0;
3504    for(size_t i=0 ; i<s_mesh.numberOfCells() ; ++i) {    for(size_t i=0 ; i<s_mesh.numberOfCells() ; ++i) {
3505        const Triangle& T=s_mesh.cell(i);      const Triangle& T=s_mesh.cell(i);
3506        TinyVector<3,const Triangle*> V=connect(s_mesh.cellNumber(T));      TinyVector<3,const Triangle*> V=connect(s_mesh.cellNumber(T));
3507        for(size_t j=0 ; j<3 ; ++j) {      for(size_t j=0 ; j<3 ; ++j) {
3508            if(V[j]==0) {        if(V[j]==0) {
3509                //ffout(4)<<"attention un des voisins n'existe pas\n";          //ffout(4)<<"attention un des voisins n'existe pas\n";
3510                ++compt;          ++compt;
           }  
3511        }        }
3512        }
3513    }    }
3514    ffout(4)<<"****************************il manque "<<compt<<" sommets\n";*/    ffout(4)<<"****************************il manque "<<compt<<" sommets\n";
3515    /////////////////////////////////////////////////////////////////////////////    /////////////////////////////////////////////////////////////////////////////
3516    
3517    int n0=0, n1=0;    int n0=0, n1=0;
# Line 3362  void plotmedit(const size_t& nobject, Line 3608  void plotmedit(const size_t& nobject,
3608      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;
3609            
3610      // sortie<<3*i+1+8*toTreatHexahedra.size()<<" "<<3*i+2+8*toTreatHexahedra.size()<<" //"<<3*i+3+8*toTreatHexahedra.size();//<<" "<<1<<"\n";      // sortie<<3*i+1+8*toTreatHexahedra.size()<<" "<<3*i+2+8*toTreatHexahedra.size()<<" //"<<3*i+3+8*toTreatHexahedra.size();//<<" "<<1<<"\n";
3611  //     if(i<n0)      //     if(i<n0)
3612  //       sortie<<" 1"<<"\n";      //       sortie<<" 1"<<"\n";
3613  //     else if(i<n1)      //     else if(i<n1)
3614  //       sortie<<" 2"<<"\n";      //       sortie<<" 2"<<"\n";
3615  //     else      //     else
3616  //       sortie<<" 4"<<"\n";      //       sortie<<" 4"<<"\n";
3617      sortie << " " << T.reference() << '\n';      sortie << " " << T.reference() << '\n';
3618    }    }
3619    ffout(4)<<"fin ecriture triangles\n";    ffout(4)<<"fin ecriture triangles\n";

Legend:
Removed from v.1.51  
changed lines
  Added in v.1.52

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