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

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