ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/aster/src/mgopt_posttraitement.cpp
Revision: 1104
Committed: Fri Sep 16 19:46:33 2022 UTC (2 years, 11 months ago) by francois
File size: 129587 byte(s)
Log Message:
Generalisation du calcul du Jacobien en 2D et 3D

File Contents

# Content
1 #include "gestionversion.h"
2 #include "mgopt_posttraitement.h"
3 #include "fem_solution.h"
4 #include "fem_maillage.h"
5 #include "mg_maillage.h"
6 #include "mg_gestionnaire.h"
7 #include "mg_triangle_peau.h"
8 #include "mg_geometrie_outils.h"
9 #include "fem_maillage_outils.h"
10 #include "tpl_fonctions_generiques.h"
11 #include "ot_algorithme_geometrique.h"
12 #include "mg_maillage_outils.h"
13 #include "mailleur_analyse.h"
14 #include "mg_arete_element.h"
15 #include "mg_face_element.h"
16 #include "mg_volume_element.h"
17 #include "mg_sommet_noeud.h"
18 #include "mg_export.h"
19 #include <stdio.h>
20 #include <stdlib.h>
21 #include <time.h>
22 #include <string.h>
23 #include <math.h>
24 #include <complex>
25 //ajout pour pouvoir utiliser les fonction venat de mailleur_analyse.cpp
26 #include "mailleur_analyse.h"
27 #include "m3d_triangle.h"
28 #include "fct_taille.h"
29 #include "mg_segment.h"
30 //---------------------------------------------------------------------------
31
32
33
34 MGOPT_POSTTRAITEMENT::MGOPT_POSTTRAITEMENT():affichageactif(0)
35 {
36 for (int i=0;i<9;i++) etat[i]=0;
37 params.ajouter("decoupage",0.,OT_PARAMETRES::DOUBLE,"0. découpage en conservant les mailles entieres 1. découpage par isodensite");
38 params.ajouter("seuil",0.5,OT_PARAMETRES::DOUBLE,"Valeur du seuil de densité pour les 2 découpages");
39 params.ajouter("nomdensextra","densiteextra.sol",OT_PARAMETRES::STRING,"Nom du fichier comportant la solution de la densité extrapolé au noeud");
40 params.ajouter("consimpose",1.,OT_PARAMETRES::DOUBLE,"0. Maillage non design non inclus dans le decoupage par mailles entieres 1. Maillage design inclus");
41 params.ajouter("consoptimise",1.,OT_PARAMETRES::DOUBLE,"0. Maillage optimise non inclus dans le decoupage par mailles entieres 1. Maillage optimise inclus");
42 params.ajouter("consnooptimise",0.,OT_PARAMETRES::DOUBLE,"0. Maillage non optimise non inclus dans le decoupage par mailles entieres 1. Maillage non optimise inclus");
43 params.ajouter("reactivation",1.,OT_PARAMETRES::DOUBLE,"1. Reactivation de mailles isolées sur le long des aretes 0. non reactivation");
44 params.ajouter("liss_aretes2017_nbiter",0.,OT_PARAMETRES::DOUBLE,"Nombre d'itérations du lissage des arêtes 0. méthode non utilisee");
45 params.ajouter("liss_aretes2017_lambda",0.8,OT_PARAMETRES::DOUBLE,"Paramète de relaxation lambda du lissage Laplacien des arêtes, la valeur doit être positive et moins de 1");
46 params.ajouter("chen2005_debut",0.,OT_PARAMETRES::DOUBLE,"Numero d'itération du debut d'application du lissage de chen2005 0. méthode non utilisee");
47 params.ajouter("chen2005_itermax",10.,OT_PARAMETRES::DOUBLE,"Nombre d'itérations maximales du lissage de chen2005");
48 params.ajouter("chen2005_epsilon",0.01,OT_PARAMETRES::DOUBLE,"Critere de convergence du lissage de chen2005");
49 params.ajouter("chen2005_sigma",0.1,OT_PARAMETRES::DOUBLE,"Écart type du filtre de lissage de chen2005");
50 params.ajouter("chen2005_filtre",1.,OT_PARAMETRES::DOUBLE,"Choix du filtre de lissage de chen2005 1. Gaussian 2.Laplacien 3.elfallahford");
51 params.ajouter("chen2008_debut",0.,OT_PARAMETRES::DOUBLE,"Numero d'itération du debut d'application deu lissage de chen2008 0. méthode non utilisee");
52 params.ajouter("chen2008_itermax",10.,OT_PARAMETRES::DOUBLE,"Nombre d'itérations maximales du lissage de chen2008");
53 params.ajouter("chen2008_epsilon",0.01,OT_PARAMETRES::DOUBLE,"Critere de convergence du lissage de chen2008");
54 params.ajouter("chen2008_sigma",0.1,OT_PARAMETRES::DOUBLE,"Écart type du filtre de lissage de chen2008");
55 params.ajouter("chen2008_gamma",1.,OT_PARAMETRES::DOUBLE,"Vitesse d'avancement de déplacement du lissage de chen2008");
56 params.ajouter("chen2008_filtre",1.,OT_PARAMETRES::DOUBLE,"Choix du filtre de lissage de chen2008 1. Gaussian 2.Laplacien 3.elfallahford");
57 params.ajouter("jiao2012_debut",0.,OT_PARAMETRES::DOUBLE,"Numero d'itération du debut d'application deu lissage de jiao2012 0. méthode non utilisee");
58 params.ajouter("jiao2012_itermax",10.,OT_PARAMETRES::DOUBLE,"Nombre d'itérations maximales du lissage de jioa2012");
59 params.ajouter("McKenzie2016_debut",0.,OT_PARAMETRES::DOUBLE,"Numero d'itération du debut d'application du lissage de McKenzie2016 0. méthode non utilisee");
60 params.ajouter("McKenzie2016_lambda",0.33,OT_PARAMETRES::DOUBLE,"Paramète lambda, la valeur doit être positive");
61 params.ajouter("McKenzie2016_nu",-0.34,OT_PARAMETRES::DOUBLE,"Paramète nu, la valeur doit être negative et respecter 0<lambda<(-nu)");
62 params.ajouter("McKenzie2016_epsilon", 0.01000, OT_PARAMETRES::DOUBLE,"Critere de convergence du lissage de chen2008");
63 params.ajouter("McKenzie2016_sigma", 0.10000, OT_PARAMETRES::DOUBLE,"Écart type du filtre de lissage de chen2008");
64 params.ajouter("McKenzie2016_gamma_", 1.0, OT_PARAMETRES::DOUBLE,"Vitesse d'avancement de déplacement du lissage de chen2008");
65 params.ajouter("McKenzie2016_filtre", 1.0, OT_PARAMETRES::DOUBLE,"Choix du filtre de lissage de chen2008 1. Gaussian 2.Laplacien 3.elfallahford");
66 params.ajouter("McKenzie2016_itertaubmax",100.,OT_PARAMETRES::DOUBLE,"Nombre d'itérations maximales de la partie taubin du lissage de McKenzie2016");
67 params.ajouter("McKenzie2016_itermax",100.,OT_PARAMETRES::DOUBLE,"Nombre d'itérations maximales de la partie chen 2008 du lissage de McKenzie2016");
68 params.ajouter("Taubin1995_debut",0.,OT_PARAMETRES::DOUBLE,"Numero d'itération du debut d'application du lissage de Taubin1995 0. méthode non utilisee");
69 params.ajouter("Taubin1995_lambda",0.33,OT_PARAMETRES::DOUBLE,"Paramète lambda, la valeur doit être positive et moins de 1");
70 params.ajouter("Taubin1995_nu",-0.34,OT_PARAMETRES::DOUBLE,"Paramète nu, la valeur doit être negative et respecter 0<lambda<(-nu)");
71 params.ajouter("Taubin1995_itermax",100.,OT_PARAMETRES::DOUBLE,"Nombre d'itérations maximales du lissage Taubin1995");
72 params.ajouter("cardinalite_debut",0.,OT_PARAMETRES::DOUBLE,"Numero d'itération du debut d'application de l'étude de cardinalite des noeuds");
73 params.ajouter("cardinalite_iter",2.,OT_PARAMETRES::DOUBLE,"Nombre d'itérations maximales de l'étude de cardinalite des noeuds");
74 params.ajouter("rmimpose_debut",0.,OT_PARAMETRES::DOUBLE,"0. Pas de suppression 1. Suppression du non_design du resultat optimal.");
75 params.ajouter("nomfichpeau","peausansnondesign.magic",OT_PARAMETRES::STRING,"Nom en .magic du fichier de sortie apres suppression du non design");
76 }
77
78
79 MGOPT_POSTTRAITEMENT::~MGOPT_POSTTRAITEMENT()
80 {
81 for (int i=0;i<lst_peau.size();i++)
82 delete lst_peau[i];
83 }
84
85 void MGOPT_POSTTRAITEMENT::active_affichage(void (*fonc)(char*))
86 {
87 affiche=fonc;
88 affichageactif=1;
89 }
90
91
92
93 void MGOPT_POSTTRAITEMENT::posttraite(char *fichier)
94 {
95 affiche((char*)"");
96 affiche((char*)"***********************************************************");
97 affiche((char*)"Post-traitement d'un calcul d'optimisation de topologie");
98 affiche((char*)"***********************************************************");
99 affiche((char*)"");
100 affiche((char*)"");
101 affiche((char*)"");
102 affiche((char*)"Écriture d'un fichier de parametres par défaut");
103 params.enregistrer(fichier);
104 }
105
106
107
108 void MGOPT_POSTTRAITEMENT::conserve(int origine)
109 {
110 etat[(origine-1000)/10]=1;
111 }
112
113
114 void MGOPT_POSTTRAITEMENT::copieorigine(FEM_MAILLAGE* mai)
115 {
116 LISTE_FEM_ELEMENT3::iterator it_tetra;
117 for (FEM_ELEMENT3* tet=mai->get_premier_element3(it_tetra);tet!=NULL;tet=mai->get_suivant_element3(it_tetra))
118 {
119 MG_TETRA* mgtet=(MG_TETRA*)tet->get_mg_element_maillage();
120 mgtet->change_origine(mgtet->get_origine());
121 }
122 }
123
124
125 void MGOPT_POSTTRAITEMENT::calcul_volume_ini(MG_MAILLAGE* mg_mai, double *volumetot,double *volumenondesign)
126 {
127 LISTE_MG_TETRA::iterator it_tetra;
128 double volume_initial = 0.;
129 double volume_final = 0.;
130 for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
131 {
132 int originetet=mgtet->get_origine();
133 MG_NOEUD* noeud1 = mgtet->get_noeud1();
134 double x1 = noeud1->get_x();
135 double y1 = noeud1->get_y();
136 double z1 = noeud1->get_z();
137 MG_NOEUD* noeud2 = mgtet->get_noeud2();
138 double x2 = noeud2->get_x();
139 double y2 = noeud2->get_y();
140 double z2 = noeud2->get_z();
141 MG_NOEUD* noeud3 = mgtet->get_noeud3();
142 double x3 = noeud3->get_x();
143 double y3 = noeud3->get_y();
144 double z3 = noeud3->get_z();
145 MG_NOEUD* noeud4 = mgtet->get_noeud4();
146 double x4 = noeud4->get_x();
147 double y4 = noeud4->get_y();
148 double z4 = noeud4->get_z();
149 double volume_tet = (x2-x1)*((y3-y1)*(z4-z1)-(z3-z1)*(y4-y1))-(y2-y1)*((x3-x1)*(z4-z1)-(z3-z1)*(x4-x1)) + (z2-z1)*((x3-x1)*(y4-y1)-(y3-y1)*(x4-x1));
150 volume_tet = volume_tet/6.;
151 volume_initial = volume_initial + volume_tet;
152 if (originetet == MAGIC::ORIGINE::IMPOSE)
153 {
154 volume_final = volume_final + volume_tet;
155 }
156 }
157 *volumetot=volume_initial;
158 *volumenondesign=volume_final;
159 }
160
161 void MGOPT_POSTTRAITEMENT::reactivation(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2)
162 {
163 std::cout << " Reactivation d'elements manquants sur aretes vives et faces planes" << std::endl;
164 LISTE_MG_TETRA::iterator it_tetra;
165 for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
166 {
167 int originetet=mgtet->get_origine();
168 int compteur = 0;
169 if (originetet == MAGIC::ORIGINE::MAILLEUR_AUTO)
170 {
171 MG_TRIANGLE* tri1 = mgtet->get_triangle1();
172 MG_TRIANGLE* tri2 = mgtet->get_triangle2();
173 MG_TRIANGLE* tri3 = mgtet->get_triangle3();
174 MG_TRIANGLE* tri4 = mgtet->get_triangle4();
175 int nb_tet1 = tri1->get_lien_tetra()->get_nb();
176 for (int k = 0;k<nb_tet1;k++)
177 {
178 MG_TETRA* tet1 = tri1->get_lien_tetra()->get(k);
179 int ori1 = tet1->get_origine();
180 if ((tet1 != mgtet) && (ori1 == MAGIC::ORIGINE::MAILLEUR_AUTO)) compteur++;
181 }
182 int nb_tet2 = tri2->get_lien_tetra()->get_nb();
183 for (int k = 0;k<nb_tet2;k++)
184 {
185 MG_TETRA* tet2 = tri2->get_lien_tetra()->get(k);
186 int ori2 = tet2->get_origine();
187 if ((tet2 != mgtet) && (ori2 == MAGIC::ORIGINE::MAILLEUR_AUTO)) compteur++;
188 }
189 int nb_tet3 = tri3->get_lien_tetra()->get_nb();
190 for (int k = 0;k<nb_tet3;k++)
191 {
192 MG_TETRA* tet3 = tri3->get_lien_tetra()->get(k);
193 int ori3 = tet3->get_origine();
194 if ((tet3 != mgtet) && (ori3 == MAGIC::ORIGINE::MAILLEUR_AUTO)) compteur++;
195 }
196 int nb_tet4 = tri4->get_lien_tetra()->get_nb();
197 for (int k = 0;k<nb_tet4;k++)
198 {
199 MG_TETRA* tet4 = tri4->get_lien_tetra()->get(k);
200 int ori4 = tet4->get_origine();
201 if ((tet4 != mgtet) && (ori4 == MAGIC::ORIGINE::MAILLEUR_AUTO)) compteur++;
202 }
203 if (compteur == 0) mgtet->change_origine(MAGIC::ORIGINE::OPTIMISE);
204 }
205 }
206 }
207
208
209
210
211 MG_MAILLAGE* MGOPT_POSTTRAITEMENT::extract_skin_par_decoupage(FEM_SOLUTION* sol,double limit,MG_GESTIONNAIRE& gest2,std::string nom)
212 {
213 affiche((char*)"Extraction de l'enveloppe par découpage");
214 sol->active_solution(0);
215 FEM_MAILLAGE* mai=sol->get_maillage();
216 affiche((char*)" Extrapolation de la densité aux noeuds");
217 LISTE_FEM_NOEUD::iterator it;
218 for (FEM_NOEUD* no=mai->get_premier_noeud(it);no!=NULL;no=mai->get_suivant_noeud(it))
219 {
220 double nume=0.;
221 double deno=0.;
222 int nb=no->get_lien_element3()->get_nb();
223 int passe=0;
224 for (int i=0;i<nb;i++)
225 {
226 FEM_ELEMENT3* tet=no->get_lien_element3()->get(i);
227 if (tet->get_mg_element_maillage()->get_origine()!=MAGIC::ORIGINE::IMPOSE)
228 {
229 double jac[9];
230 double volume=tet->get_jacobien(jac,jac,1.);
231 passe=1;
232 nume=nume+tet->get_solution()*volume;
233 deno=deno+volume;
234 }
235 }
236 if (passe==1) no->change_solution(nume/deno); else no->change_solution(1.);
237 }
238 if (nom!="")
239 {
240 affiche((char*)" Enregistrement de la densité aux noeuds");
241 MG_GESTIONNAIRE *gest=mai->get_mg_maillage()->get_gestionnaire();
242 std::string chemin=nom+".sol";
243 FEM_SOLUTION* femdens=new FEM_SOLUTION(mai,1,(char*)chemin.c_str(),mai->get_nb_fem_element3(),"Extrapolation aux noeuds",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT3_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
244 gest->ajouter_fem_solution(femdens);
245 femdens->change_legende(0,"Densite");
246 LISTE_FEM_ELEMENT3::iterator it3;
247 int i=0;
248 for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
249 {
250 MG_TETRA* tet2=(MG_TETRA*)tet->get_mg_element_maillage();
251 if (tet2->get_origine()==MAGIC::ORIGINE::IMPOSE)
252 {
253 femdens->ecrire(1.,i,0,0);
254 femdens->ecrire(1.,i,0,1);
255 femdens->ecrire(1.,i,0,2);
256 femdens->ecrire(1.,i,0,3);
257 }
258 else
259 {
260 femdens->ecrire(tet->get_fem_noeud(0)->get_solution(),i,0,0);
261 femdens->ecrire(tet->get_fem_noeud(1)->get_solution(),i,0,1);
262 femdens->ecrire(tet->get_fem_noeud(2)->get_solution(),i,0,2);
263 femdens->ecrire(tet->get_fem_noeud(3)->get_solution(),i,0,3);
264 }
265 i++;
266 }
267 MG_EXPORT exp;
268 exp.gmsh(mai,nom);
269 gest->supprimer_fem_solution(gest->get_nb_fem_solution()-1);
270 }
271 gest2.vide();
272 MG_GEOMETRIE* geo=new MG_GEOMETRIE((char*)"VIRTUEL",(char*)"VIRTUEL");
273 geo->change_valeur_unite(mai->get_mg_geometrie()->get_valeur_unite());
274 gest2.ajouter_mg_geometrie(geo);
275 MG_VOLUME_ELEMENT* vol=new MG_VOLUME_ELEMENT();
276 geo->ajouter_mg_volume(vol);
277 vol->copie_ccf(*(mai->get_mg_geometrie()->get_mg_volume(0)));
278 std::map<unsigned long,unsigned long> correspondidsom;
279 std::map<unsigned long,unsigned long> correspondid;
280 LISTE_MG_SOMMET::iterator its;
281 for (MG_SOMMET* som=mai->get_mg_geometrie()->get_premier_sommet(its);som!=NULL;som=mai->get_mg_geometrie()->get_suivant_sommet(its))
282 {
283 if (som->get_nb_ccf()>0)
284 {
285 MG_SOMMET_NOEUD* somno=new MG_SOMMET_NOEUD();
286 geo->ajouter_mg_sommet(somno);
287 somno->copie_ccf(*som);
288 correspondidsom[som->get_id()]=somno->get_id();
289 }
290 }
291 LISTE_MG_ARETE::iterator ita;
292 for (MG_ARETE* are=mai->get_mg_geometrie()->get_premier_arete(ita);are!=NULL;are=mai->get_mg_geometrie()->get_suivant_arete(ita))
293 {
294 if (are->get_nb_ccf()>0)
295 {
296 MG_ARETE_ELEMENT* areel=new MG_ARETE_ELEMENT();
297 geo->ajouter_mg_arete(areel);
298 areel->copie_ccf(*are);
299 correspondid[are->get_id()]=areel->get_id();
300 }
301 }
302 LISTE_MG_FACE::iterator itf;
303 for (MG_FACE* face=mai->get_mg_geometrie()->get_premier_face(itf);face!=NULL;face=mai->get_mg_geometrie()->get_suivant_face(itf))
304 {
305 if (face->get_nb_ccf()>0)
306 {
307 MG_FACE_ELEMENT* facel=new MG_FACE_ELEMENT();
308 geo->ajouter_mg_face(facel);
309 facel->copie_ccf(*face);
310 correspondid[face->get_id()]=facel->get_id();
311 }
312 }
313 MG_FACE_ELEMENT* faceex=new MG_FACE_ELEMENT();
314 geo->ajouter_mg_face(faceex);
315 MG_FACE_ELEMENT* facenoex=new MG_FACE_ELEMENT();
316 geo->ajouter_mg_face(facenoex);
317 MG_MAILLAGE* mgmai=new MG_MAILLAGE(geo);
318 gest2.ajouter_mg_maillage(mgmai);
319
320 MG_MAILLAGE* mgfem=mai->get_mg_maillage();
321 LISTE_MG_TRIANGLE::iterator itmg;
322 for (MG_TRIANGLE* tri=mgfem->get_premier_triangle(itmg);tri!=NULL;tri=mgfem->get_suivant_triangle(itmg))
323 tri->change_nouveau_numero(0);
324 affiche((char*)" Decoupage des tetra optimises");
325 LISTE_FEM_ELEMENT3::iterator it3;
326 for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
327 {
328 MG_TETRA* tet2=(MG_TETRA*)tet->get_mg_element_maillage();
329 if (tet2->get_origine()==MAGIC::ORIGINE::IMPOSE)
330 {
331 tet2->get_triangle1()->change_origine(MAGIC::ORIGINE::IMPOSE);
332 tet2->get_triangle2()->change_origine(MAGIC::ORIGINE::IMPOSE);
333 tet2->get_triangle3()->change_origine(MAGIC::ORIGINE::IMPOSE);
334 tet2->get_triangle4()->change_origine(MAGIC::ORIGINE::IMPOSE);
335 tet2->get_noeud1()->change_solution(tet->get_fem_noeud(0)->get_solution());
336 tet2->get_noeud2()->change_solution(tet->get_fem_noeud(1)->get_solution());
337 tet2->get_noeud3()->change_solution(tet->get_fem_noeud(2)->get_solution());
338 tet2->get_noeud4()->change_solution(tet->get_fem_noeud(3)->get_solution());
339 tet2->get_noeud1()->change_nouveau_numero(tet->get_fem_noeud(0)->get_id());
340 tet2->get_noeud2()->change_nouveau_numero(tet->get_fem_noeud(1)->get_id());
341 tet2->get_noeud3()->change_nouveau_numero(tet->get_fem_noeud(2)->get_id());
342 tet2->get_noeud4()->change_nouveau_numero(tet->get_fem_noeud(3)->get_id());
343 if (tet2->get_triangle1()->get_nouveau_numero()==0)
344 tet2->get_triangle1()->change_nouveau_numero(tet2->get_noeud4()->get_id());
345 else
346 tet2->get_triangle1()->change_nouveau_numero(-1);
347 if (tet2->get_triangle2()->get_nouveau_numero()==0)
348 tet2->get_triangle2()->change_nouveau_numero(tet2->get_noeud3()->get_id());
349 else
350 tet2->get_triangle2()->change_nouveau_numero(-1);
351 if (tet2->get_triangle3()->get_nouveau_numero()==0)
352 tet2->get_triangle3()->change_nouveau_numero(tet2->get_noeud1()->get_id());
353 else
354 tet2->get_triangle3()->change_nouveau_numero(-1);
355 if (tet2->get_triangle4()->get_nouveau_numero()==0)
356 tet2->get_triangle4()->change_nouveau_numero(tet2->get_noeud2()->get_id());
357 else
358 tet2->get_triangle4()->change_nouveau_numero(-1);
359 }
360 }
361
362 TPL_LISTE_ENTITE<MG_TRIANGLE*> tri_impose_interne;
363 for (MG_TRIANGLE* tri=mgfem->get_premier_triangle(itmg);tri!=NULL;tri=mgfem->get_suivant_triangle(itmg))
364 {
365 if (tri->get_origine()==MAGIC::ORIGINE::IMPOSE)
366 if (tri->get_lien_topologie()->get_dimension()==3)
367 if (tri->get_nouveau_numero()>0)
368 tri_impose_interne.ajouter(tri);
369 }
370 for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
371 {
372 if (tet->get_mg_element_maillage()->get_origine()!=MAGIC::ORIGINE::IMPOSE)
373 {
374 FEM_NOEUD* no1=tet->get_fem_noeud(0);
375 FEM_NOEUD* no2=tet->get_fem_noeud(1);
376 FEM_NOEUD* no3=tet->get_fem_noeud(2);
377 FEM_NOEUD* no4=tet->get_fem_noeud(3);
378 std::vector<MG_NOEUD*> tab;
379 interpole_segment(mai,no1,no2,&tab,limit,mgmai);
380 interpole_segment(mai,no1,no3,&tab,limit,mgmai);
381 interpole_segment(mai,no1,no4,&tab,limit,mgmai);
382 interpole_segment(mai,no2,no3,&tab,limit,mgmai);
383 interpole_segment(mai,no2,no4,&tab,limit,mgmai);
384 interpole_segment(mai,no3,no4,&tab,limit,mgmai);
385 int nb=tab.size();
386 FEM_NOEUD* noext;
387 if (nb>0)
388 {
389 if (no1->get_solution()<limit) noext=no1;
390 if (no2->get_solution()<limit) noext=no2;
391 if (no3->get_solution()<limit) noext=no3;
392 if (no4->get_solution()<limit) noext=no4;
393 }
394 if (nb==3)
395 {
396 MG_TRIANGLE_PEAU* tri=insere_triangle(facenoex,tab[0],tab[1],tab[2],mgmai,MAGIC::ORIGINE::TRIANGULATION);
397 oriente_tri(tri,noext->get_coord());
398 }
399 if (nb==4)
400 {
401 if (test_du_point_milieu(tab[0],tab[1],tet)==1)
402 {
403 MG_TRIANGLE_PEAU* tri=insere_triangle(facenoex,tab[0],tab[1],tab[2],mgmai,MAGIC::ORIGINE::TRIANGULATION);
404 MG_TRIANGLE_PEAU* tri2=insere_triangle(facenoex,tab[0],tab[1],tab[3],mgmai,MAGIC::ORIGINE::TRIANGULATION);
405 oriente_tri(tri,noext->get_coord());
406 oriente_tri(tri2,noext->get_coord());
407 }
408 else if (test_du_point_milieu(tab[0],tab[2],tet)==1)
409 {
410 MG_TRIANGLE_PEAU* tri=insere_triangle(facenoex,tab[0],tab[2],tab[1],mgmai,MAGIC::ORIGINE::TRIANGULATION);
411 MG_TRIANGLE_PEAU* tri2=insere_triangle(facenoex,tab[0],tab[2],tab[3],mgmai,MAGIC::ORIGINE::TRIANGULATION);
412 oriente_tri(tri,noext->get_coord());
413 oriente_tri(tri2,noext->get_coord());}
414 else if (test_du_point_milieu(tab[0],tab[3],tet)==1)
415 {
416 MG_TRIANGLE_PEAU* tri=insere_triangle(facenoex,tab[0],tab[3],tab[1],mgmai,MAGIC::ORIGINE::TRIANGULATION);
417 MG_TRIANGLE_PEAU* tri2=insere_triangle(facenoex,tab[0],tab[3],tab[2],mgmai,MAGIC::ORIGINE::TRIANGULATION);
418 oriente_tri(tri,noext->get_coord());
419 oriente_tri(tri2,noext->get_coord());}
420 else if (test_du_point_milieu(tab[1],tab[2],tet)==1)
421 {
422 MG_TRIANGLE_PEAU* tri=insere_triangle(facenoex,tab[1],tab[2],tab[0],mgmai,MAGIC::ORIGINE::TRIANGULATION);
423 MG_TRIANGLE_PEAU* tri2=insere_triangle(facenoex,tab[1],tab[2],tab[3],mgmai,MAGIC::ORIGINE::TRIANGULATION);
424 oriente_tri(tri,noext->get_coord());
425 oriente_tri(tri2,noext->get_coord());
426 }
427 else if (test_du_point_milieu(tab[1],tab[3],tet)==1)
428 {
429 MG_TRIANGLE_PEAU* tri=insere_triangle(facenoex,tab[1],tab[3],tab[0],mgmai,MAGIC::ORIGINE::TRIANGULATION);
430 MG_TRIANGLE_PEAU* tri2=insere_triangle(facenoex,tab[1],tab[3],tab[2],mgmai,MAGIC::ORIGINE::TRIANGULATION);
431 oriente_tri(tri,noext->get_coord());
432 oriente_tri(tri2,noext->get_coord());}
433 else
434 {
435 MG_TRIANGLE_PEAU* tri=insere_triangle(facenoex,tab[2],tab[3],tab[0],mgmai,MAGIC::ORIGINE::TRIANGULATION);
436 MG_TRIANGLE_PEAU* tri2=insere_triangle(facenoex,tab[2],tab[3],tab[1],mgmai,MAGIC::ORIGINE::TRIANGULATION);
437 oriente_tri(tri,noext->get_coord());
438 oriente_tri(tri2,noext->get_coord());
439 }
440
441 }
442 }
443 }
444 affiche((char*)" Decoupage des triangles non design interne au volume");
445 for (int i=0;i<tri_impose_interne.get_nb();i++)
446 {
447 MG_TRIANGLE* tri=tri_impose_interne.get(i);
448 MG_NOEUD* no1=tri->get_noeud1();
449 MG_NOEUD* no2=tri->get_noeud2();
450 MG_NOEUD* no3=tri->get_noeud3();
451 std::vector<MG_NOEUD*> tab;
452 int num1=-1;
453 int num2=-1;
454 int num3=-1;
455 if (no1->get_solution()<limit)
456 {
457 MG_NOEUD* no=get_noeud_peau(mai->get_fem_noeudid(no1->get_nouveau_numero()),mgmai,correspondidsom,geo);
458 tab.push_back(no);
459 num1=1;
460 }
461 if (no2->get_solution()<limit)
462 {
463 MG_NOEUD* no=get_noeud_peau(mai->get_fem_noeudid(no2->get_nouveau_numero()),mgmai,correspondidsom,geo);
464 tab.push_back(no);
465 num2=1;
466 }
467 if (no3->get_solution()<limit)
468 {
469 MG_NOEUD* no=get_noeud_peau(mai->get_fem_noeudid(no3->get_nouveau_numero()),mgmai,correspondidsom,geo);
470 tab.push_back(no);
471 num3=1;
472 }
473 if (num1*num2<0) interpole_segment(mai,mai->get_fem_noeudid(no1->get_nouveau_numero()),mai->get_fem_noeudid(no2->get_nouveau_numero()),&tab,limit,mgmai);
474 if (num1*num3<0) interpole_segment(mai,mai->get_fem_noeudid(no1->get_nouveau_numero()),mai->get_fem_noeudid(no3->get_nouveau_numero()),&tab,limit,mgmai);
475 if (num2*num3<0) interpole_segment(mai,mai->get_fem_noeudid(no2->get_nouveau_numero()),mai->get_fem_noeudid(no3->get_nouveau_numero()),&tab,limit,mgmai);
476 MG_NOEUD* noint=mgfem->get_mg_noeudid(tri->get_nouveau_numero());
477 int nb=tab.size();
478 if (nb==3)
479 {
480 MG_TRIANGLE_PEAU* tri=insere_triangle(facenoex,tab[0],tab[1],tab[2],mgmai,MAGIC::ORIGINE::IMPOSE);
481 oriente_tri(tri,noint->get_coord());
482 tri->inverse_sens();
483 }
484 if (nb==4)
485 {
486 MG_TRIANGLE_PEAU* tri=insere_triangle(facenoex,tab[0],tab[1],tab[2],mgmai,MAGIC::ORIGINE::IMPOSE);
487 MG_TRIANGLE_PEAU* tri2=insere_triangle(facenoex,tab[1],tab[2],tab[3],mgmai,MAGIC::ORIGINE::IMPOSE);
488 oriente_tri(tri,noint->get_coord());
489 tri->inverse_sens();
490 oriente_tri(tri2,noint->get_coord());
491 tri2->inverse_sens();
492 }
493 }
494 affiche((char*)" Decoupage des triangles sur la frontiere du volume");
495 LISTE_FEM_ELEMENT2::iterator it2;
496 for (FEM_ELEMENT2* tri=mai->get_premier_element2(it2);tri!=NULL;tri=mai->get_suivant_element2(it2))
497 {
498 std::vector<MG_NOEUD*> tab;
499 FEM_NOEUD *no1=tri->get_fem_noeud(0);
500 FEM_NOEUD *no2=tri->get_fem_noeud(1);
501 FEM_NOEUD *no3=tri->get_fem_noeud(2);
502 int ori=((MG_FACE*)tri->get_lien_topologie())->get_mg_coface(0)->get_orientation();
503 OT_VECTEUR_3D n1n2(no1->get_coord(),no2->get_coord());
504 OT_VECTEUR_3D n1n3(no1->get_coord(),no3->get_coord());
505 OT_VECTEUR_3D nor=n1n3&n1n2;
506 double xyzext[3];
507 xyzext[0]=no1->get_x()+nor.get_x()*ori;
508 xyzext[1]=no1->get_y()+nor.get_y()*ori;
509 xyzext[2]=no1->get_z()+nor.get_z()*ori;
510 if (tri->get_mg_element_maillage()->get_origine()==MAGIC::ORIGINE::IMPOSE)
511 {
512 MG_NOEUD *mgno1=get_noeud_peau(no1,mgmai,correspondidsom,geo);
513 MG_NOEUD *mgno2=get_noeud_peau(no2,mgmai,correspondidsom,geo);
514 MG_NOEUD *mgno3=get_noeud_peau(no3,mgmai,correspondidsom,geo);
515 int num1=interpole_segment(mai,no1,no2,&tab,limit,mgmai,0);
516 int num2=interpole_segment(mai,no1,no3,&tab,limit,mgmai,0);
517 int num3=interpole_segment(mai,no2,no3,&tab,limit,mgmai,0);
518 int nb=tab.size();
519 MG_MAILLAGE* maiori=mai->get_mg_maillage();
520 MG_NOEUD* n1ori=(MG_NOEUD*)no1->get_mg_element_maillage();
521 MG_NOEUD* n2ori=(MG_NOEUD*)no2->get_mg_element_maillage();
522 MG_NOEUD* n3ori=(MG_NOEUD*)no3->get_mg_element_maillage();
523 MG_FACE* face=faceex;
524 std::map<unsigned long,unsigned long>::iterator it=correspondid.find(tri->get_mg_element_maillage()->get_lien_topologie()->get_id());
525 if (it!=correspondid.end())
526 face=geo->get_mg_faceid(it->second);
527 if (nb==0)
528 {
529 MG_TRIANGLE_PEAU* tri=insere_triangle(geo,maiori,n1ori,n2ori,n3ori,correspondid,face,mgno1,mgno2,mgno3,mgmai,MAGIC::ORIGINE::IMPOSE);
530 oriente_tri(tri,xyzext);
531 }
532 if (nb==1)
533 {
534 if (num1==1)
535 {
536 MG_TRIANGLE_PEAU* tri=insere_triangle(geo,maiori,n1ori,n2ori,n3ori,correspondid,face,mgno1,mgno3,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
537 MG_TRIANGLE_PEAU* tri2=insere_triangle(geo,maiori,n1ori,n2ori,n3ori,correspondid,face,mgno2,mgno3,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
538 oriente_tri(tri,xyzext);
539 oriente_tri(tri2,xyzext);
540 }
541 if (num2==1)
542 {
543 MG_TRIANGLE_PEAU* tri=insere_triangle(geo,maiori,n1ori,n2ori,n3ori,correspondid,face,mgno1,mgno2,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
544 MG_TRIANGLE_PEAU* tri2=insere_triangle(geo,maiori,n1ori,n2ori,n3ori,correspondid,face,mgno2,mgno3,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
545 oriente_tri(tri,xyzext);
546 oriente_tri(tri2,xyzext);
547 }
548 if (num3==1)
549 {
550 MG_TRIANGLE_PEAU* tri=insere_triangle(geo,maiori,n1ori,n2ori,n3ori,correspondid,face,mgno1,mgno2,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
551 MG_TRIANGLE_PEAU* tri2=insere_triangle(geo,maiori,n1ori,n2ori,n3ori,correspondid,face,mgno1,mgno3,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
552 oriente_tri(tri,xyzext);
553 oriente_tri(tri2,xyzext);
554 }
555 }
556 if (nb==2)
557 {
558 if (num1==0)
559 {
560 MG_TRIANGLE_PEAU* tri=insere_triangle(geo,maiori,n1ori,n2ori,n3ori,correspondid,face,mgno1,mgno2,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
561 MG_TRIANGLE_PEAU* tri2=insere_triangle(geo,maiori,n1ori,n2ori,n3ori,correspondid,face,mgno2,tab[0],tab[1],mgmai,MAGIC::ORIGINE::IMPOSE);
562 MG_TRIANGLE_PEAU* tri3=insere_triangle(geo,maiori,n1ori,n2ori,n3ori,correspondid,face,mgno3,tab[0],tab[1],mgmai,MAGIC::ORIGINE::IMPOSE);
563 oriente_tri(tri,xyzext);
564 oriente_tri(tri2,xyzext);
565 oriente_tri(tri3,xyzext);
566 }
567 if (num2==0)
568 {
569 MG_TRIANGLE_PEAU* tri=insere_triangle(geo,maiori,n1ori,n2ori,n3ori,correspondid,face,mgno1,mgno3,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
570 MG_TRIANGLE_PEAU* tri2=insere_triangle(geo,maiori,n1ori,n2ori,n3ori,correspondid,face,mgno3,tab[0],tab[1],mgmai,MAGIC::ORIGINE::IMPOSE);
571 MG_TRIANGLE_PEAU* tri3=insere_triangle(geo,maiori,n1ori,n2ori,n3ori,correspondid,face,mgno2,tab[0],tab[1],mgmai,MAGIC::ORIGINE::IMPOSE);
572 oriente_tri(tri,xyzext);
573 oriente_tri(tri2,xyzext);
574 oriente_tri(tri3,xyzext);
575 }
576 if (num2==0)
577 {
578 MG_TRIANGLE_PEAU* tri=insere_triangle(geo,maiori,n1ori,n2ori,n3ori,correspondid,face,mgno2,mgno3,tab[0],mgmai,MAGIC::ORIGINE::IMPOSE);
579 MG_TRIANGLE_PEAU* tri2=insere_triangle(geo,maiori,n1ori,n2ori,n3ori,correspondid,face,mgno3,tab[0],tab[1],mgmai,MAGIC::ORIGINE::IMPOSE);
580 MG_TRIANGLE_PEAU* tri3=insere_triangle(geo,maiori,n1ori,n2ori,n3ori,correspondid,face,mgno1,tab[0],tab[1],mgmai,MAGIC::ORIGINE::IMPOSE);
581 oriente_tri(tri,xyzext);
582 oriente_tri(tri2,xyzext);
583 oriente_tri(tri3,xyzext);
584 }
585 }
586 }
587 else if ((no1->get_solution()>=limit) && (no2->get_solution()>=limit) && (no3->get_solution()>=limit))
588 {
589 MG_NOEUD *mgno1=get_noeud_peau(no1,mgmai,correspondidsom,geo);
590 MG_NOEUD *mgno2=get_noeud_peau(no2,mgmai,correspondidsom,geo);
591 MG_NOEUD *mgno3=get_noeud_peau(no3,mgmai,correspondidsom,geo);
592 MG_TRIANGLE_PEAU* tri=insere_triangle(faceex,mgno1,mgno2,mgno3,mgmai,MAGIC::ORIGINE::TRIANGULATION);
593 oriente_tri(tri,xyzext);
594 }
595 else if (!((no1->get_solution()<limit) && (no2->get_solution()<limit) && (no3->get_solution()<limit)))
596 {
597 std::vector<MG_NOEUD*> tab;
598 int num1=-1;
599 int num2=-1;
600 int num3=-1;
601 if (no1->get_solution()>=limit)
602 {
603 MG_NOEUD* no=get_noeud_peau(no1,mgmai,correspondidsom,geo);
604 tab.push_back(no);
605 num1=1;
606 }
607 if (no2->get_solution()>=limit)
608 {
609 MG_NOEUD* no=get_noeud_peau(no2,mgmai,correspondidsom,geo);
610 tab.push_back(no);
611 num2=1;
612 }
613 if (no3->get_solution()>=limit)
614 {
615 MG_NOEUD* no=get_noeud_peau(no3,mgmai,correspondidsom,geo);
616 tab.push_back(no);
617 num3=1;
618 }
619 if (num1*num2<0) interpole_segment(mai,no1,no2,&tab,limit,mgmai);
620 if (num1*num3<0) interpole_segment(mai,no1,no3,&tab,limit,mgmai);
621 if (num2*num3<0) interpole_segment(mai,no2,no3,&tab,limit,mgmai);
622 int nb=tab.size();
623 if (nb==3)
624 {
625 MG_TRIANGLE_PEAU* tri=insere_triangle(faceex,tab[0],tab[1],tab[2],mgmai,MAGIC::ORIGINE::TRIANGULATION);
626 oriente_tri(tri,xyzext);
627 }
628 if (nb==4)
629 {
630 MG_TRIANGLE_PEAU* tri=insere_triangle(faceex,tab[0],tab[1],tab[2],mgmai,MAGIC::ORIGINE::TRIANGULATION);
631 MG_TRIANGLE_PEAU* tri2=insere_triangle(faceex,tab[1],tab[2],tab[3],mgmai,MAGIC::ORIGINE::TRIANGULATION);
632 oriente_tri(tri,xyzext);
633 oriente_tri(tri2,xyzext);
634 }
635 }
636 }
637 LISTE_MG_TRIANGLE::iterator it_tri;
638 for (MG_TRIANGLE* mgtri=mgmai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mgmai->get_suivant_triangle(it_tri))
639 {
640 MG_TRIANGLE_PEAU* tripeau=(MG_TRIANGLE_PEAU*)mgtri;
641 MG_TRIANGLE_PEAU* voisin1=recherche_voisin(tripeau->get_noeud1(),tripeau->get_noeud2(),tripeau);
642 MG_TRIANGLE_PEAU* voisin2=recherche_voisin(tripeau->get_noeud2(),tripeau->get_noeud3(),tripeau);
643 MG_TRIANGLE_PEAU* voisin3=recherche_voisin(tripeau->get_noeud3(),tripeau->get_noeud1(),tripeau);
644 tripeau->change_voisin1(voisin1);
645 tripeau->change_voisin2(voisin2);
646 tripeau->change_voisin3(voisin3);
647 tripeau->change_nouveau_numero(0);
648 }
649 int fin;
650 do
651 {
652 fin=1;
653 for (MG_TRIANGLE* mgtri=mgmai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mgmai->get_suivant_triangle(it_tri))
654 {
655 MG_TRIANGLE_PEAU *tripeau=(MG_TRIANGLE_PEAU*)mgtri;
656 if (tripeau->get_nouveau_numero()==0)
657 {
658 fin=0;
659 std::vector<MG_TRIANGLE_PEAU*> *peau=new std::vector<MG_TRIANGLE_PEAU*>;
660 lst_peau.push_back(peau);
661 tripeau->change_nouveau_numero(1);
662 peau->push_back(tripeau);
663 determine_peau(peau);
664 }
665 }
666 }
667 while (fin==0);
668 return mgmai;
669 }
670
671
672 void MGOPT_POSTTRAITEMENT::suppression_peaux_isoles(MG_MAILLAGE* mgmai)
673 {
674 affiche((char*)"Suppression des peaux isolées");
675
676 char message[500];
677 for (int cas=0;cas<2;cas++)
678 {
679 if (cas==0) affiche((char*)" Analyse initiale des peaux");
680 if (cas==1) affiche((char*)" Analyse finale des peaux");
681 int nbisole=0;
682 for (int i=0;i<lst_peau.size();i++)
683 {
684 int isole=1;
685 for (int j=0;j<lst_peau[i]->size();j++)
686 if ((*lst_peau[i])[j]->get_origine()==MAGIC::ORIGINE::IMPOSE) {isole=0;break;}
687 if (isole==1)
688 {
689 nbisole++;
690 for (int j=0;j<lst_peau[i]->size();j++)
691 {
692 mgmai->supprimer_mg_triangleid((*lst_peau[i])[j]->get_id());
693 }
694 lst_peau[i]->clear();
695 }
696 }
697 sprintf(message," %d peaux, %d non isoles, %d isoles",(int)lst_peau.size(),(int)lst_peau.size()-nbisole,nbisole);
698 affiche(message);
699 for (std::vector<std::vector<MG_TRIANGLE_PEAU*> *>::iterator j=lst_peau.end()-1;j>lst_peau.begin();j--)
700 if ((*j)->size()==0)
701 lst_peau.erase(j);
702 }
703 LISTE_MG_TRIANGLE::iterator itmg;
704 for (MG_TRIANGLE* tri=mgmai->get_premier_triangle(itmg);tri!=NULL;tri=mgmai->get_suivant_triangle(itmg))
705 {
706 if (tri->get_origine()==MAGIC::ORIGINE::IMPOSE)
707 {
708 tri->get_noeud1()->change_origine(MAGIC::ORIGINE::IMPOSE);
709 tri->get_noeud2()->change_origine(MAGIC::ORIGINE::IMPOSE);
710 tri->get_noeud3()->change_origine(MAGIC::ORIGINE::IMPOSE);
711 }
712 }
713
714 }
715
716
717 void MGOPT_POSTTRAITEMENT::oriente_tri(MG_TRIANGLE_PEAU* tri,double *xyz)
718 {
719 MG_NOEUD* no1=tri->get_noeud1();
720 MG_NOEUD* no2=tri->get_noeud2();
721 MG_NOEUD* no3=tri->get_noeud3();
722 OT_VECTEUR_3D n1n2(no1->get_coord(),no2->get_coord());
723 OT_VECTEUR_3D n1n3(no1->get_coord(),no3->get_coord());
724 OT_VECTEUR_3D normal=n1n3&n1n2;
725 OT_VECTEUR_3D dir(no1->get_coord(),xyz);
726 dir.norme();
727 normal.norme();
728 double ps=normal*dir;
729 if (ps<0) tri->inverse_sens();
730 }
731
732
733
734 MG_NOEUD* MGOPT_POSTTRAITEMENT::get_noeud_peau(FEM_NOEUD* no,MG_MAILLAGE* mai,std::map<unsigned long,unsigned long> &correspondid,MG_GEOMETRIE* geo)
735 {
736 static std::map<std::string,MG_NOEUD*> correspond;
737 unsigned long id=no->get_id();
738 char message[255];
739 sprintf(message,"_%lu_",id);
740 std::string key=message;
741 MG_NOEUD* mgno=correspond[key];
742 if (mgno==NULL)
743 {
744 int origine=MAGIC::ORIGINE::TRIANGULATION;
745 if (((MG_NOEUD*)no->get_mg_element_maillage())->get_lien_topologie()->get_dimension()==1) origine=MAGIC::ORIGINE::TRIANGULATION_ARETE;
746 MG_ELEMENT_TOPOLOGIQUE* ele=NULL;
747 std::map<unsigned long,unsigned long>::iterator it=correspondid.find(no->get_mg_element_maillage()->get_lien_topologie()->get_id());
748 if (it!=correspondid.end())
749 ele=geo->get_mg_sommetid(it->second);
750 mgno=new MG_NOEUD(NULL,no->get_x(),no->get_y(),no->get_z(),origine);
751 if (it!=correspondid.end())
752 ((MG_SOMMET_NOEUD*)ele)->change_mg_noeud(mgno);
753 mai->ajouter_mg_noeud(mgno);
754 correspond[key]=mgno;
755 }
756 return mgno;
757 }
758
759
760 int MGOPT_POSTTRAITEMENT::test_du_point_milieu(MG_NOEUD* no1,MG_NOEUD* no2,FEM_ELEMENT3* tet)
761 {
762 double *xyz1=no1->get_coord();
763 double *xyz2=no2->get_coord();
764 double xyz[3];
765 xyz[0]=0.5*(xyz1[0]+xyz2[0]);
766 xyz[1]=0.5*(xyz1[1]+xyz2[1]);
767 xyz[2]=0.5*(xyz1[2]+xyz2[2]);
768 FEM_MAILLAGE_OUTILS ot;
769 if (((ot.estdansletetra(tet,xyz[0],xyz[1],xyz[2])>>1)&1)==1) return 1;
770 //if (((MG_MAILLAGE_OUTILS::estdansletetra(tet,xyz[0],xyz[1],xyz[2])>>1)&1)==1) return 1;
771 return 0;
772 }
773
774
775
776 int MGOPT_POSTTRAITEMENT::interpole_segment(FEM_MAILLAGE* fem,FEM_NOEUD* no1,FEM_NOEUD* no2,std::vector<MG_NOEUD*> *tab,double limit,MG_MAILLAGE* mai,int creation)
777 {
778 static std::map<std::string,MG_NOEUD*> correspond;
779 unsigned long id1=no1->get_id();
780 unsigned long id2=no2->get_id();
781 char message[255];
782 sprintf(message,"_%lu_%lu_",id1,id2);
783 std::string key1=message;
784 sprintf(message,"_%lu_%lu_",id2,id1);
785 std::string key2=message;
786 double sol1=no1->get_solution();
787 double sol2=no2->get_solution();
788 MG_NOEUD* n1=(MG_NOEUD*)no1->get_mg_element_maillage();
789 MG_NOEUD* n2=(MG_NOEUD*)no2->get_mg_element_maillage();
790 MG_SEGMENT *seg=fem->get_mg_maillage()->get_mg_segment(n1->get_id(),n2->get_id());
791 int origine=MAGIC::ORIGINE::TRIANGULATION;
792 if (seg->get_lien_topologie()->get_dimension()==2) origine=MAGIC::ORIGINE::TRIANGULATION_ARETE;
793 if (seg->get_lien_topologie()->get_dimension()==1) origine=MAGIC::ORIGINE::TRIANGULATION_ARETE;
794
795 if (fabs(sol1-limit)<0.0000001)
796 {
797 sprintf(message,"%lu",id1);
798 std::string key=message;
799 MG_NOEUD* no=correspond[key];
800 if (no==NULL)
801 {
802 double *xyz1=no1->get_coord();
803 no=new MG_NOEUD(NULL,xyz1[0],xyz1[1],xyz1[2],origine);
804 mai->ajouter_mg_noeud(no);
805 correspond[key]=no;
806 }
807 int present=0;
808 for (int i=0;i<tab->size();i++)
809 if ((*tab)[i]==no) present=1;
810 if (present==0)
811 {
812 tab->push_back(no);
813 return 1;
814 }
815 else return 0;
816 }
817 if (fabs(sol2-limit)<0.0000001)
818 {
819 sprintf(message,"%lu",id2);
820 std::string key=message;
821 MG_NOEUD* no=correspond[key];
822 if (no==NULL)
823 {
824 double *xyz2=no2->get_coord();
825 no=new MG_NOEUD(NULL,xyz2[0],xyz2[1],xyz2[2],origine);
826 mai->ajouter_mg_noeud(no);
827 correspond[key]=no;
828 }
829 int present=0;
830 for (int i=0;i<tab->size();i++)
831 if ((*tab)[i]==no) present=1;
832 if (present==0)
833 {
834 tab->push_back(no);
835 return 1;
836 }
837 else return 0;
838
839 }
840
841 if (((sol1<limit) && (sol2>limit))||((sol1>limit) && (sol2<limit)))
842 {
843 MG_NOEUD* no=correspond[key1];
844 if (no==NULL) no=correspond[key2];
845 if ((no==NULL) && (creation==1))
846 {
847 double t=(limit-sol1)/(sol2-sol1);
848 double xyz[3];
849 double *xyz1=no1->get_coord();
850 double *xyz2=no2->get_coord();
851 xyz[0]=xyz1[0]+t*(xyz2[0]-xyz1[0]);
852 xyz[1]=xyz1[1]+t*(xyz2[1]-xyz1[1]);
853 xyz[2]=xyz1[2]+t*(xyz2[2]-xyz1[2]);
854 no=new MG_NOEUD(NULL,xyz[0],xyz[1],xyz[2],origine);
855 mai->ajouter_mg_noeud(no);
856 correspond[key1]=no;
857 }
858 if (no!=NULL)
859 {
860 tab->push_back(no);
861 return 1;
862 }
863 }
864 return 0;
865 }
866
867 void MGOPT_POSTTRAITEMENT::lire_params(char *fichier)
868 {
869 params.lire(fichier);
870 }
871
872 void MGOPT_POSTTRAITEMENT::posttraite(FEM_SOLUTION* sol, MG_GESTIONNAIRE& gest2,char *nomparam,char *nomout)
873 {
874 strcat(nomout,".volume");
875 FILE *in=fopen(nomout,"wt");
876 fprintf(in,"Suivi du volume pendant le posttraitement d'une optimisation de topologie\n");
877 fprintf(in,"**************************************************************************\n\n\n\n");
878 if (nomparam!=NULL) lire_params(nomparam);
879 FEM_MAILLAGE* fem=sol->get_maillage();
880 affiche((char*)"Calcul du volume de depart");
881 double voltot=0.;
882 double volnon=0.;
883 calcul_volume_ini(fem->get_mg_maillage(),&voltot,&volnon);
884 fprintf(in,"Statistique initiale\n");
885 fprintf(in,"Volume total=%.14le\nVolume non design=%.14le\nVolume à optimiser=%.14le\n\n\n\n",voltot,volnon,voltot-volnon);
886 affiche((char*)"Extraction du maillage de surface");
887 int decoup=params.get_valeur("decoupage");
888 MG_MAILLAGE* mgmai;
889 if (decoup==0)
890 {
891 affiche((char*)" Découpage par mailles 3D entieres");
892 adapte_seuil(fem,sol);
893 mgmai=extract_skin_maille_entiere(fem,gest2);
894 }
895 if (decoup==1)
896 {
897 affiche((char*)" Découpage par isodensité");
898 double seuil=params.get_valeur("seuil");
899 std::string nomfich=params.get_nom("nomdensextra");
900 mgmai=extract_skin_par_decoupage(sol,seuil,gest2,nomfich.c_str());
901 }
902 suppression_peaux_isoles(mgmai);
903 affiche((char*)"Calcul du volume après extraction de la peau");
904 MG_MAILLAGE_OUTILS ot;
905 double v=ot.get_volume(mgmai);
906 fprintf(in,"Statistique après extraction de la peau\n");
907 fprintf(in,"Volume total=%.14le\nVolume non design=%.14le\nVolume design=%.14le\nf reel=%.14le\n\n\n\n",v,volnon,v-volnon,(v-volnon)/(voltot-volnon));
908 affiche((char*)"Procedure de lissage");
909
910 int liss_aretes2017_nbiter=params.get_valeur("liss_aretes2017_nbiter");
911 if(liss_aretes2017_nbiter!=0)
912 {
913 double lambda=params.get_valeur("liss_aretes2017_lambda");
914 lissage_aretes2017(mgmai,gest2,lambda,liss_aretes2017_nbiter);
915 }
916 else
917 {
918 LISTE_MG_NOEUD::iterator it_vec;
919 for(MG_NOEUD* no=mgmai->get_premier_noeud(it_vec);no!=NULL;no=mgmai->get_suivant_noeud(it_vec))
920 if((no->get_origine()==MAGIC::ORIGINE::TRIANGULATION_ARETE)||(no->get_origine()==MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE)) no->change_origine(MAGIC::ORIGINE::TRIANGULATION);
921 }
922
923 int chen2005_debut=params.get_valeur("chen2005_debut");
924 int chen2008_debut=params.get_valeur("chen2008_debut");
925 int jiao2012_debut=params.get_valeur("jiao2012_debut");
926 int McKenzie2016_debut=params.get_valeur("McKenzie2016_debut");
927 int Taubin1995_debut=params.get_valeur("Taubin1995_debut");
928 int cardinalite_debut=params.get_valeur("cardinalite_debut");
929 int cardinalite_iter=params.get_valeur("cardinalite_iter");
930 int rmimpose_debut=params.get_valeur("rmimpose_debut");
931 int rmimpose_iter=params.get_valeur("rmimpose_iter");
932 int chen2005fait=0;
933 if (chen2005_debut==0) chen2005fait=1;
934 int chen2008fait=0;
935 if (chen2008_debut==0) chen2008fait=1;
936 int jiao2012fait=0;
937 if (jiao2012_debut==0) jiao2012fait=1;
938 int McKenzie2016fait=0;
939 if (McKenzie2016_debut==0) McKenzie2016fait=1;
940 int Taubin1995fait=0;
941 if (Taubin1995_debut==0) Taubin1995fait=1;
942 int cardinalitefait=0;
943 if (cardinalite_debut==0) cardinalitefait=1;
944
945 int fin=0;
946 int iteration=1;
947 do
948 {
949 char nomiteration[500];
950 strcpy(nomiteration,"aucun lissage");
951 int next=2000000000;
952 if (chen2005fait==0) next=std::min(next,chen2005_debut);
953 if (chen2008fait==0) next=std::min(next,chen2008_debut);
954 if (jiao2012fait==0) next=std::min(next,jiao2012_debut);
955 if (McKenzie2016fait==0) next=std::min(next,McKenzie2016_debut);
956 if (Taubin1995fait==0) next=std::min(next,Taubin1995_debut);
957 if (cardinalitefait==0) next=std::min(next,cardinalite_debut);
958 if ((iteration>=chen2005_debut) && (chen2005_debut==next))
959 {
960 char message[300];
961 sprintf(message, " Iteration %d : lissage de Chen 2005",iteration);
962 affiche(message);
963 double epsilon=params.get_valeur("chen2005_epsilon");
964 double sigma=params.get_valeur("chen2005_sigma");
965 int iter_max=params.get_valeur("chen2005_itermax");
966 int nbiter=lissage_chen2005(mgmai,gest2,epsilon,sigma,iter_max);
967 if (nbiter==iter_max) sprintf(message," Arrêt de la procédure de lissage de Chen2005 après %d itérations",nbiter);
968 else sprintf(message," Convergence de la procédure de lissage de Chen2005 après %d itérations",nbiter);
969 affiche(message);
970 chen2005fait=1;
971 iteration=iteration+nbiter;
972 strcpy(nomiteration,"lissage de Chen2005");
973 }
974 else if ((iteration>=chen2008_debut) && (chen2008_debut==next))
975 {
976 char message[300];
977 sprintf(message, " Iteration %d : lissage de Chen 2008",iteration);
978 affiche(message);
979 double epsilon=params.get_valeur("chen2008_epsilon");
980 double sigma=params.get_valeur("chen2008_sigma");
981 double gamma=params.get_valeur("chen2008_gamma");
982 int iter_max=params.get_valeur("chen2008_itermax");
983 int nbiter=lissage_chen2008(mgmai,gest2,sigma,gamma,epsilon,iter_max);
984 if (nbiter==iter_max) sprintf(message," Arrêt de la procédure de lissage de Chen2008 après %d itérations",nbiter);
985 else sprintf(message," Convergence de la procédure de lissage de Chen2008 après %d itérations",nbiter);
986 affiche(message);
987 chen2008fait=1;
988 iteration=iteration+nbiter;
989 strcpy(nomiteration,"lissage de Chen2008");
990 }
991 else if ((iteration>=jiao2012_debut) && (jiao2012_debut==next))
992 {
993 char message[300];
994 sprintf(message, " Iteration %d : lissage de Jiao 2012",iteration);
995 affiche(message);
996 int iter_max=params.get_valeur("jiao2012_itermax");
997 int nbiter=lissage_jiao2012(mgmai,gest2,iter_max);
998 if (nbiter==iter_max) sprintf(message," Arrêt de la procédure de lissage de Jiao2012 après %d itérations",nbiter);
999 else sprintf(message," Convergence de la procédure de lissage de Jiao2012 après %d itérations",nbiter);
1000 affiche(message);
1001 jiao2012fait=1;
1002 iteration=iteration+nbiter;
1003 strcpy(nomiteration,"lissage de Jiao2012");
1004 }
1005 else if ((iteration>=McKenzie2016_debut) && (McKenzie2016_debut==next))
1006 {
1007 char message[300];
1008 sprintf(message, "Lissage de McKenzie2016");
1009 affiche(message);
1010 double lambda=params.get_valeur("McKenzie2016_lambda");
1011 double nu=params.get_valeur("McKenzie2016_nu");
1012 double epsilon=params.get_valeur("McKenzie2016_epsilon");
1013 double sigma=params.get_valeur("McKenzie2016_sigma");
1014 double gamma_ =params.get_valeur("McKenzie2016_sigma");
1015 int filtre=params.get_valeur("McKenzie2016_filtre");
1016 int iter_max=params.get_valeur("McKenzie2016_itermax");
1017 int itr_taubin=params.get_valeur("McKenzie2016_itertaubmax");
1018 int nbiter=lissage_McKenzie2016(mgmai,gest2,lambda,nu,epsilon,sigma,gamma_,filtre, iter_max, itr_taubin);
1019 if (nbiter==iter_max) sprintf(message," Arrêt de la procédure de lissage de McKenzie2016 après %d itérations",nbiter);
1020 else sprintf(message," Convergence de la procédure de lissage de McKenzie2016 après %d itérations",nbiter);
1021 affiche(message);
1022 McKenzie2016fait=1;
1023 iteration=iteration+nbiter;
1024 strcpy(nomiteration,"lissage de McKenzie2016");
1025 }
1026
1027 else if ((iteration>=Taubin1995_debut) && (Taubin1995_debut==next))
1028 {
1029 char message[300];
1030 sprintf(message, "Lissage de Taubin1995");
1031 affiche(message);
1032 double lambda=params.get_valeur("Taubin1995_lambda");
1033 double nu=params.get_valeur("Taubin1995_nu");
1034 int iter_max=params.get_valeur("Taubin1995_itermax");
1035 int nbiter=lissage_Taubin1995(mgmai,gest2,lambda,nu, iter_max);
1036 if (nbiter==iter_max) sprintf(message," Arrêt de la procédure de lissage de Taubin1995 après %d itérations",nbiter);
1037 else sprintf(message," Convergence de la procédure de lissage de Taubin1995 après %d itérations",nbiter);
1038 affiche(message);
1039 Taubin1995fait=1;
1040 iteration=iteration+nbiter;
1041 strcpy(nomiteration,"lissage de Taubin1995");
1042 }
1043 else if ((iteration>=cardinalite_debut) && (cardinalite_debut==next))
1044 {
1045 char message[300];
1046 sprintf(message, "Étude de la cardinalite des noeuds");
1047 affiche(message);
1048 int nb=cardinalite(mgmai,cardinalite_iter);
1049 if (nb!=cardinalite_iter)
1050 sprintf(message, "Cardinalité des noeuds convergé");
1051 else
1052 sprintf(message, "Cardinalité des noeuds non convergé");
1053 affiche(message);
1054 strcpy(nomiteration,"Cardinalité");
1055 iteration=iteration+nb;
1056 cardinalitefait=1;
1057 }
1058 else
1059 {
1060 iteration++;
1061 }
1062 if (chen2005fait==1)
1063 if (chen2008fait==1)
1064 if (jiao2012fait==1)
1065 if (McKenzie2016fait==1)
1066 if (Taubin1995fait==1)
1067 if (cardinalitefait==1)
1068 fin=1;
1069 if (strcmp(nomiteration,"aucun lissage")!=0)
1070 {
1071 MG_MAILLAGE_OUTILS ot;
1072 double v=ot.get_volume(mgmai);
1073 fprintf(in,"Statistique après itération %d %s\n",iteration-1,nomiteration);
1074 fprintf(in,"Volume total=%.14le\nVolume non design=%.14le\nVolume design=%.14le\nf reel=%.14le\n\n\n\n",v,volnon,v-volnon,(v-volnon)/(voltot-volnon));
1075 }
1076
1077 }
1078 while (fin==0);
1079
1080 fclose(in);
1081
1082 if (rmimpose_debut != 0)
1083 {
1084 char message[300];
1085 sprintf(message, "Suppression du non_design du resultat optimal");
1086 affiche(message);
1087 std::string nomfichier=params.get_nom("nomfichpeau");
1088 rmimpose(mgmai,(char*)nomfichier.c_str());
1089 }
1090 }
1091
1092
1093 MG_MAILLAGE* MGOPT_POSTTRAITEMENT::extract_skin_maille_entiere(FEM_MAILLAGE* mai,MG_GESTIONNAIRE &gest2)
1094 {
1095 //Menu de la methode de lissage
1096 affiche((char*)"Extraction du maillage de surface");
1097 int coderesu = 0;
1098 int mai2_id;
1099 int imp=params.get_valeur("consimpose");
1100 int opti=params.get_valeur("consoptimise");
1101 int m_auto=params.get_valeur("consnooptimise");
1102 if (opti==1) conserve(MAGIC::ORIGINE::OPTIMISE);
1103 if (imp==1) conserve(MAGIC::ORIGINE::IMPOSE);
1104 if (m_auto==1) conserve(MAGIC::ORIGINE::MAILLEUR_AUTO);
1105 copieorigine(mai);
1106 MG_MAILLAGE* mg_mai = (MG_MAILLAGE*)mai->get_mg_maillage();
1107 //gain_poids(mg_mai,gest2);
1108 int reactiv=params.get_valeur("reactivation");
1109 if (reactiv == 1)
1110 {
1111 reactivation(mg_mai,gest2);
1112 }
1113 affiche((char*)" Analyse des cas non manifold");
1114 do
1115 {
1116 int nbmaniare,nbmanino,nbpeau;
1117 coderesu = extract_skin(mg_mai,gest2,nbpeau,nbmaniare,nbmanino,&mai2_id);
1118 char message[500];
1119 sprintf(message," %d peaux. %d manifold par arete. %d manifold par noeud",nbpeau,nbmaniare,nbmanino);
1120 affiche(message);
1121 // gain_poids(mg_mai,gest2);
1122 }
1123 while (coderesu == 0);
1124
1125 MG_MAILLAGE* mg_mai2=gest2.get_mg_maillageid(mai2_id);
1126
1127 return mg_mai2;
1128
1129 }
1130
1131
1132 void MGOPT_POSTTRAITEMENT::adapte_seuil(class FEM_MAILLAGE* fem,FEM_SOLUTION* solution)
1133 {
1134 double seuil=params.get_valeur("seuil");
1135 solution->active_solution(0);
1136 LISTE_FEM_ELEMENT3::iterator it;
1137 for (FEM_ELEMENT3 *tet=fem->get_premier_element3(it);tet!=NULL;tet=fem->get_suivant_element3(it))
1138 {
1139 if (((MG_TETRA*)tet->get_mg_element_maillage())->get_origine()!=MAGIC::ORIGINE::IMPOSE)
1140 if (tet->get_solution()>seuil)
1141 ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(MAGIC::ORIGINE::OPTIMISE);
1142 else
1143 ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(MAGIC::ORIGINE::MAILLEUR_AUTO);
1144
1145 }
1146 }
1147
1148
1149
1150 //---------------------------- Lissage des arêtes 2017 Abdennour -------------------------------------------------------------------------
1151 void MGOPT_POSTTRAITEMENT::lissage_aretes2017(MG_MAILLAGE* mg_mai, MG_GESTIONNAIRE& gest2, double lambda, int nb_iter)
1152 {
1153 //Etape 1 - Recherche de tous les noeuds des arêtes
1154 std::vector<MG_NOEUD*> lst_no_1090_1091;
1155 LISTE_MG_NOEUD::iterator it_vec;
1156 MG_NOEUD* no=NULL;
1157 for(no=mg_mai->get_premier_noeud(it_vec);no!=NULL;no=mg_mai->get_suivant_noeud(it_vec))
1158 if((no->get_origine()==MAGIC::ORIGINE::TRIANGULATION_ARETE)||(no->get_origine()==MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE)) lst_no_1090_1091.push_back(no);
1159
1160 //Etape 2 - Déterminantion des noeuds extrêmes de chaque arête ouverte
1161 std::vector<MG_NOEUD*> lst_no_extr;
1162 MG_NOEUD* no_voisin=NULL;
1163 MG_SEGMENT* seg=NULL;
1164 for (int i=0;i<lst_no_1090_1091.size();i++)
1165 {
1166 no=lst_no_1090_1091[i];
1167 for(int j=0;j<no->get_lien_segment()->get_nb();j++)
1168 {
1169 seg=no->get_lien_segment()->get(j);
1170 bool seg_bord=seg->get_lien_triangle()->get(0)->get_lien_topologie()!=seg->get_lien_triangle()->get(1)->get_lien_topologie();
1171 if(seg->get_noeud1()!=no) no_voisin=seg->get_noeud1();
1172 else no_voisin=seg->get_noeud2();
1173 if((no_voisin->get_origine()==MAGIC::ORIGINE::IMPOSE)&&seg_bord)
1174 lst_no_extr.push_back(no_voisin);
1175 }
1176 }
1177
1178 //Etape 3 - Déterminantion des arêtes ouvertes
1179 std::vector<std::vector<MG_NOEUD*>> lst_aretes;
1180 std::map<int,MG_NOEUD*> lst_no_voisins;
1181 TPL_LISTE_ENTITE<long unsigned int> lst_idno_nonmanifold;
1182 MG_GEOMETRIE* geo=gest2.get_mg_geometrie(0);
1183 void recherche_voisins(MG_NOEUD* no, std::vector<MG_NOEUD*>* lst_no_1090_1091,std::vector<MG_NOEUD*>* arete, std::map<int,MG_NOEUD*>* lst_no_voisins, bool arete_ouverte, std::vector<MG_NOEUD*>* lst_no_extr);
1184 int passage_non_manifold(MG_NOEUD** no,std::map<int,MG_NOEUD*>* lst_no_voisins, MG_GEOMETRIE* geo, std::vector<MG_NOEUD*>* arete);
1185 bool arete_ouverte=true;
1186
1187 for (int i=0;i<lst_no_extr.size();i++)
1188 {
1189 std::vector<MG_NOEUD*> arete;
1190 no=lst_no_extr[i];
1191 arete.push_back(no);
1192 do
1193 {
1194 lst_no_voisins.clear();
1195 recherche_voisins(no,&lst_no_1090_1091,&arete,&lst_no_voisins,arete_ouverte,&lst_no_extr);
1196 if((lst_no_voisins.size()==1)||(lst_no_voisins.size()==2))
1197 {
1198 arete.push_back(lst_no_voisins.begin()->second);
1199 lst_no_1090_1091.erase(lst_no_1090_1091.begin()+lst_no_voisins.begin()->first);
1200 no=lst_no_voisins.begin()->second;
1201 }
1202 /***traitement des cas "non-manifold"***/
1203 if(lst_no_voisins.size()==3)
1204 {
1205 lst_idno_nonmanifold.ajouter(no->get_id());
1206 lst_no_1090_1091.push_back(no);
1207 lst_no_1090_1091.erase(lst_no_1090_1091.begin()+passage_non_manifold(&no,&lst_no_voisins,geo,&arete));
1208 }
1209 }
1210 while(lst_no_voisins.size()!=0);
1211 if(arete.size()>6) lst_aretes.push_back(arete);
1212 else for(int i=1;i<arete.size()-1;i++)
1213 if((arete[i]->get_origine()!=MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE)&&(lst_idno_nonmanifold.est_dans_la_liste(arete[i]->get_id())==0))
1214 arete[i]->change_origine(MAGIC::ORIGINE::TRIANGULATION);
1215 }
1216
1217 //Etape 4 - Déterminantion des arêtes fermées
1218 arete_ouverte=false;
1219 for (int i=0;i<lst_no_1090_1091.size();i++)
1220 if(lst_no_1090_1091[i]->get_origine()!=MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE)
1221 {
1222 std::vector<MG_NOEUD*> arete;
1223 no=lst_no_1090_1091[i];
1224 arete.push_back(no);
1225 do
1226 {
1227 lst_no_voisins.clear();
1228 recherche_voisins(no,&lst_no_1090_1091,&arete,&lst_no_voisins,arete_ouverte,&lst_no_extr);
1229 if((lst_no_voisins.size()==1)||(lst_no_voisins.size()==2))
1230 {
1231 arete.push_back(lst_no_voisins.begin()->second);
1232 lst_no_1090_1091.erase(lst_no_1090_1091.begin()+lst_no_voisins.begin()->first);
1233 no=lst_no_voisins.begin()->second;
1234 }
1235 /***traitement des cas "non-manifold"***/
1236 if(lst_no_voisins.size()==3)
1237 {
1238 lst_idno_nonmanifold.ajouter(no->get_id());
1239 lst_no_1090_1091.push_back(no);
1240 lst_no_1090_1091.erase(lst_no_1090_1091.begin()+passage_non_manifold(&no,&lst_no_voisins,geo,&arete));
1241 }
1242 }
1243 while(lst_no_voisins.size()!=0);
1244 if(arete.size()>6) lst_aretes.push_back(arete);
1245 else
1246 for(int i=0;i<arete.size()-1;i++)
1247 if((arete[i]->get_origine()!=MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE)&&(lst_idno_nonmanifold.est_dans_la_liste(arete[i]->get_id())==0))
1248 arete[i]->change_origine(MAGIC::ORIGINE::TRIANGULATION);
1249 }
1250
1251 //Etape 5 - Lissage des arêtes
1252 for(int iter=0;iter<nb_iter;iter++)
1253 for(int k=0;k<lst_aretes.size();k++)
1254 for(int i=0;i<lst_aretes[k].size()-1;i++)
1255 if((lst_aretes[k][i]->get_origine()==MAGIC::ORIGINE::TRIANGULATION_ARETE)&&
1256 (lst_idno_nonmanifold.est_dans_la_liste(lst_aretes[k][i]->get_id())==0))
1257 {
1258 MG_NOEUD* noeud_i=lst_aretes[k][i];
1259
1260 int i_prec=i-1;
1261 if(i==0) i_prec=lst_aretes[k].size()-2;
1262
1263 double x_i=lst_aretes[k][i]->get_x(), x_i_prec=lst_aretes[k][i_prec]->get_x(), x_i_suiv=lst_aretes[k][i+1]->get_x();
1264 double y_i=lst_aretes[k][i]->get_y(), y_i_prec=lst_aretes[k][i_prec]->get_y(), y_i_suiv=lst_aretes[k][i+1]->get_y();
1265 double z_i=lst_aretes[k][i]->get_z(), z_i_prec=lst_aretes[k][i_prec]->get_z(), z_i_suiv=lst_aretes[k][i+1]->get_z();
1266
1267 double delta_x=0.5*(x_i_prec+x_i_suiv)-x_i;
1268 double delta_y=0.5*(y_i_prec+y_i_suiv)-y_i;
1269 double delta_z=0.5*(z_i_prec+z_i_suiv)-z_i;
1270
1271 double x_i_prime=x_i+lambda*delta_x;
1272 double y_i_prime=y_i+lambda*delta_y;
1273 double z_i_prime=z_i+lambda*delta_z;
1274
1275 noeud_i->change_x(x_i_prime);
1276 noeud_i->change_y(y_i_prime);
1277 noeud_i->change_z(z_i_prime);
1278 }
1279 }
1280
1281
1282
1283 //---------------------------- Lissage Chen 2005 Gilles Philippe -------------------------------------------------------------------------
1284 int MGOPT_POSTTRAITEMENT::lissage_chen2005(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double epsilon, double sigma, int iter_max)
1285 {
1286 //Methode de lissage inspiree de Chen(2005) par Gilles Philipe Picher Martel
1287 double un_sur_pi = 1./M_PI;
1288 int compteur = 0;
1289 int fin = 0;
1290
1291
1292 do
1293 {
1294 TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
1295 TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales2;
1296 TPL_LISTE_ENTITE<double> liste_wij;
1297 LISTE_MG_TRIANGLE::iterator it_tri;
1298 int k = 0; //pour identifier les triangles pour liste_normales et liste_wij
1299 for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
1300 {
1301 MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
1302 OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
1303 normal_f_i.norme();
1304 liste_normales2.ajouter(normal_f_i);
1305 //Remplissage de la liste des voisins du triangle i
1306 TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
1307 MG_NOEUD* noeud1 = mgtri_i->get_noeud1();
1308 double nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
1309 for (int j = 0;j<nb_voisins1;j++)
1310 {
1311 MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
1312 liste_voisins.ajouter(mgtri_1);
1313 }
1314 MG_NOEUD* noeud2 = mgtri_i->get_noeud2();
1315 double nb_voisins2 = noeud2->get_lien_triangle()->get_nb();
1316 for (int j = 0;j<nb_voisins2;j++)
1317 {
1318 MG_TRIANGLE_PEAU* mgtri_2 = (MG_TRIANGLE_PEAU*) noeud2->get_lien_triangle()->get(j);
1319 liste_voisins.ajouter(mgtri_2);
1320 }
1321 MG_NOEUD* noeud3 = mgtri_i->get_noeud3();
1322 double nb_voisins3 = noeud3->get_lien_triangle()->get_nb();
1323 for (int j = 0;j<nb_voisins3;j++)
1324 {
1325 MG_TRIANGLE_PEAU* mgtri_3 = (MG_TRIANGLE_PEAU*) noeud3->get_lien_triangle()->get(j);
1326 liste_voisins.ajouter(mgtri_3);
1327 }
1328 liste_voisins.supprimer(mgtri_i);
1329 int nb_voisins = liste_voisins.get_nb();
1330 double w_ij = 1./nb_voisins;
1331 double phi_i_min = 10.;
1332 double s_i = 0.0;
1333 double phi_im = 0.0;
1334 double *phi_ij = new double[nb_voisins];
1335 OT_VECTEUR_3D normal_f_i_mean(0.,0.,0.);
1336 OT_VECTEUR_3D eta_i(0.,0.,0.);
1337 TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
1338 int j = 0;
1339 for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
1340 {
1341 OT_VECTEUR_3D normal_f_j = mgtri_j->calcul_normal();
1342 //1-Calculer la normale moyenne pour chaque triangle
1343 normal_f_i_mean = normal_f_i_mean + w_ij*normal_f_j;
1344 //2.1-On calcul l'angle entre normal_f_i et normal_f_j pour j allant de 1 a Nb_voisins
1345 double prod_scalaire = normal_f_i*normal_f_j;
1346 if (prod_scalaire > 1.)
1347 {
1348 prod_scalaire = 1.;
1349 }
1350 if (prod_scalaire < -1.)
1351 {
1352 prod_scalaire = -1.;
1353 }
1354 phi_ij[j] = acos(prod_scalaire)*un_sur_pi;
1355 //2.2-On trouve le plus petit des angles et la normale heta_i correspondante
1356 if (phi_ij[j] < phi_i_min)
1357 {
1358 phi_i_min = phi_ij[j];
1359 eta_i = normal_f_j;
1360 }
1361 //3.1-On calcul l'angle moyen phi_im
1362 phi_im = phi_im + w_ij*phi_ij[j];
1363 j++;
1364 }
1365 normal_f_i_mean.norme();
1366 j = 0;
1367 for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
1368 {
1369 //3.2-Calcul de s_i selon la variance
1370 s_i = s_i + w_ij*pow((phi_ij[j] - phi_im),2);
1371 j++;
1372 }
1373 delete[] phi_ij;
1374 //4-On calcule une nouvelle normale pour chaque triangle
1375 double pond;
1376 int num=params.get_valeur("chen2005_filtre");
1377 if (num==1) pond=ponderation_gaussian(s_i,sigma);
1378 else if (num==2) pond=ponderation_laplacian(s_i,sigma);
1379 else if (num==3) pond=ponderation_elfallahford(s_i,sigma);
1380 OT_VECTEUR_3D normal_f_i_new = pond*normal_f_i_mean + (1. - pond)*eta_i;
1381 normal_f_i_new.norme();
1382 liste_normales.ajouter(normal_f_i_new);
1383 liste_wij.ajouter(w_ij);
1384 mgtri->change_nouveau_numero(k);
1385 k++;
1386 }
1387
1388 LISTE_MG_NOEUD::iterator it_no;
1389 for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1390 {
1391 int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
1392 double w_ij_prime = 0.0;
1393 OT_VECTEUR_3D v_temp(0.,0.,0.);
1394 OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
1395 for(int j=0;j<nb_voisins_j;j++)
1396 {
1397 MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
1398 //On calcule le centroide cj du triangle mgtri_j
1399 MG_NOEUD* n1 = mgtri_j->get_noeud1();
1400 MG_NOEUD* n2 = mgtri_j->get_noeud2();
1401 MG_NOEUD* n3 = mgtri_j->get_noeud3();
1402 double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
1403 double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
1404 double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
1405 //On forme le vecteur vi_cj
1406 OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
1407 OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_j->get_nouveau_numero());
1408 // w_ij_prime correspond a la somme des aires des triangles voisins du noeuds
1409 OT_VECTEUR_3D AB(n2->get_x() - n1->get_x(),n2->get_y() - n1->get_y(),n2->get_z() - n1->get_z());
1410 OT_VECTEUR_3D AC(n3->get_x() - n1->get_x(),n3->get_y() - n1->get_y(),n3->get_z() - n1->get_z());
1411 OT_VECTEUR_3D prodvect = AB&AC;
1412 double w_ij = 0.5*prodvect.get_longueur();
1413 w_ij_prime = w_ij_prime + w_ij;
1414 v_temp = v_temp + w_ij*(vi_cj*normal_f_i_new)*normal_f_i_new;
1415 }
1416 //5-On met a jour la position des noeuds
1417 v_i = v_i + v_temp/w_ij_prime;
1418 int origine = noeud_i->get_origine();
1419 if (origine == MAGIC::ORIGINE::TRIANGULATION)
1420 {
1421 noeud_i->change_x(v_i.get_x());
1422 noeud_i->change_y(v_i.get_y());
1423 noeud_i->change_z(v_i.get_z());
1424 }
1425 }
1426 //Critere d'arret de l'algorithme
1427 int l=0;
1428 int nb_tri = mg_mai->get_nb_mg_triangle();
1429 for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
1430 {
1431 MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
1432 OT_VECTEUR_3D normal_f_i = liste_normales2.get(mgtri_i->get_nouveau_numero());
1433 OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_i->get_nouveau_numero());
1434 double critere = 1. - normal_f_i*normal_f_i_new;
1435 if (critere <= epsilon) l++;
1436 }
1437 double tolerance = 0.01*nb_tri;
1438 if (nb_tri - l <= 0) fin = 1;
1439 compteur++;
1440 }
1441 while ((fin == 0) && (compteur < iter_max));
1442
1443 return compteur;
1444 }
1445
1446 //---------------------------- Lissage Jiao 2012 -------------------------------------------------------------------------
1447 int MGOPT_POSTTRAITEMENT::lissage_jiao2012(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, int iter_max)
1448 {
1449 // Lissage global avec methode de lissage inspiree de JIAO IMR2012 par Jean Christophe sept 2013
1450 int compteur = 0;
1451
1452 do
1453 {
1454 std::vector<double> nouv_position_x;
1455 std::vector<double> nouv_position_y;
1456 std::vector<double> nouv_position_z;
1457 LISTE_MG_NOEUD::iterator it_no;
1458 for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1459 {
1460 int origine = noeud_i->get_origine();
1461 if (origine == MAGIC::ORIGINE::TRIANGULATION)
1462 {
1463 int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
1464 double wij = 0.0;
1465 double wij_sommej = 0.0;
1466 OT_VECTEUR_3D v_temp(0.,0.,0.);
1467 OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
1468 for(int j=0;j<nb_voisins_j;j++)
1469 {
1470 MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
1471 //On calcule le centroide cj du triangle mgtri_j
1472 MG_NOEUD* n1 = mgtri_j->get_noeud1();
1473 MG_NOEUD* n2 = mgtri_j->get_noeud2();
1474 MG_NOEUD* n3 = mgtri_j->get_noeud3();
1475 double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
1476 double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
1477 double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
1478 //On forme le vecteur vi_cj et le vecteur cj
1479 OT_VECTEUR_3D cj(cj_x ,cj_y,cj_z);
1480 OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
1481 wij=vi_cj.get_longueur();
1482 wij_sommej=wij_sommej+wij;
1483 v_temp = v_temp + (wij*cj);
1484 }
1485 //5-On calcule la nouvelle position des noeuds et on affecte au tableau temporaire
1486 v_i =v_temp/wij_sommej;
1487 nouv_position_x.push_back(v_i.get_x());
1488 nouv_position_y.push_back(v_i.get_y());
1489 nouv_position_z.push_back(v_i.get_z());
1490 }
1491 }
1492 //On actualise la position des noeuds
1493 int ind_noeud=0;
1494 for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1495 {
1496 int origine = noeud_i->get_origine();
1497 if (origine == MAGIC::ORIGINE::TRIANGULATION)
1498 {
1499 noeud_i->change_x(nouv_position_x[ind_noeud]);
1500 noeud_i->change_y(nouv_position_y[ind_noeud]);
1501 noeud_i->change_z(nouv_position_z[ind_noeud]);
1502 ind_noeud++;
1503 }
1504 }
1505 //Critere d'arret de l'algorithme
1506 compteur++;
1507 }
1508 while (compteur < iter_max);
1509
1510 return compteur;
1511 }
1512
1513 //---------------------------- Lissage Chen 2008 -------------------------------------------------------------------------
1514 int MGOPT_POSTTRAITEMENT::lissage_chen2008(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double sigma, double gamma, double epsilon, int iter_max)
1515 {
1516 //Lissage global inspire de Chen(2008) corrige par Jean Christophe sept 2013
1517 double un_sur_pi = 1./M_PI;
1518 int compteur = 0;
1519 int fin = 0;
1520 compteur = 0;
1521
1522 do
1523 {
1524 std::vector<double> nouv_position_x;
1525 std::vector<double> nouv_position_y;
1526 std::vector<double> nouv_position_z;
1527 TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
1528 TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales2;
1529 TPL_LISTE_ENTITE<double> liste_wij;
1530 LISTE_MG_TRIANGLE::iterator it_tri;
1531 int k = 0; //pour identifier les triangles pour liste_normales et liste_wij
1532 for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
1533 {
1534 MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
1535 OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
1536 normal_f_i.norme();
1537 liste_normales2.ajouter(normal_f_i);
1538 //Remplissage de la liste des voisins du triangle i
1539 TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
1540 MG_NOEUD* noeud1 = mgtri_i->get_noeud1();
1541 double nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
1542 for (int j = 0;j<nb_voisins1;j++)
1543 {
1544 MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
1545 liste_voisins.ajouter(mgtri_1);
1546 }
1547 MG_NOEUD* noeud2 = mgtri_i->get_noeud2();
1548 double nb_voisins2 = noeud2->get_lien_triangle()->get_nb();
1549 for (int j = 0;j<nb_voisins2;j++)
1550 {
1551 MG_TRIANGLE_PEAU* mgtri_2 = (MG_TRIANGLE_PEAU*) noeud2->get_lien_triangle()->get(j);
1552 liste_voisins.ajouter(mgtri_2);
1553 }
1554 MG_NOEUD* noeud3 = mgtri_i->get_noeud3();
1555 double nb_voisins3 = noeud3->get_lien_triangle()->get_nb();
1556 for (int j = 0;j<nb_voisins3;j++)
1557 {
1558 MG_TRIANGLE_PEAU* mgtri_3 = (MG_TRIANGLE_PEAU*) noeud3->get_lien_triangle()->get(j);
1559 liste_voisins.ajouter(mgtri_3);
1560 }
1561 liste_voisins.supprimer(mgtri_i);
1562 int nb_voisins = liste_voisins.get_nb();
1563 double w_ij = 1./nb_voisins;
1564 double phi_i_min = 10.;
1565 double s_i = 0.0;
1566 double phi_im = 0.0;
1567 double *phi_ij = new double[nb_voisins];
1568 OT_VECTEUR_3D normal_f_i_mean(0.,0.,0.);
1569 normal_f_i_mean = normal_f_i;
1570 OT_VECTEUR_3D eta_i(0.,0.,0.);
1571 TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
1572 int j = 0;
1573 for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
1574 {
1575 OT_VECTEUR_3D normal_f_j = mgtri_j->calcul_normal();
1576 //1-Calculer la normale moyenne pour chaque triangle
1577 normal_f_i_mean = normal_f_i_mean + normal_f_j;
1578 //2.1-On calcule l'angle entre normal_f_i et normal_f_j pour j allant de 1 a Nb_voisins
1579 double prod_scalaire = normal_f_i*normal_f_j;
1580 if (prod_scalaire > 1.)
1581 {
1582 prod_scalaire = 1.;
1583 }
1584 if (prod_scalaire < -1.)
1585 {
1586 prod_scalaire = -1.;
1587 }
1588 phi_ij[j] = acos(prod_scalaire)*un_sur_pi;
1589 //2.2-On trouve le plus petit des angles et la normale heta_i correspondante
1590 if (phi_ij[j] < phi_i_min)
1591 {
1592 phi_i_min = phi_ij[j];
1593 eta_i = normal_f_j;
1594 }
1595 //3.1-On calcule l'angle moyen phi_im
1596 phi_im = phi_im + w_ij*phi_ij[j];
1597 j++;
1598 }
1599 normal_f_i_mean.norme();
1600 j = 0;
1601 for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
1602 {
1603 //3.2-Calcul de s_i selon la variance
1604 s_i = s_i + w_ij*pow((phi_ij[j] - phi_im),2);
1605 j++;
1606 }
1607 delete[] phi_ij;
1608 //4-On calcule une nouvelle normale pour chaque triangle
1609 double pond;
1610 int num=params.get_valeur("chen2008_filtre");
1611 if (num==1) pond=ponderation_gaussian(s_i,sigma);
1612 else if (num==2) pond=ponderation_laplacian(s_i,sigma);
1613 else if (num==3) pond=ponderation_elfallahford(s_i,sigma);
1614 OT_VECTEUR_3D normal_f_i_new = pond*normal_f_i_mean + (1. - pond)*eta_i;
1615 normal_f_i_new.norme();
1616 liste_normales.ajouter(normal_f_i_new);
1617 liste_wij.ajouter(w_ij);
1618 mgtri->change_nouveau_numero(k);
1619 k++;
1620 }
1621
1622 LISTE_MG_NOEUD::iterator it_no;
1623 for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1624 {
1625 int origine = noeud_i->get_origine();
1626 if (origine == MAGIC::ORIGINE::TRIANGULATION)
1627 {
1628 int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
1629 double w_ij_prime = 0.0;
1630 OT_VECTEUR_3D v_temp(0.,0.,0.);
1631 OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
1632 for(int j=0;j<nb_voisins_j;j++)
1633 {
1634 MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
1635 //On calcule le centroide cj du triangle mgtri_j
1636 MG_NOEUD* n1 = mgtri_j->get_noeud1();
1637 MG_NOEUD* n2 = mgtri_j->get_noeud2();
1638 MG_NOEUD* n3 = mgtri_j->get_noeud3();
1639 double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
1640 double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
1641 double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
1642 //On forme le vecteur vi_cj
1643 OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
1644 OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_j->get_nouveau_numero());
1645 v_temp = v_temp + (vi_cj*normal_f_i_new)*normal_f_i_new;
1646 }
1647 //5-On met a jour la position des noeuds
1648 v_i = v_i + (gamma/(2*nb_voisins_j))*v_temp;
1649 nouv_position_x.push_back(v_i.get_x());
1650 nouv_position_y.push_back(v_i.get_y());
1651 nouv_position_z.push_back(v_i.get_z());
1652 }
1653 }
1654
1655 //On actualise la position des noeuds
1656 int ind_noeud=0;
1657 for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1658 {
1659 int origine = noeud_i->get_origine();
1660 if (origine == MAGIC::ORIGINE::TRIANGULATION)
1661 {
1662 noeud_i->change_x(nouv_position_x[ind_noeud]);
1663 noeud_i->change_y(nouv_position_y[ind_noeud]);
1664 noeud_i->change_z(nouv_position_z[ind_noeud]);
1665 ind_noeud++;
1666 }
1667 }
1668
1669 //Critere d'arret de l'algorithme
1670 int l=0;
1671 int nb_tri = mg_mai->get_nb_mg_triangle();
1672 for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
1673 {
1674 MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
1675 OT_VECTEUR_3D normal_f_i = liste_normales2.get(mgtri_i->get_nouveau_numero());
1676 OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_i->get_nouveau_numero());
1677 double critere = 1. - normal_f_i*normal_f_i_new;
1678 if (critere <= epsilon) l++;
1679 }
1680 double tolerance = 0.01*nb_tri;
1681 if (nb_tri - l <= 0) fin = 1;
1682 compteur++;
1683 }
1684 while ((fin == 0) && (compteur < iter_max));
1685
1686 return compteur;
1687
1688 }
1689
1690
1691
1692 //---------------------------- Cycle Taubin 1995 par Bryan McKnezie 2016-------------------------------------------------------------------------
1693 int MGOPT_POSTTRAITEMENT::cycle_taubin1995(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double lambda, double nu)
1694 {
1695
1696 /////partie lambda du lissage
1697 std::vector<double> nouv_position_x;
1698 std::vector<double> nouv_position_y;
1699 std::vector<double> nouv_position_z;
1700 LISTE_MG_NOEUD::iterator it_no;
1701 for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1702 {
1703 int origine = noeud_i->get_origine();
1704 if ((origine == MAGIC::ORIGINE::TRIANGULATION)&& (noeud_i->get_lien_topologie()==NULL))//((origine != MAGIC::ORIGINE::IMPOSE)&& (origine != MAGIC::ORIGINE::TRIANGULATION_ARETE)&&(noeud_i->get_lien_topologie()==NULL))///////////////////////////////////////////////////////////////////////////////
1705 {
1706 TPL_MAP_ENTITE<MG_NOEUD*> liste_voisin;
1707 MG_MAILLAGE_OUTILS ot;
1708 ot.get_noeud_voisin(noeud_i,liste_voisin,2);
1709
1710 int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
1711 double delta_x = 0.0;
1712 double delta_y = 0.0;
1713 double delta_z = 0.0;
1714 double v_x = noeud_i->get_x();//v_1->get_x();
1715 double v_y = noeud_i->get_y();//v_1->get_y();
1716 double v_z = noeud_i->get_z();//v_1->get_z();
1717 double multiplicateur = (1.0/nb_voisins_j);
1718
1719 for(int j=0;j<nb_voisins_j;j++)
1720 {
1721 double v_i_x=liste_voisin.get(j)->get_x();
1722 double v_i_y=liste_voisin.get(j)->get_y();
1723 double v_i_z=liste_voisin.get(j)->get_z();
1724 delta_x = delta_x+(multiplicateur*(v_i_x-v_x));
1725 delta_y = delta_y+(multiplicateur*(v_i_y-v_y));
1726 delta_z = delta_z+(multiplicateur*(v_i_z-v_z));
1727
1728 }
1729 double x_prime = (v_x)+(lambda*delta_x);
1730 double y_prime = (v_y)+(lambda*delta_y);
1731 double z_prime = (v_z)+(lambda*delta_z);
1732
1733 //On calcule la nouvelle position des noeuds et on affecte au tableau temporaire
1734 nouv_position_x.push_back(x_prime);
1735 nouv_position_y.push_back(y_prime);
1736 nouv_position_z.push_back(z_prime);
1737 }
1738 }
1739 //On actualise la position des noeuds
1740 int ind_noeud=0;
1741 for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1742 {
1743 int origine = noeud_i->get_origine();
1744 if ((origine == MAGIC::ORIGINE::TRIANGULATION)&& (noeud_i->get_lien_topologie()==NULL))//((origine != MAGIC::ORIGINE::IMPOSE)&& (origine != MAGIC::ORIGINE::TRIANGULATION_ARETE)&&(noeud_i->get_lien_topologie()==NULL))///////////////////////////////////////////////////////////////////////////////
1745 {
1746 noeud_i->change_x(nouv_position_x[ind_noeud]);
1747 noeud_i->change_y(nouv_position_y[ind_noeud]);
1748 noeud_i->change_z(nouv_position_z[ind_noeud]);
1749 ind_noeud++;
1750 }
1751 }
1752
1753
1754 //partie nu du lissage
1755 nouv_position_x.clear();
1756 nouv_position_y.clear();
1757 nouv_position_z.clear();
1758 for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1759 {
1760 int origine = noeud_i->get_origine();
1761 if ((origine == MAGIC::ORIGINE::TRIANGULATION)&& (noeud_i->get_lien_topologie()==NULL))//((origine != MAGIC::ORIGINE::IMPOSE)&& (origine != MAGIC::ORIGINE::TRIANGULATION_ARETE)&&(noeud_i2->get_lien_topologie()==NULL))///////////////////////////////////////////////////////////////////////////////
1762 {
1763 TPL_MAP_ENTITE<MG_NOEUD*> liste_voisin;
1764 MG_MAILLAGE_OUTILS ot;
1765 ot.get_noeud_voisin(noeud_i,liste_voisin,2);
1766
1767 int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
1768 double delta_x = 0.0;
1769 double delta_y = 0.0;
1770 double delta_z = 0.0;
1771 double v_x = noeud_i->get_x();//v_1->get_x();
1772 double v_y = noeud_i->get_y();//v_1->get_y();
1773 double v_z = noeud_i->get_z();//v_1->get_z();
1774 double multiplicateur = (1.0/nb_voisins_j);
1775
1776 for(int j=0;j<nb_voisins_j;j++)
1777 {
1778 double v_i_x=liste_voisin.get(j)->get_x();
1779 double v_i_y=liste_voisin.get(j)->get_y();
1780 double v_i_z=liste_voisin.get(j)->get_z();
1781 delta_x = delta_x+(multiplicateur*(v_i_x-v_x));
1782 delta_y = delta_y+(multiplicateur*(v_i_y-v_y));
1783 delta_z = delta_z+(multiplicateur*(v_i_z-v_z));
1784
1785 }
1786 double x_prime = (v_x)+(nu*delta_x);
1787 double y_prime = (v_y)+(nu*delta_y);
1788 double z_prime = (v_z)+(nu*delta_z);
1789
1790 //On calcule la nouvelle position des noeuds et on affecte au tableau temporaire
1791 nouv_position_x.push_back(x_prime);
1792 nouv_position_y.push_back(y_prime);
1793 nouv_position_z.push_back(z_prime);
1794 }
1795 }
1796 //On actualise la position des noeuds
1797 ind_noeud=0;
1798 for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1799 {
1800 int origine = noeud_i->get_origine();
1801 if ((origine == MAGIC::ORIGINE::TRIANGULATION)&& ( noeud_i->get_lien_topologie()==NULL))//((origine != MAGIC::ORIGINE::IMPOSE)&& (origine != MAGIC::ORIGINE::TRIANGULATION_ARETE)&&( noeud_i2->get_lien_topologie()==NULL))/////////////////////////////////////////////////////////////////////////////
1802 {
1803 noeud_i->change_x(nouv_position_x[ind_noeud]);
1804 noeud_i->change_y(nouv_position_y[ind_noeud]);
1805 noeud_i->change_z(nouv_position_z[ind_noeud]);
1806 ind_noeud++;
1807 }
1808 }
1809 return 1;
1810
1811 }
1812
1813 //---------------------------- Lissage Taubin 1995 par Bryan McKenzie-------------------------------------------------------------------------
1814 int MGOPT_POSTTRAITEMENT::lissage_Taubin1995(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double lambda, double nu, int iter_max)
1815 {
1816 int compt = 0;
1817 int ok = 0;
1818 char mess[200];
1819 do
1820 {
1821 int taubin=cycle_taubin1995(mg_mai, gest2, lambda, nu);
1822 sprintf(mess, "cycle taubin faite"); affiche(mess);
1823 compt = compt+1;
1824 sprintf(mess, "incrémentation compt %d",compt); affiche(mess);
1825 }while (compt < iter_max);
1826 return compt;
1827 sprintf(mess, "fin Taubin1995"); affiche(mess);
1828 }
1829
1830
1831 //---------------------------- Varience/Ecart type/ pour McKenzie2016-------------------------------------------------------------------------
1832 void MGOPT_POSTTRAITEMENT::varience_McKenzie2016(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double *vari0,double *vari1,double *vari2)
1833 {
1834 char mess[300];
1835 int ecart_type = 0;
1836 double s_i = 0.0;
1837 double phi_im = 0.0;
1838 double ecart_type_i = 0.;
1839 ecart_type = 0;
1840 do
1841 {
1842
1843 TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
1844 TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales2;
1845 TPL_LISTE_ENTITE<double> liste_wij;
1846 LISTE_MG_TRIANGLE::iterator it_tri;
1847 int k = 0; //pour identifier les triangles pour liste_normales et liste_wij
1848 for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
1849 {
1850 MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
1851 OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
1852 normal_f_i.norme();
1853 liste_normales2.ajouter(normal_f_i);
1854 //Remplissage de la liste des voisins du triangle i
1855 TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
1856 MG_NOEUD* noeud1 = mgtri_i->get_noeud1();
1857 double nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
1858 for (int j = 0;j<nb_voisins1;j++)
1859 {
1860 MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
1861 liste_voisins.ajouter(mgtri_1);
1862 }
1863 MG_NOEUD* noeud2 = mgtri_i->get_noeud2();
1864 double nb_voisins2 = noeud2->get_lien_triangle()->get_nb();
1865 for (int j = 0;j<nb_voisins2;j++)
1866 {
1867 MG_TRIANGLE_PEAU* mgtri_2 = (MG_TRIANGLE_PEAU*) noeud2->get_lien_triangle()->get(j);
1868 liste_voisins.ajouter(mgtri_2);
1869 }
1870 MG_NOEUD* noeud3 = mgtri_i->get_noeud3();
1871 double nb_voisins3 = noeud3->get_lien_triangle()->get_nb();
1872 for (int j = 0;j<nb_voisins3;j++)
1873 {
1874 MG_TRIANGLE_PEAU* mgtri_3 = (MG_TRIANGLE_PEAU*) noeud3->get_lien_triangle()->get(j);
1875 liste_voisins.ajouter(mgtri_3);
1876 }
1877 liste_voisins.supprimer(mgtri_i);
1878 int nb_voisins = liste_voisins.get_nb();
1879 double w_ij = 1./nb_voisins;
1880 double phi_i_min = 10.;
1881 s_i = 0.0;
1882 phi_im = 0.0;
1883 double *phi_ij = new double[nb_voisins];
1884 OT_VECTEUR_3D normal_f_i_mean(0.,0.,0.);
1885 normal_f_i_mean = normal_f_i;
1886 OT_VECTEUR_3D eta_i(0.,0.,0.);
1887 TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
1888 int j = 0;
1889 double un_sur_pi = 1./M_PI;
1890 for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
1891 {
1892 OT_VECTEUR_3D normal_f_j = mgtri_j->calcul_normal();
1893 //1-Calculer la normale moyenne pour chaque triangle
1894 normal_f_i_mean = normal_f_i_mean + normal_f_j;
1895 //2.1-On calcule l'angle entre normal_f_i et normal_f_j pour j allant de 1 a Nb_voisins
1896 double prod_scalaire = normal_f_i*normal_f_j;
1897 if (prod_scalaire > 1.)
1898 {
1899 prod_scalaire = 1.;
1900 }
1901 if (prod_scalaire < -1.)
1902 {
1903 prod_scalaire = -1.;
1904 }
1905 phi_ij[j] = acos(prod_scalaire)*un_sur_pi;
1906 //2.2-On trouve le plus petit des angles et la normale heta_i correspondante
1907 if (phi_ij[j] < phi_i_min)
1908 {
1909 phi_i_min = phi_ij[j];
1910 eta_i = normal_f_j;
1911 }
1912 //3.1-On calcule l'angle moyen phi_im
1913 phi_im = phi_im + w_ij*phi_ij[j];
1914 j++;
1915 }
1916 normal_f_i_mean.norme();
1917 j = 0;
1918 for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
1919 {
1920 //3.2-Calcul de s_i selon la variance
1921 s_i = s_i + w_ij*pow((phi_ij[j] - phi_im),2);
1922 j++;
1923 }
1924 delete[] phi_ij;
1925 }
1926 ecart_type_i = sqrt(s_i);
1927 ecart_type = 1;
1928 }
1929 while (ecart_type == 0);
1930 ecart_type = 0;//reinitialisation d'ecart_type pour la prochainepasse
1931 *vari0 = phi_im;
1932 *vari1 = s_i;
1933 *vari2 = ecart_type_i;
1934 //return &varience[0];
1935
1936 }
1937
1938
1939 //---------------------------- Lissage McKenzie2016 par Bryan Mckenzie 2016-------------------------------------------------------------------------
1940 int MGOPT_POSTTRAITEMENT::lissage_McKenzie2016(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double lambda, double nu, double epsilon,double sigma,double gamma_,int filtre, int iter_max, int itr_taubin)
1941 {
1942 // Lissage global avec methode de lissage inspiree de Taubin 1995 et Chen 2008 par Bryan McKenzie 2016
1943
1944
1945 //sauvegarde des differents position de noeuds pour conservation des position
1946 std::vector<double> coord_x;
1947 std::vector<double> coord_y;
1948 std::vector<double> coord_z;
1949 std::vector<double> vecteur_original;
1950 LISTE_MG_NOEUD::iterator it_no;
1951 int ii=0;
1952 for (MG_NOEUD* noeud=mg_mai->get_premier_noeud(it_no);noeud!=NULL;noeud=mg_mai->get_suivant_noeud(it_no))
1953 {
1954 MG_NOEUD* tri=(MG_NOEUD*)noeud;
1955 if(tri->get_origine()!=MAGIC::ORIGINE::MAILLEUR_AUTO)
1956 {
1957 coord_x.push_back(noeud->get_x());
1958 coord_y.push_back(noeud->get_y());
1959 coord_z.push_back(noeud->get_z());
1960 vecteur_original.push_back(sqrt(coord_x[ii]*coord_x[ii] + coord_y[ii]*coord_y[ii] + coord_z[ii]*coord_z[ii]));
1961 ii++;
1962 }
1963
1964 }
1965
1966
1967
1968 char mess[300];
1969 sprintf(mess, "Depart Taubin1995"); affiche(mess);
1970 double v_initial = 0.0;
1971 double v_final = 0.0;
1972 double v_taubin = 0.0;
1973 double v_taubin_prime = 0.0;
1974 double d_initial = 0.0;
1975 double d_final = 0.0;
1976 double d_taubin = 0.0;
1977 double nu_ = nu;
1978 double sigma_ = sigma;//variable utiliser pour sassuré que la pièce ne grossie pas dan taubin_chen 2
1979 int ecart_type = 0;
1980 double s_i = 0.0;
1981 double phi_im = 0.0;
1982 double ecart_type_i = 0.;
1983 int compteur_triangle_retourne =0;
1984 int compteur_total = 0;
1985 int compteur_taubin_initial = 0;
1986 int compteur_chen = 0;
1987 int compteur_taubin_chen = 0;
1988 int compteur_taubin_final = 0;
1989 int compteur_temp = 0;
1990 int compteur_taubin2D = 0;
1991 int triangle_final_retourne =0;
1992 int iter_taubin = itr_taubin;
1993
1994
1995
1996 int activationPI = 1; // variable qui active l'adjustement du volume dans taubin1995
1997 int activationPI2 = 1; // variable qui active l'adjustement du volume dans taubin_chen
1998 double condition_de_varience = 1000.0; //variable de sortie pour les prenières iteration taubin pour s'assuré que la surface a une certaine surface lisse
1999 int taubin_premiere_passe = 0;
2000
2001 /// Calcul du volume initial
2002 MG_MAILLAGE_OUTILS ot;
2003 v_initial=ot.get_volume(mg_mai);
2004
2005
2006
2007 //******************parti taubin pur***********************************//
2008
2009
2010
2011 //calcule de la varience et l'écart type pour la surface du maillaige au début
2012 double varience_ori0 = 10.0;
2013 double varience_ori1 = 10.0;
2014 double varience_ori2 = 10.0;
2015 varience_McKenzie2016(mg_mai,gest2, &varience_ori0, &varience_ori1, &varience_ori2);
2016 sprintf(mess," Moyenne des normes des triangles de la surface : %f",varience_ori0); affiche(mess);
2017 sprintf(mess," Varience des normes des triangles de la surface : %f",varience_ori1); affiche(mess);
2018 sprintf(mess," Ecart type des normes des triangles de la surface : %f",varience_ori2); affiche(mess);
2019 v_taubin=v_initial; //ajout vincent
2020 do
2021 {
2022 v_taubin_prime = v_taubin;//assignation du volume de l'itération précédente pour fair la loi de controle
2023 int taubin=cycle_taubin1995(mg_mai, gest2, lambda, nu);
2024
2025 /// Calcul du volume
2026 MG_MAILLAGE_OUTILS ot;
2027 v_taubin=ot.get_volume(mg_mai);
2028
2029 //implemantion d'un PI por controler le volume
2030 if (activationPI==1)
2031 {
2032 if (((fabs(v_taubin)-fabs(v_taubin_prime))/fabs(v_taubin))<-0.1)//fine adjustment makes part bigger if new iteration of part is smaller than the preceeding pass and volume is less than 10% bigger than original
2033 {
2034 nu_ =nu_-0.00005;
2035 }
2036 if (((fabs(v_taubin)-fabs(v_taubin_prime))/fabs(v_taubin))>0.1)//fine adjustment makes part smaller if part is much bigger than preceeding pass and the volume is larger than original
2037 {
2038 nu_ = nu_+0.00005;
2039 }
2040 if (fabs(v_taubin)>(fabs(v_initial)*1.015))//gross adjustment makes part smaller if part is 3% bigger than original make smaller
2041 {
2042 nu_ =nu_+0.0005;
2043 }
2044 if (fabs(v_taubin)<(fabs(v_initial)*0.985))//gross adjustment makes part bigger if it is smaller than the original
2045 {
2046 nu_ =nu_-0.0005;
2047 }
2048 }
2049 //double un_sur_pi = 1./M_PI;
2050
2051
2052
2053
2054 //code qui vérifie la qualité des different trianges
2055 MAILLEUR_ANALYSE m3d(mg_mai);
2056 m3d.change_eps_angle_retourne(0.03000);
2057 m3d.change_borne(0.1,0.2,0.5);
2058 double qualmin,qualmax,qualmoy;
2059 int tab[7];
2060 m3d.analyse_qualite_maillage_2D(NULL,qualmin,qualmax,qualmoy,tab);
2061 int nbretourne=tab[0];
2062 int itri=0;
2063
2064
2065
2066 if (tab[0]>0)
2067 {
2068 compteur_triangle_retourne++;
2069 compteur_temp--;
2070 sprintf(mess, "Triangle retourné");
2071 affiche(mess);
2072 }
2073 compteur_taubin_initial++;
2074 compteur_total++;
2075 compteur_temp++;
2076
2077
2078 //calcule de la varience et l'écart type pour la surface du maillaige au début
2079 double varience0 = 0.0;
2080 double varience1 =0.0;
2081 double varience2 =0.0;
2082 varience_McKenzie2016(mg_mai,gest2, &varience0, &varience1, &varience2);
2083
2084 if (((compteur_temp >= iter_taubin)&&(varience1<0.0275))||(compteur_taubin_initial >(2*iter_taubin)))
2085 {
2086 taubin_premiere_passe = 1;
2087 }
2088
2089
2090 }
2091 while (taubin_premiere_passe == 0);
2092
2093
2094
2095
2096
2097
2098
2099
2100
2101 //code d'utilisation de chen pour garder la forme et de taubin pour retirer les triangles retournés
2102 sprintf(mess, "Début de Chen pour %d iterations avec Taubin qui résoud les triangles retournés",iter_max); affiche(mess);
2103
2104 do
2105 {
2106
2107 int nbiter=lissage_chen2008(mg_mai,gest2,sigma,gamma_,epsilon,1);
2108
2109 // code de vérifiaction de triangle retourner
2110 MAILLEUR_ANALYSE m3d(mg_mai);
2111 m3d.change_eps_angle_retourne(0.03000);
2112 m3d.change_borne(0.1,0.2,0.5);
2113 double qualmin,qualmax,qualmoy;
2114 int tab[7];
2115 m3d.analyse_qualite_maillage_2D(NULL,qualmin,qualmax,qualmoy,tab);
2116 int triangleretourner =tab[0];
2117
2118 //code de retourenage de triangle inverser avec taubin 1995
2119 if (triangleretourner >0)
2120 {
2121 do
2122 {
2123 v_taubin_prime = v_taubin;//assignation du volume de l'itération précédente pour fair la loi de controle
2124 int taubin=cycle_taubin1995(mg_mai, gest2, lambda, nu);
2125
2126 /// Calcul du volume
2127 MG_MAILLAGE_OUTILS ot;
2128 v_taubin=ot.get_volume(mg_mai);
2129
2130 //implemantion d'un PI por controler le volume
2131 if (activationPI2 = 1)
2132 {
2133
2134 if (((fabs(v_taubin)-fabs(v_taubin_prime))/fabs(v_taubin))<-0.1)//fine adjustment makes part bigger if new iteration of part is smaller than the preceeding pass and volume is less than 10% bigger than original
2135 {
2136 nu_ =nu_-0.00005;
2137 }
2138
2139 if (((fabs(v_taubin)-fabs(v_taubin_prime))/fabs(v_taubin))>0.1)//fine adjustment makes part smaller if part is much bigger than preceeding pass and the volume is larger than original
2140 {
2141 nu_ = nu_+0.00005;
2142 }
2143
2144 if (fabs(v_taubin)>(fabs(v_initial)*1.015))//gross adjustment makes part smaller if part is 3% bigger than original make smaller
2145 {
2146 nu_ =nu_+0.0005;
2147 }
2148 if (fabs(v_taubin)<(fabs(v_initial)*0.985))//gross adjustment makes part bigger if it is smaller than the original
2149 {
2150 nu_ =nu_-0.0005;
2151 }
2152 }
2153
2154 //code qui vérifie la qualité des different trianges
2155 m3d.analyse_qualite_maillage_2D(NULL,qualmin,qualmax,qualmoy,tab);
2156 if (tab[0]>0)
2157 {
2158 triangleretourner = 1;
2159 sprintf(mess, "Trianlge retourner pas régler"); affiche(mess);
2160
2161 }
2162 else
2163 {
2164 triangleretourner = 0;
2165 sprintf(mess, "Trianlge retourner régler"); affiche(mess);
2166 }
2167 compteur_triangle_retourne++;
2168 compteur_taubin_chen++;
2169 compteur_total++;
2170 }
2171
2172 while ((triangleretourner > 0));
2173
2174 }
2175 compteur_total++;
2176 compteur_chen++;
2177
2178
2179 }
2180 while (compteur_chen < iter_max);
2181
2182 sprintf(mess, "Début du retour de volume avec Taubin1995 qui résoud les triangles retournés"); affiche(mess);
2183 // passe final de taubin pour revenir au volume initial
2184 v_taubin=ot.get_volume(mg_mai);
2185
2186 if ((fabs(fabs(v_initial)-fabs(v_taubin)))>(0.02*(fabs(v_initial))))
2187 {
2188 do
2189 {
2190 v_taubin_prime = v_taubin;
2191 int taubin=cycle_taubin1995(mg_mai, gest2, lambda, nu);
2192
2193 /// Calcul du volume
2194 v_taubin=ot.get_volume(mg_mai);
2195
2196 //implemantion d'un PDF pour controler le volume
2197
2198 if (((fabs(v_taubin)-fabs(v_taubin_prime))/fabs(v_taubin))<-0.1)//fine adjustment makes part bigger if new iteration of part is smaller than the preceeding pass and volume is less than 10% bigger than original
2199 {
2200 nu_ =nu_-0.00002;
2201 sprintf(mess, "ajustement mineur 1"); affiche(mess);
2202 }
2203
2204 if (((fabs(v_taubin)-fabs(v_taubin_prime))/fabs(v_taubin))>0.1)//fine adjustment makes part smaller if part is much bigger than preceeding pass and the volume is larger than original
2205 {
2206 nu_ = nu_+0.00002;
2207 sprintf(mess, "ajustement mineur 2"); affiche(mess);
2208 }
2209
2210 if (fabs(v_taubin)>(fabs(v_initial)*1.015))//gross adjustment makes part smaller if part is 1.5% bigger than original make smaller
2211 {
2212 nu_ =nu_+0.0005;
2213 sprintf(mess, "ajustement majeur 1"); affiche(mess);
2214 }
2215 if (fabs(v_taubin)<(fabs(v_initial)*0.985))//gross adjustment makes part bigger if it is more than 1.5% smaller than the original
2216 {
2217 nu_ =nu_-0.0005;
2218 sprintf(mess, "ajustement majeur 2"); affiche(mess);
2219 }
2220 sprintf(mess, "Volume: %f\n old vol: %f\n nu_: %f\n",v_taubin,v_taubin_prime, nu_); affiche(mess);
2221
2222
2223
2224
2225 //code qui vérifie la qualité des different trianges
2226 MAILLEUR_ANALYSE m3d(mg_mai);
2227 m3d.change_eps_angle_retourne(0.03000);
2228 m3d.change_borne(0.1,0.2,0.5);
2229 double qualmin,qualmax,qualmoy;
2230 int tab[7];
2231 m3d.analyse_qualite_maillage_2D(NULL,qualmin,qualmax,qualmoy,tab);
2232 int triangleretourner =tab[0];
2233 if (tab[0]>0)
2234 {
2235 compteur_triangle_retourne++;
2236 sprintf(mess, "Triangle retourné");
2237 affiche(mess);
2238 triangle_final_retourne =1;
2239 }
2240 else{triangle_final_retourne =0;}
2241
2242 compteur_taubin_final++;
2243 compteur_total++;
2244
2245 }
2246
2247 while ((triangle_final_retourne ==0)&&(fabs(v_initial-v_taubin)>(fabs(0.02*v_initial))));//((compteur_ < 1000));
2248
2249 }
2250
2251 ////////////////////////////calcul du mouvement des noeuds par rapport a l'original
2252 /******************code qui calcule la moyenne de pourcentage de déplacement des differents noeuds*/
2253 //LISTE_MG_NOEUD::iterator it_no;
2254 ii=0;
2255 int numdepla=0;
2256 double distot=0;
2257 double distotx=0;
2258 double distoty=0;
2259 double distotz=0;
2260 double dismin=1e300,dismax=-1e300;
2261 double disminx=1e300,dismaxx=-1e300;
2262 double disminy=1e300,dismaxy=-1e300;
2263 double disminz=1e300,dismaxz=-1e300;
2264 double xmin=1e300,xmax=-1e300;
2265 double ymin=1e300,ymax=-1e300;
2266 double zmin=1e300,zmax=-1e300;
2267 for (MG_NOEUD* noeud=mg_mai->get_premier_noeud(it_no);noeud!=NULL;noeud=mg_mai->get_suivant_noeud(it_no))
2268 {
2269 if(noeud->get_origine()!=MAGIC::ORIGINE::IMPOSE)////////
2270 {
2271 double x=noeud->get_x();
2272 double y=noeud->get_y();
2273 double z=noeud->get_z();
2274 if (x<xmin) xmin=x;
2275 if (y<ymin) ymin=y;
2276 if (z<zmin) zmin=z;
2277 if (x>xmax) xmax=x;
2278 if (y>ymax) ymax=y;
2279 if (z>zmax) zmax=z;
2280 double dx=fabs(noeud->get_x()-coord_x[ii]);
2281 double dy=fabs(noeud->get_y()-coord_y[ii]);
2282 double dz=fabs(noeud->get_z()-coord_z[ii]);
2283 double dis=sqrt(dx*dx+dy*dy+dz*dz);
2284 if (dx>dismaxx) dismaxx=dx;
2285 if (dx<disminx) disminx=dx;
2286 if (dy>dismaxy) dismaxy=dy;
2287 if (dy<disminy) disminy=dy;
2288 if (dz>dismaxz) dismaxz=dz;
2289 if (dz<disminz) disminz=dz;
2290 if (dis>dismax) dismax=dis;
2291 if (dis<dismin) dismin=dis;
2292 distot=distot+dis;
2293 distotx=distoty+dx;
2294 distoty=distoty+dy;
2295 distotz=distotz+dz;
2296 numdepla++;
2297 }
2298 ii++;
2299 }
2300 affiche((char*)"");
2301 sprintf(mess, "Déplacement absolu minimum en X : %f en Y : %f en Z : %f en amplitude : %f",disminx,disminy,disminz,dismin); affiche(mess);
2302 sprintf(mess, "Déplacement absolu maximum en X : %f en Y : %f en Z : %f en amplitude : %f",dismaxx,dismaxy,dismaxz,dismax); affiche(mess);
2303 sprintf(mess, "Déplacement absolu moyen (tous les noeuds) en X : %f en Y : %f en Z : %f en amplitude : %f",distotx/ii,distoty/ii,distotz/ii,distot/ii); affiche(mess);
2304 sprintf(mess, "Déplacement absolu moyen (les noeuds du non design) en X : %f en Y : %f en Z : %f en amplitude : %f",distotx/numdepla,distoty/numdepla,distotz/numdepla,distot/numdepla); affiche(mess);
2305 double dx=xmax-xmin;
2306 double dy=ymax-ymin;
2307 double dz=zmax-zmin;
2308 double dmax=dx;
2309 if (dy>dmax) dmax=dy;
2310 if (dz>dmax) dmax=dz;
2311 sprintf(mess, "Déplacement relatif moyen (tous les noeuds) en X : %f%% en Y : %f%% en Z : %f%% en amplitude : %f%%",distotx/ii/dx*100.,distoty/ii/dy*100.,distotz/ii/dz*100.,distot/ii/dmax*100.); affiche(mess);
2312 sprintf(mess, "Déplacement relatif moyen (les noeuds du non design) en X : %f%% en Y : %f%% en Z : %f%% en amplitude : %f%%",distotx/numdepla/dx*100.,distoty/numdepla/dy*100.,distotz/numdepla/dz*100.,distot/numdepla/dmax*100.); affiche(mess);
2313
2314
2315
2316 affiche((char*)"");
2317
2318 affiche((char*)"");
2319
2320
2321
2322 /// Calcul du volume final
2323 v_final=ot.get_volume(mg_mai);
2324
2325 double delta_vol = v_final-v_initial;
2326 double vol_pourcent = (v_final/v_initial)*100;
2327 double delta_vol_pourcent = vol_pourcent-100;
2328 sprintf(mess, "Volume initial: %f",v_initial);affiche(mess);
2329 sprintf(mess, "Volume final: %f",v_final);affiche(mess);
2330 sprintf(mess, "Volume final pourcentage: %f%%",vol_pourcent);affiche(mess);
2331 sprintf(mess, "Delta Volume: %f",delta_vol);affiche(mess);
2332 sprintf(mess, "Delta Volume pourcentage: %f%%",delta_vol_pourcent);affiche(mess);
2333 sprintf(mess, " Iteration inital Taubin1995: %i",compteur_taubin_initial); affiche(mess);
2334 sprintf(mess, " Iteration Chen2008: %i",compteur_chen); affiche(mess);
2335 sprintf(mess, " Iteration Taubin1995 dans Chen2008: %i",compteur_taubin_chen); affiche(mess);
2336 sprintf(mess, " epsilon: %f",epsilon); affiche(mess);
2337 sprintf(mess, " sigma: %f",sigma); affiche(mess);
2338 sprintf(mess, " gamma_: %f",gamma_); affiche(mess);
2339 sprintf(mess, " filtre: %i",filtre); affiche(mess);
2340 sprintf(mess, " Iteration final Taubin1995: %i",compteur_taubin_final); affiche(mess);
2341 sprintf(mess, " lambda(depart): %f",lambda); affiche(mess);
2342 sprintf(mess, " nu(depart):%f",nu); affiche(mess);
2343 sprintf(mess, " nu(fin):%f",nu_); affiche(mess);
2344 sprintf(mess, " Iteration triangle retourné: %i",compteur_triangle_retourne); affiche(mess);
2345 sprintf(mess, " Iteration Total: %i",compteur_total); affiche(mess);
2346
2347
2348
2349 //verification de nombre de triange retournée
2350 MAILLEUR_ANALYSE m3d(mg_mai);
2351 double eps_angle_retourne=0.03;
2352 m3d.change_eps_angle_retourne(eps_angle_retourne);
2353 double borne1=0.1;
2354 double borne2=0.2;
2355 double borne3=0.5;
2356 m3d.change_borne(borne1,borne2,borne3);
2357 double qualmin,qualmax,qualmoy;
2358 int tab[7];
2359 m3d.analyse_qualite_maillage_2D(NULL,qualmin,qualmax,qualmoy,tab);
2360 sprintf(mess,"\n\nANALYSEUR DE MG MAILLAGE 2D"); affiche(mess);
2361 sprintf(mess," Nombre de triangles avec des noeuds fusionnés : %d",tab[6]); affiche(mess);
2362 sprintf(mess," Nombre de triangles avec 1 voisin manquant : %d",tab[5]); affiche(mess);
2363 sprintf(mess," critere d'angle limite pour l'inversion : %lf",eps_angle_retourne); affiche(mess);
2364 sprintf(mess," qualite moyenne des triangles de frontiere : %lf",qualmoy); affiche(mess);
2365 sprintf(mess," qualite min des triangles de frontiere : %lf",qualmin); affiche(mess);
2366 sprintf(mess," qualite max des triangles de frontiere : %lf",qualmax); affiche(mess);
2367 sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne3,tab[4],tab[4]*100./mg_mai->get_nb_mg_triangle()); affiche(mess);
2368 sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne2,tab[3],tab[3]*100./mg_mai->get_nb_mg_triangle()); affiche(mess);
2369 sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne1,tab[2],tab[2]*100./mg_mai->get_nb_mg_triangle()); affiche(mess);
2370 sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",0.,tab[1],tab[1]*100./mg_mai->get_nb_mg_triangle()); affiche(mess);
2371 sprintf(mess," nombre de triangles de frontiere retournes : %d (%.2lf%%)",tab[0],tab[0]*100./mg_mai->get_nb_mg_triangle()); affiche(mess);
2372
2373 //calcule de la varience et l'écart type pour la surface du maillaige
2374 double varience00 =0.0;
2375 double varience11 =0.0;
2376 double varience22 =0.0;
2377 varience_McKenzie2016(mg_mai,gest2, &varience00, &varience11, &varience22);
2378 sprintf(mess,"\n Moyenne des normes des triangles de la surface : %f",varience00); affiche(mess);
2379 sprintf(mess," Varience des normes des triangles de la surface : %f",varience11); affiche(mess);
2380 sprintf(mess," Ecart type des normes des triangles de la surface : %f",varience22); affiche(mess);
2381
2382 return compteur_total;
2383 }
2384
2385
2386
2387
2388 //---------------------------- Extract skin -------------------------------------------------------------------------
2389 int MGOPT_POSTTRAITEMENT::extract_skin(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2,int &nbpeau,int &nbmaniare,int &nbmanino, int *mai2_id)
2390 {
2391
2392 //etape 0 - On commence par mettre a zero tous les nouveau_numero des triangles et des noeuds du maillage
2393 LISTE_MG_TRIANGLE::iterator it_tri;
2394 for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
2395 {
2396 mgtri->change_nouveau_numero(0);
2397 mgtri->change_origine(MAGIC::ORIGINE::MAILLEUR_AUTO);
2398 }
2399 LISTE_MG_NOEUD::iterator it_noeud;
2400 for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
2401 {
2402 mgnoeud->change_nouveau_numero(0);
2403 }
2404
2405 //etape 1 - On boucle ensuite tous les tetraedres pour ajouter un compteur du nombre de fois qu'un triangle appartient a 1 tetra
2406 LISTE_MG_TETRA::iterator it_tetra;
2407 for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
2408 {
2409 int origine = mgtet->get_origine();
2410 if (origine==MAGIC::ORIGINE::IMPOSE)
2411 {
2412 mgtet->get_triangle1()->change_origine(MAGIC::ORIGINE::IMPOSE);
2413 mgtet->get_triangle2()->change_origine(MAGIC::ORIGINE::IMPOSE);
2414 mgtet->get_triangle3()->change_origine(MAGIC::ORIGINE::IMPOSE);
2415 mgtet->get_triangle4()->change_origine(MAGIC::ORIGINE::IMPOSE);
2416 }
2417 if (origine==MAGIC::ORIGINE::OPTIMISE)
2418 {
2419 if (mgtet->get_triangle1()->get_origine()!=MAGIC::ORIGINE::IMPOSE) mgtet->get_triangle1()->change_origine(MAGIC::ORIGINE::OPTIMISE);
2420 if (mgtet->get_triangle2()->get_origine()!=MAGIC::ORIGINE::IMPOSE) mgtet->get_triangle2()->change_origine(MAGIC::ORIGINE::OPTIMISE);
2421 if (mgtet->get_triangle3()->get_origine()!=MAGIC::ORIGINE::IMPOSE) mgtet->get_triangle3()->change_origine(MAGIC::ORIGINE::OPTIMISE);
2422 if (mgtet->get_triangle4()->get_origine()!=MAGIC::ORIGINE::IMPOSE) mgtet->get_triangle4()->change_origine(MAGIC::ORIGINE::OPTIMISE);
2423 }
2424
2425 if (((origine == MAGIC::ORIGINE::OPTIMISE) && (etat[(MAGIC::ORIGINE::OPTIMISE-1000)/10]==1) ) ||
2426 ((origine == MAGIC::ORIGINE::IMPOSE) && (etat[(MAGIC::ORIGINE::IMPOSE-1000)/10]==1) ) ||
2427 ((origine == MAGIC::ORIGINE::MAILLEUR_AUTO) && (etat[(MAGIC::ORIGINE::MAILLEUR_AUTO-1000)/10]==1) ) ||
2428 ((origine == MAGIC::ORIGINE::TRIANGULATION) && (etat[(MAGIC::ORIGINE::TRIANGULATION-1000)/10]==1) ) ||
2429 ((origine == MAGIC::ORIGINE::MODIFICATION) && (etat[(MAGIC::ORIGINE::MODIFICATION-1000)/10]==1) ) ||
2430 ((origine == MAGIC::ORIGINE::DUPLIQUER) && (etat[(MAGIC::ORIGINE::DUPLIQUER-1000)/10]==1) ) )
2431
2432 {
2433 int num1 = mgtet->get_triangle1()->get_nouveau_numero();
2434 int num2 = mgtet->get_triangle2()->get_nouveau_numero();
2435 int num3 = mgtet->get_triangle3()->get_nouveau_numero();
2436 int num4 = mgtet->get_triangle4()->get_nouveau_numero();
2437 num1++;
2438 num2++;
2439 num3++;
2440 num4++;
2441 mgtet->get_triangle1()->change_nouveau_numero(num1);
2442 mgtet->get_triangle2()->change_nouveau_numero(num2);
2443 mgtet->get_triangle3()->change_nouveau_numero(num3);
2444 mgtet->get_triangle4()->change_nouveau_numero(num4);
2445 }
2446 }
2447
2448 //etape 2 - On boucle l'ensemble des triangles identifies a l'etape 1 pour identifier les noeuds leur appartenant
2449 for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
2450 {
2451 int num = mgtri->get_nouveau_numero();
2452 if (num == 1)
2453 {
2454 MG_NOEUD* noeud1 = mgtri->get_noeud1();
2455 MG_NOEUD* noeud2 = mgtri->get_noeud2();
2456 MG_NOEUD* noeud3 = mgtri->get_noeud3();
2457 noeud1->change_nouveau_numero(1);
2458 noeud2->change_nouveau_numero(1);
2459 noeud3->change_nouveau_numero(1);
2460 if (mgtri->get_origine()==MAGIC::ORIGINE::IMPOSE)
2461 {
2462 noeud1->change_origine(MAGIC::ORIGINE::IMPOSE);
2463 noeud2->change_origine(MAGIC::ORIGINE::IMPOSE);
2464 noeud3->change_origine(MAGIC::ORIGINE::IMPOSE);
2465 }
2466 }
2467 }
2468
2469 //etape 3 - On cree un nouveau maillage pour la peau
2470 gest2.vide();
2471 MG_GEOMETRIE* geo=new MG_GEOMETRIE((char*)"VIRTUEL",(char*)"VIRTUEL");
2472 geo->change_valeur_unite(mg_mai->get_mg_geometrie()->get_valeur_unite());
2473 gest2.ajouter_mg_geometrie(geo);
2474 MG_VOLUME_ELEMENT* vol=new MG_VOLUME_ELEMENT();
2475 geo->ajouter_mg_volume(vol);
2476 vol->copie_ccf(*(mg_mai->get_mg_geometrie()->get_mg_volume(0)));
2477 std::map<unsigned long,unsigned long> correspondidsom;
2478 std::map<unsigned long,unsigned long> correspondid;
2479 LISTE_MG_SOMMET::iterator its;
2480 for (MG_SOMMET* som=mg_mai->get_mg_geometrie()->get_premier_sommet(its);som!=NULL;som=mg_mai->get_mg_geometrie()->get_suivant_sommet(its))
2481 {
2482 if (som->get_nb_ccf()>0)
2483 {
2484 MG_SOMMET_NOEUD* somno=new MG_SOMMET_NOEUD();
2485 geo->ajouter_mg_sommet(somno);
2486 somno->copie_ccf(*som);
2487 correspondidsom[som->get_id()]=somno->get_id();
2488 }
2489 }
2490 LISTE_MG_ARETE::iterator ita;
2491 for (MG_ARETE* are=mg_mai->get_mg_geometrie()->get_premier_arete(ita);are!=NULL;are=mg_mai->get_mg_geometrie()->get_suivant_arete(ita))
2492 {
2493 if (are->get_nb_ccf()>0)
2494 {
2495 MG_ARETE_ELEMENT* areel=new MG_ARETE_ELEMENT();
2496 geo->ajouter_mg_arete(areel);
2497 areel->copie_ccf(*are);
2498 correspondid[are->get_id()]=areel->get_id();
2499 }
2500 }
2501 LISTE_MG_FACE::iterator itf;
2502 for (MG_FACE* face=mg_mai->get_mg_geometrie()->get_premier_face(itf);face!=NULL;face=mg_mai->get_mg_geometrie()->get_suivant_face(itf))
2503 {
2504 if (face->get_nb_ccf()>0)
2505 {
2506 MG_FACE_ELEMENT* facel=new MG_FACE_ELEMENT();
2507 geo->ajouter_mg_face(facel);
2508 facel->copie_ccf(*face);
2509 correspondid[face->get_id()]=facel->get_id();
2510 }
2511 }
2512 MG_FACE_ELEMENT* faceex=new MG_FACE_ELEMENT();
2513 geo->ajouter_mg_face(faceex);
2514 MG_FACE_ELEMENT* facenoex=new MG_FACE_ELEMENT();
2515 geo->ajouter_mg_face(facenoex);
2516 MG_MAILLAGE* mai2 = new MG_MAILLAGE(geo);
2517 //MG_MAILLAGE* mai2 = new MG_MAILLAGE(NULL);
2518 gest2.ajouter_mg_maillage(mai2);
2519
2520 //etape 4 - On boucle l'ensemble des noeuds identifies a l'etape 2 pour les recreer dans le second maillage
2521 for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
2522 {
2523 int num = mgnoeud->get_nouveau_numero();
2524 if (num == 1)
2525 {
2526 double x = mgnoeud->get_x();
2527 double y = mgnoeud->get_y();
2528 double z = mgnoeud->get_z();
2529 int origine=MAGIC::ORIGINE::TRIANGULATION;
2530 int nb_tri=mgnoeud->get_lien_triangle()->get_nb();
2531 bool tri2d=false,tri3d=false;
2532 for (int i=0;i<nb_tri;i++)
2533 {
2534 if (mgnoeud->get_lien_triangle()->get(i)->get_nouveau_numero()==1)
2535 {
2536 if (mgnoeud->get_lien_triangle()->get(i)->get_lien_topologie()->get_dimension()==2) tri2d=true;
2537 if (mgnoeud->get_lien_triangle()->get(i)->get_lien_topologie()->get_dimension()==3) tri3d=true;
2538 }
2539 }
2540 if ((tri2d)&&(tri3d)) origine=MAGIC::ORIGINE::TRIANGULATION_ARETE;
2541 if (mgnoeud->get_lien_topologie()->get_dimension()<2) origine=MAGIC::ORIGINE::TRIANGULATION_ARETEORIGINE;
2542 if (mgnoeud->get_origine()==MAGIC::ORIGINE::IMPOSE) origine=MAGIC::ORIGINE::IMPOSE;
2543 MG_ELEMENT_TOPOLOGIQUE* topo=NULL;
2544 std::map<unsigned long,unsigned long>::iterator it=correspondidsom.find(mgnoeud->get_lien_topologie()->get_id());
2545 if (it!=correspondidsom.end())
2546 topo=geo->get_mg_sommetid(it->second);
2547 MG_NOEUD* noeud1 = new MG_NOEUD(topo,x,y,z,origine);
2548 if (it!=correspondidsom.end())
2549 ((MG_SOMMET_NOEUD*)topo)->change_mg_noeud(noeud1);
2550 mai2->ajouter_mg_noeud(noeud1);
2551 mgnoeud->change_nouveau_numero(noeud1->get_id());
2552 noeud1->change_nouveau_numero(mgnoeud->get_id());
2553 }
2554 }
2555
2556 // etape 5 - On boucle l'ensemble des triangles identifies a l'etape 1 pour les recreer dans le maillage 2
2557 for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
2558 {
2559 int originetet=mgtet->get_origine();
2560 if (((originetet == MAGIC::ORIGINE::OPTIMISE) && (etat[(MAGIC::ORIGINE::OPTIMISE-1000)/10]==1) ) ||
2561 ((originetet == MAGIC::ORIGINE::IMPOSE) && (etat[(MAGIC::ORIGINE::IMPOSE-1000)/10]==1) ) ||
2562 ((originetet == MAGIC::ORIGINE::MAILLEUR_AUTO) && (etat[(MAGIC::ORIGINE::MAILLEUR_AUTO-1000)/10]==1) ) ||
2563 ((originetet == MAGIC::ORIGINE::TRIANGULATION) && (etat[(MAGIC::ORIGINE::TRIANGULATION-1000)/10]==1) ) ||
2564 ((originetet == MAGIC::ORIGINE::MODIFICATION) && (etat[(MAGIC::ORIGINE::MODIFICATION-1000)/10]==1) ) ||
2565 ((originetet == MAGIC::ORIGINE::DUPLIQUER) && (etat[(MAGIC::ORIGINE::DUPLIQUER-1000)/10]==1) ) )
2566 {
2567 MG_NOEUD* noeud1 = mgtet->get_noeud1();
2568 MG_NOEUD* noeud2 = mgtet->get_noeud2();
2569 MG_NOEUD* noeud3 = mgtet->get_noeud3();
2570 MG_NOEUD* noeud4 = mgtet->get_noeud4();
2571 MG_NOEUD* node1 = mai2->get_mg_noeudid(noeud1->get_nouveau_numero());
2572 MG_NOEUD* node2 = mai2->get_mg_noeudid(noeud2->get_nouveau_numero());
2573 MG_NOEUD* node3 = mai2->get_mg_noeudid(noeud3->get_nouveau_numero());
2574 MG_NOEUD* node4 = mai2->get_mg_noeudid(noeud4->get_nouveau_numero());
2575 MG_TRIANGLE* tri1=mgtet->get_triangle1();
2576 MG_TRIANGLE* tri2=mgtet->get_triangle2();
2577 MG_TRIANGLE* tri3=mgtet->get_triangle3();
2578 MG_TRIANGLE* tri4=mgtet->get_triangle4();
2579 if (tri1->get_nouveau_numero()==1)
2580 {
2581 int origine=MAGIC::ORIGINE::TRIANGULATION;
2582 if (tri1->get_origine()==MAGIC::ORIGINE::IMPOSE) origine=MAGIC::ORIGINE::IMPOSE;
2583 MG_FACE* topo=facenoex;
2584 if (tri1->get_lien_topologie()->get_dimension()==2) topo=faceex;
2585 std::map<unsigned long,unsigned long>::iterator it=correspondid.find(tri1->get_lien_topologie()->get_id());
2586 if (it!=correspondid.end())
2587 topo=geo->get_mg_faceid(it->second);
2588 MG_TRIANGLE_PEAU* tripeau = insere_triangle(geo,mg_mai,noeud1,noeud2,noeud3,correspondid,topo,node1,node2,node3,mai2,origine);
2589 tripeau->change_nouveau_numero(0);
2590 tri1->change_nouveau_numero(0);
2591 }
2592 if (tri2->get_nouveau_numero()==1)
2593 {
2594 int origine=MAGIC::ORIGINE::TRIANGULATION;
2595 if (tri2->get_origine()==MAGIC::ORIGINE::IMPOSE) origine=MAGIC::ORIGINE::IMPOSE;
2596 MG_FACE* topo=facenoex;
2597 if (tri2->get_lien_topologie()->get_dimension()==2) topo=faceex;
2598 std::map<unsigned long,unsigned long>::iterator it=correspondid.find(tri2->get_lien_topologie()->get_id());
2599 if (it!=correspondid.end())
2600 topo=geo->get_mg_faceid(it->second);
2601 MG_TRIANGLE_PEAU* tripeau = insere_triangle(geo,mg_mai,noeud1,noeud4,noeud2,correspondid,topo,node1,node4,node2,mai2,origine);
2602 tripeau->change_nouveau_numero(0);
2603 tri2->change_nouveau_numero(0);
2604 }
2605 if (tri3->get_nouveau_numero()==1)
2606 {
2607 int origine=MAGIC::ORIGINE::TRIANGULATION;
2608 MG_FACE* topo=facenoex;
2609 if (tri3->get_lien_topologie()->get_dimension()==2) topo=faceex;
2610 if (tri3->get_origine()==MAGIC::ORIGINE::IMPOSE) origine=MAGIC::ORIGINE::IMPOSE;
2611 std::map<unsigned long,unsigned long>::iterator it=correspondid.find(tri3->get_lien_topologie()->get_id());
2612 if (it!=correspondid.end())
2613 topo=geo->get_mg_faceid(it->second);
2614 MG_TRIANGLE_PEAU* tripeau = insere_triangle(geo,mg_mai,noeud2,noeud4,noeud3,correspondid,topo,node2,node4,node3,mai2,origine);
2615 tripeau->change_nouveau_numero(0);
2616 tri3->change_nouveau_numero(0);
2617 }
2618 if (tri4->get_nouveau_numero()==1)
2619 {
2620 int origine=MAGIC::ORIGINE::TRIANGULATION;
2621 MG_FACE* topo=facenoex;
2622 if (tri4->get_lien_topologie()->get_dimension()==2) topo=faceex;
2623 if (tri4->get_origine()==MAGIC::ORIGINE::IMPOSE) origine=MAGIC::ORIGINE::IMPOSE;
2624 std::map<unsigned long,unsigned long>::iterator it=correspondid.find(tri4->get_lien_topologie()->get_id());
2625 if (it!=correspondid.end())
2626 topo=geo->get_mg_faceid(it->second);
2627 MG_TRIANGLE_PEAU* tripeau = insere_triangle(geo,mg_mai,noeud1,noeud3,noeud4,correspondid,topo,node1,node3,node4,mai2,origine);
2628 tripeau->change_nouveau_numero(0);
2629 tri4->change_nouveau_numero(0);
2630 }
2631 }
2632 }
2633 // etape 6 recherche des vosins
2634 for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
2635 {
2636 MG_TRIANGLE_PEAU* tripeau=(MG_TRIANGLE_PEAU*)mgtri;
2637 MG_TRIANGLE_PEAU* voisin1=recherche_voisin(tripeau->get_noeud1(),tripeau->get_noeud2(),tripeau);
2638 MG_TRIANGLE_PEAU* voisin2=recherche_voisin(tripeau->get_noeud2(),tripeau->get_noeud3(),tripeau);
2639 MG_TRIANGLE_PEAU* voisin3=recherche_voisin(tripeau->get_noeud3(),tripeau->get_noeud1(),tripeau);
2640 tripeau->change_voisin1(voisin1);
2641 tripeau->change_voisin2(voisin2);
2642 tripeau->change_voisin3(voisin3);
2643 }
2644 // etape 7 - On recherche les peaux
2645 int fin;
2646 do
2647 {
2648 fin=1;
2649 for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
2650 {
2651 MG_TRIANGLE_PEAU *tripeau=(MG_TRIANGLE_PEAU*)mgtri;
2652 if (tripeau->get_nouveau_numero()==0)
2653 {
2654 fin=0;
2655 std::vector<MG_TRIANGLE_PEAU*> *peau=new std::vector<MG_TRIANGLE_PEAU*>;
2656 lst_peau.push_back(peau);
2657 tripeau->change_nouveau_numero(1);
2658 peau->push_back(tripeau);
2659 determine_peau(peau);
2660 }
2661 }
2662 }
2663 while (fin==0);
2664 nbpeau=lst_peau.size();
2665 //test de manifold
2666 LISTE_MG_SEGMENT::iterator itseg;
2667 for (MG_SEGMENT* seg=mai2->get_premier_segment(itseg);seg!=NULL;seg=mai2->get_suivant_segment(itseg))
2668 seg->change_nouveau_numero(0);
2669 TPL_MAP_ENTITE<MG_SEGMENT*> nonmanifoldarete;
2670 TPL_MAP_ENTITE<MG_NOEUD*> nonmanifoldnoeud;
2671 TPL_MAP_ENTITE<MG_NOEUD*> nonmanifoldnoeuddepuisarete;
2672 nbpeau=lst_peau.size();
2673 for (int i=0;i<nbpeau;i++)
2674 {
2675 int nbtri=lst_peau[i]->size();
2676 for (int j=0;j<nbtri;j++)
2677 {
2678 MG_TRIANGLE_PEAU* tripeau=(*lst_peau[i])[j];
2679 tripeau->get_segment1()->change_nouveau_numero(tripeau->get_segment1()->get_nouveau_numero()+1);
2680 tripeau->get_segment2()->change_nouveau_numero(tripeau->get_segment2()->get_nouveau_numero()+1);
2681 tripeau->get_segment3()->change_nouveau_numero(tripeau->get_segment3()->get_nouveau_numero()+1);
2682 if (tripeau->get_segment1()->get_nouveau_numero()>2)
2683 nonmanifoldarete.ajouter(tripeau->get_segment1());
2684 if (tripeau->get_segment2()->get_nouveau_numero()>2)
2685 nonmanifoldarete.ajouter(tripeau->get_segment2());
2686 if (tripeau->get_segment3()->get_nouveau_numero()>2)
2687 nonmanifoldarete.ajouter(tripeau->get_segment3());
2688 }
2689 }
2690
2691 int nbnonmanifoldarete=nonmanifoldarete.get_nb();
2692 for (int i=0;i<nbnonmanifoldarete;i++)
2693 {
2694 MG_SEGMENT* seg=nonmanifoldarete.get(i);
2695 MG_NOEUD* n1=seg->get_noeud1();
2696 MG_NOEUD* n2=seg->get_noeud2();
2697 nonmanifoldnoeuddepuisarete.ajouter(n1);
2698 nonmanifoldnoeuddepuisarete.ajouter(n2);
2699 }
2700 LISTE_MG_NOEUD::iterator itnoeud;
2701 for (MG_NOEUD* no=mai2->get_premier_noeud(itnoeud);no!=NULL;no=mai2->get_suivant_noeud(itnoeud))
2702 {
2703 if (nonmanifoldnoeuddepuisarete.existe(no)) continue;
2704 if (est_non_manifold(no)) nonmanifoldnoeud.ajouter(no);
2705 }
2706
2707 nbmaniare = nonmanifoldarete.get_nb();
2708 nbmanino = nonmanifoldnoeud.get_nb();
2709
2710 //etape 8 - Traitement des cas de non manifold
2711 //non manifold par arete
2712 for (int i=0;i<nbnonmanifoldarete;i++)
2713 {
2714 MG_SEGMENT* segment=nonmanifoldarete.get(i);
2715 MG_NOEUD* noeud1 = mg_mai->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
2716 MG_NOEUD* noeud2 = mg_mai->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
2717 int nb_tetra = noeud1->get_lien_tetra()->get_nb();
2718 for (int j=0;j<nb_tetra;j++)
2719 {
2720 MG_TETRA* mgtet =noeud1->get_lien_tetra()->get(j);
2721 int originetet=mgtet->get_origine();
2722 if (originetet == MAGIC::ORIGINE::MAILLEUR_AUTO)
2723 {
2724 //On reactive le tetra si l'autre noeud du segment lui appartient aussi
2725 MG_NOEUD* no1 = mgtet->get_noeud1();
2726 MG_NOEUD* no2 = mgtet->get_noeud2();
2727 MG_NOEUD* no3 = mgtet->get_noeud3();
2728 MG_NOEUD* no4 = mgtet->get_noeud4();
2729 if((no1 == noeud2) ||(no2 == noeud2) || (no3 == noeud2) || (no4 == noeud2))
2730 mgtet->change_origine(MAGIC::ORIGINE::OPTIMISE);
2731 }
2732 }
2733 }
2734
2735 //non manifold par noeud
2736 int nbnonmanifoldnoeud=nonmanifoldnoeud.get_nb();
2737 for (int i=0;i<nbnonmanifoldnoeud;i++)
2738 {
2739 MG_NOEUD* noeud=mg_mai->get_mg_noeudid(nonmanifoldnoeud.get(i)->get_nouveau_numero());
2740 int nb_tetra = noeud->get_lien_tetra()->get_nb();
2741 for (int j = 0;j<nb_tetra;j++)
2742 {
2743 MG_TETRA* mgtet =noeud->get_lien_tetra()->get(j);
2744 int originetet=mgtet->get_origine();
2745 if (originetet == MAGIC::ORIGINE::MAILLEUR_AUTO)
2746 {
2747 mgtet->change_origine(MAGIC::ORIGINE::OPTIMISE);
2748 }
2749 }
2750 }
2751 *mai2_id = mai2->get_id();
2752 if ((nbnonmanifoldarete != 0) || (nbnonmanifoldnoeud != 0))
2753 {
2754 for (int i=0;i<lst_peau.size();i++)
2755 {
2756 delete lst_peau[i];
2757 }
2758 lst_peau.clear();
2759 gest2.supprimer_mg_maillageid(*mai2_id);
2760 return 0;
2761 }
2762
2763 return 1;
2764 }
2765
2766
2767 int MGOPT_POSTTRAITEMENT::cardinalite(MG_MAILLAGE* mgmai, int cardinalite_iter)
2768 {
2769 for (int i=0;i<cardinalite_iter;i++)
2770 {
2771 TPL_MAP_ENTITE<MG_NOEUD*> lst;
2772 LISTE_MG_NOEUD::iterator it;
2773 for (MG_NOEUD* no=mgmai->get_premier_noeud(it);no!=NULL;no=mgmai->get_suivant_noeud(it))
2774 if (no->get_lien_triangle()->get_nb()==3)
2775 {
2776 if (no->get_lien_triangle()->get(0)->get_origine()!=MAGIC::ORIGINE::IMPOSE)
2777 if (no->get_lien_triangle()->get(1)->get_origine()!=MAGIC::ORIGINE::IMPOSE)
2778 if (no->get_lien_triangle()->get(2)->get_origine()!=MAGIC::ORIGINE::IMPOSE)
2779 lst.ajouter(no);
2780 }
2781 int nbcas=lst.get_nb();
2782 char message[255];
2783 sprintf(message," Itération %d, nombre de noeuds de cardinalite 3 : %d",i+1,nbcas);
2784 affiche(message);
2785 if (nbcas==0) return i;
2786 for (int j=0;j<nbcas;j++)
2787 {
2788 MG_NOEUD* no=lst.get(j);
2789 TPL_MAP_ENTITE<MG_NOEUD*> lstno;
2790 lstno.ajouter(no->get_lien_triangle()->get(0)->get_noeud1());
2791 lstno.ajouter(no->get_lien_triangle()->get(0)->get_noeud2());
2792 lstno.ajouter(no->get_lien_triangle()->get(0)->get_noeud3());
2793 lstno.ajouter(no->get_lien_triangle()->get(1)->get_noeud1());
2794 lstno.ajouter(no->get_lien_triangle()->get(1)->get_noeud2());
2795 lstno.ajouter(no->get_lien_triangle()->get(1)->get_noeud3());
2796 lstno.ajouter(no->get_lien_triangle()->get(2)->get_noeud1());
2797 lstno.ajouter(no->get_lien_triangle()->get(2)->get_noeud2());
2798 lstno.ajouter(no->get_lien_triangle()->get(2)->get_noeud3());
2799 lstno.supprimer(no);
2800 MG_ELEMENT_TOPOLOGIQUE* topo=no->get_lien_triangle()->get(0)->get_lien_topologie();
2801 int origine=no->get_lien_triangle()->get(0)->get_origine();
2802 OT_VECTEUR_3D normal;
2803 for (int k=0;k<3;k++)
2804 {
2805 MG_TRIANGLE* tri=no->get_lien_triangle()->get(k);
2806 MG_NOEUD *no1=tri->get_noeud1();
2807 MG_NOEUD *no2=tri->get_noeud2();
2808 MG_NOEUD *no3=tri->get_noeud3();
2809 OT_VECTEUR_3D vec1(no1->get_coord(),no2->get_coord());
2810 OT_VECTEUR_3D vec2(no1->get_coord(),no3->get_coord());
2811 OT_VECTEUR_3D vec3=vec1&vec2;
2812 vec3.norme();
2813 normal=normal+vec3;
2814 }
2815 normal.norme();
2816 mgmai->supprimer_mg_triangleid(no->get_lien_triangle()->get(2)->get_id());
2817 mgmai->supprimer_mg_triangleid(no->get_lien_triangle()->get(1)->get_id());
2818 mgmai->supprimer_mg_triangleid(no->get_lien_triangle()->get(0)->get_id());
2819 mgmai->supprimer_mg_noeudid(no->get_id());
2820 MG_NOEUD *no1=lstno.get(0);
2821 MG_NOEUD *no2=lstno.get(1);
2822 MG_NOEUD *no3=lstno.get(2);
2823 OT_VECTEUR_3D vec1(no1->get_coord(),no2->get_coord());
2824 OT_VECTEUR_3D vec2(no1->get_coord(),no3->get_coord());
2825 OT_VECTEUR_3D vec3=vec1&vec2;
2826 vec3.norme();
2827 if (vec3*normal>0) insere_triangle(topo,no1,no2,no3,mgmai,origine);
2828 else insere_triangle(topo,no2,no1,no3,mgmai,origine);
2829 }
2830
2831 }
2832
2833 return cardinalite_iter;
2834
2835 }
2836
2837
2838
2839 void MGOPT_POSTTRAITEMENT::rmimpose(MG_MAILLAGE* mgmai, char* nom)
2840 {
2841 MG_GESTIONNAIRE gest2;
2842 MG_MAILLAGE* mgmai3=new MG_MAILLAGE(NULL);
2843 gest2.ajouter_mg_maillage(mgmai3);
2844 TPL_MAP_ENTITE<MG_TRIANGLE*> lsttri;
2845 TPL_MAP_ENTITE<MG_TRIANGLE*>::ITERATEUR it_lsttri;
2846 TPL_MAP_ENTITE<MG_NOEUD*> lstno;
2847 TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR it_lstno;
2848
2849 LISTE_MG_TRIANGLE::iterator itt;
2850 LISTE_MG_NOEUD::iterator itn;
2851 for (MG_TRIANGLE *tri=mgmai->get_premier_triangle(itt);tri!=NULL;tri=mgmai->get_suivant_triangle(itt))
2852 {
2853 if (tri->get_origine() != MAGIC::ORIGINE::IMPOSE)
2854 {
2855 lsttri.ajouter(tri);
2856 lstno.ajouter(tri->get_noeud1());
2857 lstno.ajouter(tri->get_noeud2());
2858 lstno.ajouter(tri->get_noeud3());
2859 }
2860 }
2861 // Creation des noeuds du non_design dans le nouveau maillage mg_mai4
2862 for(MG_NOEUD* mgnoeud=lstno.get_premier(it_lstno);mgnoeud;mgnoeud=lstno.get_suivant(it_lstno))
2863 {
2864 double x = mgnoeud->get_x();
2865 double y = mgnoeud->get_y();
2866 double z = mgnoeud->get_z();
2867 MG_NOEUD* node = new MG_NOEUD(NULL,x,y,z,MAGIC::ORIGINE::TRIANGULATION);
2868 node->change_nouveau_numero(mgnoeud->get_id());
2869 mgmai3->ajouter_mg_noeud(node);
2870 }
2871 for(MG_TRIANGLE* mgtri=lsttri.get_premier(it_lsttri);mgtri;mgtri=lsttri.get_suivant(it_lsttri))
2872 {
2873 MG_NOEUD* noeud1 = mgtri->get_noeud1();
2874 MG_NOEUD* noeud2 = mgtri->get_noeud2();
2875 MG_NOEUD* noeud3 = mgtri->get_noeud3();
2876 MG_NOEUD* no1 = NULL;
2877 MG_NOEUD* no2 = NULL;
2878 MG_NOEUD* no3 = NULL;
2879 for (MG_NOEUD* mgnoeud=mgmai3->get_premier_noeud(itn);mgnoeud!=NULL;mgnoeud=mgmai3->get_suivant_noeud(itn))
2880 {
2881 if (mgnoeud->get_nouveau_numero() == noeud1->get_id()) no1 = mgnoeud;
2882 if (mgnoeud->get_nouveau_numero() == noeud2->get_id()) no2 = mgnoeud;
2883 if (mgnoeud->get_nouveau_numero() == noeud3->get_id()) no3 = mgnoeud;
2884 }
2885 MG_TRIANGLE* tri = mgmai3->ajouter_mg_triangle(NULL,no1,no2,no3,mgtri->get_origine());
2886 }
2887
2888 gest2.enregistrer(nom);
2889 }
2890
2891 //---------------------------- Determine peau -------------------------------------------------------------------------
2892 void MGOPT_POSTTRAITEMENT::determine_peau(std::vector<MG_TRIANGLE_PEAU*> * peau)
2893 {
2894 int num=0;
2895 while (num!=peau->size())
2896 {
2897 MG_TRIANGLE_PEAU* p=(*peau)[num];
2898 if (p->get_voisin1()->get_nouveau_numero()==0)
2899 {
2900 peau->push_back(p->get_voisin1());
2901 p->get_voisin1()->change_nouveau_numero(1);
2902 }
2903 if (p->get_voisin2()->get_nouveau_numero()==0)
2904 {
2905 peau->push_back(p->get_voisin2());
2906 p->get_voisin2()->change_nouveau_numero(1);
2907 }
2908 if (p->get_voisin3()->get_nouveau_numero()==0)
2909 {
2910 peau->push_back(p->get_voisin3());
2911 p->get_voisin3()->change_nouveau_numero(1);
2912 }
2913 num++;
2914 }
2915 }
2916
2917 //---------------------------- Insere triangle -------------------------------------------------------------------------
2918 MG_TRIANGLE_PEAU* MGOPT_POSTTRAITEMENT::insere_triangle(MG_GEOMETRIE* geo,MG_MAILLAGE* maiori,MG_NOEUD* n1ori,MG_NOEUD* n2ori,MG_NOEUD* n3ori,std::map<unsigned long,unsigned long> &correspondid,MG_ELEMENT_TOPOLOGIQUE* topo,class MG_NOEUD *mgnoeud1,class MG_NOEUD *mgnoeud2,class MG_NOEUD *mgnoeud3,MG_MAILLAGE* mg_maillage,int origine)
2919 {
2920 //MG_TRIANGLE_PEAU* voisin1=NULL,*voisin2=NULL,*voisin3=NULL;
2921 MG_SEGMENT* mgsegment1=mg_maillage->get_mg_segment(mgnoeud1->get_id(),mgnoeud2->get_id());
2922 MG_SEGMENT* mgsegment2=mg_maillage->get_mg_segment(mgnoeud2->get_id(),mgnoeud3->get_id());
2923 MG_SEGMENT* mgsegment3=mg_maillage->get_mg_segment(mgnoeud3->get_id(),mgnoeud1->get_id());
2924 if (mgsegment1==NULL)
2925 {
2926 MG_SEGMENT* seg=maiori->get_mg_segment(n1ori->get_id(),n2ori->get_id());
2927 MG_ELEMENT_TOPOLOGIQUE* toposeg=seg->get_lien_topologie();
2928 MG_ELEMENT_TOPOLOGIQUE* toposegn=topo;
2929 std::map<unsigned long,unsigned long>::iterator it=correspondid.find(toposeg->get_id());
2930 if (it!=correspondid.end())
2931 toposegn=geo->get_mg_areteid(it->second);
2932 mgsegment1=mg_maillage->ajouter_mg_segment(toposegn,mgnoeud1,mgnoeud2,origine);
2933 }
2934 if (mgsegment2==NULL)
2935 {
2936 MG_SEGMENT* seg=maiori->get_mg_segment(n2ori->get_id(),n3ori->get_id());
2937 MG_ELEMENT_TOPOLOGIQUE* toposeg=seg->get_lien_topologie();
2938 MG_ELEMENT_TOPOLOGIQUE* toposegn=topo;
2939 std::map<unsigned long,unsigned long>::iterator it=correspondid.find(toposeg->get_id());
2940 if (it!=correspondid.end())
2941 toposegn=geo->get_mg_areteid(it->second);
2942 mgsegment2=mg_maillage->ajouter_mg_segment(toposegn,mgnoeud2,mgnoeud3,origine);
2943 }
2944 if (mgsegment3==NULL)
2945 {
2946 MG_SEGMENT* seg=maiori->get_mg_segment(n3ori->get_id(),n1ori->get_id());
2947 MG_ELEMENT_TOPOLOGIQUE* toposeg=seg->get_lien_topologie();
2948 MG_ELEMENT_TOPOLOGIQUE* toposegn=topo;
2949 std::map<unsigned long,unsigned long>::iterator it=correspondid.find(toposeg->get_id());
2950 if (it!=correspondid.end())
2951 toposegn=geo->get_mg_areteid(it->second);
2952 mgsegment3=mg_maillage->ajouter_mg_segment(toposegn,mgnoeud3,mgnoeud1,origine);
2953 }
2954 MG_TRIANGLE_PEAU* mtriangle=new MG_TRIANGLE_PEAU(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,origine);
2955 mg_maillage->ajouter_mg_triangle(mtriangle);
2956 return mtriangle;
2957 }
2958
2959 MG_TRIANGLE_PEAU* MGOPT_POSTTRAITEMENT::insere_triangle(MG_ELEMENT_TOPOLOGIQUE* topo,class MG_NOEUD *mgnoeud1,class MG_NOEUD *mgnoeud2,class MG_NOEUD *mgnoeud3,MG_MAILLAGE* mg_maillage,int origine)
2960 {
2961 //MG_TRIANGLE_PEAU* voisin1=NULL,*voisin2=NULL,*voisin3=NULL;
2962 MG_SEGMENT* mgsegment1=mg_maillage->get_mg_segment(mgnoeud1->get_id(),mgnoeud2->get_id());
2963 MG_SEGMENT* mgsegment2=mg_maillage->get_mg_segment(mgnoeud2->get_id(),mgnoeud3->get_id());
2964 MG_SEGMENT* mgsegment3=mg_maillage->get_mg_segment(mgnoeud3->get_id(),mgnoeud1->get_id());
2965 if (mgsegment1==NULL)
2966 mgsegment1=mg_maillage->ajouter_mg_segment(topo,mgnoeud1,mgnoeud2,origine);
2967 if (mgsegment2==NULL)
2968 mgsegment2=mg_maillage->ajouter_mg_segment(topo,mgnoeud2,mgnoeud3,origine);
2969 if (mgsegment3==NULL)
2970 mgsegment3=mg_maillage->ajouter_mg_segment(topo,mgnoeud3,mgnoeud1,origine);
2971 MG_TRIANGLE_PEAU* mtriangle=new MG_TRIANGLE_PEAU(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,origine);
2972 mg_maillage->ajouter_mg_triangle(mtriangle);
2973 return mtriangle;
2974 }
2975 //---------------------------- Recherche voisin -------------------------------------------------------------------------
2976 MG_TRIANGLE_PEAU* MGOPT_POSTTRAITEMENT::recherche_voisin(MG_NOEUD* mg_noeud1,MG_NOEUD* mg_noeud2,MG_TRIANGLE_PEAU* triref)
2977 {
2978 MG_TRIANGLE_PEAU* trisol=NULL;
2979 double angleref=4.*M_PI;
2980 int nb1=mg_noeud1->get_lien_triangle()->get_nb();
2981 int nb2=mg_noeud2->get_lien_triangle()->get_nb();
2982 for (int i=0;i<nb1;i++)
2983 for (int j=0;j<nb2;j++)
2984 if (mg_noeud1->get_lien_triangle()->get(i)==mg_noeud2->get_lien_triangle()->get(j))
2985 if ((MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i)!=triref)
2986 {
2987 double angle=calcul_angle(triref,(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i));
2988 if (angle<angleref)
2989 {
2990 angleref=angle;
2991 trisol=(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i);
2992 }
2993 }
2994 return trisol;
2995 }
2996
2997 int MGOPT_POSTTRAITEMENT::est_non_manifold(MG_NOEUD* no)
2998 {
2999 static int compteur=0;
3000 compteur++;
3001 int nb_tri=no->get_lien_triangle()->get_nb();
3002 TPL_MAP_ENTITE<MG_TRIANGLE*> lst;
3003 for (int i=0;i<nb_tri;i++)
3004 lst.ajouter(no->get_lien_triangle()->get(i));
3005 MG_TRIANGLE_PEAU* premier_triangle=(MG_TRIANGLE_PEAU*)lst.get(0);
3006 lst.supprimer(premier_triangle);
3007 MG_TRIANGLE_PEAU* tricour=(MG_TRIANGLE_PEAU*)premier_triangle;
3008 int i=0;
3009 do
3010 {
3011 if (lst.existe(tricour->get_voisin1()) || ((tricour->get_voisin1()==premier_triangle)&& (i>1)))
3012 {
3013 tricour=tricour->get_voisin1();
3014 lst.supprimer(tricour);
3015 }
3016 else if (lst.existe(tricour->get_voisin2()) || ((tricour->get_voisin2()==premier_triangle)&& (i>1)))
3017 {
3018 tricour=tricour->get_voisin2();
3019 lst.supprimer(tricour);
3020 }
3021 else if (lst.existe(tricour->get_voisin3()) || ((tricour->get_voisin3()==premier_triangle)&& (i>1)))
3022 {
3023 tricour=tricour->get_voisin3();
3024 lst.supprimer(tricour);
3025 }
3026 i++;
3027 }
3028 while (tricour!=premier_triangle);
3029 if (lst.get_nb()>0)
3030 return 1;
3031 return 0;
3032 }
3033
3034 //---------------------------- Calcule angle -------------------------------------------------------------------------
3035 double MGOPT_POSTTRAITEMENT::calcul_angle(MG_TRIANGLE_PEAU* ft1,MG_TRIANGLE_PEAU* ft2)
3036 {
3037 MG_NOEUD* noeud1=ft1->get_noeud1();
3038 MG_NOEUD* noeud2=ft1->get_noeud2();
3039 MG_NOEUD* noeud3=ft1->get_noeud3();
3040 MG_NOEUD* noeud4=ft2->get_noeud1();
3041 MG_NOEUD* noeud5=ft2->get_noeud2();
3042 MG_NOEUD* noeud6=ft2->get_noeud3();
3043 double angle=get_angle_par_noeud<MG_NOEUD*>(noeud1,noeud2,noeud3,noeud4,noeud5,noeud6);
3044 return angle;
3045 }
3046
3047 //---------------------------- Ponderation gaussian -------------------------------------------------------------------------
3048 double MGOPT_POSTTRAITEMENT::ponderation_gaussian(double s,double sigma)
3049 {
3050 double w_s;
3051 w_s = exp(-(s*s)/(2.*sigma*sigma));
3052 return w_s;
3053 }
3054
3055 //---------------------------- Ponderation Laplacian -------------------------------------------------------------------------
3056 double MGOPT_POSTTRAITEMENT::ponderation_laplacian(double s,double sigma)
3057 {
3058 double w_s;
3059 w_s = exp(-(sqrt(2.)*fabs(s))/sigma);
3060 return w_s;
3061 }
3062
3063 //---------------------------- Ponderation elfallahford -------------------------------------------------------------------------
3064 double MGOPT_POSTTRAITEMENT::ponderation_elfallahford(double s,double sigma)
3065 {
3066 double w_s;
3067 w_s = 1./sqrt(1+pow((s/sigma),2));
3068 return w_s;
3069 }
3070
3071
3072
3073 int passage_non_manifold(MG_NOEUD** no,std::map<int,MG_NOEUD*>* lst_no_voisins, MG_GEOMETRIE* geo, std::vector<MG_NOEUD*>* arete)
3074 {
3075 MG_SEGMENT* seg=NULL;
3076 MG_NOEUD* no_voisin=NULL;
3077 std::map<int,MG_NOEUD*>::iterator it_map;
3078
3079 for(int i=0;i<(*no)->get_lien_segment()->get_nb();i++)
3080 {
3081 seg=(*no)->get_lien_segment()->get(i);
3082 if((seg->get_lien_triangle()->get(0)->get_lien_topologie()==geo->get_mg_face(geo->get_nb_mg_face()-2))&&
3083 (seg->get_lien_triangle()->get(1)->get_lien_topologie()==geo->get_mg_face(geo->get_nb_mg_face()-2)))
3084 {
3085 if(seg->get_noeud1()!=*no) (*lst_no_voisins)[i]=seg->get_noeud1();
3086 else (*lst_no_voisins)[i]=seg->get_noeud2();
3087 }
3088 }
3089 no_voisin=(*arete)[arete->size()-2];
3090 int k,j=0;
3091 for(int i=0;i<lst_no_voisins->size();i++)
3092 {
3093 bool find=false;
3094 while((j<no_voisin->get_lien_segment()->get_nb())&&!find)
3095 {
3096 seg=no_voisin->get_lien_segment()->get(j);
3097 bool seg_valide=(seg->get_lien_triangle()->get(0)->get_lien_topologie()==geo->get_mg_face(geo->get_nb_mg_face()-2))||
3098 (seg->get_lien_triangle()->get(1)->get_lien_topologie()==geo->get_mg_face(geo->get_nb_mg_face()-2));
3099 j++;
3100 for(it_map=lst_no_voisins->begin();it_map!=lst_no_voisins->end();it_map++)
3101 if((it_map->second->get_lien_segment()->est_dans_la_liste(seg))&&seg_valide)
3102 {
3103 no_voisin=it_map->second;
3104 k=it_map->first;
3105 lst_no_voisins->erase(it_map);
3106 j=0;
3107 find=true;
3108 break;
3109 }
3110 }
3111 }
3112 arete->push_back(no_voisin);
3113 *no=no_voisin;
3114 return k;
3115 }
3116
3117
3118
3119 void recherche_voisins(MG_NOEUD* no,std::vector<MG_NOEUD*>* lst_no_1090_1091,std::vector<MG_NOEUD*>* arete,std::map<int,MG_NOEUD*>* lst_no_voisins,bool arete_ouverte,std::vector<MG_NOEUD*>* lst_no_extr)
3120 {
3121 MG_SEGMENT* seg=NULL;
3122 MG_NOEUD* no_voisin=NULL;
3123
3124 for(int j=0;j<no->get_lien_segment()->get_nb();j++)
3125 {
3126 seg=no->get_lien_segment()->get(j);
3127 bool seg_bord=seg->get_lien_triangle()->get(0)->get_lien_topologie()!=seg->get_lien_triangle()->get(1)->get_lien_topologie();
3128 if(seg->get_noeud1()!=no) no_voisin=seg->get_noeud1();
3129 else no_voisin=seg->get_noeud2();
3130 for(int k=0;k<lst_no_1090_1091->size();k++)
3131 if(((*lst_no_1090_1091)[k]==no_voisin)&&seg_bord&&(no_voisin!=(*arete)[arete->size()-2]))
3132 {
3133 (*lst_no_voisins)[k]=no_voisin;
3134 break;
3135 }
3136 if(arete_ouverte&&(no_voisin->get_origine()==MAGIC::ORIGINE::IMPOSE)&&(arete->size()>2)&&seg_bord)
3137 for(int k=0;k<lst_no_extr->size();k++)
3138 if(((*lst_no_extr)[k]==no_voisin))
3139 {
3140 arete->push_back(no_voisin);
3141 lst_no_extr->erase(lst_no_extr->begin()+k);
3142 break;
3143 }
3144 }
3145 }