ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/aster/src/mgopt_posttraitement.cpp
Revision: 831
Committed: Sat Oct 1 23:33:03 2016 UTC (8 years, 7 months ago) by francois
File size: 112553 byte(s)
Log Message: 
ajout de la statistique  f reel dans l analyse des volumes

File Contents

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