ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur_analyse.cpp
Revision: 426
Committed: Wed Sep 25 20:32:34 2013 UTC (11 years, 7 months ago) by francois
File size: 13980 byte(s)
Log Message:
parametrage de l'angle limite pour l'analyse des triangles retournés

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 francois 426 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 francois 425 {
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 francois 426 void MAILLEUR_ANALYSE::change_borne(double val1,double val2,double val3)
35     {
36     borne1=val1;
37     borne2=val2;
38     borne3=val3;
39     }
40 francois 425
41 francois 426 void MAILLEUR_ANALYSE::get_borne(double &val1,double &val2,double &val3)
42     {
43     val1=borne1;
44     val2=borne2;
45     val3=borne3;
46     }
47 francois 425
48 francois 426 void MAILLEUR_ANALYSE::change_eps_angle_retourne(double val)
49     {
50     eps_angle_retourne=val;
51     }
52 francois 425
53 francois 426 double MAILLEUR_ANALYSE::get_eps_angle_retourne(void)
54     {
55     return eps_angle_retourne;
56     }
57 francois 425
58 francois 426
59 francois 425 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[5]={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_noeud1();
175     MG_NOEUD * no3=tri->get_noeud1();
176     int nb1=no1->get_lien_triangle()->get_nb();
177     int nb2=no2->get_lien_triangle()->get_nb();
178     int nb3=no3->get_lien_triangle()->get_nb();
179     OT_VECTEUR_3D vec1(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord());
180     OT_VECTEUR_3D vec2(tri->get_noeud1()->get_coord(),tri->get_noeud3()->get_coord());
181     OT_VECTEUR_3D normal=vec1&vec2;
182     normal.norme();
183     int nbretourne=0;
184     for (int i=0;i<nb1;i++)
185     for (int j=0;j<nb2;j++)
186     for (int k=0;k<nb3;k++)
187     {
188     MG_TRIANGLE* tri1=no1->get_lien_triangle()->get(i);
189     MG_TRIANGLE* tri2=no2->get_lien_triangle()->get(j);
190     MG_TRIANGLE* tri3=no3->get_lien_triangle()->get(k);
191     if ((tri1==tri2) && (tri1==tri3) && (tri1!=tri) && (tri1->get_nouveau_numero()==1))
192     {
193     OT_VECTEUR_3D vec1tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud2()->get_coord());
194     OT_VECTEUR_3D vec2tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud3()->get_coord());
195     OT_VECTEUR_3D normaltmp=vec1tmp&vec2tmp;
196     normaltmp.norme();
197     double psca=normal*normaltmp;
198 francois 426 if (psca<-cos(eps_angle_retourne)) nbretourne++;
199 francois 425 }
200     }
201     if (nbretourne>1) qual=-qual;
202     qualmoy=qualmoy+qual;
203     if (qual<qualmin) qualmin=qual;
204     if (qual>qualmax) qualmax=qual;
205     if (qual<0.) tab2[0]++;
206     else if (qual<borne1) tab2[1]++;
207     else if (qual<borne2) tab2[2]++;
208     else if (qual<borne3) tab2[3]++;
209     else tab2[4]++;
210     if (sol!=NULL)
211     if (mai->get_nb_mg_tetra()==0) sol->ecrire(qual,itri,0);
212     itri++;
213     }
214     if (nbtrifront!=0)
215     {
216 francois 426 sprintf(mess," critere d'angle limite pour l'inversion : %lf",eps_angle_retourne); affiche(mess);
217     sprintf(mess," qualite moyenne des triangles de frontiere : %lf",qualmoy/nbtrifront); affiche(mess);
218     sprintf(mess," qualite min des triangles de frontiere : %lf",qualmin); affiche(mess);
219     sprintf(mess," qualite max des triangles de frontiere : %lf",qualmax); affiche(mess);
220 francois 425 sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne3,tab2[4],tab2[4]*100./nbtrifront); affiche(mess);
221     sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne2,tab2[3],tab2[3]*100./nbtrifront); affiche(mess);
222     sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne1,tab2[2],tab2[2]*100./nbtrifront); affiche(mess);
223     sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",0.,tab2[1],tab2[1]*100./nbtrifront); affiche(mess);
224     sprintf(mess," nombre de triangles de frontiere retournes : %d (%.2lf%%)",tab2[0],tab2[0]*100./nbtrifront); affiche(mess); affiche("");
225     }
226    
227    
228    
229     qualmin=1e300,qualmax=-1e300,qualmoy=0.;
230     int tab3[5]={0,0,0,0,0};
231     int itet=0;
232     for (MG_TETRA* tet=mai->get_premier_tetra(itt);tet!=NULL;tet=mai->get_suivant_tetra(itt))
233     {
234     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());
235     OT_VECTEUR_3D vec1(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord());
236     OT_VECTEUR_3D vec2(tet->get_noeud1()->get_coord(),tet->get_noeud3()->get_coord());
237     OT_VECTEUR_3D vec3(tet->get_noeud1()->get_coord(),tet->get_noeud4()->get_coord());
238     OT_VECTEUR_3D pvec=vec1&vec2;
239     double psca=pvec*vec3;
240     if (psca<0.) qual=-qual;
241     qualmoy=qualmoy+qual;
242     if (qual<qualmin) qualmin=qual;
243     if (qual>qualmax) qualmax=qual;
244     if (qual<0.) tab3[0]++;
245     else if (qual<borne1) tab3[1]++;
246     else if (qual<borne2) tab3[2]++;
247     else if (qual<borne3) tab3[3]++;
248     else tab3[4]++;
249     if (sol!=NULL)
250     sol->ecrire(qual,itet,0);
251     itet++;
252     }
253     if (mai->get_nb_mg_tetra()!=0)
254     {
255 francois 426 sprintf(mess," qualite moyenne des tetra : %lf",qualmoy/mai->get_nb_mg_tetra()); affiche(mess);
256     sprintf(mess," qualite min des tetra : %lf",qualmin); affiche(mess);
257     sprintf(mess," qualite max des tetra : %lf",qualmax); affiche(mess);
258 francois 425 sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne3,tab3[4],tab3[4]*100./mai->get_nb_mg_tetra()); affiche(mess);
259     sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne2,tab3[3],tab3[3]*100./mai->get_nb_mg_tetra()); affiche(mess);
260     sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne1,tab3[2],tab3[2]*100./mai->get_nb_mg_tetra()); affiche(mess);
261     sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",0.,tab3[1],tab3[1]*100./mai->get_nb_mg_tetra()); affiche(mess);
262     sprintf(mess," nombre de tetra retournes : %d (%.2lf%%)",tab3[0],tab3[0]*100./mai->get_nb_mg_tetra()); affiche(mess);
263     int nbsegcorrect=0;
264     int nbsegincorrect=0;
265     for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
266     if (seg->get_nouveau_numero()!=2)
267     {
268     MG_NOEUD* no1=seg->get_noeud1();
269     MG_NOEUD* no2=seg->get_noeud2();
270     int nb1=no1->get_lien_tetra()->get_nb();
271     int nb2=no2->get_lien_tetra()->get_nb();
272     double angletot=0;
273     for (int i=0;i<nb1;i++)
274     for (int j=0;j<nb2;j++)
275     {
276     MG_TETRA* tet1=no1->get_lien_tetra()->get(i);
277     MG_TETRA* tet2=no2->get_lien_tetra()->get(j);
278     if (tet1==tet2)
279     {
280     MG_NOEUD *nn1,*nn2,*nn3,*nn4;
281     if (tet1->get_noeud1()==no1) nn1=tet1->get_noeud1();
282     if (tet1->get_noeud2()==no1) nn1=tet1->get_noeud2();
283     if (tet1->get_noeud3()==no1) nn1=tet1->get_noeud3();
284     if (tet1->get_noeud4()==no1) nn1=tet1->get_noeud4();
285     if (tet1->get_noeud1()==no2) nn2=tet1->get_noeud1();
286     if (tet1->get_noeud2()==no2) nn2=tet1->get_noeud2();
287     if (tet1->get_noeud3()==no2) nn2=tet1->get_noeud3();
288     if (tet1->get_noeud4()==no2) nn2=tet1->get_noeud4();
289     if (tet1->get_noeud1()!=nn1)
290     if (tet1->get_noeud1()!=nn2) nn3=tet1->get_noeud1();
291     if (tet1->get_noeud2()!=nn1)
292     if (tet1->get_noeud2()!=nn2) nn3=tet1->get_noeud2();
293     if (tet1->get_noeud3()!=nn1)
294     if (tet1->get_noeud3()!=nn2) nn3=tet1->get_noeud3();
295     if (tet1->get_noeud4()!=nn1)
296     if (tet1->get_noeud4()!=nn2) nn3=tet1->get_noeud4();
297     if (tet1->get_noeud1()!=nn1)
298     if (tet1->get_noeud1()!=nn2)
299     if (tet1->get_noeud1()!=nn3) nn4=tet1->get_noeud1();
300     if (tet1->get_noeud2()!=nn1)
301     if (tet1->get_noeud2()!=nn2)
302     if (tet1->get_noeud2()!=nn3) nn4=tet1->get_noeud2();
303     if (tet1->get_noeud3()!=nn1)
304     if (tet1->get_noeud3()!=nn2)
305     if (tet1->get_noeud3()!=nn3) nn4=tet1->get_noeud3();
306     if (tet1->get_noeud4()!=nn1)
307     if (tet1->get_noeud4()!=nn2)
308     if (tet1->get_noeud4()!=nn3) nn4=tet1->get_noeud4();
309     double angle=get_angle_par_noeud<MG_NOEUD*>(nn1,nn2,nn3,nn1,nn2,nn4);
310     if (angle>M_PI) angle=2.*M_PI-angle;
311     angletot=angletot+angle;
312    
313     }
314     }
315     if ((angletot<2*M_PI-0.001) || (angletot>2*M_PI+0.001))
316     {
317     sprintf(mess," segment %lu angle matiere autour de %.2lf",seg->get_id(),angletot); affiche(mess);
318     nbsegincorrect++;
319     }
320     else nbsegcorrect++;
321     }
322     sprintf(mess," segment correct %d(%.2lf%%) segment incorrect %d(%.2lf%%) ",nbsegcorrect,nbsegcorrect*100./(nbsegcorrect+nbsegincorrect),nbsegincorrect,nbsegincorrect*100./(nbsegcorrect+nbsegincorrect)); affiche(mess);
323     }
324    
325     if (nom!=NULL)
326     {
327     std::string nomf=nom;
328     nomf=nomf+".magic";
329     mai->get_gestionnaire()->enregistrer((char*)nomf.c_str());
330     }
331     }