ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur_stl.cpp
Revision: 1189
Committed: Tue Feb 4 17:26:49 2025 UTC (3 months ago) by francois
File size: 32739 byte(s)
Log Message:
Version 5.0 de MAGIC. Integration de ALGLIB pour faire de l'optimisation. ALGLIB se download automatiquement en executant un script dans le repertoire config update_magic.bash


File Contents

# Content
1 //####//------------------------------------------------------------
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 #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 #include "ot_boite_3d.h"
30 #include "ot_mathematique.h"
31
32
33
34 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 {
36 fichierpoint[0]=0;
37 if (nomfichierpoint!=NULL) (fichierpoint,nomfichierpoint);
38 }
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 affiche((char*)" Vérification de la conformité du maillage");
68 verification_conformite();
69 affiche((char*)" Adaptation de la carte de taille");
70 adapte_carte();
71 affiche((char*)" Determination des arêtes");
72 determine_arete();
73 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
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 if (seg->get_lien_triangle()->get_nb()!=2) continue;
112 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 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 if ((nbtri>2) || (nbtri<1))
284 {
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 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 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 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 void MAILLEUR_STL::calcul_nouveau_noeud(MG_SEGMENT* seg,double& x, double& y, double& z,OT_VECTEUR_3D &normale)
357 {
358 MG_TRIANGLE *tri1=seg->get_lien_triangle()->get(0);
359 MG_TRIANGLE *tri2=NULL;
360 if (seg->get_lien_triangle()->get_nb()>1) tri2=seg->get_lien_triangle()->get(1);
361 double xyz1[3],xyz2[3];
362 tri1->normal(xyz1);
363 if (tri2!=NULL) tri2->normal(xyz2);
364 OT_VECTEUR_3D nor1(xyz1);
365 if (tri2!=NULL)
366 {
367 OT_VECTEUR_3D nor2(xyz2);
368 normale=0.5*nor1+0.5*nor2;
369 }
370 else normale=nor1;
371 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 }
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 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 /*
435 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 affiche((char*)"etoile non conforme");
490 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 while (convergecoquille==false);*/
510 }
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 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 {
587 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 MG_NOEUD* nobase1=NULL;
594 MG_NOEUD* nobase2=NULL;
595 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 {
603 if (ptmp->second.second==2)
604 supprimer_seg_a_discretiser(ptmp->second.first);
605 if (ptmp->second.second==1)
606 if (ptmp->second.first==segbase)
607 supprimer_seg_a_discretiser(ptmp->second.first);
608 }
609 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 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 {
620 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 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 if (mgsegment1==NULL)
641 {//mormalement ne devrai jamais passer ici
642 mgsegment1=mai->ajouter_mg_segment(NULL,no1,no2,MAGIC::ORIGINE::TRIANGULATION);
643 ajouter_seg_a_discretiser(mgsegment1);
644 }
645 if (mgsegment2==NULL)
646 {
647 int origine=MAGIC::ORIGINE::TRIANGULATION;
648 if ((no2==nobase1)|| (no2==nobase2))
649 origine=MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE;
650 mgsegment2=mai->ajouter_mg_segment(NULL,no2,no,origine);
651 ajouter_seg_a_discretiser(mgsegment2);
652 }
653 if (mgsegment3==NULL)
654 {
655 int origine=MAGIC::ORIGINE::TRIANGULATION;
656 if ((no1==nobase1)|| (no1==nobase2))
657 origine=MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE;
658 mgsegment3=mai->ajouter_mg_segment(NULL,no,no1,origine);
659 ajouter_seg_a_discretiser(mgsegment3);
660 }
661 MG_TRIANGLE* nvtri=new MG_TRIANGLE(NULL,no1,no2,no,mgsegment1,mgsegment2,mgsegment3,MAGIC::ORIGINE::TRIANGULATION);
662 mai->ajouter_mg_triangle(nvtri);
663 listnouvtri.ajouter(nvtri);
664 /*double xyz[3];
665 nvtri->normal(xyz);
666 OT_VECTEUR_3D nvnor(xyz);
667 if (nvnor*normale<0)
668 nvtri->inverse_sens();*/
669 }
670
671 }
672 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 }
679
680
681 void MAILLEUR_STL::cree_maillage()
682 {
683 int ok=0;
684 int compteur_iter=0;
685 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
763
764
765 fclose(in);
766 }
767 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 TPL_LISTE_ENTITE<MG_TRIANGLE*> listtrietoile,listnvtri;
776 std::map<unsigned long,std::pair<MG_SEGMENT*,int> > lstsegmentetoile;
777 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 }
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
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 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 }