From 782e559d4bf66d96b6ebcd82a262d64c4c80ca29 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sun, 9 Nov 2025 17:09:11 +0200 Subject: [PATCH] refactor(basis): Non-parametric AbstractBasis for dynamic topology dimensions MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Change AbstractBasis{dim} to AbstractBasis (remove dimension type parameter) - Enable Lagrange{T,P} <: AbstractBasis inheritance (T=topology, P=polynomial degree) - Replace interface: length/size → nnodes/ndims - Remove allocating wrappers: eval_basis(), eval_dbasis() - Add nnodes() for both Lagrange instances and types - Implement nnodes formulas for all topologies: * Segment: P+1 * Triangle: (P+1)(P+2)/2 * Quadrilateral: (P+1)² * Tetrahedron: (P+1)(P+2)(P+3)/6 * Hexahedron: (P+1)³ * Pyramid: hardcoded (5, 13, 29) * Wedge: (P+1)²(P+2)/2 - Add nnodes() for old topology names (Tri3, Quad4, etc.) for backwards compatibility - BREAKING: All AbstractBasis{dim} code incompatible - Rationale: Lagrange dimension comes from topology at runtime, not compile-time constant --- src/basis/abstract.jl | 89 +++++++++++++++++++++++-------------------- 1 file changed, 48 insertions(+), 41 deletions(-) diff --git a/src/basis/abstract.jl b/src/basis/abstract.jl index af73f5a..5a2aebc 100644 --- a/src/basis/abstract.jl +++ b/src/basis/abstract.jl @@ -15,61 +15,36 @@ export Vec const Vecish{N,T} = Union{NTuple{N,T},Vec{N,T}} """ - AbstractBasis{dim} + AbstractBasis Abstract base type for all finite element basis functions. -# Type parameter -- `dim`: Dimensionality of the reference element (1, 2, or 3) - # Interface requirements Concrete basis types must implement: -- `Base.length(::Type{<:AbstractBasis})` - Number of basis functions -- `Base.size(::Type{<:AbstractBasis})` - (dim, n_basis) -- `get_reference_element_coordinates(::Type{<:AbstractBasis})` - Reference coordinates +- `nnodes(::Type{<:AbstractBasis})` - Number of basis functions +- `Base.ndims(::Type{<:AbstractBasis})` - Spatial dimension - `eval_basis!(::Type{<:AbstractBasis}, N, xi)` - Evaluate basis functions - `eval_dbasis!(::Type{<:AbstractBasis}, dN, xi)` - Evaluate basis derivatives # Example ```julia -struct Seg2 <: AbstractBasis{1} end -length(Seg2) == 2 -size(Seg2) == (1, 2) +struct Lagrange{Triangle,1} <: AbstractBasis end +nnodes(Lagrange{Triangle,1}) == 3 +ndims(Lagrange{Triangle,1}) == 2 ``` + +See also: [`Lagrange`](@ref), [`Serendipity`](@ref) """ -abstract type AbstractBasis{dim} end +abstract type AbstractBasis end # Forward methods on instances to types -# This allows calling methods on both Seg2 and Seg2() -Base.length(B::T) where {T<:AbstractBasis} = length(T) -Base.size(B::T) where {T<:AbstractBasis} = size(T) -# Updated signatures: eval_basis! and eval_dbasis! now return tuples -eval_basis!(B::T, ::Type{U}, xi) where {T<:AbstractBasis,U} = eval_basis!(T, U, xi) -eval_dbasis!(B::T, xi) where {T<:AbstractBasis} = eval_dbasis!(T, xi) +# This allows calling methods on both Lagrange{Triangle,1} and Lagrange{Triangle,1}() +Base.ndims(B::T) where {T<:AbstractBasis} = ndims(T) +nnodes(B::T) where {T<:AbstractBasis} = nnodes(T) -# Allocating versions (convenience wrappers) - now they just call and collect -""" - eval_basis(basis::AbstractBasis{dim}, xi) -> Vector{Float64} - -Evaluate basis functions at point `xi`, allocating return vector. - -See also: [`eval_basis!`](@ref) for non-allocating version that returns tuple. -""" -eval_basis(B::AbstractBasis{dim}, ::Type{T}, xi) where {dim,T} = collect(eval_basis!(B, T, xi)) - -""" - eval_dbasis(basis::AbstractBasis{dim}, xi) -> Vector{Vec{dim, Float64}} - -Evaluate basis function derivatives at point `xi`, allocating return vector. - -See also: [`eval_dbasis!`](@ref) for non-allocating version that returns tuple. -""" -eval_dbasis(B::AbstractBasis{dim}, xi) where {dim} = collect(eval_dbasis!(B, xi)) - -# Declare interface functions (will be implemented by basis generator) -function get_reference_element_coordinates end +# Declare interface functions (will be implemented by basis generator or specific basis types) function eval_basis! end function eval_dbasis! end @@ -134,9 +109,7 @@ nnodes(basis) # → 9 See also: [`AbstractBasis`](@ref), [`Serendipity`](@ref), [`Nedelec`](@ref) """ -# Lagrange inherits from AbstractBasis but dimension is determined by topology -# We cannot compute dim(T()) at type definition time, so we use methods instead -struct Lagrange{T<:AbstractTopology,P} end +struct Lagrange{T<:AbstractTopology,P} <: AbstractBasis end # Define interface methods for Lagrange # Dimension comes from topology @@ -148,32 +121,66 @@ Base.ndims(::Lagrange{T,P}) where {T,P} = dim(T()) """ nnodes(::Lagrange{T, P}) where {T, P} + nnodes(::Type{Lagrange{T, P}}) where {T, P} Compute number of nodes for Lagrange basis of degree P on topology T. +Works with both instances and types. """ # 1D: Segment nnodes(::Lagrange{Segment,P}) where {P} = P + 1 +nnodes(::Type{Lagrange{Segment,P}}) where {P} = P + 1 # 2D: Triangle (simplex) nnodes(::Lagrange{Triangle,P}) where {P} = div((P + 1) * (P + 2), 2) +nnodes(::Type{Lagrange{Triangle,P}}) where {P} = div((P + 1) * (P + 2), 2) # 2D: Quadrilateral (tensor product) nnodes(::Lagrange{Quadrilateral,P}) where {P} = (P + 1)^2 +nnodes(::Type{Lagrange{Quadrilateral,P}}) where {P} = (P + 1)^2 # 3D: Tetrahedron (simplex) nnodes(::Lagrange{Tetrahedron,P}) where {P} = div((P + 1) * (P + 2) * (P + 3), 6) +nnodes(::Type{Lagrange{Tetrahedron,P}}) where {P} = div((P + 1) * (P + 2) * (P + 3), 6) # 3D: Hexahedron (tensor product) nnodes(::Lagrange{Hexahedron,P}) where {P} = (P + 1)^3 +nnodes(::Type{Lagrange{Hexahedron,P}}) where {P} = (P + 1)^3 # 3D: Pyramid (mixed) # Pyramids don't follow a simple formula, so hardcode for known degrees nnodes(::Lagrange{Pyramid,1}) = 5 +nnodes(::Type{Lagrange{Pyramid,1}}) = 5 nnodes(::Lagrange{Pyramid,2}) = 13 +nnodes(::Type{Lagrange{Pyramid,2}}) = 13 nnodes(::Lagrange{Pyramid,3}) = 29 +nnodes(::Type{Lagrange{Pyramid,3}}) = 29 # 3D: Wedge/Prism (triangle × segment tensor product) nnodes(::Lagrange{Wedge,P}) where {P} = div((P + 1)^2 * (P + 2), 2) +nnodes(::Type{Lagrange{Wedge,P}}) where {P} = div((P + 1)^2 * (P + 2), 2) + +# Also need nnodes for the higher-order topology types themselves (Tet10, Tri6, etc.) +# These forward to the topology's nnodes() method +nnodes(::Type{T}) where {T<:AbstractTopology} = nnodes(T()) + +# For backwards compatibility, support old topology type names as if they were basis types +# This handles cases where code uses Tri3, Quad4, etc. as basis types +nnodes(::Type{Tri3}) = 3 +nnodes(::Type{Tri6}) = 6 +nnodes(::Type{Tri7}) = 7 +nnodes(::Type{Quad4}) = 4 +nnodes(::Type{Quad8}) = 8 +nnodes(::Type{Quad9}) = 9 +nnodes(::Type{Seg2}) = 2 +nnodes(::Type{Seg3}) = 3 +nnodes(::Type{Tet4}) = 4 +nnodes(::Type{Tet10}) = 10 +nnodes(::Type{Hex8}) = 8 +nnodes(::Type{Hex20}) = 20 +nnodes(::Type{Hex27}) = 27 +nnodes(::Type{Pyr5}) = 5 +nnodes(::Type{Wedge6}) = 6 +nnodes(::Type{Wedge15}) = 15 # Export the new parametric type and node count function export Lagrange, nnodes