feat(src): add loads.jl

src/assemblers/matrix_free/loads.jl | 716 ++++++++++++++++++++++++++++++++++++ 1 file changed, 716 insertions(+)
2026-09-22 10:48:32 +00:00 · 2026-05-09 16:30:35 +03:00
parent 8ac04798a0
commit 20b405429a
1 changed files with 716 additions and 0 deletions
@@ -0,0 +1,716 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+"""
+Neumann (load) boundary conditions for the DOF-based assembler.
+
+Counterpart of `dirichlet.jl` for the right-hand-side: declarative
+load types that accumulate into `f` (or any user-provided vector) via a
+single `apply_load!(f, load, cache, asm, kernel, mesh)` entry point.
+Each concrete subtype implements the integration appropriate to it,
+reusing the SoA `N_data` / `detJ_w` batches the DOF-based assembler
+already builds in `_prepare_caches!`.
+
+Load types provided include:
+
+  * `NodalForce(dofs, values)` — point loads on known global DOF
+    indices. Allocation-free and trivially zero-cost; just an indexed
+    `f[d] += value` loop.
+
+  * `UniformBodyForce(b)` — constant body force per unit volume over
+    the entire mesh. Computes `f[i, α] += ∫_Ω N_i · b_α dV` per element
+    using the same batched geometry the matrix-free `apply_M!` /
+    `apply_K!` rely on. Vector-valued for elasticity (`Vec{3,Float64}`),
+    scalar-valued for heat (`Float64` heat source per unit volume).
+
+  * `SurfaceLoad(faces, traction)` — distributed traction (or heat
+    flux) over a list of mesh faces. Computes
+    `f[i, α] += ∫_Γ N_i · t_α dS` per face using a face-element
+    Gauss quadrature. Vector-valued for elasticity, scalar-valued
+    for heat. The natural complement of `UniformBodyForce`: body
+    force lives in the volume integral, surface traction lives in
+    the surface integral.
+
+  * [`UniformMixedDarcySource`](@ref) — adds ``f \\cdot |K_e|`` to each cell
+    pressure DOF for mixed RT₀–P₀ Darcy (second equation); works with
+    [`AbstractDarcyMixedRT0P0Kernel`](@ref) (Tet4 or Hex8).
+
+  * [`MixedDarcyTet4BoundaryNormalFluxLoad`](@ref) — Tet4: uniform normal flux density
+    ``g`` on panels `(elem_id, local_face)` (triangle quadrature).
+  * [`MixedDarcyHex8BoundaryNormalFluxLoad`](@ref) — Hex8: same weak form on quad faces
+    (four Gauss points per face).
+
+`apply_load!` is *additive* on `f`, so multiple loads compose
+naturally:
+
+```julia
+apply_load!(f, body,    cache, asm, kernel, mesh)
+apply_load!(f, surface, cache, asm, kernel, mesh)
+apply_load!(f, point,   cache, asm, kernel, mesh)
+```
+
+Combined with `EliminatedDirichlet`'s `apply_constraint!(K, b, c)`
+"lift" of the RHS and the matrix-free `op` returned by
+`matrix_free_op`, this gives a complete declarative BC story for both
+direct and Krylov solves.
+"""
+
+abstract type AbstractNeumannLoad end
+
+# ----------------------------------------------------------------------------
+# NodalForce
+# ----------------------------------------------------------------------------
+
+"""
+    NodalForce(dofs, values)
+
+Concentrated point loads on a list of global DOF indices. Stores the
+DOF index list and the corresponding force values; `apply_load!` does
+`f[d] += value` for each pair.
+
+Useful for prescribed reactions, manually assembled surface tractions,
+or for verifying correctness of the body-force / matrix-free path
+against a reference assembled solve.
+"""
+struct NodalForce{IT<:AbstractVector{<:Integer},
+                  VT<:AbstractVector{Float64}} <: AbstractNeumannLoad
+    dofs::IT
+    values::VT
+
+    function NodalForce(dofs::IT, values::VT) where {
+            IT<:AbstractVector{<:Integer}, VT<:AbstractVector{Float64}}
+        @assert length(dofs) == length(values) (
+            "NodalForce: dofs ($(length(dofs))) and values ($(length(values))) " *
+            "must have the same length")
+        return new{IT, VT}(dofs, values)
+    end
+end
+
+NodalForce(dofs::AbstractVector{<:Integer},
+           values::AbstractVector{<:Real}) =
+    NodalForce(dofs, Float64.(values))
+
+"""
+    apply_load!(f::AbstractVector{Float64}, load::NodalForce,
+                cache::DOFBasedCOOCache, asm, kernel, mesh) -> f
+
+Accumulate point loads into `f`: `f[load.dofs[k]] += load.values[k]`.
+Allocation-free; ignores `cache`/`asm`/`kernel`/`mesh` (kept in the
+signature for API symmetry with body / surface loads).
+"""
+@inline function apply_load!(f::AbstractVector{Float64},
+                             load::NodalForce,
+                             cache::DOFBasedCOOCache,
+                             asm::DOFBasedCOOAssembler,
+                             kernel::AbstractKernel,
+                             mesh::AbstractMesh)
+    @inbounds for k in eachindex(load.dofs)
+        f[load.dofs[k]] += load.values[k]
+    end
+    return f
+end
+
+# ----------------------------------------------------------------------------
+# UniformBodyForce
+# ----------------------------------------------------------------------------
+
+"""
+    UniformBodyForce(b)
+
+Constant body force per unit volume over the whole mesh.
+
+* For elasticity (`ContinuumKernel`, 3 DOFs/node), `b` is a
+  `Vec{3,Float64}` — body force per unit volume, e.g. gravity
+  `Vec{3}((0.0, 0.0, -ρ * 9.81))`.
+* For heat (`HeatKernel`, 1 DOF/node), `b` is a `Float64` — heat
+  source per unit volume `[W/m³]`.
+
+`apply_load!` integrates `f[i, α] += ∫ N_i b_α dV` element by element,
+component by component, using the SoA `N_data` and `detJ_w` batches
+from the cache. Allocation-free after `_prepare_caches!` is warm.
+
+Variable / spatially-dependent body forces drop in by introducing a new
+type and overriding `_body_component(load, X, comp)` — the integration
+loop in `apply_load!` is identical.
+"""
+struct UniformBodyForce{V} <: AbstractNeumannLoad
+    b::V
+end
+
+# Component accessors — make the same integration loop work for both
+# vector (elasticity) and scalar (heat) loads. `comp` is the component
+# index from `DOFLayoutEntry` (always 1 for a scalar field).
+@inline _body_component(b::Vec{N,Float64}, comp::Int) where {N} = b[comp]
+@inline _body_component(b::Real, comp::Int) = Float64(b)
+
+"""
+    apply_load!(f::AbstractVector{Float64}, load::UniformBodyForce,
+                cache::DOFBasedCOOCache, asm, kernel, mesh) -> f
+
+Element-by-element assembly of `f[i, α] += ∫ N_i b_α dV` for every
+local DOF. Reuses the cache's SoA `N_data` and `detJ_w` batches and the
+compile-time `local_dof_layout` table, so the inner loop is fully
+type-stable and allocation-free after warmup.
+"""
+function apply_load!(f::AbstractVector{Float64},
+                     load::UniformBodyForce,
+                     cache::DOFBasedCOOCache{T,B,IPS,E,GC,Buf,FieldType,StateType},
+                     asm::DOFBasedCOOAssembler,
+                     kernel::AbstractKernel,
+                     mesh::AbstractMesh) where {T,B,IPS,E<:AbstractElement,
+                                                GC,Buf,FieldType,StateType}
+    @assert length(f) == cache.ndofs (
+        "apply_load!: f has length $(length(f)); expected $(cache.ndofs)")
+
+    # We reuse the same Pass 1 `_prepare_caches!` as `apply_K!` /
+    # `assemble_M!` so the geometry / N_data / detJ_w are populated
+    # exactly the same way. Cheap to call repeatedly because Pass 1
+    # is a fixed-cost sweep over the elements.
+    _prepare_caches!(cache, kernel, mesh)
+
+    elements         = cache.elements
+    element_caches   = cache.element_caches
+    geometry_caches  = cache.geometry_caches
+
+    layout     = local_dof_layout(E)
+    ndofs_elem = length(layout)
+
+    @inbounds for elem_idx in 1:length(elements)
+        ec = element_caches[elem_idx]
+        gc = geometry_caches[elem_idx]
+        n_ips     = length(gc.detJ_w)
+        dofs_elem = ec.dofs
+
+        @inbounds for li in 1:ndofs_elem
+            entry  = layout[li]
+            node_i = entity_local(entry)
+            comp_i = component(entry)
+            bα     = _body_component(load.b, comp_i)
+            if bα == 0.0
+                continue
+            end
+
+            sum_q = 0.0
+            @inbounds for q in 1:n_ips
+                N_i   = gc.N_data[q, node_i]
+                detJw = gc.detJ_w[q]
+                sum_q += N_i * detJw
+            end
+
+            f[Int(dofs_elem[li])] += bα * sum_q
+        end
+    end
+
+    return f
+end
+
+# ----------------------------------------------------------------------------
+# UniformMixedDarcySource — ∫ q f dΩ on P₀ pressure test functions
+# ----------------------------------------------------------------------------
+
+@inline function _mixed_darcy_flux_li(layout, iface::Int)
+    @inbounds for li in 1:length(layout)
+        e = layout[li]
+        if Int(field_idx(e)) == 1 && Int(entity_local(e)) == iface && Int(component(e)) == 1
+            return li
+        end
+    end
+    return 0
+end
+
+"""
+    UniformMixedDarcySource(f)
+
+Uniform volumetric source for the **pressure test equation** of mixed RT₀–P₀ Darcy:
+
+``\\int_\\Omega q \\,(\\nabla\\!\\cdot u)\\,\\mathrm{d}\\Omega = \\int_\\Omega q \\,f\\,\\mathrm{d}\\Omega``
+
+with piecewise constant `q` (one DOF per cell). Each cell pressure unknown receives
+``f \\cdot |K_e|`` where `|K_e|` is the element volume from ``\\sum_q \\det J \\, w``.
+
+Implemented for [`AbstractDarcyMixedRT0P0Kernel`](@ref) (Tet4 and Hex8 mixed kernels).
+The element template must place the scalar cell pressure *after* face flux fields in
+[`local_dof_layout`](@ref) (the constructor discovers the pressure local DOF by `field_idx == 2`).
+
+Calls `_prepare_caches!` once (same as [`UniformBodyForce`](@ref)).
+"""
+struct UniformMixedDarcySource <: AbstractNeumannLoad
+    f::Float64
+end
+
+function apply_load!(
+    fvec::AbstractVector{Float64},
+    load::UniformMixedDarcySource,
+    cache::DOFBasedCOOCache{T, B, IPS, E, GC, Buf, FT, ST},
+    asm::DOFBasedCOOAssembler,
+    kernel::AbstractDarcyMixedRT0P0Kernel,
+    mesh::AbstractMesh,
+) where {T, B, IPS, E <: AbstractElement, GC, Buf, FT, ST}
+    @assert length(fvec) == cache.ndofs (
+        "apply_load!: f has length $(length(fvec)); expected $(cache.ndofs)")
+
+    _prepare_caches!(cache, kernel, mesh)
+
+    layout = local_dof_layout(E)
+    p_li = 0
+    @inbounds for li in 1:length(layout)
+        if Int(field_idx(layout[li])) == 2
+            p_li = li
+            break
+        end
+    end
+    p_li > 0 || error(
+        "UniformMixedDarcySource: no local DOF with field_idx == 2 in Element{$E}; " *
+            "use @DOFSet{σ::DOF{RT0FaceFlux, Face}, p::DOF{Float64, Cell}} with σ first.",
+    )
+
+    element_caches = cache.element_caches
+    geometry_caches = cache.geometry_caches
+    nelem = length(element_caches)
+
+    src = load.f
+    @inbounds for eid in 1:nelem
+        gc = geometry_caches[eid]
+        ec = element_caches[eid]
+        vol = 0.0
+        n_ips = length(gc.detJ_w)
+        for q in 1:n_ips
+            vol += gc.detJ_w[q]
+        end
+        gdof = Int(ec.dofs[p_li])
+        fvec[gdof] += src * vol
+    end
+
+    return fvec
+end
+
+# ----------------------------------------------------------------------------
+# MixedDarcyTet4BoundaryNormalFluxLoad — ∫ g φ·n dS on RT₀ flux test functions (Tet4)
+# ----------------------------------------------------------------------------
+
+# Symmetric order-2 rule on the reference triangle ξ ≥ 0, η ≥ 0, ξ + η ≤ 1 (area 1/2).
+const _MIXED_DARCY_TRI3_AB =
+    ((1.0 / 6.0, 1.0 / 6.0), (2.0 / 3.0, 1.0 / 6.0), (1.0 / 6.0, 2.0 / 3.0))
+const _MIXED_DARCY_TRI3_W = (1.0 / 6.0, 1.0 / 6.0, 1.0 / 6.0)
+
+"""
+    MixedDarcyTet4BoundaryNormalFluxLoad(panels, g)
+
+Boundary contribution ``\\int_\\Gamma g\\, \\mathbf{\\phi}_i \\cdot \\mathbf{n}\\,\\mathrm{d}S`` to the
+**flux test equation** of mixed RT₀–P₀ Darcy on Tet4, with uniform normal flux density ``g``
+(`[m/s]`, outward positive relative to the element). Each entry of `panels` is
+`(elem_id, local_face)` with `local_face ∈ 1:4` ([`faces(::Tet4)`](@ref)).
+
+Uses three quadrature points per triangle (exact if ``g \\,\\mathbf{\\phi}\\!\\cdot\\!\\mathbf{n}`` is linear on the face).
+Requires [`DarcyMixedRT0P0Kernel`](@ref) and [`Mesh{4, Tet4}`](@ref).
+"""
+struct MixedDarcyTet4BoundaryNormalFluxLoad <: AbstractNeumannLoad
+    panels::Vector{Tuple{Int, Int}}
+    g::Float64
+end
+
+function apply_load!(
+    fvec::AbstractVector{Float64},
+    load::MixedDarcyTet4BoundaryNormalFluxLoad,
+    cache::DOFBasedCOOCache{T, B, IPS, E, GC, Buf, FT, ST},
+    asm::DOFBasedCOOAssembler,
+    kernel::DarcyMixedRT0P0Kernel,
+    mesh::Mesh{4, Tet4},
+) where {T, B, IPS, E <: AbstractElement, GC, Buf, FT, ST}
+    @assert length(fvec) == cache.ndofs (
+        "apply_load!: f has length $(length(fvec)); expected $(cache.ndofs)")
+
+    _prepare_caches!(cache, kernel, mesh)
+
+    layout = local_dof_layout(E)
+    element_caches = cache.element_caches
+    geometry_caches = cache.geometry_caches
+    nelem = length(element_caches)
+
+    @inbounds for pid in eachindex(load.panels)
+        eid, lf = load.panels[pid]
+        (1 ≤ eid ≤ nelem) || throw(ArgumentError("MixedDarcyTet4BoundaryNormalFluxLoad: elem_id $eid out of range 1:$nelem"))
+        (1 ≤ lf ≤ 4) || throw(ArgumentError("MixedDarcyTet4BoundaryNormalFluxLoad: local_face $lf out of range 1:4"))
+
+        gc = geometry_caches[eid]
+        ec = element_caches[eid]
+        X = gc.X
+
+        fc = faces(Tet4())[lf]
+        @inbounds p1 = X[Int(fc.vertices[1])]
+        @inbounds p2 = X[Int(fc.vertices[2])]
+        @inbounds p3 = X[Int(fc.vertices[3])]
+
+        Vphys = 0.0
+        n_ips = length(gc.detJ_w)
+        for q in 1:n_ips
+            Vphys += gc.detJ_w[q]
+        end
+
+        cross_vec = (p2 - p1) × (p3 - p1)
+        jac_face = norm(cross_vec)
+        jac_face > 0.0 || continue
+
+        orient = Float64(tet_face_outward_sign(X, lf))
+        n_unit = orient * cross_vec / jac_face
+
+        for k in 1:3
+            ξq, ηq = _MIXED_DARCY_TRI3_AB[k]
+            wq = _MIXED_DARCY_TRI3_W[k]
+            λ1 = 1.0 - ξq - ηq
+            xq = λ1 * p1 + ξq * p2 + ηq * p3
+            base = load.g * wq * jac_face
+            for ifi in 1:4
+                φ = _rt0_phi_tet4(X, Vphys, ifi, xq)
+                li = _mixed_darcy_flux_li(layout, ifi)
+                li > 0 || error("MixedDarcyTet4BoundaryNormalFluxLoad: no flux DOF for local face $ifi")
+                gdof = Int(ec.dofs[li])
+                fvec[gdof] += base * dot(φ, n_unit)
+            end
+        end
+    end
+
+    return fvec
+end
+
+# ----------------------------------------------------------------------------
+# SurfaceLoad — distributed traction / heat flux integrated on faces
+# ----------------------------------------------------------------------------
+
+"""
+    SurfaceLoad(faces, traction)
+
+Distributed surface load (traction `t`, heat flux `q`, …) integrated as
+`∫_Γ N_i · t dS` over a list of mesh faces. The natural complement of
+`UniformBodyForce`: body force lives in the volume integral, surface
+traction lives in the surface integral.
+
+# Arguments
+
+- `faces::Vector{NTuple{NN,Int}}` — each face is the tuple of *global*
+  mesh-node IDs of its corners. `NN = 4` for a quadrilateral face
+  (e.g. one face of a Hex8 element); `NN = 3` for a triangular face
+  (one face of a Tet4) — both are supported.
+
+- `traction` — one of
+    * `Vec{3,Float64}` — uniform vector traction on every face (3D
+      elasticity), in units of force per unit area.
+    * `Float64`        — uniform scalar flux on every face (heat
+      surface flux `[W/m²]`, or normal **Darcy flux** `[m/s]` for primal
+      potential — same `∫_Γ N_i q · dS` assembly via [`SurfaceLoad`](@ref)).
+    * `Vector{Vec{3,Float64}}` — per-face vector traction.
+    * `Vector{Float64}` — per-face scalar flux.
+
+# Quadrature
+
+* Quadrilateral face (`NN == 4`) → 2 × 2 Gauss (4 points), exact for
+  bilinear basis on a flat quad with constant traction.
+* Triangular face (`NN == 3`) → 1-point centroid (exact for linear
+  basis with constant traction; bumped to 3-point if needed by a
+  later non-flat / higher-order extension).
+
+# Surface Jacobian
+
+For a face `Γ` with corner coordinates `xᵢ`, the surface measure is
+
+    dS = ‖∂x/∂ξ × ∂x/∂η‖ dξ dη
+
+(quad) or `‖e₁ × e₂‖ / 2` (tri, where `eₖ` are the in-plane edge
+vectors). Both forms are reduced to a single `_face_metric` call below
+so the per-face loop stays type-stable.
+
+# Notes
+
+* `apply_load!` for `SurfaceLoad` does not call `_prepare_caches!`
+  — it works directly off `mesh.nodes` for the face geometry and the
+  `DOFHandler.field_starts` table for DOF lookup. There is no need
+  to materialise face elements, build a face cache, or call any
+  volume-element machinery. This keeps the surface path compositional
+  (any subset of faces from any mesh works) and disjoint from the
+  ContinuumKernel evaluation path.
+* The implementation is allocation-free after warmup. The face-corner
+  coordinate buffer is stack-allocated (`SVector` for `NN ≤ 8`).
+"""
+struct SurfaceLoad{NN, V} <: AbstractNeumannLoad
+    faces::Vector{NTuple{NN,Int}}
+    traction::V
+
+    function SurfaceLoad(faces::Vector{NTuple{NN,Int}}, traction::V) where {NN, V}
+        # `Vec{3,Float64}` is itself an `AbstractVector` (it subtypes
+        # `StaticArray`), so distinguish *per-face* tractions
+        # (heap-stored `Vector{...}`) from a uniform tensor traction
+        # by the concrete `Vector` type, not by `AbstractVector`.
+        if traction isa Vector{<:Vec{3,Float64}} || traction isa Vector{Float64}
+            @assert length(traction) == length(faces) (
+                "SurfaceLoad: per-face traction has $(length(traction)) entries " *
+                "but $(length(faces)) faces were given")
+        end
+        if !(NN == 3 || NN == 4)
+            error("SurfaceLoad: only NN=3 (triangular) and NN=4 (quadrilateral) " *
+                  "faces are supported (got NN=$NN).")
+        end
+        return new{NN, V}(faces, traction)
+    end
+end
+
+# Per-face traction accessor — dispatches on the traction-storage type
+# itself: a single value is uniform across all faces, a `Vector{...}`
+# is one entry per face. Both vector- (elasticity) and scalar- (heat)
+# valued tractions are supported by the same two methods.
+@inline _face_traction(t::Vec{3,Float64}, f::Int) = t
+@inline _face_traction(t::Float64,        f::Int) = t
+@inline _face_traction(t::Vector{<:Vec{3,Float64}}, f::Int) = @inbounds t[f]
+@inline _face_traction(t::Vector{Float64},          f::Int) = @inbounds t[f]
+
+@inline _t_component(t::Vec{3,Float64}, comp::Int) = t[comp]
+@inline _t_component(t::Float64,        comp::Int) = t
+
+# Component count of a single-face traction value (not the whole load).
+@inline _t_comp_count(::Vec{3,Float64}) = 3
+@inline _t_comp_count(::Float64)        = 1
+
+# ----- face geometry helpers ------------------------------------------------
+
+# Quad4 reference basis at (ξ, η) ∈ [-1, 1]² and its derivatives.
+@inline function _quad4_basis(ξ::Float64, η::Float64)
+    n1 = 0.25 * (1 - ξ) * (1 - η)
+    n2 = 0.25 * (1 + ξ) * (1 - η)
+    n3 = 0.25 * (1 + ξ) * (1 + η)
+    n4 = 0.25 * (1 - ξ) * (1 + η)
+    return (n1, n2, n3, n4)
+end
+
+@inline function _quad4_basis_derivs(ξ::Float64, η::Float64)
+    # ∂N/∂ξ
+    dξ = (-0.25 * (1 - η),  0.25 * (1 - η),  0.25 * (1 + η), -0.25 * (1 + η))
+    # ∂N/∂η
+    dη = (-0.25 * (1 - ξ), -0.25 * (1 + ξ),  0.25 * (1 + ξ),  0.25 * (1 - ξ))
+    return dξ, dη
+end
+
+# Tri3 reference basis at barycentric (L1, L2) on the unit triangle
+# ξ ∈ [0,1], η ∈ [0,1-ξ].
+@inline _tri3_basis(ξ::Float64, η::Float64) = (1.0 - ξ - η, ξ, η)
+@inline _tri3_basis_derivs(::Float64, ::Float64) = ((-1.0, 1.0, 0.0), (-1.0, 0.0, 1.0))
+
+# 2 × 2 Gauss for quad: 4 points, weight 1 each.
+const _GAUSS_2X2 = ((-1.0/√3, -1.0/√3, 1.0),
+                    ( 1.0/√3, -1.0/√3, 1.0),
+                    ( 1.0/√3,  1.0/√3, 1.0),
+                    (-1.0/√3,  1.0/√3, 1.0))
+
+# ----------------------------------------------------------------------------
+# MixedDarcyHex8BoundaryNormalFluxLoad — ∫ g φ·n dS on RT₀ flux test functions (Hex8)
+# ----------------------------------------------------------------------------
+
+@inline function _hex8_face_volume_ref_coords(ξf::Float64, ηf::Float64, lf::Int)
+    fc = faces(Hex8())[lf]
+    vs = fc.vertices
+    ref_c = reference_coordinates(Hex8())
+    N = _quad4_basis(ξf, ηf)
+    @inbounds return N[1] * ref_c[vs[1]] +
+        N[2] * ref_c[vs[2]] +
+        N[3] * ref_c[vs[3]] +
+        N[4] * ref_c[vs[4]]
+end
+
+"""
+    MixedDarcyHex8BoundaryNormalFluxLoad(panels, g)
+
+Boundary contribution ``\\int_\\Gamma g\\, \\mathbf{\\phi}_i \\cdot \\mathbf{n}\\,\\mathrm{d}S`` to the
+**flux test equation** of mixed RT₀–P₀ Darcy on Hex8, with uniform normal flux density ``g``
+(`[m/s]`, outward positive relative to the element). Each entry of `panels` is
+`(elem_id, local_face)` with `local_face ∈ 1:6` ([`faces(::Hex8)`](@ref)).
+
+Contravariant Piola push-forward of [`rt0_hex8_reference_basis`](@ref) matches
+[`DarcyMixedHex8RT0P0Kernel`](@ref). Uses four Gauss points per quadrilateral face.
+
+Requires [`DarcyMixedHex8RT0P0Kernel`](@ref) and [`Mesh{8, Hex8}`](@ref).
+
+Note: entries scale with ``\\int_{\\partial\\Omega} \\mathbf{\\phi}_i\\!\\cdot\\!\\mathbf{n}\\,\\mathrm{d}S`` (Piola field),
+which for a single trilinear brick need not match the lumped ``\\pm 1`` divergence coupling used in
+[`DarcyMixedHex8RT0P0Kernel`](@ref) (e.g. ``4`` on one ``[0,1]^3`` element from ``[-1,1]^3`` reference).
+"""
+struct MixedDarcyHex8BoundaryNormalFluxLoad <: AbstractNeumannLoad
+    panels::Vector{Tuple{Int, Int}}
+    g::Float64
+end
+
+function apply_load!(
+    fvec::AbstractVector{Float64},
+    load::MixedDarcyHex8BoundaryNormalFluxLoad,
+    cache::DOFBasedCOOCache{T, B, IPS, E, GC, Buf, FT, ST},
+    asm::DOFBasedCOOAssembler,
+    kernel::DarcyMixedHex8RT0P0Kernel,
+    mesh::Mesh{8, Hex8},
+) where {T, B, IPS, E <: AbstractElement, GC, Buf, FT, ST}
+    @assert length(fvec) == cache.ndofs (
+        "apply_load!: f has length $(length(fvec)); expected $(cache.ndofs)")
+
+    _prepare_caches!(cache, kernel, mesh)
+
+    layout = local_dof_layout(E)
+    element_caches = cache.element_caches
+    geometry_caches = cache.geometry_caches
+    nelem = length(element_caches)
+
+    @inbounds for pid in eachindex(load.panels)
+        eid, lf = load.panels[pid]
+        (1 ≤ eid ≤ nelem) || throw(ArgumentError("MixedDarcyHex8BoundaryNormalFluxLoad: elem_id $eid out of range 1:$nelem"))
+        (1 ≤ lf ≤ 6) || throw(ArgumentError("MixedDarcyHex8BoundaryNormalFluxLoad: local_face $lf out of range 1:6"))
+
+        gc = geometry_caches[eid]
+        ec = element_caches[eid]
+        X = gc.X
+        length(X) == 8 || continue
+
+        fc = faces(Hex8())[lf]
+        vs = fc.vertices
+        @inbounds p1 = X[Int(vs[1])]
+        @inbounds p2 = X[Int(vs[2])]
+        @inbounds p3 = X[Int(vs[3])]
+        @inbounds p4 = X[Int(vs[4])]
+
+        orient = Float64(hex8_face_outward_sign(X, lf))
+
+        @inbounds for gpt in _GAUSS_2X2
+            ξf = gpt[1]
+            ηf = gpt[2]
+            wf = gpt[3]
+            ξv = _hex8_face_volume_ref_coords(ξf, ηf, lf)
+            dN_dξ = get_basis_derivatives(Hex8(), Lagrange{1}(), ξv)
+            J = X[1] ⊗ dN_dξ[1]
+            for a in 2:8
+                J += X[a] ⊗ dN_dξ[a]
+            end
+            detJ = det(J)
+            detJ == 0.0 && continue
+
+            dξq, dηq = _quad4_basis_derivs(ξf, ηf)
+            tξ = dξq[1] * p1 + dξq[2] * p2 + dξq[3] * p3 + dξq[4] * p4
+            tη = dηq[1] * p1 + dηq[2] * p2 + dηq[3] * p3 + dηq[4] * p4
+            jac_met = norm(tξ × tη)
+            jac_met == 0.0 && continue
+            # Normal direction must match [`hex8_face_outward_sign`](@ref) (first triangle on
+            # `faces(::Hex8)`); ‖tξ×tη‖ is still the correct surface Jacobian for the quad map.
+            cross_face = (p2 - p1) × (p3 - p1)
+            nf = norm(cross_face)
+            nf == 0.0 && continue
+            n_unit = orient * cross_face / nf
+            dS_w = jac_met * wf
+
+            base = load.g * dS_w
+            @inbounds for ifi in 1:6
+                ψ = rt0_hex8_reference_basis(ifi, ξv)
+                φ = piola_contravariant(J, ψ)
+                li = _mixed_darcy_flux_li(layout, ifi)
+                li > 0 || error("MixedDarcyHex8BoundaryNormalFluxLoad: no flux DOF for local face $ifi")
+                gdof = Int(ec.dofs[li])
+                fvec[gdof] += base * dot(φ, n_unit)
+            end
+        end
+    end
+
+    return fvec
+end
+
+# 1-point centroid for tri: weight = 1/2 (unit triangle area).
+const _GAUSS_TRI1 = ((1.0/3.0, 1.0/3.0, 0.5),)
+
+# ----- core face integration ------------------------------------------------
+
+# Look up the global DOF index for (node_id, component) given the
+# *concrete* field-starts vector for field 1. Hoisted outside the hot
+# loops so we only pay one type-assertion per `apply_load!` call,
+# never per face / per component.
+@inline _node_field_dof(field_starts1::Vector{Int}, node_id::Int, comp::Int) =
+    @inbounds field_starts1[node_id] + (comp - 1)
+
+# One-face integration: writes f[d] += traction · N_i · dS · w.
+# Specialised on NN (3 or 4) so the inner loops fully unroll. Takes the
+# pre-resolved `field_starts1::Vector{Int}` so the inner DOF lookup is
+# fully concrete-typed (no allocations regardless of the cache's
+# `DOFHandler` abstract field).
+@inline function _integrate_face!(f::AbstractVector{Float64},
+                                  field_starts1::Vector{Int},
+                                  nodes_xyz,
+                                  face_nodes::NTuple{NN,Int},
+                                  t,
+                                  ::Val{NN}) where {NN}
+    n_comp = _t_comp_count(t)
+
+    # Pull face corner coordinates onto the stack as a small NTuple.
+    # Avoids any heap allocation in the hot loop.
+    X = ntuple(i -> nodes_xyz[face_nodes[i]], NN)
+
+    if NN == 4
+        @inbounds for (ξ, η, w) in _GAUSS_2X2
+            N      = _quad4_basis(ξ, η)
+            dξ, dη = _quad4_basis_derivs(ξ, η)
+            tξ = dξ[1] * X[1] + dξ[2] * X[2] + dξ[3] * X[3] + dξ[4] * X[4]
+            tη = dη[1] * X[1] + dη[2] * X[2] + dη[3] * X[3] + dη[4] * X[4]
+            dS_w = norm(tξ × tη) * w
+            for li in 1:NN
+                Ni = N[li] * dS_w
+                node = face_nodes[li]
+                for comp in 1:n_comp
+                    d = _node_field_dof(field_starts1, node, comp)
+                    f[d] += _t_component(t, comp) * Ni
+                end
+            end
+        end
+    else  # NN == 3
+        @inbounds for (ξ, η, w) in _GAUSS_TRI1
+            N      = _tri3_basis(ξ, η)
+            dξ, dη = _tri3_basis_derivs(ξ, η)
+            tξ = dξ[1] * X[1] + dξ[2] * X[2] + dξ[3] * X[3]
+            tη = dη[1] * X[1] + dη[2] * X[2] + dη[3] * X[3]
+            dS_w = norm(tξ × tη) * w
+            for li in 1:NN
+                Ni = N[li] * dS_w
+                node = face_nodes[li]
+                for comp in 1:n_comp
+                    d = _node_field_dof(field_starts1, node, comp)
+                    f[d] += _t_component(t, comp) * Ni
+                end
+            end
+        end
+    end
+    return nothing
+end
+
+"""
+    apply_load!(f::AbstractVector{Float64}, load::SurfaceLoad,
+                cache::DOFBasedCOOCache, asm, kernel, mesh) -> f
+
+Surface-traction / heat-flux integration: `f[i, α] += ∫_Γ N_i · t_α dS`
+for every face in `load.faces`. Uses the cache's `dof_connectivity`
+parent `DOFHandler` to map `(face_node, component)` to global DOF
+indices, and the mesh's nodal coordinates to compute the surface
+Jacobian on the fly. Allocation-free after warmup.
+"""
+function apply_load!(f::AbstractVector{Float64},
+                     load::SurfaceLoad{NN,V},
+                     cache::DOFBasedCOOCache,
+                     asm::DOFBasedCOOAssembler,
+                     kernel::AbstractKernel,
+                     mesh::AbstractMesh) where {NN,V}
+    @assert length(f) == cache.ndofs (
+        "apply_load!: f has length $(length(f)); expected $(cache.ndofs)")
+
+    # The cache pre-resolves `dof_handler.field_starts[1]` into a
+    # concrete `Vector{Int}` field at construction so the hot loop
+    # below stays allocation-free even though `cache.dof_handler` is
+    # stored at the abstract `DOFHandler` type.
+    field_starts1 = cache.field_starts1
+    nodes_xyz = mesh.nodes
+    nfaces = length(load.faces)
+
+    @inbounds for fi in 1:nfaces
+        face_nodes = load.faces[fi]
+        t = _face_traction(load.traction, fi)
+        _integrate_face!(f, field_starts1, nodes_xyz, face_nodes, t, Val(NN))
+    end
+
+    return f
+end