Continuum Problem Types

Description

TopOpt.jl provides standard continuum topology optimization problem domains for testing and comparing algorithms. This tutorial covers 2D and 3D problems including cantilever beams, MBB beams, L-beams, tie-beams, and INP file import.

Continuum problems model structures as continuous solid domains discretized into finite elements. Design variables are element densities (0 = void, 1 = solid).

Setup

using TopOpt

2D and 3D Point Load Cantilever

The point load cantilever is a standard benchmark: a beam fixed at one end with a point load at the free end.

E = 1.0     # Young's modulus in MPa
ν = 0.3     # Poisson's ratio
f = 1.0     # downward force in N (negative is upward)
nels = (160, 40)        # 2D: 160×40 elements
elsizes = (1.0, 1.0)    # element size in mm
problem_2d = PointLoadCantilever(nels, elsizes, E, ν, f; celltype=:Linear)

The celltype keyword specifies shape function order:

:Linear — bilinear (2D) or trilinear (3D) elements
:Quadratic — biquadratic (2D) or triquadratic (3D) elements

For 3D problems:

nels_3d = (160, 40, 40)     # 3D: 160×40×40 elements
elsizes_3d = (1.0, 1.0, 2.0) # element size in mm
problem_3d = PointLoadCantilever(nels_3d, elsizes_3d, E, ν, f; celltype=:Linear)

2D and 3D Half MBB Beam

The Half MBB (Messerschmitt–Bölkow–Blohm) beam is a simply-supported beam with a central point load. Only half the beam is modeled (symmetry):

nels = (60, 20)
elsizes = (1.0, 1.0)
problem = HalfMBB(nels, elsizes, E, ν, f; celltype=:Quadratic)

Boundary conditions:

Left edge: roller support (vertical displacement fixed)
Bottom right: pin support (both DOFs fixed)
Top center: downward point load

The 3D variant uses 3-tuples for nels and elsizes.

2D L-Beam Problem

The L-beam shows stress concentration effects at the re-entrant corner:

problem = LBeam(
    ; celltype=:Quadratic,
    length=100,
    height=100,
    upperslab=50,
    lowerslab=50,
    E=1.0,
    ν=0.3,
    force=1.0
)

Geometry:

        upperslab   
       ............
       .          .
       .          .
       .          . 
height .          .                     
       .          ......................
       .                               .
       .                               . lowerslab
       .                               .
       .................................
                    length

The load is applied at the midpoint of the “lowerslab” vertical edge.

2D Tie-Beam Problem

The tie-beam has distributed loading on specified elements:

problem = TieBeam(; celltype=:Quadratic)

Fixed supports at both ends
Distributed downward load on specified elements
2D only

Reading INP (Abaqus) Files

For complex geometries, import from .inp (Abaqus) files:

filename = joinpath(@__DIR__, "problem.inp")
problem = InpStiffness(filename)

Workflow:

Define geometry in CAD software (FreeCAD, SolidWorks)
Mesh and export as .inp format
Import into TopOpt.jl for optimization

The .inp file contains nodes, elements, materials, BCs, and loads.

Defining a custom problem type

Beyond the built-in benchmarks, you can define your own continuum problem by subtyping StiffnessTopOptProblem and providing a grid, a ConstraintHandler, and a Metadata. The example below defines a cantilever fixed on the left and loaded at the top-center — a modification of PointLoadCantilever that illustrates each step: build a RectilinearGrid, mark node sets with addnodeset!, create the displacement DofHandler and ConstraintHandler, and assemble the Metadata and load dictionary.

const TTP = TopOpt.TopOptProblems
const Ferrite = TTP.Ferrite

struct CustomCantilever{dim,T,N,M,Tr,Tc,Tf,Tm} <: TTP.StiffnessTopOptProblem{dim,T}
    rect_grid::Tr
    E::T
    ν::T
    ch::Tc
    force::T
    force_dof::Tf
    metadata::Tm
end

