ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/CAD4FE/src/CAD4FE_OptimizeEdgeSwap.cpp
Revision: 569
Committed: Thu Oct 16 14:36:31 2014 UTC (10 years, 7 months ago) by foucault
File size: 17158 byte(s)
Log Message:
Mise à jour pour CAD4FE (Gilles) : operation 1 (tentative)

File Contents

# User Rev Content
1 foucault 27 //---------------------------------------------------------------------------
2     #include <set>
3    
4     #pragma hdrstop
5    
6     #include "gestionversion.h"
7     #include <math.h>
8     #include "CAD4FE_OptimizeEdgeSwap.h"
9     #include "ot_mathematique.h"
10 foucault 569 #include "ot_boite_2d.h"
11 foucault 27 #include "tpl_map_entite.h"
12     #include "CAD4FE_m3d_MCTriangle.h"
13     #include "mg_maillage.h"
14    
15     #include "CAD4FE_Intersection_MCSegment_MCSegment.h"
16     #include "CAD4FE_MCSegment_Middle.h"
17     #include "CAD4FE_MCSegment.h"
18     #include "CAD4FE_MCNode.h"
19 foucault 569 #include "ot_algorithme_geometrique.h"
20     #include "CAD4FE_Geometric_Tools.h"
21 foucault 27
22     //---------------------------------------------------------------------------
23    
24     #pragma package(smart_init)
25    
26     using namespace CAD4FE;
27    
28     OptimizeEdgeSwap::OptimizeEdgeSwap(MG_MAILLAGE* __mesh)
29     {
30     _mesh=__mesh;
31     maxDiedralAngle = 45*M_PI/180;
32     }
33    
34     MG_SEGMENT* OptimizeEdgeSwap::GetSegmentAfter(MG_TRIANGLE *__triangle, MG_SEGMENT *__segment)
35     {
36     if (__segment == __triangle->get_segment1())
37     return __triangle->get_segment2();
38     if (__segment == __triangle->get_segment2())
39     return __triangle->get_segment3();
40     if (__segment == __triangle->get_segment3())
41     return __triangle->get_segment1();
42     return NULL;
43     }
44    
45     MG_SEGMENT* OptimizeEdgeSwap::GetSegmentBefore(MG_TRIANGLE *__triangle, MG_SEGMENT *__segment)
46     {
47     if (__segment == __triangle->get_segment1())
48     return __triangle->get_segment3();
49     if (__segment == __triangle->get_segment2())
50     return __triangle->get_segment1();
51     if (__segment == __triangle->get_segment3())
52     return __triangle->get_segment2();
53     return NULL;
54     }
55    
56     // Return: -1 = triangle is degenerate (a segment or point)
57     // 0 = triangles do not overlap (no projection)
58     // 1 = triangles do overlap but are coplanar
59     // 2 = triangles do overlap AND are coplanar
60     int OptimizeEdgeSwap::CheckOverlapTriangle(MG_NOEUD * __tri[2][3])
61     {
62     for (int i=0;i<2;i++)
63     for (int j=0;j<3;j++)
64     {
65     MG_NOEUD ** t = __tri[(i+1)%2];
66     MG_NOEUD * v = __tri[i][j];
67    
68     if (t[0]==v||t[1]==v||t[2]==v) continue;
69    
70     double I[3],IT;
71     int res = OT_ALGORITHME_GEOMETRIQUE::project_PointTriangle(v->get_coord(),t[0]->get_coord(),t[1]->get_coord(),t[2]->get_coord(),I,&IT);
72     if (res != 0)
73     return res;
74     }
75     return 0;
76     }
77    
78     double OptimizeEdgeSwap::SwapScore(MG_SEGMENT *__segment)
79     {
80     return Swap(__segment,false);
81     }
82    
83     double OptimizeEdgeSwap::SwapSegment(MG_SEGMENT *__segment)
84     {
85     MG_TRIANGLE * triangle[2]={__segment->get_lien_triangle()->get(0),__segment->get_lien_triangle()->get(1)};
86     MG_NOEUD * no[4]; // no[]={A,C,B,D} where AB is the segment shared by both triangles AB, C and D nodes opposed to segment AB
87     no[2] = GetStartNode(triangle[0],__segment);
88     no[1] = GetNodeBefore(triangle[0],no[2]);
89     no[0] = GetNodeBefore(triangle[0],no[1]);
90     no[3] = GetNodeBefore(triangle[1],no[0]);
91     MG_SEGMENT * seg[4];
92     seg[0] = _mesh->get_mg_segment(no[0]->get_id(),no[1]->get_id()); // AC
93     seg[1] = _mesh->get_mg_segment(no[1]->get_id(),no[2]->get_id()); // CB
94     seg[2] = _mesh->get_mg_segment(no[2]->get_id(),no[3]->get_id()); // BD
95     seg[3] = _mesh->get_mg_segment(no[3]->get_id(),no[0]->get_id()); // DA
96    
97     double crit21 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[3]->get_coord(),no[0]->get_coord());
98     double crit22 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[2]->get_coord(),no[3]->get_coord());
99     MCSegment * swapSegment = new MCSegment(triangle[0]->get_lien_topologie(), (MCNode*)no[1], (MCNode*)no[3] );
100     if (swapSegment->GetRefFaceMapping().size() == 0)
101     {
102     delete swapSegment;
103     return 0;
104     }
105     _mesh->ajouter_mg_segment(swapSegment);
106     M3D_MCTriangle * newTriangle[2];
107     newTriangle[0] = new M3D_MCTriangle(triangle[0]->get_lien_topologie(), (MCNode*)no[1], (MCNode*)no[2], (MCNode*)no[3], (MCSegment*)seg[1], (MCSegment*)seg[2], swapSegment);
108     newTriangle[0]->change_qualite(crit21);
109     newTriangle[1] = new M3D_MCTriangle(triangle[0]->get_lien_topologie(), (MCNode*)no[3], (MCNode*)no[0], (MCNode*)no[1], (MCSegment*)seg[3], (MCSegment*)seg[0], swapSegment);
110     newTriangle[1]->change_qualite(crit22);
111     for (int i=0; i<2; i++)_mesh->ajouter_mg_triangle(newTriangle[i]);
112     for (int i=0; i<2; i++)_mesh->supprimer_mg_triangleid(triangle[i]->get_id());
113     _mesh->supprimer_mg_segmentid(__segment->get_id());
114    
115     return 1;
116     }
117    
118     /*double OptimizeEdgeSwap::Swap(MG_SEGMENT *__segment, bool __execute)
119     {
120     double score1;
121     if (__segment->get_lien_triangle()->get_nb() != 2)
122     {
123     return -1;
124     }
125     if (__segment->get_lien_triangle()->get(0)->get_lien_topologie() != __segment->get_lien_triangle()->get(1)->get_lien_topologie())
126     {
127     return -2;
128     }
129     if (__segment->get_lien_topologie()&&__segment->get_lien_topologie()->get_dimension()!=2)
130     {
131     return -3;
132     }
133    
134     MG_TRIANGLE * triangle[2]={__segment->get_lien_triangle()->get(0),__segment->get_lien_triangle()->get(1)};
135     MG_NOEUD * no[4]; // no[]={A,C,B,D} where AB is the segment shared by both triangles AB, C and D nodes opposed to segment AB
136     no[2] = GetStartNode(triangle[0],__segment);
137     no[1] = GetNodeBefore(triangle[0],no[2]);
138     no[0] = GetNodeBefore(triangle[0],no[1]);
139     no[3] = GetNodeBefore(triangle[1],no[0]);
140     MG_SEGMENT * seg[4];
141     seg[0] = _mesh->get_mg_segment(no[0]->get_id(),no[1]->get_id()); // AC
142     seg[1] = _mesh->get_mg_segment(no[1]->get_id(),no[2]->get_id()); // CB
143     seg[2] = _mesh->get_mg_segment(no[2]->get_id(),no[3]->get_id()); // BD
144     seg[3] = _mesh->get_mg_segment(no[3]->get_id(),no[0]->get_id()); // DA
145    
146     double crit11,crit12,crit21,crit22;
147     crit11 = OPERATEUR::qualite_triangle(no[0]->get_coord(),no[1]->get_coord(),no[2]->get_coord());
148     crit12 = OPERATEUR::qualite_triangle(no[2]->get_coord(),no[3]->get_coord(),no[0]->get_coord());
149     crit21 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[3]->get_coord(),no[0]->get_coord());
150     crit22 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[2]->get_coord(),no[3]->get_coord());
151    
152     double crit1=std::min(crit11,crit12);
153     double crit2=std::min(crit21,crit22);
154    
155     score1 = (crit1 > .5) ? -1 : (crit2-crit1);
156     if (score1 < 0) return score1;
157    
158     MG_NOEUD* tris2[2][3]={{no[3],no[1],no[2]},{no[1],no[3],no[0]}};
159     int testOverlapTriangles = CheckOverlapTriangle(tris2);
160     score1 = ( testOverlapTriangles != 0 ) ? -1 : score1;
161     if (score1 < 0) return score1;
162    
163     double diedralAngle2 = GeometricTools::ComputeDiedralAngle(tris2);
164     score1 = ( diedralAngle2 > maxDiedralAngle ) ? -1 : score1;
165     if (score1 < 0) return score1;
166    
167     if (__execute && score1 > 0)
168     {
169     double crit21 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[3]->get_coord(),no[0]->get_coord());
170     double crit22 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[2]->get_coord(),no[3]->get_coord());
171     MG_ELEMENT_TOPOLOGIQUE * topo = __segment->get_lien_triangle()->get(0)->get_lien_topologie() ;
172     MCSegment * swapSegment = new MCSegment(topo, (MCNode*)no[1], (MCNode*)no[3] );
173     _mesh->ajouter_mg_segment(swapSegment);
174     M3D_MCTriangle * newTriangle[2];
175     newTriangle[0] = new M3D_MCTriangle(topo, (MCNode*)no[1], (MCNode*)no[2], (MCNode*)no[3], (MCSegment*)seg[1], (MCSegment*)seg[2], swapSegment);
176     newTriangle[0]->change_qualite(crit21);
177     newTriangle[1] = new M3D_MCTriangle(topo, (MCNode*)no[3], (MCNode*)no[0], (MCNode*)no[1], (MCSegment*)seg[3], (MCSegment*)seg[0], swapSegment);
178     newTriangle[1]->change_qualite(crit22);
179     for (int i=0; i<2; i++)_mesh->ajouter_mg_triangle(newTriangle[i]);
180     for (int i=0; i<2; i++)_mesh->supprimer_mg_triangleid(triangle[i]->get_id());
181     _mesh->supprimer_mg_segmentid(__segment->get_id());
182     }
183    
184     return score1;
185     } */
186    
187    
188     double OptimizeEdgeSwap::Swap(MG_SEGMENT *__segment, bool __execute)
189     {
190     double score1;
191     if (__segment->get_lien_triangle()->get_nb() != 2)
192     {
193     return -1;
194     }
195     if (__segment->get_lien_triangle()->get(0)->get_lien_topologie() != __segment->get_lien_triangle()->get(1)->get_lien_topologie())
196     {
197     return -2;
198     }
199     if (__segment->get_lien_topologie()&&__segment->get_lien_topologie()->get_dimension()!=2)
200     {
201     return -3;
202     }
203    
204     MG_TRIANGLE * triangle[2]={__segment->get_lien_triangle()->get(0),__segment->get_lien_triangle()->get(1)};
205     MG_NOEUD * no[4]; // no[]={A,C,B,D} where AB is the segment shared by both triangles AB, C and D nodes opposed to segment AB
206     no[2] = GetStartNode(triangle[0],__segment);
207     no[1] = GetNodeBefore(triangle[0],no[2]);
208     no[0] = GetNodeBefore(triangle[0],no[1]);
209     no[3] = GetNodeBefore(triangle[1],no[0]);
210     MG_SEGMENT * seg[4];
211     seg[0] = _mesh->get_mg_segment(no[0]->get_id(),no[1]->get_id()); // AC
212     seg[1] = _mesh->get_mg_segment(no[1]->get_id(),no[2]->get_id()); // CB
213     seg[2] = _mesh->get_mg_segment(no[2]->get_id(),no[3]->get_id()); // BD
214     seg[3] = _mesh->get_mg_segment(no[3]->get_id(),no[0]->get_id()); // DA
215    
216     MG_NOEUD* tris[2][2][3]={{{no[2],no[3],no[0]},{no[0],no[1],no[2]}}, // triangles of the current configuration
217     {{no[3],no[1],no[2]},{no[1],no[3],no[0]}}}; // triangles of the candidate configuration (after segment swapping)
218    
219    
220     int testOverlapTriangles = CheckOverlapTriangle(tris[1]);
221     if (testOverlapTriangles) return -1;
222    
223     // projection des triangles sur le plan tangent
224     // Critère de qualité dans le plan tangent
225     double crit[2][2];
226     double area[2][2];
227     for (int k=0;k<2;k++)
228     for (int i=0;i<2;i++)
229     {
230     OT_VECTEUR_3D normal[4];
231     MCFace * face = (MCFace *)triangle[i]->get_lien_topologie();
232     for (int j=0;j<3;j++)
233     ((MCNode*)tris[k][i][j])->NormalMCFace(face,normal[j]);
234     normal[3] = .3333*(normal[0]+normal[1]+normal[2]);
235     // frame (x,y,n)
236     OT_MATRICE_3D tangentPlane = GeometricTools::GetPlaneFrame(normal[3]);
237     OT_MATRICE_3D trf; tangentPlane.transpose(trf);
238     OT_VECTEUR_3D tpProjNodes[3];
239     for (int j=0;j<3;j++)
240     {
241     OT_VECTEUR_3D p(tris[k][i][j]->get_coord());
242     tpProjNodes[j] = trf*p;
243     // projection --> n component = 0
244     tpProjNodes[j][2] = 0;
245     }
246     crit[k][i]=OPERATEUR::qualite_triangle(tpProjNodes[0],tpProjNodes[1],tpProjNodes[2]);
247     area[k][i]=OT_ALGORITHME_GEOMETRIQUE::Area_Triangle(tpProjNodes[0],tpProjNodes[1],tpProjNodes[2]);
248     }
249     double crit1=std::min(crit[0][0],crit[0][1]);
250     double crit2=std::min(crit[1][0],crit[1][1]);
251    
252     // critere de précision géométrique
253     // dist_P[i] = la distance entre la diagonale de la ieme configuration et la MCFace
254     // voir page 110 de these_jfl.pdf
255     double critGeomAccuracy[2];
256     {
257     MG_ELEMENT_TOPOLOGIQUE * topo = __segment->get_lien_triangle()->get(0)->get_lien_topologie() ;
258     MCSegment * AB = (MCSegment*)__segment;
259     MCNode Pgeo;
260     // MCSegment * CD = new MCSegment(topo, (MCNode*)no[1], (MCNode*)no[3] );
261     // CD->evaluer_geo(((area[0][0]*0+area[0][1]*1)/(area[0][0]+area[0][1])),P[0]);
262     // P1 = pt d'intersection géodésique entre les diagonales AB et CD
263     {
264     double wA, wB;
265     wA=(AB->get_noeud1() == no[0])?0:1;
266     if(AB->get_orientation() == -1) wA=(wA==0)?1:0;
267     wB=(wA==0)?1:0;
268     AB->evaluer_geo(((area[1][0]*wA+area[1][1]*wB)/(area[1][0]+area[1][1])),&Pgeo);
269     }
270    
271     double P_Euclid[2][3];
272     for (int i=0;i<3;i++) {
273     // P' = ( |T1|*C + |T2|*D ) / ( |T1| + |T2| ) où |T| = aire triangle
274     P_Euclid[1][i]=((area[0][0]*no[1]->get_coord()[i]+area[0][1]*no[3]->get_coord()[i])/(area[0][0]+area[0][1]));
275     // P = ( |T1'|*A + |T2'|*B ) / ( |T1'| + |T2'| ) où |T| = aire triangle
276     P_Euclid[0][i]=((area[1][0]*no[0]->get_coord()[i]+area[1][1]*no[2]->get_coord()[i])/(area[1][0]+area[1][1]));
277     }
278     double distP[2],segLength[2];
279     for (int i=0;i<2;i++)
280     distP[i]=OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(P_Euclid[i],Pgeo.get_coord());
281     for (int i=0;i<2;i++)
282     segLength[i]=OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(no[i]->get_coord(),no[i+2]->get_coord());
283     double maxDistP=.5*std::max(segLength[0],segLength[1]);
284     for (int i=0;i<2;i++)
285     if ( distP[i] < maxDistP )
286     critGeomAccuracy[i] = 1-distP[i]/maxDistP;
287     else
288     critGeomAccuracy[i] = 0;
289     }
290    
291     crit1 += critGeomAccuracy[0];
292     crit2 += critGeomAccuracy[1];
293    
294     if (__execute && crit2 > crit1 )
295     {
296     SwapSegment(__segment);
297     }
298    
299     return (crit2>crit1)?crit2:-1;
300     }
301    
302    
303     MG_NOEUD * OptimizeEdgeSwap::GetStartNode(MG_TRIANGLE *__t, MG_SEGMENT *__s)
304     {
305     MG_NOEUD *TN[3];
306     TN[0] = __t->get_noeud1();
307     TN[1] = __t->get_noeud2();
308     TN[2] = __t->get_noeud3();
309     for (int i=0; i<3; i++)
310     {
311     if (__s->get_noeud1() == TN[i])
312     {
313     if (__s->get_noeud2() == TN[(i+1)%3])
314     return __s->get_noeud1();
315     else
316     return __s->get_noeud2();
317     }
318     }
319     return NULL;
320     }
321    
322     MG_NOEUD * OptimizeEdgeSwap::GetNodeAfter(MG_TRIANGLE *__t, MG_NOEUD *__n)
323     {
324     MG_NOEUD *TN[3];
325     TN[0] = __t->get_noeud1();
326     TN[1] = __t->get_noeud2();
327     TN[2] = __t->get_noeud3();
328     for (int i=0; i<3; i++)
329     {
330     if (TN[i]==__n)
331     return TN[(i+1)%3];
332     }
333     return NULL;
334     }
335    
336     MG_NOEUD * OptimizeEdgeSwap::GetNodeBefore(MG_TRIANGLE *__t, MG_NOEUD *__n)
337     {
338     MG_NOEUD *TN[3];
339     TN[0] = __t->get_noeud1();
340     TN[1] = __t->get_noeud2();
341     TN[2] = __t->get_noeud3();
342     for (int i=0; i<3; i++)
343     {
344     if (TN[i]==__n)
345     return TN[(i+2)%3];
346     }
347     return NULL;
348     }
349    
350     MG_NOEUD * OptimizeEdgeSwap::GetOppositeNode(MG_TRIANGLE *__t, MG_SEGMENT *__s)
351     {
352     MG_NOEUD *TN[3];
353     TN[0] = __t->get_noeud1();
354     TN[1] = __t->get_noeud2();
355     TN[2] = __t->get_noeud3();
356     for (int i=0; i<3; i++)
357     if (TN[i] != __s->get_noeud1() && TN[i] != __s->get_noeud2())
358     return TN[i];
359     return NULL;
360     }
361    
362     void OptimizeEdgeSwap::RemoveVolumeMesh()
363     {
364     LISTE_MG_SEGMENT::iterator itSeg;
365     LISTE_MG_NOEUD::iterator itNo;
366     LISTE_MG_TRIANGLE::iterator itTri;
367     _mesh->supprimer_tout_mg_tetra();
368     for (MG_TRIANGLE * tri = _mesh->get_premier_triangle(itTri);tri; tri = _mesh->get_suivant_triangle(itTri))
369     {
370     if (tri->get_lien_topologie()->get_dimension()==3)
371     _mesh->supprimer_mg_triangleid(tri->get_id());
372     }
373     for (MG_SEGMENT * seg = _mesh->get_premier_segment(itSeg);seg; seg = _mesh->get_suivant_segment(itSeg))
374     {
375     if (seg->get_lien_topologie()->get_dimension()==3)
376     _mesh->supprimer_mg_segmentid(seg->get_id());
377     }
378     for (MG_NOEUD * no = _mesh->get_premier_noeud(itNo);no; no = _mesh->get_suivant_noeud(itNo))
379     {
380     if (no->get_lien_topologie()->get_dimension()==3)
381     _mesh->supprimer_mg_noeudid(no->get_id());
382     }
383     }
384    
385     int OptimizeEdgeSwap::Optimize(MG_FACE * __face)
386     {
387     int nbSwap = 0;
388     LISTE_MG_SEGMENT::iterator itSeg;
389     MG_SEGMENT * seg = _mesh->get_premier_segment(itSeg);
390     std::map<double, MG_SEGMENT*> lst;
391    
392     for (int j=0;j<3;j++)
393     {
394     lst.clear();
395     int tmpnbSwap = nbSwap;
396     for (seg = _mesh->get_premier_segment(itSeg);seg; seg = _mesh->get_suivant_segment(itSeg))
397     {
398     MG_ELEMENT_TOPOLOGIQUE * topo = seg->get_lien_topologie();
399     if (__face != NULL && topo != __face)
400     continue;
401     if (topo->get_dimension()!=2)
402     continue;
403     lst[SwapScore(seg)]=seg;
404     }
405    
406     for (std::map<double, MG_SEGMENT*>::reverse_iterator i=lst.rbegin();i!=lst.rend();i++)
407     {
408     double score1 = i->first;
409     seg=i->second;
410    
411     if (score1<0)
412     break;
413    
414     double score2 = SwapSegment(seg);
415     if (score2 > 0)
416     {
417     nbSwap++;
418     }
419     }
420     if (nbSwap-tmpnbSwap==0)
421     break;
422     }
423    
424     return nbSwap;
425     }
426    
427     int OptimizeEdgeSwap::OptimizeAllFaces()
428     {
429     return Optimize(0);
430     }
431