ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mgopt_posttraitement.cpp
Revision: 383
Committed: Mon Jan 14 22:32:40 2013 UTC (12 years, 4 months ago) by francois
Original Path: magic/lib/optimisation/src/mg_lissage.cpp
File size: 75069 byte(s)
Log Message: 
Correction de bug dans les changements de decembre sur la nouvelle gestion des solutions

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