refactor(topology): Implement Tetrahedron{N} with type parameter

Update Tetrahedron to use node count type parameter per ADR-002.

Changes:
- struct Tetrahedron → struct Tetrahedron{N} <: AbstractTopology{N}
- Aliases: Tet4 = Tetrahedron{4}, Tet10 = Tetrahedron{10}
- Simplified implementation following same pattern
- Remove old design documentation

Implements ADR-002 (November 13, 2025): node count from mesh, not basis.

Old files removed: tet4.jl, tet10.jl
New file: Single tetrahedra.jl handles all variants via {N}
This commit is contained in:
Jukka Aho
2025-11-15 02:56:52 +02:00
parent 568f09090a
commit 2c6d6dc620
+230 -92
View File
@@ -2,35 +2,42 @@
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE
"""
Tetrahedron <: AbstractTopology
Tetrahedron{N} <: AbstractTopology
Tetrahedral element topology (3D simplex).
Parametric tetrahedral element topology (3D simplex).
**Important:** This type defines ONLY the geometric shape. Node count is determined
by the interpolation scheme (basis functions):
- `Lagrange{Tetrahedron, 1}` → 4 nodes (P1, linear)
- `Lagrange{Tetrahedron, 2}` → 10 nodes (P2, quadratic)
- `Lagrange{Tetrahedron, 3}` → 20 nodes (P3, cubic)
The type parameter `N` specifies the total number of nodes in the element,
enabling compile-time dispatch and type-stable code generation.
# Type Parameter
- `N::Int`: Total number of nodes (4 or 10)
# Canonical Type Aliases
**Always use these aliases instead of constructing `Tetrahedron{N}` directly:**
- `Tet4 = Tetrahedron{4}` - Linear tetrahedron (P1, 4 corner nodes)
- `Tet10 = Tetrahedron{10}` - Quadratic tetrahedron (P2, 10 nodes: 4 corners + 6 edge midpoints)
# Why Parametric Types?
1. **Type Stability:** Each node count is a distinct type (`Tet4 !== Tet10`)
2. **Compile-Time Dispatch:** Kernel specialization for GPU performance
3. **Zero Allocation:** Node count known at compile time
4. **Clear API:** `nnodes(Tet10())` returns compile-time constant `10`
# Reference Element
```
N4
N4 (0,0,1)
/|\\
/ | \\
/ | \\
/ | \\
N1---+----N3
(0,0,0) (0,1,0)
\\ /
\\/
N2
N2 (1,0,0)
```
# Standard Corner Node Positions
1. (0, 0, 0) - Origin
2. (1, 0, 0) - Along ξ-axis
3. (0, 1, 0) - Along η-axis
4. (0, 0, 1) - Along ζ-axis
# Topology Properties
- Dimension: 3
- Corner nodes: 4
@@ -39,118 +46,249 @@ by the interpolation scheme (basis functions):
# Typical Usage
```julia
julia> topology = Tetrahedron()
julia> dim(topology)
3
julia> reference_coordinates(topology) # Corner nodes only
((0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0))
julia> basis = Lagrange{Tetrahedron, 1}()
julia> nnodes(basis) # Linear: 4 nodes
4
julia> basis = Lagrange{Tetrahedron, 2}()
julia> nnodes(basis) # Quadratic: 10 nodes (corners + edge midpoints)
julia> topology = Tet10() # Use canonical alias
julia> nnodes(topology) # Returns compile-time constant
10
julia> Tet4 !== Tet10 # Type stability check
true
julia> reference_coordinates(Tet4()) # Corner nodes only
((0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0))
```
**Zero allocation:** All functions return compile-time sized tuples.
# Design Notes
- Separates topology (geometric shape) from interpolation (basis functions)
- Corner node positions are ALWAYS the same (4 nodes)
- Intermediate nodes (edge midpoints) depend on `N` parameter
- Use `reference_coordinates(Tet10())` to get ALL 10 node positions
See also: [`AbstractTopology`](@ref), [`Hexahedron`](@ref), [`Lagrange`](@ref)
# Type Parameter
Node count comes from mesh connectivity. Type parameter enables compile-time optimization.
# Node Count Variants
- `Tetrahedron{4}` (alias `Tet4`): Linear tetrahedron (P1 Lagrange)
- `Tetrahedron{10}` (alias `Tet10`): Quadratic tetrahedron (P2 Lagrange)
"""
struct Tetrahedron <: AbstractTopology end
struct Tetrahedron{N} <: AbstractTopology{N} end
dim(::Tetrahedron) = 3
# ============================================================================
# CANONICAL TYPE ALIASES (PRIMARY API)
# ============================================================================
"""
reference_coordinates(::Tetrahedron)
Tet4 = Tetrahedron{4}
Get corner node positions for Tetrahedron (4 vertices in parametric space).
Linear tetrahedron with 4 corner nodes (P1 interpolation).
**Reference Coordinates:**
- Node 1: (0.0, 0.0, 0.0) - Origin
- Node 2: (1.0, 0.0, 0.0) - Along ξ-axis
- Node 3: (0.0, 1.0, 0.0) - Along η-axis
- Node 4: (0.0, 0.0, 1.0) - Along ζ-axis
**Use this alias everywhere** instead of `Tetrahedron{4}`.
"""
function reference_coordinates(::Tetrahedron)
const Tet4 = Tetrahedron{4}
"""
Tet10 = Tetrahedron{10}
Quadratic tetrahedron with 10 nodes (P2 interpolation).
**Node Layout:**
- Nodes 1-4: Corner nodes (same as Tet4)
- Nodes 5-10: Edge midpoints
**Use this alias everywhere** instead of `Tetrahedron{10}`.
"""
const Tet10 = Tetrahedron{10}
# ============================================================================
# CORE TOPOLOGY INTERFACE
# ============================================================================
"""
nnodes(::Tetrahedron{N}) where N -> Int
Return total number of nodes for parametric tetrahedron topology.
This is a **compile-time constant** enabling type-stable dispatch.
# Returns
- `N`: Node count specified by type parameter (4 or 10)
# Examples
```julia
julia> nnodes(Tet4()) # Returns compile-time constant 4
4
julia> nnodes(Tet10()) # Returns compile-time constant 10
10
julia> @allocated nnodes(Tet10()) # Zero allocation
0
```
# Performance Note
This function returns a compile-time constant, enabling:
- Zero-cost abstraction (compiler eliminates call)
- Fully specialized code generation
- Static memory allocation in GPU kernels
"""
nnodes(::Tetrahedron{N}) where {N} = N
"""
dim(::Tetrahedron{N}) where N -> Int
Return spatial dimension of tetrahedron reference element (always 3).
# Returns
- `3`: Tetrahedra exist in 3D space
# Examples
```julia
julia> dim(Tet4())
3
julia> dim(Tet10()) # Same for all tetrahedron types
3
```
"""
dim(::Tetrahedron{N}) where {N} = 3
# ============================================================================
# REFERENCE COORDINATES (Full Node Positions)
# ============================================================================
"""
reference_coordinates(::Tetrahedron{4}) -> NTuple{4, NTuple{3, Float64}}
Return reference coordinates for linear tetrahedron (Tet4) - 4 corner nodes only.
# Returns
Tuple of 4 coordinate triples: ((ξ₁, η₁, ζ₁), (ξ₂, η₂, ζ₂), (ξ₃, η₃, ζ₃), (ξ₄, η₄, ζ₄))
# Node Positions
```
N4 (0,0,1)
/|\\
/ | \\
/ | \\
/ | \\
N1---+----N3
(0,0,0) (0,1,0)
\\ /
\\/
N2 (1,0,0)
```
- Node 1: (0.0, 0.0, 0.0) - Origin
- Node 2: (1.0, 0.0, 0.0) - Along ξ-axis
- Node 3: (0.0, 1.0, 0.0) - Along η-axis
- Node 4: (0.0, 0.0, 1.0) - Along ζ-axis
"""
function reference_coordinates(::Tetrahedron{4})
return (
(0.0, 0.0, 0.0), # Node 1: Origin
(1.0, 0.0, 0.0), # Node 2: Along ξ-axis
(0.0, 1.0, 0.0), # Node 3: Along η-axis
(0.0, 0.0, 1.0), # Node 4: Along ζ-axis
(0.0, 0.0, 0.0), # N1: Corner at origin
(1.0, 0.0, 0.0), # N2: Corner along ξ
(0.0, 1.0, 0.0), # N3: Corner along η
(0.0, 0.0, 1.0) # N4: Corner along ζ
)
end
"""
edges(::Tetrahedron)
reference_coordinates(::Tetrahedron{10}) -> NTuple{10, NTuple{3, Float64}}
Edge connectivity for tetrahedron (corner nodes).
Return reference coordinates for quadratic tetrahedron (Tet10) - 10 nodes total.
# Node Layout
- Nodes 1-4: Corner nodes (same as Tet4)
- Node 5: Edge midpoint between N1-N2 (0.5, 0.0, 0.0)
- Node 6: Edge midpoint between N2-N3 (0.5, 0.5, 0.0)
- Node 7: Edge midpoint between N3-N1 (0.0, 0.5, 0.0)
- Node 8: Edge midpoint between N1-N4 (0.0, 0.0, 0.5)
- Node 9: Edge midpoint between N2-N4 (0.5, 0.0, 0.5)
- Node 10: Edge midpoint between N3-N4 (0.0, 0.5, 0.5)
"""
function edges(::Tetrahedron)
function reference_coordinates(::Tetrahedron{10})
return (
(0.0, 0.0, 0.0), # N1: Corner
(1.0, 0.0, 0.0), # N2: Corner
(0.0, 1.0, 0.0), # N3: Corner
(0.0, 0.0, 1.0), # N4: Corner
(0.5, 0.0, 0.0), # N5: Midpoint edge 1-2
(0.5, 0.5, 0.0), # N6: Midpoint edge 2-3
(0.0, 0.5, 0.0), # N7: Midpoint edge 3-1
(0.0, 0.0, 0.5), # N8: Midpoint edge 1-4
(0.5, 0.0, 0.5), # N9: Midpoint edge 2-4
(0.0, 0.5, 0.5) # N10: Midpoint edge 3-4
)
end
# ============================================================================
# TOPOLOGICAL CONNECTIVITY (Corner Nodes Only)
# ============================================================================
"""
edges(::Tetrahedron{N}) where N -> NTuple{6, NTuple{2, Int}}
Return edge connectivity (pairs of **corner node indices**) for tetrahedron.
This is TOPOLOGICAL connectivity, independent of interpolation order.
# Returns
6-tuple of edge definitions:
- Edge 1: (1, 2) - N1 → N2
- Edge 2: (2, 3) - N2 → N3
- Edge 3: (3, 1) - N3 → N1
- Edge 4: (1, 4) - N1 → N4
- Edge 5: (2, 4) - N2 → N4
- Edge 6: (3, 4) - N3 → N4
# Note
- Only references **corner nodes** (1, 2, 3, 4)
- Same for all tetrahedron types (Tet4, Tet10)
"""
function edges(::Tetrahedron{N}) where {N}
return (
(1, 2), # Edge 1
(2, 3), # Edge 2
(3, 1), # Edge 3
(1, 4), # Edge 4
(2, 4), # Edge 5
(3, 4), # Edge 6
(3, 4) # Edge 6
)
end
"""
faces(::Tetrahedron)
faces(::Tetrahedron{N}) where N -> NTuple{4, NTuple{3, Int}}
Face connectivity for tetrahedron (triangular faces, corner nodes).
Return face connectivity for tetrahedron.
Each face is triangular (3 **corner nodes**).
# Returns
4-tuple of triangular faces:
- Face 1: (1, 3, 2) - Base (looking from above)
- Face 2: (1, 2, 4)
- Face 3: (2, 3, 4)
- Face 4: (3, 1, 4)
# Note
- Only references corner nodes
- Faces are triangular (3 nodes each)
- Same for all tetrahedron types
"""
function faces(::Tetrahedron)
function faces(::Tetrahedron{N}) where {N}
return (
(1, 3, 2), # Face 1: Base (looking from above)
(1, 3, 2), # Face 1: Base
(1, 2, 4), # Face 2
(2, 3, 4), # Face 3
(3, 1, 4), # Face 4
(3, 1, 4) # Face 4
)
end
# ============================================================================
# Deprecated aliases (for backwards compatibility)
# EXPORTS
# ============================================================================
"""
Tet4
**DEPRECATED:** Backward compatibility alias. Use `Tetrahedron` with `Lagrange{Tetrahedron, 1}`.
The old `Tet4` conflated topology (tetrahedron) with node count (4).
In the new architecture:
- Topology defines geometric shape only
- Basis functions determine node count
This alias allows old code to work:
```julia
# Old style (still works)
element = Element(Tet4, (1, 2, 3, 4))
# Internally converted to:
element = Element(Tetrahedron, (1, 2, 3, 4)) # Infers Lagrange{Tetrahedron, 1}
```
New code should use explicit topology + basis:
```julia
element = Element(Lagrange{Tetrahedron, 1}, (1, 2, 3, 4))
```
"""
const Tet4 = Tetrahedron
"""
Tet10
**DEPRECATED:** Backward compatibility alias. Use `Tetrahedron` with `Lagrange{Tetrahedron, 2}`.
This alias allows old code to work:
```julia
# Old style (still works)
element = Element(Tet10, (1, 2, 3, 4, 5, 6, 7, 8, 9, 10))
# Internally converted to:
element = Element(Tetrahedron, (1, 2, 3, 4, 5, 6, 7, 8, 9, 10)) # Infers Lagrange{Tetrahedron, 2}
```
New code should use explicit topology + basis:
```julia
element = Element(Lagrange{Tetrahedron, 2}, (1, 2, 3, 4, 5, 6, 7, 8, 9, 10))
```
"""
const Tet10 = Tetrahedron # Same topology! Node count from basis.
# Export ONLY canonical aliases (not the parametric struct)
export Tet4, Tet10