ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/aster/src/mgopt_simp.cpp
Revision: 242
Committed: Thu Mar 18 15:56:10 2010 UTC (15 years, 1 month ago) by francois
Original Path: magic/lib/optimisation/optimisation/src/mgopt_simp.cpp
File size: 13549 byte(s)
Log Message:
modif pour maillage non structure et non parfait

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