ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mg_lissage.cpp
Revision: 457
Committed: Thu Nov 7 00:47:01 2013 UTC (11 years, 6 months ago) by francois
File size: 65996 byte(s)
Log Message:
Nouvelle interface pour le post traitement des resultats de topo_optis (anciennement appele lissage et maintenant post traitement) 

File Contents

# User Rev Content
1 francois 222 #include "gestionversion.h"
2     #include "mg_lissage.h"
3     #include "fem_solution.h"
4     #include "fem_maillage.h"
5     #include "mg_maillage.h"
6     #include "mg_gestionnaire.h"
7     #include "mg_triangle_peau.h"
8 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 343 #include "mg_export.h"
13 francois 222 #include <stdio.h>
14 picher 230 #include <stdlib.h>
15     #include <time.h>
16 francois 222 #include <string.h>
17 francois 224 #include <math.h>
18 francois 222 //---------------------------------------------------------------------------
19    
20    
21 francois 224
22 francois 232 MG_LISSAGE::MG_LISSAGE():affichageactif(0)
23 francois 222 {
24 francois 224 for (int i=0;i<9;i++) etat[i]=0;
25 francois 457 params.ajouter("decoupage",0.,OT_PARAMETRES::DOUBLE,"0. découpage en conservant les mailles entieres 1. découpage par isodensite");
26     params.ajouter("seuil",0.8,OT_PARAMETRES::DOUBLE,"Valeur du seuil de densité pour les 2 découpages");
27     params.ajouter("nomdensextra","densiteextra.sol",OT_PARAMETRES::STRING,"Nom du fichier comportant la solution de la densité extrapolé au noeud");
28     params.ajouter("consimpose",1.,OT_PARAMETRES::DOUBLE,"0. Maillage non design non inclus dans le decoupage par mailles entieres 1. Maillage design inclus");
29     params.ajouter("consoptimise",1.,OT_PARAMETRES::DOUBLE,"0. Maillage optimise non inclus dans le decoupage par mailles entieres 1. Maillage optimise inclus");
30     params.ajouter("consnooptimise",0.,OT_PARAMETRES::DOUBLE,"0. Maillage non optimise non inclus dans le decoupage par mailles entieres 1. Maillage non optimise inclus");
31     params.ajouter("reactivation",1.,OT_PARAMETRES::DOUBLE,"1. Reactivation de mailles isolées sur le long des aretes 0. non reactivation");
32     params.ajouter("chen2005_debut",0.,OT_PARAMETRES::DOUBLE,"Numero d'itération du debut d'application deu lissage de chen2005 0. méthode non utilisee");
33     params.ajouter("chen2005_itermax",10.,OT_PARAMETRES::DOUBLE,"Nombre d'itérations maxiamles du lissage de chen2005");
34     params.ajouter("chen2005_epsilon",0.01,OT_PARAMETRES::DOUBLE,"Critere de convergence du lissage de chen2005");
35     params.ajouter("chen2005_sigma",0.1,OT_PARAMETRES::DOUBLE,"Écart type du filtre gaussien du lissage de chen2005");
36     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");
37     params.ajouter("chen2008_itermax",10.,OT_PARAMETRES::DOUBLE,"Nombre d'itérations maxiamles du lissage de chen2008");
38     params.ajouter("chen2008_epsilon",0.01,OT_PARAMETRES::DOUBLE,"Critere de convergence du lissage de chen2008");
39     params.ajouter("chen2008_sigma",0.1,OT_PARAMETRES::DOUBLE,"Écart type du filtre gaussien du lissage de chen2008");
40     params.ajouter("chen2008_gamma",1.,OT_PARAMETRES::DOUBLE,"Vitesse d'avancement de déplacement du lissage de chen2008");
41     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");
42     params.ajouter("jiao2012_itermax",10.,OT_PARAMETRES::DOUBLE,"Nombre d'itérations maxiamles du lissage de jioa2012");
43 francois 222 }
44    
45    
46     MG_LISSAGE::~MG_LISSAGE()
47     {
48 francois 224 for (int i=0;i<lst_peau.size();i++)
49     delete lst_peau[i];
50 francois 222 }
51 picher 233
52 francois 232 void MG_LISSAGE::active_affichage(void (*fonc)(char*))
53     {
54     affiche=fonc;
55     affichageactif=1;
56     }
57 francois 222
58 picher 233
59 francois 457
60     void MG_LISSAGE::lisse(char *fichier)
61     {
62     affiche((char*)"");
63     affiche((char*)"***********************************************************");
64     affiche((char*)"Post-traitement d'un calcul d'optimisation de topologie");
65     affiche((char*)"***********************************************************");
66     affiche((char*)"");
67     affiche((char*)"");
68     affiche((char*)"");
69     affiche((char*)"Écriture d'un fichier de parametres par défaut");
70     params.enregistrer(fichier);
71     }
72    
73    
74    
75 francois 224 void MG_LISSAGE::conserve(int origine)
76     {
77     etat[(origine-1000)/10]=1;
78     }
79    
80 picher 233
81 francois 457 void MG_LISSAGE::copieorigine(FEM_MAILLAGE* mai)
82 picher 231 {
83 francois 309 LISTE_FEM_ELEMENT3::iterator it_tetra;
84     for (FEM_ELEMENT3* tet=mai->get_premier_element3(it_tetra);tet!=NULL;tet=mai->get_suivant_element3(it_tetra))
85 picher 233 {
86     MG_TETRA* mgtet=(MG_TETRA*)tet->get_mg_element_maillage();
87     mgtet->change_origine(mgtet->get_origine());
88     }
89     }
90 francois 457
91    
92 picher 233 void MG_LISSAGE::gain_poids(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2)
93     {
94     LISTE_MG_TETRA::iterator it_tetra;
95     double volume_initial = 0.;
96     double volume_final = 0.;
97     for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
98     {
99     int originetet=mgtet->get_origine();
100     MG_NOEUD* noeud1 = mgtet->get_noeud1();
101     double x1 = noeud1->get_x();
102     double y1 = noeud1->get_y();
103     double z1 = noeud1->get_z();
104     MG_NOEUD* noeud2 = mgtet->get_noeud2();
105     double x2 = noeud2->get_x();
106     double y2 = noeud2->get_y();
107     double z2 = noeud2->get_z();
108     MG_NOEUD* noeud3 = mgtet->get_noeud3();
109     double x3 = noeud3->get_x();
110     double y3 = noeud3->get_y();
111     double z3 = noeud3->get_z();
112     MG_NOEUD* noeud4 = mgtet->get_noeud4();
113     double x4 = noeud4->get_x();
114     double y4 = noeud4->get_y();
115     double z4 = noeud4->get_z();
116     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));
117     volume_tet = volume_tet/6.;
118     volume_initial = volume_initial + volume_tet;
119     if (originetet != MAILLEUR_AUTO)
120     {
121     volume_final = volume_final + volume_tet;
122     }
123     }
124     }
125 francois 457
126 picher 233 void MG_LISSAGE::reactivation(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2)
127     {
128 francois 457 cout << " Reactivation d'elements manquants sur aretes vives et faces planes" << endl;
129 picher 233 LISTE_MG_TETRA::iterator it_tetra;
130     for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
131     {
132     int originetet=mgtet->get_origine();
133     int compteur = 0;
134     if (originetet == MAILLEUR_AUTO)
135     {
136     MG_TRIANGLE* tri1 = mgtet->get_triangle1();
137     MG_TRIANGLE* tri2 = mgtet->get_triangle2();
138     MG_TRIANGLE* tri3 = mgtet->get_triangle3();
139     MG_TRIANGLE* tri4 = mgtet->get_triangle4();
140     int nb_tet1 = tri1->get_lien_tetra()->get_nb();
141     for (int k = 0;k<nb_tet1;k++)
142     {
143     MG_TETRA* tet1 = tri1->get_lien_tetra()->get(k);
144     int ori1 = tet1->get_origine();
145     if ((tet1 != mgtet) && (ori1 == MAILLEUR_AUTO)) compteur++;
146     }
147     int nb_tet2 = tri2->get_lien_tetra()->get_nb();
148     for (int k = 0;k<nb_tet2;k++)
149     {
150     MG_TETRA* tet2 = tri2->get_lien_tetra()->get(k);
151     int ori2 = tet2->get_origine();
152     if ((tet2 != mgtet) && (ori2 == MAILLEUR_AUTO)) compteur++;
153     }
154     int nb_tet3 = tri3->get_lien_tetra()->get_nb();
155     for (int k = 0;k<nb_tet3;k++)
156     {
157     MG_TETRA* tet3 = tri3->get_lien_tetra()->get(k);
158     int ori3 = tet3->get_origine();
159     if ((tet3 != mgtet) && (ori3 == MAILLEUR_AUTO)) compteur++;
160     }
161     int nb_tet4 = tri4->get_lien_tetra()->get_nb();
162     for (int k = 0;k<nb_tet4;k++)
163     {
164     MG_TETRA* tet4 = tri4->get_lien_tetra()->get(k);
165     int ori4 = tet4->get_origine();
166     if ((tet4 != mgtet) && (ori4 == MAILLEUR_AUTO)) compteur++;
167     }
168 picher 248 if (compteur == 0) mgtet->change_origine(OPTIMISE);
169 picher 233 }
170     }
171     }
172    
173    
174 francois 343
175    
176 francois 457 MG_MAILLAGE* MG_LISSAGE::extract_skin_par_decoupage(FEM_SOLUTION* sol,double limit,MG_GESTIONNAIRE& gest2,std::string nom)
177 francois 343 {
178     affiche((char*)"Extraction de l'enveloppe par découpage");
179     sol->active_solution(0);
180     FEM_MAILLAGE* mai=sol->get_maillage();
181     affiche((char*)" Extrapolation de la densité aux noeuds");
182     LISTE_FEM_NOEUD::iterator it;
183     for (FEM_NOEUD* no=mai->get_premier_noeud(it);no!=NULL;no=mai->get_suivant_noeud(it))
184     {
185     double nume=0.;
186     double deno=0.;
187     int nb=no->get_lien_element3()->get_nb();
188     int passe=0;
189     for (int i=0;i<nb;i++)
190     {
191     FEM_ELEMENT3* tet=no->get_lien_element3()->get(i);
192     if (tet->get_mg_element_maillage()->get_origine()!=IMPOSE)
193     {
194     double jac[9];
195     int col,ligne;
196     double volume=tet->get_jacobien(jac,jac,ligne,col,1.);
197     passe=1;
198     nume=nume+tet->get_solution()*volume;
199     deno=deno+volume;
200     }
201     }
202     if (passe==1) no->change_solution(nume/deno); else no->change_solution(1.);
203     }
204     if (nom!="")
205     {
206     affiche((char*)" Enregistrement de la densité aux noeuds");
207     MG_GESTIONNAIRE *gest=mai->get_mg_maillage()->get_gestionnaire();
208     std::string chemin=nom+".sol";
209 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);
210 francois 343 gest->ajouter_fem_solution(femdens);
211     femdens->change_legende(0,"Densite");
212     LISTE_FEM_ELEMENT3::iterator it3;
213     int i=0;
214     for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
215     {
216     MG_TETRA* tet2=(MG_TETRA*)tet->get_mg_element_maillage();
217     if (tet2->get_origine()==IMPOSE)
218     {
219 francois 375 femdens->ecrire(1.,i,0,0);
220     femdens->ecrire(1.,i,0,1);
221     femdens->ecrire(1.,i,0,2);
222     femdens->ecrire(1.,i,0,3);
223 francois 343 }
224     else
225     {
226 francois 375 femdens->ecrire(tet->get_fem_noeud(0)->get_solution(),i,0,0);
227     femdens->ecrire(tet->get_fem_noeud(1)->get_solution(),i,0,1);
228     femdens->ecrire(tet->get_fem_noeud(2)->get_solution(),i,0,2);
229     femdens->ecrire(tet->get_fem_noeud(3)->get_solution(),i,0,3);
230 francois 343 }
231     i++;
232     }
233     MG_EXPORT exp;
234     exp.gmsh(mai,nom);
235     gest->supprimer_fem_solution(gest->get_nb_fem_solution()-1);
236     }
237     MG_MAILLAGE* mgmai=new MG_MAILLAGE(NULL);
238     gest2.ajouter_mg_maillage(mgmai);
239    
240     MG_MAILLAGE* mgfem=mai->get_mg_maillage();
241     LISTE_MG_TRIANGLE::iterator itmg;
242     for (MG_TRIANGLE* tri=mgfem->get_premier_triangle(itmg);tri!=NULL;tri=mgfem->get_suivant_triangle(itmg))
243     tri->change_nouveau_numero(0);
244     affiche((char*)" Decoupage des tetra optimises");
245     LISTE_FEM_ELEMENT3::iterator it3;
246     for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
247     {
248     MG_TETRA* tet2=(MG_TETRA*)tet->get_mg_element_maillage();
249     if (tet2->get_origine()==IMPOSE)
250     {
251     tet2->get_triangle1()->change_origine(IMPOSE);
252     tet2->get_triangle2()->change_origine(IMPOSE);
253     tet2->get_triangle3()->change_origine(IMPOSE);
254     tet2->get_triangle4()->change_origine(IMPOSE);
255     tet2->get_noeud1()->change_solution(tet->get_fem_noeud(0)->get_solution());
256     tet2->get_noeud2()->change_solution(tet->get_fem_noeud(1)->get_solution());
257     tet2->get_noeud3()->change_solution(tet->get_fem_noeud(2)->get_solution());
258     tet2->get_noeud4()->change_solution(tet->get_fem_noeud(3)->get_solution());
259     tet2->get_noeud1()->change_nouveau_numero(tet->get_fem_noeud(0)->get_id());
260     tet2->get_noeud2()->change_nouveau_numero(tet->get_fem_noeud(1)->get_id());
261     tet2->get_noeud3()->change_nouveau_numero(tet->get_fem_noeud(2)->get_id());
262     tet2->get_noeud4()->change_nouveau_numero(tet->get_fem_noeud(3)->get_id());
263     if (tet2->get_triangle1()->get_nouveau_numero()==0)
264     tet2->get_triangle1()->change_nouveau_numero(tet2->get_noeud4()->get_id());
265     else
266     tet2->get_triangle1()->change_nouveau_numero(-1);
267     if (tet2->get_triangle2()->get_nouveau_numero()==0)
268     tet2->get_triangle2()->change_nouveau_numero(tet2->get_noeud3()->get_id());
269     else
270     tet2->get_triangle2()->change_nouveau_numero(-1);
271     if (tet2->get_triangle3()->get_nouveau_numero()==0)
272     tet2->get_triangle3()->change_nouveau_numero(tet2->get_noeud1()->get_id());
273     else
274     tet2->get_triangle3()->change_nouveau_numero(-1);
275     if (tet2->get_triangle4()->get_nouveau_numero()==0)
276     tet2->get_triangle4()->change_nouveau_numero(tet2->get_noeud2()->get_id());
277     else
278     tet2->get_triangle4()->change_nouveau_numero(-1);
279     }
280     }
281    
282     TPL_LISTE_ENTITE<MG_TRIANGLE*> tri_impose_interne;
283     for (MG_TRIANGLE* tri=mgfem->get_premier_triangle(itmg);tri!=NULL;tri=mgfem->get_suivant_triangle(itmg))
284     {
285     if (tri->get_origine()==IMPOSE)
286     if (tri->get_lien_topologie()->get_dimension()==3)
287     if (tri->get_nouveau_numero()>0)
288     tri_impose_interne.ajouter(tri);
289     }
290     for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
291     {
292     if (tet->get_mg_element_maillage()->get_origine()!=IMPOSE)
293     {
294     FEM_NOEUD* no1=tet->get_fem_noeud(0);
295     FEM_NOEUD* no2=tet->get_fem_noeud(1);
296     FEM_NOEUD* no3=tet->get_fem_noeud(2);
297     FEM_NOEUD* no4=tet->get_fem_noeud(3);
298     std::vector<MG_NOEUD*> tab;
299     interpole_segment(no1,no2,&tab,limit,mgmai);
300     interpole_segment(no1,no3,&tab,limit,mgmai);
301     interpole_segment(no1,no4,&tab,limit,mgmai);
302     interpole_segment(no2,no3,&tab,limit,mgmai);
303     interpole_segment(no2,no4,&tab,limit,mgmai);
304     interpole_segment(no3,no4,&tab,limit,mgmai);
305     int nb=tab.size();
306     FEM_NOEUD* noext;
307     if (nb>0)
308     {
309     if (no1->get_solution()<limit) noext=no1;
310     if (no2->get_solution()<limit) noext=no2;
311     if (no3->get_solution()<limit) noext=no3;
312     if (no4->get_solution()<limit) noext=no4;
313     }
314     if (nb==3)
315     {
316     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[1],tab[2],mgmai,TRIANGULATION);
317     oriente_tri(tri,noext->get_coord());
318     }
319     if (nb==4)
320     {
321     if (test_du_point_milieu(tab[0],tab[1],tet)==1)
322     {
323     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[1],tab[2],mgmai,TRIANGULATION);
324     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[0],tab[1],tab[3],mgmai,TRIANGULATION);
325     oriente_tri(tri,noext->get_coord());
326     oriente_tri(tri2,noext->get_coord());
327     }
328     else if (test_du_point_milieu(tab[0],tab[2],tet)==1)
329     {
330     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[2],tab[1],mgmai,TRIANGULATION);
331     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[0],tab[2],tab[3],mgmai,TRIANGULATION);
332     oriente_tri(tri,noext->get_coord());
333     oriente_tri(tri2,noext->get_coord());}
334     else if (test_du_point_milieu(tab[0],tab[3],tet)==1)
335     {
336     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[3],tab[1],mgmai,TRIANGULATION);
337     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[0],tab[3],tab[2],mgmai,TRIANGULATION);
338     oriente_tri(tri,noext->get_coord());
339     oriente_tri(tri2,noext->get_coord());}
340     else if (test_du_point_milieu(tab[1],tab[2],tet)==1)
341     {
342     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[1],tab[2],tab[0],mgmai,TRIANGULATION);
343     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[1],tab[2],tab[3],mgmai,TRIANGULATION);
344     oriente_tri(tri,noext->get_coord());
345     oriente_tri(tri2,noext->get_coord());
346     }
347     else if (test_du_point_milieu(tab[1],tab[3],tet)==1)
348     {
349     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[1],tab[3],tab[0],mgmai,TRIANGULATION);
350     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[1],tab[3],tab[2],mgmai,TRIANGULATION);
351     oriente_tri(tri,noext->get_coord());
352     oriente_tri(tri2,noext->get_coord());}
353     else
354     {
355     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[2],tab[3],tab[0],mgmai,TRIANGULATION);
356     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[2],tab[3],tab[1],mgmai,TRIANGULATION);
357     oriente_tri(tri,noext->get_coord());
358     oriente_tri(tri2,noext->get_coord());
359     }
360    
361     }
362     }
363     }
364     affiche((char*)" Decoupage des triangles non design interne au volume");
365     for (int i=0;i<tri_impose_interne.get_nb();i++)
366     {
367     MG_TRIANGLE* tri=tri_impose_interne.get(i);
368     MG_NOEUD* no1=tri->get_noeud1();
369     MG_NOEUD* no2=tri->get_noeud2();
370     MG_NOEUD* no3=tri->get_noeud3();
371     std::vector<MG_NOEUD*> tab;
372     int num1=-1;
373     int num2=-1;
374     int num3=-1;
375     if (no1->get_solution()<limit)
376     {
377     MG_NOEUD* no=get_noeud_peau(mai->get_fem_noeudid(no1->get_nouveau_numero()),mgmai);
378     tab.push_back(no);
379     num1=1;
380     }
381     if (no2->get_solution()<limit)
382     {
383     MG_NOEUD* no=get_noeud_peau(mai->get_fem_noeudid(no2->get_nouveau_numero()),mgmai);
384     tab.push_back(no);
385     num2=1;
386     }
387     if (no3->get_solution()<limit)
388     {
389     MG_NOEUD* no=get_noeud_peau(mai->get_fem_noeudid(no3->get_nouveau_numero()),mgmai);
390     tab.push_back(no);
391     num3=1;
392     }
393     if (num1*num2<0) interpole_segment(mai->get_fem_noeudid(no1->get_nouveau_numero()),mai->get_fem_noeudid(no2->get_nouveau_numero()),&tab,limit,mgmai);
394     if (num1*num3<0) interpole_segment(mai->get_fem_noeudid(no1->get_nouveau_numero()),mai->get_fem_noeudid(no3->get_nouveau_numero()),&tab,limit,mgmai);
395     if (num2*num3<0) interpole_segment(mai->get_fem_noeudid(no2->get_nouveau_numero()),mai->get_fem_noeudid(no3->get_nouveau_numero()),&tab,limit,mgmai);
396     MG_NOEUD* noint=mgfem->get_mg_noeudid(tri->get_nouveau_numero());
397     int nb=tab.size();
398     if (nb==3)
399     {
400     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[1],tab[2],mgmai,IMPOSE);
401     oriente_tri(tri,noint->get_coord());
402     tri->inverse_sens();
403     }
404     if (nb==4)
405     {
406     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[1],tab[2],mgmai,IMPOSE);
407     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[1],tab[2],tab[3],mgmai,IMPOSE);
408     oriente_tri(tri,noint->get_coord());
409     tri->inverse_sens();
410     oriente_tri(tri2,noint->get_coord());
411     tri2->inverse_sens();
412     }
413     }
414 francois 345 affiche((char*)" Decoupage des triangles non design sur la frontiere du volume");
415 francois 343 LISTE_FEM_ELEMENT2::iterator it2;
416     for (FEM_ELEMENT2* tri=mai->get_premier_element2(it2);tri!=NULL;tri=mai->get_suivant_element2(it2))
417     {
418     std::vector<MG_NOEUD*> tab;
419     FEM_NOEUD *no1=tri->get_fem_noeud(0);
420     FEM_NOEUD *no2=tri->get_fem_noeud(1);
421     FEM_NOEUD *no3=tri->get_fem_noeud(2);
422     int ori=((MG_FACE*)tri->get_lien_topologie())->get_mg_coface(0)->get_orientation();
423     OT_VECTEUR_3D n1n2(no1->get_coord(),no2->get_coord());
424     OT_VECTEUR_3D n1n3(no1->get_coord(),no3->get_coord());
425     OT_VECTEUR_3D nor=n1n3&n1n2;
426     double xyzext[3];
427     xyzext[0]=no1->get_x()+nor.get_x()*ori;
428     xyzext[1]=no1->get_y()+nor.get_y()*ori;
429     xyzext[2]=no1->get_z()+nor.get_z()*ori;
430     if (tri->get_mg_element_maillage()->get_origine()==IMPOSE)
431     {
432     MG_NOEUD *mgno1=get_noeud_peau(no1,mgmai);
433     MG_NOEUD *mgno2=get_noeud_peau(no2,mgmai);
434     MG_NOEUD *mgno3=get_noeud_peau(no3,mgmai);
435     int num1=interpole_segment(no1,no2,&tab,limit,mgmai,0);
436     int num2=interpole_segment(no1,no3,&tab,limit,mgmai,0);
437     int num3=interpole_segment(no2,no3,&tab,limit,mgmai,0);
438     int nb=tab.size();
439     if (nb==0)
440     {
441     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,mgno1,mgno2,mgno3,mgmai,IMPOSE);
442     oriente_tri(tri,xyzext);
443     }
444     if (nb==1)
445     {
446     if (num1==1)
447     {
448     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,mgno1,mgno3,tab[0],mgmai,IMPOSE);
449     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,mgno2,mgno3,tab[0],mgmai,IMPOSE);
450     oriente_tri(tri,xyzext);
451     oriente_tri(tri2,xyzext);
452     }
453     if (num2==1)
454     {
455     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,mgno1,mgno2,tab[0],mgmai,IMPOSE);
456     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,mgno2,mgno3,tab[0],mgmai,IMPOSE);
457     oriente_tri(tri,xyzext);
458     oriente_tri(tri2,xyzext);
459     }
460     if (num3==1)
461     {
462     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,mgno1,mgno2,tab[0],mgmai,IMPOSE);
463     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,mgno1,mgno3,tab[0],mgmai,IMPOSE);
464     oriente_tri(tri,xyzext);
465     oriente_tri(tri2,xyzext);
466     }
467     }
468     if (nb==2)
469     {
470     if (num1==0)
471     {
472     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,mgno1,mgno2,tab[0],mgmai,IMPOSE);
473     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,mgno2,tab[0],tab[1],mgmai,IMPOSE);
474     MG_TRIANGLE_PEAU* tri3=insere_triangle(NULL,mgno3,tab[0],tab[1],mgmai,IMPOSE);
475     oriente_tri(tri,xyzext);
476     oriente_tri(tri2,xyzext);
477     oriente_tri(tri3,xyzext);
478     }
479     if (num2==0)
480     {
481     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,mgno1,mgno3,tab[0],mgmai,IMPOSE);
482     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,mgno3,tab[0],tab[1],mgmai,IMPOSE);
483     MG_TRIANGLE_PEAU* tri3=insere_triangle(NULL,mgno2,tab[0],tab[1],mgmai,IMPOSE);
484     oriente_tri(tri,xyzext);
485     oriente_tri(tri2,xyzext);
486     oriente_tri(tri3,xyzext);
487     }
488     if (num2==0)
489     {
490     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,mgno2,mgno3,tab[0],mgmai,IMPOSE);
491     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,mgno3,tab[0],tab[1],mgmai,IMPOSE);
492     MG_TRIANGLE_PEAU* tri3=insere_triangle(NULL,mgno1,tab[0],tab[1],mgmai,IMPOSE);
493     oriente_tri(tri,xyzext);
494     oriente_tri(tri2,xyzext);
495     oriente_tri(tri3,xyzext);
496     }
497     }
498     }
499     else if ((no1->get_solution()>=limit) && (no2->get_solution()>=limit) && (no3->get_solution()>=limit))
500     {
501     MG_NOEUD *mgno1=get_noeud_peau(no1,mgmai);
502     MG_NOEUD *mgno2=get_noeud_peau(no2,mgmai);
503     MG_NOEUD *mgno3=get_noeud_peau(no3,mgmai);
504     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,mgno1,mgno2,mgno3,mgmai,TRIANGULATION);
505     oriente_tri(tri,xyzext);
506     }
507     else if (!((no1->get_solution()<limit) && (no2->get_solution()<limit) && (no3->get_solution()<limit)))
508     {
509     std::vector<MG_NOEUD*> tab;
510     int num1=-1;
511     int num2=-1;
512     int num3=-1;
513     if (no1->get_solution()>=limit)
514     {
515     MG_NOEUD* no=get_noeud_peau(no1,mgmai);
516     tab.push_back(no);
517     num1=1;
518     }
519     if (no2->get_solution()>=limit)
520     {
521     MG_NOEUD* no=get_noeud_peau(no2,mgmai);
522     tab.push_back(no);
523     num2=1;
524     }
525     if (no3->get_solution()>=limit)
526     {
527     MG_NOEUD* no=get_noeud_peau(no3,mgmai);
528     tab.push_back(no);
529     num3=1;
530     }
531     if (num1*num2<0) interpole_segment(no1,no2,&tab,limit,mgmai);
532     if (num1*num3<0) interpole_segment(no1,no3,&tab,limit,mgmai);
533     if (num2*num3<0) interpole_segment(no2,no3,&tab,limit,mgmai);
534     int nb=tab.size();
535     if (nb==3)
536     {
537     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[1],tab[2],mgmai,TRIANGULATION);
538     oriente_tri(tri,xyzext);
539     }
540     if (nb==4)
541     {
542     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[1],tab[2],mgmai,TRIANGULATION);
543     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[1],tab[2],tab[3],mgmai,TRIANGULATION);
544     oriente_tri(tri,xyzext);
545     oriente_tri(tri2,xyzext);
546     }
547     }
548     }
549     LISTE_MG_TRIANGLE::iterator it_tri;
550     for (MG_TRIANGLE* mgtri=mgmai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mgmai->get_suivant_triangle(it_tri))
551     {
552     MG_TRIANGLE_PEAU* tripeau=(MG_TRIANGLE_PEAU*)mgtri;
553     MG_TRIANGLE_PEAU* voisin1=recherche_voisin(tripeau->get_noeud1(),tripeau->get_noeud2(),tripeau);
554     MG_TRIANGLE_PEAU* voisin2=recherche_voisin(tripeau->get_noeud2(),tripeau->get_noeud3(),tripeau);
555     MG_TRIANGLE_PEAU* voisin3=recherche_voisin(tripeau->get_noeud3(),tripeau->get_noeud1(),tripeau);
556     tripeau->change_voisin1(voisin1);
557     tripeau->change_voisin2(voisin2);
558     tripeau->change_voisin3(voisin3);
559     tripeau->change_nouveau_numero(0);
560     }
561     int fin;
562     do
563     {
564     fin=1;
565     for (MG_TRIANGLE* mgtri=mgmai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mgmai->get_suivant_triangle(it_tri))
566     {
567     MG_TRIANGLE_PEAU *tripeau=(MG_TRIANGLE_PEAU*)mgtri;
568     if (tripeau->get_nouveau_numero()==0)
569     {
570     fin=0;
571     std::vector<MG_TRIANGLE_PEAU*> *peau=new std::vector<MG_TRIANGLE_PEAU*>;
572     lst_peau.push_back(peau);
573     tripeau->change_nouveau_numero(1);
574     peau->push_back(tripeau);
575     determine_peau(peau);
576     }
577     }
578     }
579     while (fin==0);
580 francois 457 return mgmai;
581     }
582     void MG_LISSAGE::suppression_peaux_isoles(MG_MAILLAGE* mgmai)
583     {
584     affiche((char*)"Suppression des peaux isolées");
585 francois 343
586     char message[500];
587     for (int cas=0;cas<2;cas++)
588     {
589 francois 457 if (cas==0) affiche((char*)" Analyse initiale des peaux");
590     if (cas==1) affiche((char*)" Analyse finale des peaux");
591     int nbisole=0;
592 francois 343 for (int i=0;i<lst_peau.size();i++)
593     {
594     int isole=1;
595     for (int j=0;j<lst_peau[i]->size();j++)
596     if ((*lst_peau[i])[j]->get_origine()==IMPOSE) {isole=0;break;}
597     if (isole==1)
598     {
599 francois 457 nbisole++;
600 francois 343 for (int j=0;j<lst_peau[i]->size();j++)
601     {
602     mgmai->supprimer_mg_triangleid((*lst_peau[i])[j]->get_id());
603     }
604     lst_peau[i]->clear();
605     }
606     }
607 francois 457 sprintf(message," %d peaux, %d non isoles, %d isoles",(int)lst_peau.size(),(int)lst_peau.size()-nbisole,nbisole);
608     affiche(message);
609 francois 343 for (std::vector<std::vector<MG_TRIANGLE_PEAU*> *>::iterator j=lst_peau.end()-1;j>lst_peau.begin();j--)
610     if ((*j)->size()==0)
611     lst_peau.erase(j);
612     }
613 francois 457 LISTE_MG_TRIANGLE::iterator itmg;
614 francois 343 for (MG_TRIANGLE* tri=mgmai->get_premier_triangle(itmg);tri!=NULL;tri=mgmai->get_suivant_triangle(itmg))
615     {
616     if (tri->get_origine()==IMPOSE)
617     {
618     tri->get_noeud1()->change_origine(IMPOSE);
619     tri->get_noeud2()->change_origine(IMPOSE);
620     tri->get_noeud3()->change_origine(IMPOSE);
621     }
622     }
623 cuillier 438
624 francois 343 }
625    
626     void MG_LISSAGE::oriente_tri(MG_TRIANGLE_PEAU* tri,double *xyz)
627     {
628     MG_NOEUD* no1=tri->get_noeud1();
629     MG_NOEUD* no2=tri->get_noeud2();
630     MG_NOEUD* no3=tri->get_noeud3();
631     OT_VECTEUR_3D n1n2(no1->get_coord(),no2->get_coord());
632     OT_VECTEUR_3D n1n3(no1->get_coord(),no3->get_coord());
633     OT_VECTEUR_3D normal=n1n2&n1n3;
634     OT_VECTEUR_3D dir(no1->get_coord(),xyz);
635     dir.norme();
636     normal.norme();
637     double ps=normal*dir;
638     if (ps<0) tri->inverse_sens();
639     }
640    
641    
642    
643     MG_NOEUD* MG_LISSAGE::get_noeud_peau(FEM_NOEUD* no,MG_MAILLAGE* mai)
644     {
645     static std::map<std::string,MG_NOEUD*> correspond;
646     unsigned long id=no->get_id();
647     char message[255];
648     sprintf(message,"_%lu_",id);
649     std::string key=message;
650     MG_NOEUD* mgno=correspond[key];
651     if (mgno==NULL)
652     {
653     mgno=new MG_NOEUD(NULL,no->get_x(),no->get_y(),no->get_z(),TRIANGULATION);
654     mai->ajouter_mg_noeud(mgno);
655     correspond[key]=mgno;
656     }
657     return mgno;
658     }
659    
660    
661     int MG_LISSAGE::test_du_point_milieu(MG_NOEUD* no1,MG_NOEUD* no2,FEM_ELEMENT3* tet)
662     {
663     double *xyz1=no1->get_coord();
664     double *xyz2=no2->get_coord();
665     double xyz[3];
666     xyz[0]=0.5*(xyz1[0]+xyz2[0]);
667     xyz[1]=0.5*(xyz1[1]+xyz2[1]);
668     xyz[2]=0.5*(xyz1[2]+xyz2[2]);
669 francois 412 FEM_MAILLAGE_OUTILS ot;
670 francois 388 if (((ot.estdansletetra(tet,xyz[0],xyz[1],xyz[2])>>1)&1)==1) return 1;
671     //if (((MG_MAILLAGE_OUTILS::estdansletetra(tet,xyz[0],xyz[1],xyz[2])>>1)&1)==1) return 1;
672 francois 343 return 0;
673     }
674    
675    
676    
677     int MG_LISSAGE::interpole_segment(FEM_NOEUD* no1,FEM_NOEUD* no2,std::vector<MG_NOEUD*> *tab,double limit,MG_MAILLAGE* mai,int creation)
678     {
679     static std::map<std::string,MG_NOEUD*> correspond;
680     unsigned long id1=no1->get_id();
681     unsigned long id2=no2->get_id();
682     char message[255];
683     sprintf(message,"_%lu_%lu_",id1,id2);
684     std::string key1=message;
685     sprintf(message,"_%lu_%lu_",id2,id1);
686     std::string key2=message;
687     double sol1=no1->get_solution();
688     double sol2=no2->get_solution();
689 francois 345 if (fabs(sol1-limit)<0.0000001)
690     {
691     sprintf(message,"%lu",id1);
692     std::string key=message;
693     MG_NOEUD* no=correspond[key];
694     if (no==NULL)
695     {
696     double *xyz1=no1->get_coord();
697     no=new MG_NOEUD(NULL,xyz1[0],xyz1[1],xyz1[2],TRIANGULATION);
698     mai->ajouter_mg_noeud(no);
699     correspond[key]=no;
700     }
701     int present=0;
702     for (int i=0;i<tab->size();i++)
703     if ((*tab)[i]==no) present=1;
704     if (present==0)
705     {
706     tab->push_back(no);
707     return 1;
708     }
709     else return 0;
710     }
711     if (fabs(sol2-limit)<0.0000001)
712     {
713     sprintf(message,"%lu",id2);
714     std::string key=message;
715     MG_NOEUD* no=correspond[key];
716     if (no==NULL)
717     {
718     double *xyz2=no2->get_coord();
719     no=new MG_NOEUD(NULL,xyz2[0],xyz2[1],xyz2[2],TRIANGULATION);
720     mai->ajouter_mg_noeud(no);
721     correspond[key]=no;
722     }
723     int present=0;
724     for (int i=0;i<tab->size();i++)
725     if ((*tab)[i]==no) present=1;
726     if (present==0)
727     {
728     tab->push_back(no);
729     return 1;
730     }
731     else return 0;
732    
733     }
734    
735     if (((sol1<limit) && (sol2>limit))||((sol1>limit) && (sol2<limit)))
736 francois 343 {
737     MG_NOEUD* no=correspond[key1];
738     if (no==NULL) no=correspond[key2];
739     if ((no==NULL) && (creation==1))
740     {
741     double t=(limit-sol1)/(sol2-sol1);
742     double xyz[3];
743     double *xyz1=no1->get_coord();
744     double *xyz2=no2->get_coord();
745     xyz[0]=xyz1[0]+t*(xyz2[0]-xyz1[0]);
746     xyz[1]=xyz1[1]+t*(xyz2[1]-xyz1[1]);
747     xyz[2]=xyz1[2]+t*(xyz2[2]-xyz1[2]);
748     no=new MG_NOEUD(NULL,xyz[0],xyz[1],xyz[2],TRIANGULATION);
749     mai->ajouter_mg_noeud(no);
750     correspond[key1]=no;
751     }
752     if (no!=NULL)
753     {
754     tab->push_back(no);
755     return 1;
756     }
757     }
758     return 0;
759     }
760    
761 francois 457 void MG_LISSAGE::lire_params(char *fichier)
762     {
763     params.vide();
764     params.lire(fichier);
765     }
766 francois 343
767 francois 457 void MG_LISSAGE::lisse(FEM_SOLUTION* sol, MG_GESTIONNAIRE& gest2,char *nomparam)
768 picher 233 {
769 francois 457 if (nomparam!=NULL) lire_params(nomparam);
770 francois 272 affiche((char*)"Extraction du maillage de surface");
771 francois 457 FEM_MAILLAGE* fem=sol->get_maillage();
772     int decoup=params.get_valeur("decoupage");
773     MG_MAILLAGE* mgmai;
774     if (decoup==0)
775     {
776     affiche((char*)" Découpage par mailles 3D entieres");
777     adapte_seuil(fem,sol);
778     mgmai=extract_skin_maille_entiere(fem,gest2);
779     }
780     if (decoup==1)
781     {
782     affiche((char*)" Découpage par isodensité");
783     double seuil=params.get_valeur("seuil");
784     std::string nomfich=params.get_nom("nomdensextra");
785     mgmai=extract_skin_par_decoupage(sol,seuil,gest2,nomfich.c_str());
786     }
787     suppression_peaux_isoles(mgmai);
788     affiche((char*)"Procedure de lissage");
789     int chen2005_debut=params.get_valeur("chen2005_debut");
790     int chen2008_debut=params.get_valeur("chen2008_debut");
791     int jiao2012_debut=params.get_valeur("jiao2012_debut");
792     int chen2005fait=0;
793     if (chen2005_debut==0) chen2005fait=1;
794     int chen2008fait=0;
795     if (chen2008_debut==0) chen2008fait=1;
796     int jiao2012fait=0;
797     if (jiao2012_debut==0) jiao2012fait=1;
798     int fin=0;
799     int iteration=1;
800     do
801     {
802     if ((iteration>=chen2005_debut) && (chen2005fait==0))
803     {
804     char message[300];
805     sprintf(message, " Iteration %d : lissage de Chen 2005",iteration);
806     affiche(message);
807     double epsilon=params.get_valeur("chen2005_epsilon");
808     double sigma=params.get_valeur("chen2005_sigma");
809     int iter_max=params.get_valeur("chen2005_itermax");
810     int nbiter=lissage_chen2005(mgmai,gest2,epsilon,sigma,iter_max);
811     if (nbiter==iter_max) sprintf(message," Arrêt de la procédure de lissage de Chen2005 après %d itérations",nbiter);
812     else sprintf(message," Convergence de la procédure de lissage de Chen2005 après %d itérations",nbiter);
813     affiche(message);
814     chen2005fait=1;
815     iteration=iteration+nbiter;
816     }
817     else if ((iteration>=chen2008_debut) && (chen2008fait==0))
818     {
819     char message[300];
820     sprintf(message, " Iteration %d : lissage de Chen 2008",iteration);
821     affiche(message);
822     double epsilon=params.get_valeur("chen2008_epsilon");
823     double sigma=params.get_valeur("chen2008_sigma");
824     double gamma=params.get_valeur("chen2008_gamma");
825     int iter_max=params.get_valeur("chen2008_itermax");
826     int nbiter=lissage_chen2008(mgmai,gest2,sigma,gamma,epsilon,iter_max);
827     if (nbiter==iter_max) sprintf(message," Arrêt de la procédure de lissage de Chen2008 après %d itérations",nbiter);
828     else sprintf(message," Convergence de la procédure de lissage de Chen2008 après %d itérations",nbiter);
829     affiche(message);
830     chen2008fait=1;
831     iteration=iteration+nbiter;
832     }
833     else if ((iteration>=jiao2012_debut) && (jiao2012fait==0))
834     {
835     char message[300];
836     sprintf(message, " Iteration %d : lissage de Jiao 2012",iteration);
837     affiche(message);
838     int iter_max=params.get_valeur("jiao2012_itermax");
839     int nbiter=lissage_jiao2012(mgmai,gest2,iter_max);
840     if (nbiter==iter_max) sprintf(message," Arrêt de la procédure de lissage de Jiao2012 après %d itérations",nbiter);
841     else sprintf(message," Convergence de la procédure de lissage de Jiao2012 après %d itérations",nbiter);
842     affiche(message);
843     jiao2012fait=1;
844     iteration=iteration+nbiter;
845     }
846     else iteration++;
847     if (chen2005fait==1)
848     if (chen2008fait==1)
849     if (jiao2012fait==1)
850     fin=1;
851     }
852     while (fin==0);
853     }
854    
855    
856     MG_MAILLAGE* MG_LISSAGE::extract_skin_maille_entiere(FEM_MAILLAGE* mai,MG_GESTIONNAIRE &gest2)
857     {
858     //Menu de la methode de lissage
859     affiche((char*)"Extraction du maillage de surface");
860 picher 231 int coderesu = 0;
861 picher 233 int mai2_id;
862 francois 457 int imp=params.get_valeur("consimpose");
863     int opti=params.get_valeur("consoptimise");
864     int m_auto=params.get_valeur("consnooptimise");
865 francois 232 if (opti==1) conserve(OPTIMISE);
866     if (imp==1) conserve(IMPOSE);
867 picher 233 if (m_auto==1) conserve(MAILLEUR_AUTO);
868 francois 457 copieorigine(mai);
869 picher 233 MG_MAILLAGE* mg_mai = (MG_MAILLAGE*)mai->get_mg_maillage();
870 francois 457 //gain_poids(mg_mai,gest2);
871     int reactiv=params.get_valeur("reactivation");
872 picher 233 if (reactiv == 1)
873     {
874     reactivation(mg_mai,gest2);
875     }
876 francois 457 affiche((char*)" Analyse des cas non manifold");
877 picher 231 do
878     {
879 francois 457 int nbmaniare,nbmanino,nbpeau;
880     coderesu = extract_skin(mg_mai,gest2,nbpeau,nbmaniare,nbmanino,&mai2_id);
881     char message[500];
882     sprintf(message," %d peaux. %d manifold par arete. %d manifold par noeud",nbpeau,nbmaniare,nbmanino);
883     affiche(message);
884     // gain_poids(mg_mai,gest2);
885 picher 231 }
886     while (coderesu == 0);
887 picher 233
888     MG_MAILLAGE* mg_mai2=gest2.get_mg_maillageid(mai2_id);
889 francois 457
890     return mg_mai2;
891    
892 picher 231 }
893    
894 francois 457 void MG_LISSAGE::adapte_seuil(class FEM_MAILLAGE* fem,FEM_SOLUTION* solution)
895 picher 230 {
896 francois 457 double seuil=params.get_valeur("seuil");
897     solution->active_solution(0);
898     LISTE_FEM_ELEMENT3::iterator it;
899     for (FEM_ELEMENT3 *tet=fem->get_premier_element3(it);tet!=NULL;tet=fem->get_suivant_element3(it))
900 picher 230 {
901 francois 457 if (((MG_TETRA*)tet->get_mg_element_maillage())->get_origine()!=IMPOSE)
902     if (tet->get_solution()>seuil)
903     ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(OPTIMISE);
904     else
905     ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(MAILLEUR_AUTO);
906    
907 picher 230 }
908     }
909 francois 224
910 francois 457
911    
912    
913     //---------------------------- Lissage Chen 2005 Gilles Philippe -------------------------------------------------------------------------
914     int MG_LISSAGE::lissage_chen2005(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double epsilon, double sigma, int iter_max)
915 picher 230 {
916 francois 457 //Methode de lissage inspiree de Chen(2005) par Gilles Philipe Picher Martel
917 picher 230 double un_sur_pi = 1./M_PI;
918     int compteur = 0;
919     int fin = 0;
920 francois 224
921 francois 457
922 picher 230 do
923     {
924     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
925 picher 233 TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales2;
926 picher 230 TPL_LISTE_ENTITE<double> liste_wij;
927     LISTE_MG_TRIANGLE::iterator it_tri;
928     int k = 0; //pour identifier les triangles pour liste_normales et liste_wij
929     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
930     {
931     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
932     OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
933 picher 233 normal_f_i.norme();
934     liste_normales2.ajouter(normal_f_i);
935 picher 230 //Remplissage de la liste des voisins du triangle i
936     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
937     MG_NOEUD* noeud1 = mgtri_i->get_noeud1();
938     double nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
939     for (int j = 0;j<nb_voisins1;j++)
940     {
941     MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
942     liste_voisins.ajouter(mgtri_1);
943     }
944     MG_NOEUD* noeud2 = mgtri_i->get_noeud2();
945     double nb_voisins2 = noeud2->get_lien_triangle()->get_nb();
946     for (int j = 0;j<nb_voisins2;j++)
947     {
948     MG_TRIANGLE_PEAU* mgtri_2 = (MG_TRIANGLE_PEAU*) noeud2->get_lien_triangle()->get(j);
949     liste_voisins.ajouter(mgtri_2);
950     }
951     MG_NOEUD* noeud3 = mgtri_i->get_noeud3();
952     double nb_voisins3 = noeud3->get_lien_triangle()->get_nb();
953     for (int j = 0;j<nb_voisins3;j++)
954     {
955     MG_TRIANGLE_PEAU* mgtri_3 = (MG_TRIANGLE_PEAU*) noeud3->get_lien_triangle()->get(j);
956     liste_voisins.ajouter(mgtri_3);
957     }
958     liste_voisins.supprimer(mgtri_i);
959     int nb_voisins = liste_voisins.get_nb();
960     double w_ij = 1./nb_voisins;
961     double phi_i_min = 10.;
962     double s_i = 0.0;
963     double phi_im = 0.0;
964 picher 231 double *phi_ij = new double[nb_voisins];
965 picher 230 OT_VECTEUR_3D normal_f_i_mean(0.,0.,0.);
966     OT_VECTEUR_3D eta_i(0.,0.,0.);
967     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
968     int j = 0;
969     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
970     {
971     OT_VECTEUR_3D normal_f_j = mgtri_j->calcul_normal();
972     //1-Calculer la normale moyenne pour chaque triangle
973     normal_f_i_mean = normal_f_i_mean + w_ij*normal_f_j;
974 francois 457 //2.1-On calcul l'angle entre normal_f_i et normal_f_j pour j allant de 1 a Nb_voisins
975 picher 230 double prod_scalaire = normal_f_i*normal_f_j;
976 picher 233 if (prod_scalaire > 1.)
977     {
978     prod_scalaire = 1.;
979     }
980     if (prod_scalaire < -1.)
981     {
982     prod_scalaire = -1.;
983     }
984 picher 230 phi_ij[j] = acos(prod_scalaire)*un_sur_pi;
985     //2.2-On trouve le plus petit des angles et la normale heta_i correspondante
986     if (phi_ij[j] < phi_i_min)
987     {
988     phi_i_min = phi_ij[j];
989     eta_i = normal_f_j;
990     }
991     //3.1-On calcul l'angle moyen phi_im
992     phi_im = phi_im + w_ij*phi_ij[j];
993     j++;
994     }
995     normal_f_i_mean.norme();
996     j = 0;
997     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
998     {
999     //3.2-Calcul de s_i selon la variance
1000     s_i = s_i + w_ij*pow((phi_ij[j] - phi_im),2);
1001     j++;
1002     }
1003 picher 231 delete[] phi_ij;
1004 picher 230 //4-On calcule une nouvelle normale pour chaque triangle
1005     OT_VECTEUR_3D normal_f_i_new = ponderation_gaussian(s_i,sigma)*normal_f_i_mean + (1. - ponderation_gaussian(s_i,sigma))*eta_i;
1006     normal_f_i_new.norme();
1007     liste_normales.ajouter(normal_f_i_new);
1008     liste_wij.ajouter(w_ij);
1009     mgtri->change_nouveau_numero(k);
1010     k++;
1011     }
1012    
1013     LISTE_MG_NOEUD::iterator it_no;
1014     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1015     {
1016     int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
1017     double w_ij_prime = 0.0;
1018     OT_VECTEUR_3D v_temp(0.,0.,0.);
1019     OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
1020     for(int j=0;j<nb_voisins_j;j++)
1021     {
1022     MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
1023 francois 457 //On calcule le centroide cj du triangle mgtri_j
1024 picher 230 MG_NOEUD* n1 = mgtri_j->get_noeud1();
1025     MG_NOEUD* n2 = mgtri_j->get_noeud2();
1026     MG_NOEUD* n3 = mgtri_j->get_noeud3();
1027     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
1028     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
1029     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
1030     //On forme le vecteur vi_cj
1031     OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
1032     OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_j->get_nouveau_numero());
1033 francois 457 // w_ij_prime correspond a la somme des aires des triangles voisins du noeuds
1034 picher 233 OT_VECTEUR_3D AB(n2->get_x() - n1->get_x(),n2->get_y() - n1->get_y(),n2->get_z() - n1->get_z());
1035     OT_VECTEUR_3D AC(n3->get_x() - n1->get_x(),n3->get_y() - n1->get_y(),n3->get_z() - n1->get_z());
1036     OT_VECTEUR_3D prodvect = AB&AC;
1037     double w_ij = 0.5*prodvect.get_longueur();
1038 picher 230 w_ij_prime = w_ij_prime + w_ij;
1039     v_temp = v_temp + w_ij*(vi_cj*normal_f_i_new)*normal_f_i_new;
1040     }
1041 francois 457 //5-On met a jour la position des noeuds
1042 picher 230 v_i = v_i + v_temp/w_ij_prime;
1043 picher 233 int origine = noeud_i->get_origine();
1044     if (origine != IMPOSE)
1045     {
1046     noeud_i->change_x(v_i.get_x());
1047     noeud_i->change_y(v_i.get_y());
1048     noeud_i->change_z(v_i.get_z());
1049     }
1050 picher 230 }
1051 francois 457 //Critere d'arret de l'algorithme
1052 picher 230 int l=0;
1053     int nb_tri = mg_mai->get_nb_mg_triangle();
1054     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
1055     {
1056     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
1057 picher 233 OT_VECTEUR_3D normal_f_i = liste_normales2.get(mgtri_i->get_nouveau_numero());
1058 picher 230 OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_i->get_nouveau_numero());
1059 picher 248 double critere = 1. - normal_f_i*normal_f_i_new;
1060 picher 230 if (critere <= epsilon) l++;
1061     }
1062 picher 248 double tolerance = 0.01*nb_tri;
1063     if (nb_tri - l <= 0) fin = 1;
1064 picher 230 compteur++;
1065     }
1066     while ((fin == 0) && (compteur < iter_max));
1067    
1068 francois 457 return compteur;
1069     }
1070 picher 230
1071 francois 457 //---------------------------- Lissage Jiao 2012 -------------------------------------------------------------------------
1072     int MG_LISSAGE::lissage_jiao2012(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, int iter_max)
1073 picher 248 {
1074 francois 457 // Lissage global avec methode de lissage inspiree de JIAO IMR2012 par Jean Christophe sept 2013
1075 picher 248 double un_sur_pi = 1./M_PI;
1076     int compteur = 0;
1077    
1078 cuillier 438 do
1079     {
1080     vector<double> nouv_position_x;
1081     vector<double> nouv_position_y;
1082     vector<double> nouv_position_z;
1083     LISTE_MG_NOEUD::iterator it_no;
1084     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1085     {
1086     int origine = noeud_i->get_origine();
1087     if (origine != IMPOSE)
1088     {
1089     int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
1090     double wij = 0.0;
1091     double wij_sommej = 0.0;
1092     OT_VECTEUR_3D v_temp(0.,0.,0.);
1093     OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
1094     for(int j=0;j<nb_voisins_j;j++)
1095     {
1096     MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
1097 francois 457 //On calcule le centroide cj du triangle mgtri_j
1098 cuillier 438 MG_NOEUD* n1 = mgtri_j->get_noeud1();
1099     MG_NOEUD* n2 = mgtri_j->get_noeud2();
1100     MG_NOEUD* n3 = mgtri_j->get_noeud3();
1101     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
1102     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
1103     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
1104     //On forme le vecteur vi_cj et le vecteur cj
1105     OT_VECTEUR_3D cj(cj_x ,cj_y,cj_z);
1106     OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
1107     wij=vi_cj.get_longueur();
1108     wij_sommej=wij_sommej+wij;
1109     v_temp = v_temp + (wij*cj);
1110     }
1111     //5-On calcule la nouvelle position des noeuds et on affecte au tableau temporaire
1112     v_i =v_temp/wij_sommej;
1113     nouv_position_x.push_back(v_i.get_x());
1114     nouv_position_y.push_back(v_i.get_y());
1115     nouv_position_z.push_back(v_i.get_z());
1116     }
1117     }
1118     //On actualise la position des noeuds
1119     int ind_noeud=0;
1120     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1121     {
1122     int origine = noeud_i->get_origine();
1123     if (origine != IMPOSE)
1124     {
1125     noeud_i->change_x(nouv_position_x[ind_noeud]);
1126     noeud_i->change_y(nouv_position_y[ind_noeud]);
1127     noeud_i->change_z(nouv_position_z[ind_noeud]);
1128     ind_noeud++;
1129     }
1130     }
1131     //Critere d'arret de l'algorithme
1132     compteur++;
1133     }
1134     while (compteur < iter_max);
1135    
1136 francois 457 return compteur;
1137 cuillier 438 }
1138    
1139 francois 457 //---------------------------- Lissage Chen 2008 -------------------------------------------------------------------------
1140     int MG_LISSAGE::lissage_chen2008(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double sigma, double gamma, double epsilon, int iter_max)
1141 cuillier 438 {
1142 francois 457 //Lissage global inspire de Chen(2008) corrige par Jean Christophe sept 2013
1143 cuillier 438 double un_sur_pi = 1./M_PI;
1144     int compteur = 0;
1145     int fin = 0;
1146 francois 457 compteur = 0;
1147 cuillier 438
1148     do
1149     {
1150     vector<double> nouv_position_x;
1151     vector<double> nouv_position_y;
1152     vector<double> nouv_position_z;
1153     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
1154     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales2;
1155     TPL_LISTE_ENTITE<double> liste_wij;
1156     LISTE_MG_TRIANGLE::iterator it_tri;
1157     int k = 0; //pour identifier les triangles pour liste_normales et liste_wij
1158     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
1159     {
1160     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
1161     OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
1162     normal_f_i.norme();
1163     liste_normales2.ajouter(normal_f_i);
1164     //Remplissage de la liste des voisins du triangle i
1165     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
1166     MG_NOEUD* noeud1 = mgtri_i->get_noeud1();
1167     double nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
1168     for (int j = 0;j<nb_voisins1;j++)
1169     {
1170     MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
1171     liste_voisins.ajouter(mgtri_1);
1172     }
1173     MG_NOEUD* noeud2 = mgtri_i->get_noeud2();
1174     double nb_voisins2 = noeud2->get_lien_triangle()->get_nb();
1175     for (int j = 0;j<nb_voisins2;j++)
1176     {
1177     MG_TRIANGLE_PEAU* mgtri_2 = (MG_TRIANGLE_PEAU*) noeud2->get_lien_triangle()->get(j);
1178     liste_voisins.ajouter(mgtri_2);
1179     }
1180     MG_NOEUD* noeud3 = mgtri_i->get_noeud3();
1181     double nb_voisins3 = noeud3->get_lien_triangle()->get_nb();
1182     for (int j = 0;j<nb_voisins3;j++)
1183     {
1184     MG_TRIANGLE_PEAU* mgtri_3 = (MG_TRIANGLE_PEAU*) noeud3->get_lien_triangle()->get(j);
1185     liste_voisins.ajouter(mgtri_3);
1186     }
1187     liste_voisins.supprimer(mgtri_i);
1188     int nb_voisins = liste_voisins.get_nb();
1189     double w_ij = 1./nb_voisins;
1190     double phi_i_min = 10.;
1191     double s_i = 0.0;
1192     double phi_im = 0.0;
1193     double *phi_ij = new double[nb_voisins];
1194     OT_VECTEUR_3D normal_f_i_mean(0.,0.,0.);
1195     normal_f_i_mean = normal_f_i;
1196     OT_VECTEUR_3D eta_i(0.,0.,0.);
1197     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
1198     int j = 0;
1199     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
1200     {
1201     OT_VECTEUR_3D normal_f_j = mgtri_j->calcul_normal();
1202     //1-Calculer la normale moyenne pour chaque triangle
1203     normal_f_i_mean = normal_f_i_mean + normal_f_j;
1204 francois 457 //2.1-On calcule l'angle entre normal_f_i et normal_f_j pour j allant de 1 a Nb_voisins
1205 cuillier 438 double prod_scalaire = normal_f_i*normal_f_j;
1206     if (prod_scalaire > 1.)
1207     {
1208     prod_scalaire = 1.;
1209     }
1210     if (prod_scalaire < -1.)
1211     {
1212     prod_scalaire = -1.;
1213     }
1214     phi_ij[j] = acos(prod_scalaire)*un_sur_pi;
1215     //2.2-On trouve le plus petit des angles et la normale heta_i correspondante
1216     if (phi_ij[j] < phi_i_min)
1217     {
1218     phi_i_min = phi_ij[j];
1219     eta_i = normal_f_j;
1220     }
1221 francois 457 //3.1-On calcule l'angle moyen phi_im
1222 cuillier 438 phi_im = phi_im + w_ij*phi_ij[j];
1223     j++;
1224     }
1225     normal_f_i_mean.norme();
1226     j = 0;
1227     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
1228     {
1229     //3.2-Calcul de s_i selon la variance
1230     s_i = s_i + w_ij*pow((phi_ij[j] - phi_im),2);
1231     j++;
1232     }
1233     delete[] phi_ij;
1234     //4-On calcule une nouvelle normale pour chaque triangle
1235     OT_VECTEUR_3D normal_f_i_new = ponderation_gaussian(s_i,sigma)*normal_f_i_mean + (1. - ponderation_gaussian(s_i,sigma))*eta_i;
1236     normal_f_i_new.norme();
1237     liste_normales.ajouter(normal_f_i_new);
1238     liste_wij.ajouter(w_ij);
1239     mgtri->change_nouveau_numero(k);
1240     k++;
1241     }
1242    
1243     LISTE_MG_NOEUD::iterator it_no;
1244     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1245     {
1246     int origine = noeud_i->get_origine();
1247     if (origine != IMPOSE)
1248     {
1249     int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
1250     double w_ij_prime = 0.0;
1251     OT_VECTEUR_3D v_temp(0.,0.,0.);
1252     OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
1253     for(int j=0;j<nb_voisins_j;j++)
1254     {
1255     MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
1256 francois 457 //On calcule le centroide cj du triangle mgtri_j
1257 cuillier 438 MG_NOEUD* n1 = mgtri_j->get_noeud1();
1258     MG_NOEUD* n2 = mgtri_j->get_noeud2();
1259     MG_NOEUD* n3 = mgtri_j->get_noeud3();
1260     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
1261     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
1262     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
1263     //On forme le vecteur vi_cj
1264     OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
1265     OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_j->get_nouveau_numero());
1266     v_temp = v_temp + (vi_cj*normal_f_i_new)*normal_f_i_new;
1267     }
1268 francois 457 //5-On met a jour la position des noeuds
1269 cuillier 438 v_i = v_i + (gamma/(2*nb_voisins_j))*v_temp;
1270     nouv_position_x.push_back(v_i.get_x());
1271     nouv_position_y.push_back(v_i.get_y());
1272     nouv_position_z.push_back(v_i.get_z());
1273     }
1274     }
1275    
1276     //On actualise la position des noeuds
1277     int ind_noeud=0;
1278     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1279     {
1280     int origine = noeud_i->get_origine();
1281     if (origine != IMPOSE)
1282     {
1283     noeud_i->change_x(nouv_position_x[ind_noeud]);
1284     noeud_i->change_y(nouv_position_y[ind_noeud]);
1285     noeud_i->change_z(nouv_position_z[ind_noeud]);
1286     ind_noeud++;
1287     }
1288     }
1289    
1290     //Critere d'arret de l'algorithme
1291     int l=0;
1292     int nb_tri = mg_mai->get_nb_mg_triangle();
1293     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
1294     {
1295     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
1296     OT_VECTEUR_3D normal_f_i = liste_normales2.get(mgtri_i->get_nouveau_numero());
1297     OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_i->get_nouveau_numero());
1298     double critere = 1. - normal_f_i*normal_f_i_new;
1299 picher 248 if (critere <= epsilon) l++;
1300     }
1301     double tolerance = 0.01*nb_tri;
1302     if (nb_tri - l <= 0) fin = 1;
1303     compteur++;
1304     }
1305     while ((fin == 0) && (compteur < iter_max));
1306    
1307 francois 457 return compteur;
1308    
1309 picher 248 }
1310    
1311 cuillier 438
1312 francois 457 //---------------------------- Extract skin -------------------------------------------------------------------------
1313     int MG_LISSAGE::extract_skin(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2,int &nbpeau,int &nbmaniare,int &nbmanino, int *mai2_id)
1314 cuillier 438 {
1315    
1316 francois 457 //etape 0 - On commence par mettre a zero tous les nouveau_numero des triangles et des noeuds du maillage
1317 francois 222 LISTE_MG_TRIANGLE::iterator it_tri;
1318     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
1319     {
1320     mgtri->change_nouveau_numero(0);
1321 francois 224 mgtri->change_origine(MAILLEUR_AUTO);
1322 francois 222 }
1323     LISTE_MG_NOEUD::iterator it_noeud;
1324     for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
1325     {
1326     mgnoeud->change_nouveau_numero(0);
1327     }
1328    
1329 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
1330 picher 233 LISTE_MG_TETRA::iterator it_tetra;
1331     for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
1332 francois 222 {
1333     int origine = mgtet->get_origine();
1334 francois 224 if (origine==IMPOSE)
1335 francois 222 {
1336 francois 224 mgtet->get_triangle1()->change_origine(IMPOSE);
1337     mgtet->get_triangle2()->change_origine(IMPOSE);
1338     mgtet->get_triangle3()->change_origine(IMPOSE);
1339     mgtet->get_triangle4()->change_origine(IMPOSE);
1340     }
1341     if (origine==OPTIMISE)
1342     {
1343     if (mgtet->get_triangle1()->get_origine()!=IMPOSE) mgtet->get_triangle1()->change_origine(OPTIMISE);
1344     if (mgtet->get_triangle2()->get_origine()!=IMPOSE) mgtet->get_triangle2()->change_origine(OPTIMISE);
1345     if (mgtet->get_triangle3()->get_origine()!=IMPOSE) mgtet->get_triangle3()->change_origine(OPTIMISE);
1346     if (mgtet->get_triangle4()->get_origine()!=IMPOSE) mgtet->get_triangle4()->change_origine(OPTIMISE);
1347     }
1348    
1349     if (((origine == OPTIMISE) && (etat[(OPTIMISE-1000)/10]==1) ) ||
1350     ((origine == IMPOSE) && (etat[(IMPOSE-1000)/10]==1) ) ||
1351     ((origine == MAILLEUR_AUTO) && (etat[(MAILLEUR_AUTO-1000)/10]==1) ) ||
1352     ((origine == TRIANGULATION) && (etat[(TRIANGULATION-1000)/10]==1) ) ||
1353     ((origine == MODIFICATION) && (etat[(MODIFICATION-1000)/10]==1) ) ||
1354     ((origine == DUPLIQUER) && (etat[(DUPLIQUER-1000)/10]==1) ) )
1355    
1356     {
1357 francois 222 int num1 = mgtet->get_triangle1()->get_nouveau_numero();
1358     int num2 = mgtet->get_triangle2()->get_nouveau_numero();
1359     int num3 = mgtet->get_triangle3()->get_nouveau_numero();
1360     int num4 = mgtet->get_triangle4()->get_nouveau_numero();
1361     num1++;
1362     num2++;
1363     num3++;
1364     num4++;
1365     mgtet->get_triangle1()->change_nouveau_numero(num1);
1366     mgtet->get_triangle2()->change_nouveau_numero(num2);
1367     mgtet->get_triangle3()->change_nouveau_numero(num3);
1368     mgtet->get_triangle4()->change_nouveau_numero(num4);
1369     }
1370     }
1371    
1372 francois 457 //etape 2 - On boucle l'ensemble des triangles identifies a l'etape 1 pour identifier les noeuds leur appartenant
1373 francois 222 for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
1374     {
1375     int num = mgtri->get_nouveau_numero();
1376     if (num == 1)
1377     {
1378     MG_NOEUD* noeud1 = mgtri->get_noeud1();
1379     MG_NOEUD* noeud2 = mgtri->get_noeud2();
1380     MG_NOEUD* noeud3 = mgtri->get_noeud3();
1381     noeud1->change_nouveau_numero(1);
1382     noeud2->change_nouveau_numero(1);
1383     noeud3->change_nouveau_numero(1);
1384 picher 248 if (mgtri->get_origine()==IMPOSE)
1385     {
1386     noeud1->change_origine(IMPOSE);
1387     noeud2->change_origine(IMPOSE);
1388     noeud3->change_origine(IMPOSE);
1389     }
1390 francois 222 }
1391     }
1392    
1393 francois 457 //etape 3 - On cree un nouveau maillage pour la peau
1394 francois 222 MG_MAILLAGE* mai2 = new MG_MAILLAGE(NULL);
1395     gest2.ajouter_mg_maillage(mai2);
1396    
1397 francois 457 //etape 4 - On boucle l'ensemble des noeuds identifies a l'etape 2 pour les recreer dans le second maillage
1398 francois 222 for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
1399     {
1400     int num = mgnoeud->get_nouveau_numero();
1401     if (num == 1)
1402     {
1403     double x = mgnoeud->get_x();
1404     double y = mgnoeud->get_y();
1405     double z = mgnoeud->get_z();
1406 francois 224 int origine=TRIANGULATION;
1407     if (mgnoeud->get_origine()==IMPOSE) origine=IMPOSE;
1408     MG_NOEUD* noeud1 = new MG_NOEUD(NULL,x,y,z,origine);
1409 francois 222 mai2->ajouter_mg_noeud(noeud1);
1410     mgnoeud->change_nouveau_numero(noeud1->get_id());
1411 francois 234 noeud1->change_nouveau_numero(mgnoeud->get_id());
1412 francois 222 }
1413     }
1414    
1415 francois 457 // etape 5 - On boucle l'ensemble des triangles identifies a l'etape 1 pour les recreer dans le maillage 2
1416 picher 233 for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
1417 francois 222 {
1418 francois 224 int originetet=mgtet->get_origine();
1419     if (((originetet == OPTIMISE) && (etat[(OPTIMISE-1000)/10]==1) ) ||
1420     ((originetet == IMPOSE) && (etat[(IMPOSE-1000)/10]==1) ) ||
1421     ((originetet == MAILLEUR_AUTO) && (etat[(MAILLEUR_AUTO-1000)/10]==1) ) ||
1422     ((originetet == TRIANGULATION) && (etat[(TRIANGULATION-1000)/10]==1) ) ||
1423     ((originetet == MODIFICATION) && (etat[(MODIFICATION-1000)/10]==1) ) ||
1424     ((originetet == DUPLIQUER) && (etat[(DUPLIQUER-1000)/10]==1) ) )
1425 francois 222 {
1426 francois 224 MG_NOEUD* noeud1 = mgtet->get_noeud1();
1427     MG_NOEUD* noeud2 = mgtet->get_noeud2();
1428     MG_NOEUD* noeud3 = mgtet->get_noeud3();
1429     MG_NOEUD* noeud4 = mgtet->get_noeud4();
1430 francois 222 MG_NOEUD* node1 = mai2->get_mg_noeudid(noeud1->get_nouveau_numero());
1431     MG_NOEUD* node2 = mai2->get_mg_noeudid(noeud2->get_nouveau_numero());
1432     MG_NOEUD* node3 = mai2->get_mg_noeudid(noeud3->get_nouveau_numero());
1433 francois 224 MG_NOEUD* node4 = mai2->get_mg_noeudid(noeud4->get_nouveau_numero());
1434     MG_TRIANGLE* tri1=mgtet->get_triangle1();
1435     MG_TRIANGLE* tri2=mgtet->get_triangle2();
1436     MG_TRIANGLE* tri3=mgtet->get_triangle3();
1437     MG_TRIANGLE* tri4=mgtet->get_triangle4();
1438     if (tri1->get_nouveau_numero()==1)
1439     {
1440     int origine=TRIANGULATION;
1441     if (tri1->get_origine()==IMPOSE) origine=IMPOSE;
1442     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node2,node3,mai2,origine);
1443     tripeau->change_nouveau_numero(0);
1444     tri1->change_nouveau_numero(0);
1445     }
1446     if (tri2->get_nouveau_numero()==1)
1447     {
1448     int origine=TRIANGULATION;
1449     if (tri2->get_origine()==IMPOSE) origine=IMPOSE;
1450     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node4,node2,mai2,origine);
1451     tripeau->change_nouveau_numero(0);
1452     tri2->change_nouveau_numero(0);
1453     }
1454     if (tri3->get_nouveau_numero()==1)
1455     {
1456     int origine=TRIANGULATION;
1457     if (tri3->get_origine()==IMPOSE) origine=IMPOSE;
1458     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node2,node4,node3,mai2,origine);
1459     tripeau->change_nouveau_numero(0);
1460     tri3->change_nouveau_numero(0);
1461     }
1462     if (tri4->get_nouveau_numero()==1)
1463     {
1464     int origine=TRIANGULATION;
1465     if (tri4->get_origine()==IMPOSE) origine=IMPOSE;
1466     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node3,node4,mai2,origine);
1467     tripeau->change_nouveau_numero(0);
1468     tri4->change_nouveau_numero(0);
1469     }
1470 francois 222 }
1471 francois 224 }
1472     // etape 6 recherche des vosins
1473     for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
1474     {
1475     MG_TRIANGLE_PEAU* tripeau=(MG_TRIANGLE_PEAU*)mgtri;
1476     MG_TRIANGLE_PEAU* voisin1=recherche_voisin(tripeau->get_noeud1(),tripeau->get_noeud2(),tripeau);
1477     MG_TRIANGLE_PEAU* voisin2=recherche_voisin(tripeau->get_noeud2(),tripeau->get_noeud3(),tripeau);
1478     MG_TRIANGLE_PEAU* voisin3=recherche_voisin(tripeau->get_noeud3(),tripeau->get_noeud1(),tripeau);
1479     tripeau->change_voisin1(voisin1);
1480     tripeau->change_voisin2(voisin2);
1481     tripeau->change_voisin3(voisin3);
1482     }
1483 francois 457 // etape 7 - On recherche les peaux
1484 francois 224 int fin;
1485     do
1486     {
1487     fin=1;
1488     for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
1489     {
1490     MG_TRIANGLE_PEAU *tripeau=(MG_TRIANGLE_PEAU*)mgtri;
1491     if (tripeau->get_nouveau_numero()==0)
1492     {
1493     fin=0;
1494     std::vector<MG_TRIANGLE_PEAU*> *peau=new std::vector<MG_TRIANGLE_PEAU*>;
1495     lst_peau.push_back(peau);
1496     tripeau->change_nouveau_numero(1);
1497     peau->push_back(tripeau);
1498     determine_peau(peau);
1499     }
1500     }
1501     }
1502     while (fin==0);
1503 francois 457 nbpeau=lst_peau.size();
1504 francois 224 //test de manifold
1505     LISTE_MG_SEGMENT::iterator itseg;
1506     for (MG_SEGMENT* seg=mai2->get_premier_segment(itseg);seg!=NULL;seg=mai2->get_suivant_segment(itseg))
1507     seg->change_nouveau_numero(0);
1508     TPL_MAP_ENTITE<MG_SEGMENT*> nonmanifoldarete;
1509     TPL_MAP_ENTITE<MG_NOEUD*> nonmanifoldnoeud;
1510 francois 234 TPL_MAP_ENTITE<MG_NOEUD*> nonmanifoldnoeuddepuisarete;
1511 francois 457 nbpeau=lst_peau.size();
1512 francois 224 for (int i=0;i<nbpeau;i++)
1513     {
1514     int nbtri=lst_peau[i]->size();
1515     for (int j=0;j<nbtri;j++)
1516     {
1517     MG_TRIANGLE_PEAU* tripeau=(*lst_peau[i])[j];
1518     tripeau->get_segment1()->change_nouveau_numero(tripeau->get_segment1()->get_nouveau_numero()+1);
1519     tripeau->get_segment2()->change_nouveau_numero(tripeau->get_segment2()->get_nouveau_numero()+1);
1520     tripeau->get_segment3()->change_nouveau_numero(tripeau->get_segment3()->get_nouveau_numero()+1);
1521     if (tripeau->get_segment1()->get_nouveau_numero()>2)
1522     nonmanifoldarete.ajouter(tripeau->get_segment1());
1523     if (tripeau->get_segment2()->get_nouveau_numero()>2)
1524     nonmanifoldarete.ajouter(tripeau->get_segment2());
1525     if (tripeau->get_segment3()->get_nouveau_numero()>2)
1526     nonmanifoldarete.ajouter(tripeau->get_segment3());
1527     }
1528     }
1529 picher 233
1530 francois 224 int nbnonmanifoldarete=nonmanifoldarete.get_nb();
1531     for (int i=0;i<nbnonmanifoldarete;i++)
1532     {
1533     MG_SEGMENT* seg=nonmanifoldarete.get(i);
1534     MG_NOEUD* n1=seg->get_noeud1();
1535     MG_NOEUD* n2=seg->get_noeud2();
1536 francois 234 nonmanifoldnoeuddepuisarete.ajouter(n1);
1537     nonmanifoldnoeuddepuisarete.ajouter(n2);
1538 francois 224 }
1539 francois 234 LISTE_MG_NOEUD::iterator itnoeud;
1540     for (MG_NOEUD* no=mai2->get_premier_noeud(itnoeud);no!=NULL;no=mai2->get_suivant_noeud(itnoeud))
1541     {
1542     if (nonmanifoldnoeuddepuisarete.existe(no)) continue;
1543     if (est_non_manifold(no)) nonmanifoldnoeud.ajouter(no);
1544     }
1545 francois 224
1546 francois 457 nbmaniare = nonmanifoldarete.get_nb();
1547     nbmanino = nonmanifoldnoeud.get_nb();
1548 picher 230
1549 francois 457 //etape 8 - Traitement des cas de non manifold
1550     //non manifold par arete
1551 picher 233 for (int i=0;i<nbnonmanifoldarete;i++)
1552     {
1553     MG_SEGMENT* segment=nonmanifoldarete.get(i);
1554 francois 234 MG_NOEUD* noeud1 = mg_mai->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
1555     MG_NOEUD* noeud2 = mg_mai->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
1556     int nb_tetra = noeud1->get_lien_tetra()->get_nb();
1557     for (int j=0;j<nb_tetra;j++)
1558 picher 233 {
1559 francois 234 MG_TETRA* mgtet =noeud1->get_lien_tetra()->get(j);
1560     int originetet=mgtet->get_origine();
1561     if (originetet == MAILLEUR_AUTO)
1562 picher 233 {
1563 francois 457 //On reactive le tetra si l'autre noeud du segment lui appartient aussi
1564 francois 234 MG_NOEUD* no1 = mgtet->get_noeud1();
1565     MG_NOEUD* no2 = mgtet->get_noeud2();
1566     MG_NOEUD* no3 = mgtet->get_noeud3();
1567     MG_NOEUD* no4 = mgtet->get_noeud4();
1568     if((no1 == noeud2) ||(no2 == noeud2) || (no3 == noeud2) || (no4 == noeud2))
1569     mgtet->change_origine(OPTIMISE);
1570 picher 233 }
1571     }
1572     }
1573    
1574     //non manifold par noeud
1575 picher 231 int nbnonmanifoldnoeud=nonmanifoldnoeud.get_nb();
1576     for (int i=0;i<nbnonmanifoldnoeud;i++)
1577     {
1578 francois 234 MG_NOEUD* noeud=mg_mai->get_mg_noeudid(nonmanifoldnoeud.get(i)->get_nouveau_numero());
1579     int nb_tetra = noeud->get_lien_tetra()->get_nb();
1580     for (int j = 0;j<nb_tetra;j++)
1581 francois 345 {
1582 francois 234 MG_TETRA* mgtet =noeud->get_lien_tetra()->get(j);
1583     int originetet=mgtet->get_origine();
1584     if (originetet == MAILLEUR_AUTO)
1585 picher 231 {
1586 francois 234 mgtet->change_origine(OPTIMISE);
1587 picher 231 }
1588     }
1589     }
1590 francois 234 *mai2_id = mai2->get_id();
1591     if ((nbnonmanifoldarete != 0) || (nbnonmanifoldnoeud != 0))
1592     {
1593     for (int i=0;i<lst_peau.size();i++)
1594     {
1595     delete lst_peau[i];
1596     }
1597     lst_peau.clear();
1598     gest2.supprimer_mg_maillageid(*mai2_id);
1599     return 0;
1600     }
1601 picher 230
1602 francois 234 return 1;
1603 francois 222 }
1604    
1605 francois 457 //---------------------------- Determine peau -------------------------------------------------------------------------
1606 francois 224 void MG_LISSAGE::determine_peau(std::vector<MG_TRIANGLE_PEAU*> * peau)
1607 francois 222 {
1608 francois 224 int num=0;
1609     while (num!=peau->size())
1610     {
1611     MG_TRIANGLE_PEAU* p=(*peau)[num];
1612     if (p->get_voisin1()->get_nouveau_numero()==0)
1613     {
1614     peau->push_back(p->get_voisin1());
1615     p->get_voisin1()->change_nouveau_numero(1);
1616     }
1617     if (p->get_voisin2()->get_nouveau_numero()==0)
1618     {
1619     peau->push_back(p->get_voisin2());
1620     p->get_voisin2()->change_nouveau_numero(1);
1621     }
1622     if (p->get_voisin3()->get_nouveau_numero()==0)
1623     {
1624     peau->push_back(p->get_voisin3());
1625     p->get_voisin3()->change_nouveau_numero(1);
1626     }
1627     num++;
1628     }
1629     }
1630    
1631 francois 457 //---------------------------- Insere triangle -------------------------------------------------------------------------
1632 francois 224 MG_TRIANGLE_PEAU* MG_LISSAGE::insere_triangle(MG_ELEMENT_TOPOLOGIQUE* topo,class MG_NOEUD *mgnoeud1,class MG_NOEUD *mgnoeud2,class MG_NOEUD *mgnoeud3,MG_MAILLAGE* mg_maillage,int origine)
1633     {
1634 francois 222 MG_TRIANGLE_PEAU* voisin1=NULL,*voisin2=NULL,*voisin3=NULL;
1635     MG_SEGMENT* mgsegment1=mg_maillage->get_mg_segment(mgnoeud1->get_id(),mgnoeud2->get_id());
1636     MG_SEGMENT* mgsegment2=mg_maillage->get_mg_segment(mgnoeud2->get_id(),mgnoeud3->get_id());
1637     MG_SEGMENT* mgsegment3=mg_maillage->get_mg_segment(mgnoeud3->get_id(),mgnoeud1->get_id());
1638     if (mgsegment1==NULL)
1639 francois 224 mgsegment1=mg_maillage->ajouter_mg_segment(topo,mgnoeud1,mgnoeud2,origine);
1640 francois 222 if (mgsegment2==NULL)
1641 francois 224 mgsegment2=mg_maillage->ajouter_mg_segment(topo,mgnoeud2,mgnoeud3,origine);
1642 francois 222 if (mgsegment3==NULL)
1643 francois 224 mgsegment3=mg_maillage->ajouter_mg_segment(topo,mgnoeud3,mgnoeud1,origine);
1644     MG_TRIANGLE_PEAU* mtriangle=new MG_TRIANGLE_PEAU(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,origine);
1645 francois 222 mg_maillage->ajouter_mg_triangle(mtriangle);
1646     return mtriangle;
1647     }
1648    
1649 francois 457 //---------------------------- Recherche voisin -------------------------------------------------------------------------
1650 francois 224 MG_TRIANGLE_PEAU* MG_LISSAGE::recherche_voisin(MG_NOEUD* mg_noeud1,MG_NOEUD* mg_noeud2,MG_TRIANGLE_PEAU* triref)
1651 francois 222 {
1652 francois 224 MG_TRIANGLE_PEAU* trisol=NULL;
1653     double angleref=4.*M_PI;
1654 francois 222 int nb1=mg_noeud1->get_lien_triangle()->get_nb();
1655     int nb2=mg_noeud2->get_lien_triangle()->get_nb();
1656     for (int i=0;i<nb1;i++)
1657     for (int j=0;j<nb2;j++)
1658 francois 224 if (mg_noeud1->get_lien_triangle()->get(i)==mg_noeud2->get_lien_triangle()->get(j))
1659     if ((MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i)!=triref)
1660     {
1661     double angle=calcul_angle(triref,(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i));
1662     if (angle<angleref)
1663     {
1664     angleref=angle;
1665     trisol=(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i);
1666     }
1667     }
1668     return trisol;
1669 francois 222 }
1670 francois 223
1671 francois 234 int MG_LISSAGE::est_non_manifold(MG_NOEUD* no)
1672     {
1673     static int compteur=0;
1674     compteur++;
1675     int nb_tri=no->get_lien_triangle()->get_nb();
1676     TPL_MAP_ENTITE<MG_TRIANGLE*> lst;
1677     for (int i=0;i<nb_tri;i++)
1678     lst.ajouter(no->get_lien_triangle()->get(i));
1679     MG_TRIANGLE_PEAU* premier_triangle=(MG_TRIANGLE_PEAU*)lst.get(0);
1680     lst.supprimer(premier_triangle);
1681     MG_TRIANGLE_PEAU* tricour=(MG_TRIANGLE_PEAU*)premier_triangle;
1682     int i=0;
1683     do
1684     {
1685     if (lst.existe(tricour->get_voisin1()) || ((tricour->get_voisin1()==premier_triangle)&& (i>1)))
1686     {
1687     tricour=tricour->get_voisin1();
1688     lst.supprimer(tricour);
1689     }
1690     else if (lst.existe(tricour->get_voisin2()) || ((tricour->get_voisin2()==premier_triangle)&& (i>1)))
1691     {
1692     tricour=tricour->get_voisin2();
1693     lst.supprimer(tricour);
1694     }
1695     else if (lst.existe(tricour->get_voisin3()) || ((tricour->get_voisin3()==premier_triangle)&& (i>1)))
1696     {
1697     tricour=tricour->get_voisin3();
1698     lst.supprimer(tricour);
1699     }
1700     i++;
1701     }
1702     while (tricour!=premier_triangle);
1703     if (lst.get_nb()>0)
1704     return 1;
1705     return 0;
1706     }
1707 francois 224
1708 francois 457 //---------------------------- Calcule angle -------------------------------------------------------------------------
1709 francois 224 double MG_LISSAGE::calcul_angle(MG_TRIANGLE_PEAU* ft1,MG_TRIANGLE_PEAU* ft2)
1710 francois 223 {
1711 francois 388 MG_NOEUD* noeud1=ft1->get_noeud1();
1712     MG_NOEUD* noeud2=ft1->get_noeud2();
1713     MG_NOEUD* noeud3=ft1->get_noeud3();
1714     MG_NOEUD* noeud4=ft2->get_noeud1();
1715     MG_NOEUD* noeud5=ft2->get_noeud2();
1716     MG_NOEUD* noeud6=ft2->get_noeud3();
1717     double angle=get_angle_par_noeud<MG_NOEUD*>(noeud1,noeud2,noeud3,noeud4,noeud5,noeud6);
1718     return angle;
1719 francois 223 }
1720 picher 230
1721 francois 457 //---------------------------- Ponderation gaussian -------------------------------------------------------------------------
1722 picher 230 double MG_LISSAGE::ponderation_gaussian(double s,double sigma)
1723     {
1724     double w_s;
1725     w_s = exp(-(s*s)/(2.*sigma*sigma));
1726     return w_s;
1727     }
1728    
1729 francois 457 //---------------------------- Ponderation Laplacian -------------------------------------------------------------------------
1730 picher 230 double MG_LISSAGE::ponderation_laplacian(double s,double sigma)
1731     {
1732     double w_s;
1733     w_s = exp(-(sqrt(2.)*fabs(s))/sigma);
1734     return w_s;
1735     }
1736    
1737 francois 457 //---------------------------- Ponderation elfallahford -------------------------------------------------------------------------
1738 picher 230 double MG_LISSAGE::ponderation_elfallahford(double s,double sigma)
1739     {
1740     double w_s;
1741     w_s = 1./sqrt(1+pow((s/sigma),2));
1742     return w_s;
1743     }