Sets and Boundary Conditions

Overview

Boundary conditions and loads are the difference between a bare mesh and a solvable model. Earlier tutorials assembled them by hand (looping over node coordinates and filling forces vectors). Rasmah instead provides a small vocabulary of sets — PointSet/EdgeSet/FaceSet/ElementSet — plus helpers that turn a set into degrees of freedom (dofs), nodal loads (point_load/face_load), or named boundary faces (box_boundary_sets, boundary_sets). The same machinery has an SDF-side counterpart (point_set/point_sets/boundary_points/discontinuity_points) that labels points by which region of a shape they lie on.

This tutorial builds a plate in tension and drives the entire FEM setup — fixed face, load, solve, compliance — from sets, then tours the SDF labeling side.

The mesh and its boundary

boundary_faces returns the faces shared by exactly one tetrahedron (the mesh skin); boundary_edges returns the edges shared by exactly one tet (the skin’s wireframe):

using Rasmah
using LinearAlgebra

mesh = tetrahedralize_box(2.0, 1.0, 0.5, 4, 2, 1)

bf = boundary_faces(mesh)
be = boundary_edges(mesh)

println(size(mesh.nodes, 2), " nodes, ", size(mesh.elements, 2), " tetrahedra")
println(size(bf.faces, 2), " boundary faces, ", size(be.edges, 2), " boundary edges")

30 nodes, 48 tetrahedra
56 boundary faces, 14 boundary edges

The set types

Four concrete set types tag different mesh entities. They are thin wrappers around id vectors/matrices:

p = PointSet([1, 2, 3])          # node ids
el = ElementSet([4, 5])          # element ids
f = FaceSet(bf.faces[:, 1:2])    # 3 × n node-id matrix
e = EdgeSet(be.edges[:, 1:2])    # 2 × n node-id matrix

println(typeof(p), "  ", typeof(el))
println(typeof(f), "  ", typeof(e))

PointSet  ElementSet
FaceSet  EdgeSet

You rarely build a FaceSet/EdgeSet by slicing — the point is that boundary_faces, boundary_edges, box_boundary_sets, and boundary_sets produce them, and the rest of the API consumes them.

Selecting sets by a geometric predicate

nodes_where and faces_where select entities whose node coordinates satisfy a predicate — the clean way to pick “everything on the left face”:

xmin_nodes = nodes_where(mesh, n -> n[1] ≈ 0.0)
xmax_nodes = nodes_where(mesh, n -> n[1] ≈ 2.0)
println("x = 0 nodes: ", length(xmin_nodes.ids), ", x = 2 nodes: ", length(xmax_nodes.ids))

xmax_faces = faces_where(mesh, bf, n -> n[1] ≈ 2.0)
println("x = 2 faces: ", size(xmax_faces.faces, 2))

x = 0 nodes: 6, x = 2 nodes: 6
x = 2 faces: 4

Bounding-box regions (filter + box)

The most common way to place a boundary condition is “everything on this face / in this region”. BoundingBox(lo, hi) (or a (lo, hi) tuple) selects nodes or boundary faces whose coordinates fall inside the box, with an optional filter predicate to refine the selection:

# the x = 0 plane: a thin slab
left  = nodes_in_box(mesh, [0.0, 0.0, 0.0], [0.0, 1.0, 0.5])
right = faces_in_box(mesh, [2.0, 0.0, 0.0], [2.0, 1.0, 0.5])
println("left nodes: ", length(left.ids), ", right faces: ", size(right.faces, 2))

# a filter narrows the box (here: the x = 0 nodes on the z = 0 edge)
edge = nodes_in_box(mesh, BoundingBox([0.0, 0.0, 0.0], [0.0, 1.0, 0.5]);
                    filter = n -> n[3] ≈ 0.0)
println("edge nodes: ", length(edge.ids))

left nodes: 6, right faces: 4
edge nodes: 3

These regions feed straight into the solve keyword arguments — no dof loops:

u = solve(Elasticity(), mesh, steel;
          fixed = left, g = [0.0, 0.0, 0.0],
          loads = right => [1.0, 0.0, 0.0])   # a traction over the faces
println("max |u| = ", round(maximum(abs, u); digits = 6))

max |u| = 0.0 m

