ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/microstructure/src/mstruct_outils.cpp
Revision: 934
Committed: Mon May 28 19:32:31 2018 UTC (6 years, 11 months ago) by couturad
File size: 31333 byte(s)
Log Message:
Correction bugs dans MICROSTRUCTURE
Ajout d'une fonction de vérification du decalage entre deux faces pour la detection des couches minces

File Contents

# User Rev Content
1 couturad 919 #include "mstruct_outils.h"
2 couturad 926 #include "mg_cg_groupe_forme.h"
3 couturad 919 #include "ot_geometrie.h"
4     #include "mg_maillage.h"
5     #include "fem_maillage.h"
6     #include "mg_cg_forme.h"
7     #include "mg_cg_forme_multi_volume.h"
8     #include "mg_cg_forme_volume.h"
9 couturad 926 #include "ot_tenseur.h"
10     #include "mg_cg_info.h"
11     #include "fem_maillage_quadratique_outils.h"
12 couturad 919 #include <math.h>
13    
14     using namespace MICROSTRUCTURE;
15    
16 couturad 926 int OUTILS::statistiques_champ_volumique(FEM_SOLUTION* fem_sol,
17     MG_CG_GROUPE_FORME* groupe_forme,
18     BOITE_3D* boite_analyse,
19     double* moyenne,
20     double* ecart_type,
21     double* min,
22     double* max,
23     OT_HISTOGRAMME** tab_histogramme)
24 couturad 919 {
25 couturad 926 int nb_champ = fem_sol->get_nb_champ();
26     double volume_total=0.0;
27     for(int i=0;i<nb_champ;i++)
28 couturad 919 {
29 couturad 926 moyenne[i]=0.0;
30     ecart_type[i]=0.0;
31     min[i]=numeric_limits< double >::max();
32     max[i]=numeric_limits< double >::min();
33     }
34     std::vector<double*> vector_champ;
35     std::vector<double> vector_volume;
36     LISTE_FEM_ELEMENT3::iterator it;
37     FEM_MAILLAGE* fem_mail = fem_sol->get_maillage();
38     long k=0;
39     for (FEM_ELEMENT3* ele=fem_mail->get_premier_element3(it);ele!=NULL;ele=fem_mail->get_suivant_element3(it))
40     {
41     bool ok=false;
42     if(groupe_forme!=NULL)
43 couturad 919 {
44 couturad 926 TPL_MAP_ENTITE<MG_VOLUME*> tpl_map_volume = groupe_forme->get_tpl_map_volume();
45     TPL_MAP_ENTITE<MG_VOLUME*>::ITERATEUR it_volume;
46     MG_VOLUME* mg_vol=tpl_map_volume.get_premier(it_volume);
47     while(ok==false && mg_vol!=NULL)
48     {
49     if(ele->get_lien_topologie()==mg_vol) ok=true;
50     else mg_vol=tpl_map_volume.get_suivant(it_volume);
51     }
52     } else ok=true;
53     if(boite_analyse!=NULL && ok==true)
54     {
55     BOITE_3D b3d = ele->get_boite_3D();
56     if(!boite_analyse->contient(b3d.get_xcentre(),b3d.get_ycentre(),b3d.get_zcentre())) ok=false;
57 couturad 919 }
58 couturad 926 if(ok==true)
59     {
60     double *champ=new double[nb_champ];
61     double volume;
62     for(int i=0;i<nb_champ;i++)
63     {
64     double unite=fem_mail->get_mg_geometrie()->get_valeur_unite();
65     double val=0.;
66     volume=0;
67     int nbgauss=fem_sol->get_nb_gauss(k);
68     for (int r=0;r<nbgauss;r++)
69 couturad 934 {
70     int degre=ele->get_degre_gauss(nbgauss);
71     double w;
72     double uvw[3];
73     ele->get_pt_gauss(degre,r,w,uvw);
74     double jac[9];
75     int l,c;
76     double det=ele->get_jacobien(jac,uvw,l,c,unite);
77     double valeur=fem_sol->lire(k,i,0,r)*fabs(det);
78     val=val+w*valeur;
79     volume=volume+w*fabs(det);
80     }
81 couturad 926 champ[i] = val;
82     }
83     vector_champ.push_back(champ);
84     vector_volume.push_back(volume);
85     volume_total+=volume;
86     }
87     k++;
88     }
89     long nb_echantillions=vector_champ.size();
90     std::vector<double*>::iterator it_vector_champ;
91     std::vector<double>::iterator it_vector_volume;
92     for(it_vector_champ=vector_champ.begin(),it_vector_volume=vector_volume.begin();it_vector_champ!=vector_champ.end();it_vector_champ++,it_vector_volume++)
93     {
94     double *champ = *it_vector_champ;
95     double volume = *it_vector_volume;
96     for(int i=0;i<nb_champ;i++)
97     {
98     double champ_volumique = champ[i]/volume;
99     moyenne[i]+=champ_volumique;
100     if(champ_volumique<min[i]) min[i]=champ_volumique;
101     if(champ_volumique>max[i]) max[i]=champ_volumique;
102     tab_histogramme[i]->ajouter_valeur(champ[i],volume/volume_total);
103     }
104 couturad 919 }
105 couturad 926 for(int i=0;i<nb_champ;i++)
106 couturad 919 {
107 couturad 926 moyenne[i]=moyenne[i]/nb_echantillions;
108     }
109     for(it_vector_champ=vector_champ.begin(),it_vector_volume=vector_volume.begin();it_vector_champ!=vector_champ.end();it_vector_champ++,it_vector_volume++)
110     {
111     double *champ = *it_vector_champ;
112     double volume = *it_vector_volume;
113     for(int i=0;i<nb_champ;i++)
114 couturad 919 {
115 couturad 926 double champ_volumique = champ[i]/volume;
116     ecart_type[i]+=(champ_volumique-moyenne[i])*(champ_volumique-moyenne[i]);
117 couturad 919 }
118     }
119 couturad 926 for(it_vector_champ=vector_champ.begin();it_vector_champ!=vector_champ.end();it_vector_champ++) delete [] *it_vector_champ;
120     for(int i=0;i<nb_champ;i++)
121     {
122     ecart_type[i]=ecart_type[i]*(1.0/(nb_echantillions-1.0));
123     ecart_type[i]=sqrt(ecart_type[i]);
124     }
125     return OK;
126 couturad 919 }
127    
128 couturad 926 int OUTILS::statistiques_tenseur_orientation(MG_CG_GROUPE_FORME* groupe_forme,
129     BOITE_3D* boite_analyse,
130     double* moyenne,
131     double* ecart_type,
132     double* min,
133     double* max)
134 couturad 919 {
135 couturad 926 for(int i=0;i<6;i++)
136 couturad 919 {
137 couturad 926 moyenne[i]=0.0;
138     ecart_type[i]=0.0;
139     min[i]=numeric_limits< double >::max();
140     max[i]=numeric_limits< double >::min();
141 couturad 919 }
142 couturad 926 std::vector<double*> vector_tenseur;
143     std::map<long,MG_CG_FORME*>::iterator it_forme;
144     for(MG_CG_FORME* forme=groupe_forme->get_premiere_mgcg_forme(it_forme);forme!=NULL;forme=groupe_forme->get_suivante_mgcg_forme(it_forme))
145 couturad 919 {
146 couturad 926 BOITE_3D boite_forme = forme->get_boite_3D();
147     bool ok=true;
148     if(boite_analyse!=NULL)
149 couturad 919 {
150 couturad 926 if(boite_analyse->contient(boite_forme.get_xmin(),boite_forme.get_ymin(),boite_forme.get_zmin())) ok=true;
151     else if(boite_analyse->contient(boite_forme.get_xmax(),boite_forme.get_ymin(),boite_forme.get_zmin())) ok=true;
152     else if(boite_analyse->contient(boite_forme.get_xmin(),boite_forme.get_ymax(),boite_forme.get_zmin())) ok=true;
153     else if(boite_analyse->contient(boite_forme.get_xmax(),boite_forme.get_ymax(),boite_forme.get_zmin())) ok=true;
154     else if(boite_analyse->contient(boite_forme.get_xmin(),boite_forme.get_ymin(),boite_forme.get_zmax())) ok=true;
155     else if(boite_analyse->contient(boite_forme.get_xmax(),boite_forme.get_ymin(),boite_forme.get_zmax())) ok=true;
156     else if(boite_analyse->contient(boite_forme.get_xmin(),boite_forme.get_ymax(),boite_forme.get_zmax())) ok=true;
157     else if(boite_analyse->contient(boite_forme.get_xmax(),boite_forme.get_ymax(),boite_forme.get_zmax())) ok=true;
158 couturad 930 else if(boite_forme.contient(boite_analyse->get_xmin(),boite_analyse->get_ymin(),boite_analyse->get_zmin())) ok=true;
159     else if(boite_forme.contient(boite_analyse->get_xmax(),boite_analyse->get_ymin(),boite_analyse->get_zmin())) ok=true;
160     else if(boite_forme.contient(boite_analyse->get_xmin(),boite_analyse->get_ymax(),boite_analyse->get_zmin())) ok=true;
161     else if(boite_forme.contient(boite_analyse->get_xmax(),boite_analyse->get_ymax(),boite_analyse->get_zmin())) ok=true;
162     else if(boite_forme.contient(boite_analyse->get_xmin(),boite_analyse->get_ymin(),boite_analyse->get_zmax())) ok=true;
163     else if(boite_forme.contient(boite_analyse->get_xmax(),boite_analyse->get_ymin(),boite_analyse->get_zmax())) ok=true;
164     else if(boite_forme.contient(boite_analyse->get_xmin(),boite_analyse->get_ymax(),boite_analyse->get_zmax())) ok=true;
165     else if(boite_forme.contient(boite_analyse->get_xmax(),boite_analyse->get_ymax(),boite_analyse->get_zmax())) ok=true;
166 couturad 926 else ok=false;
167 couturad 919 }
168 couturad 926 if(ok)
169     {
170     double axe_forme[3];
171     MG_CG_INFO_VCT_DOUBLE* info_axe = (MG_CG_INFO_VCT_DOUBLE*)forme->get_mgcg_info("AXE");
172     std::vector<double> vct_double = info_axe->get_vct_valeur();
173     axe_forme[0]=vct_double[0];
174     axe_forme[1]=vct_double[1];
175     axe_forme[2]=vct_double[2];
176     double *tenseur = new double[6];
177     tenseur[0] = axe_forme[0]*axe_forme[0];
178     tenseur[1] = axe_forme[1]*axe_forme[1];
179     tenseur[2] = axe_forme[2]*axe_forme[2];
180     tenseur[3] = axe_forme[0]*axe_forme[1];
181     tenseur[4] = axe_forme[0]*axe_forme[2];
182     tenseur[5] = axe_forme[1]*axe_forme[2];
183     vector_tenseur.push_back(tenseur);
184     }
185 couturad 919 }
186 couturad 926 long nb_forme=vector_tenseur.size();
187     std::vector<double*>::iterator it_tenseur;
188     for(it_tenseur=vector_tenseur.begin();it_tenseur!=vector_tenseur.end();it_tenseur++)
189 couturad 919 {
190 couturad 926 double* tenseur = *it_tenseur;
191     for(int i=0;i<6;i++)
192 couturad 919 {
193 couturad 926 moyenne[i]+=tenseur[i];
194     if(tenseur[i]<min[i]) min[i]=tenseur[i];
195     if(tenseur[i]>max[i]) max[i]=tenseur[i];
196 couturad 919 }
197     }
198 couturad 926 for(int i=0;i<6;i++)
199 couturad 919 {
200 couturad 926 moyenne[i]=moyenne[i]/nb_forme;
201     }
202     for(it_tenseur=vector_tenseur.begin();it_tenseur!=vector_tenseur.end();it_tenseur++)
203     {
204     double* tenseur = *it_tenseur;
205     for(int i=0;i<6;i++)
206 couturad 919 {
207 couturad 926 ecart_type[i]+=(tenseur[i]-moyenne[i])*(tenseur[i]-moyenne[i]);
208 couturad 919 }
209 couturad 926 delete [] *it_tenseur;
210 couturad 919 }
211 couturad 926 for(int i=0;i<6;i++)
212     {
213     ecart_type[i]=ecart_type[i]*(1.0/(nb_forme-1.0));
214     ecart_type[i]=sqrt(ecart_type[i]);
215     }
216     return OK;
217 couturad 919 }
218    
219 couturad 926 int OUTILS::statistiques_CAO(BOITE_3D boite_VES,
220     TPL_MAP_ENTITE< MG_CG_FORME* >& tpl_map_forme,
221     long int& nb_forme,
222     long int& nb_volume,
223     double& volume_forme_moyenne,
224     double& volume_forme_ecart_type,
225     double& volume_forme_min,
226     double& volume_forme_max,
227     OT_HISTOGRAMME& volume_forme_histogramme,
228     double& volume_total,
229     double& fraction_volumique,
230     double precision)
231 couturad 919 {
232 couturad 926 double volume_ves = boite_VES.get_volume();
233     nb_forme=0;
234     nb_volume=0;
235     volume_forme_moyenne=0;
236     volume_forme_ecart_type=0;
237     volume_forme_min=numeric_limits< double >::max();
238     volume_forme_max=numeric_limits< double >::min();
239     volume_total=0;
240     std::vector<double> vector_volume;
241    
242     TPL_MAP_ENTITE<MG_CG_FORME*>::ITERATEUR it_forme;
243     for(MG_CG_FORME* forme=tpl_map_forme.get_premier(it_forme);forme!=NULL;forme=tpl_map_forme.get_suivant(it_forme))
244 couturad 919 {
245 couturad 926 TPL_MAP_ENTITE<MG_VOLUME*> tpl_map_volume;
246     if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::VOLUME)
247 couturad 919 {
248 couturad 926 MG_CG_FORME_VOLUME* forme_volume = (MG_CG_FORME_VOLUME*)forme;
249     tpl_map_volume.ajouter(forme_volume->get_mg_volume());
250 couturad 919 }
251 couturad 926 else if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::MULTI_VOLUME)
252     {
253     MG_CG_FORME_MULTI_VOLUME* forme_multi_volume = (MG_CG_FORME_MULTI_VOLUME*)forme;
254     std::map<long,MG_VOLUME*>::iterator it_volume;
255     for(MG_VOLUME* vol=forme_multi_volume->get_premier_mg_volume(it_volume);vol!=NULL;vol=forme_multi_volume->get_suivant_mg_volume(it_volume))
256     {
257     tpl_map_volume.ajouter(vol);
258     }
259     }
260     double volume_forme=0;
261     nb_volume += tpl_map_volume.get_nb();
262     TPL_MAP_ENTITE<MG_VOLUME*>::ITERATEUR it_mg_volume;
263     for(MG_VOLUME* vol=tpl_map_volume.get_premier(it_mg_volume);vol!=NULL;vol=tpl_map_volume.get_suivant(it_mg_volume))
264     {
265     volume_forme += OT_GEOMETRIE::get_volume(vol,precision);
266     }
267     vector_volume.push_back(volume_forme);
268 couturad 919 }
269 couturad 926 std::vector<double>::iterator it_data;
270     for(it_data=vector_volume.begin();it_data!=vector_volume.end();it_data++)
271 couturad 919 {
272 couturad 926 double volume = *it_data;
273     volume_total+=volume;
274     if(volume<volume_forme_min) volume_forme_min=volume;
275     if(volume>volume_forme_max) volume_forme_max=volume;
276 couturad 919 }
277 couturad 926 for(it_data=vector_volume.begin();it_data!=vector_volume.end();it_data++)
278     {
279     double volume = *it_data;
280     volume_forme_histogramme.ajouter_valeur(volume,volume/volume_total);
281     }
282     nb_forme=vector_volume.size();
283     volume_forme_moyenne=volume_total/nb_forme;
284     if(nb_forme>1)
285     {
286     for(it_data=vector_volume.begin();it_data!=vector_volume.end();it_data++)
287     {
288     double volume = *it_data;
289     volume_forme_ecart_type+=(volume-volume_forme_moyenne)*(volume-volume_forme_moyenne);
290     }
291     volume_forme_ecart_type=volume_forme_ecart_type*(1.0/(nb_forme-1.0));
292     volume_forme_ecart_type=sqrt(volume_forme_ecart_type);
293     }
294     else volume_forme_ecart_type = nan("");
295     fraction_volumique = volume_total/volume_ves;
296 couturad 919 }
297    
298 couturad 926 int OUTILS::statistiques_mg_maillage(BOITE_3D boite_VES,
299     MG_MAILLAGE* mai,
300     TPL_MAP_ENTITE<MG_CG_FORME*> &tpl_map_forme,
301     long &nb_element_2D,long &nb_element_3D,
302     double& qualite_moyenne_2D, double& qualite_ecart_type_2D, double& qualite_min_2D, double& qualite_max_2D, OT_HISTOGRAMME& histogramme_qualite_2D,
303     double& qualite_moyenne_3D, double& qualite_ecart_type_3D, double& qualite_min_3D, double& qualite_max_3D, OT_HISTOGRAMME& histogramme_qualite_3D,
304     double& taille_moyenne_2D, double& taille_ecart_type_2D, double& taille_min_2D, double& taille_max_2D, OT_HISTOGRAMME& histogramme_taille_2D,
305     double& taille_moyenne_3D, double& taille_ecart_type_3D, double& taille_min_3D, double& taille_max_3D, OT_HISTOGRAMME& histogramme_taille_3D,
306     double& volume, double& fraction_volumique)
307 couturad 919 {
308 couturad 926 double volume_ves = boite_VES.get_volume();
309     nb_element_2D=0;nb_element_3D=0;
310     qualite_moyenne_2D=0;qualite_ecart_type_2D=0;qualite_min_2D=numeric_limits< double >::max();;qualite_max_2D=numeric_limits< double >::min();
311     qualite_moyenne_3D=0;qualite_ecart_type_3D=0;qualite_min_3D=numeric_limits< double >::max();;qualite_max_3D=numeric_limits< double >::min();
312     taille_moyenne_2D=0;taille_ecart_type_2D=0;taille_min_2D=numeric_limits< double >::max();;taille_max_2D=numeric_limits< double >::min();
313     taille_moyenne_3D=0;taille_ecart_type_3D=0;taille_min_3D=numeric_limits< double >::max();;taille_max_3D=numeric_limits< double >::min();
314     volume=0;fraction_volumique=0;
315     TPL_MAP_ENTITE<MG_ELEMENT_MAILLAGE*> tpl_map_element_maill_3D;
316     TPL_MAP_ENTITE<MG_ELEMENT_MAILLAGE*> tpl_map_element_maill_2D;
317     TPL_MAP_ENTITE<MG_CG_FORME*>::ITERATEUR it_forme;
318     for(MG_CG_FORME* forme=tpl_map_forme.get_premier(it_forme);forme!=NULL;forme=tpl_map_forme.get_suivant(it_forme))
319 couturad 919 {
320 couturad 926 if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::VOLUME)
321 couturad 919 {
322 couturad 926 MG_CG_FORME_VOLUME* forme_volume = (MG_CG_FORME_VOLUME*)forme;
323     TPL_SET<MG_ELEMENT_MAILLAGE*> *tmp_maill = forme_volume->get_mg_volume()->get_lien_maillage();
324     TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it_maill;
325     for(MG_ELEMENT_MAILLAGE* ele=tmp_maill->get_premier(it_maill);ele!=NULL;ele=tmp_maill->get_suivant(it_maill)) tpl_map_element_maill_3D.ajouter(ele);
326     TPL_MAP_ENTITE<MG_FACE*> tpl_map_face;
327     OT_GEOMETRIE::get_map_mg_face_sous_jacent(forme_volume->get_mg_volume(),tpl_map_face);
328     TPL_MAP_ENTITE<MG_FACE*>::ITERATEUR it_face;
329     for(MG_FACE* face=tpl_map_face.get_premier(it_face);face!=NULL;face=tpl_map_face.get_suivant(it_face))
330     {
331     tmp_maill = face->get_lien_maillage();
332     for(MG_ELEMENT_MAILLAGE* ele=tmp_maill->get_premier(it_maill);ele!=NULL;ele=tmp_maill->get_suivant(it_maill)) tpl_map_element_maill_2D.ajouter(ele);
333     }
334     double volume_i;
335     OT_VECTEUR_3D cdm;
336     OT_MATRICE_3D inertieglobale;
337     OT_MATRICE_3D inertiecdm;
338     forme_volume->get_mg_volume()->get_propriete_massique(mai,volume_i,cdm,inertieglobale,inertiecdm);
339     volume+=volume_i;
340 couturad 919 }
341 couturad 926 else if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::MULTI_VOLUME)
342     {
343     MG_CG_FORME_MULTI_VOLUME* forme_multi_volume = (MG_CG_FORME_MULTI_VOLUME*)forme;
344     std::map<long,MG_VOLUME*>::iterator it_volume;
345     for(MG_VOLUME* vol=forme_multi_volume->get_premier_mg_volume(it_volume);vol!=NULL;vol=forme_multi_volume->get_suivant_mg_volume(it_volume))
346     {
347     TPL_SET<MG_ELEMENT_MAILLAGE*> *tmp_maill = vol->get_lien_maillage();
348     TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it_maill;
349     for(MG_ELEMENT_MAILLAGE* ele=tmp_maill->get_premier(it_maill);ele!=NULL;ele=tmp_maill->get_suivant(it_maill)) tpl_map_element_maill_3D.ajouter(ele);
350     TPL_MAP_ENTITE<MG_FACE*> tpl_map_face;
351     OT_GEOMETRIE::get_map_mg_face_sous_jacent(vol,tpl_map_face);
352     TPL_MAP_ENTITE<MG_FACE*>::ITERATEUR it_face;
353     for(MG_FACE* face=tpl_map_face.get_premier(it_face);face!=NULL;face=tpl_map_face.get_suivant(it_face))
354     {
355     tmp_maill = face->get_lien_maillage();
356     for(MG_ELEMENT_MAILLAGE* ele=tmp_maill->get_premier(it_maill);ele!=NULL;ele=tmp_maill->get_suivant(it_maill)) tpl_map_element_maill_2D.ajouter(ele);
357     }
358     double volume_i;
359     OT_VECTEUR_3D cdm;
360     OT_MATRICE_3D inertieglobale;
361     OT_MATRICE_3D inertiecdm;
362     vol->get_propriete_massique(mai,volume_i,cdm,inertieglobale,inertiecdm);
363     volume+=volume_i;
364     }
365     }
366 couturad 919 }
367 couturad 926 fraction_volumique = volume/volume_ves;
368     std::vector<double> vector_qualite_3D;
369     std::vector<double> vector_taille_3D;
370     LISTE_MG_TETRA::iterator it_tetra;
371     for(MG_TETRA* tet=mai->get_premier_tetra(it_tetra);tet!=NULL;tet=mai->get_suivant_tetra(it_tetra))
372 couturad 919 {
373 couturad 926 if(tpl_map_element_maill_3D.existe(tet))
374     {
375     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());
376     vector_qualite_3D.push_back(qual);
377     if(qual<qualite_min_3D) qualite_min_3D=qual;
378     if(qual>qualite_max_3D) qualite_max_3D=qual;
379     qualite_moyenne_3D+=qual;
380     double taille=OPERATEUR::taille_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord());
381     vector_taille_3D.push_back(taille);
382     if(taille<taille_min_3D) taille_min_3D=taille;
383     if(taille>taille_max_3D) taille_max_3D=taille;
384     taille_moyenne_3D+=taille;
385     nb_element_3D++;
386     }
387     }
388     qualite_moyenne_3D = qualite_moyenne_3D/nb_element_3D;
389     taille_moyenne_3D = taille_moyenne_3D/nb_element_3D;
390     std::vector<double>::iterator it_qual_3D;
391    
392     for(it_qual_3D=vector_qualite_3D.begin();it_qual_3D!=vector_qualite_3D.end();it_qual_3D++)
393     {
394     double qual = *it_qual_3D;
395     qualite_ecart_type_3D+=(qual-qualite_moyenne_3D)*(qual-qualite_moyenne_3D);
396     histogramme_qualite_3D.ajouter_valeur(qual,1.0/nb_element_3D);
397 couturad 919 }
398 couturad 926 qualite_ecart_type_3D=qualite_ecart_type_3D*(1.0/(nb_element_3D-1.0));
399     qualite_ecart_type_3D=sqrt(qualite_ecart_type_3D);
400    
401     std::vector<double>::iterator it_taille_3D;
402     for(it_taille_3D=vector_taille_3D.begin();it_taille_3D!=vector_taille_3D.end();it_taille_3D++)
403 couturad 919 {
404 couturad 926 double taille = *it_taille_3D;
405     taille_ecart_type_3D=(taille-taille_moyenne_3D)*(taille-taille_moyenne_3D);
406     histogramme_taille_3D.ajouter_valeur(taille,1.0/nb_element_3D);
407     }
408     taille_ecart_type_3D=taille_ecart_type_3D*(1.0/(nb_element_3D-1.0));
409     taille_ecart_type_3D=sqrt(taille_ecart_type_3D);
410    
411    
412    
413     std::vector<double> vector_qualite_2D;
414     std::vector<double> vector_taille_2D;
415     LISTE_MG_TRIANGLE::iterator it_tri;
416     for(MG_TRIANGLE* tri=mai->get_premier_triangle(it_tri);tri!=NULL;tri=mai->get_suivant_triangle(it_tri))
417     {
418     if(tpl_map_element_maill_2D.existe(tri))
419 couturad 919 {
420 couturad 926 double qual=OPERATEUR::qualite_triangle(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord(),tri->get_noeud3()->get_coord());
421     vector_qualite_2D.push_back(qual);
422     if(qual<qualite_min_2D) qualite_min_2D=qual;
423     if(qual>qualite_max_2D) qualite_max_2D=qual;
424     qualite_moyenne_2D+=qual;
425     double taille=OPERATEUR::taille_triangle(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord(),tri->get_noeud3()->get_coord());
426     vector_taille_2D.push_back(taille);
427     if(taille<taille_min_2D) taille_min_2D=taille;
428     if(taille>taille_max_2D) taille_max_2D=taille;
429     taille_moyenne_2D+=taille;
430     nb_element_2D++;
431 couturad 919 }
432 couturad 926 }
433     qualite_moyenne_2D = qualite_moyenne_2D/nb_element_2D;
434     taille_moyenne_2D = taille_moyenne_2D/nb_element_2D;
435     std::vector<double>::iterator it_qual_2D;
436    
437     for(it_qual_2D=vector_qualite_2D.begin();it_qual_2D!=vector_qualite_2D.end();it_qual_2D++)
438     {
439     double qual = *it_qual_2D;
440     qualite_ecart_type_2D+=(qual-qualite_moyenne_2D)*(qual-qualite_moyenne_2D);
441     histogramme_qualite_2D.ajouter_valeur(qual,1.0/nb_element_2D);
442 couturad 919 }
443 couturad 926 qualite_ecart_type_2D=qualite_ecart_type_2D*(1.0/(nb_element_2D-1.0));
444     qualite_ecart_type_2D=sqrt(qualite_ecart_type_2D);
445    
446     std::vector<double>::iterator it_taille_2D;
447     for(it_taille_2D=vector_taille_2D.begin();it_taille_2D!=vector_taille_2D.end();it_taille_2D++)
448 couturad 919 {
449 couturad 926 double taille = *it_taille_2D;
450     taille_ecart_type_2D=(taille-taille_moyenne_2D)*(taille-taille_moyenne_2D);
451     histogramme_taille_2D.ajouter_valeur(taille,1.0/nb_element_2D);
452 couturad 919 }
453 couturad 926 taille_ecart_type_2D=taille_ecart_type_2D*(1.0/(nb_element_2D-1.0));
454     taille_ecart_type_2D=sqrt(taille_ecart_type_2D);
455    
456 couturad 919 }
457    
458    
459 couturad 931 int OUTILS::statistiques_fem_maillage(FEM_MAILLAGE* fem,
460 couturad 926 BOITE_3D *boite_analyse,
461     TPL_MAP_ENTITE< MG_CG_FORME* >& tpl_map_forme,
462     long int& nb_element_2D, long int& nb_element_3D,
463     double& jacobien_min_moyenne_2D, double& jacobien_min_ecart_type_2D, double& jacobien_min_min_2D, double& jacobien_min_max_2D, OT_HISTOGRAMME& histogramme_jacobien_min_2D,
464     double& jacobien_max_moyenne_2D, double& jacobien_max_ecart_type_2D, double& jacobien_max_min_2D, double& jacobien_max_max_2D, OT_HISTOGRAMME& histogramme_jacobien_max_2D,
465     double& jacobien_min_moyenne_3D, double& jacobien_min_ecart_type_3D, double& jacobien_min_min_3D, double& jacobien_min_max_3D, OT_HISTOGRAMME& histogramme_jacobien_min_3D,
466     double& jacobien_max_moyenne_3D, double& jacobien_max_ecart_type_3D, double& jacobien_max_min_3D, double& jacobien_max_max_3D, OT_HISTOGRAMME& histogramme_jacobien_max_3D,
467     double& distortion_moyenne_2D, double& distortion_ecart_type_2D, double& distortion_min_2D, double& distortion_max_2D, OT_HISTOGRAMME& histogramme_distortion_2D,
468     double& distortion_moyenne_3D, double& distortion_ecart_type_3D, double& distortion_min_3D, double& distortion_max_3D, OT_HISTOGRAMME& histogramme_distortion_3D,
469     double& volume,
470     double& fraction_volumique)
471 couturad 919 {
472 couturad 931 double volume_ves=0;
473 couturad 926 nb_element_2D=0;nb_element_3D=0;
474     jacobien_min_moyenne_2D=0.0; jacobien_min_ecart_type_2D=0.0; jacobien_min_min_2D=numeric_limits< double >::max(); jacobien_min_max_2D=numeric_limits< double >::min();
475     jacobien_max_moyenne_2D=0.0; jacobien_max_ecart_type_2D=0.0; jacobien_max_min_2D=numeric_limits< double >::max(); jacobien_max_max_2D=numeric_limits< double >::min();
476     jacobien_min_moyenne_3D=0.0; jacobien_min_ecart_type_3D=0.0; jacobien_min_min_3D=numeric_limits< double >::max(); jacobien_min_max_3D=numeric_limits< double >::min();
477     jacobien_max_moyenne_3D=0.0; jacobien_max_ecart_type_3D=0.0; jacobien_max_min_3D=numeric_limits< double >::max(); jacobien_max_max_3D=numeric_limits< double >::min();
478     distortion_moyenne_2D=0.0; distortion_ecart_type_2D=0.0; distortion_min_2D=numeric_limits< double >::max(); distortion_max_2D=numeric_limits< double >::min();
479     distortion_moyenne_3D=0.0; distortion_ecart_type_3D=0.0; distortion_min_3D=numeric_limits< double >::max(); distortion_max_3D=numeric_limits< double >::min();
480     volume=0;
481     fraction_volumique=0;
482    
483 couturad 931 LISTE_FEM_ELEMENT3::iterator it_element3;
484     for(FEM_ELEMENT3* ele3=fem->get_premier_element3(it_element3);ele3!=NULL;ele3=fem->get_suivant_element3(it_element3))
485     {
486     BOITE_3D b3d_ele = ele3->get_boite_3D();
487     if(boite_analyse!=NULL)
488     if(!boite_analyse->contient(b3d_ele.get_xcentre(),b3d_ele.get_ycentre(),b3d_ele.get_zcentre())) continue;
489     int nb_pt_gauss=ele3->get_nb_pt_gauss(ele3->get_degremax_fonction_interpolation());
490     for (int r=0;r<nb_pt_gauss;r++)
491     {
492     int degre=ele3->get_degre_gauss(nb_pt_gauss);
493     double w;
494     double uvw[3];
495     ele3->get_pt_gauss(degre,r,w,uvw);
496     double jac[9];
497     int l,c;
498     double det=ele3->get_jacobien(jac,uvw,l,c);
499     volume_ves+=w*det;
500     }
501     }
502 couturad 926 TPL_MAP_ENTITE<FEM_ELEMENT3*> tpl_map_element_maill_3D;
503     TPL_MAP_ENTITE<FEM_ELEMENT2*> tpl_map_element_maill_2D;
504     TPL_MAP_ENTITE<MG_CG_FORME*>::ITERATEUR it_forme;
505     for(MG_CG_FORME* forme=tpl_map_forme.get_premier(it_forme);forme!=NULL;forme=tpl_map_forme.get_suivant(it_forme))
506 couturad 919 {
507 couturad 926 if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::VOLUME)
508 couturad 919 {
509 couturad 926 MG_CG_FORME_VOLUME* forme_volume = (MG_CG_FORME_VOLUME*)forme;
510     TPL_LISTE_ENTITE<FEM_ELEMENT_MAILLAGE*> *tmp_maill = forme_volume->get_mg_volume()->get_lien_fem_maillage();
511     long tmp_nb = tmp_maill->get_nb();
512     for(long i=0;i<tmp_nb;i++)
513 couturad 919 {
514 couturad 926 FEM_ELEMENT3* ele3 = (FEM_ELEMENT3*)tmp_maill->get(i);
515     BOITE_3D b3d_ele = ele3->get_boite_3D();
516     if(boite_analyse!=NULL)
517     if(!boite_analyse->contient(b3d_ele.get_xcentre(),b3d_ele.get_ycentre(),b3d_ele.get_zcentre())) continue;
518     tpl_map_element_maill_3D.ajouter(ele3);
519 couturad 919 }
520 couturad 926 TPL_MAP_ENTITE<MG_FACE*> tpl_map_face;
521     OT_GEOMETRIE::get_map_mg_face_sous_jacent(forme_volume->get_mg_volume(),tpl_map_face);
522     TPL_MAP_ENTITE<MG_FACE*>::ITERATEUR it_face;
523     for(MG_FACE* face=tpl_map_face.get_premier(it_face);face!=NULL;face=tpl_map_face.get_suivant(it_face))
524 couturad 919 {
525 couturad 926 tmp_maill = face->get_lien_fem_maillage();
526     tmp_nb = tmp_maill->get_nb();
527     for(long i=0;i<tmp_nb;i++)
528 couturad 919 {
529 couturad 926 FEM_ELEMENT2* ele2 = (FEM_ELEMENT2*)tmp_maill->get(i);
530     BOITE_3D b3d_ele = ele2->get_boite_3D();
531     if(boite_analyse!=NULL)
532     if(!boite_analyse->contient(b3d_ele.get_xcentre(),b3d_ele.get_ycentre(),b3d_ele.get_zcentre())) continue;
533     tpl_map_element_maill_2D.ajouter(ele2);
534 couturad 919 }
535     }
536     }
537 couturad 926 else if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::MULTI_VOLUME)
538     {
539     MG_CG_FORME_MULTI_VOLUME* forme_multi_volume = (MG_CG_FORME_MULTI_VOLUME*)forme;
540     std::map<long,MG_VOLUME*>::iterator it_volume;
541     for(MG_VOLUME* vol=forme_multi_volume->get_premier_mg_volume(it_volume);vol!=NULL;vol=forme_multi_volume->get_suivant_mg_volume(it_volume))
542     {
543     TPL_LISTE_ENTITE<FEM_ELEMENT_MAILLAGE*> *tmp_maill = vol->get_lien_fem_maillage();
544     long tmp_nb = tmp_maill->get_nb();
545     for(long i=0;i<tmp_nb;i++)
546     {
547     FEM_ELEMENT3* ele3 = (FEM_ELEMENT3*)tmp_maill->get(i);
548     BOITE_3D b3d_ele = ele3->get_boite_3D();
549     if(boite_analyse!=NULL)
550     if(!boite_analyse->contient(b3d_ele.get_xcentre(),b3d_ele.get_ycentre(),b3d_ele.get_zcentre())) continue;
551     tpl_map_element_maill_3D.ajouter(ele3);
552     }
553     TPL_MAP_ENTITE<MG_FACE*> tpl_map_face;
554     OT_GEOMETRIE::get_map_mg_face_sous_jacent(vol,tpl_map_face);
555     TPL_MAP_ENTITE<MG_FACE*>::ITERATEUR it_face;
556     for(MG_FACE* face=tpl_map_face.get_premier(it_face);face!=NULL;face=tpl_map_face.get_suivant(it_face))
557     {
558     tmp_maill = face->get_lien_fem_maillage();
559     tmp_nb = tmp_maill->get_nb();
560     for(long i=0;i<tmp_nb;i++)
561     {
562     FEM_ELEMENT2* ele2 = (FEM_ELEMENT2*)tmp_maill->get(i);
563     BOITE_3D b3d_ele = ele2->get_boite_3D();
564     if(boite_analyse!=NULL)
565     if(!boite_analyse->contient(b3d_ele.get_xcentre(),b3d_ele.get_ycentre(),b3d_ele.get_zcentre())) continue;
566     tpl_map_element_maill_2D.ajouter(ele2);
567     }
568     }
569     }
570     }
571 couturad 919 }
572 couturad 926 nb_element_3D=tpl_map_element_maill_3D.get_nb();
573     nb_element_2D=tpl_map_element_maill_2D.get_nb();
574     std::vector<double> vector_jacobien_min_2D;
575     std::vector<double> vector_jacobien_max_2D;
576     std::vector<double> vector_jacobien_min_3D;
577     std::vector<double> vector_jacobien_max_3D;
578     std::vector<double> vector_distortion_2D;
579     std::vector<double> vector_distortion_3D;
580     FEM_MAILLAGE_QUADRATIQUE_OUTILS ot;
581     TPL_MAP_ENTITE<FEM_ELEMENT3*>::ITERATEUR it_ele3;
582     for(FEM_ELEMENT3* ele3=tpl_map_element_maill_3D.get_premier(it_ele3);ele3!=NULL;ele3=tpl_map_element_maill_3D.get_suivant(it_ele3))
583 couturad 919 {
584 couturad 926 int nb_pt_gauss=ele3->get_nb_pt_gauss(ele3->get_degremax_fonction_interpolation());
585     for (int r=0;r<nb_pt_gauss;r++)
586     {
587     int degre=ele3->get_degre_gauss(nb_pt_gauss);
588     double w;
589     double uvw[3];
590     ele3->get_pt_gauss(degre,r,w,uvw);
591     double jac[9];
592     int l,c;
593     double det=ele3->get_jacobien(jac,uvw,l,c);
594     volume+=w*det;
595     }
596     double jmin=ot.get_jmin(ele3);
597     double jmax=ot.get_jmin(ele3);
598     double distortion = ot.get_distorsion2(ele3);
599     vector_jacobien_min_3D.push_back(jmin);
600     vector_jacobien_max_3D.push_back(jmax);
601     vector_distortion_3D.push_back(distortion);
602     jacobien_min_moyenne_3D+=jmin;
603     jacobien_max_moyenne_3D+=jmax;
604     distortion_moyenne_3D+=distortion;
605 couturad 919 }
606 couturad 926 fraction_volumique=volume/volume_ves;
607     jacobien_min_moyenne_3D=jacobien_min_moyenne_3D/nb_element_3D;
608     jacobien_max_moyenne_3D=jacobien_max_moyenne_3D/nb_element_3D;
609     distortion_moyenne_3D=distortion_moyenne_3D/nb_element_3D;
610     std::vector<double>::iterator it_jac_min_3D=vector_jacobien_min_3D.begin();
611     std::vector<double>::iterator it_jac_max_3D=vector_jacobien_max_3D.begin();
612     std::vector<double>::iterator it_dist_3D=vector_distortion_3D.begin();
613     for(long i=0;i<nb_element_3D;i++)
614 couturad 919 {
615 couturad 926 double jmin=*it_jac_min_3D;
616     double jmax=*it_jac_max_3D;
617     double dist=*it_dist_3D;
618     if(jmin<jacobien_min_min_3D)jacobien_min_min_3D=jmin;
619     if(jmin>jacobien_min_max_3D)jacobien_min_max_3D=jmin;
620     if(jmax<jacobien_max_min_3D)jacobien_max_min_3D=jmax;
621     if(jmax>jacobien_max_max_3D)jacobien_max_max_3D=jmax;
622     if(dist<distortion_min_3D)distortion_min_3D=dist;
623     if(dist>distortion_max_3D)distortion_max_3D=dist;
624     jacobien_min_ecart_type_3D=(jmin-jacobien_min_moyenne_3D)*(jmin-jacobien_min_moyenne_3D);
625     jacobien_max_ecart_type_3D=(jmax-jacobien_max_moyenne_3D)*(jmax-jacobien_max_moyenne_3D);
626     distortion_ecart_type_3D=(dist-distortion_moyenne_3D)*(dist-distortion_moyenne_3D);
627     histogramme_jacobien_min_3D.ajouter_valeur(jmin,1.0/nb_element_3D);
628     histogramme_jacobien_max_3D.ajouter_valeur(jmax,1.0/nb_element_3D);
629     histogramme_distortion_3D.ajouter_valeur(dist,1.0/nb_element_3D);
630     it_jac_min_3D++;
631     it_jac_max_3D++;
632     it_dist_3D++;
633 couturad 919 }
634 couturad 926 jacobien_min_ecart_type_3D=sqrt(jacobien_min_ecart_type_3D*(1.0/(nb_element_3D-1.0)));
635     jacobien_max_ecart_type_3D=sqrt(jacobien_max_ecart_type_3D*(1.0/(nb_element_3D-1.0)));
636     distortion_ecart_type_3D=sqrt(distortion_ecart_type_3D*(1.0/(nb_element_3D-1.0)));
637    
638     TPL_MAP_ENTITE<FEM_ELEMENT2*>::ITERATEUR it_ele2;
639     for(FEM_ELEMENT2* ele2=tpl_map_element_maill_2D.get_premier(it_ele2);ele2!=NULL;ele2=tpl_map_element_maill_2D.get_suivant(it_ele2))
640 couturad 919 {
641 couturad 926 double jmin=ot.get_jmin(ele2);
642     double jmax=ot.get_jmin(ele2);
643     double distortion = ot.get_distorsion2(ele2);
644     vector_jacobien_min_2D.push_back(jmin);
645     vector_jacobien_max_2D.push_back(jmax);
646     vector_distortion_2D.push_back(distortion);
647     jacobien_min_moyenne_2D+=jmin;
648     jacobien_max_moyenne_2D+=jmax;
649     distortion_moyenne_2D+=distortion;
650 couturad 919 }
651 couturad 926 jacobien_min_moyenne_2D=jacobien_min_moyenne_2D/nb_element_2D;
652     jacobien_max_moyenne_2D=jacobien_max_moyenne_2D/nb_element_2D;
653     distortion_moyenne_2D=distortion_moyenne_2D/nb_element_2D;
654     std::vector<double>::iterator it_jac_min_2D=vector_jacobien_min_2D.begin();
655     std::vector<double>::iterator it_jac_max_2D=vector_jacobien_max_2D.begin();
656     std::vector<double>::iterator it_dist_2D=vector_distortion_2D.begin();
657     for(long i=0;i<nb_element_2D;i++)
658 couturad 919 {
659 couturad 926 double jmin=*it_jac_min_2D;
660     double jmax=*it_jac_max_2D;
661     double dist=*it_dist_2D;
662     if(jmin<jacobien_min_min_2D)jacobien_min_min_2D=jmin;
663     if(jmin>jacobien_min_max_2D)jacobien_min_max_2D=jmin;
664     if(jmax<jacobien_max_min_2D)jacobien_max_min_2D=jmax;
665     if(jmax>jacobien_max_max_2D)jacobien_max_max_2D=jmax;
666     if(dist<distortion_min_2D)distortion_min_2D=dist;
667     if(dist>distortion_max_2D)distortion_max_2D=dist;
668     jacobien_min_ecart_type_2D=(jmin-jacobien_min_moyenne_2D)*(jmin-jacobien_min_moyenne_2D);
669     jacobien_max_ecart_type_2D=(jmax-jacobien_max_moyenne_2D)*(jmax-jacobien_max_moyenne_2D);
670     distortion_ecart_type_2D=(dist-distortion_moyenne_2D)*(dist-distortion_moyenne_2D);
671     histogramme_jacobien_min_2D.ajouter_valeur(jmin,1.0/nb_element_2D);
672     histogramme_jacobien_max_2D.ajouter_valeur(jmax,1.0/nb_element_2D);
673     histogramme_distortion_2D.ajouter_valeur(dist,1.0/nb_element_2D);
674     it_jac_min_2D++;
675     it_jac_max_2D++;
676     it_dist_2D++;
677 couturad 919 }
678 couturad 926 jacobien_min_ecart_type_2D=sqrt(jacobien_min_ecart_type_2D*(1.0/(nb_element_2D-1.0)));
679     jacobien_max_ecart_type_2D=sqrt(jacobien_max_ecart_type_2D*(1.0/(nb_element_2D-1.0)));
680     distortion_ecart_type_2D=sqrt(distortion_ecart_type_2D*(1.0/(nb_element_2D-1.0)));
681 couturad 919 }
682    
683    
684    
685    
686    
687    
688    
689    
690    
691