ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mgopt_simp.cpp
Revision: 240
Committed: Mon Mar 15 18:23:39 2010 UTC (15 years, 2 months ago) by francois
Original Path: magic/lib/optimisation/optimisation/src/mgopt_simp.cpp
File size: 11923 byte(s)
Log Message:
Ajustement de la derniere version en parametrant le chnagement de seuil

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