diff --git a/src/elements/elements.jl b/src/elements/elements.jl
index aedef99..37ab11e 100644
--- a/src/elements/elements.jl
+++ b/src/elements/elements.jl
@@ -1,963 +1,370 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
 
-"""
-    AbstractFieldSet{N<:Int}
-
-Abstract supertype for all field sets, where `N` is the length of the discrete
-fields (typically is the number of the nodes in element).
-"""
-abstract type AbstractFieldSet{N} end
-
-"""
-    EmptyFieldSet{N} <: AbstractFieldSet{N}
-
-Empty field set used as a default for all elements.
-"""
-struct EmptyFieldSet{N} <: AbstractFieldSet{N}
-end
-
-const DefaultFieldSet = EmptyFieldSet
-
-"""
-    AbstractElement{F, B}
-
-Abstract supertype for all elements.
-
-# Type Parameters
-- `F`: Field type (any type-stable container: NamedTuple, struct, etc.)
-- `B`: Basis function type (e.g., Lagrange{Triangle,1})
-"""
-abstract type AbstractElement{F,B} end
-
-"""
-    Element{N,NIP,F,B} <: AbstractElement{F,B}
-
-Immutable finite element with type-stable fields.
-
-# Type Parameters
-- `N`: Number of nodes (compile-time constant)
-- `NIP`: Number of integration points (compile-time constant)
-- `F`: Field container type (NamedTuple, struct, etc.) - **must be type-stable!**
-- `B`: Basis function type (e.g., Lagrange{Triangle,1})
-
-# Fields
-- `id`: Element identifier
-- `connectivity`: Node IDs as NTuple (zero-cost, immutable)
-- `integration_points`: Integration points as NTuple (zero-cost, immutable)
-- `fields`: Type-stable field container (can be empty tuple `()` if no fields)
-- `basis`: Basis function type instance
-
-# Design Philosophy
-- **Type-stable fields**: `F` parameter ensures compile-time type knowledge
-- **Immutable**: All fields immutable (GPU-compatible, thread-safe)
-- **Zero-allocation**: NTuple for connectivity and IPs
-- **GPU-ready**: Immutable fields can be transferred to GPU without copying
-
-# Examples
-```julia
-# With fields (NamedTuple):
-fields = (E = 210e3, ν = 0.3)
-el = Element(UInt(1), (1,2,3), (), fields, Lagrange{Triangle,1}())
-
-# Without fields (empty tuple):
-el = Element(UInt(1), (1,2,3), (), (), Lagrange{Triangle,1}())
-
-# Access fields (type-stable!):
-E = el.fields.E  # Float64, known at compile time
-```
-"""
-struct Element{N,NIP,F,B} <: AbstractElement{F,B}
-    id::UInt
-    connectivity::NTuple{N,UInt}  # Tuple for zero-cost, compile-time known size
-    integration_points::NTuple{NIP,IP}  # Tuple for zero-cost
-    fields::F  # Type-stable field container (NamedTuple, struct, or ())
-    basis::B
-end
-
-"""
-    Element(basis_type, connectivity; fields=(), id=UInt(0))
-
-Construct element from basis type with optional fields.
-
-# Arguments
-- `basis_type`: Basis function type (e.g., Lagrange{Triangle,1})
-- `connectivity`: Node IDs as tuple or vector
-- `fields`: Optional field container (NamedTuple, struct, or empty tuple)
-- `id`: Element identifier (default: 0)
-
-# Examples
-```julia
-# Simple element without fields:
-element = Element(Lagrange{Triangle,1}, (1, 2, 3))
-
-# Element with material properties:
-element = Element(Lagrange{Triangle,1}, (1, 2, 3), 
-                  fields=(E=210e3, ν=0.3))
-
-# Element with complete specification:
-element = Element(Lagrange{Quadrilateral,1}, (1,2,3,4),
-                  fields=(E=210e3, ν=0.3, thickness=0.01),
-                  id=UInt(42))
-```
-"""
-function Element(::Type{B}, connectivity::NTuple{N,<:Integer};
-    fields::F=(), id::UInt=UInt(0)) where {N,B<:AbstractBasis,F}
-    connectivity_uint = UInt.(connectivity)
-    integration_points = ntuple(i -> IP(UInt(0), 0.0, ()), 0)  # Empty initially
-    basis = B()
-    return Element{N,0,F,B}(id, connectivity_uint, integration_points, fields, basis)
-end
-
-function Element(::Type{B}, connectivity::Vector{<:Integer}; kwargs...) where B<:AbstractBasis
-    return Element(B, (connectivity...,); kwargs...)
-end
-
 # ============================================================================
