/[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.7 by dobj, Fri Mar 14 15:25:59 2003 UTC revision 1.8 by dobj, Thu Mar 20 15:47:15 2003 UTC
# Line 88  private: Line 88  private:
88        
89    //! list of cut edges and intersection vertex.    //! list of cut edges and intersection vertex.
90    std::map<Edge::Pair, size_t> edgesRef;    std::map<Edge::Pair, size_t> edgesRef;
91    
92  std::vector<ReferenceCounting<EdgePairList> > triangleListes;  std::vector<ReferenceCounting<EdgePairList> > triangleListes;
93    
94    EdgePairList localTriangleList;    EdgePairList localTriangleList;
# Line 237  __createTriangles(const Shape& S, Line 238  __createTriangles(const Shape& S,
238    
239          //! Stores the Edges on which the vertex will be computed.          //! Stores the Edges on which the vertex will be computed.
240          for (size_t i=0; i<cutEdges.size(); ++i)          for (size_t i=0; i<cutEdges.size(); ++i)
241          //  verticesList[*cutEdges[i]] =  new Vertex(Internals::__splitEdge(*cutEdges[i], ref, S));          //verticesList[*cutEdges[i]] =  new Vertex(Internals::__splitEdge(*cutEdges[i], ref, S));
242              (*verticesListes[objectNumber])[*cutEdges[i]] =  new Vertex(Internals::__splitEdge(*cutEdges[i], ref, S));              (*verticesListes[objectNumber])[*cutEdges[i]] =  new Vertex(Internals::__splitEdge(*cutEdges[i], ref, S));
243    
244          break;          break;
# Line 281  __createTriangles(const Shape& S, Line 282  __createTriangles(const Shape& S,
282            }            }
283    
284          for (size_t i=0; i<4; ++i) {          for (size_t i=0; i<4; ++i) {
285  //        verticesList[*cutEdges[i]] =  new Vertex(__splitEdge(*cutEdges[i], ref, S));            //verticesList[*cutEdges[i]] =  new Vertex(__splitEdge(*cutEdges[i], ref, S));
286              (*verticesListes[objectNumber])[*cutEdges[i]] =  new Vertex(__splitEdge(*cutEdges[i], ref, S));              (*verticesListes[objectNumber])[*cutEdges[i]] =  new Vertex(__splitEdge(*cutEdges[i], ref, S));
287    
288          }          }
# Line 544  void plot2(const SurfaceMeshOfTriangles& Line 545  void plot2(const SurfaceMeshOfTriangles&
545  void plotmedit(const SurfaceMeshOfTriangles& s_mesh,  void plotmedit(const SurfaceMeshOfTriangles& s_mesh,
546            std::set<const Cell*>& toTreatHexahedra,size_t nobj);            std::set<const Cell*>& toTreatHexahedra,size_t nobj);
547  void Bary(const Vertex& P,TinyVector<3>& P0,TinyVector<3>& V12,TinyVector<3>& V02,size_t& ncompt);  void Bary(const Vertex& P,TinyVector<3>& P0,TinyVector<3>& V12,TinyVector<3>& V02,size_t& ncompt);
548    void inout(size_t iter,const size_t nn,TinyVector<3,int>  ttn,size_t& re0,size_t& re1,size_t& re2,size_t& ii0,size_t& ii1,size_t& ii2);
549    
550    /*! Generates a surfacic mesh using only characteristic function of an object    /*! Generates a surfacic mesh using only characteristic function of an object
551      and the Structured 3D Mesh.      and the Structured 3D Mesh.
# Line 599  generateSurfacicMesh(const Object& O, Line 600  generateSurfacicMesh(const Object& O,
600    (*__internals).triangleListes.resize(nobject);    (*__internals).triangleListes.resize(nobject);
601    
602    
603    for(int no=0 ; no<nobject ; ++no)    for(size_t no=0 ; no<nobject ; ++no)
604    {    {
605        (*__internals).triangleListes[no]= new SurfaceMeshGenerator::Internals::EdgePairList;        (*__internals).triangleListes[no]= new SurfaceMeshGenerator::Internals::EdgePairList;
606        (*__internals).verticesListes[no]= new SurfaceMeshGenerator::Internals::VerticesList;        (*__internals).verticesListes[no]= new SurfaceMeshGenerator::Internals::VerticesList;
# Line 625  generateSurfacicMesh(const Object& O, Line 626  generateSurfacicMesh(const Object& O,
626    /// a triangulation has been generated inside the hexahedron    /// a triangulation has been generated inside the hexahedron
627    std::set<const Cell*> meshedHexahedra;    std::set<const Cell*> meshedHexahedra;
628    
629   // if ((*__internals).verticesList.size()==0) {  //  if ((*__internals).verticesList.size()==0) {
630  if ((*(*__internals).verticesListes[0]).size()==0) {  if ((*(*__internals).verticesListes[0]).size()==0) {
631    
632  fferr(2) << __FILE__ << ':' << __LINE__ << ": "  fferr(2) << __FILE__ << ':' << __LINE__ << ": "
# Line 843  fferr(2) << __FILE__ << ':' << __LINE__ Line 844  fferr(2) << __FILE__ << ':' << __LINE__
844                  //(*__internals).verticesList[*cutEdges[i]]                  //(*__internals).verticesList[*cutEdges[i]]
845                    (*(*__internals).verticesListes[0])[*cutEdges[i]]                    (*(*__internals).verticesListes[0])[*cutEdges[i]]
846                    =  new Vertex((*__internals).__splitEdge(*cutEdges[i], ref, S));                    =  new Vertex((*__internals).__splitEdge(*cutEdges[i], ref, S));
847                  (*__internals).edgesRef[*cutEdges[i]] = 1;                    (*__internals).edgesRef[*cutEdges[i]] = 1;
848                }                }
849                break;                break;
850              }              }
# Line 890  fferr(2) << __FILE__ << ':' << __LINE__ Line 891  fferr(2) << __FILE__ << ':' << __LINE__
891                  //(*__internals).verticesList[*cutEdges[i]]                  //(*__internals).verticesList[*cutEdges[i]]
892                    (*(*__internals).verticesListes[0])[*cutEdges[i]]                    (*(*__internals).verticesListes[0])[*cutEdges[i]]
893                    =  new Vertex((*__internals).__splitEdge(*cutEdges[i], ref, S));                    =  new Vertex((*__internals).__splitEdge(*cutEdges[i], ref, S));
894                  (*__internals).edgesRef[*cutEdges[i]] = 1;                          (*__internals).edgesRef[*cutEdges[i]] = 1;
895                }                }
896    
897                //! Addes a triangle to the list.                //! Addes a triangle to the list.
# Line 915  fferr(2) << __FILE__ << ':' << __LINE__ Line 916  fferr(2) << __FILE__ << ':' << __LINE__
916              }              }
917              case 1:              case 1:
918              case 2: {              case 2: {
919                jobDone = false;                  jobDone = false;
920                break;                break;
921              }              }
922              }              }
# Line 927  fferr(2) << __FILE__ << ':' << __LINE__ Line 928  fferr(2) << __FILE__ << ':' << __LINE__
928                     = (*__internals).localTriangleList.begin();                     = (*__internals).localTriangleList.begin();
929                   t != (*__internals).localTriangleList.end(); ++t) {                   t != (*__internals).localTriangleList.end(); ++t) {
930                  (*(*__internals).triangleListes[0]).push_front(*t);                  (*(*__internals).triangleListes[0]).push_front(*t);
931    
932              }              }
933                (*__internals).hexaTetraAssociation[&currentCell]
934                    = std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,
935                    ReferenceCounting<std::vector<Edge> > >
936                    (pT,pAddedEdges);
937                
938            } else {            } else {
939              (*__internals).hexaTetraAssociation[&currentCell]              (*__internals).hexaTetraAssociation[&currentCell]
940                = std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,                = std::pair<ReferenceCounting<std::vector<TetrahedronByEdges> >,
941                            ReferenceCounting<std::vector<Edge> > >                            ReferenceCounting<std::vector<Edge> > >
942                (pT,pAddedEdges);                (pT,pAddedEdges);
943                
944            }            }
945    
946            (*__internals).localTriangleList.clear();            (*__internals).localTriangleList.clear();
# Line 942  fferr(2) << __FILE__ << ':' << __LINE__ Line 950  fferr(2) << __FILE__ << ':' << __LINE__
950  //(*(*__internals).triangleListes[0]).clear();  //(*(*__internals).triangleListes[0]).clear();
951  //(*(*__internals).verticesListes[0]).clear();  //(*(*__internals).verticesListes[0]).clear();
952    
953  //  (*__internals).verticesList.clear();  ////  (*__internals).verticesList.clear();
954    
955    
956  //   (*__internals).__createTriangles(S,ref,0);  //(*__internals).__createTriangles(S,ref,0);
957    
958  //size_t nobj= (*(*__internals).triangleListes[0]).size();  //size_t nobj= (*(*__internals).triangleListes[0]).size();
959  //std::cout<<"nb tr objet1 "<<nobj<<" \n";  //std::cout<<"nb tr objet1 "<<nobj<<" \n";
960    
961  //int num=0;  //int num=0;
962  //   for(ObjectReferences::iterator i = otherReferences.begin();   //  for(ObjectReferences::iterator i = otherReferences.begin();
963   //      i != otherReferences.end(); ++i) {    //     i != otherReferences.end(); ++i) {
964  //       num++;     //    num++;
965  //     for(size_t n = 0; n<scene.nbObjects(); ++n) {      // for(size_t n = 0; n<scene.nbObjects(); ++n) {
966  //       Reference::Representation otherRef;       //  Reference::Representation otherRef;
967  //       otherRef.resize(scene.nbObjects());        // otherRef.resize(scene.nbObjects());
968  //  
969  //       if (scene.object(n).reference() == ((*i).first)) {   //      if (scene.object(n).reference() == ((*i).first)) {
970  //      otherRef[n] = true;   //     otherRef[n] = true;
971  //      (*__internals).__createTriangles(*(scene.object(n).shape()),   //     (*__internals).__createTriangles(*(scene.object(n).shape()),
972  //                                       otherRef,num);   //                                      otherRef,num);
973   //      }     //    }
974  //     }      // }
975  //  }    // }
976  //size_t nobjother=(*(*__internals).triangleListes[1]).size();  //size_t nobjother=(*(*__internals).triangleListes[1]).size();
977  //std::cout<<"nb tr objet2 "<<nobjother<<" \n";  //std::cout<<"nb tr objet2 "<<nobjother<<" \n";
978    
979    /*  // *******************************************************************************
980    //std::map<Edge::Pair, Vertex*>  listofVertex;
981    //std::map<Edge::Pair,int> triangleedge;
982    std::vector< TinyVector<3,int> > tt;
983    std::vector< std::vector<Vertex> > vv;
984       vv.resize((*(*__internals).triangleListes[0]).size());
985       tt.resize((*(*__internals).triangleListes[0]).size());
986       for(size_t kk=0 ;kk<nobj;kk++)
987           vv[kk].resize(3);
988       int nn=0;
989      
990      for(SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle1
991          = (*(*__internals).triangleListes[0]).begin();
992          currentTriangle1 != (*(*__internals).triangleListes[0]).end();
993          ++currentTriangle1)
994      {
995          const Edge::Pair e0=(*currentTriangle1)[0];
996          const Edge::Pair e1=(*currentTriangle1)[1];
997          const Edge::Pair e2=(*currentTriangle1)[2];
998          const Vertex& V0 = *((*(*__internals).verticesListes[0])[e0]);
999          const Vertex& V1 = *((*(*__internals).verticesListes[0])[e1]);
1000          const Vertex& V2 = *((*(*__internals).verticesListes[0])[e2]);
1001    
1002          //equation of the plan
1003          TinyVector<3>  V01,V12,V02;
1004          V01=V1.Vector3()-V0.Vector3();
1005          V12=V2.Vector3()-V1.Vector3();
1006          V02=V2.Vector3()-V0.Vector3();
1007        //  const TinyVector<3> coeff=V01^V02;
1008        //  const real_t d=-(coeff[0]*V0[0]+coeff[1]*V0[1]+coeff[2]*V0[2]);  
1009    //      std::cout<<"--------------------------------------------------------------\n";
1010    
1011          tt[nn][0]=0;
1012          tt[nn][1]=0;
1013          tt[nn][2]=0;
1014          for(SurfaceMeshGenerator::Internals::EdgePairList::iterator currentTriangle
1015              = (*(*__internals).triangleListes[1]).begin();
1016              currentTriangle != (*(*__internals).triangleListes[1]).end();
1017              ++currentTriangle)
1018          {
1019              //std::cout<<"*************************\n";
1020              const Vertex& Vo0 = *((*(*__internals).verticesListes[1])[(*currentTriangle)[0]]);
1021              const Vertex& Vo1 = *((*(*__internals).verticesListes[1])[(*currentTriangle)[1]]);
1022              const Vertex& Vo2 = *((*(*__internals).verticesListes[1])[(*currentTriangle)[2]]);
1023    
1024              //equation of plan
1025              TinyVector<3>  Vo01,Vo12,Vo02;
1026              Vo01=Vo1.Vector3()-Vo0.Vector3();
1027              Vo12=Vo2.Vector3()-Vo1.Vector3();
1028              Vo02=Vo2.Vector3()-Vo0.Vector3();
1029              
1030              TinyVector<3> coeffo;
1031              coeffo=Vo01^Vo02;
1032              real_t d2=-(coeffo[0]*Vo0[0]+coeffo[1]*Vo0[1]+coeffo[2]*Vo0[2]);
1033    //std::cout<<"verif "<<coeffo[0]*Vo0.Vector3()[0]+coeffo[1]*Vo0.Vector3()[1]+coeffo[2]*Vo0.Vector3()[2]+d2<<"\n";
1034    //std::cout<<"verif "<<coeffo[0]*Vo1.Vector3()[0]+coeffo[1]*Vo1.Vector3()[1]+coeffo[2]*Vo1.Vector3()[2]+d2<<"\n";
1035    //std::cout<<"verif "<<coeffo[0]*Vo2.Vector3()[0]+coeffo[1]*Vo2.Vector3()[1]+coeffo[2]*Vo2.Vector3()[2]+d2<<"\n";
1036              real_t pscalar,k;
1037              size_t num;
1038              {
1039                  Vertex P;
1040                  real_t norm1=Norm(V01);
1041                  pscalar=coeffo*V01;
1042                  num=0;
1043                  if(pscalar>1e-13 or pscalar<-1e-13)
1044                  {
1045                      k=-(coeffo[0]*V0[0]+coeffo[1]*V0[1]+coeffo[2]*V0[2]+d2)/pscalar;
1046                      P=V0+k*V01;
1047                      TinyVector<3> P0=P.Vector3()-Vo2.Vector3(),P1=Vo12*(-1.),P2=Vo02*(-1.);
1048                      real_t norm2=Norm(P.Vector3()-V0.Vector3());
1049                      real_t norm3=Norm(P.Vector3()-V1.Vector3());
1050                      if(norm2<=norm1 and norm3<=norm1)
1051                      {
1052                          Bary(P,P0,P1,P2,num);
1053                          if(num==3)
1054                          {
1055    //std::cout<<"droite1 "<<P.Vector3()[0]<<" "<<P.Vector3()[1]<<" "<<P.Vector3()[2]<<"  \n";
1056                              tt[nn][0]=1;
1057                              vv[nn][0]=P;
1058    //std::cout<<"P "<<nn<<" "<< vv[0][0].Vector3()[0]<<" "<<(vv[0][0]).Vector3()[1]<<" "<<(vv[0][0]).Vector3()[2]<<"\n";
1059    //std::cout<<"verif "<<coeffo[0]*P.Vector3()[0]+coeffo[1]*P.Vector3()[1]+coeffo[2]*P.Vector3()[2]+d2<<"\n";
1060    //std::cout<<"verif "<<coeff[0]*P.Vector3()[0]+coeff[1]*P.Vector3()[1]+coeff[2]*P.Vector3()[2]+d<<"\n";
1061                          }
1062                      }
1063                  }
1064              }
1065              {
1066                  Vertex P;
1067                  pscalar=coeffo*V02;
1068                  real_t norm1=Norm(V02);
1069                  if(pscalar>1e-13 or pscalar<-1e-13)
1070                  {
1071                      k=-(coeffo[0]*V0[0]+coeffo[1]*V0[1]+coeffo[2]*V0[2]+d2)/pscalar;
1072                      P=V0+k*V02;
1073                      TinyVector<3> P0=P.Vector3()-Vo2.Vector3(),P1=Vo12*(-1.),P2=Vo02*(-1.);
1074                      real_t norm2=Norm(P.Vector3()-V0.Vector3());
1075                      real_t norm3=Norm(P.Vector3()-V2.Vector3());
1076                      //std::cout<<"norm "<<norm3<<" "<<norm2<<" "<<norm1<<"\n";
1077                      if(norm2<=norm1 and norm3<=norm1)
1078                      {
1079                          Bary(P,P0,P1,P2,num);
1080                          if(num==3)
1081                          {
1082    //std::cout<<"droite2 "<<P.Vector3()[0]<<" "<<P.Vector3()[1]<<" "<<P.Vector3()[2]<<"  \n";
1083                              tt[nn][1]=1;
1084                              (vv[nn])[1]=P;
1085                          }
1086                      }
1087                  }
1088              }
1089              {
1090                  Vertex P;
1091                  pscalar=coeffo*V12;
1092                  real_t norm1=Norm(V12);
1093                  if(pscalar>1e-13 or pscalar<-1e-13)
1094                  {
1095                      k=-(coeffo[0]*V1[0]+coeffo[1]*V1[1]+coeffo[2]*V1[2]+d2)/pscalar;
1096                      P=V1+k*V12;
1097                      TinyVector<3> P0=P.Vector3()-Vo2.Vector3(),P1=Vo12*(-1.),P2=Vo02*(-1.);
1098                      real_t norm2=Norm(P.Vector3()-V1.Vector3());
1099                      real_t norm3=Norm(P.Vector3()-V2.Vector3());
1100                      
1101                      if(norm2<=norm1 and norm3<=norm1)
1102                      {
1103                          Bary(P,P0,P1,P2,num);
1104                          if(num==3)
1105                          {
1106                              tt[nn][2]=1;
1107                              (vv[nn])[2]=P;
1108    
1109                          }
1110    
1111                      }
1112                  }
1113              }
1114          }
1115          nn++;
1116      }
1117      //std::cout<<"Papres "<<nn<<" "<<nobj<<" "<<(vv[0][0]).Vector3()[0]<<" "<<(vv[0][0]).Vector3()[1]<<" "<<(vv[0][0]).Vector3()[2]<<"\n";
1118    
1119      // create triangles
1120    
1121    std::vector<Vertex> vertexplot;
1122    {    nn=0;
1123        size_t ii0,ii1,ii2;
1124        for(SurfaceMeshGenerator::Internals::EdgePairList::const_iterator currentTriangle
1125            = (*(*__internals).triangleListes[0]).begin();
1126            currentTriangle != (*(*__internals).triangleListes[0]).end();
1127            ++currentTriangle)
1128        {
1129            const Edge& e0=(*currentTriangle)[0];
1130            const Edge& e1=(*currentTriangle)[1];
1131            const Edge& e2=(*currentTriangle)[2];
1132    //      std::cout<<"ref "<<nn<<" "<<(*__internals).edgesRef[e0]<<" "<<(*__internals).edgesRef[e1]<<" "<<(*__internals).edgesRef[e2]<<"\n";
1133    //std::cout<<"tt "<<tt[nn][0]<<" "<<tt[nn][1]<<" "<<tt[nn][2]<<"\n";
1134    
1135    
1136                   inout(1,nn,tt[nn],(*__internals).edgesRef[e0],(*__internals).edgesRef[e1],(*__internals).edgesRef[e2],ii0,ii1,ii2);
1137    
1138    
1139    //std::cout<<"ii "<<ii0<<" "<<ii1<<" "<<ii2<<"\n";
1140    //             std::cout<<"apres ref "<<(*__internals).edgesRef[e0]<<" "<<(*__internals).edgesRef[e1]<<" "<<(*__internals).edgesRef[e2]<<"\n";
1141                   nn++;
1142        }
1143        nn=0;
1144            for(SurfaceMeshGenerator::Internals::EdgePairList::const_iterator currentTriangle
1145         = (*(*__internals).triangleListes[0]).begin();
1146         currentTriangle != (*(*__internals).triangleListes[0]).end();
1147         ++currentTriangle)
1148            {
1149                const Edge& e0=(*currentTriangle)[0];
1150                const Edge& e1=(*currentTriangle)[1];
1151                const Edge& e2=(*currentTriangle)[2];
1152    
1153    
1154                inout(2,nn,tt[nn],(*__internals).edgesRef[e0],(*__internals).edgesRef[e1],(*__internals).edgesRef[e2],ii0,ii1,ii2);
1155    
1156    
1157                nn++;
1158            }
1159            nn=0;
1160            for(SurfaceMeshGenerator::Internals::EdgePairList::const_iterator currentTriangle
1161         = (*(*__internals).triangleListes[0]).begin();
1162         currentTriangle != (*(*__internals).triangleListes[0]).end();
1163         ++currentTriangle)
1164            {
1165                const Edge& e0=(*currentTriangle)[0];
1166                const Edge& e1=(*currentTriangle)[1];
1167                const Edge& e2=(*currentTriangle)[2];
1168    
1169                inout(2,nn,tt[nn],(*__internals).edgesRef[e0],(*__internals).edgesRef[e1],(*__internals).edgesRef[e2],ii0,ii1,ii2);
1170    
1171    //std::cout<<ii0<<" "<<ii1<<" "<<ii2<<"\n";
1172                
1173                nn++;
1174            }
1175            
1176        nn=0;
1177        for(SurfaceMeshGenerator::Internals::EdgePairList::const_iterator currentTriangle
1178            = (*(*__internals).triangleListes[0]).begin();
1179            currentTriangle != (*(*__internals).triangleListes[0]).end();
1180            ++currentTriangle)
1181        {
1182            //if the triangle is intersected
1183    //std::cout<<"**************************************************"<<"\n";
1184             const Edge& e0=(*currentTriangle)[0];
1185             const Edge& e1=(*currentTriangle)[1];
1186           const Edge& e2=(*currentTriangle)[2];
1187    //std::cout<<"  avt i123 "<<nn<<" "<<ii0<<" "<<ii1<<" "<<ii2<<" \n";
1188           if(tt[nn][0]+tt[nn][1]+tt[nn][2]!=0)
1189           {
1190    //std::cout<<"-------------------------------------------------"<<"\n";
1191    //std::cout<<"tt "<<nn<<" "<<tt[nn][0]<<" "<<tt[nn][1]<<" "<<tt[nn][2]<<"\n";
1192           }
1193    
1194    //std::cout<<"ref "<<(*__internals).edgesRef[e0]<<" "<<(*__internals).edgesRef[e1]<<" "<<(*__internals).edgesRef[e2]<<"\n";
1195    //std::cout<<"inter "<<nn<<" "<<(vv[0][0]).Vector3()[0]<<" "<<(vv[0][0]).Vector3()[1]<<" "<<(vv[0][0]).Vector3()[2]<<" \n";
1196             // if(triangleedge[(*currentTriangle)[0]]+triangleedge[(*currentTriangle)[1]]+ triangleedge[(*currentTriangle)[2]]==2)
1197    //std::cout<<"avt "<<(*__internals).edgesRef[e0]<<" "<<(*__internals).edgesRef[e1]<<" "<<(*__internals).edgesRef[e2]<<"\n";
1198    
1199      //     if(nn==1254 or nn==1248 or nn==1252 or nn==1255 or nn==1249 )
1200           {
1201    //std::cout<<"-------------------------------------------------"<<"\n";
1202    
1203    //std::cout<<"avt "<<nn<<" "<<(*__internals).edgesRef[e0]<<" "<<(*__internals).edgesRef[e1]<<" "<<(*__internals).edgesRef[e2]<<"\n";
1204           }
1205           inout(3,nn,tt[nn],(*__internals).edgesRef[e0],(*__internals).edgesRef[e1],(*__internals).edgesRef[e2],ii0,ii1,ii2);
1206    
1207        //   if(nn==1254 or nn==1248 or nn==1252 or nn==1255 or nn==1249 )
1208           {
1209    //std::cout<<"apres "<<nn<<" "<<(*__internals).edgesRef[e0]<<" "<<(*__internals).edgesRef[e1]<<" "<<(*__internals).edgesRef[e2]<<"\n";
1210    //std::cout<<"ii "<<ii0<<" "<<ii1<<" "<<ii2<<"\n";
1211           }
1212    //std::cout<<"apres "<<(*__internals).edgesRef[e0]<<" "<<(*__internals).edgesRef[e1]<<" "<<(*__internals).edgesRef[e2]<<"\n";
1213    //std::cout<<" i123 "<<nn<<" "<<ii0<<" "<<ii1<<" "<<ii2<<" \n";
1214    
1215    //std::cout<<e0.Ref()<<" "<<e1.Ref()<<" "<<e2.Ref()<<"\n";
1216            if(tt[nn][0]+tt[nn][1]+tt[nn][2]==2)
1217            {
1218    //std::cout<<"-------------------------------------------------"<<"\n";
1219                Vertex& V0 = *((*(*__internals).verticesListes[0])[(*currentTriangle)[0]]);
1220                Vertex& V1 = *((*(*__internals).verticesListes[0])[(*currentTriangle)[1]]);
1221                Vertex& V2 = *((*(*__internals).verticesListes[0])[(*currentTriangle)[2]]);
1222    //std::cout<<"VO "<<V0.Vector3()[0]<<"  "<<V0.Vector3()[1]<<" "<<V0.Vector3()[2]<<"\n";
1223    //std::cout<<"V1 "<<V1.Vector3()[0]<<"  "<<V1.Vector3()[1]<<" "<<V1.Vector3()[2]<<"\n";
1224    //std::cout<<"V2 "<<V2.Vector3()[0]<<"  "<<V2.Vector3()[1]<<" "<<V2.Vector3()[2]<<"\n";
1225    //std::cout<<" "<<"\n";
1226    //     **************** DELETE************************************************
1227                   bool& i1=V0.Ref().insideDomain();
1228                   bool& i2=V1.Ref().insideDomain();
1229                   bool& i3=V2.Ref().insideDomain();
1230                   if(ii0==1)
1231                       i1=true;
1232                   else
1233                       i1=false;
1234                   if(ii1==1)
1235                       i2=true;
1236                   else
1237                       i2=false;
1238                   if(ii2==1)
1239                       i3=true;
1240                   else
1241                       i3=false;
1242    
1243    //std::cout<<"inside? "<<i1<<" "<<i2<<" "<<i3<<"\n";
1244    //     ****************************************************************
1245    
1246                  //if(triangleedge[(*currentTriangle)[0]]+triangleedge[(*currentTriangle)[1]]==2)
1247                  if(tt[nn][0]+tt[nn][1]==2)
1248                  {
1249                      if(V0.Ref().insideDomain())
1250                      {
1251                          Vertex A=V0;
1252                         // Vertex B=(listofVertex[(*currentTriangle)[0]]);
1253                          //Vertex C=(listofVertex[(*currentTriangle)[1]]);
1254                          Vertex B=(vv[nn][0]);
1255                          Vertex C=(vv[nn][1]);
1256                          vertexplot.push_back(A);
1257                          vertexplot.push_back(B);
1258                          vertexplot.push_back(C);
1259    //std::cout<<"A "<<A.Vector3()[0]<<"  "<<A.Vector3()[1]<<" "<<A.Vector3()[2]<<"\n";
1260    //std::cout<<"B "<<B.Vector3()[0]<<"  "<<B.Vector3()[1]<<" "<<B.Vector3()[2]<<"\n";
1261    //std::cout<<"C "<<C.Vector3()[0]<<"  "<<C.Vector3()[1]<<" "<<C.Vector3()[2]<<"\n";
1262                      }
1263                      else if(V1.Ref().insideDomain() and V2.Ref().insideDomain())
1264                      {
1265                          Vertex A=V1;
1266                          Vertex B=V2;
1267                         // Vertex C=(listofVertex[(*currentTriangle)[0]]);
1268                          //Vertex D=(listofVertex[(*currentTriangle)[1]]);
1269                          Vertex C=(vv[nn][0]);
1270                          Vertex D=(vv[nn][1]);
1271                          vertexplot.push_back(A);
1272                          vertexplot.push_back(C);
1273                          vertexplot.push_back(D);
1274                          vertexplot.push_back(A);
1275                          vertexplot.push_back(B);
1276                          vertexplot.push_back(D);
1277                      }
1278                      else
1279                      {
1280    std::cout<<"erreur1\n";
1281                      }
1282                  }
1283                 // else if(triangleedge[(*currentTriangle)[0]]+triangleedge[(*currentTriangle)[2]]==2)
1284                  else if(tt[nn][0]+tt[nn][2]==2)
1285                  {
1286                      if(V1.Ref().insideDomain())
1287                      {
1288                          Vertex A=V1;
1289                         // Vertex B=(listofVertex[(*currentTriangle)[0]]);
1290                         // Vertex C=(listofVertex[(*currentTriangle)[2]]);
1291                          Vertex B=(vv[nn][0]);
1292                          Vertex C=(vv[nn][2]);
1293                          vertexplot.push_back(A);
1294                          vertexplot.push_back(B);
1295                          vertexplot.push_back(C);
1296    //std::cout<<"A "<<A.Vector3()[0]<<"  "<<A.Vector3()[1]<<" "<<A.Vector3()[2]<<"\n";
1297    //std::cout<<"B "<<B.Vector3()[0]<<"  "<<B.Vector3()[1]<<" "<<B.Vector3()[2]<<"\n";
1298    //std::cout<<"C "<<C.Vector3()[0]<<"  "<<C.Vector3()[1]<<" "<<C.Vector3()[2]<<"\n";
1299                      }
1300                      else if(V0.Ref().insideDomain() and V2.Ref().insideDomain())
1301                      {
1302                          Vertex A=V0;
1303                          Vertex B=V2;
1304                       //   Vertex C=(listofVertex[(*currentTriangle)[0]]);
1305                        //  Vertex D=(listofVertex[(*currentTriangle)[2]]);
1306                          Vertex C=(vv[nn][0]);
1307                          Vertex D=(vv[nn][2]);
1308                          vertexplot.push_back(A);
1309                          vertexplot.push_back(C);
1310                          vertexplot.push_back(D);
1311                          vertexplot.push_back(A);
1312                          vertexplot.push_back(B);
1313                          vertexplot.push_back(D);
1314    
1315                      }
1316                      else
1317                      {
1318    std::cout<<"erreur2\n";
1319                      }
1320                  }
1321                  else
1322                  {
1323                      if(V2.Ref().insideDomain())
1324                      {
1325                          Vertex A=V2;
1326                       //   Vertex B=(listofVertex[(*currentTriangle)[1]]);
1327                         // Vertex C=(listofVertex[(*currentTriangle)[2]]);
1328                          Vertex B=(vv[nn][1]);
1329                          Vertex C=(vv[nn][2]);
1330                          vertexplot.push_back(A);
1331                          vertexplot.push_back(B);
1332                          vertexplot.push_back(C);
1333    //std::cout<<"A "<<A.Vector3()[0]<<"  "<<A.Vector3()[1]<<" "<<A.Vector3()[2]<<"\n";
1334    //std::cout<<"B "<<B.Vector3()[0]<<"  "<<B.Vector3()[1]<<" "<<B.Vector3()[2]<<"\n";
1335    //std::cout<<"C "<<C.Vector3()[0]<<"  "<<C.Vector3()[1]<<" "<<C.Vector3()[2]<<"\n";
1336                      }
1337                      else if(V1.Ref().insideDomain() and V0.Ref().insideDomain())
1338                      {
1339                          Vertex A=V0;
1340                          Vertex B=V1;
1341                      //    Vertex C=(listofVertex[(*currentTriangle)[1]]);
1342                        //  Vertex D=(listofVertex[(*currentTriangle)[2]]);
1343                          Vertex C=(vv[nn][1]);
1344                          Vertex D=(vv[nn][2]);
1345                          vertexplot.push_back(A);
1346                          vertexplot.push_back(C);
1347                          vertexplot.push_back(D);
1348                          vertexplot.push_back(A);
1349                          vertexplot.push_back(B);
1350                          vertexplot.push_back(D);
1351                      }
1352                      else
1353                      {
1354    std::cout<<"erreur3\n";
1355                      }
1356                  }
1357            }
1358            if(tt[nn][0]+tt[nn][1]+tt[nn][2]==3)
1359            {
1360                Vertex& V0 = *((*(*__internals).verticesListes[0])[(*currentTriangle)[0]]);
1361                Vertex& V1 = *((*(*__internals).verticesListes[0])[(*currentTriangle)[1]]);
1362                Vertex& V2 = *((*(*__internals).verticesListes[0])[(*currentTriangle)[2]]);
1363                bool& i1=V0.Ref().insideDomain();
1364                bool& i2=V1.Ref().insideDomain();
1365                bool& i3=V2.Ref().insideDomain();
1366                if(ii0==1)
1367                    i1=true;
1368                else
1369                    i1=false;
1370                if(ii1==1)
1371                    i2=true;
1372                else
1373                    i2=false;
1374                if(ii2==1)
1375                    i3=true;
1376                else
1377                    i3=false;
1378    //std::cout<<"be careful triangle \n";
1379    //std::cout<<"inside? "<<V0.Ref().insideDomain()<<" "<<V1.Ref().insideDomain()<<" "<<V2.Ref().insideDomain()<<"\n";
1380                if(V0.Ref().insideDomain())
1381                {
1382    std::cout<<"be careful triangle 0\n";
1383                    Vertex A=V0;
1384                    Vertex B=(vv[nn][0]);
1385                    Vertex C=(vv[nn][1]);
1386                    Vertex D=(vv[nn][2]);
1387                    vertexplot.push_back(A);
1388                    vertexplot.push_back(B);
1389                    vertexplot.push_back(D);
1390                    vertexplot.push_back(A);
1391                    vertexplot.push_back(C);
1392                    vertexplot.push_back(D);
1393    //std::cout<<"A "<<A.Vector3()[0]<<"  "<<A.Vector3()[1]<<" "<<A.Vector3()[2]<<"\n";
1394    //std::cout<<"B "<<B.Vector3()[0]<<"  "<<B.Vector3()[1]<<" "<<B.Vector3()[2]<<"\n";
1395    //std::cout<<"C "<<C.Vector3()[0]<<"  "<<C.Vector3()[1]<<" "<<C.Vector3()[2]<<"\n";
1396    //std::cout<<"D "<<D.Vector3()[0]<<"  "<<D.Vector3()[1]<<" "<<D.Vector3()[2]<<"\n";
1397    
1398                }
1399                else if(V1.Ref().insideDomain())
1400                {
1401    std::cout<<"be careful triangle 1 \n";
1402                    Vertex A=V1;
1403                    Vertex B=(vv[nn][0]);
1404                    Vertex C=(vv[nn][1]);
1405                    Vertex D=(vv[nn][2]);
1406                    vertexplot.push_back(A);
1407                    vertexplot.push_back(C);
1408                    vertexplot.push_back(D);
1409                    vertexplot.push_back(A);
1410                    vertexplot.push_back(C);
1411                    vertexplot.push_back(B);
1412                    
1413                }
1414                else if(V2.Ref().insideDomain())
1415                {
1416    std::cout<<"be careful triangle 2 \n";
1417                    Vertex A=V2;
1418                    Vertex B=(vv[nn][0]);
1419                    Vertex C=(vv[nn][1]);
1420                    Vertex D=(vv[nn][2]);
1421                    vertexplot.push_back(A);
1422                    vertexplot.push_back(C);
1423                    vertexplot.push_back(D);
1424                    vertexplot.push_back(A);
1425                    vertexplot.push_back(B);
1426                    vertexplot.push_back(D);
1427    //std::cout<<"A "<<A.Vector3()[0]<<"  "<<A.Vector3()[1]<<" "<<A.Vector3()[2]<<"\n";
1428    //std::cout<<"B "<<B.Vector3()[0]<<"  "<<B.Vector3()[1]<<" "<<B.Vector3()[2]<<"\n";
1429    //std::cout<<"C "<<C.Vector3()[0]<<"  "<<C.Vector3()[1]<<" "<<C.Vector3()[2]<<"\n";
1430    //std::cout<<"D "<<D.Vector3()[0]<<"  "<<D.Vector3()[1]<<" "<<D.Vector3()[2]<<"\n";
1431                    
1432                    
1433                }
1434            }
1435           nn++;
1436        }
1437        
1438    }
1439    */
1440  size_t ntr=0;  size_t ntr=0;
1441  for(size_t n = 0; n<scene.nbObjects(); ++n)  for(size_t n = 0; n<scene.nbObjects(); ++n)
1442  ntr+=(*(*__internals).verticesListes[n]).size();  ntr+=(*(*__internals).verticesListes[n]).size();
1443    
1444  s_mesh.setNumberOfVertices(ntr);  s_mesh.setNumberOfVertices(ntr);
1445    
1446     //   s_mesh.setNumberOfVertices((*__internals).verticesList.size());  // //  s_mesh.setNumberOfVertices((*__internals).verticesList.size());
1447    
1448  //std::cout<<"nb tot "<<s_mesh.numberOfVertices()<<"\n";  //std::cout<<"nb tot "<<s_mesh.numberOfVertices()<<"\n";
1449    
# Line 983  s_mesh.setNumberOfVertices(ntr); Line 1453  s_mesh.setNumberOfVertices(ntr);
1453  //      =  (*__internals).verticesList.begin();  //      =  (*__internals).verticesList.begin();
1454  //      splitVertex != (*__internals).verticesList.end(); ++splitVertex) {  //      splitVertex != (*__internals).verticesList.end(); ++splitVertex) {
1455  //  //
1456    //  s_mesh.vertex(i) = (*(*splitVertex).second);  //  s_mesh.vertex(i) = (*(*splitVertex).second);
1457      //delete (*splitVertex).second;  //delete (*splitVertex).second;
1458      //(*splitVertex).second = &s_mesh.vertex(i);  // (*splitVertex).second = &s_mesh.vertex(i);
1459      //i++;  // i++;
1460    //}  //}
1461  for(size_t n = 0; n<scene.nbObjects(); ++n)  for(size_t n = 0; n<scene.nbObjects(); ++n)
1462  {  {
1463      for(SurfaceMeshGenerator::Internals::VerticesList::iterator splitVertex      for(SurfaceMeshGenerator::Internals::VerticesList::iterator splitVertex
# Line 1049  for(size_t n = 0; n<scene.nbObjects(); + Line 1519  for(size_t n = 0; n<scene.nbObjects(); +
1519  }  }
1520    
1521    //   plot2(s_mesh, toTreatHexahedra, SMesh);    //   plot2(s_mesh, toTreatHexahedra, SMesh);
1522    //   plotmedit(s_mesh, toTreatHexahedra,nobj);     //  plotmedit(s_mesh, toTreatHexahedra,nobj);
1523    
1524    /*{
1525    std::ofstream sortie("plot.mesh");
1526    
1527        sortie<<"MeshVersionFormatted 1"<<" \n";
1528        sortie<<" \n";
1529        sortie<<"Dimension"<<" \n";
1530        sortie<<"3"<<" \n";
1531        sortie<<" \n";
1532        sortie<<"Vertices"<<" \n";
1533        sortie<<vertexplot.size()<<"\n";
1534        for(size_t i=0;i<vertexplot.size();++i)
1535        {
1536            sortie<<(*(&vertexplot[i])).Vector3()[0]<<" "<<(*(&vertexplot[i])).Vector3()[1]<<" "<<(*(&vertexplot[i])).Vector3()[2]<<" 0"<<" \n";
1537        }
1538    
1539        sortie<<" \n";
1540        sortie<<"Triangles"<<"\n";
1541        sortie<<vertexplot.size()/3<<"\n";
1542    
1543        for (size_t i=0; i<vertexplot.size()/3; ++i)
1544        {
1545            sortie<<3*i+1<<" "<<3*i+2<<" "<<3*i+3<<" "<<1<<"\n";
1546        }
1547        sortie<<" \n";
1548        sortie<<"End"<<"\n";
1549    }*/
1550    
1551    
1552    
1553  }  }
1554    
# Line 1067  SurfaceMeshGenerator::~SurfaceMeshGenera Line 1566  SurfaceMeshGenerator::~SurfaceMeshGenera
1566  }  }
1567    
1568    
1569    void inout(size_t iter,const size_t nn,TinyVector<3,int> ttn,size_t& re0,size_t& re1,size_t& re2,size_t& ii0,size_t& ii1,size_t& ii2)
1570    {
1571    
1572           if(re0==1 and re1==1)
1573           {
1574               if(ttn[0]!=0 and iter>1)
1575    std::cout<<"-- cas indetermine\n";
1576           }
1577           else if((re0+re1==1))
1578           {
1579               if(ttn[0]==0)
1580               {
1581    //std::cout<<"-- 1 0\n";
1582                   re0=0;
1583                   re1=0;
1584               }
1585               else
1586               {
1587    //std::cout<<"-- 1 in 0 out on sait traite\n";
1588                   ii0=re0;
1589                   ii1=re1;
1590               }
1591           }
1592           else if((re0==2 and re1==0) or
1593                   (re0==0 and re1==2))
1594           {
1595               if(ttn[0]==0)
1596               {
1597                   if(ttn[0]+ttn[1]+ttn[2]==0)
1598                   {
1599                       re0=2;
1600                       re1=2;
1601                   }
1602                   else
1603                   {
1604    std::cout<<"-- 0 2 pas coupe\n";
1605    std::cout<<ttn[0]<<" "<<ttn[1]<<" "<<ttn[2]<<"\n";
1606    std::cout<<re0<<" "<<re1<<"\n";
1607                   }
1608             /*      if(re0==2 and (ttn[1]==1 or ttn[2]==1))
1609                   {
1610                       re0=-nn;
1611                       ii0=-nn;
1612                   }
1613                   else
1614                   {
1615                       re1=-nn;
1616                       ii1=-nn;
1617                   }*/
1618                 //  re0=0;
1619                   //re1=0;
1620               }
1621               else
1622               {
1623    //std::cout<<"-- 0 2 on sait traite\n";
1624                   ii0=re0;
1625                   ii1=re1;
1626               }
1627           }
1628           else if((re0+re1==3))
1629           {
1630               if(ttn[0]==0)
1631               {
1632    //std::cout<<"-- 1 2 pas coupe\n";
1633                   re0=2;
1634                   re1=2;
1635               }
1636               else
1637               {
1638    //std::cout<<"-- 1 2 on sait traite\n";
1639                   if(re0==1)
1640                   {
1641                       ii0=0;
1642                       ii1=1;
1643                       re0=0;
1644                   }
1645                   else
1646                   {
1647                       ii0=1;
1648                       ii1=0;
1649                       re1=0;
1650                   }
1651               }
1652           }
1653      //     else
1654    //std::cout<<"--0 out 0 out rien a faire\n";
1655    
1656           if(re0==1 and re2==1)
1657           {
1658               if(ttn[1]!=0 and iter>1)
1659    std::cout<<"--* cas indetermine\n";
1660           }
1661           else if((re0==1 and re2==0) or
1662                   (re0==0 and re2==1))
1663           {
1664               if(ttn[1]==0)
1665               {
1666    //std::cout<<"--* 1 O non inter\n";
1667                   re0=0;
1668                   re2=0;
1669               }
1670               else
1671               {
1672    //std::cout<<"--* 1 in 0 out on sait traite\n";
1673                   ii0=re0;
1674                   ii2=re2;
1675               }
1676           }
1677           else if((re0==2 and re2==0) or
1678                   (re0==0 and re2==2))
1679           {
1680               if(ttn[1]==0)
1681               {
1682                   if(ttn[0]+ttn[1]+ttn[2]==0)
1683                   {
1684                       re0=2;
1685                       re2=2;
1686                   }
1687                   else
1688                   {
1689                      
1690    std::cout<<"--* 0 2 pas coupe\n";
1691    std::cout<<ttn[0]<<" "<<ttn[1]<<" "<<ttn[2]<<"\n";
1692    std::cout<<re0<<" "<<re2<<"\n";
1693                   }
1694                /*   if(re0==2 and (ttn[0]==1 or ttn[2]==1))
1695                   {
1696                       re0=-nn;
1697                       ii0=-nn;
1698                   }
1699                   else
1700                   {
1701                       re2=-nn;
1702                       ii2=-nn;
1703                   }*/
1704    
1705             //      re0=0;
1706               //    re2=0;
1707               }
1708               else
1709               {
1710    //std::cout<<"--* 0 2 on sait traite\n";
1711                   ii0=re0;
1712                   ii2=re2;
1713               }
1714           }
1715           else if((re0+re2==3))
1716           {
1717               if(ttn[1]==0)
1718               {
1719    //std::cout<<"--* 1 2 pas coupe\n";
1720                   re0=2;
1721                   re2=2;
1722               }
1723               else
1724               {
1725    //std::cout<<"--* 1 2 on sait traite\n";
1726                   if(re0==1)
1727                   {
1728                       ii0=0;
1729                       ii2=1;
1730                       re0=0;
1731                   }
1732                   else
1733                   {
1734                       ii0=1;
1735                       ii2=0;
1736                       re2=0;
1737                   }
1738               }
1739           }      
1740    //       else
1741    //std::cout<<"--* 0 out 0 out rien a faire\n";
1742    
1743           if(re2==1 and re1==1)
1744           {
1745               if(ttn[2]!=0 and iter>1)
1746    std::cout<<"--+ cas indetermine\n";
1747           }
1748           else if((re2+re1==1))
1749           {
1750               if(ttn[2]==0)
1751               {
1752    //std::cout<<"--+ 1 0 non inter\n";
1753                   re2=0;
1754                   re1=0;
1755               }
1756               else
1757               {
1758    //std::cout<<"--+ 1 in 0 out on sait traite\n";
1759                   ii1=re1;
1760                   ii2=re2;
1761               }
1762           }
1763           else if((re2==2 and re1==0) or
1764                   (re2==0 and re1==2))
1765           {
1766               if(ttn[2]==0)
1767               {
1768                   if(ttn[0]+ttn[1]+ttn[2]==0)
1769                   {
1770                       re1=2;
1771                       re2=2;
1772                   }
1773                   else
1774                   {
1775                      
1776    std::cout<<"--+ 0 2 pas coupe\n";
1777    std::cout<<ttn[0]<<" "<<ttn[1]<<" "<<ttn[2]<<"\n";
1778    std::cout<<re1<<" "<<re2<<"\n";
1779                   }
1780                /*   if(re1==2 and (ttn[1]==1 or ttn[0]==1))
1781                   {
1782                       re1=-nn;
1783                       ii1=-nn;
1784                   }
1785                   else
1786                   {
1787                       re2=-nn;
1788                       ii2=-nn;
1789                   }*/
1790                   //re1=0;
1791                   //re2=0;
1792               }
1793               else
1794               {
1795    //std::cout<<"--+ 0 2 on sait traite\n";
1796                   ii1=re1;
1797                   ii2=re2;
1798               }
1799           }
1800           else if((re2+re1==3))
1801           {
1802               if(ttn[2]==0)
1803               {
1804    //std::cout<<"--+ 1 2 pas coupe\n";
1805                   re2=2;
1806                   re1=2;
1807               }
1808               else
1809               {
1810    //std::cout<<"--+ 1 2 on sait traite\n";
1811                   if(re1==1)
1812                   {
1813                       ii1=0;
1814                       ii2=1;
1815                       re1=0;
1816                   }
1817                   else
1818                   {
1819                       ii1=1;
1820                       ii2=0;
1821                       re2=0;
1822                   }
1823               }
1824           }
1825    //       else
1826    //std::cout<<"--+ 0 out 0 out rien a faire\n";
1827    
1828        if(iter>1 and (((re2==1 and re1==1 and re0==1) or re0+re1+re2==6)
1829                       and ttn[0]+ttn[1]+ttn[2]==0))
1830        {
1831            re0=2;
1832            re1=2;
1833            re2=2;
1834            ii0=1;
1835            ii1=1;
1836            ii2=1;
1837        }
1838        if(iter>2 and re0+re1+re2==0 and ttn[0]+ttn[1]+ttn[2]!=0)
1839        {
1840    std::cout<<"pbs de precision?? tr out cut\n";
1841            ii0=0;
1842            ii1=0;
1843            ii2=0;
1844        }
1845        if(iter==1)
1846        {
1847            if(ttn[0]+ttn[1]+ttn[2]==0 and ((re0==0 and re1==1 and re2==1) or (re0==1 and re1==0 and re2==1)
1848                                     or (re0==1 and re1==1 and re2==0) or (re0==0 and re1==0 and re2==1)
1849                                     or (re0==0 and re1==1 and re2==0) or (re0==1 and re1==0 and re2==0)))
1850            {
1851                re0=0;
1852                re1=0;
1853                re2=0;
1854            }
1855        }
1856    
1857        if(ii0==2)
1858            ii0=1;
1859        if(ii1==2)
1860            ii1=1;
1861        if(ii2==2)
1862            ii2=1;
1863    
1864    }
1865    
1866    void Bary(const Vertex& P,TinyVector<3>& P0,TinyVector<3>& V12,TinyVector<3>& V02,size_t& compt)
1867    {
1868        //std::cout<<" tr vert "<<V02[0]<<" "<<V02[1]<<" "<<V02[2]<<"\n";
1869        //std::cout<<" tr vert "<<V12[0]<<" "<<V12[1]<<" "<<V12[2]<<"\n";
1870        //std::cout<<" P0 "<<P0[0]<<" "<<P0[1]<<" "<<P0[2]<<"\n";
1871        //std::cout<<" P "<<P[0]<<" "<<P[1]<<" "<<P[2]<<"\n";
1872    
1873    
1874    
1875        compt=0;
1876        real_t det,eps=-1e-10;
1877        det=V02[0]*V12[1]-V02[1]*V12[0];
1878        if(det<eps or det>-eps)
1879        {
1880            //std::cout<<"det "<<det<<"\n";
1881            TinyMatrix<2,2> mat;
1882            mat(0,0)=V02[0];
1883            mat(1,0)=V02[1];
1884            mat(0,1)=V12[0];
1885            mat(1,1)=V12[1];
1886            TinyVector<2> bary,secondmb;
1887            secondmb[0]=P0[0];
1888            secondmb[1]=P0[1];
1889            bary=mat.invert()*secondmb;
1890            //std::cout<<mat.invert()(0,0)<<" "<<mat.invert()(1,0)<<" "<<mat.invert()(0,1)<<" "<<mat.invert()(1,1)<<"\n";
1891     //std::cout<<secondmb[0]<<" "<<secondmb[1]<<"\n";
1892     //std::cout<<bary[0]<<" "<<bary[1]<<" "<<1.-bary[0]-bary[1]<<"\n";
1893     //std::cout<<bary[0]>=eps<<" "<<bary[1]>=eps<<" "<<(1.-bary[0]-bary[1])>=eps<<"\n";
1894            real_t temp=(1.-bary[0]-bary[1]);
1895            if(bary[0]>=eps && bary[1]>=eps && temp>=eps)
1896            {
1897                compt+=1;
1898              //  std::cout<<bary[0]<<" "<<bary[1]<<" "<<1.-bary[0]-bary[1]<<" "<<compt<<"\n";
1899         //std::cout<<verif[0]<<" "<<verif[1]<<" \n";
1900         //std::cout<<P0[0]<<" "<<P0[2]<<" \n";
1901            }
1902        }
1903        else
1904            compt+=1;
1905        {
1906            det=V02[0]*V12[2]-V02[2]*V12[0];
1907            if(det<eps or det>-eps)
1908            {
1909                TinyMatrix<2,2> mat;
1910                mat(0,0)=V02[0];
1911                mat(1,0)=V02[2];
1912                mat(1,1)=V12[2];
1913                mat(0,1)=V12[0];
1914                TinyVector<2> bary,secondmb;
1915                secondmb[0]=P0[0];
1916                secondmb[1]=P0[2];
1917                bary=mat.invert()*secondmb;
1918                //std::cout<<"------------------\n";
1919         //std::cout<<mat.invert()(0,0)<<" "<<mat.invert()(1,0)<<" "<<mat.invert()(0,1)<<" "<<mat.invert()(1,1)<<"\n";
1920         //std::cout<<"2 "<<bary[0]<<" "<<bary[1]<<" "<<1.-bary[0]-bary[1]<<"\n";
1921    
1922                TinyVector<2> verif=mat*bary;
1923                if(bary[0]>=eps && bary[1]>=eps && (1-bary[0]-bary[1])>=eps)
1924                {
1925                    compt+=1;
1926                    //std::cout<<bary[0]<<" "<<bary[1]<<" "<<1.-bary[0]-bary[1]<<"\n";
1927      //std::cout<<verif[0]<<" "<<verif[1]<<" \n";
1928      //std::cout<<P0[0]<<" "<<P0[2]<<" \n";
1929    
1930    
1931                }
1932            }
1933            else
1934                compt+=1;
1935        }
1936        {
1937            det=V02[1]*V12[2]-V02[2]*V12[1];
1938            if(det<eps or det>-eps)
1939            {
1940                TinyMatrix<2,2> mat;
1941                mat(0,0)=V02[1];
1942                mat(1,0)=V02[2];
1943                mat(1,1)=V12[2];
1944                mat(0,1)=V12[1];
1945                TinyVector<2> bary,secondmb;
1946                secondmb[0]=P0[1];
1947                secondmb[1]=P0[2];
1948                bary=mat.invert()*secondmb;
1949                //std::cout<<"------------------\n";
1950      //std::cou    t<<mat.invert()(0,0)<<" "<<mat.invert()(1,0)<<" "<<mat.invert()(0,1)<<" "<<mat.invert()(1,1)<<"\n";
1951      //std::cout<<"3 "<<bary[0]<<" "<<bary[1]<<" "<<1.-bary[0]-bary[1]<<"\n";
1952    
1953                TinyVector<2> verif=mat*bary;
1954    
1955                if(bary[0]>=eps && bary[1]>=eps && (1-bary[0]-bary[1])>=eps)
1956                {
1957                    compt+=1;
1958                    //std::cout<<bary[0]<<" "<<bary[1]<<" "<<1.-bary[0]-bary[1]<<" "<<compt<<"\n";
1959          //std::cout<<verif[0]<<" "<<verif[1]<<" \n";
1960          //std::cout<<P0[0]<<" "<<P0[1]<<" \n";
1961    
1962                }
1963            }
1964            else
1965                compt+=1;
1966        }
1967        
1968        //    if(compt==3)
1969        //  {
1970        //std::cout<<P.Vector3()[0]<<" "<<P.Vector3()[1]<<" "<<P.Vector3()[2]<<" 3 point intersection \n";
1971        //}
1972    }
1973    
1974    
1975    
1976  void plotmedit(const SurfaceMeshOfTriangles& s_mesh,std::set<const Cell*>& toTreatHexahedra,size_t nobj)  void plotmedit(const SurfaceMeshOfTriangles& s_mesh,std::set<const Cell*>& toTreatHexahedra,size_t nobj)
1977  {  {
1978    
1979      //**********************************************************************      // **********************************************************************
1980            
1981    
1982        std::ofstream sortie("intersection.mesh");        std::ofstream sortie("intersection.mesh");
# Line 1091  void plotmedit(const SurfaceMeshOfTriang Line 1997  void plotmedit(const SurfaceMeshOfTriang
1997              sortie<<H(j)[0]<<" "<<H(j)[1]<<" "<<H(j)[2]<<" 0"<<" \n";              sortie<<H(j)[0]<<" "<<H(j)[1]<<" "<<H(j)[2]<<" 0"<<" \n";
1998          }          }
1999      }      }
2000      for(int i=0;i<s_mesh.numberOfVertices();++i)      for(size_t i=0;i<s_mesh.numberOfVertices();++i)
2001      {      {
2002          sortie<<s_mesh.vertex(i)[0]<<" "<<s_mesh.vertex(i)[1]<<" "<<s_mesh.vertex(i)[2]<<" 0"<<" \n";          sortie<<s_mesh.vertex(i)[0]<<" "<<s_mesh.vertex(i)[1]<<" "<<s_mesh.vertex(i)[2]<<" 0"<<" \n";
2003      }      }
# Line 1128  void plotmedit(const SurfaceMeshOfTriang Line 2034  void plotmedit(const SurfaceMeshOfTriang
2034          for(size_t j=0; j<8; ++j)          for(size_t j=0; j<8; ++j)
2035              sortie<<SMesh.vertexNumber(H(j))<<" \n";              sortie<<SMesh.vertexNumber(H(j))<<" \n";
2036  }*/  }*/
2037  for(int i=0;i<toTreatHexahedra.size();i++)  for(size_t i=0;i<toTreatHexahedra.size();i++)
2038  sortie<<8*i+1<<" "<<8*i+2<<" "<<8*i+3<<" "<<8*i+4<<" "<<8*i+5<<" "<<8*i+6<<" "<<8*i+7<<" "<<8*i+8<<" 3"<<"\n";  sortie<<8*i+1<<" "<<8*i+2<<" "<<8*i+3<<" "<<8*i+4<<" "<<8*i+5<<" "<<8*i+6<<" "<<8*i+7<<" "<<8*i+8<<" 3"<<"\n";
2039      sortie<<"End"<<" \n";      sortie<<"End"<<" \n";
2040    

Legend:
Removed from v.1.7  
changed lines
  Added in v.1.8

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