ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mgopt_simp.cpp
Revision: 468
Committed: Wed Nov 27 03:09:54 2013 UTC (11 years, 5 months ago) by francois
File size: 23570 byte(s)
Log Message:
adaptation de l optimisation de topologie a la nouvelle version de carte

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