mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-19 09:54:55 +00:00
chore: removed files from docs
Removed files from docs directory: - docs/badges.jl - docs/book/nodal_assembly_concept.md - docs/build_api.jl - docs/build_doctests.jl - docs/build_lint.jl - docs/build_notebooks.jl - docs/build_unittests.jl - docs/src/contributor/coding_standards.md - docs/src/contributor/guides/gpu_elasticity_quickstart.md - docs/src/contributor/guides/gpu_nodal_assembly_quickstart.md - docs/src/contributor/guides/quick_reference_gpu.md - docs/src/contributor/README.md - docs/src/contributor/status.md - docs/src/contributor/test_fixes_needed.md - docs/src/contributor/testing_philosophy.md - docs/src/examples.md - docs/src/features.md - docs/src/links.md - docs/src/user/elasticity_quickstart.md - docs/src/user/README.md - docs/src/user/system_architecture.md
This commit is contained in:
Jukka Aho
@@ -1,53 +0,0 @@
|
||||
using Requests
|
||||
|
||||
|
||||
function make_badge_simple(subject::String, status::String, color::String, outfile::String)
|
||||
subject = replace(subject, " ", "%20")
|
||||
url = "https://img.shields.io/badge/$subject-$status-$color.svg"
|
||||
res = get(url)
|
||||
open(outfile, "w") do fid
|
||||
write(fid, res.data)
|
||||
end
|
||||
end
|
||||
|
||||
"""
|
||||
Make badge based on two numbers
|
||||
"""
|
||||
function make_badge(subject::String, a::Int, b::Int, outfile::String)
|
||||
limits = Float64[100, 80, 60, 40, 20, 0]
|
||||
colors = ["brightgreen", "green", "yellowgreen", "yellow", "orange", "red"]
|
||||
picked_color = "lightgray"
|
||||
perce = a/b*100
|
||||
println("perce: $perce")
|
||||
for (i, limit) in enumerate(limits)
|
||||
if perce >= limit
|
||||
picked_color = colors[i]
|
||||
break
|
||||
end
|
||||
end
|
||||
println("color: $picked_color")
|
||||
message = "$a/$b"
|
||||
make_badge_simple(subject, message, picked_color, outfile)
|
||||
end
|
||||
|
||||
"""
|
||||
Make badge based on one number
|
||||
"""
|
||||
function make_badge(subject::String, a::Int, outfile::ASCIIString)
|
||||
limits = Float64[0, 20, 40, 60, 80, 100]
|
||||
colors = ["brightgreen", "green", "yellowgreen", "yellow", "orange", "red"]
|
||||
picked_color = "red"
|
||||
for (i, limit) in enumerate(limits)
|
||||
if a <= limit
|
||||
picked_color = colors[i]
|
||||
break
|
||||
end
|
||||
end
|
||||
println("color: $picked_color")
|
||||
a = 100 - a
|
||||
if a < 0
|
||||
a = 0
|
||||
end
|
||||
message = "$a"
|
||||
make_badge_simple(subject, message, picked_color, outfile)
|
||||
end
|
||||
@@ -1,307 +0,0 @@
|
||||
---
|
||||
title: "Nodal Assembly: Concept and Data Structures"
|
||||
date: 2025-11-11
|
||||
author: "JuliaFEM Team"
|
||||
status: "Experimental"
|
||||
last_updated: 2025-11-11
|
||||
tags: ["assembly", "nodal", "gpu", "architecture"]
|
||||
---
|
||||
|
||||
## Introduction
|
||||
|
||||
This document describes the **nodal assembly** concept - an alternative to traditional element-based assembly that is naturally suited for:
|
||||
|
||||
- GPU parallelization (no atomic operations needed)
|
||||
- Matrix-free methods (Krylov solvers)
|
||||
- Contact mechanics (contact is inherently nodal)
|
||||
- Domain decomposition (nodes have clear ownership)
|
||||
|
||||
**Status:** Experimental concept with working prototype. See `src/nodal_assembly_structures.jl` and tests.
|
||||
|
||||
## The Problem with Element Assembly
|
||||
|
||||
Traditional FEM assembles **element by element**:
|
||||
|
||||
```julia
|
||||
# Traditional element assembly
|
||||
for element in elements
|
||||
K_local = compute_element_stiffness(element) # 30×30 for Tet10
|
||||
|
||||
# Scatter to global (requires atomic operations on GPU!)
|
||||
for i in 1:ndofs_local, j in 1:ndofs_local
|
||||
K_global[gdof[i], gdof[j]] += K_local[i,j] # Race condition!
|
||||
end
|
||||
end
|
||||
```
|
||||
|
||||
**Problems:**
|
||||
|
||||
1. **GPU:** Multiple elements write to same global DOF → need atomics → slow
|
||||
2. **Contact:** Contact forces are nodal, but assembly is elemental → mismatch
|
||||
3. **Matrix-free:** Hard to compute K*v without forming K
|
||||
|
||||
## Nodal Assembly Solution
|
||||
|
||||
Assemble **node by node** instead:
|
||||
|
||||
```julia
|
||||
# Nodal assembly
|
||||
for node_i in nodes
|
||||
# Compute contributions FROM all elements touching node_i
|
||||
K_blocks, f_int = compute_nodal_contribution(node_i, elements_touching_i)
|
||||
|
||||
# Each thread owns its node → no atomics needed!
|
||||
w[3*(node_i-1)+1:3*node_i] = matvec_nodal(K_blocks, u)
|
||||
end
|
||||
```
|
||||
|
||||
**Advantages:**
|
||||
|
||||
1. **GPU:** One thread per node, no conflicts, no atomics
|
||||
2. **Contact:** Natural fit (contact forces already nodal)
|
||||
3. **Matrix-free:** Direct K*v computation without forming global K
|
||||
|
||||
## The "Spider" Pattern
|
||||
|
||||
For node $i$, we only compute stiffness blocks for nodes it couples with:
|
||||
|
||||
```text
|
||||
j₃
|
||||
/\
|
||||
/ \
|
||||
/ \
|
||||
j₂------i------j₄ ← Node i's "spider"
|
||||
\ /
|
||||
\ /
|
||||
\/
|
||||
j₁
|
||||
```
|
||||
|
||||
**Key insight:** Most nodes couple with only ~10-30 neighbors (not all N nodes!)
|
||||
|
||||
- **Corner node:** 8 neighbors (1 element touches it)
|
||||
- **Interior node:** 27 neighbors (8 elements touch it)
|
||||
- **Face node:** 12 neighbors (intermediate)
|
||||
|
||||
**Efficiency:** Sparse connectivity preserved without storing full matrix!
|
||||
|
||||
## Data Structures
|
||||
|
||||
### 1. Inverse Mapping: Node → Elements
|
||||
|
||||
```julia
|
||||
struct ElementNodeInfo
|
||||
element_id::Int # Which element
|
||||
local_node_idx::Int # Which local node index (1-10 for Tet10)
|
||||
end
|
||||
|
||||
struct NodeToElementsMap
|
||||
node_to_elements::Vector{Vector{ElementNodeInfo}}
|
||||
nnodes::Int
|
||||
nelements::Int
|
||||
end
|
||||
|
||||
# Usage
|
||||
map = NodeToElementsMap(connectivity)
|
||||
for elem_info in map.node_to_elements[node_i]
|
||||
println("Node $node_i is local node $(elem_info.local_node_idx) ",
|
||||
"in element $(elem_info.element_id)")
|
||||
end
|
||||
```
|
||||
|
||||
**Purpose:** Given node, find all elements touching it (needed for nodal loop).
|
||||
|
||||
### 2. Spider Nodes
|
||||
|
||||
```julia
|
||||
function get_node_spider(map::NodeToElementsMap, node_id::Int,
|
||||
connectivity) -> Vector{Int}
|
||||
spider = Set{Int}()
|
||||
|
||||
# Union of all nodes in elements touching node_id
|
||||
for elem_info in map.node_to_elements[node_id]
|
||||
for node in connectivity[elem_info.element_id]
|
||||
push!(spider, node)
|
||||
end
|
||||
end
|
||||
|
||||
return sort(collect(spider))
|
||||
end
|
||||
```
|
||||
|
||||
**Purpose:** Find all nodes that couple with `node_id` (non-zero stiffness blocks).
|
||||
|
||||
### 3. Nodal Stiffness Contribution
|
||||
|
||||
```julia
|
||||
struct NodalStiffnessContribution{T}
|
||||
node_id::Int
|
||||
spider_nodes::Vector{Int} # Nodes that couple
|
||||
K_blocks::Vector{Tensor{2,3,T}} # 3×3 blocks (one per spider node)
|
||||
f_int::Vec{3,T} # Internal force at this node
|
||||
f_ext::Vec{3,T} # External force at this node
|
||||
end
|
||||
```
|
||||
|
||||
**Purpose:** Storage for nodal assembly. `K_blocks[k]` is the 3×3 coupling between `node_id` and `spider_nodes[k]`.
|
||||
|
||||
**Zero-allocation:** All quantities use `Tensors.jl` types (immutable, stack-allocated).
|
||||
|
||||
## Matrix-Free Matvec
|
||||
|
||||
Given nodal contributions, compute $\mathbf{w} = \mathbf{K} \mathbf{u}$ without forming $\mathbf{K}$:
|
||||
|
||||
```julia
|
||||
function matrix_vector_product_nodal(contrib::NodalStiffnessContribution,
|
||||
u::Vector{Vec{3}}) -> Vec{3}
|
||||
w = zero(Vec{3})
|
||||
|
||||
# Loop over spider nodes (only non-zero columns!)
|
||||
for (k, node_j) in enumerate(contrib.spider_nodes)
|
||||
K_ij = contrib.K_blocks[k] # 3×3 block
|
||||
u_j = u[node_j] # Displacement at node j
|
||||
|
||||
w += K_ij ⋅ u_j # Block matvec
|
||||
end
|
||||
|
||||
return w
|
||||
end
|
||||
```
|
||||
|
||||
**Performance:**
|
||||
|
||||
- Only computes non-zero contributions (sparse spider)
|
||||
- Zero allocations (Tensors.jl)
|
||||
- GPU-friendly (parallel over nodes)
|
||||
|
||||
## Example: 2 Tet4 Elements
|
||||
|
||||
```text
|
||||
Mesh:
|
||||
Element 1: nodes (1,2,3,4)
|
||||
Element 2: nodes (2,3,4,5)
|
||||
|
||||
Nodes 2,3,4 shared between elements
|
||||
```
|
||||
|
||||
**Node 1 (corner):**
|
||||
|
||||
- Touches: 1 element
|
||||
- Spider: [1, 2, 3, 4] (4 nodes)
|
||||
- Needs: 4 × 3×3 blocks
|
||||
|
||||
**Node 2 (interior):**
|
||||
|
||||
- Touches: 2 elements
|
||||
- Spider: [1, 2, 3, 4, 5] (5 nodes = union of both elements)
|
||||
- Needs: 5 × 3×3 blocks
|
||||
|
||||
**Node 5 (corner):**
|
||||
|
||||
- Touches: 1 element
|
||||
- Spider: [2, 3, 4, 5] (4 nodes)
|
||||
- Needs: 4 × 3×3 blocks
|
||||
|
||||
## Assembly Algorithm
|
||||
|
||||
```julia
|
||||
# 1. Build inverse mapping (once, at mesh creation)
|
||||
map = NodeToElementsMap(connectivity)
|
||||
|
||||
# 2. For each node (parallel on GPU)
|
||||
for node_i in 1:nnodes
|
||||
# Find spider
|
||||
spider = get_node_spider(map, node_i, connectivity)
|
||||
|
||||
# Allocate storage
|
||||
contrib = NodalStiffnessContribution(node_i, spider)
|
||||
|
||||
# Loop over elements touching this node
|
||||
for elem_info in map.node_to_elements[node_i]
|
||||
elem = elements[elem_info.element_id]
|
||||
local_idx = elem_info.local_node_idx
|
||||
|
||||
# Compute element contribution to node_i
|
||||
# (loop over integration points inside)
|
||||
compute_element_contribution!(contrib, elem, local_idx, u, time)
|
||||
end
|
||||
|
||||
# Matrix-free matvec: w_i = K_i * u
|
||||
w[node_i] = matrix_vector_product_nodal(contrib, u)
|
||||
end
|
||||
```
|
||||
|
||||
## Comparison to Element Assembly
|
||||
|
||||
| Aspect | Element Assembly | Nodal Assembly |
|
||||
|--------|------------------|----------------|
|
||||
| **Outer loop** | Elements | Nodes |
|
||||
| **Parallelization** | Element → atomics | Node → no atomics |
|
||||
| **Storage** | Full K matrix (sparse) | 3×3 blocks per spider |
|
||||
| **Matrix-free** | Difficult | Natural |
|
||||
| **Contact** | Mismatch | Natural fit |
|
||||
| **GPU** | Slow (atomics) | Fast (no atomics) |
|
||||
|
||||
## Connection to Golden Standard
|
||||
|
||||
This implements the architecture from `docs/src/book/multigpu_nodal_assembly.md`:
|
||||
|
||||
1. ✅ **Nodal assembly** (not element assembly)
|
||||
2. ✅ **3×3 blocks** using `Tensor{2,3}` from Tensors.jl
|
||||
3. ✅ **Matrix-free** matvec with spider pattern
|
||||
4. ✅ **Zero allocations** (immutable Tensor types)
|
||||
|
||||
**Next steps:**
|
||||
|
||||
- Implement `compute_element_contribution!()` for real elements
|
||||
- Integration with material models (already done: `compute_stress()` returns `SymmetricTensor{2,3}`)
|
||||
- GPU kernels for nodal loop
|
||||
- Contact mechanics integration
|
||||
|
||||
## Performance Implications
|
||||
|
||||
**2×2×2 Hex8 mesh (27 nodes, 81 DOFs):**
|
||||
|
||||
- **Element assembly:** 8 elements, each writes to overlapping DOFs → atomics
|
||||
- **Nodal assembly:** 27 nodes, independent writes → no atomics
|
||||
|
||||
**Spider statistics:**
|
||||
|
||||
- Corner node: 8 couplings → compute 8 × 3×3 = 72 entries
|
||||
- Interior node: 27 couplings → compute 27 × 3×3 = 243 entries (all nodes!)
|
||||
- Average node: ~12 couplings → compute 12 × 3×3 = 108 entries
|
||||
|
||||
**Memory:** No global K matrix, only local K_blocks per thread (reused).
|
||||
|
||||
## Testing
|
||||
|
||||
See `test/test_nodal_assembly_structures.jl` for working examples:
|
||||
|
||||
```bash
|
||||
cd /home/juajukka/dev/JuliaFEM.jl
|
||||
julia --project=. test/test_nodal_assembly_structures.jl
|
||||
```
|
||||
|
||||
**Tests:**
|
||||
|
||||
- ✅ Inverse mapping construction
|
||||
- ✅ Spider computation
|
||||
- ✅ Nodal contribution storage
|
||||
- ✅ Matrix-free matvec
|
||||
- ✅ Efficiency analysis (hex mesh)
|
||||
|
||||
## References
|
||||
|
||||
1. **Golden standard:** `docs/src/book/multigpu_nodal_assembly.md`
|
||||
2. **ARCHITECTURE.md:** Nodal assembly motivation
|
||||
3. **TECHNICAL_VISION.md:** Why matrix-free iterative solvers
|
||||
|
||||
## Status
|
||||
|
||||
- **Implementation:** Prototype complete ✅
|
||||
- **Testing:** Basic tests passing ✅
|
||||
- **Integration:** Not yet integrated with main JuliaFEM
|
||||
- **Performance:** Not yet benchmarked
|
||||
- **GPU:** Not yet implemented (but designed for it)
|
||||
|
||||
This is the foundation for the modern JuliaFEM architecture!
|
||||
@@ -1,71 +0,0 @@
|
||||
using Docile, Lexicon, JuliaFEM
|
||||
|
||||
const api_directory = "api"
|
||||
|
||||
"""
|
||||
Searches recursively all the modules from packages. As documentation grows, it's a bit
|
||||
troublesome to add all the new modules manually, so this function searches all the modules
|
||||
automatically.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
module_: Module
|
||||
Module where we want to search modules inside
|
||||
append_list: Array{Module, 1}
|
||||
Array, where we append Modules as we find them
|
||||
|
||||
Returns
|
||||
-------
|
||||
None. Void function, which manipulates the append_list
|
||||
"""
|
||||
function search_modules!(module_::Module, append_list::Array{Module, 1})
|
||||
all_names = names(module_, true)
|
||||
for each in all_names
|
||||
inner_module = module_.(each)
|
||||
if (typeof(inner_module) == Module) && !(inner_module in append_list)
|
||||
push!(append_list, inner_module)
|
||||
search_modules!(inner_module, append_list)
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
append_list = Array(Module, 0)
|
||||
search_modules!(JuliaFEM, append_list)
|
||||
|
||||
const modules = append_list
|
||||
|
||||
# main_folder = dirname(dirname(@__FILE__))
|
||||
# this_folder = dirname(@__FILE__)
|
||||
|
||||
# file_ = "README.md"
|
||||
# run(`cp $main_folder/$file_ $this_folder`)
|
||||
|
||||
cd(dirname(@__FILE__)) do
|
||||
# Run the doctests *before* we start to generate *any* documentation.
|
||||
# for m in modules
|
||||
# failures = failed(doctest(m))
|
||||
# if !isempty(failures.results)
|
||||
# println("\nDoctests failed, aborting commit.\n")
|
||||
# display(failures)
|
||||
# exit(1) # Bail when doctests fail.
|
||||
# end
|
||||
# end
|
||||
# Generate and save the contents of docstrings as markdown files.
|
||||
index = Index()
|
||||
for mod in modules
|
||||
Lexicon.update!(index, save(joinpath(api_directory, "$(mod).rst"), mod))
|
||||
end
|
||||
save(joinpath(api_directory, "index.rst"), index)
|
||||
|
||||
# Add a reminder not to edit the generated files.
|
||||
# open(joinpath(api_directory, "README.md"), "w") do f
|
||||
# print(f, """
|
||||
# Files in this directory are generated using the `build.jl` script. Make
|
||||
# all changes to the originating docstrings/files rather than these ones.
|
||||
# """)
|
||||
# end
|
||||
# save(joinpath(api_directory, "index.rst"), index; md_subheader = :category)
|
||||
# info("Adding all documentation changes in $(api_directory) to this commit.")
|
||||
# success(`git add $(api_directory)`) || exit(1)
|
||||
|
||||
end
|
||||
@@ -1,52 +0,0 @@
|
||||
using JuliaFEM
|
||||
using Lexicon
|
||||
|
||||
include("badges.jl")
|
||||
|
||||
"""
|
||||
Searches recursively all the modules from packages. As documentation grows, it's a bit
|
||||
troublesome to add all the new modules manually, so this function searches all the modules
|
||||
automatically.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
module_: Module
|
||||
Module where we want to search modules inside
|
||||
append_list: Array{Module, 1}
|
||||
Array, where we append Modules as we find them
|
||||
|
||||
Returns
|
||||
-------
|
||||
None. Void function, which manipulates the append_list
|
||||
"""
|
||||
function search_modules!(module_::Module, append_list::Array{Module, 1})
|
||||
all_names = names(module_, true)
|
||||
for each in all_names
|
||||
inner_module = module_.(each)
|
||||
if (typeof(inner_module) == Module) && !(inner_module in append_list)
|
||||
push!(append_list, inner_module)
|
||||
search_modules!(inner_module, append_list)
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
append_list = Array(Module, 0)
|
||||
search_modules!(JuliaFEM, append_list)
|
||||
|
||||
const modules = append_list
|
||||
|
||||
cd(dirname(@__FILE__)) do
|
||||
npassed = 0
|
||||
nfailed = 0
|
||||
nskipped = 0
|
||||
|
||||
for m in modules
|
||||
s = doctest(m)
|
||||
lnpassed, lnfailed, lnskipped = map(length, (passed(s), failed(s), skipped(s)))
|
||||
npassed += lnpassed
|
||||
nfailed += lnfailed
|
||||
nskipped += lnskipped
|
||||
end
|
||||
println("""DOCTEST: {"failed": $nfailed, "skipped": $nskipped, "passed": $npassed}""")
|
||||
make_badge("doctests", npassed, (npassed+nfailed+nskipped), "badges/doctests-status.svg")
|
||||
end
|
||||
@@ -1,17 +0,0 @@
|
||||
using Lint
|
||||
|
||||
include("badges.jl")
|
||||
|
||||
d = lintpkg("JuliaFEM", returnMsgs=true)
|
||||
cd(dirname(@__FILE__)) do
|
||||
k = 0
|
||||
open("quality/lint_report.rst", "w") do fid
|
||||
for i in d
|
||||
i = string(i)
|
||||
i = i[search(i, "JuliaFEM")[2]-1:end]
|
||||
write(fid, "| "*i*"\n")
|
||||
k += 1
|
||||
end
|
||||
end
|
||||
make_badge("code quality", k, "badges/lint-status.svg")
|
||||
end
|
||||
@@ -1,231 +0,0 @@
|
||||
# Automatically run notebooks and generate rst table for results
|
||||
|
||||
include("badges.jl")
|
||||
|
||||
""" rst file parser to find title, author and factcheck status.
|
||||
"""
|
||||
function parse_rst(filename)
|
||||
res = Dict("author" => "unknown", "status" => 0, "title" => "", "abstract" => "")
|
||||
title_found = false
|
||||
open(filename) do fid
|
||||
data = readlines(fid)
|
||||
for line in data
|
||||
line = strip(line)
|
||||
if line == ""
|
||||
continue
|
||||
end
|
||||
if !title_found
|
||||
res["title"] = line
|
||||
title_found = true
|
||||
end
|
||||
if startswith(line, "Author(s):")
|
||||
auth = split(line[12:end], ' ')
|
||||
auth = filter(s -> !('@' in s), auth)
|
||||
auth = join(auth, ' ')
|
||||
res["author"] = auth
|
||||
end
|
||||
if startswith(line, "Abstract:") # not working, todo, fixme, ...
|
||||
res["abstract"] = line[11:end]
|
||||
end
|
||||
if startswith(line, "Failed:")
|
||||
res["status"] = 1
|
||||
end
|
||||
if startswith(line, "Failure")
|
||||
res["status"] = 1
|
||||
end
|
||||
m = matchall(r"[-0-9.]+", line)
|
||||
end
|
||||
end
|
||||
return res
|
||||
end
|
||||
|
||||
|
||||
function run_notebooks()
|
||||
k = 0
|
||||
results = Dict[]
|
||||
cd(dirname(@__FILE__)) do
|
||||
for ipynb in readdir("tutorials")
|
||||
tic()
|
||||
if !endswith(ipynb, "ipynb")
|
||||
continue
|
||||
end
|
||||
runtime = 0
|
||||
status = 1
|
||||
println("Running notebook tutorials/$ipynb")
|
||||
# port = 34211+k # we're having some weird port issue with zmq
|
||||
# k += 1
|
||||
try
|
||||
run(`timeout 180 runipy -o tutorials/$ipynb --kernel=julia-0.4`)
|
||||
status = 0
|
||||
catch error
|
||||
warn("running notebook failed")
|
||||
Base.showerror(Base.STDOUT, error)
|
||||
end
|
||||
runtime = toc()
|
||||
bn = "tutorials/$(ipynb[1:end-6])"
|
||||
#try
|
||||
# run(`ipython nbconvert tutorials/$ipynb --to rst --output=$bn`)
|
||||
#catch error
|
||||
# warn("unable to convert notebook to rst format")
|
||||
# Base.showerror(Base.STDOUT, error)
|
||||
#end
|
||||
|
||||
try
|
||||
run(`ipython nbconvert tutorials/$ipynb --to html --output=$bn`)
|
||||
catch error
|
||||
warn("unable to convert notebook to html format")
|
||||
Base.showerror(Base.STDOUT, error)
|
||||
end
|
||||
|
||||
try
|
||||
run(`ipython nbconvert tutorials/$ipynb --to latex --output=$bn`)
|
||||
catch error
|
||||
warn("unable to convert notebook to tex format")
|
||||
Base.showerror(Base.STDOUT, error)
|
||||
end
|
||||
|
||||
try
|
||||
run(`lualatex --output-directory=tutorials $bn.tex`)
|
||||
catch error
|
||||
warn("unable to convert notebook from tex to pdf")
|
||||
Base.showerror(Base.STDOUT, error)
|
||||
end
|
||||
|
||||
data = Dict("author" => "unknown", "status" => status, "runtime" => runtime,
|
||||
"filename" => ipynb, "last_run" => time(), "description"=>"")
|
||||
res = parse_rst("$bn.rst")
|
||||
data["description"] = res["title"]
|
||||
data["author"] = res["author"]
|
||||
# there is two possibilities, either notbook does not run for syntax error (status = 1)
|
||||
# or notebook is fine but factcheck returns failed. This checks the latter one.
|
||||
if status == 0
|
||||
data["status"] = res["status"]
|
||||
end
|
||||
push!(results, data)
|
||||
end
|
||||
end
|
||||
return results
|
||||
end
|
||||
|
||||
"""
|
||||
Make rows from results ready to tabular form
|
||||
"""
|
||||
function makerows(results)
|
||||
rows = ["id" "author" "description" "status" "ipynb" "pdf" "last run" "runtime (s)"]
|
||||
for (i, data) in enumerate(results)
|
||||
row = ["na" "unknown" "unknown" "unknown" "unknown" "unknown" "na" "na"]
|
||||
println(i)
|
||||
row[1] = string(i)
|
||||
desc = data["description"]
|
||||
if desc == ""
|
||||
desc = data["filename"]
|
||||
end
|
||||
row[2] = data["author"]
|
||||
fn = data["filename"][1:end-6]
|
||||
row[3] = ":doc:`$desc <$fn>`"
|
||||
|
||||
if data["status"] == 0
|
||||
row[4] = ".. image:: /badges/notebook-passing.svg"
|
||||
else
|
||||
row[4] = ".. image:: /badges/notebook-failing.svg"
|
||||
end
|
||||
|
||||
fn = data["filename"]
|
||||
row[5] = ":download:`ipynb <$fn>`"
|
||||
fn = data["filename"][1:end-6]*".pdf"
|
||||
row[6] = ":download:`pdf <$fn>`"
|
||||
|
||||
row[7] = Libc.strftime("%Y-%m-%d %H:%M:%S", data["last_run"])
|
||||
row[8] = string(round(data["runtime"], 2))
|
||||
rows = vcat(rows, row)
|
||||
end
|
||||
return rows
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
write rst table. rows is array of arrays.
|
||||
|
||||
Expected output
|
||||
|
||||
+----+-----------+------------------------------+------------+---------+--------+------+-----+
|
||||
| id | author | description | last run | runtime | status | html | pdf |
|
||||
+====+===========+==============================+============+=========+========+======+=====+
|
||||
| 1 | Jukka Aho | This is placeholder notebook | 2015-07-05 | 123.45 | |pass| | html | pdf |
|
||||
+----+-----------+------------------------------+------------+---------+--------+------+-----+
|
||||
|
||||
|
||||
"""
|
||||
function write_rst_table(rows)
|
||||
# determine cell lenghts
|
||||
lenghts = zeros(Int, size(rows))
|
||||
for i=1:size(rows, 1)
|
||||
for j=1:size(rows, 2)
|
||||
lenghts[i,j] = length(rows[i,j])+3
|
||||
end
|
||||
end
|
||||
#println(lenghts)
|
||||
cl = maximum(lenghts, 1) # maximum column lenghts
|
||||
tl = sum(cl)+1 # total table width
|
||||
mp = [1 cumsum(cl, 2)+1] # marker points "+"
|
||||
println(mp)
|
||||
#println(cl)
|
||||
#println(tl)
|
||||
|
||||
function sep(l, m; s="+", d="-")
|
||||
p = ""
|
||||
for j=1:l
|
||||
if j in m
|
||||
p *= s
|
||||
else
|
||||
p *= d
|
||||
end
|
||||
end
|
||||
p *= "\n"
|
||||
return p
|
||||
end
|
||||
|
||||
s = ""
|
||||
# write rows
|
||||
s *= sep(tl, mp)
|
||||
for i = 1:size(rows, 1)
|
||||
s *= "|"
|
||||
for j=1:size(rows, 2)
|
||||
s *= " " * rows[i, j] * " "
|
||||
s *= repeat(" ", cl[j]-lenghts[i,j])
|
||||
s *= "|"
|
||||
end
|
||||
s *= "\n"
|
||||
if i == 1
|
||||
s *= sep(tl, mp; d="=")
|
||||
else
|
||||
s *= sep(tl, mp)
|
||||
end
|
||||
end
|
||||
return s
|
||||
end
|
||||
|
||||
|
||||
function main()
|
||||
results = run_notebooks()
|
||||
notebooks_total = 0
|
||||
notebooks_passing = 0
|
||||
for result in results
|
||||
notebooks_total += 1
|
||||
if result["status"] == 0
|
||||
notebooks_passing += 1
|
||||
end
|
||||
end
|
||||
rows = makerows(results)
|
||||
table = write_rst_table(rows)
|
||||
println(table)
|
||||
cd(dirname(@__FILE__)) do
|
||||
open("tutorials/notebooks.rst", "w") do fid
|
||||
write(fid, table)
|
||||
end
|
||||
make_badge("notebooks", notebooks_passing, notebooks_total, "badges/notebooks-status.svg")
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
main()
|
||||
@@ -1,10 +0,0 @@
|
||||
import FactCheck
|
||||
import Logging
|
||||
include("badges.jl")
|
||||
|
||||
cd(dirname(@__FILE__)) do
|
||||
include("../test/runtests.jl")
|
||||
s = FactCheck.getstats()
|
||||
make_badge("unittests", s["nSuccesses"], sum(values(s)), "badges/unittests-status.svg")
|
||||
end
|
||||
|
||||
@@ -1,64 +0,0 @@
|
||||
---
|
||||
title: "JuliaFEM Contributor Manual"
|
||||
description: "Technical guide for developers and contributors"
|
||||
date: 2025-11-09
|
||||
author: "Jukka Aho"
|
||||
categories: ["development", "contributor guide"]
|
||||
keywords: ["juliafem", "development", "architecture", "testing", "performance"]
|
||||
audience: "developers"
|
||||
level: "advanced"
|
||||
type: "manual"
|
||||
---
|
||||
|
||||
**Audience:** Developers, contributors, advanced users who want to extend or modify JuliaFEM.
|
||||
|
||||
This manual is **technical and detailed** - it explains HOW the code works and WHY we made certain design choices.
|
||||
|
||||
## What's Here
|
||||
|
||||
- **Testing Philosophy:** How and why we test
|
||||
- **Coding Standards:** Required conventions for all contributions (variable names, types, performance)
|
||||
- **Architecture:** Module structure, data flow, key abstractions
|
||||
- **Performance:** Zero-allocation design, profiling, benchmarking
|
||||
- **Adding Elements:** How to implement new element types
|
||||
- **CI/CD:** Continuous integration, releases, versioning
|
||||
- **Git Workflow:** Branching, commits, pull requests
|
||||
|
||||
## What's NOT Here
|
||||
|
||||
- User tutorials (see `docs/user/` for that)
|
||||
- Deep mathematical theory (see `docs/book/` for that)
|
||||
- "How do I solve problem X?" (that's user docs)
|
||||
|
||||
## Philosophy
|
||||
|
||||
**"Show me the code AND tell me why."**
|
||||
|
||||
We assume you:
|
||||
|
||||
- Know Julia reasonably well
|
||||
- Understand FEM basics
|
||||
- Want to add features or fix bugs
|
||||
- Care about performance and correctness
|
||||
- Need to understand design rationale
|
||||
|
||||
## Before Contributing
|
||||
|
||||
1. Read [Testing Philosophy](testing_philosophy.md)
|
||||
2. Follow [Coding Standards](coding_standards.md) - **REQUIRED** for all contributions
|
||||
3. Understand [Architecture](architecture.md)
|
||||
4. Check [Performance Guidelines](performance.md)
|
||||
5. Review [Git Workflow](git_workflow.md)
|
||||
|
||||
## Key Principles
|
||||
|
||||
- **Type stability:** No `Any`, no `Dict` without types
|
||||
- **Zero allocations:** Hot paths should allocate nothing
|
||||
- **Immutability:** Prefer `struct` over `mutable struct`
|
||||
- **Composition:** Use tuples and free functions, not OOP hierarchies
|
||||
- **Explicit:** No magic, user knows what happens
|
||||
- **Test first:** Write tests before fixing bugs
|
||||
|
||||
---
|
||||
|
||||
**Start here:** [Testing Philosophy](testing_philosophy.md) | [Architecture Overview](architecture.md)
|
||||
@@ -1,499 +0,0 @@
|
||||
---
|
||||
title: "JuliaFEM Coding Standards"
|
||||
description: "Required coding conventions and style guide for all contributors"
|
||||
date: 2025-11-09
|
||||
author: "Jukka Aho"
|
||||
categories: ["development", "standards", "style guide"]
|
||||
keywords: ["juliafem", "coding standards", "style", "conventions", "best practices"]
|
||||
audience: "developers"
|
||||
level: "required"
|
||||
type: "standards"
|
||||
status: "active"
|
||||
---
|
||||
|
||||
This document defines the coding standards and conventions for JuliaFEM development.
|
||||
|
||||
**Last Updated:** November 9, 2025
|
||||
**Status:** Active - all new code must follow these standards
|
||||
|
||||
---
|
||||
|
||||
## Core Principles
|
||||
|
||||
1. **Readability over cleverness** - Code should be understandable by FEM practitioners
|
||||
2. **Type stability first** - Performance depends on it (100x difference)
|
||||
3. **Zero allocations in hot paths** - Profiling required
|
||||
4. **Explicit over implicit** - No magic, show what happens
|
||||
5. **Composition over inheritance** - Structs and free functions, not OOP
|
||||
|
||||
---
|
||||
|
||||
## Variable Naming Conventions
|
||||
|
||||
### No Greek Letters in Code (Critical!)
|
||||
|
||||
**Rule:** Never use Greek letters (ξ, η, ζ, α, β, γ, etc.) in code.
|
||||
|
||||
**Rationale:**
|
||||
|
||||
- **Keyboard accessibility** - Not all keyboards support Greek input
|
||||
- **Editor compatibility** - Some editors struggle with Unicode math symbols
|
||||
- **Copy-paste issues** - Greek letters cause encoding problems
|
||||
- **Search/replace problems** - Text tools may not handle Unicode correctly
|
||||
- **Terminal rendering** - SSH sessions may not display correctly
|
||||
- **Internationalization** - Non-Western keyboards make editing difficult
|
||||
- **Git diffs** - Unicode can cause merge conflicts or display issues
|
||||
- **Accessibility** - Screen readers struggle with Greek letters
|
||||
|
||||
**Correct:**
|
||||
|
||||
```julia
|
||||
# Reference element coordinates
|
||||
function eval_basis(u::Float64, v::Float64, w::Float64)
|
||||
# Natural coordinates u, v, w ∈ [-1, 1]
|
||||
N1 = (1 - u) * (1 - v) * (1 - w) / 8
|
||||
return N1
|
||||
end
|
||||
|
||||
# Physical coordinates
|
||||
function map_to_physical(x::Vec, y::Vec, z::Vec, u::Float64, v::Float64, w::Float64)
|
||||
# Map from natural (u,v,w) to physical (x,y,z)
|
||||
end
|
||||
```
|
||||
|
||||
**Incorrect:**
|
||||
|
||||
```julia
|
||||
# ❌ DON'T DO THIS
|
||||
function eval_basis(ξ::Float64, η::Float64, ζ::Float64)
|
||||
N1 = (1 - ξ) * (1 - η) * (1 - ζ) / 8
|
||||
return N1
|
||||
end
|
||||
```
|
||||
|
||||
**Exception:** Greek letters are acceptable in:
|
||||
|
||||
- **Comments** - Mathematical notation for clarity: `# Shape function: Nᵢ(ξ)`
|
||||
- **Documentation** - Latex math blocks: `$\\xi \\in [-1, 1]$`
|
||||
- **String literals** - Plot labels: `xlabel="ξ coordinate"`
|
||||
- **Error messages** - User-facing text: `"Invalid ξ coordinate"`
|
||||
|
||||
**Standard variable names:**
|
||||
|
||||
- Reference coordinates: `u`, `v`, `w` (not ξ, η, ζ)
|
||||
- Physical coordinates: `x`, `y`, `z`
|
||||
- Derivatives: `du`, `dv`, `dw` or `dudx`, `dudy`, etc.
|
||||
- Jacobian: `J` or `jac` (not ∂)
|
||||
- Determinant: `detJ` (not |J|)
|
||||
- Inverse: `invJ` or `Jinv` (not J⁻¹)
|
||||
|
||||
---
|
||||
|
||||
## Type Naming
|
||||
|
||||
### Structs and Types
|
||||
|
||||
- **PascalCase** - All type names: `Element`, `Problem`, `Material`
|
||||
- **No abbreviations** - `Quadrilateral` not `Quad` (unless established convention)
|
||||
- **Descriptive suffixes** - Purpose clear from name
|
||||
|
||||
**Basis types** - Append "Basis" suffix to distinguish from topology:
|
||||
|
||||
```julia
|
||||
# ✅ Correct - No name collision
|
||||
struct Tri3Basis <: AbstractBasis{2} end # Interpolation scheme
|
||||
struct Tri3 <: AbstractTopology end # Element geometry
|
||||
|
||||
# ❌ Wrong - Name collision!
|
||||
struct Tri3 <: AbstractBasis{2} end
|
||||
struct Tri3 <: AbstractTopology end # ERROR: type Tri3 already defined
|
||||
```
|
||||
|
||||
### Functions
|
||||
|
||||
- **snake_case** - All function names: `assemble_element`, `solve_static`
|
||||
- **Verb-noun pattern** - Action clear: `compute_stiffness`, `evaluate_basis`
|
||||
- **Boolean predicates** - `is_*` or `has_*`: `is_converged`, `has_contact`
|
||||
|
||||
### Constants
|
||||
|
||||
- **UPPER_CASE** - Module-level constants: `MAX_ITERATIONS`, `TOLERANCE`
|
||||
- **Type-stable** - Always specify type: `const MAX_ITER::Int = 100`
|
||||
|
||||
### Internal/Private
|
||||
|
||||
- **Underscore prefix** - Not exported: `_compute_internal_forces`
|
||||
- **Not API** - Can change between versions
|
||||
|
||||
---
|
||||
|
||||
## Code Organization
|
||||
|
||||
### File Structure
|
||||
|
||||
```julia
|
||||
# Standard file header
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE
|
||||
|
||||
"""
|
||||
Brief one-line description of file purpose.
|
||||
|
||||
Extended description if needed.
|
||||
"""
|
||||
|
||||
# Imports (grouped)
|
||||
using LinearAlgebra
|
||||
using SparseArrays
|
||||
|
||||
# Local imports
|
||||
using ..JuliaFEM: AbstractElement, AbstractBasis
|
||||
|
||||
# Type definitions
|
||||
struct MyType
|
||||
# ...
|
||||
end
|
||||
|
||||
# Function implementations
|
||||
function my_function(args)
|
||||
# ...
|
||||
end
|
||||
```
|
||||
|
||||
### Import Style
|
||||
|
||||
```julia
|
||||
# ✅ Correct - Explicit imports
|
||||
using LinearAlgebra: norm, dot, cross
|
||||
using SparseArrays: sparse, spzeros
|
||||
|
||||
# ❌ Avoid - Blanket imports (pollutes namespace)
|
||||
using LinearAlgebra
|
||||
using SparseArrays
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Performance Guidelines
|
||||
|
||||
### Type Stability
|
||||
|
||||
```julia
|
||||
# ✅ Type-stable - Return type inferrable
|
||||
function compute_mass(element::Quad4, density::Float64)
|
||||
m::Float64 = 0.0
|
||||
# ...
|
||||
return m
|
||||
end
|
||||
|
||||
# ❌ Type-unstable - Return type changes!
|
||||
function compute_mass(element, density)
|
||||
if density > 0
|
||||
return 1.0 # Float64
|
||||
else
|
||||
return nothing # Nothing - type-unstable!
|
||||
end
|
||||
end
|
||||
```
|
||||
|
||||
### Zero Allocations
|
||||
|
||||
```julia
|
||||
# ✅ Pre-allocated cache
|
||||
struct AssemblyCache
|
||||
K_local::Matrix{Float64}
|
||||
f_local::Vector{Float64}
|
||||
end
|
||||
|
||||
function assemble!(cache::AssemblyCache, element)
|
||||
fill!(cache.K_local, 0.0)
|
||||
# Reuse cache.K_local - no allocations
|
||||
end
|
||||
|
||||
# ❌ Allocates every call
|
||||
function assemble(element)
|
||||
K_local = zeros(8, 8) # Allocates!
|
||||
return K_local
|
||||
end
|
||||
```
|
||||
|
||||
### Tuple Returns (Zero Allocation)
|
||||
|
||||
```julia
|
||||
# ✅ Tuple return - no allocation
|
||||
function eval_basis(element::Tri3, u::Float64, v::Float64)
|
||||
N1 = 1 - u - v
|
||||
N2 = u
|
||||
N3 = v
|
||||
return (N1, N2, N3) # NTuple{3,Float64} - stack allocated
|
||||
end
|
||||
|
||||
# ❌ Vector return - allocates!
|
||||
function eval_basis(element::Tri3, u::Float64, v::Float64)
|
||||
return [1 - u - v, u, v] # Vector{Float64} - heap allocation
|
||||
end
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Documentation Style
|
||||
|
||||
### Docstrings
|
||||
|
||||
Use Julia's docstring format with standard sections:
|
||||
|
||||
`````markdown
|
||||
"""
|
||||
assemble_element(element::Quad4, u::Vector{Float64}) -> Matrix{Float64}
|
||||
|
||||
Assemble element stiffness matrix for 4-node quadrilateral element.
|
||||
|
||||
# Arguments
|
||||
- `element::Quad4`: Quadrilateral element with nodal connectivity
|
||||
- `u::Vector{Float64}`: Nodal displacement vector (8 DOFs: u1,v1,u2,v2,...)
|
||||
|
||||
# Returns
|
||||
- `K::Matrix{Float64}`: 8×8 element stiffness matrix in global DOFs
|
||||
|
||||
# Theory
|
||||
Uses 2×2 Gauss quadrature with bilinear shape functions:
|
||||
```math
|
||||
K_{ij} = \\int_{\\Omega_e} B_i^T D B_j \\, d\\Omega
|
||||
```
|
||||
|
||||
# Example
|
||||
```julia
|
||||
nodes = [Node(0,0), Node(1,0), Node(1,1), Node(0,1)]
|
||||
element = Quad4(nodes)
|
||||
u = zeros(8)
|
||||
K = assemble_element(element, u)
|
||||
```
|
||||
|
||||
# Performance
|
||||
This function allocates. For zero-allocation assembly, use `assemble_element!`
|
||||
with pre-allocated cache.
|
||||
"""
|
||||
function assemble_element(element::Quad4, u::Vector{Float64})
|
||||
# Implementation
|
||||
end
|
||||
`````
|
||||
|
||||
### Comments
|
||||
|
||||
```julia
|
||||
# ✅ Good comments - Why, not what
|
||||
# Use RCM ordering to minimize bandwidth (10x faster solve)
|
||||
perm = rcm_permutation(mesh)
|
||||
|
||||
# Check convergence: ||Δu|| < ε||u||
|
||||
# Relative norm prevents scale-dependent tolerance
|
||||
if norm(Δu) < tol * norm(u)
|
||||
break
|
||||
end
|
||||
|
||||
# ❌ Bad comments - Obvious from code
|
||||
# Loop over elements
|
||||
for element in elements
|
||||
# Add to global matrix
|
||||
K_global += K_local
|
||||
end
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Testing Standards
|
||||
|
||||
### Test Organization
|
||||
|
||||
```julia
|
||||
@testset "Quad4 element" begin
|
||||
@testset "Stiffness matrix" begin
|
||||
# Unit square, E=1, ν=0.3
|
||||
element = Quad4([Node(0,0), Node(1,0), Node(1,1), Node(0,1)])
|
||||
K = stiffness_matrix(element, E=1.0, ν=0.3)
|
||||
|
||||
# Symmetry
|
||||
@test issymmetric(K)
|
||||
|
||||
# Positive definite (after BC)
|
||||
@test all(eigvals(K[3:end, 3:end]) .> 0)
|
||||
end
|
||||
|
||||
@testset "Patch test" begin
|
||||
# Linear displacement field must be exact
|
||||
# ... validation test ...
|
||||
end
|
||||
end
|
||||
```
|
||||
|
||||
### Floating Point Comparisons
|
||||
|
||||
```julia
|
||||
# ✅ Use tolerances
|
||||
@test result ≈ expected atol=1e-10
|
||||
@test isapprox(result, expected, rtol=1e-6)
|
||||
|
||||
# ❌ Never exact equality for floats
|
||||
@test result == expected # Fragile!
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Anti-Patterns (Don't Do This!)
|
||||
|
||||
### 1. Dict Without Type Parameters
|
||||
|
||||
```julia
|
||||
# ❌ Type-unstable Dict (100x slower!)
|
||||
fields = Dict("displacement" => u, "velocity" => v)
|
||||
|
||||
# ✅ Type-stable alternative
|
||||
fields = (displacement=u, velocity=v) # NamedTuple
|
||||
# or
|
||||
struct Fields
|
||||
displacement::Vector{Float64}
|
||||
velocity::Vector{Float64}
|
||||
end
|
||||
```
|
||||
|
||||
### 2. Abstract Types in Structs
|
||||
|
||||
```julia
|
||||
# ❌ Type-unstable struct
|
||||
struct Element
|
||||
nodes::AbstractVector # Type-unstable!
|
||||
end
|
||||
|
||||
# ✅ Parametric struct
|
||||
struct Element{N}
|
||||
nodes::NTuple{N, Node} # Type-stable, zero-allocation
|
||||
end
|
||||
```
|
||||
|
||||
### 3. Global Variables
|
||||
|
||||
```julia
|
||||
# ❌ Global mutable state
|
||||
global_stiffness = zeros(1000, 1000)
|
||||
|
||||
function assemble!(element)
|
||||
global_stiffness .+= K_local # Spooky action at a distance!
|
||||
end
|
||||
|
||||
# ✅ Explicit parameters
|
||||
function assemble!(K_global::Matrix, element)
|
||||
K_global .+= K_local # Clear data flow
|
||||
end
|
||||
```
|
||||
|
||||
### 4. Type Piracy
|
||||
|
||||
```julia
|
||||
# ❌ Extending methods on types you don't own
|
||||
Base.+(a::Vector, b::Matrix) = ... # DON'T!
|
||||
|
||||
# ✅ Wrapper type or different function
|
||||
struct MyVector
|
||||
data::Vector
|
||||
end
|
||||
Base.+(a::MyVector, b::Matrix) = ... # OK - our type
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Git Commit Style
|
||||
|
||||
### Commit Messages
|
||||
|
||||
```text
|
||||
feat(topology): Add Pyr5 pyramid element topology
|
||||
|
||||
- Implement 5-node pyramid reference element
|
||||
- Add connectivity information (faces, edges)
|
||||
- Zero-allocation tuple interface
|
||||
- Tests: reference coordinates, edge/face queries
|
||||
|
||||
Closes #123
|
||||
```
|
||||
|
||||
**Format:** `<type>(<scope>): <subject>`
|
||||
|
||||
**Types:**
|
||||
|
||||
- `feat`: New feature
|
||||
- `fix`: Bug fix
|
||||
- `docs`: Documentation only
|
||||
- `refactor`: Code restructuring (no behavior change)
|
||||
- `perf`: Performance improvement
|
||||
- `test`: Adding tests
|
||||
- `chore`: Tooling, dependencies
|
||||
|
||||
**Scope:** Module or component (topology, assembly, solver, etc.)
|
||||
|
||||
**Subject:**
|
||||
|
||||
- Imperative mood: "Add feature" not "Added feature"
|
||||
- No period at end
|
||||
- Max 72 characters
|
||||
|
||||
---
|
||||
|
||||
## Editor Configuration
|
||||
|
||||
### Recommended Settings
|
||||
|
||||
```julia
|
||||
# .editorconfig
|
||||
[*.jl]
|
||||
charset = utf-8
|
||||
end_of_line = lf
|
||||
indent_style = space
|
||||
indent_size = 4
|
||||
trim_trailing_whitespace = true
|
||||
insert_final_newline = true
|
||||
```
|
||||
|
||||
### JuliaFormatter.jl
|
||||
|
||||
```julia
|
||||
# .JuliaFormatter.toml
|
||||
indent = 4
|
||||
margin = 92
|
||||
always_for_in = true
|
||||
whitespace_typedefs = true
|
||||
whitespace_ops_in_indices = true
|
||||
remove_extra_newlines = true
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Summary Checklist
|
||||
|
||||
Before submitting code, verify:
|
||||
|
||||
- [ ] No Greek letters in variable names (use u, v, w)
|
||||
- [ ] All types are PascalCase
|
||||
- [ ] All functions are snake_case
|
||||
- [ ] Type-stable (check with `@code_warntype`)
|
||||
- [ ] Zero allocations in hot paths (check with `@allocations` or `@btime`)
|
||||
- [ ] Docstrings for exported functions
|
||||
- [ ] Tests pass locally
|
||||
- [ ] Comments explain "why" not "what"
|
||||
- [ ] No type piracy
|
||||
- [ ] No global mutable state
|
||||
|
||||
---
|
||||
|
||||
## References
|
||||
|
||||
- **Performance tips:** https://docs.julialang.org/en/v1/manual/performance-tips/
|
||||
- **Style guide:** https://docs.julialang.org/en/v1/manual/style-guide/
|
||||
- **JuliaFEM vision:** `llm/VISION_2.0.md`
|
||||
- **Architecture:** `llm/ARCHITECTURE.md`
|
||||
- **Technical lessons:** `llm/TECHNICAL_VISION.md`
|
||||
|
||||
---
|
||||
|
||||
**Enforcement:** These standards are enforced through code review. All PRs must follow these conventions.
|
||||
|
||||
**Evolution:** This document evolves with the project. Propose changes via pull request.
|
||||
@@ -1,119 +0,0 @@
|
||||
---
|
||||
title: "Quick Reference: GPU Elasticity Solver"
|
||||
date: 2025-11-10
|
||||
author: "JuliaFEM Team"
|
||||
status: "Authoritative"
|
||||
last_updated: 2025-11-10
|
||||
tags: ["gpu", "elasticity", "quickstart", "guide"]
|
||||
---
|
||||
|
||||
**Ready to use!** Complete implementation with tests.
|
||||
|
||||
---
|
||||
|
||||
## 🚀 Quick Start
|
||||
|
||||
### 1. Run Demo
|
||||
```bash
|
||||
cd /home/juajukka/dev/JuliaFEM.jl
|
||||
julia --project=. demos/cantilever_beam_demo.jl
|
||||
```
|
||||
|
||||
This will:
|
||||
- Generate cantilever mesh (10×1×1 beam)
|
||||
- Solve on GPU
|
||||
- Compare with analytical solution
|
||||
|
||||
### 2. Run Tests
|
||||
```bash
|
||||
cd test
|
||||
julia --project=.. test_gpu_elasticity.jl
|
||||
```
|
||||
|
||||
This validates:
|
||||
- Fixed boundary conditions
|
||||
- Deflection pattern
|
||||
- Analytical comparison
|
||||
|
||||
---
|
||||
|
||||
## 📁 Key Files
|
||||
|
||||
**Solver:** `src/gpu_elasticity.jl` (550 lines)
|
||||
- Main module with GPU kernels
|
||||
- CG solver
|
||||
- BC handling
|
||||
|
||||
**Mesh Generator:** `scripts/generate_cantilever_mesh.jl`
|
||||
- Creates test geometry with Gmsh
|
||||
|
||||
**Demo:** `demos/cantilever_beam_demo.jl`
|
||||
- Complete workflow example
|
||||
|
||||
**Tests:** `test/test_gpu_elasticity.jl`
|
||||
- Validation suite
|
||||
|
||||
**Docs:** `docs/design/GPU_ELASTICITY_IMPLEMENTATION.md`
|
||||
- Complete guide
|
||||
|
||||
---
|
||||
|
||||
## 🎯 What We Built
|
||||
|
||||
✅ **Complete GPU solver** - Two-phase nodal assembly
|
||||
✅ **Tensors.jl on GPU** - Natural tensor operations
|
||||
✅ **No atomics** - Node-parallel, no race conditions
|
||||
✅ **Matrix-free** - Lower memory, recompute geometry
|
||||
✅ **Test suite** - Cantilever beam validation
|
||||
✅ **Gmsh integration** - Automated mesh generation
|
||||
|
||||
---
|
||||
|
||||
## 📊 Expected Results
|
||||
|
||||
**Cantilever Beam (10×1×1 m, Steel, 1 MPa pressure):**
|
||||
- Max displacement: ~1e-4 m at free end
|
||||
- CG iterations: 50-100 (no preconditioning)
|
||||
- Analytical match: within 10-30%
|
||||
|
||||
---
|
||||
|
||||
## 🔧 Usage Example
|
||||
|
||||
```julia
|
||||
using GPUElasticity
|
||||
|
||||
# Read mesh
|
||||
mesh = read_gmsh_mesh("cantilever_beam.msh")
|
||||
|
||||
# Material (steel)
|
||||
material = ElasticMaterial(210e9, 0.3)
|
||||
|
||||
# Boundary conditions
|
||||
fixed = get_surface_nodes(mesh, "FixedEnd")
|
||||
pressure = get_surface_nodes(mesh, "PressureSurface")
|
||||
|
||||
# Solve
|
||||
problem = ElasticityProblem(mesh, material, fixed, pressure, 1e6)
|
||||
u = solve_elasticity_gpu(problem)
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## 🎯 Next Steps
|
||||
|
||||
1. **Test on GPU** - Run demo and tests
|
||||
2. **Add preconditioning** - Target 10-20 CG iters
|
||||
3. **Extend to nonlinear** - Plasticity + Newton-Krylov
|
||||
|
||||
---
|
||||
|
||||
## 📚 Documentation
|
||||
|
||||
- `docs/design/GPU_ELASTICITY_IMPLEMENTATION.md` - Full guide
|
||||
- `docs/design/gpu_nodal_assembly_architecture.md` - Architecture
|
||||
- `llm/sessions/2025-11-10_gpu_elasticity_implementation.md` - Session notes
|
||||
|
||||
---
|
||||
|
||||
**Everything is ready to test! 🚀**
|
||||
@@ -1,262 +0,0 @@
|
||||
---
|
||||
title: "GPU Nodal Assembly - Quick Start Guide"
|
||||
date: 2025-11-10
|
||||
author: "JuliaFEM Team"
|
||||
status: "Authoritative"
|
||||
last_updated: 2025-11-10
|
||||
tags: ["gpu", "nodal-assembly", "nonlinear", "quickstart", "guide"]
|
||||
---
|
||||
|
||||
**Status:** CPU ✅ Working | GPU 🔄 Ready to Test
|
||||
|
||||
---
|
||||
|
||||
## What is This?
|
||||
|
||||
A **complete GPU-resident nonlinear FEM solver** using:
|
||||
- **Nodal assembly** (matrix-free, no atomics)
|
||||
- **Tensors.jl** (natural tensor operations on GPU)
|
||||
- **Two-phase pipeline** (GP data → nodal assembly)
|
||||
- **Perfect plasticity** (von Mises with return mapping)
|
||||
|
||||
---
|
||||
|
||||
## Quick Test (CPU Reference)
|
||||
|
||||
```bash
|
||||
cd /home/juajukka/dev/JuliaFEM.jl
|
||||
julia demos/nodal_assembly_cpu.jl
|
||||
```
|
||||
|
||||
**Expected output:**
|
||||
```
|
||||
Residual norm: 727.2081516082287
|
||||
Material States: Plastic (α = 5.634921e-03) at all GPs
|
||||
Force Balance: ✅ PASSED
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Quick Test (GPU)
|
||||
|
||||
```bash
|
||||
cd /home/juajukka/dev/JuliaFEM.jl
|
||||
julia --project=. demos/nodal_assembly_gpu.jl
|
||||
```
|
||||
|
||||
**Requirements:**
|
||||
- CUDA-capable GPU
|
||||
- CUDA.jl installed
|
||||
|
||||
**Expected output:**
|
||||
- Residual norm should match CPU: ~727.2
|
||||
- All material states plastic
|
||||
- Force balance passed
|
||||
|
||||
---
|
||||
|
||||
## Architecture Overview
|
||||
|
||||
### Two-Phase Pipeline
|
||||
|
||||
```
|
||||
Phase 1: Integration Point Data (GP Kernel)
|
||||
Input: u, nodes, elements, states_old
|
||||
Output: σ_gp (stresses), states_new
|
||||
Parallelism: One thread per GP
|
||||
|
||||
↓
|
||||
|
||||
Phase 2: Nodal Assembly (Node Kernel)
|
||||
Input: σ_gp, nodes, elements, node_to_elems (CSR)
|
||||
Output: r (residual vector)
|
||||
Parallelism: One thread per node
|
||||
NO ATOMICS NEEDED!
|
||||
```
|
||||
|
||||
### Key Data Structures
|
||||
|
||||
```julia
|
||||
# Stresses (Tensors.jl on GPU!)
|
||||
σ_gp = CuArray{SymmetricTensor{2,3,Float64,6}, 1}
|
||||
|
||||
# Material states
|
||||
states = CuArray{PlasticState, 1}
|
||||
|
||||
# CSR map (which elements touch each node)
|
||||
struct NodeToElementsMap
|
||||
ptr::CuArray{Int32, 1}
|
||||
data::CuArray{Int32, 1}
|
||||
end
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Files to Know
|
||||
|
||||
### Documentation
|
||||
- **`docs/design/gpu_nodal_assembly_architecture.md`** - Complete architecture (500+ lines)
|
||||
- **`llm/sessions/2025-11-10_gpu_nodal_assembly_complete.md`** - Session summary
|
||||
|
||||
### Implementation
|
||||
- **`demos/nodal_assembly_cpu.jl`** - CPU reference (400+ lines, ✅ working)
|
||||
- **`demos/nodal_assembly_gpu.jl`** - GPU version (450+ lines, ready to test)
|
||||
|
||||
### Background
|
||||
- **`demos/newton_krylov_anderson_cpu.jl`** - Complete Newton-Krylov solver
|
||||
- **`docs/design/gpu_solver_strategy_expert_validated.md`** - Expert-validated strategy
|
||||
|
||||
---
|
||||
|
||||
## Why Nodal Assembly?
|
||||
|
||||
### ❌ Element-Based (Standard GPU FEM)
|
||||
```julia
|
||||
for elem in elements
|
||||
compute element forces
|
||||
CUDA.@atomic r[node] += f_elem[i] # ATOMIC - CONTENTION!
|
||||
end
|
||||
```
|
||||
|
||||
### ✅ Node-Based (Our Approach)
|
||||
```julia
|
||||
for node in nodes # Each thread owns ONE node
|
||||
for elem in elements_touching_node
|
||||
f_node += contribution from elem
|
||||
end
|
||||
r[node] = f_node # DIRECT WRITE - NO ATOMICS!
|
||||
end
|
||||
```
|
||||
|
||||
**Benefits:**
|
||||
- No atomic operations (faster!)
|
||||
- Matrix-free (lower memory)
|
||||
- Contact-ready (contact is nodal)
|
||||
- Scalable (perfect parallelism)
|
||||
|
||||
---
|
||||
|
||||
## Next Steps (Prioritized)
|
||||
|
||||
### 1. Test GPU Implementation 🔄 IMMEDIATE
|
||||
```bash
|
||||
julia --project=. demos/nodal_assembly_gpu.jl
|
||||
```
|
||||
Verify results match CPU reference.
|
||||
|
||||
### 2. Add Line Search ⚠️ CRITICAL
|
||||
Current Newton solver diverges. Need backtracking line search.
|
||||
|
||||
### 3. Add Preconditioning 🎯 PERFORMANCE
|
||||
Chebyshev-Jacobi → GMG. Expert says: "THE critical factor."
|
||||
|
||||
### 4. Integrate with Newton-Krylov
|
||||
Replace element assembly with GPU nodal assembly.
|
||||
|
||||
---
|
||||
|
||||
## Performance Expectations
|
||||
|
||||
### Phase 1 (GP Data)
|
||||
- **Compute-bound** (plasticity return mapping)
|
||||
- 1M GPs: ~10-100ms on modern GPU
|
||||
- Scales linearly with GP count
|
||||
|
||||
### Phase 2 (Nodal Assembly)
|
||||
- **Memory-bound** (CSR traversal, stress reads)
|
||||
- 100K nodes: ~5-50ms on modern GPU
|
||||
- Depends on node connectivity
|
||||
|
||||
### Overall
|
||||
- Small meshes (<1K elements): GPU overhead dominates
|
||||
- Medium meshes (~10K elements): Breakeven point
|
||||
- Large meshes (100K+ elements): 10-100× speedup expected
|
||||
|
||||
---
|
||||
|
||||
## Troubleshooting
|
||||
|
||||
### GPU Kernel Doesn't Compile
|
||||
- Check CUDA.jl is installed: `using CUDA; CUDA.functional()`
|
||||
- Check Tensors.jl version compatible with CUDA.jl
|
||||
- Simplify kernel (remove plasticity, test with elastic only)
|
||||
|
||||
### Results Don't Match CPU
|
||||
- Check thread indexing (1-based in Julia!)
|
||||
- Check CSR map built correctly
|
||||
- Compare GP-by-GP (print intermediate values)
|
||||
|
||||
### Force Balance Fails
|
||||
- Check Gauss weights (should sum to element volume)
|
||||
- Check detJ computation (should be positive)
|
||||
- Check node ordering (right-hand rule)
|
||||
|
||||
---
|
||||
|
||||
## The Grand Vision
|
||||
|
||||
**Goal:** Complete GPU-resident nonlinear FEM solver
|
||||
|
||||
**Pipeline:**
|
||||
```
|
||||
Augmented Lagrangian (for contact)
|
||||
↓ Anderson acceleration HERE
|
||||
Newton Loop
|
||||
↓ Line search for globalization
|
||||
GMRES (preconditioned)
|
||||
↓ GMG preconditioner
|
||||
↓ Eisenstat-Walker forcing
|
||||
Matrix-vector product:
|
||||
↓ Phase 1: compute_gp_data_kernel!()
|
||||
↓ Phase 2: nodal_assembly_kernel!()
|
||||
|
||||
ALL ON GPU - NO CPU TRANSFERS!
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Success Criteria
|
||||
|
||||
### ✅ Achieved (November 10, 2025)
|
||||
- [x] Architecture documented
|
||||
- [x] CPU reference working
|
||||
- [x] GPU implementation complete
|
||||
- [x] Tensors.jl validated
|
||||
|
||||
### 🔄 Next Session
|
||||
- [ ] GPU kernels tested on hardware
|
||||
- [ ] Results match CPU reference
|
||||
- [ ] Force balance passes on GPU
|
||||
|
||||
### 🎯 Near-Term Goals
|
||||
- [ ] Newton solver converges
|
||||
- [ ] GMRES preconditioned
|
||||
- [ ] GPU-resident solver working
|
||||
|
||||
---
|
||||
|
||||
## Quick Reference
|
||||
|
||||
**Test CPU:**
|
||||
```bash
|
||||
julia demos/nodal_assembly_cpu.jl
|
||||
```
|
||||
|
||||
**Test GPU:**
|
||||
```bash
|
||||
julia --project=. demos/nodal_assembly_gpu.jl
|
||||
```
|
||||
|
||||
**Check architecture:**
|
||||
```bash
|
||||
cat docs/design/gpu_nodal_assembly_architecture.md
|
||||
```
|
||||
|
||||
**Check session notes:**
|
||||
```bash
|
||||
cat llm/sessions/2025-11-10_gpu_nodal_assembly_complete.md
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
**Ready to test the beast! 🚀**
|
||||
@@ -1,87 +0,0 @@
|
||||
---
|
||||
title: "GPU Architecture Quick Reference"
|
||||
date: 2025-11-10
|
||||
status: "Reference Card"
|
||||
---
|
||||
|
||||
## Design Decisions (One Page Summary)
|
||||
|
||||
### Q1: Which State Management Strategy?
|
||||
|
||||
**Answer:** Strategy 2 - Separate Mutable State (SoA)
|
||||
|
||||
**Why:** 10× better memory bandwidth (800-900 GB/s vs 50-100 GB/s)
|
||||
|
||||
### Q2: How to Eliminate Nested Newton + GMRES Loops?
|
||||
|
||||
**Answer:** Three-tier optimization
|
||||
|
||||
1. **Eisenstat-Walker** (now): Adaptive tolerance → 3× speedup
|
||||
2. **Matrix-Free NK** (Month 2): No assembly → 4× speedup
|
||||
3. **Anderson** (Month 3): Superlinear → 2.5× speedup
|
||||
|
||||
**Total: 9.8× speedup demonstrated!**
|
||||
|
||||
### Q3: How to Store Data for GPU?
|
||||
|
||||
**Answer:** Structure of Arrays (SoA) with reinterpret trick
|
||||
|
||||
```julia
|
||||
# Flat storage (GPU kernel)
|
||||
u_flat = zeros(3 * N_nodes)
|
||||
|
||||
# Physical semantics (high-level)
|
||||
u_vec3 = reinterpret(Vec{3,Float64}, u_flat)
|
||||
|
||||
# Access: u_vec3[5] returns Vec{3}
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Data Layout
|
||||
|
||||
```julia
|
||||
# Hot (mutable)
|
||||
mutable struct AssemblyState{T}
|
||||
u::Vector{T}
|
||||
material_states::Vector{State} # Flat!
|
||||
end
|
||||
|
||||
# Cold (immutable)
|
||||
struct ElementGeometry
|
||||
connectivity::Matrix{Int32}
|
||||
node_coords::Matrix{Float64}
|
||||
end
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Performance Targets
|
||||
|
||||
| Metric | Current | Target | Achieved |
|
||||
|--------|---------|--------|----------|
|
||||
| Time/iter | 8.2s | 2.1s | ✅ |
|
||||
| Memory | 12GB | 1.2GB | ✅ |
|
||||
| DOF size | 10K | 1M | ✅ |
|
||||
| Speedup | 1× | 10× | **9.8×** ✅ |
|
||||
|
||||
---
|
||||
|
||||
## Documents
|
||||
|
||||
1. `STATE_MANAGEMENT_DECISION.md` - Executive summary
|
||||
2. `GPU_ARCHITECTURE_COMPLETE.md` - Full summary
|
||||
3. `gpu_state_management.md` - Technical deep dive
|
||||
4. `matrix_free_newton_krylov.md` - Tutorial + code
|
||||
5. `reinterpret_trick.md` - Data patterns
|
||||
6. `state_implementation_roadmap.md` - Week-by-week plan
|
||||
|
||||
**Total: ~88KB documentation**
|
||||
|
||||
---
|
||||
|
||||
## Next Steps
|
||||
|
||||
**Week 1:** Create `src/assembly/state.jl`
|
||||
|
||||
**Status:** READY TO IMPLEMENT! 🚀
|
||||
@@ -1,109 +0,0 @@
|
||||
---
|
||||
title: "JuliaFEM Project Status"
|
||||
description: "Current state of the revival project as of November 2025"
|
||||
date: 2025-11-08
|
||||
updated: 2025-11-08
|
||||
author: "Jukka Aho"
|
||||
categories: ["status", "progress"]
|
||||
keywords: ["status", "progress", "revival", "roadmap"]
|
||||
audience: "contributors"
|
||||
level: "intermediate"
|
||||
type: "status report"
|
||||
series: "Contributor Manual"
|
||||
---
|
||||
|
||||
# JuliaFEM Status
|
||||
|
||||
## ✅ SUCCESS: Package Loads!
|
||||
|
||||
JuliaFEM now loads successfully on Julia 1.12.1:
|
||||
```bash
|
||||
julia> using JuliaFEM
|
||||
✓ JuliaFEM loads successfully
|
||||
Exported names: 171
|
||||
```
|
||||
|
||||
## Fixed Issues
|
||||
|
||||
### 1. Element Type Signature Errors (CRITICAL)
|
||||
**Problem:** Element type changed from `Element{Basis}` to `Element{M, Basis} where M`
|
||||
**Fixed in:**
|
||||
- `vendor/FEMBase.jl/src/FEMBase.jl` - Added AbstractBasis import
|
||||
- `vendor/FEMBase.jl/src/elements_lagrange.jl` - Fixed Poi1 subtyping
|
||||
- `vendor/FEMBeam.jl/src/beam3d.jl` - 3 function signatures
|
||||
- `vendor/MortarContact2D.jl/src/mortar2d.jl` - 2 functions
|
||||
- `vendor/MortarContact2D.jl/src/contact2d.jl` - 3 functions
|
||||
- `vendor/MortarContact2DAD.jl/src/mortar2dad.jl` - 1 function
|
||||
- `vendor/MortarContact2DAD.jl/src/contact2dad.jl` - 1 function
|
||||
- `src/problems_mortar_3d.jl` - 2 functions (M renamed to FS to avoid conflict)
|
||||
- `src/problems_contact_3d.jl` - 3 functions (M renamed to FS)
|
||||
- `src/io.jl` - 15 dispatch functions
|
||||
|
||||
### 2. Merge Conflicts (Issue #250 from 2019)
|
||||
**Fixed in:**
|
||||
- `test/runtests.jl` - Removed conflict markers
|
||||
- `src/problems_elasticity.jl` - Resolved and simplified
|
||||
|
||||
### 3. Missing Package Dependencies
|
||||
**Fixed:**
|
||||
- Created `vendor/MortarContact2DAD.jl/Project.toml`
|
||||
- Updated `Manifest.toml` to use local vendor packages
|
||||
|
||||
### 4. Parallel Assembly Code
|
||||
**Fixed:**
|
||||
- Removed references to non-existent `problem.assemble_parallel` field
|
||||
- Simplified to use non-threaded assembly (threading can be added back later)
|
||||
|
||||
## Test Status
|
||||
|
||||
**Test Suite:** 5 passed, 51 errored (but package loads!)
|
||||
|
||||
The errors are due to deeper API incompatibilities with Julia 1.12:
|
||||
- Method signature mismatches (e.g., `jacobian` function)
|
||||
- Some tests expect features from incomplete multithreading branch
|
||||
- API evolution over 6+ years (Julia 0.6 → 1.12)
|
||||
|
||||
## What Works
|
||||
|
||||
✅ Package installation and loading
|
||||
✅ All vendor packages compile
|
||||
✅ No type signature errors
|
||||
✅ Core data structures intact
|
||||
✅ 171 symbols exported
|
||||
✅ Basic FEM infrastructure present
|
||||
|
||||
## Next Steps for Full Revival
|
||||
|
||||
1. **Fix jacobian/geometry method mismatches** - Update vendor/FEMBasis for Julia 1.12
|
||||
2. **Fix remaining test errors** - Systematic fixes for API changes
|
||||
3. **Add threading infrastructure** - Properly implement parallel assembly
|
||||
4. **Update documentation** - Reflect Julia 1.12 compatibility
|
||||
5. **Benchmark performance** - Establish baseline vs old version
|
||||
|
||||
## Key Learnings
|
||||
|
||||
- Multi-package ecosystems are maintenance nightmares (see llm/TECHNICAL_VISION.md)
|
||||
- Type stability critical: Dict-based fields caused 100× slowdown
|
||||
- Git history cleanup successful: 99MB → 9.8MB (90% reduction)
|
||||
- Vendor packages approach works for development
|
||||
|
||||
## Files Modified
|
||||
|
||||
**Critical fixes (this session):**
|
||||
- 10 source files with Element type fixes
|
||||
- 2 merge conflict resolutions
|
||||
- 2 dependency files (Project.toml, Manifest.toml)
|
||||
- 1 assembly simplification
|
||||
|
||||
**Scripts created:**
|
||||
- `test.sh` - Test runner
|
||||
- `fix_src_element_types.py` - Automated type fixing
|
||||
|
||||
## Conclusion
|
||||
|
||||
**Mission accomplished:** JuliaFEM loads on modern Julia!
|
||||
|
||||
While tests have errors, the **fundamental blocker (type signatures) is resolved**.
|
||||
The package is now in a state where systematic fixing of remaining issues can proceed.
|
||||
|
||||
The 51 test errors are fixable - they're API evolution issues, not architectural problems.
|
||||
@@ -1,134 +0,0 @@
|
||||
---
|
||||
title: "Test Fixes Needed"
|
||||
description: "Known test failures and fixes required for full test suite passing"
|
||||
date: 2025-11-09
|
||||
updated: 2025-11-09
|
||||
author: "Jukka Aho"
|
||||
categories: ["testing", "todo"]
|
||||
keywords: ["tests", "failures", "fixes", "todo"]
|
||||
audience: "contributors"
|
||||
level: "intermediate"
|
||||
type: "technical note"
|
||||
series: "Contributor Manual"
|
||||
status: "active maintenance"
|
||||
---
|
||||
|
||||
**Date:** November 8, 2025
|
||||
**Status:** 5 passing, 49 failing (infrastructure now in place)
|
||||
|
||||
## Summary
|
||||
|
||||
Tests are failing due to API evolution between Julia 0.6/1.0 (2018) and Julia 1.12 (2025), not fundamental architectural problems. Package loads successfully and core functionality works.
|
||||
|
||||
## Main Issues
|
||||
|
||||
### 1. Missing `aster_read_mesh` (14 tests)
|
||||
|
||||
**Problem:** Tests use `aster_read_mesh()` from IO submodule, but it requires HDF5
|
||||
**Files affected:** Most 3D elasticity tests, med file tests
|
||||
**Fix options:**
|
||||
|
||||
- A) Add HDF5 as optional dependency (Julia 1.9+ package extensions)
|
||||
- B) Skip tests that need .med files for now
|
||||
- C) Convert test meshes to .inp format (ABAQUS, which we support)
|
||||
|
||||
**Recommendation:** Option C - convert test meshes to .inp format
|
||||
|
||||
### 2. `eval_basis!` Signature Mismatch (2 tests)
|
||||
|
||||
**Problem:** `eval_basis!(::Type{Seg2}, ::Matrix, ::Tuple{Float64})`
|
||||
**Current:** `eval_basis!(::Seg2, ::Vector, ::Tuple{Float64}, time::Float64)`
|
||||
**Location:** `vendor/FEMBasis.jl`
|
||||
|
||||
**Fix:** Update signature in FEMBasis or fix call sites
|
||||
|
||||
### 3. `jacobian` Signature Mismatch (~20 tests)
|
||||
|
||||
**Problem:** Tests call `jacobian(element_type, X, xi)` with old signatures
|
||||
**Current API:** Different parameter order or types
|
||||
**Location:** `vendor/FEMBasis.jl/src/jacobian.jl`
|
||||
|
||||
**Fix:** Consolidate FEMBasis into src/basis/ with modern API
|
||||
|
||||
### 4. `allocate_buffer` Missing (2 tests)
|
||||
|
||||
**Problem:** `allocate_buffer(::Problem{Elasticity}, ::Vector{Element})`
|
||||
**Status:** Method doesn't exist in current codebase
|
||||
**Fix:** Either restore method or update tests to not need it
|
||||
|
||||
### 5. `Analysis` Missing (5 tests) - ✅ FIXED
|
||||
|
||||
**Status:** Now exported, these tests should pass
|
||||
|
||||
### 6. Statistics Package Missing (1 test) - ✅ FIXED
|
||||
|
||||
**Status:** Now in test dependencies
|
||||
|
||||
## Test Categories
|
||||
|
||||
### ✅ Passing (5 tests)
|
||||
|
||||
- Virtual work test
|
||||
- Contact 2D/3D tests
|
||||
- Mortar 2D tests
|
||||
- Heat transfer (basic)
|
||||
|
||||
### ❌ Failing - Missing HDF5 (~14 tests)
|
||||
|
||||
- test_elasticity_2d_nonlinear_with_surface_load.jl
|
||||
- test_elasticity_3d_unit_block.jl
|
||||
- test_elasticity_med_pyr5_point_load.jl
|
||||
- test_elasticity_plane_strain.jl
|
||||
- test_elasticity_pyr5_point_load.jl
|
||||
- Many more...
|
||||
|
||||
### ❌ Failing - API Mismatches (~30 tests)
|
||||
|
||||
- eval_basis! signature (2)
|
||||
- jacobian signature (~20)
|
||||
- allocate_buffer missing (2)
|
||||
- Various others (6)
|
||||
|
||||
## Action Plan
|
||||
|
||||
### Phase 1: Low-Hanging Fruit (1-2 hours)
|
||||
|
||||
1. ✅ Export Analysis types
|
||||
2. ✅ Add Statistics to test deps
|
||||
3. ⏳ Skip/comment out HDF5-dependent tests temporarily
|
||||
4. ⏳ Re-run tests, see how many pass
|
||||
|
||||
### Phase 2: API Fixes (4-6 hours)
|
||||
|
||||
1. Fix `eval_basis!` signature in FEMBasis
|
||||
2. Fix `jacobian` signature in FEMBasis
|
||||
3. Either restore `allocate_buffer` or update tests
|
||||
4. Fix any remaining signature mismatches
|
||||
|
||||
### Phase 3: Mesh Conversion (2-4 hours)
|
||||
|
||||
1. Find all .med test meshes
|
||||
2. Convert to .inp format using Code Aster or similar
|
||||
3. Update test files to use .inp instead of .med
|
||||
4. Re-run tests
|
||||
|
||||
### Phase 4: Verify All Pass (1 hour)
|
||||
|
||||
1. Run full test suite
|
||||
2. Fix any remaining issues
|
||||
3. Update CI to run tests automatically
|
||||
4. Celebrate! 🎉
|
||||
|
||||
## Expected Outcome
|
||||
|
||||
After these fixes:
|
||||
|
||||
- ~40+ tests should pass (out of 56 total)
|
||||
- CI will catch regressions automatically
|
||||
- Good foundation for further consolidation work
|
||||
|
||||
## Notes
|
||||
|
||||
The fact that package loads and 5 tests pass is actually very good news - it means the core architecture is sound. These are just API compatibility issues that accumulated over 6 years of Julia evolution.
|
||||
|
||||
Most fixes are mechanical (update signatures) rather than requiring deep understanding of the algorithms.
|
||||
@@ -1,802 +0,0 @@
|
||||
---
|
||||
title: "Testing Philosophy"
|
||||
description: "How and why we test in JuliaFEM"
|
||||
date: 2025-11-08
|
||||
author: "Jukka Aho"
|
||||
categories: ["testing", "quality assurance"]
|
||||
keywords: ["testing", "unit tests", "verification", "validation"]
|
||||
audience: "contributors"
|
||||
level: "intermediate"
|
||||
type: "guide"
|
||||
series: "Contributor Manual"
|
||||
---
|
||||
|
||||
# Testing Philosophy
|
||||
|
||||
**Date:** November 9, 2025
|
||||
**Goal:** 99% code coverage with educational, fast, well-structured tests
|
||||
**Tool:** Literate.jl for test-as-documentation
|
||||
|
||||
---
|
||||
|
||||
## Core Principles
|
||||
|
||||
### 1. Tests as Teaching Material
|
||||
|
||||
**Every test should teach something.**
|
||||
|
||||
Tests are not just validation - they're the **primary way users learn JuliaFEM**. When someone asks "How do I solve an elasticity problem?", the answer should be: "Look at `test/tutorials/elasticity_basics.jl`"
|
||||
|
||||
**Benefits:**
|
||||
|
||||
- Users learn by example (better than API docs)
|
||||
- Tests stay current (if API changes, tests must update)
|
||||
- Documentation never lies (it's tested code!)
|
||||
- Newcomers can contribute tests (learning exercise)
|
||||
|
||||
### 2. Literate.jl for Test-Driven Documentation
|
||||
|
||||
Use Literate.jl to write tests as narrative documents:
|
||||
|
||||
```julia
|
||||
# # Solving Your First Elasticity Problem
|
||||
#
|
||||
# This tutorial shows how to solve a simple 2D elasticity problem.
|
||||
# We'll create a square block, apply boundary conditions, and solve.
|
||||
|
||||
using JuliaFEM
|
||||
using Test
|
||||
|
||||
# ## Step 1: Create the Geometry
|
||||
#
|
||||
# First, define the nodes of a unit square:
|
||||
|
||||
X = Dict(
|
||||
1 => [0.0, 0.0],
|
||||
2 => [1.0, 0.0],
|
||||
3 => [1.0, 1.0],
|
||||
4 => [0.0, 1.0]
|
||||
)
|
||||
|
||||
# ## Step 2: Create Elements
|
||||
#
|
||||
# Create a single Quad4 element:
|
||||
|
||||
element = Element(Quad4, (1, 2, 3, 4))
|
||||
update!(element, "geometry", X)
|
||||
|
||||
# ... and so on
|
||||
```
|
||||
|
||||
**Output:** Same file generates both test (runs in CI) and documentation (builds HTML).
|
||||
|
||||
### 3. Structured Test Hierarchy
|
||||
|
||||
Organize tests to match learning progression:
|
||||
|
||||
```text
|
||||
test/
|
||||
├── tutorials/ # Literate.jl files (test + docs)
|
||||
│ ├── 01_fundamentals/
|
||||
│ │ ├── creating_elements.jl
|
||||
│ │ ├── fields_and_updates.jl
|
||||
│ │ ├── reading_meshes.jl
|
||||
│ │ └── basis_functions.jl
|
||||
│ ├── 02_linear_problems/
|
||||
│ │ ├── elasticity_1d.jl
|
||||
│ │ ├── elasticity_2d.jl
|
||||
│ │ ├── heat_transfer.jl
|
||||
│ │ └── boundary_conditions.jl
|
||||
│ ├── 03_nonlinear_problems/
|
||||
│ │ ├── large_deformation.jl
|
||||
│ │ ├── plasticity.jl
|
||||
│ │ └── contact_basics.jl
|
||||
│ ├── 04_advanced/
|
||||
│ │ ├── contact_2d.jl
|
||||
│ │ ├── contact_3d.jl
|
||||
│ │ ├── mortar_methods.jl
|
||||
│ │ └── friction.jl
|
||||
│ └── 05_parallel/
|
||||
│ ├── threading.jl
|
||||
│ ├── gpu_assembly.jl
|
||||
│ └── distributed.jl
|
||||
├── unit/ # Fast unit tests (not Literate)
|
||||
│ ├── basis/
|
||||
│ ├── assembly/
|
||||
│ └── solvers/
|
||||
├── verification/ # Known analytical solutions
|
||||
│ ├── timoshenko_beam.jl
|
||||
│ ├── hertz_contact.jl
|
||||
│ └── cook_membrane.jl
|
||||
└── runtests.jl # Test runner
|
||||
```
|
||||
|
||||
### 4. Fast Tests First
|
||||
|
||||
**Test pyramid:**
|
||||
|
||||
```text
|
||||
/\
|
||||
/ \ Integration tests (slow, few)
|
||||
/____\
|
||||
/ \ Tutorial tests (medium, some)
|
||||
/________\
|
||||
/ \ Unit tests (fast, many)
|
||||
/__________\
|
||||
```
|
||||
|
||||
**Timing targets:**
|
||||
|
||||
- Unit tests: < 5 minutes (run during development)
|
||||
- Tutorial tests: < 15 minutes (run before commits)
|
||||
- Full suite: < 30 minutes (run in CI)
|
||||
|
||||
**Strategy:**
|
||||
|
||||
- Use realistic meshes in tutorials (~10 elements: not too trivial, not too slow)
|
||||
- Test correctness with analytical solutions, not big problems
|
||||
- Profile and optimize slow tests
|
||||
- Mark slow tests with `@testset "slow: hertz_contact"` (skip during dev)
|
||||
- Co-locate mesh files with tests (no shared `/test/meshes/` directory)
|
||||
- Include mesh generation recipes (show how mesh was created for reproducibility)
|
||||
|
||||
### 5. Coverage-Driven Development
|
||||
|
||||
**Target: 99% code coverage**
|
||||
|
||||
Every function should have:
|
||||
|
||||
1. **Happy path test** - normal usage
|
||||
2. **Edge case tests** - empty input, single element, etc.
|
||||
3. **Error tests** - what happens with bad input?
|
||||
|
||||
**Process:**
|
||||
|
||||
1. Write tutorial (covers main API)
|
||||
2. Check coverage report
|
||||
3. Add unit tests for uncovered lines
|
||||
4. Repeat until 99%+
|
||||
|
||||
**Tools:**
|
||||
|
||||
- Coverage.jl (built into Julia)
|
||||
- LocalCoverage.jl (for local checks)
|
||||
- Codecov (in CI, we set this up yesterday)
|
||||
|
||||
---
|
||||
|
||||
## Test Structure Details
|
||||
|
||||
### Tutorial Tests (Literate.jl)
|
||||
|
||||
**Template:**
|
||||
|
||||
```julia
|
||||
# # Tutorial Title
|
||||
#
|
||||
# Brief description of what this tutorial teaches.
|
||||
# Prerequisites: what the reader should know first.
|
||||
|
||||
using JuliaFEM
|
||||
using Test
|
||||
|
||||
# ## Section 1: Concept Explanation
|
||||
#
|
||||
# Explain the concept in prose, with equations if needed:
|
||||
#
|
||||
# The strain-displacement relationship is:
|
||||
# ```math
|
||||
# ε = \frac{1}{2}(∇u + ∇u^T)
|
||||
# ```
|
||||
|
||||
# Code demonstrating the concept
|
||||
element = Element(Quad4, (1,2,3,4))
|
||||
|
||||
# ## Section 2: Building Up
|
||||
#
|
||||
# Step-by-step construction of a working example
|
||||
|
||||
# ... code ...
|
||||
|
||||
# ## Section 3: Validation
|
||||
#
|
||||
# Test that the result is correct (this is still a test!)
|
||||
|
||||
@testset "Elasticity 2D" begin
|
||||
@test isapprox(u_computed, u_analytical, rtol=1e-6)
|
||||
end
|
||||
|
||||
# ## Discussion
|
||||
#
|
||||
# What did we learn? What can we do next?
|
||||
# Links to related tutorials.
|
||||
```
|
||||
|
||||
**Generation:**
|
||||
|
||||
```julia
|
||||
using Literate
|
||||
Literate.markdown("test/tutorials/01_fundamentals/creating_elements.jl",
|
||||
"docs/src/tutorials/")
|
||||
Literate.notebook("test/tutorials/01_fundamentals/creating_elements.jl",
|
||||
"docs/notebooks/")
|
||||
```
|
||||
|
||||
### Unit Tests (Fast Validation)
|
||||
|
||||
**Purpose:** Test individual functions in isolation
|
||||
**Style:** Concise, no narrative
|
||||
**Location:** `test/unit/`
|
||||
|
||||
```julia
|
||||
@testset "Basis Functions - Quad4" begin
|
||||
@testset "Evaluation at ξ=0, η=0" begin
|
||||
N = eval_basis(Quad4, (0.0, 0.0))
|
||||
@test N ≈ [0.25, 0.25, 0.25, 0.25]
|
||||
end
|
||||
|
||||
@testset "Derivatives" begin
|
||||
dN = eval_dbasis(Quad4, (0.0, 0.0))
|
||||
@test size(dN) == (2, 4)
|
||||
end
|
||||
|
||||
@testset "Edge case: extreme ξ" begin
|
||||
N = eval_basis(Quad4, (1.0, 1.0))
|
||||
@test N[3] ≈ 1.0
|
||||
@test sum(N) ≈ 1.0
|
||||
end
|
||||
end
|
||||
```
|
||||
|
||||
### Verification Tests (Known Solutions)
|
||||
|
||||
**Purpose:** Validate against analytical solutions or published results
|
||||
**Style:** Brief explanation + reference
|
||||
**Location:** `test/verification/`
|
||||
|
||||
```julia
|
||||
# Timoshenko Beam - Verification Test
|
||||
#
|
||||
# Reference: Timoshenko & Goodier, "Theory of Elasticity", 3rd Ed.
|
||||
# Problem: Cantilever beam with end load
|
||||
# Analytical solution available for tip displacement
|
||||
|
||||
using JuliaFEM, Test
|
||||
|
||||
# Problem parameters from reference
|
||||
L = 10.0 # Length
|
||||
h = 1.0 # Height
|
||||
E = 200e3 # Young's modulus
|
||||
ν = 0.3 # Poisson's ratio
|
||||
P = 100.0 # End load
|
||||
|
||||
# ... setup and solve ...
|
||||
|
||||
# Analytical solution
|
||||
u_tip_analytical = P*L^3 / (3*E*I)
|
||||
|
||||
@testset "Timoshenko Beam Verification" begin
|
||||
@test isapprox(u_tip_computed, u_tip_analytical, rtol=0.01)
|
||||
end
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Implementation Roadmap
|
||||
|
||||
### Phase 1: Infrastructure (Week 1) ✅ IN PROGRESS
|
||||
|
||||
**Goal:** Set up Literate.jl integration and test structure
|
||||
|
||||
**Tasks:**
|
||||
|
||||
1. ✅ Create `test/tutorials/` directory structure
|
||||
2. ✅ Create new test runner (`test/runtests_new.jl`) with environment control
|
||||
3. ✅ Add Gmsh.jl to test dependencies
|
||||
4. ✅ Tutorial 1: Creating elements (5 tests passing)
|
||||
5. ✅ Tutorial 2: Reading Gmsh meshes (72 tests passing, with recipe and co-located .msh)
|
||||
6. ⏳ Tutorial 4: 1-element validation (priority for Issue #265)
|
||||
7. ⏳ Tutorial 3: Basis functions
|
||||
8. ⏳ Add Literate.jl to `docs/Project.toml`
|
||||
9. ⏳ Update `docs/make.jl` to process tutorials
|
||||
10. ⏳ Add coverage tools (LocalCoverage.jl)
|
||||
|
||||
**Status:** 77/77 tests passing (Tutorial 1: 5, Tutorial 2: 72)
|
||||
|
||||
**Deliverable:** Running CI that executes tutorials and reports coverage
|
||||
|
||||
### Phase 2: Core Tutorials (Week 2-3)
|
||||
|
||||
**Goal:** Write 10-15 fundamental tutorials covering main API
|
||||
|
||||
**Priority order:**
|
||||
|
||||
1. ✅ Creating elements and updating fields ⭐ (Done: 5 tests passing)
|
||||
2. ✅ Reading meshes (Gmsh .msh format) ⭐ (Done: 72 tests passing)
|
||||
3. **1-element validation tests** ⭐⭐ (Next: helps others validate FEM software, see Issue #265)
|
||||
4. Basis functions and integration ⭐
|
||||
5. Simple 2D elasticity (10 elements) ⭐
|
||||
6. Boundary conditions (Dirichlet)
|
||||
7. Surface loads and tractions
|
||||
8. Heat transfer basics
|
||||
9. Assembly process
|
||||
10. Solving linear systems
|
||||
|
||||
**Note on mesh format:** We use **Gmsh** (via Gmsh.jl) instead of ABAQUS because:
|
||||
|
||||
- No license required (accessible to everyone)
|
||||
- Julia native package (Gmsh.jl)
|
||||
- Modern, actively developed
|
||||
- Programmatic mesh generation (reproducible)
|
||||
- Co-located mesh files with test files (self-contained tests)
|
||||
|
||||
**Success metric:** 60%+ code coverage from tutorials alone
|
||||
|
||||
### Phase 3: Advanced Tutorials (Week 4-5)
|
||||
|
||||
**Goal:** Cover advanced features
|
||||
|
||||
**Topics:**
|
||||
|
||||
1. Nonlinear elasticity (large deformation)
|
||||
2. Contact mechanics basics (2D)
|
||||
3. Mortar methods
|
||||
4. 3D contact
|
||||
5. Material models
|
||||
|
||||
**Success metric:** 80%+ coverage
|
||||
|
||||
### Phase 4: Unit Tests (Week 6)
|
||||
|
||||
**Goal:** Fill coverage gaps with fast unit tests
|
||||
|
||||
**Process:**
|
||||
|
||||
1. Generate coverage report: `julia --code-coverage=user test/runtests.jl`
|
||||
2. Analyze: `using Coverage; LCOV.writefile("coverage.info", process_folder())`
|
||||
3. Find uncovered lines
|
||||
4. Write unit tests for each uncovered function
|
||||
5. Repeat until 99%+
|
||||
|
||||
**Success metric:** 99% coverage, < 5 min unit test runtime
|
||||
|
||||
### Phase 5: Verification (Week 7)
|
||||
|
||||
**Goal:** Validate against known solutions
|
||||
|
||||
**Tests:**
|
||||
|
||||
1. Timoshenko beam (bending)
|
||||
2. Hertz contact (2D)
|
||||
3. Cook's membrane (stress concentration)
|
||||
4. Patch tests (element validation)
|
||||
5. Manufactured solutions
|
||||
|
||||
**Success metric:** All verification tests pass with < 1% error
|
||||
|
||||
### Phase 6: Documentation (Week 8)
|
||||
|
||||
**Goal:** Polish and publish
|
||||
|
||||
**Tasks:**
|
||||
|
||||
1. Review all tutorial narratives
|
||||
2. Add figures and visualizations
|
||||
3. Cross-link tutorials
|
||||
4. Build documentation locally
|
||||
5. Deploy to GitHub Pages
|
||||
6. Write README.md guide to tutorials
|
||||
|
||||
**Success metric:** Beautiful, usable documentation website
|
||||
|
||||
---
|
||||
|
||||
## Test Runner Design
|
||||
|
||||
### `test/runtests.jl`
|
||||
|
||||
```julia
|
||||
using Test, JuliaFEM
|
||||
|
||||
# Determine what to run based on environment
|
||||
const RUN_UNIT = get(ENV, "JULIAFEM_TEST_UNIT", "true") == "true"
|
||||
const RUN_TUTORIALS = get(ENV, "JULIAFEM_TEST_TUTORIALS", "true") == "true"
|
||||
const RUN_SLOW = get(ENV, "JULIAFEM_TEST_SLOW", "false") == "true"
|
||||
const RUN_VERIFICATION = get(ENV, "JULIAFEM_TEST_VERIFICATION", "true") == "true"
|
||||
|
||||
# Fast unit tests (always run)
|
||||
if RUN_UNIT
|
||||
@testset "Unit Tests" begin
|
||||
include("unit/basis/test_quad4.jl")
|
||||
include("unit/basis/test_seg2.jl")
|
||||
# ... more unit tests
|
||||
end
|
||||
end
|
||||
|
||||
# Tutorial tests (run in CI, optional locally)
|
||||
if RUN_TUTORIALS
|
||||
@testset "Tutorials" begin
|
||||
# These are also documentation!
|
||||
include("tutorials/01_fundamentals/creating_elements.jl")
|
||||
include("tutorials/01_fundamentals/reading_meshes.jl")
|
||||
include("tutorials/02_linear_problems/elasticity_2d.jl")
|
||||
# ... more tutorials
|
||||
end
|
||||
end
|
||||
|
||||
# Slow integration tests (CI only by default)
|
||||
if RUN_SLOW
|
||||
@testset "Slow Tests" begin
|
||||
include("verification/hertz_contact.jl")
|
||||
# Large mesh tests
|
||||
end
|
||||
end
|
||||
|
||||
# Verification tests
|
||||
if RUN_VERIFICATION
|
||||
@testset "Verification" begin
|
||||
include("verification/timoshenko_beam.jl")
|
||||
include("verification/cook_membrane.jl")
|
||||
end
|
||||
end
|
||||
```
|
||||
|
||||
**Usage:**
|
||||
|
||||
```bash
|
||||
# During development (fast, < 5 min)
|
||||
julia --project=. -e 'using Pkg; Pkg.test()'
|
||||
|
||||
# Before commit (< 15 min)
|
||||
JULIAFEM_TEST_TUTORIALS=true julia --project=. test/runtests.jl
|
||||
|
||||
# Full CI run (< 30 min)
|
||||
JULIAFEM_TEST_SLOW=true julia --project=. test/runtests.jl
|
||||
|
||||
# Only unit tests (< 2 min)
|
||||
JULIAFEM_TEST_TUTORIALS=false julia --project=. test/runtests.jl
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Coverage Workflow
|
||||
|
||||
### Local Development
|
||||
|
||||
```bash
|
||||
# 1. Run tests with coverage
|
||||
julia --project=. --code-coverage=user test/runtests.jl
|
||||
|
||||
# 2. Generate coverage report
|
||||
julia --project=. -e '
|
||||
using Coverage
|
||||
coverage = process_folder()
|
||||
covered = length(filter(c -> c.coverage > 0, coverage))
|
||||
total = length(coverage)
|
||||
println("Coverage: $(round(100*covered/total, digits=2))%")
|
||||
LCOV.writefile("coverage.info", coverage)
|
||||
'
|
||||
|
||||
# 3. View in browser (requires genhtml from lcov package)
|
||||
genhtml coverage.info -o coverage/
|
||||
firefox coverage/index.html
|
||||
```
|
||||
|
||||
### GitHub Actions CI
|
||||
|
||||
Already set up yesterday! Codecov will:
|
||||
|
||||
- Track coverage over time
|
||||
- Comment on PRs with coverage changes
|
||||
- Show which lines are uncovered
|
||||
- Badge in README.md
|
||||
|
||||
---
|
||||
|
||||
## Writing Style Guide
|
||||
|
||||
### For Tutorials (Literate.jl)
|
||||
|
||||
**Do:**
|
||||
|
||||
- ✅ Explain **why**, not just **what**
|
||||
- ✅ Use realistic examples (~10 elements, not too simple or too complex)
|
||||
- ✅ Include mathematical notation where helpful
|
||||
- ✅ Show output/results
|
||||
- ✅ Link to related tutorials
|
||||
- ✅ Test the actual result (still a test!)
|
||||
- ✅ Co-locate mesh files with test files (self-contained)
|
||||
- ✅ Include mesh generation recipe (show how mesh was created)
|
||||
|
||||
**Don't:**
|
||||
|
||||
- ❌ Assume prior knowledge (explain or link)
|
||||
- ❌ Use large meshes (slow tests)
|
||||
- ❌ Skip explanation (this is documentation!)
|
||||
- ❌ Test implementation details (test behavior)
|
||||
|
||||
### For Unit Tests
|
||||
|
||||
**Do:**
|
||||
|
||||
- ✅ Test one thing per `@testset`
|
||||
- ✅ Use descriptive test names
|
||||
- ✅ Test edge cases (empty, single, many)
|
||||
- ✅ Test error conditions (`@test_throws`)
|
||||
- ✅ Be concise (no prose needed)
|
||||
|
||||
**Don't:**
|
||||
|
||||
- ❌ Mix multiple concepts in one test
|
||||
- ❌ Use real-world complex examples
|
||||
- ❌ Duplicate tutorial content
|
||||
|
||||
---
|
||||
|
||||
## Migration Plan for Existing Tests
|
||||
|
||||
### Current State (56 test files)
|
||||
|
||||
Many old tests are:
|
||||
|
||||
- Poorly documented
|
||||
- Using outdated API
|
||||
- Slow (large meshes)
|
||||
- Not structured for learning
|
||||
|
||||
### Migration Strategy
|
||||
|
||||
**Don't delete old tests immediately!** Instead:
|
||||
|
||||
1. **Categorize:** Is it tutorial material, unit test, or verification?
|
||||
2. **Rewrite:** Create new version following philosophy
|
||||
3. **Verify:** Ensure new test covers same functionality
|
||||
4. **Archive:** Move old test to `test/archive/` with note
|
||||
5. **Delete:** After new tests run successfully in CI
|
||||
|
||||
**Example:**
|
||||
|
||||
```text
|
||||
test/test_elasticity_2d_linear_with_surface_load.jl (old)
|
||||
→ test/tutorials/02_linear_problems/elasticity_2d.jl (new, Literate)
|
||||
→ test/archive/test_elasticity_2d_linear_with_surface_load.jl.old (keep for reference)
|
||||
→ delete after 1 month if no issues
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Success Metrics
|
||||
|
||||
### Quantitative
|
||||
|
||||
- **Coverage:** 99%+ by end of Phase 4
|
||||
- **Speed:** < 5 min unit tests, < 30 min full suite
|
||||
- **Count:** 15+ tutorials, 100+ unit tests
|
||||
|
||||
### Qualitative
|
||||
|
||||
- **Can a new user learn JuliaFEM from tutorials alone?** (Ask someone!)
|
||||
- **Are tutorials referenced in issue discussions?** ("See tutorial X")
|
||||
- **Do contributors write tests first?** (TDD culture)
|
||||
- **Is documentation always up-to-date?** (Literate.jl ensures it)
|
||||
|
||||
---
|
||||
|
||||
## Special: 1-Element Validation Tests
|
||||
|
||||
### Motivation (Issue #265)
|
||||
|
||||
In 2019, JuliaFEM was used to **validate another FEM software**. A user computed a reference solution with JuliaFEM and compared it against their own implementation. This is a powerful use case we should embrace!
|
||||
|
||||
### Design Philosophy
|
||||
|
||||
**Goal:** Create 1-element tests that can be:
|
||||
|
||||
1. **Hand-calculated** - Simple enough to verify by hand
|
||||
2. **Exact** - Integer or simple fractional results (no floating-point ambiguity)
|
||||
3. **Reference-quality** - Others can use to validate their FEM code
|
||||
4. **Educational** - Show the math, not just code
|
||||
|
||||
### Example Pattern
|
||||
|
||||
```julia
|
||||
# # 1-Element Validation: Quad4 Elasticity
|
||||
#
|
||||
# This tutorial computes the stiffness matrix for a single Quad4 element
|
||||
# under plane stress conditions. The solution can be verified by hand.
|
||||
#
|
||||
# **Use case:** Validating FEM implementations (see Issue #265)
|
||||
#
|
||||
# ## Problem Setup
|
||||
#
|
||||
# Geometry: Unit square [0,1] × [0,1]
|
||||
# Material: E = 100.0, ν = 0.3 (simple values)
|
||||
# Element: Single Quad4 with nodes at corners
|
||||
#
|
||||
# ## Hand Calculation
|
||||
#
|
||||
# For Quad4 under plane stress, the stiffness matrix has structure:
|
||||
# ... detailed derivation ...
|
||||
#
|
||||
# Expected K[1,1] = ... (show calculation)
|
||||
|
||||
using JuliaFEM, Test
|
||||
|
||||
# Define nodes (unit square)
|
||||
nodes = Dict(
|
||||
1 => [0.0, 0.0],
|
||||
2 => [1.0, 0.0],
|
||||
3 => [1.0, 1.0],
|
||||
4 => [0.0, 1.0]
|
||||
)
|
||||
|
||||
# Create element
|
||||
element = Element(Quad4, (1, 2, 3, 4))
|
||||
update!(element, "geometry", nodes)
|
||||
update!(element, "youngs modulus", 100.0)
|
||||
update!(element, "poissons ratio", 0.3)
|
||||
|
||||
# Assemble stiffness matrix
|
||||
K = assemble_stiffness_matrix(element)
|
||||
|
||||
# Validate against hand calculation
|
||||
@testset "1-Element Validation: Quad4 Stiffness" begin
|
||||
# Check specific entries that we calculated by hand
|
||||
@test K[1,1] ≈ 42.5 # Hand calculated value
|
||||
@test K[1,2] ≈ 10.0 # Hand calculated value
|
||||
# ... more validations
|
||||
|
||||
# Symmetry check
|
||||
@test issymmetric(K)
|
||||
|
||||
# Positive definite check (all eigenvalues > 0, after BC)
|
||||
# ... (need to apply constraints first)
|
||||
end
|
||||
```
|
||||
|
||||
### Benefits
|
||||
|
||||
1. **Validation tool** - Others can use JuliaFEM as reference
|
||||
2. **Debugging aid** - If basic case fails, know where to look
|
||||
3. **Educational** - Shows the math behind FEM
|
||||
4. **Confidence** - Proves implementation is correct
|
||||
5. **Regression test** - Any refactoring must pass these
|
||||
|
||||
### Implementation Priority
|
||||
|
||||
**High priority** - These tests are foundational. Should be in Tutorial 4 (before basis functions).
|
||||
|
||||
**Coverage:** Create 1-element tests for each element type:
|
||||
|
||||
- Seg2 (1D bar)
|
||||
- Tri3 (2D triangle)
|
||||
- Quad4 (2D quadrilateral)
|
||||
- Tet4 (3D tetrahedron)
|
||||
- Hex8 (3D hexahedron)
|
||||
|
||||
---
|
||||
|
||||
## Questions to Resolve
|
||||
|
||||
### 1. Literate.jl Execution Strategy
|
||||
|
||||
**Option A:** Tutorials are in `test/tutorials/`, executed by `Pkg.test()`
|
||||
|
||||
- Pro: Ensures tutorials always work
|
||||
- Con: Slows down test suite
|
||||
|
||||
**Option B:** Tutorials are in `docs/tutorials/`, executed by docs build
|
||||
|
||||
- Pro: Fast test suite
|
||||
- Con: Tutorials might break without noticing
|
||||
|
||||
**Recommendation:** Option A (tutorials in test/), with `RUN_TUTORIALS` env var
|
||||
|
||||
### 2. Visualization in Tutorials
|
||||
|
||||
Should tutorials include plots/visualizations?
|
||||
|
||||
**Option A:** Yes, using Makie.jl or similar
|
||||
|
||||
- Pro: Very educational, shows results
|
||||
- Con: Heavy dependency, slow to precompile
|
||||
|
||||
**Option B:** No plots in tests, but show code in docs
|
||||
|
||||
- Pro: Fast tests
|
||||
- Con: Less visual learning
|
||||
|
||||
**Recommendation:** Option B initially, add plots later as optional
|
||||
|
||||
### 3. Mesh Files ✅ DECIDED
|
||||
|
||||
**Decision:** Use Gmsh.jl for mesh generation (Option B)
|
||||
|
||||
**Rationale:**
|
||||
|
||||
- No license required (ABAQUS needs license, not accessible)
|
||||
- No heavy dependencies (HDF5/MED format requires HDF5.jl)
|
||||
- Programmatic generation (reproducible, can show recipe)
|
||||
- Modern and well-maintained
|
||||
- Julia native integration
|
||||
|
||||
**Implementation Pattern:**
|
||||
|
||||
1. Create `*_recipe.jl` script showing how mesh was generated
|
||||
2. Run recipe to generate `*.msh` file
|
||||
3. Commit both recipe and .msh file to repository
|
||||
4. Test reads the .msh file (fast, reliable)
|
||||
5. Users can see recipe to understand mesh structure
|
||||
|
||||
**Example:**
|
||||
|
||||
```text
|
||||
test/tutorials/01_fundamentals/
|
||||
├── reading_gmsh_meshes.jl # The actual test
|
||||
├── reading_gmsh_meshes.msh # Pre-generated mesh (committed)
|
||||
└── reading_gmsh_meshes_recipe.jl # Shows how mesh was created
|
||||
```
|
||||
|
||||
**Benefits:**
|
||||
|
||||
- Self-contained tests (mesh next to test file)
|
||||
- Reproducible (recipe shows exactly how to recreate)
|
||||
- Fast (don't generate mesh in CI, just read it)
|
||||
- Educational (recipe teaches mesh generation)
|
||||
- Accessible (no external software needed to run tests)
|
||||
|
||||
---
|
||||
|
||||
## Next Steps (This Week)
|
||||
|
||||
### Immediate (November 9, 2025)
|
||||
|
||||
1. ✅ Create this document (done!)
|
||||
2. ✅ Create tutorial directory structure
|
||||
3. ✅ Write first tutorial (creating elements - 5 tests passing)
|
||||
4. ✅ Write second tutorial (reading Gmsh meshes - 72 tests passing)
|
||||
5. ✅ Create new test runner (`test/runtests_new.jl`)
|
||||
6. 🔄 **Now: Tutorial 4 - 1-element validation** (priority for Issue #265)
|
||||
7. ⏳ Tutorial 3 - Basis functions
|
||||
8. ⏳ Add Literate.jl to docs dependencies
|
||||
9. ⏳ Update `docs/make.jl` for HTML generation
|
||||
|
||||
### This Week (Week 1 - November 9-15, 2025)
|
||||
|
||||
**Completed:**
|
||||
|
||||
1. ✅ Tutorial 1: Creating elements and fields (5 tests)
|
||||
2. ✅ Tutorial 2: Reading Gmsh meshes (72 tests)
|
||||
|
||||
**In Progress:**
|
||||
3. 🔄 Tutorial 4: 1-element validation (priority - helps validate other FEM software)
|
||||
4. 🔄 Tutorial 3: Basis functions and integration
|
||||
|
||||
**Planned:**
|
||||
5. Tutorial 5: Simple 2D elasticity (10 elements, use Tutorial 2 mesh)
|
||||
6. Set up coverage workflow locally
|
||||
7. Test Literate.jl → HTML generation
|
||||
8. Update CI to run new test structure
|
||||
|
||||
**Target:** 5 tutorials by end of week (currently 2/5)
|
||||
|
||||
### Next Week
|
||||
|
||||
Continue writing tutorials, aim for 10 total by end of week.
|
||||
|
||||
---
|
||||
|
||||
## Conclusion
|
||||
|
||||
**This is a complete overhaul of testing strategy** - from "fix 49 failing tests" to "build educational test suite that reaches 99% coverage."
|
||||
|
||||
**Timeline:** ~8 weeks for full implementation
|
||||
**Effort:** Significant, but creates lasting value
|
||||
**Benefit:** Tests become documentation, documentation is always tested
|
||||
|
||||
**Philosophy:** Tests are not a chore - they're the best way to teach users how to use JuliaFEM.
|
||||
|
||||
---
|
||||
|
||||
**Status:** 📋 Planning complete, ready to implement
|
||||
**Next:** Get your feedback, then start Phase 1
|
||||
@@ -1,64 +0,0 @@
|
||||
# Simple usage examples
|
||||
|
||||
A simple example demonstrating the basic usage of package. Calculate a simple
|
||||
one element model. Add pressure load on top and support block symmetrically.
|
||||
|
||||
```@example 1
|
||||
using JuliaFEM # hide
|
||||
X = Dict(
|
||||
1 => [0.0, 0.0],
|
||||
2 => [1.0, 0.0],
|
||||
3 => [1.0, 1.0],
|
||||
4 => [0.0, 1.0])
|
||||
```
|
||||
|
||||
```@example 1
|
||||
element = Element(Quad4, [1, 2, 3, 4])
|
||||
update!(element, "geometry", X)
|
||||
update!(element, "youngs modulus", 288.0)
|
||||
update!(element, "poissons ratio", 1/3)
|
||||
```
|
||||
|
||||
First define a field problem and add element to it
|
||||
```@example 1
|
||||
body = Problem(Elasticity, "test problem", 2)
|
||||
update!(body.properties,
|
||||
"formulation" => "plane_stress",
|
||||
"finite_strain" => "false",
|
||||
"geometric_stiffness" => "false")
|
||||
body.elements = [element]
|
||||
```
|
||||
|
||||
Then create element to carry on pressure
|
||||
```@example 1
|
||||
tr_el = Element(Seg2, [3, 4])
|
||||
update!(tr_el, "geometry", X)
|
||||
update!(tr_el, "displacement traction force 2", 288.0)
|
||||
traction = Problem(Elasticity, "pressure on top of block", 2)
|
||||
update!(traction.properties,
|
||||
"formulation" => "plane_stress",
|
||||
"finite_strain" => "false",
|
||||
"geometric_stiffness" => "false")
|
||||
traction.elements = [tr_el]
|
||||
```
|
||||
|
||||
Create boundary condition to support block at bottom and left
|
||||
```@example 1
|
||||
bc_el_1 = Element(Seg2, [1, 2])
|
||||
bc_el_2 = Element(Seg2, [4, 1])
|
||||
update!(bc_el_1, "displacement 2", 0.0)
|
||||
update!(bc_el_2, "displacement 1", 0.0)
|
||||
bc = Problem(Dirichlet, "add symmetry bc", 2, "displacement")
|
||||
bc.elements = [bc_el_1, bc_el_2]
|
||||
```
|
||||
|
||||
Last thing is to create a solver, push problem to solver and solve:
|
||||
```@example 1
|
||||
solver = Solver(Linear, body, traction, bc)
|
||||
solver()
|
||||
```
|
||||
|
||||
Displacement in node 3 is
|
||||
```@example 1
|
||||
solver("displacement", 0.0)[3]
|
||||
```
|
||||
@@ -1,29 +0,0 @@
|
||||
===============
|
||||
FEATURES / TODO
|
||||
===============
|
||||
|
||||
- Parallel design
|
||||
- Sparse matrices
|
||||
- Discontinuous Galerkin?
|
||||
- 100% Tested
|
||||
- Automatic coverage (Travis CI)
|
||||
- Doctest examples and tutorials
|
||||
- Intuitive to use
|
||||
- Marketing
|
||||
- Engineering porn
|
||||
- Field functions
|
||||
- Multiphysics platform
|
||||
- Elasticity
|
||||
- Thermal implemented
|
||||
- Transient and steady solvers
|
||||
- Implicit dynamics
|
||||
- Nonlinearities
|
||||
- Geometrical
|
||||
- Material
|
||||
- Contact
|
||||
- Mortar
|
||||
- Modular design
|
||||
- E.g. contact formulation can be altered by user
|
||||
- Mesh and results format, Xdmf?
|
||||
- No scalars but field variables, e.g. no constant 210GPa for steel, we interpolate variable from nodes.
|
||||
- (Constant is a special case of field variable).
|
||||
@@ -1,31 +0,0 @@
|
||||
Links
|
||||
=====
|
||||
|
||||
Discretization
|
||||
--------------
|
||||
http://code.activestate.com/recipes/579021-delaunay-triangulation/
|
||||
|
||||
|
||||
Interpolation
|
||||
-------------
|
||||
- http://www.cs.rpi.edu/~flaherje/pdf/fea4.pdf
|
||||
- http://www.sd.ruhr-uni-bochum.de/downloads/Shape_funct.pdf
|
||||
- http://what-when-how.com/the-finite-element-method/fem-for-3d-solids-finite-element-method-part-1/
|
||||
- http://www.uni-tuebingen.de/zag/teaching/environmental_modeling/S4B_FiniteElements.pdf
|
||||
- http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch10.d/AFEM.Ch10.pdf
|
||||
- http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch10.d/AFEM.Ch10.Slides.d/AFEM.Ch10.Slides.pdf
|
||||
- http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.AppI.d/AFEM.AppI.pdf
|
||||
- http://www.code-aster.org/V2/doc/default/en/man_r/r3/r3.01.01.pdf
|
||||
- http://www.researchgate.net/publication/267082822_Unified_isoparametric_3D_Lagrange_finite_elements
|
||||
|
||||
Integration
|
||||
-----------
|
||||
- http://arxiv.org/pdf/1411.1341.pdf
|
||||
|
||||
Hierarchial shape functions
|
||||
---------------------------
|
||||
- https://www.math.vt.edu/people/adjerids/research/papers/basis.pdf
|
||||
|
||||
Solvers
|
||||
-------
|
||||
- https://github.com/JuliaSparse/MultiFrontalCholesky.jl
|
||||
@@ -1,52 +0,0 @@
|
||||
---
|
||||
title: "JuliaFEM User Manual"
|
||||
description: "Quick start guide and practical tutorials for end users"
|
||||
date: 2025-11-09
|
||||
author: "Jukka Aho"
|
||||
categories: ["user guide", "tutorial"]
|
||||
keywords: ["juliafem", "quick start", "tutorial", "examples"]
|
||||
audience: "end users"
|
||||
level: "beginner"
|
||||
type: "manual"
|
||||
---
|
||||
|
||||
# JuliaFEM User Manual
|
||||
|
||||
**Audience:** End users, engineers, students who want to run simulations and get results.
|
||||
|
||||
This manual is designed to be **simple and practical** - get you from zero to running simulations as quickly as possible.
|
||||
|
||||
## What's Here
|
||||
|
||||
- **Quick Start:** Installation and first simulation in 5 minutes
|
||||
- **Tutorials:** Step-by-step guides for common problems
|
||||
- **Examples:** Pre-built simulations you can run and modify
|
||||
- **API Reference:** Function documentation (what does this do?)
|
||||
- **Troubleshooting:** Common errors and how to fix them
|
||||
|
||||
## What's NOT Here
|
||||
|
||||
- Deep theory (see `docs/book/` for that)
|
||||
- How to contribute code (see `docs/contributor/` for that)
|
||||
- Internal architecture details
|
||||
|
||||
## Philosophy
|
||||
|
||||
**"Just show me how to solve my problem."**
|
||||
|
||||
We assume you:
|
||||
- Have a problem to solve (heat transfer, elasticity, contact)
|
||||
- Want working code, not lectures
|
||||
- Will read theory when YOU need it, not when WE think you should
|
||||
|
||||
## Getting Help
|
||||
|
||||
1. Start with the Quick Start
|
||||
2. Find an example similar to your problem
|
||||
3. Modify it to fit your needs
|
||||
4. If stuck, check Troubleshooting
|
||||
5. Still stuck? Ask on GitHub Discussions
|
||||
|
||||
---
|
||||
|
||||
**Next:** Start with [Quick Start](quickstart.md) or browse [Examples](../examples/)
|
||||
@@ -1,361 +0,0 @@
|
||||
---
|
||||
title: "Quick Start: Linear Elasticity"
|
||||
date: 2025-11-10
|
||||
author: "JuliaFEM Team"
|
||||
status: "Draft"
|
||||
tags: ["user-guide", "elasticity", "quickstart", "tutorial"]
|
||||
---
|
||||
|
||||
## Overview
|
||||
|
||||
This guide shows how to solve a simple linear elasticity problem using JuliaFEM's GPU-accelerated solver.
|
||||
|
||||
**What you'll learn:**
|
||||
|
||||
- How to create a mesh with Gmsh
|
||||
- How to define materials and boundary conditions
|
||||
- How to solve and visualize results
|
||||
|
||||
**Example problem:** Cantilever beam with fixed end and pressure load.
|
||||
|
||||
## Step 1: Create the Mesh
|
||||
|
||||
First, create a 3D mesh using Gmsh:
|
||||
|
||||
```julia
|
||||
using Gmsh
|
||||
|
||||
# Initialize Gmsh
|
||||
gmsh.initialize()
|
||||
gmsh.model.add("cantilever")
|
||||
|
||||
# Geometry: 10m × 1m × 1m beam
|
||||
L, W, H = 10.0, 1.0, 1.0
|
||||
box = gmsh.model.occ.addBox(0, 0, 0, L, W, H)
|
||||
gmsh.model.occ.synchronize()
|
||||
|
||||
# Define physical groups for boundary conditions
|
||||
surfaces = gmsh.model.getBoundary([(3, box)], false, false, true)
|
||||
|
||||
for surf in surfaces
|
||||
_, surf_id = surf
|
||||
# Get surface center to identify it
|
||||
com = gmsh.model.occ.getCenterOfMass(2, surf_id)
|
||||
|
||||
if abs(com[1]) < 1e-6 # X = 0 (fixed end)
|
||||
gmsh.model.addPhysicalGroup(2, [surf_id], -1, "FixedEnd")
|
||||
elseif abs(com[3] - H) < 1e-6 # Z = H (top surface for pressure)
|
||||
gmsh.model.addPhysicalGroup(2, [surf_id], -1, "PressureSurface")
|
||||
end
|
||||
end
|
||||
|
||||
# Add volume physical group
|
||||
gmsh.model.addPhysicalGroup(3, [box], -1, "Volume")
|
||||
|
||||
# Generate 3D tetrahedral mesh
|
||||
gmsh.option.setNumber("Mesh.MeshSizeMax", 0.5)
|
||||
gmsh.model.mesh.generate(3)
|
||||
|
||||
# Save mesh
|
||||
gmsh.write("cantilever_beam.msh")
|
||||
gmsh.finalize()
|
||||
```
|
||||
|
||||
**Key concepts:**
|
||||
|
||||
- Physical groups label surfaces/volumes for boundary conditions
|
||||
- `"FixedEnd"` - nodes that will be constrained (Dirichlet BC)
|
||||
- `"PressureSurface"` - nodes where pressure is applied (Neumann BC)
|
||||
- `"Volume"` - elements for assembly
|
||||
|
||||
## Step 2: Load the Mesh
|
||||
|
||||
```julia
|
||||
include("src/gpu_elasticity.jl")
|
||||
using .GPUElasticity
|
||||
|
||||
# Read the mesh file
|
||||
mesh = read_gmsh_mesh("cantilever_beam.msh")
|
||||
|
||||
println("Mesh info:")
|
||||
println(" Nodes: $(size(mesh.nodes, 2))")
|
||||
println(" Elements: $(size(mesh.elements, 2))")
|
||||
```
|
||||
|
||||
**What you get:**
|
||||
|
||||
- `mesh.nodes` - 3×n_nodes matrix of coordinates
|
||||
- `mesh.elements` - 4×n_elements matrix of connectivity (Tet4)
|
||||
- `mesh.physical_groups` - Dictionary mapping names to node/element IDs
|
||||
|
||||
## Step 3: Extract Boundary Condition Nodes
|
||||
|
||||
```julia
|
||||
using .GPUElasticity.GmshReader: get_surface_nodes
|
||||
|
||||
# Get nodes for each boundary condition
|
||||
fixed_nodes = get_surface_nodes(mesh, "FixedEnd")
|
||||
pressure_nodes = get_surface_nodes(mesh, "PressureSurface")
|
||||
|
||||
println("Boundary conditions:")
|
||||
println(" Fixed nodes: $(length(fixed_nodes))")
|
||||
println(" Pressure nodes: $(length(pressure_nodes))")
|
||||
```
|
||||
|
||||
**Boundary condition types:**
|
||||
|
||||
1. **Dirichlet (fixed_nodes):** Zero displacement constraint
|
||||
- Nodes cannot move (u = 0)
|
||||
- Models supports, clamps, symmetry
|
||||
|
||||
2. **Neumann (pressure_nodes):** Applied force/pressure
|
||||
- External load on surface
|
||||
- Models traction, pressure, point forces
|
||||
|
||||
## Step 4: Define Material
|
||||
|
||||
```julia
|
||||
# Create material (steel)
|
||||
material = ElasticMaterial(
|
||||
210e9, # E - Young's modulus [Pa]
|
||||
0.3 # ν - Poisson's ratio [-]
|
||||
)
|
||||
```
|
||||
|
||||
**Common materials:**
|
||||
|
||||
| Material | E (GPa) | ν |
|
||||
|----------|---------|---|
|
||||
| Steel | 200-210 | 0.27-0.30 |
|
||||
| Aluminum | 69 | 0.33 |
|
||||
| Concrete | 30-40 | 0.15-0.20 |
|
||||
| Rubber | 0.01-0.1 | 0.48-0.50 |
|
||||
|
||||
## Step 5: Create Physics Problem
|
||||
|
||||
```julia
|
||||
# Define the complete problem
|
||||
physics = ElasticityPhysics(
|
||||
mesh, # Mesh with geometry
|
||||
material, # Material properties
|
||||
fixed_nodes, # Dirichlet BC nodes
|
||||
pressure_nodes, # Neumann BC nodes
|
||||
10e6 # Pressure magnitude [Pa] = 10 MPa
|
||||
)
|
||||
```
|
||||
|
||||
**What ElasticityPhysics contains:**
|
||||
|
||||
- Mesh (nodes, elements, connectivity)
|
||||
- Material (E, ν for isotropic linear elasticity)
|
||||
- Fixed nodes (where displacement = 0)
|
||||
- Pressure nodes (where external load is applied)
|
||||
- Pressure value (load magnitude)
|
||||
|
||||
## Step 6: Solve
|
||||
|
||||
```julia
|
||||
# Solve on GPU with iterative solver
|
||||
result = solve_elasticity_gpu(physics, tol=1e-6, max_iter=1000)
|
||||
|
||||
println("\nSolution converged!")
|
||||
println(" CG iterations: $(result.iterations)")
|
||||
println(" Final residual: $(result.residual)")
|
||||
```
|
||||
|
||||
**Solver parameters:**
|
||||
|
||||
- `tol` - Convergence tolerance (default: 1e-6)
|
||||
- `max_iter` - Maximum CG iterations (default: 1000)
|
||||
|
||||
**What you get:**
|
||||
|
||||
- `result.u` - Displacement field (3n_nodes vector)
|
||||
- `result.iterations` - Number of CG iterations
|
||||
- `result.residual` - Final residual norm
|
||||
|
||||
## Step 7: Post-Process Results
|
||||
|
||||
```julia
|
||||
# Extract displacement components
|
||||
n_nodes = size(mesh.nodes, 2)
|
||||
u_x = result.u[1:3:end]
|
||||
u_y = result.u[2:3:end]
|
||||
u_z = result.u[3:3:end]
|
||||
|
||||
# Find maximum displacement
|
||||
u_magnitude = sqrt.(u_x.^2 + u_y.^2 + u_z.^2)
|
||||
max_disp = maximum(u_magnitude)
|
||||
max_node = argmax(u_magnitude)
|
||||
|
||||
println("\nResults:")
|
||||
println(" Max displacement: $(max_disp * 1000) mm")
|
||||
println(" At node: $(max_node)")
|
||||
println(" Location: $(mesh.nodes[:, max_node])")
|
||||
|
||||
# Compute stresses (requires element loop - TODO)
|
||||
```
|
||||
|
||||
## Complete Example
|
||||
|
||||
Here's the full script combining all steps:
|
||||
|
||||
```julia
|
||||
using Gmsh
|
||||
include("src/gpu_elasticity.jl")
|
||||
using .GPUElasticity
|
||||
using .GPUElasticity.GmshReader: get_surface_nodes
|
||||
|
||||
# 1. Generate mesh
|
||||
gmsh.initialize()
|
||||
gmsh.model.add("cantilever")
|
||||
L, W, H = 10.0, 1.0, 1.0
|
||||
box = gmsh.model.occ.addBox(0, 0, 0, L, W, H)
|
||||
gmsh.model.occ.synchronize()
|
||||
|
||||
# Label surfaces
|
||||
surfaces = gmsh.model.getBoundary([(3, box)], false, false, true)
|
||||
for surf in surfaces
|
||||
_, surf_id = surf
|
||||
com = gmsh.model.occ.getCenterOfMass(2, surf_id)
|
||||
if abs(com[1]) < 1e-6
|
||||
gmsh.model.addPhysicalGroup(2, [surf_id], -1, "FixedEnd")
|
||||
elseif abs(com[3] - H) < 1e-6
|
||||
gmsh.model.addPhysicalGroup(2, [surf_id], -1, "PressureSurface")
|
||||
end
|
||||
end
|
||||
gmsh.model.addPhysicalGroup(3, [box], -1, "Volume")
|
||||
|
||||
# Generate and save
|
||||
gmsh.option.setNumber("Mesh.MeshSizeMax", 0.5)
|
||||
gmsh.model.mesh.generate(3)
|
||||
gmsh.write("cantilever.msh")
|
||||
gmsh.finalize()
|
||||
|
||||
# 2. Load mesh
|
||||
mesh = read_gmsh_mesh("cantilever.msh")
|
||||
|
||||
# 3. Define boundary conditions
|
||||
fixed_nodes = get_surface_nodes(mesh, "FixedEnd")
|
||||
pressure_nodes = get_surface_nodes(mesh, "PressureSurface")
|
||||
|
||||
# 4. Define material (steel)
|
||||
material = ElasticMaterial(210e9, 0.3)
|
||||
|
||||
# 5. Create physics
|
||||
physics = ElasticityPhysics(
|
||||
mesh,
|
||||
material,
|
||||
fixed_nodes,
|
||||
pressure_nodes,
|
||||
10e6 # 10 MPa pressure
|
||||
)
|
||||
|
||||
# 6. Solve
|
||||
result = solve_elasticity_gpu(physics)
|
||||
|
||||
# 7. Results
|
||||
n_nodes = size(mesh.nodes, 2)
|
||||
u_mag = sqrt.(
|
||||
result.u[1:3:end].^2 +
|
||||
result.u[2:3:end].^2 +
|
||||
result.u[3:3:end].^2
|
||||
)
|
||||
|
||||
println("Max displacement: $(maximum(u_mag) * 1000) mm")
|
||||
```
|
||||
|
||||
## Current Limitations (Linear Elasticity)
|
||||
|
||||
The current implementation (`gpu_elasticity.jl`) supports:
|
||||
|
||||
✅ **Working:**
|
||||
|
||||
- Linear elastic material (Hooke's law)
|
||||
- Isotropic materials (E, ν constant)
|
||||
- Small strain assumption
|
||||
- Dirichlet BC (fixed displacement)
|
||||
- Neumann BC (pressure on surfaces)
|
||||
- Matrix-free CG solver
|
||||
- GPU acceleration
|
||||
|
||||
❌ **Not yet implemented:**
|
||||
|
||||
- Nonlinear materials (plasticity, hyperelasticity)
|
||||
- Large deformations (geometric nonlinearity)
|
||||
- Material state variables (plastic strain, damage)
|
||||
- Contact mechanics
|
||||
- Dynamic analysis (time integration)
|
||||
- Point forces (only surface pressure)
|
||||
|
||||
## Next Steps
|
||||
|
||||
**To extend to nonlinear elasticity**, we need:
|
||||
|
||||
1. **Material state at integration points**
|
||||
- Store plastic strain εₚ, hardening α, etc.
|
||||
- Update state during Newton iterations
|
||||
|
||||
2. **Newton-Raphson solver**
|
||||
- Replace CG with Newton loop
|
||||
- Compute tangent stiffness and residual
|
||||
- Line search for globalization
|
||||
|
||||
3. **Stress update algorithms**
|
||||
- Radial return for plasticity
|
||||
- Hyperelastic stress from strain energy
|
||||
- State management (old vs new state)
|
||||
|
||||
4. **Boundary condition updates**
|
||||
- Prescribed displacement (not just zero)
|
||||
- Follower forces (load direction changes)
|
||||
- Contact constraints
|
||||
|
||||
See `docs/src/book/` for design documents on these extensions.
|
||||
|
||||
## Troubleshooting
|
||||
|
||||
### Gmsh not found
|
||||
|
||||
```julia
|
||||
using Pkg
|
||||
Pkg.add("Gmsh")
|
||||
```
|
||||
|
||||
### No CUDA device
|
||||
|
||||
CPU-only version coming soon. For now, requires NVIDIA GPU with CUDA.
|
||||
|
||||
### CG doesn't converge
|
||||
|
||||
- Increase `max_iter` parameter
|
||||
- Check boundary conditions (mesh must be constrained)
|
||||
- Add preconditioner (future work)
|
||||
|
||||
### Out of GPU memory
|
||||
|
||||
- Reduce mesh size (fewer elements)
|
||||
- Use coarser mesh (`Mesh.MeshSizeMax` larger)
|
||||
- Future: Distributed multi-GPU solver
|
||||
|
||||
## Where to Learn More
|
||||
|
||||
- **Theory:** `docs/src/book/elasticity_theory.md`
|
||||
- **Implementation:** `src/gpu_elasticity.jl` (477 lines, well-commented)
|
||||
- **Test:** `test/test_gpu_elasticity.jl`
|
||||
- **Demo:** `demos/cantilever_beam_demo.jl`
|
||||
- **Design:** `docs/src/book/design/gpu_elasticity_implementation.md`
|
||||
|
||||
## Summary
|
||||
|
||||
**Workflow:**
|
||||
|
||||
1. Generate mesh with Gmsh (label surfaces for BCs)
|
||||
2. Load mesh into JuliaFEM
|
||||
3. Extract boundary condition nodes
|
||||
4. Define material properties
|
||||
5. Create `ElasticityPhysics` struct
|
||||
6. Solve with `solve_elasticity_gpu()`
|
||||
7. Post-process displacement field
|
||||
|
||||
**Current status:** Linear elasticity works. Nonlinear extensions in progress.
|
||||
@@ -1,574 +0,0 @@
|
||||
---
|
||||
title: "JuliaFEM System Architecture: Core Concepts"
|
||||
date: 2025-11-10
|
||||
author: "JuliaFEM Team"
|
||||
status: "Authoritative"
|
||||
last_updated: 2025-11-10
|
||||
tags: ["architecture", "design", "concepts", "user-guide"]
|
||||
---
|
||||
|
||||
## Introduction
|
||||
|
||||
This document explains **why** JuliaFEM's core data structures exist and **how**
|
||||
they work together. Every struct has a purpose, every abstraction has a reason.
|
||||
|
||||
**Target audience:** Users who want to understand the system deeply,
|
||||
contributors implementing new physics, anyone asking "why is it designed this
|
||||
way?"
|
||||
|
||||
---
|
||||
|
||||
## The Big Picture: Method of Lines
|
||||
|
||||
JuliaFEM follows the **Method of Lines** approach to FEM:
|
||||
|
||||
1. **Spatial discretization** (FEM) → System of ODEs
|
||||
2. **Time discretization** (if transient) → Nonlinear algebraic equations
|
||||
3. **Linearization** (Newton) → Linear system solve
|
||||
4. **Repeat** until convergence
|
||||
|
||||
**Key insight:** Separate concerns cleanly at each level!
|
||||
|
||||
```text
|
||||
Mesh + Physics → Element Assembly → Global System → Solver → Solution
|
||||
↓ ↓ ↓ ↓ ↓
|
||||
Geometry Material Models Sparse Matrix Krylov u(x,t)
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Core Abstractions
|
||||
|
||||
### 1. `Element` - The Geometric Container
|
||||
|
||||
**Role:** Holds geometric and connectivity information for a single finite element.
|
||||
|
||||
**What it knows:**
|
||||
|
||||
- Element type (Tri3, Quad4, Tet10, etc.)
|
||||
- Node connectivity
|
||||
- Basis functions (shape functions)
|
||||
- Integration points
|
||||
|
||||
**What it does NOT know:**
|
||||
|
||||
- Physics equations
|
||||
- Material properties
|
||||
- Boundary conditions
|
||||
|
||||
**Why this design?**
|
||||
|
||||
- ✅ Element is **reusable** across different physics (same Tet10 for elasticity, heat, fluid)
|
||||
- ✅ Geometry is **immutable** (connectivity doesn't change during analysis)
|
||||
- ✅ **Type-stable dispatch** on element type enables compiler optimizations
|
||||
|
||||
**Example:**
|
||||
|
||||
```julia
|
||||
# Create a 10-node tetrahedral element
|
||||
nodes = [1, 5, 12, 23, 14, 8, 19, 27, 31, 16]
|
||||
element = Element(Tet10, nodes)
|
||||
|
||||
# Element knows its topology
|
||||
@assert nnodes(element) == 10
|
||||
@assert dim(element) == 3
|
||||
|
||||
# But element doesn't know about stress, temperature, etc.
|
||||
# That's the job of Physics!
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
### 2. `Physics` - The Equation Selector (RENAMED from "Problem")
|
||||
|
||||
**Role:** Multiple dispatch tag + configuration for physical equations.
|
||||
|
||||
**What it is:**
|
||||
|
||||
- A **type** that selects assembly methods via multiple dispatch
|
||||
- A **struct** that holds physics-specific configuration
|
||||
- A **name provider** for field names ("displacement", "temperature", etc.)
|
||||
|
||||
**What it is NOT:**
|
||||
|
||||
- ❌ Not the mesh (that's separate)
|
||||
- ❌ Not the material models (those are parameters)
|
||||
- ❌ Not the solver (that's a different layer)
|
||||
|
||||
**Why "Physics" instead of "Problem"?**
|
||||
|
||||
- ✅ **Positive connotation** ("solve physics" vs "solve problem")
|
||||
- ✅ **Accurate description** (selecting physical equations)
|
||||
- ✅ **Clear role** (what physics are we simulating?)
|
||||
|
||||
**Types of Physics:**
|
||||
|
||||
```julia
|
||||
# Elasticity: Solves ∇⋅σ = ρü + b
|
||||
struct ElasticityPhysics <: AbstractPhysics
|
||||
formulation::Symbol # :plane_stress, :plane_strain, :continuum
|
||||
finite_strain::Bool # Geometric nonlinearity
|
||||
geometric_stiffness::Bool # σ-dependent stiffness for buckling
|
||||
store_fields::Vector{Symbol} # Output fields to save
|
||||
end
|
||||
|
||||
# Heat transfer: Solves ∇⋅(k∇T) = ρcₚ∂T/∂t + Q
|
||||
struct HeatPhysics <: AbstractPhysics
|
||||
formulation::Symbol # :steady_state, :transient
|
||||
nonlinear::Bool # Temperature-dependent properties
|
||||
store_fields::Vector{Symbol}
|
||||
end
|
||||
|
||||
# Contact mechanics: Solves contact constraints
|
||||
struct ContactPhysics <: AbstractPhysics
|
||||
algorithm::Symbol # :penalty, :lagrange, :augmented_lagrange
|
||||
friction_model::Symbol # :coulomb, :frictionless
|
||||
# ...
|
||||
end
|
||||
```
|
||||
|
||||
**Multiple Dispatch in Action:**
|
||||
|
||||
```julia
|
||||
# Compiler selects correct assembly method based on physics type!
|
||||
|
||||
function assemble!(assembly, physics::ElasticityPhysics, elements, time)
|
||||
# Elasticity-specific assembly:
|
||||
# - Compute strain from displacement
|
||||
# - Call material model: ε → (σ, 𝔻)
|
||||
# - Build stiffness matrix K and force vector f
|
||||
end
|
||||
|
||||
function assemble!(assembly, physics::HeatPhysics, elements, time)
|
||||
# Heat transfer-specific assembly:
|
||||
# - Compute temperature gradient
|
||||
# - Call thermal conductivity: ∇T → q
|
||||
# - Build capacity matrix C and conductivity matrix K
|
||||
end
|
||||
|
||||
function assemble!(assembly, physics::ContactPhysics, elements, time)
|
||||
# Contact-specific assembly:
|
||||
# - Detect penetration
|
||||
# - Compute contact forces
|
||||
# - Build constraint equations
|
||||
end
|
||||
```
|
||||
|
||||
**This is Julia's superpower!** No runtime type checks, no vtables, just fast compiled code for each physics type.
|
||||
|
||||
---
|
||||
|
||||
### 3. `Material` - The Constitutive Model
|
||||
|
||||
**Role:** Maps kinematic quantities to stress/flux/response.
|
||||
|
||||
**Interface (Elasticity example):**
|
||||
|
||||
```julia
|
||||
abstract type AbstractMaterial end
|
||||
abstract type AbstractMaterialState end
|
||||
|
||||
# Material computes: (ε, state_old, Δt) → (σ, 𝔻, state_new)
|
||||
function compute_stress(
|
||||
material::AbstractMaterial,
|
||||
ε::SymmetricTensor{2,3},
|
||||
state_old::AbstractMaterialState,
|
||||
Δt::Float64
|
||||
) -> Tuple{SymmetricTensor{2,3}, SymmetricTensor{4,3}, AbstractMaterialState}
|
||||
# Returns: (stress, tangent, state_new)
|
||||
end
|
||||
```
|
||||
|
||||
**Why separate from Physics?**
|
||||
|
||||
- ✅ **Modularity:** Change material without touching assembly code
|
||||
- ✅ **Testability:** Unit-test materials independently
|
||||
- ✅ **Performance:** Compiler specializes on material type
|
||||
- ✅ **Clarity:** Material logic isolated from kinematics
|
||||
|
||||
**Material examples:**
|
||||
|
||||
```julia
|
||||
# Stateless material (no history)
|
||||
struct LinearElastic <: AbstractMaterial
|
||||
λ::Float64 # Lamé parameter
|
||||
μ::Float64 # Shear modulus
|
||||
end
|
||||
|
||||
# Stateful material (history-dependent)
|
||||
struct PerfectPlasticity <: AbstractMaterial
|
||||
E::Float64 # Young's modulus
|
||||
ν::Float64 # Poisson's ratio
|
||||
σ_y::Float64 # Yield stress
|
||||
end
|
||||
|
||||
struct PlasticityState{T} <: AbstractMaterialState
|
||||
εₚ::SymmetricTensor{2,3,T} # Plastic strain
|
||||
α::T # Equivalent plastic strain
|
||||
end
|
||||
```
|
||||
|
||||
**Performance:** 20-70 ns per material evaluation (validated Nov 10, 2025)
|
||||
|
||||
---
|
||||
|
||||
### 4. `Assembly` - The Global System Builder
|
||||
|
||||
**Role:** Accumulate element contributions into global matrices/vectors.
|
||||
|
||||
**What it holds:**
|
||||
|
||||
- Global stiffness matrix `K` (sparse)
|
||||
- Global force vector `f`
|
||||
- (Optional) Mass matrix `M`, damping `C`, geometric stiffness `Kg`
|
||||
|
||||
**What it does:**
|
||||
|
||||
- Pre-allocates sparse matrix structure
|
||||
- Accumulates element contributions: `K += Kₑ`, `f += fₑ`
|
||||
- Handles DOF mapping: local element DOFs → global system DOFs
|
||||
|
||||
**Why separate from Physics?**
|
||||
|
||||
- ✅ **Reusability:** Same Assembly struct for all physics types
|
||||
- ✅ **Optimization:** Pre-allocated structure, efficient COO→CSC conversion
|
||||
- ✅ **Parallelism:** (Future) Thread-safe assembly with color-based locking
|
||||
|
||||
**Example:**
|
||||
|
||||
```julia
|
||||
# Create assembly
|
||||
assembly = Assembly()
|
||||
|
||||
# Loop over elements
|
||||
for element in elements
|
||||
# Compute element stiffness and force
|
||||
Kₑ, fₑ = assemble_element(physics, element, time)
|
||||
|
||||
# Get global DOF indices
|
||||
gdofs = get_gdofs(element)
|
||||
|
||||
# Add to global system
|
||||
add!(assembly.K, gdofs, gdofs, Kₑ)
|
||||
add!(assembly.f, gdofs, fₑ)
|
||||
end
|
||||
|
||||
# Solve global system
|
||||
u = assembly.K \ assembly.f
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
### 5. `IntegrationPoint` - The Quadrature Point
|
||||
|
||||
**Role:** Location and weight for numerical integration.
|
||||
|
||||
**What it knows:**
|
||||
|
||||
- Position in reference element: `ξ ∈ [-1,1]ᵈⁱᵐ`
|
||||
- Integration weight: `w`
|
||||
- Index/ID for state storage
|
||||
|
||||
**What it does NOT know:**
|
||||
|
||||
- Material state (that's stored per-element per-IP)
|
||||
- Stress/strain (that's computed on-the-fly)
|
||||
|
||||
**Why this design?**
|
||||
|
||||
- ✅ **Immutable:** Integration points never change
|
||||
- ✅ **Topology-specific:** Different rules for Tri3 vs Quad4
|
||||
- ✅ **Pre-computed:** Created once, reused forever
|
||||
|
||||
**Example:**
|
||||
|
||||
```julia
|
||||
# Get integration points for element type
|
||||
ips = integration_points(Gauss{2}, Quad4)
|
||||
|
||||
# Each IP knows position and weight
|
||||
for ip in ips
|
||||
ξ = ip.ξ # Position: NTuple{2, Float64}
|
||||
w = ip.weight # Weight: Float64
|
||||
|
||||
# Evaluate basis functions at this point
|
||||
N, ∇N_ref = evaluate_basis(basis, ξ)
|
||||
|
||||
# Do integration: ∫f dΩ ≈ ∑ᵢ f(ξᵢ)⋅w(ξᵢ)
|
||||
end
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
### 6. `BasisInfo` - The Shape Function Cache
|
||||
|
||||
**Role:** Pre-allocated workspace for basis function evaluation.
|
||||
|
||||
**What it caches:**
|
||||
|
||||
- Basis function values `N`
|
||||
- Gradients in reference config `∇N_ref`
|
||||
- Gradients in current config `∇N`
|
||||
- Jacobian `J`, determinant `detJ`
|
||||
|
||||
**Why cache?**
|
||||
|
||||
- ✅ **Zero allocation:** Reuse same arrays for every element
|
||||
- ✅ **Type stability:** All sizes known at compile time
|
||||
- ✅ **Performance:** Avoid repeated memory allocation
|
||||
|
||||
**Example:**
|
||||
|
||||
```julia
|
||||
# Create cache for Tet10 elements
|
||||
bi = BasisInfo(Tet10)
|
||||
|
||||
# Reuse for every element
|
||||
for element in elements
|
||||
for ip in integration_points(element)
|
||||
# Evaluate into pre-allocated cache
|
||||
eval_basis!(bi, element.X, ip)
|
||||
|
||||
# Access cached results
|
||||
N = bi.N # Shape functions
|
||||
∇N = bi.grad # Gradients ∂N/∂x
|
||||
w = ip.weight * bi.detJ # Integration weight
|
||||
end
|
||||
end
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Data Flow: From Mesh to Solution
|
||||
|
||||
### Step 1: Problem Setup
|
||||
|
||||
```julia
|
||||
# Define physics
|
||||
physics = ElasticityPhysics(
|
||||
formulation = :continuum,
|
||||
finite_strain = false,
|
||||
store_fields = [:stress, :strain]
|
||||
)
|
||||
|
||||
# Define material
|
||||
material = LinearElastic(E=200e9, ν=0.3)
|
||||
|
||||
# Create elements with material
|
||||
elements = [Element(Tet10, conn) for conn in connectivity]
|
||||
for el in elements
|
||||
el.material = material
|
||||
el.states_old = [NoState() for _ in 1:n_integration_points]
|
||||
el.states_new = [NoState() for _ in 1:n_integration_points]
|
||||
end
|
||||
```
|
||||
|
||||
### Step 2: Assembly Loop
|
||||
|
||||
```julia
|
||||
assembly = Assembly()
|
||||
|
||||
for element in elements
|
||||
# Get element data
|
||||
X = element.geometry # Nodal coordinates
|
||||
u = element.displacement # Nodal displacements
|
||||
|
||||
# Initialize element matrices
|
||||
Kₑ = zeros(ndofs, ndofs)
|
||||
fₑ = zeros(ndofs)
|
||||
|
||||
# Integration point loop
|
||||
for (ip_idx, ip) in enumerate(element.integration_points)
|
||||
|
||||
# 1. KINEMATICS: u → ε
|
||||
∇N = shape_function_gradients(element, ip)
|
||||
ε = compute_strain_from_gradients(∇N, u)
|
||||
|
||||
# 2. MATERIAL: ε → (σ, 𝔻, state)
|
||||
state_old = element.states_old[ip_idx]
|
||||
σ, 𝔻, state_trial = compute_stress(material, ε, state_old, Δt)
|
||||
|
||||
# 3. ASSEMBLY: (∇N, σ, 𝔻) → (Kₑ, fₑ)
|
||||
w = integration_weight(ip)
|
||||
accumulate_stiffness!(Kₑ, ∇N, 𝔻, w)
|
||||
accumulate_internal_forces!(fₑ, ∇N, σ, w)
|
||||
|
||||
# DON'T update states yet (Newton iterations!)
|
||||
end
|
||||
|
||||
# 4. GLOBAL: Kₑ → K, fₑ → f
|
||||
gdofs = get_gdofs(element)
|
||||
add!(assembly.K, gdofs, gdofs, Kₑ)
|
||||
add!(assembly.f, gdofs, fₑ)
|
||||
end
|
||||
```
|
||||
|
||||
### Step 3: Solve
|
||||
|
||||
```julia
|
||||
# Linear solve: K⋅Δu = f
|
||||
Δu = assembly.K \ assembly.f
|
||||
|
||||
# Newton iteration (if nonlinear)
|
||||
while norm(residual) > tolerance
|
||||
# Re-assemble with trial displacement
|
||||
u_trial = u_old + Δu
|
||||
|
||||
# Solve linearized system
|
||||
Δu = assembly.K \ assembly.f
|
||||
|
||||
# Update
|
||||
u_trial += Δu
|
||||
end
|
||||
|
||||
# Converged! Commit state
|
||||
for element in elements
|
||||
element.states_old .= element.states_new
|
||||
end
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## Design Principles
|
||||
|
||||
### 1. Separation of Concerns
|
||||
|
||||
**Each struct has ONE job:**
|
||||
|
||||
- `Element` → Geometry
|
||||
- `Physics` → Equation selection
|
||||
- `Material` → Constitutive model
|
||||
- `Assembly` → Global system
|
||||
- `Solver` → Linear algebra
|
||||
|
||||
**Benefits:**
|
||||
|
||||
- ✅ Easy to test (unit test each component)
|
||||
- ✅ Easy to extend (add new material without touching assembly)
|
||||
- ✅ Easy to optimize (profile each layer independently)
|
||||
|
||||
### 2. Type Stability
|
||||
|
||||
**Every function has concrete return type:**
|
||||
|
||||
```julia
|
||||
# ✅ GOOD: Compiler knows return type
|
||||
function compute_stress(m::LinearElastic, ε) -> Tuple{SymmetricTensor{2,3}, SymmetricTensor{4,3}, NoState}
|
||||
# ...
|
||||
end
|
||||
|
||||
# ❌ BAD: Compiler doesn't know (Dict lookup)
|
||||
function compute_stress(element, ip)
|
||||
stress = element.fields["stress"] # Unknown type!
|
||||
# ...
|
||||
end
|
||||
```
|
||||
|
||||
**Performance impact:** 10-100× speedup from type stability alone!
|
||||
|
||||
### 3. Zero Allocation
|
||||
|
||||
**Hot paths allocate NOTHING:**
|
||||
|
||||
```julia
|
||||
# Pre-allocate once
|
||||
bi = BasisInfo(Tet10)
|
||||
Kₑ = zeros(30, 30)
|
||||
|
||||
# Reuse in loop (zero allocations!)
|
||||
for element in elements
|
||||
fill!(Kₑ, 0.0)
|
||||
for ip in integration_points(element)
|
||||
eval_basis!(bi, X, ip) # Fills cache, no allocation
|
||||
# ... assembly logic
|
||||
end
|
||||
end
|
||||
```
|
||||
|
||||
**Validated:** All material models achieve 0 bytes allocation (Nov 10, 2025)
|
||||
|
||||
### 4. Multiple Dispatch
|
||||
|
||||
**Use Julia's type system:**
|
||||
|
||||
```julia
|
||||
# Same function name, different implementations
|
||||
assemble!(assembly, ::ElasticityPhysics, elements, time)
|
||||
assemble!(assembly, ::HeatPhysics, elements, time)
|
||||
assemble!(assembly, ::ContactPhysics, elements, time)
|
||||
|
||||
# Compiler generates specialized code for each!
|
||||
```
|
||||
|
||||
**No runtime overhead, no vtables, just fast native code.**
|
||||
|
||||
---
|
||||
|
||||
## Common Questions
|
||||
|
||||
### Q: Why not use classes with methods?
|
||||
|
||||
**A:** Julia's multiple dispatch is more powerful than OOP:
|
||||
|
||||
```julia
|
||||
# OOP way (single dispatch on first argument)
|
||||
element.assemble(physics) # Only element type matters
|
||||
|
||||
# Julia way (multiple dispatch on ALL arguments)
|
||||
assemble!(assembly, physics, element, time) # All types matter!
|
||||
```
|
||||
|
||||
This enables:
|
||||
|
||||
- Compiler specialization on ALL argument types
|
||||
- Adding new methods without modifying existing types
|
||||
- True separation of concerns (no "god objects")
|
||||
|
||||
### Q: Why immutable structs?
|
||||
|
||||
**A:** Performance and safety:
|
||||
|
||||
- Structs with all concrete types are stack-allocated
|
||||
- Immutability enables compiler optimizations
|
||||
- No accidental mutation bugs
|
||||
|
||||
**Rule:** Use immutable structs unless you NEED mutability (like Assembly accumulation)
|
||||
|
||||
### Q: Why Tensors.jl instead of matrices?
|
||||
|
||||
**A:** Performance and clarity:
|
||||
|
||||
- `SymmetricTensor{2,3,Float64,6}` is **stack-allocated** (48 bytes on stack)
|
||||
- Regular `Matrix{Float64}` is **heap-allocated** (pointer + malloc)
|
||||
- Code looks like math: `σ = λ⋅tr(ε)⋅I + 2μ⋅ε`
|
||||
- Type stability: Compiler knows exact size at compile time
|
||||
|
||||
---
|
||||
|
||||
## Summary: The JuliaFEM Way
|
||||
|
||||
**Core philosophy:**
|
||||
|
||||
1. **Separate concerns** - Each struct has ONE job
|
||||
2. **Type stability** - Compiler knows ALL types
|
||||
3. **Zero allocation** - Hot paths reuse memory
|
||||
4. **Multiple dispatch** - Compiler specializes for each case
|
||||
5. **Immutability** - Stack allocation + safety
|
||||
6. **Tensors.jl** - Mathematical clarity + performance
|
||||
|
||||
**Result:** Fast, maintainable, extensible FEM code that looks like the math it implements.
|
||||
|
||||
---
|
||||
|
||||
## Further Reading
|
||||
|
||||
- **`docs/book/material_modeling.md`** - Material model implementation guide
|
||||
- **`docs/book/elasticity_refactoring_plan.md`** - Elasticity system design
|
||||
- **`docs/book/element_architecture.md`** - Element composition philosophy
|
||||
- **`llm/ARCHITECTURE.md`** - Eight-layer system architecture
|
||||
- **`llm/VISION_2.0.md`** - Project vision and philosophy
|
||||
|
||||
---
|
||||
|
||||
**Last Updated:** November 10, 2025
|
||||
**Status:** Authoritative (reflects current design decisions and implementation)
|
||||
Reference in New Issue
Block a user