fixed accepts a region (or raw dofs) and g a scalar / per-node vector / per-dof vector / position function; loads accepts a region => value pair (a face traction, or a nodal force when the region is a PointSet). The same regions work on a FESpace (order-aware, so P2 mid-edge nodes are included) via dofs(sp, BoundingBox(...)).

Degrees of freedom

dofs maps any set type to the global displacement DOF indices of its nodes (each node owns 3 DOFs). It dispatches on the set type, so a FaceSet’s corner nodes are collected automatically:

fixed_dofs = dofs(mesh, xmin_nodes)          # from a PointSet
fixed_dofs2 = dofs(mesh, xmax_faces)         # from a FaceSet
println("fixed DOFs: ", length(fixed_dofs), ", loaded-face DOFs: ", length(fixed_dofs2))

fixed DOFs: 18, loaded-face DOFs: 18

This is the set-based replacement for the manual 3i-2, 3i-1, 3i loops in the linear-elasticity tutorial.

Loads

point_load distributes a force vector over a set of nodes; face_load integrates a constant traction over a set of faces (splitting each face’s load equally over its three corners). face_area gives the area of one face:

nodal = point_load(mesh, xmax_nodes, [1.0, 0.0, 0.0])     # 1 N total, split evenly
traction = face_load(mesh, xmax_faces, [2.0, 0.0, 0.0])   # 2 Pa over the face
println("point-load  nonzeros: ", count(!=(0), nodal))
println("face-load   nonzeros: ", count(!=(0), traction))
println("area of one loaded face: ", face_area(mesh, xmax_faces.faces[:, 1]))

point-load  nonzeros: 6
face-load   nonzeros: 6
area of one loaded face: 0.125

face_load is the physically correct way to apply a pressure/traction: it weights each node by the area of the faces it touches, whereas point_load is a lumped nodal force.

For a scalar field (temperature, potential, pressure) the same distributors exist one-dof-per-node: scalar_point_load and scalar_face_load build a 1-component-per-node load vector — the thermal/flux analogues of point_load / face_load (which build the 3-components-per-node displacement vector):

q = scalar_point_load(mesh, xmax_nodes, 1.0)       # 1 unit per node (a heat source)
qface = scalar_face_load(mesh, xmax_faces, 2.0)    # 2 units per area (a heat flux)
println("scalar point-load nonzeros: ", count(!=(0), q))
println("scalar face-load nonzeros:  ", count(!=(0), qface))
println("one dof per node:          ", length(q) == size(mesh.nodes, 2))

scalar point-load nonzeros: 6
scalar face-load nonzeros:  6
one dof per node:          true

The six named faces of a box

For the standard unit-cube load case, box_boundary_sets returns the six faces keyed by name:

sets = box_boundary_sets(mesh, 2.0, 1.0, 0.5)
for (name, fs) in sort(collect(pairs(sets)); by=p -> String(p[1]))
    println("  ", name, " -> ", size(fs.faces, 2), " faces")
end

  xmax -> 4 faces
  xmin -> 4 faces
  ymax -> 8 faces
  ymin -> 8 faces
  zmax -> 16 faces
  zmin -> 16 faces

A solve driven entirely by sets

Putting it together: fix the :xmin face, pull the :xmax face with a unit traction, and read off the compliance — no hand-written loops anywhere:

E, ν = 1e3, 0.3
steel = Material(youngs_modulus=E, poisson_ratio=ν)

fixed = dofs(mesh, sets[:xmin])
forces = face_load(mesh, sets[:xmax], [1.0, 0.0, 0.0])

u = solve(Elasticity(), mesh, steel; fixed=fixed, loads=forces)
println("max |u|        = ", round(maximum(abs, ustrip(u)); digits=6))
println("compliance     = ", round(compliance(mesh, steel, fixed, forces); digits=6))

vtk(mesh; field=ustrip(u)[1:3:end], fieldname="u_x")

max |u|        = 0.002024
compliance     = 0.000967

SDF point classification

Sets also exist on the implicit side. point_set(g, x) labels a point against an SDF. The generic case returns :interior/:exterior/:surface; a box shape additionally returns the face labels :xmin…:zmax — and, for points on an edge or corner, several labels at once:

g = to_sdf(box(2.0, 2.0, 2.0))

println(point_set(g, [0.0, 0.0, 0.0]))    # inside
println(point_set(g, [1.0, 0.0, 0.0]))    # on the +x face
println(point_set(g, [1.0, 1.0, 0.0]))    # on an edge → two labels
println(point_set(g, [3.0, 0.0, 0.0]))    # outside

Set([:interior])
Set([:xmax])
Set([:xmax, :ymax])
Set([:exterior])

point_sets maps over a list; boundary_points keeps only the points that are on a boundary (neither interior nor exterior); and discontinuity_points keeps the multi-label points — the edges/corners that a mesher must respect:

xs = [[0.0, 0, 0], [1.0, 0, 0], [1.0, 1.0, 0], [1.0, 1.0, 1.0], [3.0, 0, 0]]

println("sets:              ", point_sets(g, xs))
println("boundary points:   ", boundary_points(g, xs))
println("discontinuity pts: ", discontinuity_points(g, xs))

sets:              Set{Symbol}[Set([:interior]), Set([:xmax]), Set([:xmax, :ymax]), Set([:zmax, :xmax, :ymax]), Set([:exterior])]
boundary points:   [[1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [1.0, 1.0, 1.0]]
discontinuity pts: [[1.0, 1.0, 0.0], [1.0, 1.0, 1.0]]

Generic boundary sets from an SDF

boundary_sets(mesh, g) generalizes box_boundary_sets to the SDF-driven path: it classifies each boundary face by the point_set label of its centroid. The named labels (:xmin…:zmax) come from the SDF’s own point_set, so a box reproduces the six box faces without any hard-coded coordinates.

Note that a box shape is centred at the origin (box(2,2,2) spans [-1,1]³), whereas tetrahedralize_box spans the positive octant [0,w]³ — so we shift the mesh to align it with the SDF first:

g = to_sdf(box(2.0, 2.0, 2.0))                      # spans [-1,1]³
m = tetrahedralize_box(2.0, 2.0, 2.0, 2, 2, 2)          # spans [0,2]³
m = TetMesh(m.nodes .- 1.0, m.elements)                 # shift to match the SDF

bsets = boundary_sets(m, g)
for (name, fs) in sort(collect(pairs(bsets)); by=p -> String(p[1]))
    println("  ", name, " -> ", size(fs.faces, 2), " faces")
end

  xmax -> 8 faces
  xmin -> 8 faces
  ymax -> 8 faces
  ymin -> 8 faces
  zmax -> 8 faces
  zmin -> 8 faces

A sphere or cylinder has no face labels (its point_set returns only :interior/:exterior/:surface), so boundary_sets is empty for them — the named sets are a property of the labeled shape, and boundary_sets is how that labeling flows onto a mesh (and, as shown in the BRep tutorials, onto a BRep or surface mesh too).

Named mesh data

A mesh can carry arbitrary named fields on its nodes, cells, or faces — the substrate for field visualization and for storing external-solver results. The set_*_data! functions attach a field (returning a new mesh with the data attached) and point_data/cell_data/face_data read the dictionaries back:

m = tetrahedralize_box(1.0, 1.0, 1.0, 2, 2, 2)

m = set_point_data!(m, "temperature", 300.0 .+ m.nodes[1, :] .* 100.0)
m = set_cell_data!(m, "group", fill(2, size(m.elements, 2)))

println("data container type: ", typeof(m.data))   # MeshData
println(point_data(m)["temperature"])
println("cell field 'group' length: ", length(cell_data(m)["group"]))
println("face data (empty): ", face_data(m))

vtk(m; field=point_data(m)["temperature"], fieldname="temperature")

data container type: MeshData
[300.0, 350.0, 400.0, 300.0, 350.0, 400.0, 300.0, 350.0, 400.0, 300.0  …  400.0, 300.0, 350.0, 400.0, 300.0, 350.0, 400.0, 300.0, 350.0, 400.0]
cell field 'group' length: 48
face data (empty): Dict{String, Any}()

Next steps

• Linear elasticity — the FEM solve these sets drive.
• Differentiable simulation and optimization — differentiate the compliance assembled here.
• Boundary representation — the same boundary_sets concept at the exact BRep level.