ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mgopt_simp.cpp
Revision: 246
Committed: Thu Apr 1 22:00:04 2010 UTC (15 years, 1 month ago) by francois
Original Path: magic/lib/optimisation/optimisation/src/mgopt_simp.cpp
File size: 16362 byte(s)
Log Message:
Adaptation de la distance pour les voisins dans la methode SIMP

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    
7    
8    
9    
10    
11     class SIMP_TETRA
12     {
13     public:
14     FEM_TETRA* tet;
15     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     double centre[3];
26     BOITE_3D get_boite_3D(void)
27     {
28     return tet->get_boite_3D();
29     };
30     vector<SIMP_TETRA*> voisin;
31     unsigned long get_id(void) {return tet->get_id();};
32     double distance(SIMP_TETRA* tet2)
33     {
34     double dx=centre[0]-tet2->centre[0];
35     double dy=centre[1]-tet2->centre[1];
36     double dz=centre[2]-tet2->centre[2];
37     return sqrt(dx*dx+dy*dy+dz*dz);
38     };
39     };
40    
41    
42    
43    
44    
45    
46    
47    
48    
49    
50     MGOPT_SIMP::MGOPT_SIMP():MGOPT()
51     {
52 francois 243 params.ajouter("frac",0.3,"Fraction volumique de la methode SIMP");
53 francois 246 params.ajouter("rmin",0.,"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");
54     params.ajouter("coefvoisin",1.0,"Coefficient de multiplication de rmin");
55 francois 243 params.ajouter("nbrniveau",150.,"Nombre de matériaux utilisés");
56     params.ajouter("penal",3.,"Coefficient de penalite du module d'Young");
57     params.ajouter("densite_min",0.001,"Densite minimale consideree comme nulle");
58     params.ajouter("m",0.2,"Limite d'evolution de densité entre deux itérations");
59     params.ajouter("eta",0.5,"Coefficent d'affectation du beta");
60     params.ajouter("type_lissage",1.,"0. Pas de lissage 1.Lissage classique de la litterature");
61     params.ajouter("relaxation",0.,"0. Lissage sans relaxation 1. Lissage avec relaxation");
62     params.ajouter("k",1.0,"Ponderation de la distance dans le lissage");
63     params.ajouter("convergence_lagrange",0.1,"Critere de convergence de la recherche du multiplicateur de Lagrange en \%");
64     params.ajouter("iter_max",50.,"Nombre d'iteration maximale dans le processus SIMP");
65     params.ajouter("critere_convergence",0.5,"Critere de convergence de la méthode SIMP en \%");
66     params.ajouter("seuil",0.8,"Seuil de conservation des éléments");
67 francois 239 }
68    
69     MGOPT_SIMP::MGOPT_SIMP(MGOPT_SIMP &mdd):MGOPT(mdd)
70     {
71     }
72    
73    
74     MGOPT_SIMP::~MGOPT_SIMP()
75     {
76     }
77    
78    
79     void MGOPT_SIMP::optimisation(FEM_MAILLAGE* fem)
80     {
81     affiche(" Initialisation");
82     vector<SIMP_TETRA*> lsttet;
83     double frac=params.get_valeur("frac");
84 francois 246 double rmin=params.get_valeur("rmin");
85     if ((rmin<1e-16) && (carte==NULL))
86     {
87     affiche(" Lecture de la carte de taille");
88     carte=new FCT_GENERATEUR_3D<4>;
89     std::string fichiercarte=params.get_nom("fichiercarte");
90     carte->lire((char *)fichiercarte.c_str());
91     }
92 francois 239 double volume_tot;
93     double volume_design=0;
94     double volume_non_design=0;
95     double unite=fem->get_mg_maillage()->get_mg_geometrie()->get_valeur_unite();
96     LISTE_FEM_TETRA::iterator it;
97     int ntet=0;
98     for (FEM_TETRA* tet=fem->get_premier_tetra(it);tet!=NULL;tet=fem->get_suivant_tetra(it))
99     {
100     SIMP_TETRA* tet2=new SIMP_TETRA;
101     tet2->tet=tet;
102     tet2->tet->change_numero(ntet);
103     ntet++;
104     tet2->centre[0]=0.;
105     tet2->centre[1]=0.;
106     tet2->centre[2]=0.;
107     int nbnoeud=tet->get_nb_fem_noeud();
108     for (int i=0;i<nbnoeud;i++)
109     {
110 francois 242 tet2->centre[0]=tet2->centre[0]+tet->get_fem_noeud(i)->get_x()*unite;
111     tet2->centre[1]=tet2->centre[1]+tet->get_fem_noeud(i)->get_y()*unite;
112     tet2->centre[2]=tet2->centre[2]+tet->get_fem_noeud(i)->get_z()*unite;
113 francois 239 }
114     tet2->centre[0]=tet2->centre[0]/nbnoeud;
115     tet2->centre[1]=tet2->centre[1]/nbnoeud;
116     tet2->centre[2]=tet2->centre[2]/nbnoeud;
117     int lig,col;
118     double jac[9],uv[3]={1./4.,1./4.,1./4.};
119 francois 242 tet2->volume=1./6.*tet->get_jacobien(jac,uv,lig,col,unite);
120 francois 246 tet2->distance_ref=rmin;
121     double xyz[3]={tet2->centre[0],tet2->centre[1],tet2->centre[2]};
122     double val[9];
123     if (tet2->distance_ref<1e-16)
124 francois 242 {
125 francois 246 carte->evaluer(xyz,val);
126     tet2->distance_ref=1./sqrt(val[0]);
127 francois 242 }
128     // tet2->distance_ref=params.get_valeur("coefvoisin")*pow(12.*tet2->volume/sqrt(2),1./3.);
129 francois 246 tet2->distance_ref=params.get_valeur("coefvoisin")*tet2->distance_ref*unite;
130 francois 239 if (((MG_TETRA*)tet->get_mg_element_maillage())->get_origine()!=IMPOSE)
131     {
132     tet2->design=1;
133     tet2->densite=frac;
134     volume_design=volume_design+tet2->volume;
135     }
136     else
137     {
138     tet2->design=0;
139     tet2->densite=1.;
140     volume_non_design=volume_non_design+tet2->volume;
141     }
142     lsttet.insert(lsttet.end(),tet2);
143     }
144     affiche(" Recherche de voisins");
145     int nbsimptet=lsttet.size();
146     for (int i=0;i<nbsimptet;i++)
147     {
148     SIMP_TETRA* tet=lsttet[i];
149     if (tet->design==0) continue;
150     ajouter_voisin(i,tet,lsttet);
151     }
152     int nbiteration=1;
153     int ok=0;
154     int niveaumax=(int)params.get_valeur("nbrniveau");
155     double penal=params.get_valeur("penal");
156     vector<double> Ctotiter;
157     while (ok==0)
158     {
159     char message[255];
160     sprintf(message," Iteration %d",nbiteration);
161     affiche(message);
162     vector<FEM_TETRA*> *lstniveau=new vector<FEM_TETRA*>[niveaumax+1];
163     double densmin = params.get_valeur("densite_min");
164     for (int i=0;i<nbsimptet;i++)
165     {
166     SIMP_TETRA* tet=lsttet[i];
167     tet->maille_niveau=(int)((tet->densite-densmin)*niveaumax/(1.-densmin))+1;
168     if (tet->maille_niveau>niveaumax) tet->maille_niveau=niveaumax;
169     lstniveau[tet->maille_niveau].insert(lstniveau[tet->maille_niveau].end(),tet->tet);
170     }
171     MG_EXPORT exp;
172     exp.aster(fem,nometude,2,"00000001",penal,niveaumax,lstniveau);
173     delete [] lstniveau;
174     affiche(" Calcul aster");
175     char nomfichiertmp[255];
176     sprintf(nomfichiertmp,"%s/as_run %s.export 1>aster.log 2>&1",getenv("PATHASTER"),nometude);
177 francois 242 int code=system(nomfichiertmp);
178     if (code!=0)
179     {
180     sprintf(nomfichiertmp," Code de sortie aster : %d",code);
181     affiche(nomfichiertmp);
182     }
183 francois 239 affiche(" Analyse resultat");
184     recupere_energie(lsttet);
185     double Ctot=0;
186     for (int i=0;i<nbsimptet;i++)
187     {
188     SIMP_TETRA* tet=lsttet[i];
189     tet->denergie= -penal*2.*tet->energie/tet->densite;
190     Ctot=Ctot+2.*tet->energie;
191     }
192     Ctotiter.insert(Ctotiter.end(),Ctot);
193     sprintf(message," Compliance %le",Ctot);
194     affiche(message);
195     if (nbiteration!=1)
196     {
197     int nb=Ctotiter.size();
198     double c1=Ctotiter[nb-2];
199     double c2=Ctotiter[nb-1];
200     double ecart=fabs((c2-c1)/c1)*100.;
201     if (ecart<params.get_valeur("critere_convergence")) ok=1;
202     if (nbiteration>params.get_valeur("iter_max")) ok=1;
203     sprintf(message," Ecart %lf%%",ecart);
204     affiche(message);
205     }
206 francois 242 // lissage
207     double k=params.get_valeur("k");
208 francois 243 int relaxation =(int)params.get_valeur("relaxation");
209     int type_lissage =(int)params.get_valeur("type_lissage");
210 francois 242 for (int i=0;i<nbsimptet;i++)
211     {
212     SIMP_TETRA* tet=lsttet[i];
213 francois 243 if (type_lissage==0)
214     tet->new_denergie=tet->denergie;
215     if (type_lissage==1)
216 francois 242 {
217 francois 243 if (tet->design == 1)
218 francois 242 {
219 francois 243 double somme_Hv_fois_derivecompliance=0.;
220     double somme_Hv = 0.;
221     if (relaxation==1)
222 francois 242 {
223 francois 243 somme_Hv_fois_derivecompliance = tet->densite*pow(tet->distance_ref,k)*tet->denergie;
224     somme_Hv = pow(tet->distance_ref,k);
225     }
226     int nbvoisin=tet->voisin.size();
227     for (int j = 0 ; j<nbvoisin ; j++)
228     {
229     SIMP_TETRA* tet2=tet->voisin[j];
230     if (tet2->design == 1)
231     {
232     double m_distance=tet->distance(tet2);
233     double Hv =fabs(pow(tet->distance_ref - m_distance,k));
234     somme_Hv_fois_derivecompliance = somme_Hv_fois_derivecompliance + tet2->densite*Hv*tet2->denergie;
235     somme_Hv = somme_Hv + Hv;
236     }
237 francois 242 }
238 francois 243 tet->new_denergie = somme_Hv_fois_derivecompliance/somme_Hv/tet->densite;
239     }
240     else tet->new_denergie=tet->denergie;
241 francois 242 }
242     }
243 francois 239 // application formule pas encore connue
244     double l1= 1e-8 ;
245     double l2= 1e8;
246     double m=params.get_valeur("m");
247 francois 241 double eta=params.get_valeur("eta");
248 francois 242 double convergence_lagrange=params.get_valeur("convergence_lagrange");
249     double critere_densite = 0.0;
250     int oklagrange=0;
251     int compteurlagrange=0;
252     while (oklagrange==0)
253 francois 239 {
254 francois 242 compteurlagrange++;
255     double lmid = 0.5*(l2+l1);
256     critere_densite = 0.0;
257     for (int i=0;i<nbsimptet;i++)
258 francois 239 {
259     SIMP_TETRA* tet=lsttet[i];
260     if (tet->design==1)
261     {
262     double x1 = tet->densite+m;
263 francois 242 double beta=-tet->new_denergie/lmid/tet->volume;
264 francois 241 double x2 = tet->densite*pow(beta,eta);
265     double x3 = min(x1,x2);
266     double x4 = 1.;
267     double x5 = min(x3,x4);
268     double x6 = tet->densite-m;
269     double x7 = max(x5,x6);
270     double x8 = densmin;
271 francois 239 tet->new_densite = max(x7,x8);
272 francois 242 critere_densite = critere_densite + tet->new_densite*tet->volume;
273     }
274 francois 239 else
275     {
276     tet->new_densite=1.;
277     }
278 francois 242 }
279     if (critere_densite - frac*volume_design > 0.)
280 francois 239 l1=lmid;
281     else
282     l2=lmid;
283 francois 242 if (100.*fabs(critere_densite- frac*volume_design )/(frac*volume_design)<convergence_lagrange) oklagrange=1;
284 francois 243 if (compteurlagrange>500) oklagrange=2;
285 francois 242 }
286     if (oklagrange==1) sprintf(message," Convergence multiplicateur lagrange %le%%",100.*(critere_densite- frac*volume_design )/(frac*volume_design));
287     if (oklagrange==2) sprintf(message," Divergence multiplicateur lagrange %le %le",l1,l2);
288     affiche(message);
289 francois 239 for (int i=0;i<nbsimptet;i++)
290 francois 242 {
291 francois 239 SIMP_TETRA* tet=lsttet[i];
292 francois 242 tet->densite=tet->new_densite;
293     }
294 francois 239 nbiteration++;
295     }
296 francois 245 double seuil=params.get_valeur("seuil");
297     affiche(" Ecriture des donnees finales");
298 francois 239 char message[255];
299 francois 244 sprintf(message,"%s.compliance",nometudesortie);
300 francois 239 FILE *out=fopen(message,"wt");
301     for (int i=0;i<Ctotiter.size();i++)
302     fprintf(out,"%le\n",Ctotiter[i]);
303 francois 245 double crit1=0.;
304     double crit2=0.;
305     double crit3=0.;
306     double critv1[10]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
307     double critv2[10]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
308     for (int i=0;i<nbsimptet;i++)
309     {
310     SIMP_TETRA* tet=lsttet[i];
311     if (tet->design==1)
312     {
313     crit1=crit1+tet->densite*tet->volume;
314     if (tet->densite>seuil)
315     {
316     crit2=crit2+tet->densite*tet->volume;
317     crit3=crit3+tet->volume;
318     }
319     for (int j=1;j<10;j++)
320     if (tet->densite>0.1*j)
321     {
322     critv1[j]=critv1[j]+tet->densite*tet->volume;
323     critv2[j]=critv2[j]+tet->volume;
324     }
325     }
326     }
327     fprintf(out,"\n\n\n**************************************************************\n");
328     fprintf(out,"Volume du design : %le \n",volume_design);
329     fprintf(out,"Objectif du volume de design : %le \n",volume_design*frac);
330     fprintf(out,"Volume de design obtenu : %le \n",crit1);
331     fprintf(out,"Volume de design obtenu avec le seuil: %le \n",crit2);
332     fprintf(out,"Volume reel de design obtenu avec le seuil: %le \n\n",crit3);
333     fprintf(out," : Volume : Volume reel : Objectif\n");
334     fprintf(out," : %le : %le : %le\n",volume_design,volume_design,volume_design);
335     for (int j=1;j<10;j++)
336     fprintf(out,"Volume de design obtenu avec le seuil de %.1f : %le : %le : %le \n",j*0.1,critv1[j],critv2[j],volume_design*frac);
337     fprintf(out,"*************************************************************\n");
338     fclose(out);
339 francois 239 LISTE_FEM_NOEUD::iterator itnoeud;
340     int nbfemnoeud=fem->get_nb_fem_noeud();
341 francois 242 /*double *nume=new double[nbfemnoeud];
342 francois 239 double *deno=new double[nbfemnoeud];
343     int cpt=0;
344     for (FEM_NOEUD *nd=fem->get_premier_noeud(itnoeud);nd!=NULL;nd=fem->get_suivant_noeud(itnoeud))
345     {
346     nd->change_numero(cpt);
347     nume[cpt]=0.;
348     deno[cpt]=0.;
349     cpt++;
350 francois 242 }*/
351 francois 243 sprintf(message,"%s_densite1.sol",nometudesortie);
352 francois 241 int nbsolution=gestd->get_nb_fem_solution();
353 francois 242 for (int i=nbsolution;i>0;i--)
354 francois 243 gestd->supprimer_fem_solution_du_gestionnaire(i-1);
355 francois 239 FEM_SOLUTION* solution=new FEM_SOLUTION(fem,1,message,fem->get_nb_fem_tetra(),"Optimisation",ENTITE_TETRA);
356     gestd->ajouter_fem_solution(solution);
357     solution->change_legende(0,"Densite");
358     for (int i=0;i<nbsimptet;i++)
359     {
360     SIMP_TETRA* tet=lsttet[i];
361     if (tet->design==1)
362     if (tet->densite>seuil)
363     ((MG_TETRA*)tet->tet->get_mg_element_maillage())->change_origine(OPTIMISE);
364     else
365     ((MG_TETRA*)tet->tet->get_mg_element_maillage())->change_origine(MAILLEUR_AUTO);
366     solution->ecrire(i,0,tet->densite);
367 francois 242 /*for (int j=0;j<tet->tet->get_nb_fem_noeud();j++)
368 francois 239 {
369     FEM_NOEUD* noeud=tet->tet->get_fem_noeud(j);
370     nume[noeud->get_numero()]=nume[noeud->get_numero()]+tet->volume*tet->densite;
371     deno[noeud->get_numero()]=deno[noeud->get_numero()]+tet->volume;
372 francois 242 } */
373 francois 239 }
374 francois 243 /*sprintf(message,"%s_densite2.sol",nometudesortie);
375 francois 239 FEM_SOLUTION* solution2=new FEM_SOLUTION(fem,1,message,fem->get_nb_fem_noeud(),"Optimisation",ENTITE_NOEUD);
376     gestd->ajouter_fem_solution(solution2);
377     solution2->change_legende(0,"Densite");
378     cpt=0;
379     for (FEM_NOEUD *nd=fem->get_premier_noeud(itnoeud);nd!=NULL;nd=fem->get_suivant_noeud(itnoeud))
380     {
381     solution2->ecrire(cpt,0,nume[cpt]/deno[cpt]);
382     cpt++;
383     }
384     delete [] deno;
385 francois 242 delete [] nume;*/
386 francois 239 int nb=lsttet.size();
387     for (int i=0;i<nb;i++) delete lsttet[i];
388     }
389    
390    
391 francois 240 void MGOPT_SIMP::adapte_resultat(char *nomgestd,char *nomparam)
392     {
393     if (nomparam!=NULL) lire_params(nomparam);
394     affiche("");
395     affiche("*************************");
396     affiche("Optimisation de topologie");
397     affiche("*************************");
398     affiche("");
399     affiche("");
400     affiche("Changement du seuil dans les resultats");
401     double seuil=params.get_valeur("seuil");
402     gestd=new MG_FILE(nomgestd);
403     FEM_MAILLAGE* fem=gestd->get_fem_maillage(0);
404     FEM_SOLUTION* solution=gestd->get_fem_solution(0);
405     solution->active_solution(0);
406     LISTE_FEM_TETRA::iterator it;
407     for (FEM_TETRA *tet=fem->get_premier_tetra(it);tet!=NULL;tet=fem->get_suivant_tetra(it))
408     {
409     if (((MG_TETRA*)tet->get_mg_element_maillage())->get_origine()!=IMPOSE)
410     if (tet->get_solution()>seuil)
411     ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(OPTIMISE);
412     else
413     ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(MAILLEUR_AUTO);
414    
415     }
416     affiche("Enregistrement");
417     gestd->enregistrer(nomgestd);
418 francois 245 affiche("Enregistrement sous GMSH");
419     char *p=strchr(nomgestd,'.');
420     strncpy(nometude,nomgestd,p-nomgestd);
421     nometude[p-nomgestd]=0;
422     MG_EXPORT exp;
423     char nomfichier[500];
424     sprintf(nomfichier,"%s_mg",nometude);
425     exp.gmsh(fem->get_mg_maillage(),nomfichier);
426     sprintf(nomfichier,"%s_fem",nometude);
427     exp.gmsh(fem,nomfichier);
428 francois 240 affiche("Fin");
429     }
430 francois 239
431    
432     void MGOPT_SIMP::recupere_energie(vector<class SIMP_TETRA*> lsttet)
433     {
434     char message[750];
435     sprintf(message,"%s.resu",nometude);
436     FILE* in=fopen(message,"rt");
437     int fin=0;
438     do
439     {
440     fgets(message,750,in);
441     if (feof(in)) fin=1;
442     char mot1[100];
443     char mot2[100];
444     char mot3[100];
445     char mot4[100];
446     char mot5[100];
447     char mot6[100];
448     char mot7[100];
449     char mot8[100];
450     char mot9[100];
451     char mot10[100];
452     int numlu=sscanf(message,"%s %s %s %s %s %s %s %s %s %s",mot1,mot2,mot3,mot4,mot5,mot6,mot7,mot8,mot9,mot10);
453     if (numlu>9)
454     if (strcmp(mot1,"CHAMP")==0)
455     if (strcmp(mot2,"PAR")==0)
456     if (strcmp(mot3,"ELEMENT")==0)
457     if (strcmp(mot4,"CONSTANT")==0)
458     if (strcmp(mot5,"SUR")==0)
459     if (strcmp(mot6,"L'ELEMENT")==0)
460     if (strcmp(mot7,"DE")==0)
461     if (strcmp(mot8,"NOM")==0)
462     if (strcmp(mot9,"SYMBOLIQUE")==0)
463     if (strcmp(mot10,"EPOT_ELEM_DEPL")==0)
464     {
465     int fin2=0;
466     int passe=0;
467     int nbelement=0;
468     do
469     {
470     char message[750];
471     fgets(message,750,in);
472     char mot1[500];
473     char mot2[500];
474     int numlu=sscanf(message,"%s %s",mot1,mot2);
475     int decalage;
476     if ((numlu==2) && (strcmp(mot2,"TOTALE")==0))
477     {
478     char *p=strchr(mot1,'M')+1;
479     int num=atoi(p);
480     if (passe==0) {passe=1;decalage=num;}
481     fgets(message,750,in);
482     double val;
483     sscanf(message,"%lf",&val);
484     lsttet[num-decalage]->energie=val;
485     nbelement++;
486     }
487     if (nbelement == lsttet.size()) {fin2=1;fin=0;}
488     }
489    
490     while (fin2==0);
491     }
492     }
493     while (fin==0);
494     fclose(in);
495     }
496    
497    
498    
499     void MGOPT_SIMP::ajouter_voisin(int i,SIMP_TETRA* tet,vector<SIMP_TETRA*> &lst)
500     {
501     tet->indice=i;
502     int nbnoeud=tet->tet->get_nb_fem_noeud();
503     int correspondance[4];
504     if (nbnoeud==4)
505     {
506     correspondance[0]=0;
507     correspondance[1]=1;
508     correspondance[2]=2;
509     correspondance[3]=3;
510     }
511     if (nbnoeud==10)
512     {
513     correspondance[0]=0;
514     correspondance[1]=2;
515     correspondance[2]=4;
516     correspondance[3]=9;
517     }
518     SIMP_TETRA* tetcour=tet;
519     int ok=0;
520     int compteur=0;
521     while (ok==0)
522     {
523     for (int j=0;j<4;j++)
524     {
525     int num=correspondance[j];
526     FEM_NOEUD *noeud=tetcour->tet->get_fem_noeud(num);
527     int nbtetra=noeud->get_lien_tetra()->get_nb();
528     for (int k=0;k<nbtetra;k++)
529     {
530     FEM_TETRA* ftet=noeud->get_lien_tetra()->get(k);
531     SIMP_TETRA* stet=lst[ftet->get_numero()];
532     if (stet->indice!=i)
533     {
534     stet->indice=i;
535     //if (stet->design==1)
536     if (tet->distance(stet)<tet->distance_ref)
537     tet->voisin.insert(tet->voisin.end(),stet);
538     }
539     }
540     }
541     if (compteur>=tet->voisin.size()) ok=1;
542     else
543     {
544     tetcour=tet->voisin[compteur];
545     compteur++;
546     }
547    
548     }
549    
550     }