-# Topology → Basis constructors: Accept topology types, create Lagrange basis
+# Type-Level DOF Count Computation
 # ============================================================================
 
 """
-    Element(::Type{<:AbstractTopology}, connectivity; fields=(), id=UInt(0))
+    ndofs(::Type{K}, ::Type{S}) → Int
 
-Construct element from topology type. Creates linear Lagrange basis automatically.
-
-# Examples
-```julia
-element = Element(Triangle, (1, 2, 3))        # → Lagrange{Triangle,1}
-element = Element(Quadrilateral, (1,2,3,4),   # → Lagrange{Quadrilateral,1}
-                  fields=(E=210e3, ν=0.3))
-element = Element(Tet10, (1,2,3,4,5,6,7,8,9,10))  # → Lagrange{Tet10,2}
-```
-"""
-function Element(::Type{T}, connectivity::NTuple{N,<:Integer}; kwargs...) where {N,T<:AbstractTopology}
-    # Map to base topology for basis functions (Seg3 → Segment, Tri6 → Triangle, etc.)
-    base_topo = get_base_topology(T)
-    # Determine order from number of nodes
-    order = infer_lagrange_order(T, N)
-    # Use NEW API: basis order only (topology passed separately)
-    BasisType = Lagrange{order}
-    return Element(BasisType, connectivity; kwargs...)
-end
-
-function Element(::Type{T}, connectivity::Vector{<:Integer}; kwargs...) where T<:AbstractTopology
-    return Element(T, (connectivity...,); kwargs...)
-end
-
-"""Infer Lagrange order from topology type and number of nodes"""
-function infer_lagrange_order(::Type{T}, n::Int) where T<:AbstractTopology
-    # Get the type name for pattern matching
-    tname = string(nameof(T))
-
-    # Check explicit higher-order topologies first (Tet10, Hex20, Tri6, etc.)
-    if tname == "Tet10" && n == 10
-        return 2
-    elseif tname == "Tri6" && n == 6
-        return 2
-    elseif tname == "Tri7" && n == 7
-        return 3
-    elseif tname == "Quad8" && n == 8
-        return 2  # Serendipity
-    elseif tname == "Quad9" && n == 9
-        return 2  # Full quadratic
-    elseif tname == "Hex20" && n == 20
-        return 2  # Serendipity
-    elseif tname == "Hex27" && n == 27
-        return 2  # Full quadratic
-    elseif tname == "Wedge15" && n == 15
-        return 2
-    elseif tname == "Seg3" && n == 3
-        return 2
-    end
-
-    # For base topologies (and their aliases), infer order from number of nodes
-    if T === Segment || T === Seg2
-        n == 2 && return 1
-        n == 3 && return 2
-    elseif T === Triangle || T === Tri3
-        n == 3 && return 1
-        n == 6 && return 2
-        n == 7 && return 3
-    elseif T === Quadrilateral || T === Quad4
-        n == 4 && return 1
-        n == 8 && return 2  # Serendipity
-        n == 9 && return 2  # Full
-    elseif T === Tetrahedron || T === Tet4
-        n == 4 && return 1
-        n == 10 && return 2
-    elseif T === Hexahedron || T === Hex8
-        n == 8 && return 1
-        n == 20 && return 2  # Serendipity
-        n == 27 && return 2  # Full
-    elseif T === Pyramid || T === Pyr5
-        n == 5 && return 1
-    elseif T === Wedge || T === Wedge6
-        n == 6 && return 1
-        n == 15 && return 2
-    end
-
-    # Default to order 1 (linear)
-    return 1
-end
-# Commenting out until new architecture is fully implemented.
-
-# # Helper: Map topology type to basis type
-# topology_to_basis(::Type{Poi1}) = Poi1Basis
-# topology_to_basis(::Type{Seg2}) = Seg2Basis
-# topology_to_basis(::Type{Seg3}) = Seg3Basis
-# topology_to_basis(::Type{Tri3}) = Tri3Basis
-# topology_to_basis(::Type{Tri6}) = Tri6Basis
-# topology_to_basis(::Type{Tri7}) = Tri7Basis
-# topology_to_basis(::Type{Quad4}) = Quad4Basis
-# topology_to_basis(::Type{Quad8}) = Quad8Basis
-# topology_to_basis(::Type{Quad9}) = Quad9Basis
-# topology_to_basis(::Type{Tet4}) = Tet4Basis
-# topology_to_basis(::Type{Tet10}) = Tet10Basis
-# topology_to_basis(::Type{Hex8}) = Hex8Basis
-# topology_to_basis(::Type{Hex20}) = Hex20Basis
-# topology_to_basis(::Type{Hex27}) = Hex27Basis
-# topology_to_basis(::Type{Pyr5}) = Pyr5Basis
-# topology_to_basis(::Type{Wedge6}) = Wedge6Basis
-# topology_to_basis(::Type{Wedge15}) = Wedge15Basis
-
-# # Shim constructors that forward to basis-based constructors
-# function Element(::Type{T}, connectivity::NTuple{N,<:Integer}) where {N,T<:AbstractTopology}
-#     return Element(topology_to_basis(T), connectivity)
-# end
-
-# function Element(::Type{T}, ::Type{M}, connectivity::NTuple{N,<:Integer}) where {N,M<:AbstractFieldSet,T<:AbstractTopology}
-#     return Element(topology_to_basis(T), M, connectivity)
-# end
-
-# function Element(::Type{T}, connectivity::Vector{<:Integer}) where T<:AbstractTopology
-#     return Element(topology_to_basis(T), connectivity)
-# end
-
-# ============================================================================
-# Immutable update: Returns new element with updated fields
-# ============================================================================
-
-"""
-    update(element::Element, pairs::Pair...) -> Element
-    update(element::Element; kwargs...) -> Element
-
-Create a new element with updated fields (immutable operation).
-
-Since elements are immutable, this function returns a **new** element with merged fields.
-The original element is unchanged.
-
-# Arguments
-- `element`: Original element
-- `pairs`: Field updates as `Pair` objects (`:field => value`)
-- `kwargs`: Field updates as keyword arguments
-
-# Returns
-New `Element` with updated fields (same connectivity, basis, id)
-
-# Examples
-```julia
-# Create element with initial fields:
-el = Element(Lagrange{Triangle,1}, (1,2,3), fields=(E=210e3, ν=0.3))
-
-# Update single field (returns NEW element):
-el2 = update(el, :E => 200e3)  # el.fields.E still 210e3, el2.fields.E is 200e3
-
-# Update multiple fields:
-el3 = update(el, :E => 200e3, :ν => 0.35)
-
-# Keyword syntax:
-el4 = update(el; E=200e3, ν=0.35)
-
-# Add new field:
-el5 = update(el, :temperature => 293.15)
-```
-
-# Performance
-Zero-allocation when fields are type-stable NamedTuples with same field types.
-May allocate if adding new fields with different types.
-
-# See Also
-- Old API: `update!(element, "field", value)` (DEPRECATED - mutates element)
-- New API: `update(element, :field => value)` (returns new element)
-"""
-function update(element::Element{N,NIP,F,B}, pairs::Pair{Symbol,<:Any}...) where {N,NIP,F,B}
-    # Convert pairs to NamedTuple
-    new_fields_dict = Dict{Symbol,Any}(pairs)
-
-    # Merge with existing fields
-    if F === Tuple{} || element.fields === ()
-        # No existing fields, create new NamedTuple
-        merged = NamedTuple(new_fields_dict)
-    else
-        # Merge existing fields with new ones
-        old_fields = element.fields
-        for field_name in fieldnames(typeof(old_fields))
-            if !haskey(new_fields_dict, field_name)
-                new_fields_dict[field_name] = getfield(old_fields, field_name)
-            end
-        end
-        merged = NamedTuple(new_fields_dict)
-    end
-
-    # Create new element with updated fields
-    return Element{N,NIP,typeof(merged),B}(
-        element.id,
-        element.connectivity,
-        element.integration_points,
-        merged,
-        element.basis
-    )
-end
-
-# Keyword argument version
-function update(element::Element{N,NIP,F,B}; kwargs...) where {N,NIP,F,B}
-    pairs = [k => v for (k, v) in kwargs]
-    return update(element, pairs...)
-end
-
-function get_element_id(element::AbstractElement)
-    return element.id
-end
-
-function get_element_type(::AbstractElement{M,T}) where {M,T}
-    return T
-end
-
-function is_element_type(::AbstractElement{M,T}, element_type) where {M,T}
-    return T === element_type
-end
-
-function filter_by_element_type(element_type, elements)
-    return Iterators.filter(element -> is_element_type(element, element_type), elements)
-end
-
-function get_connectivity(element::AbstractElement)
-    return element.connectivity
-end
-
-"""
-    group_by_element_type(elements)
-
-Given a vector of elements, group elements by element type to several vectors.
-Returns a dictionary, where key is the element type and value is a vector
-containing all elements of type `element_type`.
-"""
-function group_by_element_type(elements)
-    eltypes = map(T -> typeof(T), elements)
-    elgroups = Dict(T => T[] for T in eltypes)
-    for element in elements
-        T = typeof(element)
-        push!(elgroups[T], element)
-    end
-    return elgroups
-end
-
-### dfields - dynamically defined fields
-
-# ============================================================================
-# COMPATIBILITY SHIM: Old dfields API → New fields API
-# ============================================================================
-# NOTE: This is temporary until field system is redesigned (Phase 3)
-# Old API: element.dfields (mutable Dict for dynamic fields)
-# New API: element.fields (type-stable NamedTuple/struct, immutable)
-#
-# For now, we map dfields to the same place as sfields (element.fields)
-# since the new API doesn't distinguish between static and dynamic fields.
-
-function has_dfield(element, field_name)
-    # In new API, dfields don't exist as separate from sfields
-    # Check if field_name is in element.fields
-    F = typeof(element.fields)
-    if F === Tuple{}
-        return false  # Empty fields
-    elseif F <: NamedTuple
-        return haskey(element.fields, field_name)
-    else
-        return isdefined(element.fields, field_name)
-    end
-end
-
-function get_dfield(element, field_name)
-    # Return the field value from element.fields
-    F = typeof(element.fields)
-    if F === Tuple{}
-        error("Element has no fields (empty tuple)")
-    end
-    return getfield(element.fields, field_name)
-end
-
-function create_dfield!(element, field_name, field_::AbstractField)
-    # WARNING: In new API, fields are IMMUTABLE!
-    # This function cannot actually create the field in-place
-    # TODO Phase 3: Redesign field system to support mutable time-varying fields
-    @warn "create_dfield!() called but fields are immutable in new API. Field not created." maxlog = 1
-    return
-end
-
-function create_dfield!(element, field_name, field_data)
-    # WARNING: In new API, fields are IMMUTABLE!
-    @warn "create_dfield!() called but fields are immutable in new API. Field not created." maxlog = 1
-    return
-end
-
-function update_dfield!(element, field_name, field_data)
-    # WARNING: In new API, fields are IMMUTABLE!
-    # This function cannot actually update the field in-place
-    # TODO Phase 3: Redesign field system to support mutable time-varying fields
-    @warn "update_dfield!() called but fields are immutable in new API. Field not updated." maxlog = 1
-    return
-end
-
-# A helper function to pick element data from dictionary
-function pick_data_(element, field_data)
-    connectivity = get_connectivity(element)
-    N = length(connectivity)
-    picked_data = ntuple(i -> getindex(field_data, connectivity[i]), N)
-    return picked_data
-end
-
-function update_dfield!(element, field_name, (time, field_data)::Pair{Float64,Dict{Int,V}}) where V
-    # WARNING: In new API, fields are IMMUTABLE!
-    @warn "update_dfield!() called but fields are immutable in new API. Field not updated." maxlog = 1
-    return
-end
-
-function update_dfield!(element, field_name, field_data::Dict{Int,V}) where V
-    # WARNING: In new API, fields are IMMUTABLE!
-    @warn "update_dfield!() called but fields are immutable in new API. Field not updated." maxlog = 1
-    return
-end
-
-function update_dfield!(element, field_name, field_data::Function)
-    # WARNING: In new API, fields are IMMUTABLE!
-    @warn "update_dfield!() called but fields are immutable in new API. Field not updated." maxlog = 1
-    return
-end
-
-function interpolate_dfield(element, field_name, time)
-    # In new API, fields are static (no time variation)
-    # Just return the field value
-    return get_dfield(element, field_name)
-end
-
-### sfields statically defined fields
-
-# ============================================================================
-# COMPATIBILITY SHIM: Old sfields/dfields API → New fields API
-# ============================================================================
-# NOTE: This is temporary until field system is redesigned (Phase 3)
-# Old API: element.sfields (static fields) and element.dfields (dynamic fields)
-# New API: element.fields (type-stable NamedTuple/struct, immutable)
-#
-# Problem: Tests use update_field!() to mutate fields, but new fields are immutable
-# Solution: Make has_sfield/get_sfield work with new API, but warn about mutations
-
-function has_sfield(element, field_name)
-    # In new API, sfields don't exist - but we can check if field_name is in element.fields
-    F = typeof(element.fields)
-    if F === Tuple{}
-        return false  # Empty fields
-    elseif F <: NamedTuple
-        return haskey(element.fields, field_name)
-    else
-        return isdefined(element.fields, field_name)
-    end
-end
-
-function get_sfield(element, field_name)
-    # Return the field value from element.fields
-    F = typeof(element.fields)
-    if F === Tuple{}
-        error("Element has no fields (empty tuple)")
-    end
-    return getfield(element.fields, field_name)
-end
-
-function update_sfield!(element, field_name, field_data)
-    # WARNING: In new API, fields are IMMUTABLE!
-    # This function cannot actually update the field in-place
-    # TODO Phase 3: Redesign field system to support mutable time-varying fields
-    @warn "update_sfield!() called but fields are immutable in new API. Field not updated." maxlog = 1
-    return nothing
-end
-
-function interpolate_sfield(element, field_name, time)
-    # In new API, fields are static (no time variation)
-    # Just return the field value
-    return get_sfield(element, field_name)
-end
-
-### dfield & sfield -- common routines
-
-function has_field(element, field_name)
-    return has_sfield(element, field_name) || has_dfield(element, field_name)
-end
-
-function get_field(element, field_name)
-    if has_sfield(element, field_name)
-        return get_sfield(element, field_name)
-    else
-        return get_dfield(element, field_name)
-    end
-end
-
-function update_field!(element, field_name, field_data)
-    # In new API: fields are immutable, so we can't actually update them
-    # For now, just check if field exists and warn
-    if has_sfield(element, field_name)
-        update_sfield!(element, field_name, field_data)
-    elseif has_dfield(element, field_name)
-        update_dfield!(element, field_name, field_data)
-    else
-        # Field doesn't exist - this is OK in new API where fields are optional
-        # Silently ignore (many tests add fields dynamically that we don't need)
-    end
-end
-
-function interpolate_field(element, field_name::Symbol, time)
-    if has_sfield(element, field_name)
-        return interpolate_sfield(element, field_name, time)
-    elseif has_dfield(element, field_name)
-        return interpolate_dfield(element, field_name, time)
-    else
-        error("Cannot interpolate from field $field_name: no such field.")
-    end
-end
-
-function interpolate(element::AbstractElement, field_name, time)
-    return interpolate_field(element, field_name, time)
-end
-
-function update_field!(elements::Vector{Element}, field_name, field_data)
-    for element in elements
-        update_field!(element, field_name, field_data)
-    end
-end
-
-# Update fields when given a dictionary or time => dictionary:
-# pick data from dictionary diven by the connectivity information of element
-#=
-function update_field!(element::AbstractElement, field::F,
-                       data::Dict{T,V}) where {F<:DVTI,T,V}
-    connectivity = get_connectivity(element)
-    N = length(connectivity)
-    picked_data = ntuple(i -> data[connectivity[i]], N)
-    update_field!(field, picked_data)
-end
-
-function update_field!(element::AbstractElement, field::F,
-                       ddata::Pair{Float64, Dict{T,V}}) where {F<:DVTV,T,V}
-    time, data = ddata
-    connectivity = get_connectivity(element)
-    N = length(connectivity)
-    picked_data = ntuple(i -> data[connectivity[i]], N)
-    update_field!(field, time => picked_data)
-end
-=#
-
-"""
-    interpolate(element, field_name, time)
-
-Interpolate field `field_name` from element at given `time`.
+Compute total number of DOFs for field spec S on topology K.
 
 # Example
-```
-element = Element(Seg2, [1, 2])
-data1 = Dict(1 => 1.0, 2 => 2.0)
-data2 = Dict(1 => 2.0, 2 => 3.0)
-update!(element, "my field", 0.0 => data1)
-update!(element, "my field", 1.0 => data2)
-interpolate(element, "my field", 0.5)
-
-# output
-
-(1.5, 2.5)
-
+```julia
+S = @DOFSet{u::DOF{Displacement{3}, Vertex}}
+ndofs(Tetrahedron{4}, S)  # → 12 (4 nodes × 3 components)
 ```
 """
