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

Les paramètres dans ~/.magic.

File Contents

# User Rev Content
1 francois 756 //------------------------------------------------------------
2     //------------------------------------------------------------
3     // MAGiC
4     // Jean Christophe Cuilli�re et Vincent FRANCOIS
5     // D�partement de G�nie M�canique - UQTR
6     //------------------------------------------------------------
7     // Le projet MAGIC est un projet de recherche du d�partement
8     // de g�nie m�canique de l'Universit� du Qu�bec �
9     // Trois Rivi�res
10     // Les librairies ne peuvent �tre utilis�es sans l'accord
11     // des auteurs (contact : francois@uqtr.ca)
12     //------------------------------------------------------------
13     //------------------------------------------------------------
14     //
15     // mg_ves.cpp
16     //
17     //------------------------------------------------------------
18     //------------------------------------------------------------
19     // COPYRIGHT 2000
20     // Version du 02/03/2006 � 11H23
21     //------------------------------------------------------------
22     //------------------------------------------------------------
23    
24 couturad 801 #include "ve_ves.h"
25     #include "ve_definition.h"
26 francois 756 #include "gestionversion.h"
27 couturad 801 #include "mg_gestionnaire.h"
28     #include "mg_arbre.h"
29     #include "mg_assemblage.h"
30     #include "mg_primitive.h"
31     #include "mg_primitive_boite.h"
32     #include "mg_primitive_sphere.h"
33     #include "mg_primitive_cylindre.h"
34     #include "mg_primitive_cone.h"
35     #include "mg_primitive_tore.h"
36     #include "mg_primitive_complexe.h"
37     #include "mg_operateur_boolean.h"
38     #include "mg_operateur_boolean_difference.h"
39     #include "mg_operateur_boolean_intersection.h"
40     #include "mg_operateur_boolean_union.h"
41     #include "mg_sommet_noeud.h"
42     #include "fct_taille.h"
43     #include "fct_generateur_frontiere.h"
44     #include "fct_generateur_constante.h"
45     #include "fct_generateur_fichier.h"
46     #include "fct_taille_fem_solution_generateur_constant.h"
47     #include "fct_taille_fem_solution_generateur_echantillon_frontiere.h"
48     #include "fct_taille_fem_solution_generateur_fichier.h"
49     #include "mailleur.h"
50     #include "mailleur0d.h"
51     #include "mailleur1d.h"
52     #include "mailleur2d.h"
53     #include "mailleur3d.h"
54 francois 883 #include "mailleur_fem.h"
55 couturad 801 #include "tpl_map_entite.h"
56     #include "ot_mathematique.h"
57     #include "fem_maillage.h"
58     #include "mgaster.h"
59     #include "mg_export.h"
60     #include "sphere_grille.h"
61     #include <string.h>
62     #include <map>
63     #include <random>
64     #include <stdio.h>
65     #include <iostream>
66     #include <chrono>
67 francois 756
68 couturad 801 VE_VES::VE_VES()
69     {
70     gest = new MG_GESTIONNAIRE;
71     arbre = new MG_ARBRE("Arbre");
72     assemblage = new MG_ASSEMBLAGE("Assemblage");
73     param_ver = new OT_PARAMETRES();
74     gest->ajouter_mg_arbre(arbre);
75     arbre->ajouter_mg_assemblage(assemblage);
76     ini_param_VER(param_ver);
77     dim_VES[0] = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_VER_Dim_longueur_x]);
78     dim_VES[1] = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_VER_Dim_longueur_y]);
79     dim_VES[2] = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_VER_Dim_longueur_z]);
80     }
81 francois 756
82 couturad 801 VE_VES::VE_VES(VE_VES& mdd)
83     {
84 francois 756
85 couturad 801 }
86    
87     VE_VES::~VE_VES()
88 francois 756 {
89 couturad 801 delete assemblage;
90     delete param_ver;
91     delete gest;
92     }
93 francois 756
94 couturad 801 int VE_VES::construire(void)
95     {
96     int type_ver = (int)param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_VER_Type]);
97     switch(type_ver)
98     {
99     case VE_TYPE_VER::HOMOGENE:
100     {
101     generer_VES_homogene();
102     break;
103     }
104     case VE_TYPE_VER::UNE_INCLUSION:
105     {
106     generer_VES_une_inclusion();
107     break;
108     }
109     case VE_TYPE_VER::ALEATOIRE:
110     {
111     generer_VES_aleatoire();
112     break;
113     }
114     }
115 couturad 814 affiche((char*)"Evaluation de la geometrie :\n");
116 couturad 801 evaluer_geometrie();
117     creer_carte();
118     creer_mg_maillage();
119     appliquer_conditions_limites(param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_VER_CL_CHARGEMENT_type]));
120     appliquer_materiaux();
121     creer_fem_maillage();
122     calcul_aster();
123     homogeneisation();
124 couturad 814 // char nom_fichier_magic[1000];
125     // sprintf(nom_fichier_magic,"%s",param_ver->get_nom(VE_PARAM_CHAR[VE_PARAM::VE_IO_nom_fichier_MAGiC]).c_str());
126     // gest->enregistrer(nom_fichier_magic);
127 couturad 801 return 0;
128 francois 756 }
129    
130 couturad 801 MG_GESTIONNAIRE* VE_VES::get_gestionnaire(void)
131     {
132     return gest;
133     }
134 francois 756
135 couturad 801 MG_ARBRE* VE_VES::get_arbre(void)
136     {
137     return arbre;
138     }
139 francois 756
140 couturad 801 MG_ASSEMBLAGE* VE_VES::get_assemblage(void)
141 francois 756 {
142 couturad 801 return assemblage;
143     }
144 francois 756
145 couturad 801 void VE_VES::active_affichage(fonction_affiche* fonc)
146     {
147     fonc_affiche = fonc;
148     affichageactif = 1;
149 francois 756 }
150    
151 couturad 801 void VE_VES::affiche(char* message)
152     {
153     if(affichageactif==1) fonc_affiche(message);
154     }
155 francois 756
156 couturad 801 void VE_VES::lire_param_VES(char* fichier)
157     {
158     param_ver->lire(fichier);
159     dim_VES[0] = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_VER_Dim_longueur_x]);
160     dim_VES[1] = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_VER_Dim_longueur_y]);
161     dim_VES[2] = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_VER_Dim_longueur_z]);
162     }
163 francois 756
164 couturad 801 void VE_VES::ecrire_param_VES(char* fichier)
165 francois 756 {
166 couturad 801 param_ver->enregistrer(fichier);
167     }
168 francois 756
169 couturad 801 int VE_VES::generer_VES_homogene(void)
170     {
171     affiche((char*)"Generation VER Homogene :\n");
172     MG_PRIMITIVE_BOITE *cube = new MG_PRIMITIVE_BOITE(0,0,0,dim_VES[0],dim_VES[1],dim_VES[2]);
173     arbre->ajouter_mg_primitive(cube);
174     cube->construit();
175     assemblage->ajouter_mg_primitive(cube);
176     return 0;
177 francois 756 }
178 couturad 801
179     int VE_VES::generer_VES_une_inclusion(void)
180     {
181     affiche((char*)"Generation VER Une Inclusion :\n");
182     MG_PRIMITIVE_BOITE *cube = new MG_PRIMITIVE_BOITE(0,0,0,dim_VES[0],dim_VES[1],dim_VES[2]);
183     arbre->ajouter_mg_primitive(cube);
184     cube->construit();
185    
186     MG_PRIMITIVE *matrice;
187     MG_PRIMITIVE *inclusion;
188     MG_PRIMITIVE *inclusion_outil;
189    
190    
191     int type_inclusion = (int)param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_INCLU_Type_INCLUSION]);
192     switch(type_inclusion)
193     {
194     case VE_TYPE_INCLUSION::BOITE:
195     {
196     double xmin,ymin,zmin,xmax,ymax,zmax;
197     double longueur_x,longueur_y,longueur_z;
198     longueur_x = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_INCLU_Dim_mu_longueur_x]);
199     longueur_y = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_INCLU_Dim_mu_longueur_y]);
200     longueur_z = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_INCLU_Dim_mu_longueur_z]);
201     xmin = dim_VES[0]/2.-(longueur_x/2.);
202     ymin = dim_VES[1]/2.-(longueur_y/2.);
203     zmin = dim_VES[2]/2.-(longueur_z/2.);
204     xmax = dim_VES[0]/2.+(longueur_x/2.);
205     ymax = dim_VES[1]/2.+(longueur_y/2.);
206     zmax = dim_VES[2]/2.+(longueur_z/2.);
207     inclusion = new MG_PRIMITIVE_BOITE(xmin,ymin,zmin,xmax,ymax,zmax);
208     inclusion_outil = new MG_PRIMITIVE_BOITE(xmin,ymin,zmin,xmax,ymax,zmax);
209     break;
210     }
211     case VE_TYPE_INCLUSION::SPHERE:
212     {
213     double rayon;
214     rayon = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_INCLU_Dim_mu_rayon]);
215     inclusion = new MG_PRIMITIVE_SPHERE(dim_VES[0]/2.,dim_VES[1]/2.,dim_VES[2]/2.,rayon);
216     inclusion_outil = new MG_PRIMITIVE_SPHERE(dim_VES[0]/2.,dim_VES[1]/2.,dim_VES[2]/2.,rayon);
217     break;
218     }
219     case VE_TYPE_INCLUSION::CYLINDRE:
220     {
221     break;
222     }
223     }
224     arbre->ajouter_mg_primitive(inclusion_outil);
225     arbre->ajouter_mg_primitive(inclusion);
226     inclusion_outil->construit();
227     inclusion->construit();
228     MG_OPERATEUR_BOOLEAN_DIFFERENCE *op_difference=new MG_OPERATEUR_BOOLEAN_DIFFERENCE(MAGIC::SEMANTIQUECSG::SOUSTRACTION,cube,inclusion_outil);
229     matrice=op_difference->construit();
230     arbre->ajouter_mg_operateur_boolean(op_difference);
231     arbre->ajouter_mg_primitive(matrice);
232     assemblage->ajouter_mg_primitive(inclusion);
233     assemblage->ajouter_mg_primitive(matrice);
234     return 0;
235     }
236    
237     int VE_VES::generer_VES_aleatoire(void)
238     {
239 couturad 814 affiche((char*)"Generation VES Aleatoire :\n");
240     MG_PRIMITIVE_BOITE *cube = new MG_PRIMITIVE_BOITE(0,0,0,dim_VES[0],dim_VES[1],dim_VES[2]);
241     arbre->ajouter_mg_primitive(cube);
242     cube->construit();
243     TPL_GRILLE<SPHERE_GRILLE*> grille;
244     double cible_frac_volumique = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_INCLU_Dim_pourcentage_volumique]);
245     long nb_max_inclusion = (long)param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_INCLU_Dim_nb_max]);
246     int type_inclusion = (int)param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_INCLU_Type_INCLUSION]);
247     double facteur_distance=1.1;
248     switch(type_inclusion)
249     {
250     case VE_TYPE_INCLUSION::BOITE:
251     {
252     double mu_longueur_x = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_INCLU_Dim_mu_longueur_x]);
253     double sigma_longueur_x = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_INCLU_Dim_sigma_longueur_x]);
254     double mu_longueur_y = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_INCLU_Dim_mu_longueur_y]);
255     double sigma_longueur_y = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_INCLU_Dim_sigma_longueur_y]);
256     double mu_longueur_z = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_INCLU_Dim_mu_longueur_z]);
257     double sigma_longueur_z = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_INCLU_Dim_sigma_longueur_z]);
258     std::random_device seed;
259     std::mt19937_64 generateur(seed());
260     std::normal_distribution<double> distribution_longueur_x(mu_longueur_x,sigma_longueur_x);
261     std::normal_distribution<double> distribution_longueur_y(mu_longueur_y,sigma_longueur_y);
262     std::normal_distribution<double> distribution_longueur_z(mu_longueur_z,sigma_longueur_z);
263    
264     std::uniform_real_distribution<double> distribution_position_x(-mu_longueur_x,dim_VES[0]+mu_longueur_x);
265     std::uniform_real_distribution<double> distribution_position_y(-mu_longueur_y,dim_VES[1]+mu_longueur_y);
266     std::uniform_real_distribution<double> distribution_position_z(-mu_longueur_z,dim_VES[2]+mu_longueur_z);
267     double frac_volumique=0;
268     double volume_inclusion=0;
269     long nb_inclusion=0;
270     long it_max=10000;
271     long it=0;
272     double nb_pas_grille;
273     nb_pas_grille = std::max(dim_VES[0]/(mu_longueur_x-3*sigma_longueur_x),dim_VES[1]/(mu_longueur_y-3*sigma_longueur_y));
274     nb_pas_grille = std::max(nb_pas_grille,dim_VES[2]/(mu_longueur_z-3*sigma_longueur_z));
275     grille.initialiser(0,0,0,dim_VES[0],dim_VES[1],dim_VES[2],nb_pas_grille,nb_pas_grille,nb_pas_grille);
276     long id_sphere_grille=0;
277     std::vector<SPHERE_GRILLE*> vector_sphere_grille;
278     MG_PRIMITIVE* matrice_ves = cube;
279     while(frac_volumique< cible_frac_volumique && nb_inclusion<nb_max_inclusion && it<it_max)
280     {
281     double position[3];
282     position[0] = distribution_position_x(generateur);
283     position[1] = distribution_position_y(generateur);
284     position[2] = distribution_position_z(generateur);
285     double longueur[3];
286     longueur[0] = distribution_longueur_x(generateur);
287     longueur[1] = distribution_longueur_y(generateur);
288     longueur[2] = distribution_longueur_z(generateur);
289     double rayon_sphere_grille = sqrt(0.5*longueur[0]*0.5*longueur[0]+0.5*longueur[1]*0.5*longueur[1]+0.5*longueur[2]*0.5*longueur[2]);
290     bool ok_position=true;
291     // if((position[0])>0.0 && (position[0])<0.2) ok_position=false;
292     // else if((position[1])>0.0 && (position[1])<0.2) ok_position=false;
293     // else if((position[2])>0.0 && (position[2])<0.2) ok_position=false;
294     // else if((position[0]+longueur[0])<dim_VES[0] && (position[0]+longueur[0])>(dim_VES[0]-0.2)) ok_position=false;
295     // else if((position[1]+longueur[1])<dim_VES[1] && (position[1]+longueur[1])>(dim_VES[1]-0.2)) ok_position=false;
296     // else if((position[2]+longueur[2])<dim_VES[2] && (position[2]+longueur[2])>(dim_VES[2]-0.2)) ok_position=false;
297     if((position[0])<0.05) ok_position=false;
298     else if((position[1])<0.05) ok_position=false;
299     else if((position[2])<0.05) ok_position=false;
300     else if((position[0]+longueur[0])>(dim_VES[0]-0.05)) ok_position=false;
301     else if((position[1]+longueur[1])>(dim_VES[1]-0.05)) ok_position=false;
302     else if((position[2]+longueur[2])>(dim_VES[2]-0.05)) ok_position=false;
303     if(ok_position)
304     {
305     double rayon_recherche = rayon_sphere_grille*facteur_distance;
306     TPL_MAP_ENTITE<SPHERE_GRILLE*> map_recherche;
307     grille.rechercher(position[0]+0.5*longueur[0],position[1]+0.5*longueur[1],position[2]+0.5*longueur[2],rayon_recherche,map_recherche);
308     bool ok_insert=true;
309     if(map_recherche.get_nb()>0)
310     {
311     TPL_MAP_ENTITE<SPHERE_GRILLE*>::ITERATEUR it_map_recherche;
312     for(SPHERE_GRILLE* sph = map_recherche.get_premier(it_map_recherche);sph!=NULL;sph=map_recherche.get_suivant(it_map_recherche))
313     {
314     double xyz_sph[3];
315     sph->get_xyz(xyz_sph);
316     if(sqrt((xyz_sph[0]-(position[0]+0.5*longueur[0]))*(xyz_sph[0]-(position[0]+0.5*longueur[0]))+
317     (xyz_sph[1]-(position[1]+0.5*longueur[1]))*(xyz_sph[1]-(position[1]+0.5*longueur[1]))+
318     (xyz_sph[2]-(position[2]+0.5*longueur[2]))*(xyz_sph[2]-(position[2]+0.5*longueur[2])))<(1.1*(sph->get_rayon()+rayon_sphere_grille)))
319     ok_insert=false;
320     }
321     }
322     else if(map_recherche.get_nb()==0) ok_insert=true;
323     if(ok_insert)
324     {
325     MG_PRIMITIVE_BOITE *boite = new MG_PRIMITIVE_BOITE(position[0],position[1],position[2],position[0]+longueur[0],position[1]+longueur[1],position[2]+longueur[2]);
326     boite->construit();
327     MG_PRIMITIVE_BOITE *boite_outil = new MG_PRIMITIVE_BOITE(position[0],position[1],position[2],position[0]+longueur[0],position[1]+longueur[1],position[2]+longueur[2]);
328     boite_outil->construit();
329     if(((position[0] > 0.0) && ((position[0]+longueur[0]) < dim_VES[0])) &&
330     ((position[1] > 0.0) && ((position[1]+longueur[1]) < dim_VES[1])) &&
331     ((position[2] > 0.0) && ((position[2]+longueur[2]) < dim_VES[0])))
332     {
333     char message[1000];
334     sprintf(message,"Inclusion #%li",id_sphere_grille);
335     affiche(message);
336     nb_inclusion++;
337     arbre->ajouter_mg_primitive(boite);
338     arbre->ajouter_mg_primitive(boite_outil);
339     MG_OPERATEUR_BOOLEAN_DIFFERENCE *op_difference=new MG_OPERATEUR_BOOLEAN_DIFFERENCE(MAGIC::SEMANTIQUECSG::SANS,matrice_ves,boite_outil);
340     matrice_ves=op_difference->construit();
341     arbre->ajouter_mg_operateur_boolean(op_difference);
342     assemblage->ajouter_mg_primitive(boite);
343     volume_inclusion+=boite->get_volume();
344     SPHERE_GRILLE *sp = new SPHERE_GRILLE(position[0]+0.5*longueur[0],position[1]+0.5*longueur[1],position[2]+0.5*longueur[2],rayon_sphere_grille,id_sphere_grille);
345     vector_sphere_grille.push_back(sp);
346     grille.inserer(sp);
347     id_sphere_grille++;
348     it=0;
349     }
350     else
351     {
352    
353     MG_PRIMITIVE_BOITE *cube_inter = new MG_PRIMITIVE_BOITE(0,0,0,dim_VES[0],dim_VES[1],dim_VES[2]);
354     cube_inter->construit();
355     MG_OPERATEUR_BOOLEAN_INTERSECTION *op_intersection = new MG_OPERATEUR_BOOLEAN_INTERSECTION(MAGIC::SEMANTIQUECSG::SANS,boite,cube_inter);
356     //MG_OPERATEUR_BOOLEAN_INTERSECTION *op_intersection = new MG_OPERATEUR_BOOLEAN_INTERSECTION(MAGIC::SEMANTIQUECSG::SANS,cube_inter,boite);
357     MG_PRIMITIVE *boite_inter = op_intersection->construit();
358     if(boite_inter->get_volume()> 0.25*(boite->get_volume()))
359     {
360     volume_inclusion+=boite_inter->get_volume();
361     char message[1000];
362     sprintf(message,"Inclusion #%li : COUPÉE",id_sphere_grille);
363     affiche(message);
364     nb_inclusion++;
365     MG_OPERATEUR_BOOLEAN_DIFFERENCE *op_difference=new MG_OPERATEUR_BOOLEAN_DIFFERENCE(MAGIC::SEMANTIQUECSG::SANS,matrice_ves,boite_outil);
366     matrice_ves=op_difference->construit();
367     arbre->ajouter_mg_primitive(boite);
368     arbre->ajouter_mg_primitive(boite_outil);
369     arbre->ajouter_mg_operateur_boolean(op_difference);
370     arbre->ajouter_mg_primitive(boite_inter);
371     arbre->ajouter_mg_primitive(cube_inter);
372     assemblage->ajouter_mg_primitive(boite_inter);
373     SPHERE_GRILLE *sp = new SPHERE_GRILLE(position[0]+0.5*longueur[0],position[1]+0.5*longueur[1],position[2]+0.5*longueur[2],rayon_sphere_grille,id_sphere_grille);
374     vector_sphere_grille.push_back(sp);
375     grille.inserer(sp);
376     id_sphere_grille++;
377     it=0;
378     }
379     else
380     {
381     delete boite;
382     delete boite_outil;
383     delete boite_inter;
384     delete cube_inter;
385     delete op_intersection;
386     }
387     }
388     frac_volumique = volume_inclusion/(cube->get_volume());
389     }
390     it++;
391     }
392     }
393     arbre->ajouter_mg_primitive(matrice_ves);
394     assemblage->ajouter_mg_primitive(matrice_ves);
395     char message[1000];
396     sprintf(message,"echo -n %li:%lf: >> resultats.txt",nb_inclusion,frac_volumique);
397     system(message);
398     break;
399     }
400     }
401     return 0;
402 couturad 801 }
403    
404     int VE_VES::evaluer_geometrie(void)
405     {
406     char nom_fichier[1000];
407     sprintf(nom_fichier,"%s",param_ver->get_nom(VE_PARAM_CHAR[VE_PARAM::VE_IO_nom_fichier_STEP_BREP]).c_str());
408     bool export_STEP, avec_fusion,import_STL;
409     export_STEP = (int)param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_IO_Bool_export_STEP]);
410     avec_fusion = (int)param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_IO_Bool_fusion]);
411     import_STL = (int)param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_IO_Bool_import_STL]);
412     double precision,eps;
413     precision = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_Dim_precision]);
414     eps = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_IO_STL_eps]);
415     arbre->evaluer_geometrie(assemblage,nom_fichier,avec_fusion,export_STEP,import_STL,precision,eps);
416     return 0;
417     }
418    
419     int VE_VES::creer_mg_maillage(void)
420     {
421     affiche((char*)"Maillage du VER :");
422     MG_GEOMETRIE *geo = gest->get_mg_geometrie(0);
423     MG_MAILLAGE *maillage = new MG_MAILLAGE(geo);
424     gest->ajouter_mg_maillage(maillage);
425     MAILLEUR0D m0d(maillage,geo);
426     affiche((char*)" Maillage 0D :");
427     m0d.maille();
428     MAILLEUR1D m1d(maillage,geo,metrique);
429     affiche((char*)" Maillage 1D :");
430     //m1d.active_affichage(fonc_affiche);
431     m1d.maille();
432     MAILLEUR2D m2d(maillage,geo,metrique);
433     affiche((char*)" Maillage 2D :");
434     //m2d.active_affichage(fonc_affiche);
435     m2d.maille();
436     MAILLEUR3D m3d(maillage,geo,metrique,false);
437     affiche((char*)" Maillage 3D :");
438     m3d.active_affichage(fonc_affiche);
439     m3d.maille();
440     return 0;
441     }
442    
443     int VE_VES::creer_carte(void)
444     {
445     int type_carte = (int)param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_CARTE_Type]);
446     switch(type_carte)
447     {
448     case VE_TYPE_CARTE::CONSTANTE:
449     {
450     affiche((char*)"Generation de carte de taille constante V1");
451     double ecart_nodal = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_CARTE_Dim_ecart_nodal]);
452     double fechantillonnage = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_CARTE_Dim_fechantillonnage]);
453     double festimation = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_CARTE_Dim_festimation]);
454     FCT_GENERATEUR_CONSTANTE *carte = new FCT_GENERATEUR_CONSTANTE(*gest,ecart_nodal);
455     affiche((char*)"Construction de la carte");
456     carte->construit(fechantillonnage,festimation);
457     affiche((char*)"Enregistrement de la carte");
458     char nom_fichier[1000];
459     sprintf(nom_fichier,"%s",param_ver->get_nom(VE_PARAM_CHAR[VE_PARAM::VE_IO_nom_fichier_CARTE]).c_str());
460     carte->enregistrer(nom_fichier);
461     metrique = carte;
462     affiche((char*)"Fin");
463     break;
464     }
465     }
466     return 0;
467     }
468    
469     int VE_VES::appliquer_conditions_limites_plan(TPL_MAP_ENTITE< MG_FACE* >* plan, char* condition, double valeur,bool topo_sous_jacente)
470     {
471     TPL_MAP_ENTITE<MG_FACE*>::ITERATEUR it_face_plan;
472     for(MG_FACE *face = plan->get_premier(it_face_plan);face!=NULL;face=plan->get_suivant(it_face_plan))
473     {
474     face->ajouter_ccf(condition,valeur);
475     if(topo_sous_jacente==true)
476     {
477     TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> map_topo;
478     face->get_topologie_sousjacente(&map_topo);
479     TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*>::ITERATEUR it_topo;
480     for(MG_ELEMENT_TOPOLOGIQUE *topo=map_topo.get_premier(it_topo);topo!=NULL;topo=map_topo.get_suivant(it_topo))
481     {
482     topo->ajouter_ccf(condition,valeur);
483     }
484     }
485     }
486     }
487    
488     int VE_VES::appliquer_conditions_limites_plan(TPL_MAP_ENTITE< MG_FACE* >* plan, char* condition, std::string formule,std::vector<std::string> &listvariable,bool topo_sous_jacente)
489     {
490     TPL_MAP_ENTITE<MG_FACE*>::ITERATEUR it_face_plan;
491     for(MG_FACE *face = plan->get_premier(it_face_plan);face!=NULL;face=plan->get_suivant(it_face_plan))
492     {
493     face->ajouter_ccf(condition,formule,listvariable);
494     if(topo_sous_jacente==true)
495     {
496     TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> map_topo;
497     face->get_topologie_sousjacente(&map_topo);
498     TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*>::ITERATEUR it_topo;
499     for(MG_ELEMENT_TOPOLOGIQUE *topo=map_topo.get_premier(it_topo);topo!=NULL;topo=map_topo.get_suivant(it_topo))
500     {
501     topo->ajouter_ccf(condition,formule,listvariable);
502     }
503     }
504     }
505     }
506    
507     int VE_VES::appliquer_conditions_limites(int type_chargement)
508     {
509     TPL_MAP_ENTITE<MG_FACE*> plan_xy_z0;
510     TPL_MAP_ENTITE<MG_FACE*> plan_xy_z1;
511     TPL_MAP_ENTITE<MG_FACE*> plan_yz_x0;
512     TPL_MAP_ENTITE<MG_FACE*> plan_yz_x1;
513     TPL_MAP_ENTITE<MG_FACE*> plan_xz_y0;
514     TPL_MAP_ENTITE<MG_FACE*> plan_xz_y1;
515     MG_SOMMET* sommet_origine;
516     MG_SOMMET* sommet_x1y0z0;
517     MG_SOMMET* sommet_x0y1z0;
518     MG_SOMMET* sommet_x0y0z1;
519     double eps = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_Dim_precision]);
520    
521     MG_GEOMETRIE *geo = gest->get_mg_geometrie(0);
522     std::map<unsigned long,MG_FACE*,std::less<unsigned long>>::iterator it_face;
523     for(MG_FACE *face = geo->get_premier_face(it_face);face!=NULL;face=geo->get_suivant_face(it_face))
524     {
525     TPL_SET<MG_ELEMENT_MAILLAGE*> *liste_element_maillage = face->get_lien_maillage();
526     MG_TRIANGLE *tri = (MG_TRIANGLE*)liste_element_maillage->get(0);
527     MG_NOEUD *nd[3];
528     nd[0] = tri->get_noeud1();
529     nd[1] = tri->get_noeud2();
530     nd[2] = tri->get_noeud3();
531    
532     if(OPERATEUR::egal(nd[0]->get_z(),0,eps) && OPERATEUR::egal(nd[1]->get_z(),0,eps) && OPERATEUR::egal(nd[2]->get_z(),0,eps))
533     {
534     plan_xy_z0.ajouter(face);
535     }
536     else if(OPERATEUR::egal(nd[0]->get_z(),dim_VES[2],eps) && OPERATEUR::egal(nd[1]->get_z(),dim_VES[2],eps) && OPERATEUR::egal(nd[2]->get_z(),dim_VES[2],eps))
537     {
538     plan_xy_z1.ajouter(face);
539     }
540     else if(OPERATEUR::egal(nd[0]->get_x(),0,eps) && OPERATEUR::egal(nd[1]->get_x(),0,eps) && OPERATEUR::egal(nd[2]->get_x(),0,eps))
541     {
542     plan_yz_x0.ajouter(face);
543     }
544     else if(OPERATEUR::egal(nd[0]->get_x(),dim_VES[0],eps) && OPERATEUR::egal(nd[1]->get_x(),dim_VES[0],eps) && OPERATEUR::egal(nd[2]->get_x(),dim_VES[0],eps))
545     {
546     plan_yz_x1.ajouter(face);
547     }
548     else if(OPERATEUR::egal(nd[0]->get_y(),0,eps) && OPERATEUR::egal(nd[1]->get_y(),0,eps) && OPERATEUR::egal(nd[2]->get_y(),0,eps))
549     {
550     plan_xz_y0.ajouter(face);
551     }
552     else if(OPERATEUR::egal(nd[0]->get_y(),dim_VES[1],eps) && OPERATEUR::egal(nd[1]->get_y(),dim_VES[1],eps) && OPERATEUR::egal(nd[2]->get_y(),dim_VES[1],eps))
553     {
554     plan_xz_y1.ajouter(face);
555     }
556     }
557    
558     std::map<unsigned long,MG_SOMMET*,std::less<unsigned long>>::iterator it_sommet;
559     for(MG_SOMMET* som=geo->get_premier_sommet(it_sommet);som!=NULL;som=geo->get_suivant_sommet(it_sommet))
560     {
561     TPL_SET<MG_ELEMENT_MAILLAGE*> *liste_element_maillage = som->get_lien_maillage();
562     MG_NOEUD* nd = (MG_NOEUD*)liste_element_maillage->get(0);
563     if(OPERATEUR::egal(nd->get_x(),0,eps) && OPERATEUR::egal(nd->get_y(),0,eps) && OPERATEUR::egal(nd->get_z(),0,eps))
564     {
565     sommet_origine=som;
566     }
567     else if(OPERATEUR::egal(nd->get_x(),dim_VES[0],eps) && OPERATEUR::egal(nd->get_y(),0,eps) && OPERATEUR::egal(nd->get_z(),0,eps))
568     {
569     sommet_x1y0z0=som;
570     }
571     else if(OPERATEUR::egal(nd->get_x(),0,eps) && OPERATEUR::egal(nd->get_y(),dim_VES[1],eps) && OPERATEUR::egal(nd->get_z(),0,eps))
572     {
573     sommet_x0y1z0=som;
574     }
575     else if(OPERATEUR::egal(nd->get_x(),0,eps) && OPERATEUR::egal(nd->get_y(),0,eps) && OPERATEUR::egal(nd->get_z(),dim_VES[2],eps))
576     {
577     sommet_x0y0z1=som;
578     }
579     }
580    
581     int type_cl = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_VER_CL_Type]);
582     double valeur_cl = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_VER_CL_Dim_valeur]);
583    
584     switch(type_chargement)
585     {
586     case VE_TYPE_CHARGEMENT::SPHERIQUE:
587     {
588     switch(type_cl)
589     {
590     case VE_TYPE_CL::DEFORMATION_HOMOGENE:
591     {
592     appliquer_conditions_limites_plan(&plan_xy_z0,(char*)"Dz",0.0,true);
593     appliquer_conditions_limites_plan(&plan_xy_z1,(char*)"Dz",valeur_cl,true);
594     appliquer_conditions_limites_plan(&plan_yz_x0,(char*)"Dx",0.0,true);
595     appliquer_conditions_limites_plan(&plan_yz_x1,(char*)"Dx",valeur_cl,true);
596     appliquer_conditions_limites_plan(&plan_xz_y0,(char*)"Dy",0.0,true);
597     appliquer_conditions_limites_plan(&plan_xz_y1,(char*)"Dy",valeur_cl,true);
598     break;
599     }
600     case VE_TYPE_CL::CONTRAINTE_HOMOGENE:
601     {
602     appliquer_conditions_limites_plan(&plan_xy_z0,(char*)"Dz",0.0,true);
603     appliquer_conditions_limites_plan(&plan_xy_z1,(char*)"Pz",valeur_cl,false);
604     appliquer_conditions_limites_plan(&plan_yz_x0,(char*)"Dx",0.0,true);
605     appliquer_conditions_limites_plan(&plan_yz_x1,(char*)"Px",valeur_cl,false);
606     appliquer_conditions_limites_plan(&plan_xz_y0,(char*)"Dy",0.0,true);
607     appliquer_conditions_limites_plan(&plan_xz_y1,(char*)"Py",valeur_cl,false);
608     break;
609     }
610     case VE_TYPE_CL::FLUX_CHALEUR_HOMOGENE:
611     {
612     break;
613     }
614     case VE_TYPE_CL::GRADIENT_TEMPERATURE_HOMOGENE:
615     {
616     break;
617     }
618     }
619     break;
620     }
621    
622     case VE_TYPE_CHARGEMENT::DEVIATORIQUE:
623     {
624     switch(type_cl)
625     {
626     case VE_TYPE_CL::DEFORMATION_HOMOGENE:
627     {
628     char chr_valeur_cl[1000];
629     sprintf(chr_valeur_cl,"%lf",valeur_cl);
630     std::string str_valeur_cl = chr_valeur_cl;
631     std::string formule_Dx = str_valeur_cl+ "*Y";
632     std::string formule_Dy = str_valeur_cl+ "*Z";
633     std::string formule_Dz = str_valeur_cl+ "*X";
634     std::vector<std::string> liste_variable_formule_Dx;
635     std::vector<std::string> liste_variable_formule_Dy;
636     std::vector<std::string> liste_variable_formule_Dz;
637     liste_variable_formule_Dx.push_back((std::string)"Y");
638     liste_variable_formule_Dy.push_back((std::string)"Z");
639     liste_variable_formule_Dz.push_back((std::string)"X");
640     appliquer_conditions_limites_plan(&plan_xy_z0,(char*)"Dx",formule_Dx,liste_variable_formule_Dx,true);
641     appliquer_conditions_limites_plan(&plan_xy_z0,(char*)"Dy",formule_Dy,liste_variable_formule_Dy,true);
642     appliquer_conditions_limites_plan(&plan_xy_z0,(char*)"Dz",formule_Dz,liste_variable_formule_Dz,true);
643     appliquer_conditions_limites_plan(&plan_xy_z1,(char*)"Dx",formule_Dx,liste_variable_formule_Dx,true);
644     appliquer_conditions_limites_plan(&plan_xy_z1,(char*)"Dy",formule_Dy,liste_variable_formule_Dy,true);
645     appliquer_conditions_limites_plan(&plan_xy_z1,(char*)"Dz",formule_Dz,liste_variable_formule_Dz,true);
646     appliquer_conditions_limites_plan(&plan_yz_x0,(char*)"Dx",formule_Dx,liste_variable_formule_Dx,true);
647     appliquer_conditions_limites_plan(&plan_yz_x0,(char*)"Dy",formule_Dy,liste_variable_formule_Dy,true);
648     appliquer_conditions_limites_plan(&plan_yz_x0,(char*)"Dz",formule_Dz,liste_variable_formule_Dz,true);
649     appliquer_conditions_limites_plan(&plan_yz_x1,(char*)"Dx",formule_Dx,liste_variable_formule_Dx,true);
650     appliquer_conditions_limites_plan(&plan_yz_x1,(char*)"Dy",formule_Dy,liste_variable_formule_Dy,true);
651     appliquer_conditions_limites_plan(&plan_yz_x1,(char*)"Dz",formule_Dz,liste_variable_formule_Dz,true);
652     appliquer_conditions_limites_plan(&plan_xz_y0,(char*)"Dx",formule_Dx,liste_variable_formule_Dx,true);
653     appliquer_conditions_limites_plan(&plan_xz_y0,(char*)"Dy",formule_Dy,liste_variable_formule_Dy,true);
654     appliquer_conditions_limites_plan(&plan_xz_y0,(char*)"Dz",formule_Dz,liste_variable_formule_Dz,true);
655     appliquer_conditions_limites_plan(&plan_xz_y1,(char*)"Dx",formule_Dx,liste_variable_formule_Dx,true);
656     appliquer_conditions_limites_plan(&plan_xz_y1,(char*)"Dy",formule_Dy,liste_variable_formule_Dy,true);
657     appliquer_conditions_limites_plan(&plan_xz_y1,(char*)"Dz",formule_Dz,liste_variable_formule_Dz,true);
658     break;
659     }
660    
661     case VE_TYPE_CL::CONTRAINTE_HOMOGENE:
662     {
663     sommet_origine->ajouter_ccf((char*)"Dt",0.0);
664     sommet_x1y0z0->ajouter_ccf((char*)"Dy",0.0);
665     sommet_x1y0z0->ajouter_ccf((char*)"Dz",0.0);
666     sommet_x0y1z0->ajouter_ccf((char*)"Dz",0.0);
667     appliquer_conditions_limites_plan(&plan_xy_z0,(char*)"Px",-valeur_cl,true);
668     appliquer_conditions_limites_plan(&plan_xy_z0,(char*)"Py",-valeur_cl,true);
669     appliquer_conditions_limites_plan(&plan_xy_z1,(char*)"Px",valeur_cl,true);
670     appliquer_conditions_limites_plan(&plan_xy_z1,(char*)"Py",valeur_cl,true);
671     appliquer_conditions_limites_plan(&plan_xz_y0,(char*)"Px",-valeur_cl,true);
672     appliquer_conditions_limites_plan(&plan_xz_y0,(char*)"Pz",-valeur_cl,true);
673     appliquer_conditions_limites_plan(&plan_xz_y1,(char*)"Px",valeur_cl,true);
674     appliquer_conditions_limites_plan(&plan_xz_y1,(char*)"Pz",valeur_cl,true);
675     appliquer_conditions_limites_plan(&plan_yz_x0,(char*)"Py",-valeur_cl,true);
676     appliquer_conditions_limites_plan(&plan_yz_x0,(char*)"Pz",-valeur_cl,true);
677     appliquer_conditions_limites_plan(&plan_yz_x1,(char*)"Py",valeur_cl,true);
678     appliquer_conditions_limites_plan(&plan_yz_x1,(char*)"Pz",valeur_cl,true);
679     break;
680     }
681     case VE_TYPE_CL::FLUX_CHALEUR_HOMOGENE:
682     {
683     break;
684     }
685     case VE_TYPE_CL::GRADIENT_TEMPERATURE_HOMOGENE:
686     {
687     break;
688     }
689     }
690     break;
691     }
692     }
693     return 0;
694     }
695    
696     int VE_VES::appliquer_materiaux(void)
697     {
698     double nu_inclu = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_INCLU_Dim_mu_nu]);
699     double nu_matrice = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_MATRICE_Dim_mu_nu]);
700     double Em_inclu = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_INCLU_Dim_mu_Em]);
701     double Em_matrice = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_MATRICE_Dim_mu_Em]);
702     MG_GEOMETRIE *geo = gest->get_mg_geometrie(0);
703     MG_VOLUME *matrice = geo->get_mg_volume(geo->get_nb_mg_volume()-1);
704     std::map<unsigned long,MG_VOLUME*,std::less<unsigned long>>::iterator it_volume;
705     for(MG_VOLUME *volume = geo->get_premier_volume(it_volume);volume!=matrice;volume=geo->get_suivant_volume(it_volume))
706     {
707     volume->ajouter_ccf((char*)"nu",nu_inclu);
708     volume->ajouter_ccf((char*)"Em",Em_inclu);
709     }
710     matrice->ajouter_ccf((char*)"nu",nu_matrice);
711     matrice->ajouter_ccf((char*)"Em",Em_matrice);
712     return 0;
713     }
714    
715     int VE_VES::creer_fem_maillage(void)
716     {
717     MG_GEOMETRIE *geo = gest->get_mg_geometrie(0);
718     MG_MAILLAGE *mai = gest->get_mg_maillage(gest->get_nb_mg_maillage()-1);
719     int degre_mail = (int)param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_MAIL_Dim_degre_maillage]);
720     FEM_MAILLAGE *fem = new FEM_MAILLAGE(geo,mai,degre_mail);
721     gest->ajouter_fem_maillage(fem);
722 francois 883 MAILLEUR_FEM m;
723     m.maille(fem);
724 couturad 801 return 0;
725     }
726    
727     int VE_VES::calcul_aster(void)
728     {
729     char coderesu[500];
730     char nomparam[500];
731     char nometude[500];
732     sprintf(nometude,"%s",param_ver->get_nom(VE_PARAM_CHAR[VE_PARAM::VE_IO_nom_fichier_ETUDE]).c_str());
733     sprintf(nomparam,"%s",param_ver->get_nom(VE_PARAM_CHAR[VE_PARAM::VE_IO_nom_fichier_param_ASTER]).c_str());
734     sprintf(coderesu,"%s",param_ver->get_nom(VE_PARAM_CHAR[VE_PARAM::VE_IO_code_resu]).c_str());
735     MGASTER mgaster;
736     mgaster.active_affichage(fonc_affiche);
737     mgaster.calcule(nomparam,gest->get_fem_maillage(0),nometude,MAGIC::CALCUL_ASTER::ELASTIQUE,coderesu);
738     }
739    
740     int VE_VES::homogeneisation(void)
741     {
742     FEM_SOLUTION* sol=gest->get_fem_solution(1);
743     FEM_SOLUTION* sol2=gest->get_fem_solution(2);
744     double exx=sol->get_moyenne_volumique_champs(NULL,0);
745     double eyy=sol->get_moyenne_volumique_champs(NULL,1);
746     double ezz=sol->get_moyenne_volumique_champs(NULL,2);
747     double exy=sol->get_moyenne_volumique_champs(NULL,3);
748     double exz=sol->get_moyenne_volumique_champs(NULL,4);
749     double eyz=sol->get_moyenne_volumique_champs(NULL,5);
750     double sxx=sol2->get_moyenne_volumique_champs(NULL,0);
751     double syy=sol2->get_moyenne_volumique_champs(NULL,1);
752     double szz=sol2->get_moyenne_volumique_champs(NULL,2);
753     double sxy=sol2->get_moyenne_volumique_champs(NULL,3);
754     double sxz=sol2->get_moyenne_volumique_champs(NULL,4);
755     double syz=sol2->get_moyenne_volumique_champs(NULL,5);
756     double traces=sxx+syy+szz;
757     double tracee=exx+eyy+ezz;
758     char message[1000];
759     affiche((char*)"");
760     sprintf(message, "Tenseur des deformations :"); affiche(message);
761     sprintf(message, " | exx=% 3.2le exy=% 3.2le exz=% 3.2le |",exx,exy,exz); affiche(message);
762     sprintf(message, " | exy=% 3.2le eyy=% 3.2le eyz=% 3.2le |",exy,eyy,eyz); affiche(message);
763     sprintf(message, " | exz=% 3.2le eyz=% 3.2le ezz=% 3.2le |",exz,eyz,ezz); affiche(message);
764     affiche((char*)"");
765     sprintf(message, " Trace e = % 3.2le",tracee); affiche(message);
766     affiche((char*)"");
767     sprintf(message, "Tenseur des contraintes :"); affiche(message);
768     sprintf(message, " | sxx=% 3.2le sxy=% 3.2le sxz=% 3.2le |",sxx,sxy,sxz); affiche(message);
769     sprintf(message, " | sxy=% 3.2le syy=% 3.2le syz=% 3.2le |",sxy,syy,syz); affiche(message);
770     sprintf(message, " | sxz=% 3.2le syz=% 3.2le szz=% 3.2le |",sxz,syz,szz); affiche(message);
771     affiche((char*)"");
772     sprintf(message, " Trace s = % 3.2le",traces); affiche(message);
773     affiche((char*)"");
774     int type_chargement = param_ver->get_valeur(VE_PARAM_CHAR[VE_PARAM::VE_VER_CL_CHARGEMENT_type]);
775     if(type_chargement==VE_TYPE_CHARGEMENT::SPHERIQUE)
776     {
777     double K = traces/3./tracee;
778     sprintf(message, "Module de compressibilite apparent :"); affiche(message);
779     sprintf(message, " Kapp = % 3.5le",K); affiche(message);
780 couturad 814 char messagesortie[1000];
781     sprintf(messagesortie,"echo %lf >> resultats.txt",K);
782     system(messagesortie);
783 couturad 801 }
784     else if(type_chargement==VE_TYPE_CHARGEMENT::DEVIATORIQUE)
785     {
786     double G = (1./3.)*((syz/(2.*eyz))+(sxz/(2.*exz))+(sxy/(2.*exy)));
787     sprintf(message, "Module de cisaillement apparent :"); affiche(message);
788     sprintf(message, " Gapp = % 3.5le",G); affiche(message);
789 couturad 814 char messagesortie[1000];
790     sprintf(messagesortie,"echo %lf >> resultats.txt",G);
791     system(messagesortie);
792 couturad 801 }
793 couturad 814
794 couturad 801 }