ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur_analyse.cpp
Revision: 425
Committed: Tue Sep 24 22:01:46 2013 UTC (11 years, 11 months ago) by francois
File size: 13407 byte(s)
Log Message:
ajout d'une procedure d'analyse de qualité de maillage + amelioration de la sortie sur terminal des informations dans les mailleurs

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)
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
35
36
37
38 void MAILLEUR_ANALYSE::analyse(char *nom)
39 {
40 char mess[500];
41 sprintf(mess,"ANALYSEUR DE MG MAILLAGE 3D");
42 affiche(mess);
43 sprintf(mess,"---------------------------");
44 affiche(mess);
45 sprintf(mess,"");
46 affiche(mess);
47
48 LISTE_MG_NOEUD::iterator itn;
49 for (MG_NOEUD* no=mai->get_premier_noeud(itn);no!=NULL;no=mai->get_suivant_noeud(itn))
50 no->change_nouveau_numero(0);
51
52 LISTE_MG_SEGMENT::iterator its;
53 for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
54 seg->change_nouveau_numero(0);
55
56
57 LISTE_MG_TRIANGLE::iterator ittr;
58 for (MG_TRIANGLE* tri=mai->get_premier_triangle(ittr);tri!=NULL;tri=mai->get_suivant_triangle(ittr))
59 tri->change_nouveau_numero(0);
60
61 LISTE_MG_QUADRANGLE::iterator itq;
62 for (MG_QUADRANGLE* quad=mai->get_premier_quadrangle(itq);quad!=NULL;quad=mai->get_suivant_quadrangle(itq))
63 quad->change_nouveau_numero(0);
64
65 LISTE_MG_TETRA::iterator itt;
66 for (MG_TETRA* tet=mai->get_premier_tetra(itt);tet!=NULL;tet=mai->get_suivant_tetra(itt))
67 {
68 tet->get_triangle1()->change_nouveau_numero(tet->get_triangle1()->get_nouveau_numero()+1);
69 tet->get_triangle2()->change_nouveau_numero(tet->get_triangle2()->get_nouveau_numero()+1);
70 tet->get_triangle3()->change_nouveau_numero(tet->get_triangle3()->get_nouveau_numero()+1);
71 tet->get_triangle4()->change_nouveau_numero(tet->get_triangle4()->get_nouveau_numero()+1);
72 }
73
74 LISTE_MG_HEXA::iterator ith;
75 for (MG_HEXA* hex=mai->get_premier_hexa(ith);hex!=NULL;hex=mai->get_suivant_hexa(ith))
76 {
77 hex->get_quadrangle1()->change_nouveau_numero(hex->get_quadrangle1()->get_nouveau_numero()+1);
78 hex->get_quadrangle2()->change_nouveau_numero(hex->get_quadrangle2()->get_nouveau_numero()+1);
79 hex->get_quadrangle3()->change_nouveau_numero(hex->get_quadrangle3()->get_nouveau_numero()+1);
80 hex->get_quadrangle4()->change_nouveau_numero(hex->get_quadrangle4()->get_nouveau_numero()+1);
81 hex->get_quadrangle5()->change_nouveau_numero(hex->get_quadrangle5()->get_nouveau_numero()+1);
82 hex->get_quadrangle6()->change_nouveau_numero(hex->get_quadrangle6()->get_nouveau_numero()+1);
83 }
84 int nbtrifront=0;
85 for (MG_TRIANGLE* tri=mai->get_premier_triangle(ittr);tri!=NULL;tri=mai->get_suivant_triangle(ittr))
86 {
87 if (mai->get_nb_mg_tetra()==0)
88 tri->change_nouveau_numero(1);
89 if (tri->get_nouveau_numero()==1)
90 {
91 nbtrifront++;
92 tri->get_segment1()->change_nouveau_numero(tri->get_segment1()->get_nouveau_numero()+1);
93 tri->get_segment2()->change_nouveau_numero(tri->get_segment2()->get_nouveau_numero()+1);
94 tri->get_segment3()->change_nouveau_numero(tri->get_segment3()->get_nouveau_numero()+1);
95 tri->get_noeud1()->change_nouveau_numero(1);
96 tri->get_noeud2()->change_nouveau_numero(1);
97 tri->get_noeud3()->change_nouveau_numero(1);
98
99 }
100 }
101 int nbquadfront=0;
102 for (MG_QUADRANGLE* quad=mai->get_premier_quadrangle(itq);quad!=NULL;quad=mai->get_suivant_quadrangle(itq))
103 {
104 if (mai->get_nb_mg_hexa()==0)
105 quad->change_nouveau_numero(1);
106 if (quad->get_nouveau_numero()==1)
107 {
108 nbquadfront++;
109 quad->get_segment1()->change_nouveau_numero(quad->get_segment1()->get_nouveau_numero()+1);
110 quad->get_segment2()->change_nouveau_numero(quad->get_segment2()->get_nouveau_numero()+1);
111 quad->get_segment3()->change_nouveau_numero(quad->get_segment3()->get_nouveau_numero()+1);
112 quad->get_segment4()->change_nouveau_numero(quad->get_segment4()->get_nouveau_numero()+1);
113 quad->get_noeud1()->change_nouveau_numero(1);
114 quad->get_noeud2()->change_nouveau_numero(1);
115 quad->get_noeud3()->change_nouveau_numero(1);
116 quad->get_noeud4()->change_nouveau_numero(1);
117 }
118 }
119 int nbsegfront=0;
120 for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
121 if (seg->get_nouveau_numero()==2) nbsegfront++;
122 int nbnofront=0;
123 for (MG_NOEUD* no=mai->get_premier_noeud(itn);no!=NULL;no=mai->get_suivant_noeud(itn))
124 if (no->get_nouveau_numero()==1) nbnofront++;
125 sprintf(mess," Maillage numero : %lu",mai->get_id());affiche(mess);
126 sprintf(mess," %d noeuds",mai->get_nb_mg_noeud()); affiche(mess);
127 sprintf(mess," %d segments",mai->get_nb_mg_segment()); affiche(mess);
128 sprintf(mess," %d triangles",mai->get_nb_mg_triangle()); affiche(mess);
129 sprintf(mess," %d quadrangles",mai->get_nb_mg_quadrangle()); affiche(mess);
130 sprintf(mess," %d tetras",mai->get_nb_mg_tetra()); affiche(mess);
131 sprintf(mess," %d hexas",mai->get_nb_mg_hexa()); affiche(mess); affiche("");
132 sprintf(mess," %d noeuds frontiere",nbnofront); affiche(mess);
133 sprintf(mess," %d segments frontiere",nbsegfront); affiche(mess);
134 sprintf(mess," %d triangles frontiere",nbtrifront); affiche(mess);
135 sprintf(mess," %d quadrangles frontiere",nbquadfront); affiche(mess);affiche("");
136
137 MG_SOLUTION *sol=NULL;
138 if (nom!=NULL)
139 {
140 std::string nomsol=nom;nomsol=nomsol+".sol";
141 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);
142 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);
143 mai->get_gestionnaire()->ajouter_mg_solution(sol);
144 }
145 double qualmin=1e300,qualmax=-1e300,qualmoy=0.;
146 int tab2[5]={0,0,0,0,0};
147 int itri=0;
148 for (MG_TRIANGLE* tri=mai->get_premier_triangle(ittr);tri!=NULL;tri=mai->get_suivant_triangle(ittr))
149 if (tri->get_nouveau_numero()==1)
150 {
151 double qual=OPERATEUR::qualite_triangle(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord(),tri->get_noeud3()->get_coord());
152 MG_NOEUD * no1=tri->get_noeud1();
153 MG_NOEUD * no2=tri->get_noeud1();
154 MG_NOEUD * no3=tri->get_noeud1();
155 int nb1=no1->get_lien_triangle()->get_nb();
156 int nb2=no2->get_lien_triangle()->get_nb();
157 int nb3=no3->get_lien_triangle()->get_nb();
158 OT_VECTEUR_3D vec1(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord());
159 OT_VECTEUR_3D vec2(tri->get_noeud1()->get_coord(),tri->get_noeud3()->get_coord());
160 OT_VECTEUR_3D normal=vec1&vec2;
161 normal.norme();
162 int nbretourne=0;
163 for (int i=0;i<nb1;i++)
164 for (int j=0;j<nb2;j++)
165 for (int k=0;k<nb3;k++)
166 {
167 MG_TRIANGLE* tri1=no1->get_lien_triangle()->get(i);
168 MG_TRIANGLE* tri2=no2->get_lien_triangle()->get(j);
169 MG_TRIANGLE* tri3=no3->get_lien_triangle()->get(k);
170 if ((tri1==tri2) && (tri1==tri3) && (tri1!=tri) && (tri1->get_nouveau_numero()==1))
171 {
172 OT_VECTEUR_3D vec1tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud2()->get_coord());
173 OT_VECTEUR_3D vec2tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud3()->get_coord());
174 OT_VECTEUR_3D normaltmp=vec1tmp&vec2tmp;
175 normaltmp.norme();
176 double psca=normal*normaltmp;
177 if (psca<-0.999) nbretourne++;
178 }
179 }
180 if (nbretourne>1) qual=-qual;
181 qualmoy=qualmoy+qual;
182 if (qual<qualmin) qualmin=qual;
183 if (qual>qualmax) qualmax=qual;
184 if (qual<0.) tab2[0]++;
185 else if (qual<borne1) tab2[1]++;
186 else if (qual<borne2) tab2[2]++;
187 else if (qual<borne3) tab2[3]++;
188 else tab2[4]++;
189 if (sol!=NULL)
190 if (mai->get_nb_mg_tetra()==0) sol->ecrire(qual,itri,0);
191 itri++;
192 }
193 if (nbtrifront!=0)
194 {
195 sprintf(mess," qualite moyenne des triangles de frontiere %lf",qualmoy/nbtrifront); affiche(mess);
196 sprintf(mess," qualite min des triangles de frontiere %lf",qualmin); affiche(mess);
197 sprintf(mess," qualite max des triangles de frontiere %lf",qualmax); affiche(mess);
198 sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne3,tab2[4],tab2[4]*100./nbtrifront); affiche(mess);
199 sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne2,tab2[3],tab2[3]*100./nbtrifront); affiche(mess);
200 sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne1,tab2[2],tab2[2]*100./nbtrifront); affiche(mess);
201 sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",0.,tab2[1],tab2[1]*100./nbtrifront); affiche(mess);
202 sprintf(mess," nombre de triangles de frontiere retournes : %d (%.2lf%%)",tab2[0],tab2[0]*100./nbtrifront); affiche(mess); affiche("");
203 }
204
205
206
207 qualmin=1e300,qualmax=-1e300,qualmoy=0.;
208 int tab3[5]={0,0,0,0,0};
209 int itet=0;
210 for (MG_TETRA* tet=mai->get_premier_tetra(itt);tet!=NULL;tet=mai->get_suivant_tetra(itt))
211 {
212 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());
213 OT_VECTEUR_3D vec1(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord());
214 OT_VECTEUR_3D vec2(tet->get_noeud1()->get_coord(),tet->get_noeud3()->get_coord());
215 OT_VECTEUR_3D vec3(tet->get_noeud1()->get_coord(),tet->get_noeud4()->get_coord());
216 OT_VECTEUR_3D pvec=vec1&vec2;
217 double psca=pvec*vec3;
218 if (psca<0.) qual=-qual;
219 qualmoy=qualmoy+qual;
220 if (qual<qualmin) qualmin=qual;
221 if (qual>qualmax) qualmax=qual;
222 if (qual<0.) tab3[0]++;
223 else if (qual<borne1) tab3[1]++;
224 else if (qual<borne2) tab3[2]++;
225 else if (qual<borne3) tab3[3]++;
226 else tab3[4]++;
227 if (sol!=NULL)
228 sol->ecrire(qual,itet,0);
229 itet++;
230 }
231 if (mai->get_nb_mg_tetra()!=0)
232 {
233 sprintf(mess," qualite moyenne des tetra %lf",qualmoy/mai->get_nb_mg_tetra()); affiche(mess);
234 sprintf(mess," qualite min des tetra %lf",qualmin); affiche(mess);
235 sprintf(mess," qualite max des tetra %lf",qualmax); affiche(mess);
236 sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne3,tab3[4],tab3[4]*100./mai->get_nb_mg_tetra()); affiche(mess);
237 sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne2,tab3[3],tab3[3]*100./mai->get_nb_mg_tetra()); affiche(mess);
238 sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne1,tab3[2],tab3[2]*100./mai->get_nb_mg_tetra()); affiche(mess);
239 sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",0.,tab3[1],tab3[1]*100./mai->get_nb_mg_tetra()); affiche(mess);
240 sprintf(mess," nombre de tetra retournes : %d (%.2lf%%)",tab3[0],tab3[0]*100./mai->get_nb_mg_tetra()); affiche(mess);
241 int nbsegcorrect=0;
242 int nbsegincorrect=0;
243 for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
244 if (seg->get_nouveau_numero()!=2)
245 {
246 MG_NOEUD* no1=seg->get_noeud1();
247 MG_NOEUD* no2=seg->get_noeud2();
248 int nb1=no1->get_lien_tetra()->get_nb();
249 int nb2=no2->get_lien_tetra()->get_nb();
250 double angletot=0;
251 for (int i=0;i<nb1;i++)
252 for (int j=0;j<nb2;j++)
253 {
254 MG_TETRA* tet1=no1->get_lien_tetra()->get(i);
255 MG_TETRA* tet2=no2->get_lien_tetra()->get(j);
256 if (tet1==tet2)
257 {
258 MG_NOEUD *nn1,*nn2,*nn3,*nn4;
259 if (tet1->get_noeud1()==no1) nn1=tet1->get_noeud1();
260 if (tet1->get_noeud2()==no1) nn1=tet1->get_noeud2();
261 if (tet1->get_noeud3()==no1) nn1=tet1->get_noeud3();
262 if (tet1->get_noeud4()==no1) nn1=tet1->get_noeud4();
263 if (tet1->get_noeud1()==no2) nn2=tet1->get_noeud1();
264 if (tet1->get_noeud2()==no2) nn2=tet1->get_noeud2();
265 if (tet1->get_noeud3()==no2) nn2=tet1->get_noeud3();
266 if (tet1->get_noeud4()==no2) nn2=tet1->get_noeud4();
267 if (tet1->get_noeud1()!=nn1)
268 if (tet1->get_noeud1()!=nn2) nn3=tet1->get_noeud1();
269 if (tet1->get_noeud2()!=nn1)
270 if (tet1->get_noeud2()!=nn2) nn3=tet1->get_noeud2();
271 if (tet1->get_noeud3()!=nn1)
272 if (tet1->get_noeud3()!=nn2) nn3=tet1->get_noeud3();
273 if (tet1->get_noeud4()!=nn1)
274 if (tet1->get_noeud4()!=nn2) nn3=tet1->get_noeud4();
275 if (tet1->get_noeud1()!=nn1)
276 if (tet1->get_noeud1()!=nn2)
277 if (tet1->get_noeud1()!=nn3) nn4=tet1->get_noeud1();
278 if (tet1->get_noeud2()!=nn1)
279 if (tet1->get_noeud2()!=nn2)
280 if (tet1->get_noeud2()!=nn3) nn4=tet1->get_noeud2();
281 if (tet1->get_noeud3()!=nn1)
282 if (tet1->get_noeud3()!=nn2)
283 if (tet1->get_noeud3()!=nn3) nn4=tet1->get_noeud3();
284 if (tet1->get_noeud4()!=nn1)
285 if (tet1->get_noeud4()!=nn2)
286 if (tet1->get_noeud4()!=nn3) nn4=tet1->get_noeud4();
287 double angle=get_angle_par_noeud<MG_NOEUD*>(nn1,nn2,nn3,nn1,nn2,nn4);
288 if (angle>M_PI) angle=2.*M_PI-angle;
289 angletot=angletot+angle;
290
291 }
292 }
293 if ((angletot<2*M_PI-0.001) || (angletot>2*M_PI+0.001))
294 {
295 sprintf(mess," segment %lu angle matiere autour de %.2lf",seg->get_id(),angletot); affiche(mess);
296 nbsegincorrect++;
297 }
298 else nbsegcorrect++;
299 }
300 sprintf(mess," segment correct %d(%.2lf%%) segment incorrect %d(%.2lf%%) ",nbsegcorrect,nbsegcorrect*100./(nbsegcorrect+nbsegincorrect),nbsegincorrect,nbsegincorrect*100./(nbsegcorrect+nbsegincorrect)); affiche(mess);
301 }
302
303 if (nom!=NULL)
304 {
305 std::string nomf=nom;
306 nomf=nomf+".magic";
307 mai->get_gestionnaire()->enregistrer((char*)nomf.c_str());
308 }
309 }