Note

Go to the end to download the full example code.

Full robust formulation

Robust formulation, where each of the eroded, nominal, and dilated designs are evaluated.

This is the ‘expensive’ version of the robust formulation for compliance problems. The poor man’s version can be found in Robust formulation, which only evaluates the eroded version.

In this particular example, Python functions are used to generate the repeated sequence of modules used for geometry projection and performance evaluation of each design.

References:

Wang, F., Lazarov, B. S., & Sigmund, O. (2011). On projection methods, convergence and robust formulations in topology optimization. Structural and multidisciplinary optimization, 43, 767-784. DOI: https://doi.org/10.1007/s00158-010-0602-y

 18 import numpy as np
 19 import pymoto as pym
 20
 21
 22 class Continuation(pym.Module):
 23     """ Module that generates a continuated value """
 24
 25     def __init__(self, start=0.0, stop=1.0, nsteps=80, stepstart=10, interval=10):
 26         self.startval = start
 27         self.endval = stop
 28         self.interval = interval
 29         self.dval = (stop - start) / nsteps
 30         self.nstart = stepstart
 31         self.iter = -1
 32         self.val = self.startval
 33
 34     def __call__(self):
 35         maxval = max(self.startval, self.endval)
 36         minval = min(self.startval, self.endval)
 37         if self.iter % self.interval == 0:  # Only update value every `interval` iterations
 38             self.val = np.clip(self.startval + self.dval * (self.iter - self.nstart), minval, maxval)
 39
 40         self.iter += 1
 41         return self.val
 42
 43
 44 if __name__ == '__main__':
 45     print(__doc__)
 46
 47     # Define the network as a normal python function
 48     def my_function(sx_in, s_beta, f, domain, eta=0.5, xmin=1e-9, bc=None):
 49         """Generate Heaviside projected eigenvalue calculation
 50         """
 51         # Heaviside projection
 52         heaviside = "(tanh(inp1 * {0}) + tanh(inp1 * (inp0 - {0}))) / (tanh(inp1 * {0}) + tanh(inp1 * (1 - {0})))"
 53         sx_projected = pym.MathExpression(heaviside.format(eta))(sx_in, s_beta)
 54
 55         # SIMP material interpolation
 56         sx_SIMP = pym.MathExpression(f"{xmin} + {1-xmin}*inp0^3")(sx_projected)
 57
 58         # Stiffness matrix assembly
 59         sK = pym.AssembleStiffness(domain, bc=bc)(sx_SIMP)
 60
 61         # Eigenvalue calculation
 62         s_u = pym.LinSolve()(sK, f)
 63
 64         # Compliance calculation
 65         s_compl = pym.EinSum('i,i->')(s_u, f)
 66
 67         return s_compl, sx_projected
 68
 69     # Now that the function is defined, it can be used as a 'recipe' for various inputs. The function will add the
 70     # modules that are requested and return the signal for 'frequency' in this case.
 71
 72     # Setup a domain
 73     domain_2d = pym.VoxelDomain(100, 50)
 74     bc = domain_2d.get_dofnumber(domain_2d.nodes[0, :], ndof=2)
 75     f = np.zeros(domain_2d.nnodes*2)
 76     f[domain_2d.nodes[-1, 10]*2 + 1] = 1
 77
 78     # Setup design vector
 79     sx = pym.Signal('x', np.ones(domain_2d.nel) * 0.5)
 80     sxfilt = pym.DensityFilter(domain_2d, radius=2)(sx)
 81
 82     # Continuation parameter `beta`
 83     sbeta = Continuation(start=1e-3, stop=40.0, stepstart=10, nsteps=80, interval=5)()
 84     sbeta.tag = "beta"
 85
 86     # Generate the three designs and evaluation for those
 87     s_compl_ero, sx_ero = my_function(sxfilt, sbeta, f, domain_2d, eta=0.7, bc=bc)
 88     s_compl_nom, sx_nom = my_function(sxfilt, sbeta, f, domain_2d, eta=0.5, bc=bc)
 89     s_compl_dil, sx_dil = my_function(sxfilt, sbeta, f, domain_2d, eta=0.3, bc=bc)
 90
 91     # After this you can mix and match any of the outputs of the functions, to use for an optimization.
 92
 93     # Plot the nominal domain
 94     pym.PlotDomain(domain_2d)(sx_nom)
 95
 96     # In this example the root mean square value of the compliances is minimized.
 97     s_compl_avg = pym.MathExpression("sqrt(inp0^2 + inp1^2 + inp2^2)")(s_compl_ero, s_compl_nom, s_compl_dil)
 98     s_obj = pym.Scaling(scaling=100.0)(s_compl_avg)
 99
100     # Perform the optimization
101     pym.minimize_oc(sx, s_obj, maxit=100)
102
103     # You can also evaluate the compliance for a range of eta values using the same function again
104     eta_array = np.linspace(0.2, 0.8, 30)
105     compliance_array = np.zeros(eta_array.size)
106     for i, eta in enumerate(eta_array):
107         s_compl, _ = my_function(sxfilt, sbeta, f, domain_2d, eta=eta, bc=bc)
108         compliance_array[i] = s_compl.state
109
110     # Show the compliance vs eta parameter. From this plot can clearly be seen that high eta (dilation) is always worse
111     # for performance in the case of compliance optimization. For other types of optimization (e.g. eigenfrequency) this
112     # is not the case. Try it out for yourself if you can see this behavior by changing the function to do an
113     # eigensolve to calculate eigenfrequencies.
114     import matplotlib.pyplot as plt
115     plt.figure()
116     plt.plot(eta_array, compliance_array)
117     plt.xlabel('eta')
118     plt.ylabel('Compliance')
119     plt.show()

Gallery generated by Sphinx-Gallery