ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mgopt_posttraitement.cpp
Revision: 232
Committed: Fri Jan 15 20:42:24 2010 UTC (15 years, 4 months ago) by francois
Original Path: magic/lib/optimisation/optimisation/src/mg_lissage.cpp
File size: 25838 byte(s)
Log Message:
Correction de bud apparu dans le lissage : destruction d'une entité non attaché a la géométrie + version d'affichage d'info dans le lissage + routine de gestion du lissage

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