diff --git a/docs/book/element_architecture.md b/docs/book/element_architecture.md
index 689fe06..f096ae1 100644
--- a/docs/book/element_architecture.md
+++ b/docs/book/element_architecture.md
@@ -21,19 +21,21 @@ What is a finite element? This seemingly simple question has profound implicatio
 
 A finite element is fundamentally composed of **four independent concerns**:
 
-### 1. Topology (Connectivity/Graph Theory)
+### 1. Topology (Reference Element Geometry)
 
-**What it is:** The combinatorial structure of how nodes connect to form an element.
+**What it is:** The **geometric shape** of the reference element in parametric space.
 
-- **Examples:** `Tri3` (3-node triangle), `Quad4` (4-node quadrilateral), `Tet10` (10-node tetrahedron)
-- **Properties:** Number of nodes, edges, faces; reference element geometry
-- **Mathematics:** Graph theory, combinatorics
-- **Rarely changes:** Topology is a mathematical object, not implementation-dependent
+- **Examples:** `Triangle`, `Quadrilateral`, `Tetrahedron`, `Hexahedron`, `Pyramid`, `Wedge`
+- **Properties:** Dimension, edges, faces; parametric domain
+- **Mathematics:** Differential geometry, topology
+- **Key insight:** Topology is **pure geometry**, independent of node count or DOF placement
 
-**Reference element:** The element in parametric coordinates $\xi \in [-1, 1]^d$
+**Critical:** Topology does NOT specify node count! That's determined by the interpolation scheme.
+
+**Reference element:** The element in parametric coordinates $\xi \in \Omega_{ref}$
 
 ```text
-Tri3 reference element:
+Triangle reference element (parametric domain):
      η
      ^
      |
@@ -43,28 +45,54 @@ Tri3 reference element:
      |      \
      +---------> ξ
   (0,0)    (1,0)
+
+Same topology works for:
+- 3 nodes (linear, P1)
+- 6 nodes (quadratic, P2) 
+- 10 nodes (cubic, P3)
+- Edge DOFs (Nédélec)
+- Face DOFs (Raviart-Thomas)
 ```
 
-### 2. Interpolation (Basis Functions)
+### 2. Interpolation (Basis Functions + DOF Placement)
 
-**What it is:** How to interpolate field values between nodes.
+**What it is:** How to interpolate field values AND where degrees of freedom live.
 
-- **Examples:** Lagrange polynomials, hierarchical polynomials, NURBS
-- **Properties:** Polynomial order, continuity, partition of unity
+- **Examples:** Lagrange (nodal DOFs), Nédélec (edge DOFs), Raviart-Thomas (face DOFs)
+- **Properties:** Polynomial order, continuity, DOF location, partition of unity
 - **Mathematics:** Approximation theory, functional analysis
 - **Can vary:** Same topology with different interpolation schemes
 
 **Interpolation formula:** $u(\xi) = \sum_{i=1}^n N_i(\xi) u_i$
 
-where $N_i(\xi)$ are basis functions and $u_i$ are nodal values.
+where $N_i(\xi)$ are basis functions and $u_i$ are DOF values (not necessarily at nodes!).
 
-**Key property:** Basis functions are **independent of topology** (mostly).
+**Critical insight:** Interpolation determines BOTH polynomial degree AND node count:
+- `Lagrange{Triangle, 1}` → P1 → 3 nodes (vertices only)
+- `Lagrange{Triangle, 2}` → P2 → 6 nodes (vertices + edge midpoints)
+- `Lagrange{Triangle, 3}` → P3 → 10 nodes (vertices + edges + interior)
+- `Nedelec{Triangle, 1}` → Edge elements → 3 DOFs on edges, NOT at nodes!
 
-- Linear Lagrange on `Tri3`: $N_1 = 1 - \xi - \eta$, $N_2 = \xi$, $N_3 = \eta$
-- Linear Lagrange on `Quad4`: $N_1 = (1-\xi)(1-\eta)/4$, ...
-- Hierarchical on `Tri3`: $N_1 = 1 - \xi - \eta$, $N_2 = \xi(1-\xi-\eta)$, ...
+**Node count vs DOF count:**
+- **Nodes:** Geometric points for element connectivity (graph structure)
+- **DOFs:** Where unknowns live (can be at nodes, edges, faces, interior)
 
