Note

Go to the end to download the full example code.

Compliance minimization

This example demonstrates how to minimize the compliance of a 2D or 3D structure using pymoto

Two different physics are implemented: a static mechanical analysis (= stiffness maximization) and a static thermal analysis (= conductivity maximization). The example contains the essential modules needed to run a topology optimization problem:

And next to that some important functions that operate on the pymoto.Network:

 26 import numpy as np
 27
 28 import pymoto as pym
 29
 30 nx, ny, nz = 120, 40, 0  # Set nz to zero for the 2D problem, nz > 0 runs a 3D problem
 31 xmin = 1e-9
 32 filter_radius = 2.0
 33 volfrac = 0.5
 34 thermal = False  # True = Static thermal analysis; False = Static mechanical analysis will be done
 35
 36 if __name__ == "__main__":
 37     print(__doc__)
 38
 39     if nz == 0:  # 2D analysis
 40         # Generate a grid
 41         domain = pym.VoxelDomain(nx, ny)
 42
 43         if thermal:
 44             ndof = 1  # Number of dofs per node
 45             # Get dof numbers at the boundary
 46             boundary_dofs = domain.nodes[0, int(ny/4):-int(ny/4)].flatten()
 47
 48             # Make a force vector
 49             force_dofs = domain.nodes[1:,:].flatten()
 50
 51         else:  # Mechanical
 52             ndof = 2
 53             # Calculate boundary dof indices
 54             boundary_dofs = domain.get_dofnumber(domain.nodes[0, :], np.arange(ndof), ndof=ndof).flatten()
 55
 56             # Which dofs to put a force on? The 1 is added for a force in y-direction (x-direction would be zero)
 57             force_dofs = ndof * domain.nodes[nx, int(ny/2)] + 1
 58
 59     else:
 60         domain = pym.VoxelDomain(nx, ny, nz)
 61
 62         if thermal:
 63             ndof = 1
 64             boundary_dofs = domain.nodes[0, int(3*ny/8):-int(3*ny/8), int(3*nz/8):-int(3*nz/8)].flatten()
 65             force_dofs = domain.nodes[int(nx/4):, :, :].flatten()
 66         else:  # Mechanical
 67             ndof = 3
 68             boundary_dofs = domain.get_dofnumber(domain.nodes[0, :, :], np.arange(ndof), ndof=ndof).flatten()
 69             force_dofs = ndof * domain.nodes[nx, :, int(nz/2)] + 2  # Z-direction
 70
 71     # Generate a force vector
 72     f = np.zeros(domain.nnodes * ndof)
 73     f[force_dofs] = 1.0  # Uniform force of 1.0 at all selected dofs
 74
 75     # Make signal for design vector, and fill with initial values
 76     sx = pym.Signal('x', state=np.ones(domain.nel) * volfrac)
 77
 78     # Start building the modular network
 79     # Filter
 80     sxfilt = pym.FilterConv(domain, radius=filter_radius, ymin_bc=0, ymax_bc=0, zmin_bc=0, zmax_bc=0)(sx)
 81     sxfilt.tag = "Filtered design"
 82     sx_analysis = sxfilt
 83
 84     # Show the design on the screen as it optimizes
 85     if domain.dim == 2:
 86         pym.PlotDomain(domain, saveto="out/design", clim=[0, 1])(sx_analysis)
 87
 88     # SIMP material interpolation
 89     sSIMP = pym.MathExpression(f"{xmin} + {1.0 - xmin}*inp0^3")(sx_analysis)
 90
 91     # System matrix assembly module
 92     if thermal:
 93         sK = pym.AssemblePoisson(domain, bc=boundary_dofs)(sSIMP)
 94     else:
 95         sK = pym.AssembleStiffness(domain, bc=boundary_dofs)(sSIMP)
 96
 97     # Linear system solver. The linear solver can be chosen by uncommenting any of the following lines.
 98     # The matrix properties 'symmetric' and 'positive_definite' are passed to the LinSolve module, which are optional
 99     # but can improve the solver choice. This may increase the speed of solving the linear system.
100     solver = None  # Default (solver is automatically chosen based on matrix properties)
101     # solver = pym.solvers.SolverSparsePardiso()  # Requires Intel MKL installed
102     # solver = pym.solvers.SolverSparseCholeskyCVXOPT()  # Requires cvxopt installed
103     # solver = pym.solvers.SolverSparseCholeskyScikit()  # Requires scikit installed
104     su = pym.LinSolve(symmetric=True, positive_definite=True, solver=solver)(sK, f)
105
106     # # Output the design, deformation, and force field to a Paraview file
107     pym.WriteToVTI(domain, saveto='out/dat.vti')(sx_analysis, su, f)
108
109     # Compliance calculation c = f^T u
110     scompl = pym.EinSum('i, i->')(su, f)
111     scompl.tag = 'compliance'
112
113     # MMA needs correct scaling of the objective
114     sg0 = pym.Scaling(scaling=100.0)(scompl)
115     sg0.tag = "objective"
116
117     # Calculate the volume of the domain by adding all design densities together
118     svol = pym.EinSum('i->')(sx_analysis)
119     svol.tag = 'volume'
120
121     # Volume constraint
122     sg1 = pym.MathExpression(f'10*(inp0/{domain.nel} - {volfrac})')(svol)
123     sg1.tag = "volume constraint"
124
125     # Maybe you want to check the design-sensitivities?
126     do_finite_difference = False
127     if do_finite_difference:
128         pym.finite_difference(sx, [sg0, sg1], dx=1e-4)
129         exit()
130
131     pym.PlotIter()(sg0, sg1)  # Plot iteration history
132
133     # Do the optimization with MMA
134     # pym.minimize_mma(sx, [sg0, sg1])
135
136     # Do the optimization with OC
137     pym.minimize_oc(sx, sg0)
138
139     # Here you can do some post processing
140     print("The optimization has finished!")
141     print(f"The final compliance value obtained is {scompl.state}")
142     print(f"The maximum {'temperature' if thermal else 'displacement'} is {max(np.absolute(su.state))}")

Gallery generated by Sphinx-Gallery