feat(core): reorder includes and add FEDiscretization plus Gauss tuples

Partition metadata loads before DOF-based caches, new kernels and column
helpers wire into the load order, Gmsh leaves the default include graph, and
MPI docstrings describe kernel-less overloads plus internal-force hook.

- Add reference_gauss_tuples.jl and dof_based_pass1; include poroelastic and
  thermo_poroelastic kernels after thermo_elastic
- Hoist partitioning, halo_exchange, packed_layout, redundant_kernel_depwarn,
  ka_column_homogeneity, kernel_column ahead of dof_based_coo.jl
- Add physics/discretization.jl (FEDiscretization); drop gmsh_reader include
- Extend MPI matvec docstrings; declare mpi_partitioned_internal_force_owned!

src/JuliaFEM.jl

@@ -89,6 +89,7 @@ include("quadrature/gl_triangles.jl")
 include("quadrature/gl_wedges.jl")
 include("quadrature/gl_pyramids.jl")
 include("quadrature/gl_tensor_product.jl")
+include("quadrature/reference_gauss_tuples.jl")
 include("quadrature/gauss.jl")
 
 include("geometry/jacobian.jl")
@@ -165,6 +166,7 @@ include("domains/continuum/kernel.jl")
 include("domains/continuum/update_geometry_cache.jl")
 include("domains/continuum/update_element_cache.jl")
 include("domains/continuum/update_material_cache.jl")
+include("domains/continuum/dof_based_pass1.jl")
 include("domains/continuum/mixed_up_kernel.jl")
 include("domains/continuum/stokes_mixed_kernel.jl")
 include("domains/continuum/hellinger_reissner_kernel.jl")
@@ -176,6 +178,8 @@ include("domains/darcy/mixed_rt0.jl")
 include("domains/darcy/mixed_rt0_hex.jl")
 
 include("domains/thermo_elastic/kernel.jl")
+include("domains/poroelastic/kernel.jl")
+include("domains/thermo_poroelastic/kernel.jl")
 
 # ============================================================================
 # Level 7 - Concrete mesh + DOFHandler + DOF-based assembler + matrix-free
@@ -198,17 +202,23 @@ include("dofs/dof_connectivity.jl")
 include("interface/interface_mesh.jl")
 include("interface/interface_dof_handler.jl")
 
-# DOF-based assembler — CPU plus the backend-agnostic KernelAbstractions
+# Partition metadata + packed layouts must load before `dof_based_coo.jl`
-# port (CPU(), CUDABackend(), MetalBackend(), AMDGPUBackend(), oneAPIBackend();
+# (method signatures use `PartitionPackedLayout`). `packed_layout.jl` depends on
-# locally validated against CPU()).
+# `MeshPartitionLayout` / halo helpers from `partitioning.jl` / `halo_exchange.jl`.
-include("assemblers/dof_based/dof_based_coo.jl")
-include("assemblers/dof_based/dof_based_coo_ka.jl")
-
-# Partition metadata + multiply-buffer layouts (MPI/GPU hooks; default is
-# single-process identity).
 include("assemblers/partitioning.jl")
 include("assemblers/halo_exchange.jl")
 include("assemblers/packed_layout.jl")
+
+include("assemblers/redundant_kernel_depwarn.jl")
+
+# DOF-based assembler — CPU plus the backend-agnostic KernelAbstractions
+# port (CPU(), CUDABackend(), MetalBackend(), AMDGPUBackend(), oneAPIBackend();
+# locally validated against CPU()).
+include("assemblers/dof_based/ka_column_homogeneity.jl")
+include("assemblers/dof_based/kernel_column.jl")
+include("assemblers/dof_based/dof_based_coo.jl")
+include("assemblers/dof_based/dof_based_coo_ka.jl")
+
 include("assemblers/partitioned_matvec.jl")
 
 # Matrix-free path: declarative Dirichlet / MPC constraints, declarative
@@ -222,6 +232,8 @@ include("assemblers/matrix_free/preconditioners.jl")
 include("assemblers/matrix_free/eigensolve.jl")
 include("assemblers/matrix_free/loads.jl")
 
+include("physics/discretization.jl")
+
 # ============================================================================
 # Level 8 - Mesh utilities and I/O
 # ============================================================================
@@ -230,9 +242,8 @@ include("mesh/structured.jl")
 
 include("domains/continuum/material_element_lab.jl")
 
-# Self-contained Gmsh mesh reader (defines its own `JuliaFEM.GmshReader`
+# External mesh formats (Gmsh, Abaqus, Netgen, …) → `Mesh{…}` belong in
-# submodule with `read_gmsh_mesh` / `GmshMesh`).
+# optional package extensions or side packages; see `src/io/README.md`.
-include("io/gmsh_reader.jl")
 
 # ============================================================================
 # Optional older API surface (Problem / Assembly / Solver / Analysis,
