diff --git a/src/domains/continuum/abstract.jl b/src/domains/continuum/abstract.jl
index 4a84903..0f18cb8 100644
--- a/src/domains/continuum/abstract.jl
+++ b/src/domains/continuum/abstract.jl
@@ -1,5 +1,5 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+# SPDX-FileCopyrightText: 2015-2026 Jukka Aho
+# SPDX-License-Identifier: MIT
 
 """
 Abstract types for continuum mechanics domain.
@@ -24,7 +24,7 @@ DESIGN PHILOSOPHY: Formulations are DOMAIN-AGNOSTIC dimensionality concepts.
 Formulation (domain-agnostic):
 - Describes DIMENSIONALITY and geometric simplifications
 - Used by multiple physics domains
-- Examples: FullThreeD, Axisymmetric
+- Examples: ThreeDimensional, Axisymmetric
 - Can be reused across continuum, heat, acoustics, etc.
 
 Theory (domain-specific):
@@ -34,8 +34,8 @@ Theory (domain-specific):
 
 # Why Separate Them?
 
-Problem: Heat transfer needs FullThreeD and Axisymmetric, just like continuum!
-If FullThreeD is defined in domains/continuum/, heat can't use it without duplication.
+Problem: Heat transfer needs ThreeDimensional and Axisymmetric, just like continuum!
+If ThreeDimensional is defined in domains/continuum/, heat can't use it without duplication.
 
 Solution: Formulations are dimensionality (shared), theories are physics (domain-specific).
 
@@ -43,7 +43,7 @@ Solution: Formulations are dimensionality (shared), theories are physics (domain
 
 ```julia
 # Domain-agnostic formulations
-FullThreeD()              # Used by: continuum, heat, poisson, acoustics
+ThreeDimensional()              # Used by: continuum, heat, poisson, acoustics
 Axisymmetric()            # Used by: continuum, heat, etc.
 
 # Domain-specific theories (in domains/*/types.jl)
@@ -75,7 +75,7 @@ PlaneStrain (ε_xx, ε_yy, ε_xy, ε_zz = 0):
 - Out-of-plane strain ε_zz = 0
 - Examples: Dams, tunnels, retaining walls, long cylinders
 
-FullThreeD:
+ThreeDimensional:
 - No simplifications, all six stress/strain components
 - Most accurate but most expensive
 
@@ -97,6 +97,6 @@ Plane Strain (thick section):
 - Constitutive: 3×3 reduced stiffness matrix (different from plane stress!)
 
 # See Also
-- Concrete theories: `continuum/types.jl` (FullThreeD, PlaneStress, PlaneStrain, Axisymmetric)
+- Concrete theories: `continuum/types.jl` (ThreeDimensional, PlaneStress, PlaneStrain, Axisymmetric)
 """
 abstract type AbstractContinuumTheory end
