Note
Go to the end to download the full example code.
Showcase for 3D MMB beam compliance
This example showcases the topology optimization of a 3D MBB beam (i.e. bridge-like) structure
The structure is subjected to 3-point bending and the optimization seeks to maximize the stiffness. As the problem is 2-fold symmetric, only a quarter of the domain needs analysis and optimization. With post-processing the design can be mirrored to restore the complete domain (e.g. using Paraview).
This example is based on Multigrid preconditioned CG solver.
12 import numpy as np
13 import pymoto as pym
14
15 n = 32 # Resolution (# elements). Use powers of 2 for good GMG performance (e.g. 4, 8, 16, ...)
16 nx, ny, nz = 4*n, n, 2*n # Domain sizes (# elements)
17 Lx = 10.0 # Length in (m), other dimensions are deduced from `ny` and `nz`
18 xmin = 1e-9 # Minimum density value
19 filter_radius = 2.0 # Filter radius
20 volfrac = 0.1 # Target volume fraction
21
22 if __name__ == "__main__":
23 print(__doc__)
24
25 # Setup domain
26 h = Lx/nx # Element size
27 domain = pym.VoxelDomain(nx, ny, nz, unitx=h, unity=h, unitz=h) # Setup domain object
28 ndof = 3 # Number of dofs per node
29
30 # Symmetry BC in x-direction
31 nidx_xmin = domain.nodes[0, :, :].flatten()
32 bc_x = domain.get_dofnumber(nidx_xmin, dof_idx=0, ndof=ndof)
33
34 # Symmetry BC in y-direction
35 nidx_ymin = domain.nodes[:, 0, :].flatten()
36 bc_y = domain.get_dofnumber(nidx_ymin, dof_idx=1, ndof=ndof)
37
38 # Simple support in z-direction
39 nidx_supp = domain.nodes[-1, :, 0].flatten()
40 bc_z = domain.get_dofnumber(nidx_supp, dof_idx=2, ndof=ndof)
41
42 # Combine all boundary conditions
43 boundary_dofs = np.concatenate([bc_x, bc_y, bc_z])
44
45 # Load vector on top in z-direction
46 delem = int(np.ceil(filter_radius))
47 eidx_force = domain.elements[:delem, :, -delem:].flatten()
48 nidx_force = domain.conn[eidx_force, :].flatten()
49 didx_force = domain.get_dofnumber(nidx_force, dof_idx=2, ndof=ndof)
50 f = np.zeros(domain.nnodes * ndof)
51 np.add.at(f, didx_force, 1.0)
52 f *= -1 / (4*np.sum(f)) # Normalize to -1/4N unit force (-1N for the total symmetric domain)
53
54 # Make design vector and fill with initial values
55 sx = pym.Signal('x', state=np.ones(domain.nel) * volfrac)
56
57 # Padding elements for density filter (set densities to 0)
58 delem0 = delem - 1
59 eidx_xmax = domain.elements[-delem0:, :, :].flatten()
60 eidx_ymax = domain.elements[:, -delem0:, :].flatten()
61 eidx_zmin = domain.elements[:, :, :delem0].flatten()
62 eidx_zmax = domain.elements[:, :, -delem0:].flatten()
63 eidx_0 = np.unique(np.concatenate([eidx_xmax, eidx_ymax, eidx_zmin, eidx_zmax]))
64
65 # Solid elements at force and simple-support location (set densities to 1)
66 eidx_bc = domain.elements[-delem:, :, :delem].flatten()
67 eidx_1 = np.unique(np.concatenate([eidx_force, eidx_bc]))
68
69 # Start building the modular network
70 with pym.Network() as func:
71 # Set elements at the edge to zero to prevent the design 'sticking' to the boundary
72 sx0 = pym.SetValue(indices=eidx_0, value=0.0)(sx)
73
74 # Density filter
75 sxfilt = pym.FilterConv(domain, radius=filter_radius, xmax_bc=0, ymax_bc=0, zmin_bc=0, zmax_bc=0)(sx0)
76
77 # Set elements to one at location where the force acts and where the simple support is
78 sx1 = pym.SetValue(indices=eidx_1, value=1.0)(sxfilt)
79 sx_analysis = sx1
80 sx_analysis.tag = 'design'
81
82 # SIMP material interpolation
83 sSIMP = pym.MathExpression(f"{xmin} + {1.0 - xmin}*inp0^3")(sx_analysis)
84
85 # System matrix assembly module
86 sK = pym.AssembleStiffness(domain, bc=boundary_dofs, e_modulus=70e+9, poisson_ratio=0.3)(sSIMP)
87
88 # Linear system solver. The linear solver can be chosen by uncommenting any of the following lines.
89 solver = None # Default (solver is automatically chosen based on matrix properties)
90 # solver = pym.solvers.SolverSparsePardiso() # Requires Intel MKL installed
91 # solver = pym.solvers.SolverSparseCholeskyCVXOPT() # Requires cvxopt installed
92 # solver = pym.solvers.SolverSparseCholeskyScikit() # Requires scikit installed
93
94 # Iterative solver: CG with geometric multigrid preconditioning
95 # Set up fine level multigrid operator
96 mg1 = pym.solvers.GeometricMultigrid(domain)
97 mgs = [mg1] # List of multigrid levels, starting with the finest grid
98 while True: # Setup multiple sub-levels
99 # Stop if any of the new dimensions is odd
100 if any([nn % 2 != 0 for nn in mgs[-1].sub_domain.size]):
101 break
102
103 # Stop if any of the coarse dimensions is smaller than 8
104 if any(mgs[-1].sub_domain.size < 8):
105 break
106
107 # Create a new coarser grid
108 mgs.append(pym.solvers.GeometricMultigrid(mgs[-1].sub_domain))
109 mgs[-2].inner_level = mgs[-1] # Set the inner level to the next coarser grid
110
111 print(f"Fine grid size: {domain.nelx} x {domain.nely} x {domain.nelz}")
112 print(f"Multigrid levels: {len(mgs)}")
113 print(f"Coarsest grid size: {mgs[-1].sub_domain.nelx} x {mgs[-1].sub_domain.nely} x {mgs[-1].sub_domain.nelz}")
114
115 # Set up the solver (comment out to use the default factorization, try this to see the difference in time)
116 solver = pym.solvers.CG(preconditioner=mg1, verbosity=0, tol=1e-5)
117
118 su = pym.LinSolve(hermitian=True, solver=solver)(sK, f)
119
120 # Output the design, deformation, and force field to a Paraview file
121 pym.WriteToVTI(domain, saveto='out/dat.vti')(sx_analysis, su)
122
123 # Compliance calculation c = f^T u
124 scompl = pym.EinSum('i,i->')(su, f)
125 scompl.tag = 'compliance'
126
127 # MMA needs correct scaling of the objective
128 sg0 = pym.Scaling(scaling=100.0)(scompl)
129 sg0.tag = "objective"
130
131 # Calculate the volume of the domain by adding all design densities together
132 svol = pym.EinSum('i->')(sx_analysis)
133 svol.tag = 'volume'
134
135 # Volume constraint
136 sg1 = pym.Scaling(scaling=10.0, maxval=volfrac*domain.nel)(svol)
137 sg1.tag = "volume constraint"
138
139 pym.PlotIter()(sg0, sg1) # Plot iteration history
140
141 # Do the optimization with MMA
142 # pym.minimize_mma([sx], [sg0, sg1], function=func)
143
144 # Do the optimization with OC
145 pym.minimize_oc(sx, sg0, function=func)
146
147 print("The optimization has finished!")
148 print(f"The final compliance value obtained is {scompl.state}")
149 print(f"The maximum displacement is {max(np.absolute(su.state))}")