-function interpolate(element::AbstractElement, field_name::String, time::Float64)
-    field = element[field_name]
-    # If field is an AbstractField, interpolate it
-    if field isa AbstractField
-        result = interpolate(field, time)
-    else
-        # For raw data (Dict, NamedTuple, etc.), it's time-invariant
-        result = field
+@generated function ndofs(::Type{K}, ::Type{S}) where {K, S}
+    field_names = fieldnames(S)
+    total = 0
+    
+    for fname in field_names
+        field_spec = fieldtype(S, fname)  # Tuple{Displacement{3}, Vertex}
+        field_type = field_spec.parameters[1]  # Displacement{3}
+        entity_type = field_spec.parameters[2]  # Vertex
+        
+        # Extract quantity type via trait
+        Q = quantity_type(field_spec)  # Vec{3}
+        
+        # Count entities
+        n_entities = if entity_type === Vertex
+            nnodes(K())
+        elseif entity_type === Edge
+            nedges(K())
+        elseif entity_type === Face
+            nfaces(K())
+        else
+            error("Unsupported entity type: $entity_type")
+        end
+        
+        # Count components per entity
+        n_components = if Q === Float64
+            1
+        elseif Q isa UnionAll && Q.body <: Tensor && Q.body.parameters[1] == 1
+            Q.body.parameters[2]
+        else
+            error("Unsupported quantity type: $Q")
+        end
+        
+        total += n_entities * n_components
     end
