ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mgopt_posttraitement.cpp
Revision: 388
Committed: Tue Feb 12 15:43:42 2013 UTC (12 years, 3 months ago) by francois
Original Path: magic/lib/optimisation/src/mg_lissage.cpp
File size: 75496 byte(s)
Log Message: 
Meilleuire structuration du calcul de l'angle matière entre deux triangles

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