From f0353699b0a2b64f54e105a8e08b74a264c2e196 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Fri, 12 Dec 2025 22:01:57 +0200 Subject: [PATCH] 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 --- docs/badges.jl | 53 -- docs/book/nodal_assembly_concept.md | 307 ------- docs/build_api.jl | 71 -- docs/build_doctests.jl | 52 -- docs/build_lint.jl | 17 - docs/build_notebooks.jl | 231 ----- docs/build_unittests.jl | 10 - docs/src/contributor/README.md | 64 -- docs/src/contributor/coding_standards.md | 499 ----------- .../guides/gpu_elasticity_quickstart.md | 119 --- .../guides/gpu_nodal_assembly_quickstart.md | 262 ------ .../contributor/guides/quick_reference_gpu.md | 87 -- docs/src/contributor/status.md | 109 --- docs/src/contributor/test_fixes_needed.md | 134 --- docs/src/contributor/testing_philosophy.md | 802 ------------------ docs/src/examples.md | 64 -- docs/src/features.md | 29 - docs/src/links.md | 31 - docs/src/user/README.md | 52 -- docs/src/user/elasticity_quickstart.md | 361 -------- docs/src/user/system_architecture.md | 574 ------------- 21 files changed, 3928 deletions(-) delete mode 100644 docs/badges.jl delete mode 100644 docs/book/nodal_assembly_concept.md delete mode 100644 docs/build_api.jl delete mode 100644 docs/build_doctests.jl delete mode 100644 docs/build_lint.jl delete mode 100644 docs/build_notebooks.jl delete mode 100644 docs/build_unittests.jl delete mode 100644 docs/src/contributor/README.md delete mode 100644 docs/src/contributor/coding_standards.md delete mode 100644 docs/src/contributor/guides/gpu_elasticity_quickstart.md delete mode 100644 docs/src/contributor/guides/gpu_nodal_assembly_quickstart.md delete mode 100644 docs/src/contributor/guides/quick_reference_gpu.md delete mode 100644 docs/src/contributor/status.md delete mode 100644 docs/src/contributor/test_fixes_needed.md delete mode 100644 docs/src/contributor/testing_philosophy.md delete mode 100644 docs/src/examples.md delete mode 100644 docs/src/features.md delete mode 100644 docs/src/links.md delete mode 100644 docs/src/user/README.md delete mode 100644 docs/src/user/elasticity_quickstart.md delete mode 100644 docs/src/user/system_architecture.md diff --git a/docs/badges.jl b/docs/badges.jl deleted file mode 100644 index 2a1f774..0000000 --- a/docs/badges.jl +++ /dev/null @@ -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 diff --git a/docs/book/nodal_assembly_concept.md b/docs/book/nodal_assembly_concept.md deleted file mode 100644 index 8200849..0000000 --- a/docs/book/nodal_assembly_concept.md +++ /dev/null @@ -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! diff --git a/docs/build_api.jl b/docs/build_api.jl deleted file mode 100644 index c1fcf84..0000000 --- a/docs/build_api.jl +++ /dev/null @@ -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 diff --git a/docs/build_doctests.jl b/docs/build_doctests.jl deleted file mode 100644 index 345977c..0000000 --- a/docs/build_doctests.jl +++ /dev/null @@ -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 diff --git a/docs/build_lint.jl b/docs/build_lint.jl deleted file mode 100644 index f6bb08b..0000000 --- a/docs/build_lint.jl +++ /dev/null @@ -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 diff --git a/docs/build_notebooks.jl b/docs/build_notebooks.jl deleted file mode 100644 index 6d7c427..0000000 --- a/docs/build_notebooks.jl +++ /dev/null @@ -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() diff --git a/docs/build_unittests.jl b/docs/build_unittests.jl deleted file mode 100644 index e61502d..0000000 --- a/docs/build_unittests.jl +++ /dev/null @@ -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 - diff --git a/docs/src/contributor/README.md b/docs/src/contributor/README.md deleted file mode 100644 index ae1d0d8..0000000 --- a/docs/src/contributor/README.md +++ /dev/null @@ -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) diff --git a/docs/src/contributor/coding_standards.md b/docs/src/contributor/coding_standards.md deleted file mode 100644 index c7c4428..0000000 --- a/docs/src/contributor/coding_standards.md +++ /dev/null @@ -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:** `(): ` - -**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. diff --git a/docs/src/contributor/guides/gpu_elasticity_quickstart.md b/docs/src/contributor/guides/gpu_elasticity_quickstart.md deleted file mode 100644 index 6fbf08f..0000000 --- a/docs/src/contributor/guides/gpu_elasticity_quickstart.md +++ /dev/null @@ -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! 🚀** diff --git a/docs/src/contributor/guides/gpu_nodal_assembly_quickstart.md b/docs/src/contributor/guides/gpu_nodal_assembly_quickstart.md deleted file mode 100644 index 23ff8da..0000000 --- a/docs/src/contributor/guides/gpu_nodal_assembly_quickstart.md +++ /dev/null @@ -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! 🚀** diff --git a/docs/src/contributor/guides/quick_reference_gpu.md b/docs/src/contributor/guides/quick_reference_gpu.md deleted file mode 100644 index bf9de6a..0000000 --- a/docs/src/contributor/guides/quick_reference_gpu.md +++ /dev/null @@ -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! 🚀 diff --git a/docs/src/contributor/status.md b/docs/src/contributor/status.md deleted file mode 100644 index c14765c..0000000 --- a/docs/src/contributor/status.md +++ /dev/null @@ -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. diff --git a/docs/src/contributor/test_fixes_needed.md b/docs/src/contributor/test_fixes_needed.md deleted file mode 100644 index 939d079..0000000 --- a/docs/src/contributor/test_fixes_needed.md +++ /dev/null @@ -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. diff --git a/docs/src/contributor/testing_philosophy.md b/docs/src/contributor/testing_philosophy.md deleted file mode 100644 index 2e756c3..0000000 --- a/docs/src/contributor/testing_philosophy.md +++ /dev/null @@ -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 diff --git a/docs/src/examples.md b/docs/src/examples.md deleted file mode 100644 index 0af8fd1..0000000 --- a/docs/src/examples.md +++ /dev/null @@ -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] -``` diff --git a/docs/src/features.md b/docs/src/features.md deleted file mode 100644 index bfb80d1..0000000 --- a/docs/src/features.md +++ /dev/null @@ -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). diff --git a/docs/src/links.md b/docs/src/links.md deleted file mode 100644 index 2063fb9..0000000 --- a/docs/src/links.md +++ /dev/null @@ -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 diff --git a/docs/src/user/README.md b/docs/src/user/README.md deleted file mode 100644 index db12394..0000000 --- a/docs/src/user/README.md +++ /dev/null @@ -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/) diff --git a/docs/src/user/elasticity_quickstart.md b/docs/src/user/elasticity_quickstart.md deleted file mode 100644 index e99456c..0000000 --- a/docs/src/user/elasticity_quickstart.md +++ /dev/null @@ -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. diff --git a/docs/src/user/system_architecture.md b/docs/src/user/system_architecture.md deleted file mode 100644 index e6f99a8..0000000 --- a/docs/src/user/system_architecture.md +++ /dev/null @@ -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)