diff --git a/src/domains/continuum/dof_based_pass1.jl b/src/domains/continuum/dof_based_pass1.jl
new file mode 100644
index 0000000..ec6cc04
--- /dev/null
+++ b/src/domains/continuum/dof_based_pass1.jl
@@ -0,0 +1,172 @@
+# SPDX-FileCopyrightText: 2015-2026 Jukka Aho
+# SPDX-License-Identifier: MIT
+
+"""
+Pass~1 helpers for DOF-based assembly with [`ContinuumKernel`](@ref):
+vertex displacement scatter and material-workspace preparation.
+"""
+
+using Tensors
+using Tensors: basevec, symmetric, Tensor
+
+using ..JuliaFEM: AbstractMaterial, compute_stress, continuum_kinematics,
+    GreenLagrangeKinematics, SmallStrainKinematics
+
+"""
+    scatter_vertex_displacements_from_global!(
+        u_buffer, dofs_storage, configuration, ::Type{E},
+    ) -> Nothing
+
+Scatter global displacement DOFs from `configuration` into
+`u_buffer` (vertex-major `Vec{3}` per topology node) using compile-time
+[`local_dof_layout`](@ref)`(E)` and the element's global DOF list in
+`dofs_storage`. Only entries with `1 ≤ component ≤ 3` and valid
+`entity_local` vertex indices participate.
+
+No heap allocation in the loop.
+"""
+@inline function scatter_vertex_displacements_from_global!(
+    u_buffer::Vector{Vec{3,Float64}},
+    dofs_storage::Vector{Int},
+    configuration::AbstractVector{Float64},
+    ::Type{E},
+) where {E<:AbstractElement}
+    layout = local_dof_layout(E)
+    nbuf = length(u_buffer)
+    @inbounds for v in 1:nbuf
+        u_buffer[v] = zero(Vec{3,Float64})
+    end
+    @inbounds for li in eachindex(layout)
+        ent = entity_local(layout[li])
+        (1 ≤ ent ≤ nbuf) || continue
+        comp = component(layout[li])
+        (1 ≤ comp ≤ 3) || continue
+        d = dofs_storage[li]
+        u_buffer[ent] = u_buffer[ent] + basevec(Vec{3,Float64}, comp) * configuration[d]
+    end
+    return nothing
+end
+
+"""
+    _continuum_fill_workspace_stress_from_displacement!(
+        material_workspace, geometry_cache, element_cache, material, Δt, empty_state,
+    ) -> Nothing
+
+Recompute `(σ, 𝔻)` at every IP from nodal displacements in `element_cache.u_buffer`
+using [`continuum_kinematics`](@ref)`(material)` (small strain or Green–Lagrange).
+
+Used when `material_behavior(material) isa StatelessConstantTangent` but a global
+`configuration` is supplied so Cauchy stress tracks the current displacement while
+the tangent remains the constitutive tangent returned by [`compute_stress`](@ref).
+
+Zero-allocation in the IP loop.
+"""
+@inline function _continuum_fill_workspace_stress_from_displacement!(
+    material_workspace::AssemblyMaterialWorkspace{FieldType, StateType},
+    geometry_cache::GeometryCache,
+    element_cache::ElementCache,
+    material::AbstractMaterial,
+    Δt::Float64,
+    empty_state::NamedTuple,
+) where {FieldType, StateType}
+    fields_mw = getfield(material_workspace, 1)
+    states_mw = getfield(material_workspace, 2)
+    nips = length(element_cache.ips)
+    nnodes = length(geometry_cache.X)
+    kin = continuum_kinematics(material)
+    I = one(Tensor{2,3,Float64,9})
+    @inbounds for q in 1:nips
+        strain_measure = if kin isa SmallStrainKinematics
+            ε = zero(SymmetricTensor{2,3,Float64,6})
+            for k in 1:nnodes
+                u_k = element_cache.u_buffer[k]
+                ∇N_k_q = geometry_cache.∇N_data[q, k]
+                ε += symmetric(u_k ⊗ ∇N_k_q)
+            end
+            ε
+        elseif kin isa GreenLagrangeKinematics
+            F = I
+            for k in 1:nnodes
+                u_k = element_cache.u_buffer[k]
+                ∇N_k_q = geometry_cache.∇N_data[q, k]
+                F += u_k ⊗ ∇N_k_q
+            end
+            C_tensor = symmetric(F' ⋅ F)
+            SymmetricTensor{2,3}(0.5 * (C_tensor - I))
+        else
+            throw(ArgumentError("unknown continuum kinematics $(typeof(kin))"))
+        end
+        σ, 𝔻, _ = compute_stress(material, strain_measure, NamedTuple(), Δt)
+        fields_mw[q] = (σ=σ, 𝔻=𝔻)
+        states_mw[q] = empty_state
+    end
+    return nothing
+end
+
+@inline function prepare_dof_based_material_workspace!(
+    k_e::ContinuumKernel,
+    material_workspace::AssemblyMaterialWorkspace,
+    geometry_cache::GeometryCache,
+    element_cache::ElementCache,
+    eid::Int,
+    configuration::Union{Nothing,AbstractVector{Float64}},
+    global_material_cache::Union{Nothing,GlobalMaterialCache},
+    Δt::Float64,
+    ::Type{E},
+) where {E<:AbstractElement}
+    mat = k_e.material
+    beh = material_behavior(mat)
+    fields_mw = getfield(material_workspace, 1)
+    states_mw = getfield(material_workspace, 2)
+    ips_ec = getfield(element_cache, :ips)
+    nips = length(ips_ec)
+    if beh isa StatelessConstantTangent
+        fields_ref_e, empty_state_e = reference_fields(k_e)
+        if configuration !== nothing
+            scatter_vertex_displacements_from_global!(
+                element_cache.u_buffer, element_cache.dofs, configuration, E,
+            )
+            _continuum_fill_workspace_stress_from_displacement!(
+                material_workspace, geometry_cache, element_cache, mat, Δt, empty_state_e,
+            )
+        else
+            @inbounds for q in 1:nips
+                fields_mw[q] = fields_ref_e
+                states_mw[q] = empty_state_e
+            end
+        end
+    elseif beh isa StatelessStrainDependent
+        if configuration !== nothing
+            scatter_vertex_displacements_from_global!(
+                element_cache.u_buffer, element_cache.dofs, configuration, E,
+            )
+        end
+        update_material_cache_stateless_strain!(
+            material_workspace, geometry_cache, mat, element_cache, Δt,
+        )
+    elseif beh isa StatefulStrainDependent
+        global_material_cache === nothing && throw(ArgumentError(
+            "DOF-based Pass 1: StatefulStrainDependent material requires keyword " *
+            "`global_material_cache=create_global_material_cache(mat; n_ips, n_elems)`",
+        ))
+        if configuration !== nothing
+            scatter_vertex_displacements_from_global!(
+                element_cache.u_buffer, element_cache.dofs, configuration, E,
+            )
+        end
+        update_material_cache!(
+            material_workspace,
+            geometry_cache,
+            mat,
+            element_cache,
+            global_material_cache,
+            eid,
+            Δt,
+        )
+    else
+        throw(ArgumentError(
+            "unsupported material behavior $(typeof(beh)) for ContinuumKernel in DOF-based Pass 1",
+        ))
+    end
+    return nothing
+end
diff --git a/src/domains/continuum/hellinger_reissner_kernel.jl b/src/domains/continuum/hellinger_reissner_kernel.jl
index 0519812..664d172 100644
--- a/src/domains/continuum/hellinger_reissner_kernel.jl
+++ b/src/domains/continuum/hellinger_reissner_kernel.jl
@@ -78,7 +78,7 @@ the discrete `σ`–`σ` block uses `G M⁻¹ G` in the Voigt component basis.
 ```julia
 S = @DOFSet{u::DOF{Displacement{3}, Vertex}, σ::DOF{SymmetricTensor{2,3}, Cell}}
 kernel = HellingerReissnerKernel(
-    ContinuumFormulation{FullThreeD}(),
+    ContinuumFormulation{ThreeDimensional}(),
     LinearElastic(E = 210e9, ν = 0.3),
 )
 ```
diff --git a/src/domains/continuum/hu_washizu_kernel.jl b/src/domains/continuum/hu_washizu_kernel.jl
index cb94edd..255ec49 100644
--- a/src/domains/continuum/hu_washizu_kernel.jl
+++ b/src/domains/continuum/hu_washizu_kernel.jl
@@ -101,7 +101,7 @@ S = @DOFSet{
     sig::DOF{SymmetricTensor{2,3}, Cell},
 }
 kernel = HuWashizuKernel(
-    ContinuumFormulation{FullThreeD}(),
+    ContinuumFormulation{ThreeDimensional}(),
     LinearElastic(E = 210e9, ν = 0.3),
 )
 ```
diff --git a/src/domains/continuum/kernel.jl b/src/domains/continuum/kernel.jl
index ab2cb85..593c340 100644
--- a/src/domains/continuum/kernel.jl
+++ b/src/domains/continuum/kernel.jl
@@ -1,5 +1,5 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+# SPDX-FileCopyrightText: 2015-2026 Jukka Aho
+# SPDX-License-Identifier: MIT
 
 """
 Continuum mechanics kernel - defines the weak form only.
@@ -10,9 +10,10 @@ This module defines:
 3. Block builder: compute_stiffness_block (builds D×D blocks)
 
 The DOF-based / matrix-free assembler microkernel surface
-(`qpoint_buffer_eltype`, `update_qpoint_buffer!`, `evaluate_entry`,
-`evaluate_mass_entry`, `reference_fields`) is implemented further
-down in this file. Everything else (geometry preprocessing,
+(`qpoint_buffer_eltype`, `prepare_dof_based_material_workspace!`,
+`update_qpoint_buffer!`, `evaluate_entry`, `evaluate_mass_entry`,
+`reference_fields`) is implemented in this file and in
+`dof_based_pass1.jl`. Everything else (geometry preprocessing,
 integration, assembly, DOF mapping) belongs elsewhere.
 """
 
@@ -29,7 +30,7 @@ optional density (carried on the kernel rather than the material so
 existing material structs stay untouched).
 
 # Type Parameters
-- `Theory`: Continuum theory (FullThreeD, PlaneStress, PlaneStrain, Axisymmetric)
+- `Theory`: Continuum theory (ThreeDimensional, PlaneStress, PlaneStrain, Axisymmetric)
 - `Mat`: Material model (LinearElastic, NeoHookean, etc.)
 
 # Fields
@@ -45,7 +46,7 @@ existing material structs stay untouched).
 
 ```julia
 kernel = ContinuumKernel(
-    ContinuumFormulation{FullThreeD}(),
+    ContinuumFormulation{ThreeDimensional}(),
     LinearElastic(E=210e9, ν=0.3),
     Displacement{3}();
     density = 7850.0,        # for mass matrix; omit for static-only
@@ -223,6 +224,35 @@ keeps the entire chain inside the symmetric-tensor methods of
     return dcontract(B_k_α, dcontract(C, B_l_β))
 end
 
+"""
+    compute_internal_force_value(grad_i::Vec{3,F}, σ::SymmetricTensor{2,3,F}, α::Int) where {F}
+
+Scalar factor for the Galerkin internal-force row of a displacement test function
+associated with shape function ``N_i`` (gradient ``\\nabla N_i``) and Cartesian
+component ``\\alpha``:
+
+``\\sigma_{j\\alpha} \\, \\partial N_i / \\partial x_j``
+
+(sum over ``j = 1\\ldots 3``). The caller multiplies by ``\\det J \\cdot w`` per
+quadrature point and accumulates over IPs and elements.
+
+Cauchy stress ``\\sigma`` is the value stored in the material workspace at the IP
+(small-strain or finite-strain model, depending on the constitutive update).
+
+Zero allocation.
+"""
+@inline function compute_internal_force_value(
+    grad_i::Vec{3,F},
+    σ::SymmetricTensor{2,3,F},
+    α::Int,
+) where {F<:AbstractFloat}
+    s = zero(F)
+    @inbounds for j in 1:3
+        s += grad_i[j] * σ[j, α]
+    end
+    return s
+end
+
 """
     compute_stiffness_block(
         grad_k::Vec{D},
diff --git a/src/domains/continuum/material_element_lab.jl b/src/domains/continuum/material_element_lab.jl
index 940315e..ac00e68 100644
--- a/src/domains/continuum/material_element_lab.jl
+++ b/src/domains/continuum/material_element_lab.jl
@@ -101,7 +101,7 @@ function hex8_symmetric_uniaxial_eliminated_dirichlet(
 end
 
 """
-    material_lab_linear_elastic_uniaxial_solve(mesh, handler, elements, E, ν, δx; formulation = ContinuumFormulation{FullThreeD}())
+    material_lab_linear_elastic_uniaxial_solve(mesh, handler, elements, E, ν, δx; formulation = ContinuumFormulation{ThreeDimensional}())
 
 Assemble `K`, apply [`hex8_symmetric_uniaxial_eliminated_dirichlet`](@ref), solve `K u = 0`
 with elimination lift, and return `u`.
@@ -115,7 +115,7 @@ function material_lab_linear_elastic_uniaxial_solve(
     E::Float64,
     ν::Float64,
     δx::Float64;
-    formulation = ContinuumFormulation{FullThreeD}(),
+    formulation = ContinuumFormulation{ThreeDimensional}(),
 )
     material = LinearElastic(E = E, ν = ν)
     kernel = ContinuumKernel(formulation, material, Displacement{3}())
diff --git a/src/domains/continuum/mixed_up_kernel.jl b/src/domains/continuum/mixed_up_kernel.jl
index f65db3c..42b6636 100644
--- a/src/domains/continuum/mixed_up_kernel.jl
+++ b/src/domains/continuum/mixed_up_kernel.jl
@@ -47,7 +47,7 @@ This is the first mixed kernel: `evaluate_entry` dispatches on
 using Tensors
 
 using ..JuliaFEM: AbstractKernel, AbstractFormulation
-using ..JuliaFEM: ContinuumFormulation, FullThreeD, AbstractContinuumTheory
+using ..JuliaFEM: ContinuumFormulation, ThreeDimensional, AbstractContinuumTheory
 using ..JuliaFEM: AbstractMaterial, Displacement
 using ..JuliaFEM: AssemblyMaterialWorkspace, compute_stress
 import ..JuliaFEM: qpoint_buffer_eltype, update_qpoint_buffer!, evaluate_entry,
@@ -63,7 +63,7 @@ Mixed `u`–`p` kernel: 3D vertex displacement (field 1) + scalar cell pressure
 (field 2). See the file-level docstring for the weak form.
 
 # Fields
-- `formulation::ContinuumFormulation{Theory}` — geometric driver (`FullThreeD`, …)
+- `formulation::ContinuumFormulation{Theory}` — geometric driver (`ThreeDimensional`, …)
 - `material::Mat` — mechanical material (`LinearElastic`, …)
 - `inv_bulk::Float64` — `1/κ` for the `−κ⁻¹ ∫ p q dΩ` term (`0` = incompressible limit)
 
@@ -72,7 +72,7 @@ Mixed `u`–`p` kernel: 3D vertex displacement (field 1) + scalar cell pressure
 ```julia
 S = @DOFSet{u::DOF{Displacement{3}, Vertex}, p::DOF{Float64, Cell}}
 kernel = MixedUPKernel(
-    ContinuumFormulation{FullThreeD}(),
+    ContinuumFormulation{ThreeDimensional}(),
     LinearElastic(E = 210e9, ν = 0.3),
     inv_bulk = 1.0 / (210e9 / 3),  # order-of-magnitude compressible term
 )
diff --git a/src/domains/continuum/stokes_mixed_kernel.jl b/src/domains/continuum/stokes_mixed_kernel.jl
index 53b3985..07a7bfb 100644
--- a/src/domains/continuum/stokes_mixed_kernel.jl
+++ b/src/domains/continuum/stokes_mixed_kernel.jl
@@ -42,7 +42,7 @@ Newtonian Stokes mixed kernel: `u` (vertex, three components) + scalar
 `p` on `Cell`. See the file-level docstring for the weak form.
 
 # Fields
-- `formulation::ContinuumFormulation{Theory}` — geometric driver (`FullThreeD`, …)
+- `formulation::ContinuumFormulation{Theory}` — geometric driver (`ThreeDimensional`, …)
 - `μ::Float64` — dynamic viscosity (Stokes: `σ = 2μ ε(u)` with symmetric gradient `ε`)
 - `inv_bulk::Float64` — `1/κ` for the `−κ⁻¹ ∫ p q dΩ` term (`0` = incompressible)
 
@@ -50,7 +50,7 @@ Newtonian Stokes mixed kernel: `u` (vertex, three components) + scalar
 
 ```julia
 S = @DOFSet{u::DOF{Displacement{3}, Vertex}, p::DOF{Float64, Cell}}
-kernel = StokesMixedKernel(ContinuumFormulation{FullThreeD}(); μ = 1.0e-3, inv_bulk = 0.0)
+kernel = StokesMixedKernel(ContinuumFormulation{ThreeDimensional}(); μ = 1.0e-3, inv_bulk = 0.0)
 ```
 """
 struct StokesMixedKernel{Theory<:AbstractContinuumTheory} <: AbstractKernel
diff --git a/src/domains/continuum/types.jl b/src/domains/continuum/types.jl
index b59f214..1915cf0 100644
--- a/src/domains/continuum/types.jl
+++ b/src/domains/continuum/types.jl
@@ -1,12 +1,13 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+# SPDX-FileCopyrightText: 2015-2026 Jukka Aho
+# SPDX-License-Identifier: MIT
 
 """
 Concrete types for continuum mechanics formulations and theories.
 
-Abstract types are in abstract.jl, implementations are in formulations.jl.
+Abstract types are in `abstract.jl`; concrete theory structs and
+`ContinuumFormulation` live in this file.
 
-Must be included after abstract.jl.
+Must be included after `abstract.jl`.
 """
 
 # ============================================================================
@@ -14,13 +15,15 @@ Must be included after abstract.jl.
 # ============================================================================
 
 """
-    FullThreeD <: AbstractContinuumTheory
+    ThreeDimensional <: AbstractContinuumTheory
 
-Full 3D analysis with no simplifications. All six stress / flux components
-are carried; no geometric simplifications. Domain-agnostic — used by both
-`ContinuumKernel` (solid mechanics) and `HeatKernel` (heat conduction).
+Bulk three-dimensional model: no in-plane or axisymmetric reduction. All
+independent tensor components are retained at the quadrature point (six for
+symmetric mechanical stress / strain; three for isotropic flux, etc.).
+Domain-agnostic tag shared by `ContinuumKernel`, `HeatKernel`, Darcy-style
+kernels, and other `ContinuumFormulation{…}` drivers on 3D meshes.
 """
-struct FullThreeD <: AbstractContinuumTheory end
+struct ThreeDimensional <: AbstractContinuumTheory end
 
 """
     PlaneStress <: AbstractContinuumTheory
@@ -60,7 +63,7 @@ Used as a type tag inside `ContinuumKernel{Theory, Material, Field}`.
 # Examples
 
 ```julia
-ContinuumFormulation{FullThreeD}()
+ContinuumFormulation{ThreeDimensional}()
 ContinuumFormulation{PlaneStress}()
 ContinuumFormulation{PlaneStrain}()
 ContinuumFormulation{Axisymmetric}()
diff --git a/src/domains/continuum/update_geometry_cache.jl b/src/domains/continuum/update_geometry_cache.jl
index d6f0d25..4781e93 100644
--- a/src/domains/continuum/update_geometry_cache.jl
+++ b/src/domains/continuum/update_geometry_cache.jl
@@ -1,5 +1,5 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+# SPDX-FileCopyrightText: 2015-2026 Jukka Aho
+# SPDX-License-Identifier: MIT
 
 """
 Geometry cache update functions for continuum elements.
@@ -17,14 +17,24 @@ fields written are:
 
 - `geometry_cache.X` — node coordinates (one entry per element node)
 - `geometry_cache.N_data[ip, k]` — basis values
-- `geometry_cache.∇N_data[ip, k]` — physical gradients `∇N`
-- `geometry_cache.detJ_w[ip]` — `det(J) * w` for integration
+- `geometry_cache.∇N_data[ip, k]` — physical gradients `∇N` as `Vec{3}` (tangent to the
+  embedded surface for `dim(topology)==2`; all three components may be nonzero)
+- `geometry_cache.detJ_w[ip]` — `det(J) * w` for `D==3`, or `√(det G) · w` for `D==2`
 
 The function is allocation-free; it reads topology, basis, and the
 integration points from `element_cache` and writes back into the
-pre-allocated `geometry_cache` arrays. For each integration point the
-Jacobian `J = X ⊗ ∇_ξ N` is built on the fly, then physical gradients
-are obtained via `J^{-T} · ∇_ξ N`.
+pre-allocated `geometry_cache` arrays.
+
+For **`dim(topology) == 3`**, the Jacobian `J = Σ_k X_k ⊗ ∇_ξ N_k` is a `Tensor{2,3}`,
+inverted in the usual way, and `∇N` uses the full `Vec{3}` chain rule.
+
+For **`dim(topology) == 2`**, node coordinates are `Vec{3}` but the isoparametric map
+`x(ξ, η) ∈ ℝ³` is only two-parameter. Let `v_α = ∂x/∂ξ^α = Σ_k X_k ∂N_k/∂ξ^α` for
+`α ∈ {1,2}` (columns of the `3 × 2` Jacobian), `G_{αβ} = v_α · v_β` (Gram matrix),
+`detJ_w = √(det G) · w` (surface area measure on the embedded patch), and
+`∇N_k = v_1 (G^{-1} ∂_ξ N_k)_1 + v_2 (G^{-1} ∂_ξ N_k)_2`. This reduces to the former
+`(x, y)` / `J_2^{-T}` formula when the element lies in the global `xy` plane.
+Degenerate `det(G) ≤ 0` throws `ArgumentError`.
 """
 @inline function update_geometry_cache!(
     geometry_cache::GeometryCache,
@@ -35,8 +45,6 @@ are obtained via `J^{-T} · ∇_ξ N`.
     conn = mesh.connectivity[elem_id]
     nnodes = length(conn)
 
-    # Extract node coordinates (mesh.nodes already contains Vec{3}).
-    # Indexed loop avoids the iterator allocation that `enumerate` introduces.
     @inbounds for i in 1:nnodes
         node = conn[i]
         geometry_cache.X[i] = mesh.nodes[node]
@@ -44,28 +52,75 @@ are obtained via `J^{-T} · ∇_ξ N`.
 
     ips = element_cache.ips
     nips = length(ips)
+    D = dim(element_cache.topology)
 
-    @inbounds for ip_idx in 1:nips
-        ip = ips[ip_idx]
-        ξ = ip.coords
+    if D == 3
+        @inbounds for ip_idx in 1:nips
+            ip = ips[ip_idx]
+            ξ = ip.coords
 
-        N_vals = get_basis_functions(  element_cache.topology, element_cache.basis, ξ)
-        dN_dξ  = get_basis_derivatives(element_cache.topology, element_cache.basis, ξ)
+            N_vals = get_basis_functions(  element_cache.topology, element_cache.basis, ξ)
+            dN_dξ  = get_basis_derivatives(element_cache.topology, element_cache.basis, ξ)
 
-        # Jacobian J = X ⊗ ∇_ξ N
-        J = geometry_cache.X[1] ⊗ dN_dξ[1]
-        for i in 2:nnodes
-            J += geometry_cache.X[i] ⊗ dN_dξ[i]
+            J = geometry_cache.X[1] ⊗ dN_dξ[1]
+            for i in 2:nnodes
+                J += geometry_cache.X[i] ⊗ dN_dξ[i]
+            end
+
+            J_inv_T = transpose(inv(J))
+
+            for k in 1:nnodes
+                geometry_cache.N_data[ip_idx, k]  = N_vals[k]
+                geometry_cache.∇N_data[ip_idx, k] = J_inv_T ⋅ dN_dξ[k]
+            end
+
+            geometry_cache.detJ_w[ip_idx] = det(J) * ip.weight
         end
+    elseif D == 2
+        F = geometry_eltype(geometry_cache)
+        @inbounds for ip_idx in 1:nips
+            ip = ips[ip_idx]
+            ξ = ip.coords
 
-        J_inv_T = transpose(inv(J))
+            N_vals = get_basis_functions(  element_cache.topology, element_cache.basis, ξ)
+            dN_dξ  = get_basis_derivatives(element_cache.topology, element_cache.basis, ξ)
 
-        for k in 1:nnodes
-            geometry_cache.N_data[ip_idx, k]  = N_vals[k]
-            geometry_cache.∇N_data[ip_idx, k] = J_inv_T ⋅ dN_dξ[k]
+            v1 = zero(Vec{3,F})
+            v2 = zero(Vec{3,F})
+            for i in 1:nnodes
+                xk = geometry_cache.X[i]
+                di = dN_dξ[i]
+                v1 += xk * di[1]
+                v2 += xk * di[2]
+            end
+            G11 = v1 ⋅ v1
+            G12 = v1 ⋅ v2
+            G22 = v2 ⋅ v2
+            G = Tensor{2,2,F,4}((G11, G12, G12, G22))
+            detG = det(G)
+            if !(detG > zero(F))
+                throw(ArgumentError(
+                    "update_geometry_cache!: singular or non-right-handed 2D element " *
+                    "(det(G) = $detG); check node ordering and geometry",
+                ))
+            end
+            Ginv = inv(G)
+            detJ_w = sqrt(detG) * ip.weight
+
+            for k in 1:nnodes
+                geometry_cache.N_data[ip_idx, k] = N_vals[k]
+                h = Ginv ⋅ dN_dξ[k]
+                g = v1 * h[1] + v2 * h[2]
+                geometry_cache.∇N_data[ip_idx, k] = Vec{3,F}((g[1], g[2], g[3]))
+            end
+
+            geometry_cache.detJ_w[ip_idx] = detJ_w
         end
-
-        geometry_cache.detJ_w[ip_idx] = det(J) * ip.weight
+    else
+        throw(ArgumentError(
+            "update_geometry_cache!: topology spatial dimension $D is not supported " *
+            "(expected 2 or 3 for continuum assembly)",
+        ))
     end
 
     return nothing
diff --git a/src/domains/continuum/update_material_cache.jl b/src/domains/continuum/update_material_cache.jl
index 2d20a09..5d871da 100644
--- a/src/domains/continuum/update_material_cache.jl
+++ b/src/domains/continuum/update_material_cache.jl
@@ -1,5 +1,5 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+# SPDX-FileCopyrightText: 2015-2026 Jukka Aho
+# SPDX-License-Identifier: MIT
 
 """
 Material cache update functions for continuum elements.
@@ -10,6 +10,7 @@ Computes stress, tangent modulus, and internal state at integration points.
 using Tensors
 using ..JuliaFEM: GlobalMaterialCache, get_old_state, set_state!
 using ..JuliaFEM: continuum_kinematics, SmallStrainKinematics, GreenLagrangeKinematics
+using ..JuliaFEM: material_behavior, StatelessStrainDependent
 
 # ============================================================================
 # GLOBAL MATERIAL CACHE — behavior-dispatched implementations
@@ -126,6 +127,55 @@ end
     return nothing
 end
 
+"""
+    update_material_cache_stateless_strain!(
+        material_workspace, geometry_cache, material, element_cache, Δt,
+    ) -> Nothing
+
+Finite-strain / hyperelastic branch with **no** persistent integration-point
+state in [`GlobalMaterialCache`](@ref). Requires
+`material_behavior(material) isa StatelessStrainDependent`.
+
+`element_cache.u_buffer` must hold the current nodal displacements (vertex
+ordering matching `geometry_cache.∇N_data`). Used by the DOF-based Pass 1
+when the configuration vector is supplied (or zero displacement when it
+is not).
+
+Zero-allocation in the integration loop.
+"""
+@inline function update_material_cache_stateless_strain!(
+    material_workspace::AssemblyMaterialWorkspace,
+    geometry_cache::GeometryCache,
+    material::AbstractMaterial,
+    element_cache::ElementCache,
+    Δt::Float64,
+)
+    material_behavior(material) isa StatelessStrainDependent ||
+        throw(ArgumentError("update_material_cache_stateless_strain! requires StatelessStrainDependent material"))
+    nips = length(element_cache.ips)
+    nnodes = length(geometry_cache.X)
+    I = one(Tensor{2,3,Float64,9})
+
+    @inbounds for q in 1:nips
+        F = I
+        for k in 1:nnodes
+            u_k = element_cache.u_buffer[k]
+            ∇N_k_q = geometry_cache.∇N_data[q, k]
+            F += u_k ⊗ ∇N_k_q
+        end
+
+        C_tensor = symmetric(F' ⋅ F)
+        E = SymmetricTensor{2,3}(0.5 * (C_tensor - I))
+
+        σ, 𝔻, _ = compute_stress(material, E, NamedTuple(), Δt)
+
+        @inbounds material_workspace.fields[q] = (σ=σ, 𝔻=𝔻)
+        material_workspace.states[q] = NamedTuple()
+    end
+
+    return nothing
+end
+
 # StatelessStrainDependent — strain at each IP, no persistent state.
 @inline function update_material_cache!(
     material_workspace::AssemblyMaterialWorkspace,
@@ -137,30 +187,9 @@ end
     elem_id::Int,
     Δt::Float64,
 )
-    nips = length(element_cache.ips)
-    nnodes = length(geometry_cache.X)
-    I = one(Tensor{2,3,Float64,9})
-
-    @inbounds for q in 1:nips
-        # Deformation gradient F = I + ∇u
-        F = I
-        for k in 1:nnodes
-            u_k = element_cache.u_buffer[k]
-            ∇N_k_q = geometry_cache.∇N_data[q, k]
-            F += u_k ⊗ ∇N_k_q
-        end
-
-        # Green–Lagrange strain E = ½(F'F − I)
-        C_tensor = symmetric(F' ⋅ F)
-        E = SymmetricTensor{2,3}(0.5 * (C_tensor - I))
-
-        σ, 𝔻, _ = compute_stress(material, E, NamedTuple(), 0.0)
-
-        @inbounds material_workspace.fields[q] = (σ=σ, 𝔻=𝔻)
-        material_workspace.states[q] = NamedTuple()
-    end
-
-    return nothing
+    return update_material_cache_stateless_strain!(
+        material_workspace, geometry_cache, material, element_cache, Δt,
+    )
 end
 
 # StatefulStrainDependent — read old state from `global_cache`, compute the