ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/microstructure/src/mstruct_outils.cpp
Revision: 926
Committed: Tue May 1 20:38:42 2018 UTC (7 years ago) by couturad
File size: 29690 byte(s)
Log Message:
Ajout de nouvelles fonctionnalités de representation et d'analyse des VES
*Modification du CMakeLists.txt de microstructure

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     {
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,0,i,r)*fabs(det);
78     val=val+w*valeur;
79     volume=volume+w*fabs(det);
80     }
81     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     else ok=false;
159 couturad 919 }
160 couturad 926 if(ok)
161     {
162     double axe_forme[3];
163     MG_CG_INFO_VCT_DOUBLE* info_axe = (MG_CG_INFO_VCT_DOUBLE*)forme->get_mgcg_info("AXE");
164     std::vector<double> vct_double = info_axe->get_vct_valeur();
165     axe_forme[0]=vct_double[0];
166     axe_forme[1]=vct_double[1];
167     axe_forme[2]=vct_double[2];
168     double *tenseur = new double[6];
169     tenseur[0] = axe_forme[0]*axe_forme[0];
170     tenseur[1] = axe_forme[1]*axe_forme[1];
171     tenseur[2] = axe_forme[2]*axe_forme[2];
172     tenseur[3] = axe_forme[0]*axe_forme[1];
173     tenseur[4] = axe_forme[0]*axe_forme[2];
174     tenseur[5] = axe_forme[1]*axe_forme[2];
175     vector_tenseur.push_back(tenseur);
176     }
177 couturad 919 }
178 couturad 926 long nb_forme=vector_tenseur.size();
179     std::vector<double*>::iterator it_tenseur;
180     for(it_tenseur=vector_tenseur.begin();it_tenseur!=vector_tenseur.end();it_tenseur++)
181 couturad 919 {
182 couturad 926 double* tenseur = *it_tenseur;
183     for(int i=0;i<6;i++)
184 couturad 919 {
185 couturad 926 moyenne[i]+=tenseur[i];
186     if(tenseur[i]<min[i]) min[i]=tenseur[i];
187     if(tenseur[i]>max[i]) max[i]=tenseur[i];
188 couturad 919 }
189     }
190 couturad 926 for(int i=0;i<6;i++)
191 couturad 919 {
192 couturad 926 moyenne[i]=moyenne[i]/nb_forme;
193     }
194     for(it_tenseur=vector_tenseur.begin();it_tenseur!=vector_tenseur.end();it_tenseur++)
195     {
196     double* tenseur = *it_tenseur;
197     for(int i=0;i<6;i++)
198 couturad 919 {
199 couturad 926 ecart_type[i]+=(tenseur[i]-moyenne[i])*(tenseur[i]-moyenne[i]);
200 couturad 919 }
201 couturad 926 delete [] *it_tenseur;
202 couturad 919 }
203 couturad 926 for(int i=0;i<6;i++)
204     {
205     ecart_type[i]=ecart_type[i]*(1.0/(nb_forme-1.0));
206     ecart_type[i]=sqrt(ecart_type[i]);
207     }
208     return OK;
209 couturad 919 }
210    
211 couturad 926 int OUTILS::statistiques_CAO(BOITE_3D boite_VES,
212     TPL_MAP_ENTITE< MG_CG_FORME* >& tpl_map_forme,
213     long int& nb_forme,
214     long int& nb_volume,
215     double& volume_forme_moyenne,
216     double& volume_forme_ecart_type,
217     double& volume_forme_min,
218     double& volume_forme_max,
219     OT_HISTOGRAMME& volume_forme_histogramme,
220     double& volume_total,
221     double& fraction_volumique,
222     double precision)
223 couturad 919 {
224 couturad 926 double volume_ves = boite_VES.get_volume();
225     nb_forme=0;
226     nb_volume=0;
227     volume_forme_moyenne=0;
228     volume_forme_ecart_type=0;
229     volume_forme_min=numeric_limits< double >::max();
230     volume_forme_max=numeric_limits< double >::min();
231     volume_total=0;
232     std::vector<double> vector_volume;
233    
234     TPL_MAP_ENTITE<MG_CG_FORME*>::ITERATEUR it_forme;
235     for(MG_CG_FORME* forme=tpl_map_forme.get_premier(it_forme);forme!=NULL;forme=tpl_map_forme.get_suivant(it_forme))
236 couturad 919 {
237 couturad 926 TPL_MAP_ENTITE<MG_VOLUME*> tpl_map_volume;
238     if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::VOLUME)
239 couturad 919 {
240 couturad 926 MG_CG_FORME_VOLUME* forme_volume = (MG_CG_FORME_VOLUME*)forme;
241     tpl_map_volume.ajouter(forme_volume->get_mg_volume());
242 couturad 919 }
243 couturad 926 else if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::MULTI_VOLUME)
244     {
245     MG_CG_FORME_MULTI_VOLUME* forme_multi_volume = (MG_CG_FORME_MULTI_VOLUME*)forme;
246     std::map<long,MG_VOLUME*>::iterator it_volume;
247     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))
248     {
249     tpl_map_volume.ajouter(vol);
250     }
251     }
252     double volume_forme=0;
253     nb_volume += tpl_map_volume.get_nb();
254     TPL_MAP_ENTITE<MG_VOLUME*>::ITERATEUR it_mg_volume;
255     for(MG_VOLUME* vol=tpl_map_volume.get_premier(it_mg_volume);vol!=NULL;vol=tpl_map_volume.get_suivant(it_mg_volume))
256     {
257     volume_forme += OT_GEOMETRIE::get_volume(vol,precision);
258     }
259     vector_volume.push_back(volume_forme);
260 couturad 919 }
261 couturad 926 std::vector<double>::iterator it_data;
262     for(it_data=vector_volume.begin();it_data!=vector_volume.end();it_data++)
263 couturad 919 {
264 couturad 926 double volume = *it_data;
265     volume_total+=volume;
266     if(volume<volume_forme_min) volume_forme_min=volume;
267     if(volume>volume_forme_max) volume_forme_max=volume;
268 couturad 919 }
269 couturad 926 for(it_data=vector_volume.begin();it_data!=vector_volume.end();it_data++)
270     {
271     double volume = *it_data;
272     volume_forme_histogramme.ajouter_valeur(volume,volume/volume_total);
273     }
274     nb_forme=vector_volume.size();
275     volume_forme_moyenne=volume_total/nb_forme;
276     if(nb_forme>1)
277     {
278     for(it_data=vector_volume.begin();it_data!=vector_volume.end();it_data++)
279     {
280     double volume = *it_data;
281     volume_forme_ecart_type+=(volume-volume_forme_moyenne)*(volume-volume_forme_moyenne);
282     }
283     volume_forme_ecart_type=volume_forme_ecart_type*(1.0/(nb_forme-1.0));
284     volume_forme_ecart_type=sqrt(volume_forme_ecart_type);
285     }
286     else volume_forme_ecart_type = nan("");
287     fraction_volumique = volume_total/volume_ves;
288 couturad 919 }
289    
290 couturad 926 int OUTILS::statistiques_mg_maillage(BOITE_3D boite_VES,
291     MG_MAILLAGE* mai,
292     TPL_MAP_ENTITE<MG_CG_FORME*> &tpl_map_forme,
293     long &nb_element_2D,long &nb_element_3D,
294     double& qualite_moyenne_2D, double& qualite_ecart_type_2D, double& qualite_min_2D, double& qualite_max_2D, OT_HISTOGRAMME& histogramme_qualite_2D,
295     double& qualite_moyenne_3D, double& qualite_ecart_type_3D, double& qualite_min_3D, double& qualite_max_3D, OT_HISTOGRAMME& histogramme_qualite_3D,
296     double& taille_moyenne_2D, double& taille_ecart_type_2D, double& taille_min_2D, double& taille_max_2D, OT_HISTOGRAMME& histogramme_taille_2D,
297     double& taille_moyenne_3D, double& taille_ecart_type_3D, double& taille_min_3D, double& taille_max_3D, OT_HISTOGRAMME& histogramme_taille_3D,
298     double& volume, double& fraction_volumique)
299 couturad 919 {
300 couturad 926 double volume_ves = boite_VES.get_volume();
301     nb_element_2D=0;nb_element_3D=0;
302     qualite_moyenne_2D=0;qualite_ecart_type_2D=0;qualite_min_2D=numeric_limits< double >::max();;qualite_max_2D=numeric_limits< double >::min();
303     qualite_moyenne_3D=0;qualite_ecart_type_3D=0;qualite_min_3D=numeric_limits< double >::max();;qualite_max_3D=numeric_limits< double >::min();
304     taille_moyenne_2D=0;taille_ecart_type_2D=0;taille_min_2D=numeric_limits< double >::max();;taille_max_2D=numeric_limits< double >::min();
305     taille_moyenne_3D=0;taille_ecart_type_3D=0;taille_min_3D=numeric_limits< double >::max();;taille_max_3D=numeric_limits< double >::min();
306     volume=0;fraction_volumique=0;
307    
308     TPL_MAP_ENTITE<MG_ELEMENT_MAILLAGE*> tpl_map_element_maill_3D;
309     TPL_MAP_ENTITE<MG_ELEMENT_MAILLAGE*> tpl_map_element_maill_2D;
310     TPL_MAP_ENTITE<MG_CG_FORME*>::ITERATEUR it_forme;
311     for(MG_CG_FORME* forme=tpl_map_forme.get_premier(it_forme);forme!=NULL;forme=tpl_map_forme.get_suivant(it_forme))
312 couturad 919 {
313 couturad 926 if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::VOLUME)
314 couturad 919 {
315 couturad 926 MG_CG_FORME_VOLUME* forme_volume = (MG_CG_FORME_VOLUME*)forme;
316     TPL_SET<MG_ELEMENT_MAILLAGE*> *tmp_maill = forme_volume->get_mg_volume()->get_lien_maillage();
317     TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it_maill;
318     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);
319     TPL_MAP_ENTITE<MG_FACE*> tpl_map_face;
320     OT_GEOMETRIE::get_map_mg_face_sous_jacent(forme_volume->get_mg_volume(),tpl_map_face);
321     TPL_MAP_ENTITE<MG_FACE*>::ITERATEUR it_face;
322     for(MG_FACE* face=tpl_map_face.get_premier(it_face);face!=NULL;face=tpl_map_face.get_suivant(it_face))
323     {
324     tmp_maill = face->get_lien_maillage();
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_2D.ajouter(ele);
326     }
327     double volume_i;
328     OT_VECTEUR_3D cdm;
329     OT_MATRICE_3D inertieglobale;
330     OT_MATRICE_3D inertiecdm;
331     forme_volume->get_mg_volume()->get_propriete_massique(mai,volume_i,cdm,inertieglobale,inertiecdm);
332     volume+=volume_i;
333 couturad 919 }
334 couturad 926 else if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::MULTI_VOLUME)
335     {
336     MG_CG_FORME_MULTI_VOLUME* forme_multi_volume = (MG_CG_FORME_MULTI_VOLUME*)forme;
337     std::map<long,MG_VOLUME*>::iterator it_volume;
338     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))
339     {
340     TPL_SET<MG_ELEMENT_MAILLAGE*> *tmp_maill = vol->get_lien_maillage();
341     TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it_maill;
342     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);
343     TPL_MAP_ENTITE<MG_FACE*> tpl_map_face;
344     OT_GEOMETRIE::get_map_mg_face_sous_jacent(vol,tpl_map_face);
345     TPL_MAP_ENTITE<MG_FACE*>::ITERATEUR it_face;
346     for(MG_FACE* face=tpl_map_face.get_premier(it_face);face!=NULL;face=tpl_map_face.get_suivant(it_face))
347     {
348     tmp_maill = face->get_lien_maillage();
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_2D.ajouter(ele);
350     }
351     double volume_i;
352     OT_VECTEUR_3D cdm;
353     OT_MATRICE_3D inertieglobale;
354     OT_MATRICE_3D inertiecdm;
355     vol->get_propriete_massique(mai,volume_i,cdm,inertieglobale,inertiecdm);
356     volume+=volume_i;
357     }
358     }
359 couturad 919 }
360 couturad 926 fraction_volumique = volume/volume_ves;
361     std::vector<double> vector_qualite_3D;
362     std::vector<double> vector_taille_3D;
363     LISTE_MG_TETRA::iterator it_tetra;
364     for(MG_TETRA* tet=mai->get_premier_tetra(it_tetra);tet!=NULL;tet=mai->get_suivant_tetra(it_tetra))
365 couturad 919 {
366 couturad 926 if(tpl_map_element_maill_3D.existe(tet))
367     {
368     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());
369     vector_qualite_3D.push_back(qual);
370     if(qual<qualite_min_3D) qualite_min_3D=qual;
371     if(qual>qualite_max_3D) qualite_max_3D=qual;
372     qualite_moyenne_3D+=qual;
373     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());
374     vector_taille_3D.push_back(taille);
375     if(taille<taille_min_3D) taille_min_3D=taille;
376     if(taille>taille_max_3D) taille_max_3D=taille;
377     taille_moyenne_3D+=taille;
378     nb_element_3D++;
379     }
380     }
381     qualite_moyenne_3D = qualite_moyenne_3D/nb_element_3D;
382     taille_moyenne_3D = taille_moyenne_3D/nb_element_3D;
383     std::vector<double>::iterator it_qual_3D;
384    
385     for(it_qual_3D=vector_qualite_3D.begin();it_qual_3D!=vector_qualite_3D.end();it_qual_3D++)
386     {
387     double qual = *it_qual_3D;
388     qualite_ecart_type_3D+=(qual-qualite_moyenne_3D)*(qual-qualite_moyenne_3D);
389     histogramme_qualite_3D.ajouter_valeur(qual,1.0/nb_element_3D);
390 couturad 919 }
391 couturad 926 qualite_ecart_type_3D=qualite_ecart_type_3D*(1.0/(nb_element_3D-1.0));
392     qualite_ecart_type_3D=sqrt(qualite_ecart_type_3D);
393    
394     std::vector<double>::iterator it_taille_3D;
395     for(it_taille_3D=vector_taille_3D.begin();it_taille_3D!=vector_taille_3D.end();it_taille_3D++)
396 couturad 919 {
397 couturad 926 double taille = *it_taille_3D;
398     taille_ecart_type_3D=(taille-taille_moyenne_3D)*(taille-taille_moyenne_3D);
399     histogramme_taille_3D.ajouter_valeur(taille,1.0/nb_element_3D);
400     }
401     taille_ecart_type_3D=taille_ecart_type_3D*(1.0/(nb_element_3D-1.0));
402     taille_ecart_type_3D=sqrt(taille_ecart_type_3D);
403    
404    
405    
406     std::vector<double> vector_qualite_2D;
407     std::vector<double> vector_taille_2D;
408     LISTE_MG_TRIANGLE::iterator it_tri;
409     for(MG_TRIANGLE* tri=mai->get_premier_triangle(it_tri);tri!=NULL;tri=mai->get_suivant_triangle(it_tri))
410     {
411     if(tpl_map_element_maill_2D.existe(tri))
412 couturad 919 {
413 couturad 926 double qual=OPERATEUR::qualite_triangle(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord(),tri->get_noeud3()->get_coord());
414     vector_qualite_2D.push_back(qual);
415     if(qual<qualite_min_2D) qualite_min_2D=qual;
416     if(qual>qualite_max_2D) qualite_max_2D=qual;
417     qualite_moyenne_2D+=qual;
418     double taille=OPERATEUR::taille_triangle(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord(),tri->get_noeud3()->get_coord());
419     vector_taille_2D.push_back(taille);
420     if(taille<taille_min_2D) taille_min_2D=taille;
421     if(taille>taille_max_2D) taille_max_2D=taille;
422     taille_moyenne_2D+=taille;
423     nb_element_2D++;
424 couturad 919 }
425 couturad 926 }
426     qualite_moyenne_2D = qualite_moyenne_2D/nb_element_2D;
427     taille_moyenne_2D = taille_moyenne_2D/nb_element_2D;
428     std::vector<double>::iterator it_qual_2D;
429    
430     for(it_qual_2D=vector_qualite_2D.begin();it_qual_2D!=vector_qualite_2D.end();it_qual_2D++)
431     {
432     double qual = *it_qual_2D;
433     qualite_ecart_type_2D+=(qual-qualite_moyenne_2D)*(qual-qualite_moyenne_2D);
434     histogramme_qualite_2D.ajouter_valeur(qual,1.0/nb_element_2D);
435 couturad 919 }
436 couturad 926 qualite_ecart_type_2D=qualite_ecart_type_2D*(1.0/(nb_element_2D-1.0));
437     qualite_ecart_type_2D=sqrt(qualite_ecart_type_2D);
438    
439     std::vector<double>::iterator it_taille_2D;
440     for(it_taille_2D=vector_taille_2D.begin();it_taille_2D!=vector_taille_2D.end();it_taille_2D++)
441 couturad 919 {
442 couturad 926 double taille = *it_taille_2D;
443     taille_ecart_type_2D=(taille-taille_moyenne_2D)*(taille-taille_moyenne_2D);
444     histogramme_taille_2D.ajouter_valeur(taille,1.0/nb_element_2D);
445 couturad 919 }
446 couturad 926 taille_ecart_type_2D=taille_ecart_type_2D*(1.0/(nb_element_2D-1.0));
447     taille_ecart_type_2D=sqrt(taille_ecart_type_2D);
448    
449 couturad 919 }
450    
451    
452 couturad 926 int OUTILS::statistiques_fem_maillage(BOITE_3D boite_VES,
453     FEM_MAILLAGE* fem,
454     BOITE_3D *boite_analyse,
455     TPL_MAP_ENTITE< MG_CG_FORME* >& tpl_map_forme,
456     long int& nb_element_2D, long int& nb_element_3D,
457     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,
458     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,
459     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,
460     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,
461     double& distortion_moyenne_2D, double& distortion_ecart_type_2D, double& distortion_min_2D, double& distortion_max_2D, OT_HISTOGRAMME& histogramme_distortion_2D,
462     double& distortion_moyenne_3D, double& distortion_ecart_type_3D, double& distortion_min_3D, double& distortion_max_3D, OT_HISTOGRAMME& histogramme_distortion_3D,
463     double& volume,
464     double& fraction_volumique)
465 couturad 919 {
466 couturad 926 double volume_ves=boite_VES.get_volume();
467     nb_element_2D=0;nb_element_3D=0;
468     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();
469     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();
470     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();
471     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();
472     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();
473     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();
474     volume=0;
475     fraction_volumique=0;
476    
477     TPL_MAP_ENTITE<FEM_ELEMENT3*> tpl_map_element_maill_3D;
478     TPL_MAP_ENTITE<FEM_ELEMENT2*> tpl_map_element_maill_2D;
479     TPL_MAP_ENTITE<MG_CG_FORME*>::ITERATEUR it_forme;
480     for(MG_CG_FORME* forme=tpl_map_forme.get_premier(it_forme);forme!=NULL;forme=tpl_map_forme.get_suivant(it_forme))
481 couturad 919 {
482 couturad 926 if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::VOLUME)
483 couturad 919 {
484 couturad 926 MG_CG_FORME_VOLUME* forme_volume = (MG_CG_FORME_VOLUME*)forme;
485     TPL_LISTE_ENTITE<FEM_ELEMENT_MAILLAGE*> *tmp_maill = forme_volume->get_mg_volume()->get_lien_fem_maillage();
486     long tmp_nb = tmp_maill->get_nb();
487     for(long i=0;i<tmp_nb;i++)
488 couturad 919 {
489 couturad 926 FEM_ELEMENT3* ele3 = (FEM_ELEMENT3*)tmp_maill->get(i);
490     BOITE_3D b3d_ele = ele3->get_boite_3D();
491     if(boite_analyse!=NULL)
492     if(!boite_analyse->contient(b3d_ele.get_xcentre(),b3d_ele.get_ycentre(),b3d_ele.get_zcentre())) continue;
493     tpl_map_element_maill_3D.ajouter(ele3);
494 couturad 919 }
495 couturad 926 TPL_MAP_ENTITE<MG_FACE*> tpl_map_face;
496     OT_GEOMETRIE::get_map_mg_face_sous_jacent(forme_volume->get_mg_volume(),tpl_map_face);
497     TPL_MAP_ENTITE<MG_FACE*>::ITERATEUR it_face;
498     for(MG_FACE* face=tpl_map_face.get_premier(it_face);face!=NULL;face=tpl_map_face.get_suivant(it_face))
499 couturad 919 {
500 couturad 926 tmp_maill = face->get_lien_fem_maillage();
501     tmp_nb = tmp_maill->get_nb();
502     for(long i=0;i<tmp_nb;i++)
503 couturad 919 {
504 couturad 926 FEM_ELEMENT2* ele2 = (FEM_ELEMENT2*)tmp_maill->get(i);
505     BOITE_3D b3d_ele = ele2->get_boite_3D();
506     if(boite_analyse!=NULL)
507     if(!boite_analyse->contient(b3d_ele.get_xcentre(),b3d_ele.get_ycentre(),b3d_ele.get_zcentre())) continue;
508     tpl_map_element_maill_2D.ajouter(ele2);
509 couturad 919 }
510     }
511     }
512 couturad 926 else if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::MULTI_VOLUME)
513     {
514     MG_CG_FORME_MULTI_VOLUME* forme_multi_volume = (MG_CG_FORME_MULTI_VOLUME*)forme;
515     std::map<long,MG_VOLUME*>::iterator it_volume;
516     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))
517     {
518     TPL_LISTE_ENTITE<FEM_ELEMENT_MAILLAGE*> *tmp_maill = vol->get_lien_fem_maillage();
519     long tmp_nb = tmp_maill->get_nb();
520     for(long i=0;i<tmp_nb;i++)
521     {
522     FEM_ELEMENT3* ele3 = (FEM_ELEMENT3*)tmp_maill->get(i);
523     BOITE_3D b3d_ele = ele3->get_boite_3D();
524     if(boite_analyse!=NULL)
525     if(!boite_analyse->contient(b3d_ele.get_xcentre(),b3d_ele.get_ycentre(),b3d_ele.get_zcentre())) continue;
526     tpl_map_element_maill_3D.ajouter(ele3);
527     }
528     TPL_MAP_ENTITE<MG_FACE*> tpl_map_face;
529     OT_GEOMETRIE::get_map_mg_face_sous_jacent(vol,tpl_map_face);
530     TPL_MAP_ENTITE<MG_FACE*>::ITERATEUR it_face;
531     for(MG_FACE* face=tpl_map_face.get_premier(it_face);face!=NULL;face=tpl_map_face.get_suivant(it_face))
532     {
533     tmp_maill = face->get_lien_fem_maillage();
534     tmp_nb = tmp_maill->get_nb();
535     for(long i=0;i<tmp_nb;i++)
536     {
537     FEM_ELEMENT2* ele2 = (FEM_ELEMENT2*)tmp_maill->get(i);
538     BOITE_3D b3d_ele = ele2->get_boite_3D();
539     if(boite_analyse!=NULL)
540     if(!boite_analyse->contient(b3d_ele.get_xcentre(),b3d_ele.get_ycentre(),b3d_ele.get_zcentre())) continue;
541     tpl_map_element_maill_2D.ajouter(ele2);
542     }
543     }
544     }
545     }
546 couturad 919 }
547 couturad 926 nb_element_3D=tpl_map_element_maill_3D.get_nb();
548     nb_element_2D=tpl_map_element_maill_2D.get_nb();
549     std::vector<double> vector_jacobien_min_2D;
550     std::vector<double> vector_jacobien_max_2D;
551     std::vector<double> vector_jacobien_min_3D;
552     std::vector<double> vector_jacobien_max_3D;
553     std::vector<double> vector_distortion_2D;
554     std::vector<double> vector_distortion_3D;
555     FEM_MAILLAGE_QUADRATIQUE_OUTILS ot;
556     TPL_MAP_ENTITE<FEM_ELEMENT3*>::ITERATEUR it_ele3;
557     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))
558 couturad 919 {
559 couturad 926 int nb_pt_gauss=ele3->get_nb_pt_gauss(ele3->get_degremax_fonction_interpolation());
560     for (int r=0;r<nb_pt_gauss;r++)
561     {
562     int degre=ele3->get_degre_gauss(nb_pt_gauss);
563     double w;
564     double uvw[3];
565     ele3->get_pt_gauss(degre,r,w,uvw);
566     double jac[9];
567     int l,c;
568     double det=ele3->get_jacobien(jac,uvw,l,c);
569     volume+=w*det;
570     }
571     double jmin=ot.get_jmin(ele3);
572     double jmax=ot.get_jmin(ele3);
573     double distortion = ot.get_distorsion2(ele3);
574     vector_jacobien_min_3D.push_back(jmin);
575     vector_jacobien_max_3D.push_back(jmax);
576     vector_distortion_3D.push_back(distortion);
577     jacobien_min_moyenne_3D+=jmin;
578     jacobien_max_moyenne_3D+=jmax;
579     distortion_moyenne_3D+=distortion;
580 couturad 919 }
581 couturad 926 fraction_volumique=volume/volume_ves;
582     jacobien_min_moyenne_3D=jacobien_min_moyenne_3D/nb_element_3D;
583     jacobien_max_moyenne_3D=jacobien_max_moyenne_3D/nb_element_3D;
584     distortion_moyenne_3D=distortion_moyenne_3D/nb_element_3D;
585     std::vector<double>::iterator it_jac_min_3D=vector_jacobien_min_3D.begin();
586     std::vector<double>::iterator it_jac_max_3D=vector_jacobien_max_3D.begin();
587     std::vector<double>::iterator it_dist_3D=vector_distortion_3D.begin();
588     for(long i=0;i<nb_element_3D;i++)
589 couturad 919 {
590 couturad 926 double jmin=*it_jac_min_3D;
591     double jmax=*it_jac_max_3D;
592     double dist=*it_dist_3D;
593     if(jmin<jacobien_min_min_3D)jacobien_min_min_3D=jmin;
594     if(jmin>jacobien_min_max_3D)jacobien_min_max_3D=jmin;
595     if(jmax<jacobien_max_min_3D)jacobien_max_min_3D=jmax;
596     if(jmax>jacobien_max_max_3D)jacobien_max_max_3D=jmax;
597     if(dist<distortion_min_3D)distortion_min_3D=dist;
598     if(dist>distortion_max_3D)distortion_max_3D=dist;
599     jacobien_min_ecart_type_3D=(jmin-jacobien_min_moyenne_3D)*(jmin-jacobien_min_moyenne_3D);
600     jacobien_max_ecart_type_3D=(jmax-jacobien_max_moyenne_3D)*(jmax-jacobien_max_moyenne_3D);
601     distortion_ecart_type_3D=(dist-distortion_moyenne_3D)*(dist-distortion_moyenne_3D);
602     histogramme_jacobien_min_3D.ajouter_valeur(jmin,1.0/nb_element_3D);
603     histogramme_jacobien_max_3D.ajouter_valeur(jmax,1.0/nb_element_3D);
604     histogramme_distortion_3D.ajouter_valeur(dist,1.0/nb_element_3D);
605     it_jac_min_3D++;
606     it_jac_max_3D++;
607     it_dist_3D++;
608 couturad 919 }
609 couturad 926 jacobien_min_ecart_type_3D=sqrt(jacobien_min_ecart_type_3D*(1.0/(nb_element_3D-1.0)));
610     jacobien_max_ecart_type_3D=sqrt(jacobien_max_ecart_type_3D*(1.0/(nb_element_3D-1.0)));
611     distortion_ecart_type_3D=sqrt(distortion_ecart_type_3D*(1.0/(nb_element_3D-1.0)));
612    
613     TPL_MAP_ENTITE<FEM_ELEMENT2*>::ITERATEUR it_ele2;
614     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))
615 couturad 919 {
616 couturad 926 double jmin=ot.get_jmin(ele2);
617     double jmax=ot.get_jmin(ele2);
618     double distortion = ot.get_distorsion2(ele2);
619     vector_jacobien_min_2D.push_back(jmin);
620     vector_jacobien_max_2D.push_back(jmax);
621     vector_distortion_2D.push_back(distortion);
622     jacobien_min_moyenne_2D+=jmin;
623     jacobien_max_moyenne_2D+=jmax;
624     distortion_moyenne_2D+=distortion;
625 couturad 919 }
626 couturad 926 jacobien_min_moyenne_2D=jacobien_min_moyenne_2D/nb_element_2D;
627     jacobien_max_moyenne_2D=jacobien_max_moyenne_2D/nb_element_2D;
628     distortion_moyenne_2D=distortion_moyenne_2D/nb_element_2D;
629     std::vector<double>::iterator it_jac_min_2D=vector_jacobien_min_2D.begin();
630     std::vector<double>::iterator it_jac_max_2D=vector_jacobien_max_2D.begin();
631     std::vector<double>::iterator it_dist_2D=vector_distortion_2D.begin();
632     for(long i=0;i<nb_element_2D;i++)
633 couturad 919 {
634 couturad 926 double jmin=*it_jac_min_2D;
635     double jmax=*it_jac_max_2D;
636     double dist=*it_dist_2D;
637     if(jmin<jacobien_min_min_2D)jacobien_min_min_2D=jmin;
638     if(jmin>jacobien_min_max_2D)jacobien_min_max_2D=jmin;
639     if(jmax<jacobien_max_min_2D)jacobien_max_min_2D=jmax;
640     if(jmax>jacobien_max_max_2D)jacobien_max_max_2D=jmax;
641     if(dist<distortion_min_2D)distortion_min_2D=dist;
642     if(dist>distortion_max_2D)distortion_max_2D=dist;
643     jacobien_min_ecart_type_2D=(jmin-jacobien_min_moyenne_2D)*(jmin-jacobien_min_moyenne_2D);
644     jacobien_max_ecart_type_2D=(jmax-jacobien_max_moyenne_2D)*(jmax-jacobien_max_moyenne_2D);
645     distortion_ecart_type_2D=(dist-distortion_moyenne_2D)*(dist-distortion_moyenne_2D);
646     histogramme_jacobien_min_2D.ajouter_valeur(jmin,1.0/nb_element_2D);
647     histogramme_jacobien_max_2D.ajouter_valeur(jmax,1.0/nb_element_2D);
648     histogramme_distortion_2D.ajouter_valeur(dist,1.0/nb_element_2D);
649     it_jac_min_2D++;
650     it_jac_max_2D++;
651     it_dist_2D++;
652 couturad 919 }
653 couturad 926 jacobien_min_ecart_type_2D=sqrt(jacobien_min_ecart_type_2D*(1.0/(nb_element_2D-1.0)));
654     jacobien_max_ecart_type_2D=sqrt(jacobien_max_ecart_type_2D*(1.0/(nb_element_2D-1.0)));
655     distortion_ecart_type_2D=sqrt(distortion_ecart_type_2D*(1.0/(nb_element_2D-1.0)));
656 couturad 919 }
657    
658    
659    
660    
661    
662    
663    
664    
665    
666