diff --git a/src/basis/abstract.jl b/src/basis/abstract.jl
deleted file mode 100644
index e356e95..0000000
--- a/src/basis/abstract.jl
+++ /dev/null
@@ -1,374 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE
-
-# AbstractBasis type and interface
-# Consolidated from jl package
-
-using Tensors
-using LinearAlgebra
-# import Calculus  # Only needed for symbolic basis generation (create_basis.jl)
-
-# Re-export Vec for convenience (from Tensors.jl)
-export Vec
-
-# Type alias for coordinate inputs (tuples or Vec)
-const Vecish{N,T} = Union{NTuple{N,T},Vec{N,T}}
-
-"""
-    AbstractBasis
-
-Abstract base type for all finite element basis functions.
-
-# Interface requirements
-
-Concrete basis types must implement:
-- `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 Lagrange{Triangle,1} <: AbstractBasis end
-nnodes(Lagrange{Triangle,1}) == 3
-ndims(Lagrange{Triangle,1}) == 2
-```
-
-See also: [`Lagrange`](@ref), [`Serendipity`](@ref)
-"""
-abstract type AbstractBasis end
-
-# Forward methods on instances to types
-# 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)
-
-# ============================================================================
-# DEPRECATED API (Remove after migration)
-# ============================================================================
-# These functions are DEPRECATED and should not be used in new code.
-# They remain for backward compatibility during migration.
-#
-# OLD API: eval_basis!(basis_type, T, xi) → Returns NTuple{N, T}
-# NEW API: get_basis_functions(topology, basis, xi) → Returns NTuple{N, Float64}
-#
-# Migration guide:
-# OLD: N = eval_basis!(Lagrange{Triangle,1}, Float64, xi)
-# NEW: N = get_basis_functions(Triangle(), Lagrange{1}(), xi)
-#
-# OLD: dN = eval_dbasis!(Lagrange{Triangle,1}, xi)
-# NEW: dN = get_basis_derivatives(Triangle(), Lagrange{1}(), xi)
-#
-# Why deprecated:
-# - Topology and basis are separate concerns (should be passed separately)
-# - Naming is unclear (eval_basis! suggests mutation, but returns value)
-# - Type parameter T is redundant (always Float64 in practice)
-
-"""
-    eval_basis!(B::Type{<:AbstractBasis}, T::Type, xi::Vec) -> NTuple{N, T}
-
-**DEPRECATED:** Use `get_basis_functions(topology, basis, xi)` instead.
-
-This function will be removed in a future release. The new API separates
-topology and basis, has clearer naming, and follows consistent conventions.
-
-# Migration Example
-
-```julia
-# OLD (deprecated):
-N = eval_basis!(Lagrange{Triangle,1}, Float64, xi)
-
-# NEW (recommended):
-N = get_basis_functions(Triangle(), Lagrange{1}(), xi)
-```
-
-See: [`get_basis_functions`](@ref)
-"""
-# Note: eval_basis! stub is defined in basis_api.jl
-
-"""
-    eval_dbasis!(B::Type{<:AbstractBasis}, xi::Vec) -> NTuple{N, Vec{D}}
-
-**DEPRECATED:** Use `get_basis_derivatives(topology, basis, xi)` instead.
-
-This function will be removed in a future release. The new API separates
-topology and basis, has clearer naming, and follows consistent conventions.
-
-# Migration Example
-
-```julia
-# OLD (deprecated):
-dN = eval_dbasis!(Lagrange{Triangle,1}, xi)
-
-# NEW (recommended):
-dN = get_basis_derivatives(Triangle(), Lagrange{1}(), xi)
-```
-
-See: [`get_basis_derivatives`](@ref)
-"""
-# Note: eval_dbasis! stub is defined in basis_api.jl
-
-# ============================================================================
-# NEW API (Recommended)
-# ============================================================================
-# These are the recommended functions for all new code.
-# See basis_api.jl for full documentation and examples.
-
-"""
-    get_basis_functions(topology, basis, xi) -> NTuple{N, Float64}
-
-Evaluate all basis functions at parametric point `xi`.
-
-This is the **recommended API** for evaluating basis functions. See `basis_api.jl`
-for full documentation and performance benchmarks.
-
-# Quick Example
-
-```julia
-topology = Triangle()
-basis = Lagrange{1}()  # Linear
-xi = Vec(0.25, 0.25)
-
-N = get_basis_functions(topology, basis, xi)
-# Returns: (N1, N2, N3) as NTuple{3, Float64}
-```
-
-See also: [`get_basis_derivatives`](@ref)
-"""
-# Note: get_basis_functions stub is defined in basis_api.jl
-
-"""
-    get_basis_derivatives(topology, basis, xi) -> NTuple{N, Vec{D, Float64}}
-
-Evaluate all basis function derivatives at parametric point `xi`.
-
-This is the **recommended API** for evaluating basis derivatives. See `basis_api.jl`
-for full documentation and performance benchmarks.
-
-# Quick Example
-
-```julia
-topology = Triangle()
-basis = Lagrange{1}()  # Linear
-xi = Vec(0.25, 0.25)
-
-dN = get_basis_derivatives(topology, basis, xi)
-# Returns: (∇N1, ∇N2, ∇N3) as NTuple{3, Vec{2, Float64}}
-```
-
-See also: [`get_basis_functions`](@ref)
-"""
-# Note: get_basis_derivatives stub is defined in basis_api.jl
-
-# ============================================================================
-# Parametric Lagrange Basis Type (OLD - with topology parameter)
-# ============================================================================
-# NOTE: This is the OLD API with topology in the type parameter.
-# New code should use LagrangeP{P} below (topology passed separately).
-# This is kept for backward compatibility during migration.
-
-"""
-    Lagrange{T<:AbstractTopology, P} <: AbstractBasis
-
-Parametric Lagrange basis functions for topology `T` with polynomial degree `P`.
-
-The node count is automatically derived from the topology and polynomial degree:
-- `Lagrange{Triangle, 1}`: P1 → 3 nodes (vertices)
-- `Lagrange{Triangle, 2}`: P2 → 6 nodes (vertices + edge midpoints)
-- `Lagrange{Quadrilateral, 1}`: Q1 → 4 nodes (corners)
-- `Lagrange{Quadrilateral, 2}`: Q2 → 9 nodes (full tensor product)
-- `Lagrange{Tetrahedron, 1}`: P1 → 4 nodes (vertices)
-- `Lagrange{Tetrahedron, 2}`: P2 → 10 nodes (vertices + edge midpoints)
-- `Lagrange{Hexahedron, 1}`: Q1 → 8 nodes (corners)
-- `Lagrange{Hexahedron, 2}`: Q2 → 27 nodes (full tensor product)
-
-# Type Parameters
-- `T`: Topology type (Triangle, Quadrilateral, Tetrahedron, Hexahedron, etc.)
-- `P`: Polynomial degree (1 = linear, 2 = quadratic, 3 = cubic, ...)
-
-# Mathematical Background
-
-Lagrange basis functions satisfy the cardinal property:
-```
-Nᵢ(xⱼ) = δᵢⱼ  (Kronecker delta)
-```
-
-where `xⱼ` are the interpolation nodes.
-
-For polynomial degree P:
-- 1D: P+1 nodes
-- Triangle: (P+1)(P+2)/2 nodes
-- Quadrilateral: (P+1)² nodes (tensor product)
-- Tetrahedron: (P+1)(P+2)(P+3)/6 nodes
-- Hexahedron: (P+1)³ nodes (tensor product)
-
-# Example
-
-```julia
-# Linear triangle (P1)
-basis = Lagrange{Triangle, 1}()
-nnodes(basis)  # → 3
-
-# Quadratic triangle (P2)
-basis = Lagrange{Triangle, 2}()
-nnodes(basis)  # → 6
-
-# Bilinear quadrilateral (Q1)
-basis = Lagrange{Quadrilateral, 1}()
-nnodes(basis)  # → 4
-
-# Biquadratic quadrilateral (Q2)
-basis = Lagrange{Quadrilateral, 2}()
-nnodes(basis)  # → 9
-```
-
-See also: [`AbstractBasis`](@ref), [`Serendipity`](@ref), [`Nedelec`](@ref)
-"""
-struct Lagrange{T<:AbstractTopology,P} <: AbstractBasis end
-
-"""
-    Serendipity{T<:AbstractTopology, P} <: AbstractBasis
-
-Serendipity basis family for quadrilateral and hexahedral elements.
-
-Serendipity elements use a reduced set of nodes compared to full tensor-product Lagrange
-elements by omitting interior nodes while maintaining the polynomial order on element edges.
-
-# Common examples:
-- `Serendipity{Quadrilateral, 2}`: 8-node quadrilateral (no center node)
-- `Serendipity{Hexahedron, 2}`: 20-node hexahedron (no interior nodes)
-
-# Comparison with Lagrange:
-- Quad8 (Serendipity): 8 nodes (4 corners + 4 edge midpoints, no center)
-- Quad9 (Lagrange): 9 nodes (4 corners + 4 edge midpoints + center)
-
-See also: [`AbstractBasis`](@ref), [`Lagrange`](@ref)
-"""
-struct Serendipity{T<:AbstractTopology,P} <: AbstractBasis end
-
-# Define interface methods for Lagrange
-# Dimension comes from topology
-Base.ndims(::Type{Lagrange{T,P}}) where {T,P} = dim(T())
-Base.ndims(::Lagrange{T,P}) where {T,P} = dim(T())
-
-# Define interface methods for Serendipity
-Base.ndims(::Type{Serendipity{T,P}}) where {T,P} = dim(T())
-Base.ndims(::Serendipity{T,P}) where {T,P} = dim(T())
-
-# Node count formulas for different topologies and polynomial degrees
-# These replace the hardcoded node counts in old Tri3, Quad4, etc. types
-
-"""
-    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)
-
-"""
-    nnodes(::Serendipity{T, P}) where {T, P}
-    nnodes(::Type{Serendipity{T, P}}) where {T, P}
-
-Compute number of nodes for Serendipity basis of degree P on topology T.
-
-Serendipity elements omit interior nodes, using only edge and corner nodes:
-- Quadrilateral P=2: 8 nodes (4 corners + 4 edge midpoints)
-- Hexahedron P=2: 20 nodes (8 corners + 12 edge midpoints)
-"""
-# 2D: Quadrilateral (Serendipity)
-# P=2: 8 nodes (4 corners + 4 edge midpoints, no center)
-nnodes(::Serendipity{Quadrilateral,2}) = 8
-nnodes(::Type{Serendipity{Quadrilateral,2}}) = 8
-
-# 3D: Hexahedron (Serendipity)
-# P=2: 20 nodes (8 corners + 12 edge midpoints, no face/interior nodes)
-nnodes(::Serendipity{Hexahedron,2}) = 20
-nnodes(::Type{Serendipity{Hexahedron,2}}) = 20
-
-# 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
-
-# ============================================================================
-# Additional Interface Methods
-# ============================================================================
-
-"""
-    ndofs(basis::AbstractBasis) -> Int
-    ndofs(::Type{<:AbstractBasis}) -> Int
-
-Number of degrees of freedom for this basis.
-
-For standard nodal elements (Lagrange), ndofs == nnodes.
-For plate elements and others with multiple DOFs per node, ndofs > nnodes.
-
-# Example
-
-```julia
-ndofs(Lagrange{Triangle, 1}())  # → 3 (same as nnodes)
-ndofs(DKT())                     # → 9 (3 nodes × 3 DOFs/node)
-```
-
-See also: [`nnodes`](@ref), [`AbstractPlateBasis`](@ref)
-"""
-function ndofs end
-
-# Default: For standard elements, ndofs == nnodes
-ndofs(B::AbstractBasis) = nnodes(B)
-ndofs(B::Type{<:AbstractBasis}) = nnodes(B)
-
-# Export the new parametric type and interface functions
-export Lagrange, nnodes, ndofs
diff --git a/src/basis/api.jl b/src/basis/api.jl
new file mode 100644
index 0000000..5610a39
--- /dev/null
+++ b/src/basis/api.jl
@@ -0,0 +1,164 @@
+# Basis API - Definitions and Interfaces
+# This file combines the abstract basis type and the basis evaluation API.
+
+using Tensors
+using LinearAlgebra
+
+# Re-export Vec for convenience (from Tensors.jl)
+export Vec
+
+"""
+    AbstractBasisDescription
+
+Abstract description of how to construct or generate a basis for a given topology.
+Different description types can encode symbolic (Vandermonde) generation, external
+rules, or manually provided functions.
+"""
+abstract type AbstractBasisDescription end
+
+"""
+    AbstractBasis
+
+Abstract base type for all finite element basis families.
+
+Basis types describe interpolation schemes (H¹ nodal, H(curl), H(div), plate/shell, etc.).
+Topology is passed separately to evaluation routines—basis never owns geometry.
+"""
+abstract type AbstractBasis end
+
+"""
+    Lagrange{P} <: AbstractBasis
+
+Standard nodal Lagrange basis of polynomial order `P`.
+Topology is supplied separately to evaluation routines.
+"""
+struct Lagrange{P} <: AbstractBasis end
+
+"""
+    Serendipity{P} <: AbstractBasis
+
+Serendipity basis family (reduced tensor product) of order `P` for quads/hexes.
+Topology is supplied separately to evaluation routines.
+"""
+struct Serendipity{P} <: AbstractBasis end
+
+"""
+    VandermondeBasisDescription{Topo,Basis} <: AbstractBasisDescription
+
+Description for bases generated from a polynomial ansatz via a Vandermonde system.
+Stores the topology, basis family/order, and the polynomial terms used to build
+shape functions.
+"""
+struct VandermondeBasisDescription{Topo<:AbstractTopology,Basis<:AbstractBasis} <: AbstractBasisDescription
+    name::String
+    description::String
+    topology::Type{Topo}
+    family::Type{Basis}
+    ansatz::Tuple
+    function VandermondeBasisDescription(; name, description, topology::Type{Topo}, family::Type{Basis}, ansatz::Tuple) where {Topo<:AbstractTopology,Basis<:AbstractBasis}
+        new{Topo,Basis}(name, description, topology, family, ansatz)
+    end
+end
+
+basis_family(desc::VandermondeBasisDescription) = desc.family
+basis_topology(desc::VandermondeBasisDescription) = desc.topology
+basis_order(::VandermondeBasisDescription{<:AbstractTopology,<:Lagrange{P}}) where {P} = P
+basis_order(::VandermondeBasisDescription{<:AbstractTopology,<:Serendipity{P}}) where {P} = P
+basis_order(desc::VandermondeBasisDescription) = error("basis order not defined for $(desc.family)")
+
+reference_coordinates(desc::VandermondeBasisDescription) = reference_coordinates(desc.topology())
+
+# ---------------------------------------------------------------------------
+# Basis evaluation API (new design: topology passed separately)
+# ---------------------------------------------------------------------------
+
+"""
+    get_basis_functions(topology::AbstractTopology, basis::AbstractBasis, xi::Vec)
+
+Evaluate all basis functions at parametric point `xi`.
+
+Returns `NTuple{N, Float64}` where `N` is the number of basis functions for the
+given topology–basis combination. Implementations live in generated code
+(e.g., lagrange_generated.jl) or custom basis modules.
+"""
+function get_basis_functions end
+
+"""
+    get_basis_derivatives(topology::AbstractTopology, basis::AbstractBasis, xi::Vec)
+
+Evaluate all basis function derivatives at parametric point `xi`.
+
+Returns `NTuple{N, Vec{D, Float64}}` where `D` is the parametric dimension.
+Implementations live in generated code (e.g., lagrange_generated.jl) or
+custom basis modules.
+"""
+function get_basis_derivatives end
+
+"""
+    get_basis_function(topology, basis, xi, i::Int)
+
+Convenience accessor for a single basis function value.
+Equivalent to `get_basis_functions(topology, basis, xi)[i]`.
+"""
+@inline function get_basis_function(topology::AbstractTopology,
+    basis::AbstractBasis,
+    xi::Vec,
+    i::Int)
+    return get_basis_functions(topology, basis, xi)[i]
+end
+
+"""
+    get_basis_derivative(topology, basis, xi, i::Int)
+
+Convenience accessor for a single basis function derivative.
+Equivalent to `get_basis_derivatives(topology, basis, xi)[i]`.
+"""
+@inline function get_basis_derivative(topology::AbstractTopology,
+    basis::AbstractBasis,
+    xi::Vec,
+    i::Int)
+    return get_basis_derivatives(topology, basis, xi)[i]
+end
+
+# ---------------------------------------------------------------------------
+# Deprecated bridge (old API names)
+# ---------------------------------------------------------------------------
+
+"""
+    eval_basis!(basis_type, xi) (DEPRECATED)
+
+Use `get_basis_functions(topology, basis, xi)` instead.
+Provided temporarily for migration.
+"""
+function eval_basis! end
+
+"""
+    eval_dbasis!(basis_type, xi) (DEPRECATED)
+
+Use `get_basis_derivatives(topology, basis, xi)` instead.
+Provided temporarily for migration.
+"""
+function eval_dbasis! end
+
+# ---------------------------------------------------------------------------
+# Degrees of freedom utility
+# ---------------------------------------------------------------------------
+
+"""
+    ndofs(basis::AbstractBasis)
+    ndofs(::Type{<:AbstractBasis})
+
+Total degrees of freedom for this basis. For standard nodal bases, this is
+usually equal to the number of basis functions; specialized bases (e.g., plates)
+can override to return multiple DOFs per node.
+"""
+function ndofs end
+
+# Default: basis implementations should override; no assumption about nnodes here.
+ndofs(::AbstractBasis) = error("ndofs not implemented for this basis")
+ndofs(::Type{<:AbstractBasis}) = error("ndofs not implemented for this basis type")
+
+export AbstractBasis, Lagrange, Serendipity
+export ndofs
+export get_basis_functions, get_basis_derivatives, get_basis_function, get_basis_derivative
+export eval_basis!, eval_dbasis!
diff --git a/src/basis/basis_api.jl b/src/basis/basis_api.jl
deleted file mode 100644
index 7811190..0000000
--- a/src/basis/basis_api.jl
+++ /dev/null
@@ -1,210 +0,0 @@
-# Basis Function API - New Design (Nov 2025)
-#
-# This file implements the new basis function API based on comprehensive
-# benchmarking results (see docs/book/adr-003-basis-function-api.md).
-#
-# Key design decisions:
-# 1. Topology passed separately: get_basis_functions(topology, basis, xi)
-# 2. Return tuples (zero allocation, type-stable)
-# 3. Use simple runtime indexing for single access (fastest!)
-# 4. Clear naming: get_basis_functions, get_basis_derivatives
-#
-# Performance: 6.5 ns for Tet10 derivatives, zero allocations
-
-"""
-    get_basis_functions(topology::AbstractTopology, basis::AbstractBasis, xi::Vec)
-
-Evaluate all basis functions at parametric point `xi`.
-
-Returns `NTuple{N, Float64}` where N is the number of nodes/DOFs for the given
-topology and basis combination.
-
-# Arguments
-- `topology`: Element topology (e.g., `Tetrahedron()`, `Triangle()`)
-- `basis`: Interpolation scheme (e.g., `Lagrange{1}()`, `Lagrange{2}()`)
-- `xi`: Parametric coordinates as `Vec{D, T}` where D is spatial dimension
-
-# Returns
-- `NTuple{N, Float64}`: All N basis function values
-
-# Examples
-
-```julia
-# Tet10: 10-node quadratic tetrahedron
-topology = Tetrahedron()
-basis = Lagrange{2}()
-xi = Vec(0.25, 0.25, 0.2)
-
-N_all = get_basis_functions(topology, basis, xi)
-# Returns: (N1, N2, ..., N10) as NTuple{10, Float64}
-
-# Access single basis function (simple runtime indexing)
-N_5 = N_all[5]
-
-# Verify partition of unity
-@assert abs(sum(N_all) - 1.0) < 1e-10
-```
-
-# Performance
-- Tet10 (10 nodes): ~3.6 ns, zero allocations
-- Triangle P1 (3 nodes): ~2.5 ns, zero allocations
-
-See also: [`get_basis_derivatives`](@ref), [`get_basis_function`](@ref)
-"""
-function get_basis_functions end
-
-"""
-    get_basis_derivatives(topology::AbstractTopology, basis::AbstractBasis, xi::Vec)
-
-Evaluate all basis function derivatives at parametric point `xi`.
-
-Returns `NTuple{N, Vec{D, Float64}}` where:
-- N is the number of nodes/DOFs
-- D is the spatial dimension
-
-Each derivative is ∇N_i = (∂N_i/∂ξ₁, ∂N_i/∂ξ₂, ..., ∂N_i/∂ξ_D)
-
-# Arguments
-- `topology`: Element topology (e.g., `Tetrahedron()`, `Triangle()`)
-- `basis`: Interpolation scheme (e.g., `Lagrange{1}()`, `Lagrange{2}()`)
-- `xi`: Parametric coordinates as `Vec{D, T}` where D is spatial dimension
-
-# Returns
-- `NTuple{N, Vec{D, Float64}}`: All N basis function gradients
-
-# Examples
-
-```julia
-# Tet10: 10-node quadratic tetrahedron
-topology = Tetrahedron()
-basis = Lagrange{2}()
-xi = Vec(0.25, 0.25, 0.2)
-
-dN_all = get_basis_derivatives(topology, basis, xi)
-# Returns: (∇N1, ∇N2, ..., ∇N10) as NTuple{10, Vec{3, Float64}}
-
-# Access single derivative (simple runtime indexing)
-dN_5 = dN_all[5]  # Vec{3, Float64} gradient for node 5
-
-# Use in assembly (typical pattern)
-for i in 1:10
-    for j in 1:10
-        dNi = dN_all[i]
-        dNj = dN_all[j]
-        K_local[i,j] += dot(dNi, dNj) * detJ  # Simplified
-    end
-end
-```
-
-# Performance
-- Tet10 (10 nodes, 3D): ~6.5 ns, zero allocations ← HOT PATH!
-- Triangle P2 (6 nodes, 2D): ~4.0 ns, zero allocations
-
-This is the critical path for stiffness matrix assembly!
-
-See also: [`get_basis_functions`](@ref), [`get_basis_derivative`](@ref)
-"""
-function get_basis_derivatives end
-
-"""
-    get_basis_function(topology, basis, xi, i::Int)
-
-Convenience function to get a single basis function value.
-
-Equivalent to `get_basis_functions(topology, basis, xi)[i]`.
-
-# Note
-Simple runtime tuple indexing is fastest (~1 ns overhead). No need for
-Val dispatch or @generated functions (benchmarks showed those are 25-300× slower!).
-
-# Examples
-
-```julia
-N_5 = get_basis_function(Tetrahedron(), Lagrange{2}(), xi, 5)
-# Equivalent to:
-N_5 = get_basis_functions(Tetrahedron(), Lagrange{2}(), xi)[5]
-```
-"""
-@inline function get_basis_function(topology::AbstractTopology,
-    basis::AbstractBasis,
-    xi::Vec,
-    i::Int)
-    return get_basis_functions(topology, basis, xi)[i]
-end
-
-"""
-    get_basis_derivative(topology, basis, xi, i::Int)
-
-Convenience function to get a single basis function derivative.
-
-Equivalent to `get_basis_derivatives(topology, basis, xi)[i]`.
-
-# Note
-Simple runtime tuple indexing is fastest (~1 ns overhead).
-
-# Examples
-
-```julia
-dN_5 = get_basis_derivative(Tetrahedron(), Lagrange{2}(), xi, 5)
-# Equivalent to:
-dN_5 = get_basis_derivatives(Tetrahedron(), Lagrange{2}(), xi)[5]
-```
-"""
-@inline function get_basis_derivative(topology::AbstractTopology,
-    basis::AbstractBasis,
-    xi::Vec,
-    i::Int)
-    return get_basis_derivatives(topology, basis, xi)[i]
-end
-
-# ============================================================================
-# Implementation for Lagrange Basis
-# ============================================================================
-# 
-# These implementations are auto-generated for all topology/degree combinations.
-# See src/basis/lagrange_generator.jl for the code generation.
-#
-# Pattern for each topology:
-#
-# @inline function get_basis_functions(::Topology, ::Lagrange{P}, xi::Vec{D,T}) where T
-#     # Compute basis functions
-#     return (N1, N2, ..., Nn)  # NTuple{N, Float64}
-# end
-#
-# @inline function get_basis_derivatives(::Topology, ::Lagrange{P}, xi::Vec{D,T}) where T
-#     # Compute derivatives
-#     return (dN1, dN2, ..., dNn)  # NTuple{N, Vec{D, Float64}}
-# end
-
-# Include auto-generated implementations
-# TODO: Update lagrange_generator.jl to generate get_basis_* functions
-# For now, implementations will be added to lagrange_generated.jl
-
-# ============================================================================
-# Backward Compatibility Bridge (Temporary)
-# ============================================================================
-#
-# These functions bridge to the new API from old API.
-# Will be deprecated after full migration.
-
-"""
-    eval_basis!(basis_type, xi) (DEPRECATED)
-
-**DEPRECATED**: Use `get_basis_functions(topology, basis, xi)` instead.
-
-This function is provided for backward compatibility during migration.
-"""
-function eval_basis! end
-
-"""
-    eval_dbasis!(basis_type, xi) (DEPRECATED)
-
-**DEPRECATED**: Use `get_basis_derivatives(topology, basis, xi)` instead.
-
-This function is provided for backward compatibility during migration.
-"""
-function eval_dbasis! end
-
-# TODO: Add deprecation warnings after new API is implemented
-# @deprecate eval_basis!(args...) get_basis_functions(args...)
-# @deprecate eval_dbasis!(args...) get_basis_derivatives(args...)