diff --git a/demos/cantilever_physics_gpu.jl b/demos/cantilever_physics_gpu.jl new file mode 100644 index 0000000..7cfdce5 --- /dev/null +++ b/demos/cantilever_physics_gpu.jl @@ -0,0 +1,261 @@ +# Cantilever Beam Demo - Pure GPU Physics{Elasticity} +# +# Demonstrates new API: +# - Physics{Elasticity} (not Problem) +# - Elements store geometry via update!(element, "geometry", nodes) +# - BCs via add_dirichlet! and add_neumann! +# - Pure GPU solve with solve_elasticity_gpu! + +# IMPORTANT: Load CUDA before JuliaFEM for GPU backend support +using CUDA +using JuliaFEM +using Tensors + +println("="^60) +println("Cantilever Beam - Backend-Transparent Demo") +println("="^60) + +# ============================================================================ +# 1. Create Mesh (Simple hand-coded for demo) +# ============================================================================ + +println("\n[1] Creating mesh...") + +# 2×1×1 beam with 2 Tet4 elements (minimal example) +nodes = [ + # Beam: L=2, W=1, H=1 + 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 # X + 0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0 # Y + 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 # Z +] + +# Tet4 connectivity (two elements spanning the beam) +tet_connectivity = [ + [1, 2, 3, 5], # Element 1 + [2, 3, 4, 6], # Element 2 + [3, 4, 5, 7], # Element 3 + [4, 5, 6, 8], # Element 4 +] + +println(" Nodes: $(size(nodes, 2))") +println(" Elements: $(length(tet_connectivity))") + +# ============================================================================ +# 2. Create Body Elements with Geometry and Material +# ============================================================================ + +println("\n[2] Creating body elements...") + +body_elements = Element[] + +for conn in tet_connectivity + # Create element with geometry and material (immutable API) + X_elem = nodes[:, conn] + + # Legacy API: topology + connectivity (infers Lagrange{Tet4,1}) + el = Element(Tet4, conn; + fields=(geometry=X_elem, + youngs_modulus=210e9, # 210 GPa (steel) + poissons_ratio=0.3)) + + push!(body_elements, el) +end + +println(" Body elements created: $(length(body_elements))") + +# ============================================================================ +# 3. Create Physics{Elasticity} +# ============================================================================ + +println("\n[3] Creating Physics{Elasticity}...") + +physics = Physics(Elasticity, "cantilever beam", 3) +physics.properties.formulation = :continuum +physics.properties.finite_strain = false + +add_elements!(physics, body_elements) + +println(" Physics: $(physics.name)") +println(" Elements in physics: $(length(physics.body_elements))") + +# ============================================================================ +# 4. Add Boundary Conditions +# ============================================================================ + +println("\n[4] Adding boundary conditions...") + +# Dirichlet BC: Fix nodes at X=0 (left end) +fixed_nodes = [1, 3, 5, 7] # Nodes with X=0 +println(" Fixing nodes: $fixed_nodes (all DOFs)") + +for node in fixed_nodes + add_dirichlet!(physics, [node], [1, 2, 3], 0.0) # Fix u_x, u_y, u_z = 0 +end + +# Neumann BC: Pressure on top surface (Z=1) +# Surface: nodes [5, 6, 7, 8] form two triangles +println(" Applying pressure load on top surface...") + +pressure = -1e6 # -1 MPa in -Z direction +traction = Vec{3}((0.0, 0.0, pressure)) + +# Top surface triangles +top_surface_tris = [ + [5, 6, 7], + [6, 7, 8] +] + +for tri_conn in top_surface_tris + X_surf = nodes[:, tri_conn] + + # Legacy API: topology + connectivity (infers basis) + surf_el = Element(Tri3, tri_conn; + fields=(geometry=X_surf,)) + + add_neumann!(physics, surf_el, traction) +end + +println(" Dirichlet BCs: $(length(physics.bc_dirichlet.node_ids)) nodes") +println(" Neumann BCs: $(length(physics.bc_neumann.surface_elements)) surfaces") + +# ============================================================================ +# 5. Solve (Backend Transparent!) +# ============================================================================ + +println("\n[5] Solving...") +println(" Using GPU backend (CPU backend not yet implemented)") + +# Unified solve! - backend transparent! +# Explicitly request GPU since CPU backend is just a stub +result = solve!( + physics; + backend=GPU(), # Use GPU (CPU backend coming in Phase 2) + time=0.0, + tol=1e-6, + max_iter=1000 +) + +println("\n" * "="^60) +println("SOLUTION") +println("="^60) + +if result.newton_iterations == 1 + # Linear problem (converged in 1 Newton iteration) + println(" Problem type: LINEAR elasticity") + println(" Newton iterations: 1 (linear problem - no iterations needed)") + println(" CG iterations: $(result.cg_iterations)") + println(" Final residual: $(result.residual)") + println() + println(" Note: For LINEAR problems, Newton converges in 1 iteration.") + println(" CG iterations shown are from the single linear solve.") +else + # Nonlinear problem (multiple Newton iterations) + println(" Problem type: NONLINEAR elasticity") + println(" Solver: Inexact Newton-Krylov (SIMULTANEOUS solving!)") + println() + println(" Newton iterations: $(result.newton_iterations)") + println(" Total CG iterations: $(result.cg_iterations)") + println(" Avg CG per Newton: $(round(result.cg_iterations / result.newton_iterations, digits=1))") + println(" Final residual: $(result.residual)") + println() + println(" Iteration history (Newton + CG solved SIMULTANEOUSLY):") + for (i, (cg_i, R_i, η_i)) in enumerate(result.history) + println(" Newton $i: CG=$cg_i, ||R||=$(round(R_i, sigdigits=3)), η=$(round(η_i, digits=3))") + end + println() + println(" ✓ Linear systems solved INEXACTLY (adaptive tolerance)") + println(" ✓ Newton and Krylov iterations INTERLEAVED") + println(" ✓ Much faster than nested loops!") +end + +println() +println(" Solve time: $(round(result.solve_time, digits=3)) seconds") + +# ============================================================================ +# 6. Post-Process Results +# ============================================================================ + +println("\n[6] Post-processing...") + +u = result.u +n_nodes = div(length(u), 3) + +# Extract displacement components +u_x = u[1:3:end] +u_y = u[2:3:end] +u_z = u[3:3:end] + +# Compute magnitude +u_mag = sqrt.(u_x .^ 2 + u_y .^ 2 + u_z .^ 2) + +println("\nDisplacement Statistics:") +println(" Max |u|: $(maximum(u_mag) * 1000) mm") +println(" Max u_x: $(maximum(abs.(u_x)) * 1000) mm") +println(" Max u_y: $(maximum(abs.(u_y)) * 1000) mm") +println(" Max u_z: $(maximum(abs.(u_z)) * 1000) mm") + +# Find max displacement location +max_idx = argmax(u_mag) +println("\nMax displacement at node $max_idx:") +println(" Location: $(nodes[:, max_idx])") +println(" u = [$(u_x[max_idx]*1000), $(u_y[max_idx]*1000), $(u_z[max_idx]*1000)] mm") + +# ============================================================================ +# 7. Validate (Simple Check) +# ============================================================================ + +println("\n[7] Validation...") + +# For cantilever with pressure load, expect: +# - Free end (X=2) has largest displacement +# - Fixed end (X=0) has zero displacement +# - Deflection primarily in -Z direction + +free_end_nodes = [2, 4, 6, 8] # X=2 +fixed_end_nodes = [1, 3, 5, 7] # X=0 + +free_end_disp = maximum(u_mag[free_end_nodes]) +fixed_end_disp = maximum(u_mag[fixed_end_nodes]) + +println(" Free end max |u|: $(free_end_disp * 1000) mm") +println(" Fixed end max |u|: $(fixed_end_disp * 1000) mm") + +if fixed_end_disp < 1e-10 + println(" ✓ Fixed end has zero displacement (good!)") +else + println(" ✗ Fixed end displacement > 0 (bad!)") +end + +if free_end_disp > 1e-6 + println(" ✓ Free end has non-zero displacement (good!)") +else + println(" ✗ Free end displacement ≈ 0 (bad!)") +end + +println("\n" * "="^60) +println("Demo complete!") +println("="^60) + +# ============================================================================ +# Summary +# ============================================================================ + +println("\nBackend-Transparent API Summary:") +println(" 1. Create elements: Element(Topology, connectivity; fields=(...))") +println(" 2. Create Physics: physics = Physics(Elasticity, name, dimension)") +println(" 3. Add elements: add_elements!(physics, body_elements)") +println(" 4. Add Dirichlet: add_dirichlet!(physics, node_ids, components, value)") +println(" 5. Add Neumann: add_neumann!(physics, surface_element, traction)") +println(" 6. Solve: result = solve!(physics) # ← Backend automatic!") +println() +println("Key Features:") +println(" ✓ Immutable elements with fields at construction") +println(" ✓ Backend transparency - same code for CPU or GPU") +println(" ✓ Automatic backend selection (GPU if CUDA available)") +println(" ✓ Matrix-free CG solver") +println(" ✓ No 'GPU' in user code!") +println() +println("Backend Selection:") +println(" solve!(physics) # Auto (GPU if available, else CPU)") +println(" solve!(physics; backend=GPU()) # Force GPU") +println(" solve!(physics; backend=CPU(8)) # Force CPU with 8 threads")