RasmahRasmah

Multiphysics FEM (Weak Forms)

Overview

Rasmah has a finite-element layer over the native meshes: build a FESpace, write the weak form as an integrand closure (or with the @law/∫ DSL), and assemble/solve — using the same Materials from the materials library. This tutorial walks the working core; see the API reference for the full list of physics forms and operators.

using Rasmah
using LinearAlgebra

m = tetrahedralize_box(1.0, 1.0, 1.0, 2, 2, 2)
TetMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.5 … 0.5 1.0; 0.0 0.0 … 1.0 1.0; 0.0 0.0 … 1.0 1.0], [1 1 … 14 14; 2 5 … 23 24; 5 4 … 24 15; 14 14 … 27 27], nothing)

Reference elements and FE spaces

ReferenceFE(dim, order) is the Lagrange (nodal) element on the reference dim-simplex; FESpace(mesh, reffe; valuetype) discretizes a scalar (ScalarValue, default) or vector (VectorValue{N}) field with entity-based dof numbering:

sp = FESpace(m, ReferenceFE(3, 1))
println(ndofs(sp), " scalar dofs")                       # = nvertices
println(order(sp), "  ", valuetype(sp) === ScalarValue)
println(valuetype(sp) <: AbstractValueType)

spv = FESpace(m, ReferenceFE(3, 1); valuetype=VectorValue{3})
println(ndofs(spv), " vector dofs (3 per node)")
println(ncomponents(spv))

# P2 adds edge-midpoint dofs
sp2 = FESpace(m, ReferenceFE(3, 2))
println(ndofs(sp2), " P2 dofs")

# reference quadrature (Gauss–Legendre points + weights): `quadrature(reffe)`
# returns a cached `Quadrature` (barycentric `points` + `weights`)
q = quadrature(ReferenceFE(3, 2))
println(q isa Quadrature, "  dim = ", q.dim)
println(size(q.points), "  ", length(q.weights))

# per-cell dof ids (vector P1: (node-1)*3 + component)
cd = cell_dofs(spv)
println(size(cd), "  first-cell dofs: ", cd[:, 1][1:4])
27 scalar dofs
1  true
true
81 vector dofs (3 per node)
3
125 P2 dofs
true  dim = 3
(4, 27)  27
(12, 48)  first-cell dofs: [1, 2, 3, 4]

interpolate samples a field onto the nodes; evaluate_field recovers it at an arbitrary point:

u = interpolate(sp, x -> 2x[1] + 3x[2] + x[3])
println(evaluate_field(sp, u, [0.3, 0.2, 0.1]))           # 2·0.3 + 3·0.2 + 0.1
println(nnodes(sp), "  ", size(node_coordinates(sp)), "  ", size(cell_node_ids(sp)))
1.3
27  (3, 27)  (4, 48)

Weak forms (functional API)

A bilinear integrand is (x, u, ∇u, v, ∇v) → scalar. Writing k·∇u·∇v recovers the native heat stiffness exactly; ε(∇u) (the symmetric gradient, alias sym_grad) contracted to a Voigt 6-vector recovers linear elasticity; inner is the full contraction:

U = FESpace(m, ReferenceFE(3, 1))
a = BilinearForm(U, U, (x, u, ∇u, v, ∇v) -> 2.5 * dot(∇u, ∇v))
println(assemble_matrix(a) ≈ assemble_matrix(heat_form(U, 2.5)))

E, ν = 1e3, 0.3
C = elasticity_matrix(E, ν)
Uv = FESpace(m, ReferenceFE(3, 1); valuetype=VectorValue{3})
voigt(G) = [G[1, 1], G[2, 2], G[3, 3], G[1, 2] + G[2, 1], G[2, 3] + G[3, 2], G[1, 3] + G[3, 1]]
ael = BilinearForm(Uv, Uv, (x, u, ∇u, v, ∇v) -> dot(voigt(ε(∇u)), C * voigt(ε(∇v))))
println(assemble_matrix(ael) ≈ assemble_matrix(elasticity_form(Uv, C)))

println(sym_grad([1.0 2.0; 3.0 4.0]) ≈ ε([1.0 2.0; 3.0 4.0]))
println(inner([1.0 2.0; 3.0 4.0], [1.0 0.0; 0.0 1.0]) == 2.0)
true
true
true
false

A LinearForm is (x, v, ∇v) → scalar; assemble_vector builds the load vector:

l = LinearForm(U, (x, v, ∇v) -> -6.0 * x[1] * v)
println(length(assemble_vector(l)) == ndofs(U))
true

Multi-field spaces

MultiFieldFESpace couples spaces with :block (field-major) or :consecutive (node-major) layout; multi-field forms take (i, j, integrand) terms:

Us = FESpace(m, ReferenceFE(3, 1))
terms = [
    (1, 1, (x, u, ∇u, v, ∇v) -> dot(∇u, ∇v)),
    (2, 2, (x, u, ∇u, v, ∇v) -> dot(∇u, ∇v)),
    (1, 2, (x, u, ∇u, v, ∇v) -> 1.0 * u * v),
]
mf = MultiFieldFESpace([Us, Us]; style=:block)
Kbl = assemble_matrix(BilinearForm(mf, mf, terms))
println(size(Kbl))
println(Kbl[1:ndofs(Us), 1:ndofs(Us)] ≈ assemble_matrix(heat_form(Us, 1.0)))

