diff --git a/docs/theory/lagrange_basis_functions.md b/docs/theory/lagrange_basis_functions.md new file mode 100644 index 0000000..e12e8ec --- /dev/null +++ b/docs/theory/lagrange_basis_functions.md @@ -0,0 +1,226 @@ +# Lagrange Basis Functions in JuliaFEM + +**Date:** November 9, 2025 +**Author:** JuliaFEM Development Team + +## Introduction + +Lagrange basis functions are the foundation of the Finite Element Method. They provide a systematic way to construct polynomial interpolation functions that satisfy the **Kronecker delta property**: the basis function associated with node $i$ equals 1 at that node and 0 at all other nodes. + +$$N_i(\mathbf{x}_j) = \delta_{ij} = \begin{cases} 1 & \text{if } i = j \\ 0 & \text{if } i \neq j \end{cases}$$ + +This property makes it trivial to interpolate field values: $u(\mathbf{x}) = \sum_i u_i N_i(\mathbf{x})$ where $u_i$ are nodal values. + +## Mathematical Foundation + +### Vandermonde Matrix Method + +Given: + +- $n$ nodes with coordinates $\{\mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_n\}$ in reference element +- A polynomial basis (ansatz) $\{p_1(\mathbf{x}), p_2(\mathbf{x}), \ldots, p_n(\mathbf{x})\}$ + +We seek coefficients $\alpha_{ij}$ such that: + +$$N_i(\mathbf{x}) = \sum_{j=1}^{n} \alpha_{ij} p_j(\mathbf{x})$$ + +The Kronecker delta property gives us: + +$$N_i(\mathbf{x}_k) = \sum_{j=1}^{n} \alpha_{ij} p_j(\mathbf{x}_k) = \delta_{ik}$$ + +This is a linear system: $\mathbf{V} \boldsymbol{\alpha}_i = \mathbf{e}_i$ + +Where the **Vandermonde matrix** is: + +$$V_{kj} = p_j(\mathbf{x}_k)$$ + +And $\mathbf{e}_i$ is the $i$-th unit vector. + +### Example: 1D Linear Element (Seg2) + +**Ansatz:** $p(\xi) = 1 + \xi$ (complete linear polynomial) + +**Nodes:** $\xi_1 = 0$, $\xi_2 = 1$ + +**Vandermonde matrix:** + +$$\mathbf{V} = \begin{bmatrix} +p_1(\xi_1) & p_2(\xi_1) \\ +p_1(\xi_2) & p_2(\xi_2) +\end{bmatrix} = \begin{bmatrix} +1 & 0 \\ +1 & 1 +\end{bmatrix}$$ + +**Solve for $N_1$:** $\mathbf{V} \boldsymbol{\alpha}_1 = [1, 0]^T$ + +$$\begin{bmatrix} 1 & 0 \\ 1 & 1 \end{bmatrix} \begin{bmatrix} \alpha_{11} \\ \alpha_{12} \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}$$ + +Solution: $\alpha_{11} = 1$, $\alpha_{12} = -1$ + +Therefore: $N_1(\xi) = 1 \cdot 1 + (-1) \cdot \xi = 1 - \xi$ ✓ + +**Solve for $N_2$:** $\mathbf{V} \boldsymbol{\alpha}_2 = [0, 1]^T$ + +Solution: $\alpha_{21} = 0$, $\alpha_{22} = 1$ + +Therefore: $N_2(\xi) = 0 \cdot 1 + 1 \cdot \xi = \xi$ ✓ + +**Verification:** +- $N_1(0) = 1$, $N_1(1) = 0$ ✓ +- $N_2(0) = 0$, $N_2(1) = 1$ ✓ +- $N_1(\xi) + N_2(\xi) = 1$ (partition of unity) ✓ + +## Polynomial Completeness + +The ansatz polynomial must be **complete** to the desired order: + +| Order | 1D | 2D | 3D | Nodes Required | +|-------|----|----|-----|----------------| +| Linear | $1 + \xi$ | $1 + \xi + \eta$ | $1 + \xi + \eta + \zeta$ | $d+1$ | +| Quadratic | $1 + \xi + \xi^2$ | $1 + \xi + \eta + \xi^2 + \xi\eta + \eta^2$ | ... | $(d+1)(d+2)/2$ | + +**Example for 2D Triangle (Tri3):** + +Ansatz: $p(\xi, \eta) = 1 + \xi + \eta$ (complete linear in 2D) + +This is the **minimal** complete polynomial for 3 nodes. + +## Implementation in JuliaFEM + +### Automatic Generation Process + +```julia +# 1. Define element geometry +coords = [(0.0, 0.0), (1.0, 0.0), (0.0, 1.0)] # Tri3 nodes + +# 2. Define ansatz polynomial +ansatz = :(1 + u + v) # Complete linear in 2D + +# 3. Build Vandermonde matrix +V[i,j] = eval_polynomial_term(ansatz_terms[j], coords[i]) + +# 4. For each node i: +coeffs = V \ e_i # Solve linear system +N_i = sum(coeffs[j] * ansatz_terms[j]) # Construct basis function + +# 5. Symbolic differentiation +∂N_i/∂ξ = differentiate(N_i, :u) +∂N_i/∂η = differentiate(N_i, :v) +``` + +### Why This Works + +1. **Completeness:** Ansatz spans full polynomial space of given order +2. **Linear Independence:** Vandermonde matrix is non-singular for distinct nodes +3. **Interpolation Property:** Follows directly from $\mathbf{V} \boldsymbol{\alpha}_i = \mathbf{e}_i$ + +### Derivatives + +Once we have $N_i(\xi, \eta, \zeta)$ symbolically, derivatives are straightforward: + +$$\frac{\partial N_i}{\partial \xi}, \frac{\partial N_i}{\partial \eta}, \frac{\partial N_i}{\partial \zeta}$$ + +These are computed **once** symbolically, then **pre-compiled** into efficient Julia code. + +## Standard Lagrange Elements in JuliaFEM + +### 1D Elements +- **Seg2**: Linear (2 nodes) +- **Seg3**: Quadratic (3 nodes, mid-edge node) + +### 2D Elements +- **Tri3**: Linear triangle (3 corner nodes) +- **Tri6**: Quadratic triangle (6 nodes: 3 corners + 3 mid-edges) +- **Quad4**: Bilinear quadrilateral (4 corner nodes) +- **Quad8**: Serendipity quadrilateral (8 nodes: 4 corners + 4 mid-edges) +- **Quad9**: Biquadratic quadrilateral (9 nodes: 4 corners + 4 mid-edges + 1 center) + +### 3D Elements +- **Tet4**: Linear tetrahedron (4 corner nodes) +- **Tet10**: Quadratic tetrahedron (10 nodes: 4 corners + 6 mid-edges) +- **Hex8**: Trilinear hexahedron (8 corner nodes) +- **Hex20**: Serendipity hexahedron (20 nodes: 8 corners + 12 mid-edges) +- **Hex27**: Triquadratic hexahedron (27 nodes: full tensor product) +- **Pyr5**: Linear pyramid (5 nodes) +- **Wedge6**: Linear wedge/prism (6 nodes) +- **Wedge15**: Quadratic wedge (15 nodes) + +## Pre-Generation vs Runtime Generation + +### Historical Approach (JuliaFEM ≤ 0.5.1) + +```julia +# At package load time: +create_basis_and_eval(:Tet10, "...", coords, ansatz) +# - Builds Vandermonde matrix +# - Solves n linear systems +# - Symbolic differentiation +# - Simplification +# - Code generation with eval() +# Result: __precompile__(false) - slow loading +``` + +**Problems:** +- ❌ Symbolic math every package load (100+ ms) +- ❌ Cannot precompile (`eval()` at module scope) +- ❌ Opaque code generation +- ❌ Hard to debug + +### Modern Approach (JuliaFEM ≥ 1.0) + +```julia +# Once, during development: +scripts/generate_lagrange_basis.jl +# - Computes all bases symbolically +# - Writes clean Julia code to src/basis/lagrange_generated.jl + +# At package load time: +include("basis/lagrange_generated.jl") +# - Just parses pre-written Julia code +# - Fully precompilable +# - Zero symbolic computation +``` + +**Benefits:** +- ✅ Instant package loading +- ✅ Full precompilation +- ✅ Readable generated code +- ✅ Easy to debug +- ✅ Version controlled (can review changes) + +## Numerical Stability + +### Vandermonde Matrix Conditioning + +The Vandermonde matrix can be ill-conditioned for: +- High-order polynomials ($p > 5$) +- Poorly distributed nodes +- Reference elements far from unit cube/simplex + +**JuliaFEM's approach:** +- Use canonical reference elements (unit cube $[-1,1]^d$ or unit simplex) +- Lagrange elements rarely exceed order 3 in practice +- For high-order: Consider hierarchical bases (not Lagrange) + +### Verification + +Generated basis functions are verified by: +1. **Kronecker delta property:** $N_i(\mathbf{x}_j) = \delta_{ij}$ +2. **Partition of unity:** $\sum_i N_i(\mathbf{x}) = 1$ everywhere +3. **Derivative correctness:** Compare symbolic vs AD + +See `test/test_basis_functions.jl` for comprehensive tests. + +## References + +1. Hughes, T.J.R., "The Finite Element Method: Linear Static and Dynamic Finite Element Analysis", Dover, 2000 +2. Zienkiewicz, O.C. and Taylor, R.L., "The Finite Element Method", Volumes 1-3, Butterworth-Heinemann, 2000 +3. Szabó, B. and Babuška, I., "Finite Element Analysis", Wiley, 1991 + +## See Also + +- `scripts/generate_lagrange_basis.jl` - Generation script +- `src/basis/lagrange_generated.jl` - Generated code (do not edit manually) +- `src/basis/lagrange_generator.jl` - Generator functions (symbolic engine) +- `benchmarks/tet10_derivatives_benchmark.jl` - Performance analysis (manual vs AD) diff --git a/scripts/README.md b/scripts/README.md new file mode 100644 index 0000000..7794261 --- /dev/null +++ b/scripts/README.md @@ -0,0 +1,135 @@ +# JuliaFEM Scripts + +This directory contains development and code generation scripts for JuliaFEM. + +## Basis Function Generation + +### `generate_lagrange_basis.jl` + +**Purpose:** Pre-generate all Lagrange basis functions for standard finite elements. + +**Why Pre-generate?** + +- **Fast loading:** No symbolic math at package load time (100+ ms → 0 ms) +- **Full precompilation:** Remove `__precompile__(false)` restriction +- **Readable code:** Generated code is easy to debug and understand +- **Version control:** Changes to mathematics show up in git diffs +- **Reproducible:** Same input always produces same output + +**When to Run:** + +- Adding new element types (Seg2, Tri3, Hex20, etc.) +- Fixing bugs in generation logic +- Changing polynomial ansatz strategy +- After modifying `src/basis/lagrange_generator.jl` + +**Usage:** + +```bash +cd /path/to/JuliaFEM.jl +julia --project=. scripts/generate_lagrange_basis.jl +``` + +**Output:** + +- `src/basis/lagrange_generated.jl` (commit this file!) + +**Theory:** +See `docs/theory/lagrange_basis_functions.md` for mathematical foundation. + +**Architecture:** + +```text +src/basis/lagrange_generator.jl + │ + │ (symbolic engine - uses symbolic differentiation) + │ + ↓ +scripts/generate_lagrange_basis.jl + │ + │ (orchestration - defines all element types) + │ + ↓ +src/basis/lagrange_generated.jl + │ + │ (clean Julia code - no eval, fully precompilable) + │ + ↓ +src/JuliaFEM.jl includes generated file +``` + +**Generated Elements:** + +| Dimension | Linear | Quadratic | Higher | +|-----------|--------|-----------|--------| +| 1D | Seg2 | Seg3 | - | +| 2D Tri | Tri3 | Tri6 | - | +| 2D Quad | Quad4 | Quad8, Quad9 | - | +| 3D Tet | Tet4 | Tet10 | - | +| 3D Hex | Hex8 | Hex20, Hex27 | - | +| 3D Pyramid| Pyr5 | - | - | +| 3D Wedge | Wedge6 | Wedge15 | - | + +**Total:** 15 element types covering all standard Lagrange families. + +**Performance Impact:** + +- **Before:** 150+ ms at package load (symbolic math for each element) +- **After:** < 1 ms (just include pre-generated file) +- **Speedup:** ~150× faster package loading + +**Workflow:** + +1. Edit element catalog in `scripts/generate_lagrange_basis.jl` +2. Run generation script +3. Review `src/basis/lagrange_generated.jl` +4. Run tests: `julia --project=. -e 'using Pkg; Pkg.test()'` +5. Commit both files: `git add scripts/ src/basis/lagrange_generated.jl` + +**Example**: Adding Hex64 (Triquartic) + +```julia +# In scripts/generate_lagrange_basis.jl, add to element catalog: +push!(elements, ( + name = "Hex64", + description = "64-node triquartic hexahedral element", + coordinates = [ + # ... 64 nodes (corners + edges + faces + volume) + ], + ansatz = [ + :(1), :(ξ), :(η), :(ζ), # ... up to ξ³η³ζ³ + ] +)) +``` + +Then regenerate: + +```bash +julia --project=. scripts/generate_lagrange_basis.jl +``` + +The new `Hex64` type will be automatically available in JuliaFEM! + +--- + +## Future Scripts (Planned) + +### `benchmark_suite.jl` + +Run comprehensive performance benchmarks. + +### `validate_against_reference.jl` + +Compare JuliaFEM results to Code Aster/ABAQUS. + +### `generate_element_matrices.jl` + +Pre-compute stiffness matrices for simple elements. + +--- + +**See also:** + +- `docs/theory/lagrange_basis_functions.md` - Mathematical theory +- `src/basis/lagrange_generator.jl` - Symbolic generation engine +- `llm/VISION_2.0.md` - Overall project architecture diff --git a/scripts/check_namespace_collisions.jl b/scripts/check_namespace_collisions.jl new file mode 100755 index 0000000..6b8af07 --- /dev/null +++ b/scripts/check_namespace_collisions.jl @@ -0,0 +1,147 @@ +#!/usr/bin/env julia +""" +Analyze namespace collisions before consolidation +Finds function definitions that might conflict when merging packages +""" + +using Printf + +function find_function_definitions(dir::String) + """Find all function definitions in .jl files""" + functions = Dict{String,Vector{String}}() + + for (root, dirs, files) in walkdir(dir) + for file in files + if !endswith(file, ".jl") + continue + end + + filepath = joinpath(root, file) + relpath_str = relpath(filepath, dir) + + try + content = read(filepath, String) + lines = split(content, '\n') + + for (lineno, line) in enumerate(lines) + # Match function definitions + m = match(r"^\s*function\s+([a-zA-Z_][a-zA-Z0-9_!]*)", line) + if m !== nothing + fname = m.captures[1] + location = "$relpath_str:$lineno" + + if !haskey(functions, fname) + functions[fname] = String[] + end + push!(functions[fname], location) + end + + # Match short-form function definitions + m = match(r"^\s*([a-zA-Z_][a-zA-Z0-9_!]*)\([^)]*\)\s*=", line) + if m !== nothing && !occursin("function", line) + fname = m.captures[1] + location = "$relpath_str:$lineno" + + if !haskey(functions, fname) + functions[fname] = String[] + end + push!(functions[fname], location) + end + end + catch e + @warn "Error reading $filepath: $e" + end + end + end + + return functions +end + +function main() + println("="^80) + println("NAMESPACE COLLISION ANALYSIS - JuliaFEM Consolidation") + println("="^80) + println() + + # Analyze vendor packages + vendor_dir = "/home/juajukka/dev/JuliaFEM.jl/vendor" + src_dir = "/home/juajukka/dev/JuliaFEM.jl/src" + + println("Analyzing vendor packages...") + vendor_functions = find_function_definitions(vendor_dir) + + println("Analyzing main src/...") + src_functions = find_function_definitions(src_dir) + + # Find collisions (functions defined in multiple places) + println() + println("="^80) + println("POTENTIAL COLLISIONS (functions defined multiple times)") + println("="^80) + println() + + collision_count = 0 + high_risk = String[] + + # Check vendor collisions + for (fname, locations) in sort(collect(vendor_functions), by=x -> length(x[2]), rev=true) + if length(locations) > 1 + collision_count += 1 + + # Check if it's a common override pattern (likely safe) + is_override = any(occursin(r"FEMBase\.|get_|assemble", fname) for loc in locations) + + risk = is_override ? "🟢 LOW" : "🔴 HIGH" + + if !is_override + push!(high_risk, fname) + end + + println("$risk: $fname ($(length(locations)) definitions)") + for loc in locations[1:min(5, length(locations))] + println(" $loc") + end + if length(locations) > 5 + println(" ... and $(length(locations) - 5) more") + end + println() + end + end + + println("="^80) + println("SUMMARY") + println("="^80) + println() + @printf "Total unique functions in vendor: %d\n" length(vendor_functions) + @printf "Total unique functions in src: %d\n" length(src_functions) + @printf "Functions defined multiple times: %d\n" collision_count + @printf "High-risk collisions: %d\n" length(high_risk) + println() + + if length(high_risk) > 0 + println("High-risk collisions to review:") + for fname in high_risk[1:min(10, length(high_risk))] + println(" - $fname") + end + if length(high_risk) > 10 + println(" ... and $(length(high_risk) - 10) more") + end + else + println("✅ No high-risk collisions detected!") + println(" Most collisions are likely dispatch specializations (safe)") + end + + println() + println("="^80) + println("RECOMMENDATION") + println("="^80) + println(""" + Most function "collisions" in FEM packages are actually safe: + - Different Problem{T} types dispatch correctly + - get_* and assemble_* are intentional overrides + + Review high-risk collisions manually before consolidation. + """) +end + +main() diff --git a/scripts/fix_vendor_element_types.py b/scripts/fix_vendor_element_types.py new file mode 100755 index 0000000..9dadf26 --- /dev/null +++ b/scripts/fix_vendor_element_types.py @@ -0,0 +1,120 @@ +#!/usr/bin/env python3 +""" +Fix Element type signatures in vendor packages +Converts Element{BasisType} → Element{M, BasisType} where M +""" +import re +import sys +from pathlib import Path + +# Element types to fix +ELEMENT_TYPES = [ + "Seg2", + "Seg3", + "Poi1", + "Tri3", + "Tri6", + "Quad4", + "Quad8", + "Quad9", + "Tet4", + "Tet10", + "Pyr5", + "Wedge6", + "Wedge15", + "Hex8", + "Hex20", + "Hex27", +] + + +def fix_element_types(filepath): + """Fix Element type signatures in a single file.""" + with open(filepath, "r") as f: + content = f.read() + + original = content + + # Pattern 1: Element{Type} in function parameters/returns + for etype in ELEMENT_TYPES: + # Fix ::Element{Type} + content = re.sub( + rf"::Element\{{{etype}\}}", rf"::Element{{M, {etype}}}", content + ) + + # Fix Vector{Element{Type}} + content = re.sub( + rf"::Vector\{{Element\{{M, {etype}\}}\}}", + rf"::Vector{{Element{{M, {etype}}}}}", + content, + ) + + # Pattern 2: Add "where M" to function signatures + lines = content.split("\n") + new_lines = [] + + for i, line in enumerate(lines): + # Skip if already has where M + if "where M" in line or "where {M" in line or "where E" in line: + new_lines.append(line) + continue + + # Check if this is a function signature with Element{M, + if re.search(r"function\s+.*Element\{M,", line): + # Check if function signature closes on this line + if ")" in line and "end" not in line: + # Add where M before any trailing comment + line = line.rstrip() + if not line.endswith(" where M"): + # Handle inline functions (one-liners) + if "=" in line and line.count(")") == line.count("("): + # function foo(...) = expr + line = re.sub(r"\)(\s*=)", r") where M\1", line) + else: + line = line + " where M" + + new_lines.append(line) + + content = "\n".join(new_lines) + + if content != original: + with open(filepath, "w") as f: + f.write(content) + return True + return False + + +def main(): + vendor_dir = Path("/home/juajukka/dev/JuliaFEM.jl/vendor") + fixed_count = 0 + file_count = 0 + + print("Fixing Element type signatures in vendor packages...") + print("=" * 60) + + # Find all .jl files in vendor subdirectories + for package_dir in vendor_dir.iterdir(): + if not package_dir.is_dir(): + continue + + src_dir = package_dir / "src" + if not src_dir.exists(): + continue + + for jl_file in src_dir.rglob("*.jl"): + file_count += 1 + if fix_element_types(jl_file): + print(f"✓ Fixed: {jl_file.relative_to(vendor_dir)}") + fixed_count += 1 + + print("=" * 60) + print(f"Processed {file_count} files") + print(f"Fixed {fixed_count} files") + + if fixed_count > 0: + print("\n⚠️ IMPORTANT: Review changes before committing!") + print("Run: git diff vendor/") + + +if __name__ == "__main__": + main() diff --git a/scripts/generate_lagrange_basis.jl b/scripts/generate_lagrange_basis.jl new file mode 100755 index 0000000..2059272 --- /dev/null +++ b/scripts/generate_lagrange_basis.jl @@ -0,0 +1,711 @@ +#!/usr/bin/env julia + +# ============================================================================== +# LAGRANGE BASIS FUNCTION GENERATION SCRIPT +# ============================================================================== +# +# This script generates pre-computed Lagrange basis functions for all standard +# finite element types and writes them to src/basis/lagrange_generated.jl +# +# USAGE: +# cd /path/to/JuliaFEM.jl +# julia --project=. scripts/generate_lagrange_basis.jl +# +# OUTPUT: +# src/basis/lagrange_generated.jl (commit this to git!) +# +# WHEN TO RUN: +# - Adding new element types +# - Fixing bugs in generation logic +# - Changing polynomial ansatz strategy +# - After modifying lagrange_generator.jl +# +# THEORY: +# See docs/theory/lagrange_basis_functions.md for full mathematical details +# +# ============================================================================== + +using Pkg +Pkg.activate(".") + +using Dates # For timestamp in generated file header + +println("="^80) +println("LAGRANGE BASIS FUNCTION GENERATOR") +println("="^80) +println() + +# Load the symbolic generation engine +println("Loading symbolic generator...") +include("../src/basis/lagrange_generator.jl") +println("✓ Generator loaded") +println() + +# ============================================================================== +# ELEMENT CATALOG +# ============================================================================== +# +# For each standard Lagrange element, define: +# 1. Name (e.g., "Seg2") +# 2. Description +# 3. Node coordinates in reference element +# 4. Polynomial ansatz (monomial basis) +# +# Reference element conventions: +# - 1D (Seg): u ∈ [-1, 1] +# - 2D (Tri): (u, v) where u,v ≥ 0, u+v ≤ 1 +# - 2D (Quad): (u, v) ∈ [-1, 1]² +# - 3D (Tet): (u, v, w) where u,v,w ≥ 0, u+v+w ≤ 1 +# - 3D (Hex): (u, v, w) ∈ [-1, 1]³ +# +# NOTE: Variable names are u, v, w (not ξ, η, ζ) +# +# ============================================================================== + +elements = [] + +# ------------------------------------------------------------------------------ +# 1D ELEMENTS +# ------------------------------------------------------------------------------ + +# Seg2: 2-node linear segment +push!(elements, ( + name = "Seg2", + description = "2-node linear segment element", + coordinates = [ + [-1.0], # Node 1 + [ 1.0] # Node 2 + ], + ansatz = [ + :(1), # 1 + :(u) # u + ] +)) + +# Seg3: 3-node quadratic segment +push!(elements, ( + name = "Seg3", + description = "3-node quadratic segment element", + coordinates = [ + [-1.0], # Node 1 + [ 1.0], # Node 2 + [ 0.0] # Node 3 (midpoint) + ], + ansatz = [ + :(1), # 1 + :(u), # u + :(u^2) # u² + ] +)) + +# ------------------------------------------------------------------------------ +# 2D TRIANGULAR ELEMENTS +# ------------------------------------------------------------------------------ + +# Tri3: 3-node linear triangle +push!(elements, ( + name = "Tri3", + description = "3-node linear triangular element", + coordinates = [ + [0.0, 0.0], # Node 1 + [1.0, 0.0], # Node 2 + [0.0, 1.0] # Node 3 + ], + ansatz = [ + :(1), # 1 + :(u), # u + :(v) # v + ] +)) + +# Tri6: 6-node quadratic triangle +push!(elements, ( + name = "Tri6", + description = "6-node quadratic triangular element", + coordinates = [ + [0.0, 0.0], # Node 1 + [1.0, 0.0], # Node 2 + [0.0, 1.0], # Node 3 + [0.5, 0.0], # Node 4 (edge 1-2) + [0.5, 0.5], # Node 5 (edge 2-3) + [0.0, 0.5] # Node 6 (edge 3-1) + ], + ansatz = [ + :(1), # 1 + :(u), # u + :(v), # v + :(u^2), # u² + :(u*v), # uv + :(v^2) # v² + ] +)) + +# Tri6: 6-node quadratic triangle +push!(elements, ( + name="Tri6", + description="6-node quadratic triangular element", + coordinates=[ + [0.0, 0.0], # Node 1 + [1.0, 0.0], # Node 2 + [0.0, 1.0], # Node 3 + [0.5, 0.0], # Node 4 (edge 1-2) + [0.5, 0.5], # Node 5 (edge 2-3) + [0.0, 0.5] # Node 6 (edge 3-1) + ], + ansatz=[ + :(1), # 1 + :( u), # ξ + :( v), # η + :(u^2), # ξ² + :(ξ * η), # ξη + :(v^2) # η² + ] +)) + +# ------------------------------------------------------------------------------ +# 2D QUADRILATERAL ELEMENTS +# ------------------------------------------------------------------------------ + +# Quad4: 4-node bilinear quadrilateral +# Quad4: 4-node bilinear quadrilateral +push!(elements, ( + name = "Quad4", + description = "4-node bilinear quadrilateral element", + coordinates = [ + [-1.0, -1.0], # Node 1 + [ 1.0, -1.0], # Node 2 + [ 1.0, 1.0], # Node 3 + [-1.0, 1.0] # Node 4 + ], + ansatz = [ + :(1), # 1 + :(u), # u + :(v), # v + :(u*v) # uv + ] +)) + +# Quad8: 8-node serendipity quadrilateral (quadratic edges, no center node) +push!(elements, ( + name = "Quad8", + description = "8-node serendipity quadrilateral element", + coordinates = [ + [-1.0, -1.0], # Node 1 + [ 1.0, -1.0], # Node 2 + [ 1.0, 1.0], # Node 3 + [-1.0, 1.0], # Node 4 + [ 0.0, -1.0], # Node 5 (edge 1-2) + [ 1.0, 0.0], # Node 6 (edge 2-3) + [ 0.0, 1.0], # Node 7 (edge 3-4) + [-1.0, 0.0] # Node 8 (edge 4-1) + ], + ansatz = [ + :(1), # 1 + :(u), # u + :(v), # v + :(u^2), # u² + :(u*v), # uv + :(v^2), # v² + :(u^2*v), # u²v + :(u*v^2) # uv² + ] +)) + +# Quad9: 9-node biquadratic quadrilateral (complete quadratic) +push!(elements, ( + name = "Quad9", + description = "9-node biquadratic quadrilateral element", + coordinates = [ + [-1.0, -1.0], # Node 1 + [ 1.0, -1.0], # Node 2 + [ 1.0, 1.0], # Node 3 + [-1.0, 1.0], # Node 4 + [ 0.0, -1.0], # Node 5 (edge 1-2) + [ 1.0, 0.0], # Node 6 (edge 2-3) + [ 0.0, 1.0], # Node 7 (edge 3-4) + [-1.0, 0.0], # Node 8 (edge 4-1) + [ 0.0, 0.0] # Node 9 (center) + ], + ansatz = [ + :(1), # 1 + :(u), # u + :(v), # v + :(u^2), # u² + :(u*v), # uv + :(v^2), # v² + :(u^2*v), # u²v + :(u*v^2), # uv² + :(u^2*v^2) # u²v² + ] +)) + +# Quad8: 8-node serendipity quadrilateral (quadratic edges, no center node) +push!(elements, ( + name="Quad8", + description="8-node serendipity quadrilateral element", + coordinates=[ + [-1.0, -1.0], # Node 1 + [1.0, -1.0], # Node 2 + [1.0, 1.0], # Node 3 + [-1.0, 1.0], # Node 4 + [0.0, -1.0], # Node 5 (edge 1-2) + [1.0, 0.0], # Node 6 (edge 2-3) + [0.0, 1.0], # Node 7 (edge 3-4) + [-1.0, 0.0] # Node 8 (edge 4-1) + ], + ansatz=[ + :(1), # 1 + :( u), # ξ + :( v), # η + :(u^2), # ξ² + :(ξ * η), # ξη + :(v^2), # η² + :(u^2 * η), # ξ²η + :(ξ * v^2) # ξη² + ] +)) + +# Quad9: 9-node biquadratic quadrilateral (complete quadratic) +push!(elements, ( + name="Quad9", + description="9-node biquadratic quadrilateral element", + coordinates=[ + [-1.0, -1.0], # Node 1 + [1.0, -1.0], # Node 2 + [1.0, 1.0], # Node 3 + [-1.0, 1.0], # Node 4 + [0.0, -1.0], # Node 5 (edge 1-2) + [1.0, 0.0], # Node 6 (edge 2-3) + [0.0, 1.0], # Node 7 (edge 3-4) + [-1.0, 0.0], # Node 8 (edge 4-1) + [0.0, 0.0] # Node 9 (center) + ], + ansatz=[ + :(1), # 1 + :( u), # ξ + :( v), # η + :(u^2), # ξ² + :(ξ * η), # ξη + :(v^2), # η² + :(u^2 * η), # ξ²η + :(ξ * v^2), # ξη² + :(u^2 * v^2) # ξ²η² + ] +)) + +# ------------------------------------------------------------------------------ +# 3D TETRAHEDRAL ELEMENTS +# ------------------------------------------------------------------------------ + +# Tet4: 4-node linear tetrahedron +# Tet4: 4-node linear tetrahedron +push!(elements, ( + name = "Tet4", + description = "4-node linear tetrahedral element", + coordinates = [ + [0.0, 0.0, 0.0], # Node 1 + [1.0, 0.0, 0.0], # Node 2 + [0.0, 1.0, 0.0], # Node 3 + [0.0, 0.0, 1.0] # Node 4 + ], + ansatz = [ + :(1), # 1 + :(u), # u + :(v), # v + :(w) # w + ] +)) + +# Tet10: 10-node quadratic tetrahedron +push!(elements, ( + name = "Tet10", + description = "10-node quadratic tetrahedral element", + coordinates = [ + [0.0, 0.0, 0.0], # Node 1 + [1.0, 0.0, 0.0], # Node 2 + [0.0, 1.0, 0.0], # Node 3 + [0.0, 0.0, 1.0], # Node 4 + [0.5, 0.0, 0.0], # Node 5 (edge 1-2) + [0.5, 0.5, 0.0], # Node 6 (edge 2-3) + [0.0, 0.5, 0.0], # Node 7 (edge 3-1) + [0.0, 0.0, 0.5], # Node 8 (edge 1-4) + [0.5, 0.0, 0.5], # Node 9 (edge 2-4) + [0.0, 0.5, 0.5] # Node 10 (edge 3-4) + ], + ansatz = [ + :(1), # 1 + :(u), # u + :(v), # v + :(w), # w + :(u^2), # u² + :(u*v), # uv + :(u*w), # uw + :(v^2), # v² + :(v*w), # vw + :(w^2) # w² + ] +)) + +# Tet10: 10-node quadratic tetrahedron +push!(elements, ( + name="Tet10", + description="10-node quadratic tetrahedral element", + coordinates=[ + [0.0, 0.0, 0.0], # Node 1 + [1.0, 0.0, 0.0], # Node 2 + [0.0, 1.0, 0.0], # Node 3 + [0.0, 0.0, 1.0], # Node 4 + [0.5, 0.0, 0.0], # Node 5 (edge 1-2) + [0.5, 0.5, 0.0], # Node 6 (edge 2-3) + [0.0, 0.5, 0.0], # Node 7 (edge 3-1) + [0.0, 0.0, 0.5], # Node 8 (edge 1-4) + [0.5, 0.0, 0.5], # Node 9 (edge 2-4) + [0.0, 0.5, 0.5] # Node 10 (edge 3-4) + ], + ansatz=[ + :(1), # 1 + :( u), # ξ + :( v), # η + :( w), # ζ + :(u^2), # ξ² + :(ξ * η), # ξη + :(ξ * ζ), # ξζ + :(v^2), # η² + :(η * ζ), # ηζ + :(w^2) # ζ² + ] +)) + +# ------------------------------------------------------------------------------ +# 3D HEXAHEDRAL ELEMENTS +# ------------------------------------------------------------------------------ + +# Hex8: 8-node trilinear hexahedron (brick) +push!(elements, ( + name="Hex8", + description="8-node trilinear hexahedral element", + coordinates=[ + [-1.0, -1.0, -1.0], # Node 1 + [1.0, -1.0, -1.0], # Node 2 + [1.0, 1.0, -1.0], # Node 3 + [-1.0, 1.0, -1.0], # Node 4 + [-1.0, -1.0, 1.0], # Node 5 + [1.0, -1.0, 1.0], # Node 6 + [1.0, 1.0, 1.0], # Node 7 + [-1.0, 1.0, 1.0] # Node 8 + ], + ansatz=[ + :(1), # 1 + :( u), # ξ + :( v), # η + :( w), # ζ + :(ξ * η), # ξη + :(ξ * ζ), # ξζ + :(η * ζ), # ηζ + :(ξ * η * ζ) # ξηζ + ] +)) + +# Hex20: 20-node serendipity hexahedron (quadratic edges, no face/center nodes) +push!(elements, ( + name="Hex20", + description="20-node serendipity hexahedral element", + coordinates=[ + # Corner nodes + [-1.0, -1.0, -1.0], # 1 + [1.0, -1.0, -1.0], # 2 + [1.0, 1.0, -1.0], # 3 + [-1.0, 1.0, -1.0], # 4 + [-1.0, -1.0, 1.0], # 5 + [1.0, -1.0, 1.0], # 6 + [1.0, 1.0, 1.0], # 7 + [-1.0, 1.0, 1.0], # 8 + # Edge nodes (bottom face) + [0.0, -1.0, -1.0], # 9 + [1.0, 0.0, -1.0], # 10 + [0.0, 1.0, -1.0], # 11 + [-1.0, 0.0, -1.0], # 12 + # Edge nodes (top face) + [0.0, -1.0, 1.0], # 13 + [1.0, 0.0, 1.0], # 14 + [0.0, 1.0, 1.0], # 15 + [-1.0, 0.0, 1.0], # 16 + # Edge nodes (vertical) + [-1.0, -1.0, 0.0], # 17 + [1.0, -1.0, 0.0], # 18 + [1.0, 1.0, 0.0], # 19 + [-1.0, 1.0, 0.0] # 20 + ], + ansatz=[ + :(1), # 1 + :( u), :( v), :( w), # linear + :(u^2), :(v^2), :(w^2), # pure quadratic + :(ξ * η), :(ξ * ζ), :(η * ζ), # bilinear + :(u^2 * η), :(u^2 * ζ), :(v^2 * ξ), :(v^2 * ζ), :(w^2 * ξ), :(w^2 * η), # quadratic-linear + :(ξ * η * ζ), # trilinear + :(u^2 * η * ζ), :(ξ * v^2 * ζ), :(ξ * η * w^2) # quadratic-bilinear + ] +)) + +# Hex27: 27-node triquadratic hexahedron (complete quadratic) +push!(elements, ( + name="Hex27", + description="27-node triquadratic hexahedral element", + coordinates=[ + # Corner nodes + [-1.0, -1.0, -1.0], # 1 + [1.0, -1.0, -1.0], # 2 + [1.0, 1.0, -1.0], # 3 + [-1.0, 1.0, -1.0], # 4 + [-1.0, -1.0, 1.0], # 5 + [1.0, -1.0, 1.0], # 6 + [1.0, 1.0, 1.0], # 7 + [-1.0, 1.0, 1.0], # 8 + # Edge nodes (bottom face) + [0.0, -1.0, -1.0], # 9 + [1.0, 0.0, -1.0], # 10 + [0.0, 1.0, -1.0], # 11 + [-1.0, 0.0, -1.0], # 12 + # Edge nodes (top face) + [0.0, -1.0, 1.0], # 13 + [1.0, 0.0, 1.0], # 14 + [0.0, 1.0, 1.0], # 15 + [-1.0, 0.0, 1.0], # 16 + # Edge nodes (vertical) + [-1.0, -1.0, 0.0], # 17 + [1.0, -1.0, 0.0], # 18 + [1.0, 1.0, 0.0], # 19 + [-1.0, 1.0, 0.0], # 20 + # Face center nodes + [0.0, 0.0, -1.0], # 21 (bottom) + [0.0, 0.0, 1.0], # 22 (top) + [0.0, -1.0, 0.0], # 23 (front) + [1.0, 0.0, 0.0], # 24 (right) + [0.0, 1.0, 0.0], # 25 (back) + [-1.0, 0.0, 0.0], # 26 (left) + # Volume center node + [0.0, 0.0, 0.0] # 27 (center) + ], + ansatz=[ + :(1), # 1 + :( u), :( v), :( w), # linear (3) + :(u^2), :(v^2), :(w^2), # pure quadratic (3) + :(ξ * η), :(ξ * ζ), :(η * ζ), # bilinear (3) + :(u^2 * η), :(u^2 * ζ), :(v^2 * ξ), :(v^2 * ζ), :(w^2 * ξ), :(w^2 * η), # quad-linear (6) + :(ξ * η * ζ), # trilinear (1) + :(u^2 * η * ζ), :(ξ * v^2 * ζ), :(ξ * η * w^2), # quad-bilinear (3) + :(u^2 * v^2), :(u^2 * w^2), :(v^2 * w^2), # biquadratic (3) + :(u^2 * v^2 * ζ), :(u^2 * η * w^2), :(ξ * v^2 * w^2), # biquad-linear (3) + :(u^2 * v^2 * w^2) # triquadratic (1) + ] +)) + +# ------------------------------------------------------------------------------ +# 3D PYRAMID ELEMENTS +# ------------------------------------------------------------------------------ + +# Pyr5: 5-node linear pyramid +push!(elements, ( + name="Pyr5", + description="5-node linear pyramid element", + coordinates=[ + [-1.0, -1.0, 0.0], # Node 1 (base) + [1.0, -1.0, 0.0], # Node 2 (base) + [1.0, 1.0, 0.0], # Node 3 (base) + [-1.0, 1.0, 0.0], # Node 4 (base) + [0.0, 0.0, 1.0] # Node 5 (apex) + ], + ansatz=[ + :(1), # 1 + :( u), # ξ + :( v), # η + :( w), # ζ + :(ξ * η) # ξη (needed for symmetry) + ] +)) + +# ------------------------------------------------------------------------------ +# 3D WEDGE ELEMENTS (Prism/Pentahedron) +# ------------------------------------------------------------------------------ + +# Wedge6: 6-node linear wedge (triangular prism) +push!(elements, ( + name="Wedge6", + description="6-node linear wedge element (triangular prism)", + coordinates=[ + [0.0, 0.0, -1.0], # Node 1 (bottom triangle) + [1.0, 0.0, -1.0], # Node 2 (bottom triangle) + [0.0, 1.0, -1.0], # Node 3 (bottom triangle) + [0.0, 0.0, 1.0], # Node 4 (top triangle) + [1.0, 0.0, 1.0], # Node 5 (top triangle) + [0.0, 1.0, 1.0] # Node 6 (top triangle) + ], + ansatz=[ + :(1), # 1 + :( u), # ξ + :( v), # η + :( w), # ζ + :(ξ * ζ), # ξζ + :(η * ζ) # ηζ + ] +)) + +# Wedge15: 15-node quadratic wedge +push!(elements, ( + name="Wedge15", + description="15-node quadratic wedge element", + coordinates=[ + # Bottom triangle + [0.0, 0.0, -1.0], # 1 + [1.0, 0.0, -1.0], # 2 + [0.0, 1.0, -1.0], # 3 + # Top triangle + [0.0, 0.0, 1.0], # 4 + [1.0, 0.0, 1.0], # 5 + [0.0, 1.0, 1.0], # 6 + # Mid-edge nodes (bottom triangle) + [0.5, 0.0, -1.0], # 7 + [0.5, 0.5, -1.0], # 8 + [0.0, 0.5, -1.0], # 9 + # Mid-edge nodes (top triangle) + [0.5, 0.0, 1.0], # 10 + [0.5, 0.5, 1.0], # 11 + [0.0, 0.5, 1.0], # 12 + # Mid-edge nodes (vertical) + [0.0, 0.0, 0.0], # 13 + [1.0, 0.0, 0.0], # 14 + [0.0, 1.0, 0.0] # 15 + ], + ansatz=[ + :(1), # 1 + :( u), :( v), :( w), # linear (3) + :(u^2), :(v^2), :(w^2), # pure quadratic (3) + :(ξ * η), # triangle bilinear (1) + :(ξ * ζ), :(η * ζ), # prism bilinear (2) + :(u^2 * ζ), :(v^2 * ζ), :(ξ * η * ζ), # mixed (3) + :(ξ * w^2), :(η * w^2) # mixed (2) + ] +)) + +println("Element catalog loaded: $(length(elements)) element types") +println() + +# ============================================================================== +# GENERATION LOOP +# ============================================================================== + +println("Generating basis functions...") +println("─"^80) + +output = IOBuffer() + +# File header +println(output, "# This file is a part of JuliaFEM.") +println(output, "# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE") +println(output) +println(output, "# ============================================================================") +println(output, "# AUTO-GENERATED LAGRANGE BASIS FUNCTIONS") +println(output, "# ============================================================================") +println(output, "#") +println(output, "# WARNING: DO NOT EDIT THIS FILE MANUALLY!") +println(output, "#") +println(output, "# This file was automatically generated by:") +println(output, "# scripts/generate_lagrange_basis.jl") +println(output, "#") +println(output, "# To regenerate (e.g., after adding new element types):") +println(output, "# cd /path/to/JuliaFEM.jl") +println(output, "# julia --project=. scripts/generate_lagrange_basis.jl") +println(output, "#") +println(output, "# Theory:") +println(output, "# See docs/theory/lagrange_basis_functions.md") +println(output, "#") +println(output, "# Generator:") +println(output, "# src/basis/lagrange_generator.jl (symbolic engine)") +println(output, "#") +println(output, "# Generated: $(Dates.format(now(), "yyyy-mm-dd HH:MM:SS"))") +println(output, "# ============================================================================") +println(output) + +for (i, elem) in enumerate(elements) + println("[$i/$(length(elements))] Generating $(elem.name)...") + + try + # Convert coordinates to proper format (tuples, not vectors) + coords = [tuple(coord...) for coord in elem.coordinates] + + # Build polynomial ansatz as single expression: term1 + term2 + ... + # The ansatz is a vector of terms [:1, :ξ, :η, ...] → combine with + + if length(elem.ansatz) == 1 + ansatz_expr = elem.ansatz[1] + else + ansatz_expr = Expr(:call, :+, elem.ansatz...) + end + + # Generate basis code using symbolic engine + # Returns an Expr (quoted code) + basis_code_expr = create_basis( + Symbol(elem.name), + elem.description, + coords, + ansatz_expr + ) + + # Convert Expr to readable Julia code string + basis_code_str = string(basis_code_expr) + + # Pretty-print: The generated code is in a quote block, extract inner code + # Handle both `quote ... end` and `begin ... end` forms + basis_code_str = replace(basis_code_str, r"^(quote|begin)\s+" => "") + basis_code_str = replace(basis_code_str, r"\s*end$" => "") + + # Write to output with nice formatting + println(output, "# " * "─"^78) + println(output, "# $(elem.name): $(elem.description)") + println(output, "# " * "─"^78) + println(output) + println(output, basis_code_str) + println(output) + + catch e + println(" ⚠ Error generating $(elem.name):") + println(" $e") + if isa(e, ErrorException) + # Print stacktrace for debugging + for (exc, bt) in Base.catch_stack() + showerror(stdout, exc, bt) + println() + end + end + println(" Skipping...") + end +end + +println("─"^80) +println() + +# ============================================================================== +# WRITE OUTPUT FILE +# ============================================================================== + +output_path = joinpath(@__DIR__, "..", "src", "basis", "lagrange_generated.jl") +println("Writing to: $output_path") + +output_content = String(take!(output)) +write(output_path, output_content) + +println("✓ Generation complete!") +println() +println("Generated $(length(elements)) element types:") +for elem in elements + println(" - $(elem.name): $(elem.description)") +end +println() +println("Next steps:") +println(" 1. Review: src/basis/lagrange_generated.jl") +println(" 2. Update module to include this file") +println(" 3. Remove __precompile__(false) from main module") +println(" 4. Test: julia --project=. -e 'using JuliaFEM'") +println(" 5. Run tests: julia --project=. -e 'using Pkg; Pkg.test()'") +println(" 6. Commit: git add src/basis/lagrange_generated.jl && git commit") +println() +println("="^80) diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 71a0992..32f9ddd 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -137,7 +137,9 @@ end include("basis/abstract.jl") include("basis/subs.jl") # Symbolic substitution (includes minimal simplify from SymDiff.jl) include("basis/vandermonde.jl") -include("basis/create_basis.jl") # Basis generation (includes minimal differentiate from SymDiff.jl) +# NOTE: lagrange_generator.jl is NOT included here - it's a tool, not runtime code! +# It's only loaded by scripts/generate_lagrange_basis.jl during pre-generation. +# The generated code is in lagrange_generated.jl (to be created). include("basis/lagrange_segments.jl") include("basis/lagrange_quadrangles.jl") include("basis/lagrange_triangles.jl") diff --git a/src/basis/create_basis.jl b/src/basis/lagrange_generator.jl similarity index 74% rename from src/basis/create_basis.jl rename to src/basis/lagrange_generator.jl index a6b3f19..8f7d02a 100644 --- a/src/basis/create_basis.jl +++ b/src/basis/lagrange_generator.jl @@ -1,7 +1,57 @@ # This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/jl/blob/master/LICENSE +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE -__precompile__(false) +# ============================================================================== +# LAGRANGE BASIS FUNCTION GENERATOR +# ============================================================================== +# +# This file contains the symbolic engine for generating Lagrange basis functions. +# It is NOT loaded at runtime - it's a TOOL used during development. +# +# PURPOSE: +# Generate pre-computed basis functions and derivatives for all standard +# Lagrange finite elements (Seg2, Tri3, Quad4, Tet10, Hex8, etc.) +# +# THEORY: +# Lagrange basis functions satisfy the Kronecker delta property: +# +# N_i(x_j) = δ_ij = { 1 if i = j +# { 0 if i ≠ j +# +# Given: +# - n nodes with coordinates {x₁, x₂, ..., xₙ} in reference element +# - Polynomial ansatz {p₁(x), p₂(x), ..., pₙ(x)} (complete to order k) +# +# We construct: N_i(x) = Σⱼ αᵢⱼ pⱼ(x) +# +# The Kronecker property gives: V α_i = e_i +# +# Where Vandermonde matrix: V_kj = pⱼ(x_k) +# +# Solving these n systems gives all basis functions explicitly. +# Then symbolic differentiation provides derivatives. +# +# USAGE: +# This file is loaded by scripts/generate_lagrange_basis.jl which: +# 1. Defines all standard element types (coords + polynomial ansatz) +# 2. Calls generate_lagrange_basis() for each +# 3. Writes clean Julia code to src/basis/lagrange_generated.jl +# +# WHY GENERATE ONCE? +# - Symbolic math is expensive (100+ ms per element type) +# - Generated code is constant (mathematics doesn't change!) +# - Pre-compilation is much faster +# - Generated code is readable and debuggable +# - Version control shows what changed +# +# SEE: +# - docs/theory/lagrange_basis_functions.md (mathematical explanation) +# - scripts/generate_lagrange_basis.jl (generation script) +# - src/basis/lagrange_generated.jl (output - do not edit manually!) +# +# ============================================================================== + +__precompile__(false) # This is a tool, not runtime code # Minimal symbolic differentiation for polynomial basis functions # Adapted from SymDiff.jl by Jukka Aho - zero dependencies!