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