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