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: |
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 |
|
|
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 |
|
|
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, |
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, |
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, |
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 |
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; |
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); |
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; |
789 |
|
|
790 |
} |
} |
791 |
// } |
// } |
792 |
//listOfVertexMesh.clear(); |
//listOfVertexMesh.clear(); |
793 |
} |
} |
794 |
|
|
795 |
void SurfaceMeshGenerator::Internals::__dataStructureConvertion(ObjectToTreat& O, |
void SurfaceMeshGenerator::Internals::__dataStructureConvertion(ObjectToTreat& O, |
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 |
|
|
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, |
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]); |
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; |
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) { |
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 |
} |
} |
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 |
|
|
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"; |
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); |
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 |
} |
} |
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 |
|
|
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: { |
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; |
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; |
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) { |
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 |
} |
} |
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); |
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; |
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) { |
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 |
|
|
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); |
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"; |
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, |
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); |
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 |
|
|
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() ; |
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 |
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; |
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:: |
3084 |
currentCell, |
currentCell, |
3085 |
cellList); |
cellList); |
3086 |
if(!(jobDone2)) { |
if(!(jobDone2)) { |
3087 |
jobDone=jobDone2; |
jobDone=jobDone2; |
3088 |
} |
} |
3089 |
} |
} |
3090 |
|
|
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(); |
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 |
} |
} |
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; |
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"; |