-    if isa(result, Dict)
-        connectivity = get_connectivity(element)
-        return tuple((result[i] for i in connectivity)...)
-    else
-        return result
-    end
-end
-
-function info_update_field(elements, field_name, data)
-    nelements = length(elements)
-    @info("Updating field `$field_name` for $nelements elements.")
-end
-
-function info_update_field(elements, field_name, data::Float64)
-    nelements = length(elements)
-    @info("Updating field `$field_name` => $data for $nelements elements.")
+    
+    return total
 end
 
 """
-    update!(elements, field_name, data)
+    AbstractElement{K, P, S, N}
 
-Given a list of elements, field name and data, update field to elements. Data
-is passed directly to the `field`-function.
+Abstract supertype for finite elements following Ciarlet's triple (K, P, Σ).
+
+# Type Parameters
+- `K <: AbstractTopology`: Reference domain
+- `P <: AbstractBasis`: Polynomial space
+- `S`: Field specification (determines Σ functionals)
+- `N::Int`: Total number of DOFs (inferred from S and K)
+
+See `src/elements/README.md` for complete documentation.
+"""
+abstract type AbstractElement{K<:AbstractTopology, P<:AbstractBasis, S<:DOFSet, N} end
+
+"""
+    Element{K, P, S, N}
+
+Finite element implementing Ciarlet's triple (K, P, Σ).
+
+# Type Parameters
+- `K`: Topology (Triangle{3}, Tetrahedron{4}, ...)
+- `P`: Basis (Lagrange{1}, Lagrange{2}, ...)
+- `S`: Field spec with quantity types and entity locations
+- `N::Int`: Total DOF count (automatically inferred from S and K)
+
+# Fields
+- `id::UInt`: Element identifier
+- `dof_indices::NTuple{N,UInt64}`: Flat tuple of global DOF indices
 
 # Examples
-
-Create two elements with topology `Seg2`, one is connecting to nodes (1, 2) and
-the other is connecting to (2, 3). Some examples of updating fields:
-
 ```julia
-elements = [Element(Seg2, [1, 2]), Element(Seg2, [2, 3])]
-X = Dict(1 => 0.0, 2 => 1.0, 3 => 2.0)
-u = Dict(1 => 0.0, 2 => 0.0, 3 => 0.0)
-update!(elements, "geometry", X)
-update!(elements, "displacement", 0.0 => u)
-update!(elements, "youngs modulus", 210.0e9)
-update!(elements, "time-dependent force", 0.0 => 0.0)
-update!(elements, "time-dependent force", 1.0 => 100.0)
+# Single field: 3D displacement (12 DOFs = 4 nodes × 3 components)
+S = @DOFSet{u::DOF{Displacement{3}, Vertex}}
+Element{Tetrahedron{4}, Lagrange{1}, S}(UInt(1), (1,2,3,4,5,6,7,8,9,10,11,12))
+
+# Multi-field: Thermo-mechanical (16 DOFs = 4 T + 12 u)
+S = @DOFSet{T::DOF{Temperature,Vertex}, u::DOF{Displacement{3},Vertex}}
+Element{Tetrahedron{4}, Lagrange{1}, S}(UInt(1), (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16))
 ```
 
-When using dictionaries in definition of fields, key of dictionary corresponds
-to node id, that is, updating field `geometry` in the example above is updating
-values `(0.0, 1.0)` for the first elements and values `(1.0, 2.0)` to the second
-element. For time dependent field, syntax `time => data` is used. If field is
-initialized without time-dependency, it cannot be changed to be time-dependent
-afterwards. If unsure, it's better to initialize field with time dependency.
+See `src/elements/README.md` for comprehensive documentation.
+"""
+struct Element{K<:AbstractTopology, P<:AbstractBasis, S<:DOFSet, N} <: AbstractElement{K,P,S,N}
+    id::UInt
+    dof_indices::NTuple{N,UInt64}
+    
+    # Inner constructor validates N matches spec
+    function Element{K,P,S,N}(id::UInt, dof_indices::NTuple{N,UInt64}) where {K,P,S,N}
+        expected = ndofs(K, S)
+        if N != expected
+            error("Element{$K,$P,$S,$N}: Expected $expected DOFs (from spec), got $N")
+        end
+        return new{K,P,S,N}(id, dof_indices)
+    end
+end
+
+# Outer constructor infers N from tuple length
+function Element{K,P,S}(id::UInt, dof_indices::NTuple{N,UInt64}) where {K,P,S,N}
+    return Element{K,P,S,N}(id, dof_indices)
+end
+
+# Convenience constructor from varargs or vector
+function Element{K,P,S}(id::UInt, dof_indices::UInt64...) where {K,P,S}
+    return Element{K,P,S}(id, dof_indices)
+end
+
+function Element{K,P,S}(id::UInt, dof_indices::AbstractVector{<:Integer}) where {K,P,S}
+    return Element{K,P,S}(id, tuple((UInt64(i) for i in dof_indices)...))
+end
+
+# ============================================================================
+# Type-Level Queries
+# ============================================================================
 
 """
-function update!(elements, field_name, data)
-    info_update_field(elements, field_name, data)
-    for element in elements
-        update!(element, field_name, data)
-    end
+    topology_type(::Element{K,P,S,N}) → Type{K}
+
+Extract topology type K from element.
+"""
+topology_type(::Element{K,P,S,N}) where {K,P,S,N} = K
+topology_type(::Type{Element{K,P,S,N}}) where {K,P,S,N} = K
+
+"""
+    basis_type(::Element{K,P,S,N}) → Type{P}
+
+Extract basis type P from element.
+"""
+basis_type(::Element{K,P,S,N}) where {K,P,S,N} = P
+basis_type(::Type{Element{K,P,S,N}}) where {K,P,S,N} = P
+
+"""
+    dof_type(::Element{K,P,S,N}) → Type{S}
+
+Extract DOF specification type S from element.
+"""
+dof_type(::Element{K,P,S,N}) where {K,P,S,N} = S
+dof_type(::Type{Element{K,P,S,N}}) where {K,P,S,N} = S
+
+# ============================================================================
+# Local-Global DOF Mapping for Coupled Assembly
+# ============================================================================
+
+"""
+    local_dof_count(elem::Element) → Int
+
+Total number of local DOFs for this element (sum over all fields).
+"""
+@inline function local_dof_count(elem::Element{K,P,S,N}) where {K,P,S,N}
+    return N  # Now directly available as type parameter!
 end
 
