ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mg_lissage.cpp
Revision: 343
Committed: Mon Jun 18 17:53:32 2012 UTC (13 years, 2 months ago) by francois
File size: 74079 byte(s)
Log Message: 
ajout d'une procedure de lissage dans l'optimisation de topologie basée sur une peau extraite par iso densité

File Contents

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