ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur_stl.cpp
Revision: 1085
Committed: Fri Jan 7 19:52:56 2022 UTC (3 years, 4 months ago) by francois
Original Path: magic/lib/mailleur_auto/src/mailleur_stl.cpp
File size: 26178 byte(s)
Log Message:
mailleur_stl qui conserve la conformité topologique du maillage

File Contents

# User Rev Content
1 francois 1082 #include "mailleur_stl.h"
2     #include "mg_maillage.h"
3     #include "tpl_fonctions_generiques.h"
4     #include "fct_taille.h"
5     #include <tpl_map_entite.h>
6     #include <mg_maillage_outils.h>
7     #include <fct_taille_fem_solution_generateur_gradient.h>
8    
9    
10    
11     MAILLEUR_STL::MAILLEUR_STL(MG_MAILLAGE* maiori, MG_GESTIONNAIRE* gt, FCT_TAILLE* carte,double limite,double angle1,double angle2):MAILLEUR(false),maiin(maiori),gest(gt),metrique(carte),limite_discretisation(limite),angle_arete(angle1),angle_limite(angle2),pas(4),nx(25),ny(25),nz(25),qualmaxaopt(0.5),nbcoucheopt(1)
12     {
13     }
14    
15     MAILLEUR_STL::MAILLEUR_STL(MAILLEUR_STL& mdd):MAILLEUR(mdd),maiin(mdd.maiin),mai(mdd.mai),metrique(mdd.metrique),pas(mdd.pas),angle_arete(mdd.angle_arete),angle_limite(mdd.angle_limite),limite_discretisation(mdd.limite_discretisation),nx(mdd.nx),ny(mdd.ny),nz(mdd.nz),qualmaxaopt(mdd.qualmaxaopt),nbcoucheopt(mdd.nbcoucheopt)
16     {
17     }
18    
19    
20    
21     MAILLEUR_STL::~MAILLEUR_STL()
22     {
23     delete gesttmp;
24     }
25    
26    
27    
28    
29     int MAILLEUR_STL::maille(class MG_GROUPE_TOPOLOGIQUE* mggt)
30     {
31    
32     affiche((char*)" Parametre du mailleur");
33     char mess[255];
34     sprintf(mess,(char*)" Ratio limite de discrétisation = %f",limite_discretisation);
35     affiche(mess);
36     sprintf(mess,(char*)" Ratio limite angle arete = %f",angle_arete);
37     affiche(mess);
38     sprintf(mess,(char*)" Ratio limite angle de la cavité de Delaunay = %f",angle_limite);
39     affiche(mess);
40     affiche((char*)" Duplication du maillage original");
41     //mai=maiin->dupliquer(gest);
42     mai=maiin->detacher(gest);
43 francois 1083 affiche((char*)" Vérification de la conformité du maillage");
44     verification_conformite();
45 francois 1082 affiche((char*)" Adaptation de la carte de taille");
46     adapte_carte();
47 francois 1083 affiche((char*)" Determination des arêtes");
48     determine_arete();
49 francois 1082 affiche((char*)" Prepartion de la projection");
50     prepare_projection();
51     affiche((char*)" Maillage");
52     cree_maillage();
53     affiche((char*)" Optimisation de maillage");
54     int oldnb=mai->get_nb_mg_triangle();
55     for (int i=0;i<nbcoucheopt;i++)
56     {
57     sprintf(mess," Passe numero %d",i+1);
58     affiche(mess);
59     int nb=optimise_maillage();
60     if ((nb==1)||(nb==oldnb)) i=nbcoucheopt;
61     oldnb=nb;
62     }
63     return OK;
64     }
65    
66 francois 1083
67    
68     void MAILLEUR_STL::verification_conformite()
69     {
70     TPL_MAP_ENTITE<MG_TRIANGLE*> listtriverifie;
71     TPL_LISTE_ENTITE<MG_TRIANGLE*> listtriverifieseq;
72     MG_TRIANGLE* tri=mai->get_mg_triangle(0);
73     listtriverifie.ajouter(tri);
74     listtriverifieseq.ajouter(tri);
75     int i=0;
76     int nbcorrection=0;
77     while (listtriverifie.get_nb()!=mai->get_nb_mg_triangle())
78     {
79     MG_TRIANGLE* tri=listtriverifieseq.get(i);
80     MG_SEGMENT *lstseg[3];
81     lstseg[0]=tri->get_segment1();
82     lstseg[1]=tri->get_segment2();
83     lstseg[2]=tri->get_segment3();
84     for (int j=0;j<3;j++)
85     {
86     MG_SEGMENT* seg=lstseg[j];
87     MG_TRIANGLE *tri1=seg->get_lien_triangle()->get(0);
88     MG_TRIANGLE *tri2=seg->get_lien_triangle()->get(1);
89     if (tri1!=tri)
90     {
91     tri2=tri1;
92     tri1=tri;
93     }
94     MG_NOEUD* no1=seg->get_noeud1();
95     MG_NOEUD* no2=seg->get_noeud2();
96     int sens1=-1;
97     if ((no1==tri1->get_noeud1()) && (no2==tri1->get_noeud2())) sens1=1;
98     if ((no1==tri1->get_noeud2()) && (no2==tri1->get_noeud3())) sens1=1;
99     if ((no1==tri1->get_noeud3()) && (no2==tri1->get_noeud1())) sens1=1;
100     int sens2=-1;
101     if ((no1==tri2->get_noeud1()) && (no2==tri2->get_noeud2())) sens2=1;
102     if ((no1==tri2->get_noeud2()) && (no2==tri2->get_noeud3())) sens2=1;
103     if ((no1==tri2->get_noeud3()) && (no2==tri2->get_noeud1())) sens2=1;
104     if (sens1*sens2>0)
105     {
106     tri2->inverse_sens();
107     nbcorrection++;
108     }
109     int nblistavant=listtriverifie.get_nb();
110     listtriverifie.ajouter(tri2);
111     if (listtriverifie.get_nb()!=nblistavant)
112     listtriverifieseq.ajouter(tri2);
113     }
114     i++;
115     }
116     char mess[200];
117     if (nbcorrection!=0) sprintf(mess," Inversion de triangle %d",nbcorrection);
118     else sprintf(mess," Maillage conforme");
119     affiche(mess);
120    
121     }
122    
123    
124 francois 1082 int MAILLEUR_STL::optimise_maillage(void)
125     {
126     std::multimap<double,MG_TRIANGLE*> lsttriaopt;
127     LISTE_MG_TRIANGLE::iterator it;
128     for (MG_TRIANGLE* tri=mai->get_premier_triangle(it);tri!=NULL;tri=mai->get_suivant_triangle(it))
129     {
130     double val=OPERATEUR::qualite_triangle(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord(),tri->get_noeud3()->get_coord());
131     lsttriaopt.insert(std::pair<double,MG_TRIANGLE*>(val,tri));
132     tri->change_solution(val);
133     }
134     int ok=0;
135     int compteur_iter=0;
136     do
137     {
138     compteur_iter++;
139     if (lsttriaopt.begin()->first<qualmaxaopt)
140     {
141     MG_TRIANGLE *tri=lsttriaopt.begin()->second;
142     double crit1,crit2,crit3;
143     double x1,y1,z1;
144     double x2,y2,z2;
145     double x3,y3,z3;
146     int res1=bouge_point(tri->get_noeud1(),crit1,x1,y1,z1);
147     int res2=bouge_point(tri->get_noeud2(),crit2,x2,y2,z2);
148     int res3=bouge_point(tri->get_noeud3(),crit3,x3,y3,z3);
149     double critmax=0;
150     double nvx,nvy,nvz;
151     MG_NOEUD* nno;
152     if (res1==1)
153     if (crit1>critmax)
154     {critmax=crit1;nvx=x1;nvy=y1;nvz=z1;nno=tri->get_noeud1();}
155     if (res2==1)
156     if (crit2>critmax)
157     {critmax=crit2;nvx=x2;nvy=y2;nvz=z2;nno=tri->get_noeud2();}
158     if (res3==1)
159     if (crit3>critmax)
160     {critmax=crit3;nvx=x3;nvy=y3;nvz=z3;nno=tri->get_noeud3();}
161     if (critmax>lsttriaopt.begin()->first)
162     {
163     projete_point_sur_triangulation(nvx,nvy,nvz,nvx,nvy,nvz);
164     nno->change_x(nvx);
165     nno->change_y(nvy);
166     nno->change_z(nvz);
167     for (int i=0;i<nno->get_lien_triangle()->get_nb();i++)
168     {
169     MG_TRIANGLE* tri=nno->get_lien_triangle()->get(i);
170     double *xyz1=tri->get_noeud1()->get_coord();
171     double *xyz2=tri->get_noeud2()->get_coord();
172     double *xyz3=tri->get_noeud3()->get_coord();
173     double val=OPERATEUR::qualite_triangle(xyz1,xyz2,xyz3);
174     std::multimap<double,MG_TRIANGLE*>:: iterator it=lsttriaopt.lower_bound(tri->get_solution());
175     int trouve=0;
176     while (trouve==0)
177     if (it->first>tri->get_solution()) trouve=-1;
178     else if (it->second==tri) trouve=1;
179     else it++;
180     if (trouve==1) lsttriaopt.erase(it);
181     //lsttriaopt.insert(std::pair<double,MG_TRIANGLE*>(val,tri));
182     //tri->change_solution(val);
183     }
184    
185     }
186     else lsttriaopt.erase(lsttriaopt.begin());
187    
188    
189     }
190     else
191     ok=1;
192     }
193     while (ok==0);
194     char chaine[500];
195     sprintf(chaine," %d iterations effectuées",compteur_iter);
196     affiche(chaine);
197     return(compteur_iter);
198     }
199    
200    
201     void MAILLEUR_STL::adapte_carte(void)
202     {
203     LISTE_MG_SEGMENT::iterator it;
204     std::vector<double> lstraf;
205     for (MG_SEGMENT* seg=mai->get_premier_segment(it);seg!=NULL;seg=mai->get_suivant_segment(it))
206     {
207     double val=metrique->calcul_distance_metrique(seg,pas);
208     if (val<1.)
209     {
210     double l=seg->get_longueur()/limite_discretisation/1.01;
211     double *xyz=seg->get_noeud1()->get_coord();
212     lstraf.push_back(xyz[0]);
213     lstraf.push_back(xyz[1]);
214     lstraf.push_back(xyz[2]);
215     lstraf.push_back(l);
216     xyz=seg->get_noeud2()->get_coord();
217     lstraf.push_back(xyz[0]);
218     lstraf.push_back(xyz[1]);
219     lstraf.push_back(xyz[2]);
220     lstraf.push_back(l);
221     }
222     }
223     gesttmp=new MG_GESTIONNAIRE();
224     FCT_TAILLE_FEM_SOLUTION_GENERATEUR_GRADIENT* metriquetmp=new FCT_TAILLE_FEM_SOLUTION_GENERATEUR_GRADIENT(gesttmp,mai,metrique,2,nx,ny,nz,(char*)"carteadapte.sol",lstraf);
225     metriquetmp->construit();
226     metrique2=metrique;
227     metrique=metriquetmp;
228     }
229    
230    
231    
232     void MAILLEUR_STL::prepare_projection()
233     {
234     LISTE_MG_NOEUD::iterator it;
235     MG_NOEUD* no=maiin->get_premier_noeud(it);
236     BOITE_3D boite(no->get_boite_3D());
237     TPL_LISTE_ENTITE<MG_NOEUD*> lstini;
238     for (MG_NOEUD* no=maiin->get_premier_noeud(it);no!=NULL;no=maiin->get_suivant_noeud(it))
239     {
240     boite=no->get_boite_3D()+boite;
241     lstini.ajouter(no);
242     }
243     boite.change_grosseur(1.2);
244     octree.initialiser(&lstini,1,boite.get_xmin(),boite.get_ymin(),boite.get_zmin(),boite.get_xmax(),boite.get_ymax(),boite.get_zmax());
245     LISTE_MG_TRIANGLE::iterator itt;
246     for (MG_TRIANGLE* tri=maiin->get_premier_triangle(itt);tri!=NULL;tri=maiin->get_suivant_triangle(itt))
247     octree.inserer(tri);
248     }
249    
250    
251     void MAILLEUR_STL::determine_arete(void)
252     {
253     LISTE_MG_SEGMENT::iterator it;
254     int nbsegmentarete=0;
255     for (MG_SEGMENT* seg=mai->get_premier_segment(it);seg!=NULL;seg=mai->get_suivant_segment(it))
256     {
257     int nbtri=seg->get_lien_triangle()->get_nb();
258     if (nbtri!=2)
259     {
260     char chaine[500];
261     sprintf(chaine," Maillage STL invalide : segment à %d voisin : noeud 1 %lu noeud 2 %lu ",nbtri,seg->get_noeud1()->get_id(),seg->get_noeud2()->get_id());
262     affiche(chaine);
263     }
264     if (nbtri==2)
265     {
266     MG_TRIANGLE *tri1=seg->get_lien_triangle()->get(0);
267     MG_TRIANGLE *tri2=seg->get_lien_triangle()->get(1);
268     double angle=get_angle_par_noeud<MG_NOEUD*>(tri1->get_noeud1(),tri1->get_noeud2(),tri1->get_noeud3(),tri2->get_noeud1(),tri2->get_noeud2(),tri2->get_noeud3());
269     if ((angle<M_PI-angle_arete) || (angle > M_PI+angle_arete))
270     {
271     seg->change_origine(MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE);
272     seg->get_noeud1()->change_origine(MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE);
273     seg->get_noeud2()->change_origine(MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE);
274     nbsegmentarete++;
275     }
276     ajouter_seg_a_discretiser(seg);
277     }
278     }
279     char mess[255];
280     sprintf(mess,(char*)" Nombre de segments aretes = %d",nbsegmentarete);
281     affiche(mess);
282     }
283    
284    
285 francois 1083 void MAILLEUR_STL::valide_maillage(void)
286     {
287     LISTE_MG_SEGMENT::iterator it;
288     for (MG_SEGMENT* seg=mai->get_premier_segment(it);seg!=NULL;seg=mai->get_suivant_segment(it))
289     {
290     int nbtri=seg->get_lien_triangle()->get_nb();
291     if (nbtri!=2)
292     affiche((char*)" **** erreur **** un triangle n'a pas deux voisins");
293     MG_TRIANGLE *tri1=seg->get_lien_triangle()->get(0);
294     MG_TRIANGLE *tri2=seg->get_lien_triangle()->get(1);
295     MG_NOEUD* no1=seg->get_noeud1();
296     MG_NOEUD* no2=seg->get_noeud2();
297     int sens1=-1;
298     if ((no1==tri1->get_noeud1()) && (no2==tri1->get_noeud2())) sens1=1;
299     if ((no1==tri1->get_noeud2()) && (no2==tri1->get_noeud3())) sens1=1;
300     if ((no1==tri1->get_noeud3()) && (no2==tri1->get_noeud1())) sens1=1;
301     int sens2=-1;
302     if ((no1==tri2->get_noeud1()) && (no2==tri2->get_noeud2())) sens2=1;
303     if ((no1==tri2->get_noeud2()) && (no2==tri2->get_noeud3())) sens2=1;
304     if ((no1==tri2->get_noeud3()) && (no2==tri2->get_noeud1())) sens2=1;
305     if (sens1*sens2>0)
306     affiche((char*)" **** erreur **** un segment n'est pas pris dans chaque sens pour deux mailles vosines");
307     double xyz1[3];
308     double xyz2[3];
309     tri1->normal(xyz1);
310     tri2->normal(xyz2);
311     OT_VECTEUR_3D vec1(xyz1);
312     OT_VECTEUR_3D vec2(xyz2);
313     if (vec1*vec2<0)
314     if (seg->get_origine()!=MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE)
315     affiche((char*)" **** erreur **** deux triangles voisins n'ont pas la normal dans le meme sens");
316    
317    
318     }
319     }
320    
321    
322    
323 francois 1082 void MAILLEUR_STL::calcul_nouveau_noeud_et_etoile(MG_SEGMENT* seg,double& x, double& y, double& z,OT_VECTEUR_3D &normale, TPL_LISTE_ENTITE<MG_TRIANGLE*> &listuniquetri,std::map<unsigned long,std::pair<MG_SEGMENT*,int> > &lstsegment)
324     {
325     MG_TRIANGLE *tri1=seg->get_lien_triangle()->get(0);
326     MG_TRIANGLE *tri2=seg->get_lien_triangle()->get(1);
327     double xyz1[3],xyz2[3];
328     tri1->normal(xyz1);
329     tri2->normal(xyz2);
330     OT_VECTEUR_3D nor1(xyz1);
331     OT_VECTEUR_3D nor2(xyz2);
332     normale=0.5*nor1+0.5*nor2;
333     MG_NOEUD* no1=seg->get_noeud1();
334     MG_NOEUD* no2=seg->get_noeud2();
335     x=0.5*(no1->get_x()+no2->get_x());
336     y=0.5*(no1->get_y()+no2->get_y());
337     z=0.5*(no1->get_z()+no2->get_z());
338     listuniquetri.ajouter(tri1);
339     listuniquetri.ajouter(tri2);
340     for (int i=0;i<listuniquetri.get_nb();i++)
341     {
342     MG_TRIANGLE* tri=listuniquetri.get(i);
343     for (int cote=1;cote<4;cote++)
344     {
345     MG_TRIANGLE* triv=get_voisin(tri,cote);
346     if (!(listuniquetri.est_dans_la_liste(triv)))
347     if (triv!=NULL)
348     {
349     double xyzv[3];
350     triv->normal(xyzv);
351     OT_VECTEUR_3D norv(xyzv);
352     if (norv*normale>cos(angle_limite))
353     {
354     bool res=respecte_delaunay(x,y,z,triv);
355     if (!res)
356     listuniquetri.ajouter(triv);
357     }
358     }
359     }
360    
361    
362    
363     }
364     for (int i=0;i<listuniquetri.get_nb();i++)
365     {
366     MG_TRIANGLE* tri=listuniquetri.get(i);
367     std::map<unsigned long,std::pair<MG_SEGMENT*,int> >::iterator ptmp;
368     MG_SEGMENT* seg=tri->get_segment1();
369     ptmp=lstsegment.find(seg->get_id());
370     if (ptmp==lstsegment.end()) lstsegment[seg->get_id()]=std::pair<MG_SEGMENT*,int>(seg,1);
371     else ptmp->second.second=ptmp->second.second+1;
372     seg=tri->get_segment2();
373     ptmp=lstsegment.find(seg->get_id());
374     if (ptmp==lstsegment.end()) lstsegment[seg->get_id()]=std::pair<MG_SEGMENT*,int>(seg,1);
375     else ptmp->second.second=ptmp->second.second+1;
376     seg=tri->get_segment3();
377     ptmp=lstsegment.find(seg->get_id());
378     if (ptmp==lstsegment.end()) lstsegment[seg->get_id()]=std::pair<MG_SEGMENT*,int>(seg,1);
379     else ptmp->second.second=ptmp->second.second+1;
380     }
381    
382    
383    
384     projete_point_sur_triangulation(x,y,z,x,y,z);
385    
386    
387     bool convergecoquille;
388    
389     do
390     {
391     convergecoquille=true;
392     TPL_LISTE_ENTITE<MG_SEGMENT*> contour;
393     for (std::map<unsigned long,std::pair<MG_SEGMENT*,int> >::iterator ptmp=lstsegment.begin();ptmp!=lstsegment.end();ptmp++)
394     if (ptmp->second.second==1)
395     {
396     MG_SEGMENT* seg=ptmp->second.first;
397     contour.ajouter(seg);
398     }
399    
400     for (int i=0;i<contour.get_nb();i++)
401     {
402     MG_SEGMENT* seg=contour.get(i);
403     MG_TRIANGLE *triin,*triext;
404     if (listuniquetri.est_dans_la_liste(seg->get_lien_triangle()->get(0)))
405     {
406     triin=seg->get_lien_triangle()->get(0);
407     triext=seg->get_lien_triangle()->get(1);
408     }
409     else
410     {
411     triin=seg->get_lien_triangle()->get(1);
412     triext=seg->get_lien_triangle()->get(0);
413     }
414     MG_NOEUD* no1=seg->get_noeud1();
415     MG_NOEUD* no2=seg->get_noeud2();
416     int sens=-1;
417     if ((triin->get_noeud1()==no1) && (triin->get_noeud2()==no2)) sens=1;
418     if ((triin->get_noeud2()==no1) && (triin->get_noeud3()==no2)) sens=1;
419     if ((triin->get_noeud3()==no1) && (triin->get_noeud1()==no2)) sens=1;
420     double *xyz1,*xyz2;
421     if (sens==-1)
422     {
423     xyz1=no2->get_coord();
424     xyz2=no1->get_coord();
425     }
426     else
427     {
428     xyz1=no1->get_coord();
429     xyz2=no2->get_coord();
430     }
431     double xyz3[3]={x,y,z};
432     OT_VECTEUR_3D vec1(xyz1,xyz2);
433     OT_VECTEUR_3D vec2(xyz1,xyz3);
434     OT_VECTEUR_3D nor=vec1&vec2;
435     nor.norme();
436     double xyzin[3];
437     triin->normal(xyzin);
438     OT_VECTEUR_3D norin(xyzin);
439     if (nor*norin<0.)
440     {
441 francois 1083 affiche((char*)"etoile non conforme");
442 francois 1082 listuniquetri.ajouter(triext);
443     std::map<unsigned long,std::pair<MG_SEGMENT*,int> >::iterator ptmp;
444     MG_SEGMENT* seg=triext->get_segment1();
445     ptmp=lstsegment.find(seg->get_id());
446     ptmp=lstsegment.find(seg->get_id());
447     if (ptmp==lstsegment.end()) lstsegment[seg->get_id()]=std::pair<MG_SEGMENT*,int>(seg,1);
448     else ptmp->second.second=ptmp->second.second+1;
449     seg=triext->get_segment2();
450     ptmp=lstsegment.find(seg->get_id());
451     if (ptmp==lstsegment.end()) lstsegment[seg->get_id()]=std::pair<MG_SEGMENT*,int>(seg,1);
452     else ptmp->second.second=ptmp->second.second+1;
453     seg=triext->get_segment3();
454     ptmp=lstsegment.find(seg->get_id());
455     if (ptmp==lstsegment.end()) lstsegment[seg->get_id()]=std::pair<MG_SEGMENT*,int>(seg,1);
456     else ptmp->second.second=ptmp->second.second+1;
457     convergecoquille=false;
458     }
459     }
460     }
461     while (convergecoquille==false);
462     }
463    
464     int MAILLEUR_STL::bouge_point(MG_NOEUD* mg_noeud,double& crit,double& x,double& y, double& z)
465     {
466     if (mg_noeud->get_origine()==MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE)
467     return 0;
468     TPL_MAP_ENTITE<MG_NOEUD*> lstnoeud;
469     double critdep=1.;
470     for (int i=0;i<mg_noeud->get_lien_triangle()->get_nb();i++)
471     {
472     MG_TRIANGLE* tri=mg_noeud->get_lien_triangle()->get(i);
473     double val=OPERATEUR::qualite_triangle(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord(),tri->get_noeud3()->get_coord());
474     if (val<critdep) critdep=val;
475     if (tri->get_noeud1()!=mg_noeud) lstnoeud.ajouter(tri->get_noeud1());
476     if (tri->get_noeud2()!=mg_noeud) lstnoeud.ajouter(tri->get_noeud2());
477     if (tri->get_noeud3()!=mg_noeud) lstnoeud.ajouter(tri->get_noeud3());
478     }
479     TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR it;
480     double xi=0.;
481     double yi=0.;
482     double zi=0.;
483     for (MG_NOEUD* no=lstnoeud.get_premier(it);no!=NULL;no=lstnoeud.get_suivant(it))
484     {
485     xi=xi+no->get_x();
486     yi=yi+no->get_y();
487     zi=zi+no->get_z();
488     }
489     xi=xi/lstnoeud.get_nb();
490     yi=yi/lstnoeud.get_nb();
491     zi=zi/lstnoeud.get_nb();
492     OT_VECTEUR_3D dir(xi-mg_noeud->get_x(),yi-mg_noeud->get_y(),zi-mg_noeud->get_z());
493     if (dir.get_longueur()<1e-10) return 0;
494     double alpha=0.1;
495     bool converge=false;
496     double critencours=critdep;
497     int passe=0;
498     while (converge==false)
499     {
500     double xyz[3];
501     xyz[0]=mg_noeud->get_x()+alpha*dir.get_x();
502     xyz[1]=mg_noeud->get_y()+alpha*dir.get_y();
503     xyz[2]=mg_noeud->get_z()+alpha*dir.get_z();
504     double critcourant=1.;
505     for (int i=0;i<mg_noeud->get_lien_triangle()->get_nb();i++)
506     {
507     MG_TRIANGLE* tri=mg_noeud->get_lien_triangle()->get(i);
508     double *xyz1=tri->get_noeud1()->get_coord();
509     double *xyz2=tri->get_noeud2()->get_coord();
510     double *xyz3=tri->get_noeud3()->get_coord();
511     if (tri->get_noeud1()==mg_noeud) xyz1=xyz;
512     if (tri->get_noeud2()==mg_noeud) xyz2=xyz;
513     if (tri->get_noeud3()==mg_noeud) xyz3=xyz;
514     double val=OPERATEUR::qualite_triangle(xyz1,xyz2,xyz3);
515     if (val<critcourant) critcourant=val;
516     }
517     if (critencours>critcourant)
518     {
519     converge=true;
520     return passe;
521     }
522     else
523     {
524     alpha=alpha+0.1;
525     x=xyz[0];
526     y=xyz[1];
527     z=xyz[2];
528     critencours=critcourant;
529     crit=critencours;
530     passe=1;
531     }
532    
533     }
534     return 0;
535     }
536    
537     void MAILLEUR_STL::maillage_et_segadiscretiser_update(MG_SEGMENT* segbase,double x, double y, double z, OT_VECTEUR_3D normale,TPL_LISTE_ENTITE<MG_TRIANGLE *>& listuniquetri, std::map<unsigned long, std::pair<MG_SEGMENT *, int> >& lstsegment)
538     {
539 francois 1085 //projete_point_sur_triangulation(x,y,z,x,y,z);
540 francois 1082 MG_NOEUD *no=mai->ajouter_mg_noeud(NULL,x,y,z,segbase->get_origine());
541     MG_NOEUD* nobase1=NULL;
542     MG_NOEUD* nobase2=NULL;
543     if (segbase->get_origine()==MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE)
544     {
545     nobase1=segbase->get_noeud1();
546     nobase2=segbase->get_noeud2();
547     }
548     for (std::map<unsigned long,std::pair<MG_SEGMENT*,int> >::iterator ptmp=lstsegment.begin();ptmp!=lstsegment.end();ptmp++)
549     if (ptmp->second.second==2)
550     supprimer_seg_a_discretiser(ptmp->second.first);
551     for (int i=0;i<listuniquetri.get_nb();i++)
552     {
553     MG_TRIANGLE* tri=listuniquetri.get(i);
554     mai->supprimer_mg_triangleid(tri->get_id());
555     }
556     for (std::map<unsigned long,std::pair<MG_SEGMENT*,int> >::iterator ptmp=lstsegment.begin();ptmp!=lstsegment.end();ptmp++)
557     if (ptmp->second.second==1)
558     {
559     MG_SEGMENT* seg=ptmp->second.first;
560 francois 1085 MG_TRIANGLE* triext=seg->get_lien_triangle()->get(0);
561     int sens=1;
562     if ((triext->get_noeud1()==seg->get_noeud1()) && (triext->get_noeud2()==seg->get_noeud2()) ) sens=-1;
563     if ((triext->get_noeud2()==seg->get_noeud1()) && (triext->get_noeud3()==seg->get_noeud2()) ) sens=-1;
564     if ((triext->get_noeud3()==seg->get_noeud1()) && (triext->get_noeud1()==seg->get_noeud2()) ) sens=-1;
565     MG_NOEUD *no1,*no2;
566     if (sens==1)
567     {
568     no1=seg->get_noeud1();
569     no2=seg->get_noeud2();
570     }
571     else
572     {
573     no2=seg->get_noeud1();
574     no1=seg->get_noeud2();
575     }
576    
577    
578     MG_SEGMENT* mgsegment1=mai->get_mg_segment(no1->get_id(),no2->get_id());
579     MG_SEGMENT* mgsegment2=mai->get_mg_segment(no2->get_id(),no->get_id());
580     MG_SEGMENT* mgsegment3=mai->get_mg_segment(no->get_id(),no1->get_id());
581 francois 1082 if (mgsegment1==NULL)
582 francois 1083 {//mormalement ne devrai jamais passer ici
583 francois 1085 mgsegment1=mai->ajouter_mg_segment(NULL,no1,no2,MAGIC::ORIGINE::TRIANGULATION);
584 francois 1082 ajouter_seg_a_discretiser(mgsegment1);
585     }
586     if (mgsegment2==NULL)
587     {
588     int origine=MAGIC::ORIGINE::TRIANGULATION;
589 francois 1085 if ((no2==nobase1)|| (no2==nobase2))
590 francois 1082 origine=MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE;
591 francois 1085 mgsegment2=mai->ajouter_mg_segment(NULL,no2,no,origine);
592 francois 1082 ajouter_seg_a_discretiser(mgsegment2);
593     }
594     if (mgsegment3==NULL)
595     {
596     int origine=MAGIC::ORIGINE::TRIANGULATION;
597 francois 1085 if ((no1==nobase1)|| (no1==nobase2))
598 francois 1082 origine=MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE;
599 francois 1085 mgsegment3=mai->ajouter_mg_segment(NULL,no,no1,origine);
600 francois 1082 ajouter_seg_a_discretiser(mgsegment3);
601     }
602 francois 1085 MG_TRIANGLE* nvtri=new MG_TRIANGLE(NULL,no1,no2,no,mgsegment1,mgsegment2,mgsegment3,MAGIC::ORIGINE::TRIANGULATION);
603 francois 1082 mai->ajouter_mg_triangle(nvtri);
604 francois 1085 /*double xyz[3];
605 francois 1082 nvtri->normal(xyz);
606     OT_VECTEUR_3D nvnor(xyz);
607     if (nvnor*normale<0)
608 francois 1085 nvtri->inverse_sens();*/
609 francois 1082
610 francois 1083
611 francois 1082 }
612     }
613    
614    
615     void MAILLEUR_STL::cree_maillage()
616     {
617     int ok=0;
618     int compteur_iter=0;
619 francois 1083 //valide_maillage();
620    
621 francois 1082 do
622     {
623     compteur_iter++;
624     if (lstsegadiscretiser.begin()->first>limite_discretisation)
625     {
626     MG_SEGMENT* seg=lstsegadiscretiser.begin()->second;
627     double x,y,z;
628     OT_VECTEUR_3D normale;
629     TPL_LISTE_ENTITE<MG_TRIANGLE*> listtrietoile;
630     std::map<unsigned long,std::pair<MG_SEGMENT*,int> > lstsegmentetoile;
631     calcul_nouveau_noeud_et_etoile(seg,x,y,z,normale,listtrietoile,lstsegmentetoile);
632     maillage_et_segadiscretiser_update(seg,x,y,z,normale,listtrietoile,lstsegmentetoile);
633     }
634     else
635     ok=1;
636     }
637     while (ok==0);
638     char chaine[500];
639     sprintf(chaine," %d iterations effectuées",compteur_iter);
640     affiche(chaine);
641     }
642    
643    
644    
645    
646     bool MAILLEUR_STL::respecte_delaunay(double x, double y, double z, MG_TRIANGLE* tri)
647     {
648     OT_VECTEUR_3D vec1(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord());
649     OT_VECTEUR_3D vec2(tri->get_noeud1()->get_coord(),tri->get_noeud3()->get_coord());
650     OT_VECTEUR_3D vec3=vec1&vec2;
651     vec1.norme();
652     vec3.norme();
653     vec2=vec3&vec1;
654     OT_MATRICE_3D base(vec1,vec2,vec3);
655     OT_MATRICE_3D baset=base.transpose();;
656     OT_VECTEUR_3D pt1(tri->get_noeud1()->get_coord());
657     OT_VECTEUR_3D pt2(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord());
658     OT_VECTEUR_3D pt3(tri->get_noeud1()->get_coord(),tri->get_noeud3()->get_coord());
659     OT_VECTEUR_3D pt1l(0.,0.,0.);
660     OT_VECTEUR_3D pt2l=baset*pt2;
661     OT_VECTEUR_3D pt3l=baset*pt3;
662    
663     OT_VECTEUR_3D colx(pt1l.get_x(),pt2l.get_x(),pt3l.get_x());
664     OT_VECTEUR_3D coly(pt1l.get_y(),pt2l.get_y(),pt3l.get_y());
665     OT_VECTEUR_3D colz(1.,1.,1.);
666     OT_VECTEUR_3D colx2y2(pt1l.get_x()*pt1l.get_x()+pt1l.get_y()*pt1l.get_y(),pt2l.get_x()*pt2l.get_x()+pt2l.get_y()*pt2l.get_y(),pt3l.get_x()*pt3l.get_x()+pt3l.get_y()*pt3l.get_y());
667     OT_MATRICE_3D sys1(colx,coly,colz);
668     double det=2.*sys1.get_determinant();
669     OT_MATRICE_3D sys2(colx2y2,coly,colz);
670     double xcl=sys2.get_determinant()/det;
671     OT_MATRICE_3D sys3(colx2y2,colx,colz);
672     double ycl=-sys3.get_determinant()/det;
673     OT_VECTEUR_3D centrel(xcl,ycl,0.);
674     OT_VECTEUR_3D centre=pt1+base*centrel;
675    
676     OT_VECTEUR_3D r1(centre.get_xyz(),tri->get_noeud1()->get_coord());
677     double rayon=r1.get_longueur();
678     OT_VECTEUR_3D dist(x,y,z);
679     dist=dist-centre;
680     bool rep=true;
681     if (dist.get_longueur()<rayon)
682     rep=false;
683     return rep;
684     }
685    
686    
687    
688     void MAILLEUR_STL::projete_point_sur_triangulation(double x, double y, double z, double& nvx, double& nvy, double& nvz)
689     {
690     double mat[9];
691     double xyz[3]={x,y,z};
692     metrique->evaluer(xyz,mat);
693     double dref=1./sqrt(mat[0]);
694     TPL_MAP_ENTITE<MG_TRIANGLE*> lst;
695     while (lst.get_nb()==0)
696     {
697     octree.rechercher(x,y,z,dref,lst);
698     dref=dref*1.05;
699     }
700     MG_MAILLAGE_OUTILS ot;
701     TPL_MAP_ENTITE<MG_TRIANGLE*>::ITERATEUR it;
702     double dis=1e308;
703     for (MG_TRIANGLE* tri=lst.get_premier(it);tri!=NULL;tri=lst.get_suivant(it))
704     {
705     double xt,yt,zt;
706     ot.projetedansletriangle(tri,x,y,z,xt,yt,zt);
707     OT_VECTEUR_3D vec(x-xt,y-yt,z-zt);
708     double l=vec.get_longueur();
709     if (l<dis)
710     {
711     dis=l;
712     nvx=xt;
713     nvy=yt;
714     nvz=zt;
715     }
716     }
717     }
718    
719    
720    
721     void MAILLEUR_STL::change_angle_arete(double val)
722     {
723     angle_arete=val;
724     }
725    
726    
727     double MAILLEUR_STL::get_angle_arete()
728     {
729     return angle_arete;
730     }
731    
732     MG_TRIANGLE * MAILLEUR_STL::get_voisin(MG_TRIANGLE* tri, int cote)
733     {
734     MG_SEGMENT* seg;
735     if (cote==1) seg=tri->get_segment1();
736     if (cote==2) seg=tri->get_segment2();
737     if (cote==3) seg=tri->get_segment3();
738     if (seg->get_origine()==MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE) return NULL;
739     if (seg->get_lien_triangle()->get_nb()!=2) return NULL;
740     MG_TRIANGLE* tri1=seg->get_lien_triangle()->get(0);
741     MG_TRIANGLE* tri2=seg->get_lien_triangle()->get(1);
742     if (tri1==tri) return tri2;
743     return tri1;
744     }
745    
746     void MAILLEUR_STL::ajouter_seg_a_discretiser(MG_SEGMENT* seg)
747     {
748     double l=metrique->calcul_distance_metrique(seg,pas);
749     seg->change_distance_metrique(l);
750     lstsegadiscretiser.insert(std::pair<double,MG_SEGMENT*>(l,seg));
751     }
752    
753     void MAILLEUR_STL::supprimer_seg_a_discretiser(MG_SEGMENT* seg)
754     {
755     double l;
756     seg->get_distance_metrique(&l);
757     std::multimap<double,MG_SEGMENT*,std::greater<double > >:: iterator it=lstsegadiscretiser.lower_bound(l);
758     while (it->second!=seg)
759     it++;
760     lstsegadiscretiser.erase(it);
761     }
762    
763     void MAILLEUR_STL::change_pas(int val)
764     {
765     pas=val;
766     }
767    
768     void MAILLEUR_STL::change_nbpasseoptimisation(int nb)
769     {
770     nbcoucheopt=nb;
771     }
772    
773    
774     void MAILLEUR_STL::change_nx(int val)
775     {
776     nx=val;
777     }
778    
779    
780     void MAILLEUR_STL::change_ny(int val)
781     {
782     ny=val;
783     }
784    
785    
786     void MAILLEUR_STL::change_nz(int val)
787     {
788     nz=val;
789     }
790    
791     void MAILLEUR_STL::change_nxnynz(int nvx, int nvy, int nvz)
792     {
793     nx=nvx;
794     ny=nvy;
795     nz=nvz;
796     }
797    
798    
799     void MAILLEUR_STL::change_qualmaxaoptimiser(double val)
800     {
801     qualmaxaopt=val;
802     }