ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mgopt_simp.cpp
Revision: 433
Committed: Thu Oct 17 14:04:40 2013 UTC (11 years, 6 months ago) by francois
File size: 23277 byte(s)
Log Message:
Magic V4. Nouvelle approche pour les exe (suite)

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