ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur_analyse.cpp
Revision: 429
Committed: Fri Sep 27 21:45:39 2013 UTC (11 years, 11 months ago) by francois
File size: 15481 byte(s)
Log Message:
Correction dans l'importation et l'analyse des triangulations. Attention un STL est ecrit en float et non en double

File Contents

# Content
1 #include "gestionversion.h"
2 #include "mailleur_analyse.h"
3 #include "mg_gestionnaire.h"
4 #include "m3d_triangle.h"
5 #include "tpl_fonctions_generiques.h"
6
7
8
9
10
11
12
13 MAILLEUR_ANALYSE::MAILLEUR_ANALYSE(MG_MAILLAGE* m,OT_CPU* comp):MAILLEUR(comp),mai(m),borne1(0.1),borne2(0.2),borne3(0.5),eps_angle_retourne(0.03)
14 {
15
16 }
17
18
19
20
21
22 MAILLEUR_ANALYSE::MAILLEUR_ANALYSE(MAILLEUR_ANALYSE &mdd):MAILLEUR(mdd),mai(mdd.mai),borne1(mdd.borne1),borne2(mdd.borne2),borne3(mdd.borne3)
23 {
24
25 }
26
27
28
29 MAILLEUR_ANALYSE::~MAILLEUR_ANALYSE()
30 {
31
32 }
33
34 void MAILLEUR_ANALYSE::change_borne(double val1,double val2,double val3)
35 {
36 borne1=val1;
37 borne2=val2;
38 borne3=val3;
39 }
40
41 void MAILLEUR_ANALYSE::get_borne(double &val1,double &val2,double &val3)
42 {
43 val1=borne1;
44 val2=borne2;
45 val3=borne3;
46 }
47
48 void MAILLEUR_ANALYSE::change_eps_angle_retourne(double val)
49 {
50 eps_angle_retourne=val;
51 }
52
53 double MAILLEUR_ANALYSE::get_eps_angle_retourne(void)
54 {
55 return eps_angle_retourne;
56 }
57
58
59 void MAILLEUR_ANALYSE::analyse(char *nom)
60 {
61 char mess[500];
62 sprintf(mess,"ANALYSEUR DE MG MAILLAGE 3D");
63 affiche(mess);
64 sprintf(mess,"---------------------------");
65 affiche(mess);
66 sprintf(mess,"");
67 affiche(mess);
68
69 LISTE_MG_NOEUD::iterator itn;
70 for (MG_NOEUD* no=mai->get_premier_noeud(itn);no!=NULL;no=mai->get_suivant_noeud(itn))
71 no->change_nouveau_numero(0);
72
73 LISTE_MG_SEGMENT::iterator its;
74 for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
75 seg->change_nouveau_numero(0);
76
77
78 LISTE_MG_TRIANGLE::iterator ittr;
79 for (MG_TRIANGLE* tri=mai->get_premier_triangle(ittr);tri!=NULL;tri=mai->get_suivant_triangle(ittr))
80 tri->change_nouveau_numero(0);
81
82 LISTE_MG_QUADRANGLE::iterator itq;
83 for (MG_QUADRANGLE* quad=mai->get_premier_quadrangle(itq);quad!=NULL;quad=mai->get_suivant_quadrangle(itq))
84 quad->change_nouveau_numero(0);
85
86 LISTE_MG_TETRA::iterator itt;
87 for (MG_TETRA* tet=mai->get_premier_tetra(itt);tet!=NULL;tet=mai->get_suivant_tetra(itt))
88 {
89 tet->get_triangle1()->change_nouveau_numero(tet->get_triangle1()->get_nouveau_numero()+1);
90 tet->get_triangle2()->change_nouveau_numero(tet->get_triangle2()->get_nouveau_numero()+1);
91 tet->get_triangle3()->change_nouveau_numero(tet->get_triangle3()->get_nouveau_numero()+1);
92 tet->get_triangle4()->change_nouveau_numero(tet->get_triangle4()->get_nouveau_numero()+1);
93 }
94
95 LISTE_MG_HEXA::iterator ith;
96 for (MG_HEXA* hex=mai->get_premier_hexa(ith);hex!=NULL;hex=mai->get_suivant_hexa(ith))
97 {
98 hex->get_quadrangle1()->change_nouveau_numero(hex->get_quadrangle1()->get_nouveau_numero()+1);
99 hex->get_quadrangle2()->change_nouveau_numero(hex->get_quadrangle2()->get_nouveau_numero()+1);
100 hex->get_quadrangle3()->change_nouveau_numero(hex->get_quadrangle3()->get_nouveau_numero()+1);
101 hex->get_quadrangle4()->change_nouveau_numero(hex->get_quadrangle4()->get_nouveau_numero()+1);
102 hex->get_quadrangle5()->change_nouveau_numero(hex->get_quadrangle5()->get_nouveau_numero()+1);
103 hex->get_quadrangle6()->change_nouveau_numero(hex->get_quadrangle6()->get_nouveau_numero()+1);
104 }
105 int nbtrifront=0;
106 for (MG_TRIANGLE* tri=mai->get_premier_triangle(ittr);tri!=NULL;tri=mai->get_suivant_triangle(ittr))
107 {
108 if (mai->get_nb_mg_tetra()==0)
109 tri->change_nouveau_numero(1);
110 if (tri->get_nouveau_numero()==1)
111 {
112 nbtrifront++;
113 tri->get_segment1()->change_nouveau_numero(tri->get_segment1()->get_nouveau_numero()+1);
114 tri->get_segment2()->change_nouveau_numero(tri->get_segment2()->get_nouveau_numero()+1);
115 tri->get_segment3()->change_nouveau_numero(tri->get_segment3()->get_nouveau_numero()+1);
116 tri->get_noeud1()->change_nouveau_numero(1);
117 tri->get_noeud2()->change_nouveau_numero(1);
118 tri->get_noeud3()->change_nouveau_numero(1);
119
120 }
121 }
122 int nbquadfront=0;
123 for (MG_QUADRANGLE* quad=mai->get_premier_quadrangle(itq);quad!=NULL;quad=mai->get_suivant_quadrangle(itq))
124 {
125 if (mai->get_nb_mg_hexa()==0)
126 quad->change_nouveau_numero(1);
127 if (quad->get_nouveau_numero()==1)
128 {
129 nbquadfront++;
130 quad->get_segment1()->change_nouveau_numero(quad->get_segment1()->get_nouveau_numero()+1);
131 quad->get_segment2()->change_nouveau_numero(quad->get_segment2()->get_nouveau_numero()+1);
132 quad->get_segment3()->change_nouveau_numero(quad->get_segment3()->get_nouveau_numero()+1);
133 quad->get_segment4()->change_nouveau_numero(quad->get_segment4()->get_nouveau_numero()+1);
134 quad->get_noeud1()->change_nouveau_numero(1);
135 quad->get_noeud2()->change_nouveau_numero(1);
136 quad->get_noeud3()->change_nouveau_numero(1);
137 quad->get_noeud4()->change_nouveau_numero(1);
138 }
139 }
140 int nbsegfront=0;
141 for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
142 if (seg->get_nouveau_numero()==2) nbsegfront++;
143 int nbnofront=0;
144 for (MG_NOEUD* no=mai->get_premier_noeud(itn);no!=NULL;no=mai->get_suivant_noeud(itn))
145 if (no->get_nouveau_numero()==1) nbnofront++;
146 sprintf(mess," Maillage numero : %lu",mai->get_id());affiche(mess);
147 sprintf(mess," %d noeuds",mai->get_nb_mg_noeud()); affiche(mess);
148 sprintf(mess," %d segments",mai->get_nb_mg_segment()); affiche(mess);
149 sprintf(mess," %d triangles",mai->get_nb_mg_triangle()); affiche(mess);
150 sprintf(mess," %d quadrangles",mai->get_nb_mg_quadrangle()); affiche(mess);
151 sprintf(mess," %d tetras",mai->get_nb_mg_tetra()); affiche(mess);
152 sprintf(mess," %d hexas",mai->get_nb_mg_hexa()); affiche(mess); affiche("");
153 sprintf(mess," %d noeuds frontiere",nbnofront); affiche(mess);
154 sprintf(mess," %d segments frontiere",nbsegfront); affiche(mess);
155 sprintf(mess," %d triangles frontiere",nbtrifront); affiche(mess);
156 sprintf(mess," %d quadrangles frontiere",nbquadfront); affiche(mess);affiche("");
157
158 MG_SOLUTION *sol=NULL;
159 if (nom!=NULL)
160 {
161 std::string nomsol=nom;nomsol=nomsol+".sol";
162 if (mai->get_nb_mg_tetra()==0) sol=new MG_SOLUTION(mai,1,(char*)nomsol.c_str(),mai->get_nb_mg_triangle(),"Qualité",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2,MAGIC::TYPE_SOLUTION::SCALAIRE);
163 if (mai->get_nb_mg_tetra()!=0) sol=new MG_SOLUTION(mai,1,(char*)nomsol.c_str(),mai->get_nb_mg_tetra(),"Qualité",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT3,MAGIC::TYPE_SOLUTION::SCALAIRE);
164 mai->get_gestionnaire()->ajouter_mg_solution(sol);
165 }
166 double qualmin=1e300,qualmax=-1e300,qualmoy=0.;
167 int tab2[7]={0,0,0,0,0,0,0};
168 int itri=0;
169 for (MG_TRIANGLE* tri=mai->get_premier_triangle(ittr);tri!=NULL;tri=mai->get_suivant_triangle(ittr))
170 if (tri->get_nouveau_numero()==1)
171 {
172 double qual=OPERATEUR::qualite_triangle(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord(),tri->get_noeud3()->get_coord());
173 MG_NOEUD * no1=tri->get_noeud1();
174 MG_NOEUD * no2=tri->get_noeud2();
175 MG_NOEUD * no3=tri->get_noeud3();
176 if (no1==no2)
177 if (no1==no3)
178 tab2[6]++;
179 int nb1=no1->get_lien_triangle()->get_nb();
180 int nb2=no2->get_lien_triangle()->get_nb();
181 int nb3=no3->get_lien_triangle()->get_nb();
182 OT_VECTEUR_3D vec1(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord());
183 OT_VECTEUR_3D vec2(tri->get_noeud1()->get_coord(),tri->get_noeud3()->get_coord());
184 OT_VECTEUR_3D normal=vec1&vec2;
185 normal.norme();
186 int nbretourne=0;
187 int nbvoisin=0;
188 for (int i=0;i<nb1;i++)
189 for (int j=0;j<nb2;j++)
190 {
191 MG_TRIANGLE* tri1=no1->get_lien_triangle()->get(i);
192 MG_TRIANGLE* tri2=no2->get_lien_triangle()->get(j);
193 if ( (tri1==tri2) && (tri1!=tri) && (tri1->get_nouveau_numero()==1))
194 {
195 OT_VECTEUR_3D vec1tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud2()->get_coord());
196 OT_VECTEUR_3D vec2tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud3()->get_coord());
197 OT_VECTEUR_3D normaltmp=vec1tmp&vec2tmp;
198 normaltmp.norme();
199 double psca=normal*normaltmp;
200 if (psca<-cos(eps_angle_retourne)) nbretourne++;
201 nbvoisin++;
202 }
203 }
204 for (int i=0;i<nb1;i++)
205 for (int j=0;j<nb3;j++)
206 {
207 MG_TRIANGLE* tri1=no1->get_lien_triangle()->get(i);
208 MG_TRIANGLE* tri2=no3->get_lien_triangle()->get(j);
209 if ( (tri1==tri2) && (tri1!=tri) && (tri1->get_nouveau_numero()==1))
210 {
211 OT_VECTEUR_3D vec1tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud2()->get_coord());
212 OT_VECTEUR_3D vec2tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud3()->get_coord());
213 OT_VECTEUR_3D normaltmp=vec1tmp&vec2tmp;
214 normaltmp.norme();
215 double psca=normal*normaltmp;
216 if (psca<-cos(eps_angle_retourne)) nbretourne++;
217 nbvoisin++;
218 }
219 }
220 for (int i=0;i<nb2;i++)
221 for (int j=0;j<nb3;j++)
222 {
223 MG_TRIANGLE* tri1=no2->get_lien_triangle()->get(i);
224 MG_TRIANGLE* tri2=no3->get_lien_triangle()->get(j);
225 if ( (tri1==tri2) && (tri1!=tri) && (tri1->get_nouveau_numero()==1))
226 {
227 OT_VECTEUR_3D vec1tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud2()->get_coord());
228 OT_VECTEUR_3D vec2tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud3()->get_coord());
229 OT_VECTEUR_3D normaltmp=vec1tmp&vec2tmp;
230 normaltmp.norme();
231 double psca=normal*normaltmp;
232 if (psca<-cos(eps_angle_retourne)) nbretourne++;
233 nbvoisin++;
234 }
235 }
236 if (nbvoisin!=3)
237 tab2[5]++;
238 if (nbretourne>1) qual=-qual;
239 qualmoy=qualmoy+qual;
240 if (qual<qualmin) qualmin=qual;
241 if (qual>qualmax) qualmax=qual;
242 if (qual<0.) tab2[0]++;
243 else if (qual<borne1) tab2[1]++;
244 else if (qual<borne2) tab2[2]++;
245 else if (qual<borne3) tab2[3]++;
246 else tab2[4]++;
247 if (sol!=NULL)
248 if (mai->get_nb_mg_tetra()==0) sol->ecrire(qual,itri,0);
249 itri++;
250 }
251 if (nbtrifront!=0)
252 {
253 sprintf(mess," Nombre de triangles avec des noeuds fusionnés : %d",tab2[6]); affiche(mess);
254 sprintf(mess," Nombre de triangles avec 1 voisin manquant : %d",tab2[5]); affiche(mess);
255 sprintf(mess," critere d'angle limite pour l'inversion : %lf",eps_angle_retourne); affiche(mess);
256 sprintf(mess," qualite moyenne des triangles de frontiere : %lf",qualmoy/nbtrifront); affiche(mess);
257 sprintf(mess," qualite min des triangles de frontiere : %lf",qualmin); affiche(mess);
258 sprintf(mess," qualite max des triangles de frontiere : %lf",qualmax); affiche(mess);
259 sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne3,tab2[4],tab2[4]*100./nbtrifront); affiche(mess);
260 sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne2,tab2[3],tab2[3]*100./nbtrifront); affiche(mess);
261 sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne1,tab2[2],tab2[2]*100./nbtrifront); affiche(mess);
262 sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",0.,tab2[1],tab2[1]*100./nbtrifront); affiche(mess);
263 sprintf(mess," nombre de triangles de frontiere retournes : %d (%.2lf%%)",tab2[0],tab2[0]*100./nbtrifront); affiche(mess); affiche("");
264 }
265
266
267
268 qualmin=1e300,qualmax=-1e300,qualmoy=0.;
269 int tab3[5]={0,0,0,0,0};
270 int itet=0;
271 for (MG_TETRA* tet=mai->get_premier_tetra(itt);tet!=NULL;tet=mai->get_suivant_tetra(itt))
272 {
273 double qual=OPERATEUR::qualite_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord());
274 OT_VECTEUR_3D vec1(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord());
275 OT_VECTEUR_3D vec2(tet->get_noeud1()->get_coord(),tet->get_noeud3()->get_coord());
276 OT_VECTEUR_3D vec3(tet->get_noeud1()->get_coord(),tet->get_noeud4()->get_coord());
277 OT_VECTEUR_3D pvec=vec1&vec2;
278 double psca=pvec*vec3;
279 if (psca<0.) qual=-qual;
280 qualmoy=qualmoy+qual;
281 if (qual<qualmin) qualmin=qual;
282 if (qual>qualmax) qualmax=qual;
283 if (qual<0.) tab3[0]++;
284 else if (qual<borne1) tab3[1]++;
285 else if (qual<borne2) tab3[2]++;
286 else if (qual<borne3) tab3[3]++;
287 else tab3[4]++;
288 if (sol!=NULL)
289 sol->ecrire(qual,itet,0);
290 itet++;
291 }
292 if (mai->get_nb_mg_tetra()!=0)
293 {
294 sprintf(mess," qualite moyenne des tetra : %lf",qualmoy/mai->get_nb_mg_tetra()); affiche(mess);
295 sprintf(mess," qualite min des tetra : %lf",qualmin); affiche(mess);
296 sprintf(mess," qualite max des tetra : %lf",qualmax); affiche(mess);
297 sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne3,tab3[4],tab3[4]*100./mai->get_nb_mg_tetra()); affiche(mess);
298 sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne2,tab3[3],tab3[3]*100./mai->get_nb_mg_tetra()); affiche(mess);
299 sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne1,tab3[2],tab3[2]*100./mai->get_nb_mg_tetra()); affiche(mess);
300 sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",0.,tab3[1],tab3[1]*100./mai->get_nb_mg_tetra()); affiche(mess);
301 sprintf(mess," nombre de tetra retournes : %d (%.2lf%%)",tab3[0],tab3[0]*100./mai->get_nb_mg_tetra()); affiche(mess);
302 int nbsegcorrect=0;
303 int nbsegincorrect=0;
304 for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
305 if (seg->get_nouveau_numero()!=2)
306 {
307 MG_NOEUD* no1=seg->get_noeud1();
308 MG_NOEUD* no2=seg->get_noeud2();
309 int nb1=no1->get_lien_tetra()->get_nb();
310 int nb2=no2->get_lien_tetra()->get_nb();
311 double angletot=0;
312 for (int i=0;i<nb1;i++)
313 for (int j=0;j<nb2;j++)
314 {
315 MG_TETRA* tet1=no1->get_lien_tetra()->get(i);
316 MG_TETRA* tet2=no2->get_lien_tetra()->get(j);
317 if (tet1==tet2)
318 {
319 MG_NOEUD *nn1,*nn2,*nn3,*nn4;
320 if (tet1->get_noeud1()==no1) nn1=tet1->get_noeud1();
321 if (tet1->get_noeud2()==no1) nn1=tet1->get_noeud2();
322 if (tet1->get_noeud3()==no1) nn1=tet1->get_noeud3();
323 if (tet1->get_noeud4()==no1) nn1=tet1->get_noeud4();
324 if (tet1->get_noeud1()==no2) nn2=tet1->get_noeud1();
325 if (tet1->get_noeud2()==no2) nn2=tet1->get_noeud2();
326 if (tet1->get_noeud3()==no2) nn2=tet1->get_noeud3();
327 if (tet1->get_noeud4()==no2) nn2=tet1->get_noeud4();
328 if (tet1->get_noeud1()!=nn1)
329 if (tet1->get_noeud1()!=nn2) nn3=tet1->get_noeud1();
330 if (tet1->get_noeud2()!=nn1)
331 if (tet1->get_noeud2()!=nn2) nn3=tet1->get_noeud2();
332 if (tet1->get_noeud3()!=nn1)
333 if (tet1->get_noeud3()!=nn2) nn3=tet1->get_noeud3();
334 if (tet1->get_noeud4()!=nn1)
335 if (tet1->get_noeud4()!=nn2) nn3=tet1->get_noeud4();
336 if (tet1->get_noeud1()!=nn1)
337 if (tet1->get_noeud1()!=nn2)
338 if (tet1->get_noeud1()!=nn3) nn4=tet1->get_noeud1();
339 if (tet1->get_noeud2()!=nn1)
340 if (tet1->get_noeud2()!=nn2)
341 if (tet1->get_noeud2()!=nn3) nn4=tet1->get_noeud2();
342 if (tet1->get_noeud3()!=nn1)
343 if (tet1->get_noeud3()!=nn2)
344 if (tet1->get_noeud3()!=nn3) nn4=tet1->get_noeud3();
345 if (tet1->get_noeud4()!=nn1)
346 if (tet1->get_noeud4()!=nn2)
347 if (tet1->get_noeud4()!=nn3) nn4=tet1->get_noeud4();
348 double angle=get_angle_par_noeud<MG_NOEUD*>(nn1,nn2,nn3,nn1,nn2,nn4);
349 if (angle>M_PI) angle=2.*M_PI-angle;
350 angletot=angletot+angle;
351
352 }
353 }
354 if ((angletot<2*M_PI-0.001) || (angletot>2*M_PI+0.001))
355 {
356 sprintf(mess," segment %lu angle matiere autour de %.2lf",seg->get_id(),angletot); affiche(mess);
357 nbsegincorrect++;
358 }
359 else nbsegcorrect++;
360 }
361 sprintf(mess," segment correct %d(%.2lf%%) segment incorrect %d(%.2lf%%) ",nbsegcorrect,nbsegcorrect*100./(nbsegcorrect+nbsegincorrect),nbsegincorrect,nbsegincorrect*100./(nbsegcorrect+nbsegincorrect)); affiche(mess);
362 }
363
364 if (nom!=NULL)
365 {
366 std::string nomf=nom;
367 nomf=nomf+".magic";
368 mai->get_gestionnaire()->enregistrer((char*)nomf.c_str());
369 }
370 }