+"""
+    global_dof_indices(elem::Element) → Vector{UInt64}
 
-## Interpolate fields in spatial direction
+Flattened vector of global DOF indices for this element.
 
-const ConstantField = Union{DCTI,DCTV}
-const VariableFields = Union{DVTV,DVTI}
-const DictionaryFields = Union{DVTVd,DVTId}
-
-function interpolate_field(::AbstractElement, field::ConstantField, ip, time)
-    return interpolate_field(field, time)
+See `src/elements/README.md` for assembly patterns.
+"""
+function global_dof_indices(elem::Element)
+    return collect(elem.dof_indices)  # NTuple → Vector
 end
 
-function interpolate_field(element::AbstractElement, field::VariableFields, ip, time)
-    data = interpolate_field(field, time)
-    basis = get_basis(element, ip, time)
-    N = length(basis)
-    return sum(data[i] * basis[i] for i = 1:N)
+"""
+    local_to_global_map(elem::Element) → NTuple{N,UInt64}
+
+Mapping from local DOF index to global DOF index.
+`global_dof = map[local_dof]` where `local_dof ∈ 1:N`.
+
+Returns tuple (not Vector) for type stability and compiler optimization.
+Used for coupled assembly. See `src/elements/README.md`.
+"""
+@inline function local_to_global_map(elem::Element{K,P,S,N}) where {K,P,S,N}
+    return elem.dof_indices  # Already flat!
 end
 
