Note
Go to the end to download the full example code.
TopologyOptimisation1#
An educational implementation of topology optimization inspired by Week 10- Topology Optimisation — A Step-by-Step Tutorial created by (Dr Wei Tan, Queen Mary University of London), which in turn builds upon the seminal 88-line topology optimization MATLAB code by Ole Sigmund (2001), published in Structural and Multidisciplinary Optimization, 21(2), pp. 120–127.
Iteration 01, compliance = 6.865e-01, volume fraction = 0.400, err = 3.530e-01
Iteration 02, compliance = 4.338e-01, volume fraction = 0.400, err = 2.860e-01
Iteration 03, compliance = 3.075e-01, volume fraction = 0.400, err = 2.197e-01
Iteration 04, compliance = 2.481e-01, volume fraction = 0.400, err = 1.499e-01
Iteration 05, compliance = 2.252e-01, volume fraction = 0.400, err = 1.268e-01
Iteration 06, compliance = 2.048e-01, volume fraction = 0.400, err = 1.131e-01
Iteration 07, compliance = 1.911e-01, volume fraction = 0.400, err = 1.029e-01
Iteration 08, compliance = 1.762e-01, volume fraction = 0.400, err = 9.478e-02
Iteration 09, compliance = 1.655e-01, volume fraction = 0.400, err = 9.155e-02
Iteration 10, compliance = 1.530e-01, volume fraction = 0.400, err = 8.394e-02
Iteration 11, compliance = 1.424e-01, volume fraction = 0.400, err = 8.602e-02
Iteration 12, compliance = 1.299e-01, volume fraction = 0.400, err = 7.927e-02
Iteration 13, compliance = 1.190e-01, volume fraction = 0.400, err = 7.109e-02
Iteration 14, compliance = 1.102e-01, volume fraction = 0.400, err = 5.719e-02
Iteration 15, compliance = 1.050e-01, volume fraction = 0.400, err = 4.511e-02
Iteration 16, compliance = 1.020e-01, volume fraction = 0.400, err = 3.708e-02
Iteration 17, compliance = 1.001e-01, volume fraction = 0.400, err = 3.229e-02
Iteration 18, compliance = 9.868e-02, volume fraction = 0.400, err = 2.715e-02
Iteration 19, compliance = 9.765e-02, volume fraction = 0.400, err = 2.329e-02
Iteration 20, compliance = 9.682e-02, volume fraction = 0.400, err = 2.167e-02
Iteration 21, compliance = 9.612e-02, volume fraction = 0.400, err = 2.040e-02
Iteration 22, compliance = 9.551e-02, volume fraction = 0.400, err = 1.928e-02
Iteration 23, compliance = 9.497e-02, volume fraction = 0.400, err = 1.557e-02
Iteration 24, compliance = 9.447e-02, volume fraction = 0.400, err = 1.441e-02
Iteration 25, compliance = 9.403e-02, volume fraction = 0.400, err = 1.327e-02
Iteration 26, compliance = 9.364e-02, volume fraction = 0.400, err = 1.223e-02
Iteration 27, compliance = 9.330e-02, volume fraction = 0.400, err = 1.156e-02
Iteration 28, compliance = 9.300e-02, volume fraction = 0.400, err = 1.094e-02
Iteration 29, compliance = 9.272e-02, volume fraction = 0.400, err = 1.042e-02
Iteration 30, compliance = 9.247e-02, volume fraction = 0.400, err = 9.755e-03
Iteration 31, compliance = 9.224e-02, volume fraction = 0.400, err = 7.533e-03
Iteration 32, compliance = 9.202e-02, volume fraction = 0.400, err = 6.743e-03
Iteration 33, compliance = 9.180e-02, volume fraction = 0.400, err = 6.453e-03
Iteration 34, compliance = 9.161e-02, volume fraction = 0.400, err = 6.228e-03
Iteration 35, compliance = 9.143e-02, volume fraction = 0.400, err = 6.105e-03
Iteration 36, compliance = 9.126e-02, volume fraction = 0.400, err = 5.983e-03
Iteration 37, compliance = 9.110e-02, volume fraction = 0.400, err = 5.893e-03
Iteration 38, compliance = 9.095e-02, volume fraction = 0.400, err = 5.830e-03
Iteration 39, compliance = 9.081e-02, volume fraction = 0.400, err = 5.795e-03
Iteration 40, compliance = 9.067e-02, volume fraction = 0.400, err = 5.681e-03
Iteration 41, compliance = 9.055e-02, volume fraction = 0.400, err = 5.482e-03
Iteration 42, compliance = 9.042e-02, volume fraction = 0.400, err = 5.380e-03
Iteration 43, compliance = 9.031e-02, volume fraction = 0.400, err = 5.259e-03
Iteration 44, compliance = 9.020e-02, volume fraction = 0.400, err = 4.837e-03
Generate movie 01/44 (2.27 %) 4.85 s
Generate movie 02/44 (4.55 %) 4.00 s
Generate movie 03/44 (6.82 %) 3.80 s
Generate movie 04/44 (9.09 %) 3.73 s
Generate movie 05/44 (11.36 %) 3.71 s
Generate movie 06/44 (13.64 %) 3.66 s
Generate movie 07/44 (15.91 %) 3.52 s
Generate movie 08/44 (18.18 %) 3.45 s
Generate movie 09/44 (20.45 %) 3.42 s
Generate movie 10/44 (22.73 %) 3.35 s
Generate movie 11/44 (25.00 %) 3.18 s
Generate movie 12/44 (27.27 %) 3.06 s
Generate movie 13/44 (29.55 %) 2.97 s
Generate movie 14/44 (31.82 %) 2.86 s
Generate movie 15/44 (34.09 %) 2.75 s
Generate movie 16/44 (36.36 %) 2.65 s
Generate movie 17/44 (38.64 %) 2.67 s
Generate movie 18/44 (40.91 %) 2.53 s
Generate movie 19/44 (43.18 %) 2.41 s
Generate movie 20/44 (45.45 %) 2.28 s
Generate movie 21/44 (47.73 %) 2.22 s
Generate movie 22/44 (50.00 %) 2.16 s
Generate movie 23/44 (52.27 %) 2.05 s
Generate movie 24/44 (54.55 %) 1.95 s
Generate movie 25/44 (56.82 %) 1.83 s
Generate movie 26/44 (59.09 %) 1.75 s
Generate movie 27/44 (61.36 %) 1.64 s
Generate movie 28/44 (63.64 %) 1.54 s
Generate movie 29/44 (65.91 %) 1.44 s
Generate movie 30/44 (68.18 %) 1.35 s
Generate movie 31/44 (70.45 %) 1.26 s
Generate movie 32/44 (72.73 %) 1.14 s
Generate movie 33/44 (75.00 %) 1.11 s
Generate movie 34/44 (77.27 %) 960.90 ms
Generate movie 35/44 (79.55 %) 867.41 ms
Generate movie 36/44 (81.82 %) 800.11 ms
Generate movie 37/44 (84.09 %) 684.50 ms
Generate movie 38/44 (86.36 %) 577.27 ms
Generate movie 39/44 (88.64 %) 481.04 ms
Generate movie 40/44 (90.91 %) 382.94 ms
Generate movie 41/44 (93.18 %) 286.94 ms
Generate movie 42/44 (95.45 %) 192.91 ms
Generate movie 43/44 (97.73 %) 99.35 ms
Generate movie 44/44 (100.00 %) 0.00 µs
14 import matplotlib.pyplot as plt
15 import numpy as np
16
17 from EasyFEA import Display, Folder, PyVista, ElemType, Models, Simulations
18 from EasyFEA.FEM import FeArray, Field, BiLinearForm, Sym_Grad, Trace
19 from EasyFEA.Geoms import Domain
20
21 if __name__ == "__main__":
22
23 Display.Clear()
24
25 # ----------------------------------------------
26 # Configuration
27 # ----------------------------------------------
28
29 dim = 2
30
31 L, H = 60, 30
32 # L, H = 120, 60
33
34 # optim topo
35 iterMax = 60
36 volFrac = 0.4
37 penal = 3
38 rMin = 3
39
40 # outputs
41 generateMovie = True
42 folder = Folder.Results_Dir()
43
44 # ----------------------------------------------
45 # Mesh
46 # ----------------------------------------------
47
48 meshSize = 1 if dim == 2 else H / 10
49 contour = Domain((0, 0), (L, H), meshSize)
50 assert H / meshSize % 2 == 0
51
52 if dim == 2:
53 mesh = contour.Mesh_2D([], ElemType.QUAD4, isOrganised=True)
54 else:
55 mesh = contour.Mesh_Extrude(
56 [], [0, 0, H], [H / meshSize], ElemType.HEXA8, isOrganised=True
57 )
58
59 nodesX0 = mesh.Nodes_Conditions(lambda x, y, z: x == 0)
60
61 zMean = 0 if dim == 2 else H / 2
62 nodesLoad = mesh.Nodes_Point((L, H / 2, zMean))
63 # nodesLoad = mesh.Nodes_Conditions(lambda x, y, z: y == H)
64
65 # ----------------------------------------------
66 # Mesh-Independence Sensitivity Filter (Sigmund, 1998)
67 # ----------------------------------------------
68
69 # get the coordinates of each elements
70 coord_e = mesh.coord[mesh.connect].mean(1)
71
72 # compute Hij
73 elements = mesh.groupElem.elements
74 Hij = np.array(
75 [
76 np.maximum(0, rMin - np.linalg.norm(coord_e[i] - coord_e, axis=-1) + 1e-12)
77 for i in range(mesh.Ne)
78 ]
79 )
80
81 # # plot neighbor elements
82 # elem = 100
83 # ax = Display.Plot_Mesh(mesh, alpha=0)
84 # Display.Plot_Elements(
85 # mesh, mesh.connect[Hij[elem] != 0].ravel(), 2, color="blue", ax=ax
86 # )
87 # Display.Plot_Elements(mesh, mesh.connect[elem].ravel(), 2, ax=ax)
88
89 # ----------------------------------------------
90 # Formulations
91 # ----------------------------------------------
92
93 elastic = Models.Elastic.Isotropic(dim, E=1, v=0.3, planeStress=True)
94 mu = elastic.get_mu()
95 lmbda = elastic.get_lambda()
96
97 def S(u: Field) -> FeArray:
98 Eps = Sym_Grad(u)
99 return 2 * mu * Eps + lmbda * Trace(Eps) * np.eye(dim)
100
101 p_e = np.ones(mesh.Ne, dtype=float) * volFrac
102
103 @BiLinearForm
104 def ComputeK(u: Field, v: Field):
105 Sig = S(u)
106 Eps = Sym_Grad(v)
107 return Sig.ddot(Eps)
108
109 @BiLinearForm
110 def ComputePenalizedK(u: Field, v: Field):
111 simpScaling = FeArray.asfearray(np.reshape(p_e**penal, (-1, 1)))
112 return simpScaling * ComputeK(u, v)
113
114 field = Field(mesh.groupElem, dim)
115 model = Models.WeakForms(field, ComputePenalizedK, thickness=H)
116
117 # ----------------------------------------------
118 # Simulation
119 # ----------------------------------------------
120
121 simu = Simulations.WeakForms(mesh, model)
122
123 simu.add_dirichlet(nodesX0, [0] * dim, simu.Get_unknowns())
124 simu.add_neumann(nodesLoad, [-1], ["y"])
125
126 # ----------------------------------------------
127 # Optim topo
128 # ----------------------------------------------
129
130 err = 1.0
131 list_compliance: list[float] = []
132 list_p_e: list[np.ndarray] = []
133 iter = 0
134
135 while err > 0.005 and iter < iterMax:
136 iter += 1
137 pOld_e = p_e.copy()
138
139 # solve u
140 simu.Need_Update()
141 u = simu.Solve()
142 simu.Save_Iter()
143
144 # compute compliance for elements
145 u_e = field.groupElem.Locates_sol_e(u, dim)
146 K_e = ComputeK.Integrate_e(field)
147 uKu_e = np.einsum("ei,eij,ej->e", u_e, K_e, u_e, optimize="optimal")
148 c = (p_e**penal * uKu_e).sum()
149
150 # compute sensitivity for elements
151 dCdP_e = -(penal * (p_e ** (penal - 1.0)) * uKu_e)
152
153 # use sensitivity filter
154 dCdP_e = np.einsum("ij,j,j", Hij, p_e, dCdP_e) / (
155 np.einsum("i,ij", p_e, Hij) + 1e-12
156 )
157
158 # OC update (enforce volume)
159 lmin, lmax = 0.0, 1e5
160 pmin, pmax = 0.001, 1.0
161 move = 0.2
162 pNew_e = p_e.copy()
163
164 while (lmax - lmin) > 1e-4 * (lmax + lmin + 1e-16):
165 lmid = 0.5 * (lmax + lmin)
166 candidate = p_e * np.sqrt(np.maximum(-dCdP_e / lmid, 1e-9))
167 # Apply move limits and physical bounds [pmin, pmax]
168 pNew_e = np.maximum(
169 pmin,
170 np.maximum(
171 p_e - move,
172 np.minimum(pmax, np.minimum(p_e + move, candidate)),
173 ),
174 )
175 # update lambda to fit volume fraction
176 if pNew_e.sum() - volFrac * mesh.Ne > 0:
177 lmin = lmid
178 else:
179 lmax = lmid
180
181 # get updated density and compliance
182 p_e = pNew_e
183 list_p_e.append(p_e)
184 list_compliance.append(c)
185
186 # compute relative error : || p_e - pOld_e || / || pOld_e ||
187 err = np.linalg.norm(p_e - pOld_e) / np.linalg.norm(pOld_e)
188
189 Display.MyPrint(
190 f"Iteration {str(iter).zfill(len(str(iterMax)))}, compliance = {c:.3e}, volume fraction = {p_e.mean():.3f}, err = {err:.3e}",
191 end="\r",
192 )
193
194 # ----------------------------------------------
195 # Results
196 # ----------------------------------------------
197
198 axC = Display.Init_Axes()
199 axC.plot(range(len(list_compliance)), list_compliance, ls="-", marker=".")
200 axC.set_xlabel("Iteration")
201 axC.set_ylabel("Compliance")
202 plt.show()
203
204 PyVista.Plot_BoundaryConditions(simu).show()
205
206 def get_thresh(p_e: np.ndarray, min=0.5, max=1.0):
207 grid = PyVista._pvMesh(mesh, p_e, nodeValues=False)
208 for result in simu.Results_Available():
209 grid[result] = simu.Result(result).reshape(mesh.Nn, -1)
210 thresh = grid.threshold((min, max))
211 return thresh
212
213 thresh = get_thresh(p_e)
214 PyVista.Plot(thresh, color="k").show()
215 PyVista.Plot(thresh, "uy").show()
216
217 if generateMovie:
218
219 def Func(plotter: PyVista.pv.Plotter, iter):
220 simu.Set_Iter(iter)
221 thresh = get_thresh(list_p_e[iter])
222 plotter.add_title(
223 f"{str(iter+1).zfill(len(str(iterMax)))}/{len(list_compliance)}"
224 )
225 PyVista.Plot(thresh, color="k", plotter=plotter)
226
227 PyVista.Movie_func(Func, len(list_compliance), folder, "optim.gif")
Total running time of the script: (0 minutes 11.795 seconds)