ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur_analyse.cpp
Revision: 1008
Committed: Mon Mar 25 16:37:35 2019 UTC (6 years, 3 months ago) by francois
Original Path: magic/lib/mailleur_auto/src/mailleur_analyse.cpp
File size: 25952 byte(s)
Log Message:
parametrisation de la methode de calcul des polycristaux

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 francois 1008 MAILLEUR_ANALYSE::MAILLEUR_ANALYSE(FEM_MAILLAGE* m,OT_CPU* comp):MAILLEUR(comp),fem(m),mai(NULL),borne1(0.1),borne2(0.2),borne3(0.5),eps_angle_retourne(0.03)
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 francois 876 int nbtet4=0,nbtet10=0,nbhex8=0,nbhex20=0,nbpenta5=0,nbpenta15=0;
117 francois 828 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 francois 876 if (ele->get_nb_fem_noeud()==5) nbpenta5++;
124     if (ele->get_nb_fem_noeud()==15) nbpenta15++;
125 francois 828 }
126     if (degre==1) {sprintf(mess," %d tetra4",nbtet4); affiche(mess); }
127     if (degre==2) {sprintf(mess," %d tetra10",nbtet10); affiche(mess); }
128     if (degre==1) {sprintf(mess," %d hex8",nbhex8); affiche(mess); }
129     if (degre==2) {sprintf(mess," %d hex20",nbhex20); affiche(mess); }
130 francois 876 if (degre==1) {sprintf(mess," %d penta5",nbpenta5); affiche(mess); }
131     if (degre==2) {sprintf(mess," %d penta15",nbpenta15); affiche(mess); }
132 francois 828 affiche((char*)"");
133     sprintf(mess,"Qualité du maillage");
134     affiche(mess);
135     if (degre==1)
136     {
137     mai=fem->get_mg_maillage();
138     analyse_mg(nom);
139     }
140     if (degre==2)
141     {
142 francois 882 FEM_SOLUTION *sol=NULL;
143     if (nom!=NULL)
144     {
145     std::string nomsol=nom;nomsol=nomsol+".sol";
146 francois 883 if (fem->get_nb_fem_element3()==0) sol=new FEM_SOLUTION(fem,4,(char*)nomsol.c_str(),fem->get_nb_fem_element2(),"Qualité",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2,MAGIC::TYPE_SOLUTION::SCALAIRE);
147     if (fem->get_nb_fem_element3()!=0) sol=new FEM_SOLUTION(fem,4,(char*)nomsol.c_str(),fem->get_nb_fem_element3(),"Qualité",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT3,MAGIC::TYPE_SOLUTION::SCALAIRE);
148 francois 882 fem->get_mg_maillage()->get_gestionnaire()->ajouter_fem_solution(sol);
149     sol->change_legende(0,"jmin");
150     sol->change_legende(1,"jmax");
151     sol->change_legende(2,"distorsion");
152 francois 883 sol->change_legende(3,"Qlin");
153 francois 882 }
154 francois 828 if (nbtri6>0)
155     {
156     sprintf(mess," Qualité des triangle6");
157     affiche(mess);
158     FEM_MAILLAGE_QUADRATIQUE_OUTILS ot;
159     double jminmin=1e300,jmaxmax=1e-300,dismoy=0,dismin=1e300,dismax=-1e300;
160     int nb=0;
161 francois 882 int i=0;
162 francois 828 for (FEM_ELEMENT2 *ele=fem->get_premier_element2(it2);ele!=NULL;ele=fem->get_suivant_element2(it2))
163     {
164     if (ele->get_nb_fem_noeud()==6)
165     {
166     double jmin=ot.get_jmin(ele);
167     double jmax=ot.get_jmax(ele);
168     double dis=ot.get_distorsion2(ele);
169     nb++;
170     dismoy=dismoy+dis;
171     if (jmin<jminmin) jminmin=jmin;
172     if (jmax>jmaxmax) jmaxmax=jmax;
173     if (dis<dismin) dismin=dis;
174     if (dis>dismax) dismax=dis;
175 francois 882 if (sol!=NULL)
176     if (sol->get_entite_solution()==MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2)
177     {
178     sol->ecrire(jmin,i,0);
179     sol->ecrire(jmax,i,1);
180     sol->ecrire(dis,i,2);
181     }
182 francois 828 }
183 francois 882 i++;
184 francois 828 }
185 francois 915 sprintf(mess," Jacobien min %le",jminmin);
186 francois 828 affiche(mess);
187 francois 915 sprintf(mess," Jacobien max %le",jmaxmax);
188 francois 828 affiche(mess);
189     sprintf(mess," Distortion min %lf",dismin);
190     affiche(mess);
191     sprintf(mess," Distortion moy %lf",dismoy/nb);
192     affiche(mess);
193     sprintf(mess," Distortion max %lf",dismax);
194     affiche(mess);
195    
196     }
197     if (nbtet10>0)
198     {
199     sprintf(mess," Qualité des tetra10");
200     affiche(mess);
201     FEM_MAILLAGE_QUADRATIQUE_OUTILS ot;
202     double jminmin=1e300,jmaxmax=1e-300,dismoy=0,dismin=1e300,dismax=-1e300;
203     int nb=0;
204 francois 882 int i=0;
205 francois 828 for (FEM_ELEMENT3 *ele=fem->get_premier_element3(it3);ele!=NULL;ele=fem->get_suivant_element3(it3))
206     {
207     if (ele->get_nb_fem_noeud()==10)
208     {
209     double jmin=ot.get_jmin(ele);
210     double jmax=ot.get_jmax(ele);
211     double dis=ot.get_distorsion2(ele);
212     nb++;
213     dismoy=dismoy+dis;
214     if (jmin<jminmin) jminmin=jmin;
215     if (jmax>jmaxmax) jmaxmax=jmax;
216     if (dis<dismin) dismin=dis;
217     if (dis>dismax) dismax=dis;
218 francois 882 if (sol!=NULL)
219     if (sol->get_entite_solution()==MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT3)
220     {
221     sol->ecrire(jmin,i,0);
222     sol->ecrire(jmax,i,1);
223     sol->ecrire(dis,i,2);
224 francois 883 MG_TETRA* tet=(MG_TETRA*)ele->get_mg_element_maillage();
225     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());
226     sol->ecrire(qual,i,3);
227     if (jmin<0.)
228     {
229     sprintf(mess," Erreur maille numero %d identificateur %lu\n Jmin=%lf Jmax=%lf Distorsion=%lf Qualite lineaire=%lf\n",i,ele->get_id(),jmin,jmax,dis,qual);
230     affiche(mess);
231     }
232 francois 882 }
233 francois 828 }
234 francois 882 i++;
235 francois 828 }
236 francois 915 sprintf(mess," Jacobien min %le",jminmin);
237 francois 828 affiche(mess);
238 francois 915 sprintf(mess," Jacobien max %le",jmaxmax);
239 francois 828 affiche(mess);
240     sprintf(mess," Distortion min %lf",dismin);
241     affiche(mess);
242     sprintf(mess," Distortion moy %lf",dismoy/nb);
243     affiche(mess);
244     sprintf(mess," Distortion max %lf",dismax);
245     affiche(mess);
246    
247     }
248 francois 882 if (nom!=NULL)
249     {
250     std::string nomf=nom;
251     nomf=nomf+".magic";
252     fem->get_mg_maillage()->get_gestionnaire()->enregistrer((char*)nomf.c_str());
253     }
254    
255 francois 828 }
256     }
257    
258    
259     void MAILLEUR_ANALYSE::analyse_mg(char *nom)
260     {
261     char mess[500];
262 francois 425 sprintf(mess,"ANALYSEUR DE MG MAILLAGE 3D");
263     affiche(mess);
264     sprintf(mess,"---------------------------");
265     affiche(mess);
266 francois 660 affiche((char*)"");
267 francois 828 sprintf(mess,"Constitution du maillage");
268 francois 805 affiche(mess);
269     affiche((char*)"");
270     int nbnofront,nbsegfront,nbtrifront,nbquadfront;
271     denombre_maillage(nbnofront,nbsegfront,nbtrifront,nbquadfront);
272     sprintf(mess," Maillage numero : %lu",mai->get_id());affiche(mess);
273     sprintf(mess," %d noeuds",mai->get_nb_mg_noeud()); affiche(mess);
274     sprintf(mess," %d segments",mai->get_nb_mg_segment()); affiche(mess);
275     sprintf(mess," %d triangles",mai->get_nb_mg_triangle()); affiche(mess);
276     sprintf(mess," %d quadrangles",mai->get_nb_mg_quadrangle()); affiche(mess);
277     sprintf(mess," %d tetras",mai->get_nb_mg_tetra()); affiche(mess);
278     sprintf(mess," %d hexas",mai->get_nb_mg_hexa()); affiche(mess); affiche((char*)"");
279     sprintf(mess," %d noeuds frontiere",nbnofront); affiche(mess);
280     sprintf(mess," %d segments frontiere",nbsegfront); affiche(mess);
281     sprintf(mess," %d triangles frontiere",nbtrifront); affiche(mess);
282     sprintf(mess," %d quadrangles frontiere",nbquadfront); affiche(mess);affiche((char*)"");
283 francois 425
284 francois 805
285     MG_SOLUTION *sol=NULL;
286     if (nom!=NULL)
287     {
288     std::string nomsol=nom;nomsol=nomsol+".sol";
289     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);
290     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);
291     mai->get_gestionnaire()->ajouter_mg_solution(sol);
292     }
293    
294    
295     if (nbtrifront!=0)
296     {
297     affiche((char*)"");
298     sprintf(mess,"Analyse maillage triangulaire 2D");
299     affiche(mess);
300     affiche((char*)"");
301     double qualmin,qualmax,qualmoy;
302     int tab[7];
303     analyse_qualite_maillage_2D(sol,qualmin,qualmax,qualmoy,tab);
304     sprintf(mess," Nombre de triangles avec des noeuds fusionnés : %d",tab[6]); affiche(mess);
305     sprintf(mess," Nombre de triangles avec 1 voisin manquant : %d",tab[5]); affiche(mess);
306     sprintf(mess," critere d'angle limite pour l'inversion : %lf",eps_angle_retourne); affiche(mess);
307     sprintf(mess," qualite moyenne des triangles de frontiere : %lf",qualmoy); affiche(mess);
308     sprintf(mess," qualite min des triangles de frontiere : %lf",qualmin); affiche(mess);
309     sprintf(mess," qualite max des triangles de frontiere : %lf",qualmax); affiche(mess);
310     sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne3,tab[4],tab[4]*100./nbtrifront); affiche(mess);
311     sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne2,tab[3],tab[3]*100./nbtrifront); affiche(mess);
312     sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",borne1,tab[2],tab[2]*100./nbtrifront); affiche(mess);
313     sprintf(mess," nombre de triangles de frontiere >%lf : %d (%.2lf%%)",0.,tab[1],tab[1]*100./nbtrifront); affiche(mess);
314     sprintf(mess," nombre de triangles de frontiere retournes : %d (%.2lf%%)",tab[0],tab[0]*100./nbtrifront); affiche(mess); affiche((char*)"");
315     }
316     if (mai->get_nb_mg_tetra()!=0)
317     {
318     affiche((char*)"");
319 francois 830 sprintf(mess,"Analyse maillage tetraedrique 3D");
320 francois 805 affiche(mess);
321     affiche((char*)"");
322     double qualmin,qualmax,qualmoy;
323     int tab[5];
324     analyse_qualite_maillage_3D(sol,qualmin,qualmax,qualmoy,tab);
325     sprintf(mess," qualite moyenne des tetra : %lf",qualmoy); affiche(mess);
326 francois 1008 if (qualmin>1e-6)
327     {sprintf(mess," qualite min des tetra : %lf",qualmin); affiche(mess); }
328     else
329     {sprintf(mess," qualite min des tetra : %le",qualmin); affiche(mess); }
330 francois 805 sprintf(mess," qualite max des tetra : %lf",qualmax); affiche(mess);
331     sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne3,tab[4],tab[4]*100./mai->get_nb_mg_tetra()); affiche(mess);
332     sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne2,tab[3],tab[3]*100./mai->get_nb_mg_tetra()); affiche(mess);
333     sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne1,tab[2],tab[2]*100./mai->get_nb_mg_tetra()); affiche(mess);
334     sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",0.,tab[1],tab[1]*100./mai->get_nb_mg_tetra()); affiche(mess);
335     sprintf(mess," nombre de tetra retournes : %d (%.2lf%%)",tab[0],tab[0]*100./mai->get_nb_mg_tetra()); affiche(mess);
336     affiche((char*)"");
337 francois 830 sprintf(mess,"Validite maillage tétraedrique 3D");
338 francois 805 affiche(mess);
339     affiche((char*)"");
340     int nbsegcorrect,nbsegincorrect;
341     analyse_validite_maillage_3D(nbsegcorrect,nbsegincorrect);
342    
343     sprintf(mess," segment correct %d(%.2lf%%) segment incorrect %d(%.2lf%%) ",nbsegcorrect,nbsegcorrect*100./(nbsegcorrect+nbsegincorrect),nbsegincorrect,nbsegincorrect*100./(nbsegcorrect+nbsegincorrect)); affiche(mess);
344    
345     }
346    
347     if (nom!=NULL)
348     {
349     std::string nomf=nom;
350     nomf=nomf+".magic";
351     mai->get_gestionnaire()->enregistrer((char*)nomf.c_str());
352     }
353    
354     }
355    
356     void MAILLEUR_ANALYSE::denombre_maillage(int &nbnofront,int &nbsegfront,int &nbtrifront,int &nbquadfront)
357     {
358 francois 425 LISTE_MG_NOEUD::iterator itn;
359     for (MG_NOEUD* no=mai->get_premier_noeud(itn);no!=NULL;no=mai->get_suivant_noeud(itn))
360     no->change_nouveau_numero(0);
361    
362     LISTE_MG_SEGMENT::iterator its;
363     for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
364     seg->change_nouveau_numero(0);
365    
366    
367     LISTE_MG_TRIANGLE::iterator ittr;
368     for (MG_TRIANGLE* tri=mai->get_premier_triangle(ittr);tri!=NULL;tri=mai->get_suivant_triangle(ittr))
369     tri->change_nouveau_numero(0);
370    
371     LISTE_MG_QUADRANGLE::iterator itq;
372     for (MG_QUADRANGLE* quad=mai->get_premier_quadrangle(itq);quad!=NULL;quad=mai->get_suivant_quadrangle(itq))
373     quad->change_nouveau_numero(0);
374    
375     LISTE_MG_TETRA::iterator itt;
376     for (MG_TETRA* tet=mai->get_premier_tetra(itt);tet!=NULL;tet=mai->get_suivant_tetra(itt))
377     {
378     tet->get_triangle1()->change_nouveau_numero(tet->get_triangle1()->get_nouveau_numero()+1);
379     tet->get_triangle2()->change_nouveau_numero(tet->get_triangle2()->get_nouveau_numero()+1);
380     tet->get_triangle3()->change_nouveau_numero(tet->get_triangle3()->get_nouveau_numero()+1);
381     tet->get_triangle4()->change_nouveau_numero(tet->get_triangle4()->get_nouveau_numero()+1);
382     }
383    
384     LISTE_MG_HEXA::iterator ith;
385     for (MG_HEXA* hex=mai->get_premier_hexa(ith);hex!=NULL;hex=mai->get_suivant_hexa(ith))
386     {
387     hex->get_quadrangle1()->change_nouveau_numero(hex->get_quadrangle1()->get_nouveau_numero()+1);
388     hex->get_quadrangle2()->change_nouveau_numero(hex->get_quadrangle2()->get_nouveau_numero()+1);
389     hex->get_quadrangle3()->change_nouveau_numero(hex->get_quadrangle3()->get_nouveau_numero()+1);
390     hex->get_quadrangle4()->change_nouveau_numero(hex->get_quadrangle4()->get_nouveau_numero()+1);
391     hex->get_quadrangle5()->change_nouveau_numero(hex->get_quadrangle5()->get_nouveau_numero()+1);
392     hex->get_quadrangle6()->change_nouveau_numero(hex->get_quadrangle6()->get_nouveau_numero()+1);
393     }
394 francois 805 nbtrifront=0;
395 francois 425 for (MG_TRIANGLE* tri=mai->get_premier_triangle(ittr);tri!=NULL;tri=mai->get_suivant_triangle(ittr))
396     {
397     if (mai->get_nb_mg_tetra()==0)
398     tri->change_nouveau_numero(1);
399     if (tri->get_nouveau_numero()==1)
400     {
401     nbtrifront++;
402     tri->get_segment1()->change_nouveau_numero(tri->get_segment1()->get_nouveau_numero()+1);
403     tri->get_segment2()->change_nouveau_numero(tri->get_segment2()->get_nouveau_numero()+1);
404     tri->get_segment3()->change_nouveau_numero(tri->get_segment3()->get_nouveau_numero()+1);
405     tri->get_noeud1()->change_nouveau_numero(1);
406     tri->get_noeud2()->change_nouveau_numero(1);
407     tri->get_noeud3()->change_nouveau_numero(1);
408    
409     }
410     }
411 francois 805 nbquadfront=0;
412 francois 425 for (MG_QUADRANGLE* quad=mai->get_premier_quadrangle(itq);quad!=NULL;quad=mai->get_suivant_quadrangle(itq))
413     {
414     if (mai->get_nb_mg_hexa()==0)
415     quad->change_nouveau_numero(1);
416     if (quad->get_nouveau_numero()==1)
417     {
418     nbquadfront++;
419     quad->get_segment1()->change_nouveau_numero(quad->get_segment1()->get_nouveau_numero()+1);
420     quad->get_segment2()->change_nouveau_numero(quad->get_segment2()->get_nouveau_numero()+1);
421     quad->get_segment3()->change_nouveau_numero(quad->get_segment3()->get_nouveau_numero()+1);
422     quad->get_segment4()->change_nouveau_numero(quad->get_segment4()->get_nouveau_numero()+1);
423     quad->get_noeud1()->change_nouveau_numero(1);
424     quad->get_noeud2()->change_nouveau_numero(1);
425     quad->get_noeud3()->change_nouveau_numero(1);
426     quad->get_noeud4()->change_nouveau_numero(1);
427     }
428     }
429 francois 805 nbsegfront=0;
430 francois 425 for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
431     if (seg->get_nouveau_numero()==2) nbsegfront++;
432 francois 805 nbnofront=0;
433 francois 425 for (MG_NOEUD* no=mai->get_premier_noeud(itn);no!=NULL;no=mai->get_suivant_noeud(itn))
434     if (no->get_nouveau_numero()==1) nbnofront++;
435 francois 805 }
436 francois 425
437 francois 805 void MAILLEUR_ANALYSE::analyse_qualite_maillage_2D(MG_SOLUTION *sol,double &qualmin,double &qualmax,double &qualmoy,int *tab)
438     {
439     int nbnofront,nbsegfront,nbtrifront,nbquadfront;
440     denombre_maillage(nbnofront,nbsegfront,nbtrifront,nbquadfront);
441     qualmin=1e300,qualmax=-1e300,qualmoy=0.;
442     tab[0]=0;tab[1]=0;tab[2]=0;tab[3]=0;tab[4]=0;tab[5]=0;tab[6]=0;
443 francois 425 int itri=0;
444 francois 805 LISTE_MG_TRIANGLE::iterator ittr;
445 francois 425 for (MG_TRIANGLE* tri=mai->get_premier_triangle(ittr);tri!=NULL;tri=mai->get_suivant_triangle(ittr))
446     if (tri->get_nouveau_numero()==1)
447     {
448     double qual=OPERATEUR::qualite_triangle(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord(),tri->get_noeud3()->get_coord());
449     MG_NOEUD * no1=tri->get_noeud1();
450 francois 429 MG_NOEUD * no2=tri->get_noeud2();
451     MG_NOEUD * no3=tri->get_noeud3();
452     if (no1==no2)
453     if (no1==no3)
454 francois 805 tab[6]++;
455 francois 425 int nb1=no1->get_lien_triangle()->get_nb();
456     int nb2=no2->get_lien_triangle()->get_nb();
457     int nb3=no3->get_lien_triangle()->get_nb();
458     OT_VECTEUR_3D vec1(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord());
459     OT_VECTEUR_3D vec2(tri->get_noeud1()->get_coord(),tri->get_noeud3()->get_coord());
460     OT_VECTEUR_3D normal=vec1&vec2;
461     normal.norme();
462     int nbretourne=0;
463 francois 429 int nbvoisin=0;
464 francois 425 for (int i=0;i<nb1;i++)
465     for (int j=0;j<nb2;j++)
466     {
467     MG_TRIANGLE* tri1=no1->get_lien_triangle()->get(i);
468     MG_TRIANGLE* tri2=no2->get_lien_triangle()->get(j);
469 francois 429 if ( (tri1==tri2) && (tri1!=tri) && (tri1->get_nouveau_numero()==1))
470 francois 425 {
471     OT_VECTEUR_3D vec1tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud2()->get_coord());
472     OT_VECTEUR_3D vec2tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud3()->get_coord());
473     OT_VECTEUR_3D normaltmp=vec1tmp&vec2tmp;
474     normaltmp.norme();
475     double psca=normal*normaltmp;
476 francois 426 if (psca<-cos(eps_angle_retourne)) nbretourne++;
477 francois 429 nbvoisin++;
478 francois 425 }
479     }
480 francois 429 for (int i=0;i<nb1;i++)
481     for (int j=0;j<nb3;j++)
482     {
483     MG_TRIANGLE* tri1=no1->get_lien_triangle()->get(i);
484     MG_TRIANGLE* tri2=no3->get_lien_triangle()->get(j);
485     if ( (tri1==tri2) && (tri1!=tri) && (tri1->get_nouveau_numero()==1))
486     {
487     OT_VECTEUR_3D vec1tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud2()->get_coord());
488     OT_VECTEUR_3D vec2tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud3()->get_coord());
489     OT_VECTEUR_3D normaltmp=vec1tmp&vec2tmp;
490     normaltmp.norme();
491     double psca=normal*normaltmp;
492     if (psca<-cos(eps_angle_retourne)) nbretourne++;
493     nbvoisin++;
494     }
495     }
496     for (int i=0;i<nb2;i++)
497     for (int j=0;j<nb3;j++)
498     {
499     MG_TRIANGLE* tri1=no2->get_lien_triangle()->get(i);
500     MG_TRIANGLE* tri2=no3->get_lien_triangle()->get(j);
501     if ( (tri1==tri2) && (tri1!=tri) && (tri1->get_nouveau_numero()==1))
502     {
503     OT_VECTEUR_3D vec1tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud2()->get_coord());
504     OT_VECTEUR_3D vec2tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud3()->get_coord());
505     OT_VECTEUR_3D normaltmp=vec1tmp&vec2tmp;
506     normaltmp.norme();
507     double psca=normal*normaltmp;
508     if (psca<-cos(eps_angle_retourne)) nbretourne++;
509     nbvoisin++;
510     }
511     }
512     if (nbvoisin!=3)
513 francois 805 tab[5]++;
514 francois 425 if (nbretourne>1) qual=-qual;
515     qualmoy=qualmoy+qual;
516     if (qual<qualmin) qualmin=qual;
517     if (qual>qualmax) qualmax=qual;
518 francois 805 if (qual<0.) tab[0]++;
519     else if (qual<borne1) tab[1]++;
520     else if (qual<borne2) tab[2]++;
521     else if (qual<borne3) tab[3]++;
522     else tab[4]++;
523 francois 425 if (sol!=NULL)
524     if (mai->get_nb_mg_tetra()==0) sol->ecrire(qual,itri,0);
525     itri++;
526     }
527 francois 805 qualmoy=qualmoy/nbtrifront;
528     }
529    
530     void MAILLEUR_ANALYSE::analyse_qualite_maillage_3D(MG_SOLUTION *sol,double &qualmin,double &qualmax,double &qualmoy,int *tab)
531 francois 425 {
532     qualmin=1e300,qualmax=-1e300,qualmoy=0.;
533 francois 805 tab[0]=0;tab[1]=0;tab[2]=0;tab[3]=0;tab[4]=0;
534     qualmin=1e300,qualmax=-1e300,qualmoy=0.;
535 francois 425 int itet=0;
536 francois 805 LISTE_MG_TETRA::iterator itt;
537 francois 425 for (MG_TETRA* tet=mai->get_premier_tetra(itt);tet!=NULL;tet=mai->get_suivant_tetra(itt))
538     {
539     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());
540     OT_VECTEUR_3D vec1(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord());
541     OT_VECTEUR_3D vec2(tet->get_noeud1()->get_coord(),tet->get_noeud3()->get_coord());
542     OT_VECTEUR_3D vec3(tet->get_noeud1()->get_coord(),tet->get_noeud4()->get_coord());
543     OT_VECTEUR_3D pvec=vec1&vec2;
544     double psca=pvec*vec3;
545     if (psca<0.) qual=-qual;
546     qualmoy=qualmoy+qual;
547     if (qual<qualmin) qualmin=qual;
548     if (qual>qualmax) qualmax=qual;
549 francois 805 if (qual<0.) tab[0]++;
550     else if (qual<borne1) tab[1]++;
551     else if (qual<borne2) tab[2]++;
552     else if (qual<borne3) tab[3]++;
553     else tab[4]++;
554 francois 425 if (sol!=NULL)
555     sol->ecrire(qual,itet,0);
556     itet++;
557     }
558 francois 805 qualmoy=qualmoy/mai->get_nb_mg_tetra();
559     }
560    
561     void MAILLEUR_ANALYSE::analyse_validite_maillage_3D(int &nbsegcorrect,int &nbsegincorrect)
562     {
563    
564     nbsegcorrect=0;
565     nbsegincorrect=0;
566     LISTE_MG_SEGMENT::iterator its;
567 francois 425 for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
568     if (seg->get_nouveau_numero()!=2)
569     {
570     MG_NOEUD* no1=seg->get_noeud1();
571     MG_NOEUD* no2=seg->get_noeud2();
572     int nb1=no1->get_lien_tetra()->get_nb();
573     int nb2=no2->get_lien_tetra()->get_nb();
574     double angletot=0;
575     for (int i=0;i<nb1;i++)
576     for (int j=0;j<nb2;j++)
577     {
578     MG_TETRA* tet1=no1->get_lien_tetra()->get(i);
579     MG_TETRA* tet2=no2->get_lien_tetra()->get(j);
580     if (tet1==tet2)
581     {
582     MG_NOEUD *nn1,*nn2,*nn3,*nn4;
583     if (tet1->get_noeud1()==no1) nn1=tet1->get_noeud1();
584     if (tet1->get_noeud2()==no1) nn1=tet1->get_noeud2();
585     if (tet1->get_noeud3()==no1) nn1=tet1->get_noeud3();
586     if (tet1->get_noeud4()==no1) nn1=tet1->get_noeud4();
587     if (tet1->get_noeud1()==no2) nn2=tet1->get_noeud1();
588     if (tet1->get_noeud2()==no2) nn2=tet1->get_noeud2();
589     if (tet1->get_noeud3()==no2) nn2=tet1->get_noeud3();
590     if (tet1->get_noeud4()==no2) nn2=tet1->get_noeud4();
591     if (tet1->get_noeud1()!=nn1)
592     if (tet1->get_noeud1()!=nn2) nn3=tet1->get_noeud1();
593     if (tet1->get_noeud2()!=nn1)
594     if (tet1->get_noeud2()!=nn2) nn3=tet1->get_noeud2();
595     if (tet1->get_noeud3()!=nn1)
596     if (tet1->get_noeud3()!=nn2) nn3=tet1->get_noeud3();
597     if (tet1->get_noeud4()!=nn1)
598     if (tet1->get_noeud4()!=nn2) nn3=tet1->get_noeud4();
599     if (tet1->get_noeud1()!=nn1)
600     if (tet1->get_noeud1()!=nn2)
601     if (tet1->get_noeud1()!=nn3) nn4=tet1->get_noeud1();
602     if (tet1->get_noeud2()!=nn1)
603     if (tet1->get_noeud2()!=nn2)
604     if (tet1->get_noeud2()!=nn3) nn4=tet1->get_noeud2();
605     if (tet1->get_noeud3()!=nn1)
606     if (tet1->get_noeud3()!=nn2)
607     if (tet1->get_noeud3()!=nn3) nn4=tet1->get_noeud3();
608     if (tet1->get_noeud4()!=nn1)
609     if (tet1->get_noeud4()!=nn2)
610     if (tet1->get_noeud4()!=nn3) nn4=tet1->get_noeud4();
611     double angle=get_angle_par_noeud<MG_NOEUD*>(nn1,nn2,nn3,nn1,nn2,nn4);
612     if (angle>M_PI) angle=2.*M_PI-angle;
613     angletot=angletot+angle;
614    
615     }
616     }
617     if ((angletot<2*M_PI-0.001) || (angletot>2*M_PI+0.001))
618     {
619 francois 805 char mess[255];
620 francois 425 sprintf(mess," segment %lu angle matiere autour de %.2lf",seg->get_id(),angletot); affiche(mess);
621     nbsegincorrect++;
622     }
623     else nbsegcorrect++;
624     }
625 francois 805 }
626 francois 514
627    
628 francois 577 void MAILLEUR_ANALYSE::compare_maillage_carte_isotrope(FCT_TAILLE* carte,char *nom,double *tab)
629 francois 514 {
630     MG_GESTIONNAIRE *gest=mai->get_gestionnaire();
631 francois 577 LISTE_MG_TETRA::iterator ittet;
632     double vol=0;
633     for (MG_TETRA* tet=mai->get_premier_tetra(ittet);tet!=NULL;tet=mai->get_suivant_tetra(ittet))
634     vol=vol+carte->calcul_volume_tetra_metrique(tet);
635     tab[0]=vol;
636     tab[1]=mai->get_nb_mg_tetra();
637     tab[2]=(tab[1]-tab[0])/tab[0]*100;
638 francois 514 MG_SOLUTION *sol=new MG_SOLUTION(mai,4,nom,1,"Comparaison_taille",MAGIC::ENTITE_SOLUTION::ENTITE_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
639     gest->ajouter_mg_solution(sol);
640     sol->change_legende(0,(char*)"taille_reelle");
641     sol->change_legende(1,(char*)"taille_carte");
642     sol->change_legende(2,(char*)"erreur");
643     sol->change_legende(3,(char*)"erreur absolue");
644     LISTE_MG_NOEUD::iterator it;
645     int i=0;
646     for (MG_NOEUD* no=mai->get_premier_noeud(it);no!=NULL;no=mai->get_suivant_noeud(it))
647     {
648     double *xyz=no->get_coord();
649     double tenseur[9];
650 francois 532 carte->evaluer(xyz,tenseur);
651 francois 514 double cartetaille=1./sqrt(tenseur[0]);
652     int nbseg=no->get_lien_segment()->get_nb();
653     double taille=0.;
654 francois 532 int nbreelseg=0;
655 francois 514 for (int j=0;j<nbseg;j++)
656 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++;}
657     taille=taille/nbreelseg;
658 francois 514 double erreur=(taille-cartetaille)*100./cartetaille;
659     sol->ecrire(taille,i,0);
660     sol->ecrire(cartetaille,i,1);
661     sol->ecrire(erreur,i,2);
662     sol->ecrire(fabs(erreur),i,3);
663     i++;
664     }
665     }