mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-02 01:06:23 +00:00
refactor(basis): Non-parametric AbstractBasis for dynamic topology dimensions
- 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
This commit is contained in:
Jukka Aho
+48
-41
@@ -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
|
||||
|
||||
Reference in New Issue
Block a user