function CustomCantilever(nels, sizes, E=1.0, ν=0.3, force=1.0)
    T = float(promote_type(eltype(sizes), typeof(E), typeof(ν), typeof(force)))
    rect_grid = TTP.RectilinearGrid(nels, T.(sizes); celltype=:Linear)

    # Fixed support on the left edge
    haskey(rect_grid.grid.nodesets, "fixed_all") && pop!(rect_grid.grid.nodesets, "fixed_all")
    TTP.addnodeset!(rect_grid.grid, "fixed_all", x -> TTP.left(rect_grid, x))

    # Downward load at the top-center
    haskey(rect_grid.grid.nodesets, "top_load") && pop!(rect_grid.grid.nodesets, "top_load")
    TTP.addnodeset!(rect_grid.grid, "top_load", x -> TTP.top(rect_grid, x) && TTP.middlex(rect_grid, x))

    # Displacement field (2 DOFs/node) and Dirichlet BC
    dh = Ferrite.DofHandler(rect_grid.grid)
    refshape = Ferrite.getrefshape(eltype(rect_grid.grid.cells))
    Ferrite.add!(dh, :u, Ferrite.Lagrange{refshape,1}()^2)
    Ferrite.close!(dh)

    ch = Ferrite.ConstraintHandler(dh)
    Ferrite.add!(
        ch,
        Ferrite.Dirichlet(
            :u, Ferrite.getnodeset(rect_grid.grid, "fixed_all"), (x, t) -> zeros(T, 2), collect(1:2)
        ),
    )
    Ferrite.close!(ch)
    Ferrite.update!(ch, T(0))

    metadata = TTP.Metadata(dh)
    fnode = Tuple(Ferrite.getnodeset(rect_grid.grid, "top_load"))[1]
    force_dof = metadata.node_dofs[2, fnode]

    N = TTP.nnodespercell(rect_grid)
    M = TTP.nfacespercell(rect_grid)
    return CustomCantilever{2,T,N,M,typeof(rect_grid),typeof(ch),typeof(force_dof),typeof(metadata)}(
        rect_grid, E, ν, ch, force, force_dof, metadata
    )
end

# Trait methods the FEA assembly relies on
TTP.nnodespercell(p::CustomCantilever) = TTP.nnodespercell(p.rect_grid)
function TTP.getcloaddict(p::CustomCantilever{2,T}) where {T}
    fnode = Tuple(Ferrite.getnodeset(p.rect_grid.grid, "top_load"))[1]
    return Dict{Int,Vector{T}}(fnode => [0.0, -p.force])
end

The two trait methods are the only glue required beyond the struct itself: the remaining accessors (getdim, floattype, getE, getν, getdh, …) have defaults on StiffnessTopOptProblem. The custom problem is then used exactly like a built-in one:

problem = CustomCantilever((40, 20), (1.0, 1.0), E, ν, f)
solver = FEASolver(DirectSolver, problem; xmin=0.001, penalty=PowerPenaltyFun(3.0))
filter = DensityFilterFun(solver; rmin=2.0)
comp = ComplianceFun(solver)
volfrac = VolumeFun(solver)

obj = x -> comp(filter(PseudoDensities(x)))
constr = x -> volfrac(filter(PseudoDensities(x))) - 0.5
model = Model(obj)
addvar!(model, zeros(getncells(problem)), ones(getncells(problem)))
add_ineq_constraint!(model, constr)
r = optimize(model, MMA87(), fill(0.5, getncells(problem));
    options=MMAOptions(; maxiter=100, tol=Tolerance(; kkt=1e-4, f=1e-4)))
println("custom problem: compliance = $(round(obj(r.minimizer), digits=2)), volume = $(round(constr(r.minimizer) + 0.5, digits=3))")

