Maximum Local Compliance
This loal compliance maximization designs structures with the maximum displacement in one node. This can be used to design MEMS actuators as is discussed at Maximum Local Compliance. An example as how to use the optimization is shown in an example optimization example.py
Density Constraints
Constraints class used to specify the density constraints of the topology optimisation problem. It contains functions for minimum and maximum element density in the upcoming iteration and the magnitude of the volume constraint function itself of the current design. This version of the code is used for the compliant design, local displacement maximisation.
Bram Lagerweij Aerospace Structures and Materials Department TU Delft 2018
-
class
src_Actuator.constraints.
DensityConstraint
(nelx, nely, move, volume_frac, density_min=0.0, density_max=1.0) This object relates to the constraints used in this optimization. It can be used for the MMA update scheme to derive what the limit is for all element densities at every iteration. The class itself is not changed by the iterations.
- Parameters
nelx (int) – Number of elements in x direction.
nely (int) – Number of elements in y direction.
move (float) – Maximum change in density of an element over 1 itteration.
volume_frac (float) – Maximum volume that can be filled with material.
volume_derivative (2D array size(1, nelx*nely)) – Sensitivity of the density constraint to the density in each element.
density_min (float (optional)) – Minimum density, set at 0.0 if not specified.
density_max (float (optional)) – Maximum density, set at 0.0 if not specified.
-
nelx
Number of elements in x direction.
- Type
int
-
nely
Number of elements in y direction.
- Type
int
-
move
Maximum change in density of an element over 1 iteration.
- Type
float
-
volume_frac
Maximum volume that can be filled with material.
- Type
float
-
volume_derivative
Sensitivity of the density constraint to the density in each element.
- Type
2D array size(1, nelx*nely)
-
density_min
Minimum density, set at 0.0 if not specified.
- Type
float, optional
-
density_max
Maximum density, set at 0.0 if not specified.
- Type
float, optional
-
current_volconstrain
(x) Calculates the current magnitude of the volume constraint function:
\[V_{\text{constraint}} = \frac{\sum v_e X_e}{ V_{\max}}-1\]- Parameters
x (2D array size(nely, nelx)) – Density distribution of this iteration.
- Returns
curvol – Current value of the density constraint function.
- Return type
float
-
xmax
(x) This function calculates the maximum density value of all elements of this iteration.
- Parameters
x (2D array size(nely, nelx)) – Density distribution of this iteration.
- Returns
xmax – Maximum density values of this itteration after updating.
- Return type
2D array size(nely, nelx)
-
xmin
(x) This function calculates the minimum density value of all elements of this iteration.
- Parameters
x (2D array size(nely, nelx)) – Density distribution of this iteration.
- Returns
xmin – Minimum density values of this iteration for the update scheme.
- Return type
2D array size(nely, nelx)
Load Cases
This file contains the Load class that allows the generation of an object that contains geometric, mesh, loads and boundary conditions that belong to the load case. This version of the code is meant for local compliant maximization.
Bram Lagerweij Aerospace Structures and Materials Department TU Delft 2018
Parent Load Case
-
class
src_Actuator.loads.
Load
(nelx, nely, young, Emin, poisson, ext_stiff) Load parent class that contains the basic functions used in all load cases. This class and its children do cantain information about the load case conciderd in the optimisation. The load case consists of the mesh, the loads, and the boundaries conditions. The class is constructed such that new load cases can be generated simply by adding a child and changing the function related to the geometry, loads and boundaries.
- Parameters
nelx (int) – Number of elements in x direction.
nely (int) – Number of elements in y direction.
young (float) – Youngs modulus of the materias.
Emin (float) – Artifical Youngs modulus of the material to ensure a stable FEA. It is used in the SIMP based material model.
poisson (float) – Poisson ration of the material.
-
nelx
Number of elements in x direction.
- Type
int
-
nely
Number of elements in y direction.
- Type
int
-
young
Youngs modulus of the materias.
- Type
float
-
Emin
Artifical Youngs modulus of the material to ensure a stable FEA. It is used in the SIMP based material model.
- Type
float
-
poisson
Poisson ration of the material.
- Type
float
-
dim
Amount of dimensions conciderd in the problem, set at 2.
- Type
int
-
ext_stiff
Extra stiffness to be added to global stiffness matrix. Due to interactions with meganisms outside design domain.
- Type
float
-
alldofs
() Returns a list with all degrees of freedom.
- Returns
all – List with numbers from 0 to the maximum degree of freedom number.
- Return type
1-D list
-
displaceloc
() Returns a zero vector, there is supposed to be an 1 implemented at the index where displacment output should be maximised, such that u·l = u_out
- Returns
l – Empty for the governing class.
- Return type
1-D column array length covering all degrees of freedom
-
edof
() Generates an array with the position of the nodes of each element in the global stiffness matrix.
- Returns
edof (2-D array size(nelx*nely, 8)) – The list with all elements and their degree of freedom numbers.
x_list (1-D array len(nelx*nely*8*8)) – The list with the x indices of all ellements to be inserted into the global stiffniss matrix.
y_list (1-D array len(nelx*nely*8*8)) – The list with the y indices of all ellements to be inserted into the global stiffniss matrix.
-
fixdofs
() Returns a list with indices that are fixed by the boundary conditions.
- Returns
fix – List with all the numbers of fixed degrees of freedom. This list is empty in this parrent class.
- Return type
1-D list
-
force
() Returns an 1D array, the force vector of the loading condition.
- Returns
f – Empy force vector.
- Return type
1-D column array length covering all degrees of freedom
-
freedofs
() Returns a list of arr indices that are not fixed
- Returns
free – List containing all elemens of alldogs except those that appear in the freedofs list.
- Return type
1-D list
-
lk
(E, nu) Calculates the local siffness matrix depending on E and nu.
- Parameters
E (float) – Youngs modulus of the material.
nu (float) – Poisson ratio of the material.
- Returns
ke – Local stiffness matrix.
- Return type
2-D array size(8, 8)
-
node
(elx, ely) Calculates the topleft node number of the requested element.
- Parameters
elx (int) – X position of the conciderd element.
ely (int) – Y position of the conciderd element.
- Returns
topleft – The node number of the top left node.
- Return type
int
-
nodes
(elx, ely) Calculates all node numbers of the requested element
- Parameters
elx (int) – X position of the conciderd element.
ely (int) – Y position of the conciderd element.
- Returns
n1 (int) – The node number of the top left node.
n2 (int) – The node number of the top right node.
n3 (int) – The node number of the bottom right node.
n4 (int) – The node number of the bottom left node.
-
passive
() Retuns three lists containing the location and magnitude of fixed density values
- Returns
elx (1-D list) – X coordinates of all passive elements, empty for the parrent class.
ely (1-D list) – Y ccordinates of all passive elements, empty for the parrent class.
values (1-D list) – Density values of all passive elements, empty for the parrent class.
fix_ele (1-D list) – List with all element numbers that are allowed to change.
Child Load Cases
-
class
src_Actuator.loads.
Inverter
(nelx, nely, young, Emin, poisson, ext_stiff) Bases:
src_Actuator.loads.Load
This child of the Load class represents a top half of the symetric inverter design used for MEMS actuators. It contains an positive horizontal force at the bottom left corner which causes a negative displacement at the bottom right corner.
No methods are added compared to the parrent class. Only the force, displaceloc and fixdof equations are changed to contain the propper values for the boundary conditions and optimisation objective.
-
displaceloc
() The maximisation should occur in negative x direction at the bottom right corner. Positive movement is thus in negative x direction.
- Returns
l – Value of -1 at the index related to the bottom right node.
- Return type
1-D column array length covering all degrees of freedom
-
fixdofs
() The boundary conditions of this problem fixes the bottom of the desing space in y direction (due to symetry). While the topleft element is fixed in both x and y direction on the left side.
- Returns
fix – List with all the numbers of fixed degrees of freedom.
- Return type
1-D list
-
force
() The force vector containts a load in positive x direction at the bottom left corner node.
- Returns
f – Value of 1 at the index related to the bottom left node.
- Return type
1-D column array length covering all degrees of freedom
-
Finite Element Solvers
Finite element solvers for the displacement from stiffness matrix, force and adjoin vector. This version of the code is meant for local compliant maximization.
Bram Lagerweij Aerospace Structures and Materials Department TU Delft 2018
Parent Solver
-
class
src_Actuator.fesolvers.
FESolver
(verbose=False) This parent FEA class can only assemble the global stiffness matrix and exclude all fixed degrees of freedom from it. This stiffness csc-sparse stiffness matrix is assembled in the gk_freedof method. This class solves the FE problem with a sparse LU-solver based upon umfpack. This solver is slow and inefficient. It is however more robust.
For this local compliance (actuator) maximization this solver solves two problems, the equilibrium and the adjoint problem which will be required to compute the gradients.
- Parameters
verbose (bool, optional) – False if the FEA should not print updates
-
verbose
False if the FEA should not print updates.
- Type
bool
-
displace
(load, x, ke, kmin, penal) FE solver based upon the sparse SciPy solver that uses umfpack.
- Parameters
load (object, child of the Loads class) – The loadcase(s) considered for this optimisation problem.
x (2-D array size(nely, nelx)) – Current density distribution.
ke (2-D array size(8, 8)) – Local fully dense stiffness matrix.
kmin (2-D array size(8, 8)) – Local stiffness matrix for an empty element.
penal (float) – Material model penalisation (SIMP).
- Returns
u (1-D column array shape(max(edof), 1)) – The displacement vector.
lamba (1-D column array shape(max(edof), 1)) – Adjoint equation solution.
-
gk_freedofs
(load, x, ke, kmin, penal) Generates the global stiffness matrix with deleted fixed degrees of freedom. It generates a list with stiffness values and their x and y indices in the global stiffness matrix. Some combination of x and y appear multiple times as the degree of freedom might appear in multiple elements of the FEA. The SciPy coo_matrix function adds them up at the background. At the location of the force introduction and displacement output an external stiffness is added due to stability reasons.
- Parameters
load (object, child of the Loads class) – The loadcase(s) considered for this optimisation problem.
x (2-D array size(nely, nelx)) – Current density distribution.
ke (2-D array size(8, 8)) – Local fully dense stiffness matrix.
kmin (2-D array size(8, 8)) – Local stiffness matrix for an empty element.
penal (float) – Material model penalisation (SIMP).
- Returns
k – Global stiffness matrix without fixed degrees of freedom.
- Return type
2-D sparse csc matrix
Child Solvers
-
class
src_Actuator.fesolvers.
CvxFEA
(verbose=False) Bases:
src_Actuator.fesolvers.FESolver
This parent FEA class can assemble the global stiffness matrix and solve the FE problem with a Supernodal Sparse Cholesky Factorization. It solves for both the equilibrium and adjoin problems.
-
verbose
False if the FEA should not print updates.
- Type
bool
-
displace
(load, x, ke, kmin, penal) FE solver based upon a Supernodal Sparse Cholesky Factorization. It requires the installation of the cvx module. It solves both the FEA equilibrium and adjoint problems. 1
- Parameters
load (object, child of the Loads class) – The loadcase(s) considerd for this optimisation problem.
x (2-D array size(nely, nelx)) – Current density distribution.
ke (2-D array size(8, 8)) – Local fully dense stiffness matrix.
kmin (2-D array size(8, 8)) – Local stiffness matrix for an empty element.
penal (float) – Material model penalisation (SIMP).
- Returns
u (1-D column array shape(max(edof), 1)) – The displacement vector.
lamba (1-D column array shape(max(edof), 1)) – Adjoint equation solution.
References
- 1
Y. Chen, T. A. Davis, W. W. Hager, S. Rajamanickam, “Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization and Update/Downdate”, ACM Transactions on Mathematical Software, 35(3), 22:1-22:14, 2008.
-
-
class
src_Actuator.fesolvers.
CGFEA
(verbose=False) Bases:
src_Actuator.fesolvers.FESolver
This parent FEA class can assemble the global stiffness matrix and solve the FE problem with a sparse solver based upon a preconditioned conjugate gradient solver. The preconditioning is based upon the inverse of the diagonal of the stiffness matrix.
Recommendations
Make the tolerance change over the iterations, low accuracy is required for first iteration, more accuracy for the later ones.
Add more advanced preconditioned.
Add gpu acceleration.
-
verbose
False if the FEA should not print updates.
- Type
bool
-
ufree_old
Displacement field of previous iteration.
- Type
array len(freedofs)
-
lambafree_old
Ajoint equation result of previous iteration.
- Type
array len(freedofs)
-
displace
(load, x, ke, kmin, penal) FE solver based upon the sparse SciPy solver that uses a preconditioned conjugate gradient solver, preconditioning is based upon the inverse of the diagonal of the stiffness matrix. Currently the relative tolerance is hardcoded as 1e-5. It solves both the equilibrium and adjoint problems.
- Parameters
load (object, child of the Loads class) – The loadcase(s) considered for this optimisation problem.
x (2-D array size(nely, nelx)) – Current density distribution.
ke (2-D array size(8, 8)) – Local fully dense stiffness matrix.
kmin (2-D array size(8, 8)) – Local stiffness matrix for an empty element.
penal (float) – Material model penalisation (SIMP).
- Returns
u (1-D array len(max(edof)+1)) – Displacement of all degrees of freedom
lamba (1-D column array shape(max(edof), 1)) – Adjoint equation solution.
Optimization Module
Topology Optimization class that handles the iterations, objective functions, filters and update scheme. It requires to call upon a constraint, load case and FE solver classes. This version of the code is meant for local compliant maximization (Actuator design).
Bram Lagerweij Aerospace Structures and Materials Department TU Delft 2018
-
class
src_Actuator.topopt.
Topopt
(constraint, load, fesolver, verbose=False) This is the optimisation object itself. It contains the initialisation of the density distribution.
- Parameters
constraint (object of DensityConstraint class) – The constraints for this optimization problem.
load (object, child of the Loads class) – The loadcase(s) considered for this optimisation problem.
fesolver (object, child of the CSCStiffnessMatrix class) – The finite element solver.
verbose (bool, optional) – Printing itteration results.
-
constraint
The constraints for this optimization problem.
- Type
object of DensityConstraint class
-
load
The loadcase(s) considered for this optimisation problem.
- Type
object, child of the Loads class
-
fesolver
The finite element solver.
- Type
object, child of the CSCStiffnessMatrix class
-
verbose
Printing iteration results.
- Type
bool
-
itr
Number of iterations performed
- Type
int
-
x
Array containing the current densities of every element.
- Type
2-D array size(nely, nelx)
-
xold1
Flattened density distribution one iteration ago.
- Type
1D array len(nelx*nely)
-
xold2
Flattened density distribution two iteration ago.
- Type
1D array len(nelx*nely)
-
low
Column vector with the lower asymptotes, calculated and used in the MMA subproblem of the previous iteration.
- Type
1D array len(nelx*nely)
-
upp
Column vector with the lower asymptotes, calculated and used in the MMA subproblem of the previous iteration.
- Type
1D array len(nelx*nely)
-
densityfilt
(rmin, filt) Filters with a normalized convolution on the densities with a radius of rmin if:
>>> filt=='density'
- Parameters
rmin (float) – Filter size.
filt (str) – The filter type that is selected, either ‘sensitivity’ or ‘density’.
- Returns
xf – Filterd density distribution.
- Return type
2-D array size(nely, nelx)
-
disp
(x, u, lamba, ke, penal) This function calculates displacement of the objective node and its sensitivity to the densities.
- Parameters
x (2-D array size(nely, nelx)) – Possibly filtered density distribution.
u (1-D array size(max(edof), 1)) – Displacement of all degrees of freedom.
lamba (2-D array size(max(edof), 1)) –
ke (2-D array size(8, 8)) – Element stiffness matrix with full density.
penal (float) – Material model penalisation (SIMP).
- Returns
uout (float) – Displacement objective.
duout (2-D array size(nely, nelx)) – Displacement objective sensitivity to density changes.
-
iter
(penal, rmin, filt) This function performs one iteration of the topology optimisation problem. It
loads the constraints,
calculates the stiffness matrices,
executes the density filter,
executes the FEA solver,
calls upon the displacement objective and its sensitivity calculation,
executes the sensitivity filter,
executes the MMA update scheme,
and finally updates density distribution (design).
- Parameters
penal (float) – Material model penalisation (SIMP).
rmin (float) – Filter size.
filt (str) – The filter type that is selected, either ‘sensitivity’ or ‘density’.
- Returns
change (float) – Largest difference between the new and old density distribution.
uout (float) – Displacement at the objective node for the current design.
-
layout
(penal, rmin, delta, loopy, filt, history=False) Solves the topology optimisation problem by looping over the iter function.
- Parameters
penal (float) – Material model penalisation (SIMP).
rmin (float) – Filter size.
delta (float) – Convergence is roached when delta > change.
loopy (int) – Amount of iteration allowed.
filt (str) – The filter type that is selected, either ‘sensitivity’ or ‘density’.
history (bool, optional) – Do the intermediate results need to be stored.
- Returns
xf (array size(nely, nelx)) – Density distribution resulting from the optimisation.
xf_history (list of arrays len(iterations size(nely, nelx))) – List with the density distributions of all iterations, None when history != True.
-
mma
(m, n, itr, xval, xmin, xmax, xold1, xold2, f0val, df0dx, fval, dfdx, low, upp, a0, a, c, d) This function mmasub performs one MMA-iteration, aimed at solving the nonlinear programming problem:
\[\begin{split}\min & f_0(x) & + a_0z + \sum_{i=1}^m \left(c_iy_i + \frac{1}{2}d_iy_i^2\right) \\ &\text{s.t. }& f_i(x) - a_iz - y_i \leq 0 \hspace{1cm} & i \in \{1,2,\dots,m \} \\ & & x_{\min} \geq x_j \geq x_{\max} & j \in \{1,2,\dots,n \} \\ & & y_i \leq 0 & i \in \{1,2,\dots,m \} \\ & & z \geq 0\end{split}\]- Parameters
m (int) – The number of general constraints.
n (int) – The number of variables \(x_j\).
itr (int) – Current iteration number (=1 the first time mmasub is called).
xval (1-D array len(n)) – Vector with the current values of the variables \(x_j\).
xmin (1-D array len(n)) – Vector with the lower bounds for the variables \(x_j\).
xmax (1-D array len(n)) – Vector with the upper bounds for the variables \(x_j\).
xold1 (1-D array len (n)) – xval, one iteration ago when iter>1, zero othewise.
xold2 (1-D array len (n)) – xval, two iteration ago when iter>2, zero othewise.
f0val (float) – The value of the objective function \(f_0\) at xval.
df0dx (1-D array len(n)) – Vector with the derivatives of the objective function \(f_0\) with respect to the variables \(x_j\), calculated at xval.
fval (1-D array len(m)) – Vector with the values of the constraint functions \(f_i\), calculated at xval.
dfdx (2-D array size(m x n)) – (m x n)-matrix with the derivatives of the constraint functions \(f_i\). with respect to the variables \(x_j\), calculated at xval.
low (1-D array len(n)) – Vector with the lower asymptotes from the previous iteration (provided that iter>1).
upp (1-D array len(n)) – Vector with the upper asymptotes from the previous iteration (provided that iter>1).
a0 (float) – The constants \(a_0\) in the term \(a_0 z\).
a (1-D array len(m)) – Vector with the constants \(a_i1 in the terms :math:\).
c (1-D array len(m)) – Vector with the constants \(c_i\) in the terms \(c_i*y_i\).
d (1-D array len(m)) – Vector with the constants \(d_i\) in the terms \(0.5d_i (y_i)^2\).
- Returns
xmma (1-D array len(n)) – Column vector with the optimal values of the variables \(x_j\) in the current MMA subproblem.
low (1-D array len(n)) – Column vector with the lower asymptotes, calculated and used in the current MMA subproblem.
upp (1-D array len(n)) – Column vector with the upper asymptotes, calculated and used in the current MMA subproblem.
Version September 2007 (and a small change August 2008)
Krister Svanberg <krille@math.kth.se>
Department of Mathematics KTH, SE-10044 Stockholm, Sweden.
Translated to python 3 by A.J.J. Lagerweij TU Delft June 2018
-
sensitivityfilt
(x, rmin, duout, filt) Filters with a normalized convolution on the sensitivity with a radius of rmin if:
>>> filt=='sensitivity'
- Parameters
x (2-D array size(nely, nelx)) – Current density ditribution.
duout (2-D array size(nely, nelx) – Displacement objective sensitivity to density changes.
rmin (float) – Filter size.
filt (str) – The filter type that is selected, either ‘sensitivity’ or ‘density’.
- Returns
duoutf – Filterd sensitivity distribution.
- Return type
2-D array size(nely, nelx)
-
solvemma
(m, n, epsimin, low, upp, alfa, beta, p0, q0, P, Q, a0, a, b, c, d) This function solves the MMA subproblem with a primal-dual Newton method.
\[\begin{split}\min &\sum_{j-1}^n& \left( \frac{p_{0j}^{(k)}}{U_j^{(k)} - x_j} + \frac{q_{0j}^{(k)}}{x_j - L_j^{(k)}} \right) + a_0z + \sum_{i=1}^m \left(c_iy_i + \frac{1}{2}d_iy_i^2\right) \\ &\text{s.t. }& \sum_{j-1}^n \left(\frac{p_{ij}^{(k)}}{U_j^{(k)} - x_j} + \frac{q_{ij}^{(k)}}{x_j - L_j^{(k)}} \right) - a_iz - y_i \leq b_i \\ & & \alpha_j \geq x_j \geq \beta_j\\ & & z \geq 0\end{split}\]- Returns
x – Column vector with the optimal values of the variables x_j in the current MMA subproblem.
- Return type
1-D array len(n)
Plotting Module
Plotting the simulated TopOpt geometry with boundary conditions and loads.
Bram Lagerweij Aerospace Structures and Materials Department TU Delft 2018
-
class
src_Actuator.plotting.
Plot
(load, title=None) This class contains functions that allows the visualisation of the TopOpt algorithm. It can print the density distribution, the boundary conditions and the forces.
- Parameters
load (object, child of the Loads class) – The loadcase(s) considered for this optimisation problem.
title (str) – Title of the plot if required.
-
nelx
Number of elements in x direction.
- Type
int
-
nely
Number of elements in y direction.
- Type
int
-
fig
An empty figure of size nelx/10 and nely/10*1.2 inch.
- Type
matplotlib.pyplot figure
-
ax
The axis system that belongs to fig.
- Type
matplotlib.pyplot axis
-
images
This list contains all density distributions that need to be plotted.
- Type
1-D list with imshow objects
-
add
(x, animated=False) Adding a plot of the density distribution to the figure.
- Parameters
x (2-D array size(nely, nelx)) – The density distribution.
animated (bool, optional) – An animated figure is generated when history = True.
-
boundary
(load) Plotting the boundary conditions.
- Parameters
load (object, child of the Loads class) – The loadcase(s) considered for this optimisation problem.
-
loading
(load) Plotting the loading conditions.
- Parameters
load (object, child of the Loads class) – The loadcase(s) considered for this optimisation problem.
-
save
(filename, fps=10) Saving an plot in svg or mp4 format, depending on the length of the images list. The FasterFFMpegWriter is used when videos are generated. These videos are encoded with a hardware accelerated h264 codec with the .mp4 file format. Other codecs and encoders can be set within the function itself.
- Parameters
filename (str) – Name of the file, excluding the file extension.
fps (int, optional) – Amount of frames per second if the plots are animations.
-
show
() Showing the plot in a window.
-
class
src_Actuator.plotting.
FasterFFMpegWriter
(**kwargs) Bases:
matplotlib.animation.FFMpegWriter
FFMpeg-pipe writer bypassing figure.savefig. To improof saving speed
-
classmethod
bin_path
() Return the binary path to the commandline tool used by a specific subclass. This is a class method so that the tool can be looked for before making a particular MovieWriter subclass available.
-
cleanup
() Clean-up and collect the process used to write the movie file.
-
finish
() Finish any processing for writing the movie.
-
frame_size
A tuple
(width, height)
in pixels of a movie frame.
-
grab_frame
(**savefig_kwargs) Grab the image information from the figure and save as a movie frame.
Doesn’t use savefig to be faster: savefig_kwargs will be ignored.
-
classmethod
isAvailable
() Check to see if a MovieWriter subclass is actually available.
-
saving
(fig, outfile, dpi, *args, **kwargs) Context manager to facilitate writing the movie file.
*args, **kw
are any parameters that should be passed to setup.
-
setup
(fig, outfile, dpi=None) Perform setup for writing the movie file.
- Parameters
fig (~matplotlib.figure.Figure) – The figure object that contains the information for frames
outfile (str) – The filename of the resulting movie file
dpi (int, optional) – The DPI (or resolution) for the file. This controls the size in pixels of the resulting movie file. Default is fig.dpi.
-
classmethod