mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
docs(architecture): Complete rewrite explaining topology vs basis separation
Major documentation update (380 additions, 167 deletions):
- Explained topology is pure geometry (NO hardcoded node counts)
- Clarified basis determines BOTH polynomial degree AND node count
- Distinguished node count (connectivity) vs DOF count (unknowns)
- Added examples: Nedelec (edge DOFs), Raviart-Thomas (face DOFs)
- Documented Lagrange{Topology, P} parametric architecture
- Showed why Tri3/Quad4/Tet10 names are anti-pattern
- Updated all code examples to use new architecture
- Explained Serendipity vs full Lagrange tensor products
- Added performance implications and trade-offs
- Showed how one Triangle topology works for P1/P2/P3/Nedelec/etc
This commit is contained in:
Jukka Aho
+377
-164
@@ -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**.
|
||||
|
||||
---
|
||||
|
||||
|
||||
Reference in New Issue
Block a user