ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/aster/src/mgopt_simp.cpp
Revision: 748
Committed: Thu Oct 8 19:14:25 2015 UTC (9 years, 7 months ago) by francois
File size: 26918 byte(s)
Log Message:
Calcul multi volume plus choix du solveur de code aster via le fichier de parametres

File Contents

# User Rev Content
1 francois 239 #include "gestionversion.h"
2     #include "mgopt_simp.h"
3     #include "mg_file.h"
4     #include "mg_export.h"
5     #include <string.h>
6 francois 258 #include <math.h>
7 francois 468 #include "fct_taille_fem_solution.h"
8 francois 731 #include "fct_generateur_3d.h"
9 nana 565 #include "ot_cpu.h"
10 francois 239
11    
12     class SIMP_TETRA
13     {
14     public:
15 francois 309 FEM_ELEMENT3* tet;
16 francois 239 double densite;
17 francois 242 double new_densite;
18 francois 239 int design;
19 francois 242 double new_denergie;
20 francois 239 double energie;
21     double denergie;
22     double volume;
23     int maille_niveau;
24     int indice;
25     double distance_ref;
26 francois 247 double distance_ref2;
27 francois 239 double centre[3];
28     BOITE_3D get_boite_3D(void)
29     {
30     return tet->get_boite_3D();
31     };
32 francois 468 std::vector<SIMP_TETRA*> voisin;
33     std::vector<SIMP_TETRA*> voisin2;
34 francois 239 unsigned long get_id(void) {return tet->get_id();};
35     double distance(SIMP_TETRA* tet2)
36     {
37     double dx=centre[0]-tet2->centre[0];
38     double dy=centre[1]-tet2->centre[1];
39     double dz=centre[2]-tet2->centre[2];
40     return sqrt(dx*dx+dy*dy+dz*dz);
41     };
42     };
43    
44    
45    
46    
47    
48    
49    
50    
51    
52    
53 francois 494 MGOPT_SIMP::MGOPT_SIMP(bool save):MGOPT(save)
54 francois 239 {
55 francois 396 params.ajouter("f",0.3,OT_PARAMETRES::DOUBLE,"Fraction volumique de la methode SIMP");
56     params.ajouter("rminc",0.,OT_PARAMETRES::DOUBLE,"0. Valeur de la carte de taille (si pas de maillage le parametre fichiercarte permet d'utiliser une carte de taille) sinon valeur de rmin pour le lissage de la compliance");
57     params.ajouter("rmind",0.,OT_PARAMETRES::DOUBLE,"0. Valeur de la carte de taille (si pas de maillage le parametre fichiercarte permet d'utiliser une carte de taille) sinon valeur de rmin pour le lissage de la densite");
58     params.ajouter("coefvoisinc",1.0,OT_PARAMETRES::DOUBLE,"Coefficient de multiplication de rminc");
59     params.ajouter("coefvoisind",1.0,OT_PARAMETRES::DOUBLE,"Coefficient de multiplication de rmind");
60     params.ajouter("nbrniveau",150.,OT_PARAMETRES::DOUBLE,"Nombre de matériaux utilisés");
61     params.ajouter("p",3.,OT_PARAMETRES::DOUBLE,"Coefficient de penalite du module d'Young");
62     params.ajouter("ro_void",0.001,OT_PARAMETRES::DOUBLE,"Densite minimale consideree comme nulle");
63     params.ajouter("m",0.2,OT_PARAMETRES::DOUBLE,"Limite d'evolution de densité entre deux itérations");
64     params.ajouter("eta",0.5,OT_PARAMETRES::DOUBLE,"Coefficent d'affectation du beta");
65     params.ajouter("type_lissage",1.,OT_PARAMETRES::DOUBLE,"0. Pas de lissage 1. Lissage de la compliance 2. Lissage de la densité 3. Lissage de la compliance et de la densite");
66     params.ajouter("lissage_densite",1.,OT_PARAMETRES::DOUBLE,"0. Complet 1. Derniere iteration");
67     params.ajouter("type_lissage_densite",1.,OT_PARAMETRES::DOUBLE,"0. Distance 1. Distance pondéré 2. Gaussien 3. Gaussien pondéré");
68     params.ajouter("type_lissage_compliance",1.,OT_PARAMETRES::DOUBLE,"0. Sigmund 1997 1. Sigmund 2007");
69     params.ajouter("kc",1.0,OT_PARAMETRES::DOUBLE,"Ponderation de la distance dans le lissage de la compliance");
70     params.ajouter("kd",1.0,OT_PARAMETRES::DOUBLE,"Ponderation de la distance dans le lissage de la densité");
71     params.ajouter("convergence_lagrange",0.1,OT_PARAMETRES::DOUBLE,"Critere de convergence de la recherche du multiplicateur de Lagrange en \%");
72     params.ajouter("iter_max",50.,OT_PARAMETRES::DOUBLE,"Nombre d'iteration maximale dans le processus SIMP");
73     params.ajouter("critere_convergence",0.5,OT_PARAMETRES::DOUBLE,"Critere de convergence de la méthode SIMP en \%");
74     params.ajouter("ro_min",0.8,OT_PARAMETRES::DOUBLE,"Seuil de conservation des éléments");
75 francois 433 params.ajouter("iter_vue",0.,OT_PARAMETRES::DOUBLE,"generation d'un fichier gmsh tous les iter_vue. 0 pas de generation");
76 nana 565 params.ajouter("type_filtrage",0.,OT_PARAMETRES::DOUBLE,"0. filtrage par distance 1. filtrage par couches");
77 francois 578 params.ajouter("nb_couchesc",1.,OT_PARAMETRES::DOUBLE,"nombre de couches filtrage de la compliance");
78     params.ajouter("nb_couchesd",1.,OT_PARAMETRES::DOUBLE,"nombre de couches filtrage de la densite");
79 francois 239 }
80    
81     MGOPT_SIMP::MGOPT_SIMP(MGOPT_SIMP &mdd):MGOPT(mdd)
82     {
83     }
84    
85    
86     MGOPT_SIMP::~MGOPT_SIMP()
87     {
88     }
89    
90    
91 francois 706 void MGOPT_SIMP::optimisation(int num_adapt,FEM_MAILLAGE* fem,char *nomparamaster,int iter)
92 francois 239 {
93 francois 748 MG_GEOMETRIE* geo=fem->get_mg_geometrie();
94 francois 272 affiche((char*)" Initialisation");
95 francois 468 std::vector<SIMP_TETRA*> lsttet;
96 francois 247 double f=params.get_valeur("f");
97     double rminc=params.get_valeur("rminc");
98     double rmind=params.get_valeur("rmind");
99     if (((rminc<1e-16)||(rmind<1e-16)) && (carte==NULL))
100 francois 246 {
101 francois 272 affiche((char*)" Lecture de la carte de taille");
102 francois 468 int typecarte=(int)params.get_valeur("typecarte");
103     std::string fichiercarte=params.get_nom("fichiercarte");
104     if ((typecarte==0) &&(typecarte==1))
105     {
106     carte=new FCT_GENERATEUR_3D<4>;
107     carte->lire((char *)fichiercarte.c_str());
108     }
109     if ((typecarte==2) &&(typecarte==3)) carte=new FCT_TAILLE_FEM_SOLUTION((char *)fichiercarte.c_str());
110     }
111 francois 239 double volume_tot;
112     double volume_design=0;
113     double volume_non_design=0;
114     double unite=fem->get_mg_maillage()->get_mg_geometrie()->get_valeur_unite();
115 francois 309 LISTE_FEM_ELEMENT3::iterator it;
116 francois 239 int ntet=0;
117 francois 309 for (FEM_ELEMENT3* tet=fem->get_premier_element3(it);tet!=NULL;tet=fem->get_suivant_element3(it))
118 francois 239 {
119     SIMP_TETRA* tet2=new SIMP_TETRA;
120     tet2->tet=tet;
121     tet2->tet->change_numero(ntet);
122     ntet++;
123     tet2->centre[0]=0.;
124     tet2->centre[1]=0.;
125     tet2->centre[2]=0.;
126     int nbnoeud=tet->get_nb_fem_noeud();
127     for (int i=0;i<nbnoeud;i++)
128     {
129 francois 242 tet2->centre[0]=tet2->centre[0]+tet->get_fem_noeud(i)->get_x()*unite;
130     tet2->centre[1]=tet2->centre[1]+tet->get_fem_noeud(i)->get_y()*unite;
131     tet2->centre[2]=tet2->centre[2]+tet->get_fem_noeud(i)->get_z()*unite;
132 francois 239 }
133     tet2->centre[0]=tet2->centre[0]/nbnoeud;
134     tet2->centre[1]=tet2->centre[1]/nbnoeud;
135     tet2->centre[2]=tet2->centre[2]/nbnoeud;
136     int lig,col;
137     double jac[9],uv[3]={1./4.,1./4.,1./4.};
138 francois 242 tet2->volume=1./6.*tet->get_jacobien(jac,uv,lig,col,unite);
139 francois 247 tet2->distance_ref=rminc;
140     tet2->distance_ref2=rmind;
141 cuillier 395 double xyz[3]={tet2->centre[0]/unite,tet2->centre[1]/unite,tet2->centre[2]/unite};
142 francois 246 double val[9];
143     if (tet2->distance_ref<1e-16)
144 francois 242 {
145 francois 246 carte->evaluer(xyz,val);
146     tet2->distance_ref=1./sqrt(val[0]);
147 francois 242 }
148 francois 247 if (tet2->distance_ref2<1e-16)
149     {
150     carte->evaluer(xyz,val);
151     tet2->distance_ref2=1./sqrt(val[0]);
152     }
153     tet2->distance_ref=params.get_valeur("coefvoisinc")*tet2->distance_ref*unite;
154     tet2->distance_ref2=params.get_valeur("coefvoisind")*tet2->distance_ref2*unite;
155 francois 239 if (((MG_TETRA*)tet->get_mg_element_maillage())->get_origine()!=IMPOSE)
156     {
157     tet2->design=1;
158 francois 247 tet2->densite=f;
159 francois 239 volume_design=volume_design+tet2->volume;
160     }
161     else
162     {
163     tet2->design=0;
164     tet2->densite=1.;
165     volume_non_design=volume_non_design+tet2->volume;
166     }
167     lsttet.insert(lsttet.end(),tet2);
168     }
169 francois 268 volume_tot=volume_design+volume_non_design;
170 francois 309 double volumemoyen=volume_tot/fem->get_nb_fem_element3();
171 francois 272 affiche((char*)" Recherche de voisins");
172 francois 239 int nbsimptet=lsttet.size();
173 nana 565 if (params.get_valeur("type_filtrage")==0.)
174 francois 239 {
175 nana 565 for (int i=0;i<nbsimptet;i++)
176     {
177     SIMP_TETRA* tet=lsttet[i];
178     if (tet->design==0) continue;
179     ajouter_voisin_distance(i,tet,lsttet);
180     }
181 francois 239 }
182 nana 565 if (params.get_valeur("type_filtrage")==1.)
183     {
184     for (int i=0;i<nbsimptet;i++)
185     {
186     SIMP_TETRA* tet=lsttet[i];
187     if (tet->design==0) continue;
188 francois 578 ajouter_voisin_couche(i,tet,lsttet,params.get_valeur("nb_couchesc"),params.get_valeur("nb_couchesd"));
189 nana 565 }
190     }
191    
192 francois 239 int nbiteration=1;
193     int ok=0;
194     int niveaumax=(int)params.get_valeur("nbrniveau");
195 francois 339 int itervue=(int)params.get_valeur("iter_vue");
196 francois 247 double p=params.get_valeur("p");
197 francois 468 std::vector<double> Ctotiter;
198 francois 239 while (ok==0)
199     {
200 francois 247 if (nbiteration>params.get_valeur("iter_max")) break;
201 francois 239 char message[255];
202 francois 691 if (num_adapt==0) sprintf(message," Iteration \033[1;32m%d\033[1;33m",nbiteration);
203     else sprintf(message," Iteration \033[1;32m%d.%d\033[1;33m",num_adapt,nbiteration);
204 francois 239 affiche(message);
205 francois 468 std::vector<FEM_ELEMENT3*> *lstniveau=new std::vector<FEM_ELEMENT3*>[niveaumax+1];
206 francois 247 double densmin = params.get_valeur("ro_void");
207 francois 239 for (int i=0;i<nbsimptet;i++)
208     {
209     SIMP_TETRA* tet=lsttet[i];
210     tet->maille_niveau=(int)((tet->densite-densmin)*niveaumax/(1.-densmin))+1;
211     if (tet->maille_niveau>niveaumax) tet->maille_niveau=niveaumax;
212     lstniveau[tet->maille_niveau].insert(lstniveau[tet->maille_niveau].end(),tet->tet);
213     }
214     MG_EXPORT exp;
215 francois 512 //int version_aster=params.get_valeur("version_aster");
216 francois 706 exp.lire_params_aster(nomparamaster);
217 francois 748 double version=exp.aster(geo,fem,nometude,MAGIC::CALCUL_ASTER::OPTIMISATIONTOPOLOGIQUE,(char*)"00000001",p,niveaumax,lstniveau);
218 francois 239 delete [] lstniveau;
219 francois 581 char messagever[255];
220 francois 691 sprintf(messagever," Calcul aster \033[1;32m%.1lf\033[1;33m",version);
221 francois 581 affiche(messagever);
222 francois 239 char nomfichiertmp[255];
223 francois 725 int choix=(int)param.get_valeur("Affiche_Code_Aster");
224 francois 726 #ifdef CHOIXTERMASTER
225 francois 725 if (choix!=0) sprintf(nomfichiertmp,"xterm -e \"%s/as_run %s.export \"",getenv("PATHASTER"),nometude);
226     else sprintf(nomfichiertmp,"%s/as_run %s.export 1>aster.log 2>&1",getenv("PATHASTER"),nometude);
227 francois 726 #else
228     sprintf(nomfichiertmp,"%s/as_run %s.export 1>aster.log 2>&1",getenv("PATHASTER"),nometude);
229     #endif
230 francois 242 int code=system(nomfichiertmp);
231     if (code!=0)
232     {
233     sprintf(nomfichiertmp," Code de sortie aster : %d",code);
234     affiche(nomfichiertmp);
235     }
236 francois 272 affiche((char*)" Analyse resultat");
237 francois 239 recupere_energie(lsttet);
238     double Ctot=0;
239     for (int i=0;i<nbsimptet;i++)
240     {
241     SIMP_TETRA* tet=lsttet[i];
242 francois 247 tet->denergie= -p*2.*tet->energie/tet->densite;
243 francois 239 Ctot=Ctot+2.*tet->energie;
244     }
245     Ctotiter.insert(Ctotiter.end(),Ctot);
246 francois 691 sprintf(message," Compliance \033[1;31m%le\033[1;33m",Ctot);
247 francois 239 affiche(message);
248     if (nbiteration!=1)
249     {
250     int nb=Ctotiter.size();
251     double c1=Ctotiter[nb-2];
252     double c2=Ctotiter[nb-1];
253     double ecart=fabs((c2-c1)/c1)*100.;
254     if (ecart<params.get_valeur("critere_convergence")) ok=1;
255 francois 691 sprintf(message," Ecart \033[1;31m%lf%%\033[1;33m",ecart);
256 francois 239 affiche(message);
257     }
258 francois 247 // lissage compliance
259     double kc=params.get_valeur("kc");
260 francois 243 int type_lissage =(int)params.get_valeur("type_lissage");
261 francois 247 int lissage_densite =(int)params.get_valeur("lissage_densite");
262     int type_lissage_densite =(int)params.get_valeur("type_lissage_densite");
263 francois 267 int type_lissage_compliance=(int)params.get_valeur("type_lissage_compliance");
264 francois 242 for (int i=0;i<nbsimptet;i++)
265     {
266     SIMP_TETRA* tet=lsttet[i];
267 francois 247 if ((type_lissage==0) || (type_lissage==2))
268 francois 243 tet->new_denergie=tet->denergie;
269 francois 247 if ((type_lissage==1) || (type_lissage==3))
270 francois 242 {
271 francois 243 if (tet->design == 1)
272 francois 242 {
273 francois 247 double hi= pow(tet->distance_ref,kc);
274     double hisensi = hi*tet->densite*tet->denergie;
275 francois 267 if (type_lissage_compliance==1)
276     {
277     hi=hi/tet->volume;
278 francois 268 hisensi=hisensi/max(volumemoyen*0.00001,tet->volume);
279 francois 267 }
280     int nbvoisin=tet->voisin.size();
281 francois 243 for (int j = 0 ; j<nbvoisin ; j++)
282     {
283     SIMP_TETRA* tet2=tet->voisin[j];
284     if (tet2->design == 1)
285     {
286 francois 247 double dist=tet->distance(tet2);
287     double hj =fabs(pow(tet->distance_ref - dist,kc));
288 francois 267 if (type_lissage_compliance==0)
289     {
290     hisensi=hisensi+tet2->densite*hj*tet2->denergie;
291     hi=hi+hj;
292     }
293     if (type_lissage_compliance==1)
294     {
295 francois 268 hisensi=hisensi+tet2->densite*hj*tet2->denergie/max(volumemoyen,tet2->volume);
296 francois 267 hi=hi+hj/tet2->volume;
297     }
298     hisensi=hisensi+tet2->densite*hj*tet2->denergie;
299 francois 247 hi=hi+hj;
300 francois 243 }
301 francois 242 }
302 francois 268 tet->new_denergie = hisensi/hi/max(tet->densite,densmin);
303 francois 243 }
304     else tet->new_denergie=tet->denergie;
305 francois 242 }
306     }
307 francois 239 // application formule pas encore connue
308     double l1= 1e-8 ;
309     double l2= 1e8;
310     double m=params.get_valeur("m");
311 francois 241 double eta=params.get_valeur("eta");
312 francois 242 double convergence_lagrange=params.get_valeur("convergence_lagrange");
313     double critere_densite = 0.0;
314     int oklagrange=0;
315     int compteurlagrange=0;
316     while (oklagrange==0)
317 francois 239 {
318 francois 242 compteurlagrange++;
319     double lmid = 0.5*(l2+l1);
320     critere_densite = 0.0;
321     for (int i=0;i<nbsimptet;i++)
322 francois 239 {
323     SIMP_TETRA* tet=lsttet[i];
324     if (tet->design==1)
325     {
326     double x1 = tet->densite+m;
327 francois 242 double beta=-tet->new_denergie/lmid/tet->volume;
328 francois 241 double x2 = tet->densite*pow(beta,eta);
329     double x3 = min(x1,x2);
330     double x4 = 1.;
331     double x5 = min(x3,x4);
332     double x6 = tet->densite-m;
333     double x7 = max(x5,x6);
334     double x8 = densmin;
335 francois 239 tet->new_densite = max(x7,x8);
336 francois 242 critere_densite = critere_densite + tet->new_densite*tet->volume;
337     }
338 francois 239 else
339     {
340     tet->new_densite=1.;
341     }
342 francois 242 }
343 francois 247 if (critere_densite - f*volume_design > 0.)
344 francois 239 l1=lmid;
345     else
346     l2=lmid;
347 francois 247 if (100.*fabs(critere_densite- f*volume_design )/(f*volume_design)<convergence_lagrange) oklagrange=1;
348 francois 243 if (compteurlagrange>500) oklagrange=2;
349 francois 242 }
350 francois 691 if (oklagrange==1) sprintf(message," Convergence multiplicateur lagrange \033[1;31m%le%%\033[1;33m",100.*(critere_densite- f*volume_design )/(f*volume_design));
351     if (oklagrange==2) sprintf(message," Divergence multiplicateur lagrange \033[1;31m%le %le\033[1;33m",l1,l2);
352 francois 242 affiche(message);
353 francois 247 // lissage de densite.
354     double kd=params.get_valeur("kd");
355 francois 239 for (int i=0;i<nbsimptet;i++)
356 francois 242 {
357 francois 239 SIMP_TETRA* tet=lsttet[i];
358 francois 247 if ((type_lissage==2) || (type_lissage==3))
359     {
360     if (((lissage_densite==1)&&(ok==1)) || (lissage_densite==0))
361     {
362     if (tet->design == 1)
363     {
364     double widensite=0.;
365     double wi= 0.;
366     wi=poid_lissage(0.,tet->distance_ref2,kd,tet->volume,type_lissage_densite);
367     widensite = wi*tet->new_densite;
368     int nbvoisin=tet->voisin2.size();
369     for (int j = 0 ; j<nbvoisin ; j++)
370     {
371     SIMP_TETRA* tet2=tet->voisin2[j];
372     if (tet2->design == 1)
373     {
374     double dist=tet->distance(tet2);
375     double wj=poid_lissage(dist,tet->distance_ref2,kd,tet2->volume,type_lissage_densite);
376     wi = wi+wj;
377     widensite = widensite + tet2->new_densite*wj;
378     }
379     }
380     tet->densite = widensite/wi;
381     }
382     else tet->densite=tet->new_densite;
383     }
384     else tet->densite=tet->new_densite;
385     }
386     else tet->densite=tet->new_densite;
387 francois 242 }
388 francois 339 if (itervue!=0)
389     if (nbiteration%itervue==0)
390     {
391     affiche((char*)" Sauvegarde résultat iteration");
392 cuillier 341 sprintf(message,"%s_densite_iter%d.soltmp",nometudesortie,nbiteration);
393 francois 339 //int nbsolution=gestd->get_nb_fem_solution();
394 nana 565 //for (int i=nbsolution;i>0;i--version 2 du for)
395 francois 339 //gestd->supprimer_fem_solution_du_gestionnaire(i-1);
396 francois 342 char message2[50];
397 francois 691 static int step=0;
398     sprintf(message2,"Iteration%%%d",step);
399     step++;
400 francois 375 FEM_SOLUTION* solution=new FEM_SOLUTION(fem,1,message,fem->get_nb_fem_element3(),message2,MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT3);
401 francois 339 gestd->ajouter_fem_solution(solution);
402     sprintf(message,"Densite");
403     solution->change_legende(0,message);
404     for (int i=0;i<nbsimptet;i++)
405     {
406     SIMP_TETRA* tet=lsttet[i];
407 francois 375 solution->ecrire(tet->densite,i,0);
408 francois 339 }
409     }
410 francois 239 nbiteration++;
411     }
412 francois 339 if (itervue!=0)
413 nana 565 {
414 francois 339 affiche((char*)" Sauvegarde GMSH résultats iterations");
415     MG_EXPORT exp;
416     char nomfichier[500];
417     sprintf(nomfichier,"%s_iterations",nometudesortie);
418     exp.gmsh(fem,nomfichier);
419     }
420 francois 247 double seuil=params.get_valeur("ro_min");
421 francois 272 affiche((char*)" Ecriture des donnees finales");
422 francois 239 char message[255];
423 francois 512 if (iter==0) sprintf(message,"%s.compliance",nometudesortie);
424     else sprintf(message,"%s_iter%d.compliance",nometudesortie,iter);
425 francois 239 FILE *out=fopen(message,"wt");
426 francois 271 time_t heurefin=time(NULL);
427     struct tm tfin = *localtime(&heurefin);
428     fprintf(out,"*****************************************************\n");
429     fprintf(out," Optimisateur de topologie\n");
430     fprintf(out," ERICCA - UQTR\n");
431     fprintf(out,"*****************************************************\n");
432     fprintf(out," Etude : %s\n",nometude);
433 francois 274 fprintf(out," Date : %02u-%02u-%04u\n", tdebut.tm_mday, tdebut.tm_mon+1, 1900 + tdebut.tm_year);
434 francois 271 fprintf(out," Heure debut : %02u:%02u:%02u\n", tdebut.tm_hour, tdebut.tm_min, tdebut.tm_sec);
435     fprintf(out," Heure fin : %02u:%02u:%02u\n", tfin.tm_hour, tfin.tm_min, tfin.tm_sec);
436     fprintf(out," Parametres : \n");
437     int nbparam=params.get_nb();
438     for (int i=0;i<nbparam;i++)
439     {
440 francois 397 if (params.get_type(i)==OT_PARAMETRES::DOUBLE) fprintf(out," %s = %lf\n",params.get_nom(i).c_str(),params.get_valeur(i));
441     if (params.get_type(i)==OT_PARAMETRES::STRING) fprintf(out," %s = %s\n",params.get_nom(i).c_str(),params.get_nom(params.get_nom(i).c_str()).c_str());
442 francois 271 }
443     fprintf(out,"*****************************************************\n");
444     fprintf(out," Compliance après chaque itération\n");
445     fprintf(out,"*****************************************************\n");
446 francois 239 for (int i=0;i<Ctotiter.size();i++)
447     fprintf(out,"%le\n",Ctotiter[i]);
448 francois 355 fprintf(out,"*****************************************************\n");
449     fprintf(out," Variation de la compliance après chaque itération\n");
450     fprintf(out,"*****************************************************\n");
451     for (int i=1;i<Ctotiter.size();i++)
452     fprintf(out,"%lf%%\n",(Ctotiter[i]-Ctotiter[i-1])*100./Ctotiter[i-1]);
453 francois 245 double crit1=0.;
454     double crit2=0.;
455     double crit3=0.;
456     double critv1[10]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
457     double critv2[10]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
458     for (int i=0;i<nbsimptet;i++)
459     {
460     SIMP_TETRA* tet=lsttet[i];
461     if (tet->design==1)
462     {
463     crit1=crit1+tet->densite*tet->volume;
464     if (tet->densite>seuil)
465     {
466     crit2=crit2+tet->densite*tet->volume;
467     crit3=crit3+tet->volume;
468     }
469     for (int j=1;j<10;j++)
470     if (tet->densite>0.1*j)
471     {
472     critv1[j]=critv1[j]+tet->densite*tet->volume;
473     critv2[j]=critv2[j]+tet->volume;
474     }
475     }
476     }
477 francois 271 fprintf(out,"*****************************************************\n");
478     fprintf(out," Résultat\n");
479     fprintf(out,"*****************************************************\n");
480 francois 245 fprintf(out,"Volume du design : %le \n",volume_design);
481 francois 247 fprintf(out,"Objectif du volume de design : %le \n",volume_design*f);
482 francois 245 fprintf(out,"Volume de design obtenu : %le \n",crit1);
483     fprintf(out,"Volume de design obtenu avec le seuil: %le \n",crit2);
484     fprintf(out,"Volume reel de design obtenu avec le seuil: %le \n\n",crit3);
485     fprintf(out," : Volume : Volume reel : Objectif\n");
486     fprintf(out," : %le : %le : %le\n",volume_design,volume_design,volume_design);
487     for (int j=1;j<10;j++)
488 francois 247 fprintf(out,"Volume de design obtenu avec le seuil de %.1f : %le : %le : %le \n",j*0.1,critv1[j],critv2[j],volume_design*f);
489 francois 245 fprintf(out,"*************************************************************\n");
490     fclose(out);
491 francois 239 LISTE_FEM_NOEUD::iterator itnoeud;
492     int nbfemnoeud=fem->get_nb_fem_noeud();
493 francois 242 /*double *nume=new double[nbfemnoeud];
494 francois 239 double *deno=new double[nbfemnoeud];
495     int cpt=0;
496     for (FEM_NOEUD *nd=fem->get_premier_noeud(itnoeud);nd!=NULL;nd=fem->get_suivant_noeud(itnoeud))
497     {
498     nd->change_numero(cpt);
499     nume[cpt]=0.;
500     deno[cpt]=0.;
501     cpt++;
502 francois 242 }*/
503 francois 512 if (iter==0) sprintf(message,"%s_densite.sol",nometudesortie);
504     else sprintf(message,"%s_densite_iter%d.sol",nometudesortie,iter);
505 francois 241 int nbsolution=gestd->get_nb_fem_solution();
506 francois 242 for (int i=nbsolution;i>0;i--)
507 cuillier 341 {
508 francois 342 FEM_SOLUTION* sol=gestd->get_fem_solution(i-1);
509 cuillier 341 std::string nom=sol->get_nom_fichier();
510     char message[500],extension[500];
511     sprintf(message,"%s",nom.c_str());
512     char *p=strrchr(message,'.');
513     strncpy(extension,p,strlen(message)+message-p);
514     if (strcmp(extension,".soltmp")==0)
515     gestd->supprimer_fem_solution(i-1);
516     else gestd->supprimer_fem_solution_du_gestionnaire(i-1);
517     }
518    
519    
520 francois 375 FEM_SOLUTION* solution=new FEM_SOLUTION(fem,1,message,fem->get_nb_fem_element3(),"Optimisation",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT3);
521 francois 239 gestd->ajouter_fem_solution(solution);
522     solution->change_legende(0,"Densite");
523     for (int i=0;i<nbsimptet;i++)
524     {
525     SIMP_TETRA* tet=lsttet[i];
526     if (tet->design==1)
527     if (tet->densite>seuil)
528     ((MG_TETRA*)tet->tet->get_mg_element_maillage())->change_origine(OPTIMISE);
529     else
530     ((MG_TETRA*)tet->tet->get_mg_element_maillage())->change_origine(MAILLEUR_AUTO);
531 francois 384 tet->tet->change_solution(tet->densite);
532 francois 375 solution->ecrire(tet->densite,i,0);
533 francois 242 /*for (int j=0;j<tet->tet->get_nb_fem_noeud();j++)
534 francois 239 {
535     FEM_NOEUD* noeud=tet->tet->get_fem_noeud(j);
536     nume[noeud->get_numero()]=nume[noeud->get_numero()]+tet->volume*tet->densite;
537     deno[noeud->get_numero()]=deno[noeud->get_numero()]+tet->volume;
538 francois 242 } */
539 francois 239 }
540 francois 243 /*sprintf(message,"%s_densite2.sol",nometudesortie);
541 francois 239 FEM_SOLUTION* solution2=new FEM_SOLUTION(fem,1,message,fem->get_nb_fem_noeud(),"Optimisation",ENTITE_NOEUD);
542     gestd->ajouter_fem_solution(solution2);
543     solution2->change_legende(0,"Densite");
544     cpt=0;
545     for (FEM_NOEUD *nd=fem->get_premier_noeud(itnoeud);nd!=NULL;nd=fem->get_suivant_noeud(itnoeud))
546     {
547 francois 375 solution2->ecrire(nume[cpt]/deno[cpt],cpt,0,);
548 francois 239 cpt++;
549     }
550     delete [] deno;
551 francois 242 delete [] nume;*/
552 francois 239 int nb=lsttet.size();
553     for (int i=0;i<nb;i++) delete lsttet[i];
554     }
555    
556 francois 247 double MGOPT_SIMP::poid_lissage(double dist,double distref,double k,double volume,int type)
557     {
558     double wi;
559     if (type==0) wi=pow(distref-dist,k);
560     if (type==1) wi=pow(distref-dist,k)*volume;
561     if (type==2) wi=exp(-dist*dist/2./distref/distref/9.)/2./M_PI/(distref/3.);
562     if (type==3) wi=volume*exp(-dist*dist/2./distref/distref/9.)/2./M_PI/(distref/3.);
563     return fabs(wi);
564    
565     }
566 francois 239
567 francois 240 void MGOPT_SIMP::adapte_resultat(char *nomgestd,char *nomparam)
568     {
569     if (nomparam!=NULL) lire_params(nomparam);
570 francois 272 affiche((char*)"");
571     affiche((char*)"*************************");
572     affiche((char*)"Optimisation de topologie");
573     affiche((char*)"*************************");
574     affiche((char*)"");
575     affiche((char*)"");
576     affiche((char*)"Changement du seuil dans les resultats");
577 francois 240 double seuil=params.get_valeur("seuil");
578     gestd=new MG_FILE(nomgestd);
579     FEM_MAILLAGE* fem=gestd->get_fem_maillage(0);
580     FEM_SOLUTION* solution=gestd->get_fem_solution(0);
581     solution->active_solution(0);
582 francois 309 LISTE_FEM_ELEMENT3::iterator it;
583     for (FEM_ELEMENT3 *tet=fem->get_premier_element3(it);tet!=NULL;tet=fem->get_suivant_element3(it))
584 francois 240 {
585     if (((MG_TETRA*)tet->get_mg_element_maillage())->get_origine()!=IMPOSE)
586     if (tet->get_solution()>seuil)
587     ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(OPTIMISE);
588     else
589     ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(MAILLEUR_AUTO);
590    
591     }
592 francois 272 affiche((char*)"Enregistrement");
593 francois 240 gestd->enregistrer(nomgestd);
594 francois 272 affiche((char*)"Enregistrement sous GMSH");
595 francois 245 char *p=strchr(nomgestd,'.');
596     strncpy(nometude,nomgestd,p-nomgestd);
597     nometude[p-nomgestd]=0;
598     MG_EXPORT exp;
599     char nomfichier[500];
600     sprintf(nomfichier,"%s_mg",nometude);
601     exp.gmsh(fem->get_mg_maillage(),nomfichier);
602     sprintf(nomfichier,"%s_fem",nometude);
603     exp.gmsh(fem,nomfichier);
604 francois 272 affiche((char*)"Fin");
605 francois 240 }
606 francois 239
607    
608 francois 468 void MGOPT_SIMP::recupere_energie(std::vector<class SIMP_TETRA*> lsttet)
609 francois 239 {
610     char message[750];
611     sprintf(message,"%s.resu",nometude);
612     FILE* in=fopen(message,"rt");
613     int fin=0;
614     do
615     {
616 francois 276 char* res=fgets(message,750,in);
617 francois 239 if (feof(in)) fin=1;
618     char mot1[100];
619     char mot2[100];
620     char mot3[100];
621     char mot4[100];
622     char mot5[100];
623     char mot6[100];
624     char mot7[100];
625     char mot8[100];
626     char mot9[100];
627     char mot10[100];
628     int numlu=sscanf(message,"%s %s %s %s %s %s %s %s %s %s",mot1,mot2,mot3,mot4,mot5,mot6,mot7,mot8,mot9,mot10);
629 francois 674
630    
631     if (numlu>2)
632     if (strcmp(mot2,"TOTAL_JOB")==0)
633     {
634     double val;
635     sscanf(mot8,"%lf",&val);
636     CODE_ASTER_CPU=CODE_ASTER_CPU+val;
637 francois 727 sscanf(mot10,"%lf",&val);
638     CODE_ASTER_ECOULE=CODE_ASTER_ECOULE+val;
639 francois 674 }
640    
641    
642 francois 239 if (numlu>9)
643     if (strcmp(mot1,"CHAMP")==0)
644     if (strcmp(mot2,"PAR")==0)
645     if (strcmp(mot3,"ELEMENT")==0)
646     if (strcmp(mot4,"CONSTANT")==0)
647     if (strcmp(mot5,"SUR")==0)
648 francois 322 if (strcmp(mot6,"ELEMENT")==0)
649 francois 239 if (strcmp(mot7,"DE")==0)
650     if (strcmp(mot8,"NOM")==0)
651     if (strcmp(mot9,"SYMBOLIQUE")==0)
652 francois 322 if (strcmp(mot10,"ENERGIE")==0)
653 francois 239 {
654     int fin2=0;
655     int passe=0;
656     int nbelement=0;
657     do
658     {
659     char message[750];
660 francois 276 char* res=fgets(message,750,in);
661 francois 239 char mot1[500];
662     char mot2[500];
663     int numlu=sscanf(message,"%s %s",mot1,mot2);
664     int decalage;
665     if ((numlu==2) && (strcmp(mot2,"TOTALE")==0))
666     {
667     char *p=strchr(mot1,'M')+1;
668     int num=atoi(p);
669     if (passe==0) {passe=1;decalage=num;}
670 francois 276 char* res=fgets(message,750,in);
671 francois 239 double val;
672     sscanf(message,"%lf",&val);
673     lsttet[num-decalage]->energie=val;
674     nbelement++;
675     }
676     if (nbelement == lsttet.size()) {fin2=1;fin=0;}
677     }
678    
679     while (fin2==0);
680     }
681     }
682     while (fin==0);
683     fclose(in);
684     }
685    
686    
687    
688 nana 565 void MGOPT_SIMP::ajouter_voisin_distance(int i,SIMP_TETRA* tet,std::vector<SIMP_TETRA*> &lst)
689 francois 239 {
690     tet->indice=i;
691     int nbnoeud=tet->tet->get_nb_fem_noeud();
692     int correspondance[4];
693     if (nbnoeud==4)
694     {
695     correspondance[0]=0;
696     correspondance[1]=1;
697     correspondance[2]=2;
698     correspondance[3]=3;
699     }
700     if (nbnoeud==10)
701     {
702     correspondance[0]=0;
703     correspondance[1]=2;
704     correspondance[2]=4;
705     correspondance[3]=9;
706     }
707     SIMP_TETRA* tetcour=tet;
708     int ok=0;
709     int compteur=0;
710 francois 468 std::vector<SIMP_TETRA*>& lstvoisin=tet->voisin;
711 francois 247 if (tet->distance_ref<tet->distance_ref2)
712     lstvoisin=tet->voisin2;
713 francois 239 while (ok==0)
714     {
715     for (int j=0;j<4;j++)
716     {
717     int num=correspondance[j];
718     FEM_NOEUD *noeud=tetcour->tet->get_fem_noeud(num);
719 francois 309 int nbtetra=noeud->get_lien_element3()->get_nb();
720 francois 239 for (int k=0;k<nbtetra;k++)
721     {
722 francois 309 FEM_ELEMENT3* ftet=noeud->get_lien_element3()->get(k);
723 francois 239 SIMP_TETRA* stet=lst[ftet->get_numero()];
724     if (stet->indice!=i)
725     {
726     stet->indice=i;
727 francois 247 double dist=tet->distance(stet);
728     if (dist<tet->distance_ref)
729 francois 239 tet->voisin.insert(tet->voisin.end(),stet);
730 francois 247 if (dist<tet->distance_ref2)
731     tet->voisin2.insert(tet->voisin2.end(),stet);
732 francois 239 }
733     }
734     }
735 francois 247 if (compteur>=lstvoisin.size()) ok=1;
736 francois 239 else
737     {
738 francois 247 tetcour=lstvoisin[compteur];
739 francois 239 compteur++;
740     }
741 francois 247
742 francois 239 }
743    
744 francois 274 }
745 nana 565
746    
747 francois 578 void MGOPT_SIMP::ajouter_voisin_couche(int i,SIMP_TETRA* tet,std::vector<SIMP_TETRA*> &lst,int nb_couches,int nb_couches2)
748 nana 565 {
749     tet->indice=i;
750     int nbnoeud=tet->tet->get_nb_fem_noeud();
751     int correspondance[4];
752     if (nbnoeud==4)
753     {
754     correspondance[0]=0;
755     correspondance[1]=1;
756     correspondance[2]=2;
757     correspondance[3]=3;
758     }
759     if (nbnoeud==10)
760     {
761     correspondance[0]=0;
762     correspondance[1]=2;
763     correspondance[2]=4;
764     correspondance[3]=9;
765     }
766     std::vector<SIMP_TETRA*>& lstvoisin=tet->voisin;
767     if (tet->distance_ref<tet->distance_ref2)
768     lstvoisin=tet->voisin2;
769 francois 578 tet->voisin.insert(tet->voisin.end(),tet);
770     tet->voisin2.insert(tet->voisin2.end(),tet);
771     tet->indice=i;
772     int debut=0,fin=1;
773     for (int ok=0;ok<max(nb_couches,nb_couches2);ok++)
774 nana 565 {
775 francois 578 for (int iliste=debut;iliste<fin;iliste++)
776     {
777 francois 579 SIMP_TETRA* tetcour;
778     if (nb_couches>nb_couches2 ) tetcour=tet->voisin[iliste]; else tetcour=tet->voisin2[iliste];
779 francois 578 for (int j=0;j<4;j++)
780 nana 565 {
781     int num=correspondance[j];
782     FEM_NOEUD *noeud=tetcour->tet->get_fem_noeud(num);
783     int nbtetra=noeud->get_lien_element3()->get_nb();
784     for (int k=0;k<nbtetra;k++)
785     {
786     FEM_ELEMENT3* ftet=noeud->get_lien_element3()->get(k);
787     SIMP_TETRA* stet=lst[ftet->get_numero()];
788     if (stet->indice!=i)
789     {
790     stet->indice=i;
791 francois 578 if (ok<nb_couches) tet->voisin.insert(tet->voisin.end(),stet);
792     if (ok<nb_couches2) tet->voisin2.insert(tet->voisin2.end(),stet);
793 nana 565 }
794     }
795     }
796 francois 578 }
797     debut=fin;
798     if (nb_couches>nb_couches2 ) fin=tet->voisin.size(); else fin=tet->voisin2.size();
799 nana 565 }
800     }