-**Important:** Interpolation scheme determines polynomial order, NOT topology.
+```julia
+# Electromagnetics: DOFs on edges, not nodes
+element = Element(Triangle(), Nedelec{Triangle, 1}(), Gauss{2}(), (1,2,3))
+nnodes(element) # → 3 (vertices for connectivity)
+ndofs(element)  # → 3 (one DOF per edge)
+
+# Standard mechanics: DOFs at nodes
+element = Element(Triangle(), Lagrange{Triangle, 1}(), Gauss{1}(), (1,2,3))
+nnodes(element) # → 3
+ndofs(element)  # → 3 (coincide for nodal elements)
+
+# Quadratic: More nodes than linear
+element = Element(Triangle(), Lagrange{Triangle, 2}(), Gauss{3}(), (1,2,3,4,5,6))
+nnodes(element) # → 6 (vertices + edge midpoints)
+ndofs(element)  # → 6
+```
 
 ### 3. Integration (Quadrature Rules)
 
@@ -106,94 +134,124 @@ fields = Dict(
 
 **Key property:** Fields are **completely independent** of topology, interpolation, and integration.
 
-## The Anti-Pattern: Abaqus's Mistake
+## The Anti-Pattern: Code Aster (and Abaqus)
 
-Abaqus (and many commercial codes) conflate these concerns, leading to **combinatorial explosion**:
+Commercial FEM codes conflate topology, node count, and interpolation, leading to **hardcoded combinatorial explosion**:
 
-### Hexahedral Element Examples
+### The Hardcoded Node Count Anti-Pattern
 
-| Element Type | Topology | Interpolation | Integration | Modes |
-|--------------|----------|---------------|-------------|-------|
-| `C3D8`       | Hex8     | Linear        | Full (2×2×2)| None  |
-| `C3D8R`      | Hex8     | Linear        | Reduced (1) | None  |
-| `C3D8I`      | Hex8     | Linear        | Full (2×2×2)| Incompatible |
-| `C3D20`      | Hex20    | Quadratic     | Full (3×3×3)| None  |
-| `C3D20R`     | Hex20    | Quadratic     | Reduced (2×2×2) | None |
-| `C3D20RH`    | Hex20    | Quadratic     | Reduced     | Hybrid |
-| `C3D27`      | Hex27    | Quadratic     | Full (3×3×3)| None  |
-| `C3D27R`     | Hex27    | Quadratic     | Reduced     | None  |
+| Element Type | Topology | Node Count | Polynomial Order | Integration |
+|--------------|----------|------------|------------------|-------------|
+| `TRIA3`      | Triangle | **3** (hardcoded!) | P1 | Default |
+| `TRIA6`      | Triangle | **6** (hardcoded!) | P2 | Default |
+| `QUAD4`      | Quadrilateral | **4** | Q1 | Full |
+| `QUAD8`      | Quadrilateral | **8** | Q2 Serendipity | Full |
+| `QUAD9`      | Quadrilateral | **9** | Q2 | Full |
+| `TETRA4`     | Tetrahedron | **4** | P1 | Default |
+| `TETRA10`    | Tetrahedron | **10** | P2 | Default |
+| `HEXA8`      | Hexahedron | **8** | Q1 | Full (2×2×2) |
+| `HEXA20`     | Hexahedron | **20** | Q2 Serendipity | Full |
+| `HEXA27`     | Hexahedron | **27** | Q2 | Full (3×3×3) |
 
-**Result:** 8 different "element types" for what is fundamentally **one topology** with different choices for interpolation and integration!
+**The fundamental mistake:** Node count is **hardcoded into the type name**, when it should be a *consequence* of:
+1. Topology (geometric shape)
+2. Interpolation scheme (polynomial degree)
+
+**Result:** Cannot use Triangle with edge DOFs (Nédélec), cannot add interior DOFs for pressure, cannot use hierarchical basis with same topology.
 
 ### The Problem with Conflation
 
 ```c
-// Abaqus-style (pseudo-code)
-class C3D8 {
-    // Everything mixed together
-    Node nodes[8];
+// Code Aster style (pseudo-code)
+class TRIA3 {
+    // Topology, node count, interpolation all mixed
+    Node nodes[3];  // Hardcoded!
     void stiffness_matrix() {
-        // Hardcoded: 8 nodes, linear shape functions, 2×2×2 Gauss
+        // Hardcoded: 3 nodes, P1 shape functions, default integration
     }
 };
 
-class C3D8R {
-    // Almost identical code, but different integration
-    Node nodes[8];
+class TRIA6 {
+    // Almost identical code for same topology!
+    Node nodes[6];  // Different hardcoded count
     void stiffness_matrix() {
-        // Hardcoded: 8 nodes, linear shape functions, 1 point
+        // Hardcoded: 6 nodes, P2 shape functions, more integration points
     }
 };
 
-// Now need C3D8I, C3D20, C3D20R, ... → code duplication nightmare
+// Now need TRIA7, TRIA10 (cubic), QUAD4, QUAD8, QUAD9, ... → explosion
 ```
 
 **Issues:**
 
-- ❌ Code duplication (each element type reimplements similar logic)
-- ❌ Combinatorial explosion (n topologies × m interpolations × k integrations)
-- ❌ Maintenance nightmare (bug fix must be repeated in all variants)
-- ❌ Cannot mix-and-match (user stuck with pre-defined combinations)
-- ❌ No compile-time optimization (runtime dispatch on element type)
+- ❌ **Hardcoded node count** prevents using same topology with different basis
+- ❌ **Cannot use edge/face DOFs** (Nédélec, Raviart-Thomas for electromagnetics)
+- ❌ **Cannot add interior DOFs** (pressure in mixed formulations)
+- ❌ **Code duplication** (TRIA3 and TRIA6 nearly identical except node count)
+- ❌ **Combinatorial explosion** (n shapes × m node-counts × k integrations × ...)
+- ❌ **Maintenance nightmare** (bug fix must be repeated in all variants)
+- ❌ **No runtime dispatch** possible (everything statically hardcoded)
 
 ## JuliaFEM's Approach: Composition Over Conflation
 
 ### Separation of Concerns
 
 ```julia
-# 1. Define topology (reference element)
+# 1. Define topology (pure geometry, NO node count)
 abstract type AbstractTopology end
-struct Tri3 <: AbstractTopology 
-    nnodes::Int = 3
-    dim::Int = 2
-end
-struct Quad4 <: AbstractTopology 
-    nnodes::Int = 4
-    dim::Int = 2
-end
-struct Hex8 <: AbstractTopology 
-    nnodes::Int = 8
-    dim::Int = 3
-end
 
-# 2. Define interpolation schemes
+struct Point <: AbstractTopology end
+struct Segment <: AbstractTopology end
+struct Triangle <: AbstractTopology end
+struct Quadrilateral <: AbstractTopology end
+struct Tetrahedron <: AbstractTopology end
+struct Hexahedron <: AbstractTopology end
+struct Pyramid <: AbstractTopology end
+struct Wedge <: AbstractTopology end  # Prism
+
+# Properties come from topology itself
+dim(::Triangle) = 2
+dim(::Tetrahedron) = 3
+
+# 2. Define interpolation schemes (determines node count AND DOF placement)
 abstract type AbstractBasis end
-struct Lagrange{P} <: AbstractBasis end  # P = polynomial order
-struct Hierarchical{P} <: AbstractBasis end
-struct NURBS{P} <: AbstractBasis end
+
+# Lagrange family: nodal DOFs, polynomial degree P
+struct Lagrange{T<:AbstractTopology, P} <: AbstractBasis end
+
+# Serendipity: reduced node count (no center nodes)
+struct Serendipity{T<:AbstractTopology, P} <: AbstractBasis end
+
+# Nédélec: edge DOFs for H(curl) spaces (electromagnetics)
+struct Nedelec{T<:AbstractTopology, P} <: AbstractBasis end
+
+# Raviart-Thomas: face DOFs for H(div) spaces (fluid flow)
+struct RaviartThomas{T<:AbstractTopology, P} <: AbstractBasis end
+
+# Hermite: nodal values + derivatives
+struct Hermite{T<:AbstractTopology, P} <: AbstractBasis end
+
+# Node count is DERIVED from topology + basis:
+nnodes(::Lagrange{Triangle, 1}) = 3   # P1: vertices only
+nnodes(::Lagrange{Triangle, 2}) = 6   # P2: vertices + edge midpoints
+nnodes(::Lagrange{Triangle, 3}) = 10  # P3: vertices + edges + interior
+
+nnodes(::Lagrange{Quadrilateral, 1}) = 4   # Q1: corners
+nnodes(::Lagrange{Quadrilateral, 2}) = 9   # Q2: full tensor product
+nnodes(::Serendipity{Quadrilateral, 2}) = 8  # Q2 without center
 
 # 3. Define integration rules
 abstract type AbstractIntegration end
-struct Gauss{N} <: AbstractIntegration end  # N = number of points
+struct Gauss{N} <: AbstractIntegration end  # N = order (not point count!)
 struct Lobatto{N} <: AbstractIntegration end
 struct Reduced <: AbstractIntegration end
 
-# 4. Element composes all three
-struct Element{T <: AbstractTopology, B <: AbstractBasis, I <: AbstractIntegration, N}
+# 4. Element composes all three + connectivity
+struct Element{T<:AbstractTopology, B<:AbstractBasis, I<:AbstractIntegration, N}
     topology::T
     basis::B
     integration::I
-    connectivity::NTuple{N, UInt}
+    connectivity::NTuple{N, UInt}  # Tuple, not Vector!
     fields::Dict{Symbol, Any}  # TODO: Type-stable structure
 end
 ```
@@ -202,15 +260,33 @@ end
 
 ```julia
 # Create element by composing concerns
-topology = Tri3()
-basis = Lagrange{1}()       # Linear interpolation
-integration = Gauss{3}()    # 3-point Gauss quadrature
+topology = Triangle()
+basis = Lagrange{Triangle, 1}()     # Linear P1 → 3 nodes
+integration = Gauss{2}()             # Order 2 (3 points for triangles)
 
-element = Element(topology, basis, integration, 
-                  connectivity=(1, 2, 3))
+element = Element(topology, basis, integration, (1, 2, 3))  # Tuple!
 
 # Type-stable construction (preferred)
-element = Element{Tri3, Lagrange{1}, Gauss{3}}(...)
+element = Element{Triangle, Lagrange{Triangle,1}, Gauss{2}, 3}(
+    Triangle(),
+    Lagrange{Triangle, 1}(),
+    Gauss{2}(),
+    (1, 2, 3),  # Tuple for connectivity
+    Dict{Symbol, Any}()
+)
+
+# Same topology, different polynomial order:
+element_p2 = Element(Triangle(), Lagrange{Triangle, 2}(), Gauss{3}(), 
+                      (1, 2, 3, 4, 5, 6))  # 6 nodes for P2
+
+# Same topology, edge DOFs (electromagnetics):
+element_nedelec = Element(Triangle(), Nedelec{Triangle, 1}(), Gauss{2}(),
+                           (1, 2, 3))  # 3 nodes, but DOFs on edges!
+
+# Same topology, different integration:
+element_reduced = Element(Triangle(), Lagrange{Triangle, 1}(), Reduced(),
+                           (1, 2, 3))
+```
 ```
 
 **Benefits:**
@@ -225,30 +301,47 @@ element = Element{Tri3, Lagrange{1}, Gauss{3}}(...)
 ```text
 src/
   topology/
-    tri3.jl         # Reference triangle
-    quad4.jl        # Reference quadrilateral
-    tet10.jl        # Reference tetrahedron
-    hex8.jl         # Reference hexahedron
-    ...
+    point.jl           # 0D: Point
+    segment.jl         # 1D: Segment (line)
+    triangle.jl        # 2D: Triangle (no node count!)
+    quadrilateral.jl   # 2D: Quadrilateral
+    tetrahedron.jl     # 3D: Tetrahedron
+    hexahedron.jl      # 3D: Hexahedron (brick)
+    pyramid.jl         # 3D: Pyramid
+    wedge.jl           # 3D: Wedge/prism
+    
+    # Each file defines pure geometry:
+    # - Parametric domain
+    # - Edges, faces
+    # - Reference coordinates (for standard node placements)
+    # NO node count hardcoded!
   
   basis/
-    lagrange.jl              # Lagrange polynomial bases
-    lagrange_generated.jl    # Pre-generated for compile-time
-    hierarchical.jl          # Hierarchical/p-refinement
-    nurbs.jl                 # NURBS for isogeometric
-    ...
+    lagrange.jl              # Lagrange{T, P} implementation
+    serendipity.jl           # Serendipity{T, P} (reduced nodes)
+    nedelec.jl               # Nedelec{T, P} (edge elements)
+    raviart_thomas.jl        # RaviartThomas{T, P} (face elements)
+    hermite.jl               # Hermite{T, P} (C1 continuous)
+    hierarchical.jl          # Hierarchical{T, P} (p-refinement)
+    nurbs.jl                 # NURBS (isogeometric analysis)
+    
+    # Each basis determines:
+    # - Node count (function of topology + polynomial degree)
+    # - DOF placement (nodes, edges, faces, interior)
+    # - Basis function evaluation
   
   integration/
     gauss.jl        # Gauss-Legendre quadrature
     lobatto.jl      # Gauss-Lobatto quadrature
-    reduced.jl      # Reduced integration
-    ...
+    reduced.jl      # Reduced integration (underintegration)
+    
+    # Maps integration scheme + topology → quadrature points
   
   elements/
     element.jl      # Element type definition
     integrate.jl    # Integration loop
     assemble.jl     # Global assembly
-    ...
+```
 ```
 
 **Rationale:**
@@ -276,9 +369,11 @@ where:
 
 ```julia
 function element_stiffness(element::Element{T, B, I}) where {T, B, I}
-    K = zeros(nnodes(T) * ndofs, nnodes(T) * ndofs)
+    n = nnodes(element.basis)  # Node count from BASIS, not topology
+    ndof = ndofs_per_node(element.basis)
+    K = zeros(n * ndof, n * ndof)
     
-    # Get integration points from integration scheme
+    # Get integration points from integration scheme + topology
     ips = integration_points(element.integration, element.topology)
     
     for ip in ips
@@ -286,7 +381,7 @@ function element_stiffness(element::Element{T, B, I}) where {T, B, I}
         N = evaluate_basis(element.basis, ip.ξ)
         dN = evaluate_basis_derivatives(element.basis, ip.ξ)
         
-        # Jacobian (depends on topology + node coordinates)
+        # Jacobian (depends on topology + actual node coordinates)
         J = jacobian(element.topology, element.connectivity, dN)
         
         # Strain-displacement matrix (depends on basis derivatives)
@@ -302,9 +397,10 @@ end
 
 **Notice:** Each concern is accessed through clean interfaces:
 
-- `integration_points()` → integration scheme
-- `evaluate_basis()` → interpolation scheme
-- `jacobian()` → topology + connectivity
+- `nnodes(basis)` → basis determines node count, NOT topology!
+- `integration_points(integration, topology)` → integration scheme
+- `evaluate_basis(basis, ξ)` → interpolation scheme
+- `jacobian(topology, connectivity, dN)` → geometric mapping
 
 ### Type-Stability for Performance
 
@@ -313,12 +409,13 @@ With concrete types, the compiler can specialize:
 ```julia
 # This becomes a specialized function with no runtime overhead
 function element_stiffness(
-    element::Element{Tri3, Lagrange{1}, Gauss{3}, 3}
+    element::Element{Triangle, Lagrange{Triangle,1}, Gauss{2}, 3}
 )
     # Compiler knows at compile time:
-    # - 3 nodes (Tri3)
-    # - 3 basis functions (Lagrange{1})
-    # - 3 integration points (Gauss{3})
+    # - Triangle topology (2D, 3 edges)
+    # - 3 nodes (from Lagrange{Triangle, 1})
+    # - 3 basis functions (P1)
+    # - 3 integration points (Gauss{2} on triangle)
     # - connectivity is NTuple{3, UInt}
     
     # Generated code has:
@@ -336,57 +433,105 @@ end
 ### Adding a New Topology
 
 ```julia
-# File: src/topology/hex27.jl
-struct Hex27 <: AbstractTopology
-    nnodes::Int = 27
-    dim::Int = 3
+# File: src/topology/prism.jl
+"""
+Prism/Wedge element: triangular cross-section extruded in z-direction.
+
+Parametric domain: Triangle × [-1, 1]
+  - (ξ, η) ∈ Triangle (base)
+  - ζ ∈ [-1, 1] (height)
+
+Note: Does NOT specify node count! That comes from basis.
+"""
+struct Prism <: AbstractTopology end
+
+dim(::Prism) = 3
+
+# Parametric domain edges/faces
+function edges(::Prism)
+    # 9 edges: 3 on bottom, 3 on top, 3 vertical
+    return ((1,2), (2,3), (3,1), (4,5), (5,6), (6,4), (1,4), (2,5), (3,6))
 end
 
-# Reference element coordinates
-reference_coordinates(::Hex27) = [
-    # 8 corner nodes
-    (-1, -1, -1), (1, -1, -1), (1, 1, -1), (-1, 1, -1),
-    (-1, -1,  1), (1, -1,  1), (1, 1,  1), (-1, 1,  1),
-    # 12 mid-edge nodes
-    (0, -1, -1), (1, 0, -1), (0, 1, -1), (-1, 0, -1),
-    # ... (continue for all 27 nodes)
-]
+function faces(::Prism)
+    # 5 faces: 2 triangular (top/bottom), 3 quadrilateral (sides)
+    return ((1,2,3), (4,5,6), (1,2,5,4), (2,3,6,5), (3,1,4,6))
+end
 
-# Topology is a pure mathematical object
-# No need to implement assembly, integration, etc.
+# Standard node placements for common basis functions
+function reference_node_positions(::Prism, ::Type{Lagrange{Prism, 1}})
+    # 6 nodes for linear (P1)
+    return [(-1,0,0), (1,0,0), (0,1,0),   # Bottom triangle
+            (-1,0,1), (1,0,1), (0,1,1)]   # Top triangle
+end
+
+function reference_node_positions(::Prism, ::Type{Lagrange{Prism, 2}})
+    # 18 nodes for quadratic (P2)
+    # 6 corners + 9 edge midpoints + 3 face centers
+    return [...]  # Full list
+end
 ```
 
 **Usage:**
 
 ```julia
-element = Element{Hex27, Lagrange{2}, Gauss{3}}(...)
+element = Element(Prism(), Lagrange{Prism, 1}(), Gauss{2}(), (1,2,3,4,5,6))
 # Automatically works with existing assembly code!
+
+element_p2 = Element(Prism(), Lagrange{Prism, 2}(), Gauss{3}(), 
+                      (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18))
+# Same topology, quadratic basis → 18 nodes
 ```
 
 ### Adding a New Interpolation Scheme
 
 ```julia
 # File: src/basis/hierarchical.jl
-struct Hierarchical{P} <: AbstractBasis end
+"""
+Hierarchical polynomial basis for p-refinement.
+
+Unlike Lagrange (nodal basis), hierarchical basis has:
+- Low-order modes at vertices (vertex bubbles)
+- Higher-order modes as edge/face/volume bubbles
+- Easier adaptivity (can increase P without changing low-order modes)
+
+Node count depends on polynomial degree:
+  P=1: same as Lagrange (vertices only)
+  P=2: vertices + edge modes
+  P=3: vertices + edge modes + face modes + volume modes
+"""
+struct Hierarchical{T<:AbstractTopology, P} <: AbstractBasis end
+
+# Node count = vertices + edges*(P-1) + faces*(P-1)*(P-2)/2 + ...
+nnodes(::Hierarchical{Triangle, 1}) = 3
+nnodes(::Hierarchical{Triangle, 2}) = 3 + 3*1 = 6
+nnodes(::Hierarchical{Triangle, 3}) = 3 + 3*2 + 1 = 10
 
 # Evaluate basis functions
-function evaluate_basis(basis::Hierarchical{P}, ξ::Vec) where P
+function evaluate_basis(basis::Hierarchical{Triangle, P}, ξ::Vec) where P
     # Implement hierarchical polynomial evaluation
-    # Return NTuple{N, Float64}
+    # First 3 are vertex functions (same as Lagrange P1)
+    # Next modes are edge bubbles, then face bubbles
+    return NTuple{nnodes(basis), Float64}(...)  # Zero allocation!
 end
 
 # Evaluate basis derivatives
-function evaluate_basis_derivatives(basis::Hierarchical{P}, ξ::Vec) where P
-    # Implement derivatives
-    # Return NTuple{N, Vec}
+function evaluate_basis_derivatives(basis::Hierarchical{Triangle, P}, ξ::Vec) where P
+    # Return tuple of gradients
+    return NTuple{nnodes(basis), Vec}(...)
 end
 ```
 
 **Usage:**
 
 ```julia
-element = Element{Tri3, Hierarchical{3}, Gauss{4}}(...)
-# Use same Tri3 topology with hierarchical basis!
+# Same triangle topology, hierarchical basis instead of Lagrange
+element = Element(Triangle(), Hierarchical{Triangle, 3}(), Gauss{4}(), 
+                  (1,2,3,4,5,6,7,8,9,10))  # 10 nodes for P3
+
+# Can do p-refinement by just changing basis degree!
+element_p2 = Element(Triangle(), Hierarchical{Triangle, 2}(), Gauss{3}(),
+                      (1,2,3,4,5,6))
 ```
 
 ### Adding a New Integration Rule
@@ -395,16 +540,27 @@ element = Element{Tri3, Hierarchical{3}, Gauss{4}}(...)
 # File: src/integration/lobatto.jl
 struct Lobatto{N} <: AbstractIntegration end
 
-function integration_points(::Lobatto{N}, topology::T) where {N, T <: AbstractTopology}
+function integration_points(scheme::Lobatto{N}, topology::T) where {N, T<:AbstractTopology}
     # Return integration points and weights for Lobatto quadrature
+    # Lobatto includes endpoints (useful for spectral methods)
     # Specific to topology dimension
+    
+    if T === Segment
+        # 1D Lobatto points
+        return lobatto_1d(N)
+    elseif T === Triangle
+        # 2D Lobatto-like scheme for triangles
+        return lobatto_triangle(N)
+    # ... other topologies
+    end
 end
 ```
 
 **Usage:**
 
 ```julia
-element = Element{Quad4, Lagrange{1}, Lobatto{3}}(...)
+element = Element(Quadrilateral(), Lagrange{Quadrilateral, 1}(), Lobatto{3}(),
+                  (1,2,3,4))
 # Use Lobatto instead of Gauss for same element!
 ```
 
@@ -415,32 +571,52 @@ element = Element{Quad4, Lagrange{1}, Lobatto{3}}(...)
 The type system prevents invalid combinations:
 
 ```julia
-# ✅ Valid: Tri3 with 2D basis
-element = Element{Tri3, Lagrange{1}, Gauss{3}}(...)
+# ✅ Valid: Triangle with P1 Lagrange basis
+element = Element(Triangle(), Lagrange{Triangle, 1}(), Gauss{2}(), (1,2,3))
 
-# ❌ Compile error: Cannot use 3D topology with 2D basis (if we enforce)
-element = Element{Hex8, TriangularBasis, Gauss{3}}(...)
+# ✅ Valid: Triangle with P2 Lagrange basis (6 nodes)
+element = Element(Triangle(), Lagrange{Triangle, 2}(), Gauss{3}(), 
+                  (1,2,3,4,5,6))
 
-# ✅ Valid: Mix different integration rules
-element1 = Element{Quad4, Lagrange{1}, Gauss{4}}(...)     # Full integration
-element2 = Element{Quad4, Lagrange{1}, Reduced}(...)       # Reduced integration
-element3 = Element{Quad4, Lagrange{2}, Gauss{9}}(...)     # Quadratic + more points
+# ✅ Valid: Triangle with Nédélec edge elements
+element = Element(Triangle(), Nedelec{Triangle, 1}(), Gauss{2}(), (1,2,3))
+# Note: 3 nodes, but DOFs are on edges!
+
+# ❌ Compile error: Hexahedron basis on Triangle topology (if enforced)
+element = Element(Triangle(), Lagrange{Hexahedron, 1}(), Gauss{2}(), (1,2,3))
+# Type mismatch: basis topology must match element topology
+
+# ✅ Valid: Same topology, different integration rules
+element1 = Element(Quadrilateral(), Lagrange{Quadrilateral,1}(), Gauss{4}(), 
+                    (1,2,3,4))     # Full integration
+element2 = Element(Quadrilateral(), Lagrange{Quadrilateral,1}(), Reduced(), 
+                    (1,2,3,4))     # Reduced integration
+element3 = Element(Quadrilateral(), Lagrange{Quadrilateral,2}(), Gauss{9}(), 
+                    (1,2,3,4,5,6,7,8,9))  # Quadratic + more points
 ```
 
 ### Number of Nodes Known at Compile Time
 
 ```julia
-# Connectivity is NTuple{N, UInt} where N is known at compile time
+# Connectivity is NTuple{N, UInt} where N is determined by basis, not topology!
 struct Element{T, B, I, N}
     topology::T
     basis::B
     integration::I
-    connectivity::NTuple{N, UInt}  # N from topology
+    connectivity::NTuple{N, UInt}  # N = nnodes(B)
 end
 
+# For Lagrange{Triangle, 1}: N = 3
+# For Lagrange{Triangle, 2}: N = 6
+# For Nedelec{Triangle, 1}: N = 3 (still 3 nodes, DOFs on edges)
+
 # Compiler can unroll loops over connectivity
-for i in 1:length(element.connectivity)
-    # Loop is unrolled at compile time!
+function process_element(element::Element{T, B, I, N}) where {T, B, I, N}
+    for i in 1:N  # N known at compile time
+        # Loop is unrolled at compile time!
+        node_id = element.connectivity[i]
+        # ...
+    end
 end
 ```
 
@@ -448,32 +624,44 @@ end
 
 ## Backward Compatibility
 
-### Type Aliases for Old Code
+### Migration from Hardcoded Types
 
 ```julia
-# Old code used to write:
-# element = Element("Tri3", ...)
+# Old Code Aster style (what we're moving away from):
+# element_type = "TRIA3"  # Hardcoded node count
 
-# Provide type aliases:
-const Tri3Element = Element{Tri3, Lagrange{1}, Gauss{3}}
-const Quad4Element = Element{Quad4, Lagrange{1}, Gauss{4}}
+# JuliaFEM modern approach:
+element = Element(Triangle(), Lagrange{Triangle, 1}(), Gauss{2}(), (1,2,3))
 
-# Old code still works:
-element = Tri3Element(connectivity=(1,2,3))
+# For transition, provide type aliases:
+const TRIA3 = Lagrange{Triangle, 1}
+const TRIA6 = Lagrange{Triangle, 2}
+const QUAD4 = Lagrange{Quadrilateral, 1}
+const QUAD8 = Serendipity{Quadrilateral, 2}
+const QUAD9 = Lagrange{Quadrilateral, 2}
+
+# Old code can use aliases:
+element = Element(Triangle(), TRIA3(), Gauss{2}(), (1,2,3))
 ```
 
 ### Constructor Convenience
 
 ```julia
-# Convenience constructors for common cases
-function Element(::Type{Tri3}, connectivity::NTuple{3, UInt})
-    Element{Tri3, Lagrange{1}, Gauss{3}}(
-        Tri3(), Lagrange{1}(), Gauss{3}(), connectivity, Dict()
-    )
+# Convenience constructors for common cases (implicit defaults)
+function Element(topology::Triangle, connectivity::NTuple{3, Int})
+    # Assume: P1 Lagrange + standard Gauss quadrature
+    Element(topology, Lagrange{Triangle, 1}(), Gauss{2}(), 
+            UInt.(connectivity))
+end
+
+function Element(topology::Triangle, connectivity::NTuple{6, Int})
+    # Infer P2 from 6 nodes
+    Element(topology, Lagrange{Triangle, 2}(), Gauss{3}(), 
+            UInt.(connectivity))
 end
 
 # User can still write simple code:
-element = Element(Tri3, (1, 2, 3))
+element = Element(Triangle(), (1, 2, 3))  # Defaults to P1 + Gauss{2}
 ```
 
 ## Performance Implications
@@ -592,28 +780,53 @@ Deal.II has sophisticated separation:
 
 **Element** = Topology + Interpolation + Integration + Fields
 
-This simple equation guides JuliaFEM's architecture:
+But with critical insight:
 
-1. **Topology** defines the reference element (mathematical object)
-2. **Interpolation** defines how to interpolate (approximation theory)
-3. **Integration** defines how to integrate (numerical analysis)
-4. **Fields** define what data lives on the element (problem-specific)
+- **Topology** = geometric shape ONLY (no hardcoded node count!)
+- **Interpolation** = polynomial space + DOF placement (determines node count)
+- **Integration** = numerical quadrature (independent choice)
+- **Fields** = problem-specific data
+
+**The Code Aster/Abaqus anti-pattern we avoid:**
+- ❌ `TRIA3`, `TRIA6`, `TRIA7`, `TRIA10` → hardcoded node counts
+- ✅ `Triangle` + `Lagrange{Triangle, P}` → node count derived from P
 
 **Benefits:**
 
-- ✅ Clear separation of concerns
-- ✅ Mix-and-match flexibility
-- ✅ Type system enforcement
-- ✅ 100× performance improvement
-- ✅ Maintainable, extensible codebase
+- ✅ **One topology, infinite possibilities** (P1, P2, P3, Nédélec, Raviart-Thomas, ...)
+- ✅ **Clear separation** (geometry ≠ approximation ≠ integration)
+- ✅ **Type system enforcement** (compiler catches mismatches)
+- ✅ **Zero-cost abstractions** (100× performance improvement)
+- ✅ **Extensible** (add new basis without touching topology)
+- ✅ **Educational** (code teaches FEM mathematics correctly)
 
-**Trade-off:**
+**Key architectural decision:**
+```julia
+// ❌ WRONG (Code Aster style)
+struct TRIA3 { int nnodes = 3; }  // Hardcoded!
 
-- ⚠️ More complex initial setup
-- ⚠️ Requires understanding of Julia's type system
-- ⚠️ Documentation must be excellent
+// ✅ CORRECT (JuliaFEM style)
+struct Triangle <: AbstractTopology end  // Pure geometry
+nnodes(::Lagrange{Triangle, 1}) = 3     // Derived from basis
+nnodes(::Lagrange{Triangle, 2}) = 6     // Different basis → different count
+nnodes(::Nedelec{Triangle, 1}) = 3      // Edge DOFs, still 3 nodes
+```
 
-**Result:** A modern, high-performance, extensible FEM library that teaches good software engineering alongside finite element methods.
+**Trade-offs:**
+
+- ⚠️ More complex type system (but Julia handles it elegantly)
+- ⚠️ Requires understanding separation of concerns
+- ⚠️ Documentation must be excellent (this document!)
+
+**Result:** A modern, mathematically correct, high-performance, extensible FEM library that can handle:
+- Standard nodal FEM (Lagrange)
+- Edge elements (electromagnetics with Nédélec)
+- Face elements (fluid flow with Raviart-Thomas)
+- Mixed formulations (Taylor-Hood, MINI, ...)
+- Isogeometric analysis (NURBS)
+- hp-refinement (hierarchical basis)
+
+All with **one unified Element type** and **zero runtime overhead**.
 
 ---