ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/aster/src/mgopt_posttraitement.cpp
Revision: 905
Committed: Tue Sep 19 18:58:56 2017 UTC (7 years, 11 months ago) by amroune
File size: 129539 byte(s)
Log Message:
ajout d'une méthode de lissage des arètes dans le post traitement de TO

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