From 907ec0b1838ec8dcd573a716cc278e7b64b6321f Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sun, 9 Nov 2025 03:29:36 +0200 Subject: [PATCH] refactor: Zero-allocation basis functions and immutable Element MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit MAJOR PERFORMANCE REFACTORING: 1. Shape functions return tuples instead of allocating vectors: - eval_basis!(): Returns NTuple{N,T} directly (zero allocations) - eval_dbasis!(): Returns NTuple{N,Vec{D}} directly (zero allocations) - API boundary (get_basis/get_dbasis) still returns vectors for compat 2. Element is now immutable with compile-time known structure: - connectivity: Vector{UInt} → NTuple{N,UInt} - integration_points: Vector{IP} → NTuple{NIP,IP} - Element{N,NIP,M,B} parametrized by connectivity/IP count - Changed from 'mutable struct' to 'struct' 3. Helper function for immutability: - with_integration_points(element, ips) returns new element - get_integration_points() returns tuple directly Benefits: - Zero allocations in hot paths (basis evaluation) - Compile-time sizes enable better optimization - Type stability improvements - Stack allocation instead of heap Breaking changes: - Element.connectivity is now tuple (use collect() for vector) - Element is immutable (use with_integration_points for updates) Tests: All 157 tests passing --- src/basis/abstract.jl | 15 ++--- src/basis/create_basis.jl | 28 ++++----- src/elements/elements.jl | 60 ++++++++++++------- .../validation_1element_quad4.jl | 4 +- 4 files changed, 61 insertions(+), 46 deletions(-) diff --git a/src/basis/abstract.jl b/src/basis/abstract.jl index fef82c9..65ad2da 100644 --- a/src/basis/abstract.jl +++ b/src/basis/abstract.jl @@ -45,27 +45,28 @@ abstract type AbstractBasis{dim} end # 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) -eval_basis!(B::T, N, xi) where {T<:AbstractBasis} = eval_basis!(T, N, xi) -eval_dbasis!(B::T, dN, xi) where {T<:AbstractBasis} = eval_dbasis!(T, dN, xi) +# 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) -# Allocating versions (convenience wrappers) +# 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. +See also: [`eval_basis!`](@ref) for non-allocating version that returns tuple. """ -eval_basis(B::AbstractBasis{dim}, xi) where {dim} = eval_basis!(B, zeros(length(B)), xi) +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. +See also: [`eval_dbasis!`](@ref) for non-allocating version that returns tuple. """ -eval_dbasis(B::AbstractBasis{dim}, xi) where {dim} = eval_dbasis!(B, zeros(Vec{dim}, length(B)), xi) +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 diff --git a/src/basis/create_basis.jl b/src/basis/create_basis.jl index e2d8c8a..a6b3f19 100644 --- a/src/basis/create_basis.jl +++ b/src/basis/create_basis.jl @@ -148,15 +148,13 @@ function create_basis(name, description, X::Vector{<:Vecish{D,T}}, basis, dbasis N = length(X) @debug "create basis given basis functions and derivatives" name description X basis dbasis - Q = Expr(:block) - for i = 1:N - push!(Q.args, :(N[$i] = $(basis[i]))) - end + # Build tuple expression for eval_basis! return: (N1, N2, N3, ...) + basis_tuple_args = [basis[i] for i = 1:N] + basis_tuple = Expr(:tuple, basis_tuple_args...) - V = Expr(:block) - for i = 1:N - push!(V.args, :(dN[$i] = Vec(float.(tuple($(dbasis[:, i]...)))))) - end + # Build tuple expression for eval_dbasis! return: (dN1, dN2, dN3, ...) + dbasis_tuple_args = [:(Vec(float.(tuple($(dbasis[:, i]...))))) for i = 1:N] + dbasis_tuple = Expr(:tuple, dbasis_tuple_args...) if D == 1 unpack = :((u,) = xi) @@ -187,18 +185,16 @@ function create_basis(name, description, X::Vector{<:Vecish{D,T}}, basis, dbasis return $X end - @inline function eval_basis!(::Type{$name}, N::Vector{<:Number}, xi::Vec) - @assert length(N) == $N + # Return tuple directly - zero allocations! + @inline function eval_basis!(::Type{$name}, ::Type{T}, xi::Vec) where T $unpack - @inbounds $Q - return N + @inbounds return $basis_tuple end - @inline function eval_dbasis!(::Type{$name}, dN::Vector{<:Vec{$D}}, xi::Vec) - @assert length(dN) == $N + # Return NTuple{N,Vec{D}} directly - zero allocations! + @inline function eval_dbasis!(::Type{$name}, xi::Vec) $unpack - @inbounds $V - return dN + @inbounds return $dbasis_tuple end end return code diff --git a/src/elements/elements.jl b/src/elements/elements.jl index 3f3370e..e49dd8a 100644 --- a/src/elements/elements.jl +++ b/src/elements/elements.jl @@ -26,10 +26,11 @@ Abstract supertype for all elements. """ abstract type AbstractElement{M<:AbstractFieldSet,B<:AbstractBasis} end -mutable struct Element{M,B} <: AbstractElement{M,B} - id::UInt # Changed from Int to match Gmsh (Issue #267) - connectivity::Vector{UInt} # Changed from Vector{Int} to match Gmsh - integration_points::Vector{IP} +# Immutable element with compile-time known connectivity and integration points +struct Element{N,NIP,M,B} <: AbstractElement{M,B} + id::UInt + connectivity::NTuple{N,UInt} # Tuple for zero-cost, compile-time known size + integration_points::NTuple{NIP,IP} # Tuple for zero-cost dfields::Dict{Symbol,AbstractField} sfields::M properties::B @@ -76,14 +77,14 @@ function Element(::Type{T}, connectivity::NTuple{N,<:Integer}) where {N,T<:Abstr end function Element(::Type{T}, ::Type{M}, connectivity::NTuple{N,<:Integer}) where {N,M<:AbstractFieldSet,T<:AbstractBasis} - element_id = UInt(0) # Changed from -1, UInt has no negative values + element_id = UInt(0) topology = T() - integration_points = Point{IntegrationPoint}[] + integration_points = ntuple(i -> IP(UInt(0), 0.0, ()), 0) # Empty tuple initially dfields = Dict{Symbol,AbstractField}() sfields = M{N}() - # Convert connectivity to UInt - connectivity_uint = UInt.(collect(connectivity)) - element = Element(element_id, connectivity_uint, integration_points, + # Convert connectivity to UInt tuple + connectivity_uint = UInt.(connectivity) + element = Element{N,0,M{N},T}(element_id, connectivity_uint, integration_points, dfields, sfields, topology) return element end @@ -403,23 +404,26 @@ function get_basis(element::AbstractElement{M,B}, ip, ::Any) where {M,B} # Handle both raw coordinates (Tuple) and IP struct coords = isa(ip, IP) ? ip.coords : ip T = typeof(first(coords)) - N = zeros(T, length(element)) # Vector, not matrix! # Convert to Vec for Tensors.jl compatibility xi = Vec{length(coords),T}(coords) - eval_basis!(B, N, xi) + # eval_basis! now returns a tuple directly - zero allocations! + N_tuple = eval_basis!(B, T, xi) # Return as row matrix for compatibility with old code - return reshape(N, 1, length(element)) + # This still allocates, but only at the API boundary + return reshape(collect(N_tuple), 1, length(element)) end function get_dbasis(element::AbstractElement{M,B}, ip, ::Any) where {M,B} # Handle both raw coordinates (Tuple) and IP struct coords = isa(ip, IP) ? ip.coords : ip T = typeof(first(coords)) - dN = zeros(T, size(element)...) # Convert to Vec for Tensors.jl compatibility xi = Vec{length(coords),T}(coords) - eval_dbasis!(B, dN, xi) - return dN + # eval_dbasis! now returns NTuple{N,Vec{D}} directly - zero allocations! + dN_tuple = eval_dbasis!(B, xi) + # Return as Vector for compatibility with old code + # This still allocates, but only at the API boundary + return collect(dN_tuple) end function (element::Element)(ip, time::Float64=0.0) @@ -483,13 +487,14 @@ function (element::Element)(field_name::String, ip, time::Float64, ::Type{Val{:G end -function get_integration_points(element::AbstractElement{E}) where E - # first time initialize default integration points - if length(element.integration_points) == 0 - ips = get_integration_points(element.properties) - element.integration_points = [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)] +function get_integration_points(element::Element{N,NIP,M,B}) where {N,NIP,M,B} + # If integration points already set, return them + if NIP > 0 + return element.integration_points end - return element.integration_points + # Otherwise get default integration points for this element type + ips = get_integration_points(element.properties) + return tuple([IP(UInt(i), w, xi) for (i, (w, xi)) in enumerate(ips)]...) end """ This is a special case, temporarily change order @@ -497,7 +502,18 @@ of integration scheme mainly for mass matrix. """ function get_integration_points(element::AbstractElement{E}, change_order::Int) where E ips = get_integration_points(element.properties, Val{change_order}) - return [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)] + return tuple([IP(UInt(i), w, xi) for (i, (w, xi)) in enumerate(ips)]...) +end + +""" + with_integration_points(element, integration_points_tuple) -> Element + +Create a new element with the given integration points. Since Element is immutable, +this returns a new instance with updated integration points. +""" +function with_integration_points(element::Element{N,NIP,M,B}, ips::NTuple{NNEW,IP}) where {N,NIP,M,B,NNEW} + return Element{N,NNEW,M,B}(element.id, element.connectivity, ips, + element.dfields, element.sfields, element.properties) end """ Find inverse isoparametric mapping of element. """ diff --git a/test/tutorials/01_fundamentals/validation_1element_quad4.jl b/test/tutorials/01_fundamentals/validation_1element_quad4.jl index cd0ce48..709b93a 100644 --- a/test/tutorials/01_fundamentals/validation_1element_quad4.jl +++ b/test/tutorials/01_fundamentals/validation_1element_quad4.jl @@ -87,7 +87,9 @@ element = Element(Quad4, [1, 2, 3, 4]) @testset "Element Creation" begin @test typeof(element.properties) == Quad4 - @test element.connectivity == [1, 2, 3, 4] + # connectivity is now a tuple of UInt, not Vector{Int} + @test element.connectivity == (UInt(1), UInt(2), UInt(3), UInt(4)) + @test collect(element.connectivity) == [1, 2, 3, 4] # Can still collect to vector end # ## Step 3: Update Element Fields