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 = modelDifference{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 → insideorigin (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 2new 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
- Parametric CAD: engineering parts from scratch — more standard parts built feature by feature.
- Sketches, extrusion, revolution — build geometry from 2D sketches.
- Linear elasticity — simulate this same lug under load with the finite-element method.
