ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur_analyse.cpp
Revision: 828
Committed: Thu Sep 15 21:18:11 2016 UTC (8 years, 7 months ago) by francois
Original Path: magic/lib/mailleur_auto/src/mailleur_analyse.cpp
File size: 23678 byte(s)
Log Message:
Amekioration de la quadratisation du maillage. Meme la ou la parametrisation est disordu le maillage quadratique ne l'est pas.
Analyseur de qualité de maillage FEM

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 francois 514 #include "fct_taille.h"
7 francois 828 #include "fem_maillage_quadratique_outils.h"
8 francois 425
9    
10    
11    
12    
13 francois 828 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),fem(NULL)
14     {
15    
16     }
17 francois 425
18 francois 828
19     MAILLEUR_ANALYSE::MAILLEUR_ANALYSE(FEM_MAILLAGE* m,OT_CPU* comp):MAILLEUR(comp),fem(m),mai(NULL)
20 francois 425 {
21    
22     }
23    
24    
25    
26    
27    
28     MAILLEUR_ANALYSE::MAILLEUR_ANALYSE(MAILLEUR_ANALYSE &mdd):MAILLEUR(mdd),mai(mdd.mai),borne1(mdd.borne1),borne2(mdd.borne2),borne3(mdd.borne3)
29     {
30    
31     }
32    
33    
34    
35     MAILLEUR_ANALYSE::~MAILLEUR_ANALYSE()
36     {
37    
38     }
39    
40 francois 426 void MAILLEUR_ANALYSE::change_borne(double val1,double val2,double val3)
41     {
42     borne1=val1;
43     borne2=val2;
44     borne3=val3;
45     }
46 francois 425
47 francois 426 void MAILLEUR_ANALYSE::get_borne(double &val1,double &val2,double &val3)
48     {
49     val1=borne1;
50     val2=borne2;
51     val3=borne3;
52     }
53 francois 425
54 francois 426 void MAILLEUR_ANALYSE::change_eps_angle_retourne(double val)
55     {
56     eps_angle_retourne=val;
57     }
58 francois 425
59 francois 426 double MAILLEUR_ANALYSE::get_eps_angle_retourne(void)
60     {
61     return eps_angle_retourne;
62     }
63 francois 425
64 francois 426
65 francois 425 void MAILLEUR_ANALYSE::analyse(char *nom)
66     {
67 francois 828 if (mai!=NULL) analyse_mg(nom);
68     if (fem!=NULL) analyse_fem(nom);
69     }
70    
71     void MAILLEUR_ANALYSE::analyse_fem(char *nom)
72     {
73 francois 425 char mess[500];
74 francois 828 sprintf(mess,"ANALYSEUR DE FEM MAILLAGE 3D");
75     affiche(mess);
76     sprintf(mess,"---------------------------");
77     affiche(mess);
78     affiche((char*)"");
79     sprintf(mess,"Constitution du maillage");
80     affiche(mess);
81     affiche((char*)"");
82     int nbnoeud=fem->get_nb_fem_noeud();
83     int nbele1=fem->get_nb_fem_element1();
84     int nbele2=fem->get_nb_fem_element2();
85     int nbele3=fem->get_nb_fem_element3();
86     int degre=fem->get_degre();
87     sprintf(mess," Maillage numero : %lu degre : %d",fem->get_id(),degre);affiche(mess);
88     sprintf(mess," %d noeuds",nbnoeud); affiche(mess);
89     sprintf(mess," %d elements 1D",nbele1); affiche(mess);
90    
91     LISTE_FEM_ELEMENT1::iterator it1;
92     int nbseg2=0,nbseg3=0;
93     for (FEM_ELEMENT1 *ele=fem->get_premier_element1(it1);ele!=NULL;ele=fem->get_suivant_element1(it1))
94     {
95     if (ele->get_nb_fem_noeud()==2) nbseg2++;
96     if (ele->get_nb_fem_noeud()==3) nbseg3++;
97     }
98     if (degre==1) {sprintf(mess," %d segment2",nbseg2); affiche(mess); }
99     if (degre==2) {sprintf(mess," %d segment3",nbseg3); affiche(mess); }
100     sprintf(mess," %d elements 2D",nbele2); affiche(mess);
101     LISTE_FEM_ELEMENT2::iterator it2;
102     int nbtri3=0,nbtri6=0,nbquad4=0,nbquad8=0;
103     for (FEM_ELEMENT2 *ele=fem->get_premier_element2(it2);ele!=NULL;ele=fem->get_suivant_element2(it2))
104     {
105     if (ele->get_nb_fem_noeud()==3) nbtri3++;
106     if (ele->get_nb_fem_noeud()==6) nbtri6++;
107     if (ele->get_nb_fem_noeud()==4) nbquad4++;
108     if (ele->get_nb_fem_noeud()==8) nbquad8++;
109     }
110     if (degre==1) {sprintf(mess," %d triangle3",nbtri3); affiche(mess); }
111     if (degre==2) {sprintf(mess," %d triangle6",nbtri6); affiche(mess); }
112     if (degre==1) {sprintf(mess," %d quadrangle4",nbquad4); affiche(mess); }
113     if (degre==2) {sprintf(mess," %d quadrangle8",nbquad8); affiche(mess); }
114     sprintf(mess," %d elements 3D",nbele3); affiche(mess);
115     LISTE_FEM_ELEMENT3::iterator it3;
116     int nbtet4=0,nbtet10=0,nbhex8=0,nbhex20=0;
117     for (FEM_ELEMENT3 *ele=fem->get_premier_element3(it3);ele!=NULL;ele=fem->get_suivant_element3(it3))
118     {
119     if (ele->get_nb_fem_noeud()==4) nbtet4++;
120     if (ele->get_nb_fem_noeud()==10) nbtet10++;
121     if (ele->get_nb_fem_noeud()==8) nbhex8++;
122     if (ele->get_nb_fem_noeud()==20) nbhex20++;
123     }
124     if (degre==1) {sprintf(mess," %d tetra4",nbtet4); affiche(mess); }
125     if (degre==2) {sprintf(mess," %d tetra10",nbtet10); affiche(mess); }
126     if (degre==1) {sprintf(mess," %d hex8",nbhex8); affiche(mess); }
127     if (degre==2) {sprintf(mess," %d hex20",nbhex20); affiche(mess); }
128     affiche((char*)"");
129     sprintf(mess,"Qualité du maillage");
130     affiche(mess);
131     if (degre==1)
132     {
133     mai=fem->get_mg_maillage();
134     analyse_mg(nom);
135     }
136     if (degre==2)
137     {
138     if (nbtri6>0)
139     {
140     sprintf(mess," Qualité des triangle6");
141     affiche(mess);
142     FEM_MAILLAGE_QUADRATIQUE_OUTILS ot;
143     double jminmin=1e300,jmaxmax=1e-300,dismoy=0,dismin=1e300,dismax=-1e300;
144     int nb=0;
145     for (FEM_ELEMENT2 *ele=fem->get_premier_element2(it2);ele!=NULL;ele=fem->get_suivant_element2(it2))
146     {
147     if (ele->get_nb_fem_noeud()==6)
148     {
149     double jmin=ot.get_jmin(ele);
150     double jmax=ot.get_jmax(ele);
151     double dis=ot.get_distorsion2(ele);
152     nb++;
153     dismoy=dismoy+dis;
154     if (jmin<jminmin) jminmin=jmin;
155     if (jmax>jmaxmax) jmaxmax=jmax;
156     if (dis<dismin) dismin=dis;
157     if (dis>dismax) dismax=dis;
158     }
159     }
160     sprintf(mess," Jacobien min %lf",jminmin);
161     affiche(mess);
162     sprintf(mess," Jacobien max %lf",jmaxmax);
163     affiche(mess);
164     sprintf(mess," Distortion min %lf",dismin);
165     affiche(mess);
166     sprintf(mess," Distortion moy %lf",dismoy/nb);
167     affiche(mess);
168     sprintf(mess," Distortion max %lf",dismax);
169     affiche(mess);
170    
171     }
172     if (nbtet10>0)
173     {
174     sprintf(mess," Qualité des tetra10");
175     affiche(mess);
176     FEM_MAILLAGE_QUADRATIQUE_OUTILS ot;
177     double jminmin=1e300,jmaxmax=1e-300,dismoy=0,dismin=1e300,dismax=-1e300;
178     int nb=0;
179     for (FEM_ELEMENT3 *ele=fem->get_premier_element3(it3);ele!=NULL;ele=fem->get_suivant_element3(it3))
180     {
181     if (ele->get_nb_fem_noeud()==10)
182     {
183     double jmin=ot.get_jmin(ele);
184     double jmax=ot.get_jmax(ele);
185     double dis=ot.get_distorsion2(ele);
186     nb++;
187     dismoy=dismoy+dis;
188     if (jmin<jminmin) jminmin=jmin;
189     if (jmax>jmaxmax) jmaxmax=jmax;
190     if (dis<dismin) dismin=dis;
191     if (dis>dismax) dismax=dis;
192     }
193     }
194     sprintf(mess," Jacobien min %lf",jminmin);
195     affiche(mess);
196     sprintf(mess," Jacobien max %lf",jmaxmax);
197     affiche(mess);
198     sprintf(mess," Distortion min %lf",dismin);
199     affiche(mess);
200     sprintf(mess," Distortion moy %lf",dismoy/nb);
201     affiche(mess);
202     sprintf(mess," Distortion max %lf",dismax);
203     affiche(mess);
204    
205     }
206     }
207     }
208    
209    
210     void MAILLEUR_ANALYSE::analyse_mg(char *nom)
211     {
212     char mess[500];
213 francois 425 sprintf(mess,"ANALYSEUR DE MG MAILLAGE 3D");
214     affiche(mess);
215     sprintf(mess,"---------------------------");
216     affiche(mess);
217 francois 660 affiche((char*)"");
218 francois 828 sprintf(mess,"Constitution du maillage");
219 francois 805 affiche(mess);
220     affiche((char*)"");
221     int nbnofront,nbsegfront,nbtrifront,nbquadfront;
222     denombre_maillage(nbnofront,nbsegfront,nbtrifront,nbquadfront);
223     sprintf(mess," Maillage numero : %lu",mai->get_id());affiche(mess);
224     sprintf(mess," %d noeuds",mai->get_nb_mg_noeud()); affiche(mess);
225     sprintf(mess," %d segments",mai->get_nb_mg_segment()); affiche(mess);
226     sprintf(mess," %d triangles",mai->get_nb_mg_triangle()); affiche(mess);
227     sprintf(mess," %d quadrangles",mai->get_nb_mg_quadrangle()); affiche(mess);
228     sprintf(mess," %d tetras",mai->get_nb_mg_tetra()); affiche(mess);
229     sprintf(mess," %d hexas",mai->get_nb_mg_hexa()); affiche(mess); affiche((char*)"");
230     sprintf(mess," %d noeuds frontiere",nbnofront); affiche(mess);
231     sprintf(mess," %d segments frontiere",nbsegfront); affiche(mess);
232     sprintf(mess," %d triangles frontiere",nbtrifront); affiche(mess);
233     sprintf(mess," %d quadrangles frontiere",nbquadfront); affiche(mess);affiche((char*)"");
234 francois 425
235 francois 805
236     MG_SOLUTION *sol=NULL;
237     if (nom!=NULL)
238     {
239     std::string nomsol=nom;nomsol=nomsol+".sol";
240     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);
241     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);
242     mai->get_gestionnaire()->ajouter_mg_solution(sol);
243     }
244    
245    
246     if (nbtrifront!=0)
247     {
248     affiche((char*)"");
249     sprintf(mess,"Analyse maillage triangulaire 2D");
250     affiche(mess);
251     affiche((char*)"");
252     double qualmin,qualmax,qualmoy;
253     int tab[7];
254     analyse_qualite_maillage_2D(sol,qualmin,qualmax,qualmoy,tab);
255     sprintf(mess," Nombre de triangles avec des noeuds fusionnés : %d",tab[6]); affiche(mess);
256     sprintf(mess," Nombre de triangles avec 1 voisin manquant : %d",tab[5]); affiche(mess);
257     sprintf(mess," critere d'angle limite pour l'inversion : %lf",eps_angle_retourne); affiche(mess);
258     sprintf(mess," qualite moyenne des triangles de frontiere : %lf",qualmoy); affiche(mess);
259     sprintf(mess," qualite min des triangles de frontiere : %lf",qualmin); affiche(mess);
260     sprintf(mess," qualite max des triangles de frontiere : %lf",qualmax); affiche(mess);
261     sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne3,tab[4],tab[4]*100./nbtrifront); affiche(mess);
262     sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne2,tab[3],tab[3]*100./nbtrifront); affiche(mess);
263     sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne1,tab[2],tab[2]*100./nbtrifront); affiche(mess);
264     sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",0.,tab[1],tab[1]*100./nbtrifront); affiche(mess);
265     sprintf(mess," nombre de triangles de frontiere retournes : %d (%.2lf%%)",tab[0],tab[0]*100./nbtrifront); affiche(mess); affiche((char*)"");
266     }
267     if (mai->get_nb_mg_tetra()!=0)
268     {
269     affiche((char*)"");
270     sprintf(mess,"Analyse maillage triangulaire 3D");
271     affiche(mess);
272     affiche((char*)"");
273     double qualmin,qualmax,qualmoy;
274     int tab[5];
275     analyse_qualite_maillage_3D(sol,qualmin,qualmax,qualmoy,tab);
276     sprintf(mess," qualite moyenne des tetra : %lf",qualmoy); affiche(mess);
277     sprintf(mess," qualite min des tetra : %lf",qualmin); affiche(mess);
278     sprintf(mess," qualite max des tetra : %lf",qualmax); affiche(mess);
279     sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne3,tab[4],tab[4]*100./mai->get_nb_mg_tetra()); affiche(mess);
280     sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne2,tab[3],tab[3]*100./mai->get_nb_mg_tetra()); affiche(mess);
281     sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne1,tab[2],tab[2]*100./mai->get_nb_mg_tetra()); affiche(mess);
282     sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",0.,tab[1],tab[1]*100./mai->get_nb_mg_tetra()); affiche(mess);
283     sprintf(mess," nombre de tetra retournes : %d (%.2lf%%)",tab[0],tab[0]*100./mai->get_nb_mg_tetra()); affiche(mess);
284     affiche((char*)"");
285     sprintf(mess,"Validite maillage triangulaire 3D");
286     affiche(mess);
287     affiche((char*)"");
288     int nbsegcorrect,nbsegincorrect;
289     analyse_validite_maillage_3D(nbsegcorrect,nbsegincorrect);
290    
291     sprintf(mess," segment correct %d(%.2lf%%) segment incorrect %d(%.2lf%%) ",nbsegcorrect,nbsegcorrect*100./(nbsegcorrect+nbsegincorrect),nbsegincorrect,nbsegincorrect*100./(nbsegcorrect+nbsegincorrect)); affiche(mess);
292    
293     }
294    
295     if (nom!=NULL)
296     {
297     std::string nomf=nom;
298     nomf=nomf+".magic";
299     mai->get_gestionnaire()->enregistrer((char*)nomf.c_str());
300     }
301    
302     }
303    
304     void MAILLEUR_ANALYSE::denombre_maillage(int &nbnofront,int &nbsegfront,int &nbtrifront,int &nbquadfront)
305     {
306 francois 425 LISTE_MG_NOEUD::iterator itn;
307     for (MG_NOEUD* no=mai->get_premier_noeud(itn);no!=NULL;no=mai->get_suivant_noeud(itn))
308     no->change_nouveau_numero(0);
309    
310     LISTE_MG_SEGMENT::iterator its;
311     for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
312     seg->change_nouveau_numero(0);
313    
314    
315     LISTE_MG_TRIANGLE::iterator ittr;
316     for (MG_TRIANGLE* tri=mai->get_premier_triangle(ittr);tri!=NULL;tri=mai->get_suivant_triangle(ittr))
317     tri->change_nouveau_numero(0);
318    
319     LISTE_MG_QUADRANGLE::iterator itq;
320     for (MG_QUADRANGLE* quad=mai->get_premier_quadrangle(itq);quad!=NULL;quad=mai->get_suivant_quadrangle(itq))
321     quad->change_nouveau_numero(0);
322    
323     LISTE_MG_TETRA::iterator itt;
324     for (MG_TETRA* tet=mai->get_premier_tetra(itt);tet!=NULL;tet=mai->get_suivant_tetra(itt))
325     {
326     tet->get_triangle1()->change_nouveau_numero(tet->get_triangle1()->get_nouveau_numero()+1);
327     tet->get_triangle2()->change_nouveau_numero(tet->get_triangle2()->get_nouveau_numero()+1);
328     tet->get_triangle3()->change_nouveau_numero(tet->get_triangle3()->get_nouveau_numero()+1);
329     tet->get_triangle4()->change_nouveau_numero(tet->get_triangle4()->get_nouveau_numero()+1);
330     }
331    
332     LISTE_MG_HEXA::iterator ith;
333     for (MG_HEXA* hex=mai->get_premier_hexa(ith);hex!=NULL;hex=mai->get_suivant_hexa(ith))
334     {
335     hex->get_quadrangle1()->change_nouveau_numero(hex->get_quadrangle1()->get_nouveau_numero()+1);
336     hex->get_quadrangle2()->change_nouveau_numero(hex->get_quadrangle2()->get_nouveau_numero()+1);
337     hex->get_quadrangle3()->change_nouveau_numero(hex->get_quadrangle3()->get_nouveau_numero()+1);
338     hex->get_quadrangle4()->change_nouveau_numero(hex->get_quadrangle4()->get_nouveau_numero()+1);
339     hex->get_quadrangle5()->change_nouveau_numero(hex->get_quadrangle5()->get_nouveau_numero()+1);
340     hex->get_quadrangle6()->change_nouveau_numero(hex->get_quadrangle6()->get_nouveau_numero()+1);
341     }
342 francois 805 nbtrifront=0;
343 francois 425 for (MG_TRIANGLE* tri=mai->get_premier_triangle(ittr);tri!=NULL;tri=mai->get_suivant_triangle(ittr))
344     {
345     if (mai->get_nb_mg_tetra()==0)
346     tri->change_nouveau_numero(1);
347     if (tri->get_nouveau_numero()==1)
348     {
349     nbtrifront++;
350     tri->get_segment1()->change_nouveau_numero(tri->get_segment1()->get_nouveau_numero()+1);
351     tri->get_segment2()->change_nouveau_numero(tri->get_segment2()->get_nouveau_numero()+1);
352     tri->get_segment3()->change_nouveau_numero(tri->get_segment3()->get_nouveau_numero()+1);
353     tri->get_noeud1()->change_nouveau_numero(1);
354     tri->get_noeud2()->change_nouveau_numero(1);
355     tri->get_noeud3()->change_nouveau_numero(1);
356    
357     }
358     }
359 francois 805 nbquadfront=0;
360 francois 425 for (MG_QUADRANGLE* quad=mai->get_premier_quadrangle(itq);quad!=NULL;quad=mai->get_suivant_quadrangle(itq))
361     {
362     if (mai->get_nb_mg_hexa()==0)
363     quad->change_nouveau_numero(1);
364     if (quad->get_nouveau_numero()==1)
365     {
366     nbquadfront++;
367     quad->get_segment1()->change_nouveau_numero(quad->get_segment1()->get_nouveau_numero()+1);
368     quad->get_segment2()->change_nouveau_numero(quad->get_segment2()->get_nouveau_numero()+1);
369     quad->get_segment3()->change_nouveau_numero(quad->get_segment3()->get_nouveau_numero()+1);
370     quad->get_segment4()->change_nouveau_numero(quad->get_segment4()->get_nouveau_numero()+1);
371     quad->get_noeud1()->change_nouveau_numero(1);
372     quad->get_noeud2()->change_nouveau_numero(1);
373     quad->get_noeud3()->change_nouveau_numero(1);
374     quad->get_noeud4()->change_nouveau_numero(1);
375     }
376     }
377 francois 805 nbsegfront=0;
378 francois 425 for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
379     if (seg->get_nouveau_numero()==2) nbsegfront++;
380 francois 805 nbnofront=0;
381 francois 425 for (MG_NOEUD* no=mai->get_premier_noeud(itn);no!=NULL;no=mai->get_suivant_noeud(itn))
382     if (no->get_nouveau_numero()==1) nbnofront++;
383 francois 805 }
384 francois 425
385 francois 805 void MAILLEUR_ANALYSE::analyse_qualite_maillage_2D(MG_SOLUTION *sol,double &qualmin,double &qualmax,double &qualmoy,int *tab)
386     {
387     int nbnofront,nbsegfront,nbtrifront,nbquadfront;
388     denombre_maillage(nbnofront,nbsegfront,nbtrifront,nbquadfront);
389     qualmin=1e300,qualmax=-1e300,qualmoy=0.;
390     tab[0]=0;tab[1]=0;tab[2]=0;tab[3]=0;tab[4]=0;tab[5]=0;tab[6]=0;
391 francois 425 int itri=0;
392 francois 805 LISTE_MG_TRIANGLE::iterator ittr;
393 francois 425 for (MG_TRIANGLE* tri=mai->get_premier_triangle(ittr);tri!=NULL;tri=mai->get_suivant_triangle(ittr))
394     if (tri->get_nouveau_numero()==1)
395     {
396     double qual=OPERATEUR::qualite_triangle(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord(),tri->get_noeud3()->get_coord());
397     MG_NOEUD * no1=tri->get_noeud1();
398 francois 429 MG_NOEUD * no2=tri->get_noeud2();
399     MG_NOEUD * no3=tri->get_noeud3();
400     if (no1==no2)
401     if (no1==no3)
402 francois 805 tab[6]++;
403 francois 425 int nb1=no1->get_lien_triangle()->get_nb();
404     int nb2=no2->get_lien_triangle()->get_nb();
405     int nb3=no3->get_lien_triangle()->get_nb();
406     OT_VECTEUR_3D vec1(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord());
407     OT_VECTEUR_3D vec2(tri->get_noeud1()->get_coord(),tri->get_noeud3()->get_coord());
408     OT_VECTEUR_3D normal=vec1&vec2;
409     normal.norme();
410     int nbretourne=0;
411 francois 429 int nbvoisin=0;
412 francois 425 for (int i=0;i<nb1;i++)
413     for (int j=0;j<nb2;j++)
414     {
415     MG_TRIANGLE* tri1=no1->get_lien_triangle()->get(i);
416     MG_TRIANGLE* tri2=no2->get_lien_triangle()->get(j);
417 francois 429 if ( (tri1==tri2) && (tri1!=tri) && (tri1->get_nouveau_numero()==1))
418 francois 425 {
419     OT_VECTEUR_3D vec1tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud2()->get_coord());
420     OT_VECTEUR_3D vec2tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud3()->get_coord());
421     OT_VECTEUR_3D normaltmp=vec1tmp&vec2tmp;
422     normaltmp.norme();
423     double psca=normal*normaltmp;
424 francois 426 if (psca<-cos(eps_angle_retourne)) nbretourne++;
425 francois 429 nbvoisin++;
426 francois 425 }
427     }
428 francois 429 for (int i=0;i<nb1;i++)
429     for (int j=0;j<nb3;j++)
430     {
431     MG_TRIANGLE* tri1=no1->get_lien_triangle()->get(i);
432     MG_TRIANGLE* tri2=no3->get_lien_triangle()->get(j);
433     if ( (tri1==tri2) && (tri1!=tri) && (tri1->get_nouveau_numero()==1))
434     {
435     OT_VECTEUR_3D vec1tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud2()->get_coord());
436     OT_VECTEUR_3D vec2tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud3()->get_coord());
437     OT_VECTEUR_3D normaltmp=vec1tmp&vec2tmp;
438     normaltmp.norme();
439     double psca=normal*normaltmp;
440     if (psca<-cos(eps_angle_retourne)) nbretourne++;
441     nbvoisin++;
442     }
443     }
444     for (int i=0;i<nb2;i++)
445     for (int j=0;j<nb3;j++)
446     {
447     MG_TRIANGLE* tri1=no2->get_lien_triangle()->get(i);
448     MG_TRIANGLE* tri2=no3->get_lien_triangle()->get(j);
449     if ( (tri1==tri2) && (tri1!=tri) && (tri1->get_nouveau_numero()==1))
450     {
451     OT_VECTEUR_3D vec1tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud2()->get_coord());
452     OT_VECTEUR_3D vec2tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud3()->get_coord());
453     OT_VECTEUR_3D normaltmp=vec1tmp&vec2tmp;
454     normaltmp.norme();
455     double psca=normal*normaltmp;
456     if (psca<-cos(eps_angle_retourne)) nbretourne++;
457     nbvoisin++;
458     }
459     }
460     if (nbvoisin!=3)
461 francois 805 tab[5]++;
462 francois 425 if (nbretourne>1) qual=-qual;
463     qualmoy=qualmoy+qual;
464     if (qual<qualmin) qualmin=qual;
465     if (qual>qualmax) qualmax=qual;
466 francois 805 if (qual<0.) tab[0]++;
467     else if (qual<borne1) tab[1]++;
468     else if (qual<borne2) tab[2]++;
469     else if (qual<borne3) tab[3]++;
470     else tab[4]++;
471 francois 425 if (sol!=NULL)
472     if (mai->get_nb_mg_tetra()==0) sol->ecrire(qual,itri,0);
473     itri++;
474     }
475 francois 805 qualmoy=qualmoy/nbtrifront;
476     }
477    
478     void MAILLEUR_ANALYSE::analyse_qualite_maillage_3D(MG_SOLUTION *sol,double &qualmin,double &qualmax,double &qualmoy,int *tab)
479 francois 425 {
480     qualmin=1e300,qualmax=-1e300,qualmoy=0.;
481 francois 805 tab[0]=0;tab[1]=0;tab[2]=0;tab[3]=0;tab[4]=0;
482     qualmin=1e300,qualmax=-1e300,qualmoy=0.;
483 francois 425 int itet=0;
484 francois 805 LISTE_MG_TETRA::iterator itt;
485 francois 425 for (MG_TETRA* tet=mai->get_premier_tetra(itt);tet!=NULL;tet=mai->get_suivant_tetra(itt))
486     {
487     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());
488     OT_VECTEUR_3D vec1(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord());
489     OT_VECTEUR_3D vec2(tet->get_noeud1()->get_coord(),tet->get_noeud3()->get_coord());
490     OT_VECTEUR_3D vec3(tet->get_noeud1()->get_coord(),tet->get_noeud4()->get_coord());
491     OT_VECTEUR_3D pvec=vec1&vec2;
492     double psca=pvec*vec3;
493     if (psca<0.) qual=-qual;
494     qualmoy=qualmoy+qual;
495     if (qual<qualmin) qualmin=qual;
496     if (qual>qualmax) qualmax=qual;
497 francois 805 if (qual<0.) tab[0]++;
498     else if (qual<borne1) tab[1]++;
499     else if (qual<borne2) tab[2]++;
500     else if (qual<borne3) tab[3]++;
501     else tab[4]++;
502 francois 425 if (sol!=NULL)
503     sol->ecrire(qual,itet,0);
504     itet++;
505     }
506 francois 805 qualmoy=qualmoy/mai->get_nb_mg_tetra();
507     }
508    
509     void MAILLEUR_ANALYSE::analyse_validite_maillage_3D(int &nbsegcorrect,int &nbsegincorrect)
510     {
511    
512     nbsegcorrect=0;
513     nbsegincorrect=0;
514     LISTE_MG_SEGMENT::iterator its;
515 francois 425 for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
516     if (seg->get_nouveau_numero()!=2)
517     {
518     MG_NOEUD* no1=seg->get_noeud1();
519     MG_NOEUD* no2=seg->get_noeud2();
520     int nb1=no1->get_lien_tetra()->get_nb();
521     int nb2=no2->get_lien_tetra()->get_nb();
522     double angletot=0;
523     for (int i=0;i<nb1;i++)
524     for (int j=0;j<nb2;j++)
525     {
526     MG_TETRA* tet1=no1->get_lien_tetra()->get(i);
527     MG_TETRA* tet2=no2->get_lien_tetra()->get(j);
528     if (tet1==tet2)
529     {
530     MG_NOEUD *nn1,*nn2,*nn3,*nn4;
531     if (tet1->get_noeud1()==no1) nn1=tet1->get_noeud1();
532     if (tet1->get_noeud2()==no1) nn1=tet1->get_noeud2();
533     if (tet1->get_noeud3()==no1) nn1=tet1->get_noeud3();
534     if (tet1->get_noeud4()==no1) nn1=tet1->get_noeud4();
535     if (tet1->get_noeud1()==no2) nn2=tet1->get_noeud1();
536     if (tet1->get_noeud2()==no2) nn2=tet1->get_noeud2();
537     if (tet1->get_noeud3()==no2) nn2=tet1->get_noeud3();
538     if (tet1->get_noeud4()==no2) nn2=tet1->get_noeud4();
539     if (tet1->get_noeud1()!=nn1)
540     if (tet1->get_noeud1()!=nn2) nn3=tet1->get_noeud1();
541     if (tet1->get_noeud2()!=nn1)
542     if (tet1->get_noeud2()!=nn2) nn3=tet1->get_noeud2();
543     if (tet1->get_noeud3()!=nn1)
544     if (tet1->get_noeud3()!=nn2) nn3=tet1->get_noeud3();
545     if (tet1->get_noeud4()!=nn1)
546     if (tet1->get_noeud4()!=nn2) nn3=tet1->get_noeud4();
547     if (tet1->get_noeud1()!=nn1)
548     if (tet1->get_noeud1()!=nn2)
549     if (tet1->get_noeud1()!=nn3) nn4=tet1->get_noeud1();
550     if (tet1->get_noeud2()!=nn1)
551     if (tet1->get_noeud2()!=nn2)
552     if (tet1->get_noeud2()!=nn3) nn4=tet1->get_noeud2();
553     if (tet1->get_noeud3()!=nn1)
554     if (tet1->get_noeud3()!=nn2)
555     if (tet1->get_noeud3()!=nn3) nn4=tet1->get_noeud3();
556     if (tet1->get_noeud4()!=nn1)
557     if (tet1->get_noeud4()!=nn2)
558     if (tet1->get_noeud4()!=nn3) nn4=tet1->get_noeud4();
559     double angle=get_angle_par_noeud<MG_NOEUD*>(nn1,nn2,nn3,nn1,nn2,nn4);
560     if (angle>M_PI) angle=2.*M_PI-angle;
561     angletot=angletot+angle;
562    
563     }
564     }
565     if ((angletot<2*M_PI-0.001) || (angletot>2*M_PI+0.001))
566     {
567 francois 805 char mess[255];
568 francois 425 sprintf(mess," segment %lu angle matiere autour de %.2lf",seg->get_id(),angletot); affiche(mess);
569     nbsegincorrect++;
570     }
571     else nbsegcorrect++;
572     }
573 francois 805 }
574 francois 514
575    
576 francois 577 void MAILLEUR_ANALYSE::compare_maillage_carte_isotrope(FCT_TAILLE* carte,char *nom,double *tab)
577 francois 514 {
578     MG_GESTIONNAIRE *gest=mai->get_gestionnaire();
579 francois 577 LISTE_MG_TETRA::iterator ittet;
580     double vol=0;
581     for (MG_TETRA* tet=mai->get_premier_tetra(ittet);tet!=NULL;tet=mai->get_suivant_tetra(ittet))
582     vol=vol+carte->calcul_volume_tetra_metrique(tet);
583     tab[0]=vol;
584     tab[1]=mai->get_nb_mg_tetra();
585     tab[2]=(tab[1]-tab[0])/tab[0]*100;
586 francois 514 MG_SOLUTION *sol=new MG_SOLUTION(mai,4,nom,1,"Comparaison_taille",MAGIC::ENTITE_SOLUTION::ENTITE_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
587     gest->ajouter_mg_solution(sol);
588     sol->change_legende(0,(char*)"taille_reelle");
589     sol->change_legende(1,(char*)"taille_carte");
590     sol->change_legende(2,(char*)"erreur");
591     sol->change_legende(3,(char*)"erreur absolue");
592     LISTE_MG_NOEUD::iterator it;
593     int i=0;
594     for (MG_NOEUD* no=mai->get_premier_noeud(it);no!=NULL;no=mai->get_suivant_noeud(it))
595     {
596     double *xyz=no->get_coord();
597     double tenseur[9];
598 francois 532 carte->evaluer(xyz,tenseur);
599 francois 514 double cartetaille=1./sqrt(tenseur[0]);
600     int nbseg=no->get_lien_segment()->get_nb();
601     double taille=0.;
602 francois 532 int nbreelseg=0;
603 francois 514 for (int j=0;j<nbseg;j++)
604 francois 532 if (mai->get_mg_segmentid(no->get_lien_segment()->get(j)->get_id())!=NULL) {taille=taille+no->get_lien_segment()->get(j)->get_longueur();nbreelseg++;}
605     taille=taille/nbreelseg;
606 francois 514 double erreur=(taille-cartetaille)*100./cartetaille;
607     sol->ecrire(taille,i,0);
608     sol->ecrire(cartetaille,i,1);
609     sol->ecrire(erreur,i,2);
610     sol->ecrire(fabs(erreur),i,3);
611     i++;
612     }
613     }