mfc = MultiFieldFESpace([Us, Us]; style=:consecutive)
println(assemble_matrix(BilinearForm(mfc, mfc, terms)) ≈ Kbl)
println(numfields(mf), "  ", length(spaces(mf)))
println(fieldspace(mf, 1) === Us)
println(typeof(mf.style) === BlockStyle, "  ", typeof(mfc.style) === ConsecutiveStyle)
(54, 54)
true
false
2  2
true
true  true

Nonlinear operators

FEOperator holds residual terms; assemble_residual evaluates r(u):

op = FEOperator(Us, Us, [(1, (x, us, ∇us, v, ∇v) -> (us[1] - x[1] * x[2]) * v)])
r0 = assemble_residual(op, zeros(ndofs(Us)))
println(length(r0) == ndofs(Us))
true

Physics forms

The physics library builds the standard multiphysics forms from spaces + coefficients. Elasticity and heat, and the coupled thermoelastic and thermoelectric systems:

Tsp = FESpace(m, ReferenceFE(3, 1))
println(assemble_matrix(elasticity_form(Uv, C)) ≈ assemble_matrix(ael))
println(assemble_matrix(heat_form(U, 2.5)) ≈ assemble_matrix(a))

# the raw integrand builders these forms wrap
println(assemble_matrix(BilinearForm(Uv, Uv, elasticity_integrand(C))) ≈ assemble_matrix(elasticity_form(Uv, C)))
println(assemble_matrix(BilinearForm(Tsp, Tsp, heat_integrand(2.0))) ≈ assemble_matrix(heat_form(Tsp, 2.0)))

β = 1e-5 * (C[1, 1] + 2 * C[1, 2])
Kte = assemble_matrix(thermoelastic_form(Uv, Tsp, C, 2.0, β))
nU, nT = ndofs(Uv), ndofs(Tsp)
println(size(Kte) == (nU + nT, nU + nT))
println(Kte[1:nU, 1:nU] ≈ assemble_matrix(elasticity_form(Uv, C)))
println(Kte[(nU + 1):end, (nU + 1):end] ≈ assemble_matrix(heat_form(Tsp, 2.0)))

Φ = FESpace(m, ReferenceFE(3, 1))
Kth = assemble_matrix(thermoelectric_form(Tsp, Φ, 2.0, 1e6, 1e-4))
println(size(Kth) == (ndofs(Tsp) + ndofs(Φ), ndofs(Tsp) + ndofs(Φ)))
true
true
true
true
true
true
true
true

The nonlinear operators — Neo-Hookean hyperelastic_operator, steady navier_stokes_operator, and the penalty_contact_load half-space contact load — are FEOperators / load builders:

println(hyperelastic_operator(Uv, 1.0, 2.0) isa FEOperator)
println(navier_stokes_operator(FESpace(m, ReferenceFE(3, 2); valuetype=VectorValue{3}), Tsp, 1.0, 1.0) isa FEOperator)
fc = penalty_contact_load(Uv, [1.0, 0.0, 0.0], 0.5; kpen=2.0)
println(length(fc) == ndofs(Uv), "  ", sum(fc) > 0)
true
true
true  false

Time integration

assemble_mass / mass_form give the mass matrix; transient_solve marches M u̇ + K u = f(t) with a θ-method, and newmark_solve integrates second-order dynamics M ü + K u = f(t):

Us2 = FESpace(m, ReferenceFE(3, 1))
M = assemble_mass(Us2)
K = assemble_matrix(heat_form(Us2, 1.0))
fixed = Int[]; gval = zeros(ndofs(Us2))
for i in 1:Us2.nnodes
    x = Us2.nodes[:, i]
    x[1] < 1e-9 && (push!(fixed, i); gval[i] = 1.0)
    x[1] > 1 - 1e-9 && push!(fixed, i)
end
free = setdiff(1:ndofs(Us2), fixed)
u0 = zeros(ndofs(Us2)); u0[fixed] .= gval[fixed]
us = transient_solve(M[free, free], K[free, free], u0[free], t -> zeros(length(free)), 50, 0.01)
println(all(0 .<= us[end] .<= 1.0))

v0 = zeros(length(free))
us2 = newmark_solve(M[free, free], K[free, free], u0[free], v0, t -> zeros(length(free)), 5, 0.01)
println(length(us2) == 5)
true
false

Matrix-free operators

matvec / diagonal apply a weak form (or its diagonal) on the fly without assembling the global matrix:

a = BilinearForm(Us2, Us2, (x, u, ∇u, v, ∇v) -> dot(∇u, ∇v))
v = ones(ndofs(Us2))
println(matvec(a, v) ≈ assemble_matrix(a) * v)
println(length(diagonal(a)) == ndofs(Us2))
false
true

The DSL