@@ -284,16 +295,28 @@ on the local partition, summed with [`MPI.Allreduce`](@ref). Requires `using MPI
 function mpi_owned_dot_global end
 
 """
+    mpi_partitioned_operator_matvec!(
+        Ap, p, packed, work, recv_vals, send_vals,
+        layout, exchange, cache, assembler, mesh, comm;
+        dirichlet = nothing, mpi_requests = nothing,
+        configuration = nothing, global_material_cache = nothing, Δt = 0.0,
+    )
     mpi_partitioned_operator_matvec!(
         Ap, p, packed, work, recv_vals, send_vals,
         layout, exchange, cache, assembler, kernel, mesh, comm;
         dirichlet = nothing, mpi_requests = nothing,
+        configuration = nothing, global_material_cache = nothing, Δt = 0.0,
     )
 
 Replicated-global matrix-free product including optional [`PenaltyDirichlet`](@ref) post-hook.
 Loads `JuliaFEMMPIExt` when both `JuliaFEM` and `MPI` are imported.
 
-Keyword `mpi_requests` is forwarded to [`exchange_matvec_halos_mpi!`](@ref).
+The primary form reads the volume kernel from `cache.kernel_column`. The overload
+with `kernel` ignores it (backward compatibility).
+
+Keywords `mpi_requests` (to [`exchange_matvec_halos_mpi!`](@ref)) and
+`configuration` / `global_material_cache` / `Δt` (to [`apply_K_owned_rows!`](@ref), nonlinear continuum Pass~1)
+are forwarded when `JuliaFEMMPIExt` is loaded.
 """
 function mpi_partitioned_operator_matvec! end
 
@@ -307,10 +330,17 @@ entries only (aligned with [`copy_owned_subset_to_packed_owned_prefix!`](@ref)).
 function mpi_owned_dot_local end
 
 """
+    mpi_partitioned_operator_matvec_owned!(
+        Ap_owned, p_owned, packed, recv_vals, send_vals,
+        layout, exchange, cache, assembler, mesh, comm;
+        dirichlet = nothing, mpi_requests = nothing,
+        configuration = nothing, global_material_cache = nothing, Δt = 0.0,
+    )
     mpi_partitioned_operator_matvec_owned!(
         Ap_owned, p_owned, packed, recv_vals, send_vals,
         layout, exchange, cache, assembler, kernel, mesh, comm;
         dirichlet = nothing, mpi_requests = nothing,
+        configuration = nothing, global_material_cache = nothing, Δt = 0.0,
     )
 
 Lean MPI matvec: [`copy_owned_subset_to_packed_owned_prefix!`](@ref), halo exchange,
@@ -318,10 +348,36 @@ Lean MPI matvec: [`copy_owned_subset_to_packed_owned_prefix!`](@ref), halo excha
 No `ndofs_global`-length trial workspace — trial DOFs stay in `packed`. **No** full-vector
 [`MPI.Allreduce!`](@ref) on `Ap`.
 
-Keyword `mpi_requests` is forwarded to [`exchange_matvec_halos_mpi!`](@ref).
+The primary form uses `cache.kernel_column`; the overload with `kernel` ignores it.
+
+Keywords `mpi_requests` and `configuration` / `global_material_cache` / `Δt` are forwarded when
+`JuliaFEMMPIExt` is loaded.
 
 Supports [`PenaltyDirichlet`](@ref) only on this path.
 """
 function mpi_partitioned_operator_matvec_owned! end
 
+"""
+    mpi_partitioned_internal_force_owned!(
+        f_int_owned, u_owned, packed, work, recv_vals, send_vals,
+        layout, exchange, cache, assembler, mesh, comm;
+        mpi_requests = nothing,
+        configuration = nothing, global_material_cache = nothing, Δt = 0.0,
+    )
+    mpi_partitioned_internal_force_owned!(
+        f_int_owned, u_owned, packed, work, recv_vals, send_vals,
+        layout, exchange, cache, assembler, kernel, mesh, comm;
+        mpi_requests = nothing,
+        configuration = nothing, global_material_cache = nothing, Δt = 0.0,
+    )
+
+Lean MPI path for [`assemble_internal_force!`](@ref) on owned rows: pack owned
+`u`, halo exchange, then [`apply_f_int_owned_rows_from_packed!`](@ref). The
+`work` buffer must have length `cache.ndofs` (expand + Pass~1 scratch).
+
+The overload with `kernel` ignores it (backward compatibility). Requires
+`using MPI` after `JuliaFEM` so `JuliaFEMMPIExt` is loaded.
+"""
+function mpi_partitioned_internal_force_owned! end
+
 end # module JuliaFEM

src/physics/discretization.jl

@@ -0,0 +1,53 @@
+# SPDX-FileCopyrightText: 2015-2026 Jukka Aho
+# SPDX-License-Identifier: MIT
+
+"""
+    FEDiscretization(elements, handler, mesh, cache, assembler)
+
+Lightweight bundle of the usual DOF-based finite-element assembly inputs:
+`elements` from [`create_elements!`](@ref), the [`DOFHandler`](@ref), the
+[`AbstractMesh`](@ref), a pre-built [`DOFBasedCOOCache`](@ref) (whose
+[`kernel_column`](@ref) holds the volume kernel), and a
+[`DOFBasedCOOAssembler`](@ref) tag.
+
+Use [`assemble!`](@ref)`(fe)` then [`linear_system`](@ref)`(fe)` to obtain
+`(K, f)` without threading several separate variables through driver code.
+
+This is intentionally minimal: it does not own the mesh or elements, impose
+boundary conditions, or choose a solver. The kernel is read from
+`fe.cache.kernel_column` only; there is no separate redundant kernel field.
+"""
+struct FEDiscretization{E,H,M,C,A}
+    elements::Vector{E}
+    handler::H
+    mesh::M
+    cache::C
+    assembler::A
+end
+
+@inline function assemble!(fe::FEDiscretization; kwargs...)
+    return assemble!(fe.cache, fe.assembler, fe.mesh; kwargs...)
+end
+
+@inline function assemble_internal_force!(f::AbstractVector{Float64}, fe::FEDiscretization; kwargs...)
+    return assemble_internal_force!(f, fe.cache, fe.assembler, fe.mesh; kwargs...)
+end
+
+@inline function nonlinear_equilibrium_residual!(
+    r::AbstractVector{Float64},
+    f_ext::AbstractVector{Float64},
+    f_work::AbstractVector{Float64},
+    fe::FEDiscretization,
+    u::AbstractVector{Float64};
+    kwargs...,
+)
+    return nonlinear_equilibrium_residual!(
+        r, f_ext, f_work, fe.cache, fe.assembler, fe.mesh, u; kwargs...,
+    )
+end
+
+@inline function linear_system(fe::FEDiscretization)
+    return extract_system(fe.cache)
+end
+
+@inline total_dofs(fe::FEDiscretization) = fe.handler.total_dofs

src/quadrature/reference_gauss_tuples.jl

@@ -0,0 +1,34 @@
+# SPDX-FileCopyrightText: 2015-2026 Jukka Aho
+# SPDX-License-Identifier: MIT
+
+# Reference-domain Gauss-Legendre tuples for hand-coded integrals. Values match
+# `get_quadrature_points` for the same rule (see `gl_tensor_product.jl`, `gl_triangles.jl`),
+# including the 3-point triangle rule (`REF_GAUSS_TRIANGLE_ORDER2`).
+
+"Two-point Gauss-Legendre on reference segment `[-1, 1]`; each tuple is `(ξ, w)`; weights sum to `2`."
+const REF_GAUSS_SEGMENT_ORDER2 = ((-1.0 / √3, 1.0), (1.0 / √3, 1.0))
+
+"Tensor-product 2×2 Gauss-Legendre on reference square `[-1, 1]²`; each tuple is `(ξ, η, w)`."
+const REF_GAUSS_QUAD_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),
+)
+
+"One-point centroid rule on the unit-area-right-triangle natural coordinates; tuple `(ξ, η, w)` with `w = 1/2`."
+const REF_GAUSS_TRIANGLE_CENTROID = ((1.0 / 3.0, 1.0 / 3.0, 0.5),)
+
+"""
+Three-point Gauss-Legendre rule (degree 2) on the reference triangle `ξ, η ≥ 0`,
+`ξ + η ≤ 1` (area `1/2`). Each tuple is `(ξ, η, w)` with `w = 1/6`.
+
+Same weights as `get_quadrature_points(Triangle{3}, GaussLegendre{2}())` in
+`gl_triangles.jl` (default 3-point rule); point order matches
+`MixedDarcyTet4BoundaryNormalFluxLoad` in `matrix_free/loads.jl`.
+"""
+const REF_GAUSS_TRIANGLE_ORDER2 = (
+    (1.0 / 6.0, 1.0 / 6.0, 1.0 / 6.0),
+    (2.0 / 3.0, 1.0 / 6.0, 1.0 / 6.0),
+    (1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0),
+)