[ Info:   iter       obj      Δobj  violation  kkt_residual  
[ Info:      0   6.8e+01       Inf   0.0e+00   9.6e+00
[ Info:      1   2.9e+01   3.8e+01   0.0e+00   2.4e+01
[ Info:      2   1.4e+01   1.5e+01   0.0e+00   3.0e+00
[ Info:      3   1.1e+01   3.4e+00   0.0e+00   5.2e-01
[ Info:      4   9.7e+00   8.7e-01   0.0e+00   1.1e-01
[ Info:      5   9.5e+00   2.5e-01   0.0e+00   4.4e-02
[ Info:      6   9.4e+00   1.2e-01   0.0e+00   2.5e-02
[ Info:      7   9.3e+00   6.6e-02   0.0e+00   1.5e-02
[ Info:      8   9.3e+00   3.6e-02   0.0e+00   1.2e-02
[ Info:      9   9.3e+00   1.9e-02   0.0e+00   1.1e-02
[ Info:     10   9.3e+00   8.4e-03   0.0e+00   9.1e-03
[ Info:     11   9.2e+00   5.0e-03   0.0e+00   8.1e-03
[ Info:     12   9.2e+00   4.2e-03   0.0e+00   7.3e-03
[ Info:     13   9.2e+00   3.8e-03   0.0e+00   4.6e-03
[ Info:     14   9.2e+00   3.1e-03   0.0e+00   4.5e-03
[ Info:     15   9.2e+00   1.8e-03   0.0e+00   3.5e-03
[ Info:     16   9.2e+00   1.3e-03   0.0e+00   3.6e-03
[ Info:     17   9.2e+00   8.1e-04   0.0e+00   3.2e-03
[ Info:     18   9.2e+00   3.7e-04   0.0e+00   2.4e-03
[ Info:     19   9.2e+00   3.1e-04   0.0e+00   1.1e-03
[ Info:     20   9.2e+00   1.3e-04   0.0e+00   1.4e-03
[ Info:     21   9.2e+00   9.4e-05   0.0e+00   1.4e-03
[ Info:     22   9.2e+00   1.1e-04   0.0e+00   1.4e-03
[ Info:     23   9.2e+00   1.3e-04   0.0e+00   1.4e-03
[ Info:     24   9.2e+00   1.6e-04   0.0e+00   1.4e-03
[ Info:     25   9.2e+00   1.9e-04   0.0e+00   1.4e-03
[ Info:     26   9.2e+00   2.2e-04   0.0e+00   1.4e-03
[ Info:     27   9.2e+00   2.5e-04   0.0e+00   1.1e-03
[ Info:     28   9.2e+00   1.6e-04   0.0e+00   2.0e-04
[ Info:     29   9.2e+00   2.6e-05   0.0e+00   2.2e-04
[ Info:     30   9.2e+00   8.4e-06   0.0e+00   2.1e-04
[ Info:     31   9.2e+00   9.2e-06   0.0e+00   2.0e-04
[ Info:     32   9.2e+00   9.6e-06   0.0e+00   1.9e-04
[ Info:     33   9.2e+00   9.6e-06   0.0e+00   1.9e-04
[ Info:     34   9.2e+00   9.5e-06   0.0e+00   1.9e-04
[ Info:     35   9.2e+00   9.1e-06   0.0e+00   1.8e-04
[ Info:     36   9.2e+00   8.7e-06   0.0e+00   1.7e-04
[ Info:     37   9.2e+00   8.3e-06   0.0e+00   1.7e-04
[ Info:     38   9.2e+00   7.8e-06   0.0e+00   1.6e-04
[ Info:     39   9.2e+00   6.5e-06   0.0e+00   1.5e-04
[ Info:     40   9.2e+00   4.9e-06   0.0e+00   1.3e-04
[ Info:     41   9.2e+00   3.9e-06   0.0e+00   1.2e-04
[ Info:     42   9.2e+00   3.1e-06   0.0e+00   1.1e-04
[ Info:     43   9.2e+00   2.4e-06   0.0e+00   9.7e-05
custom problem: compliance = 9.23, volume = 0.5