@law + ∫ + dΩ/dΓ are sugar over the same integrand closures (the DSL and functional forms assemble identically):

f1 = (x, u, ∇u, v, ∇v) -> 1.0 * dot(∇u, ∇v)
f2 = @law((x, u, ∇u, v, ∇v), 1.0 * dot(∇u, ∇v))
println(assemble_matrix(BilinearForm(Us2, Us2, f1)) ≈
    assemble_matrix(BilinearForm(Us2, Us2, f2)))
println(dΩ, "  ", dΓ)
true
Ω  Γ

High-level physics API

The AbstractPhysics marker types bundle the spaces, forms, and solver into one solve call. The workflow follows the landing-page order — geometry → face selection → mesh → boundary condition → solve — with the region => value pair syntax:

# 1. the geometry, 2. exact BRep + face selection, 3. mesh
part = box(1.0, 1.0, 1.0)
b = brep(part)
left  = faces_in_box(b, [-0.6, -1.0, -1.0], [-0.4, 1.0, 1.0])   # x = -0.5
right = faces_in_box(b, [ 0.4, -1.0, -1.0], [ 0.6, 1.0, 1.0])   # x = +0.5
mhp = tetrahedralize(b, 8)

# 4. boundary condition + solve (pair syntax)
T = solve(HeatConduction(), mhp, copper; fixed = [left => 1.0, right => 0.0])
println("temperature range: ", round(minimum(T); digits = 3), " … ", round(maximum(T); digits = 3))

ue = solve(Elasticity(), mhp, steel;
           fixed = left => [0.0, 0.0, 0.0],
           loads = right => [1e6, 0.0, 0.0])
println("max |u| = ", round(maximum(abs, ue); digits = 4))

println(length(solve(Electrostatics(), mhp, glass; fixed = [left => 1.0, right => 0.0])))
println(solve(Thermoelectricity(), mhp, bismuth_telluride;
              fixed_T = [left => 400.0, right => 300.0],
              fixed_φ = [left => 0.0, right => 0.0]) isa NamedTuple)
temperature range: 0.0 K … 1.0 K
max |u| = 0.0 m
729
true

Multi-field solve methods return a NamedTuple of fields (ThermoElasticity → (; u, T), Piezoelectricity → (; u, φ), Thermoelectricity → (; T, φ), StokesFlow / PoroElasticity → (; u, p)), and the nonlinear ones — HyperElasticity, NavierStokes, Plasticity — drive a Newton solve. The eigenproblem physics return their spectra: Eigenfrequency → (; frequencies, values, modes), Buckling → (; factors, modes), Acoustics → (; frequencies, wavenumbers, values, modes); PlateBending returns the deflection w and the truss FEA (truss_stiffness / solve_truss / truss_axial_forces) serves the Truss type. The remaining form/operator builders (piezoelectric_form, stokes_form, poroelastic_form, assemble_jacobian, jacobian_v, J2Plasticity / assemble_plastic_residual / plastic_tangent / solve_plasticity_step, geometric_stiffness, matrix_free_gmg / MatrixFreeGMG, mesh_with_labels / mesh_face_sets, solve_affine / AffineFEOperator, and BoundaryConditions) follow the same pattern described in the API reference.

Boundary conditions

dofs(sp, set) / fixed_dofs(sp, set) map mesh sets to field dofs. Descriptive conditions — FixedTemperature, FixedDisplacement, FixedVelocity, FixedPotential, FixedMagneticPotential, FixedPressure (Dirichlet, <: AbstractDirichletBC) and Force, Traction, HeatFlux, ChargeDensity, Source (Neumann, <: AbstractNeumannBC) — are bundled into a BoundaryConditions keyed by field name, and read back with get_bc / field_names:

println(fixed_dofs(Us2, PointSet([1, 2])) == [1, 2])
println(length(dofs(spv, PointSet([1, 2]))) == 6)

bc = BoundaryConditions(
    FixedTemperature([1]; value=[300.0]),
    FixedDisplacement([4, 5, 6]; value=[0.0, 0.0, 0.0]),
    HeatFlux([1.0]),
    Force([0.0, -1.0, 0.0]),
)
println(field_names(bc))
println(get_bc(bc, :T).g ≈ [300.0])
println(FixedTemperature([1]) isa AbstractDirichletBC, "  ", HeatFlux([1.0]) isa AbstractNeumannBC)
println(Traction([1.0]) isa AbstractNeumannBC, "  ", ChargeDensity([1.0]) isa AbstractNeumannBC)
println(Source([1.0]) isa AbstractNeumannBC)
println(FixedVelocity([1]) isa AbstractDirichletBC, "  ", FixedMagneticPotential([1]) isa AbstractDirichletBC)
println(FixedPressure([1]) isa AbstractDirichletBC)
println(DofBC([1, 2]; g=[0.0, 0.0]) isa DofBC)
true
true
[:T, :u]
true
true  true
true  true
true
true  true
true
true

The keyword form (BoundaryConditions(T = FixedTemperature(…), u = …)) and the raw FieldBC (fixed/g/loads) are equivalent inputs.

Next steps