Topology optimization
Keywords: Topology optimization, weak and strong form, non-linear problem, Laplacian, grid topology
Figure 1: Optimization of the bending beam. Evolution of the density for fixed total mass.
This example is also available as a Jupyter notebook: topology_optimization.ipynb
.
Introduction
Topology optimization is the task of finding structures that are mechanically ideal. In this example we cover the bending beam, where we specify a load, boundary conditions and the total mass. Then, our objective is to find the most suitable geometry within the design space minimizing the compliance (i.e. the inverse stiffness) of the structure. We shortly introduce our simplified model for regular meshes. A detailed derivation of the method and advanced techniques can be found in [14] and [15].
We start by introducing the local, elementwise density $\chi \in [\chi_{\text{min}}, 1]$ of the material, where we choose $\chi_{\text{min}}$ slightly above zero to prevent numerical instabilities. Here, $\chi = \chi_{\text{min}}$ means void and $\chi=1$ means bulk material. Then, we use a SIMP ansatz (solid isotropic material with penalization) for the stiffness tensor $C(\chi) = \chi^p C_0$, where $C_0$ is the stiffness of the bulk material. The SIMP exponent $p>1$ ensures that the model prefers the density values void and bulk before the intermediate values. The variational formulation then yields the modified Gibbs energy
\[G = \int_{\Omega} \frac{1}{2} \chi^p \varepsilon : C : \varepsilon \; \text{d}V - \int_{\Omega} \boldsymbol{f} \cdot \boldsymbol{u} \; \text{d}V - \int_{\partial\Omega} \boldsymbol{t} \cdot \boldsymbol{u} \; \text{d}A.\]
Furthermore, we receive the evolution equation of the density and the additional Neumann boundary condition in the strong form
\[p_\chi + \eta \dot{\chi} + \lambda + \gamma - \beta \nabla^2 \chi \ni 0 \quad \forall \textbf{x} \in \Omega,\]
\[\beta \nabla \chi \cdot \textbf{n} = 0 \quad \forall \textbf{x} \in \partial \Omega,\]
with the thermodynamic driving force
\[p_\chi = \frac{1}{2} p \chi^{p-1} \varepsilon : C : \varepsilon.\]
We obtain the mechanical displacement field by applying the Finite Element Method to the weak form of the Gibbs energy using Ferrite. In contrast, we use the evolution equation (i.e. the strong form) to calculate the value of the density field $\chi$. The advantage of this "split" approach is the very high computation speed. The evolution equation consists of the driving force, the damping parameter $\eta$, the regularization parameter $\beta$ times the Laplacian, which is necessary to avoid numerical issues like mesh dependence or checkerboarding, and the constraint parameters $\lambda$, to keep the mass constant, and $\gamma$, to avoid leaving the set $[\chi_{\text{min}}, 1]$. By including gradient regularization, it becomes necessary to calculate the Laplacian. The Finite Difference Method for square meshes with the edge length $\Delta h$ approximates the Laplacian as follows:
\[\nabla^2 \chi_p = \frac{1}{(\Delta h)^2} (\chi_n + \chi_s + \chi_w + \chi_e - 4 \chi_p)\]
Here, the indices refer to the different cardinal directions. Boundary element do not have neighbors in each direction. However, we can calculate the central difference to fulfill Neumann boundary condition. For example, if the element is on the left boundary, we have to fulfill
\[\nabla \chi_p \cdot \textbf{n} = \frac{1}{\Delta h} (\chi_w - \chi_e) = 0\]
from which follows $\chi_w = \chi_e$. Thus for boundary elements we can replace the value for the missing neighbor by the value of the opposite neighbor. In order to find the corresponding neighbor elements, we will make use of Ferrites grid topology funcionalities.
Commented Program
We now solve the problem in Ferrite. What follows is a program spliced with comments. The full program, without comments, can be found in the next section.
First we load all necessary packages.
using Ferrite, SparseArrays, LinearAlgebra, Tensors, Printf
Next, we create a simple square grid of the size 2x1. We apply a fixed Dirichlet boundary condition to the left facet set, called clamped
. On the right facet, we create a small set traction
, where we will later apply a force in negative y-direction.
function create_grid(n)
corners = [
Vec{2}((0.0, 0.0)),
Vec{2}((2.0, 0.0)),
Vec{2}((2.0, 1.0)),
Vec{2}((0.0, 1.0)),
]
grid = generate_grid(Quadrilateral, (2 * n, n), corners)
# node-/facesets for boundary conditions
addnodeset!(grid, "clamped", x -> x[1] ≈ 0.0)
addfacetset!(grid, "traction", x -> x[1] ≈ 2.0 && norm(x[2] - 0.5) <= 0.05)
return grid
end
Next, we create the FE values, the DofHandler and the Dirichlet boundary condition.
function create_values()
# quadrature rules
qr = QuadratureRule{RefQuadrilateral}(2)
facet_qr = FacetQuadratureRule{RefQuadrilateral}(2)
# cell and facetvalues for u
ip = Lagrange{RefQuadrilateral, 1}()^2
cellvalues = CellValues(qr, ip)
facetvalues = FacetValues(facet_qr, ip)
return cellvalues, facetvalues
end
function create_dofhandler(grid)
dh = DofHandler(grid)
add!(dh, :u, Lagrange{RefQuadrilateral, 1}()^2) # displacement
close!(dh)
return dh
end
function create_bc(dh)
dbc = ConstraintHandler(dh)
add!(dbc, Dirichlet(:u, getnodeset(dh.grid, "clamped"), (x, t) -> zero(Vec{2}), [1, 2]))
close!(dbc)
t = 0.0
update!(dbc, t)
return dbc
end
Now, we define a struct to store all necessary material parameters (stiffness tensor of the bulk material and the parameters for topology optimization) and add a constructor to the struct to initialize it by using the common material parameters Young's modulus and Poisson number.
struct MaterialParameters{T, S <: SymmetricTensor{4, 2, T}}
C::S
χ_min::T
p::T
β::T
η::T
end
function MaterialParameters(E, ν, χ_min, p, β, η)
δ(i, j) = i == j ? 1.0 : 0.0 # helper function
G = E / 2(1 + ν) # =μ
λ = E * ν / (1 - ν^2) # correction for plane stress included
C = SymmetricTensor{4, 2}((i, j, k, l) -> λ * δ(i, j) * δ(k, l) + G * (δ(i, k) * δ(j, l) + δ(i, l) * δ(j, k)))
return MaterialParameters(C, χ_min, p, β, η)
end
To store the density and the strain required to calculate the driving forces, we create the struct MaterialState
. We add a constructor to initialize the struct. The function update_material_states!
updates the density values once we calculated the new values.
mutable struct MaterialState{T, S <: AbstractArray{SymmetricTensor{2, 2, T, 3}, 1}}
χ::T # density
ε::S # strain in each quadrature point
end
function MaterialState(ρ, n_qp)
return MaterialState(ρ, Array{SymmetricTensor{2, 2, Float64, 3}, 1}(undef, n_qp))
end
function update_material_states!(χn1, states, dh)
for (element, state) in zip(CellIterator(dh), states)
state.χ = χn1[cellid(element)]
end
return
end
Next, we define a function to calculate the driving forces for all elements. For this purpose, we iterate through all elements and calculate the average strain in each element. Then, we compute the driving force from the formula introduced at the beginning. We create a second function to collect the density in each element.
function compute_driving_forces(states, mp, dh, χn)
pΨ = zeros(length(states))
for (element, state) in zip(CellIterator(dh), states)
i = cellid(element)
ε = sum(state.ε) / length(state.ε) # average element strain
pΨ[i] = 1 / 2 * mp.p * χn[i]^(mp.p - 1) * (ε ⊡ mp.C ⊡ ε)
end
return pΨ
end
function compute_densities(states, dh)
χn = zeros(length(states))
for (element, state) in zip(CellIterator(dh), states)
i = cellid(element)
χn[i] = state.χ
end
return χn
end
For the Laplacian we need some neighboorhood information which is constant throughout the analysis so we compute it once and cache it. We iterate through each facet of each element, obtaining the neighboring element by using the getneighborhood
function. For boundary facets, the function call will return an empty object. In that case we use the dictionary to instead find the opposite facet, as discussed in the introduction.
function cache_neighborhood(dh, topology)
nbgs = Vector{Vector{Int}}(undef, getncells(dh.grid))
_nfacets = nfacets(dh.grid.cells[1])
opp = Dict(1 => 3, 2 => 4, 3 => 1, 4 => 2)
for element in CellIterator(dh)
nbg = zeros(Int, _nfacets)
i = cellid(element)
for j in 1:_nfacets
nbg_cellid = getneighborhood(topology, dh.grid, FacetIndex(i, j))
if !isempty(nbg_cellid)
nbg[j] = first(nbg_cellid)[1] # assuming only one facet neighbor per cell
else # boundary facet
nbg[j] = first(getneighborhood(topology, dh.grid, FacetIndex(i, opp[j])))[1]
end
end
nbgs[i] = nbg
end
return nbgs
end
Now we calculate the Laplacian using the previously cached neighboorhood information.
function approximate_laplacian(nbgs, χn, Δh)
∇²χ = zeros(length(nbgs))
for i in 1:length(nbgs)
nbg = nbgs[i]
∇²χ[i] = (χn[nbg[1]] + χn[nbg[2]] + χn[nbg[3]] + χn[nbg[4]] - 4 * χn[i]) / (Δh^2)
end
return ∇²χ
end
For the iterative computation of the solution, a function is needed to update the densities in each element. To ensure that the mass is kept constant, we have to calculate the constraint parameter $\lambda$, which we do via the bisection method. We repeat the calculation until the difference between the average density (calculated from the element-wise trial densities) and the target density nearly vanishes. By using the extremal values of $\Delta \chi$ as the starting interval, we guarantee that the method converges eventually.
function compute_χn1(χn, Δχ, ρ, ηs, χ_min)
n_el = length(χn)
χ_trial = zeros(n_el)
ρ_trial = 0.0
λ_lower = minimum(Δχ) - ηs
λ_upper = maximum(Δχ) + ηs
λ_trial = 0.0
while abs(ρ - ρ_trial) > 1.0e-7
for i in 1:n_el
Δχt = 1 / ηs * (Δχ[i] - λ_trial)
χ_trial[i] = max(χ_min, min(1.0, χn[i] + Δχt))
end
ρ_trial = 0.0
for i in 1:n_el
ρ_trial += χ_trial[i] / n_el
end
if ρ_trial > ρ
λ_lower = λ_trial
elseif ρ_trial < ρ
λ_upper = λ_trial
end
λ_trial = 1 / 2 * (λ_upper + λ_lower)
end
return χ_trial
end
Lastly, we use the following helper function to compute the average driving force, which is later used to normalize the driving forces. This makes the used material parameters and numerical parameters independent of the problem.
function compute_average_driving_force(mp, pΨ, χn)
n = length(pΨ)
w = zeros(n)
for i in 1:n
w[i] = (χn[i] - mp.χ_min) * (1 - χn[i])
end
p_Ω = sum(w .* pΨ) / sum(w) # average driving force
return p_Ω
end
Finally, we put everything together to update the density. The loop ensures the stability of the updated solution.
function update_density(dh, states, mp, ρ, neighboorhoods, Δh)
n_j = Int(ceil(6 * mp.β / (mp.η * Δh^2))) # iterations needed for stability
χn = compute_densities(states, dh) # old density field
χn1 = zeros(length(χn))
for j in 1:n_j
∇²χ = approximate_laplacian(neighboorhoods, χn, Δh) # Laplacian
pΨ = compute_driving_forces(states, mp, dh, χn) # driving forces
p_Ω = compute_average_driving_force(mp, pΨ, χn) # average driving force
Δχ = pΨ / p_Ω + mp.β * ∇²χ
χn1 = compute_χn1(χn, Δχ, ρ, mp.η, mp.χ_min)
if j < n_j
χn[:] = χn1[:]
end
end
return χn1
end
Now, we move on to the Finite Element part of the program. We use the following function to assemble our linear system.
function doassemble!(cellvalues::CellValues, facetvalues::FacetValues, K::SparseMatrixCSC, grid::Grid, dh::DofHandler, mp::MaterialParameters, u, states)
r = zeros(ndofs(dh))
assembler = start_assemble(K, r)
nu = getnbasefunctions(cellvalues)
re = zeros(nu) # local residual vector
Ke = zeros(nu, nu) # local stiffness matrix
for (element, state) in zip(CellIterator(dh), states)
fill!(Ke, 0)
fill!(re, 0)
eldofs = celldofs(element)
ue = u[eldofs]
elmt!(Ke, re, element, cellvalues, facetvalues, grid, mp, ue, state)
assemble!(assembler, celldofs(element), Ke, re)
end
return K, r
end
The element routine is used to calculate the elementwise stiffness matrix and the residual. In contrast to a purely elastomechanic problem, for topology optimization we additionally use our material state to receive the density value of the element and to store the strain at each quadrature point.
function elmt!(Ke, re, element, cellvalues, facetvalues, grid, mp, ue, state)
n_basefuncs = getnbasefunctions(cellvalues)
reinit!(cellvalues, element)
χ = state.χ
# We only assemble lower half triangle of the stiffness matrix and then symmetrize it.
@inbounds for q_point in 1:getnquadpoints(cellvalues)
dΩ = getdetJdV(cellvalues, q_point)
state.ε[q_point] = function_symmetric_gradient(cellvalues, q_point, ue)
for i in 1:n_basefuncs
δεi = shape_symmetric_gradient(cellvalues, q_point, i)
for j in 1:i
δεj = shape_symmetric_gradient(cellvalues, q_point, j)
Ke[i, j] += (χ)^(mp.p) * (δεi ⊡ mp.C ⊡ δεj) * dΩ
end
re[i] += (-δεi ⊡ ((χ)^(mp.p) * mp.C ⊡ state.ε[q_point])) * dΩ
end
end
symmetrize_lower!(Ke)
@inbounds for facet in 1:nfacets(getcells(grid, cellid(element)))
if (cellid(element), facet) ∈ getfacetset(grid, "traction")
reinit!(facetvalues, element, facet)
t = Vec((0.0, -1.0)) # force pointing downwards
for q_point in 1:getnquadpoints(facetvalues)
dΓ = getdetJdV(facetvalues, q_point)
for i in 1:n_basefuncs
δu = shape_value(facetvalues, q_point, i)
re[i] += (δu ⋅ t) * dΓ
end
end
end
end
return
end
function symmetrize_lower!(K)
for i in 1:size(K, 1)
for j in (i + 1):size(K, 1)
K[i, j] = K[j, i]
end
end
return
end
We put everything together in the main function. Here the user may choose the radius parameter, which is related to the regularization parameter as $\beta = ra^2$, the starting density, the number of elements in vertical direction and finally the name of the output. Additionally, the user may choose whether only the final design (default) or every iteration step is saved.
First, we compute the material parameters and create the grid, DofHandler, boundary condition and FE values. Then we prepare the iterative Newton-Raphson method by pre-allocating all important vectors. Furthermore, we create material states for each element and construct the topology of the grid.
During each iteration step, first we solve our FE problem in the Newton-Raphson loop. With the solution of the elastomechanic problem, we check for convergence of our topology design. The criteria has to be fulfilled twice in a row to avoid oscillations. If no convergence is reached yet, we update our design and prepare the next iteration step. Finally, we output the results in paraview and calculate the relative stiffness of the final design, i.e. how much how the stiffness increased compared to the starting point.
function topopt(ra, ρ, n, filename; output = :false)
# material
mp = MaterialParameters(210.0e3, 0.3, 1.0e-3, 3.0, ra^2, 15.0)
# grid, dofhandler, boundary condition
grid = create_grid(n)
dh = create_dofhandler(grid)
Δh = 1 / n # element edge length
dbc = create_bc(dh)
# cellvalues
cellvalues, facetvalues = create_values()
# Pre-allocate solution vectors, etc.
n_dofs = ndofs(dh) # total number of dofs
u = zeros(n_dofs) # solution vector
un = zeros(n_dofs) # previous solution vector
Δu = zeros(n_dofs) # previous displacement correction
ΔΔu = zeros(n_dofs) # new displacement correction
# create material states
states = [MaterialState(ρ, getnquadpoints(cellvalues)) for _ in 1:getncells(dh.grid)]
χ = zeros(getncells(dh.grid))
r = zeros(n_dofs) # residual
K = allocate_matrix(dh) # stiffness matrix
i_max = 300 ## maximum number of iteration steps
tol = 1.0e-4
compliance = 0.0
compliance_0 = 0.0
compliance_n = 0.0
conv = :false
topology = ExclusiveTopology(grid)
neighboorhoods = cache_neighborhood(dh, topology)
# Newton-Raphson loop
NEWTON_TOL = 1.0e-8
print("\n Starting Newton iterations\n")
for it in 1:i_max
apply_zero!(u, dbc)
newton_itr = -1
while true
newton_itr += 1
if newton_itr > 10
error("Reached maximum Newton iterations, aborting")
break
end
# current guess
u .= un .+ Δu
K, r = doassemble!(cellvalues, facetvalues, K, grid, dh, mp, u, states)
norm_r = norm(r[Ferrite.free_dofs(dbc)])
if (norm_r) < NEWTON_TOL
break
end
apply_zero!(K, r, dbc)
ΔΔu = Symmetric(K) \ r
apply_zero!(ΔΔu, dbc)
Δu .+= ΔΔu
end # of loop while NR-Iteration
# calculate compliance
compliance = 1 / 2 * u' * K * u
if it == 1
compliance_0 = compliance
end
# check convergence criterium (twice!)
if abs(compliance - compliance_n) / compliance < tol
if conv
println("Converged at iteration number: ", it)
break
else
conv = :true
end
else
conv = :false
end
# update density
χ = update_density(dh, states, mp, ρ, neighboorhoods, Δh)
# update old displacement, density and compliance
un .= u
Δu .= 0.0
update_material_states!(χ, states, dh)
compliance_n = compliance
# output during calculation
if output
i = @sprintf("%3.3i", it)
filename_it = string(filename, "_", i)
VTKGridFile(filename_it, grid) do vtk
write_cell_data(vtk, χ, "density")
end
end
end
# export converged results
if !output
VTKGridFile(filename, grid) do vtk
write_cell_data(vtk, χ, "density")
end
end
@printf "Rel. stiffness: %.4f \n" compliance^(-1) / compliance_0^(-1)
return
end
Lastly, we call our main function and compare the results. To create the complete output with all iteration steps, it is possible to set the output parameter to true
.
grid, χ =topopt(0.02, 0.5, 60, "small_radius"; output=:false);
@time topopt(0.03, 0.5, 60, "large_radius"; output = :false);
#topopt(0.02, 0.5, 60, "topopt_animation"; output=:true); # can be used to create animations
Starting Newton iterations
Converged at iteration number: 65
Rel. stiffness: 4.8466
4.338780 seconds (2.24 M allocations: 2.005 GiB, 2.39% gc time, 1.87% compilation time)
We observe, that the stiffness for the lower value of $ra$ is higher, but also requires more iterations until convergence and finer structures to be manufactured, as can be seen in Figure 2:
Figure 2: Optimization results of the bending beam for smaller (left) and larger (right) value of the regularization parameter $\beta$.
To prove mesh independence, the user could vary the mesh resolution and compare the results.
References
- [14]
- D. R. Jantos, K. Hackl and P. Junker. An accurate and fast regularization approach to thermodynamic topology optimization. International Journal for Numerical Methods in Engineering 117, 991–1017 (2019).
- [15]
- M. Blaszczyk, D. R. Jantos and P. Junker. Application of Taylor series combined with the weighted least square method to thermodynamic topology optimization. Computer Methods in Applied Mechanics and Engineering 393, 114698 (2022).
Plain program
Here follows a version of the program without any comments. The file is also available here: topology_optimization.jl
.
using Ferrite, SparseArrays, LinearAlgebra, Tensors, Printf
function create_grid(n)
corners = [
Vec{2}((0.0, 0.0)),
Vec{2}((2.0, 0.0)),
Vec{2}((2.0, 1.0)),
Vec{2}((0.0, 1.0)),
]
grid = generate_grid(Quadrilateral, (2 * n, n), corners)
# node-/facesets for boundary conditions
addnodeset!(grid, "clamped", x -> x[1] ≈ 0.0)
addfacetset!(grid, "traction", x -> x[1] ≈ 2.0 && norm(x[2] - 0.5) <= 0.05)
return grid
end
function create_values()
# quadrature rules
qr = QuadratureRule{RefQuadrilateral}(2)
facet_qr = FacetQuadratureRule{RefQuadrilateral}(2)
# cell and facetvalues for u
ip = Lagrange{RefQuadrilateral, 1}()^2
cellvalues = CellValues(qr, ip)
facetvalues = FacetValues(facet_qr, ip)
return cellvalues, facetvalues
end
function create_dofhandler(grid)
dh = DofHandler(grid)
add!(dh, :u, Lagrange{RefQuadrilateral, 1}()^2) # displacement
close!(dh)
return dh
end
function create_bc(dh)
dbc = ConstraintHandler(dh)
add!(dbc, Dirichlet(:u, getnodeset(dh.grid, "clamped"), (x, t) -> zero(Vec{2}), [1, 2]))
close!(dbc)
t = 0.0
update!(dbc, t)
return dbc
end
struct MaterialParameters{T, S <: SymmetricTensor{4, 2, T}}
C::S
χ_min::T
p::T
β::T
η::T
end
function MaterialParameters(E, ν, χ_min, p, β, η)
δ(i, j) = i == j ? 1.0 : 0.0 # helper function
G = E / 2(1 + ν) # =μ
λ = E * ν / (1 - ν^2) # correction for plane stress included
C = SymmetricTensor{4, 2}((i, j, k, l) -> λ * δ(i, j) * δ(k, l) + G * (δ(i, k) * δ(j, l) + δ(i, l) * δ(j, k)))
return MaterialParameters(C, χ_min, p, β, η)
end
mutable struct MaterialState{T, S <: AbstractArray{SymmetricTensor{2, 2, T, 3}, 1}}
χ::T # density
ε::S # strain in each quadrature point
end
function MaterialState(ρ, n_qp)
return MaterialState(ρ, Array{SymmetricTensor{2, 2, Float64, 3}, 1}(undef, n_qp))
end
function update_material_states!(χn1, states, dh)
for (element, state) in zip(CellIterator(dh), states)
state.χ = χn1[cellid(element)]
end
return
end
function compute_driving_forces(states, mp, dh, χn)
pΨ = zeros(length(states))
for (element, state) in zip(CellIterator(dh), states)
i = cellid(element)
ε = sum(state.ε) / length(state.ε) # average element strain
pΨ[i] = 1 / 2 * mp.p * χn[i]^(mp.p - 1) * (ε ⊡ mp.C ⊡ ε)
end
return pΨ
end
function compute_densities(states, dh)
χn = zeros(length(states))
for (element, state) in zip(CellIterator(dh), states)
i = cellid(element)
χn[i] = state.χ
end
return χn
end
function cache_neighborhood(dh, topology)
nbgs = Vector{Vector{Int}}(undef, getncells(dh.grid))
_nfacets = nfacets(dh.grid.cells[1])
opp = Dict(1 => 3, 2 => 4, 3 => 1, 4 => 2)
for element in CellIterator(dh)
nbg = zeros(Int, _nfacets)
i = cellid(element)
for j in 1:_nfacets
nbg_cellid = getneighborhood(topology, dh.grid, FacetIndex(i, j))
if !isempty(nbg_cellid)
nbg[j] = first(nbg_cellid)[1] # assuming only one facet neighbor per cell
else # boundary facet
nbg[j] = first(getneighborhood(topology, dh.grid, FacetIndex(i, opp[j])))[1]
end
end
nbgs[i] = nbg
end
return nbgs
end
function approximate_laplacian(nbgs, χn, Δh)
∇²χ = zeros(length(nbgs))
for i in 1:length(nbgs)
nbg = nbgs[i]
∇²χ[i] = (χn[nbg[1]] + χn[nbg[2]] + χn[nbg[3]] + χn[nbg[4]] - 4 * χn[i]) / (Δh^2)
end
return ∇²χ
end
function compute_χn1(χn, Δχ, ρ, ηs, χ_min)
n_el = length(χn)
χ_trial = zeros(n_el)
ρ_trial = 0.0
λ_lower = minimum(Δχ) - ηs
λ_upper = maximum(Δχ) + ηs
λ_trial = 0.0
while abs(ρ - ρ_trial) > 1.0e-7
for i in 1:n_el
Δχt = 1 / ηs * (Δχ[i] - λ_trial)
χ_trial[i] = max(χ_min, min(1.0, χn[i] + Δχt))
end
ρ_trial = 0.0
for i in 1:n_el
ρ_trial += χ_trial[i] / n_el
end
if ρ_trial > ρ
λ_lower = λ_trial
elseif ρ_trial < ρ
λ_upper = λ_trial
end
λ_trial = 1 / 2 * (λ_upper + λ_lower)
end
return χ_trial
end
function compute_average_driving_force(mp, pΨ, χn)
n = length(pΨ)
w = zeros(n)
for i in 1:n
w[i] = (χn[i] - mp.χ_min) * (1 - χn[i])
end
p_Ω = sum(w .* pΨ) / sum(w) # average driving force
return p_Ω
end
function update_density(dh, states, mp, ρ, neighboorhoods, Δh)
n_j = Int(ceil(6 * mp.β / (mp.η * Δh^2))) # iterations needed for stability
χn = compute_densities(states, dh) # old density field
χn1 = zeros(length(χn))
for j in 1:n_j
∇²χ = approximate_laplacian(neighboorhoods, χn, Δh) # Laplacian
pΨ = compute_driving_forces(states, mp, dh, χn) # driving forces
p_Ω = compute_average_driving_force(mp, pΨ, χn) # average driving force
Δχ = pΨ / p_Ω + mp.β * ∇²χ
χn1 = compute_χn1(χn, Δχ, ρ, mp.η, mp.χ_min)
if j < n_j
χn[:] = χn1[:]
end
end
return χn1
end
function doassemble!(cellvalues::CellValues, facetvalues::FacetValues, K::SparseMatrixCSC, grid::Grid, dh::DofHandler, mp::MaterialParameters, u, states)
r = zeros(ndofs(dh))
assembler = start_assemble(K, r)
nu = getnbasefunctions(cellvalues)
re = zeros(nu) # local residual vector
Ke = zeros(nu, nu) # local stiffness matrix
for (element, state) in zip(CellIterator(dh), states)
fill!(Ke, 0)
fill!(re, 0)
eldofs = celldofs(element)
ue = u[eldofs]
elmt!(Ke, re, element, cellvalues, facetvalues, grid, mp, ue, state)
assemble!(assembler, celldofs(element), Ke, re)
end
return K, r
end
function elmt!(Ke, re, element, cellvalues, facetvalues, grid, mp, ue, state)
n_basefuncs = getnbasefunctions(cellvalues)
reinit!(cellvalues, element)
χ = state.χ
# We only assemble lower half triangle of the stiffness matrix and then symmetrize it.
@inbounds for q_point in 1:getnquadpoints(cellvalues)
dΩ = getdetJdV(cellvalues, q_point)
state.ε[q_point] = function_symmetric_gradient(cellvalues, q_point, ue)
for i in 1:n_basefuncs
δεi = shape_symmetric_gradient(cellvalues, q_point, i)
for j in 1:i
δεj = shape_symmetric_gradient(cellvalues, q_point, j)
Ke[i, j] += (χ)^(mp.p) * (δεi ⊡ mp.C ⊡ δεj) * dΩ
end
re[i] += (-δεi ⊡ ((χ)^(mp.p) * mp.C ⊡ state.ε[q_point])) * dΩ
end
end
symmetrize_lower!(Ke)
@inbounds for facet in 1:nfacets(getcells(grid, cellid(element)))
if (cellid(element), facet) ∈ getfacetset(grid, "traction")
reinit!(facetvalues, element, facet)
t = Vec((0.0, -1.0)) # force pointing downwards
for q_point in 1:getnquadpoints(facetvalues)
dΓ = getdetJdV(facetvalues, q_point)
for i in 1:n_basefuncs
δu = shape_value(facetvalues, q_point, i)
re[i] += (δu ⋅ t) * dΓ
end
end
end
end
return
end
function symmetrize_lower!(K)
for i in 1:size(K, 1)
for j in (i + 1):size(K, 1)
K[i, j] = K[j, i]
end
end
return
end
function topopt(ra, ρ, n, filename; output = :false)
# material
mp = MaterialParameters(210.0e3, 0.3, 1.0e-3, 3.0, ra^2, 15.0)
# grid, dofhandler, boundary condition
grid = create_grid(n)
dh = create_dofhandler(grid)
Δh = 1 / n # element edge length
dbc = create_bc(dh)
# cellvalues
cellvalues, facetvalues = create_values()
# Pre-allocate solution vectors, etc.
n_dofs = ndofs(dh) # total number of dofs
u = zeros(n_dofs) # solution vector
un = zeros(n_dofs) # previous solution vector
Δu = zeros(n_dofs) # previous displacement correction
ΔΔu = zeros(n_dofs) # new displacement correction
# create material states
states = [MaterialState(ρ, getnquadpoints(cellvalues)) for _ in 1:getncells(dh.grid)]
χ = zeros(getncells(dh.grid))
r = zeros(n_dofs) # residual
K = allocate_matrix(dh) # stiffness matrix
i_max = 300 ## maximum number of iteration steps
tol = 1.0e-4
compliance = 0.0
compliance_0 = 0.0
compliance_n = 0.0
conv = :false
topology = ExclusiveTopology(grid)
neighboorhoods = cache_neighborhood(dh, topology)
# Newton-Raphson loop
NEWTON_TOL = 1.0e-8
print("\n Starting Newton iterations\n")
for it in 1:i_max
apply_zero!(u, dbc)
newton_itr = -1
while true
newton_itr += 1
if newton_itr > 10
error("Reached maximum Newton iterations, aborting")
break
end
# current guess
u .= un .+ Δu
K, r = doassemble!(cellvalues, facetvalues, K, grid, dh, mp, u, states)
norm_r = norm(r[Ferrite.free_dofs(dbc)])
if (norm_r) < NEWTON_TOL
break
end
apply_zero!(K, r, dbc)
ΔΔu = Symmetric(K) \ r
apply_zero!(ΔΔu, dbc)
Δu .+= ΔΔu
end # of loop while NR-Iteration
# calculate compliance
compliance = 1 / 2 * u' * K * u
if it == 1
compliance_0 = compliance
end
# check convergence criterium (twice!)
if abs(compliance - compliance_n) / compliance < tol
if conv
println("Converged at iteration number: ", it)
break
else
conv = :true
end
else
conv = :false
end
# update density
χ = update_density(dh, states, mp, ρ, neighboorhoods, Δh)
# update old displacement, density and compliance
un .= u
Δu .= 0.0
update_material_states!(χ, states, dh)
compliance_n = compliance
# output during calculation
if output
i = @sprintf("%3.3i", it)
filename_it = string(filename, "_", i)
VTKGridFile(filename_it, grid) do vtk
write_cell_data(vtk, χ, "density")
end
end
end
# export converged results
if !output
VTKGridFile(filename, grid) do vtk
write_cell_data(vtk, χ, "density")
end
end
@printf "Rel. stiffness: %.4f \n" compliance^(-1) / compliance_0^(-1)
return
end
@time topopt(0.03, 0.5, 60, "large_radius"; output = :false);
#topopt(0.02, 0.5, 60, "topopt_animation"; output=:true); # can be used to create animations
This page was generated using Literate.jl.