Jeux de données
top_oc.dgibi
Script Cast3M top_oc.dgibi
1************************************************************************
2* Example of topological optimization *
3* Density method + penalization (SIMP) *
4* *
5* This test case is a reproduction of the reference code: *
6* Andreassen, Clausen, Schevenels, Lazarov, & Sigmund (2011) *
7* "Efficient topology optimization in MATLAB using 88 lines of code" *
8* Structural and Multidisciplinary Optimization, 43(1), 1-16 *
9* https://doi.org/10.1007/s00158-010-0594-7 *
10* *
11* Application to a beam under bending in 2D plane stresses *
12* *
13* Optimization algorithm: Optimality Criteria *
14* *
15* Problem: minimize the compliance with a constraint *
16* on the volume fraction *
17* *
18* min C(x) = u^T.F *
19* s.t. G(x) = vf(x) = volfrac *
20* *
21* | *
22* | Force *
23* v *
24* O>+------------------------------------------------------+ *
25* | p1 | *
26* | | *
27* O>| | *
28* | | *
29* | | *
30* O>+------------------------------------------------------+ p2 *
31* Ʌ *
32* O *
33* *
34************************************************************************
35
36
37** Filtering procedure
38DEBP HFILT cham1*'MCHAML' mo*'MMODEL' mf*'RIGIDITE' mpt*'MAILLAGE' ;
39 chp1 = MANU 'CHPO' mpt 'SCAL' (EXTR cham1 'VALE' 'SCAL') ;
40 chp2 = mf * chp1 ;
41 cham2 = MANU 'CHML' mo 'REPA' 'SCAL' (EXTR chp2 'VALE' mpt) 'TYPE' 'SCALAIRE' 'GRAVITE' ;
42FINP cham2 ;
43
44** Global parameters
45itrac = VRAI ;
46OPTI 'DIME' 2 'MODE' 'PLAN' 'CONT' 'ELEM' 'QUA4' 'ECHO' ;
47
48** Geometrical parameters (width and height)
49l = 60. ;
50h = 20. ;
51e = 1. ;
52
53** Mesh parameters
54nelx = 60 ;
55nely = 20 ;
56
57** Material parameters (isotropic elasticity)
58e0 = 1. ;
59emin = e0 / 1.E9 ;
60nu = 0.3 ;
61
62** Topology optimization parameters
63* penal : SIMP penalization coefficient
64* volfrac : minimal volume fraction
65* ft : = 1 apply filter on the compliance sensitivity field
66* = 2 apply spatial filter on the density field
67* = 3 apply spatial filter on the density field + thresholding by heaviside function
68* rmin : filter radius (in m)
69* beta : initial value of the slope of the heaviside function (only for ft = 3)
70* will be doubled every 50 iterations
71* move : limit of the maximal increment of density
72* changmax : optimization stop criterion
73* xmin/xmax : min and max density bounds
74penal = 3. ;
75volfrac = 0.5 ;
76ft = 2 ;
77rmin = 0.04 * l ;
78beta = 1 ;
79move = 0.2 ;
80changmax = 0.01 ;
81xmin = 0. ;
82xmax = 1. ;
83
84** Mesh
85p0 = 0. 0. ;
86p1 = 0. h ;
87ll = DROI nely p0 p1 ;
88mesh = ll TRAN nelx (l 0.) ;
89con = CONT mesh ;
90p2 = con POIN 'PROC' (l 0.) ;
91
92** Mechanical model
93mod = MODE mesh 'MECANIQUE' ;
94
95** Material properties field with a unit Young modulus
96maun = MATE mod 'YOUN' 1. 'NU' nu 'DIM3' e ;
97
98** Boundary conditions
99blo = (BLOQ 'UX' ll) ET (BLOQ 'UY' p2) ;
100
101** Load (local force)
102f = FORC (0. -1.) p1 ;
103
104** Volume of each element (ve) and total volume (vtot)
105un = MANU 'CHML' mod 'SCAL' 1. 'GRAVITE' ;
106ve = INTG 'ELEM' mod un maun ;
107vtot = INTG mod un maun ;
108
109** Initial density field
110xini = volfrac ;
111x = MANU 'CHML' mod 'SCAL' xini 'GRAVITE' ;
112change = 1. ;
113
114** Physical density field
115SI ((ft EGA 1) OU (ft EGA 2)) ;
116 xphys = x ;
117FINSI ;
118SI (ft EGA 3) ;
119 xtilde = x ;
120 xphys = 1. - (EXP (-1. * beta * xtilde)) + (xtilde * (EXP (-1. * beta))) ;
121FINSI ;
122
123** Initial volume fraction
124vfx = (INTG mod xphys maun) / vtot ;
125
126** Gravity centers and filtering matrix
127* the weight depends on the volume of each element
128ptg = un POIN 'SUPERIEUR' 0. ;
129wg = MANU 'CHPO' ptg 'SCAL' (EXTR ve 'VALE' 'SCAL') ;
130kfil = MFIL wg rmin 1. 0. ;
131ung = MANU 'CHPO' ptg 1 'SCAL' 1. ;
132ks = kfil * ung ;
133kfil = NFIL kfil ks ;
134
135** Lets's start a timer
136TEMP 'ZERO' ;
137
138** Topology optimization loop
139liso = PROG 0. 'PAS' 0.05 1. ;
140loopbeta = 0 ;
141liter = PROG ;
142lobj = PROG ;
143lvf = PROG ;
144lchange = PROG ;
145REPE b1 500 ;
146 loopbeta = loopbeta + 1 ;
147* penalization of the stiffness matrix (modified SIMP)
148 ep = emin + ((xphys ** penal) * (e0 - emin)) ;
149 map = MATE mod 'YOUN' ep 'NU' nu 'DIM3' e ;
150 k = RIGI mod map ;
151* resolution of the FE problem
152 kbc = k ET blo ;
153 u = RESO kbc f ;
154* objective function: compliance = u^T.F
155* value
156 c = MAXI (RESU (PSCA u f (MOTS 'UX' 'UY') (MOTS 'FX' 'FY'))) ;
157* sensitivity (gradient with respect to the physical variable)
158 eps = EPSI 'LINE' mod u ;
159 sigun = ELAS mod maun eps ;
160 eneun = ENER mod eps sigun ;
161 eneun = INTG 'ELEM' mod eneun ;
162 dc = -1. * penal * (xphys ** (penal - 1.)) * (e0 - emin) * eneun ;
163* constraint function: g = vf(x) - volfrac
164* value
165 g = vfx - volfrac ;
166* sensitivity (gradient with respect to the physical variable)
167 dg = ve / vtot ;
168* ft = 1 --> filtering the complicance sensitivity
169 SI (ft EGA 1) ;
170 dc = (HFILT (x * dc) mod kfil ptg) / (BORN x 'MINIMUM' 1.E-3) ;
171 FINSI ;
172* ft = 2,3 --> updating the sensitivies (to become gradients with respect to the design variable)
173 SI (ft EGA 2) ;
174 dc = HFILT dc mod kfil ptg ;
175 dg = HFILT dg mod kfil ptg ;
176 FINSI ;
177 SI (ft EGA 3) ;
178 dx = (beta * (EXP (-1. * beta * xtilde))) + (EXP (-1. * beta)) ;
179 dctilde = dc * dx ;
180 dc = HFILT dctilde mod kfil ptg ;
181 dgtilde = dg * dx ;
182 dg = HFILT dgtilde mod kfil ptg ;
183 FINSI ;
184* information about the current topology
185 info = CHAI 'It:' (&b1 - 1) / 5 'Obj:' / 10 c > 1 'Vol. frac:' > 4 vfx > 1 'Change:' > 4 change > 1 'Beta:' > 4 beta > 2 ;
186 MESS info ;
187 SI itrac ;
188 def1 = DEFO mesh u 0.05 ;
189 TRAC xphys mod con def1 liso 'TITR' info 'NCLK' ;
190 FINSI ;
191* update the topology by optimization
192 l1 = 0. ;
193 l2 = 1.E9 ;
194 REPE b2 200 ;
195 SI (((l2 - l1) / (l1 + l2)) < 0.001) ;
196 QUIT b2 ;
197 FINSI ;
198 lmid = 0.5 * (l1 + l2) ;
199 b = -1. * dc / (lmid * dg) ;
200 xinf = BORN (x - move) 'MINIMUM' xmin ;
201 xsup = BORN (x + move) 'MAXIMUM' xmax ;
202 xnew = x * (b ** 0.5) ;
203 minf = (xnew - xinf) MASQ 'INFERIEUR' 0. ;
204 mmil = ((xnew - xinf) MASQ 'SUPERIEUR' 0.) * ((xnew - xsup) MASQ 'INFERIEUR' 0.) ;
205 msup = (xnew - xsup) MASQ 'SUPERIEUR' 0. ;
206 xnew = (xinf * minf) + (xnew * mmil) + (xsup * msup) ;
207* update the physical density (by filtering and thresholding)
208 SI (ft EGA 1) ;
209 xphys = xnew ;
210 FINSI ;
211 SI (ft EGA 2) ;
212 xphys = HFILT xnew mod kfil ptg ;
213 FINSI ;
214 SI (ft EGA 3) ;
215 xtilde = HFILT xnew mod kfil ptg ;
216 xphys = 1. - (EXP (-1. * beta * xtilde)) + (xtilde * (EXP (-1. * beta))) ;
217 FINSI ;
218* updating the volume fraction and the constraint function
219 vfx = (INTG mod xphys maun) / vtot ;
220 g = vfx - volfrac ;
221 SI (g > 0.) ;
222 l1 = lmid ;
223 SINON ;
224 l2 = lmid ;
225 FINSI ;
226 FIN b2 ;
227* summary of the current iteration
228 change = MAXI 'ABS' (xnew - x) ;
229 liter = liter ET &b1 ;
230 lobj = lobj ET c ;
231 lvf = lvf ET vfx ;
232 lchange = lchange ET change ;
233* preparing the next iteration
234 x = xnew ;
235* ft = 3 --> updatind the slope of the thresholding function
236 SI (ft EGA 3) ;
237 SI ((beta < 512) ET ((loopbeta >EG 50) OU (change <EG changmax))) ;
238 beta = 2 * beta ;
239 loopbeta = 0 ;
240 change = 1. ;
241 FINSI ;
242 FINSI ;
243* stop criterion
244 SI (change < changmax) ;
245 info = CHAI 'It:' &b1 / 5 'Obj:' / 10 c > 1 'Vol. frac:' > 4 vfx > 1 'Change:' > 4 change > 1 'Beta:' > 4 beta > 2 ;
246 MESS info ;
247 QUIT b1 ;
248 FINSI ;
249FIN b1 ;
250
251fin;
252** Elapsed time
253TEMP 'IMPR' 'SOMM' 'HORL' ;
254
255** Plotting the final topology
256evobj = EVOL 'ROUG' 'MANU' 'Iterations' liter 'Compliance' lobj ;
257evvf = EVOL 'ORAN' 'MANU' 'Iterations' liter 'Frac. vol.' lvf ;
258evchange = EVOL 'VERT' 'MANU' 'Iterations' liter 'Max. change' lchange ;
259SI itrac ;
260 info = CHAI '[Final topology]' ' ' info ;
261 TRAC xphys mod con liso 'TITR' info ;
262 DESS evobj 'TITR' 'Objective function' ;
263 DESS evvf 'TITR' 'Volume fraction' ;
264 DESS evchange 'TITR' 'Max. change' ;
265FINSI ;
266
267** Mesh of the final topology
268tab1 = TABL ;
269tab1 . 'EPAISSEUR' = e ;
270tab1 . 'MODELE' = mod ;
271tab1 . 'TOPOLOGIE' = xphys ;
272meshf = (TOPOSURF tab1) COUL 'GRIS' ;
273edge = ARET meshf ;
274SI itrac ;
275 TRAC 'FACE' meshf 'ARET' edge 'TITR' 'Mesh of the final topology' ;
276FINSI ;
277
278FIN ;
top_mma.dgibi
Script Cast3M top_mma.dgibi
1************************************************************************
2* Example of topological optimization *
3* Density method + penalization (SIMP) *
4* *
5* This test case is a reproduction of the reference code: *
6* Andreassen, Clausen, Schevenels, Lazarov, & Sigmund (2011) *
7* "Efficient topology optimization in MATLAB using 88 lines of code" *
8* Structural and Multidisciplinary Optimization, 43(1), 1-16 *
9* https://doi.org/10.1007/s00158-010-0594-7 *
10* *
11* Application to a beam under bending in 2D plane stresses *
12* *
13* Optimization algorithm: Method of Moving Asymptotes *
14* *
15* Problem: minimize the compliance with a unilateral constraint *
16* on the volume fraction *
17* *
18* min C(x) = u^T.F *
19* s.t. G(x) = vf(x)/volfrac - 1 < 0 *
20* *
21* | *
22* | Force *
23* v *
24* O>+------------------------------------------------------+ *
25* | p1 | *
26* | | *
27* O>| | *
28* | | *
29* | | *
30* O>+------------------------------------------------------+ p2 *
31* Ʌ *
32* O *
33* *
34************************************************************************
35
36
37** Filtering procedure
38DEBP HFILT cham1*'MCHAML' mo*'MMODEL' mf*'RIGIDITE' mpt*'MAILLAGE' ;
39 chp1 = MANU 'CHPO' mpt 'SCAL' (EXTR cham1 'VALE' 'SCAL') ;
40 chp2 = mf * chp1 ;
41 cham2 = MANU 'CHML' mo 'REPA' 'SCAL' (EXTR chp2 'VALE' mpt) 'TYPE' 'SCALAIRE' 'GRAVITE' ;
42FINP cham2 ;
43
44** Global parameters
45itrac = VRAI ;
46OPTI 'DIME' 2 'MODE' 'PLAN' 'CONT' 'ELEM' 'QUA4' 'ECHO' ;
47
48** Geometrical parameters (width and height)
49l = 60. ;
50h = 20. ;
51e = 1. ;
52
53** Mesh parameters
54nelx = 60 ;
55nely = 20 ;
56
57** Material parameters (isotropic elasticity)
58e0 = 1. ;
59emin = e0 / 1.E9 ;
60nu = 0.3 ;
61
62** Topology optimization parameters
63* penal : SIMP penalization coefficient
64* volfrac : minimal volume fraction
65* ft : = 1 apply filter on the compliance sensitivity field
66* = 2 apply spatial filter on the density field
67* = 3 apply spatial filter on the density field + thresholding by heaviside function
68* rmin : filter radius (in m)
69* beta : initial value of the slope of the heaviside function (only for ft = 3)
70* will be doubled every 50 iterations
71* move : limit of the maximal increment of density
72* changmax : optimization stop criterion
73* xmin/xmax : min and max density bounds
74penal = 3. ;
75volfrac = 0.5 ;
76ft = 2 ;
77rmin = 0.04 * l ;
78beta = 1 ;
79move = 0.2 ;
80changmax = 0.01 ;
81xmin = 0. ;
82xmax = 1. ;
83
84** Mesh
85p0 = 0. 0. ;
86p1 = 0. h ;
87ll = DROI nely p0 p1 ;
88mesh = ll TRAN nelx (l 0.) ;
89con = CONT mesh ;
90p2 = con POIN 'PROC' (l 0.) ;
91
92** Mechanical model
93mod = MODE mesh 'MECANIQUE' ;
94
95** Material properties field with a unit Young modulus
96maun = MATE mod 'YOUN' 1. 'NU' nu 'DIM3' e ;
97
98** Boundary conditions
99blo = (BLOQ 'UX' ll) ET (BLOQ 'UY' p2) ;
100
101** Load (local force)
102f = FORC (0. -1.) p1 ;
103
104** Volume of each element (ve) and total volume (vtot)
105un = MANU 'CHML' mod 'SCAL' 1. 'GRAVITE' ;
106ve = INTG 'ELEM' mod un maun ;
107vtot = INTG mod un maun ;
108
109** Initial density field
110xini = volfrac ;
111x = MANU 'CHML' mod 'SCAL' xini 'GRAVITE' ;
112change = 1. ;
113
114** Physical density field
115SI ((ft EGA 1) OU (ft EGA 2)) ;
116 xphys = x ;
117FINSI ;
118SI (ft EGA 3) ;
119 xtilde = x ;
120 xphys = 1. - (EXP (-1. * beta * xtilde)) + (xtilde * (EXP (-1. * beta))) ;
121FINSI ;
122
123** Initial volume fraction
124vfx = (INTG mod xphys maun) / vtot ;
125
126** Gravity centers and filtering matrix
127* the weight depends on the volume of each element
128ptg = un POIN 'SUPERIEUR' 0. ;
129wg = MANU 'CHPO' ptg 'SCAL' (EXTR ve 'VALE' 'SCAL') ;
130kfil = MFIL wg rmin 1. 0. ;
131ung = MANU 'CHPO' ptg 1 'SCAL' 1. ;
132ks = kfil * ung ;
133kfil = NFIL kfil ks ;
134
135** Initialization of the table for the MMA
136nx = NBEL mesh ;
137tmma = TABL ;
138* initial values of the design variables x
139lx = EXTR x 'VALE' 'SCAL' ;
140tmma . 'X' = lx ;
141* bounds for x
142tmma . 'XMIN' = xmin ;
143tmma . 'XMAX' = xmax ;
144* asymptotes
145tmma . 'LOW' = PROG nx*1. ;
146tmma . 'UPP' = PROG nx*1. ;
147* other parameters
148tmma . 'A0' = 1. ;
149tmma . 'A' = PROG 0. ;
150tmma . 'C' = PROG 1.E4 ;
151tmma . 'D' = PROG 0. ;
152tmma . 'MOVE' = move ;
153
154** Lets's start a timer
155TEMP 'ZERO' ;
156
157** Topology optimization loop
158liso = PROG 0. 'PAS' 0.05 1. ;
159loopbeta = 0 ;
160liter = PROG ;
161lobj = PROG ;
162lvf = PROG ;
163lchange = PROG ;
164REPE b1 500 ;
165 loopbeta = loopbeta + 1 ;
166* penalization of the stiffness matrix (modified SIMP)
167 ep = emin + ((xphys ** penal) * (e0 - emin)) ;
168 map = MATE mod 'YOUN' ep 'NU' nu 'DIM3' e ;
169 k = RIGI mod map ;
170* resolution of the FE problem
171 kbc = k ET blo ;
172 u = RESO kbc f ;
173* objective function: compliance = u^T.K.u
174* value
175 c = MAXI (RESU (PSCA u f (MOTS 'UX' 'UY') (MOTS 'FX' 'FY'))) ;
176* sensitivity (gradient with respect to the physical variable)
177 eps = EPSI 'LINE' mod u ;
178 sigun = ELAS mod maun eps ;
179 eneun = ENER mod eps sigun ;
180 eneun = INTG 'ELEM' mod eneun ;
181 dc = -1. * penal * (xphys ** (penal - 1.)) * (e0 - emin) * eneun ;
182* constraint function: g = vf(x)/volfrac - 1
183* value
184 g = (vfx / volfrac) - 1. ;
185* sensitivity (gradient with respect to the physical variable)
186 dg = ve / vtot / volfrac ;
187* ft = 1 --> filtering the complicance sensitivity
188 SI (ft EGA 1) ;
189 dc = (HFILT (x * dc) mod kfil ptg) / (BORN x 'MINIMUM' 1.E-3) ;
190 FINSI ;
191* ft = 2,3 --> updating the sensitivies (to become gradients with respect to the design variable)
192 SI (ft EGA 2) ;
193 dc = HFILT dc mod kfil ptg ;
194 dg = HFILT dg mod kfil ptg ;
195 FINSI ;
196 SI (ft EGA 3) ;
197 dx = (beta * (EXP (-1. * beta * xtilde))) + (EXP (-1. * beta)) ;
198 dctilde = dc * dx ;
199 dc = HFILT dctilde mod kfil ptg ;
200 dgtilde = dg * dx ;
201 dg = HFILT dgtilde mod kfil ptg ;
202 FINSI ;
203* information about the current topology
204 info = CHAI 'It:' (&b1 - 1) / 5 'Obj:' / 10 c > 1 'Vol. frac:' > 4 vfx > 1 'Change:' > 4 change > 1 'Beta:' > 4 beta > 2 ;
205 MESS info ;
206 SI itrac ;
207 def1 = DEFO mesh u 0.05 ;
208 TRAC xphys mod con def1 liso 'TITR' info 'NCLK' ;
209 FINSI ;
210* update the topology by optimization
211 tmma . 'F0VAL' = c ;
212 tmma . 'DF0DX' = EXTR dc 'VALE' 'SCAL' ;
213 tmma . 'FVAL' = PROG g ;
214 tmma . 'DFDX' = TABL ;
215 tmma . 'DFDX' . 1 = EXTR dg 'VALE' 'SCAL' ;
216 MMA tmma ;
217 lxnew = tmma . 'X' ;
218 xnew = MANU 'CHML' mod 'REPA' 'SCAL' lxnew 'TYPE' 'SCALAIRE' 'GRAVITE' ;
219* update the physical density (by filtering and thresholding)
220 SI (ft EGA 1) ;
221 xphys = xnew ;
222 FINSI ;
223 SI (ft EGA 2) ;
224 xphys = HFILT xnew mod kfil ptg ;
225 FINSI ;
226 SI (ft EGA 3) ;
227 xtilde = HFILT xnew mod kfil ptg ;
228 xphys = 1. - (EXP (-1. * beta * xtilde)) + (xtilde * (EXP (-1. * beta))) ;
229 FINSI ;
230* updating the volume fraction
231 vfx = (INTG mod xphys maun) / vtot ;
232* summary of the current iteration
233 change = MINI (MAXI 'ABS' (lxnew - (tmma . 'XOLD1')))
234 (MAXI 'ABS' (lxnew - (tmma . 'XOLD2'))) ;
235 liter = liter ET &b1 ;
236 lobj = lobj ET c ;
237 lvf = lvf ET vfx ;
238 lchange = lchange ET change ;
239* preparing the next iteration
240 x = xnew ;
241 lx = lxnew ;
242* ft = 3 --> updatind the slope of the thresholding function
243 SI (ft EGA 3) ;
244 SI ((beta < 512) ET ((loopbeta >EG 50) OU (change <EG changmax))) ;
245 beta = 2 * beta ;
246 loopbeta = 0 ;
247 change = 1. ;
248 FINSI ;
249 FINSI ;
250* stop criterion
251 SI (change < changmax) ;
252 info = CHAI 'It:' &b1 / 5 'Obj:' / 10 c > 1 'Vol. frac:' > 4 vfx > 1 'Change:' > 4 change > 1 'Beta:' > 4 beta > 2 ;
253 MESS info ;
254 QUIT b1 ;
255 FINSI ;
256FIN b1 ;
257
258fin;
259** Elapsed time
260TEMP 'IMPR' 'SOMM' 'HORL' ;
261
262** Plotting the final topology
263evobj = EVOL 'ROUG' 'MANU' 'Iterations' liter 'Compliance' lobj ;
264evvf = EVOL 'ORAN' 'MANU' 'Iterations' liter 'Frac. vol.' lvf ;
265evchange = EVOL 'VERT' 'MANU' 'Iterations' liter 'Max. change' lchange ;
266SI itrac ;
267 info = CHAI '[Final topology]' ' ' info ;
268 TRAC xphys mod con liso 'TITR' info ;
269 DESS evobj 'TITR' 'Objective function' ;
270 DESS evvf 'TITR' 'Volume fraction' ;
271 DESS evchange 'TITR' 'Max. change' ;
272FINSI ;
273
274** Mesh of the final topology
275tab1 = TABL ;
276tab1 . 'EPAISSEUR' = e ;
277tab1 . 'MODELE' = mod ;
278tab1 . 'TOPOLOGIE' = xphys ;
279meshf = (TOPOSURF tab1) COUL 'GRIS' ;
280edge = ARET meshf ;
281SI itrac ;
282 TRAC 'FACE' meshf 'ARET' edge 'TITR' 'Mesh of the final topology' ;
283FINSI ;
284
285FIN ;