RasmahRasmah

Native BRep loads and boundary conditions

Overview

Loads and boundary conditions can be applied directly on the BRep — the CAD model — rather than on mesh indices. mesh_face_sets(b, m) classifies every boundary face of a volume mesh onto the BRep face it belongs to (nearest-surface triangle over the BRep tessellation), returning one FaceSet per BRep face. Those face sets then drive the descriptive FixedDisplacement / Force conditions — the same pattern mesh_with_labels performs in a single call.

The part is an 80 × 60 × 8 mm aluminium plate: clamped on one x-face and pulled with a 10 MPa traction on the other, then checked against the rod formula u = σL/E. Dimensions are millimetres.

using Rasmah

# an 80 × 60 × 8 plate, expressed in the same coordinates `tetrahedralize_box`
# uses (origin-based), so the mesh and the BRep tessellation align
part = translate(box(80.0, 60.0, 8.0), 40.0, 30.0, 4.0)
b = brep(part)

m = tetrahedralize_box(80.0, 60.0, 8.0, 16, 12, 3)
fs = mesh_face_sets(b, m)                # fs[i] is the FaceSet of BRep face i
println("BRep faces: ", length(fs), ", boundary faces: ",
        size(boundary_faces(m).faces, 2))
BRep faces: 6, boundary faces: 1104

Applying conditions on BRep faces

Pick the faces by their centroid (here the two x-faces) and constrain or load them through their face set:

U = FESpace(m, ReferenceFE(3, 1); valuetype = VectorValue{3})

function face_centroid(fs)
    nodes = m.nodes[:, unique(vec(fs.faces))]
    vec(sum(nodes; dims = 2)) ./ size(nodes, 2)
end
left  = argmin(i -> face_centroid(fs[i])[1], 1:length(fs))
right = argmax(i -> face_centroid(fs[i])[1], 1:length(fs))

# clamp the left face, pull the right face in +x with 10 MPa
fixed = dofs(U, fs[left])
loads = face_load(m, fs[right], [10.0, 0.0, 0.0])     # MPa = N/mm² traction
2652-element Vector{Float64}:
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  ⋮
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
 44.44444444444446
  0.0
  0.0

face_load distributes a traction over a face; dofs(sp, set) returns the field dofs of a point/edge/face set (component-aware, and for higher-order spaces it includes the edge/face midpoints lying on the boundary). Clamping a full face removes all rigid-body modes.

The solve is SI internally, so for a millimetre-and-MPa hand check the material constants are rescaled to N/mm² (E in MPa):

alu = override(aluminum; youngs_modulus = 6.89e4, poisson_ratio = 0.33)  # mm/MPa

u = solve(Elasticity(), m, alu;
          bc = BoundaryConditions(FixedDisplacement(fixed), Force(loads)))
println("max |u|      = ", round(maximum(abs, ustrip(u)); digits = 4), " mm")
println("analytic σL/E = ", round(10.0 * 80.0 / 6.89e4; digits = 4), " mm")
max |u|      = 0.0115 mm
analytic σL/E = 0.0116 mm

The end displacement lands within a percent of the rod formula σL/E — the same closed-form check a design office writes on the drawing margin. The stress field is the applied 10 MPa, with the usual boundary bumps at the clamped and loaded faces:

σ = von_mises_stress(m, ustrip(u), alu)
println("von Mises band = ", round(minimum(σ); digits = 1), " … ",
        round(maximum(σ); digits = 1), " MPa   (applied 10 MPa)")
println("plate mass (aluminum) = ", round(80.0 * 60.0 * 8.0 * 2.70e-3; digits = 1), " g")
von Mises band = 5.2 … 12.6 MPa   (applied 10 MPa)
plate mass (aluminum) = 103.7 g

Because the conditions live on the BRep faces, changing the geometry — a thicker plate, a hole through it — re-derives the face sets from the CAD model, and the same FixedDisplacement / Force block keeps working.

Note on mesh_with_labels

mesh_with_labels(b, n) is the one-call convenience that meshes a BRep and labels its boundary faces. It meshes through the SDF marching_tetrahedra path, which produces sliver tetrahedra and unreferenced grid nodes, so it is not yet suitable for FEM assembly. For FEM-grade solves, mesh with tetrahedralize_box (or any clean conforming mesher) and classify with mesh_face_sets, as above.

Next steps