From 6fd99fa32335573915fabf6c20a79e074d3d2352 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sun, 9 Nov 2025 17:29:35 +0200 Subject: [PATCH] feat(basis): Update generator for parametric Lagrange{T,P} architecture MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Changed create_basis() signature from (name, desc, X, ...) to (topology_type, poly_degree, desc, X, ...) - Generator now produces methods for Lagrange{Segment,1}, Lagrange{Triangle,1}, etc. - Added ELEMENT_TO_LAGRANGE mapping dict (old names → topology_type + poly_degree) - Fixed reference coordinates to return tuples instead of vectors - Removed struct definitions (now use parametric Lagrange{T,P} type) - Removed Base.size(), Base.length() methods (use nnodes() instead) - Fixed typo: 'antsatz' → 'ansatz' All 15 element types regenerate successfully: Segment (1,2), Triangle (1,2), Quadrilateral (1,2,2), Tetrahedron (1,2), Hexahedron (1,2,2), Pyramid (1), Wedge (1,2) Tests pass for all element types. --- src/basis/lagrange_generator.jl | 106 +++++++++++++++++++++----------- 1 file changed, 71 insertions(+), 35 deletions(-) diff --git a/src/basis/lagrange_generator.jl b/src/basis/lagrange_generator.jl index 6510408..6720a4e 100644 --- a/src/basis/lagrange_generator.jl +++ b/src/basis/lagrange_generator.jl @@ -222,23 +222,23 @@ function calculate_interpolation_polynomial_derivatives(basis, D) return dbasis end -function create_basis(name, description, X::Vector{<:Vecish{D}}, p::Expr) where D - @debug "create basis given antsatz polynomial" name description X p +function create_basis(topology_type::Symbol, polynomial_degree::Int, description, X::Vector{<:Vecish{D}}, p::Expr) where D + @debug "create basis given ansatz polynomial" topology_type polynomial_degree description X p V = vandermonde_matrix(p, X) basis = calculate_interpolation_polynomials(p, V) - return create_basis(name, description, X, basis) + return create_basis(topology_type, polynomial_degree, description, X, basis) end -function create_basis(name, description, X::Vector{<:Vecish{D}}, basis::Vector) where D +function create_basis(topology_type::Symbol, polynomial_degree::Int, description, X::Vector{<:Vecish{D}}, basis::Vector) where D @assert length(X) == length(basis) - @debug "create basis given basis functions" name description X basis + @debug "create basis given basis functions" topology_type polynomial_degree description X basis dbasis = calculate_interpolation_polynomial_derivatives(basis, D) - return create_basis(name, description, Vec.(X), basis, dbasis) + return create_basis(topology_type, polynomial_degree, description, Vec.(X), basis, dbasis) end -function create_basis(name, description, X::Vector{<:Vecish{D,T}}, basis, dbasis) where {D,T} +function create_basis(topology_type::Symbol, polynomial_degree::Int, description, X::Vector{<:Vecish{D,T}}, basis, dbasis) where {D,T} N = length(X) - @debug "create basis given basis functions and derivatives" name description X basis dbasis + @debug "create basis given basis functions and derivatives" topology_type polynomial_degree description X basis dbasis # Build tuple expression for eval_basis! return: (N1, N2, N3, ...) basis_tuple_args = [basis[i] for i = 1:N] @@ -247,6 +247,10 @@ function create_basis(name, description, X::Vector{<:Vecish{D,T}}, basis, dbasis # 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...) + + # Build tuple expression for reference coordinates: (Vec(...), Vec(...), ...) + coord_tuple_args = [:(Vec{$D,$T}(tuple($(X[i]...)))) for i = 1:N] + coord_tuple = Expr(:tuple, coord_tuple_args...) if D == 1 unpack = :((u,) = xi) @@ -256,35 +260,35 @@ function create_basis(name, description, X::Vector{<:Vecish{D,T}}, basis, dbasis unpack = :((u, v, w) = xi) end + # Generate code for Lagrange{Topology, P} instead of TopologyBasis code = quote - struct $name <: AbstractBasis{$D} + # $description + function get_reference_element_coordinates(::Type{Lagrange{$topology_type,$polynomial_degree}}) + return $coord_tuple end - - Base.@pure function Base.size(::Type{$name}) - return ($D, $N) - end - - function Base.size(::Type{$name}, j::Int) - j == 1 && return $D - j == 2 && return $N - end - - Base.@pure function Base.length(::Type{$name}) - return $N - end - - function get_reference_element_coordinates(::Type{$name}) - return $X + + function get_reference_element_coordinates(::Lagrange{$topology_type,$polynomial_degree}) + return $coord_tuple end # Return tuple directly - zero allocations! - @inline function eval_basis!(::Type{$name}, ::Type{T}, xi::Vec) where T + @inline function eval_basis!(::Type{Lagrange{$topology_type,$polynomial_degree}}, ::Type{T}, xi::Vec) where T + $unpack + @inbounds return $basis_tuple + end + + @inline function eval_basis!(::Lagrange{$topology_type,$polynomial_degree}, ::Type{T}, xi::Vec) where T $unpack @inbounds return $basis_tuple end # Return NTuple{N,Vec{D}} directly - zero allocations! - @inline function eval_dbasis!(::Type{$name}, xi::Vec) + @inline function eval_dbasis!(::Type{Lagrange{$topology_type,$polynomial_degree}}, xi::Vec) + $unpack + @inbounds return $dbasis_tuple + end + + @inline function eval_dbasis!(::Lagrange{$topology_type,$polynomial_degree}, xi::Vec) $unpack @inbounds return $dbasis_tuple end @@ -294,6 +298,26 @@ end create_basis_and_eval(args...) = eval(create_basis(args...)) +# Mapping from old element names to (topology_type, polynomial_degree) +# This allows the generator to use the new parametric Lagrange{T,P} system +const ELEMENT_TO_LAGRANGE = Dict{String, Tuple{Symbol, Int}}( + "Seg2" => (:Segment, 1), + "Seg3" => (:Segment, 2), + "Tri3" => (:Triangle, 1), + "Tri6" => (:Triangle, 2), + "Quad4" => (:Quadrilateral, 1), + "Quad8" => (:Quadrilateral, 2), # Serendipity, special case + "Quad9" => (:Quadrilateral, 2), + "Tet4" => (:Tetrahedron, 1), + "Tet10" => (:Tetrahedron, 2), + "Hex8" => (:Hexahedron, 1), + "Hex20" => (:Hexahedron, 2), # Serendipity, special case + "Hex27" => (:Hexahedron, 2), + "Pyr5" => (:Pyramid, 1), + "Wedge6" => (:Wedge, 1), + "Wedge15" => (:Wedge, 2), # Serendipity, special case +) + # ============================================================================== # GENERATION SCRIPT - Run as: julia --project=. src/basis/lagrange_generator.jl # ============================================================================== @@ -755,16 +779,27 @@ if abspath(PROGRAM_FILE) == @__FILE__ println(output, "# ============================================================================") println(output) - # Export all basis types - basis_names = [Symbol(elem.name * "Basis") for elem in elements] - println(output, "# Export all basis types") - println(output, "export ", join(string.(basis_names), ", ")) + # No exports needed - we generate methods for the parametric Lagrange{T,P} type + println(output, "# This file generates methods for Lagrange{T,P} where:") + println(output, "# T = topology type (Segment, Triangle, Quadrilateral, etc.)") + println(output, "# P = polynomial degree (1, 2, 3, ...)") + println(output, "#") + println(output, "# Example: Lagrange{Triangle, 2} is a quadratic triangular element") println(output) for (i, elem) in enumerate(elements) println("[$i/$(length(elements))] Generating $(elem.name)...") try + # Map old element name to (topology_type, polynomial_degree) + if !haskey(ELEMENT_TO_LAGRANGE, elem.name) + println(" ⚠ Warning: $(elem.name) not in ELEMENT_TO_LAGRANGE mapping") + println(" Skipping...") + continue + end + + topology_type, poly_degree = ELEMENT_TO_LAGRANGE[elem.name] + # Convert coordinates to proper format (tuples, not vectors) coords = [tuple(coord...) for coord in elem.coordinates] @@ -776,10 +811,10 @@ if abspath(PROGRAM_FILE) == @__FILE__ end # Generate basis code using symbolic engine - # APPEND "Basis" SUFFIX to resolve name conflicts with topology types - basis_type_name = Symbol(elem.name * "Basis") + # Now generates methods for Lagrange{topology_type, poly_degree} basis_code_expr = create_basis( - basis_type_name, + topology_type, + poly_degree, elem.description, coords, ansatz_expr @@ -798,7 +833,8 @@ if abspath(PROGRAM_FILE) == @__FILE__ # Write to output with nice formatting println(output, "# " * "─"^78) - println(output, "# $(elem.name): $(elem.description)") + println(output, "# Lagrange{$(topology_type), $(poly_degree)}: $(elem.description)") + println(output, "# (Old name: $(elem.name))") println(output, "# " * "─"^78) println(output) println(output, basis_code_str)