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.53 s
Generate movie 02/44 (4.55 %) 3.68 s
Generate movie 03/44 (6.82 %) 3.70 s
Generate movie 04/44 (9.09 %) 3.52 s
Generate movie 05/44 (11.36 %) 3.42 s
Generate movie 06/44 (13.64 %) 3.29 s
Generate movie 07/44 (15.91 %) 3.33 s
Generate movie 08/44 (18.18 %) 3.17 s
Generate movie 09/44 (20.45 %) 3.14 s
Generate movie 10/44 (22.73 %) 3.22 s
Generate movie 11/44 (25.00 %) 2.88 s
Generate movie 12/44 (27.27 %) 2.87 s
Generate movie 13/44 (29.55 %) 2.79 s
Generate movie 14/44 (31.82 %) 2.67 s
Generate movie 15/44 (34.09 %) 2.59 s
Generate movie 16/44 (36.36 %) 2.48 s
Generate movie 17/44 (38.64 %) 2.46 s
Generate movie 18/44 (40.91 %) 2.40 s
Generate movie 19/44 (43.18 %) 2.28 s
Generate movie 20/44 (45.45 %) 2.15 s
Generate movie 21/44 (47.73 %) 2.08 s
Generate movie 22/44 (50.00 %) 2.07 s
Generate movie 23/44 (52.27 %) 1.87 s
Generate movie 24/44 (54.55 %) 1.80 s
Generate movie 25/44 (56.82 %) 1.68 s
Generate movie 26/44 (59.09 %) 1.69 s
Generate movie 27/44 (61.36 %) 1.57 s
Generate movie 28/44 (63.64 %) 1.48 s
Generate movie 29/44 (65.91 %) 1.38 s
Generate movie 30/44 (68.18 %) 1.30 s
Generate movie 31/44 (70.45 %) 1.16 s
Generate movie 32/44 (72.73 %) 1.09 s
Generate movie 33/44 (75.00 %) 991.54 ms
Generate movie 34/44 (77.27 %) 893.75 ms
Generate movie 35/44 (79.55 %) 833.73 ms
Generate movie 36/44 (81.82 %) 742.00 ms
Generate movie 37/44 (84.09 %) 632.92 ms
Generate movie 38/44 (86.36 %) 532.19 ms
Generate movie 39/44 (88.64 %) 459.04 ms
Generate movie 40/44 (90.91 %) 372.91 ms
Generate movie 41/44 (93.18 %) 265.92 ms
Generate movie 42/44 (95.45 %) 179.90 ms
Generate movie 43/44 (97.73 %) 89.12 ms
Generate movie 44/44 (100.00 %) 0.00 µs
13 from EasyFEA import Display, Folder, PyVista, np, ElemType, Models, Simulations
14 from EasyFEA.fem import FeArray, Field, BiLinearForm, Sym_Grad, Trace
15 from EasyFEA.Geoms import Domain
16
17 if __name__ == "__main__":
18
19 Display.Clear()
20
21 # ----------------------------------------------
22 # Configuration
23 # ----------------------------------------------
24
25 dim = 2
26
27 L, H = 60, 30
28 # L, H = 120, 60
29
30 # optim topo
31 iterMax = 60
32 volFrac = 0.4
33 penal = 3
34 rMin = 3
35
36 # outputs
37 generateMovie = True
38 folder = Folder.Join(
39 Folder.RESULTS_DIR,
40 "WeakForms",
41 "TopologyOptimisation1",
42 )
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.ElasIsot(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.WeakFormSimu(mesh, model)
122 simu._Solver_Set_PETSc4Py_Options("none", "lu")
123
124 simu.add_dirichlet(nodesX0, [0] * dim, simu.Get_unknowns())
125 simu.add_neumann(nodesLoad, [-1], ["y"])
126
127 # ----------------------------------------------
128 # Optim topo
129 # ----------------------------------------------
130
131 err = 1.0
132 list_compliance: list[float] = []
133 list_p_e: list[np.ndarray] = []
134 iter = 0
135
136 while err > 0.005 and iter < iterMax:
137 iter += 1
138 pOld_e = p_e.copy()
139
140 # solve u
141 simu.Need_Update()
142 u = simu.Solve()
143 simu.Save_Iter()
144
145 # compute compliance for elements
146 u_e = field.groupElem.Locates_sol_e(u, dim)
147 K_e = ComputeK.Integrate_e(field)
148 uKu_e = np.einsum("ei,eij,ej->e", u_e, K_e, u_e, optimize="optimal")
149 c = (p_e**penal * uKu_e).sum()
150
151 # compute sensitivity for elements
152 dCdP_e = -(penal * (p_e ** (penal - 1.0)) * uKu_e)
153
154 # use sensitivity filter
155 dCdP_e = np.einsum("ij,j,j", Hij, p_e, dCdP_e) / (
156 np.einsum("i,ij", p_e, Hij) + 1e-12
157 )
158
159 # OC update (enforce volume)
160 lmin, lmax = 0.0, 1e5
161 pmin, pmax = 0.001, 1.0
162 move = 0.2
163 pNew_e = p_e.copy()
164
165 while (lmax - lmin) > 1e-4 * (lmax + lmin + 1e-16):
166 lmid = 0.5 * (lmax + lmin)
167 candidate = p_e * np.sqrt(np.maximum(-dCdP_e / lmid, 1e-9))
168 # Apply move limits and physical bounds [pmin, pmax]
169 pNew_e = np.maximum(
170 pmin,
171 np.maximum(
172 p_e - move,
173 np.minimum(pmax, np.minimum(p_e + move, candidate)),
174 ),
175 )
176 # update lambda to fit volume fraction
177 if pNew_e.sum() - volFrac * mesh.Ne > 0:
178 lmin = lmid
179 else:
180 lmax = lmid
181
182 # get updated density and compliance
183 p_e = pNew_e
184 list_p_e.append(p_e)
185 list_compliance.append(c)
186
187 # compute relative error : || p_e - pOld_e || / || pOld_e ||
188 err = np.linalg.norm(p_e - pOld_e) / np.linalg.norm(pOld_e)
189
190 Display.MyPrint(
191 f"Iteration {str(iter).zfill(len(str(iterMax)))}, compliance = {c:.3e}, volume fraction = {p_e.mean():.3f}, err = {err:.3e}",
192 end="\r",
193 )
194
195 # ----------------------------------------------
196 # Results
197 # ----------------------------------------------
198
199 axC = Display.Init_Axes()
200 axC.plot(range(len(list_compliance)), list_compliance, ls="-", marker=".")
201 axC.set_xlabel("Iteration")
202 axC.set_ylabel("Compliance")
203 Display.plt.show()
204
205 PyVista.Plot_BoundaryConditions(simu).show()
206
207 def get_thresh(p_e: np.ndarray, min=0.5, max=1.0):
208 grid = PyVista._pvMesh(mesh, p_e, nodeValues=False)
209 for result in simu.Results_Available():
210 grid[result] = simu.Result(result).reshape(mesh.Nn, -1)
211 thresh = grid.threshold((min, max))
212 return thresh
213
214 thresh = get_thresh(p_e)
215 PyVista.Plot(thresh, color="k").show()
216 PyVista.Plot(thresh, "uy").show()
217
218 if generateMovie:
219
220 def Func(plotter: PyVista.pv.Plotter, iter):
221 simu.Set_Iter(iter)
222 thresh = get_thresh(list_p_e[iter])
223 plotter.add_title(
224 f"{str(iter+1).zfill(len(str(iterMax)))}/{len(list_compliance)}"
225 )
226 PyVista.Plot(thresh, color="k", plotter=plotter)
227
228 PyVista.Movie_func(Func, len(list_compliance), folder, "optim.gif")
Total running time of the script: (0 minutes 10.710 seconds)