-function interpolate_field(element::AbstractElement, field::DictionaryFields, ip, time)
-    data = interpolate_field(field, time)
-    basis = element(ip, time)
-    N = length(element)
-    c = get_connectivity(element)
-    return sum(data[c[i]] * basis[i] for i = 1:N)
-end
+# ============================================================================
+# Compile-Time Helper Functions for @generated field_dof_range
+# ============================================================================
 
-function interpolate_field(::AbstractElement, field::CVTV, ip, time)
-    return field(ip, time)
-end
-
-function interpolate(element::AbstractElement, field_name, ip, time)
-    field = get_field(element, field_name)
-    interpolate_field(element, field, ip, time)
-end
-
-
-## Other stuff
-
-# TEMPORARY: Commented out during basis refactoring - uses old Lagrange{T,P} API
-# TODO: Rewrite to use new Lagrange{P} API with topology passed separately
-#=
-# Helper: Create topology instance from basis type
-function _create_topology_instance(::Type{Lagrange{T,P}}) where {T,P}
-    N_nodes = nnodes(Lagrange{T,P}())
-    if T <: Segment
-        return Segment{N_nodes}()
-    elseif T <: Triangle
-        return Triangle()
-    elseif T <: Quadrilateral
-        return Quadrilateral()
-    elseif T <: Tetrahedron
-        return Tetrahedron()
-    elseif T <: Hexahedron
-        return Hexahedron()
-    elseif T <: Wedge
-        return Wedge()
-    elseif T <: Pyramid
-        return Pyramid()
+# Helper: Compute ndofs at compile time
+function _compile_time_ndofs(@nospecialize(field_type), @nospecialize(topology_type))
+    # Handle DOF{FieldType, EntityType} format (new format)
+    if field_type isa DataType && field_type <: DOF && length(field_type.parameters) == 2
+        FieldType = field_type.parameters[1]  # e.g., Displacement{3}
+        E = field_type.parameters[2]          # e.g., Vertex
+        
+        # Extract quantity type via trait (handles Displacement{3} → Vec{3})
+        Q = quantity_type(field_type)
+        
+        # Number of DOFs = dof_per_entity * number_of_entities
+        return _dof_per_entity(Q) * _count_entities_compiletime(topology_type, E)
+    # Handle Tuple{FieldType, E} format (legacy format)
+    elseif field_type isa DataType && field_type <: Tuple && length(field_type.parameters) == 2
+        FieldType = field_type.parameters[1]  # Could be Displacement{3} or Vec{3}
+        E = field_type.parameters[2]
+        
+        # Extract quantity type via trait (handles both field types and quantity types)
+        Q = quantity_type(field_type)
+        
+        # Number of DOFs = dof_per_entity * number_of_entities
+        return _dof_per_entity(Q) * _count_entities_compiletime(topology_type, E)
     else
-        error("Unknown topology type: $T")
+        error("Cannot compute ndofs for field type $field_type (expected DOF{...} or Tuple{...})")
     end
 end
 
-# Jacobian computation using new API
-# Handles embedded elements (e.g., 1D element in 2D/3D space)
-function jacobian(B::Lagrange{T,P}, X::Vector{<:Vec}, xi::Vec) where {T,P}
-    topo_instance = _create_topology_instance(Lagrange{T,P})
-    dN_dξ_tuple = get_basis_derivatives(topo_instance, B, xi)
-
-    # Handle embedding: X can be Vec{D_phys} while dN_dξ is Vec{D_param}
-    # where D_phys > D_param (e.g., 1D element in 2D/3D space)
-    dim_physical = length(first(X))
-    dim_parametric = length(xi)
-
-    # Build Jacobian matrix manually for embedding case
-    # J[i,j] = ∂X_i/∂ξ_j = sum_k X_k[i] * dN_k/dξ_j
-    J_data = zeros(dim_physical, dim_parametric)
-    @inbounds for k in 1:length(X)
-        for i in 1:dim_physical
-            for j in 1:dim_parametric
-                J_data[i, j] += X[k][i] * dN_dξ_tuple[k][j]
-            end
+function _dof_per_entity(@nospecialize(Q))
+    # Use dof_size which handles all quantity types properly (Displacement{3}, Vec{3}, Float64, UnionAll, etc.)
+    # This is the most robust approach
+    try
+        return dof_size(Q)
+    catch e
+        # Fallback for specific cases if dof_size fails
+        if Q === Float64
+            return 1
+        else
+            error("Cannot determine dof_size for quantity type $Q: $e")
         end
     end
-
-    return J_data
 end
 
