ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur_stl.cpp
Revision: 1106
Committed: Thu Nov 10 14:26:25 2022 UTC (2 years, 6 months ago) by francois
Original Path: magic/lib/mailleur_auto/src/mailleur_stl.cpp
File size: 26208 byte(s)
Log Message:
correction orthographique

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