RasmahRasmah

Constructive Solid Geometry and Meshing

println("Julia ", VERSION)
using Rasmah

# Geometry is stored in SI (metres); these helpers render millimetres.
mm(v)  = round(uconvert(u"mm", v); digits = 1)
mm2(v) = round(uconvert(u"mm^2", v); digits = 0)
mm3(v) = round(uconvert(u"mm^3", v); digits = 0)
Julia 1.12.6
mm3 (generic function with 1 method)

Overview

Constructive solid geometry (CSG) builds complex shapes by combining simpler ones with boolean operations and transforming them in space. In this tutorial we assemble a real part — a steel lifting lug, the 120 × 80 × 20 mm plate with a Ø32 pin hole and two Ø10 mounting holes that gets welded to the side of a machine frame — then turn it into a surface mesh and a volume mesh. Everything remains purely geometric; no simulation yet. Dimensions are in millimetres; the code writes them as metres.

Boolean operations

Rasmah provides the three standard boolean operations as features, each with a signed-distance-field (SDF) rule:

operation SDF rule meaning
difference max(a, -b) a minus b
intersection max(a, b) overlap of a,b
union (+/∪) min(a, b) both shapes

The lug is a box with three cylinders drilled out of it. The plate is centred at the origin, so the pin hole sits in the middle with no extra transform; the two mounting holes move to a 70 mm pitch with translate:

model = difference(
    difference(
        difference(
            box(0.120, 0.080, 0.020),                # the plate: 120 × 80 × 20 mm
            cylinder(0.016, 0.020),                  # the Ø32 pin hole, drilled through
        ),
        translate(cylinder(0.005, 0.020), -0.035, 0.0, 0.0),  # Ø10 mounting hole
    ),
    translate(cylinder(0.005, 0.020), 0.035, 0.0, 0.0),       # …and its mate
)

lug = model
Difference{Difference{Difference{Box{Float64, Float64, Float64}, Cylinder{Float64, Float64}}, Translate{Cylinder{Float64, Float64}, Float64, Float64, Float64}}, Translate{Cylinder{Float64, Float64}, Float64, Float64, Float64}}(Difference{Difference{Box{Float64, Float64, Float64}, Cylinder{Float64, Float64}}, Translate{Cylinder{Float64, Float64}, Float64, Float64, Float64}}(Difference{Box{Float64, Float64, Float64}, Cylinder{Float64, Float64}}(Box{Float64, Float64, Float64}(0.12, 0.08, 0.02), Cylinder{Float64, Float64}(0.016, 0.02)), Translate{Cylinder{Float64, Float64}, Float64, Float64, Float64}(Cylinder{Float64, Float64}(0.005, 0.02), -0.035, 0.0, 0.0)), Translate{Cylinder{Float64, Float64}, Float64, Float64, Float64}(Cylinder{Float64, Float64}(0.005, 0.02), 0.035, 0.0, 0.0))

Sampling the field shows the classic CSG behaviour: the field of a difference is max(a, -b), so points inside a hole (where the cylinder’s field is negative) are outside the final shape:

println("origin (in the pin hole) : ", field(lug, [0.0, 0.0, 0.0]))   # positive → outside
println("inside the plate arm    : ", field(lug, [0.045, 0.0, 0.0]))  # negative → inside
origin (in the pin hole) : 0.01
inside the plate arm    : -0.004999999999999995

Transformations

translate, scale, and rotate wrap any geometry and shift, resize, or reorient it in space. rotate takes a 3×3 rotation matrix, and Rasmah ships helpers rotation_x/y/z for the usual axes:

using LinearAlgebra

# translate: move a unit sphere 2 units along +x
moved = translate(sphere(1.0), 2.0, 0.0, 0.0)
println("new centre (x=2)   : ", field(moved, [2.0, 0.0, 0.0]))   # -1 → inside
println("shifted boundary   : ", field(moved, [1.0, 0.0, 0.0]))   #  0 → surface at x=1
println("old centre (origin): ", field(moved, [0.0, 0.0, 0.0]))   #  1 → now outside

# rotate: lay a z-aligned cylinder onto the x axis
rotated = rotate(cylinder(0.3, 2.0), rotation_y(π / 2))
println("rotated end cap    : ", field(rotated, [1.0, 0.0, 0.0]))   #  0 → cap now at x=1

# scale: enlarge a unit sphere to radius 2
scaled = scale(sphere(1.0), 2.0)
println("scaled centre      : ", field(scaled, [0.0, 0.0, 0.0]))    # -2 → radius 2
new centre (x=2)   : -1.0
shifted boundary   : 0.0
old centre (origin): 1.0
rotated end cap    : 0.0
scaled centre      : -2.0

Smooth (blended) boolean operations

Exact booleans create sharp max/min creases, which are non-smooth (their gradient is discontinuous). For differentiable pipelines it is often useful to blend the transition with a smoothing parameter k — the SDF answer to the fillet weld that joins two plates with a rounded throat:

