ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur_stl.cpp
Revision: 1158
Committed: Thu Jun 13 22:18:49 2024 UTC (10 months, 4 weeks ago) by francois
Original Path: magic/lib/mailleur_auto/src/mailleur_stl.cpp
File size: 32739 byte(s)
Log Message:
compatibilité Ubuntu 22.04
Suppression des refeences à Windows
Ajout d'une banière

File Contents

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