-function get_basis(element::AbstractElement{M,B}, ip, ::Any) where {M,B<:Lagrange}
-    # Handle both raw coordinates (Tuple) and IP struct
-    coords = isa(ip, IP) ? ip.coords : ip
-    T = typeof(first(coords))
-    # Convert to Vec for Tensors.jl compatibility
-    xi = Vec{length(coords),T}(coords)
-
-    # Create topology instance once per basis type (will be inlined/constant folded)
-    topo_instance = _create_topology_instance(B)
-
-    # Call new API
-    N_tuple = get_basis_functions(topo_instance, B(), xi)
-    # Return as row matrix for compatibility with old code
-    return reshape(collect(N_tuple), 1, length(element))
-end
-
-function get_dbasis(element::AbstractElement{M,B}, ip, ::Any) where {M,B<:Lagrange}
-    # Handle both raw coordinates (Tuple) and IP struct
-    coords = isa(ip, IP) ? ip.coords : ip
-    T = typeof(first(coords))
-    # Convert to Vec for Tensors.jl compatibility
-    xi = Vec{length(coords),T}(coords)
-
-    # Create topology instance once per basis type (will be inlined/constant folded)
-    topo_instance = _create_topology_instance(B)
-
-    # Call new API
-    dN_tuple = get_basis_derivatives(topo_instance, B(), xi)
-    # Return as Vector for compatibility with old code
-    return collect(dN_tuple)
-end
-=#
-
-# END COMMENTED OUT SECTION
-
-function (element::Element)(ip, time::Float64=0.0)
-    return get_basis(element, ip, time)
-end
-
-#"""
-#Examples
-#julia> el = Element(Quad4, [1, 2, 3, 4]);
-#julia> el([0.0, 0.0], 0.0, 1)
-#1x4 Array{Float64,2}:
-# 0.25  0.25  0.25  0.25
-#julia> el([0.0, 0.0], 0.0, 2)
-#2x8 Array{Float64,2}:
-# 0.25  0.0   0.25  0.0   0.25  0.0   0.25  0.0
-# 0.0   0.25  0.0   0.25  0.0   0.25  0.0   0.25
-#"""
-function (element::Element)(ip, time::Float64, dim::Int)
-    dim == 1 && return get_basis(element, ip, time)
-    Ni = vec(get_basis(element, ip, time))
-    N = zeros(dim, length(element) * dim)
-    for i = 1:dim
-        N[i, i:dim:end] += Ni
-    end
-    return N
-end
-
-function (element::Element)(ip, time, ::Type{Val{:Jacobian}})
-    X_dict = element("geometry", time)
-    # Convert Dict to Vector{Vec} in connectivity order
-    # (Dict iteration is unordered, must use element.connectivity)
-    X = [Vec(X_dict[node_id]...) for node_id in element.connectivity]
-    # Convert ip to Vec - handle both Tuple and IntegrationPoint
-    if isa(ip, Tuple)
-        xi = Vec(ip)
+function _count_entities_compiletime(@nospecialize(K), @nospecialize(E))
+    # This must match count_entities(topology, entity_type) at runtime
+    # K is a TYPE (e.g., Tet4), not an instance
+    if E === Vertex
+        return nnodes(K)  # nnodes accepts Type
+    elseif E === Edge
+        return nedges(K)  # nedges accepts Type
+    elseif E === Face
+        return nfaces(K)  # nfaces accepts Type
+    elseif E === Cell
+        return 1  # One cell per element
     else
-        xi = Vec(ip.coords)
-    end
-    J = jacobian(element.basis, X, xi)
-    return J
-end
-
-function (element::Element)(ip, time::Float64, ::Type{Val{:detJ}})
-    J = element(ip, time, Val{:Jacobian})
-    n, m = size(J)
-    if n == m  # volume element
-        return det(J)
-    end
-    # For embedded elements (1D in 2D/3D, 2D in 3D):
-    # detJ = || ∂X/∂ξ || for 1D elements
-    # detJ = || ∂X/∂ξ₁ × ∂X/∂ξ₂ || for 2D elements
-    JT = transpose(J)
-    if m == 1  # 1D element (boundary of 2D or 3D), J is n×1, JT is 1×n
-        return norm(JT)
-    else  # 2D element (manifold on 3D problem), J is 3×2, JT is 2×3
-        return norm(cross(JT[:, 1], JT[:, 2]))
+        error("Unknown entity type $E")
     end
 end
 
-function (element::Element)(ip, time::Float64, ::Type{Val{:Grad}})
-    J = element(ip, time, Val{:Jacobian})
-    return inv(J) * get_dbasis(element, ip, time)
-end
+# ============================================================================
+# Local DOF Range Computation (COMPILE-TIME via @generated)
+# ============================================================================
 
-function (element::Element)(field_name::String, ip, time::Float64, ::Type{Val{:Grad}})
-    X = element("geometry", time)
-    u = element(field_name, time)
-    return grad(element.properties, u, X, ip)
-end
-
-
-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
-    # Otherwise get default integration points based on basis type
-    ips = get_integration_points_from_basis(B)
-    return tuple([IP(UInt(i), w, xi) for (i, (w, xi)) in enumerate(ips)]...)
-end
-
-""" This is a special case, temporarily change order
-of integration scheme mainly for mass matrix.
 """
-function get_integration_points(element::AbstractElement{M,B}, change_order::Int) where {M,B}
-    ips = get_integration_points_from_basis(B, change_order)
-    # Convert from new API format (weight, Vec{D}) to old IP structs
-    # Vec{D} needs to be converted to Tuple for IP constructor
-    return tuple([IP(UInt(i), w, Tuple(xi)) for (i, (w, xi)) in enumerate(ips)]...)
-end
+    field_dof_range(elem::Element, field::Symbol) → UnitRange{Int}
 
