Note

Go to the end to download the full example code.

Overhang filter

In this example the usage of an overhang filter is demonstrated to obtain a printable design

The module pymoto.OverhangFilter removes all overhanging material, which forces the optimizer to add supports to load-bearing parts of the design. Compared to the example Compliance minimization, only the overhang filter is added (see line 64-65 of this example).

References:

  • Langelaar, M. (2017). An additive manufacturing filter for topology optimization of print-ready designs. Structural and Multidisciplinary Optimization, 55(3), 871-883. doi: 10.1007/s00158-016-1522-2

  • Langelaar, M. (2016). Topology optimization of 3D self-supporting structures for additive manufacturing. Additive Manufacturing, 12, 60-70. doi: 10.1016/j.addma.2016.06.010

  • Delissen, A. et al. (2022). Realization and assessment of metal additive manufacturing and topology optimization for high-precision motion systems. Additive Manufacturing, 58, 103012. doi: 10.1016/j.addma.2022.103012

 22 import numpy as np
 23
 24 import pymoto as pym
 25
 26 nx, ny, nz = 120, 40, 0  # Set nz to zero for the 2D problem, nz > 0 runs a 3D problem
 27 xmin = 1e-9
 28 filter_radius = 2.0
 29 volfrac = 0.5
 30
 31 if __name__ == "__main__":
 32     print(__doc__)
 33
 34     if nz == 0:  # 2D analysis
 35         # Generate a grid
 36         domain = pym.VoxelDomain(nx, ny)
 37         ndof = domain.dim
 38         # Calculate boundary dof indices
 39         boundary_dofs = domain.get_dofnumber(domain.nodes[0, :], np.arange(ndof), ndof=ndof).flatten()
 40
 41         # Which dofs to put a force on? The 1 is added for a force in y-direction (x-direction would be zero)
 42         force_dofs = ndof * domain.nodes[nx, int(ny/2)] + 1
 43
 44     else:
 45         domain = pym.VoxelDomain(nx, ny, nz)
 46         ndof = domain.dim
 47         boundary_dofs = domain.get_dofnumber(domain.nodes[0, :, :], np.arange(ndof), ndof=ndof).flatten()
 48         force_dofs = ndof * domain.nodes[nx, :, int(nz/2)] + 2  # Z-direction
 49
 50     # Generate a force vector
 51     f = np.zeros(domain.nnodes * ndof)
 52     f[force_dofs] = 1.0  # Uniform force of 1.0 at all selected dofs
 53
 54     # Make signal for design vector, and fill with initial values
 55     sx = pym.Signal('x', state=np.ones(domain.nel) * volfrac)
 56
 57     # Start building the modular network
 58     # Filter
 59     sxfilt = pym.FilterConv(domain, radius=filter_radius, ymin_bc=0, ymax_bc=0, zmin_bc=0, zmax_bc=0)(sx)
 60     sxfilt.tag = "Filtered design"
 61
 62     # -----> Only the :py:class:`pymoto.OverhangFilter` module needs to be added after filtering. The rest of the
 63     # example is equal to :ref:`sphx_glr_auto_examples_topology_optimization_ex_compliance.py`
 64     sxoverhang = pym.OverhangFilter(domain, direction='+y')(sxfilt)
 65     sxoverhang.tag = "Printable design"
 66     # <-----
 67
 68     sx_analysis = sxoverhang
 69
 70     # Show the design on the screen as it optimizes
 71     if domain.dim == 2:
 72         pym.PlotDomain(domain, saveto="out/design", clim=[0, 1])(sx_analysis)
 73
 74     # SIMP material interpolation
 75     sSIMP = pym.MathExpression(f"{xmin} + {1.0 - xmin}*inp0^3")(sx_analysis)
 76
 77     # System matrix assembly module
 78     sK = pym.AssembleStiffness(domain, bc=boundary_dofs)(sSIMP)
 79
 80     # Linear system solver
 81     su = pym.LinSolve(symmetric=True, positive_definite=True)(sK, f)
 82
 83     # # Output the design, deformation, and force field to a Paraview file
 84     pym.WriteToVTI(domain, saveto='out/dat.vti')(sx_analysis, su, f)
 85
 86     # Compliance calculation c = f^T u
 87     scompl = pym.EinSum('i, i->')(su, f)
 88     scompl.tag = 'compliance'
 89
 90     # MMA needs correct scaling of the objective
 91     sg0 = pym.Scaling(scaling=100.0)(scompl)
 92     sg0.tag = "objective"
 93
 94     # Calculate the volume of the domain by adding all design densities together
 95     svol = pym.EinSum('i->')(sx_analysis)
 96     svol.tag = 'volume'
 97
 98     # Volume constraint
 99     sg1 = pym.MathExpression(f'10*(inp0/{domain.nel} - {volfrac})')(svol)
100     sg1.tag = "volume constraint"
101
102     # Maybe you want to check the design-sensitivities?
103     do_finite_difference = False
104     if do_finite_difference:
105         pym.finite_difference(sx, [sg0, sg1], dx=1e-4)
106         exit()
107
108     pym.PlotIter()(sg0, sg1)  # Plot iteration history
109
110     # Do the optimization with MMA
111     pym.minimize_mma(sx, [sg0, sg1])
112
113     # Here you can do some post processing
114     print("The optimization has finished!")
115     print(f"The final compliance value obtained is {scompl.state}")
116     print(f"The maximum displacement is {max(np.absolute(su.state))}")

Gallery generated by Sphinx-Gallery