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, 7 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

# User Rev Content
1 francois 425 #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     }