-# TEMPORARY: Commented out during basis refactoring - uses old Lagrange{T,P} API
-#=
-# Helper function to map Lagrange basis types to new integration points API
-function get_integration_points_from_basis(::Type{Lagrange{T,P}}, order::Int=0) where {T,P}
-    # Map topology type to base topology for integration points
-    # (Seg2, Seg3 → Segment; Tri3, Tri6, Tri7 → Triangle, etc.)
-    base_topology = get_base_topology(T)
-
-    # Map polynomial order to Gauss quadrature order
-    # For polynomial order P, we need at least (P+1)/2 integration order
-    # With order parameter, we can increase accuracy
-    gauss_order = max(P, P + order)
-
-    # Return integration points using new zero-allocation API
-    # Note: get_gauss_points! expects Type{Gauss{N}}, not instance
-    if gauss_order == 1
-        return get_gauss_points!(base_topology, Gauss{1})
-    elseif gauss_order == 2
-        return get_gauss_points!(base_topology, Gauss{2})
-    elseif gauss_order == 3
-        return get_gauss_points!(base_topology, Gauss{3})
-    elseif gauss_order == 4
-        return get_gauss_points!(base_topology, Gauss{4})
-    elseif gauss_order == 5
-        return get_gauss_points!(base_topology, Gauss{5})
-    else
-        error("Gauss quadrature order $gauss_order not supported for topology $T")
-    end
-end
-=#
-# END COMMENTED OUT SECTION (get_integration_points_from_basis)
-
-# Map topology types to their base forms for integration
-# NOTE: Restored from commented section - needed by Element constructor
-function get_base_topology(::Type{T}) where T
-    name = string(nameof(T))
-    # 1D elements
-    if startswith(name, "Seg") || name == "Segment"
-        return Segment
-        # 2D elements
-    elseif startswith(name, "Tri") || name == "Triangle"
-        return Triangle
-    elseif startswith(name, "Quad") || name == "Quadrilateral"
-        return Quadrilateral
-        # 3D elements
-    elseif startswith(name, "Tet") || name == "Tetrahedron"
-        return Tetrahedron
-    elseif startswith(name, "Hex") || name == "Hexahedron"
-        return Hexahedron
-    elseif startswith(name, "Wedge")
-        return Wedge
-    elseif startswith(name, "Pyr") || name == "Pyramid"
-        return Pyramid
-    else
-        error("Unknown topology type: $T")
-    end
-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.
+Local DOF range for a specific field. Computed at compile time via @generated.
+See `src/elements/README.md` for usage examples.
 """
-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. """
-function get_local_coordinates(element::AbstractElement, X::Vector, time::Float64; max_iterations=10, tolerance=1.0e-6)
-    haskey(element, "geometry") || error("element geometry not defined, cannot calculate inverse isoparametric mapping")
-    dim = size(element, 1)
-    dim == length(X) || error("manifolds not supported.")
-    xi = zeros(dim)
-    dX = element("geometry", xi, time) - X
-    for i = 1:max_iterations
-        J = element(xi, time, Val{:Jacobian})'
-        xi -= J \ dX
-        dX = element("geometry", xi, time) - X
-        norm(dX) < tolerance && return xi
+@generated function field_dof_range(::Element{K,P,S,N}, field::Symbol) where {K,P,S,N}
+    # This runs at COMPILE TIME!
+    # S is the NamedTuple type containing field specifications
+    
+    if S <: NamedTuple
+        # Multi-field case
+        field_types = S.parameters[2]  # Tuple of field types
+        field_names = fieldnames(S)
+        
+        # Compute offset for each field at compile time
+        offset = 0
+        field_ranges = Expr(:block)
+        
+        for (i, fname) in enumerate(field_names)
+            ftype = field_types.parameters[i]
+            n = _compile_time_ndofs(ftype, K)
+            range_expr = :($offset+1:$offset+$n)
+            
+            # Generate: if field === :fname return range_expr end
+            push!(field_ranges.args, quote
+                if field === $(QuoteNode(fname))
+                    return $range_expr
+                end
+            end)
+            
+            offset += n
+        end
+        
+        # Add error case
+        push!(field_ranges.args, :(error("Field ", field, " not found in element type $S")))
+        
+        return field_ranges
+    else
+        # Single-field case (S <: AbstractDOF)
+        n = _compile_time_ndofs(S, K)
+        return :(return 1:$n)
     end
-    debug("get_local_coordinates", X, dX, xi)
-    error("Unable to find inverse isoparametric mapping for element $element for X = $X")
 end
 
-""" Test is X inside element. """
-function inside(element::AbstractElement{M,B}, X, time) where {M,B}
-    xi = get_local_coordinates(element, X, time)
-    return inside(B, xi)
+# ============================================================================
+# Local-to-Global Mapping (Type-Stable Tuple Version)
+# ============================================================================
+
+# ============================================================================
+# Type Extraction (previously defined above)
+# ============================================================================
+
+# These were defined earlier but are here for reference
+# topology_type, basis_type, dof_type already defined above
+
+# ============================================================================
+# Element Queries
+# ============================================================================
+
+"""
+    element_id(elem::Element) → UInt
+
+Get element ID (index in mesh).
+"""
+element_id(elem::Element) = elem.id
+
+"""
+    element_dofs(elem::Element) → NTuple{N,UInt64}
+
+Get all global DOF indices as flat tuple.
+"""
+element_dofs(elem::Element) = elem.dof_indices
+
+"""
+    element_dofs(elem::Element, field::Symbol) → Tuple
+
+Get global DOF indices for specific field by extracting from flat tuple.
+
+# Example
+```julia
+element_dofs(elem, :T)  # Extracts T indices from flat tuple
+element_dofs(elem, :u)  # Extracts u indices from flat tuple
+```
+"""
+function element_dofs(elem::Element{K,P,S,N}, field::Symbol) where {K,P,S,N}
+    range = field_dof_range(elem, field)
+    return elem.dof_indices[range]
 end
 
-## Convenience functions
+"""
+    n_element_dofs(elem::Element) → Int
 
-# element("displacement", 0.0)
-function (element::Element)(field_name::String, time::Float64)
-    return interpolate(element, field_name, time)
+Get total number of DOFs for this element (all fields).
+"""
+n_element_dofs(elem::Element{K,P,S,N}) where {K,P,S,N} = N
+
+"""
+    nnodes(::Element{K,P,S,N}) → Int
+
+Get number of nodes from topology.
+"""
+nnodes(::Element{K,P,S,N}) where {K,P,S,N} = nnodes(K)
+nnodes(::Type{Element{K,P,S,N}}) where {K,P,S,N} = nnodes(K)
+
+# ============================================================================
+# Display
+# ============================================================================
+
+function Base.show(io::IO, elem::Element{K,P,S,N}) where {K,P,S,N}
+    print(io, "Element{$K, $P, $S}(id=$(elem.id), ndofs=$N)")
 end
-
-# element("displacement", (0.0, 0.0), 0.0)
-function (element::Element)(field_name::String, ip, time::Float64)
-    return interpolate(element, field_name, ip, time)
-end
-
-# OLD: Uses BasisInfo which was defined in math.jl (commented out)
-# function element_info!(bi::BasisInfo{T}, element::AbstractElement{M,T}, ip, time) where {M,T}
-#     X = interpolate(element, "geometry", time)
-#     eval_basis!(bi, X, ip)
-#     return bi.J, bi.detJ, bi.N, bi.grad
-# end