ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/aster/src/mgopt_posttraitement.cpp
Revision: 438
Committed: Thu Oct 17 21:10:11 2013 UTC (11 years, 6 months ago) by cuillier
Original Path: magic/lib/optimisation/src/mg_lissage.cpp
File size: 94270 byte(s)
Log Message:
Nouvel executable pour lisser le resultat optimisation topologie

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 222 }
26    
27    
28     MG_LISSAGE::~MG_LISSAGE()
29     {
30 francois 224 for (int i=0;i<lst_peau.size();i++)
31     delete lst_peau[i];
32 francois 222 }
33 picher 233
34 francois 232 void MG_LISSAGE::active_affichage(void (*fonc)(char*))
35     {
36     affiche=fonc;
37     affichageactif=1;
38     }
39 francois 222
40 picher 233
41 francois 224 void MG_LISSAGE::conserve(int origine)
42     {
43     etat[(origine-1000)/10]=1;
44     }
45    
46 picher 233
47     void MG_LISSAGE::copieorigine(FEM_MAILLAGE* mai,MG_GESTIONNAIRE& gest2)
48 picher 231 {
49 francois 309 LISTE_FEM_ELEMENT3::iterator it_tetra;
50     for (FEM_ELEMENT3* tet=mai->get_premier_element3(it_tetra);tet!=NULL;tet=mai->get_suivant_element3(it_tetra))
51 picher 233 {
52     MG_TETRA* mgtet=(MG_TETRA*)tet->get_mg_element_maillage();
53     mgtet->change_origine(mgtet->get_origine());
54     }
55     }
56     void MG_LISSAGE::gain_poids(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2)
57     {
58     LISTE_MG_TETRA::iterator it_tetra;
59     double volume_initial = 0.;
60     double volume_final = 0.;
61     for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
62     {
63     int originetet=mgtet->get_origine();
64     MG_NOEUD* noeud1 = mgtet->get_noeud1();
65     double x1 = noeud1->get_x();
66     double y1 = noeud1->get_y();
67     double z1 = noeud1->get_z();
68     MG_NOEUD* noeud2 = mgtet->get_noeud2();
69     double x2 = noeud2->get_x();
70     double y2 = noeud2->get_y();
71     double z2 = noeud2->get_z();
72     MG_NOEUD* noeud3 = mgtet->get_noeud3();
73     double x3 = noeud3->get_x();
74     double y3 = noeud3->get_y();
75     double z3 = noeud3->get_z();
76     MG_NOEUD* noeud4 = mgtet->get_noeud4();
77     double x4 = noeud4->get_x();
78     double y4 = noeud4->get_y();
79     double z4 = noeud4->get_z();
80     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));
81     volume_tet = volume_tet/6.;
82     volume_initial = volume_initial + volume_tet;
83     if (originetet != MAILLEUR_AUTO)
84     {
85     volume_final = volume_final + volume_tet;
86     }
87     }
88     cout << "volume_initial = " << volume_initial << endl;
89     cout << "volume_final = " << volume_final << endl;
90     }
91     void MG_LISSAGE::reactivation(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2)
92     {
93     cout << "Reactivation d'elements manquants sur aretes vives et faces planes" << endl;
94     LISTE_MG_TETRA::iterator it_tetra;
95     for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
96     {
97     int originetet=mgtet->get_origine();
98     int compteur = 0;
99     if (originetet == MAILLEUR_AUTO)
100     {
101     MG_TRIANGLE* tri1 = mgtet->get_triangle1();
102     MG_TRIANGLE* tri2 = mgtet->get_triangle2();
103     MG_TRIANGLE* tri3 = mgtet->get_triangle3();
104     MG_TRIANGLE* tri4 = mgtet->get_triangle4();
105     int nb_tet1 = tri1->get_lien_tetra()->get_nb();
106     for (int k = 0;k<nb_tet1;k++)
107     {
108     MG_TETRA* tet1 = tri1->get_lien_tetra()->get(k);
109     int ori1 = tet1->get_origine();
110     if ((tet1 != mgtet) && (ori1 == MAILLEUR_AUTO)) compteur++;
111     }
112     int nb_tet2 = tri2->get_lien_tetra()->get_nb();
113     for (int k = 0;k<nb_tet2;k++)
114     {
115     MG_TETRA* tet2 = tri2->get_lien_tetra()->get(k);
116     int ori2 = tet2->get_origine();
117     if ((tet2 != mgtet) && (ori2 == MAILLEUR_AUTO)) compteur++;
118     }
119     int nb_tet3 = tri3->get_lien_tetra()->get_nb();
120     for (int k = 0;k<nb_tet3;k++)
121     {
122     MG_TETRA* tet3 = tri3->get_lien_tetra()->get(k);
123     int ori3 = tet3->get_origine();
124     if ((tet3 != mgtet) && (ori3 == MAILLEUR_AUTO)) compteur++;
125     }
126     int nb_tet4 = tri4->get_lien_tetra()->get_nb();
127     for (int k = 0;k<nb_tet4;k++)
128     {
129     MG_TETRA* tet4 = tri4->get_lien_tetra()->get(k);
130     int ori4 = tet4->get_origine();
131     if ((tet4 != mgtet) && (ori4 == MAILLEUR_AUTO)) compteur++;
132     }
133 picher 248 if (compteur == 0) mgtet->change_origine(OPTIMISE);
134 picher 233 }
135     }
136     }
137     void MG_LISSAGE::visualisation2(MG_MAILLAGE* mai2,MG_GESTIONNAIRE& gest2,char *nomfichier,int id)
138     {
139     MG_MAILLAGE* mai3 = new MG_MAILLAGE(NULL);
140     MG_GESTIONNAIRE gest3;
141     gest3.ajouter_mg_maillage(mai3);
142     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
143     MG_NOEUD* noeud1 = mai2->get_mg_noeudid(id);
144     int nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
145     for (int j = 0;j<nb_voisins1;j++)
146     {
147     MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
148     liste_voisins.ajouter(mgtri_1);
149     MG_NOEUD* no1 = mgtri_1->get_noeud1();
150     int nb_vois1 = no1->get_lien_triangle()->get_nb();
151     for (int k = 0;k<nb_vois1;k++)
152     {
153     MG_TRIANGLE_PEAU* mgtri1 = (MG_TRIANGLE_PEAU*) no1->get_lien_triangle()->get(k);
154     liste_voisins.ajouter(mgtri1);
155     }
156     MG_NOEUD* no2 = mgtri_1->get_noeud2();
157     int nb_vois2 = no2->get_lien_triangle()->get_nb();
158     for (int k = 0;k<nb_vois2;k++)
159     {
160     MG_TRIANGLE_PEAU* mgtri2 = (MG_TRIANGLE_PEAU*) no2->get_lien_triangle()->get(k);
161     liste_voisins.ajouter(mgtri2);
162     }
163     MG_NOEUD* no3 = mgtri_1->get_noeud3();
164     int nb_vois3 = no3->get_lien_triangle()->get_nb();
165     for (int k = 0;k<nb_vois3;k++)
166     {
167     MG_TRIANGLE_PEAU* mgtri3 = (MG_TRIANGLE_PEAU*) no3->get_lien_triangle()->get(k);
168     liste_voisins.ajouter(mgtri3);
169     }
170     }
171     int nb_voisins = liste_voisins.get_nb();
172     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
173     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
174     {
175     MG_NOEUD* noeud1 = mgtri_j->get_noeud1();
176     MG_NOEUD* noeud2 = mgtri_j->get_noeud2();
177     MG_NOEUD* noeud3 = mgtri_j->get_noeud3();
178     mai3->ajouter_mg_noeud(noeud1);
179     mai3->ajouter_mg_noeud(noeud2);
180     mai3->ajouter_mg_noeud(noeud3);
181     mai3->ajouter_mg_triangle(mgtri_j);
182     }
183     gest3.enregistrer(nomfichier);
184     }
185    
186     void MG_LISSAGE::visualisation(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2,char *nomfichier)
187     {
188     //Routine pour visualiser le traitement des cas de non-manifold
189     MG_MAILLAGE* mai3 = new MG_MAILLAGE(NULL);
190     MG_GESTIONNAIRE gest3;
191     gest3.ajouter_mg_maillage(mai3);
192     MG_NOEUD* noeud_nm1 = mg_mai->get_mg_noeudid(41526);//placer ici le ID d'un noeud non manifold
193 francois 234 //MG_NOEUD* noeud_nm2 = mg_mai->get_mg_noeudid(12437);//placer ici le ID d'un second noeud non manifold (pour une ar�te)
194 picher 233 TPL_MAP_ENTITE<MG_TETRA*> liste_voisins;
195     int nb_voisins1 = noeud_nm1->get_lien_tetra()->get_nb();
196     for (int j = 0;j<nb_voisins1;j++)
197     {
198     MG_TETRA* mgtet_1 = (MG_TETRA*) noeud_nm1->get_lien_tetra()->get(j);
199     int origine = mgtet_1->get_origine();
200     if ((origine == OPTIMISE) || (origine == IMPOSE))
201     {
202     liste_voisins.ajouter(mgtet_1);
203     MG_NOEUD* no1 = mgtet_1->get_noeud1();
204     int nb_vois1 = no1->get_lien_tetra()->get_nb();
205     for (int k = 0;k<nb_vois1;k++)
206     {
207     MG_TETRA* mgtet1 = (MG_TETRA*) no1->get_lien_tetra()->get(k);
208     int origine1 = mgtet1->get_origine();
209     if ((origine1 == OPTIMISE) || (origine1 == IMPOSE)) liste_voisins.ajouter(mgtet1);
210     }
211     MG_NOEUD* no2 = mgtet_1->get_noeud2();
212     int nb_vois2 = no2->get_lien_tetra()->get_nb();
213     for (int k = 0;k<nb_vois2;k++)
214     {
215     MG_TETRA* mgtet2 = (MG_TETRA*) no2->get_lien_tetra()->get(k);
216     int origine2 = mgtet2->get_origine();
217     if ((origine2 == OPTIMISE) || (origine2 == IMPOSE)) liste_voisins.ajouter(mgtet2);
218     }
219     MG_NOEUD* no3 = mgtet_1->get_noeud3();
220     int nb_vois3 = no3->get_lien_tetra()->get_nb();
221     for (int k = 0;k<nb_vois3;k++)
222     {
223     MG_TETRA* mgtet3 = (MG_TETRA*) no3->get_lien_tetra()->get(k);
224     int origine3 = mgtet3->get_origine();
225     if ((origine3 == OPTIMISE) || (origine3 == IMPOSE)) liste_voisins.ajouter(mgtet3);
226     }
227     MG_NOEUD* no4 = mgtet_1->get_noeud4();
228     int nb_vois4 = no4->get_lien_tetra()->get_nb();
229     for (int k = 0;k<nb_vois4;k++)
230     {
231     MG_TETRA* mgtet4 = (MG_TETRA*) no4->get_lien_tetra()->get(k);
232     int origine4 = mgtet4->get_origine();
233     if ((origine4 == OPTIMISE) || (origine4 == IMPOSE)) liste_voisins.ajouter(mgtet4);
234     }
235     }
236     }
237     // int nb_voisins2 = noeud_nm2->get_lien_tetra()->get_nb();
238     // for (int j = 0;j<nb_voisins2;j++)
239     // {
240     // MG_TETRA* mgtet_1 = (MG_TETRA*) noeud_nm2->get_lien_tetra()->get(j);
241     // int origine = mgtet_1->get_origine();
242     // if ((origine == OPTIMISE) || (origine == IMPOSE))
243     // {
244     // liste_voisins.ajouter(mgtet_1);
245     // // MG_NOEUD* no1 = mgtet_1->get_noeud1();
246     // // int nb_vois1 = no1->get_lien_tetra()->get_nb();
247     // // for (int k = 0;k<nb_vois1;k++)
248     // // {
249     // // MG_TETRA* mgtet1 = (MG_TETRA*) no1->get_lien_tetra()->get(k);
250     // // int origine1 = mgtet1->get_origine();
251     // // if ((origine1 == OPTIMISE) || (origine1 == IMPOSE)) liste_voisins.ajouter(mgtet1);
252     // // }
253     // // MG_NOEUD* no2 = mgtet_1->get_noeud2();
254     // // int nb_vois2 = no2->get_lien_tetra()->get_nb();
255     // // for (int k = 0;k<nb_vois2;k++)
256     // // {
257     // // MG_TETRA* mgtet2 = (MG_TETRA*) no2->get_lien_tetra()->get(k);
258     // // int origine2 = mgtet2->get_origine();
259     // // if ((origine2 == OPTIMISE) || (origine2 == IMPOSE)) liste_voisins.ajouter(mgtet2);
260     // // }
261     // // MG_NOEUD* no3 = mgtet_1->get_noeud3();
262     // // int nb_vois3 = no3->get_lien_tetra()->get_nb();
263     // // for (int k = 0;k<nb_vois3;k++)
264     // // {
265     // // MG_TETRA* mgtet3 = (MG_TETRA*) no3->get_lien_tetra()->get(k);
266     // // int origine3 = mgtet3->get_origine();
267     // // if ((origine3 == OPTIMISE) || (origine3 == IMPOSE)) liste_voisins.ajouter(mgtet3);
268     // // }
269     // // MG_NOEUD* no4 = mgtet_1->get_noeud4();
270     // // int nb_vois4 = no4->get_lien_tetra()->get_nb();
271     // // for (int k = 0;k<nb_vois4;k++)
272     // // {
273     // // MG_TETRA* mgtet4 = (MG_TETRA*) no4->get_lien_tetra()->get(k);
274     // // int origine4 = mgtet4->get_origine();
275     // // if ((origine4 == OPTIMISE) || (origine4 == IMPOSE)) liste_voisins.ajouter(mgtet4);
276     // // }
277     // }
278     // }
279     int nb_voisins = liste_voisins.get_nb();
280     TPL_MAP_ENTITE<MG_TETRA*>::ITERATEUR it;
281     for(MG_TETRA* mgtet_j=liste_voisins.get_premier(it);mgtet_j!=NULL;mgtet_j=liste_voisins.get_suivant(it))
282     {
283     MG_NOEUD* noeud1 = mgtet_j->get_noeud1();
284     double xx1 = noeud1->get_x();
285     double yy1 = noeud1->get_y();
286     double zz1 = noeud1->get_z();
287     int ori_no1 = noeud1->get_origine();
288     MG_NOEUD* noeud11 = new MG_NOEUD(NULL,xx1,yy1,zz1,ori_no1);
289     MG_NOEUD* noeud2 = mgtet_j->get_noeud2();
290     double xx2 = noeud2->get_x();
291     double yy2 = noeud2->get_y();
292     double zz2 = noeud2->get_z();
293     int ori_no2 = noeud2->get_origine();
294     MG_NOEUD* noeud22 = new MG_NOEUD(NULL,xx2,yy2,zz2,ori_no2);
295     MG_NOEUD* noeud3 = mgtet_j->get_noeud3();
296     double xx3 = noeud3->get_x();
297     double yy3 = noeud3->get_y();
298     double zz3 = noeud3->get_z();
299     int ori_no3 = noeud3->get_origine();
300     MG_NOEUD* noeud33 = new MG_NOEUD(NULL,xx3,yy3,zz3,ori_no3);
301     MG_NOEUD* noeud4 = mgtet_j->get_noeud4();
302     double xx4 = noeud4->get_x();
303     double yy4 = noeud4->get_y();
304     double zz4 = noeud4->get_z();
305     int ori_no4 = noeud4->get_origine();
306     MG_NOEUD* noeud44 = new MG_NOEUD(NULL,xx4,yy4,zz4,ori_no4);
307     mai3->ajouter_mg_noeud(noeud11);
308     mai3->ajouter_mg_noeud(noeud22);
309     mai3->ajouter_mg_noeud(noeud33);
310     mai3->ajouter_mg_noeud(noeud44);
311     MG_TRIANGLE* tri1 = mgtet_j->get_triangle1();
312     MG_NOEUD* no11 = tri1->get_noeud1();
313     MG_NOEUD* no12 = tri1->get_noeud2();
314     MG_NOEUD* no13 = tri1->get_noeud3();
315     MG_SEGMENT* seg11 = tri1->get_segment1();
316     MG_SEGMENT* seg12 = tri1->get_segment2();
317     MG_SEGMENT* seg13 = tri1->get_segment3();
318     int ori_tri1 = tri1->get_origine();
319     MG_TRIANGLE* tri11 = new MG_TRIANGLE(NULL,no11,no12,no13,seg11,seg12,seg13,ori_tri1);
320     MG_TRIANGLE* tri2 = mgtet_j->get_triangle2();
321     MG_NOEUD* no21 = tri2->get_noeud1();
322     MG_NOEUD* no22 = tri2->get_noeud2();
323     MG_NOEUD* no23 = tri2->get_noeud3();
324     MG_SEGMENT* seg21 = tri2->get_segment1();
325     MG_SEGMENT* seg22 = tri2->get_segment2();
326     MG_SEGMENT* seg23 = tri2->get_segment3();
327     int ori_tri2 = tri2->get_origine();
328     MG_TRIANGLE* tri22 = new MG_TRIANGLE(NULL,no21,no22,no23,seg21,seg22,seg23,ori_tri2);
329     MG_TRIANGLE* tri3 = mgtet_j->get_triangle3();
330     MG_NOEUD* no31 = tri3->get_noeud1();
331     MG_NOEUD* no32 = tri3->get_noeud2();
332     MG_NOEUD* no33 = tri3->get_noeud3();
333     MG_SEGMENT* seg31 = tri3->get_segment1();
334     MG_SEGMENT* seg32 = tri3->get_segment2();
335     MG_SEGMENT* seg33 = tri3->get_segment3();
336     int ori_tri3 = tri3->get_origine();
337     MG_TRIANGLE* tri33 = new MG_TRIANGLE(NULL,no31,no32,no33,seg31,seg32,seg33,ori_tri3);
338     MG_TRIANGLE* tri4 = mgtet_j->get_triangle4();
339     MG_NOEUD* no41 = tri4->get_noeud1();
340     MG_NOEUD* no42 = tri4->get_noeud2();
341     MG_NOEUD* no43 = tri4->get_noeud3();
342     MG_SEGMENT* seg41 = tri4->get_segment1();
343     MG_SEGMENT* seg42 = tri4->get_segment2();
344     MG_SEGMENT* seg43 = tri4->get_segment3();
345     int ori_tri4 = tri4->get_origine();
346     MG_TRIANGLE* tri44 = new MG_TRIANGLE(NULL,no41,no42,no43,seg41,seg42,seg43,ori_tri4);
347     mai3->ajouter_mg_triangle(tri11);
348     mai3->ajouter_mg_triangle(tri22);
349     mai3->ajouter_mg_triangle(tri33);
350     mai3->ajouter_mg_triangle(tri44);
351     int ori_tet = mgtet_j->get_origine();
352     MG_TETRA* mgtet = new MG_TETRA(NULL,noeud11,noeud22,noeud33,noeud44,tri11,tri22,tri33,tri44,ori_tet);
353     mai3->ajouter_mg_tetra(mgtet);
354     }
355     gest3.enregistrer(nomfichier);
356     }
357    
358 francois 343
359    
360     void MG_LISSAGE::extract_skin_par_decoupage(FEM_SOLUTION* sol,double limit,MG_GESTIONNAIRE& gest2,int liss,double epsilon,double sigma,double gamma,int iter_max,double sigmaf,double sigmag,std::string nom)
361     {
362     affiche((char*)"Extraction de l'enveloppe par découpage");
363     sol->active_solution(0);
364     FEM_MAILLAGE* mai=sol->get_maillage();
365     affiche((char*)" Extrapolation de la densité aux noeuds");
366     LISTE_FEM_NOEUD::iterator it;
367     for (FEM_NOEUD* no=mai->get_premier_noeud(it);no!=NULL;no=mai->get_suivant_noeud(it))
368     {
369     double nume=0.;
370     double deno=0.;
371     int nb=no->get_lien_element3()->get_nb();
372     int passe=0;
373     for (int i=0;i<nb;i++)
374     {
375     FEM_ELEMENT3* tet=no->get_lien_element3()->get(i);
376     if (tet->get_mg_element_maillage()->get_origine()!=IMPOSE)
377     {
378     double jac[9];
379     int col,ligne;
380     double volume=tet->get_jacobien(jac,jac,ligne,col,1.);
381     passe=1;
382     nume=nume+tet->get_solution()*volume;
383     deno=deno+volume;
384     }
385     }
386     if (passe==1) no->change_solution(nume/deno); else no->change_solution(1.);
387     }
388     if (nom!="")
389     {
390     affiche((char*)" Enregistrement de la densité aux noeuds");
391     MG_GESTIONNAIRE *gest=mai->get_mg_maillage()->get_gestionnaire();
392     std::string chemin=nom+".sol";
393 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);
394 francois 343 gest->ajouter_fem_solution(femdens);
395     femdens->change_legende(0,"Densite");
396     LISTE_FEM_ELEMENT3::iterator it3;
397     int i=0;
398     for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
399     {
400     MG_TETRA* tet2=(MG_TETRA*)tet->get_mg_element_maillage();
401     if (tet2->get_origine()==IMPOSE)
402     {
403 francois 375 femdens->ecrire(1.,i,0,0);
404     femdens->ecrire(1.,i,0,1);
405     femdens->ecrire(1.,i,0,2);
406     femdens->ecrire(1.,i,0,3);
407 francois 343 }
408     else
409     {
410 francois 375 femdens->ecrire(tet->get_fem_noeud(0)->get_solution(),i,0,0);
411     femdens->ecrire(tet->get_fem_noeud(1)->get_solution(),i,0,1);
412     femdens->ecrire(tet->get_fem_noeud(2)->get_solution(),i,0,2);
413     femdens->ecrire(tet->get_fem_noeud(3)->get_solution(),i,0,3);
414 francois 343 }
415     i++;
416     }
417     MG_EXPORT exp;
418     exp.gmsh(mai,nom);
419     gest->supprimer_fem_solution(gest->get_nb_fem_solution()-1);
420     }
421     MG_MAILLAGE* mgmai=new MG_MAILLAGE(NULL);
422     gest2.ajouter_mg_maillage(mgmai);
423    
424     MG_MAILLAGE* mgfem=mai->get_mg_maillage();
425     LISTE_MG_TRIANGLE::iterator itmg;
426     for (MG_TRIANGLE* tri=mgfem->get_premier_triangle(itmg);tri!=NULL;tri=mgfem->get_suivant_triangle(itmg))
427     tri->change_nouveau_numero(0);
428     affiche((char*)" Decoupage des tetra optimises");
429     LISTE_FEM_ELEMENT3::iterator it3;
430     for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
431     {
432     MG_TETRA* tet2=(MG_TETRA*)tet->get_mg_element_maillage();
433     if (tet2->get_origine()==IMPOSE)
434     {
435     tet2->get_triangle1()->change_origine(IMPOSE);
436     tet2->get_triangle2()->change_origine(IMPOSE);
437     tet2->get_triangle3()->change_origine(IMPOSE);
438     tet2->get_triangle4()->change_origine(IMPOSE);
439     tet2->get_noeud1()->change_solution(tet->get_fem_noeud(0)->get_solution());
440     tet2->get_noeud2()->change_solution(tet->get_fem_noeud(1)->get_solution());
441     tet2->get_noeud3()->change_solution(tet->get_fem_noeud(2)->get_solution());
442     tet2->get_noeud4()->change_solution(tet->get_fem_noeud(3)->get_solution());
443     tet2->get_noeud1()->change_nouveau_numero(tet->get_fem_noeud(0)->get_id());
444     tet2->get_noeud2()->change_nouveau_numero(tet->get_fem_noeud(1)->get_id());
445     tet2->get_noeud3()->change_nouveau_numero(tet->get_fem_noeud(2)->get_id());
446     tet2->get_noeud4()->change_nouveau_numero(tet->get_fem_noeud(3)->get_id());
447     if (tet2->get_triangle1()->get_nouveau_numero()==0)
448     tet2->get_triangle1()->change_nouveau_numero(tet2->get_noeud4()->get_id());
449     else
450     tet2->get_triangle1()->change_nouveau_numero(-1);
451     if (tet2->get_triangle2()->get_nouveau_numero()==0)
452     tet2->get_triangle2()->change_nouveau_numero(tet2->get_noeud3()->get_id());
453     else
454     tet2->get_triangle2()->change_nouveau_numero(-1);
455     if (tet2->get_triangle3()->get_nouveau_numero()==0)
456     tet2->get_triangle3()->change_nouveau_numero(tet2->get_noeud1()->get_id());
457     else
458     tet2->get_triangle3()->change_nouveau_numero(-1);
459     if (tet2->get_triangle4()->get_nouveau_numero()==0)
460     tet2->get_triangle4()->change_nouveau_numero(tet2->get_noeud2()->get_id());
461     else
462     tet2->get_triangle4()->change_nouveau_numero(-1);
463     }
464     }
465    
466     TPL_LISTE_ENTITE<MG_TRIANGLE*> tri_impose_interne;
467     for (MG_TRIANGLE* tri=mgfem->get_premier_triangle(itmg);tri!=NULL;tri=mgfem->get_suivant_triangle(itmg))
468     {
469     if (tri->get_origine()==IMPOSE)
470     if (tri->get_lien_topologie()->get_dimension()==3)
471     if (tri->get_nouveau_numero()>0)
472     tri_impose_interne.ajouter(tri);
473     }
474     for (FEM_ELEMENT3* tet=mai->get_premier_element3(it3);tet!=NULL;tet=mai->get_suivant_element3(it3))
475     {
476     if (tet->get_mg_element_maillage()->get_origine()!=IMPOSE)
477     {
478     FEM_NOEUD* no1=tet->get_fem_noeud(0);
479     FEM_NOEUD* no2=tet->get_fem_noeud(1);
480     FEM_NOEUD* no3=tet->get_fem_noeud(2);
481     FEM_NOEUD* no4=tet->get_fem_noeud(3);
482     std::vector<MG_NOEUD*> tab;
483     interpole_segment(no1,no2,&tab,limit,mgmai);
484     interpole_segment(no1,no3,&tab,limit,mgmai);
485     interpole_segment(no1,no4,&tab,limit,mgmai);
486     interpole_segment(no2,no3,&tab,limit,mgmai);
487     interpole_segment(no2,no4,&tab,limit,mgmai);
488     interpole_segment(no3,no4,&tab,limit,mgmai);
489     int nb=tab.size();
490     FEM_NOEUD* noext;
491     if (nb>0)
492     {
493     if (no1->get_solution()<limit) noext=no1;
494     if (no2->get_solution()<limit) noext=no2;
495     if (no3->get_solution()<limit) noext=no3;
496     if (no4->get_solution()<limit) noext=no4;
497     }
498     if (nb==3)
499     {
500     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[1],tab[2],mgmai,TRIANGULATION);
501     oriente_tri(tri,noext->get_coord());
502     }
503     if (nb==4)
504     {
505     if (test_du_point_milieu(tab[0],tab[1],tet)==1)
506     {
507     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[1],tab[2],mgmai,TRIANGULATION);
508     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[0],tab[1],tab[3],mgmai,TRIANGULATION);
509     oriente_tri(tri,noext->get_coord());
510     oriente_tri(tri2,noext->get_coord());
511     }
512     else if (test_du_point_milieu(tab[0],tab[2],tet)==1)
513     {
514     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[2],tab[1],mgmai,TRIANGULATION);
515     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[0],tab[2],tab[3],mgmai,TRIANGULATION);
516     oriente_tri(tri,noext->get_coord());
517     oriente_tri(tri2,noext->get_coord());}
518     else if (test_du_point_milieu(tab[0],tab[3],tet)==1)
519     {
520     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[3],tab[1],mgmai,TRIANGULATION);
521     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[0],tab[3],tab[2],mgmai,TRIANGULATION);
522     oriente_tri(tri,noext->get_coord());
523     oriente_tri(tri2,noext->get_coord());}
524     else if (test_du_point_milieu(tab[1],tab[2],tet)==1)
525     {
526     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[1],tab[2],tab[0],mgmai,TRIANGULATION);
527     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[1],tab[2],tab[3],mgmai,TRIANGULATION);
528     oriente_tri(tri,noext->get_coord());
529     oriente_tri(tri2,noext->get_coord());
530     }
531     else if (test_du_point_milieu(tab[1],tab[3],tet)==1)
532     {
533     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[1],tab[3],tab[0],mgmai,TRIANGULATION);
534     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[1],tab[3],tab[2],mgmai,TRIANGULATION);
535     oriente_tri(tri,noext->get_coord());
536     oriente_tri(tri2,noext->get_coord());}
537     else
538     {
539     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[2],tab[3],tab[0],mgmai,TRIANGULATION);
540     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[2],tab[3],tab[1],mgmai,TRIANGULATION);
541     oriente_tri(tri,noext->get_coord());
542     oriente_tri(tri2,noext->get_coord());
543     }
544    
545     }
546     }
547     }
548     affiche((char*)" Decoupage des triangles non design interne au volume");
549     for (int i=0;i<tri_impose_interne.get_nb();i++)
550     {
551     MG_TRIANGLE* tri=tri_impose_interne.get(i);
552     MG_NOEUD* no1=tri->get_noeud1();
553     MG_NOEUD* no2=tri->get_noeud2();
554     MG_NOEUD* no3=tri->get_noeud3();
555     std::vector<MG_NOEUD*> tab;
556     int num1=-1;
557     int num2=-1;
558     int num3=-1;
559     if (no1->get_solution()<limit)
560     {
561     MG_NOEUD* no=get_noeud_peau(mai->get_fem_noeudid(no1->get_nouveau_numero()),mgmai);
562     tab.push_back(no);
563     num1=1;
564     }
565     if (no2->get_solution()<limit)
566     {
567     MG_NOEUD* no=get_noeud_peau(mai->get_fem_noeudid(no2->get_nouveau_numero()),mgmai);
568     tab.push_back(no);
569     num2=1;
570     }
571     if (no3->get_solution()<limit)
572     {
573     MG_NOEUD* no=get_noeud_peau(mai->get_fem_noeudid(no3->get_nouveau_numero()),mgmai);
574     tab.push_back(no);
575     num3=1;
576     }
577     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);
578     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);
579     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);
580     MG_NOEUD* noint=mgfem->get_mg_noeudid(tri->get_nouveau_numero());
581     int nb=tab.size();
582     if (nb==3)
583     {
584     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[1],tab[2],mgmai,IMPOSE);
585     oriente_tri(tri,noint->get_coord());
586     tri->inverse_sens();
587     }
588     if (nb==4)
589     {
590     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[1],tab[2],mgmai,IMPOSE);
591     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[1],tab[2],tab[3],mgmai,IMPOSE);
592     oriente_tri(tri,noint->get_coord());
593     tri->inverse_sens();
594     oriente_tri(tri2,noint->get_coord());
595     tri2->inverse_sens();
596     }
597     }
598 francois 345 affiche((char*)" Decoupage des triangles non design sur la frontiere du volume");
599 francois 343 LISTE_FEM_ELEMENT2::iterator it2;
600     for (FEM_ELEMENT2* tri=mai->get_premier_element2(it2);tri!=NULL;tri=mai->get_suivant_element2(it2))
601     {
602     std::vector<MG_NOEUD*> tab;
603     FEM_NOEUD *no1=tri->get_fem_noeud(0);
604     FEM_NOEUD *no2=tri->get_fem_noeud(1);
605     FEM_NOEUD *no3=tri->get_fem_noeud(2);
606     int ori=((MG_FACE*)tri->get_lien_topologie())->get_mg_coface(0)->get_orientation();
607     OT_VECTEUR_3D n1n2(no1->get_coord(),no2->get_coord());
608     OT_VECTEUR_3D n1n3(no1->get_coord(),no3->get_coord());
609     OT_VECTEUR_3D nor=n1n3&n1n2;
610     double xyzext[3];
611     xyzext[0]=no1->get_x()+nor.get_x()*ori;
612     xyzext[1]=no1->get_y()+nor.get_y()*ori;
613     xyzext[2]=no1->get_z()+nor.get_z()*ori;
614     if (tri->get_mg_element_maillage()->get_origine()==IMPOSE)
615     {
616     MG_NOEUD *mgno1=get_noeud_peau(no1,mgmai);
617     MG_NOEUD *mgno2=get_noeud_peau(no2,mgmai);
618     MG_NOEUD *mgno3=get_noeud_peau(no3,mgmai);
619     int num1=interpole_segment(no1,no2,&tab,limit,mgmai,0);
620     int num2=interpole_segment(no1,no3,&tab,limit,mgmai,0);
621     int num3=interpole_segment(no2,no3,&tab,limit,mgmai,0);
622     int nb=tab.size();
623     if (nb==0)
624     {
625     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,mgno1,mgno2,mgno3,mgmai,IMPOSE);
626     oriente_tri(tri,xyzext);
627     }
628     if (nb==1)
629     {
630     if (num1==1)
631     {
632     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,mgno1,mgno3,tab[0],mgmai,IMPOSE);
633     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,mgno2,mgno3,tab[0],mgmai,IMPOSE);
634     oriente_tri(tri,xyzext);
635     oriente_tri(tri2,xyzext);
636     }
637     if (num2==1)
638     {
639     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,mgno1,mgno2,tab[0],mgmai,IMPOSE);
640     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,mgno2,mgno3,tab[0],mgmai,IMPOSE);
641     oriente_tri(tri,xyzext);
642     oriente_tri(tri2,xyzext);
643     }
644     if (num3==1)
645     {
646     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,mgno1,mgno2,tab[0],mgmai,IMPOSE);
647     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,mgno1,mgno3,tab[0],mgmai,IMPOSE);
648     oriente_tri(tri,xyzext);
649     oriente_tri(tri2,xyzext);
650     }
651     }
652     if (nb==2)
653     {
654     if (num1==0)
655     {
656     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,mgno1,mgno2,tab[0],mgmai,IMPOSE);
657     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,mgno2,tab[0],tab[1],mgmai,IMPOSE);
658     MG_TRIANGLE_PEAU* tri3=insere_triangle(NULL,mgno3,tab[0],tab[1],mgmai,IMPOSE);
659     oriente_tri(tri,xyzext);
660     oriente_tri(tri2,xyzext);
661     oriente_tri(tri3,xyzext);
662     }
663     if (num2==0)
664     {
665     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,mgno1,mgno3,tab[0],mgmai,IMPOSE);
666     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,mgno3,tab[0],tab[1],mgmai,IMPOSE);
667     MG_TRIANGLE_PEAU* tri3=insere_triangle(NULL,mgno2,tab[0],tab[1],mgmai,IMPOSE);
668     oriente_tri(tri,xyzext);
669     oriente_tri(tri2,xyzext);
670     oriente_tri(tri3,xyzext);
671     }
672     if (num2==0)
673     {
674     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,mgno2,mgno3,tab[0],mgmai,IMPOSE);
675     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,mgno3,tab[0],tab[1],mgmai,IMPOSE);
676     MG_TRIANGLE_PEAU* tri3=insere_triangle(NULL,mgno1,tab[0],tab[1],mgmai,IMPOSE);
677     oriente_tri(tri,xyzext);
678     oriente_tri(tri2,xyzext);
679     oriente_tri(tri3,xyzext);
680     }
681     }
682     }
683     else if ((no1->get_solution()>=limit) && (no2->get_solution()>=limit) && (no3->get_solution()>=limit))
684     {
685     MG_NOEUD *mgno1=get_noeud_peau(no1,mgmai);
686     MG_NOEUD *mgno2=get_noeud_peau(no2,mgmai);
687     MG_NOEUD *mgno3=get_noeud_peau(no3,mgmai);
688     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,mgno1,mgno2,mgno3,mgmai,TRIANGULATION);
689     oriente_tri(tri,xyzext);
690     }
691     else if (!((no1->get_solution()<limit) && (no2->get_solution()<limit) && (no3->get_solution()<limit)))
692     {
693     std::vector<MG_NOEUD*> tab;
694     int num1=-1;
695     int num2=-1;
696     int num3=-1;
697     if (no1->get_solution()>=limit)
698     {
699     MG_NOEUD* no=get_noeud_peau(no1,mgmai);
700     tab.push_back(no);
701     num1=1;
702     }
703     if (no2->get_solution()>=limit)
704     {
705     MG_NOEUD* no=get_noeud_peau(no2,mgmai);
706     tab.push_back(no);
707     num2=1;
708     }
709     if (no3->get_solution()>=limit)
710     {
711     MG_NOEUD* no=get_noeud_peau(no3,mgmai);
712     tab.push_back(no);
713     num3=1;
714     }
715     if (num1*num2<0) interpole_segment(no1,no2,&tab,limit,mgmai);
716     if (num1*num3<0) interpole_segment(no1,no3,&tab,limit,mgmai);
717     if (num2*num3<0) interpole_segment(no2,no3,&tab,limit,mgmai);
718     int nb=tab.size();
719     if (nb==3)
720     {
721     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[1],tab[2],mgmai,TRIANGULATION);
722     oriente_tri(tri,xyzext);
723     }
724     if (nb==4)
725     {
726     MG_TRIANGLE_PEAU* tri=insere_triangle(NULL,tab[0],tab[1],tab[2],mgmai,TRIANGULATION);
727     MG_TRIANGLE_PEAU* tri2=insere_triangle(NULL,tab[1],tab[2],tab[3],mgmai,TRIANGULATION);
728     oriente_tri(tri,xyzext);
729     oriente_tri(tri2,xyzext);
730     }
731     }
732     }
733     affiche((char*)" Suppression des peaux isolées");
734     LISTE_MG_TRIANGLE::iterator it_tri;
735     for (MG_TRIANGLE* mgtri=mgmai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mgmai->get_suivant_triangle(it_tri))
736     {
737     MG_TRIANGLE_PEAU* tripeau=(MG_TRIANGLE_PEAU*)mgtri;
738     MG_TRIANGLE_PEAU* voisin1=recherche_voisin(tripeau->get_noeud1(),tripeau->get_noeud2(),tripeau);
739     MG_TRIANGLE_PEAU* voisin2=recherche_voisin(tripeau->get_noeud2(),tripeau->get_noeud3(),tripeau);
740     MG_TRIANGLE_PEAU* voisin3=recherche_voisin(tripeau->get_noeud3(),tripeau->get_noeud1(),tripeau);
741     tripeau->change_voisin1(voisin1);
742     tripeau->change_voisin2(voisin2);
743     tripeau->change_voisin3(voisin3);
744     tripeau->change_nouveau_numero(0);
745     }
746     // //�tape 7 - On recherche les peaux
747     int fin;
748     do
749     {
750     fin=1;
751     for (MG_TRIANGLE* mgtri=mgmai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mgmai->get_suivant_triangle(it_tri))
752     {
753     MG_TRIANGLE_PEAU *tripeau=(MG_TRIANGLE_PEAU*)mgtri;
754     if (tripeau->get_nouveau_numero()==0)
755     {
756     fin=0;
757     std::vector<MG_TRIANGLE_PEAU*> *peau=new std::vector<MG_TRIANGLE_PEAU*>;
758     lst_peau.push_back(peau);
759     tripeau->change_nouveau_numero(1);
760     peau->push_back(tripeau);
761     determine_peau(peau);
762     }
763     }
764     }
765     while (fin==0);
766    
767     char message[500];
768     for (int cas=0;cas<2;cas++)
769     {
770     if (cas==0) affiche("Analyse initiale des peaux");
771     if (cas==1) affiche("Analyse finale des peaux");
772     sprintf(message," %d peaux",lst_peau.size());
773     affiche(message);
774     for (int i=0;i<lst_peau.size();i++)
775     {
776     sprintf(message," peau %d : %d triangles ",i+1,lst_peau[i]->size());
777     affiche(message);
778     int isole=1;
779     for (int j=0;j<lst_peau[i]->size();j++)
780     if ((*lst_peau[i])[j]->get_origine()==IMPOSE) {isole=0;break;}
781     if (isole==1)
782     {
783     affiche(" peau isolee");
784     for (int j=0;j<lst_peau[i]->size();j++)
785     {
786     mgmai->supprimer_mg_triangleid((*lst_peau[i])[j]->get_id());
787     }
788     lst_peau[i]->clear();
789     }
790     else affiche(" peau non isolee");
791     }
792     for (std::vector<std::vector<MG_TRIANGLE_PEAU*> *>::iterator j=lst_peau.end()-1;j>lst_peau.begin();j--)
793     if ((*j)->size()==0)
794     lst_peau.erase(j);
795     }
796     for (MG_TRIANGLE* tri=mgmai->get_premier_triangle(itmg);tri!=NULL;tri=mgmai->get_suivant_triangle(itmg))
797     {
798     if (tri->get_origine()==IMPOSE)
799     {
800     tri->get_noeud1()->change_origine(IMPOSE);
801     tri->get_noeud2()->change_origine(IMPOSE);
802     tri->get_noeud3()->change_origine(IMPOSE);
803     }
804     }
805     if (liss == 1)
806     {
807     affiche((char*)"Procedure de lissage 1 ");
808     lissage(mgmai,gest2,epsilon,sigma,iter_max);
809     }
810     if (liss == 2)
811     {
812     affiche((char*)"Procedure de lissage 2");
813     lissage2(mgmai,gest2,sigmaf,sigmag,iter_max);
814     }
815     if (liss == 3)
816     {
817     affiche((char*)"Procedure de lissage 3");
818     lissage3(mgmai,gest2,sigma,gamma,epsilon,iter_max);
819     }
820 cuillier 438 if (liss == 4)
821     {
822     affiche((char*)"Procedure de lissage 4");
823     lissage4(mgmai,gest2,iter_max);
824     }
825     if (liss == 5)
826     {
827     affiche((char*)"Procedure de lissage 5");
828     lissage5(mgmai,gest2,iter_max);
829     }
830     if (liss == 6)
831     {
832     affiche((char*)"Procedure de lissage 6");
833     lissage6(mgmai,gest2,sigma,gamma,epsilon,iter_max);
834     }
835     if (liss == 7)
836     {
837     affiche((char*)"Procedure de lissage 7");
838     lissage7(mgmai,gest2,iter_max);
839     }
840    
841 francois 343 }
842    
843     void MG_LISSAGE::oriente_tri(MG_TRIANGLE_PEAU* tri,double *xyz)
844     {
845     MG_NOEUD* no1=tri->get_noeud1();
846     MG_NOEUD* no2=tri->get_noeud2();
847     MG_NOEUD* no3=tri->get_noeud3();
848     OT_VECTEUR_3D n1n2(no1->get_coord(),no2->get_coord());
849     OT_VECTEUR_3D n1n3(no1->get_coord(),no3->get_coord());
850     OT_VECTEUR_3D normal=n1n2&n1n3;
851     OT_VECTEUR_3D dir(no1->get_coord(),xyz);
852     dir.norme();
853     normal.norme();
854     double ps=normal*dir;
855     if (ps<0) tri->inverse_sens();
856     }
857    
858    
859    
860     MG_NOEUD* MG_LISSAGE::get_noeud_peau(FEM_NOEUD* no,MG_MAILLAGE* mai)
861     {
862     static std::map<std::string,MG_NOEUD*> correspond;
863     unsigned long id=no->get_id();
864     char message[255];
865     sprintf(message,"_%lu_",id);
866     std::string key=message;
867     MG_NOEUD* mgno=correspond[key];
868     if (mgno==NULL)
869     {
870     mgno=new MG_NOEUD(NULL,no->get_x(),no->get_y(),no->get_z(),TRIANGULATION);
871     mai->ajouter_mg_noeud(mgno);
872     correspond[key]=mgno;
873     }
874     return mgno;
875     }
876    
877    
878     int MG_LISSAGE::test_du_point_milieu(MG_NOEUD* no1,MG_NOEUD* no2,FEM_ELEMENT3* tet)
879     {
880     double *xyz1=no1->get_coord();
881     double *xyz2=no2->get_coord();
882     double xyz[3];
883     xyz[0]=0.5*(xyz1[0]+xyz2[0]);
884     xyz[1]=0.5*(xyz1[1]+xyz2[1]);
885     xyz[2]=0.5*(xyz1[2]+xyz2[2]);
886 francois 412 FEM_MAILLAGE_OUTILS ot;
887 francois 388 if (((ot.estdansletetra(tet,xyz[0],xyz[1],xyz[2])>>1)&1)==1) return 1;
888     //if (((MG_MAILLAGE_OUTILS::estdansletetra(tet,xyz[0],xyz[1],xyz[2])>>1)&1)==1) return 1;
889 francois 343 return 0;
890     }
891    
892    
893    
894     int MG_LISSAGE::interpole_segment(FEM_NOEUD* no1,FEM_NOEUD* no2,std::vector<MG_NOEUD*> *tab,double limit,MG_MAILLAGE* mai,int creation)
895     {
896     static std::map<std::string,MG_NOEUD*> correspond;
897     unsigned long id1=no1->get_id();
898     unsigned long id2=no2->get_id();
899     char message[255];
900     sprintf(message,"_%lu_%lu_",id1,id2);
901     std::string key1=message;
902     sprintf(message,"_%lu_%lu_",id2,id1);
903     std::string key2=message;
904     double sol1=no1->get_solution();
905     double sol2=no2->get_solution();
906 francois 345 if (fabs(sol1-limit)<0.0000001)
907     {
908     sprintf(message,"%lu",id1);
909     std::string key=message;
910     MG_NOEUD* no=correspond[key];
911     if (no==NULL)
912     {
913     double *xyz1=no1->get_coord();
914     no=new MG_NOEUD(NULL,xyz1[0],xyz1[1],xyz1[2],TRIANGULATION);
915     mai->ajouter_mg_noeud(no);
916     correspond[key]=no;
917     }
918     int present=0;
919     for (int i=0;i<tab->size();i++)
920     if ((*tab)[i]==no) present=1;
921     if (present==0)
922     {
923     tab->push_back(no);
924     return 1;
925     }
926     else return 0;
927     }
928     if (fabs(sol2-limit)<0.0000001)
929     {
930     sprintf(message,"%lu",id2);
931     std::string key=message;
932     MG_NOEUD* no=correspond[key];
933     if (no==NULL)
934     {
935     double *xyz2=no2->get_coord();
936     no=new MG_NOEUD(NULL,xyz2[0],xyz2[1],xyz2[2],TRIANGULATION);
937     mai->ajouter_mg_noeud(no);
938     correspond[key]=no;
939     }
940     int present=0;
941     for (int i=0;i<tab->size();i++)
942     if ((*tab)[i]==no) present=1;
943     if (present==0)
944     {
945     tab->push_back(no);
946     return 1;
947     }
948     else return 0;
949    
950     }
951    
952     if (((sol1<limit) && (sol2>limit))||((sol1>limit) && (sol2<limit)))
953 francois 343 {
954     MG_NOEUD* no=correspond[key1];
955     if (no==NULL) no=correspond[key2];
956     if ((no==NULL) && (creation==1))
957     {
958     double t=(limit-sol1)/(sol2-sol1);
959     double xyz[3];
960     double *xyz1=no1->get_coord();
961     double *xyz2=no2->get_coord();
962     xyz[0]=xyz1[0]+t*(xyz2[0]-xyz1[0]);
963     xyz[1]=xyz1[1]+t*(xyz2[1]-xyz1[1]);
964     xyz[2]=xyz1[2]+t*(xyz2[2]-xyz1[2]);
965     no=new MG_NOEUD(NULL,xyz[0],xyz[1],xyz[2],TRIANGULATION);
966     mai->ajouter_mg_noeud(no);
967     correspond[key1]=no;
968     }
969     if (no!=NULL)
970     {
971     tab->push_back(no);
972     return 1;
973     }
974     }
975     return 0;
976     }
977    
978    
979    
980 picher 248 void MG_LISSAGE::lisse(FEM_MAILLAGE* mai, MG_GESTIONNAIRE& gest2, double epsilon, double sigma, double sigmaf, double sigmag, double gamma, int iter_max, int reactiv, int bruit, int liss,int opti,int imp,int m_auto)
981 picher 233 {
982 francois 234 //Menu de la m�thode de lissage
983 francois 272 affiche((char*)"Extraction du maillage de surface");
984 picher 231 int coderesu = 0;
985 picher 233 double frac_min = 0.005;
986     int mai2_id;
987 francois 232 if (opti==1) conserve(OPTIMISE);
988     if (imp==1) conserve(IMPOSE);
989 picher 233 if (m_auto==1) conserve(MAILLEUR_AUTO);
990     copieorigine(mai,gest2);
991     MG_MAILLAGE* mg_mai = (MG_MAILLAGE*)mai->get_mg_maillage();
992     gain_poids(mg_mai,gest2);
993     if (reactiv == 1)
994     {
995     reactivation(mg_mai,gest2);
996     }
997 picher 231 do
998     {
999 picher 233 coderesu = extract_skin(mg_mai,gest2,frac_min,&mai2_id);
1000 francois 234 gain_poids(mg_mai,gest2);
1001 picher 231 }
1002     while (coderesu == 0);
1003 picher 233
1004     MG_MAILLAGE* mg_mai2=gest2.get_mg_maillageid(mai2_id);
1005 francois 232 if (bruit == 1)
1006 picher 231 {
1007 francois 272 affiche((char*)"Bruitage du maillage initial");
1008 picher 233 bruitage(mg_mai2,gest2);
1009 picher 231 }
1010 francois 309 //Choix de la m�thode de lissage (1 = chen(2005), 2 = Jones(2007), 3 = chen(2008))
1011 francois 232 if (liss == 1)
1012 picher 231 {
1013 cuillier 438 affiche((char*)"Procedure de lissage 1");
1014 picher 233 lissage(mg_mai2,gest2,epsilon,sigma,iter_max);
1015 picher 231 }
1016 picher 248 if (liss == 2)
1017     {
1018 cuillier 438 affiche((char*)"Procedure de lissage 2");
1019 picher 248 lissage2(mg_mai2,gest2,sigmaf,sigmag,iter_max);
1020     }
1021     if (liss == 3)
1022     {
1023 cuillier 438 affiche((char*)"Procedure de lissage 3");
1024 picher 248 lissage3(mg_mai2,gest2,sigma,gamma,epsilon,iter_max);
1025     }
1026 cuillier 438 if (liss == 4)
1027     {
1028     affiche((char*)"Procedure de lissage 4");
1029     lissage4(mg_mai2,gest2,iter_max);
1030     }
1031     if (liss == 5)
1032     {
1033     affiche((char*)"Procedure de lissage 5");
1034     lissage5(mg_mai2,gest2,iter_max);
1035     }
1036     if (liss == 6)
1037     {
1038     affiche((char*)"Procedure de lissage 6");
1039     lissage6(mg_mai2,gest2,sigma,gamma,epsilon,iter_max);
1040     }
1041     if (liss == 7)
1042     {
1043     affiche((char*)"Procedure de lissage 7");
1044     lissage7(mg_mai2,gest2,iter_max);
1045     }
1046 picher 231 }
1047    
1048 picher 230 void MG_LISSAGE::bruitage(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2)
1049     {
1050 francois 234 //M�thode test pour bruiter un maillage initial et tester la m�thode de lissage
1051 picher 230 srand(time(NULL));
1052     LISTE_MG_NOEUD::iterator it_no;
1053     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1054     {
1055 francois 234 //On modifie la position des noeuds l�g�rement (ajout d'un bruit au maillage initial)
1056 picher 230 int signe;
1057     int delta;
1058     double delta_x,delta_y,delta_z;
1059     signe = rand()%2 + 1;
1060     delta = rand()%2 + 1;
1061     if(signe == 1)
1062     {
1063 picher 248 delta_x = delta*0.5;
1064 picher 230 }
1065     else //(signe == 2)
1066     {
1067 picher 248 delta_x = -delta*0.5;
1068 picher 230 }
1069     signe = rand()%2 + 1;
1070     delta = rand()%2 + 1;
1071     if(signe == 1)
1072     {
1073 picher 248 delta_y = delta*0.5;
1074 picher 230 }
1075     else //(signe == 2)
1076     {
1077 picher 248 delta_y = -delta*0.5;
1078 picher 230 }
1079     signe = rand()%2 + 1;
1080     delta = rand()%2 + 1;
1081     if(signe == 1)
1082     {
1083 picher 248 delta_z = delta*0.5;
1084 picher 230 }
1085     else //(signe == 2)
1086     {
1087 picher 248 delta_z = -delta*0.5;
1088 picher 230 }
1089     double x_new = noeud_i->get_x() + delta_x;
1090     double y_new = noeud_i->get_y() + delta_y;
1091     double z_new = noeud_i->get_z() + delta_z;
1092     noeud_i->change_x(x_new);
1093     noeud_i->change_y(y_new);
1094     noeud_i->change_z(z_new);
1095     }
1096     }
1097 francois 224
1098 picher 230 void MG_LISSAGE::lissage(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double epsilon, double sigma, int iter_max)
1099     {
1100 francois 234 //M�thode de lissage inspir�e de Chen(2005)
1101 picher 230 double un_sur_pi = 1./M_PI;
1102     int compteur = 0;
1103     int fin = 0;
1104 picher 231 FILE *out = fopen("convergence.txt","wt");
1105 francois 224
1106 picher 230 do
1107     {
1108     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
1109 picher 233 TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales2;
1110 picher 230 TPL_LISTE_ENTITE<double> liste_wij;
1111     LISTE_MG_TRIANGLE::iterator it_tri;
1112     int k = 0; //pour identifier les triangles pour liste_normales et liste_wij
1113     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
1114     {
1115     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
1116     OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
1117 picher 233 normal_f_i.norme();
1118     liste_normales2.ajouter(normal_f_i);
1119 picher 230 //Remplissage de la liste des voisins du triangle i
1120     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
1121     MG_NOEUD* noeud1 = mgtri_i->get_noeud1();
1122     double nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
1123     for (int j = 0;j<nb_voisins1;j++)
1124     {
1125     MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
1126     liste_voisins.ajouter(mgtri_1);
1127     }
1128     MG_NOEUD* noeud2 = mgtri_i->get_noeud2();
1129     double nb_voisins2 = noeud2->get_lien_triangle()->get_nb();
1130     for (int j = 0;j<nb_voisins2;j++)
1131     {
1132     MG_TRIANGLE_PEAU* mgtri_2 = (MG_TRIANGLE_PEAU*) noeud2->get_lien_triangle()->get(j);
1133     liste_voisins.ajouter(mgtri_2);
1134     }
1135     MG_NOEUD* noeud3 = mgtri_i->get_noeud3();
1136     double nb_voisins3 = noeud3->get_lien_triangle()->get_nb();
1137     for (int j = 0;j<nb_voisins3;j++)
1138     {
1139     MG_TRIANGLE_PEAU* mgtri_3 = (MG_TRIANGLE_PEAU*) noeud3->get_lien_triangle()->get(j);
1140     liste_voisins.ajouter(mgtri_3);
1141     }
1142     liste_voisins.supprimer(mgtri_i);
1143     int nb_voisins = liste_voisins.get_nb();
1144     double w_ij = 1./nb_voisins;
1145     double phi_i_min = 10.;
1146     double s_i = 0.0;
1147     double phi_im = 0.0;
1148 picher 231 double *phi_ij = new double[nb_voisins];
1149 picher 230 OT_VECTEUR_3D normal_f_i_mean(0.,0.,0.);
1150     OT_VECTEUR_3D eta_i(0.,0.,0.);
1151     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
1152     int j = 0;
1153     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
1154     {
1155     OT_VECTEUR_3D normal_f_j = mgtri_j->calcul_normal();
1156     //1-Calculer la normale moyenne pour chaque triangle
1157     normal_f_i_mean = normal_f_i_mean + w_ij*normal_f_j;
1158 francois 234 //2.1-On calcul l'angle entre normal_f_i et normal_f_j pour j allant de 1 � Nb_voisins
1159 picher 230 double prod_scalaire = normal_f_i*normal_f_j;
1160 picher 233 if (prod_scalaire > 1.)
1161     {
1162     prod_scalaire = 1.;
1163     }
1164     if (prod_scalaire < -1.)
1165     {
1166     prod_scalaire = -1.;
1167     }
1168 picher 230 phi_ij[j] = acos(prod_scalaire)*un_sur_pi;
1169     //2.2-On trouve le plus petit des angles et la normale heta_i correspondante
1170     if (phi_ij[j] < phi_i_min)
1171     {
1172     phi_i_min = phi_ij[j];
1173     eta_i = normal_f_j;
1174     }
1175     //3.1-On calcul l'angle moyen phi_im
1176     phi_im = phi_im + w_ij*phi_ij[j];
1177     j++;
1178     }
1179     normal_f_i_mean.norme();
1180     j = 0;
1181     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
1182     {
1183     //3.2-Calcul de s_i selon la variance
1184     s_i = s_i + w_ij*pow((phi_ij[j] - phi_im),2);
1185     j++;
1186     }
1187 picher 231 delete[] phi_ij;
1188 picher 230 //4-On calcule une nouvelle normale pour chaque triangle
1189     OT_VECTEUR_3D normal_f_i_new = ponderation_gaussian(s_i,sigma)*normal_f_i_mean + (1. - ponderation_gaussian(s_i,sigma))*eta_i;
1190     normal_f_i_new.norme();
1191     liste_normales.ajouter(normal_f_i_new);
1192     liste_wij.ajouter(w_ij);
1193     mgtri->change_nouveau_numero(k);
1194     k++;
1195     }
1196    
1197     LISTE_MG_NOEUD::iterator it_no;
1198     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1199     {
1200     int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
1201     double w_ij_prime = 0.0;
1202     OT_VECTEUR_3D v_temp(0.,0.,0.);
1203     OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
1204     for(int j=0;j<nb_voisins_j;j++)
1205     {
1206     MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
1207 francois 234 //On calcule le centro�de cj du triangle mgtri_j
1208 picher 230 MG_NOEUD* n1 = mgtri_j->get_noeud1();
1209     MG_NOEUD* n2 = mgtri_j->get_noeud2();
1210     MG_NOEUD* n3 = mgtri_j->get_noeud3();
1211     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
1212     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
1213     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
1214     //On forme le vecteur vi_cj
1215     OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
1216     OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_j->get_nouveau_numero());
1217 francois 234 // w_ij_prime correspond � la somme des aires des triangles voisins du noeuds
1218 picher 233 OT_VECTEUR_3D AB(n2->get_x() - n1->get_x(),n2->get_y() - n1->get_y(),n2->get_z() - n1->get_z());
1219     OT_VECTEUR_3D AC(n3->get_x() - n1->get_x(),n3->get_y() - n1->get_y(),n3->get_z() - n1->get_z());
1220     OT_VECTEUR_3D prodvect = AB&AC;
1221     double w_ij = 0.5*prodvect.get_longueur();
1222 picher 230 w_ij_prime = w_ij_prime + w_ij;
1223     v_temp = v_temp + w_ij*(vi_cj*normal_f_i_new)*normal_f_i_new;
1224     }
1225 francois 234 //5-On met � jour la position des noeuds
1226 picher 230 v_i = v_i + v_temp/w_ij_prime;
1227 picher 233 int origine = noeud_i->get_origine();
1228     if (origine != IMPOSE)
1229     {
1230     noeud_i->change_x(v_i.get_x());
1231     noeud_i->change_y(v_i.get_y());
1232     noeud_i->change_z(v_i.get_z());
1233     }
1234 picher 230 }
1235    
1236 francois 234 //Crit�re d'arr�t de l'algorithme
1237 picher 230 int l=0;
1238     int nb_tri = mg_mai->get_nb_mg_triangle();
1239     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
1240     {
1241     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
1242 picher 233 OT_VECTEUR_3D normal_f_i = liste_normales2.get(mgtri_i->get_nouveau_numero());
1243 picher 230 OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_i->get_nouveau_numero());
1244 picher 248 double critere = 1. - normal_f_i*normal_f_i_new;
1245 picher 231 fprintf(out,"%lf\n",critere);
1246 picher 230 if (critere <= epsilon) l++;
1247 picher 233 if (critere >= 1.)
1248     {
1249     int no_id = mgtri->get_noeud1()->get_id();
1250     char str[50];
1251     sprintf (str, "test%d.magic", no_id);
1252     }
1253 picher 230 }
1254 picher 248 double tolerance = 0.01*nb_tri;
1255     if (nb_tri - l <= 0) fin = 1;
1256 picher 230 compteur++;
1257 picher 248 cout << "Fin de l'iteration #" << compteur << endl;
1258     fprintf(out,"\nFin de l'iteration #%d\n\n",compteur);
1259 picher 230 }
1260     while ((fin == 0) && (compteur < iter_max));
1261    
1262     if (fin == 1)
1263     {
1264 picher 248 cout << "Convergence de la methode en " << compteur << " iterations" << endl;
1265 picher 230 }
1266     else
1267     {
1268 picher 248 cout << "Arret de la procedure apres " << compteur << " iterations" << endl;
1269 picher 230 }
1270    
1271 picher 231 fclose(out);
1272 picher 230 } //Fin de lissage
1273    
1274 picher 248 void MG_LISSAGE::lissage2(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double sigmaf, double sigmag, int iter_max)
1275     {
1276     //Non-Iterative Smoothing (NIS) (Jones(2007))
1277     LISTE_MG_NOEUD::iterator it_no;
1278     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_noeuds;
1279     int l = 0;
1280     //1ere passe de NIS
1281     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1282     {
1283     double somme_num_x = 0.;
1284     double somme_num_y = 0.;
1285     double somme_num_z = 0.;
1286     double somme_denom = 0.;
1287     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
1288     int nb_voisins1 = noeud_i->get_lien_triangle()->get_nb();
1289     for (int j = 0;j<nb_voisins1;j++)
1290     {
1291     MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
1292     liste_voisins.ajouter(mgtri_1);
1293     MG_NOEUD* no1 = mgtri_1->get_noeud1();
1294     int nb_vois1 = no1->get_lien_triangle()->get_nb();
1295     for (int k = 0;k<nb_vois1;k++)
1296     {
1297     MG_TRIANGLE_PEAU* mgtri1 = (MG_TRIANGLE_PEAU*) no1->get_lien_triangle()->get(k);
1298     liste_voisins.ajouter(mgtri1);
1299     }
1300     MG_NOEUD* no2 = mgtri_1->get_noeud2();
1301     int nb_vois2 = no2->get_lien_triangle()->get_nb();
1302     for (int k = 0;k<nb_vois2;k++)
1303     {
1304     MG_TRIANGLE_PEAU* mgtri2 = (MG_TRIANGLE_PEAU*) no2->get_lien_triangle()->get(k);
1305     liste_voisins.ajouter(mgtri2);
1306     }
1307     MG_NOEUD* no3 = mgtri_1->get_noeud3();
1308     int nb_vois3 = no3->get_lien_triangle()->get_nb();
1309     for (int k = 0;k<nb_vois3;k++)
1310     {
1311     MG_TRIANGLE_PEAU* mgtri3 = (MG_TRIANGLE_PEAU*) no3->get_lien_triangle()->get(k);
1312     liste_voisins.ajouter(mgtri3);
1313     }
1314     }
1315     int nb_voisins = liste_voisins.get_nb();
1316     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
1317     for(MG_TRIANGLE_PEAU* mgtri_peau=liste_voisins.get_premier(it);mgtri_peau!=NULL;mgtri_peau=liste_voisins.get_suivant(it))
1318     {
1319     OT_VECTEUR_3D normale_j = mgtri_peau->calcul_normal();
1320     //Centroide du triangle
1321     MG_NOEUD* n1 = mgtri_peau->get_noeud1();
1322     MG_NOEUD* n2 = mgtri_peau->get_noeud2();
1323     MG_NOEUD* n3 = mgtri_peau->get_noeud3();
1324     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
1325     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
1326     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
1327     //spatial weight
1328     double px = noeud_i->get_x();
1329     double py = noeud_i->get_y();
1330     double pz = noeud_i->get_z();
1331     double spatial = sqrt(pow(cj_x - px,2) + pow(cj_y - py,2) +pow(cj_z - pz,2));
1332     double f = ponderation_gaussian(spatial,sigmaf/2.);
1333     //calcul de l'aire du triangle
1334     OT_VECTEUR_3D AB(n2->get_x() - n1->get_x(),n2->get_y() - n1->get_y(),n2->get_z() - n1->get_z());
1335     OT_VECTEUR_3D AC(n3->get_x() - n1->get_x(),n3->get_y() - n1->get_y(),n3->get_z() - n1->get_z());
1336     OT_VECTEUR_3D prodvect = AB&AC;
1337     double aire_q = 0.5*prodvect.get_longueur();
1338     somme_num_x = somme_num_x + cj_x*aire_q*f;
1339     somme_num_y = somme_num_y + cj_y*aire_q*f;
1340     somme_num_z = somme_num_z + cj_z*aire_q*f;
1341     somme_denom = somme_denom + aire_q*f;
1342     }
1343 francois 309 //On ne change pas r�ellement la position des noeuds lors de la premi�re passe, mais on conserve les normales mollifi�es !!
1344 picher 248 OT_VECTEUR_3D p_prime(somme_num_x/somme_denom,somme_num_y/somme_denom,somme_num_z/somme_denom);
1345     liste_noeuds.ajouter(p_prime);
1346     noeud_i->change_nouveau_numero(l);
1347     l++;
1348     }
1349     //Calcul des nouvelles normales
1350     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
1351     LISTE_MG_TRIANGLE::iterator it_tri;
1352     l = 0;
1353     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
1354     {
1355     MG_NOEUD* n1 = mgtri->get_noeud1();
1356     OT_VECTEUR_3D p1 = liste_noeuds.get(n1->get_nouveau_numero());
1357     MG_NOEUD* n2 = mgtri->get_noeud2();
1358     OT_VECTEUR_3D p2 = liste_noeuds.get(n2->get_nouveau_numero());
1359     MG_NOEUD* n3 = mgtri->get_noeud3();
1360     OT_VECTEUR_3D p3 = liste_noeuds.get(n3->get_nouveau_numero());
1361     OT_VECTEUR_3D AB(p2.get_x() - p1.get_x(),p2.get_y() - p1.get_y(),p2.get_z() - p1.get_z());
1362     OT_VECTEUR_3D AC(p3.get_x() - p1.get_x(),p3.get_y() - p1.get_y(),p3.get_z() - p1.get_z());
1363     OT_VECTEUR_3D normale = AB&AC;
1364     normale.norme();
1365     liste_normales.ajouter(normale);
1366     mgtri->change_nouveau_numero(l);
1367     l++;
1368     }
1369     //2e passe de NIS
1370     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_pprime;
1371     l = 0;
1372     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1373     {
1374     double somme_num_x = 0.;
1375     double somme_num_y = 0.;
1376     double somme_num_z = 0.;
1377     double somme_denom = 0.;
1378     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
1379     int nb_voisins1 = noeud_i->get_lien_triangle()->get_nb();
1380     for (int j = 0;j<nb_voisins1;j++)
1381     {
1382     MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
1383     liste_voisins.ajouter(mgtri_1);
1384     MG_NOEUD* no1 = mgtri_1->get_noeud1();
1385     int nb_vois1 = no1->get_lien_triangle()->get_nb();
1386     for (int k = 0;k<nb_vois1;k++)
1387     {
1388     MG_TRIANGLE_PEAU* mgtri1 = (MG_TRIANGLE_PEAU*) no1->get_lien_triangle()->get(k);
1389     liste_voisins.ajouter(mgtri1);
1390     }
1391     MG_NOEUD* no2 = mgtri_1->get_noeud2();
1392     int nb_vois2 = no2->get_lien_triangle()->get_nb();
1393     for (int k = 0;k<nb_vois2;k++)
1394     {
1395     MG_TRIANGLE_PEAU* mgtri2 = (MG_TRIANGLE_PEAU*) no2->get_lien_triangle()->get(k);
1396     liste_voisins.ajouter(mgtri2);
1397     }
1398     MG_NOEUD* no3 = mgtri_1->get_noeud3();
1399     int nb_vois3 = no3->get_lien_triangle()->get_nb();
1400     for (int k = 0;k<nb_vois3;k++)
1401     {
1402     MG_TRIANGLE_PEAU* mgtri3 = (MG_TRIANGLE_PEAU*) no3->get_lien_triangle()->get(k);
1403     liste_voisins.ajouter(mgtri3);
1404     }
1405     }
1406     int nb_voisins = liste_voisins.get_nb();
1407     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
1408     for(MG_TRIANGLE_PEAU* mgtri_peau=liste_voisins.get_premier(it);mgtri_peau!=NULL;mgtri_peau=liste_voisins.get_suivant(it))
1409     {
1410     OT_VECTEUR_3D normale_j = liste_normales.get(mgtri_peau->get_nouveau_numero());
1411     //Centroide du triangle
1412     MG_NOEUD* n1 = mgtri_peau->get_noeud1();
1413     MG_NOEUD* n2 = mgtri_peau->get_noeud2();
1414     MG_NOEUD* n3 = mgtri_peau->get_noeud3();
1415     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
1416     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
1417     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
1418     //spatial weight
1419     double px = noeud_i->get_x();
1420     double py = noeud_i->get_y();
1421     double pz = noeud_i->get_z();
1422     double spatial = sqrt(pow(cj_x - px,2) + pow(cj_y - py,2) +pow(cj_z - pz,2));
1423     double f = ponderation_gaussian(spatial,sigmaf);
1424     //calcul de l'aire du triangle
1425     OT_VECTEUR_3D AB(n2->get_x() - n1->get_x(),n2->get_y() - n1->get_y(),n2->get_z() - n1->get_z());
1426     OT_VECTEUR_3D AC(n3->get_x() - n1->get_x(),n3->get_y() - n1->get_y(),n3->get_z() - n1->get_z());
1427     OT_VECTEUR_3D prodvect = AB&AC;
1428     double aire_q = 0.5*prodvect.get_longueur();
1429     //Calcul de la projection du noeuds sur le plan tangent au triangle
1430     double a = normale_j.get_x();
1431     double b = normale_j.get_y();
1432     double c = normale_j.get_z();
1433     double d = -(a*n1->get_x() + b*n1->get_y() + c*n1->get_z());
1434     double k = -(a*px + b*py + c*pz + d)/(a*a + b*b + c*c);
1435     double proj_x = px + k*a;
1436     double proj_y = py + k*b;
1437     double proj_z = pz + k*c;
1438     //influence weight
1439     double influence = sqrt(pow(proj_x - px,2) + pow(proj_y - py,2) + pow(proj_z - pz,2));
1440     double g = ponderation_gaussian(influence,sigmag);
1441     somme_num_x = somme_num_x + proj_x*aire_q*f*g;
1442     somme_num_y = somme_num_y + proj_y*aire_q*f*g;
1443     somme_num_z = somme_num_z + proj_z*aire_q*f*g;
1444     somme_denom = somme_denom + aire_q*f*g;
1445     }
1446     OT_VECTEUR_3D pprime(somme_num_x/somme_denom,somme_num_y/somme_denom,somme_num_z/somme_denom);
1447     liste_pprime.ajouter(pprime);
1448     noeud_i->change_nouveau_numero(l);
1449     l++;
1450     }
1451 francois 309 //Mise � jour en bloc de la position des noeuds
1452 picher 248 for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1453     {
1454     int origine = noeud_i->get_origine();
1455     if (origine != IMPOSE)
1456     {
1457     OT_VECTEUR_3D pprime = liste_pprime.get(noeud_i->get_nouveau_numero());
1458     noeud_i->change_x(pprime.get_x());
1459     noeud_i->change_y(pprime.get_y());
1460     noeud_i->change_z(pprime.get_z());
1461     }
1462     }
1463    
1464     } //Fin de lissage2
1465    
1466     void MG_LISSAGE::lissage3(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double sigma, double gamma, double epsilon, int iter_max)
1467     {
1468     //M�thode de lissage inspir�e de Chen(2008)
1469     double un_sur_pi = 1./M_PI;
1470     int compteur = 0;
1471     int fin = 0;
1472     do
1473     {
1474     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
1475     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales2;
1476     TPL_LISTE_ENTITE<double> liste_wij;
1477     LISTE_MG_TRIANGLE::iterator it_tri;
1478     int k = 0; //pour identifier les triangles pour liste_normales et liste_wij
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 = mgtri_i->calcul_normal();
1483     normal_f_i.norme();
1484     liste_normales2.ajouter(normal_f_i);
1485     //Remplissage de la liste des voisins du triangle i
1486     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
1487     MG_NOEUD* noeud1 = mgtri_i->get_noeud1();
1488     double nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
1489     for (int j = 0;j<nb_voisins1;j++)
1490     {
1491     MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
1492     liste_voisins.ajouter(mgtri_1);
1493     }
1494     MG_NOEUD* noeud2 = mgtri_i->get_noeud2();
1495     double nb_voisins2 = noeud2->get_lien_triangle()->get_nb();
1496     for (int j = 0;j<nb_voisins2;j++)
1497     {
1498     MG_TRIANGLE_PEAU* mgtri_2 = (MG_TRIANGLE_PEAU*) noeud2->get_lien_triangle()->get(j);
1499     liste_voisins.ajouter(mgtri_2);
1500     }
1501     MG_NOEUD* noeud3 = mgtri_i->get_noeud3();
1502     double nb_voisins3 = noeud3->get_lien_triangle()->get_nb();
1503     for (int j = 0;j<nb_voisins3;j++)
1504     {
1505     MG_TRIANGLE_PEAU* mgtri_3 = (MG_TRIANGLE_PEAU*) noeud3->get_lien_triangle()->get(j);
1506     liste_voisins.ajouter(mgtri_3);
1507     }
1508     liste_voisins.supprimer(mgtri_i);
1509     int nb_voisins = liste_voisins.get_nb();
1510     double w_ij = 1./nb_voisins;
1511     double phi_i_min = 10.;
1512     double s_i = 0.0;
1513     double phi_im = 0.0;
1514     double *phi_ij = new double[nb_voisins];
1515     OT_VECTEUR_3D normal_f_i_mean(0.,0.,0.);
1516     normal_f_i_mean = normal_f_i;
1517     OT_VECTEUR_3D eta_i(0.,0.,0.);
1518     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
1519     int j = 0;
1520     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
1521     {
1522     OT_VECTEUR_3D normal_f_j = mgtri_j->calcul_normal();
1523     //1-Calculer la normale moyenne pour chaque triangle
1524     normal_f_i_mean = normal_f_i_mean + normal_f_j;
1525     //2.1-On calcul l'angle entre normal_f_i et normal_f_j pour j allant de 1 � Nb_voisins
1526     double prod_scalaire = normal_f_i*normal_f_j;
1527     if (prod_scalaire > 1.)
1528     {
1529     prod_scalaire = 1.;
1530     }
1531     if (prod_scalaire < -1.)
1532     {
1533     prod_scalaire = -1.;
1534     }
1535     phi_ij[j] = acos(prod_scalaire)*un_sur_pi;
1536     //2.2-On trouve le plus petit des angles et la normale heta_i correspondante
1537     if (phi_ij[j] < phi_i_min)
1538     {
1539     phi_i_min = phi_ij[j];
1540     eta_i = normal_f_j;
1541     }
1542     //3.1-On calcul l'angle moyen phi_im
1543     phi_im = phi_im + w_ij*phi_ij[j];
1544     j++;
1545     }
1546     normal_f_i_mean.norme();
1547     j = 0;
1548     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
1549     {
1550     //3.2-Calcul de s_i selon la variance
1551     s_i = s_i + w_ij*pow((phi_ij[j] - phi_im),2);
1552     j++;
1553     }
1554     delete[] phi_ij;
1555     //4-On calcule une nouvelle normale pour chaque triangle
1556     OT_VECTEUR_3D normal_f_i_new = ponderation_gaussian(s_i,sigma)*normal_f_i_mean + (1. - ponderation_gaussian(s_i,sigma))*eta_i;
1557     normal_f_i_new.norme();
1558     liste_normales.ajouter(normal_f_i_new);
1559     liste_wij.ajouter(w_ij);
1560     mgtri->change_nouveau_numero(k);
1561     k++;
1562     }
1563    
1564     LISTE_MG_NOEUD::iterator it_no;
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     {
1567     int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
1568     double w_ij_prime = 0.0;
1569     OT_VECTEUR_3D v_temp(0.,0.,0.);
1570     OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
1571     for(int j=0;j<nb_voisins_j;j++)
1572     {
1573     MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
1574     //On calcule le centro�de cj du triangle mgtri_j
1575     MG_NOEUD* n1 = mgtri_j->get_noeud1();
1576     MG_NOEUD* n2 = mgtri_j->get_noeud2();
1577     MG_NOEUD* n3 = mgtri_j->get_noeud3();
1578     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
1579     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
1580     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
1581     //On forme le vecteur vi_cj
1582     OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
1583     OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_j->get_nouveau_numero());
1584     v_temp = v_temp + (vi_cj*normal_f_i_new)*normal_f_i_new;
1585     }
1586     //5-On met � jour la position des noeuds
1587     v_i = v_i + (gamma/(2*nb_voisins_j))*v_temp;
1588     int origine = noeud_i->get_origine();
1589     if (origine != IMPOSE)
1590     {
1591     noeud_i->change_x(v_i.get_x());
1592     noeud_i->change_y(v_i.get_y());
1593     noeud_i->change_z(v_i.get_z());
1594     }
1595     }
1596    
1597     //Crit�re d'arr�t de l'algorithme
1598     int l=0;
1599     int nb_tri = mg_mai->get_nb_mg_triangle();
1600     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
1601     {
1602     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
1603     OT_VECTEUR_3D normal_f_i = liste_normales2.get(mgtri_i->get_nouveau_numero());
1604     OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_i->get_nouveau_numero());
1605     double critere = 1. - normal_f_i*normal_f_i_new;
1606 cuillier 438 if (critere <= epsilon) l++;
1607     if (critere >= 1.)
1608     {
1609     int no_id = mgtri->get_noeud1()->get_id();
1610     char str[50];
1611     sprintf (str, "test%d.magic", no_id);
1612     }
1613     }
1614     double tolerance = 0.01*nb_tri;
1615     if (nb_tri - l <= 0) fin = 1;
1616     compteur++;
1617     cout << "Fin de l'iteration #" << compteur << endl;
1618     }
1619     while ((fin == 0) && (compteur < iter_max));
1620    
1621     if (fin == 1)
1622     {
1623     cout << "Convergence de la methode en " << compteur << " iterations" << endl;
1624     }
1625     else
1626     {
1627     cout << "Arret de la procedure apres " << compteur << " iterations" << endl;
1628     }
1629     }
1630    
1631     void MG_LISSAGE::lissage4(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, int iter_max)
1632     {
1633     // Lissage de proche en proche avec methode de lissage inspiree de JIAO IMR2012
1634     double un_sur_pi = 1./M_PI;
1635     int compteur = 0;
1636    
1637     do
1638     {
1639     LISTE_MG_NOEUD::iterator it_no;
1640     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1641     {
1642     int origine = noeud_i->get_origine();
1643     if (origine != IMPOSE)
1644     {
1645     int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
1646     double wij = 0.0;
1647     double wij_sommej = 0.0;
1648     OT_VECTEUR_3D v_temp(0.,0.,0.);
1649     OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
1650     for(int j=0;j<nb_voisins_j;j++)
1651     {
1652     MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
1653     //On calcule le centro�de cj du triangle mgtri_j
1654     MG_NOEUD* n1 = mgtri_j->get_noeud1();
1655     MG_NOEUD* n2 = mgtri_j->get_noeud2();
1656     MG_NOEUD* n3 = mgtri_j->get_noeud3();
1657     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
1658     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
1659     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
1660     //On forme le vecteur vi_cj et le vecteur cj
1661     OT_VECTEUR_3D cj(cj_x ,cj_y,cj_z);
1662     OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
1663     wij=vi_cj.get_longueur();
1664     wij_sommej=wij_sommej+wij;
1665     v_temp = v_temp + (wij*cj);
1666     }
1667     //5-On met � jour la position des noeuds
1668     v_i = v_temp/wij_sommej;
1669     noeud_i->change_x(v_i.get_x());
1670     noeud_i->change_y(v_i.get_y());
1671     noeud_i->change_z(v_i.get_z());
1672     }
1673     }
1674     //Crit�re d'arr�t de l'algorithme
1675     compteur++;
1676     cout << "Fin de l'iteration #" << compteur << endl;
1677     }
1678     while (compteur < iter_max);
1679    
1680     cout << "Arret de la procedure lissage 4 apres " << compteur << " iterations" << endl;
1681     }
1682    
1683    
1684     void MG_LISSAGE::lissage5(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, int iter_max)
1685     {
1686     // Lissage global avec methode de lissage inspiree de JIAO IMR2012
1687     double un_sur_pi = 1./M_PI;
1688     int compteur = 0;
1689    
1690     do
1691     {
1692     vector<double> nouv_position_x;
1693     vector<double> nouv_position_y;
1694     vector<double> nouv_position_z;
1695     LISTE_MG_NOEUD::iterator it_no;
1696     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1697     {
1698     int origine = noeud_i->get_origine();
1699     if (origine != IMPOSE)
1700     {
1701     int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
1702     double wij = 0.0;
1703     double wij_sommej = 0.0;
1704     OT_VECTEUR_3D v_temp(0.,0.,0.);
1705     OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
1706     for(int j=0;j<nb_voisins_j;j++)
1707     {
1708     MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
1709     //On calcule le centro�de cj du triangle mgtri_j
1710     MG_NOEUD* n1 = mgtri_j->get_noeud1();
1711     MG_NOEUD* n2 = mgtri_j->get_noeud2();
1712     MG_NOEUD* n3 = mgtri_j->get_noeud3();
1713     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
1714     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
1715     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
1716     //On forme le vecteur vi_cj et le vecteur cj
1717     OT_VECTEUR_3D cj(cj_x ,cj_y,cj_z);
1718     OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
1719     wij=vi_cj.get_longueur();
1720     wij_sommej=wij_sommej+wij;
1721     v_temp = v_temp + (wij*cj);
1722     }
1723     //5-On calcule la nouvelle position des noeuds et on affecte au tableau temporaire
1724     v_i =v_temp/wij_sommej;
1725     nouv_position_x.push_back(v_i.get_x());
1726     nouv_position_y.push_back(v_i.get_y());
1727     nouv_position_z.push_back(v_i.get_z());
1728     }
1729     }
1730     //On actualise la position des noeuds
1731     int ind_noeud=0;
1732     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1733     {
1734     int origine = noeud_i->get_origine();
1735     if (origine != IMPOSE)
1736     {
1737     noeud_i->change_x(nouv_position_x[ind_noeud]);
1738     noeud_i->change_y(nouv_position_y[ind_noeud]);
1739     noeud_i->change_z(nouv_position_z[ind_noeud]);
1740     ind_noeud++;
1741     }
1742     }
1743     //Critere d'arret de l'algorithme
1744     compteur++;
1745     cout << "Fin de l'iteration #" << compteur << endl;
1746     }
1747     while (compteur < iter_max);
1748    
1749     cout << "Arret de la procedure lissage 5 apres " << compteur << " iterations" << endl;
1750     }
1751    
1752    
1753     void MG_LISSAGE::lissage6(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double sigma, double gamma, double epsilon, int iter_max)
1754     {
1755     //Quelques iterations type JIAO IMR2012 suivies de lissage inspiree de Chen(2008) corrigee par jcc sept 2013
1756     double un_sur_pi = 1./M_PI;
1757     int compteur = 0;
1758     int fin = 0;
1759     FILE *out = fopen("convergence.txt","wt");
1760    
1761     do
1762     {
1763     vector<double> nouv_position_x;
1764     vector<double> nouv_position_y;
1765     vector<double> nouv_position_z;
1766     LISTE_MG_NOEUD::iterator it_no;
1767     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1768     {
1769     int origine = noeud_i->get_origine();
1770     if (origine != IMPOSE)
1771     {
1772     int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
1773     double wij = 0.0;
1774     double wij_sommej = 0.0;
1775     OT_VECTEUR_3D v_temp(0.,0.,0.);
1776     OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
1777     for(int j=0;j<nb_voisins_j;j++)
1778     {
1779     MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
1780     //On calcule le centro�de cj du triangle mgtri_j
1781     MG_NOEUD* n1 = mgtri_j->get_noeud1();
1782     MG_NOEUD* n2 = mgtri_j->get_noeud2();
1783     MG_NOEUD* n3 = mgtri_j->get_noeud3();
1784     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
1785     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
1786     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
1787     //On forme le vecteur vi_cj et le vecteur cj
1788     OT_VECTEUR_3D cj(cj_x ,cj_y,cj_z);
1789     OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
1790     wij=vi_cj.get_longueur();
1791     wij_sommej=wij_sommej+wij;
1792     v_temp = v_temp + (wij*cj);
1793     }
1794     //5-On calcule la nouvelle position des noeuds et on affecte au tableau temporaire
1795     v_i =v_temp/wij_sommej;
1796     nouv_position_x.push_back(v_i.get_x());
1797     nouv_position_y.push_back(v_i.get_y());
1798     nouv_position_z.push_back(v_i.get_z());
1799     }
1800     }
1801     //On actualise la position des noeuds
1802     int ind_noeud=0;
1803     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1804     {
1805     int origine = noeud_i->get_origine();
1806     if (origine != IMPOSE)
1807     {
1808     noeud_i->change_x(nouv_position_x[ind_noeud]);
1809     noeud_i->change_y(nouv_position_y[ind_noeud]);
1810     noeud_i->change_z(nouv_position_z[ind_noeud]);
1811     ind_noeud++;
1812     }
1813     }
1814     //Critere d'arret de l'algorithme
1815     compteur++;
1816     cout << "Fin de l'iteration #" << compteur << endl;
1817     }
1818     while (compteur < 10);
1819    
1820     cout << "Arret de la procedure lissage laplacien " << endl;
1821     compteur = 0;
1822     do
1823     {
1824     vector<double> nouv_position_x;
1825     vector<double> nouv_position_y;
1826     vector<double> nouv_position_z;
1827     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
1828     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales2;
1829     TPL_LISTE_ENTITE<double> liste_wij;
1830     LISTE_MG_TRIANGLE::iterator it_tri;
1831     int k = 0; //pour identifier les triangles pour liste_normales et liste_wij
1832     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
1833     {
1834     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
1835     OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
1836     normal_f_i.norme();
1837     liste_normales2.ajouter(normal_f_i);
1838     //Remplissage de la liste des voisins du triangle i
1839     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
1840     MG_NOEUD* noeud1 = mgtri_i->get_noeud1();
1841     double nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
1842     for (int j = 0;j<nb_voisins1;j++)
1843     {
1844     MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
1845     liste_voisins.ajouter(mgtri_1);
1846     }
1847     MG_NOEUD* noeud2 = mgtri_i->get_noeud2();
1848     double nb_voisins2 = noeud2->get_lien_triangle()->get_nb();
1849     for (int j = 0;j<nb_voisins2;j++)
1850     {
1851     MG_TRIANGLE_PEAU* mgtri_2 = (MG_TRIANGLE_PEAU*) noeud2->get_lien_triangle()->get(j);
1852     liste_voisins.ajouter(mgtri_2);
1853     }
1854     MG_NOEUD* noeud3 = mgtri_i->get_noeud3();
1855     double nb_voisins3 = noeud3->get_lien_triangle()->get_nb();
1856     for (int j = 0;j<nb_voisins3;j++)
1857     {
1858     MG_TRIANGLE_PEAU* mgtri_3 = (MG_TRIANGLE_PEAU*) noeud3->get_lien_triangle()->get(j);
1859     liste_voisins.ajouter(mgtri_3);
1860     }
1861     liste_voisins.supprimer(mgtri_i);
1862     int nb_voisins = liste_voisins.get_nb();
1863     double w_ij = 1./nb_voisins;
1864     double phi_i_min = 10.;
1865     double s_i = 0.0;
1866     double phi_im = 0.0;
1867     double *phi_ij = new double[nb_voisins];
1868     OT_VECTEUR_3D normal_f_i_mean(0.,0.,0.);
1869     normal_f_i_mean = normal_f_i;
1870     OT_VECTEUR_3D eta_i(0.,0.,0.);
1871     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
1872     int j = 0;
1873     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
1874     {
1875     OT_VECTEUR_3D normal_f_j = mgtri_j->calcul_normal();
1876     //1-Calculer la normale moyenne pour chaque triangle
1877     normal_f_i_mean = normal_f_i_mean + normal_f_j;
1878     //2.1-On calcul l'angle entre normal_f_i et normal_f_j pour j allant de 1 � Nb_voisins
1879     double prod_scalaire = normal_f_i*normal_f_j;
1880     if (prod_scalaire > 1.)
1881     {
1882     prod_scalaire = 1.;
1883     }
1884     if (prod_scalaire < -1.)
1885     {
1886     prod_scalaire = -1.;
1887     }
1888     phi_ij[j] = acos(prod_scalaire)*un_sur_pi;
1889     //2.2-On trouve le plus petit des angles et la normale heta_i correspondante
1890     if (phi_ij[j] < phi_i_min)
1891     {
1892     phi_i_min = phi_ij[j];
1893     eta_i = normal_f_j;
1894     }
1895     //3.1-On calcul l'angle moyen phi_im
1896     phi_im = phi_im + w_ij*phi_ij[j];
1897     j++;
1898     }
1899     normal_f_i_mean.norme();
1900     j = 0;
1901     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
1902     {
1903     //3.2-Calcul de s_i selon la variance
1904     s_i = s_i + w_ij*pow((phi_ij[j] - phi_im),2);
1905     j++;
1906     }
1907     delete[] phi_ij;
1908     //4-On calcule une nouvelle normale pour chaque triangle
1909     OT_VECTEUR_3D normal_f_i_new = ponderation_gaussian(s_i,sigma)*normal_f_i_mean + (1. - ponderation_gaussian(s_i,sigma))*eta_i;
1910     normal_f_i_new.norme();
1911     liste_normales.ajouter(normal_f_i_new);
1912     liste_wij.ajouter(w_ij);
1913     mgtri->change_nouveau_numero(k);
1914     k++;
1915     }
1916    
1917     LISTE_MG_NOEUD::iterator it_no;
1918     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1919     {
1920     int origine = noeud_i->get_origine();
1921     if (origine != IMPOSE)
1922     {
1923     int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
1924     double w_ij_prime = 0.0;
1925     OT_VECTEUR_3D v_temp(0.,0.,0.);
1926     OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
1927     for(int j=0;j<nb_voisins_j;j++)
1928     {
1929     MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
1930     //On calcule le centro�de cj du triangle mgtri_j
1931     MG_NOEUD* n1 = mgtri_j->get_noeud1();
1932     MG_NOEUD* n2 = mgtri_j->get_noeud2();
1933     MG_NOEUD* n3 = mgtri_j->get_noeud3();
1934     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
1935     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
1936     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
1937     //On forme le vecteur vi_cj
1938     OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
1939     OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_j->get_nouveau_numero());
1940     v_temp = v_temp + (vi_cj*normal_f_i_new)*normal_f_i_new;
1941     }
1942     //5-On met � jour la position des noeuds
1943     v_i = v_i + (gamma/(2*nb_voisins_j))*v_temp;
1944     nouv_position_x.push_back(v_i.get_x());
1945     nouv_position_y.push_back(v_i.get_y());
1946     nouv_position_z.push_back(v_i.get_z());
1947     }
1948     }
1949    
1950     //On actualise la position des noeuds
1951     int ind_noeud=0;
1952     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
1953     {
1954     int origine = noeud_i->get_origine();
1955     if (origine != IMPOSE)
1956     {
1957     noeud_i->change_x(nouv_position_x[ind_noeud]);
1958     noeud_i->change_y(nouv_position_y[ind_noeud]);
1959     noeud_i->change_z(nouv_position_z[ind_noeud]);
1960     ind_noeud++;
1961     }
1962     }
1963    
1964     //Critere d'arret de l'algorithme
1965     int l=0;
1966     int nb_tri = mg_mai->get_nb_mg_triangle();
1967     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
1968     {
1969     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
1970     OT_VECTEUR_3D normal_f_i = liste_normales2.get(mgtri_i->get_nouveau_numero());
1971     OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_i->get_nouveau_numero());
1972     double critere = 1. - normal_f_i*normal_f_i_new;
1973 picher 248 fprintf(out,"%lf\n",critere);
1974     if (critere <= epsilon) l++;
1975     if (critere >= 1.)
1976     {
1977     int no_id = mgtri->get_noeud1()->get_id();
1978     char str[50];
1979     sprintf (str, "test%d.magic", no_id);
1980     }
1981     }
1982     double tolerance = 0.01*nb_tri;
1983     if (nb_tri - l <= 0) fin = 1;
1984     compteur++;
1985     cout << "Fin de l'iteration #" << compteur << endl;
1986     fprintf(out,"\nFin de l'iteration #%d\n\n",compteur);
1987     }
1988     while ((fin == 0) && (compteur < iter_max));
1989    
1990     if (fin == 1)
1991     {
1992     cout << "Convergence de la methode en " << compteur << " iterations" << endl;
1993     }
1994     else
1995     {
1996     cout << "Arret de la procedure apres " << compteur << " iterations" << endl;
1997     }
1998    
1999     fclose(out);
2000     }
2001    
2002 cuillier 438
2003    
2004    
2005     void MG_LISSAGE::lissage7(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, int iter_max)
2006     {
2007     //Methode de lissage inspire de Semenova 2004
2008     double un_sur_pi = 1./M_PI;
2009     int compteur = 0;
2010     do
2011     {
2012     LISTE_MG_NOEUD::iterator it_no;
2013     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
2014     {
2015     int origine = noeud_i->get_origine();
2016     if (origine != IMPOSE)
2017     {
2018     TPL_MAP_ENTITE<MG_SEGMENT*> lst_seg;
2019     int nb_seg_voisin_j = noeud_i->get_lien_segment()->get_nb();
2020     for(int j=0;j<nb_seg_voisin_j;j++)
2021     {
2022     MG_SEGMENT* mgseg_j = (MG_SEGMENT*) noeud_i->get_lien_segment()->get(j);
2023     lst_seg.ajouter(mgseg_j);
2024     }
2025     LISTE_MG_SEGMENT::iterator it_seg;
2026     MG_NOEUD* noeud_j;
2027     MG_NOEUD* noeud_k;
2028    
2029     OT_VECTEUR_3D p_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
2030     double p_in_x=0.0;double p_in_y=0.0;double p_in_z=0.0;
2031     do
2032     {
2033     MG_SEGMENT* segment_j=lst_seg.get_premier(it_seg);
2034     if (segment_j->get_noeud1()==noeud_i) {noeud_j=segment_j->get_noeud2();} else {noeud_j=segment_j->get_noeud1();}
2035     OT_VECTEUR_3D p_j(noeud_j->get_x(),noeud_j->get_y(),noeud_j->get_z());
2036     double angle_max=-10.0;
2037     MG_SEGMENT* seg_candidat;
2038     MG_NOEUD* noeud_c;
2039     for(int k=0;k<nb_seg_voisin_j;k++)
2040     {
2041     MG_SEGMENT* segment_k = (MG_SEGMENT*) noeud_i->get_lien_segment()->get(k);
2042     if (segment_k->get_noeud1()==noeud_i) {noeud_k=segment_k->get_noeud2();} else {noeud_k=segment_k->get_noeud1();}
2043     OT_VECTEUR_3D p_k(noeud_k->get_x(),noeud_k->get_y(),noeud_k->get_z());
2044     OT_VECTEUR_3D p_i_p_j(p_j.get_x()-p_i.get_x(),p_j.get_y()-p_i.get_y(),p_j.get_z()-p_i.get_z());
2045     OT_VECTEUR_3D p_i_p_k(p_k.get_x()-p_i.get_x(),p_k.get_y()-p_i.get_y(),p_k.get_z()-p_i.get_z());
2046     double l_pipj=p_i_p_j.get_longueur();
2047     double l_pipk=p_i_p_k.get_longueur();
2048     double norm=(p_i_p_j.get_x()*p_i_p_k.get_x()+p_i_p_j.get_y()*p_i_p_k.get_y()+p_i_p_j.get_z()*p_i_p_k.get_z())/(l_pipj*l_pipk);
2049     if (norm > 1.0)
2050     { norm = 1.0;}
2051     if (norm < -1.0)
2052     { norm = -1.0;}
2053     double angle = acos (norm);
2054     //cout << "angle" << angle << endl;
2055     if (angle >= angle_max)
2056     {
2057     angle_max=angle;
2058     seg_candidat=segment_k;
2059     }
2060     }
2061     //cout << "angle max" << angle_max << endl;
2062     if (seg_candidat->get_noeud1()==noeud_i) {noeud_c=seg_candidat->get_noeud2(); } else {noeud_c=seg_candidat->get_noeud1(); }
2063     OT_VECTEUR_3D p_c(noeud_c->get_x(),noeud_c->get_y(),noeud_c->get_z());
2064     OT_VECTEUR_3D p_m(0.5*(p_j.get_x()+p_c.get_x()),0.5*(p_j.get_y()+p_c.get_y()),0.5*(p_j.get_z()+p_c.get_z()));
2065     OT_VECTEUR_3D p_m_p_i(p_i.get_x()-p_m.get_x(),p_i.get_y()-p_m.get_y(),p_i.get_z()-p_m.get_z());
2066     OT_VECTEUR_3D p_m_p_c(p_c.get_x()-p_m.get_x(),p_c.get_y()-p_m.get_y(),p_c.get_z()-p_m.get_z());
2067     OT_VECTEUR_3D p_c_p_j(p_j.get_x()-p_c.get_x(),p_j.get_y()-p_c.get_y(),p_j.get_z()-p_c.get_z());
2068     double l_pmpc=p_m_p_c.get_longueur();
2069     double l_pcpj=p_c_p_j.get_longueur();
2070     double prod=p_m_p_i.get_x()*p_m_p_c.get_x()+p_m_p_i.get_y()*p_m_p_c.get_y()+p_m_p_i.get_z()*p_m_p_c.get_z();
2071     p_in_x=p_in_x+p_i.get_x()+(prod/(l_pcpj*l_pmpc))*p_c_p_j.get_x();
2072     p_in_y=p_in_y+p_i.get_y()+(prod/(l_pcpj*l_pmpc))*p_c_p_j.get_y();
2073     p_in_z=p_in_z+p_i.get_z()+(prod/(l_pcpj*l_pmpc))*p_c_p_j.get_z();
2074     lst_seg.supprimer(segment_j);
2075     }
2076     while(lst_seg.get_nb()>0);
2077     //5-On met � jour la position du noeud i
2078     p_in_x=p_in_x/(double)nb_seg_voisin_j;
2079     p_in_y=p_in_y/(double)nb_seg_voisin_j;
2080     p_in_z=p_in_z/(double)nb_seg_voisin_j;
2081     noeud_i->change_x(p_in_x);
2082     noeud_i->change_y(p_in_y);
2083     noeud_i->change_z(p_in_z);
2084     }
2085     }
2086     //Crit�re d'arr�t de l'algorithme
2087     compteur++;
2088     cout << "Fin de l'iteration #" << compteur << endl;
2089     }
2090     while (compteur < iter_max);
2091    
2092     cout << "Arret de la procedure apres " << compteur << " iterations" << endl;
2093     }
2094    
2095    
2096 picher 233 int MG_LISSAGE::extract_skin(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2,double frac_min, int *mai2_id)
2097 francois 222 {
2098    
2099     //�tape 0 - On commence par mettre � z�ro tous les nouveau_numero des triangles et des noeuds du maillage
2100     LISTE_MG_TRIANGLE::iterator it_tri;
2101     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
2102     {
2103     mgtri->change_nouveau_numero(0);
2104 francois 224 mgtri->change_origine(MAILLEUR_AUTO);
2105 francois 222 }
2106     LISTE_MG_NOEUD::iterator it_noeud;
2107     for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
2108     {
2109     mgnoeud->change_nouveau_numero(0);
2110     }
2111    
2112     //�tape 1 - On boucle ensuite tous les t�tra�dres pour ajouter un compteur du nombre de fois qu'un triangle appartient � 1 t�tra
2113 picher 233 LISTE_MG_TETRA::iterator it_tetra;
2114     for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
2115 francois 222 {
2116     int origine = mgtet->get_origine();
2117 francois 224 if (origine==IMPOSE)
2118 francois 222 {
2119 francois 224 mgtet->get_triangle1()->change_origine(IMPOSE);
2120     mgtet->get_triangle2()->change_origine(IMPOSE);
2121     mgtet->get_triangle3()->change_origine(IMPOSE);
2122     mgtet->get_triangle4()->change_origine(IMPOSE);
2123     }
2124     if (origine==OPTIMISE)
2125     {
2126     if (mgtet->get_triangle1()->get_origine()!=IMPOSE) mgtet->get_triangle1()->change_origine(OPTIMISE);
2127     if (mgtet->get_triangle2()->get_origine()!=IMPOSE) mgtet->get_triangle2()->change_origine(OPTIMISE);
2128     if (mgtet->get_triangle3()->get_origine()!=IMPOSE) mgtet->get_triangle3()->change_origine(OPTIMISE);
2129     if (mgtet->get_triangle4()->get_origine()!=IMPOSE) mgtet->get_triangle4()->change_origine(OPTIMISE);
2130     }
2131    
2132     if (((origine == OPTIMISE) && (etat[(OPTIMISE-1000)/10]==1) ) ||
2133     ((origine == IMPOSE) && (etat[(IMPOSE-1000)/10]==1) ) ||
2134     ((origine == MAILLEUR_AUTO) && (etat[(MAILLEUR_AUTO-1000)/10]==1) ) ||
2135     ((origine == TRIANGULATION) && (etat[(TRIANGULATION-1000)/10]==1) ) ||
2136     ((origine == MODIFICATION) && (etat[(MODIFICATION-1000)/10]==1) ) ||
2137     ((origine == DUPLIQUER) && (etat[(DUPLIQUER-1000)/10]==1) ) )
2138    
2139     {
2140 francois 222 int num1 = mgtet->get_triangle1()->get_nouveau_numero();
2141     int num2 = mgtet->get_triangle2()->get_nouveau_numero();
2142     int num3 = mgtet->get_triangle3()->get_nouveau_numero();
2143     int num4 = mgtet->get_triangle4()->get_nouveau_numero();
2144     num1++;
2145     num2++;
2146     num3++;
2147     num4++;
2148     mgtet->get_triangle1()->change_nouveau_numero(num1);
2149     mgtet->get_triangle2()->change_nouveau_numero(num2);
2150     mgtet->get_triangle3()->change_nouveau_numero(num3);
2151     mgtet->get_triangle4()->change_nouveau_numero(num4);
2152     }
2153     }
2154    
2155     //�tape 2 - On boucle l'ensemble des triangles identifi�s � l'�tape 1 pour identifier les noeuds leur appartenant
2156     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
2157     {
2158     int num = mgtri->get_nouveau_numero();
2159     if (num == 1)
2160     {
2161     MG_NOEUD* noeud1 = mgtri->get_noeud1();
2162     MG_NOEUD* noeud2 = mgtri->get_noeud2();
2163     MG_NOEUD* noeud3 = mgtri->get_noeud3();
2164     noeud1->change_nouveau_numero(1);
2165     noeud2->change_nouveau_numero(1);
2166     noeud3->change_nouveau_numero(1);
2167 picher 248 if (mgtri->get_origine()==IMPOSE)
2168     {
2169     noeud1->change_origine(IMPOSE);
2170     noeud2->change_origine(IMPOSE);
2171     noeud3->change_origine(IMPOSE);
2172     }
2173 francois 222 }
2174     }
2175    
2176 francois 234 //�tape 3 - On cr�e un nouveau maillage pour la peau
2177 francois 222 MG_MAILLAGE* mai2 = new MG_MAILLAGE(NULL);
2178     gest2.ajouter_mg_maillage(mai2);
2179    
2180     // //�tape 4 - On boucle l'ensemble des noeuds identifi�s � l'�tape 2 pour les recr�er dans le second maillage
2181     for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
2182     {
2183     int num = mgnoeud->get_nouveau_numero();
2184     if (num == 1)
2185     {
2186     double x = mgnoeud->get_x();
2187     double y = mgnoeud->get_y();
2188     double z = mgnoeud->get_z();
2189 francois 224 int origine=TRIANGULATION;
2190     if (mgnoeud->get_origine()==IMPOSE) origine=IMPOSE;
2191     MG_NOEUD* noeud1 = new MG_NOEUD(NULL,x,y,z,origine);
2192 francois 222 mai2->ajouter_mg_noeud(noeud1);
2193     mgnoeud->change_nouveau_numero(noeud1->get_id());
2194 francois 234 noeud1->change_nouveau_numero(mgnoeud->get_id());
2195 francois 222 }
2196     }
2197    
2198     // //�tape 5 - On boucle l'ensemble des triangles identifi�s � l'�tape 1 pour les recr�er dans le maillage 2
2199 picher 233 for (MG_TETRA* mgtet=mg_mai->get_premier_tetra(it_tetra);mgtet!=NULL;mgtet=mg_mai->get_suivant_tetra(it_tetra))
2200 francois 222 {
2201 francois 224 int originetet=mgtet->get_origine();
2202     if (((originetet == OPTIMISE) && (etat[(OPTIMISE-1000)/10]==1) ) ||
2203     ((originetet == IMPOSE) && (etat[(IMPOSE-1000)/10]==1) ) ||
2204     ((originetet == MAILLEUR_AUTO) && (etat[(MAILLEUR_AUTO-1000)/10]==1) ) ||
2205     ((originetet == TRIANGULATION) && (etat[(TRIANGULATION-1000)/10]==1) ) ||
2206     ((originetet == MODIFICATION) && (etat[(MODIFICATION-1000)/10]==1) ) ||
2207     ((originetet == DUPLIQUER) && (etat[(DUPLIQUER-1000)/10]==1) ) )
2208 francois 222 {
2209 francois 224 MG_NOEUD* noeud1 = mgtet->get_noeud1();
2210     MG_NOEUD* noeud2 = mgtet->get_noeud2();
2211     MG_NOEUD* noeud3 = mgtet->get_noeud3();
2212     MG_NOEUD* noeud4 = mgtet->get_noeud4();
2213 francois 222 MG_NOEUD* node1 = mai2->get_mg_noeudid(noeud1->get_nouveau_numero());
2214     MG_NOEUD* node2 = mai2->get_mg_noeudid(noeud2->get_nouveau_numero());
2215     MG_NOEUD* node3 = mai2->get_mg_noeudid(noeud3->get_nouveau_numero());
2216 francois 224 MG_NOEUD* node4 = mai2->get_mg_noeudid(noeud4->get_nouveau_numero());
2217     MG_TRIANGLE* tri1=mgtet->get_triangle1();
2218     MG_TRIANGLE* tri2=mgtet->get_triangle2();
2219     MG_TRIANGLE* tri3=mgtet->get_triangle3();
2220     MG_TRIANGLE* tri4=mgtet->get_triangle4();
2221     if (tri1->get_nouveau_numero()==1)
2222     {
2223     int origine=TRIANGULATION;
2224     if (tri1->get_origine()==IMPOSE) origine=IMPOSE;
2225     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node2,node3,mai2,origine);
2226     tripeau->change_nouveau_numero(0);
2227     tri1->change_nouveau_numero(0);
2228     }
2229     if (tri2->get_nouveau_numero()==1)
2230     {
2231     int origine=TRIANGULATION;
2232     if (tri2->get_origine()==IMPOSE) origine=IMPOSE;
2233     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node4,node2,mai2,origine);
2234     tripeau->change_nouveau_numero(0);
2235     tri2->change_nouveau_numero(0);
2236     }
2237     if (tri3->get_nouveau_numero()==1)
2238     {
2239     int origine=TRIANGULATION;
2240     if (tri3->get_origine()==IMPOSE) origine=IMPOSE;
2241     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node2,node4,node3,mai2,origine);
2242     tripeau->change_nouveau_numero(0);
2243     tri3->change_nouveau_numero(0);
2244     }
2245     if (tri4->get_nouveau_numero()==1)
2246     {
2247     int origine=TRIANGULATION;
2248     if (tri4->get_origine()==IMPOSE) origine=IMPOSE;
2249     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node3,node4,mai2,origine);
2250     tripeau->change_nouveau_numero(0);
2251     tri4->change_nouveau_numero(0);
2252     }
2253 francois 222 }
2254 francois 224 }
2255     cout << mai2->get_nb_mg_triangle()<< " triangles" <<endl;
2256     // etape 6 recherche des vosins
2257     for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
2258     {
2259     MG_TRIANGLE_PEAU* tripeau=(MG_TRIANGLE_PEAU*)mgtri;
2260     MG_TRIANGLE_PEAU* voisin1=recherche_voisin(tripeau->get_noeud1(),tripeau->get_noeud2(),tripeau);
2261     MG_TRIANGLE_PEAU* voisin2=recherche_voisin(tripeau->get_noeud2(),tripeau->get_noeud3(),tripeau);
2262     MG_TRIANGLE_PEAU* voisin3=recherche_voisin(tripeau->get_noeud3(),tripeau->get_noeud1(),tripeau);
2263     tripeau->change_voisin1(voisin1);
2264     tripeau->change_voisin2(voisin2);
2265     tripeau->change_voisin3(voisin3);
2266     }
2267     // //�tape 7 - On recherche les peaux
2268     int fin;
2269     do
2270     {
2271     fin=1;
2272     for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
2273     {
2274     MG_TRIANGLE_PEAU *tripeau=(MG_TRIANGLE_PEAU*)mgtri;
2275     if (tripeau->get_nouveau_numero()==0)
2276     {
2277     fin=0;
2278     std::vector<MG_TRIANGLE_PEAU*> *peau=new std::vector<MG_TRIANGLE_PEAU*>;
2279     lst_peau.push_back(peau);
2280     tripeau->change_nouveau_numero(1);
2281     peau->push_back(tripeau);
2282     determine_peau(peau);
2283     }
2284     }
2285     }
2286     while (fin==0);
2287 picher 233
2288 francois 224 std::cout << lst_peau.size() << " peau" << endl;
2289     for (int i=0;i<lst_peau.size();i++)
2290     std::cout << " peau " << i << " " << lst_peau[i]->size() << " triangles" << endl;
2291     //test de manifold
2292     LISTE_MG_SEGMENT::iterator itseg;
2293     for (MG_SEGMENT* seg=mai2->get_premier_segment(itseg);seg!=NULL;seg=mai2->get_suivant_segment(itseg))
2294     seg->change_nouveau_numero(0);
2295     TPL_MAP_ENTITE<MG_SEGMENT*> nonmanifoldarete;
2296     TPL_MAP_ENTITE<MG_NOEUD*> nonmanifoldnoeud;
2297 francois 234 TPL_MAP_ENTITE<MG_NOEUD*> nonmanifoldnoeuddepuisarete;
2298 francois 224 int nbpeau=lst_peau.size();
2299     for (int i=0;i<nbpeau;i++)
2300     {
2301     int nbtri=lst_peau[i]->size();
2302     for (int j=0;j<nbtri;j++)
2303     {
2304     MG_TRIANGLE_PEAU* tripeau=(*lst_peau[i])[j];
2305     tripeau->get_segment1()->change_nouveau_numero(tripeau->get_segment1()->get_nouveau_numero()+1);
2306     tripeau->get_segment2()->change_nouveau_numero(tripeau->get_segment2()->get_nouveau_numero()+1);
2307     tripeau->get_segment3()->change_nouveau_numero(tripeau->get_segment3()->get_nouveau_numero()+1);
2308     if (tripeau->get_segment1()->get_nouveau_numero()>2)
2309     nonmanifoldarete.ajouter(tripeau->get_segment1());
2310     if (tripeau->get_segment2()->get_nouveau_numero()>2)
2311     nonmanifoldarete.ajouter(tripeau->get_segment2());
2312     if (tripeau->get_segment3()->get_nouveau_numero()>2)
2313     nonmanifoldarete.ajouter(tripeau->get_segment3());
2314     }
2315     }
2316 picher 233
2317 francois 224 int nbnonmanifoldarete=nonmanifoldarete.get_nb();
2318     for (int i=0;i<nbnonmanifoldarete;i++)
2319     {
2320     MG_SEGMENT* seg=nonmanifoldarete.get(i);
2321     MG_NOEUD* n1=seg->get_noeud1();
2322     MG_NOEUD* n2=seg->get_noeud2();
2323 francois 234 nonmanifoldnoeuddepuisarete.ajouter(n1);
2324     nonmanifoldnoeuddepuisarete.ajouter(n2);
2325 francois 224 }
2326 francois 234 LISTE_MG_NOEUD::iterator itnoeud;
2327     for (MG_NOEUD* no=mai2->get_premier_noeud(itnoeud);no!=NULL;no=mai2->get_suivant_noeud(itnoeud))
2328     {
2329     if (nonmanifoldnoeuddepuisarete.existe(no)) continue;
2330     if (est_non_manifold(no)) nonmanifoldnoeud.ajouter(no);
2331     }
2332 francois 224
2333     cout << "Non manifold par arete " << nonmanifoldarete.get_nb() << endl;
2334 picher 233 cout << "Non manifold par noeud " << nonmanifoldnoeud.get_nb() << endl;
2335 picher 230
2336 francois 234 //�tape 8 - Traitement des cas de non manifold
2337     //non manifold par ar�te
2338 picher 233 for (int i=0;i<nbnonmanifoldarete;i++)
2339     {
2340     MG_SEGMENT* segment=nonmanifoldarete.get(i);
2341 francois 234 MG_NOEUD* noeud1 = mg_mai->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
2342     MG_NOEUD* noeud2 = mg_mai->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
2343     int nb_tetra = noeud1->get_lien_tetra()->get_nb();
2344     for (int j=0;j<nb_tetra;j++)
2345 picher 233 {
2346 francois 234 MG_TETRA* mgtet =noeud1->get_lien_tetra()->get(j);
2347     int originetet=mgtet->get_origine();
2348     if (originetet == MAILLEUR_AUTO)
2349 picher 233 {
2350 francois 234 //On r�active le tetra si l'autre noeud du segment lui appartient aussi
2351     MG_NOEUD* no1 = mgtet->get_noeud1();
2352     MG_NOEUD* no2 = mgtet->get_noeud2();
2353     MG_NOEUD* no3 = mgtet->get_noeud3();
2354     MG_NOEUD* no4 = mgtet->get_noeud4();
2355     if((no1 == noeud2) ||(no2 == noeud2) || (no3 == noeud2) || (no4 == noeud2))
2356     mgtet->change_origine(OPTIMISE);
2357 picher 233 }
2358     }
2359     }
2360    
2361     //non manifold par noeud
2362 picher 231 int nbnonmanifoldnoeud=nonmanifoldnoeud.get_nb();
2363     for (int i=0;i<nbnonmanifoldnoeud;i++)
2364     {
2365 francois 234 MG_NOEUD* noeud=mg_mai->get_mg_noeudid(nonmanifoldnoeud.get(i)->get_nouveau_numero());
2366     int nb_tetra = noeud->get_lien_tetra()->get_nb();
2367     for (int j = 0;j<nb_tetra;j++)
2368 francois 345 {
2369 francois 234 MG_TETRA* mgtet =noeud->get_lien_tetra()->get(j);
2370     int originetet=mgtet->get_origine();
2371     if (originetet == MAILLEUR_AUTO)
2372 picher 231 {
2373 francois 234 mgtet->change_origine(OPTIMISE);
2374 picher 231 }
2375     }
2376     }
2377 francois 234 *mai2_id = mai2->get_id();
2378     if ((nbnonmanifoldarete != 0) || (nbnonmanifoldnoeud != 0))
2379     {
2380     for (int i=0;i<lst_peau.size();i++)
2381     {
2382     delete lst_peau[i];
2383     }
2384     lst_peau.clear();
2385     gest2.supprimer_mg_maillageid(*mai2_id);
2386     return 0;
2387     }
2388 picher 230
2389 francois 345 //�tape 9 - Suppression des peaux de tr�s faible tail -gamma 0.1 -liss3 -densite densiterockerle (mati�re flottante)
2390 picher 233 cout << "Suppression des peaux de faible taille" << endl;
2391     int nbtri_total = mai2->get_nb_mg_triangle();
2392 picher 230 for (int i=0;i<nbpeau;i++)
2393     {
2394     int nbtri=lst_peau[i]->size();
2395 picher 233 if (nbtri < frac_min*nbtri_total)
2396 picher 230 {
2397     for (int j=0;j<nbtri;j++)
2398     {
2399     MG_TRIANGLE_PEAU* tripeau=(*lst_peau[i])[j];
2400 picher 233 mai2->supprimer_mg_triangleid(tripeau->get_id());
2401 picher 230 }
2402     }
2403 picher 233 }
2404     cout << mai2->get_nb_mg_triangle()<< " triangles" <<endl;
2405    
2406 francois 234 return 1;
2407 francois 222 }
2408    
2409 francois 224 void MG_LISSAGE::determine_peau(std::vector<MG_TRIANGLE_PEAU*> * peau)
2410 francois 222 {
2411 francois 224 int num=0;
2412     while (num!=peau->size())
2413     {
2414     MG_TRIANGLE_PEAU* p=(*peau)[num];
2415     if (p->get_voisin1()->get_nouveau_numero()==0)
2416     {
2417     peau->push_back(p->get_voisin1());
2418     p->get_voisin1()->change_nouveau_numero(1);
2419     }
2420     if (p->get_voisin2()->get_nouveau_numero()==0)
2421     {
2422     peau->push_back(p->get_voisin2());
2423     p->get_voisin2()->change_nouveau_numero(1);
2424     }
2425     if (p->get_voisin3()->get_nouveau_numero()==0)
2426     {
2427     peau->push_back(p->get_voisin3());
2428     p->get_voisin3()->change_nouveau_numero(1);
2429     }
2430     num++;
2431     }
2432     }
2433    
2434    
2435     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)
2436     {
2437 francois 222 MG_TRIANGLE_PEAU* voisin1=NULL,*voisin2=NULL,*voisin3=NULL;
2438     MG_SEGMENT* mgsegment1=mg_maillage->get_mg_segment(mgnoeud1->get_id(),mgnoeud2->get_id());
2439     MG_SEGMENT* mgsegment2=mg_maillage->get_mg_segment(mgnoeud2->get_id(),mgnoeud3->get_id());
2440     MG_SEGMENT* mgsegment3=mg_maillage->get_mg_segment(mgnoeud3->get_id(),mgnoeud1->get_id());
2441     if (mgsegment1==NULL)
2442 francois 224 mgsegment1=mg_maillage->ajouter_mg_segment(topo,mgnoeud1,mgnoeud2,origine);
2443 francois 222 if (mgsegment2==NULL)
2444 francois 224 mgsegment2=mg_maillage->ajouter_mg_segment(topo,mgnoeud2,mgnoeud3,origine);
2445 francois 222 if (mgsegment3==NULL)
2446 francois 224 mgsegment3=mg_maillage->ajouter_mg_segment(topo,mgnoeud3,mgnoeud1,origine);
2447     MG_TRIANGLE_PEAU* mtriangle=new MG_TRIANGLE_PEAU(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,origine);
2448 francois 222 mg_maillage->ajouter_mg_triangle(mtriangle);
2449     return mtriangle;
2450     }
2451    
2452 francois 224 MG_TRIANGLE_PEAU* MG_LISSAGE::recherche_voisin(MG_NOEUD* mg_noeud1,MG_NOEUD* mg_noeud2,MG_TRIANGLE_PEAU* triref)
2453 francois 222 {
2454 francois 224 MG_TRIANGLE_PEAU* trisol=NULL;
2455     double angleref=4.*M_PI;
2456 francois 222 int nb1=mg_noeud1->get_lien_triangle()->get_nb();
2457     int nb2=mg_noeud2->get_lien_triangle()->get_nb();
2458     for (int i=0;i<nb1;i++)
2459     for (int j=0;j<nb2;j++)
2460 francois 224 if (mg_noeud1->get_lien_triangle()->get(i)==mg_noeud2->get_lien_triangle()->get(j))
2461     if ((MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i)!=triref)
2462     {
2463     double angle=calcul_angle(triref,(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i));
2464     if (angle<angleref)
2465     {
2466     angleref=angle;
2467     trisol=(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i);
2468     }
2469     }
2470     return trisol;
2471 francois 222 }
2472 francois 223
2473 francois 234 int MG_LISSAGE::est_non_manifold(MG_NOEUD* no)
2474     {
2475     static int compteur=0;
2476     compteur++;
2477     int nb_tri=no->get_lien_triangle()->get_nb();
2478     TPL_MAP_ENTITE<MG_TRIANGLE*> lst;
2479     for (int i=0;i<nb_tri;i++)
2480     lst.ajouter(no->get_lien_triangle()->get(i));
2481     MG_TRIANGLE_PEAU* premier_triangle=(MG_TRIANGLE_PEAU*)lst.get(0);
2482     lst.supprimer(premier_triangle);
2483     MG_TRIANGLE_PEAU* tricour=(MG_TRIANGLE_PEAU*)premier_triangle;
2484     int i=0;
2485     do
2486     {
2487     if (lst.existe(tricour->get_voisin1()) || ((tricour->get_voisin1()==premier_triangle)&& (i>1)))
2488     {
2489     tricour=tricour->get_voisin1();
2490     lst.supprimer(tricour);
2491     }
2492     else if (lst.existe(tricour->get_voisin2()) || ((tricour->get_voisin2()==premier_triangle)&& (i>1)))
2493     {
2494     tricour=tricour->get_voisin2();
2495     lst.supprimer(tricour);
2496     }
2497     else if (lst.existe(tricour->get_voisin3()) || ((tricour->get_voisin3()==premier_triangle)&& (i>1)))
2498     {
2499     tricour=tricour->get_voisin3();
2500     lst.supprimer(tricour);
2501     }
2502     i++;
2503     }
2504     while (tricour!=premier_triangle);
2505     if (lst.get_nb()>0)
2506     return 1;
2507     return 0;
2508     }
2509 francois 224
2510     double MG_LISSAGE::calcul_angle(MG_TRIANGLE_PEAU* ft1,MG_TRIANGLE_PEAU* ft2)
2511 francois 223 {
2512 francois 388 MG_NOEUD* noeud1=ft1->get_noeud1();
2513     MG_NOEUD* noeud2=ft1->get_noeud2();
2514     MG_NOEUD* noeud3=ft1->get_noeud3();
2515     MG_NOEUD* noeud4=ft2->get_noeud1();
2516     MG_NOEUD* noeud5=ft2->get_noeud2();
2517     MG_NOEUD* noeud6=ft2->get_noeud3();
2518     double angle=get_angle_par_noeud<MG_NOEUD*>(noeud1,noeud2,noeud3,noeud4,noeud5,noeud6);
2519     return angle;
2520 francois 223 }
2521 picher 230
2522     double MG_LISSAGE::ponderation_gaussian(double s,double sigma)
2523     {
2524     double w_s;
2525     w_s = exp(-(s*s)/(2.*sigma*sigma));
2526     return w_s;
2527     }
2528    
2529     double MG_LISSAGE::ponderation_laplacian(double s,double sigma)
2530     {
2531     double w_s;
2532     w_s = exp(-(sqrt(2.)*fabs(s))/sigma);
2533     return w_s;
2534     }
2535    
2536     double MG_LISSAGE::ponderation_elfallahford(double s,double sigma)
2537     {
2538     double w_s;
2539     w_s = 1./sqrt(1+pow((s/sigma),2));
2540     return w_s;
2541     }