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.73 s
Generate movie 02/44 (4.55 %) 3.85 s
Generate movie 03/44 (6.82 %) 3.60 s
Generate movie 04/44 (9.09 %) 3.63 s
Generate movie 05/44 (11.36 %) 3.40 s
Generate movie 06/44 (13.64 %) 3.35 s
Generate movie 07/44 (15.91 %) 3.28 s
Generate movie 08/44 (18.18 %) 3.19 s
Generate movie 09/44 (20.45 %) 3.08 s
Generate movie 10/44 (22.73 %) 2.99 s
Generate movie 11/44 (25.00 %) 2.92 s
Generate movie 12/44 (27.27 %) 2.93 s
Generate movie 13/44 (29.55 %) 2.81 s
Generate movie 14/44 (31.82 %) 2.64 s
Generate movie 15/44 (34.09 %) 2.56 s
Generate movie 16/44 (36.36 %) 2.44 s
Generate movie 17/44 (38.64 %) 2.44 s
Generate movie 18/44 (40.91 %) 2.36 s
Generate movie 19/44 (43.18 %) 2.28 s
Generate movie 20/44 (45.45 %) 2.13 s
Generate movie 21/44 (47.73 %) 2.05 s
Generate movie 22/44 (50.00 %) 1.96 s
Generate movie 23/44 (52.27 %) 1.88 s
Generate movie 24/44 (54.55 %) 1.78 s
Generate movie 25/44 (56.82 %) 1.68 s
Generate movie 26/44 (59.09 %) 1.61 s
Generate movie 27/44 (61.36 %) 1.50 s
Generate movie 28/44 (63.64 %) 1.43 s
Generate movie 29/44 (65.91 %) 1.37 s
Generate movie 30/44 (68.18 %) 1.28 s
Generate movie 31/44 (70.45 %) 1.21 s
Generate movie 32/44 (72.73 %) 1.16 s
Generate movie 33/44 (75.00 %) 1.01 s
Generate movie 34/44 (77.27 %) 907.84 ms
Generate movie 35/44 (79.55 %) 821.38 ms
Generate movie 36/44 (81.82 %) 744.87 ms
Generate movie 37/44 (84.09 %) 641.37 ms
Generate movie 38/44 (86.36 %) 545.62 ms
Generate movie 39/44 (88.64 %) 468.38 ms
Generate movie 40/44 (90.91 %) 376.22 ms
Generate movie 41/44 (93.18 %) 280.46 ms
Generate movie 42/44 (95.45 %) 188.30 ms
Generate movie 43/44 (97.73 %) 90.52 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.Results_Dir()
39
40 # ----------------------------------------------
41 # Mesh
42 # ----------------------------------------------
43
44 meshSize = 1 if dim == 2 else H / 10
45 contour = Domain((0, 0), (L, H), meshSize)
46 assert H / meshSize % 2 == 0
47
48 if dim == 2:
49 mesh = contour.Mesh_2D([], ElemType.QUAD4, isOrganised=True)
50 else:
51 mesh = contour.Mesh_Extrude(
52 [], [0, 0, H], [H / meshSize], ElemType.HEXA8, isOrganised=True
53 )
54
55 nodesX0 = mesh.Nodes_Conditions(lambda x, y, z: x == 0)
56
57 zMean = 0 if dim == 2 else H / 2
58 nodesLoad = mesh.Nodes_Point((L, H / 2, zMean))
59 # nodesLoad = mesh.Nodes_Conditions(lambda x, y, z: y == H)
60
61 # ----------------------------------------------
62 # Mesh-Independence Sensitivity Filter (Sigmund, 1998)
63 # ----------------------------------------------
64
65 # get the coordinates of each elements
66 coord_e = mesh.coord[mesh.connect].mean(1)
67
68 # compute Hij
69 elements = mesh.groupElem.elements
70 Hij = np.array(
71 [
72 np.maximum(0, rMin - np.linalg.norm(coord_e[i] - coord_e, axis=-1) + 1e-12)
73 for i in range(mesh.Ne)
74 ]
75 )
76
77 # # plot neighbor elements
78 # elem = 100
79 # ax = Display.Plot_Mesh(mesh, alpha=0)
80 # Display.Plot_Elements(
81 # mesh, mesh.connect[Hij[elem] != 0].ravel(), 2, color="blue", ax=ax
82 # )
83 # Display.Plot_Elements(mesh, mesh.connect[elem].ravel(), 2, ax=ax)
84
85 # ----------------------------------------------
86 # Formulations
87 # ----------------------------------------------
88
89 elastic = Models.Elastic.Isotropic(dim, E=1, v=0.3, planeStress=True)
90 mu = elastic.get_mu()
91 lmbda = elastic.get_lambda()
92
93 def S(u: Field) -> FeArray:
94 Eps = Sym_Grad(u)
95 return 2 * mu * Eps + lmbda * Trace(Eps) * np.eye(dim)
96
97 p_e = np.ones(mesh.Ne, dtype=float) * volFrac
98
99 @BiLinearForm
100 def ComputeK(u: Field, v: Field):
101 Sig = S(u)
102 Eps = Sym_Grad(v)
103 return Sig.ddot(Eps)
104
105 @BiLinearForm
106 def ComputePenalizedK(u: Field, v: Field):
107 simpScaling = FeArray.asfearray(np.reshape(p_e**penal, (-1, 1)))
108 return simpScaling * ComputeK(u, v)
109
110 field = Field(mesh.groupElem, dim)
111 model = Models.WeakForms(field, ComputePenalizedK, thickness=H)
112
113 # ----------------------------------------------
114 # Simulation
115 # ----------------------------------------------
116
117 simu = Simulations.WeakForm(mesh, model)
118 simu._Solver_Set_PETSc4Py_Options("none", "lu")
119
120 simu.add_dirichlet(nodesX0, [0] * dim, simu.Get_unknowns())
121 simu.add_neumann(nodesLoad, [-1], ["y"])
122
123 # ----------------------------------------------
124 # Optim topo
125 # ----------------------------------------------
126
127 err = 1.0
128 list_compliance: list[float] = []
129 list_p_e: list[np.ndarray] = []
130 iter = 0
131
132 while err > 0.005 and iter < iterMax:
133 iter += 1
134 pOld_e = p_e.copy()
135
136 # solve u
137 simu.Need_Update()
138 u = simu.Solve()
139 simu.Save_Iter()
140
141 # compute compliance for elements
142 u_e = field.groupElem.Locates_sol_e(u, dim)
143 K_e = ComputeK.Integrate_e(field)
144 uKu_e = np.einsum("ei,eij,ej->e", u_e, K_e, u_e, optimize="optimal")
145 c = (p_e**penal * uKu_e).sum()
146
147 # compute sensitivity for elements
148 dCdP_e = -(penal * (p_e ** (penal - 1.0)) * uKu_e)
149
150 # use sensitivity filter
151 dCdP_e = np.einsum("ij,j,j", Hij, p_e, dCdP_e) / (
152 np.einsum("i,ij", p_e, Hij) + 1e-12
153 )
154
155 # OC update (enforce volume)
156 lmin, lmax = 0.0, 1e5
157 pmin, pmax = 0.001, 1.0
158 move = 0.2
159 pNew_e = p_e.copy()
160
161 while (lmax - lmin) > 1e-4 * (lmax + lmin + 1e-16):
162 lmid = 0.5 * (lmax + lmin)
163 candidate = p_e * np.sqrt(np.maximum(-dCdP_e / lmid, 1e-9))
164 # Apply move limits and physical bounds [pmin, pmax]
165 pNew_e = np.maximum(
166 pmin,
167 np.maximum(
168 p_e - move,
169 np.minimum(pmax, np.minimum(p_e + move, candidate)),
170 ),
171 )
172 # update lambda to fit volume fraction
173 if pNew_e.sum() - volFrac * mesh.Ne > 0:
174 lmin = lmid
175 else:
176 lmax = lmid
177
178 # get updated density and compliance
179 p_e = pNew_e
180 list_p_e.append(p_e)
181 list_compliance.append(c)
182
183 # compute relative error : || p_e - pOld_e || / || pOld_e ||
184 err = np.linalg.norm(p_e - pOld_e) / np.linalg.norm(pOld_e)
185
186 Display.MyPrint(
187 f"Iteration {str(iter).zfill(len(str(iterMax)))}, compliance = {c:.3e}, volume fraction = {p_e.mean():.3f}, err = {err:.3e}",
188 end="\r",
189 )
190
191 # ----------------------------------------------
192 # Results
193 # ----------------------------------------------
194
195 axC = Display.Init_Axes()
196 axC.plot(range(len(list_compliance)), list_compliance, ls="-", marker=".")
197 axC.set_xlabel("Iteration")
198 axC.set_ylabel("Compliance")
199 Display.plt.show()
200
201 PyVista.Plot_BoundaryConditions(simu).show()
202
203 def get_thresh(p_e: np.ndarray, min=0.5, max=1.0):
204 grid = PyVista._pvMesh(mesh, p_e, nodeValues=False)
205 for result in simu.Results_Available():
206 grid[result] = simu.Result(result).reshape(mesh.Nn, -1)
207 thresh = grid.threshold((min, max))
208 return thresh
209
210 thresh = get_thresh(p_e)
211 PyVista.Plot(thresh, color="k").show()
212 PyVista.Plot(thresh, "uy").show()
213
214 if generateMovie:
215
216 def Func(plotter: PyVista.pv.Plotter, iter):
217 simu.Set_Iter(iter)
218 thresh = get_thresh(list_p_e[iter])
219 plotter.add_title(
220 f"{str(iter+1).zfill(len(str(iterMax)))}/{len(list_compliance)}"
221 )
222 PyVista.Plot(thresh, color="k", plotter=plotter)
223
224 PyVista.Movie_func(Func, len(list_compliance), folder, "optim.gif")
Total running time of the script: (0 minutes 10.770 seconds)