ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/aster/src/mgopt_simp.cpp
Revision: 731
Committed: Mon Sep 21 22:31:27 2015 UTC (9 years, 7 months ago) by francois
File size: 26933 byte(s)
Log Message: 
Correction des periodicites dans opencascade
Correction des importation de surface et courbe dans opencascade en respectant tous nos besoins
Ajout d'un langage script pour construire des arbres de construction sous opencascade avec un cas test

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 319 MG_VOLUME* vol=fem->get_mg_geometrie()->get_mg_volume(0);
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 581 double version=exp.aster(vol,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     }