ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mgopt_simp.cpp
Revision: 243
Committed: Wed Mar 24 13:37:21 2010 UTC (15 years, 1 month ago) by francois
Original Path: magic/lib/optimisation/optimisation/src/mgopt_simp.cpp
File size: 14473 byte(s)
Log Message:
Parametrage lissage plus correction bug lissage plus optimisation de code dans la creation des maillages FEM

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