ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur3d_optimisation.cpp
Revision: 1189
Committed: Tue Feb 4 17:26:49 2025 UTC (3 months ago) by francois
File size: 25259 byte(s)
Log Message:
Version 5.0 de MAGIC. Integration de ALGLIB pour faire de l'optimisation. ALGLIB se download automatiquement en executant un script dans le repertoire config update_magic.bash


File Contents

# Content
1 //####//------------------------------------------------------------
2 //####//------------------------------------------------------------
3 //####// MAGiC
4 //####// Jean Christophe Cuilliere et Vincent FRANCOIS
5 //####// Departement de Genie Mecanique - UQTR
6 //####//------------------------------------------------------------
7 //####// MAGIC est un projet de recherche de l equipe ERICCA
8 //####// du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
9 //####// http://www.uqtr.ca/ericca
10 //####// http://www.uqtr.ca/
11 //####//------------------------------------------------------------
12 //####//------------------------------------------------------------
13 //####//
14 //####// mailleur3d_optimisation.cpp
15 //####//
16 //####//------------------------------------------------------------
17 //####//------------------------------------------------------------
18 //####// COPYRIGHT 2000-2024
19 //####// jeu 13 jun 2024 11:58:55 EDT
20 //####//------------------------------------------------------------
21 //####//------------------------------------------------------------
22
23
24 #include "gestionversion.h"
25 #include "mailleur3d_optimisation.h"
26 #include "ot_mathematique.h"
27 #include "m3d_tetra.h"
28 #include "m3d_triangle.h"
29 #include "m3d_noeud.h"
30 #include "mg_maillage.h"
31 #include <math.h>
32
33
34
35
36 MAILLEUR3D_OPTIMISATION::MAILLEUR3D_OPTIMISATION(MG_MAILLAGE* mgmai,int niv):MAILLEUR(false),mg_maillage(mgmai),niveau_optimisation(niv)
37 {
38 o3d_data();
39 o3d_data2();
40 LISTE_MG_TETRA::iterator it;
41 TPL_MAP_ENTITE<M3D_TETRA*> listetet3d;
42 for (MG_TETRA* tet=mg_maillage->get_premier_tetra(it);tet!=NULL;tet=mg_maillage->get_suivant_tetra(it))
43 if (tet->get_type_entite()!=MAGIC::TYPE_ENTITE::IDM3D_TETRA)
44 {
45 M3D_TETRA* tet3d=new M3D_TETRA(tet->get_id(),tet->get_lien_topologie(),tet->get_noeud1(),tet->get_noeud2(),tet->get_noeud3(),tet->get_noeud4(),tet->get_triangle1(),tet->get_triangle2(),tet->get_triangle3(),tet->get_triangle4(),tet->get_origine());
46 listetet3d.ajouter(tet3d);
47 }
48
49 TPL_MAP_ENTITE<M3D_TETRA*>::ITERATEUR it2;
50
51 for (M3D_TETRA* tet=listetet3d.get_premier(it2);tet!=NULL;tet=listetet3d.get_suivant(it2))
52 {
53 mg_maillage->supprimer_mg_tetraid(tet->get_id());
54 mg_maillage->ajouter_mg_tetra(tet);
55 }
56 }
57
58 MAILLEUR3D_OPTIMISATION::~MAILLEUR3D_OPTIMISATION()
59 {
60 }
61
62
63
64 double MAILLEUR3D_OPTIMISATION::get_volume(MG_TETRA* tet)
65 {
66 double *xyz1=tet->get_noeud1()->get_coord();
67 double *xyz2=tet->get_noeud2()->get_coord();
68 double *xyz3=tet->get_noeud3()->get_coord();
69 double *xyz4=tet->get_noeud4()->get_coord();
70 OT_VECTEUR_3D ab(xyz1,xyz2);
71 OT_VECTEUR_3D ac(xyz1,xyz3);
72 OT_VECTEUR_3D ad(xyz1,xyz4);
73 double vol=(ab&ac)*ad;
74 vol=vol/6.;
75 return vol;
76 }
77
78 double MAILLEUR3D_OPTIMISATION::get_volume(double *xyz1,double *xyz2,double *xyz3,double *xyz4)
79 {
80 OT_VECTEUR_3D ab(xyz1,xyz2);
81 OT_VECTEUR_3D ac(xyz1,xyz3);
82 OT_VECTEUR_3D ad(xyz1,xyz4);
83 double vol=(ab&ac)*ad;
84 vol=vol/6.;
85 return vol;
86 }
87
88 void MAILLEUR3D_OPTIMISATION::optimise(MG_VOLUME* mgvol)
89 {
90 char mess[255];
91 sprintf(mess," Niveau d'optimisation : %d",niveau_optimisation);
92 if (affichageactif==1) affiche(mess);
93 int nbaoptimiser;
94 int nbaoptimiserapres=mg_maillage->get_nb_mg_tetra();
95 do
96 {
97 nbaoptimiser=nbaoptimiserapres;
98 optimise(mgvol,nbaoptimiserapres);
99 if (nbaoptimiserapres==0)
100 {
101 nbaoptimiser=0;
102 sprintf(mess," Tout est optimise");
103 if (affichageactif==1) affiche(mess);
104 }
105
106 }
107 while (nbaoptimiserapres!=nbaoptimiser);
108 }
109
110
111
112 void MAILLEUR3D_OPTIMISATION::optimise(MG_VOLUME* mgvol,int& nbmauvais)
113 {
114 int passe=1;
115 //passe++;
116 int nbtet=0;
117 if (mgvol!=NULL)
118 {
119 TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
120 for (MG_TETRA* tet=(MG_TETRA*)mgvol->get_lien_maillage()->get_premier(it);tet;tet=(MG_TETRA*)mgvol->get_lien_maillage()->get_suivant(it))
121 {
122 if (mg_maillage->get_mg_tetraid(tet->get_id())==NULL) continue;
123 M3D_TETRA* mtet=(M3D_TETRA*)tet;
124 double qual=OPERATEUR::qualite_tetra(mtet->get_noeud1()->get_coord(),mtet->get_noeud2()->get_coord(),mtet->get_noeud3()->get_coord(),mtet->get_noeud4()->get_coord());
125 mtet->change_qualite(qual);
126 if (qual<0.1*niveau_optimisation) nbtet++;
127 ajouter_ordre_tetra(mtet,0);
128 }
129 }
130 else
131 {
132 LISTE_MG_TETRA::iterator it;
133 for (MG_TETRA* tet=mg_maillage->get_premier_tetra(it);tet;tet=mg_maillage->get_suivant_tetra(it))
134 {
135 if (mg_maillage->get_mg_tetraid(tet->get_id())==NULL) continue;
136 M3D_TETRA* mtet=(M3D_TETRA*)tet;
137 double qual=OPERATEUR::qualite_tetra(mtet->get_noeud1()->get_coord(),mtet->get_noeud2()->get_coord(),mtet->get_noeud3()->get_coord(),mtet->get_noeud4()->get_coord());
138 mtet->change_qualite(qual);
139 if (qual<0.1*niveau_optimisation) nbtet++;
140 ajouter_ordre_tetra(mtet,0);
141 }
142 }
143
144 if (nbtet==0) return;
145 for (int phase=0;phase<2;phase++)
146 {
147 char mess[255];
148 if (phase==0) sprintf(mess," Phase %d : %d a optimiser",2*passe-1+phase,nbtet);
149 if (phase==1) sprintf(mess," Phase %d",2*passe-1+phase);
150 if (affichageactif==1) affiche(mess);
151 while (lst_tetra[phase].size()>0)
152 {
153 ORDRE_TETRA::iterator i=lst_tetra[phase].begin();
154 M3D_TETRA* tet=(*i).second;
155 supprimer_ordre_tetra(tet);
156 double crit_avant=tet->get_qualite();
157 if (crit_avant<0.1*niveau_optimisation)
158 {
159 double crit[4]={-1.,-1.,-1.,-1.};
160 MG_NOEUD* no[4];
161 no[0]=(MG_NOEUD*)tet->get_noeud1();
162 no[1]=(MG_NOEUD*)tet->get_noeud2();
163 no[2]=(MG_NOEUD*)tet->get_noeud3();
164 no[3]=(MG_NOEUD*)tet->get_noeud4();
165 double x[4],y[4],z[4];
166 bouge_point(mgvol,no[0],crit[0],x[0],y[0],z[0]);
167 bouge_point(mgvol,no[1],crit[1],x[1],y[1],z[1]);
168 bouge_point(mgvol,no[2],crit[2],x[2],y[2],z[2]);
169 bouge_point(mgvol,no[3],crit[3],x[3],y[3],z[3]);
170 std::vector<double> vecteur_crit(crit,crit+4);
171 std::vector<double>::iterator it1 = max_element(vecteur_crit.begin(), vecteur_crit.end());
172 double crit_opt=*it1;
173 if (crit_opt>crit_avant)
174 {
175 int num=it1-vecteur_crit.begin();
176 no[num]->change_x(x[num]);
177 no[num]->change_y(y[num]);
178 no[num]->change_z(z[num]);
179 int nb=no[num]->get_lien_tetra()->get_nb();
180 for (int j=0;j<nb;j++)
181 {
182 M3D_TETRA* mtet=(M3D_TETRA*)no[num]->get_lien_tetra()->get(j);
183 double qual=OPERATEUR::qualite_tetra(mtet->get_noeud1()->get_coord(),mtet->get_noeud2()->get_coord(),mtet->get_noeud3()->get_coord(),mtet->get_noeud4()->get_coord());
184 mtet->change_qualite(qual);
185 }
186 crit_avant=crit_opt;
187 }
188 }
189
190 int coquille_non_optimise=1;
191 if (crit_avant<0.1*niveau_optimisation)
192 {
193 COQUILLE coque[6];
194 double crit[6]={-1,-1,-1,-1,-1,-1};
195 remaille_coquille(tet->get_noeud1(),tet->get_noeud2(),crit[0],coque[0]);
196 remaille_coquille(tet->get_noeud1(),tet->get_noeud3(),crit[1],coque[1]);
197 remaille_coquille(tet->get_noeud1(),tet->get_noeud4(),crit[2],coque[2]);
198 remaille_coquille(tet->get_noeud2(),tet->get_noeud3(),crit[3],coque[3]);
199 remaille_coquille(tet->get_noeud2(),tet->get_noeud3(),crit[4],coque[4]);
200 remaille_coquille(tet->get_noeud3(),tet->get_noeud4(),crit[5],coque[5]);
201 std::vector<double> vecteur_crit(crit,crit+6);
202 std::vector<double>::iterator it1 = max_element(vecteur_crit.begin(), vecteur_crit.end());
203 double crit_opt=*it1;
204 if (crit_opt>crit_avant)
205 {
206 coquille_non_optimise=0;
207 int num=it1-vecteur_crit.begin();
208 for (int i=0;i<coque[num].taille;i++)
209 {
210 M3D_TETRA *tettmp=(M3D_TETRA*)coque[num].tet[i];
211 if (tet!=tettmp) supprimer_ordre_tetra(tettmp);
212 mg_maillage->supprimer_mg_tetraid(coque[num].tet[i]->get_id());
213 }
214 for (int i=0;i<2*coque[num].taille-4;i++)
215 {
216 MG_NOEUD* noeud1=coque[num].new_tetra[4*i];
217 MG_NOEUD* noeud2=coque[num].new_tetra[4*i+1];
218 MG_NOEUD* noeud3=coque[num].new_tetra[4*i+2];
219 MG_NOEUD* noeud4=coque[num].new_tetra[4*i+3];
220 MG_TRIANGLE* triangle1=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud3->get_id(),noeud2->get_id());
221 MG_TRIANGLE* triangle2=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud2->get_id(),noeud4->get_id());
222 MG_TRIANGLE* triangle3=mg_maillage->get_mg_triangle(noeud2->get_id(),noeud3->get_id(),noeud4->get_id());
223 MG_TRIANGLE* triangle4=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud4->get_id(),noeud3->get_id());
224 if (triangle1==NULL) triangle1=insere_triangle(mgvol,noeud1,noeud3,noeud2,MAGIC::ORIGINE::MAILLEUR_AUTO);
225 if (triangle2==NULL) triangle2=insere_triangle(mgvol,noeud1,noeud2,noeud4,MAGIC::ORIGINE::MAILLEUR_AUTO);
226 if (triangle3==NULL) triangle3=insere_triangle(mgvol,noeud2,noeud3,noeud4,MAGIC::ORIGINE::MAILLEUR_AUTO);
227 if (triangle4==NULL) triangle4=insere_triangle(mgvol,noeud1,noeud4,noeud3,MAGIC::ORIGINE::MAILLEUR_AUTO);
228 M3D_TETRA* mtet=new M3D_TETRA(mgvol,noeud1,noeud2,noeud3,noeud4,triangle1,triangle2,triangle3,triangle4,MAGIC::ORIGINE::MAILLEUR_AUTO);
229 mg_maillage->ajouter_mg_tetra(mtet);
230 double qual=OPERATEUR::qualite_tetra(mtet->get_noeud1()->get_coord(),mtet->get_noeud2()->get_coord(),mtet->get_noeud3()->get_coord(),mtet->get_noeud4()->get_coord());
231 mtet->change_qualite(qual);
232 if (phase==0)
233 ajouter_ordre_tetra(mtet,1);
234 }
235
236
237 }
238
239 }
240 if ((phase==0) && (coquille_non_optimise==1))
241 ajouter_ordre_tetra(tet,1);
242 }
243
244 }
245 int nb=0;
246 //for (MG_TETRA* tet=mg_maillage->get_premier_tetra(it);tet;tet=mg_maillage->get_suivant_tetra(it))
247 if (mgvol!=NULL)
248 {
249 TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
250 for (MG_TETRA* tet=(MG_TETRA*)mgvol->get_lien_maillage()->get_premier(it);tet;tet=(MG_TETRA*)mgvol->get_lien_maillage()->get_suivant(it))
251 {
252 if (mg_maillage->get_mg_tetraid(tet->get_id())==NULL) continue;
253 M3D_TETRA* mtet=(M3D_TETRA*)tet;
254 if (mtet->get_qualite()<0.1*niveau_optimisation) nb++;
255 }
256 }
257 else
258 {
259 LISTE_MG_TETRA::iterator it;
260 for (MG_TETRA* tet=mg_maillage->get_premier_tetra(it);tet;tet=mg_maillage->get_suivant_tetra(it))
261 {
262 if (mg_maillage->get_mg_tetraid(tet->get_id())==NULL) continue;
263 M3D_TETRA* mtet=(M3D_TETRA*)tet;
264 if (mtet->get_qualite()<0.1*niveau_optimisation) nb++;
265 }
266 }
267
268 nbmauvais=nb;
269
270 }
271
272
273
274 int MAILLEUR3D_OPTIMISATION::bouge_point(MG_VOLUME* mgvol,MG_NOEUD* noeud,double& crit,double& x,double& y, double& z)
275 {
276 if (noeud->get_origine()!=MAGIC::ORIGINE::MAILLEUR_AUTO) return 0;
277 // if (noeud->get_type_entite()!=IDM3D_NOEUD) return 0;
278
279 if (mgvol!=NULL)
280 if (noeud->get_lien_topologie()!=mgvol) return 0;
281 if (mgvol==NULL)
282 if (noeud->get_type_entite()==MAGIC::TYPE_ENTITE::IDM3D_NOEUD)
283 if (((M3D_NOEUD*)noeud)->get_etat()==MAGIC::MAILLEURFRONTALETAT::INACTIF) return 0;
284 double xopt=0.;
285 double yopt=0.;
286 double zopt=0.;
287 double qual_dep=1.;
288 int nb_tet=noeud->get_lien_tetra()->get_nb();
289 double tab_coord[3000];
290 double gamma=0;
291 for (int i=0;i<nb_tet;i++)
292 {
293 M3D_TETRA* tet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
294 qual_dep=std::min(qual_dep,tet->get_qualite());
295 MG_NOEUD* no1=tet->get_noeud1();
296 MG_NOEUD* no2=tet->get_noeud2();
297 MG_NOEUD* no3=tet->get_noeud3();
298 MG_NOEUD* no4=tet->get_noeud4();
299 MG_NOEUD *noeud1,*noeud2,*noeud3;
300 if (no1==noeud)
301 {
302 noeud1=no2;
303 noeud2=no4;
304 noeud3=no3;
305 }
306 if (no2==noeud)
307 {
308 noeud1=no1;
309 noeud2=no3;
310 noeud3=no4;
311 }
312 if (no3==noeud)
313 {
314 noeud1=no1;
315 noeud2=no4;
316 noeud3=no2;
317 }
318 if (no4==noeud)
319 {
320
321 noeud1=no1;
322 noeud2=no2;
323 noeud3=no3;
324 }
325 double xyzi[3];
326 double *xyz1=noeud1->get_coord();
327 double *xyz2=noeud2->get_coord();
328 double *xyz3=noeud3->get_coord();
329 xyzi[0]=(xyz1[0]+xyz2[0]+xyz3[0])/3.;
330 xyzi[1]=(xyz1[1]+xyz2[1]+xyz3[1])/3.;
331 xyzi[2]=(xyz1[2]+xyz2[2]+xyz3[2])/3.;
332 OT_VECTEUR_3D v12(xyz1,xyz2);
333 OT_VECTEUR_3D v13(xyz1,xyz3);
334 OT_VECTEUR_3D v23(xyz2,xyz3);
335 OT_VECTEUR_3D normal=v12&v13;
336 double perimetre=v12.get_longueur()+v13.get_longueur()+v23.get_longueur();
337 double hauteur = (perimetre/3.) * 0.8 ;
338 normal.norme();
339 tab_coord[3*i]=xyzi[0]+normal.get_x()*hauteur;
340 tab_coord[3*i+1]=xyzi[1]+normal.get_y()*hauteur;
341 tab_coord[3*i+2]=xyzi[2]+normal.get_z()*hauteur;
342 if (tet->get_qualite()> 1e-10) gamma=gamma+1./tet->get_qualite()/tet->get_qualite();
343 }
344 gamma=1./gamma;
345 for (int i=0;i<nb_tet;i++)
346 {
347 M3D_TETRA* tet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
348 if (tet->get_qualite()> 1e-10)
349 {
350 double alpha=gamma/tet->get_qualite()/tet->get_qualite();
351 xopt=xopt+alpha*tab_coord[3*i];
352 yopt=yopt+alpha*tab_coord[3*i+1];
353 zopt=zopt+alpha*tab_coord[3*i+2];
354 }
355
356 }
357 double delta[3]={xopt-noeud->get_x(),yopt-noeud->get_y(),zopt-noeud->get_z()};
358 double bmin=0.,bmax=1.;
359 double vieuxx=noeud->get_x();
360 double vieuxy=noeud->get_y();
361 double vieuxz=noeud->get_z();
362 double qualini=0.;
363 for (int iteration=0;iteration<5;iteration++)
364 {
365 double alpha=(bmin+bmax)*0.5;
366 noeud->change_x(vieuxx+alpha*delta[0]);
367 noeud->change_y(vieuxy+alpha*delta[1]);
368 noeud->change_z(vieuxz+alpha*delta[2]);
369 double qualcour=1.;
370 for (int i=0;i<nb_tet;i++)
371 {
372 M3D_TETRA* mtet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
373 double qual=OPERATEUR::qualite_tetra(mtet->get_noeud1()->get_coord(),mtet->get_noeud2()->get_coord(),mtet->get_noeud3()->get_coord(),mtet->get_noeud4()->get_coord());
374 qualcour=std::min(qualcour,qual);
375 }
376 double alpha_eps=(bmin+bmax)*0.5-(bmax-bmin)/50.;
377 noeud->change_x(vieuxx+alpha_eps*delta[0]);
378 noeud->change_y(vieuxy+alpha_eps*delta[1]);
379 noeud->change_z(vieuxz+alpha_eps*delta[2]);
380 double qualcour_eps=1.;
381 for (int i=0;i<nb_tet;i++)
382 {
383 M3D_TETRA* mtet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
384 double qual=OPERATEUR::qualite_tetra(mtet->get_noeud1()->get_coord(),mtet->get_noeud2()->get_coord(),mtet->get_noeud3()->get_coord(),mtet->get_noeud4()->get_coord());
385 qualcour_eps=std::min(qualcour_eps,qual);
386 }
387 if (qualcour>qualcour_eps) bmin =alpha;
388 else bmax=alpha;
389 qualini=std::max(qualini,qualcour);
390 }
391 noeud->change_x(vieuxx);
392 noeud->change_y(vieuxy);
393 noeud->change_z(vieuxz);
394 if (qualini>qual_dep)
395 {
396 x=vieuxx+(bmin+bmax)*0.5*delta[0];
397 y=vieuxy+(bmin+bmax)*0.5*delta[1];
398 z=vieuxz+(bmin+bmax)*0.5*delta[2];
399 crit=qualini;
400 return 1;
401 }
402
403 return 0;
404 }
405
406
407
408 void MAILLEUR3D_OPTIMISATION::remaille_coquille(MG_NOEUD* noeud1,MG_NOEUD* noeud2, double& crit, COQUILLE& coque)
409 {
410 MG_NOEUD* noeud;
411 if (noeud1->get_id()>noeud2->get_id()) noeud=noeud2;
412 else noeud=noeud1;
413 int nb=noeud->get_lien_petit_segment()->get_nb();
414 for (int i=0;i<nb;i++)
415 {
416 MG_SEGMENT* seg=noeud->get_lien_petit_segment()->get(i);
417 if ( ((seg->get_noeud1()==noeud1) && (seg->get_noeud2()==noeud2)) || ((seg->get_noeud1()==noeud2) && (seg->get_noeud2()==noeud1)))
418 {
419 //if (mg_maillage->get_mg_geometrie()==NULL)
420 //{
421 if (seg->get_dimension_topo_null()==2)
422 {
423 coque.taille=0;
424 crit=0.;
425 return;
426 }
427 //}
428 //else
429 if (seg->get_lien_topologie()!=NULL) if ((seg->get_lien_topologie()->get_dimension()!=3) || (seg->get_origine()!=MAGIC::ORIGINE::MAILLEUR_AUTO))
430 {
431 coque.taille=0;
432 crit=0.;
433 return;
434 }
435 }
436 }
437 int nb1=noeud1->get_lien_tetra()->get_nb();
438 int nb2=noeud2->get_lien_tetra()->get_nb();
439 MG_TETRA* coq[100];
440 int nb_coq=0;
441 for (int i=0;i<nb1;i++)
442 for (int j=0;j<nb2;j++)
443 if (noeud1->get_lien_tetra()->get(i)==noeud2->get_lien_tetra()->get(j))
444 {
445 coq[nb_coq]=noeud1->get_lien_tetra()->get(i);
446 nb_coq++;
447 }
448 if ((nb_coq<4) || (nb_coq>10))
449 {
450 coque.taille=0;
451 crit=0.;
452 return;
453 }
454 MG_NOEUD* tet_noeud[50];
455 int nb_tet_noeud=0;
456 coque.taille=nb_coq;
457 for (int i=0;i<nb_coq;i++)
458 {
459 coque.tet[i]=coq[i];
460
461 if ((coq[i]->get_noeud1()!=noeud1) && (coq[i]->get_noeud1()!=noeud2)) tet_noeud[nb_tet_noeud++]=coq[i]->get_noeud1();
462 if ((coq[i]->get_noeud2()!=noeud1) && (coq[i]->get_noeud2()!=noeud2)) tet_noeud[nb_tet_noeud++]=coq[i]->get_noeud2();
463 if ((coq[i]->get_noeud3()!=noeud1) && (coq[i]->get_noeud3()!=noeud2)) tet_noeud[nb_tet_noeud++]=coq[i]->get_noeud3();
464 if ((coq[i]->get_noeud4()!=noeud1) && (coq[i]->get_noeud4()!=noeud2)) tet_noeud[nb_tet_noeud++]=coq[i]->get_noeud4();
465 coque.volume=coque.volume+fabs(get_volume(coq[i]));
466 }
467 MG_NOEUD* polygone[20];
468 int nb_poly=2;
469 polygone[0]=tet_noeud[0];
470 tet_noeud[0]=NULL;
471 polygone[1]=tet_noeud[1];
472 tet_noeud[1]=NULL;
473 while (nb_poly!=nb_coq+1)
474 {
475 for (int j=0;j<nb_coq;j++)
476 {
477 if (tet_noeud[2*j]==polygone[nb_poly-1])
478 {
479 polygone[nb_poly++]=tet_noeud[2*j+1];
480 tet_noeud[2*j]=NULL;
481 tet_noeud[2*j+1]=NULL;
482 }
483 if (tet_noeud[2*j+1]==polygone[nb_poly-1])
484 {
485 polygone[nb_poly++]=tet_noeud[2*j];
486 tet_noeud[2*j]=NULL;
487 tet_noeud[2*j+1]=NULL;
488 }
489
490 }
491 }
492 double vol1=get_volume(polygone[0]->get_coord(),polygone[1]->get_coord(),polygone[nb_poly-1]->get_coord(),noeud1->get_coord());
493 double vol2=get_volume(polygone[0]->get_coord(),polygone[nb_poly-1]->get_coord(),polygone[1]->get_coord(),noeud2->get_coord());
494 if (vol1*vol2<0.)
495 {
496 crit=0.;
497 coque.taille=0;
498 return;
499 }
500 MG_NOEUD *noeuda=noeud1;
501 MG_NOEUD *noeudb=noeud2;
502 if (vol1<0.)
503 {
504 noeuda=noeud2;
505 noeudb=noeud1;
506 }
507 int nb_solution,nb_triangle;
508
509 if (nb_coq==4) {
510 nb_solution=2;
511 nb_triangle=4;
512 }
513 if (nb_coq==5) {
514 nb_solution=5;
515 nb_triangle=10;
516 }
517 if (nb_coq==6) {
518 nb_solution=14;
519 nb_triangle=20;
520 }
521 if (nb_coq==7) {
522 nb_solution=42;
523 nb_triangle=35;
524 }
525 if (nb_coq==8) {
526 nb_solution=132;
527 nb_triangle=56;
528 }
529 if (nb_coq==9) {
530 nb_solution=429;
531 nb_triangle=84;
532 }
533 if (nb_coq==10) {
534 nb_solution=1430;
535 nb_triangle=120;
536 }
537 double crit_triangle[120];
538 for (int i=0;i<nb_triangle;i++)
539 {
540 double crit1=OPERATEUR::qualite_tetra(polygone[tab_face[nb_coq-4][i][0]]->get_coord(),polygone[tab_face[nb_coq-4][i][1]]->get_coord(),polygone[tab_face[nb_coq-4][i][2]]->get_coord(),noeuda->get_coord());
541 double crit2=OPERATEUR::qualite_tetra(polygone[tab_face[nb_coq-4][i][0]]->get_coord(),polygone[tab_face[nb_coq-4][i][2]]->get_coord(),polygone[tab_face[nb_coq-4][i][1]]->get_coord(),noeudb->get_coord());
542 crit_triangle[i]=std::min(crit1,crit2);
543 }
544 /* examen de chaque solution */
545 double crit_opt=0.;
546 int numero_solution= -1;
547 double crit_solution[1430];
548 for (int i=0;i<nb_solution;i++)
549 {
550 double volume=0.;
551 for (int j=0;j<nb_coq-2;j++)
552 {
553 MG_NOEUD* noa=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][i][j]][0]];
554 MG_NOEUD* nob=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][i][j]][1]];
555 MG_NOEUD* noc=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][i][j]][2]];
556 MG_NOEUD* nod=noeuda;
557 volume=volume+fabs(get_volume(noa->get_coord(),nob->get_coord(),noc->get_coord(),nod->get_coord()));
558 nod=noeudb;
559 volume=volume+fabs(get_volume(noa->get_coord(),nob->get_coord(),noc->get_coord(),nod->get_coord()));
560 }
561 double eps=0.0018*pow(volume,0.666666666);
562 if (OPERATEUR::egal(volume,coque.volume,eps))
563 {
564 crit_solution[i]=1.;
565 for (int j=0;j<nb_coq-2;j++)
566 crit_solution[i]=std::min(crit_solution[i],crit_triangle[tab_solution[nb_coq-4][i][j]]);
567 }
568 else crit_solution[i]=0.;
569 if (crit_opt<crit_solution[i])
570 {
571 crit_opt=crit_solution[i];
572 numero_solution=i;
573 }
574
575 }
576 if (numero_solution==(-1))
577 {
578 crit=0.;
579 coque.taille=0;
580 return;
581 }
582 crit=crit_opt;
583 for (int j=0;j<nb_coq-2;j++)
584 {
585 MG_NOEUD* noa=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][numero_solution][j]][0]];
586 MG_NOEUD* nob=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][numero_solution][j]][1]];
587 MG_NOEUD* noc=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][numero_solution][j]][2]];
588 coque.new_tetra[8*j]=noa;
589 coque.new_tetra[8*j+1]=nob;
590 coque.new_tetra[8*j+2]=noc;
591 coque.new_tetra[8*j+3]=noeuda;
592 coque.new_tetra[8*j+4]=noa;
593 coque.new_tetra[8*j+5]=noc;
594 coque.new_tetra[8*j+6]=nob;
595 coque.new_tetra[8*j+7]=noeudb;
596 }
597 }
598
599
600 void MAILLEUR3D_OPTIMISATION::ajouter_ordre_tetra(M3D_TETRA* tet,int num)
601 {
602 double val;
603 if (num==0) val=tet->get_qualite();
604 if (num==1) val=1./tet->get_qualite();
605 std::pair<double,M3D_TETRA*> tmp(val,tet);
606 ORDRE_TETRA::iterator p=lst_tetra[num].insert(tmp);
607 lst_tetraid[num][tet->get_id()]=p;
608 }
609
610 void MAILLEUR3D_OPTIMISATION::supprimer_ordre_tetra(M3D_TETRA* tet)
611 {
612 int num=0;
613 ORDRE_TETRA_PARID::iterator it=lst_tetraid[num].find(tet->get_id());
614 if (it==lst_tetraid[num].end())
615 {
616 num=1;
617 it=lst_tetraid[num].find(tet->get_id());
618 }
619 if (it==lst_tetraid[num].end()) return;
620 lst_tetra[num].erase(it->second);
621 lst_tetraid[num].erase(it);
622 }
623
624
625 void MAILLEUR3D_OPTIMISATION::change_niveau_optimisation(int num)
626 {
627 niveau_optimisation=num;
628 }
629
630 int MAILLEUR3D_OPTIMISATION::get_niveau_optimisation(void)
631 {
632 return niveau_optimisation;
633 }
634
635
636
637 class MG_TRIANGLE* MAILLEUR3D_OPTIMISATION::insere_triangle(class MG_VOLUME* mgvol,MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,int origine)
638 {
639 MG_TRIANGLE* triangle=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud3->get_id(),noeud2->get_id());
640 if (triangle!=NULL) return NULL;
641 MG_SEGMENT* segment1=mg_maillage->get_mg_segment(noeud1->get_id(),noeud2->get_id());
642 MG_SEGMENT* segment2=mg_maillage->get_mg_segment(noeud1->get_id(),noeud3->get_id());
643 MG_SEGMENT* segment3=mg_maillage->get_mg_segment(noeud2->get_id(),noeud3->get_id());
644 if (segment1==NULL) segment1=cree_segment(mgvol,noeud1,noeud2,origine);
645 if (segment2==NULL) segment2=cree_segment(mgvol,noeud1,noeud3,origine);
646 if (segment3==NULL) segment3=cree_segment(mgvol,noeud2,noeud3,origine);
647 MG_TRIANGLE* tri=cree_triangle(mgvol,noeud1,noeud2,noeud3,segment1,segment2,segment3,origine);
648 return tri;
649 }
650
651
652
653
654
655
656 MG_TRIANGLE* MAILLEUR3D_OPTIMISATION::cree_triangle(class MG_VOLUME* mgvol,MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_SEGMENT* segment1,MG_SEGMENT* segment2,MG_SEGMENT* segment3,int origine)
657 {
658 M3D_TRIANGLE* tri=new M3D_TRIANGLE(mgvol,noeud1,noeud2,noeud3,segment1,segment2,segment3,origine);
659 mg_maillage->ajouter_mg_triangle(tri);
660 return tri;
661 }
662
663
664
665 MG_SEGMENT* MAILLEUR3D_OPTIMISATION::cree_segment(class MG_VOLUME* mgvol,MG_NOEUD* noeud1,MG_NOEUD* noeud2,int origine)
666 {
667 MG_SEGMENT* seg=mg_maillage->ajouter_mg_segment(mgvol,noeud1,noeud2,origine);
668 seg->change_dimension_topo_null(3);
669 return seg;
670 }