SIMP example: Point Load Cantilever

Tip

This example is also available as a Jupyter notebook: simp.ipynb

Commented Program

What follows is a program spliced with comments. The full program, without comments, can be found in the next section.

using TopOpt

Define the problem

E = 1.0 # Young’s modulus
v = 0.3 # Poisson’s ratio
f = 1.0; # downward force

nels = (30, 10, 10)
problem = PointLoadCantilever(Val{:Linear}, nels, (1.0, 1.0, 1.0), E, v, f);

See also the detailed API of PointLoadCantilever:

TopOpt.TopOptProblems.PointLoadCantilever — Type

///**********************************
///*                                *
///*                                * |
///*                                * |
///********************************** v


@params struct PointLoadCantilever{dim, T, N, M} <: StiffnessTopOptProblem{dim, T}
    rect_grid::RectilinearGrid{dim, T, N, M}
    E::T
    ν::T
    ch::ConstraintHandler{<:DofHandler{dim, <:Cell{dim,N,M}, T}, T}
    force::T
    force_dof::Integer
    black::AbstractVector
    white::AbstractVector
    varind::AbstractVector{Int}
    metadata::Metadata
end

dim: dimension of the problem
T: number type for computations and coordinates
N: number of nodes in a cell of the grid
M: number of faces in a cell of the grid
rect_grid: a RectilinearGrid struct
E: Young's modulus
ν: Poisson's ration
force: force at the center right of the cantilever beam (positive is downward)
force_dof: dof number at which the force is applied
ch: a Ferrite.ConstraintHandler struct
metadata: Metadata having various cell-node-dof relationships
black: a BitVector of length equal to the number of elements where black[e] is 1 iff the e^th element must be part of the final design
white: a BitVector of length equal to the number of elements where white[e] is 1 iff the e^th element must not be part of the final design
varind: an AbstractVector{Int} of length equal to the number of elements where varind[e] gives the index of the decision variable corresponding to element e. Because some elements can be fixed to be black or white, not every element has a decision variable associated.

source

Parameter settings

V = 0.3 # volume fraction
xmin = 1e-6 # minimum density
rmin = 2.0; # density filter radius

Define a finite element solver

penalty = TopOpt.PowerPenalty(3.0)
solver = FEASolver(Direct, problem; xmin=xmin, penalty=penalty)

Define compliance objective

comp = TopOpt.Compliance(solver)
filter = DensityFilter(solver; rmin=rmin)
obj = x -> comp(filter(PseudoDensities(x)))

#1 (generic function with 1 method)

Define volume constraint

volfrac = TopOpt.Volume(solver)
constr = x -> volfrac(filter(PseudoDensities(x))) - V

#3 (generic function with 1 method)

Define subproblem optimizer

x0 = fill(V, length(solver.vars))
model = Model(obj)
addvar!(model, zeros(length(x0)), ones(length(x0)))
add_ineq_constraint!(model, constr)
alg = MMA87()
convcriteria = Nonconvex.KKTCriteria()
options = MMAOptions(;
    maxiter=3000, tol=Nonconvex.Tolerance(; x=1e-3, f=1e-3, kkt=0.001), convcriteria
)
r = optimize(model, alg, x0; options)

@show obj(r.minimizer)

42.23812907385913

(Optional) Visualize the result using Makie.jl

Need to run using Pkg; Pkg.add(["Makie", "GLMakie"]) first

using Makie, GLMakie
using TopOpt.TopOptProblems.Visualization: visualize
fig = visualize(
    problem; topology = r.minimizer,
    default_exagg_scale = 0.07, scale_range = 10.0,
    vector_linewidth = 3, vector_arrowsize = 0.5,
)
Makie.display(fig)

or convert it to a Mesh Need to run using Pkg; Pkg.add(GeometryBasics) first

import Makie, GeometryBasics
result_mesh = GeometryBasics.Mesh(problem, r.minimizer);
Makie.mesh(result_mesh)

Plain Program

Below follows a version of the program without any comments. The file is also available here: simp.jl

using TopOpt

E = 1.0 # Young’s modulus
v = 0.3 # Poisson’s ratio
f = 1.0; # downward force

nels = (30, 10, 10)
problem = PointLoadCantilever(Val{:Linear}, nels, (1.0, 1.0, 1.0), E, v, f);

V = 0.3 # volume fraction
xmin = 1e-6 # minimum density
rmin = 2.0; # density filter radius

penalty = TopOpt.PowerPenalty(3.0)
solver = FEASolver(Direct, problem; xmin=xmin, penalty=penalty)

comp = TopOpt.Compliance(solver)
filter = DensityFilter(solver; rmin=rmin)
obj = x -> comp(filter(PseudoDensities(x)))

volfrac = TopOpt.Volume(solver)
constr = x -> volfrac(filter(PseudoDensities(x))) - V

x0 = fill(V, length(solver.vars))
model = Model(obj)
addvar!(model, zeros(length(x0)), ones(length(x0)))
add_ineq_constraint!(model, constr)
alg = MMA87()
convcriteria = Nonconvex.KKTCriteria()
options = MMAOptions(;
    maxiter=3000, tol=Nonconvex.Tolerance(; x=1e-3, f=1e-3, kkt=0.001), convcriteria
)
r = optimize(model, alg, x0; options)

@show obj(r.minimizer)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

This page was generated using Literate.jl.