Introduce root-level AGENTS.md as the architecture and workflow summary for
humans and coding agents. Add docs/src/repository_layout.md as a placement
and anti-pattern reference for sources, tests, scripts, and gitignored trees.
- AGENTS.md links deep docs and states 2.0 invariants
- repository_layout.md includes .cursor/ as gitignored local tooling
Replace outdated badges and roadmap prose with the current 2.0 focus areas,
installation notes, and pointers to AGENTS.md. Rewrite CONTRIBUTING as a
concise fork/test/layer-contract checklist aligned with the repo layout.
- README: drop stale Travis/Gitter stack; document legacy module caveat
- CONTRIBUTING: remove Quarto front matter; link architecture_layers checker
Runnable 4×4×4 Hex8 linear elasticity assembly used as the Documenter/README
minimal example source of truth; returns (; ndofs, nnz_stiffness).
- Intended for @literalinclude, test/docs, and scripts/verify_docs_quickstart.jl
Reorganize API documentation with organized sections for better
navigation and discoverability.
- Add structured sections: Elements, Physics, Fields, Formulations,
Materials, Mesh, Topology, Basis Functions, DOF System, Solvers
- Include key types and functions in each section
- Change title from "API documentation" to "API Reference"
- Keep index at the end for complete reference
Simplify documentation build script to generate API docs only.
Remove legacy code for generating examples, packages, and guides.
- Replace complex build system with simple Documenter.jl setup
- Generate markdown output for Quarto integration
- Build directly to juliafem.github.io/api/ directory
- Remove functions: copy_docs, add_page!, generate_* functions
- Use DocumenterMarkdown format instead of HTML
Remove documentation files from docs/ directory as they have
been moved to juliafem.github.io/ according to the website
structure. Manual documentation now lives in the website
repository directly.
Removed files:
- docs/book/adr-004-integration-points-api.md
- docs/book/adr-005-element-integration-point-indices.md
- docs/book/deformation_gradient_implementation.md
- docs/book/design/backend_transparent_architecture.md
- docs/book/design/gmres_algorithm_gpu.md
- docs/book/design/matrix_free_dirichlet_bc.md
- docs/book/element_assembly_implementation.md
- docs/book/linear_elastic_implementation.md
- docs/book/migration-guide-basis-api.md
- docs/book/neo_hookean_implementation.md
- docs/book/perfect_plasticity_implementation.md
Remove docs/src/book directory as it has been moved to another
repository
- Remove docs/src/book/ directory and all contents
- 40 files removed including documentation and design docs
Major revision of ADR-002 based on implementation experience.
Decision change:
- OLD (Nov 9): Topology without node count (pure geometry)
- NEW (Nov 13): Topology with node count type parameter
Rationale:
Node count comes from MESH FILES, not from basis choice. When reading
Abaqus .inp, Code Aster .med, or GMSH .msh files, the mesh explicitly
specifies node count in element connectivity:
- Hex element (1,2,...,8) → 8 nodes
- Hex element (1,2,...,20) → 20 nodes
- Hex element (1,2,...,27) → 27 nodes
The mesh reader knows node count BEFORE basis functions are selected.
Therefore, topology must include node count: Hexahedron{N}.
New design:
struct Hexahedron{N} <: AbstractTopology end
const Hex8 = Hexahedron{8}
const Hex20 = Hexahedron{20}
This maintains type stability (N known at compile time) while
acknowledging that N originates from mesh data, not basis choice.
Supersedes: ADR-002 (November 9, 2025)
- Complete cantilever beam example from mesh to visualization
- Gmsh mesh creation with physical groups for BCs
- Material definition (Young's modulus and Poisson's ratio)
- Dirichlet (fixed) and Neumann (pressure) boundary conditions
- ElasticityPhysics problem setup and solve!() call
- Results visualization with stress and displacement
- 361 lines: Step-by-step user tutorial for beginners
- 4 design documents for GPU implementation (2588 lines total)
- gpu_physics_architecture.md: Physics{Elasticity} GPU-first design
- gpu_cpu_backend_architecture.md: Backend selection and dispatch
- gpu_cpu_migration_guide.md: Migration from old API to new
- gpu_elasticity_refactoring.md: Complete refactoring strategy
- Zero CPU-GPU transfer during solve, matrix-free CG
- Elements store geometry, no mesh dependency
- Breaking changes allowed for GPU performance
- Design principle: users never see CPU/GPU differences
- Three-layer architecture: User API / Backend Abstraction / Implementations
- Auto() backend selection based on hardware availability
- Physics{ElasticityPhysicsType} as single problem type
- Internal conversion between CPU arrays and GPU arrays
- solve!() with automatic dispatch to CPU or GPU backend
- 611 lines: Complete architecture design proposal
- get_basis_functions() and get_basis_derivatives() recommended
- Separation of concerns: topology and basis as separate parameters
- Side-by-side examples for all common element types
- Complete assembly example showing migration path
- Type-stable implementation with no performance penalty
- 302 lines: Comprehensive migration documentation
- Alternative to element-by-element assembly for GPU/matrix-free
- Node-by-node loop eliminates atomic operations on GPU
- Spider pattern: nodes couple with 10-30 neighbors not all N
- NodeToElementsMap: inverse connectivity (node → elements)
- get_node_spider() finds coupled nodes for sparse stiffness
- NodalStiffnessContribution: 3×3 blocks per node
- 307 lines: Experimental architecture with working prototype
- Compressible Neo-Hookean strain energy function
- Automatic differentiation for stress and tangent computation
- Dual constructor: Lamé (μ,λ) or engineering (E,ν)
- Total Lagrangian formulation with 2nd Piola-Kirchhoff stress
- Zero-allocation AD via Tensors.jl
- When to use: rubber, large deformation, contact mechanics
- 571 lines: Complete AD-based material model documentation
- Complete mathematical foundation of Hooke's law in tensor form
- Lamé parameters derived from Young's modulus and Poisson's ratio
- compute_stress() implementation achieving ~25 ns execution
- Fourth-order elasticity tensor with symmetries
- Zero-allocation SIMD-optimized implementation
- Physical constraints and thermodynamic admissibility
- 736 lines: Authoritative implementation documentation
- Reference implementation of element-by-element assembly
- ElementAssemblyData and ElementContribution data structures
- Sparse matrix assembly in COO then CSC format
- scatter_to_global!() adds local to global system
- Penalty method for Dirichlet BCs
- Matrix-vector product interface for GMRES
- 479 lines: Complete documentation with examples and tests
- Mathematical derivation of F = I + ∇u for finite strain
- Zero-allocation implementation achieving 34 ns median
- LLVM IR analysis confirms 0 heap allocations
- 92 SIMD vector operations detected
- Small strain vs finite strain formulations
- Comparison with old deprecated eval_dbasis!() API
- 502 lines: Complete performance analysis with benchmarks
Document decision to store integration point indices instead of data in Element struct.
Key rationale: Elements should store relationships (indices), not data, for memory
efficiency and consistency with node connectivity pattern. Aligns with nodal assembly
approach and GPU-friendly architecture.
Architectural Decision Record documenting design of integration points
API for high-performance finite element assembly.
Decision: Compile-time function returning tuple of (weight, Vec{D})
matching eval_basis! zero-cost abstraction pattern.
Problem context:
- OLD API: Runtime dispatch with mutable struct containing Dict
- Performance penalty: ~50× slower due to type instability
- Allocations: New struct created every query
- Impact: Millions of calls during assembly
Solution properties:
- Compile-time generation (fully inlined)
- Vec{D} from Tensors.jl for FEM math
- Zero allocation (tuples, stack-only)
- Type-stable (all types known at compile time)
- GPU compatible (no heap allocations)
API signature:
get_gauss_points!(::Type{Topology}, ::Type{Gauss{order}})
→ NTuple{N, Tuple{Float64, Vec{D}}}
Alternatives rejected:
- Plain tuples (less convenient for FEM math)
- Store in element (overhead, less flexible)
- Global constants (not composable)
- Runtime dispatch (type-unstable, slow)
Status: Accepted, implemented in src/integration/ (193 lines)
- Delete docs/blog/ directory (files moved to docs/src/book/blog/)
- Delete docs/design/ directory (files moved to docs/src/book/design/)
- Cleanup after three-tier documentation reorganization
- Old locations no longer needed after migration to docs/src/ structure
- Relocate docs/contributor/ to docs/src/contributor/
- Add three GPU quickstart guides (renamed from UPPERCASE to snake_case):
- gpu_elasticity_quickstart.md
- gpu_nodal_assembly_quickstart.md
- quick_reference_gpu.md
- Part of three-tier docs reorganization following Documenter.jl standard
- All files now under docs/src/ for automatic rendering
- Relocate docs/blog/immutability_performance.md to docs/src/book/blog/
- Comprehensive guide on immutable material models with Tensors.jl
- Covers LinearElastic, NeoHookean, PerfectPlasticity implementations
- Includes full benchmarks: 5× speedup for linear, 21× for plasticity
- Zero allocation performance validated
- Part of three-tier docs reorganization under standard docs/src/ structure
- Relocate docs/book/README.md to docs/src/book/README.md
- Follows standard Julia documentation structure where all source files live under docs/src/
- File contains YAML header and book philosophy/structure overview
- Part of three-tier documentation reorganization (user/contributor/book)
- Converted header metadata to YAML frontmatter format
- Added categories and tags for documentation site compatibility
- Preserved all existing content (only header format changed)
- Status: IMPLEMENTED, Phase: Phase 1B
- Links to benchmark: element_immutability_benchmark.jl
New 139-line quick-reference article covering:
- Side-by-side code comparisons (mutable vs immutable)
- 130x speedup summary with key metrics
- Type stability explanation with timing breakdown
- Compiler optimization differences
- Real-world impact table (2.4s → 0.02s)
- Mental model shift (1990s C++ → 2025 modern compilers)
- Quick command to run benchmark
- Links to full article for details
Implements compatibility layer to allow old test code to run with new
immutable element design (though fields won't actually update).
src/elements/elements.jl:
- Replaced has_dfield/get_dfield to work with new fields API
- Fixed get_sfield/get_dfield to handle empty Tuple{} fields
- All dfield functions now map to element.fields (immutable NamedTuple)
src/topology/*.jl (seg2, tri3, quad4, tet4, hex8):
- Added nnodes() implementation for each topology type
- Returns corner node count (backwards compatibility)
- Example: nnodes(::Triangle) = 3, nnodes(::Hexahedron) = 8
- Note: Actual node count depends on basis degree in new architecture
Test Results:
- test_topology_standalone.jl: 36/36 tests passing ✓
- Full test suite: 43 errors (same as before)
- Error breakdown:
* 40+ tests: Problem types not defined (Elasticity, Heat, Mortar)
* 2 tests: Mesh readers not defined (aster_read_mesh)
* 1 test: Tries to mutate empty element (test_elasticity_1d)
Next Steps:
- Tests that create empty elements then mutate need rewriting
- Pattern: Element(Seg2, (1,2)) + update!() → not compatible
- New pattern: Element(..., fields=(geometry=X, displacement=u))
- See docs/design/IMMUTABILITY.md for migration guide
Created comprehensive documentation and benchmark demonstrating why immutable
elements with type-stable fields are 40-130x faster than mutable Dict-based
elements.
benchmarks/element_immutability_benchmark.jl:
- Compares mutable (Dict) vs immutable (NamedTuple) implementations
- Measures field access, updates, assembly loops, large-scale meshes
- Results: 40x faster field access, 130x faster assembly, zero allocations
docs/design/IMMUTABILITY.md:
- Explains counterintuitive API change: element = update(element, ...)
- Benchmarks show 40-130x speedup despite 'copying' elements
- Key insight: Type stability >> mutation, compiler optimizes away copies
- Migration guide: old mutable API → new immutable API
- GPU/HPC rationale: Only bits types work on GPU (no pointers)
Key Results:
- Field access: 1ns vs 45ns (40x faster)
- Assembly: 9ns vs 1124ns per element (130x faster)
- Large mesh: 0.01ms vs 1.2ms for 1000 elements (120x faster)
- Memory: 0 allocations vs 70,000 allocations
- GPU: Compatible (bits types) vs Incompatible (pointers)
This documents a fundamental architectural decision for JuliaFEM 1.0.
Design for handling both nodal and element fields in nodal assembly:
Architecture:
- Nodes have geometry (immutable)
- Elements have connectivity + fields (immutable struct)
- Nodal fields: displacement, temperature, contact pressure
- Element fields: integration point data (σ, ε_plastic, α, C)
Update pattern:
- Create new field containers (NamedTuples)
- Create new elements with updated fields
- Shallow copy element vector, replace elements
- All immutable (GPU-compatible, thread-safe)
GPU kernel:
- Loops over nodes (nodal assembly)
- Accesses nodal_fields for global quantities
- Accesses element.fields for integration point data
- Gathers from connected elements (node_to_elements)
- No atomic operations (each node owns DOFs)
Material state update:
- Process elements in parallel (Threads.@threads)
- Extract nodal displacements from solution
- Compute strains at integration points
- Run material model (plasticity, damage, etc.)
- Create new elements with updated state
- Return new problem with updated fields
Newton iteration:
- Residual uses element.fields.C (current tangent)
- GMRES with matrix-free matvec (nodal assembly)
- Material update after each iteration
- All data structures immutable throughout
Benchmarks show creating new containers ~1000× faster than deepcopy
Explains why JuliaFEM uses nodal assembly instead of element assembly:
Five major advantages:
1. No atomic operations on GPU (each node writes to own DOFs)
2. Contact mechanics is natural (forces at nodes, not elements)
3. Clean domain decomposition (explicit node ownership for MPI)
4. Better cache locality (sequential node processing)
5. Adaptive refinement easier (local node operations)
Key data structure:
- NodeSet contains nodes + elements + node_to_elements connectivity
- Inverse connectivity enables gathering from connected elements
- Fields accessed via node_set.fields (type-stable)
Algorithm:
- Loop over nodes (not elements)
- Each node gathers contributions from connected elements
- Direct write to owned DOFs (no race conditions)
- Perfect for matrix-free Krylov methods
Validated with demo:
- CPU/GPU results match exactly (0.0 relative error)
- Average 3.24 elements per node (efficient gathering)
- Natural integration with contact mechanics
Compares to traditional element assembly:
- Element: scatter to nodes (atomic ops, cache misses)
- Nodal: gather from elements (no atomics, better cache)
Complete algorithm for GPU-resident FEM solver with nodal assembly:
- Data partitioning by node ownership (domain decomposition)
- GPU-resident data structures (nodes, elements, connectivity, state)
- Three GPU kernels: residual, matvec, state update
- MPI communication patterns for interface nodes
- Full Newton-GMRES loop on GPU (data stays resident)
Architecture:
- Each GPU owns subset of nodes (exclusive ownership)
- Ghost elements copied for gathering during assembly
- node_to_elements connectivity enables nodal assembly
- No atomic operations (each GPU writes to owned DOFs only)
Key features:
- Data moves to GPU once at start, back once at end
- GMRES iterations entirely on GPU (Arnoldi steps)
- Material state updates on GPU (integration points)
- MPI exchanges only for interface DOFs between iterations
- O(N) memory per GPU (matrix-free)
Handles nonlinearity:
- Element state contains σ, ε_plastic, α, C (tangent)
- Residual kernel uses current stress/tangent
- State update kernel after convergence
- Natural for contact mechanics (nodal forces)
Status: Design document for future GPU implementation
Explains the correct pattern for matrix-vector products in Krylov methods:
- Fields accessed through element_set (not passed separately)
- GPU kernel computes y=K*x (not K itself)
- O(N) memory (vs O(N²) for stored matrix)
- Type-stable field access (compile-time types)
Key insights:
- GMRES needs matvec operation, not the matrix
- ElementSet contains elements + fields together
- Zero allocations with immutable connectivity/fields
- Natural pattern for contact mechanics (nodal updates)
- Material state separate from field parameters
Compares old vs new approach:
- Old: Dict{String,Any} in element (type-unstable)
- New: NamedTuple in ElementSet (type-stable)
- Old: O(N²) matrix storage
- New: O(N) matrix-free operator
Validated with gpu_elementset_matvec_demo.jl:
- GPU/CPU results match exactly
- Fields accessed naturally through element_set
- Returns y vector (what Krylov methods need)
Analyzes field storage patterns and recommends ElementSet approach:
- Element has NO field type parameter (simpler type)
- ElementSet groups elements + shared fields
- Fields can be NamedTuple, struct, any type-stable container
- Embraces immutability (GPU-compatible, thread-safe)
- Separates mutable state from immutable parameters
Design rationale:
- Benchmarks show NamedTuple gives 9-92× speedup vs Dict
- Immutability enables GPU execution without copying
- Creating new containers ~1000× faster than deepcopy
- Matches physical reality (material properties per set)
Compares three options:
1. Fields as type parameter (type proliferation)
2. ElementSet pattern (RECOMMENDED)
3. Hybrid approach (too complex)
Addresses common concerns:
- Time-dependent fields (use interpolation)
- Material state (separate mutable arrays)
- Custom field types (any type-stable container works)
Status: Ready for implementation
Remove duplicate "JuliaFEM Documentation" header (line 11) and add blank lines
for consistent list formatting in three-manual organization document.
Changes:
- Line 11: Removed duplicate H1 title (already in frontmatter)
- Lines 22, 41, 61, 110: Added blank lines before list items for markdown clarity
- Maintains three-audience structure: Users, Contributors, Researchers
- Preserves content organization and cross-reference section
This is a formatting-only change to improve readability. No content modified.
Create 299-line focused design rationale explaining why type stability is essential
for CPU/GPU/MPI performance, without mandating specific storage patterns.
Executive summary (lines 16-26):
- v0.5.1 Dict{String,Any}: 9-92× performance degradation
- Type-unstable code CANNOT run on GPUs
- Significant MPI communication overhead
- Document does NOT prescribe storage location
- Demonstrates why type stability at access points is essential
- Key: Storage pattern matters less than type inference
Problem analysis (lines 28-67):
- Type instability definition: Runtime dispatch when type unknown at compile time
- Why it matters: 10-100× slower CPU, GPU compilation fails, MPI serialization
- Measured impact table: 9-92× speedup, 0 allocations with type stability
- Critical: Zero allocations required for GPU kernels
Design requirements (lines 69-136):
1. Type stability at access points (compiler must infer types)
- Fields could be element-local, global arrays, or arguments
- Access pattern must be type-stable regardless
2. Zero allocations in hot paths (GPU/MPI requirement)
- Assembly loop must allocate nothing
3. Contiguous memory layout (GPU/MPI optimization)
- CUDA transfers contiguous arrays directly
4. Immutable where possible (safe parallelism)
- Thread-safe reads without locks
Demonstrated solutions (lines 138-201) - EXAMPLES, not mandates:
1. NamedTuple container: Simple, type-stable, immutable
2. Struct with typed fields: Explicit, self-documenting
3. Passed as arguments: Maximum type stability, explicit dependencies
- All three achieve type stability
- Choice depends on use case, not performance
GPU and MPI rationale (lines 203-237):
- GPU execution: CUDA requires all code type-stable
- Mock demonstration in benchmarks/gpu_mpi_mock.jl
- MPI communication: Typed arrays use fast memcpy vs slow serialization
- Type stability enables identical code for CPU/GPU
Recommendations (lines 239-256):
- Use type-stable access patterns (REQUIRED)
- Prefer immutable data structures (threading/GPU)
- Pre-allocate caches (zero allocations)
- Use contiguous arrays (GPU/MPI transfer)
- Profile with @btime (verify zero allocations)
- Does NOT mandate: Storage location, container type, dynamic vs static
Validation (lines 258-275):
- benchmarks/field_storage_comparison.jl: 9-92× CPU speedup
- benchmarks/gpu_mpi_mock.jl: GPU/MPI patterns
- benchmarks/VALIDATION_RESULTS.md: Summary table
Conclusion (lines 277-299):
- Type stability is fundamental requirement, not implementation detail
- Enables: High CPU performance, GPU execution, efficient MPI, safe threading
- v1.0 must ensure type stability at access points
- Storage pattern is secondary concern (memory, cache, API)
- Next steps: Review, benchmark, choose pattern, implement, validate CUDA
Key difference from v1: Shorter (299 vs 1114 lines), focused on WHY not HOW,
explicitly states storage pattern is flexible, emphasizes GPU/MPI requirements.
Platform: Julia 1.12.1, November 9, 2025
Series: The JuliaFEM Book, Chapter 5
Status: Design rationale with validated measurements
Create 415-line blog post combining technical demonstration with philosophical
vision for JuliaFEM v1.0 nodal assembly architecture.
Content structure:
- Lines 1-32: Why Krylov+nodal is brilliant for contact mechanics
* Contact is inherently nodal (constraints at nodes, not elements)
* Krylov only needs matvec (never forms global matrix)
* Nodal assembly provides natural row-by-row interface
* O(N) memory vs O(N²) for traditional element assembly
- Lines 34-73: Controversial hypothesis about nodal material modeling
* Claims material state should be at nodes, not integration points
* Argues integration points constrain physics to numerical method
* Variational consistency, physical meaning, scalability arguments
* "I will show them they're wrong" - experimental vision
- Lines 75-95: GMRES advantage for unsymmetric systems
* Material nonlinearity, contact, large deformation all unsymmetric
* GMRES solves positive definite unsymmetric systems
* O(N·iter) time, O(N) memory vs O(N³)/O(N²) for direct solvers
- Lines 97-415: Working GMRES demonstration on 10×10 unsymmetric system
* Problem setup: Positive definite but unsymmetric matrix (lines 105-145)
* Nodal assembly pattern: get_row() interface (lines 147-193)
* Simplified GMRES implementation (lines 195-282)
* Execution and verification (lines 284-318)
* Results: Converged in 10 iterations, 6.28×10⁻¹⁶ relative error
* Key insights section explaining significance (lines 320-365)
* Development roadmap: Immediate → Near-term → Long-term → Vision (lines 367-401)
* Conclusion: Philosophical statement about nodal correctness (lines 403-415)
Technical validation:
- Matrix: 10×10, eigenvalues [9.94, 38.06], condition number 3.83
- Nodal matvec: 1.59×10⁻¹⁴ error vs direct computation
- GMRES: 10 iterations to convergence
- Solution accuracy: 1.23×10⁻¹⁴ absolute error, 6.28×10⁻¹⁶ relative
Dependencies: LinearAlgebra, Random, Printf
Format: Literate.jl (# # for section headers, # for narrative)
Target: Blog post for JuliaFEM v1.0 development documentation
Tone: Opinionated, controversial, technically rigorous
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
Jukka Aho