ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/aster/src/mgopt_simp.cpp
Revision: 565
Committed: Tue Oct 14 21:29:43 2014 UTC (10 years, 7 months ago) by nana
File size: 25364 byte(s)
Log Message:
mise  a jour SIMP adaptative etape 2

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