ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur_stl.cpp
Revision: 1111
Committed: Tue Apr 4 13:48:28 2023 UTC (2 years, 1 month ago) by francois
Original Path: magic/lib/mailleur_auto/src/mailleur_stl.cpp
File size: 31713 byte(s)
Log Message:
mailleur STL avec surface ouverte + liste de noeud contraint 

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