ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur3d_optimisation.cpp
Revision: 1152
Committed: Mon Jun 3 20:21:02 2024 UTC (11 months, 1 week ago) by ghazal
Original Path: magic/lib/mailleur_auto/src/mailleur3d_optimisation.cpp
File size: 52779 byte(s)
Log Message:
Optimisateur3D pour geometrie virtuelle
quadratisation pour geometrie virtuelle

File Contents

# User Rev Content
1 francois 283 //------------------------------------------------------------
2     //------------------------------------------------------------
3     // MAGiC
4     // Jean Christophe Cuilli�e et Vincent FRANCOIS
5     // D�artement de G�ie M�anique - UQTR
6     //------------------------------------------------------------
7     // Le projet MAGIC est un projet de recherche du d�artement
8     // de g�ie m�anique de l'Universit�du Qu�ec �
9     // Trois Rivi�es
10     // Les librairies ne peuvent �re utilis�s sans l'accord
11     // des auteurs (contact : francois@uqtr.ca)
12     //------------------------------------------------------------
13     //------------------------------------------------------------
14     //
15     // mailleur3d.cpp
16     //
17     //------------------------------------------------------------
18     //------------------------------------------------------------
19     // COPYRIGHT 2000
20     // Version du 02/03/2006 �11H23
21     //------------------------------------------------------------
22     //------------------------------------------------------------
23    
24    
25     #include "gestionversion.h"
26 francois 287 #include "mailleur3d_optimisation.h"
27 francois 283 #include "ot_mathematique.h"
28     #include "m3d_tetra.h"
29 francois 287 #include "m3d_triangle.h"
30 francois 283 #include "m3d_noeud.h"
31 francois 287 #include "mg_maillage.h"
32 francois 283 //#include "mg_gestionnaire.h"
33     #include <math.h>
34    
35     //---------------------------------------------------------------------------
36    
37     //#pragma package(smart_init)
38    
39 francois 287
40 francois 494 MAILLEUR3D_OPTIMISATION::MAILLEUR3D_OPTIMISATION(MG_MAILLAGE* mgmai,int niv):MAILLEUR(false),mg_maillage(mgmai),niveau_optimisation(niv)
41 francois 283 {
42 francois 287 o3d_data();
43     o3d_data2();
44 francois 1137 LISTE_MG_TETRA::iterator it;
45     TPL_MAP_ENTITE<M3D_TETRA*> listetet3d;
46     for (MG_TETRA* tet=mg_maillage->get_premier_tetra(it);tet!=NULL;tet=mg_maillage->get_suivant_tetra(it))
47 francois 1150 if (tet->get_type_entite()!=MAGIC::TYPE_ENTITE::IDM3D_TETRA)
48 francois 1137 {
49     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());
50     listetet3d.ajouter(tet3d);
51     }
52    
53     TPL_MAP_ENTITE<M3D_TETRA*>::ITERATEUR it2;
54    
55     for (M3D_TETRA* tet=listetet3d.get_premier(it2);tet!=NULL;tet=listetet3d.get_suivant(it2))
56     {
57     mg_maillage->supprimer_mg_tetraid(tet->get_id());
58     mg_maillage->ajouter_mg_tetra(tet);
59     }
60 francois 287 }
61    
62     MAILLEUR3D_OPTIMISATION::~MAILLEUR3D_OPTIMISATION()
63     {
64     }
65    
66    
67    
68     double MAILLEUR3D_OPTIMISATION::get_volume(MG_TETRA* tet)
69     {
70 francois 283 double *xyz1=tet->get_noeud1()->get_coord();
71     double *xyz2=tet->get_noeud2()->get_coord();
72     double *xyz3=tet->get_noeud3()->get_coord();
73     double *xyz4=tet->get_noeud4()->get_coord();
74     OT_VECTEUR_3D ab(xyz1,xyz2);
75     OT_VECTEUR_3D ac(xyz1,xyz3);
76     OT_VECTEUR_3D ad(xyz1,xyz4);
77     double vol=(ab&ac)*ad;
78     vol=vol/6.;
79     return vol;
80     }
81    
82 francois 287 double MAILLEUR3D_OPTIMISATION::get_volume(double *xyz1,double *xyz2,double *xyz3,double *xyz4)
83 francois 283 {
84     OT_VECTEUR_3D ab(xyz1,xyz2);
85     OT_VECTEUR_3D ac(xyz1,xyz3);
86     OT_VECTEUR_3D ad(xyz1,xyz4);
87     double vol=(ab&ac)*ad;
88     vol=vol/6.;
89     return vol;
90     }
91    
92 francois 288 void MAILLEUR3D_OPTIMISATION::optimise(MG_VOLUME* mgvol)
93     {
94     char mess[255];
95     sprintf(mess," Niveau d'optimisation : %d",niveau_optimisation);
96     if (affichageactif==1) affiche(mess);
97     int nbaoptimiser;
98     int nbaoptimiserapres=mg_maillage->get_nb_mg_tetra();
99     do
100     {
101     nbaoptimiser=nbaoptimiserapres;
102     optimise(mgvol,nbaoptimiserapres);
103     if (nbaoptimiserapres==0)
104     {
105     nbaoptimiser=0;
106     sprintf(mess," Tout est optimise");
107     if (affichageactif==1) affiche(mess);
108     }
109    
110     }
111     while (nbaoptimiserapres!=nbaoptimiser);
112     }
113    
114    
115    
116 francois 287 void MAILLEUR3D_OPTIMISATION::optimise(MG_VOLUME* mgvol,int& nbmauvais)
117 francois 283 {
118 francois 603 int passe=1;
119     //passe++;
120 francois 283 //preparation
121     int nbtet=0;
122 francois 420 if (mgvol!=NULL)
123 francois 283 {
124 francois 420 TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
125     for (MG_TETRA* tet=(MG_TETRA*)mgvol->get_lien_maillage()->get_premier(it);tet;tet=(MG_TETRA*)mgvol->get_lien_maillage()->get_suivant(it))
126     {
127 francois 339 if (mg_maillage->get_mg_tetraid(tet->get_id())==NULL) continue;
128     M3D_TETRA* mtet=(M3D_TETRA*)tet;
129 francois 283 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());
130     mtet->change_qualite(qual);
131     if (qual<0.1*niveau_optimisation) nbtet++;
132     ajouter_ordre_tetra(mtet,0);
133 francois 420 }
134 francois 283 }
135 francois 420 else
136     {
137     LISTE_MG_TETRA::iterator it;
138     for (MG_TETRA* tet=mg_maillage->get_premier_tetra(it);tet;tet=mg_maillage->get_suivant_tetra(it))
139     {
140     if (mg_maillage->get_mg_tetraid(tet->get_id())==NULL) continue;
141     M3D_TETRA* mtet=(M3D_TETRA*)tet;
142     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());
143     mtet->change_qualite(qual);
144     if (qual<0.1*niveau_optimisation) nbtet++;
145     ajouter_ordre_tetra(mtet,0);
146     }
147     }
148    
149 francois 283 //on commence
150 francois 339 if (nbtet==0) return;
151 francois 283 for (int phase=0;phase<2;phase++)
152     {
153     char mess[255];
154     if (phase==0) sprintf(mess," Phase %d : %d a optimiser",2*passe-1+phase,nbtet);
155     if (phase==1) sprintf(mess," Phase %d",2*passe-1+phase);
156     if (affichageactif==1) affiche(mess);
157     while (lst_tetra[phase].size()>0)
158     {
159 francois 830 ORDRE_TETRA::iterator i=lst_tetra[phase].begin();
160 francois 283 M3D_TETRA* tet=(*i).second;
161     supprimer_ordre_tetra(tet);
162     double crit_avant=tet->get_qualite();
163 francois 339 if (crit_avant<0.1*niveau_optimisation)
164 francois 283 {
165     double crit[4]={-1.,-1.,-1.,-1.};
166 ghazal 1152 MG_NOEUD* no[4];
167     no[0]=(MG_NOEUD*)tet->get_noeud1();
168     no[1]=(MG_NOEUD*)tet->get_noeud2();
169     no[2]=(MG_NOEUD*)tet->get_noeud3();
170     no[3]=(MG_NOEUD*)tet->get_noeud4();
171 francois 283 double x[4],y[4],z[4];
172     bouge_point(mgvol,no[0],crit[0],x[0],y[0],z[0]);
173     bouge_point(mgvol,no[1],crit[1],x[1],y[1],z[1]);
174     bouge_point(mgvol,no[2],crit[2],x[2],y[2],z[2]);
175     bouge_point(mgvol,no[3],crit[3],x[3],y[3],z[3]);
176     std::vector<double> vecteur_crit(crit,crit+4);
177     std::vector<double>::iterator it1 = max_element(vecteur_crit.begin(), vecteur_crit.end());
178     double crit_opt=*it1;
179     if (crit_opt>crit_avant)
180     {
181     int num=it1-vecteur_crit.begin();
182     no[num]->change_x(x[num]);
183     no[num]->change_y(y[num]);
184     no[num]->change_z(z[num]);
185     int nb=no[num]->get_lien_tetra()->get_nb();
186     for (int j=0;j<nb;j++)
187     {
188     M3D_TETRA* mtet=(M3D_TETRA*)no[num]->get_lien_tetra()->get(j);
189     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());
190     mtet->change_qualite(qual);
191     }
192     crit_avant=crit_opt;
193     }
194     }
195 ghazal 1152
196 francois 283 int coquille_non_optimise=1;
197 francois 339 if (crit_avant<0.1*niveau_optimisation)
198 francois 283 {
199     COQUILLE coque[6];
200     double crit[6]={-1,-1,-1,-1,-1,-1};
201     remaille_coquille(tet->get_noeud1(),tet->get_noeud2(),crit[0],coque[0]);
202     remaille_coquille(tet->get_noeud1(),tet->get_noeud3(),crit[1],coque[1]);
203     remaille_coquille(tet->get_noeud1(),tet->get_noeud4(),crit[2],coque[2]);
204     remaille_coquille(tet->get_noeud2(),tet->get_noeud3(),crit[3],coque[3]);
205     remaille_coquille(tet->get_noeud2(),tet->get_noeud3(),crit[4],coque[4]);
206     remaille_coquille(tet->get_noeud3(),tet->get_noeud4(),crit[5],coque[5]);
207     std::vector<double> vecteur_crit(crit,crit+6);
208     std::vector<double>::iterator it1 = max_element(vecteur_crit.begin(), vecteur_crit.end());
209     double crit_opt=*it1;
210     if (crit_opt>crit_avant)
211     {
212     coquille_non_optimise=0;
213     int num=it1-vecteur_crit.begin();
214     for (int i=0;i<coque[num].taille;i++)
215     {
216     M3D_TETRA *tettmp=(M3D_TETRA*)coque[num].tet[i];
217     if (tet!=tettmp) supprimer_ordre_tetra(tettmp);
218     mg_maillage->supprimer_mg_tetraid(coque[num].tet[i]->get_id());
219     }
220     for (int i=0;i<2*coque[num].taille-4;i++)
221     {
222     MG_NOEUD* noeud1=coque[num].new_tetra[4*i];
223     MG_NOEUD* noeud2=coque[num].new_tetra[4*i+1];
224     MG_NOEUD* noeud3=coque[num].new_tetra[4*i+2];
225     MG_NOEUD* noeud4=coque[num].new_tetra[4*i+3];
226     MG_TRIANGLE* triangle1=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud3->get_id(),noeud2->get_id());
227     MG_TRIANGLE* triangle2=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud2->get_id(),noeud4->get_id());
228     MG_TRIANGLE* triangle3=mg_maillage->get_mg_triangle(noeud2->get_id(),noeud3->get_id(),noeud4->get_id());
229     MG_TRIANGLE* triangle4=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud4->get_id(),noeud3->get_id());
230 francois 791 if (triangle1==NULL) triangle1=insere_triangle(mgvol,noeud1,noeud3,noeud2,MAGIC::ORIGINE::MAILLEUR_AUTO);
231     if (triangle2==NULL) triangle2=insere_triangle(mgvol,noeud1,noeud2,noeud4,MAGIC::ORIGINE::MAILLEUR_AUTO);
232     if (triangle3==NULL) triangle3=insere_triangle(mgvol,noeud2,noeud3,noeud4,MAGIC::ORIGINE::MAILLEUR_AUTO);
233     if (triangle4==NULL) triangle4=insere_triangle(mgvol,noeud1,noeud4,noeud3,MAGIC::ORIGINE::MAILLEUR_AUTO);
234     M3D_TETRA* mtet=new M3D_TETRA(mgvol,noeud1,noeud2,noeud3,noeud4,triangle1,triangle2,triangle3,triangle4,MAGIC::ORIGINE::MAILLEUR_AUTO);
235 francois 283 mg_maillage->ajouter_mg_tetra(mtet);
236     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());
237     mtet->change_qualite(qual);
238     if (phase==0)
239     ajouter_ordre_tetra(mtet,1);
240     }
241    
242    
243     }
244    
245     }
246     if ((phase==0) && (coquille_non_optimise==1))
247     ajouter_ordre_tetra(tet,1);
248     }
249    
250     }
251     int nb=0;
252 francois 339 //for (MG_TETRA* tet=mg_maillage->get_premier_tetra(it);tet;tet=mg_maillage->get_suivant_tetra(it))
253 francois 420 if (mgvol!=NULL)
254 francois 283 {
255 francois 420 TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
256     for (MG_TETRA* tet=(MG_TETRA*)mgvol->get_lien_maillage()->get_premier(it);tet;tet=(MG_TETRA*)mgvol->get_lien_maillage()->get_suivant(it))
257     {
258 francois 339 if (mg_maillage->get_mg_tetraid(tet->get_id())==NULL) continue;
259 francois 283 M3D_TETRA* mtet=(M3D_TETRA*)tet;
260     if (mtet->get_qualite()<0.1*niveau_optimisation) nb++;
261 francois 420 }
262 francois 283 }
263 francois 420 else
264     {
265     LISTE_MG_TETRA::iterator it;
266     for (MG_TETRA* tet=mg_maillage->get_premier_tetra(it);tet;tet=mg_maillage->get_suivant_tetra(it))
267     {
268     if (mg_maillage->get_mg_tetraid(tet->get_id())==NULL) continue;
269     M3D_TETRA* mtet=(M3D_TETRA*)tet;
270     if (mtet->get_qualite()<0.1*niveau_optimisation) nb++;
271     }
272     }
273 francois 283
274     nbmauvais=nb;
275    
276     }
277    
278    
279    
280 ghazal 1152 int MAILLEUR3D_OPTIMISATION::bouge_point(MG_VOLUME* mgvol,MG_NOEUD* noeud,double& crit,double& x,double& y, double& z)
281 francois 283 {
282 francois 791 if (noeud->get_origine()!=MAGIC::ORIGINE::MAILLEUR_AUTO) return 0;
283 francois 558 // if (noeud->get_type_entite()!=IDM3D_NOEUD) return 0;
284 ghazal 1152
285     if (mgvol!=NULL)
286 francois 288 if (noeud->get_lien_topologie()!=mgvol) return 0;
287     if (mgvol==NULL)
288 ghazal 1152 if (noeud->get_type_entite()==MAGIC::TYPE_ENTITE::IDM3D_NOEUD)
289     if (((M3D_NOEUD*)noeud)->get_etat()==MAGIC::MAILLEURFRONTALETAT::INACTIF) return 0;
290 francois 283 double xopt=0.;
291     double yopt=0.;
292     double zopt=0.;
293     double qual_dep=1.;
294     int nb_tet=noeud->get_lien_tetra()->get_nb();
295     double tab_coord[3000];
296     double gamma=0;
297     for (int i=0;i<nb_tet;i++)
298     {
299     M3D_TETRA* tet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
300     qual_dep=std::min(qual_dep,tet->get_qualite());
301     MG_NOEUD* no1=tet->get_noeud1();
302     MG_NOEUD* no2=tet->get_noeud2();
303     MG_NOEUD* no3=tet->get_noeud3();
304     MG_NOEUD* no4=tet->get_noeud4();
305     MG_NOEUD *noeud1,*noeud2,*noeud3;
306     if (no1==noeud)
307     {
308     noeud1=no2;
309     noeud2=no4;
310     noeud3=no3;
311     }
312     if (no2==noeud)
313     {
314     noeud1=no1;
315     noeud2=no3;
316     noeud3=no4;
317     }
318     if (no3==noeud)
319     {
320     noeud1=no1;
321     noeud2=no4;
322     noeud3=no2;
323     }
324     if (no4==noeud)
325     {
326 ghazal 1152
327 francois 283 noeud1=no1;
328     noeud2=no2;
329     noeud3=no3;
330     }
331     double xyzi[3];
332     double *xyz1=noeud1->get_coord();
333     double *xyz2=noeud2->get_coord();
334     double *xyz3=noeud3->get_coord();
335     xyzi[0]=(xyz1[0]+xyz2[0]+xyz3[0])/3.;
336     xyzi[1]=(xyz1[1]+xyz2[1]+xyz3[1])/3.;
337     xyzi[2]=(xyz1[2]+xyz2[2]+xyz3[2])/3.;
338     OT_VECTEUR_3D v12(xyz1,xyz2);
339     OT_VECTEUR_3D v13(xyz1,xyz3);
340     OT_VECTEUR_3D v23(xyz2,xyz3);
341     OT_VECTEUR_3D normal=v12&v13;
342     double perimetre=v12.get_longueur()+v13.get_longueur()+v23.get_longueur();
343     double hauteur = (perimetre/3.) * 0.8 ;
344     normal.norme();
345     tab_coord[3*i]=xyzi[0]+normal.get_x()*hauteur;
346     tab_coord[3*i+1]=xyzi[1]+normal.get_y()*hauteur;
347     tab_coord[3*i+2]=xyzi[2]+normal.get_z()*hauteur;
348     if (tet->get_qualite()> 1e-10) gamma=gamma+1./tet->get_qualite()/tet->get_qualite();
349     }
350     gamma=1./gamma;
351     for (int i=0;i<nb_tet;i++)
352     {
353     M3D_TETRA* tet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
354     if (tet->get_qualite()> 1e-10)
355     {
356     double alpha=gamma/tet->get_qualite()/tet->get_qualite();
357     xopt=xopt+alpha*tab_coord[3*i];
358     yopt=yopt+alpha*tab_coord[3*i+1];
359     zopt=zopt+alpha*tab_coord[3*i+2];
360     }
361    
362     }
363     double delta[3]={xopt-noeud->get_x(),yopt-noeud->get_y(),zopt-noeud->get_z()};
364     double bmin=0.,bmax=1.;
365     double vieuxx=noeud->get_x();
366     double vieuxy=noeud->get_y();
367     double vieuxz=noeud->get_z();
368     double qualini=0.;
369     for (int iteration=0;iteration<5;iteration++)
370     {
371     double alpha=(bmin+bmax)*0.5;
372     noeud->change_x(vieuxx+alpha*delta[0]);
373     noeud->change_y(vieuxy+alpha*delta[1]);
374     noeud->change_z(vieuxz+alpha*delta[2]);
375     double qualcour=1.;
376     for (int i=0;i<nb_tet;i++)
377     {
378     M3D_TETRA* mtet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
379     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());
380     qualcour=std::min(qualcour,qual);
381     }
382     double alpha_eps=(bmin+bmax)*0.5-(bmax-bmin)/50.;
383     noeud->change_x(vieuxx+alpha_eps*delta[0]);
384     noeud->change_y(vieuxy+alpha_eps*delta[1]);
385     noeud->change_z(vieuxz+alpha_eps*delta[2]);
386     double qualcour_eps=1.;
387     for (int i=0;i<nb_tet;i++)
388     {
389     M3D_TETRA* mtet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
390     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());
391     qualcour_eps=std::min(qualcour_eps,qual);
392     }
393     if (qualcour>qualcour_eps) bmin =alpha;
394     else bmax=alpha;
395     qualini=std::max(qualini,qualcour);
396     }
397     noeud->change_x(vieuxx);
398     noeud->change_y(vieuxy);
399     noeud->change_z(vieuxz);
400     if (qualini>qual_dep)
401     {
402     x=vieuxx+(bmin+bmax)*0.5*delta[0];
403     y=vieuxy+(bmin+bmax)*0.5*delta[1];
404     z=vieuxz+(bmin+bmax)*0.5*delta[2];
405     crit=qualini;
406     return 1;
407     }
408 ghazal 1152
409 francois 283 return 0;
410     }
411    
412 ghazal 1152 #ifdef TOTOVAALAPLAGE
413 francois 283
414 ghazal 1152 void MAILLEUR3D_OPTIMISATION::optimise_virtu(MG_VOLUME* mgvol)
415     {
416     char mess[255];
417     sprintf(mess," Niveau d'optimisation : %d",niveau_optimisation);
418     if (affichageactif==1) affiche(mess);
419     int nbaoptimiser;
420     int nbaoptimiserapres=mg_maillage->get_nb_mg_tetra();
421     do
422     {
423     nbaoptimiser=nbaoptimiserapres;
424     optimise_virtu(mgvol,nbaoptimiserapres);
425     if (nbaoptimiserapres==0)
426     {
427     nbaoptimiser=0;
428     sprintf(mess," Tout est optimise");
429     if (affichageactif==1) affiche(mess);
430     }
431    
432     }
433     while (nbaoptimiserapres!=nbaoptimiser);
434     }
435 francois 283
436 ghazal 1152 void MAILLEUR3D_OPTIMISATION::optimise_virtu(MG_VOLUME* mgvol,int& nbmauvais)
437     {
438     int passe=1;
439     //passe++;
440     //preparation
441     int nbtet=0;
442     if (mgvol!=NULL)
443     {
444     TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
445     for (MG_TETRA* tet=(MG_TETRA*)mgvol->get_lien_maillage()->get_premier(it);tet;tet=(MG_TETRA*)mgvol->get_lien_maillage()->get_suivant(it))
446     {
447     if (mg_maillage->get_mg_tetraid(tet->get_id())==NULL) continue;
448     M3D_TETRA* mtet=(M3D_TETRA*)tet;
449     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());
450     mtet->change_qualite(qual);
451     if (qual<0.1*niveau_optimisation) nbtet++;
452     ajouter_ordre_tetra(mtet,0);
453     }
454     }
455     else
456     {
457     LISTE_MG_TETRA::iterator it;
458     for (MG_TETRA* tet=mg_maillage->get_premier_tetra(it);tet;tet=mg_maillage->get_suivant_tetra(it))
459     {
460     if (mg_maillage->get_mg_tetraid(tet->get_id())==NULL) continue;
461     M3D_TETRA* mtet=(M3D_TETRA*)tet;
462     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());
463     mtet->change_qualite(qual);
464     if (qual<0.1*niveau_optimisation) nbtet++;
465     ajouter_ordre_tetra(mtet,0);
466     }
467     }
468    
469     //on commence
470     if (nbtet==0) return;
471     for (int phase=0;phase<2;phase++)
472     {
473     char mess[255];
474     if (phase==0) sprintf(mess," Phase %d : %d a optimiser",2*passe-1+phase,nbtet);
475     if (phase==1) sprintf(mess," Phase %d",2*passe-1+phase);
476     if (affichageactif==1) affiche(mess);
477     while (lst_tetra[phase].size()>0)
478     {
479     ORDRE_TETRA::iterator i=lst_tetra[phase].begin();
480     M3D_TETRA* tet=(*i).second;
481     supprimer_ordre_tetra(tet);
482     double crit_avant=tet->get_qualite();
483     if (crit_avant<0.1*niveau_optimisation)
484     {
485     double crit[4]={-1.,-1.,-1.,-1.};
486     MG_NOEUD* no[4];
487     no[0]=(MG_NOEUD*)tet->get_noeud1();
488     no[1]=(MG_NOEUD*)tet->get_noeud2();
489     no[2]=(MG_NOEUD*)tet->get_noeud3();
490     no[3]=(MG_NOEUD*)tet->get_noeud4();
491     double x[4],y[4],z[4];
492     bouge_point_virtu(mgvol,no[0],crit[0],x[0],y[0],z[0]);
493     bouge_point_virtu(mgvol,no[1],crit[1],x[1],y[1],z[1]);
494     bouge_point_virtu(mgvol,no[2],crit[2],x[2],y[2],z[2]);
495     bouge_point_virtu(mgvol,no[3],crit[3],x[3],y[3],z[3]);
496     std::vector<double> vecteur_crit(crit,crit+4);
497     std::vector<double>::iterator it1 = max_element(vecteur_crit.begin(), vecteur_crit.end());
498     double crit_opt=*it1;
499     if (crit_opt>crit_avant)
500     {
501     int num=it1-vecteur_crit.begin();
502     no[num]->change_x(x[num]);
503     no[num]->change_y(y[num]);
504     no[num]->change_z(z[num]);
505     int nb=no[num]->get_lien_tetra()->get_nb();
506     for (int j=0;j<nb;j++)
507     {
508     M3D_TETRA* mtet=(M3D_TETRA*)no[num]->get_lien_tetra()->get(j);
509     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());
510     // mtet->change_qualite(qual);
511     }
512     crit_avant=crit_opt;
513     }
514     }
515    
516     int coquille_non_optimise=1;
517     if (crit_avant<0.1*niveau_optimisation)
518     {
519     COQUILLE coque[6];
520     double crit[6]={-1,-1,-1,-1,-1,-1};
521     remaille_coquille(tet->get_noeud1(),tet->get_noeud2(),crit[0],coque[0]);
522     remaille_coquille(tet->get_noeud1(),tet->get_noeud3(),crit[1],coque[1]);
523     remaille_coquille(tet->get_noeud1(),tet->get_noeud4(),crit[2],coque[2]);
524     remaille_coquille(tet->get_noeud2(),tet->get_noeud3(),crit[3],coque[3]);
525     remaille_coquille(tet->get_noeud2(),tet->get_noeud3(),crit[4],coque[4]);
526     remaille_coquille(tet->get_noeud3(),tet->get_noeud4(),crit[5],coque[5]);
527     std::vector<double> vecteur_crit(crit,crit+6);
528     std::vector<double>::iterator it1 = max_element(vecteur_crit.begin(), vecteur_crit.end());
529     double crit_opt=*it1;
530     if (crit_opt>crit_avant)
531     {
532     coquille_non_optimise=0;
533     int num=it1-vecteur_crit.begin();
534     for (int i=0;i<coque[num].taille;i++)
535     {
536     M3D_TETRA *tettmp=(M3D_TETRA*)coque[num].tet[i];
537     if (tet!=tettmp) supprimer_ordre_tetra(tettmp);
538     mg_maillage->supprimer_mg_tetraid(coque[num].tet[i]->get_id());
539     }
540     for (int i=0;i<2*coque[num].taille-4;i++)
541     {
542     MG_NOEUD* noeud1=coque[num].new_tetra[4*i];
543     MG_NOEUD* noeud2=coque[num].new_tetra[4*i+1];
544     MG_NOEUD* noeud3=coque[num].new_tetra[4*i+2];
545     MG_NOEUD* noeud4=coque[num].new_tetra[4*i+3];
546     MG_TRIANGLE* triangle1=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud3->get_id(),noeud2->get_id());
547     MG_TRIANGLE* triangle2=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud2->get_id(),noeud4->get_id());
548     MG_TRIANGLE* triangle3=mg_maillage->get_mg_triangle(noeud2->get_id(),noeud3->get_id(),noeud4->get_id());
549     MG_TRIANGLE* triangle4=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud4->get_id(),noeud3->get_id());
550     if (triangle1==NULL) triangle1=insere_triangle(mgvol,noeud1,noeud3,noeud2,MAGIC::ORIGINE::MAILLEUR_AUTO);
551     if (triangle2==NULL) triangle2=insere_triangle(mgvol,noeud1,noeud2,noeud4,MAGIC::ORIGINE::MAILLEUR_AUTO);
552     if (triangle3==NULL) triangle3=insere_triangle(mgvol,noeud2,noeud3,noeud4,MAGIC::ORIGINE::MAILLEUR_AUTO);
553     if (triangle4==NULL) triangle4=insere_triangle(mgvol,noeud1,noeud4,noeud3,MAGIC::ORIGINE::MAILLEUR_AUTO);
554     M3D_TETRA* mtet=new M3D_TETRA(mgvol,noeud1,noeud2,noeud3,noeud4,triangle1,triangle2,triangle3,triangle4,MAGIC::ORIGINE::MAILLEUR_AUTO);
555     mg_maillage->ajouter_mg_tetra(mtet);
556     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());
557     mtet->change_qualite(qual);
558     if (phase==0)
559     ajouter_ordre_tetra(mtet,1);
560     }
561    
562    
563     }
564    
565     }
566     if ((phase==0) && (coquille_non_optimise==1))
567     ajouter_ordre_tetra(tet,1);
568     }
569    
570     }
571     int nb=0;
572     //for (MG_TETRA* tet=mg_maillage->get_premier_tetra(it);tet;tet=mg_maillage->get_suivant_tetra(it))
573     if (mgvol!=NULL)
574     {
575     TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
576     for (MG_TETRA* tet=(MG_TETRA*)mgvol->get_lien_maillage()->get_premier(it);tet;tet=(MG_TETRA*)mgvol->get_lien_maillage()->get_suivant(it))
577     {
578     if (mg_maillage->get_mg_tetraid(tet->get_id())==NULL) continue;
579     M3D_TETRA* mtet=(M3D_TETRA*)tet;
580     if (mtet->get_qualite()<0.1*niveau_optimisation) nb++;
581     }
582     }
583     else
584     {
585     LISTE_MG_TETRA::iterator it;
586     for (MG_TETRA* tet=mg_maillage->get_premier_tetra(it);tet;tet=mg_maillage->get_suivant_tetra(it))
587     {
588     if (mg_maillage->get_mg_tetraid(tet->get_id())==NULL) continue;
589     M3D_TETRA* mtet=(M3D_TETRA*)tet;
590     if (mtet->get_qualite()<0.1*niveau_optimisation) nb++;
591     }
592     }
593    
594     nbmauvais=nb;
595    
596     }
597     int MAILLEUR3D_OPTIMISATION::bouge_point_virtu(MG_VOLUME* mgvol,MG_NOEUD* noeud,double& crit,double& x,double& y, double& z)
598     {
599     if (noeud->get_origine()!=MAGIC::ORIGINE::MAILLEUR_AUTO) return 0;
600     // if (noeud->get_type_entite()!=IDM3D_NOEUD) return 0;
601    
602    
603     //===================GHAZAL=============================================================
604    
605     if( noeud->get_lien_topologie()->get_type()==MG_ELEMENT_TOPOLOGIQUE::TYPE_ELEMENT_TOPOLOGIQUE::FACE)
606     {
607     // printf("On face\n\n");
608    
609     double x_point=noeud->get_lien_topologie()->get_contrainte(0);
610     double y_point=noeud->get_lien_topologie()->get_contrainte(1);
611     double z_point=noeud->get_lien_topologie()->get_contrainte(2);
612    
613     double x_norm1=noeud->get_lien_topologie()->get_contrainte(3);
614     double y_norm1=noeud->get_lien_topologie()->get_contrainte(4);
615     double z_norm1=noeud->get_lien_topologie()->get_contrainte(5);
616    
617     double x_norm2=noeud->get_lien_topologie()->get_contrainte(6);
618     double y_norm2=noeud->get_lien_topologie()->get_contrainte(7);
619     double z_norm2=noeud->get_lien_topologie()->get_contrainte(8);
620    
621     OT_VECTEUR_3D n1(x_norm1,y_norm1,z_norm1);
622     OT_VECTEUR_3D n2(x_norm2,y_norm2,z_norm2);
623    
624    
625     double P1[3]={x_point,y_point,z_point};
626    
627     OT_VECTEUR_3D normal_plane=n1&n2;
628    
629     double xopt=0.;
630     double yopt=0.;
631     double zopt=0.;
632     double qual_dep=1.;
633     int nb_tet=noeud->get_lien_tetra()->get_nb();
634     double tab_coord[3000];
635     double gamma=0;
636    
637     for (int i=0;i<nb_tet;i++)
638     {
639     M3D_TETRA* tet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
640     qual_dep=std::min(qual_dep,OPERATEUR::qualite_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord()));
641    
642     // if(tet->get_id()==109676)
643     // {
644     // int ggg=0;
645     // }
646    
647     MG_NOEUD* no1=tet->get_noeud1();
648     MG_NOEUD* no2=tet->get_noeud2();
649     MG_NOEUD* no3=tet->get_noeud3();
650     MG_NOEUD* no4=tet->get_noeud4();
651     MG_NOEUD *noeud1,*noeud2,*noeud3;
652     if (no1==noeud)
653     {
654     noeud1=no2;
655     noeud2=no4;
656     noeud3=no3;
657     }
658     if (no2==noeud)
659     {
660     noeud1=no1;
661     noeud2=no3;
662     noeud3=no4;
663     }
664     if (no3==noeud)
665     {
666     noeud1=no1;
667     noeud2=no4;
668     noeud3=no2;
669     }
670     if (no4==noeud)
671     {
672     noeud1=no1;
673     noeud2=no2;
674     noeud3=no3;
675     }
676     double xyzi[3];
677     double *xyz1=noeud1->get_coord();
678     double *xyz2=noeud2->get_coord();
679     double *xyz3=noeud3->get_coord();
680     xyzi[0]=(xyz1[0]+xyz2[0]+xyz3[0])/3.;
681     xyzi[1]=(xyz1[1]+xyz2[1]+xyz3[1])/3.;
682     xyzi[2]=(xyz1[2]+xyz2[2]+xyz3[2])/3.;
683     OT_VECTEUR_3D v12(xyz1,xyz2);
684     OT_VECTEUR_3D v13(xyz1,xyz3);
685     OT_VECTEUR_3D v23(xyz2,xyz3);
686     OT_VECTEUR_3D normal=v12&v13;
687     double perimetre=v12.get_longueur()+v13.get_longueur()+v23.get_longueur();
688     double hauteur = (perimetre/3.) * 0.8 ;
689     normal.norme();
690    
691     tab_coord[3*i]=xyzi[0]+normal.get_x()*hauteur;
692     tab_coord[3*i+1]=xyzi[1]+normal.get_y()*hauteur;
693     tab_coord[3*i+2]=xyzi[2]+normal.get_z()*hauteur;
694    
695     if (OPERATEUR::qualite_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord())> 1e-10) gamma=gamma+1./OPERATEUR::qualite_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord())/OPERATEUR::qualite_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord());
696     }
697    
698     //-------------------------------------
699    
700     gamma=1./gamma;
701     for (int i=0;i<nb_tet;i++)
702     {
703     M3D_TETRA* tet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
704     if (OPERATEUR::qualite_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord())> 1e-10)
705     {
706     double alpha=gamma/OPERATEUR::qualite_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord())/OPERATEUR::qualite_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord());
707    
708     xopt=xopt+alpha*tab_coord[3*i];
709     yopt=yopt+alpha*tab_coord[3*i+1];
710     zopt=zopt+alpha*tab_coord[3*i+2];
711    
712    
713     }
714    
715     }
716    
717     //----------------Ghazal Added--------------
718    
719     double x_popt_ori=xopt-noeud->get_x();
720     double y_popt_ori=yopt-noeud->get_y();
721     double z_popt_ori=zopt-noeud->get_z();
722    
723     double dist=x_popt_ori*normal_plane(0)+y_popt_ori*normal_plane(1)+z_popt_ori*normal_plane(2);
724    
725     double xopt_planar=xopt-dist*normal_plane(0);
726     double yopt_planar=yopt-dist*normal_plane(1);
727     double zopt_planar=zopt-dist*normal_plane(2);
728    
729     xopt=xopt_planar;
730     yopt=yopt_planar;
731     zopt=zopt_planar;
732    
733     //-------------------------------------
734    
735     double delta[3]={xopt-noeud->get_x(),yopt-noeud->get_y(),zopt-noeud->get_z()};
736     double bmin=0.,bmax=1.;
737     double vieuxx=noeud->get_x();
738     double vieuxy=noeud->get_y();
739     double vieuxz=noeud->get_z();
740     double qualini=0.;
741    
742     for (int iteration=0;iteration<5;iteration++)
743     {
744     double alpha=(bmin+bmax)*0.5;
745     noeud->change_x(vieuxx+alpha*delta[0]);
746     noeud->change_y(vieuxy+alpha*delta[1]);
747     noeud->change_z(vieuxz+alpha*delta[2]);
748     double qualcour=1.;
749     for (int i=0;i<nb_tet;i++)
750     {
751     M3D_TETRA* mtet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
752     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());
753     qualcour=std::min(qualcour,qual);
754     }
755    
756    
757    
758     double alpha_eps=(bmin+bmax)*0.5-(bmax-bmin)/50.;
759     noeud->change_x(vieuxx+alpha_eps*delta[0]);
760     noeud->change_y(vieuxy+alpha_eps*delta[1]);
761     noeud->change_z(vieuxz+alpha_eps*delta[2]);
762     double qualcour_eps=1.;
763     for (int i=0;i<nb_tet;i++)
764     {
765     M3D_TETRA* mtet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
766     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());
767     qualcour_eps=std::min(qualcour_eps,qual);
768     }
769    
770    
771     if (qualcour>qualcour_eps) bmin =alpha;
772     else if(qualcour<qualcour_eps) bmax=alpha;
773     else
774     {
775     bmax=alpha/2;
776     }
777     qualini=std::max(qualini,qualcour);
778     }
779     //-------------------------------------
780    
781     noeud->change_x(vieuxx);
782     noeud->change_y(vieuxy);
783     noeud->change_z(vieuxz);
784    
785     if (qualini>qual_dep)
786     {
787     x=vieuxx+(bmin+bmax)*0.5*delta[0];
788     y=vieuxy+(bmin+bmax)*0.5*delta[1];
789     z=vieuxz+(bmin+bmax)*0.5*delta[2];
790     crit=qualini;
791     return 1;
792     }
793     //-------------------------------------
794    
795    
796    
797     }
798     else if( noeud->get_lien_topologie()->get_type()==MG_ELEMENT_TOPOLOGIQUE::TYPE_ELEMENT_TOPOLOGIQUE::ARETE)
799     {
800     // printf("On edge\n\n");
801    
802     double x_norm1=noeud->get_lien_topologie()->get_contrainte(3);
803     double y_norm1=noeud->get_lien_topologie()->get_contrainte(4);
804     double z_norm1=noeud->get_lien_topologie()->get_contrainte(5);
805    
806    
807     int nb_seg=noeud->get_lien_segment()->get_nb();
808    
809    
810     double dist_nodes;
811     double dist_nodes_max=0.0;
812     for (int i=0;i<nb_seg;i++)
813     {
814     MG_SEGMENT* segg=noeud->get_lien_segment()->get(i);
815    
816     dist_nodes=segg->get_longueur();
817    
818     dist_nodes_max=std::max(dist_nodes_max,dist_nodes);
819    
820     }
821    
822     double qual_dep=1.;
823     int nb_tet=noeud->get_lien_tetra()->get_nb();
824    
825    
826    
827     for (int i=0;i<nb_tet;i++)
828     {
829     M3D_TETRA* tet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
830     qual_dep=std::min(qual_dep,OPERATEUR::qualite_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord()));
831    
832     }
833     //-------------------------------------
834    
835    
836     double bmin=0.,bmax=1.;
837     double bmin2=0.,bmax2=1.;
838     double vieuxx=noeud->get_x();
839     double vieuxy=noeud->get_y();
840     double vieuxz=noeud->get_z();
841     double qualini1=0.;
842     double qualini2=0.;
843    
844    
845     for (int iteration=0;iteration<5;iteration++)
846     {
847     double alpha=(bmin+bmax)*0.5;
848     noeud->change_x(vieuxx+alpha*dist_nodes_max*x_norm1);
849     noeud->change_y(vieuxy+alpha*dist_nodes_max*y_norm1);
850     noeud->change_z(vieuxz+alpha*dist_nodes_max*z_norm1);
851     double qualcour=1.;
852     for (int i=0;i<nb_tet;i++)
853     {
854     M3D_TETRA* mtet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
855     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());
856     qualcour=std::min(qualcour,qual);
857     }
858    
859    
860    
861     double alpha_eps=(bmin+bmax)*0.5-(bmax-bmin)/50.;
862     noeud->change_x(vieuxx+alpha_eps*dist_nodes_max*x_norm1);
863     noeud->change_y(vieuxy+alpha_eps*dist_nodes_max*y_norm1);
864     noeud->change_z(vieuxz+alpha_eps*dist_nodes_max*z_norm1);
865     double qualcour_eps=1.;
866     for (int i=0;i<nb_tet;i++)
867     {
868     M3D_TETRA* mtet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
869     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());
870     qualcour_eps=std::min(qualcour_eps,qual);
871     }
872    
873    
874     if (qualcour>qualcour_eps) bmin =alpha;
875     else if(qualcour<qualcour_eps) bmax=alpha;
876     else
877     {
878     bmax=alpha/2;
879     }
880     qualini1=std::max(qualini1,qualcour);
881     }
882    
883     for (int iteration2=0;iteration2<5;iteration2++)
884     {
885     double alpha2=(bmin2+bmax2)*0.5;
886     noeud->change_x(vieuxx-alpha2*dist_nodes_max*x_norm1);
887     noeud->change_y(vieuxy-alpha2*dist_nodes_max*y_norm1);
888     noeud->change_z(vieuxz-alpha2*dist_nodes_max*z_norm1);
889     double qualcour=1.;
890     for (int i=0;i<nb_tet;i++)
891     {
892     M3D_TETRA* mtet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
893     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());
894     qualcour=std::min(qualcour,qual);
895     }
896    
897    
898    
899     double alpha_eps=(bmin2+bmax2)*0.5-(bmax2-bmin2)/50.;
900     noeud->change_x(vieuxx-alpha_eps*dist_nodes_max*x_norm1);
901     noeud->change_y(vieuxy-alpha_eps*dist_nodes_max*y_norm1);
902     noeud->change_z(vieuxz-alpha_eps*dist_nodes_max*z_norm1);
903     double qualcour_eps=1.;
904     for (int i=0;i<nb_tet;i++)
905     {
906     M3D_TETRA* mtet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
907     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());
908     qualcour_eps=std::min(qualcour_eps,qual);
909     }
910    
911    
912     if (qualcour>qualcour_eps) bmin2 =alpha2;
913     else if(qualcour<qualcour_eps) bmax2=alpha2;
914     else
915     {
916     bmax2=alpha2/2;
917     }
918     qualini2=std::max(qualini2,qualcour);
919     }
920     //-------------------------------------
921    
922     noeud->change_x(vieuxx);
923     noeud->change_y(vieuxy);
924     noeud->change_z(vieuxz);
925    
926    
927    
928     if ((qualini1>qual_dep)&&(qualini1>qualini2))
929     {
930     x=vieuxx+(bmin+bmax)*0.5*dist_nodes_max*x_norm1;
931     y=vieuxy+(bmin+bmax)*0.5*dist_nodes_max*y_norm1;
932     z=vieuxz+(bmin+bmax)*0.5*dist_nodes_max*z_norm1;
933     crit=qualini1;
934     return 1;
935     }
936     else if ((qualini2>qual_dep)&&(qualini2>qualini1))
937     {
938     x=vieuxx-(bmin2+bmax2)*0.5*dist_nodes_max*x_norm1;
939     y=vieuxy-(bmin2+bmax2)*0.5*dist_nodes_max*y_norm1;
940     z=vieuxz-(bmin2+bmax2)*0.5*dist_nodes_max*z_norm1;
941     crit=qualini2;
942     return 1;
943     }
944    
945     }
946     else if( noeud->get_lien_topologie()->get_type()==MG_ELEMENT_TOPOLOGIQUE::TYPE_ELEMENT_TOPOLOGIQUE::VOLUME)
947     {
948    
949     double xopt=0.;
950     double yopt=0.;
951     double zopt=0.;
952     double qual_dep=1.;
953     int nb_tet=noeud->get_lien_tetra()->get_nb();
954     double tab_coord[3000];
955     double gamma=0;
956    
957     for (int i=0;i<nb_tet;i++)
958     {
959     M3D_TETRA* tet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
960     qual_dep=std::min(qual_dep,OPERATEUR::qualite_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord()));
961    
962    
963     MG_NOEUD* no1=tet->get_noeud1();
964     MG_NOEUD* no2=tet->get_noeud2();
965     MG_NOEUD* no3=tet->get_noeud3();
966     MG_NOEUD* no4=tet->get_noeud4();
967     MG_NOEUD *noeud1,*noeud2,*noeud3;
968     if (no1==noeud)
969     {
970     noeud1=no2;
971     noeud2=no4;
972     noeud3=no3;
973     }
974     if (no2==noeud)
975     {
976     noeud1=no1;
977     noeud2=no3;
978     noeud3=no4;
979     }
980     if (no3==noeud)
981     {
982     noeud1=no1;
983     noeud2=no4;
984     noeud3=no2;
985     }
986     if (no4==noeud)
987     {
988     noeud1=no1;
989     noeud2=no2;
990     noeud3=no3;
991     }
992     double xyzi[3];
993     double *xyz1=noeud1->get_coord();
994     double *xyz2=noeud2->get_coord();
995     double *xyz3=noeud3->get_coord();
996     xyzi[0]=(xyz1[0]+xyz2[0]+xyz3[0])/3.;
997     xyzi[1]=(xyz1[1]+xyz2[1]+xyz3[1])/3.;
998     xyzi[2]=(xyz1[2]+xyz2[2]+xyz3[2])/3.;
999     OT_VECTEUR_3D v12(xyz1,xyz2);
1000     OT_VECTEUR_3D v13(xyz1,xyz3);
1001     OT_VECTEUR_3D v23(xyz2,xyz3);
1002     OT_VECTEUR_3D normal=v12&v13;
1003     double perimetre=v12.get_longueur()+v13.get_longueur()+v23.get_longueur();
1004     double hauteur = (perimetre/3.) * 0.8 ;
1005     normal.norme();
1006    
1007     tab_coord[3*i]=xyzi[0]+normal.get_x()*hauteur;
1008     tab_coord[3*i+1]=xyzi[1]+normal.get_y()*hauteur;
1009     tab_coord[3*i+2]=xyzi[2]+normal.get_z()*hauteur;
1010    
1011     if (OPERATEUR::qualite_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord())> 1e-10) gamma=gamma+1./OPERATEUR::qualite_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord())/OPERATEUR::qualite_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord());
1012     }
1013    
1014    
1015    
1016    
1017    
1018    
1019    
1020    
1021     gamma=1./gamma;
1022     for (int i=0;i<nb_tet;i++)
1023     {
1024     M3D_TETRA* tet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
1025     if (OPERATEUR::qualite_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord())> 1e-10)
1026     {
1027     double alpha=gamma/OPERATEUR::qualite_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord())/OPERATEUR::qualite_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord());
1028    
1029     xopt=xopt+alpha*tab_coord[3*i];
1030     yopt=yopt+alpha*tab_coord[3*i+1];
1031     zopt=zopt+alpha*tab_coord[3*i+2];
1032     }
1033    
1034     }
1035    
1036    
1037    
1038    
1039    
1040    
1041     double delta[3]={xopt-noeud->get_x(),yopt-noeud->get_y(),zopt-noeud->get_z()};
1042     double bmin=0.,bmax=1.;
1043     double vieuxx=noeud->get_x();
1044     double vieuxy=noeud->get_y();
1045     double vieuxz=noeud->get_z();
1046     double qualini=0.;
1047    
1048     for (int iteration=0;iteration<5;iteration++)
1049     {
1050     double alpha=(bmin+bmax)*0.5;
1051     noeud->change_x(vieuxx+alpha*delta[0]);
1052     noeud->change_y(vieuxy+alpha*delta[1]);
1053     noeud->change_z(vieuxz+alpha*delta[2]);
1054     double qualcour=1.;
1055     for (int i=0;i<nb_tet;i++)
1056     {
1057     M3D_TETRA* mtet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
1058     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());
1059     qualcour=std::min(qualcour,qual);
1060     }
1061     double alpha_eps=(bmin+bmax)*0.5-(bmax-bmin)/50.;
1062     noeud->change_x(vieuxx+alpha_eps*delta[0]);
1063     noeud->change_y(vieuxy+alpha_eps*delta[1]);
1064     noeud->change_z(vieuxz+alpha_eps*delta[2]);
1065     double qualcour_eps=1.;
1066     for (int i=0;i<nb_tet;i++)
1067     {
1068     M3D_TETRA* mtet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
1069     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());
1070     qualcour_eps=std::min(qualcour_eps,qual);
1071     }
1072     if (qualcour>qualcour_eps) bmin =alpha;
1073     else if(qualcour<qualcour_eps) bmax=alpha;
1074     else
1075     {
1076     bmax=alpha/2;
1077     }
1078     qualini=std::max(qualini,qualcour);
1079     }
1080    
1081    
1082    
1083    
1084    
1085     noeud->change_x(vieuxx);
1086     noeud->change_y(vieuxy);
1087     noeud->change_z(vieuxz);
1088    
1089    
1090    
1091    
1092    
1093    
1094     if (qualini>qual_dep)
1095     {
1096     x=vieuxx+(bmin+bmax)*0.5*delta[0];
1097     y=vieuxy+(bmin+bmax)*0.5*delta[1];
1098     z=vieuxz+(bmin+bmax)*0.5*delta[2];
1099     crit=qualini;
1100     return 1;
1101     }
1102    
1103    
1104     }
1105    
1106     //=================FIN=================================================
1107    
1108    
1109     return 0;
1110     }
1111     #endif
1112    
1113    
1114 francois 287 void MAILLEUR3D_OPTIMISATION::remaille_coquille(MG_NOEUD* noeud1,MG_NOEUD* noeud2, double& crit, COQUILLE& coque)
1115 francois 283 {
1116     //recherche du segment a supprimer
1117     MG_NOEUD* noeud;
1118     if (noeud1->get_id()>noeud2->get_id()) noeud=noeud2;
1119     else noeud=noeud1;
1120     int nb=noeud->get_lien_petit_segment()->get_nb();
1121     for (int i=0;i<nb;i++)
1122     {
1123     MG_SEGMENT* seg=noeud->get_lien_petit_segment()->get(i);
1124     if ( ((seg->get_noeud1()==noeud1) && (seg->get_noeud2()==noeud2)) || ((seg->get_noeud1()==noeud2) && (seg->get_noeud2()==noeud1)))
1125 francois 288 {
1126 francois 420 //if (mg_maillage->get_mg_geometrie()==NULL)
1127     //{
1128 francois 288 if (seg->get_dimension_topo_null()==2)
1129     {
1130     coque.taille=0;
1131     crit=0.;
1132     return;
1133     }
1134 francois 420 //}
1135     //else
1136 francois 791 if (seg->get_lien_topologie()!=NULL) if ((seg->get_lien_topologie()->get_dimension()!=3) || (seg->get_origine()!=MAGIC::ORIGINE::MAILLEUR_AUTO))
1137 francois 283 {
1138     coque.taille=0;
1139     crit=0.;
1140     return;
1141     }
1142 francois 288 }
1143 francois 283 }
1144     //recherche des tetra qui s enroule autour du segment
1145     int nb1=noeud1->get_lien_tetra()->get_nb();
1146     int nb2=noeud2->get_lien_tetra()->get_nb();
1147     MG_TETRA* coq[100];
1148     int nb_coq=0;
1149     for (int i=0;i<nb1;i++)
1150     for (int j=0;j<nb2;j++)
1151     if (noeud1->get_lien_tetra()->get(i)==noeud2->get_lien_tetra()->get(j))
1152     {
1153     coq[nb_coq]=noeud1->get_lien_tetra()->get(i);
1154     nb_coq++;
1155     }
1156     if ((nb_coq<4) || (nb_coq>10))
1157     {
1158     coque.taille=0;
1159     crit=0.;
1160     return;
1161     }
1162     //recherche des noeuds qui compose la coquille
1163     MG_NOEUD* tet_noeud[50];
1164     int nb_tet_noeud=0;
1165     coque.taille=nb_coq;
1166     for (int i=0;i<nb_coq;i++)
1167     {
1168     coque.tet[i]=coq[i];
1169 ghazal 1152
1170 francois 283 if ((coq[i]->get_noeud1()!=noeud1) && (coq[i]->get_noeud1()!=noeud2)) tet_noeud[nb_tet_noeud++]=coq[i]->get_noeud1();
1171     if ((coq[i]->get_noeud2()!=noeud1) && (coq[i]->get_noeud2()!=noeud2)) tet_noeud[nb_tet_noeud++]=coq[i]->get_noeud2();
1172     if ((coq[i]->get_noeud3()!=noeud1) && (coq[i]->get_noeud3()!=noeud2)) tet_noeud[nb_tet_noeud++]=coq[i]->get_noeud3();
1173     if ((coq[i]->get_noeud4()!=noeud1) && (coq[i]->get_noeud4()!=noeud2)) tet_noeud[nb_tet_noeud++]=coq[i]->get_noeud4();
1174     coque.volume=coque.volume+fabs(get_volume(coq[i]));
1175     }
1176     //ordonnancement des noeuds qui compose la coquille en polygone
1177     MG_NOEUD* polygone[20];
1178     int nb_poly=2;
1179     polygone[0]=tet_noeud[0];
1180     tet_noeud[0]=NULL;
1181     polygone[1]=tet_noeud[1];
1182     tet_noeud[1]=NULL;
1183     while (nb_poly!=nb_coq+1)
1184     {
1185     for (int j=0;j<nb_coq;j++)
1186     {
1187     if (tet_noeud[2*j]==polygone[nb_poly-1])
1188     {
1189     polygone[nb_poly++]=tet_noeud[2*j+1];
1190     tet_noeud[2*j]=NULL;
1191     tet_noeud[2*j+1]=NULL;
1192     }
1193     if (tet_noeud[2*j+1]==polygone[nb_poly-1])
1194     {
1195     polygone[nb_poly++]=tet_noeud[2*j];
1196     tet_noeud[2*j]=NULL;
1197     tet_noeud[2*j+1]=NULL;
1198     }
1199    
1200     }
1201     }
1202     //etude du positionnement du polygone
1203     double vol1=get_volume(polygone[0]->get_coord(),polygone[1]->get_coord(),polygone[nb_poly-1]->get_coord(),noeud1->get_coord());
1204     double vol2=get_volume(polygone[0]->get_coord(),polygone[nb_poly-1]->get_coord(),polygone[1]->get_coord(),noeud2->get_coord());
1205     if (vol1*vol2<0.)
1206     {
1207     crit=0.;
1208     coque.taille=0;
1209     return;
1210     }
1211     MG_NOEUD *noeuda=noeud1;
1212     MG_NOEUD *noeudb=noeud2;
1213     if (vol1<0.)
1214     {
1215     noeuda=noeud2;
1216     noeudb=noeud1;
1217     }
1218     //Examen des solutions
1219     int nb_solution,nb_triangle;
1220 ghazal 1152
1221 francois 283 if (nb_coq==4) {
1222     nb_solution=2;
1223     nb_triangle=4;
1224     }
1225     if (nb_coq==5) {
1226     nb_solution=5;
1227     nb_triangle=10;
1228     }
1229     if (nb_coq==6) {
1230     nb_solution=14;
1231     nb_triangle=20;
1232     }
1233     if (nb_coq==7) {
1234     nb_solution=42;
1235     nb_triangle=35;
1236     }
1237     if (nb_coq==8) {
1238     nb_solution=132;
1239     nb_triangle=56;
1240     }
1241     if (nb_coq==9) {
1242     nb_solution=429;
1243     nb_triangle=84;
1244     }
1245     if (nb_coq==10) {
1246     nb_solution=1430;
1247     nb_triangle=120;
1248     }
1249     double crit_triangle[120];
1250     for (int i=0;i<nb_triangle;i++)
1251     {
1252     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());
1253     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());
1254     crit_triangle[i]=std::min(crit1,crit2);
1255     }
1256     /* examen de chaque solution */
1257     double crit_opt=0.;
1258     int numero_solution= -1;
1259     double crit_solution[1430];
1260     for (int i=0;i<nb_solution;i++)
1261     {
1262     double volume=0.;
1263     for (int j=0;j<nb_coq-2;j++)
1264     {
1265     MG_NOEUD* noa=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][i][j]][0]];
1266     MG_NOEUD* nob=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][i][j]][1]];
1267     MG_NOEUD* noc=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][i][j]][2]];
1268     MG_NOEUD* nod=noeuda;
1269     volume=volume+fabs(get_volume(noa->get_coord(),nob->get_coord(),noc->get_coord(),nod->get_coord()));
1270     nod=noeudb;
1271     volume=volume+fabs(get_volume(noa->get_coord(),nob->get_coord(),noc->get_coord(),nod->get_coord()));
1272     }
1273     double eps=0.0018*pow(volume,0.666666666);
1274     if (OPERATEUR::egal(volume,coque.volume,eps))
1275     {
1276     crit_solution[i]=1.;
1277     for (int j=0;j<nb_coq-2;j++)
1278     crit_solution[i]=std::min(crit_solution[i],crit_triangle[tab_solution[nb_coq-4][i][j]]);
1279     }
1280     else crit_solution[i]=0.;
1281     if (crit_opt<crit_solution[i])
1282     {
1283     crit_opt=crit_solution[i];
1284     numero_solution=i;
1285     }
1286    
1287     }
1288     if (numero_solution==(-1))
1289     {
1290     crit=0.;
1291     coque.taille=0;
1292     return;
1293     }
1294     crit=crit_opt;
1295     for (int j=0;j<nb_coq-2;j++)
1296     {
1297     MG_NOEUD* noa=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][numero_solution][j]][0]];
1298     MG_NOEUD* nob=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][numero_solution][j]][1]];
1299     MG_NOEUD* noc=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][numero_solution][j]][2]];
1300     coque.new_tetra[8*j]=noa;
1301     coque.new_tetra[8*j+1]=nob;
1302     coque.new_tetra[8*j+2]=noc;
1303     coque.new_tetra[8*j+3]=noeuda;
1304     coque.new_tetra[8*j+4]=noa;
1305     coque.new_tetra[8*j+5]=noc;
1306     coque.new_tetra[8*j+6]=nob;
1307     coque.new_tetra[8*j+7]=noeudb;
1308     }
1309     }
1310    
1311    
1312 francois 287 void MAILLEUR3D_OPTIMISATION::ajouter_ordre_tetra(M3D_TETRA* tet,int num)
1313 francois 283 {
1314 francois 339 double val;
1315     if (num==0) val=tet->get_qualite();
1316     if (num==1) val=1./tet->get_qualite();
1317     std::pair<double,M3D_TETRA*> tmp(val,tet);
1318 francois 283 ORDRE_TETRA::iterator p=lst_tetra[num].insert(tmp);
1319     lst_tetraid[num][tet->get_id()]=p;
1320     }
1321    
1322 francois 287 void MAILLEUR3D_OPTIMISATION::supprimer_ordre_tetra(M3D_TETRA* tet)
1323 francois 283 {
1324     int num=0;
1325     ORDRE_TETRA_PARID::iterator it=lst_tetraid[num].find(tet->get_id());
1326     if (it==lst_tetraid[num].end())
1327     {
1328     num=1;
1329     it=lst_tetraid[num].find(tet->get_id());
1330     }
1331     if (it==lst_tetraid[num].end()) return;
1332     lst_tetra[num].erase(it->second);
1333     lst_tetraid[num].erase(it);
1334     }
1335    
1336    
1337 francois 287 void MAILLEUR3D_OPTIMISATION::change_niveau_optimisation(int num)
1338 francois 283 {
1339     niveau_optimisation=num;
1340     }
1341    
1342 francois 287 int MAILLEUR3D_OPTIMISATION::get_niveau_optimisation(void)
1343 francois 283 {
1344     return niveau_optimisation;
1345     }
1346    
1347 francois 287
1348    
1349     class MG_TRIANGLE* MAILLEUR3D_OPTIMISATION::insere_triangle(class MG_VOLUME* mgvol,MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,int origine)
1350     {
1351     MG_TRIANGLE* triangle=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud3->get_id(),noeud2->get_id());
1352     if (triangle!=NULL) return NULL;
1353     // insertion du triangle
1354     MG_SEGMENT* segment1=mg_maillage->get_mg_segment(noeud1->get_id(),noeud2->get_id());
1355     MG_SEGMENT* segment2=mg_maillage->get_mg_segment(noeud1->get_id(),noeud3->get_id());
1356     MG_SEGMENT* segment3=mg_maillage->get_mg_segment(noeud2->get_id(),noeud3->get_id());
1357     if (segment1==NULL) segment1=cree_segment(mgvol,noeud1,noeud2,origine);
1358     if (segment2==NULL) segment2=cree_segment(mgvol,noeud1,noeud3,origine);
1359     if (segment3==NULL) segment3=cree_segment(mgvol,noeud2,noeud3,origine);
1360     MG_TRIANGLE* tri=cree_triangle(mgvol,noeud1,noeud2,noeud3,segment1,segment2,segment3,origine);
1361     return tri;
1362     }
1363    
1364    
1365    
1366    
1367    
1368    
1369     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)
1370     {
1371     M3D_TRIANGLE* tri=new M3D_TRIANGLE(mgvol,noeud1,noeud2,noeud3,segment1,segment2,segment3,origine);
1372     mg_maillage->ajouter_mg_triangle(tri);
1373     return tri;
1374     }
1375    
1376    
1377    
1378     MG_SEGMENT* MAILLEUR3D_OPTIMISATION::cree_segment(class MG_VOLUME* mgvol,MG_NOEUD* noeud1,MG_NOEUD* noeud2,int origine)
1379     {
1380     MG_SEGMENT* seg=mg_maillage->ajouter_mg_segment(mgvol,noeud1,noeud2,origine);
1381 francois 288 seg->change_dimension_topo_null(3);
1382 francois 287 return seg;
1383     }