exact = difference(sphere(1.0), translate(sphere(1.0), 1.0, 0.0, 0.0))
smooth = smooth_difference(sphere(1.0), translate(sphere(1.0), 1.0, 0.0, 0.0), 0.3)

println("exact difference at (0,0,0)  : ", field(exact, [0.0, 0.0, 0.0]))
println("smooth difference at (0,0,0) : ", field(smooth, [0.0, 0.0, 0.0]))
exact difference at (0,0,0)  : -0.0
smooth difference at (0,0,0) : -0.0

k controls the radius of the blend: k = 0 recovers the exact boolean, and larger k rounds the joint more.

Marching cubes (surface mesh)

An SDF is a continuous description; to visualize or measure it we need a mesh. marching_cubes extracts a triangle mesh of the level set f(x) = 0 on a grid. We take the shape’s bounding box and pad it slightly so the surface is not flush with the grid boundary:

lo, hi = bounds(lug)
margin = 0.05 .* (hi .- lo)

mesh = marching_cubes(lug, lo .- margin, hi .+ margin, 32)

println("vertices: ", size(mesh.vertices, 2))
println("faces:    ", size(mesh.faces, 2))

vtk(mesh)
vertices: 26064
faces:    13040

After the field, the three positional arguments are (lo, hi, n): lo and hi define the sampling box, and n = 32 is the grid resolution — the box is split into 32 steps along each axis, i.e. 32 × 32 × 32 = 32768 cells. Marching cubes checks each cell and emits triangles wherever the surface crosses its edges, so a larger n gives a finer, more accurate mesh (more triangles) at a higher cost.

surface_mesh does exactly this for you: it computes the bounding box, pads it, and runs marching cubes. The keyword n is the grid resolution — the box is split into an n × n × n grid, so larger n produces a finer mesh (more triangles, more accurate area) at a higher cost. surface_area integrates the triangle areas, which converge to the analytic total — the plate’s six faces plus the three hole walls, minus the material the holes removed:

println("n = 32 → area ≈ ", mm2(surface_area(surface_mesh(lug; n=32))), " mm²")
println("n = 64 → area ≈ ", mm2(surface_area(surface_mesh(lug; n=64))), " mm²")
println("analytic       : ",
        120 * 80 * 2 + 120 * 20 * 2 + 80 * 20 * 2 - 2 * π * (16^2 + 2 * 5^2)
        + 2π * 16 * 20 + 2 * 2π * 5 * 20, " mm²")
n = 32 → area ≈ 27674.0 mm²
n = 64 → area ≈ 28463.0 mm²
analytic       : 28544.601655736435 mm²

Mesh smoothness

The resolution n is a direct cost/quality dial: a coarse grid produces a faceted approximation of the curved holes, while a fine grid converges to the smooth surface. Here is the same lug at three resolutions, with edges shown so the triangle size is visible:

vtk(marching_cubes(lug, lo .- margin, hi .+ margin, 16); edges=true)
vtk(marching_cubes(lug, lo .- margin, hi .+ margin, 32); edges=true)
vtk(marching_cubes(lug, lo .- margin, hi .+ margin, 48))

At n = 16 the facets around the pin hole are obvious; at n = 48 the surface is visually smooth (and surface_area has converged to the analytic value).

Marching tetrahedra (volume mesh)

For solid simulation we need a tetrahedral volume mesh. Rasmah has a native marching-tetrahedra implementation that fills the interior with tetrahedra whose vertices lie exactly on the implicit surface. It takes the same (lo, hi, n) arguments as marching_cubes; here n = 32 is the grid resolution (32 × 32 × 32 cells):

volmesh = marching_tetrahedra(lug, lo .- margin, hi .+ margin, 32)

println("nodes:    ", size(volmesh.nodes, 2))
println("elements: ", size(volmesh.elements, 2))
nodes:    46964
elements: 200928

For a volume mesh, vtk draws the boundary surface of the tetrahedra:

vtk(volmesh)

We use a coarser grid than for the surface because each crossed cell can spawn several tetrahedra, so volume meshing costs much more per cell.

volume on a TetMesh sums the signed tetrahedron volumes, giving a numerical volume estimate that approaches the analytic value as the grid refines. The analytic value is exact arithmetic — the plate minus one pin hole and two mounting holes:

v_numeric  = volume(volmesh)
v_analytic = volume(box(0.120, 0.080, 0.020)) - volume(cylinder(0.016, 0.020)) -
             2 * volume(cylinder(0.005, 0.020))
println("numeric volume  ≈ ", mm3(v_numeric), " mm³")
println("analytic volume = ", mm3(v_analytic), " mm³")
println("mass (steel)    ≈ ",
        round(uconvert(u"kg", 7850u"kg/m^3" * v_numeric); digits = 2), " kg")
numeric volume  ≈ 172106.0 mm³
analytic volume = 172773.0 mm³
mass (steel)    ≈ 1.35 kg

Next steps