ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mg_lissage.cpp
Revision: 231
Committed: Fri Jan 15 17:20:10 2010 UTC (15 years, 4 months ago) by picher
Original Path: magic/lib/optimisation/optimisation/src/mg_lissage.cpp
File size: 25829 byte(s)
Log Message:
Modification de la hierarchisation de la classe mg_lissage pour l'ajout d'une methode menu

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     #include <stdio.h>
9 picher 230 #include <stdlib.h>
10     #include <time.h>
11 francois 222 #include <string.h>
12 francois 224 #include <math.h>
13 francois 222 //---------------------------------------------------------------------------
14    
15    
16 francois 224
17 francois 222 MG_LISSAGE::MG_LISSAGE()
18     {
19 francois 224 for (int i=0;i<9;i++) etat[i]=0;
20 francois 222 }
21    
22    
23     MG_LISSAGE::~MG_LISSAGE()
24     {
25 francois 224 for (int i=0;i<lst_peau.size();i++)
26     delete lst_peau[i];
27 francois 222 }
28    
29 francois 224 void MG_LISSAGE::conserve(int origine)
30     {
31     etat[(origine-1000)/10]=1;
32     }
33    
34 picher 231 void MG_LISSAGE::menu(FEM_MAILLAGE* mai, MG_GESTIONNAIRE& gest2, double epsilon, double sigma, int iter_max, int bruitage, int lissage)
35     {
36     //Menu de la méthode de lissage
37    
38     cout << "Extraction du maillage de surface" << endl;
39     int coderesu = 0;
40     MG_LISSAGE liss;
41     liss.conserve(OPTIMISE);
42     liss.conserve(IMPOSE);
43     liss.conserve(MAILLEUR_AUTO); //si la pièce de départ ne provient pas d'un calcul d'optimisation il faut utiliser cette ligne
44     do
45     {
46     coderesu = liss.extract_skin(mai,gest2);
47     cout << "coderesu = " << coderesu << endl;
48     }
49     while (coderesu == 0);
50     MG_MAILLAGE* mg_mai=gest2.get_mg_maillageid(1);
51     if (bruitage == 1)
52     {
53     cout << "Bruitage du maillage initial" << endl;
54     liss.bruitage(mg_mai,gest2);
55     }
56     if (lissage == 1)
57     {
58     cout << "Procédure de lissage" << endl;
59     liss.lissage(mg_mai,gest2,epsilon,sigma,iter_max);
60     }
61    
62     }
63    
64 picher 230 void MG_LISSAGE::bruitage(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2)
65     {
66     //Méthode test pour bruiter un maillage initial et tester la méthode de lissage
67     srand(time(NULL));
68     LISTE_MG_NOEUD::iterator it_no;
69     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
70     {
71     //On modifie la position des noeuds légèrement (ajout d'un bruit au maillage initial)
72     int signe;
73     int delta;
74     double delta_x,delta_y,delta_z;
75     signe = rand()%2 + 1;
76     delta = rand()%2 + 1;
77     if(signe == 1)
78     {
79     delta_x = delta*0.5;
80     }
81     else //(signe == 2)
82     {
83     delta_x = -delta*0.5;
84     }
85     signe = rand()%2 + 1;
86     delta = rand()%2 + 1;
87     if(signe == 1)
88     {
89     delta_y = delta*0.5;
90     }
91     else //(signe == 2)
92     {
93     delta_y = -delta*0.5;
94     }
95     signe = rand()%2 + 1;
96     delta = rand()%2 + 1;
97     if(signe == 1)
98     {
99     delta_z = delta*0.5;
100     }
101     else //(signe == 2)
102     {
103     delta_z = -delta*0.5;
104     }
105     double x_new = noeud_i->get_x() + delta_x;
106     double y_new = noeud_i->get_y() + delta_y;
107     double z_new = noeud_i->get_z() + delta_z;
108     noeud_i->change_x(x_new);
109     noeud_i->change_y(y_new);
110     noeud_i->change_z(z_new);
111     }
112     }
113 francois 224
114 picher 230 void MG_LISSAGE::lissage(MG_MAILLAGE* mg_mai,MG_GESTIONNAIRE& gest2, double epsilon, double sigma, int iter_max)
115     {
116     //Méthode de lissage inspirée de Chen(2005)
117     double un_sur_pi = 1./M_PI;
118     int compteur = 0;
119     int fin = 0;
120 picher 231 FILE *out = fopen("convergence.txt","wt");
121 francois 224
122 picher 231
123 picher 230 do
124     {
125     TPL_LISTE_ENTITE<OT_VECTEUR_3D> liste_normales;
126     TPL_LISTE_ENTITE<double> liste_wij;
127     LISTE_MG_TRIANGLE::iterator it_tri;
128     int k = 0; //pour identifier les triangles pour liste_normales et liste_wij
129     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
130     {
131     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
132     OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
133     //Remplissage de la liste des voisins du triangle i
134     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*> liste_voisins;
135     MG_NOEUD* noeud1 = mgtri_i->get_noeud1();
136     double nb_voisins1 = noeud1->get_lien_triangle()->get_nb();
137     for (int j = 0;j<nb_voisins1;j++)
138     {
139     MG_TRIANGLE_PEAU* mgtri_1 = (MG_TRIANGLE_PEAU*) noeud1->get_lien_triangle()->get(j);
140     liste_voisins.ajouter(mgtri_1);
141     }
142     MG_NOEUD* noeud2 = mgtri_i->get_noeud2();
143     double nb_voisins2 = noeud2->get_lien_triangle()->get_nb();
144     for (int j = 0;j<nb_voisins2;j++)
145     {
146     MG_TRIANGLE_PEAU* mgtri_2 = (MG_TRIANGLE_PEAU*) noeud2->get_lien_triangle()->get(j);
147     liste_voisins.ajouter(mgtri_2);
148     }
149     MG_NOEUD* noeud3 = mgtri_i->get_noeud3();
150     double nb_voisins3 = noeud3->get_lien_triangle()->get_nb();
151     for (int j = 0;j<nb_voisins3;j++)
152     {
153     MG_TRIANGLE_PEAU* mgtri_3 = (MG_TRIANGLE_PEAU*) noeud3->get_lien_triangle()->get(j);
154     liste_voisins.ajouter(mgtri_3);
155     }
156     liste_voisins.supprimer(mgtri_i);
157     int nb_voisins = liste_voisins.get_nb();
158     double w_ij = 1./nb_voisins;
159     double phi_i_min = 10.;
160     double s_i = 0.0;
161     double phi_im = 0.0;
162 picher 231 double *phi_ij = new double[nb_voisins];
163 picher 230 OT_VECTEUR_3D normal_f_i_mean(0.,0.,0.);
164     OT_VECTEUR_3D eta_i(0.,0.,0.);
165     TPL_MAP_ENTITE<MG_TRIANGLE_PEAU*>::ITERATEUR it;
166     int j = 0;
167     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
168     {
169     OT_VECTEUR_3D normal_f_j = mgtri_j->calcul_normal();
170     //1-Calculer la normale moyenne pour chaque triangle
171     normal_f_i_mean = normal_f_i_mean + w_ij*normal_f_j;
172     //2.1-On calcul l'angle entre normal_f_i et normal_f_j pour j allant de 1 à Nb_voisins
173     double prod_scalaire = normal_f_i*normal_f_j;
174     if (fabs(prod_scalaire) > 1.) prod_scalaire = 1.;
175     phi_ij[j] = acos(prod_scalaire)*un_sur_pi;
176     //2.2-On trouve le plus petit des angles et la normale heta_i correspondante
177     if (phi_ij[j] < phi_i_min)
178     {
179     phi_i_min = phi_ij[j];
180     eta_i = normal_f_j;
181     }
182     //3.1-On calcul l'angle moyen phi_im
183     phi_im = phi_im + w_ij*phi_ij[j];
184     j++;
185     }
186     normal_f_i_mean.norme();
187     j = 0;
188     for(MG_TRIANGLE_PEAU* mgtri_j=liste_voisins.get_premier(it);mgtri_j!=NULL;mgtri_j=liste_voisins.get_suivant(it))
189     {
190     //3.2-Calcul de s_i selon la variance
191     s_i = s_i + w_ij*pow((phi_ij[j] - phi_im),2);
192     j++;
193     }
194 picher 231 delete[] phi_ij;
195 picher 230 //4-On calcule une nouvelle normale pour chaque triangle
196     OT_VECTEUR_3D normal_f_i_new = ponderation_gaussian(s_i,sigma)*normal_f_i_mean + (1. - ponderation_gaussian(s_i,sigma))*eta_i;
197     normal_f_i_new.norme();
198     liste_normales.ajouter(normal_f_i_new);
199     liste_wij.ajouter(w_ij);
200     mgtri->change_nouveau_numero(k);
201     k++;
202     }
203    
204     LISTE_MG_NOEUD::iterator it_no;
205     for (MG_NOEUD* noeud_i=mg_mai->get_premier_noeud(it_no);noeud_i!=NULL;noeud_i=mg_mai->get_suivant_noeud(it_no))
206     {
207     int nb_voisins_j = noeud_i->get_lien_triangle()->get_nb();
208     double w_ij_prime = 0.0;
209     OT_VECTEUR_3D v_temp(0.,0.,0.);
210     OT_VECTEUR_3D v_i(noeud_i->get_x(),noeud_i->get_y(),noeud_i->get_z());
211     for(int j=0;j<nb_voisins_j;j++)
212     {
213     MG_TRIANGLE_PEAU* mgtri_j = (MG_TRIANGLE_PEAU*) noeud_i->get_lien_triangle()->get(j);
214     //On calcule le centroïde cj du triangle mgtri_j
215     MG_NOEUD* n1 = mgtri_j->get_noeud1();
216     MG_NOEUD* n2 = mgtri_j->get_noeud2();
217     MG_NOEUD* n3 = mgtri_j->get_noeud3();
218     double cj_x = 0.333333333333333*(n1->get_x() + n2->get_x() + n3->get_x());
219     double cj_y = 0.333333333333333*(n1->get_y() + n2->get_y() + n3->get_y());
220     double cj_z = 0.333333333333333*(n1->get_z() + n2->get_z() + n3->get_z());
221     //On forme le vecteur vi_cj
222     OT_VECTEUR_3D vi_cj(cj_x - noeud_i->get_x(),cj_y - noeud_i->get_y(),cj_z - noeud_i->get_z());
223     OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_j->get_nouveau_numero());
224     double w_ij = liste_wij.get(mgtri_j->get_nouveau_numero());
225     w_ij_prime = w_ij_prime + w_ij;
226     v_temp = v_temp + w_ij*(vi_cj*normal_f_i_new)*normal_f_i_new;
227     }
228     //5-On met à jour la position des noeuds
229     v_i = v_i + v_temp/w_ij_prime;
230     noeud_i->change_x(v_i.get_x());
231     noeud_i->change_y(v_i.get_y());
232     noeud_i->change_z(v_i.get_z());
233     }
234    
235     //Critère d'arrêt de l'algorithme
236     int l=0;
237     int nb_tri = mg_mai->get_nb_mg_triangle();
238     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
239     {
240     MG_TRIANGLE_PEAU* mgtri_i = (MG_TRIANGLE_PEAU*)mgtri;
241     OT_VECTEUR_3D normal_f_i = mgtri_i->calcul_normal();
242     OT_VECTEUR_3D normal_f_i_new = liste_normales.get(mgtri_i->get_nouveau_numero());
243     double critere = 1. - normal_f_i*normal_f_i_new;
244 picher 231 fprintf(out,"%lf\n",critere);
245 picher 230 if (critere <= epsilon) l++;
246     }
247     if (l == nb_tri) fin = 1;
248     compteur++;
249     cout << "Fin de l'itération #" << compteur << endl;
250 picher 231 fprintf(out,"\nFin de l'itération #%d\n\n",compteur);
251 picher 230 }
252     while ((fin == 0) && (compteur < iter_max));
253    
254     if (fin == 1)
255     {
256     cout << "Convergence de la méthode en " << compteur << " itérations" << endl;
257     }
258     else
259     {
260     cout << "Arrêt de la procédure après " << compteur << " itérations" << endl;
261     }
262    
263 picher 231 fclose(out);
264 picher 230 } //Fin de lissage
265    
266 picher 231 int MG_LISSAGE::extract_skin(FEM_MAILLAGE* mai,MG_GESTIONNAIRE& gest2)
267 francois 222 {
268    
269     //�tape 0 - On commence par mettre � z�ro tous les nouveau_numero des triangles et des noeuds du maillage
270     MG_MAILLAGE* mg_mai = (MG_MAILLAGE*)mai->get_mg_maillage();
271     LISTE_MG_TRIANGLE::iterator it_tri;
272     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
273     {
274     mgtri->change_nouveau_numero(0);
275 francois 224 mgtri->change_origine(MAILLEUR_AUTO);
276 francois 222 }
277     LISTE_MG_NOEUD::iterator it_noeud;
278     for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
279     {
280     mgnoeud->change_nouveau_numero(0);
281     }
282    
283     //�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
284     LISTE_FEM_TETRA::iterator it_tetra;
285     for (FEM_TETRA* tet=mai->get_premier_tetra(it_tetra);tet!=NULL;tet=mai->get_suivant_tetra(it_tetra))
286     {
287     MG_TETRA* mgtet=(MG_TETRA*)tet->get_mg_element_maillage();
288     int origine = mgtet->get_origine();
289 francois 224 if (origine==IMPOSE)
290 francois 222 {
291 francois 224 mgtet->get_triangle1()->change_origine(IMPOSE);
292     mgtet->get_triangle2()->change_origine(IMPOSE);
293     mgtet->get_triangle3()->change_origine(IMPOSE);
294     mgtet->get_triangle4()->change_origine(IMPOSE);
295     }
296     if (origine==OPTIMISE)
297     {
298     if (mgtet->get_triangle1()->get_origine()!=IMPOSE) mgtet->get_triangle1()->change_origine(OPTIMISE);
299     if (mgtet->get_triangle2()->get_origine()!=IMPOSE) mgtet->get_triangle2()->change_origine(OPTIMISE);
300     if (mgtet->get_triangle3()->get_origine()!=IMPOSE) mgtet->get_triangle3()->change_origine(OPTIMISE);
301     if (mgtet->get_triangle4()->get_origine()!=IMPOSE) mgtet->get_triangle4()->change_origine(OPTIMISE);
302     }
303    
304     if (((origine == OPTIMISE) && (etat[(OPTIMISE-1000)/10]==1) ) ||
305     ((origine == IMPOSE) && (etat[(IMPOSE-1000)/10]==1) ) ||
306     ((origine == MAILLEUR_AUTO) && (etat[(MAILLEUR_AUTO-1000)/10]==1) ) ||
307     ((origine == TRIANGULATION) && (etat[(TRIANGULATION-1000)/10]==1) ) ||
308     ((origine == MODIFICATION) && (etat[(MODIFICATION-1000)/10]==1) ) ||
309     ((origine == DUPLIQUER) && (etat[(DUPLIQUER-1000)/10]==1) ) )
310    
311     {
312 francois 222 int num1 = mgtet->get_triangle1()->get_nouveau_numero();
313     int num2 = mgtet->get_triangle2()->get_nouveau_numero();
314     int num3 = mgtet->get_triangle3()->get_nouveau_numero();
315     int num4 = mgtet->get_triangle4()->get_nouveau_numero();
316     num1++;
317     num2++;
318     num3++;
319     num4++;
320     mgtet->get_triangle1()->change_nouveau_numero(num1);
321     mgtet->get_triangle2()->change_nouveau_numero(num2);
322     mgtet->get_triangle3()->change_nouveau_numero(num3);
323     mgtet->get_triangle4()->change_nouveau_numero(num4);
324     }
325     }
326    
327     //�tape 2 - On boucle l'ensemble des triangles identifi�s � l'�tape 1 pour identifier les noeuds leur appartenant
328     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
329     {
330     int num = mgtri->get_nouveau_numero();
331     if (num == 1)
332     {
333     MG_NOEUD* noeud1 = mgtri->get_noeud1();
334     MG_NOEUD* noeud2 = mgtri->get_noeud2();
335     MG_NOEUD* noeud3 = mgtri->get_noeud3();
336     noeud1->change_nouveau_numero(1);
337     noeud2->change_nouveau_numero(1);
338     noeud3->change_nouveau_numero(1);
339     }
340     }
341    
342     //�tape 3 - On cr�e un nouveau gestionnaire et un nouveau maillage
343     MG_MAILLAGE* mai2 = new MG_MAILLAGE(NULL);
344     gest2.ajouter_mg_maillage(mai2);
345    
346     // //�tape 4 - On boucle l'ensemble des noeuds identifi�s � l'�tape 2 pour les recr�er dans le second maillage
347     for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
348     {
349     int num = mgnoeud->get_nouveau_numero();
350     if (num == 1)
351     {
352     double x = mgnoeud->get_x();
353     double y = mgnoeud->get_y();
354     double z = mgnoeud->get_z();
355 francois 224 int origine=TRIANGULATION;
356     if (mgnoeud->get_origine()==IMPOSE) origine=IMPOSE;
357     MG_NOEUD* noeud1 = new MG_NOEUD(NULL,x,y,z,origine);
358 francois 222 mai2->ajouter_mg_noeud(noeud1);
359     mgnoeud->change_nouveau_numero(noeud1->get_id());
360     }
361     }
362    
363     // //�tape 5 - On boucle l'ensemble des triangles identifi�s � l'�tape 1 pour les recr�er dans le maillage 2
364 francois 224 for (FEM_TETRA* tet=mai->get_premier_tetra(it_tetra);tet!=NULL;tet=mai->get_suivant_tetra(it_tetra))
365 francois 222 {
366 francois 224 MG_TETRA* mgtet=(MG_TETRA*)tet->get_mg_element_maillage();
367     int originetet=mgtet->get_origine();
368     if (((originetet == OPTIMISE) && (etat[(OPTIMISE-1000)/10]==1) ) ||
369     ((originetet == IMPOSE) && (etat[(IMPOSE-1000)/10]==1) ) ||
370     ((originetet == MAILLEUR_AUTO) && (etat[(MAILLEUR_AUTO-1000)/10]==1) ) ||
371     ((originetet == TRIANGULATION) && (etat[(TRIANGULATION-1000)/10]==1) ) ||
372     ((originetet == MODIFICATION) && (etat[(MODIFICATION-1000)/10]==1) ) ||
373     ((originetet == DUPLIQUER) && (etat[(DUPLIQUER-1000)/10]==1) ) )
374 francois 222 {
375 francois 224 MG_NOEUD* noeud1 = mgtet->get_noeud1();
376     MG_NOEUD* noeud2 = mgtet->get_noeud2();
377     MG_NOEUD* noeud3 = mgtet->get_noeud3();
378     MG_NOEUD* noeud4 = mgtet->get_noeud4();
379 francois 222 MG_NOEUD* node1 = mai2->get_mg_noeudid(noeud1->get_nouveau_numero());
380     MG_NOEUD* node2 = mai2->get_mg_noeudid(noeud2->get_nouveau_numero());
381     MG_NOEUD* node3 = mai2->get_mg_noeudid(noeud3->get_nouveau_numero());
382 francois 224 MG_NOEUD* node4 = mai2->get_mg_noeudid(noeud4->get_nouveau_numero());
383     MG_TRIANGLE* tri1=mgtet->get_triangle1();
384     MG_TRIANGLE* tri2=mgtet->get_triangle2();
385     MG_TRIANGLE* tri3=mgtet->get_triangle3();
386     MG_TRIANGLE* tri4=mgtet->get_triangle4();
387     if (tri1->get_nouveau_numero()==1)
388     {
389     int origine=TRIANGULATION;
390     if (tri1->get_origine()==IMPOSE) origine=IMPOSE;
391     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node2,node3,mai2,origine);
392     tripeau->change_nouveau_numero(0);
393     tri1->change_nouveau_numero(0);
394     }
395     if (tri2->get_nouveau_numero()==1)
396     {
397     int origine=TRIANGULATION;
398     if (tri2->get_origine()==IMPOSE) origine=IMPOSE;
399     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node4,node2,mai2,origine);
400     tripeau->change_nouveau_numero(0);
401     tri2->change_nouveau_numero(0);
402     }
403     if (tri3->get_nouveau_numero()==1)
404     {
405     int origine=TRIANGULATION;
406     if (tri3->get_origine()==IMPOSE) origine=IMPOSE;
407     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node2,node4,node3,mai2,origine);
408     tripeau->change_nouveau_numero(0);
409     tri3->change_nouveau_numero(0);
410     }
411     if (tri4->get_nouveau_numero()==1)
412     {
413     int origine=TRIANGULATION;
414     if (tri4->get_origine()==IMPOSE) origine=IMPOSE;
415     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node3,node4,mai2,origine);
416     tripeau->change_nouveau_numero(0);
417     tri4->change_nouveau_numero(0);
418     }
419 francois 222 }
420 francois 224 }
421     cout << mai2->get_nb_mg_triangle()<< " triangles" <<endl;
422     // etape 6 recherche des vosins
423     for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
424     {
425     MG_TRIANGLE_PEAU* tripeau=(MG_TRIANGLE_PEAU*)mgtri;
426     MG_TRIANGLE_PEAU* voisin1=recherche_voisin(tripeau->get_noeud1(),tripeau->get_noeud2(),tripeau);
427     MG_TRIANGLE_PEAU* voisin2=recherche_voisin(tripeau->get_noeud2(),tripeau->get_noeud3(),tripeau);
428     MG_TRIANGLE_PEAU* voisin3=recherche_voisin(tripeau->get_noeud3(),tripeau->get_noeud1(),tripeau);
429     tripeau->change_voisin1(voisin1);
430     tripeau->change_voisin2(voisin2);
431     tripeau->change_voisin3(voisin3);
432     }
433     // //�tape 7 - On recherche les peaux
434     int fin;
435     do
436     {
437     fin=1;
438     for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
439     {
440     MG_TRIANGLE_PEAU *tripeau=(MG_TRIANGLE_PEAU*)mgtri;
441     if (tripeau->get_nouveau_numero()==0)
442     {
443     fin=0;
444     std::vector<MG_TRIANGLE_PEAU*> *peau=new std::vector<MG_TRIANGLE_PEAU*>;
445     lst_peau.push_back(peau);
446     tripeau->change_nouveau_numero(1);
447     peau->push_back(tripeau);
448     determine_peau(peau);
449     }
450     }
451     }
452     while (fin==0);
453     std::cout << lst_peau.size() << " peau" << endl;
454     for (int i=0;i<lst_peau.size();i++)
455     std::cout << " peau " << i << " " << lst_peau[i]->size() << " triangles" << endl;
456     //test de manifold
457     for (MG_NOEUD* mgnoeud=mai2->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mai2->get_suivant_noeud(it_noeud))
458     mgnoeud->change_nouveau_numero(0);
459     LISTE_MG_SEGMENT::iterator itseg;
460     for (MG_SEGMENT* seg=mai2->get_premier_segment(itseg);seg!=NULL;seg=mai2->get_suivant_segment(itseg))
461     seg->change_nouveau_numero(0);
462     TPL_MAP_ENTITE<MG_SEGMENT*> nonmanifoldarete;
463     TPL_MAP_ENTITE<MG_NOEUD*> nonmanifoldnoeud;
464     int nbpeau=lst_peau.size();
465     for (int i=0;i<nbpeau;i++)
466     {
467     int nbtri=lst_peau[i]->size();
468     for (int j=0;j<nbtri;j++)
469     {
470     MG_TRIANGLE_PEAU* tripeau=(*lst_peau[i])[j];
471     tripeau->get_segment1()->change_nouveau_numero(tripeau->get_segment1()->get_nouveau_numero()+1);
472     tripeau->get_segment2()->change_nouveau_numero(tripeau->get_segment2()->get_nouveau_numero()+1);
473     tripeau->get_segment3()->change_nouveau_numero(tripeau->get_segment3()->get_nouveau_numero()+1);
474     if (tripeau->get_segment1()->get_nouveau_numero()>2)
475     nonmanifoldarete.ajouter(tripeau->get_segment1());
476     if (tripeau->get_segment2()->get_nouveau_numero()>2)
477     nonmanifoldarete.ajouter(tripeau->get_segment2());
478     if (tripeau->get_segment3()->get_nouveau_numero()>2)
479     nonmanifoldarete.ajouter(tripeau->get_segment3());
480     if ((tripeau->get_noeud1()->get_nouveau_numero()!=0) && (tripeau->get_noeud1()->get_nouveau_numero()!=i+1))
481     nonmanifoldnoeud.ajouter(tripeau->get_noeud1());
482     if ((tripeau->get_noeud2()->get_nouveau_numero()!=0) && (tripeau->get_noeud2()->get_nouveau_numero()!=i+1))
483     nonmanifoldnoeud.ajouter(tripeau->get_noeud2());
484     if ((tripeau->get_noeud3()->get_nouveau_numero()!=0) && (tripeau->get_noeud3()->get_nouveau_numero()!=i+1))
485     nonmanifoldnoeud.ajouter(tripeau->get_noeud3());
486     tripeau->get_noeud1()->change_nouveau_numero(i+1);
487     tripeau->get_noeud2()->change_nouveau_numero(i+1);
488     tripeau->get_noeud3()->change_nouveau_numero(i+1);
489     }
490     }
491     int nbnonmanifoldarete=nonmanifoldarete.get_nb();
492     for (int i=0;i<nbnonmanifoldarete;i++)
493     {
494     MG_SEGMENT* seg=nonmanifoldarete.get(i);
495     MG_NOEUD* n1=seg->get_noeud1();
496     MG_NOEUD* n2=seg->get_noeud2();
497     nonmanifoldnoeud.supprimer(n1);
498     nonmanifoldnoeud.supprimer(n2);
499     }
500 francois 222
501 francois 224
502     cout << "Non manifold par arete " << nonmanifoldarete.get_nb() << endl;
503     cout << "Non manifold par noeud " << nonmanifoldnoeud.get_nb() << endl;;
504 picher 230
505     //Étape 8 - Traitement des cas de non manifold
506 picher 231 int nbnonmanifoldnoeud=nonmanifoldnoeud.get_nb();
507     for (int i=0;i<nbnonmanifoldnoeud;i++)
508     {
509     MG_NOEUD* no=nonmanifoldnoeud.get(i);
510     for (FEM_TETRA* tet=mai->get_premier_tetra(it_tetra);tet!=NULL;tet=mai->get_suivant_tetra(it_tetra))
511     {
512     MG_TETRA* mgtet=(MG_TETRA*)tet->get_mg_element_maillage();
513     int originetet=mgtet->get_origine();
514     if (originetet == MAILLEUR_AUTO)
515     {
516     MG_NOEUD* noeud1 = mgtet->get_noeud1();
517     MG_NOEUD* noeud2 = mgtet->get_noeud2();
518     MG_NOEUD* noeud3 = mgtet->get_noeud3();
519     MG_NOEUD* noeud4 = mgtet->get_noeud4();
520     if ((noeud1 == no) || (noeud2 == no) || (noeud3 == no) || (noeud4 == no))
521     {
522     //On réactive le tétra en question
523     mgtet->change_origine(OPTIMISE);
524     }
525     }
526     }
527     }
528 picher 230
529 picher 231 if ((nbnonmanifoldarete == 0) && (nbnonmanifoldnoeud == 0))
530     {
531     return 1;
532     }
533     else
534     {
535     return 0;
536     }
537    
538 picher 230 //Étape 9 - Suppression des peaux de très faible taille
539 picher 231 /*int nbtri_total = mai2->get_nb_mg_triangle();
540 picher 230 for (int i=0;i<nbpeau;i++)
541     {
542     int nbtri=lst_peau[i]->size();
543     if (nbtri < 0.005*nbtri_total)
544     {
545     for (int j=0;j<nbtri;j++)
546     {
547     MG_TRIANGLE_PEAU* tripeau=(*lst_peau[i])[j];
548     int id = tripeau->get_id();
549     cout << "tripeau id:" << id <<endl;
550     mai2->supprimer_mg_triangleid(id);
551     }
552     }
553 picher 231 }*/
554 francois 222 }
555    
556 francois 224
557     void MG_LISSAGE::determine_peau(std::vector<MG_TRIANGLE_PEAU*> * peau)
558 francois 222 {
559 francois 224 int num=0;
560     while (num!=peau->size())
561     {
562     MG_TRIANGLE_PEAU* p=(*peau)[num];
563     if (p->get_voisin1()->get_nouveau_numero()==0)
564     {
565     peau->push_back(p->get_voisin1());
566     p->get_voisin1()->change_nouveau_numero(1);
567     }
568     if (p->get_voisin2()->get_nouveau_numero()==0)
569     {
570     peau->push_back(p->get_voisin2());
571     p->get_voisin2()->change_nouveau_numero(1);
572     }
573     if (p->get_voisin3()->get_nouveau_numero()==0)
574     {
575     peau->push_back(p->get_voisin3());
576     p->get_voisin3()->change_nouveau_numero(1);
577     }
578     num++;
579     }
580     }
581    
582    
583     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)
584     {
585 francois 222 MG_TRIANGLE_PEAU* voisin1=NULL,*voisin2=NULL,*voisin3=NULL;
586     MG_SEGMENT* mgsegment1=mg_maillage->get_mg_segment(mgnoeud1->get_id(),mgnoeud2->get_id());
587     MG_SEGMENT* mgsegment2=mg_maillage->get_mg_segment(mgnoeud2->get_id(),mgnoeud3->get_id());
588     MG_SEGMENT* mgsegment3=mg_maillage->get_mg_segment(mgnoeud3->get_id(),mgnoeud1->get_id());
589     if (mgsegment1==NULL)
590 francois 224 mgsegment1=mg_maillage->ajouter_mg_segment(topo,mgnoeud1,mgnoeud2,origine);
591 francois 222 if (mgsegment2==NULL)
592 francois 224 mgsegment2=mg_maillage->ajouter_mg_segment(topo,mgnoeud2,mgnoeud3,origine);
593 francois 222 if (mgsegment3==NULL)
594 francois 224 mgsegment3=mg_maillage->ajouter_mg_segment(topo,mgnoeud3,mgnoeud1,origine);
595     MG_TRIANGLE_PEAU* mtriangle=new MG_TRIANGLE_PEAU(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,origine);
596 francois 222 mg_maillage->ajouter_mg_triangle(mtriangle);
597     return mtriangle;
598     }
599    
600 francois 224 MG_TRIANGLE_PEAU* MG_LISSAGE::recherche_voisin(MG_NOEUD* mg_noeud1,MG_NOEUD* mg_noeud2,MG_TRIANGLE_PEAU* triref)
601 francois 222 {
602 francois 224 MG_TRIANGLE_PEAU* trisol=NULL;
603     double angleref=4.*M_PI;
604 francois 222 int nb1=mg_noeud1->get_lien_triangle()->get_nb();
605     int nb2=mg_noeud2->get_lien_triangle()->get_nb();
606     for (int i=0;i<nb1;i++)
607     for (int j=0;j<nb2;j++)
608 francois 224 if (mg_noeud1->get_lien_triangle()->get(i)==mg_noeud2->get_lien_triangle()->get(j))
609     if ((MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i)!=triref)
610     {
611     double angle=calcul_angle(triref,(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i));
612     if (angle<angleref)
613     {
614     angleref=angle;
615     trisol=(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i);
616     }
617     }
618     return trisol;
619 francois 222 }
620 francois 223
621 francois 224
622    
623     double MG_LISSAGE::calcul_angle(MG_TRIANGLE_PEAU* ft1,MG_TRIANGLE_PEAU* ft2)
624 francois 223 {
625 francois 224 MG_NOEUD* noeud1=ft1->get_noeud1();
626     MG_NOEUD* noeud2=ft1->get_noeud2();
627     MG_NOEUD* noeud3=ft1->get_noeud3();
628     MG_NOEUD* noeud4=ft2->get_noeud1();
629     MG_NOEUD* noeud5=ft2->get_noeud2();
630     MG_NOEUD* noeud6=ft2->get_noeud3();
631    
632     MG_NOEUD* noeuda=NULL;
633     MG_NOEUD* noeudb=NULL;
634     MG_NOEUD* noeudc=NULL;
635     MG_NOEUD* noeudd=NULL;
636    
637    
638     if (!(noeud1==noeud4))
639     if (!(noeud1==noeud5))
640     if (!(noeud1==noeud6))
641     {
642     noeuda=noeud1;
643     noeudb=noeud2;
644     noeudc=noeud3;
645     }
646     if (!(noeud2==noeud4))
647     if (!(noeud2==noeud5))
648     if (!(noeud2==noeud6))
649     {
650     noeuda=noeud2;
651     noeudb=noeud1;
652     noeudc=noeud3;
653     }
654     if (!(noeud3==noeud4))
655     if (!(noeud3==noeud5))
656     if (!(noeud3==noeud6))
657     {
658     noeuda=noeud3;
659     noeudb=noeud1;
660     noeudc=noeud2;
661     }
662     if (noeuda==NULL) return 2*M_PI;
663     if (!(noeud4==noeudb))
664     if (!(noeud4==noeudc)) noeudd=noeud4;
665     if (!(noeud5==noeudb))
666     if (!(noeud5==noeudc)) noeudd=noeud5;
667     if (!(noeud6==noeudb))
668     if (!(noeud6==noeudc)) noeudd=noeud6;
669    
670     OT_VECTEUR_3D milieu(0.5*(noeudb->get_x()+noeudc->get_x()),0.5*(noeudb->get_y()+noeudc->get_y()),0.5*(noeudb->get_z()+noeudc->get_z()));
671     OT_VECTEUR_3D bc(noeudc->get_x()-noeudb->get_x(),noeudc->get_y()-noeudb->get_y(),noeudc->get_z()-noeudb->get_z());
672     bc.norme();
673     double xyz[3];
674     xyz[0]=noeuda->get_x();
675     xyz[1]=noeuda->get_y();
676     xyz[2]=noeuda->get_z();
677     OT_VECTEUR_3D n1(xyz,milieu.get_xyz());
678     n1.norme();
679     OT_VECTEUR_3D n=bc&n1;
680     n1=n&bc;
681     n1.norme();
682     xyz[0]=noeudd->get_x();
683     xyz[1]=noeudd->get_y();
684     xyz[2]=noeudd->get_z();
685     OT_VECTEUR_3D n2(xyz,milieu.get_xyz());
686     n2.norme();
687     n=bc&n2;
688     n2=n&bc;
689     n2.norme();
690    
691     double ps=n1*n2;
692     if (ps>1.) ps=1.;
693     if (ps<-1.) ps=-1.;
694     double angle=acos(ps);
695     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
696     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
697     OT_VECTEUR_3D nbnd(noeudd->get_x()-noeudb->get_x(),noeudd->get_y()-noeudb->get_y(),noeudd->get_z()-noeudb->get_z());
698     n=n1n2&n1n3;
699     n.norme();
700     nbnd.norme();
701    
702     ps=n*nbnd;
703     if (ps<-0.0001) angle=2*M_PI-angle;
704     return angle;
705    
706    
707 francois 223 }
708 picher 230
709     double MG_LISSAGE::ponderation_gaussian(double s,double sigma)
710     {
711     double w_s;
712     w_s = exp(-(s*s)/(2.*sigma*sigma));
713     return w_s;
714     }
715    
716     double MG_LISSAGE::ponderation_laplacian(double s,double sigma)
717     {
718     double w_s;
719     w_s = exp(-(sqrt(2.)*fabs(s))/sigma);
720     return w_s;
721     }
722    
723     double MG_LISSAGE::ponderation_elfallahford(double s,double sigma)
724     {
725     double w_s;
726     w_s = 1./sqrt(1+pow((s/sigma),2));
727     return w_s;
728     }