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