ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur_analyse.cpp
Revision: 883
Committed: Thu Apr 20 13:38:18 2017 UTC (8 years ago) by francois
Original Path: magic/lib/mailleur_auto/src/mailleur_analyse.cpp
File size: 25783 byte(s)
Log Message:
Creation d'un mailleur FEM pour pouvoir avoir des stratégies paramétrées de maillage. L'ancienne méthode construit disparait et est remplacée par l'utilisation d'un MAILLEUR_FEM.
Stratégie de quadratisation mise en place : déplacer les noeuds pour s'assurer que les tetras quadratiques présentent une distortion au dessu d'une borne inférieure. Mais ces noeuds quittent la géométrie.

Les paramètres dans ~/.magic.

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 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     sprintf(mess," Jacobien min %lf",jminmin);
186     affiche(mess);
187     sprintf(mess," Jacobien max %lf",jmaxmax);
188     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     sprintf(mess," Jacobien min %lf",jminmin);
237     affiche(mess);
238     sprintf(mess," Jacobien max %lf",jmaxmax);
239     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     sprintf(mess," qualite min des tetra : %lf",qualmin); affiche(mess);
327     sprintf(mess," qualite max des tetra : %lf",qualmax); affiche(mess);
328     sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne3,tab[4],tab[4]*100./mai->get_nb_mg_tetra()); affiche(mess);
329     sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne2,tab[3],tab[3]*100./mai->get_nb_mg_tetra()); affiche(mess);
330     sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",borne1,tab[2],tab[2]*100./mai->get_nb_mg_tetra()); affiche(mess);
331     sprintf(mess," nombre de tetra >%lf : %d (%.2lf%%)",0.,tab[1],tab[1]*100./mai->get_nb_mg_tetra()); affiche(mess);
332     sprintf(mess," nombre de tetra retournes : %d (%.2lf%%)",tab[0],tab[0]*100./mai->get_nb_mg_tetra()); affiche(mess);
333     affiche((char*)"");
334 francois 830 sprintf(mess,"Validite maillage tétraedrique 3D");
335 francois 805 affiche(mess);
336     affiche((char*)"");
337     int nbsegcorrect,nbsegincorrect;
338     analyse_validite_maillage_3D(nbsegcorrect,nbsegincorrect);
339    
340     sprintf(mess," segment correct %d(%.2lf%%) segment incorrect %d(%.2lf%%) ",nbsegcorrect,nbsegcorrect*100./(nbsegcorrect+nbsegincorrect),nbsegincorrect,nbsegincorrect*100./(nbsegcorrect+nbsegincorrect)); affiche(mess);
341    
342     }
343    
344     if (nom!=NULL)
345     {
346     std::string nomf=nom;
347     nomf=nomf+".magic";
348     mai->get_gestionnaire()->enregistrer((char*)nomf.c_str());
349     }
350    
351     }
352    
353     void MAILLEUR_ANALYSE::denombre_maillage(int &nbnofront,int &nbsegfront,int &nbtrifront,int &nbquadfront)
354     {
355 francois 425 LISTE_MG_NOEUD::iterator itn;
356     for (MG_NOEUD* no=mai->get_premier_noeud(itn);no!=NULL;no=mai->get_suivant_noeud(itn))
357     no->change_nouveau_numero(0);
358    
359     LISTE_MG_SEGMENT::iterator its;
360     for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
361     seg->change_nouveau_numero(0);
362    
363    
364     LISTE_MG_TRIANGLE::iterator ittr;
365     for (MG_TRIANGLE* tri=mai->get_premier_triangle(ittr);tri!=NULL;tri=mai->get_suivant_triangle(ittr))
366     tri->change_nouveau_numero(0);
367    
368     LISTE_MG_QUADRANGLE::iterator itq;
369     for (MG_QUADRANGLE* quad=mai->get_premier_quadrangle(itq);quad!=NULL;quad=mai->get_suivant_quadrangle(itq))
370     quad->change_nouveau_numero(0);
371    
372     LISTE_MG_TETRA::iterator itt;
373     for (MG_TETRA* tet=mai->get_premier_tetra(itt);tet!=NULL;tet=mai->get_suivant_tetra(itt))
374     {
375     tet->get_triangle1()->change_nouveau_numero(tet->get_triangle1()->get_nouveau_numero()+1);
376     tet->get_triangle2()->change_nouveau_numero(tet->get_triangle2()->get_nouveau_numero()+1);
377     tet->get_triangle3()->change_nouveau_numero(tet->get_triangle3()->get_nouveau_numero()+1);
378     tet->get_triangle4()->change_nouveau_numero(tet->get_triangle4()->get_nouveau_numero()+1);
379     }
380    
381     LISTE_MG_HEXA::iterator ith;
382     for (MG_HEXA* hex=mai->get_premier_hexa(ith);hex!=NULL;hex=mai->get_suivant_hexa(ith))
383     {
384     hex->get_quadrangle1()->change_nouveau_numero(hex->get_quadrangle1()->get_nouveau_numero()+1);
385     hex->get_quadrangle2()->change_nouveau_numero(hex->get_quadrangle2()->get_nouveau_numero()+1);
386     hex->get_quadrangle3()->change_nouveau_numero(hex->get_quadrangle3()->get_nouveau_numero()+1);
387     hex->get_quadrangle4()->change_nouveau_numero(hex->get_quadrangle4()->get_nouveau_numero()+1);
388     hex->get_quadrangle5()->change_nouveau_numero(hex->get_quadrangle5()->get_nouveau_numero()+1);
389     hex->get_quadrangle6()->change_nouveau_numero(hex->get_quadrangle6()->get_nouveau_numero()+1);
390     }
391 francois 805 nbtrifront=0;
392 francois 425 for (MG_TRIANGLE* tri=mai->get_premier_triangle(ittr);tri!=NULL;tri=mai->get_suivant_triangle(ittr))
393     {
394     if (mai->get_nb_mg_tetra()==0)
395     tri->change_nouveau_numero(1);
396     if (tri->get_nouveau_numero()==1)
397     {
398     nbtrifront++;
399     tri->get_segment1()->change_nouveau_numero(tri->get_segment1()->get_nouveau_numero()+1);
400     tri->get_segment2()->change_nouveau_numero(tri->get_segment2()->get_nouveau_numero()+1);
401     tri->get_segment3()->change_nouveau_numero(tri->get_segment3()->get_nouveau_numero()+1);
402     tri->get_noeud1()->change_nouveau_numero(1);
403     tri->get_noeud2()->change_nouveau_numero(1);
404     tri->get_noeud3()->change_nouveau_numero(1);
405    
406     }
407     }
408 francois 805 nbquadfront=0;
409 francois 425 for (MG_QUADRANGLE* quad=mai->get_premier_quadrangle(itq);quad!=NULL;quad=mai->get_suivant_quadrangle(itq))
410     {
411     if (mai->get_nb_mg_hexa()==0)
412     quad->change_nouveau_numero(1);
413     if (quad->get_nouveau_numero()==1)
414     {
415     nbquadfront++;
416     quad->get_segment1()->change_nouveau_numero(quad->get_segment1()->get_nouveau_numero()+1);
417     quad->get_segment2()->change_nouveau_numero(quad->get_segment2()->get_nouveau_numero()+1);
418     quad->get_segment3()->change_nouveau_numero(quad->get_segment3()->get_nouveau_numero()+1);
419     quad->get_segment4()->change_nouveau_numero(quad->get_segment4()->get_nouveau_numero()+1);
420     quad->get_noeud1()->change_nouveau_numero(1);
421     quad->get_noeud2()->change_nouveau_numero(1);
422     quad->get_noeud3()->change_nouveau_numero(1);
423     quad->get_noeud4()->change_nouveau_numero(1);
424     }
425     }
426 francois 805 nbsegfront=0;
427 francois 425 for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
428     if (seg->get_nouveau_numero()==2) nbsegfront++;
429 francois 805 nbnofront=0;
430 francois 425 for (MG_NOEUD* no=mai->get_premier_noeud(itn);no!=NULL;no=mai->get_suivant_noeud(itn))
431     if (no->get_nouveau_numero()==1) nbnofront++;
432 francois 805 }
433 francois 425
434 francois 805 void MAILLEUR_ANALYSE::analyse_qualite_maillage_2D(MG_SOLUTION *sol,double &qualmin,double &qualmax,double &qualmoy,int *tab)
435     {
436     int nbnofront,nbsegfront,nbtrifront,nbquadfront;
437     denombre_maillage(nbnofront,nbsegfront,nbtrifront,nbquadfront);
438     qualmin=1e300,qualmax=-1e300,qualmoy=0.;
439     tab[0]=0;tab[1]=0;tab[2]=0;tab[3]=0;tab[4]=0;tab[5]=0;tab[6]=0;
440 francois 425 int itri=0;
441 francois 805 LISTE_MG_TRIANGLE::iterator ittr;
442 francois 425 for (MG_TRIANGLE* tri=mai->get_premier_triangle(ittr);tri!=NULL;tri=mai->get_suivant_triangle(ittr))
443     if (tri->get_nouveau_numero()==1)
444     {
445     double qual=OPERATEUR::qualite_triangle(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord(),tri->get_noeud3()->get_coord());
446     MG_NOEUD * no1=tri->get_noeud1();
447 francois 429 MG_NOEUD * no2=tri->get_noeud2();
448     MG_NOEUD * no3=tri->get_noeud3();
449     if (no1==no2)
450     if (no1==no3)
451 francois 805 tab[6]++;
452 francois 425 int nb1=no1->get_lien_triangle()->get_nb();
453     int nb2=no2->get_lien_triangle()->get_nb();
454     int nb3=no3->get_lien_triangle()->get_nb();
455     OT_VECTEUR_3D vec1(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord());
456     OT_VECTEUR_3D vec2(tri->get_noeud1()->get_coord(),tri->get_noeud3()->get_coord());
457     OT_VECTEUR_3D normal=vec1&vec2;
458     normal.norme();
459     int nbretourne=0;
460 francois 429 int nbvoisin=0;
461 francois 425 for (int i=0;i<nb1;i++)
462     for (int j=0;j<nb2;j++)
463     {
464     MG_TRIANGLE* tri1=no1->get_lien_triangle()->get(i);
465     MG_TRIANGLE* tri2=no2->get_lien_triangle()->get(j);
466 francois 429 if ( (tri1==tri2) && (tri1!=tri) && (tri1->get_nouveau_numero()==1))
467 francois 425 {
468     OT_VECTEUR_3D vec1tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud2()->get_coord());
469     OT_VECTEUR_3D vec2tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud3()->get_coord());
470     OT_VECTEUR_3D normaltmp=vec1tmp&vec2tmp;
471     normaltmp.norme();
472     double psca=normal*normaltmp;
473 francois 426 if (psca<-cos(eps_angle_retourne)) nbretourne++;
474 francois 429 nbvoisin++;
475 francois 425 }
476     }
477 francois 429 for (int i=0;i<nb1;i++)
478     for (int j=0;j<nb3;j++)
479     {
480     MG_TRIANGLE* tri1=no1->get_lien_triangle()->get(i);
481     MG_TRIANGLE* tri2=no3->get_lien_triangle()->get(j);
482     if ( (tri1==tri2) && (tri1!=tri) && (tri1->get_nouveau_numero()==1))
483     {
484     OT_VECTEUR_3D vec1tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud2()->get_coord());
485     OT_VECTEUR_3D vec2tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud3()->get_coord());
486     OT_VECTEUR_3D normaltmp=vec1tmp&vec2tmp;
487     normaltmp.norme();
488     double psca=normal*normaltmp;
489     if (psca<-cos(eps_angle_retourne)) nbretourne++;
490     nbvoisin++;
491     }
492     }
493     for (int i=0;i<nb2;i++)
494     for (int j=0;j<nb3;j++)
495     {
496     MG_TRIANGLE* tri1=no2->get_lien_triangle()->get(i);
497     MG_TRIANGLE* tri2=no3->get_lien_triangle()->get(j);
498     if ( (tri1==tri2) && (tri1!=tri) && (tri1->get_nouveau_numero()==1))
499     {
500     OT_VECTEUR_3D vec1tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud2()->get_coord());
501     OT_VECTEUR_3D vec2tmp(tri1->get_noeud1()->get_coord(),tri1->get_noeud3()->get_coord());
502     OT_VECTEUR_3D normaltmp=vec1tmp&vec2tmp;
503     normaltmp.norme();
504     double psca=normal*normaltmp;
505     if (psca<-cos(eps_angle_retourne)) nbretourne++;
506     nbvoisin++;
507     }
508     }
509     if (nbvoisin!=3)
510 francois 805 tab[5]++;
511 francois 425 if (nbretourne>1) qual=-qual;
512     qualmoy=qualmoy+qual;
513     if (qual<qualmin) qualmin=qual;
514     if (qual>qualmax) qualmax=qual;
515 francois 805 if (qual<0.) tab[0]++;
516     else if (qual<borne1) tab[1]++;
517     else if (qual<borne2) tab[2]++;
518     else if (qual<borne3) tab[3]++;
519     else tab[4]++;
520 francois 425 if (sol!=NULL)
521     if (mai->get_nb_mg_tetra()==0) sol->ecrire(qual,itri,0);
522     itri++;
523     }
524 francois 805 qualmoy=qualmoy/nbtrifront;
525     }
526    
527     void MAILLEUR_ANALYSE::analyse_qualite_maillage_3D(MG_SOLUTION *sol,double &qualmin,double &qualmax,double &qualmoy,int *tab)
528 francois 425 {
529     qualmin=1e300,qualmax=-1e300,qualmoy=0.;
530 francois 805 tab[0]=0;tab[1]=0;tab[2]=0;tab[3]=0;tab[4]=0;
531     qualmin=1e300,qualmax=-1e300,qualmoy=0.;
532 francois 425 int itet=0;
533 francois 805 LISTE_MG_TETRA::iterator itt;
534 francois 425 for (MG_TETRA* tet=mai->get_premier_tetra(itt);tet!=NULL;tet=mai->get_suivant_tetra(itt))
535     {
536     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());
537     OT_VECTEUR_3D vec1(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord());
538     OT_VECTEUR_3D vec2(tet->get_noeud1()->get_coord(),tet->get_noeud3()->get_coord());
539     OT_VECTEUR_3D vec3(tet->get_noeud1()->get_coord(),tet->get_noeud4()->get_coord());
540     OT_VECTEUR_3D pvec=vec1&vec2;
541     double psca=pvec*vec3;
542     if (psca<0.) qual=-qual;
543     qualmoy=qualmoy+qual;
544     if (qual<qualmin) qualmin=qual;
545     if (qual>qualmax) qualmax=qual;
546 francois 805 if (qual<0.) tab[0]++;
547     else if (qual<borne1) tab[1]++;
548     else if (qual<borne2) tab[2]++;
549     else if (qual<borne3) tab[3]++;
550     else tab[4]++;
551 francois 425 if (sol!=NULL)
552     sol->ecrire(qual,itet,0);
553     itet++;
554     }
555 francois 805 qualmoy=qualmoy/mai->get_nb_mg_tetra();
556     }
557    
558     void MAILLEUR_ANALYSE::analyse_validite_maillage_3D(int &nbsegcorrect,int &nbsegincorrect)
559     {
560    
561     nbsegcorrect=0;
562     nbsegincorrect=0;
563     LISTE_MG_SEGMENT::iterator its;
564 francois 425 for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
565     if (seg->get_nouveau_numero()!=2)
566     {
567     MG_NOEUD* no1=seg->get_noeud1();
568     MG_NOEUD* no2=seg->get_noeud2();
569     int nb1=no1->get_lien_tetra()->get_nb();
570     int nb2=no2->get_lien_tetra()->get_nb();
571     double angletot=0;
572     for (int i=0;i<nb1;i++)
573     for (int j=0;j<nb2;j++)
574     {
575     MG_TETRA* tet1=no1->get_lien_tetra()->get(i);
576     MG_TETRA* tet2=no2->get_lien_tetra()->get(j);
577     if (tet1==tet2)
578     {
579     MG_NOEUD *nn1,*nn2,*nn3,*nn4;
580     if (tet1->get_noeud1()==no1) nn1=tet1->get_noeud1();
581     if (tet1->get_noeud2()==no1) nn1=tet1->get_noeud2();
582     if (tet1->get_noeud3()==no1) nn1=tet1->get_noeud3();
583     if (tet1->get_noeud4()==no1) nn1=tet1->get_noeud4();
584     if (tet1->get_noeud1()==no2) nn2=tet1->get_noeud1();
585     if (tet1->get_noeud2()==no2) nn2=tet1->get_noeud2();
586     if (tet1->get_noeud3()==no2) nn2=tet1->get_noeud3();
587     if (tet1->get_noeud4()==no2) nn2=tet1->get_noeud4();
588     if (tet1->get_noeud1()!=nn1)
589     if (tet1->get_noeud1()!=nn2) nn3=tet1->get_noeud1();
590     if (tet1->get_noeud2()!=nn1)
591     if (tet1->get_noeud2()!=nn2) nn3=tet1->get_noeud2();
592     if (tet1->get_noeud3()!=nn1)
593     if (tet1->get_noeud3()!=nn2) nn3=tet1->get_noeud3();
594     if (tet1->get_noeud4()!=nn1)
595     if (tet1->get_noeud4()!=nn2) nn3=tet1->get_noeud4();
596     if (tet1->get_noeud1()!=nn1)
597     if (tet1->get_noeud1()!=nn2)
598     if (tet1->get_noeud1()!=nn3) nn4=tet1->get_noeud1();
599     if (tet1->get_noeud2()!=nn1)
600     if (tet1->get_noeud2()!=nn2)
601     if (tet1->get_noeud2()!=nn3) nn4=tet1->get_noeud2();
602     if (tet1->get_noeud3()!=nn1)
603     if (tet1->get_noeud3()!=nn2)
604     if (tet1->get_noeud3()!=nn3) nn4=tet1->get_noeud3();
605     if (tet1->get_noeud4()!=nn1)
606     if (tet1->get_noeud4()!=nn2)
607     if (tet1->get_noeud4()!=nn3) nn4=tet1->get_noeud4();
608     double angle=get_angle_par_noeud<MG_NOEUD*>(nn1,nn2,nn3,nn1,nn2,nn4);
609     if (angle>M_PI) angle=2.*M_PI-angle;
610     angletot=angletot+angle;
611    
612     }
613     }
614     if ((angletot<2*M_PI-0.001) || (angletot>2*M_PI+0.001))
615     {
616 francois 805 char mess[255];
617 francois 425 sprintf(mess," segment %lu angle matiere autour de %.2lf",seg->get_id(),angletot); affiche(mess);
618     nbsegincorrect++;
619     }
620     else nbsegcorrect++;
621     }
622 francois 805 }
623 francois 514
624    
625 francois 577 void MAILLEUR_ANALYSE::compare_maillage_carte_isotrope(FCT_TAILLE* carte,char *nom,double *tab)
626 francois 514 {
627     MG_GESTIONNAIRE *gest=mai->get_gestionnaire();
628 francois 577 LISTE_MG_TETRA::iterator ittet;
629     double vol=0;
630     for (MG_TETRA* tet=mai->get_premier_tetra(ittet);tet!=NULL;tet=mai->get_suivant_tetra(ittet))
631     vol=vol+carte->calcul_volume_tetra_metrique(tet);
632     tab[0]=vol;
633     tab[1]=mai->get_nb_mg_tetra();
634     tab[2]=(tab[1]-tab[0])/tab[0]*100;
635 francois 514 MG_SOLUTION *sol=new MG_SOLUTION(mai,4,nom,1,"Comparaison_taille",MAGIC::ENTITE_SOLUTION::ENTITE_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
636     gest->ajouter_mg_solution(sol);
637     sol->change_legende(0,(char*)"taille_reelle");
638     sol->change_legende(1,(char*)"taille_carte");
639     sol->change_legende(2,(char*)"erreur");
640     sol->change_legende(3,(char*)"erreur absolue");
641     LISTE_MG_NOEUD::iterator it;
642     int i=0;
643     for (MG_NOEUD* no=mai->get_premier_noeud(it);no!=NULL;no=mai->get_suivant_noeud(it))
644     {
645     double *xyz=no->get_coord();
646     double tenseur[9];
647 francois 532 carte->evaluer(xyz,tenseur);
648 francois 514 double cartetaille=1./sqrt(tenseur[0]);
649     int nbseg=no->get_lien_segment()->get_nb();
650     double taille=0.;
651 francois 532 int nbreelseg=0;
652 francois 514 for (int j=0;j<nbseg;j++)
653 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++;}
654     taille=taille/nbreelseg;
655 francois 514 double erreur=(taille-cartetaille)*100./cartetaille;
656     sol->ecrire(taille,i,0);
657     sol->ecrire(cartetaille,i,1);
658     sol->ecrire(erreur,i,2);
659     sol->ecrire(fabs(erreur),i,3);
660     i++;
661     }
662     }