ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mgopt_simp.cpp
Revision: 241
Committed: Mon Mar 15 22:13:49 2010 UTC (15 years, 2 months ago) by francois
Original Path: magic/lib/optimisation/optimisation/src/mgopt_simp.cpp
File size: 12124 byte(s)
Log Message:
parametrage de eta dans la methode SIMP + correction bug affichage dans mgoperation

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     int design;
17     double new_densite;
18     double energie;
19     double denergie;
20     double volume;
21     int maille_niveau;
22     int indice;
23     double distance_ref;
24     double centre[3];
25     BOITE_3D get_boite_3D(void)
26     {
27     return tet->get_boite_3D();
28     };
29     vector<SIMP_TETRA*> voisin;
30     unsigned long get_id(void) {return tet->get_id();};
31     double distance(SIMP_TETRA* tet2)
32     {
33     double dx=centre[0]-tet2->centre[0];
34     double dy=centre[1]-tet2->centre[1];
35     double dz=centre[2]-tet2->centre[2];
36     return sqrt(dx*dx+dy*dy+dz*dz);
37     };
38     };
39    
40    
41    
42    
43    
44    
45    
46    
47    
48    
49     MGOPT_SIMP::MGOPT_SIMP():MGOPT()
50     {
51     params.ajouter("frac",0.3);
52     params.ajouter("coefvoisin",2.5);
53     params.ajouter("nbrniveau",150.);
54     params.ajouter("penal",3.);
55     params.ajouter("densite_min",0.001);
56     params.ajouter("m",0.2);
57 francois 241 params.ajouter("eta",0.5);
58 francois 239 params.ajouter("relaxation",0.);
59     params.ajouter("k",3.0);
60     params.ajouter("iter_max",50.);
61     params.ajouter("critere_convergence",0.5);
62     params.ajouter("seuil",0.8);
63     }
64    
65     MGOPT_SIMP::MGOPT_SIMP(MGOPT_SIMP &mdd):MGOPT(mdd)
66     {
67     }
68    
69    
70     MGOPT_SIMP::~MGOPT_SIMP()
71     {
72     }
73    
74    
75     void MGOPT_SIMP::optimisation(FEM_MAILLAGE* fem)
76     {
77     affiche(" Initialisation");
78     vector<SIMP_TETRA*> lsttet;
79     double frac=params.get_valeur("frac");
80     double volume_tot;
81     double volume_design=0;
82     double volume_non_design=0;
83     double unite=fem->get_mg_maillage()->get_mg_geometrie()->get_valeur_unite();
84     LISTE_FEM_TETRA::iterator it;
85     int ntet=0;
86     for (FEM_TETRA* tet=fem->get_premier_tetra(it);tet!=NULL;tet=fem->get_suivant_tetra(it))
87     {
88     SIMP_TETRA* tet2=new SIMP_TETRA;
89     tet2->tet=tet;
90     tet2->tet->change_numero(ntet);
91     ntet++;
92     tet2->centre[0]=0.;
93     tet2->centre[1]=0.;
94     tet2->centre[2]=0.;
95     int nbnoeud=tet->get_nb_fem_noeud();
96     for (int i=0;i<nbnoeud;i++)
97     {
98     tet2->centre[0]=tet2->centre[0]+tet->get_fem_noeud(i)->get_x();
99     tet2->centre[1]=tet2->centre[1]+tet->get_fem_noeud(i)->get_y();
100     tet2->centre[2]=tet2->centre[2]+tet->get_fem_noeud(i)->get_z();
101     }
102     tet2->centre[0]=tet2->centre[0]/nbnoeud;
103     tet2->centre[1]=tet2->centre[1]/nbnoeud;
104     tet2->centre[2]=tet2->centre[2]/nbnoeud;
105     int lig,col;
106     double jac[9],uv[3]={1./4.,1./4.,1./4.};
107     tet2->volume=1./6.*tet->get_jacobien(jac,uv,lig,col,1.0);
108     tet2->distance_ref=params.get_valeur("coefvoisin")*pow(12.*tet2->volume/sqrt(2),1./3.);
109     if (((MG_TETRA*)tet->get_mg_element_maillage())->get_origine()!=IMPOSE)
110     {
111     tet2->design=1;
112     tet2->densite=frac;
113     volume_design=volume_design+tet2->volume;
114     }
115     else
116     {
117     tet2->design=0;
118     tet2->densite=1.;
119     volume_non_design=volume_non_design+tet2->volume;
120     }
121     lsttet.insert(lsttet.end(),tet2);
122     }
123     affiche(" Recherche de voisins");
124     int nbsimptet=lsttet.size();
125     for (int i=0;i<nbsimptet;i++)
126     {
127     SIMP_TETRA* tet=lsttet[i];
128     if (tet->design==0) continue;
129     ajouter_voisin(i,tet,lsttet);
130     }
131     int nbiteration=1;
132     int ok=0;
133     int niveaumax=(int)params.get_valeur("nbrniveau");
134     double penal=params.get_valeur("penal");
135     vector<double> Ctotiter;
136     while (ok==0)
137     {
138     char message[255];
139     sprintf(message," Iteration %d",nbiteration);
140     affiche(message);
141     vector<FEM_TETRA*> *lstniveau=new vector<FEM_TETRA*>[niveaumax+1];
142     double densmin = params.get_valeur("densite_min");
143     for (int i=0;i<nbsimptet;i++)
144     {
145     SIMP_TETRA* tet=lsttet[i];
146     tet->maille_niveau=(int)((tet->densite-densmin)*niveaumax/(1.-densmin))+1;
147     if (tet->maille_niveau>niveaumax) tet->maille_niveau=niveaumax;
148     lstniveau[tet->maille_niveau].insert(lstniveau[tet->maille_niveau].end(),tet->tet);
149     }
150     MG_EXPORT exp;
151     exp.aster(fem,nometude,2,"00000001",penal,niveaumax,lstniveau);
152     delete [] lstniveau;
153     affiche(" Calcul aster");
154     char nomfichiertmp[255];
155     sprintf(nomfichiertmp,"%s/as_run %s.export 1>aster.log 2>&1",getenv("PATHASTER"),nometude);
156     system(nomfichiertmp);
157     affiche(" Analyse resultat");
158     recupere_energie(lsttet);
159     double Ctot=0;
160     for (int i=0;i<nbsimptet;i++)
161     {
162     SIMP_TETRA* tet=lsttet[i];
163     tet->denergie= -penal*2.*tet->energie/tet->densite;
164     Ctot=Ctot+2.*tet->energie;
165     }
166     Ctotiter.insert(Ctotiter.end(),Ctot);
167     sprintf(message," Compliance %le",Ctot);
168     affiche(message);
169     if (nbiteration!=1)
170     {
171     int nb=Ctotiter.size();
172     double c1=Ctotiter[nb-2];
173     double c2=Ctotiter[nb-1];
174     double ecart=fabs((c2-c1)/c1)*100.;
175     if (ecart<params.get_valeur("critere_convergence")) ok=1;
176     if (nbiteration>params.get_valeur("iter_max")) ok=1;
177     sprintf(message," Ecart %lf%%",ecart);
178     affiche(message);
179     }
180     // application formule pas encore connue
181     double l1= 1e-8 ;
182     double l2= 1e8;
183     double m=params.get_valeur("m");
184 francois 241 double eta=params.get_valeur("eta");
185 francois 239 while ((l2-l1) > 1e-12)
186     {
187     double lmid = 0.5*(l2+l1);
188     double critere_densite = 0.0;
189     for (int i=0;i<nbsimptet;i++)
190     {
191     SIMP_TETRA* tet=lsttet[i];
192     if (tet->design==1)
193     {
194     double x1 = tet->densite+m;
195 francois 241 double beta=-tet->denergie/lmid/tet->volume;
196     double x2 = tet->densite*pow(beta,eta);
197     double x3 = min(x1,x2);
198     double x4 = 1.;
199     double x5 = min(x3,x4);
200     double x6 = tet->densite-m;
201     double x7 = max(x5,x6);
202     double x8 = densmin;
203 francois 239 tet->new_densite = max(x7,x8);
204     }
205     else
206     {
207     tet->new_densite=1.;
208     }
209     critere_densite = critere_densite + tet->new_densite*tet->volume;
210     }
211     if (critere_densite - (frac*volume_design + volume_non_design) > 0)
212     l1=lmid;
213     else
214     l2=lmid;
215     }
216     // lissage
217     double k=params.get_valeur("k");
218     for (int i=0;i<nbsimptet;i++)
219     {
220     SIMP_TETRA* tet=lsttet[i];
221     if (tet->design == 0)
222     {
223     tet->densite = 1.;
224     }
225     else
226     {
227     double somme_Hv_fois_densite = 0.;
228     double somme_Hv = 0.;
229     if ((int)params.get_valeur("relaxation")==1)
230     {
231     somme_Hv_fois_densite = tet->volume*pow(tet->distance_ref,k)*tet->new_densite;
232     somme_Hv = tet->volume*pow(tet->distance_ref,k);
233     }
234     int nbvoisin=tet->voisin.size();
235     for (int j = 0 ; j<nbvoisin ; j++)
236     {
237     SIMP_TETRA* tet2=tet->voisin[j];
238     if (tet2->design == 1)
239     {
240     double m_distance=tet->distance(tet2);
241     double Hv =fabs(tet2->volume*pow(tet->distance_ref - m_distance,k));
242     somme_Hv_fois_densite = somme_Hv_fois_densite + Hv*tet2->new_densite;
243     somme_Hv = somme_Hv + Hv;
244     }
245     }
246     tet->densite = somme_Hv_fois_densite/somme_Hv;
247     tet->densite=max(densmin,tet->densite);
248     }
249     }
250     nbiteration++;
251     }
252     char message[255];
253     sprintf(message,"%s.compliance",nometude);
254     FILE *out=fopen(message,"wt");
255     for (int i=0;i<Ctotiter.size();i++)
256     fprintf(out,"%le\n",Ctotiter[i]);
257     affiche(" Ecriture des donnees finales");
258     LISTE_FEM_NOEUD::iterator itnoeud;
259     int nbfemnoeud=fem->get_nb_fem_noeud();
260     double *nume=new double[nbfemnoeud];
261     double *deno=new double[nbfemnoeud];
262     int cpt=0;
263     for (FEM_NOEUD *nd=fem->get_premier_noeud(itnoeud);nd!=NULL;nd=fem->get_suivant_noeud(itnoeud))
264     {
265     nd->change_numero(cpt);
266     nume[cpt]=0.;
267     deno[cpt]=0.;
268     cpt++;
269     }
270     double seuil=params.get_valeur("seuil");
271     sprintf(message,"%s_densite1.sol",nometude);
272 francois 241 int nbsolution=gestd->get_nb_fem_solution();
273     for (int i=0;i<nbsolution;i++)
274     gestd->supprimer_fem_solution(i);
275 francois 239 FEM_SOLUTION* solution=new FEM_SOLUTION(fem,1,message,fem->get_nb_fem_tetra(),"Optimisation",ENTITE_TETRA);
276     gestd->ajouter_fem_solution(solution);
277     solution->change_legende(0,"Densite");
278     for (int i=0;i<nbsimptet;i++)
279     {
280     SIMP_TETRA* tet=lsttet[i];
281     if (tet->design==1)
282     if (tet->densite>seuil)
283     ((MG_TETRA*)tet->tet->get_mg_element_maillage())->change_origine(OPTIMISE);
284     else
285     ((MG_TETRA*)tet->tet->get_mg_element_maillage())->change_origine(MAILLEUR_AUTO);
286     solution->ecrire(i,0,tet->densite);
287     for (int j=0;j<tet->tet->get_nb_fem_noeud();j++)
288     {
289     FEM_NOEUD* noeud=tet->tet->get_fem_noeud(j);
290     nume[noeud->get_numero()]=nume[noeud->get_numero()]+tet->volume*tet->densite;
291     deno[noeud->get_numero()]=deno[noeud->get_numero()]+tet->volume;
292     }
293     }
294     sprintf(message,"%s_densite2.sol",nometude);
295     FEM_SOLUTION* solution2=new FEM_SOLUTION(fem,1,message,fem->get_nb_fem_noeud(),"Optimisation",ENTITE_NOEUD);
296     gestd->ajouter_fem_solution(solution2);
297     solution2->change_legende(0,"Densite");
298     cpt=0;
299     for (FEM_NOEUD *nd=fem->get_premier_noeud(itnoeud);nd!=NULL;nd=fem->get_suivant_noeud(itnoeud))
300     {
301     solution2->ecrire(cpt,0,nume[cpt]/deno[cpt]);
302     cpt++;
303     }
304     delete [] deno;
305     delete [] nume;
306     int nb=lsttet.size();
307     for (int i=0;i<nb;i++) delete lsttet[i];
308     }
309    
310    
311 francois 240 void MGOPT_SIMP::adapte_resultat(char *nomgestd,char *nomparam)
312     {
313     if (nomparam!=NULL) lire_params(nomparam);
314     affiche("");
315     affiche("*************************");
316     affiche("Optimisation de topologie");
317     affiche("*************************");
318     affiche("");
319     affiche("");
320     affiche("Changement du seuil dans les resultats");
321     double seuil=params.get_valeur("seuil");
322     gestd=new MG_FILE(nomgestd);
323     FEM_MAILLAGE* fem=gestd->get_fem_maillage(0);
324     FEM_SOLUTION* solution=gestd->get_fem_solution(0);
325     solution->active_solution(0);
326     LISTE_FEM_TETRA::iterator it;
327     for (FEM_TETRA *tet=fem->get_premier_tetra(it);tet!=NULL;tet=fem->get_suivant_tetra(it))
328     {
329     if (((MG_TETRA*)tet->get_mg_element_maillage())->get_origine()!=IMPOSE)
330     if (tet->get_solution()>seuil)
331     ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(OPTIMISE);
332     else
333     ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(MAILLEUR_AUTO);
334    
335     }
336     affiche("Enregistrement");
337     gestd->enregistrer(nomgestd);
338     affiche("Fin");
339     }
340 francois 239
341    
342     void MGOPT_SIMP::recupere_energie(vector<class SIMP_TETRA*> lsttet)
343     {
344     char message[750];
345     sprintf(message,"%s.resu",nometude);
346     FILE* in=fopen(message,"rt");
347     int fin=0;
348     do
349     {
350     fgets(message,750,in);
351     if (feof(in)) fin=1;
352     char mot1[100];
353     char mot2[100];
354     char mot3[100];
355     char mot4[100];
356     char mot5[100];
357     char mot6[100];
358     char mot7[100];
359     char mot8[100];
360     char mot9[100];
361     char mot10[100];
362     int numlu=sscanf(message,"%s %s %s %s %s %s %s %s %s %s",mot1,mot2,mot3,mot4,mot5,mot6,mot7,mot8,mot9,mot10);
363     if (numlu>9)
364     if (strcmp(mot1,"CHAMP")==0)
365     if (strcmp(mot2,"PAR")==0)
366     if (strcmp(mot3,"ELEMENT")==0)
367     if (strcmp(mot4,"CONSTANT")==0)
368     if (strcmp(mot5,"SUR")==0)
369     if (strcmp(mot6,"L'ELEMENT")==0)
370     if (strcmp(mot7,"DE")==0)
371     if (strcmp(mot8,"NOM")==0)
372     if (strcmp(mot9,"SYMBOLIQUE")==0)
373     if (strcmp(mot10,"EPOT_ELEM_DEPL")==0)
374     {
375     int fin2=0;
376     int passe=0;
377     int nbelement=0;
378     do
379     {
380     char message[750];
381     fgets(message,750,in);
382     char mot1[500];
383     char mot2[500];
384     int numlu=sscanf(message,"%s %s",mot1,mot2);
385     int decalage;
386     if ((numlu==2) && (strcmp(mot2,"TOTALE")==0))
387     {
388     char *p=strchr(mot1,'M')+1;
389     int num=atoi(p);
390     if (passe==0) {passe=1;decalage=num;}
391     fgets(message,750,in);
392     double val;
393     sscanf(message,"%lf",&val);
394     lsttet[num-decalage]->energie=val;
395     nbelement++;
396     }
397     if (nbelement == lsttet.size()) {fin2=1;fin=0;}
398     }
399    
400     while (fin2==0);
401     }
402     }
403     while (fin==0);
404     fclose(in);
405     }
406    
407    
408    
409     void MGOPT_SIMP::ajouter_voisin(int i,SIMP_TETRA* tet,vector<SIMP_TETRA*> &lst)
410     {
411     tet->indice=i;
412     int nbnoeud=tet->tet->get_nb_fem_noeud();
413     int correspondance[4];
414     if (nbnoeud==4)
415     {
416     correspondance[0]=0;
417     correspondance[1]=1;
418     correspondance[2]=2;
419     correspondance[3]=3;
420     }
421     if (nbnoeud==10)
422     {
423     correspondance[0]=0;
424     correspondance[1]=2;
425     correspondance[2]=4;
426     correspondance[3]=9;
427     }
428     SIMP_TETRA* tetcour=tet;
429     int ok=0;
430     int compteur=0;
431     while (ok==0)
432     {
433     for (int j=0;j<4;j++)
434     {
435     int num=correspondance[j];
436     FEM_NOEUD *noeud=tetcour->tet->get_fem_noeud(num);
437     int nbtetra=noeud->get_lien_tetra()->get_nb();
438     for (int k=0;k<nbtetra;k++)
439     {
440     FEM_TETRA* ftet=noeud->get_lien_tetra()->get(k);
441     SIMP_TETRA* stet=lst[ftet->get_numero()];
442     if (stet->indice!=i)
443     {
444     stet->indice=i;
445     //if (stet->design==1)
446     if (tet->distance(stet)<tet->distance_ref)
447     tet->voisin.insert(tet->voisin.end(),stet);
448     }
449     }
450     }
451     if (compteur>=tet->voisin.size()) ok=1;
452     else
453     {
454     tetcour=tet->voisin[compteur];
455     compteur++;
456     }
457    
458     }
459    
460     }