From 212ffb0d66a9eba3444582c60d8d797e36dd9624 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Mon, 11 May 2026 02:40:28 +0300 Subject: [PATCH] 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 | 82 ++++++++++++++++++++---- src/physics/discretization.jl | 53 +++++++++++++++ src/quadrature/reference_gauss_tuples.jl | 34 ++++++++++ 3 files changed, 156 insertions(+), 13 deletions(-) create mode 100644 src/physics/discretization.jl create mode 100644 src/quadrature/reference_gauss_tuples.jl diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 6a82a81..2b4ddc5 100644 --- a/src/JuliaFEM.jl +++ b/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 -# port (CPU(), CUDABackend(), MetalBackend(), AMDGPUBackend(), oneAPIBackend(); -# locally validated against CPU()). -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). +# Partition metadata + packed layouts must load before `dof_based_coo.jl` +# (method signatures use `PartitionPackedLayout`). `packed_layout.jl` depends on +# `MeshPartitionLayout` / halo helpers from `partitioning.jl` / `halo_exchange.jl`. 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` -# submodule with `read_gmsh_mesh` / `GmshMesh`). -include("io/gmsh_reader.jl") +# External mesh formats (Gmsh, Abaqus, Netgen, …) → `Mesh{…}` belong in +# optional package extensions or side packages; see `src/io/README.md`. # ============================================================================ # 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 diff --git a/src/physics/discretization.jl b/src/physics/discretization.jl new file mode 100644 index 0000000..8953a45 --- /dev/null +++ b/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 diff --git a/src/quadrature/reference_gauss_tuples.jl b/src/quadrature/reference_gauss_tuples.jl new file mode 100644 index 0000000..4e8a57e --- /dev/null +++ b/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), +)