.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/WeakForms/TopologyOptimisation1.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_WeakForms_TopologyOptimisation1.py: 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. .. GENERATED FROM PYTHON SOURCE LINES 13-228 .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/WeakForms/images/sphx_glr_TopologyOptimisation1_001.png :alt: TopologyOptimisation1 :srcset: /examples/WeakForms/images/sphx_glr_TopologyOptimisation1_001.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.9.2/docs/examples/WeakForms/images/sphx_glr_TopologyOptimisation1_001.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/WeakForms/images/sphx_glr_TopologyOptimisation1_002.png :alt: TopologyOptimisation1 :srcset: /examples/WeakForms/images/sphx_glr_TopologyOptimisation1_002.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.9.2/docs/examples/WeakForms/images/sphx_glr_TopologyOptimisation1_002.vtksz .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/WeakForms/images/sphx_glr_TopologyOptimisation1_003.png :alt: TopologyOptimisation1 :srcset: /examples/WeakForms/images/sphx_glr_TopologyOptimisation1_003.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.9.2/docs/examples/WeakForms/images/sphx_glr_TopologyOptimisation1_003.vtksz .. image-sg:: /examples/WeakForms/images/sphx_glr_TopologyOptimisation1_004.gif :alt: TopologyOptimisation1 :srcset: /examples/WeakForms/images/sphx_glr_TopologyOptimisation1_004.gif :class: sphx-glr-single-img .. image-sg:: /examples/WeakForms/images/sphx_glr_TopologyOptimisation1_005.png :alt: TopologyOptimisation1 :srcset: /examples/WeakForms/images/sphx_glr_TopologyOptimisation1_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none 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 %) 5.83 s Generate movie 02/44 (4.55 %) 5.07 s Generate movie 03/44 (6.82 %) 4.88 s Generate movie 04/44 (9.09 %) 4.80 s Generate movie 05/44 (11.36 %) 4.73 s Generate movie 06/44 (13.64 %) 4.68 s Generate movie 07/44 (15.91 %) 4.54 s Generate movie 08/44 (18.18 %) 4.39 s Generate movie 09/44 (20.45 %) 4.24 s Generate movie 10/44 (22.73 %) 4.09 s Generate movie 11/44 (25.00 %) 4.02 s Generate movie 12/44 (27.27 %) 3.90 s Generate movie 13/44 (29.55 %) 3.79 s Generate movie 14/44 (31.82 %) 3.70 s Generate movie 15/44 (34.09 %) 3.58 s Generate movie 16/44 (36.36 %) 3.42 s Generate movie 17/44 (38.64 %) 3.33 s Generate movie 18/44 (40.91 %) 3.16 s Generate movie 19/44 (43.18 %) 3.15 s Generate movie 20/44 (45.45 %) 2.98 s Generate movie 21/44 (47.73 %) 2.84 s Generate movie 22/44 (50.00 %) 2.82 s Generate movie 23/44 (52.27 %) 2.65 s Generate movie 24/44 (54.55 %) 2.49 s Generate movie 25/44 (56.82 %) 2.38 s Generate movie 26/44 (59.09 %) 2.21 s Generate movie 27/44 (61.36 %) 2.13 s Generate movie 28/44 (63.64 %) 2.02 s Generate movie 29/44 (65.91 %) 1.88 s Generate movie 30/44 (68.18 %) 1.76 s Generate movie 31/44 (70.45 %) 1.62 s Generate movie 32/44 (72.73 %) 1.48 s Generate movie 33/44 (75.00 %) 1.38 s Generate movie 34/44 (77.27 %) 1.24 s Generate movie 35/44 (79.55 %) 1.12 s Generate movie 36/44 (81.82 %) 997.30 ms Generate movie 37/44 (84.09 %) 869.26 ms Generate movie 38/44 (86.36 %) 771.05 ms Generate movie 39/44 (88.64 %) 630.34 ms Generate movie 40/44 (90.91 %) 503.02 ms Generate movie 41/44 (93.18 %) 372.40 ms Generate movie 42/44 (95.45 %) 249.31 ms Generate movie 43/44 (97.73 %) 124.74 ms Generate movie 44/44 (100.00 %) 0.00 µs | .. code-block:: Python :lineno-start: 14 import matplotlib.pyplot as plt import numpy as np from EasyFEA import Display, Folder, PyVista, ElemType, Models, Simulations from EasyFEA.FEM import FeArray, Field, BiLinearForm, Sym_Grad, Trace from EasyFEA.Geoms import Domain if __name__ == "__main__": Display.Clear() # ---------------------------------------------- # Configuration # ---------------------------------------------- dim = 2 L, H = 60, 30 # L, H = 120, 60 # optim topo iterMax = 60 volFrac = 0.4 penal = 3 rMin = 3 # outputs generateMovie = True folder = Folder.Results_Dir() # ---------------------------------------------- # Mesh # ---------------------------------------------- meshSize = 1 if dim == 2 else H / 10 contour = Domain((0, 0), (L, H), meshSize) assert H / meshSize % 2 == 0 if dim == 2: mesh = contour.Mesh_2D([], ElemType.QUAD4, isOrganised=True) else: mesh = contour.Mesh_Extrude( [], [0, 0, H], [H / meshSize], ElemType.HEXA8, isOrganised=True ) nodesX0 = mesh.Nodes_Conditions(lambda x, y, z: x == 0) zMean = 0 if dim == 2 else H / 2 nodesLoad = mesh.Nodes_Point((L, H / 2, zMean)) # nodesLoad = mesh.Nodes_Conditions(lambda x, y, z: y == H) # ---------------------------------------------- # Mesh-Independence Sensitivity Filter (Sigmund, 1998) # ---------------------------------------------- # get the coordinates of each elements coord_e = mesh.coord[mesh.connect].mean(1) # compute Hij elements = mesh.groupElem.elements Hij = np.array( [ np.maximum(0, rMin - np.linalg.norm(coord_e[i] - coord_e, axis=-1) + 1e-12) for i in range(mesh.Ne) ] ) # # plot neighbor elements # elem = 100 # ax = Display.Plot_Mesh(mesh, alpha=0) # Display.Plot_Elements( # mesh, mesh.connect[Hij[elem] != 0].ravel(), 2, color="blue", ax=ax # ) # Display.Plot_Elements(mesh, mesh.connect[elem].ravel(), 2, ax=ax) # ---------------------------------------------- # Formulations # ---------------------------------------------- elastic = Models.Elastic.Isotropic(dim, E=1, v=0.3, planeStress=True) mu = elastic.get_mu() lmbda = elastic.get_lambda() def S(u: Field) -> FeArray: Eps = Sym_Grad(u) return 2 * mu * Eps + lmbda * Trace(Eps) * np.eye(dim) p_e = np.ones(mesh.Ne, dtype=float) * volFrac @BiLinearForm def ComputeK(u: Field, v: Field): Sig = S(u) Eps = Sym_Grad(v) return Sig.ddot(Eps) @BiLinearForm def ComputePenalizedK(u: Field, v: Field): simpScaling = FeArray.asfearray(np.reshape(p_e**penal, (-1, 1))) return simpScaling * ComputeK(u, v) field = Field(mesh.groupElem, dim) model = Models.WeakForms(field, ComputePenalizedK, thickness=H) # ---------------------------------------------- # Simulation # ---------------------------------------------- simu = Simulations.WeakForms(mesh, model) simu.add_dirichlet(nodesX0, [0] * dim, simu.Get_unknowns()) simu.add_neumann(nodesLoad, [-1], ["y"]) # ---------------------------------------------- # Optim topo # ---------------------------------------------- err = 1.0 list_compliance: list[float] = [] list_p_e: list[np.ndarray] = [] iter = 0 while err > 0.005 and iter < iterMax: iter += 1 pOld_e = p_e.copy() # solve u simu.Need_Update() u = simu.Solve() simu.Save_Iter() # compute compliance for elements u_e = field.groupElem.Locates_sol_e(u, dim) K_e = ComputeK.Integrate_e(field) uKu_e = np.einsum("ei,eij,ej->e", u_e, K_e, u_e, optimize="optimal") c = (p_e**penal * uKu_e).sum() # compute sensitivity for elements dCdP_e = -(penal * (p_e ** (penal - 1.0)) * uKu_e) # use sensitivity filter dCdP_e = np.einsum("ij,j,j", Hij, p_e, dCdP_e) / ( np.einsum("i,ij", p_e, Hij) + 1e-12 ) # OC update (enforce volume) lmin, lmax = 0.0, 1e5 pmin, pmax = 0.001, 1.0 move = 0.2 pNew_e = p_e.copy() while (lmax - lmin) > 1e-4 * (lmax + lmin + 1e-16): lmid = 0.5 * (lmax + lmin) candidate = p_e * np.sqrt(np.maximum(-dCdP_e / lmid, 1e-9)) # Apply move limits and physical bounds [pmin, pmax] pNew_e = np.maximum( pmin, np.maximum( p_e - move, np.minimum(pmax, np.minimum(p_e + move, candidate)), ), ) # update lambda to fit volume fraction if pNew_e.sum() - volFrac * mesh.Ne > 0: lmin = lmid else: lmax = lmid # get updated density and compliance p_e = pNew_e list_p_e.append(p_e) list_compliance.append(c) # compute relative error : || p_e - pOld_e || / || pOld_e || err = np.linalg.norm(p_e - pOld_e) / np.linalg.norm(pOld_e) Display.MyPrint( f"Iteration {str(iter).zfill(len(str(iterMax)))}, compliance = {c:.3e}, volume fraction = {p_e.mean():.3f}, err = {err:.3e}", end="\r", ) # ---------------------------------------------- # Results # ---------------------------------------------- axC = Display.Init_Axes() axC.plot(range(len(list_compliance)), list_compliance, ls="-", marker=".") axC.set_xlabel("Iteration") axC.set_ylabel("Compliance") plt.show() PyVista.Plot_BoundaryConditions(simu).show() def get_thresh(p_e: np.ndarray, min=0.5, max=1.0): grid = PyVista._pvMesh(mesh, p_e, nodeValues=False) for result in simu.Results_Available(): grid[result] = simu.Result(result).reshape(mesh.Nn, -1) thresh = grid.threshold((min, max)) return thresh thresh = get_thresh(p_e) PyVista.Plot(thresh, color="k").show() PyVista.Plot(thresh, "uy").show() if generateMovie: def Func(plotter: PyVista.pv.Plotter, iter): simu.Set_Iter(iter) thresh = get_thresh(list_p_e[iter]) plotter.add_title( f"{str(iter+1).zfill(len(str(iterMax)))}/{len(list_compliance)}" ) PyVista.Plot(thresh, color="k", plotter=plotter) PyVista.Movie_func(Func, len(list_compliance), folder, "optim.gif") .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 14.287 seconds) .. _sphx_glr_download_examples_WeakForms_TopologyOptimisation1.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: TopologyOptimisation1.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: TopologyOptimisation1.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: TopologyOptimisation1.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_