RasmahRasmah

Level-Set Topology Optimization

Overview

The density method (SIMP) in the topology-optimization tutorial represents the design as a per-element density field. The level-set method instead represents the design boundary implicitly: a scalar field φ is negative inside the material and positive outside, and the shape is the zero level set {φ = 0}. Moving φ moves the boundary, giving crisp (not gray) interfaces and smooth, mesh-independent shapes.

Rasmah provides the level-set machinery as a self-contained toolkit — a LevelSetGrid, signed-distance reinitialization, Hamilton–Jacobi advection, velocity extension, a radial-basis-function (RBF) parametrization, and the levelset_compliance_optimize! driver.

The level-set grid and a field

LevelSetGrid(nx, ny, dx, dy) is a uniform 2D grid; x_coords / y_coords return the node coordinates. Start from the signed-distance function of a circle:

using Rasmah

g = LevelSetGrid(41, 41, 0.025, 0.025)

φ0 = [sqrt((x_coords(g)[i] - 0.5)^2 + (y_coords(g)[j] - 0.5)^2) - 0.3
      for i in 1:g.nx, j in 1:g.ny]

println("grid: ", g.nx, "×", g.ny, ", φ at centre = ", round(φ0[21, 21]; digits = 3))
grid: 41×41, φ at centre = -0.3

Reinitialization

Advecting a level set distorts its gradient, so the field must be periodically reinitialized to a signed distance (|∇φ| ≈ 1). reinitialize does this with a fast marching-style iteration:

φr = reinitialize(φ0, g; iters = 50)
println("reinitialized field finite: ", all(isfinite, φr))
println("signed distance at (0.5, 0.5): ", round(φr[21, 21]; digits = 3))
reinitialized field finite: true
signed distance at (0.5, 0.5): -0.282

Hamilton–Jacobi advection

solve_hj(φ, v, grid, dt) evolves the interface under a normal velocity field — the level-set analogue of ∂φ/∂t + v|∇φ| = 0. A positive inward velocity expands the material region:

φmoved = solve_hj(φr, fill(1.0, g.nx, g.ny), g, 0.005; steps = 10)

# the rightmost zero-crossing of φ < 0 moves outward
r0 = first([x_coords(g)[i] for i in 1:g.nx if x_coords(g)[i] >= 0.5 && φr[i, 21] >= 0])
r1 = first([x_coords(g)[i] for i in 1:g.nx if x_coords(g)[i] >= 0.5 && φmoved[i, 21] >= 0])
println("interface moved from x = ", round(r0; digits = 3), " to x = ", round(r1; digits = 3))
interface moved from x = 0.8 to x = 0.85

extend_velocity! propagates an interface velocity outward along normals, so the front moves smoothly rather than from a single layer of cells.

An RBF level set

RBFLevelSet(centers, α; kernel) parametrizes the level set by a small set of coefficients α on scattered centers — a low-dimensional design space that is automatically differentiable (the shape is a smooth function of α). The default Wendland kernel is compactly supported, so each coefficient only deforms the interface within a radius c of its centre — a local control that stays stable as the number of centers grows:

centers = [0.2 0.8 0.5; 0.5 0.5 0.8]
rbf = RBFLevelSet(centers, [1.0, -0.5, 0.3]; kernel = :wendland, c = 2.0)
println("φ at centre (0.5, 0.5) = ", round(rbf([0.5, 0.5]); digits = 4))
φ at centre (0.5, 0.5) = 0.6682

Level-set compliance minimization

levelset_compliance_optimize! drives the level set to minimize compliance subject to a volume-fraction constraint, using a reaction–diffusion update on the grid. Removing material to meet the volume constraint raises the compliance — that is the expected trade-off, and the optimizer keeps it as low as possible for the given fraction. The grid discretizes to a TriangleMesh via grid_triangle_mesh:

E, ν = 1.0, 0.3
n = 12
gg = LevelSetGrid(n, n, 1.0 / (n - 1), 1.0 / (n - 1))
m = grid_triangle_mesh(gg)
N = size(m.vertices, 2)

fixed = Int[]; loads = zeros(2N)
for i in 1:N
    x = m.vertices[:, i]
    x[1] < 1e-12 && append!(fixed, 2i - 1, 2i)
    x[1] > 1 - 1e-12 && (loads[2i] = -1.0)
end

φinit = [-min(x_coords(gg)[i], y_coords(gg)[j], 1 - x_coords(gg)[i], 1 - y_coords(gg)[j])
         for i in 1:n, j in 1:n]

res = levelset_compliance_optimize!(gg, E, ν, fixed, loads, φinit;
                                    volfrac = 0.5, maxiter = 12)

println("compliance: ", round(res.history[1]; digits = 3), " → ", round(res.history[end]; digits = 3))
println("final volume fraction = ", round(sum(res.γ) / length(res.γ); digits = 3))

vtk(m; field = res.γ, fieldname = "level-set density")
compliance: 1122.129 → 1997.673
final volume fraction = 0.5

Next steps