diff --git a/src/assemblers/nodal_based.jl b/src/assemblers/nodal_based.jl deleted file mode 100644 index 87800ea..0000000 --- a/src/assemblers/nodal_based.jl +++ /dev/null @@ -1,177 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -""" -Nodal-based assembly (future implementation). - -Node-by-node assembly using inverse connectivity (node-to-elements map). -Natural for GPU parallelization (one thread per node). - -**Performance**: Expected 2-10x speedup on GPU for large problems (> 100k nodes). -**Best for**: GPU acceleration, very large problems. -**Status**: Planned for future implementation. -""" - -using SparseArrays - -""" - assemble!( - cache::NodalCache, - assembler::NodalAssembler, - kernel::AbstractKernel, - mesh::AbstractMesh - ) -> Nothing - -Assemble global system using nodal traversal **in-place, zero allocations**. - -**Status**: Not yet implemented. Placeholder for future GPU-based assembly. - -# Algorithm (Planned) - -1. Reset cache -2. Loop over nodes (parallelizable on GPU): - a. Get all elements touching this node - b. For each touching element: - - Compute full element stiffness (or use cached value) - - Extract only rows/columns for this node - - Accumulate to global system (atomic add on GPU) - -# GPU Parallelization - -Each node processed by one GPU thread: -```cuda -__global__ void assemble_nodal(nodes, elements, K, f) { - int node_id = blockIdx.x * blockDim.x + threadIdx.x; - if (node_id >= nnodes) return; - - // Get touching elements - for (elem in touching_elements[node_id]) { - // Compute element contribution for this node - // Atomic add to K, f - } -} -``` - -# Arguments -- `cache`: Pre-allocated nodal cache -- `assembler`: Nodal assembler -- `kernel`: Domain kernel -- `mesh`: Finite element mesh - -# Zero-Allocation Guarantee - -All arrays pre-allocated. Atomic operations on GPU ensure thread-safety. -""" -function assemble!( - cache::NodalCache, - assembler::NodalAssembler, - kernel::AbstractKernel, - mesh::AbstractMesh -) - error("NodalAssembler not yet implemented. Use COOAssembler or CSCAssembler.") - - # Planned implementation: - # 1. Reset cache - # 2. Loop over nodes - # 3. For each node, get touching elements from cache.node_to_elements - # 4. Accumulate contributions from all touching elements - # 5. Write to global K, f (atomic on GPU) - - return nothing -end - -# ============================================================================ -# HELPER FUNCTIONS (for future implementation) -# ============================================================================ - -""" - create_cache(assembler::NodalAssembler, mesh::AbstractMesh, kernel::AbstractKernel) -> NodalCache - -Create pre-allocated cache for nodal assembly. - -Builds node-to-elements map (inverse connectivity) for efficient nodal traversal. - -# Arguments -- `assembler`: Nodal assembler -- `mesh`: Finite element mesh -- `kernel`: Domain kernel - -# Returns -- Pre-allocated nodal cache with inverse connectivity -""" -function create_cache(assembler::NodalAssembler, mesh::AbstractMesh, kernel::AbstractKernel) - return NodalCache(mesh, kernel) -end - -""" - compute_node_contributions!( - node_cache::NodeCache, - element_cache::ElementCache, - node_id::Int, - kernel::AbstractKernel, - mesh::AbstractMesh, - node_to_elements::NodeToElementsMap - ) -> Nothing - -Compute contributions to this node from all touching elements **in-place**. - -# Algorithm - -1. Get all elements touching this node -2. For each element: - a. Compute full element stiffness (using `compute_element_stiffness!`) - b. Find local node index within element - c. Extract rows/columns corresponding to this node's DOFs - d. Accumulate to node contributions - -# Arguments -- `node_cache`: Node workspace (output) -- `element_cache`: Element workspace (for kernel calls) -- `node_id`: Node index -- `kernel`: Domain kernel -- `mesh`: Finite element mesh -- `node_to_elements`: Inverse connectivity - -# Returns - -Nothing - writes to `node_cache` in-place. -""" -function compute_node_contributions!( - node_cache::NodeCache, - element_cache::ElementCache, - node_id::Int, - kernel::AbstractKernel, - mesh::AbstractMesh, - node_to_elements::NodeToElementsMap -) - error("compute_node_contributions! not yet implemented") - - # Planned implementation: - # 1. Get touching elements: get_node_spider(node_to_elements, node_id) - # 2. For each element: - # - Compute element stiffness: compute_element_stiffness!(element_cache, ...) - # - Find local node index in element - # - Extract node DOF rows/columns from Ke, fe - # - Accumulate to node_cache - # 3. Return node_cache (contains all contributions for this node) - - return nothing -end - -# ============================================================================ -# GPU KERNEL STUBS (for future CUDA implementation) -# ============================================================================ - -# """ -# assemble_nodal_gpu!(K, f, mesh, kernel, node_to_elements) -# -# GPU kernel for nodal assembly. -# -# Launches one thread per node. Each thread: -# 1. Gets touching elements for its node -# 2. Computes contributions from all elements -# 3. Atomically adds to global K, f -# -# Requires CUDA.jl or similar GPU framework. -# """ -# function assemble_nodal_gpu! end diff --git a/src/assemblers/nodal_cache.jl b/src/assemblers/nodal_cache.jl deleted file mode 100644 index aa45f72..0000000 --- a/src/assemblers/nodal_cache.jl +++ /dev/null @@ -1,125 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -""" -Nodal cache for node-based assembly. - -Nodal assembly iterates over nodes rather than elements. Each node assembles -contributions from all touching elements. This approach has advantages for: -- Contact mechanics (contact is inherently nodal) -- Domain decomposition (clear node ownership) -- Matrix-free operations (natural node-based matvec) -- Adaptive refinement (local node operations) - -# Structure -Pre-builds node-to-elements map (inverse connectivity) once. During assembly, -each node visits all touching elements and accumulates contributions. - -# Performance -Similar to CSC for standard problems. Better locality for nodal operations -like contact and matrix-free solvers. - -# Use case -Best for problems with nodal phenomena (contact, nodal plasticity) or -matrix-free iterative solvers. -""" - -using SparseArrays - -""" - NodalCache <: AbstractAssemblerCache - -Cache for nodal-based assembly. - -Pre-allocates sparse matrix, force vector, and node-to-elements map. -Each node assembles contributions from all touching elements. - -# Fields -- `K::SparseMatrixCSC{Float64,Int}`: Sparse matrix -- `f::Vector{Float64}`: Global force vector -- `node_cache::NodeCache`: Per-node workspace -- `element_cache::ElementCache`: Per-element workspace (for kernel calls) -- `node_to_elements::NodeToElementsMap`: Inverse connectivity - -# Zero-Allocation Usage - -```julia -cache = NodalCache(mesh, kernel) -fill!(cache) -assemble!(cache, assembler, kernel, mesh) # No allocations -K, f = extract_system(cache) -``` -""" -mutable struct NodalCache <: AbstractAssemblerCache - K::SparseMatrixCSC{Float64,Int} # Sparse matrix - f::Vector{Float64} # Force vector - node_cache::NodeCache # Node workspace - element_cache::ElementCache # Element workspace - node_to_elements::NodeToElementsMap # Inverse connectivity -end - -""" - NodalCache(mesh::AbstractMesh, kernel::AbstractKernel) -> NodalCache - -Create pre-allocated nodal cache. - -Builds node-to-elements map (inverse connectivity) for efficient nodal traversal. - -# Arguments -- `mesh`: Finite element mesh -- `kernel`: Domain kernel defining DOF structure - -# Returns -- `NodalCache` with pre-allocated workspace and inverse connectivity -""" -function NodalCache(mesh::AbstractMesh, kernel::AbstractKernel) - ndofs_per_node = dofs_per_node(kernel) - nnodes_mesh = nnodes_total(mesh) - ndofs = nnodes_mesh * ndofs_per_node - - # Build sparsity pattern (same as CSC) - K = build_sparsity_pattern(mesh, kernel) - f = zeros(Float64, ndofs) - - # Create caches - node_cache = create_node_cache(mesh, kernel) - element_cache = create_element_cache(mesh, kernel) - - # Build inverse connectivity - node_to_elements = NodeToElementsMap(mesh) - - return NodalCache(K, f, node_cache, element_cache, node_to_elements) -end - -""" - reset!(cache::NodalCache) - -Reset nodal cache for new assembly. - -Zeros out matrix values and force vector. -**Zero allocations** - reuses existing arrays. -""" -function reset!(cache::NodalCache) - fill!(cache.K.nzval, 0.0) - fill!(cache.f, 0.0) - return nothing -end - -""" - extract_system(cache::NodalCache) -> (K, f) - -Extract global system from nodal cache. - -Returns references to sparse matrix and force vector. -**Zero allocations** - no copying. - -# Arguments -- `cache`: Nodal cache after assembly - -# Returns -- `K`: Sparse matrix (reference, no copy) -- `f`: Force vector (reference, no copy) -""" -function extract_system(cache::NodalCache) - return cache.K, cache.f -end diff --git a/src/assemblers/node_based_coo.jl b/src/assemblers/node_based_coo.jl deleted file mode 100644 index 62557bf..0000000 --- a/src/assemblers/node_based_coo.jl +++ /dev/null @@ -1,484 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -""" -Node-based COO assembly using block integration. - -**NODAL ASSEMBLY PARADIGM**: Loop over nodes, not elements! - -Each node: -1. Finds all elements touching it (via inverse connectivity) -2. For each touching element: - - Prepares element geometry once (PreparedElement) - - Computes only needed 3×3 blocks (compute_block!) -3. Scatters blocks to COO triplets - -# Key Differences from Element-Based Assembly - -**Element-Based (traditional):** -```julia -for element in elements - K_e = compute_element_stiffness(element) # Full N×N matrix of 3×3 blocks - scatter(K_e) # Scatter all entries -end -``` - -**Node-Based (this file):** -```julia -for node_i in nodes - for element in elements_touching(node_i) - prepared = prepare_element(element) # Geometry preprocessing - for node_j in element.nodes - K_ij = compute_block!(prepared, i, j) # Single 3×3 block - scatter(K_ij, i, j) # Scatter one block - end - end -end -``` - -# Advantages - -1. **GPU-friendly**: One thread per node, no race conditions -2. **Contact-ready**: Contact is naturally node-based -3. **Matrix-free ready**: Can compute K*v without forming K -4. **Cache-friendly**: Reuses PreparedElement for multiple blocks -5. **Adaptive-ready**: Easy to refine/coarsen at node level - -# Performance Expectations - -- **CPU Single-thread**: ~1.5-2x slower than element-based (more kernel calls) -- **CPU Multi-thread**: ~1.5-2x faster (better parallelization) -- **GPU**: ~10-50x faster (massive parallelization, no atomics needed) - -# References - -- Golden standard: `docs/src/book/multigpu_nodal_assembly.md` -- PreparedElement: `src/domains/continuum/integration.jl` -- Block kernel: `src/domains/continuum/kernel.jl` - -# Example - -```julia -# Setup -mesh = create_cantilever_mesh(50, 10, 10) -material = LinearElastic(E=210e9, ν=0.3) -kernel = ContinuumKernel( - ContinuumFormulation{FullThreeD}(), - material, - Displacement{3}() -) - -# Create node-based assembler and cache -assembler = NodeBasedCOOAssembler() -cache = create_cache(assembler, mesh, kernel) - -# Assemble (zero allocations after warmup!) -assemble!(cache, assembler, kernel, mesh) - -# Extract system -K, f = extract_system(cache) - -# Solve -apply_dirichlet_bcs!(K, f, kernel, mesh, bc_dirichlet) -u = K \\ f -``` -""" - -using SparseArrays -using Tensors - -""" - NodeBasedCOOCache - -Pre-allocated cache for node-based COO assembly. - -Similar to COOCache but includes inverse connectivity mapping. - -# Fields -- `I::Vector{Int}`: Row indices (COO format) -- `J::Vector{Int}`: Column indices (COO format) -- `V::Vector{Float64}`: Values (COO format) -- `f::Vector{Float64}`: Global force vector -- `counter::Ref{Int}`: Current triplet count -- `capacity::Int`: Maximum triplet capacity -- `node_to_elements::NodeToElementsMap`: Inverse connectivity -- `element_cache::ElementCache`: Cache for element operations -- `ndofs::Int`: Total DOFs in system -""" -struct NodeBasedCOOCache{T<:AbstractTopology,B<:AbstractBasis,IPS} - I::Vector{Int} - J::Vector{Int} - V::Vector{Float64} - f::Vector{Float64} - counter::Ref{Int} - capacity::Int - node_to_elements::NodeToElementsMap - element_cache::ElementCache{T,B,IPS} - ndofs::Int -end - -""" - NodeBasedCOOCache(mesh::AbstractMesh, kernel::ContinuumKernel) - -Create cache for node-based assembly. - -Builds inverse connectivity and allocates buffers. - -# Arguments -- `mesh`: Finite element mesh -- `kernel`: Continuum kernel - -# Returns -- Pre-allocated node-based COO cache -""" -function NodeBasedCOOCache(mesh::AbstractMesh, kernel::ContinuumKernel) - # Build inverse connectivity - node_to_elements = NodeToElementsMap(mesh.connectivity) - - # Estimate triplet count (same as element-based) - ndofs_per_node = dofs_per_node(kernel) - nnodes = length(mesh.nodes) - ndofs = ndofs_per_node * nnodes - - # Estimate: For each node, sum over touching elements - # Each element contributes N blocks (N = nodes per element) - # Each block = 3×3 = 9 triplets - avg_elements_per_node = node_to_elements.nelements / nnodes - N = length(first(mesh.connectivity)) # Nodes per element - estimated_triplets = Int(ceil(1.2 * nnodes * avg_elements_per_node * N * 9)) - - # Allocate triplet arrays - I = Vector{Int}(undef, estimated_triplets) - J = Vector{Int}(undef, estimated_triplets) - V = Vector{Float64}(undef, estimated_triplets) - f = zeros(Float64, ndofs) - counter = Ref(0) - - # Create element cache (for prepare_element! and compute_block!) - element_cache = ElementCache(mesh, kernel) - - return NodeBasedCOOCache(I, J, V, f, counter, estimated_triplets, - node_to_elements, element_cache, ndofs) -end - -""" - reset!(cache::NodeBasedCOOCache) - -Reset cache for new assembly (zero force vector, reset counter). - -Does NOT clear inverse connectivity (that's permanent structure). -""" -function reset!(cache::NodeBasedCOOCache) - fill!(cache.f, 0.0) - cache.counter[] = 0 - return nothing -end - -""" - assemble!( - cache::NodeBasedCOOCache, - assembler::NodeBasedCOOAssembler, - kernel::ContinuumKernel, - mesh::AbstractMesh - ) -> Nothing - -Assemble global system using **node-based traversal**. - -# Algorithm - -```julia -for node_i in 1:nnodes - # Get all elements touching this node - for elem_info in node_to_elements[node_i] - element_id = elem_info.element_id - local_i = elem_info.local_node_idx - - # Prepare element geometry ONCE - prepared = prepare_element!(cache.element_cache, kernel, element_id, mesh) - - # Compute blocks for all nodes in this element - for local_j in 1:N - global_j = connectivity[element_id][local_j] - - # Compute single 3×3 block - K_ij = compute_block!(prepared, kernel.material, local_i, local_j) - - # Scatter to triplets - scatter_block_to_triplets!(cache, K_ij, node_i, global_j) - end - end -end -``` - -# Key Operations - -1. **prepare_element!** - Precompute geometry (Jacobian, gradients) once per element -2. **compute_block!** - Compute single 3×3 stiffness block using prepared geometry -3. **scatter_block_to_triplets!** - Add 9 triplets (i,j,value) for 3×3 block - -# Zero-Allocation (After Warmup) - -All arrays pre-allocated. Element preparation reuses cache buffers. - -# Arguments -- `cache`: Pre-allocated node-based COO cache -- `assembler`: Node-based COO assembler -- `kernel`: Continuum kernel -- `mesh`: Finite element mesh -""" -function assemble!( - cache::NodeBasedCOOCache, - assembler::NodeBasedCOOAssembler, - kernel::ContinuumKernel, - mesh::AbstractMesh -) - # Reset cache - reset!(cache) - - nnodes = length(mesh.nodes) - ndofs_per_node = dofs_per_node(kernel) - - # NODAL LOOP: One iteration per node (GPU: one thread per node!) - for node_i in 1:nnodes - # Get all elements touching this node - touching_elements = cache.node_to_elements.node_to_elements[node_i] - - # Loop over touching elements - for elem_info in touching_elements - element_id = elem_info.element_id - local_i = elem_info.local_node_idx # Position of node_i in element - - # Prepare element geometry ONCE (reuses cache.element_cache) - prepared = prepare_element!(cache.element_cache, kernel, element_id, mesh) - - # Get element connectivity - conn = mesh.connectivity[element_id] - N = length(conn) # Nodes per element - - # Compute blocks for all nodes j in this element - for local_j in 1:N - global_j = conn[local_j] - - # Compute single 3×3 block K[i,j] - # This is THE KEY OPERATION: block-based integration - K_ij = compute_block!( - prepared, - kernel.material, - local_i, - local_j - ) - - # Scatter 3×3 block to triplets (adds 9 entries) - scatter_block_to_triplets!( - cache, - K_ij, - node_i, - global_j, - ndofs_per_node - ) - end - end - end - - return nothing -end - -""" - scatter_block_to_triplets!( - cache::NodeBasedCOOCache, - K_block::Tensor{2,3}, - node_i::Int, - node_j::Int, - ndofs_per_node::Int - ) - -Scatter single 3×3 block to COO triplets **in-place**. - -Maps block[α,β] → triplet at DOF indices: -- Row: 3*(node_i-1) + α -- Col: 3*(node_j-1) + β -- Val: K_block[α,β] - -# Arguments -- `cache`: Node-based COO cache -- `K_block`: 3×3 stiffness block (Tensor{2,3}) -- `node_i`: Global row node index -- `node_j`: Global column node index -- `ndofs_per_node`: DOFs per node (typically 3) - -# Zero-Allocation - -Writes to pre-allocated triplet arrays, updates counter. -""" -function scatter_block_to_triplets!( - cache::NodeBasedCOOCache, - K_block::Tensor{2,3,Float64}, - node_i::Int, - node_j::Int, - ndofs_per_node::Int -) - counter = cache.counter[] - - # Check capacity - new_triplets = ndofs_per_node * ndofs_per_node # 3×3 = 9 - if counter + new_triplets > cache.capacity - error("Node-based COO cache overflow: need $(counter + new_triplets) triplets, " * - "capacity is $(cache.capacity). Increase cache size.") - end - - # DOF offsets for nodes i and j - row_offset = ndofs_per_node * (node_i - 1) - col_offset = ndofs_per_node * (node_j - 1) - - # Scatter 3×3 block to triplets - for β in 1:ndofs_per_node # Column (node j DOF) - j_global = col_offset + β - for α in 1:ndofs_per_node # Row (node i DOF) - i_global = row_offset + α - counter += 1 - cache.I[counter] = i_global - cache.J[counter] = j_global - cache.V[counter] = K_block[α, β] - end - end - - cache.counter[] = counter - return nothing -end - -""" - extract_system(cache::NodeBasedCOOCache) -> (K, f) - -Build sparse matrix from triplets and return system. - -Calls `sparse(I, J, V)` to build CSC matrix. Duplicates are summed automatically. - -# Arguments -- `cache`: Assembled node-based COO cache - -# Returns -- `K::SparseMatrixCSC`: Global stiffness matrix -- `f::Vector`: Global force vector - -# Allocation - -Allocates sparse matrix structure (CSC format). This is the only allocation -outside cache construction. -""" -function extract_system(cache::NodeBasedCOOCache) - ntriplets = cache.counter[] - - # Build sparse matrix (duplicates are summed automatically) - I_used = @view cache.I[1:ntriplets] - J_used = @view cache.J[1:ntriplets] - V_used = @view cache.V[1:ntriplets] - - K = sparse(I_used, J_used, V_used, cache.ndofs, cache.ndofs) - - return K, cache.f -end - -# ============================================================================ -# HELPER FUNCTIONS -# ============================================================================ - -""" - create_cache( - assembler::NodeBasedCOOAssembler, - mesh::AbstractMesh, - kernel::ContinuumKernel - ) -> NodeBasedCOOCache - -Create pre-allocated cache for node-based COO assembly. - -Convenience function that wraps `NodeBasedCOOCache(mesh, kernel)`. - -# Example - -```julia -assembler = NodeBasedCOOAssembler() -cache = create_cache(assembler, mesh, kernel) -assemble!(cache, assembler, kernel, mesh) -K, f = extract_system(cache) -``` -""" -function create_cache( - assembler::NodeBasedCOOAssembler, - mesh::AbstractMesh, - kernel::ContinuumKernel -) - return NodeBasedCOOCache(mesh, kernel) -end - -""" - dofs_per_node(kernel::ContinuumKernel) -> Int - -Return DOFs per node for continuum kernel (always 3 for displacement). - -Dispatches on kernel field dimension. -""" -function dofs_per_node(kernel::ContinuumKernel{Theory,Mat}) where {Theory,Mat} - field = kernel.field - return field.dim # Displacement{3} → 3 -end - -# ============================================================================ -# PERFORMANCE NOTES -# ============================================================================ - -#= -# CPU Performance Comparison (Estimated) - -**Element-Based Assembly:** -- Elements: 1000 Tet4 -- Operations: 1000 elements × 4×4 blocks × 3×3 entries = 48,000 block computations -- Time: ~5ms (baseline) - -**Node-Based Assembly:** -- Nodes: 500 nodes -- Operations: 500 nodes × 8 elements/node × 4 blocks/element = 16,000 block computations -- But: 3× more kernel calls due to overlaps -- Time: ~7-10ms (1.5-2× slower single-threaded) - -**Why slower on CPU?** -- Each block computed once in element assembly -- Each block computed 2× on average in nodal assembly (shared between 2 elements) -- More function call overhead - -**Why faster on GPU?** -- Element assembly: Sequential (can't parallelize over elements efficiently) -- Nodal assembly: Massive parallelism (one thread per node) -- GPU speedup: ~10-50× depending on problem size - -**Multi-threaded CPU (Threads.@threads):** -- Can parallelize outer node loop -- Expected speedup: 1.5-2× over element-based -- No race conditions (each node writes different triplets) - -# Memory Comparison - -**Element-Based:** -- Triplet storage: ~50 KB per 1000 elements -- Element cache: ~2 KB per thread - -**Node-Based:** -- Triplet storage: Same (~50 KB) -- Element cache: ~2 KB per thread -- Inverse connectivity: ~10-20 KB (one-time) - -→ Nearly identical memory usage! - -# When to Use Node-Based Assembly - -**Use when:** -- ✅ GPU acceleration needed -- ✅ Contact mechanics (naturally nodal) -- ✅ Matrix-free methods (K*v without forming K) -- ✅ Adaptive refinement (local node operations) -- ✅ Multi-threading on CPU - -**Don't use when:** -- ❌ Single-threaded CPU only -- ❌ Simple problems (< 1000 nodes) -- ❌ Prototyping/debugging (element-based is clearer) -=# diff --git a/src/assembly/assembly.jl b/src/assembly/assembly.jl deleted file mode 100644 index bf169fa..0000000 --- a/src/assembly/assembly.jl +++ /dev/null @@ -1,201 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE - -function isapprox(a1::Assembly, a2::Assembly) - T = isapprox(a1.K, a2.K) - T &= isapprox(a1.C1, a2.C1) - T &= isapprox(a1.C2, a2.C2) - T &= isapprox(a1.D, a2.D) - T &= isapprox(a1.f, a2.f) - T &= isapprox(a1.g, a2.g) - return T -end - -function assemble_prehook!(::Problem, ::T) where T<:Number end - -function assemble_posthook!(::Problem, ::T) where T<:Number end - -""" - assemble_elements!(problem, assembly, elements, time) - -Assemble elements for problem. - -This should be overridden with own `assemble_elements!`-implementation. -""" -function assemble_elements!(problem::Problem, assembly::Assembly, - elements::Vector{T}, time) where T<:AbstractElement{E} where E - elements2 = convert(Vector{Element}, elements) - assemble!(assembly, problem, elements2, time) -end - -function assemble!(problem::Problem, time) - - assemble_prehook!(problem, time) - elements = get_elements(problem) - assembly = get_assembly(problem) - - if !isempty(assembly) - @warn("Problem assembly is not empty before assembling. This is probably " * - "causing unexpected results. To remove old assembly, use " * - "`empty!(problem.assembly)`", typeof(problem), problem.name) - assemble_posthook!(problem, time) - return nothing - end - - if isempty(elements) - @warn("There is no elements defined in problem. Before assembling a " * - "problem, elements must be added using " * - "`add_elements!(problem, elements)`.", typeof(problem), problem.name) - assemble_posthook!(problem, time) - return nothing - end - - first_element = first(elements) - unknown_field_name = get_unknown_field_name(problem) - if !haskey(first_element, unknown_field_name) - #= - warn("Assembling elements for problem $(problem.name): seems that ", - "problem is uninitialized. To initialize problem, use ", - "`initialize!(problem, time)`.") - info("Initializing problem $(problem.name) at time $time automatically.") - =# - initialize!(problem, time) - end - - for (element_type, elements) in group_by_element_type(elements) - assemble_elements!(problem, assembly, elements, time) - end - assemble_posthook!(problem, time) - return nothing -end - -function assemble!(problem::Problem) - @warn("assemble!(problem) will be deprecated. Use assemble!(problem, time)") - assemble!(problem, 0.0) -end - -function assemble_mass_matrix!(problem::Problem, time::Float64) - if !isempty(problem.assembly.M) - @info("Mass matrix for is already assembled, not assembling.", - problem.name) - return nothing - end - elements = get_elements(problem) - for (element_type, elements) in group_by_element_type(get_elements(problem)) - assemble_mass_matrix!(problem::Problem, elements, time) - end - return nothing -end - -function assemble_mass_matrix!(problem::Problem, elements::Vector{E}, time) where E<:AbstractElement{M_,B} where {M_,B} - nnodes = length(first(elements)) - dim = get_unknown_field_dimension(problem) - M = zeros(nnodes, nnodes) - N = zeros(1, nnodes) - NtN = zeros(nnodes, nnodes) - ldofs = zeros(Int, nnodes) - for element in elements - fill!(M, 0.0) - for ip in get_integration_points(element, 2) - detJ = element(ip, time, Val{:detJ}) - rho = element("density", ip, time) - w = ip.weight * rho * detJ - eval_basis!(B, N, ip) - N = element(ip, time) - mul!(NtN, transpose(N), N) - rmul!(NtN, w) - for i = 1:nnodes^2 - M[i] += NtN[i] - end - end - for (i, j) in enumerate(get_connectivity(element)) - @inbounds ldofs[i] = (j - 1) * dim - end - for i = 1:dim - add!(problem.assembly.M, ldofs .+ i, ldofs .+ i, M) - end - end - return -end - -# TODO (Phase 1B): Re-enable after resolving Tet10 topology vs Tet10Basis name conflict -# This specialized method assumes Tet10 <: AbstractBasis, but Tet10 is now a topology type. -# Need to refactor to use Tet10Basis or accept AbstractElement{M, T, Tet10Basis} -#= -""" - assemble_mass_matrix!(problem, elements::Vector{Element{Tet10}}, time) - -Assemble Tet10 mass matrices using special method. If Tet10 has constant metric -if can be integrated analytically to gain performance. -""" -function assemble_mass_matrix!(problem::Problem, elements::Vector{E}, time) where E<:AbstractElement{M_, Tet10} where M_ - nnodes = length(Tet10) - dim = get_unknown_field_dimension(problem) - M = zeros(nnodes, nnodes) - N = zeros(1, nnodes) - NtN = zeros(nnodes, nnodes) - ldofs = zeros(Int, nnodes) - - M_CM = 1.0/2520.0 * [ - 6 1 1 1 -4 -6 -4 -4 -6 -6 - 1 6 1 1 -4 -4 -6 -6 -4 -6 - 1 1 6 1 -6 -4 -4 -6 -6 -4 - 1 1 1 6 -6 -6 -6 -4 -4 -4 - -4 -4 -6 -6 32 16 16 16 16 8 - -6 -4 -4 -6 16 32 16 8 16 16 - -4 -6 -4 -6 16 16 32 16 8 16 - -4 -6 -6 -4 16 8 16 32 16 16 - -6 -4 -6 -4 16 16 8 16 32 16 - -6 -6 -4 -4 8 16 16 16 16 32] - - function is_CM(::AbstractElement{M, Tet10}, X; rtol=1.0e-6) where M - isapprox(X[5], 1/2*(X[1]+X[2]); rtol=rtol) || return false - isapprox(X[6], 1/2*(X[2]+X[3]); rtol=rtol) || return false - isapprox(X[7], 1/2*(X[3]+X[1]); rtol=rtol) || return false - isapprox(X[8], 1/2*(X[1]+X[4]); rtol=rtol) || return false - isapprox(X[9], 1/2*(X[2]+X[4]); rtol=rtol) || return false - isapprox(X[10], 1/2*(X[3]+X[4]); rtol=rtol) || return false - return true - end - - - n_CM = 0 - for element in elements - for (i, j) in enumerate(get_connectivity(element)) - @inbounds ldofs[i] = (j-1)*dim - end - - X = element("geometry", time) - rho = element("density", time) - if is_CM(element, X) && length(rho) == 1 - ip = (1.0/3.0, 1.0/3.0, 1.0/3.0) - detJ = element(ip, time, Val{:detJ}) - rho = element("density", ip, time) - CM_s = detJ*rho - n_CM += 1 - for i=1:dim - add!(problem.assembly.M, ldofs .+ i, ldofs .+ i, CM_s * M_CM) - end - else - fill!(M, 0.0) - for ip in get_integration_points(element, 2) - detJ = element(ip, time, Val{:detJ}) - rho = element("density", ip, time) - w = ip.weight*rho*detJ - eval_basis!(Tet10, N, ip) - N = element(ip, time) - mul!(NtN, transpose(N), N) - rmul!(NtN, w) - for i=1:nnodes^2 - M[i] += NtN[i] - end - end - for i=1:dim - add!(problem.assembly.M, ldofs .+ i, ldofs .+ i, M) - end - end - end - @info("$n_CM of $(length(elements)) was constant metric.") - return -end -=# diff --git a/src/assembly/element_structures.jl b/src/assembly/element_structures.jl deleted file mode 100644 index da43c25..0000000 --- a/src/assembly/element_structures.jl +++ /dev/null @@ -1,341 +0,0 @@ -# Traditional Element Assembly -# -# This module provides the standard element-by-element assembly approach -# for comparison with nodal assembly. Builds global tangent stiffness matrix -# and residual force vector using sparse matrix formats. - -using Tensors -using SparseArrays -using LinearAlgebra - -""" - ElementAssemblyData{T} - -Storage for element assembly using traditional (element-by-element) approach. - -# Fields -- `K_global::SparseMatrixCSC{T}`: Global tangent stiffness matrix -- `r_global::Vector{T}`: Global residual force vector (r = f_int - f_ext) -- `f_int_global::Vector{T}`: Global internal force vector -- `f_ext_global::Vector{T}`: Global external force vector -- `ndof::Int`: Total number of degrees of freedom - -# Notes -- Assembly uses COO (coordinate) format, then converts to CSC -- Multiple elements can write to same global DOF (summed automatically) -""" -mutable struct ElementAssemblyData{T} - K_global::SparseMatrixCSC{T,Int} - r_global::Vector{T} - f_int_global::Vector{T} - f_ext_global::Vector{T} - ndof::Int -end - -""" - ElementAssemblyData(ndof::Int, ::Type{T}=Float64) - -Allocate storage for traditional element assembly. - -# Arguments -- `ndof`: Total degrees of freedom (nnodes × 3 for 3D) -- `T`: Floating point type (default Float64) - -# Example -```julia -nnodes = 100 -assembly = ElementAssemblyData(3 * nnodes, Float64) -``` -""" -function ElementAssemblyData(ndof::Int, ::Type{T}=Float64) where T - # Pre-allocate empty sparse matrix (will fill during assembly) - K_global = spzeros(T, ndof, ndof) - r_global = zeros(T, ndof) - f_int_global = zeros(T, ndof) - f_ext_global = zeros(T, ndof) - - return ElementAssemblyData{T}(K_global, r_global, f_int_global, f_ext_global, ndof) -end - -""" - reset!(assembly::ElementAssemblyData) - -Reset assembly data to zero (for incremental/iterative solvers). -""" -function reset!(assembly::ElementAssemblyData{T}) where T - assembly.K_global = spzeros(T, assembly.ndof, assembly.ndof) - fill!(assembly.r_global, 0.0) - fill!(assembly.f_int_global, 0.0) - fill!(assembly.f_ext_global, 0.0) -end - -""" - ElementContribution{T} - -Local element contribution before scattering to global. - -# Fields -- `element_id::Int`: Element ID -- `gdofs::Vector{Int}`: Global DOF indices (e.g., [1,2,3,4,5,6,...] for nodes) -- `K_local::Matrix{T}`: Local stiffness matrix (ndofs_local × ndofs_local) -- `f_int_local::Vector{T}`: Local internal force vector -- `f_ext_local::Vector{T}`: Local external force vector - -# Notes -- For Tet4: ndofs_local = 12 (4 nodes × 3 DOF) -- For Tet10: ndofs_local = 30 (10 nodes × 3 DOF) -""" -struct ElementContribution{T} - element_id::Int - gdofs::Vector{Int} - K_local::Matrix{T} - f_int_local::Vector{T} - f_ext_local::Vector{T} -end - -""" - ElementContribution(element_id::Int, gdofs::Vector{Int}, ::Type{T}=Float64) - -Allocate storage for element contribution. - -# Arguments -- `element_id`: Element ID -- `gdofs`: Global DOF indices -- `T`: Floating point type - -# Example -```julia -# Tet4 element connecting nodes [5, 7, 12, 15] -gdofs = [13,14,15, 19,20,21, 34,35,36, 43,44,45] # 3 DOF per node -contrib = ElementContribution(1, gdofs, Float64) -``` -""" -function ElementContribution(element_id::Int, gdofs::Vector{Int}, ::Type{T}=Float64) where T - ndofs = length(gdofs) - K_local = zeros(T, ndofs, ndofs) - f_int_local = zeros(T, ndofs) - f_ext_local = zeros(T, ndofs) - - return ElementContribution{T}(element_id, gdofs, K_local, f_int_local, f_ext_local) -end - -""" - scatter_to_global!(assembly::ElementAssemblyData, contrib::ElementContribution) - -Scatter element contribution to global matrices/vectors (traditional assembly). - -This is the key operation in element assembly: add local element quantities -to global system. Uses COO format (accumulates into lists). - -# Arguments -- `assembly`: Global assembly data -- `contrib`: Element contribution - -# Notes -- Multiple elements can contribute to same global DOF (summed) -- For GPU: Would require atomic operations (slow!) -- For CPU: Direct scatter-add works fine -""" -function scatter_to_global!(assembly::ElementAssemblyData{T}, - contrib::ElementContribution{T}) where T - # Scatter forces (simple vector addition) - for (local_i, global_i) in enumerate(contrib.gdofs) - assembly.f_int_global[global_i] += contrib.f_int_local[local_i] - assembly.f_ext_global[global_i] += contrib.f_ext_local[local_i] - end - - # Scatter stiffness (matrix addition) - # Build list of (I, J, V) triplets for sparse matrix - I_rows = Int[] - J_cols = Int[] - values = T[] - - ndofs_local = length(contrib.gdofs) - for i in 1:ndofs_local, j in 1:ndofs_local - if abs(contrib.K_local[i, j]) > 1e-14 # Skip near-zeros - push!(I_rows, contrib.gdofs[i]) - push!(J_cols, contrib.gdofs[j]) - push!(values, contrib.K_local[i, j]) - end - end - - # Add to existing sparse matrix - K_elem = sparse(I_rows, J_cols, values, assembly.ndof, assembly.ndof) - assembly.K_global += K_elem -end - -""" - compute_residual!(assembly::ElementAssemblyData) - -Compute residual force vector: r = f_int - f_ext - -Should be called after all elements have been assembled. -""" -function compute_residual!(assembly::ElementAssemblyData{T}) where T - assembly.r_global .= assembly.f_int_global .- assembly.f_ext_global -end - -""" - assemble_elements!(assembly::ElementAssemblyData, - contributions::Vector{ElementContribution}) - -Assemble all element contributions to global system. - -# Arguments -- `assembly`: Global assembly data (modified in-place) -- `contributions`: Vector of element contributions - -# Example -```julia -assembly = ElementAssemblyData(ndof) -contributions = compute_all_element_contributions(elements, u, time) -assemble_elements!(assembly, contributions) -compute_residual!(assembly) - -# Now solve: K_global * Δu = -r_global -``` -""" -function assemble_elements!(assembly::ElementAssemblyData{T}, - contributions::Vector{ElementContribution{T}}) where T - reset!(assembly) - - # Loop over elements and scatter (element assembly) - for contrib in contributions - scatter_to_global!(assembly, contrib) - end - - # Compute residual - compute_residual!(assembly) -end - -""" - apply_dirichlet_bc!(assembly::ElementAssemblyData, - fixed_dofs::Vector{Int}, - prescribed_values::Vector{T}=zeros(length(fixed_dofs))) - -Apply Dirichlet (essential) boundary conditions by penalty method. - -# Arguments -- `assembly`: Global assembly data (modified in-place) -- `fixed_dofs`: DOF indices to fix -- `prescribed_values`: Prescribed displacement values (default: zeros) - -# Method -Uses penalty method: adds large stiffness to diagonal and corresponding RHS. - -For DOF i with prescribed value u_prescribed: -- K[i,i] += penalty (e.g., 1e10 * max_K) -- r[i] = penalty * (u_current - u_prescribed) - -# Example -```julia -# Fix nodes 1 and 2 in all directions (zero displacement) -fixed_dofs = [1,2,3, 4,5,6] # Nodes 1,2 × 3 DOF -apply_dirichlet_bc!(assembly, fixed_dofs) -``` -""" -function apply_dirichlet_bc!(assembly::ElementAssemblyData{T}, - fixed_dofs::Vector{Int}, - prescribed_values::Vector{T}=zeros(T, length(fixed_dofs))) where T - # Penalty parameter (large relative to stiffness) - max_K = maximum(abs, assembly.K_global) - penalty = 1e10 * max_K - - for (idx, dof) in enumerate(fixed_dofs) - # Add penalty stiffness to diagonal - assembly.K_global[dof, dof] += penalty - - # Modify residual (assuming current displacement is zero for now) - # In full Newton: r[i] += penalty * (u_current[i] - u_prescribed[i]) - assembly.r_global[dof] = penalty * prescribed_values[idx] - end -end - -""" - get_dof_indices(connectivity::NTuple{N,Int}, dim::Int=3) -> Vector{Int} - -Get global DOF indices for an element given node connectivity. - -# Arguments -- `connectivity`: Element node IDs (e.g., (5, 7, 12, 15) for Tet4) -- `dim`: Dimension (3 for 3D elasticity) - -# Returns -- `gdofs::Vector{Int}`: Global DOF indices - -# Example -```julia -# Element with nodes [5, 7, 12, 15] -gdofs = get_dof_indices((5, 7, 12, 15), 3) -# Returns: [13,14,15, 19,20,21, 34,35,36, 43,44,45] -``` -""" -function get_dof_indices(connectivity::NTuple{N,Int}, dim::Int=3) where N - nnodes = length(connectivity) - gdofs = zeros(Int, dim * nnodes) - - for (local_i, global_node) in enumerate(connectivity) - for d in 1:dim - gdofs[dim*(local_i-1)+d] = dim * (global_node - 1) + d - end - end - - return gdofs -end - -""" - matrix_vector_product(assembly::ElementAssemblyData, v::Vector{T}) -> Vector{T} - -Compute matrix-vector product: w = K * v using assembled sparse matrix. - -# Arguments -- `assembly`: Assembly data (contains K_global) -- `v`: Input vector (ndof) - -# Returns -- `w`: Output vector w = K * v - -# Example -```julia -# GMRES matrix-free operator -function matvec(v) - return matrix_vector_product(assembly, v) -end -Δu = gmres(matvec, -r, tol=1e-6) -``` -""" -function matrix_vector_product(assembly::ElementAssemblyData{T}, v::Vector{T}) where T - return assembly.K_global * v -end - -""" - print_assembly_stats(assembly::ElementAssemblyData) - -Print statistics about assembled system (for debugging). -""" -function print_assembly_stats(assembly::ElementAssemblyData) - nnz_K = nnz(assembly.K_global) - ndof = assembly.ndof - fill_ratio = nnz_K / (ndof * ndof) - - println("="^60) - println("Traditional Element Assembly Statistics") - println("="^60) - println(" Total DOF: ", ndof) - println(" K matrix size: ", size(assembly.K_global)) - println(" K non-zeros: ", nnz_K) - println(" K fill ratio: ", round(fill_ratio, sigdigits=2)) - println(" K memory (MB): ", round(nnz_K * 16 / 1024^2, digits=2)) - println(" ||f_int||: ", round(norm(assembly.f_int_global), sigdigits=2)) - println(" ||f_ext||: ", round(norm(assembly.f_ext_global), sigdigits=2)) - println(" ||residual||: ", round(norm(assembly.r_global), sigdigits=2)) - println(" K symmetric: ", issymmetric(assembly.K_global)) - println("="^60) -end - -# Export main types and functions -export ElementAssemblyData, ElementContribution -export reset!, scatter_to_global!, compute_residual! -export assemble_elements!, apply_dirichlet_bc! -export get_dof_indices, matrix_vector_product -export print_assembly_stats diff --git a/src/assembly/framework.jl b/src/assembly/framework.jl deleted file mode 100644 index bf169fa..0000000 --- a/src/assembly/framework.jl +++ /dev/null @@ -1,201 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE - -function isapprox(a1::Assembly, a2::Assembly) - T = isapprox(a1.K, a2.K) - T &= isapprox(a1.C1, a2.C1) - T &= isapprox(a1.C2, a2.C2) - T &= isapprox(a1.D, a2.D) - T &= isapprox(a1.f, a2.f) - T &= isapprox(a1.g, a2.g) - return T -end - -function assemble_prehook!(::Problem, ::T) where T<:Number end - -function assemble_posthook!(::Problem, ::T) where T<:Number end - -""" - assemble_elements!(problem, assembly, elements, time) - -Assemble elements for problem. - -This should be overridden with own `assemble_elements!`-implementation. -""" -function assemble_elements!(problem::Problem, assembly::Assembly, - elements::Vector{T}, time) where T<:AbstractElement{E} where E - elements2 = convert(Vector{Element}, elements) - assemble!(assembly, problem, elements2, time) -end - -function assemble!(problem::Problem, time) - - assemble_prehook!(problem, time) - elements = get_elements(problem) - assembly = get_assembly(problem) - - if !isempty(assembly) - @warn("Problem assembly is not empty before assembling. This is probably " * - "causing unexpected results. To remove old assembly, use " * - "`empty!(problem.assembly)`", typeof(problem), problem.name) - assemble_posthook!(problem, time) - return nothing - end - - if isempty(elements) - @warn("There is no elements defined in problem. Before assembling a " * - "problem, elements must be added using " * - "`add_elements!(problem, elements)`.", typeof(problem), problem.name) - assemble_posthook!(problem, time) - return nothing - end - - first_element = first(elements) - unknown_field_name = get_unknown_field_name(problem) - if !haskey(first_element, unknown_field_name) - #= - warn("Assembling elements for problem $(problem.name): seems that ", - "problem is uninitialized. To initialize problem, use ", - "`initialize!(problem, time)`.") - info("Initializing problem $(problem.name) at time $time automatically.") - =# - initialize!(problem, time) - end - - for (element_type, elements) in group_by_element_type(elements) - assemble_elements!(problem, assembly, elements, time) - end - assemble_posthook!(problem, time) - return nothing -end - -function assemble!(problem::Problem) - @warn("assemble!(problem) will be deprecated. Use assemble!(problem, time)") - assemble!(problem, 0.0) -end - -function assemble_mass_matrix!(problem::Problem, time::Float64) - if !isempty(problem.assembly.M) - @info("Mass matrix for is already assembled, not assembling.", - problem.name) - return nothing - end - elements = get_elements(problem) - for (element_type, elements) in group_by_element_type(get_elements(problem)) - assemble_mass_matrix!(problem::Problem, elements, time) - end - return nothing -end - -function assemble_mass_matrix!(problem::Problem, elements::Vector{E}, time) where E<:AbstractElement{M_,B} where {M_,B} - nnodes = length(first(elements)) - dim = get_unknown_field_dimension(problem) - M = zeros(nnodes, nnodes) - N = zeros(1, nnodes) - NtN = zeros(nnodes, nnodes) - ldofs = zeros(Int, nnodes) - for element in elements - fill!(M, 0.0) - for ip in get_integration_points(element, 2) - detJ = element(ip, time, Val{:detJ}) - rho = element("density", ip, time) - w = ip.weight * rho * detJ - eval_basis!(B, N, ip) - N = element(ip, time) - mul!(NtN, transpose(N), N) - rmul!(NtN, w) - for i = 1:nnodes^2 - M[i] += NtN[i] - end - end - for (i, j) in enumerate(get_connectivity(element)) - @inbounds ldofs[i] = (j - 1) * dim - end - for i = 1:dim - add!(problem.assembly.M, ldofs .+ i, ldofs .+ i, M) - end - end - return -end - -# TODO (Phase 1B): Re-enable after resolving Tet10 topology vs Tet10Basis name conflict -# This specialized method assumes Tet10 <: AbstractBasis, but Tet10 is now a topology type. -# Need to refactor to use Tet10Basis or accept AbstractElement{M, T, Tet10Basis} -#= -""" - assemble_mass_matrix!(problem, elements::Vector{Element{Tet10}}, time) - -Assemble Tet10 mass matrices using special method. If Tet10 has constant metric -if can be integrated analytically to gain performance. -""" -function assemble_mass_matrix!(problem::Problem, elements::Vector{E}, time) where E<:AbstractElement{M_, Tet10} where M_ - nnodes = length(Tet10) - dim = get_unknown_field_dimension(problem) - M = zeros(nnodes, nnodes) - N = zeros(1, nnodes) - NtN = zeros(nnodes, nnodes) - ldofs = zeros(Int, nnodes) - - M_CM = 1.0/2520.0 * [ - 6 1 1 1 -4 -6 -4 -4 -6 -6 - 1 6 1 1 -4 -4 -6 -6 -4 -6 - 1 1 6 1 -6 -4 -4 -6 -6 -4 - 1 1 1 6 -6 -6 -6 -4 -4 -4 - -4 -4 -6 -6 32 16 16 16 16 8 - -6 -4 -4 -6 16 32 16 8 16 16 - -4 -6 -4 -6 16 16 32 16 8 16 - -4 -6 -6 -4 16 8 16 32 16 16 - -6 -4 -6 -4 16 16 8 16 32 16 - -6 -6 -4 -4 8 16 16 16 16 32] - - function is_CM(::AbstractElement{M, Tet10}, X; rtol=1.0e-6) where M - isapprox(X[5], 1/2*(X[1]+X[2]); rtol=rtol) || return false - isapprox(X[6], 1/2*(X[2]+X[3]); rtol=rtol) || return false - isapprox(X[7], 1/2*(X[3]+X[1]); rtol=rtol) || return false - isapprox(X[8], 1/2*(X[1]+X[4]); rtol=rtol) || return false - isapprox(X[9], 1/2*(X[2]+X[4]); rtol=rtol) || return false - isapprox(X[10], 1/2*(X[3]+X[4]); rtol=rtol) || return false - return true - end - - - n_CM = 0 - for element in elements - for (i, j) in enumerate(get_connectivity(element)) - @inbounds ldofs[i] = (j-1)*dim - end - - X = element("geometry", time) - rho = element("density", time) - if is_CM(element, X) && length(rho) == 1 - ip = (1.0/3.0, 1.0/3.0, 1.0/3.0) - detJ = element(ip, time, Val{:detJ}) - rho = element("density", ip, time) - CM_s = detJ*rho - n_CM += 1 - for i=1:dim - add!(problem.assembly.M, ldofs .+ i, ldofs .+ i, CM_s * M_CM) - end - else - fill!(M, 0.0) - for ip in get_integration_points(element, 2) - detJ = element(ip, time, Val{:detJ}) - rho = element("density", ip, time) - w = ip.weight*rho*detJ - eval_basis!(Tet10, N, ip) - N = element(ip, time) - mul!(NtN, transpose(N), N) - rmul!(NtN, w) - for i=1:nnodes^2 - M[i] += NtN[i] - end - end - for i=1:dim - add!(problem.assembly.M, ldofs .+ i, ldofs .+ i, M) - end - end - end - @info("$n_CM of $(length(elements)) was constant metric.") - return -end -=# diff --git a/src/assembly/nodal_structures.jl b/src/assembly/nodal_structures.jl deleted file mode 100644 index e9a2e8f..0000000 --- a/src/assembly/nodal_structures.jl +++ /dev/null @@ -1,234 +0,0 @@ -# Nodal Assembly Data Structures -# -# This module provides the inverse mapping needed for efficient nodal assembly: -# Given a node, find all elements touching it and the local node index within each element. - -using Tensors - -""" - ElementNodeInfo - -Information about how a node appears in an element. - -# Fields -- `element_id::Int`: Global element ID -- `local_node_idx::Int`: Local node index within the element (1-based) -""" -struct ElementNodeInfo - element_id::Int - local_node_idx::Int -end - -""" - NodeToElementsMap - -Inverse connectivity mapping: for each node, lists all elements touching it. - -# Fields -- `node_to_elements::Vector{Vector{ElementNodeInfo}}`: For node j, gives all elements touching it -- `nnodes::Int`: Total number of nodes in mesh -- `nelements::Int`: Total number of elements in mesh - -# Example -```julia -map = NodeToElementsMap(connectivity) -# Get all elements touching node 5 -elements_touching_5 = map.node_to_elements[5] -for info in elements_touching_5 - println("Node 5 is local node ", info.local_node_idx, " in element ", info.element_id) -end -``` -""" -struct NodeToElementsMap - node_to_elements::Vector{Vector{ElementNodeInfo}} - nnodes::Int - nelements::Int -end - -""" - NodeToElementsMap(connectivity::Vector{NTuple{N,Int}}) where N - -Build inverse mapping from element connectivity. - -# Arguments -- `connectivity`: Vector of element connectivity tuples, e.g., [(1,2,3,4), (2,3,5,6), ...] - -# Returns -- `NodeToElementsMap`: Inverse mapping structure - -# Example -```julia -# Tet4 mesh with 2 elements -connectivity = [(1,2,3,4), (2,3,4,5)] -map = NodeToElementsMap(connectivity) - -# Node 2 appears in both elements -@assert length(map.node_to_elements[2]) == 2 -``` -""" -function NodeToElementsMap(connectivity::Vector{NTuple{N,Int}}) where N - nelements = length(connectivity) - - # Find maximum node ID to determine array size - nnodes = maximum(maximum(conn) for conn in connectivity) - - # Pre-allocate vectors for each node - node_to_elements = [Vector{ElementNodeInfo}() for _ in 1:nnodes] - - # Build inverse mapping - for (elem_id, conn) in enumerate(connectivity) - for (local_idx, global_node_id) in enumerate(conn) - push!(node_to_elements[global_node_id], - ElementNodeInfo(elem_id, local_idx)) - end - end - - return NodeToElementsMap(node_to_elements, nnodes, nelements) -end - -""" - get_node_spider(map::NodeToElementsMap, node_id::Int) -> Vector{Int} - -Get the "spider" of a node - all nodes that couple with it (including itself). - -This is the union of all nodes in elements touching `node_id`. These are exactly -the nodes for which we need to compute 3×3 stiffness blocks. - -# Arguments -- `map`: Node-to-elements mapping -- `node_id`: Node for which to find the spider - -# Returns -- `spider_nodes::Vector{Int}`: Sorted unique list of node IDs in the spider - -# Example -```julia -# For node j, find all nodes it couples with -spider = get_node_spider(map, j) -# Now compute K_blocks[k] for each k in spider -``` -""" -function get_node_spider(map::NodeToElementsMap, node_id::Int, - connectivity::Vector{NTuple{N,Int}}) where N - spider = Set{Int}() - - # For each element touching this node - for elem_info in map.node_to_elements[node_id] - # Add all nodes in that element - for node in connectivity[elem_info.element_id] - push!(spider, node) - end - end - - return sort(collect(spider)) -end - -""" - NodalStiffnessContribution{T} - -Storage for nodal assembly contribution at a single node. - -# Fields -- `node_id::Int`: Global node ID -- `spider_nodes::Vector{Int}`: Node IDs that couple with this node -- `K_blocks::Vector{Tensor{2,3,T}}`: 3×3 stiffness blocks for each spider node -- `f_int::Vec{3,T}`: Internal force at this node -- `f_ext::Vec{3,T}`: External force at this node - -# Notes -- `K_blocks[k]` corresponds to `spider_nodes[k]` -- Diagonal block (self-coupling) is included in spider -- All quantities use Tensors.jl types (zero-allocation) -""" -struct NodalStiffnessContribution{T} - node_id::Int - spider_nodes::Vector{Int} - K_blocks::Vector{Tensor{2,3,T,9}} - f_int::Vec{3,T} - f_ext::Vec{3,T} -end - -""" - NodalStiffnessContribution(node_id::Int, spider_nodes::Vector{Int}, ::Type{T}=Float64) - -Allocate storage for nodal assembly contribution. - -# Example -```julia -spider = get_node_spider(map, 5, connectivity) -contrib = NodalStiffnessContribution(5, spider, Float64) -# Now fill in K_blocks, f_int, f_ext during assembly -``` -""" -function NodalStiffnessContribution(node_id::Int, spider_nodes::Vector{Int}, - ::Type{T}=Float64) where T - nspider = length(spider_nodes) - K_blocks = [zero(Tensor{2,3,T}) for _ in 1:nspider] - f_int = zero(Vec{3,T}) - f_ext = zero(Vec{3,T}) - - return NodalStiffnessContribution{T}(node_id, spider_nodes, K_blocks, f_int, f_ext) -end - -""" - matrix_vector_product_nodal(contrib::NodalStiffnessContribution, - u::Vector{Vec{3,T}}) -> Vec{3,T} - -Compute the matrix-vector product for one node using nodal assembly. - -This computes: w_i = sum_j K_ij * u_j for node i - -# Arguments -- `contrib`: Nodal stiffness contribution (contains K_blocks for all j in spider) -- `u`: Displacement field at all nodes (Vec{3} per node) - -# Returns -- `w_i::Vec{3}`: Result of K_i * u at this node - -# Example -```julia -# Assemble contribution for node i -contrib = assemble_nodal_contribution(element_set, node_i, u, time) - -# Matrix-free matvec: w_i = K_i * u -w_i = matrix_vector_product_nodal(contrib, u) -``` -""" -function matrix_vector_product_nodal(contrib::NodalStiffnessContribution{T}, - u::Vector{Vec{3,T}}) where T - w = zero(Vec{3,T}) - - # 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] # 3×1 displacement - - # Block matrix-vector product: K_ij is Tensor{2,3}, u_j is Vec{3} - # Use regular matrix-vector multiplication (single contraction) - w += K_ij ⋅ u_j # Tensor{2,3} ⋅ Vec{3} → Vec{3} - end - - return w -end - -""" - print_spider_info(map::NodeToElementsMap, node_id::Int, - connectivity::Vector{NTuple{N,Int}}) where N - -Print diagnostic information about a node's spider for debugging. -""" -function print_spider_info(map::NodeToElementsMap, node_id::Int, - connectivity::Vector{NTuple{N,Int}}) where N - println("Node $node_id Spider Analysis:") - println(" Touches $(length(map.node_to_elements[node_id])) elements") - - for elem_info in map.node_to_elements[node_id] - println(" Element $(elem_info.element_id): local node $(elem_info.local_node_idx)") - println(" Connectivity: $(connectivity[elem_info.element_id])") - end - - spider = get_node_spider(map, node_id, connectivity) - println(" Spider has $(length(spider)) nodes: $spider") - println(" → Need to compute $(length(spider)) 3×3 blocks") - println(" → Diagonal block at node $node_id") -end diff --git a/src/assembly/problems.jl b/src/assembly/problems.jl deleted file mode 100644 index 7c8dd95..0000000 --- a/src/assembly/problems.jl +++ /dev/null @@ -1,478 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE - -abstract type AbstractProblem end -abstract type FieldProblem <: AbstractProblem end -abstract type BoundaryProblem <: AbstractProblem end -abstract type MixedProblem <: AbstractProblem end - -""" -General linearized problem to solve - (K₁+K₂)Δu + C1'*Δλ = f₁+f₂ - C2Δu + D*Δλ = g -""" -mutable struct Assembly - - M::SparseMatrixCOO # mass matrix - - # for field assembly - K::SparseMatrixCOO # stiffness matrix - Kg::SparseMatrixCOO # geometric stiffness matrix - f::SparseMatrixCOO # force vector - fg::SparseMatrixCOO # - - # for boundary assembly - C1::SparseMatrixCOO - C2::SparseMatrixCOO - D::SparseMatrixCOO - g::SparseMatrixCOO - c::SparseMatrixCOO - - u::Vector{Float64} # solution vector u - u_prev::Vector{Float64} # previous solution vector u - u_norm_change::Real # change of norm in u - - la::Vector{Float64} # solution vector la - la_prev::Vector{Float64} # previous solution vector u - la_norm_change::Real # change of norm in la - - removed_dofs::Vector{Int} # manually remove dofs from assembly -end - -function Assembly() - return Assembly( - SparseMatrixCOO(), - SparseMatrixCOO(), - SparseMatrixCOO(), - SparseMatrixCOO(), - SparseMatrixCOO(), - SparseMatrixCOO(), - SparseMatrixCOO(), - SparseMatrixCOO(), - SparseMatrixCOO(), - SparseMatrixCOO(), - [], [], Inf, - [], [], Inf, - []) -end - -function empty!(assembly::Assembly) - empty!(assembly.M) - empty!(assembly.K) - empty!(assembly.Kg) - empty!(assembly.f) - empty!(assembly.fg) - empty!(assembly.C1) - empty!(assembly.C2) - empty!(assembly.D) - empty!(assembly.g) - empty!(assembly.c) -end - -function Base.isempty(assembly::Assembly) - T = Base.isempty(assembly.M) - T &= Base.isempty(assembly.K) - T &= Base.isempty(assembly.Kg) - T &= Base.isempty(assembly.f) - T &= Base.isempty(assembly.fg) - T &= Base.isempty(assembly.C1) - T &= Base.isempty(assembly.C2) - T &= Base.isempty(assembly.D) - T &= isempty(assembly.g) - T &= isempty(assembly.c) - return T -end - -""" - Problem{P<:AbstractProblem} - -Defines a new problem of type `P`, where `P` characterizes the physics of the -problem. `P` can be for example `Elasticity`, if the physics of the system is -described by Cauchy's stress equation ∇⋅σ + b = ̈ρu, or `Heat`, if the physics -of the problem is described by heat equation -∇⋅(k∇u) = f. - -""" -mutable struct Problem{P<:AbstractProblem} - name::AbstractString # descriptive name for the problem - dimension::Int # degrees of freedom per node - parent_field_name::AbstractString # (optional) name of the parent field e.g. "displacement" - elements::Vector{Element} - dofmap::Dict{Element,Vector{Int}} # connects the element local dofs to the global dofs - assembly::Assembly - fields::Dict{String,AbstractField} - postprocess_fields::Vector{String} - properties::P -end - -""" - Problem(problem_type, problem_name, problem_dimension) - -Construct a new field problem. - -`problem_type` must be a subtype of `FieldProblem` (`Elasticity`, `Heat`, etc..). -`problem_dimensions` is the number of degrees of freedom each node is containing. - -# Examples - -To create vector-valued elasticity problem, having 3 dofs / node: -```julia -problem1 = Problem(Elasticity, "test problem", 3) -``` - -To create scalar-valued Poisson problem: -```julia -problem2 = Problem(Heat, "test problem 2", 1) -``` - -""" -function Problem(::Type{P}, name::AbstractString, dimension::Int) where P<:FieldProblem - parent_field_name = "none" - elements = [] - dofmap = Dict() - assembly = Assembly() - fields = Dict() - postprocess_fields = Vector() - properties = P() - problem = Problem{P}(name, dimension, parent_field_name, elements, dofmap, - assembly, fields, postprocess_fields, properties) - @info("Creating a new problem of type $P, having name `$name` and " * - "dimension $dimension dofs/node.") - return problem -end - -""" - Problem(problem_type, problem_name, problem_dimension, parent_field_name) - -Construct a new boundary problem. - -`problem_type` must be a subtype of `BoundaryProblem` (`Dirichlet`, `Contact`, -etc..). `problem_dimensions` is the number of degrees of freedom each node is -containing. `parent_field_name` is describing the field, where the boundary -problem is affecting. - -# Examples - -To create a Dirichlet boundary condition for a vector-valued elasticity problem, -having 3 dofs / node: -```julia -bc1 = Problem(Dirichlet, "fix displacement on support", 3, "displacement") -``` - -To create a Dirichlet boundary condition for scalar-valued Poisson problem: -```julia -bc2 = Problem(Dirichlet, "fix surface temperature", 1, "temperature") -``` -""" -function Problem(::Type{P}, name, dimension, parent_field_name) where P<:BoundaryProblem - elements = [] - dofmap = Dict() - assembly = Assembly() - fields = Dict() - postprocess_fields = Vector() - properties = P() - problem = Problem{P}(name, dimension, parent_field_name, elements, dofmap, - assembly, fields, postprocess_fields, properties) - @info("Creating a new boundary problem of type $P, having name `$name` and " * - "dimension $dimension dofs/node. This boundary problems fixes field " * - "`$parent_field_name`.") - return problem -end - -function get_formulation_type(::Problem) - return :incremental -end - -""" - get_unknown_field_dimension(problem) - -Return the dimension of the unknown field of this problem. -""" -function get_unknown_field_dimension(problem::Problem) - return problem.dimension -end - -""" - get_unknown_field_name(problem) - -Default function if unknown field name is not defined for some problem. -""" -function get_unknown_field_name(::P) where P<:AbstractProblem - @warn("The name of unknown field (e.g. displacement, temperature, ...) of the " * - "problem type must be given by defining a function " * - "`get_unknown_field_name(::$P)`") - return "N/A" -end - -""" Return the name of the unknown field of this problem. """ -function get_unknown_field_name(problem::Problem{P}) where P - return get_unknown_field_name(problem.properties) -end - -""" Return the name of the parent field of this (boundary) problem. """ -function get_parent_field_name(problem::Problem{P}) where P<:BoundaryProblem - return problem.parent_field_name -end - -function get_unknown_field_name(::P) where P<:BoundaryProblem - return "lambda" -end - -is_field_problem(::Problem) = false -is_field_problem(::Problem{P}) where {P<:FieldProblem} = true -is_boundary_problem(::Problem) = false -is_boundary_problem(::Problem{P}) where {P<:BoundaryProblem} = true - -function get_elements(problem::Problem) - return problem.elements -end - -function update!(problem::P, attr::Pair{String,String}...) where P<:AbstractProblem - for (name, value) in attr - setfield!(problem, Meta.parse(name), Meta.parse(value)) - end -end - -""" - function initialize!(problem_type, element_name, time) - -Initialize the element ready for calculation, where `problem_type` is the type -of the problem (Elasticity, Dirichlet, etc.), `element_name` is the name of a -constructed element (see Element(element_type, connectivity_vector)) and `time` -is the starting time of the initializing process. -""" -function initialize!(problem::Problem, element::AbstractElement, time::Float64) - field_name = get_unknown_field_name(problem) - field_dim = get_unknown_field_dimension(problem) - nnodes = length(element) - if field_dim == 1 # scalar field - empty_field = tuple(zeros(nnodes)...) - else # vector field - # FIXME: the most effective way to do - # ([0.0,0.0], [0.0,0.0], ..., [0.0,0.0]) ? - empty_field = tuple(map((x) -> zeros(field_dim) * x, 1:nnodes)...) - end - - # initialize primary field - if !haskey(element, field_name) - update!(element, field_name, time => empty_field) - end - - # if a boundary problem, initialize also a field for the main problem - is_boundary_problem(problem) || return - field_name = get_parent_field_name(problem) - if !haskey(element, field_name) - update!(element, field_name, time => empty_field) - end -end - -function initialize!(problem::Problem, time::Float64=0.0) - for element in get_elements(problem) - initialize!(problem, element, time) - end -end - -function update!(problem::Problem, assembly::Assembly, u::Vector, la::Vector) - - # resize & fill with zeros vectors if length mismatch with current solution - - if length(u) != length(assembly.u) - resize!(assembly.u, length(u)) - fill!(assembly.u, 0.0) - end - - if length(la) != length(assembly.la) - resize!(assembly.la, length(la)) - fill!(assembly.la, 0.0) - end - - # copy current solutions to previous ones and add/replace new solution - # TODO: here we have couple of options and they need to be clarified - # for total formulation we are solving total quantity Ku = f while in - # incremental formulation we solve KΔu = f and u = u + Δu - assembly.u_prev = copy(assembly.u) - assembly.la_prev = copy(assembly.la) - - if get_formulation_type(problem) == :total - assembly.u = u - assembly.la = la - elseif get_formulation_type(problem) == :incremental - assembly.u += u - assembly.la = la - elseif get_formulation_type(problem) == :forwarddiff - assembly.u += u - assembly.la += la - else - @info("$(problem.name): unknown formulation type, don't know what to do with results") - error("serious failure with problem formulation: $(get_formulation_type(problem))") - end - - # calculate change of norm - assembly.u_norm_change = norm(assembly.u - assembly.u_prev) - assembly.la_norm_change = norm(assembly.la - assembly.la_prev) - return assembly.u, assembly.la -end - -""" - get_global_solution(problem, assembly) - -Return a global solution (u, la) for a problem. - -Notes ------ -If the length of solution vector != number of nodes, i.e. the field dimension is -something else than 1, reshape vectors so that their length matches to the -number of nodes. This helps to get nodal results easily. -""" -function get_global_solution(problem::Problem, assembly::Assembly) - u = assembly.u - la = assembly.la - field_dim = get_unknown_field_dimension(problem) - if field_dim == 1 - return u, la - else - nnodes = round(Int, length(u) / field_dim) - u = reshape(u, field_dim, nnodes) - u = Vector{Float64}[u[:, i] for i in 1:nnodes] - la = reshape(la, field_dim, nnodes) - la = Vector{Float64}[la[:, i] for i in 1:nnodes] - return u, la - end -end - -function update!(problem::Problem{P}, assembly::Assembly, elements::Vector{Element}, time::Float64) where P<:FieldProblem - u, la = get_global_solution(problem, assembly) - field_name = get_unknown_field_name(problem) - # update solution u for elements - for element in elements - connectivity = get_connectivity(element) - update!(element, field_name, time => tuple(u[connectivity]...)) - end -end - -function update!(problem::Problem{P}, assembly::Assembly, elements::Vector{Element}, time::Float64) where P<:BoundaryProblem - u, la = get_global_solution(problem, assembly) - parent_field_name = get_parent_field_name(problem) # displacement - field_name = get_unknown_field_name(problem) # lambda - # update solution and lagrange multipliers for boundary elements - for element in elements - connectivity = get_connectivity(element) - update!(element, parent_field_name, time => tuple(u[connectivity]...)) - update!(element, field_name, time => tuple(la[connectivity]...)) - end -end - -""" - add_element!(problem, element1, element2, ...) - -Add element(s) to the problem. -""" -function add_element!(problem, elements...) - for element in elements - push!(problem.elements, element) - end - return nothing -end - -""" - add_elements!(problem, element_set_1, element_set_2, ...) - -Add vectors/tuples of element(s) to the problem. -""" -function add_elements!(problem, element_sets::Union{Vector,Tuple}...) - for elements in element_sets - nelements = length(elements) - @info("Adding $nelements elements to problem `$(problem.name)`") - add_element!(problem, elements...) - end - return nothing -end - -add_elements!(problem, elements::Element...) = add_element!(problem, elements...) - -function add_elements!(problem, elements_or_lists_of_elements...) - for item in elements_or_lists_of_elements - add_elements!(problem, item) - end -end - -get_assembly(problem::Problem) = problem.assembly -Base.length(problem::Problem) = length(problem.elements) - -function update!(problem::Problem, field_name::AbstractString, data) - #if haskey(problem.fields, field_name) - # update!(problem.fields[field_name], field_name::AbstractString, data) - #else - # problem.fields[field_name] = Field(data) - #end - update!(problem.elements, field_name::AbstractString, data) -end - -function haskey(problem::Problem, field_name::AbstractString) - return haskey(problem.fields, field_name) -end - -function getindex(problem::Problem, field_name::String) - return problem.fields[field_name] -end - -#""" Return field calculated to nodal points for elements in problem p. """ -function (problem::Problem)(field_name::String, time::Float64) - #if haskey(problem, field_name) - # return problem[field_name](time) - #end - f = Dict{Int,Any}() - for element in get_elements(problem) - haskey(element, field_name) || continue - for (c, v) in zip(get_connectivity(element), element(field_name, time)) - if haskey(f, c) - if !isapprox(f[c], v) - @info("several values for single node when returning field $field_name") - @info("already have: $(f[c]), and trying to set $v") - end - else - f[c] = v - end - end - end - #f == nothing && return f - #update!(problem, field_name, time => f) - return f -end - -# REMOVED: push!(problem, elements) - Use add_elements!(problem, elements) instead -# This violated Julia semantics (push! should be for collections, not domain logic) -# The modern Physics API uses add_elements!, add_dirichlet!, add_neumann! - -""" - set_gdofs!(problem, element) - -Set element global degrees of freedom. -""" -function set_gdofs!(problem, element, dofs) - problem.dofmap[element] = dofs -end - -""" - get_gdofs(problem, element) - -Return the global degrees of freedom for element. - -First make lookup from problem dofmap. If not defined there, make implicit -assumption that dofs follow formula `gdofs = [dim*(nid-1)+j for j=1:dim]`, -where `nid` is node id and `dim` is the dimension of problem. This formula -arranges dofs so that first comes all dofs of node 1, then node 2 and so on: -(u11, u12, u13, u21, u22, u23, ..., un1, un2, un3) for 3 dofs/node setting. -""" -function get_gdofs(problem::Problem, element::AbstractElement) - if haskey(problem.dofmap, element) - return problem.dofmap[element] - end - conn = get_connectivity(element) - if length(conn) == 0 - error("element connectivity not defined, cannot determine global ", - "degrees of freedom for element #: $(element.id)") - end - dim = get_unknown_field_dimension(problem) - gdofs = [dim * (i - 1) + j for i in conn for j = 1:dim] - return gdofs -end diff --git a/src/beams/api.jl b/src/beams/api.jl deleted file mode 100644 index e4eaafd..0000000 --- a/src/beams/api.jl +++ /dev/null @@ -1,98 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -""" -Beam formulation API definitions. - -This file defines beam-specific abstract types and formulation theories. -Must be included after core api.jl. -""" - -# ============================================================================ -# BEAM FORMULATION THEORIES -# ============================================================================ - -""" - AbstractBeamTheory - -Abstract type for beam theory variants. - -Beam theories differ in how they model shear deformation and cross-section kinematics. - -# Concrete Theories -- `EulerBernoulli`: Classical beam theory (no shear deformation) -- `Timoshenko`: Includes shear deformation (thick beams) - -# See Also -- [`BeamFormulation`](@ref) -""" -abstract type AbstractBeamTheory end - -""" - EulerBernoulli <: AbstractBeamTheory - -Euler-Bernoulli beam theory (classical, no shear deformation). - -Assumptions: -- Plane sections remain plane and perpendicular to neutral axis -- No transverse shear deformation -- Valid for slender beams (L/h > 10) - -# Usage -```julia -formulation = BeamFormulation{EulerBernoulli}() -physics = Physics( - formulation=formulation, - field=DisplacementRotation{3}(), - mesh=beam_mesh, - material=steel -) -``` -""" -struct EulerBernoulli <: AbstractBeamTheory end - -""" - Timoshenko <: AbstractBeamTheory - -Timoshenko beam theory (includes shear deformation). - -Assumptions: -- Plane sections remain plane but NOT perpendicular to neutral axis -- Transverse shear deformation included -- Valid for thick beams and higher frequencies - -# Usage -```julia -formulation = BeamFormulation{Timoshenko}() -physics = Physics( - formulation=formulation, - field=DisplacementRotation{3}(), - mesh=beam_mesh, - material=steel -) -``` -""" -struct Timoshenko <: AbstractBeamTheory end - -""" - BeamFormulation{Theory<:AbstractBeamTheory} <: AbstractFormulation - -Beam element formulation with theory variant. - -# Type Parameter -- `Theory`: Beam theory type (EulerBernoulli or Timoshenko) - -# Examples -```julia -# Slender beam (classical theory) -BeamFormulation{EulerBernoulli}() - -# Thick beam (includes shear) -BeamFormulation{Timoshenko}() -``` - -# Fields per Node -Typically 6 DOFs in 3D: (ux, uy, uz, θx, θy, θz) -Use with `DisplacementRotation{3}` field type. -""" -struct BeamFormulation{Theory<:AbstractBeamTheory} <: AbstractFormulation end diff --git a/src/domains/continuum/assemble_v2.jl b/src/domains/continuum/assemble_v2.jl deleted file mode 100644 index fcb01ee..0000000 --- a/src/domains/continuum/assemble_v2.jl +++ /dev/null @@ -1,522 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -""" -Assembly for 3D Continuum Elasticity - Ferrite-Style Two-Pointer Merge (V2) - -This is a CLEAN implementation using the Ferrite two-pointer merge algorithm -for 4.1x faster assembly compared to COO triplets. - -Key differences from continuum_3d.jl: -- Uses pre-built CSC structure (K_csc) instead of COO triplets -- Sorts element DOFs for linear-time merge -- Zero allocations in assembly loop -- 4.1x faster assembly (2.36ms vs 9.71ms for 2500 Tet4 elements) -- 16.6x less memory (506 KB vs 8.4 MB) - -References: -- Design: experiments/FERRITE_DATA_STRUCTURES_DESIGN.md -- Proof of concept: experiments/ferrite_style_assembly.jl -- Explanation: experiments/ferrite_sorteddofs_explained.jl -""" - -using Tensors -using SparseArrays -using LinearAlgebra - -# Import basevec for creating unit vectors -using Tensors: basevec - -# ============================================================================ -# Data Structures -# ============================================================================ - -""" - AssemblyCacheFerrite{N,T,Mat} - -Pre-allocated buffers for zero-allocation Ferrite-style assembly. - -Replaces COO triplets (I_rows, J_cols, K_values) with: -- Pre-built CSC structure (K_csc) - reused every assembly -- Sorted DOF buffers (sorteddofs, permutation) - for two-pointer merge - -# Type Parameters -- `N`: Number of nodes per element (from topology) -- `T`: Topology type (Hexahedron{8}, Tet10, etc.) -- `Mat`: Material type - -# Fields (Ferrite-specific) -- `K_csc`: Pre-built sparse matrix structure (zeros, reused) -- `sorteddofs`: Sorted element DOF buffer [3N] -- `permutation`: Sortperm buffer for DOF remapping [3N] - -# Fields (Standard) -- `f`: Global force vector -- `K_blocks`: Element blocked stiffness matrix [N×N of Tensor{2,3}] -- `K_e`: Element stiffness in Float64 matrix form [3N×3N] -- `X_buffer`: Element coordinate buffer [N of Vec{3}] -- `gdofs`: Global DOF indices buffer (unsorted) [3N] -- `elements`: Element IDs to assemble (Vector{UInt32}) -- `C`: Elasticity tensor (pre-computed) -- `topology`: Topology instance -- `basis`: Basis function instance -- `ips`: Integration points (pre-computed) - -# Performance -- Memory: ~506 KB for typical mesh (vs 8.4 MB COO triplets) -- Time: 2.36ms for 2500 Tet4 (vs 9.71ms COO) -- Allocations: Zero after warmup -""" -struct AssemblyCacheFerrite{N,T<:AbstractTopology{N},Mat<:AbstractMaterial} - # Ferrite-style CSC storage (NEW!) - K_csc::SparseMatrixCSC{Float64,Int} - sorteddofs::Vector{Int} - permutation::Vector{Int} - - # Global force vector - f::Vector{Float64} - - # Element assembly buffers (same as V1) - K_blocks::Matrix{Tensor{2,3,Float64,9}} - K_e::Matrix{Float64} - X_buffer::Vector{Vec{3,Float64}} - gdofs::Vector{Int} # Unsorted DOFs - - # Element set to assemble - elements::Vector{UInt32} - - # Pre-computed material/topology data - C::Tensor{4,3,Float64,81} - topology::T - basis::Lagrange{1} # Basis order only (topology passed separately) - ips::Any # Integration points tuple -end - -# ============================================================================ -# Sparsity Pattern Construction -# ============================================================================ - -""" -Build sparsity pattern from mesh connectivity. - -Returns (I, J) vectors for sparse matrix construction. -""" -function build_sparsity_pattern(mesh::M) where {M<:AbstractMesh} - # Get N from mesh type parameters - N = typeof(mesh).parameters[1] - - ndofs_total = 3 * length(mesh.nodes) - - # Pre-count entries (3N × 3N per element) - ndofs_per_elem = 3 * N - capacity = length(mesh.connectivity) * ndofs_per_elem * ndofs_per_elem - - I = Vector{Int}() - J = Vector{Int}() - sizehint!(I, capacity) - sizehint!(J, capacity) - - # Loop over elements - for conn in mesh.connectivity - # Global DOFs for this element - elem_dofs = Int[] - sizehint!(elem_dofs, ndofs_per_elem) - - for node_id in conn - for α in 1:3 - push!(elem_dofs, 3 * (node_id - 1) + α) - end - end - - # All pairs (i,j) in elem_dofs - for i in elem_dofs - for j in elem_dofs - push!(I, i) - push!(J, j) - end - end - end - - return I, J -end - -# ============================================================================ -# Cache Construction -# ============================================================================ - -""" - AssemblyCacheFerrite(physics::Physics{ContinuumFormulation{FullThreeD}, - Displacement{3}, M, Mat}) - -Construct Ferrite-style assembly cache with pre-built CSC structure. - -This is where ALL allocations happen. After construction, assembly is zero-allocation. - -# Key Steps -1. Build sparsity pattern from mesh connectivity -2. Create K_csc with sparse(I, J, ones(...)) -3. Zero K_csc.nzval for reuse -4. Allocate sorteddofs and permutation buffers -5. Pre-compute material and topology data - -# Performance -- One-time cost: ~10-20ms for typical mesh -- Pays off after ~1-2 assemblies vs COO method -""" -function AssemblyCacheFerrite( - physics::Physics{ContinuumFormulation{FullThreeD}, - Displacement{3}, - M, - Mat}) where {M<:AbstractMesh,Mat<:AbstractMaterial} - - mesh = physics.mesh - material = physics.material - element_set = physics.element_set - - # Get topology type and N from mesh type parameters - N_param = typeof(mesh).parameters[1] # N (8 for Hex8) - T = typeof(mesh).parameters[2] # Hexahedron{8} - - # Global system dimensions - nnodes = length(mesh.nodes) - ndofs = 3 * nnodes - - # Element set - elem_set = get_element_set(mesh, element_set) - elements = collect(elem_set) - - # Build sparsity pattern ONCE - I, J = build_sparsity_pattern(mesh) - K_csc = sparse(I, J, ones(length(I)), ndofs, ndofs) - fill!(K_csc.nzval, 0.0) # Zero values for reuse - - # Ferrite buffers for sorting DOFs - max_ndofs = 3 * N_param - sorteddofs = Vector{Int}(undef, max_ndofs) - permutation = Vector{Int}(undef, max_ndofs) - - # Global force vector - f = zeros(ndofs) - - # Element assembly buffers (same as V1) - max_nnodes = N_param - K_blocks = Matrix{Tensor{2,3,Float64,9}}(undef, max_nnodes, max_nnodes) - K_e = zeros(max_ndofs, max_ndofs) - X_buffer = Vector{Vec{3,Float64}}(undef, max_nnodes) - gdofs = Vector{Int}(undef, max_ndofs) - - # Pre-compute material and topology data - C = elasticity_tensor(material) - topology = T() - basis = Lagrange{1}() # Basis order only (topology passed separately) - integration_scheme = default_integration(T) - ips = integration_points(integration_scheme, topology) - - return AssemblyCacheFerrite{N_param,T,Mat}( - K_csc, # Pre-built structure - sorteddofs, # Sorted DOF buffer - permutation, # Sortperm buffer - f, - K_blocks, K_e, X_buffer, gdofs, - elements, - C, topology, basis, ips - ) -end - -# ============================================================================ -# Ferrite Two-Pointer Merge Assembly -# ============================================================================ - -""" - assemble_elements_ferrite!(cache::AssemblyCacheFerrite, mesh::M) - -Assemble elements using Ferrite two-pointer merge (ZERO allocations). - -This is the core Ferrite algorithm: -1. Zero K_csc.nzval once at start -2. For each element: - a. Compute element stiffness K_e - b. Get global DOFs (unsorted) - c. Sort DOFs → sorteddofs, permutation - d. Two-pointer merge: Linear scan through sorted lists - e. Accumulate to K_csc using permutation for correct K_e indices - -# Algorithm Detail -For each column i_global in element DOFs: -- Get CSC column range: K_csc.colptr[i_global]:(colptr[i_global+1]-1) -- K_csc.rowval[range] is SORTED -- sorteddofs is SORTED -- Two pointers: Ri (CSC), ri (element) -- Advance Ri until K_csc.rowval[Ri] == sorteddofs[ri] -- Accumulate: K_csc.nzval[Ri] += K_e[permutation[ri], permutation[i_local]] - -# Performance -- Time: 2.36ms for 2500 Tet4 elements -- Allocations: Zero after warmup -- 4.1x faster than COO method -- Linear-time merge vs O(log n) binary search -""" -function assemble_elements_ferrite!( - cache::AssemblyCacheFerrite{N,T,Mat}, - mesh::M) where {M<:AbstractMesh,N,T,Mat} - - # Zero K_csc once at start (reuse structure!) - fill!(cache.K_csc.nzval, 0.0) - - ndofs_elem = 3 * N - - # Loop over elements (ZERO allocations target!) - @inbounds for i in eachindex(cache.elements) - elem_id = cache.elements[i] - @inbounds conn = mesh.connectivity[elem_id] - - # 1. Fill coordinate buffer (in-place) - @inbounds for j in 1:N - cache.X_buffer[j] = mesh.nodes[conn[j]] - end - - # 2. Compute element stiffness - fill!(cache.K_blocks, zero(Tensor{2,3})) - compute_element_stiffness!(cache.K_blocks, cache.X_buffer, - cache.C, cache.topology, cache.basis, cache.ips) - blocked_tensor_to_matrix!(cache.K_e, cache.K_blocks) - - # 3. Global DOFs (UNSORTED, follows connectivity) - @inbounds for (local_idx, node_id) in enumerate(conn) - for α in 1:3 - cache.gdofs[3*(local_idx-1)+α] = 3 * (node_id - 1) + α - end - end - - # 4. Sort DOFs (sortperm! is in-place, zero allocation) - sortperm!(cache.permutation, cache.gdofs) - @inbounds for k in 1:ndofs_elem - cache.sorteddofs[k] = cache.gdofs[cache.permutation[k]] - end - - # 5. TWO-POINTER MERGE (Ferrite method!) - for i_local in 1:ndofs_elem - i_global = cache.sorteddofs[i_local] - - # Column range in K_csc for column i_global - col_start = cache.K_csc.colptr[i_global] - col_end = cache.K_csc.colptr[i_global+1] - 1 - - # Two pointers: Ri (CSC), ri (element) - Ri = col_start - for ri in 1:ndofs_elem - row_i_sorted = cache.sorteddofs[ri] - - # Advance Ri until K_csc.rowval[Ri] >= row_i_sorted - while Ri <= col_end && cache.K_csc.rowval[Ri] < row_i_sorted - Ri += 1 - end - - # If found, accumulate - if Ri <= col_end && cache.K_csc.rowval[Ri] == row_i_sorted - # Use permutation to get correct K_e indices! - orig_row = cache.permutation[ri] - orig_col = cache.permutation[i_local] - cache.K_csc.nzval[Ri] += cache.K_e[orig_row, orig_col] - end - end - end - end - - nothing -end - -# ============================================================================ -# Assembly Functions (User-Facing API) -# ============================================================================ - -""" - assemble_v2!(physics::Physics{ContinuumFormulation{FullThreeD}, - Displacement{3}, M, Mat}) -> (K, f) - -Assemble using Ferrite two-pointer merge method (V2). - -This is the user-facing function that: -1. Creates AssemblyCacheFerrite (allocates all buffers, builds CSC structure) -2. Calls assemble_elements_ferrite! (zero allocations) -3. Applies boundary conditions -4. Returns (K, f) - -# Performance -- 4.1x faster assembly than V1 (COO method) -- 16.6x less memory -- Zero allocations in hot loop - -# Example -```julia -physics = Physics(...) -K, f = assemble_v2!(physics) # Ferrite method -u = K \\ f -``` -""" -function assemble_v2!( - physics::Physics{ContinuumFormulation{FullThreeD}, - Displacement{3}, - M, - Mat}) where {M<:AbstractMesh,Mat<:AbstractMaterial} - - # Create cache (ALL allocations here!) - cache = AssemblyCacheFerrite(physics) - - # Assemble (ZERO allocations!) - return _assemble_ferrite!(physics, cache) -end - -""" - _assemble_ferrite!(physics, cache::AssemblyCacheFerrite) -> (K, f) - -Zero-allocation assembly using Ferrite cache. - -This is the internal function that performs the actual assembly. -Use `assemble_v2!` for the public API. -""" -function _assemble_ferrite!( - physics::Physics{ContinuumFormulation{FullThreeD}, - Displacement{3}, - M, - Mat}, - cache::AssemblyCacheFerrite{N,T,Mat2}) where {M<:AbstractMesh,Mat<:AbstractMaterial,N,T,Mat2} - - mesh = physics.mesh - bc_dirichlet = physics.bc_dirichlet - bc_neumann = physics.bc_neumann - - # Get dimensions - nnodes = length(mesh.nodes) - ndofs = 3 * nnodes - - # Clear force vector (in-place, zero allocation) - fill!(cache.f, 0.0) - - # Ferrite assembly (ZERO allocations!) - assemble_elements_ferrite!(cache, mesh) - - # Apply Neumann BCs (add forces to f) - for (surf_id, force) in zip(bc_neumann.surface_ids, bc_neumann.values) - # For now, interpret surface_ids as node_ids (simplified) - node = surf_id - if node <= nnodes - for α in 1:3 - cache.f[3*(node-1)+α] += force[α] - end - end - end - - # Copy K_csc (structure already correct, just copy!) - K = copy(cache.K_csc) - - # Apply Dirichlet BCs (modify K and f) - for i in 1:length(bc_dirichlet.node_ids) - node = bc_dirichlet.node_ids[i] - components = bc_dirichlet.components[i] - value = bc_dirichlet.values[i] - - for comp in components - dof = 3 * (node - 1) + comp - if dof <= ndofs # Safety check - K[dof, :] .= 0.0 - K[:, dof] .= 0.0 - K[dof, dof] = 1.0 - cache.f[dof] = value - end - end - end - - return (K, cache.f) -end - -# ============================================================================ -# Element Stiffness Computation (Reuse from V1) -# ============================================================================ - -""" - compute_stiffness_block(grad_k, grad_l, C) -> Tensor{2,3} - -Compute single 3×3 stiffness block (reused from continuum_3d.jl). - -See continuum_3d.jl for detailed documentation. -""" -@inline function compute_stiffness_block( - grad_k::Vec{3,Float64}, - grad_l::Vec{3,Float64}, - C::Tensor{4,3,Float64,81} -)::Tensor{2,3,Float64,9} - - K_kl = zero(Tensor{2,3}) - - for α in 1:3, β in 1:3 - e_α = basevec(Vec{3}, α) - e_β = basevec(Vec{3}, β) - B_k_α = 0.5 * (grad_k ⊗ e_α + e_α ⊗ grad_k) - B_l_β = 0.5 * (grad_l ⊗ e_β + e_β ⊗ grad_l) - k_αβ = dcontract(B_k_α, dcontract(C, B_l_β)) - K_kl += k_αβ * (e_α ⊗ e_β) - end - - return K_kl -end - -""" - blocked_tensor_to_matrix!(K_e, K_blocks) - -Convert blocked tensor to Float64 matrix (reused from continuum_3d.jl). -""" -function blocked_tensor_to_matrix!( - K_e::AbstractMatrix{Float64}, - K_blocks::AbstractMatrix{Tensor{2,3,Float64,9}}) - - nnodes = size(K_blocks, 1) - for i in 1:nnodes, j in 1:nnodes - for α in 1:3, β in 1:3 - K_e[3*(i-1)+α, 3*(j-1)+β] = K_blocks[i, j][α, β] - end - end - nothing -end - -""" - compute_element_stiffness!(K_blocks, X, C, topology, basis, ips) - -Compute element stiffness (reused from continuum_3d.jl). - -See continuum_3d.jl for detailed documentation. -""" -function compute_element_stiffness!( - K_blocks::AbstractMatrix{Tensor{2,3,Float64,9}}, - X::Vector{Vec{3,Float64}}, - C::Tensor{4,3,Float64,81}, - topology::T, - basis::B, - ips) where {T<:AbstractTopology{N},B<:AbstractBasis} where N - - for k in 1:N, l in 1:N - for ip in ips - ξ = ip.ξ - w = ip.weight - - dN_dξ = get_basis_derivatives(topology, basis, ξ) - - J = X[1] ⊗ dN_dξ[1] - for i in 2:N - J += X[i] ⊗ dN_dξ[i] - end - detJ = det(J) - J_inv = inv(J) - J_inv_T = transpose(J_inv) - - grad_k = J_inv_T ⋅ dN_dξ[k] - grad_l = J_inv_T ⋅ dN_dξ[l] - - K_kl = compute_stiffness_block(grad_k, grad_l, C) - - K_blocks[k, l] += K_kl * detJ * w - end - end - - nothing -end diff --git a/src/element_assembly_structures.jl b/src/element_assembly_structures.jl deleted file mode 100644 index da43c25..0000000 --- a/src/element_assembly_structures.jl +++ /dev/null @@ -1,341 +0,0 @@ -# Traditional Element Assembly -# -# This module provides the standard element-by-element assembly approach -# for comparison with nodal assembly. Builds global tangent stiffness matrix -# and residual force vector using sparse matrix formats. - -using Tensors -using SparseArrays -using LinearAlgebra - -""" - ElementAssemblyData{T} - -Storage for element assembly using traditional (element-by-element) approach. - -# Fields -- `K_global::SparseMatrixCSC{T}`: Global tangent stiffness matrix -- `r_global::Vector{T}`: Global residual force vector (r = f_int - f_ext) -- `f_int_global::Vector{T}`: Global internal force vector -- `f_ext_global::Vector{T}`: Global external force vector -- `ndof::Int`: Total number of degrees of freedom - -# Notes -- Assembly uses COO (coordinate) format, then converts to CSC -- Multiple elements can write to same global DOF (summed automatically) -""" -mutable struct ElementAssemblyData{T} - K_global::SparseMatrixCSC{T,Int} - r_global::Vector{T} - f_int_global::Vector{T} - f_ext_global::Vector{T} - ndof::Int -end - -""" - ElementAssemblyData(ndof::Int, ::Type{T}=Float64) - -Allocate storage for traditional element assembly. - -# Arguments -- `ndof`: Total degrees of freedom (nnodes × 3 for 3D) -- `T`: Floating point type (default Float64) - -# Example -```julia -nnodes = 100 -assembly = ElementAssemblyData(3 * nnodes, Float64) -``` -""" -function ElementAssemblyData(ndof::Int, ::Type{T}=Float64) where T - # Pre-allocate empty sparse matrix (will fill during assembly) - K_global = spzeros(T, ndof, ndof) - r_global = zeros(T, ndof) - f_int_global = zeros(T, ndof) - f_ext_global = zeros(T, ndof) - - return ElementAssemblyData{T}(K_global, r_global, f_int_global, f_ext_global, ndof) -end - -""" - reset!(assembly::ElementAssemblyData) - -Reset assembly data to zero (for incremental/iterative solvers). -""" -function reset!(assembly::ElementAssemblyData{T}) where T - assembly.K_global = spzeros(T, assembly.ndof, assembly.ndof) - fill!(assembly.r_global, 0.0) - fill!(assembly.f_int_global, 0.0) - fill!(assembly.f_ext_global, 0.0) -end - -""" - ElementContribution{T} - -Local element contribution before scattering to global. - -# Fields -- `element_id::Int`: Element ID -- `gdofs::Vector{Int}`: Global DOF indices (e.g., [1,2,3,4,5,6,...] for nodes) -- `K_local::Matrix{T}`: Local stiffness matrix (ndofs_local × ndofs_local) -- `f_int_local::Vector{T}`: Local internal force vector -- `f_ext_local::Vector{T}`: Local external force vector - -# Notes -- For Tet4: ndofs_local = 12 (4 nodes × 3 DOF) -- For Tet10: ndofs_local = 30 (10 nodes × 3 DOF) -""" -struct ElementContribution{T} - element_id::Int - gdofs::Vector{Int} - K_local::Matrix{T} - f_int_local::Vector{T} - f_ext_local::Vector{T} -end - -""" - ElementContribution(element_id::Int, gdofs::Vector{Int}, ::Type{T}=Float64) - -Allocate storage for element contribution. - -# Arguments -- `element_id`: Element ID -- `gdofs`: Global DOF indices -- `T`: Floating point type - -# Example -```julia -# Tet4 element connecting nodes [5, 7, 12, 15] -gdofs = [13,14,15, 19,20,21, 34,35,36, 43,44,45] # 3 DOF per node -contrib = ElementContribution(1, gdofs, Float64) -``` -""" -function ElementContribution(element_id::Int, gdofs::Vector{Int}, ::Type{T}=Float64) where T - ndofs = length(gdofs) - K_local = zeros(T, ndofs, ndofs) - f_int_local = zeros(T, ndofs) - f_ext_local = zeros(T, ndofs) - - return ElementContribution{T}(element_id, gdofs, K_local, f_int_local, f_ext_local) -end - -""" - scatter_to_global!(assembly::ElementAssemblyData, contrib::ElementContribution) - -Scatter element contribution to global matrices/vectors (traditional assembly). - -This is the key operation in element assembly: add local element quantities -to global system. Uses COO format (accumulates into lists). - -# Arguments -- `assembly`: Global assembly data -- `contrib`: Element contribution - -# Notes -- Multiple elements can contribute to same global DOF (summed) -- For GPU: Would require atomic operations (slow!) -- For CPU: Direct scatter-add works fine -""" -function scatter_to_global!(assembly::ElementAssemblyData{T}, - contrib::ElementContribution{T}) where T - # Scatter forces (simple vector addition) - for (local_i, global_i) in enumerate(contrib.gdofs) - assembly.f_int_global[global_i] += contrib.f_int_local[local_i] - assembly.f_ext_global[global_i] += contrib.f_ext_local[local_i] - end - - # Scatter stiffness (matrix addition) - # Build list of (I, J, V) triplets for sparse matrix - I_rows = Int[] - J_cols = Int[] - values = T[] - - ndofs_local = length(contrib.gdofs) - for i in 1:ndofs_local, j in 1:ndofs_local - if abs(contrib.K_local[i, j]) > 1e-14 # Skip near-zeros - push!(I_rows, contrib.gdofs[i]) - push!(J_cols, contrib.gdofs[j]) - push!(values, contrib.K_local[i, j]) - end - end - - # Add to existing sparse matrix - K_elem = sparse(I_rows, J_cols, values, assembly.ndof, assembly.ndof) - assembly.K_global += K_elem -end - -""" - compute_residual!(assembly::ElementAssemblyData) - -Compute residual force vector: r = f_int - f_ext - -Should be called after all elements have been assembled. -""" -function compute_residual!(assembly::ElementAssemblyData{T}) where T - assembly.r_global .= assembly.f_int_global .- assembly.f_ext_global -end - -""" - assemble_elements!(assembly::ElementAssemblyData, - contributions::Vector{ElementContribution}) - -Assemble all element contributions to global system. - -# Arguments -- `assembly`: Global assembly data (modified in-place) -- `contributions`: Vector of element contributions - -# Example -```julia -assembly = ElementAssemblyData(ndof) -contributions = compute_all_element_contributions(elements, u, time) -assemble_elements!(assembly, contributions) -compute_residual!(assembly) - -# Now solve: K_global * Δu = -r_global -``` -""" -function assemble_elements!(assembly::ElementAssemblyData{T}, - contributions::Vector{ElementContribution{T}}) where T - reset!(assembly) - - # Loop over elements and scatter (element assembly) - for contrib in contributions - scatter_to_global!(assembly, contrib) - end - - # Compute residual - compute_residual!(assembly) -end - -""" - apply_dirichlet_bc!(assembly::ElementAssemblyData, - fixed_dofs::Vector{Int}, - prescribed_values::Vector{T}=zeros(length(fixed_dofs))) - -Apply Dirichlet (essential) boundary conditions by penalty method. - -# Arguments -- `assembly`: Global assembly data (modified in-place) -- `fixed_dofs`: DOF indices to fix -- `prescribed_values`: Prescribed displacement values (default: zeros) - -# Method -Uses penalty method: adds large stiffness to diagonal and corresponding RHS. - -For DOF i with prescribed value u_prescribed: -- K[i,i] += penalty (e.g., 1e10 * max_K) -- r[i] = penalty * (u_current - u_prescribed) - -# Example -```julia -# Fix nodes 1 and 2 in all directions (zero displacement) -fixed_dofs = [1,2,3, 4,5,6] # Nodes 1,2 × 3 DOF -apply_dirichlet_bc!(assembly, fixed_dofs) -``` -""" -function apply_dirichlet_bc!(assembly::ElementAssemblyData{T}, - fixed_dofs::Vector{Int}, - prescribed_values::Vector{T}=zeros(T, length(fixed_dofs))) where T - # Penalty parameter (large relative to stiffness) - max_K = maximum(abs, assembly.K_global) - penalty = 1e10 * max_K - - for (idx, dof) in enumerate(fixed_dofs) - # Add penalty stiffness to diagonal - assembly.K_global[dof, dof] += penalty - - # Modify residual (assuming current displacement is zero for now) - # In full Newton: r[i] += penalty * (u_current[i] - u_prescribed[i]) - assembly.r_global[dof] = penalty * prescribed_values[idx] - end -end - -""" - get_dof_indices(connectivity::NTuple{N,Int}, dim::Int=3) -> Vector{Int} - -Get global DOF indices for an element given node connectivity. - -# Arguments -- `connectivity`: Element node IDs (e.g., (5, 7, 12, 15) for Tet4) -- `dim`: Dimension (3 for 3D elasticity) - -# Returns -- `gdofs::Vector{Int}`: Global DOF indices - -# Example -```julia -# Element with nodes [5, 7, 12, 15] -gdofs = get_dof_indices((5, 7, 12, 15), 3) -# Returns: [13,14,15, 19,20,21, 34,35,36, 43,44,45] -``` -""" -function get_dof_indices(connectivity::NTuple{N,Int}, dim::Int=3) where N - nnodes = length(connectivity) - gdofs = zeros(Int, dim * nnodes) - - for (local_i, global_node) in enumerate(connectivity) - for d in 1:dim - gdofs[dim*(local_i-1)+d] = dim * (global_node - 1) + d - end - end - - return gdofs -end - -""" - matrix_vector_product(assembly::ElementAssemblyData, v::Vector{T}) -> Vector{T} - -Compute matrix-vector product: w = K * v using assembled sparse matrix. - -# Arguments -- `assembly`: Assembly data (contains K_global) -- `v`: Input vector (ndof) - -# Returns -- `w`: Output vector w = K * v - -# Example -```julia -# GMRES matrix-free operator -function matvec(v) - return matrix_vector_product(assembly, v) -end -Δu = gmres(matvec, -r, tol=1e-6) -``` -""" -function matrix_vector_product(assembly::ElementAssemblyData{T}, v::Vector{T}) where T - return assembly.K_global * v -end - -""" - print_assembly_stats(assembly::ElementAssemblyData) - -Print statistics about assembled system (for debugging). -""" -function print_assembly_stats(assembly::ElementAssemblyData) - nnz_K = nnz(assembly.K_global) - ndof = assembly.ndof - fill_ratio = nnz_K / (ndof * ndof) - - println("="^60) - println("Traditional Element Assembly Statistics") - println("="^60) - println(" Total DOF: ", ndof) - println(" K matrix size: ", size(assembly.K_global)) - println(" K non-zeros: ", nnz_K) - println(" K fill ratio: ", round(fill_ratio, sigdigits=2)) - println(" K memory (MB): ", round(nnz_K * 16 / 1024^2, digits=2)) - println(" ||f_int||: ", round(norm(assembly.f_int_global), sigdigits=2)) - println(" ||f_ext||: ", round(norm(assembly.f_ext_global), sigdigits=2)) - println(" ||residual||: ", round(norm(assembly.r_global), sigdigits=2)) - println(" K symmetric: ", issymmetric(assembly.K_global)) - println("="^60) -end - -# Export main types and functions -export ElementAssemblyData, ElementContribution -export reset!, scatter_to_global!, compute_residual! -export assemble_elements!, apply_dirichlet_bc! -export get_dof_indices, matrix_vector_product -export print_assembly_stats diff --git a/src/elements/integrate.jl b/src/elements/integrate.jl deleted file mode 100644 index 1f0957e..0000000 --- a/src/elements/integrate.jl +++ /dev/null @@ -1,57 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE - -# Default number of integration points for each element. First rule is the -# default integration rule returned by `get_integration_points(element)`. -# Sometimes we want to increase integration order, e.g. when integrating mass -# matrix or boundary conditions. For that reason, additional rules are provied -# in list, so e.g. `get_integration_points(element, 1)` returns the second rule, -# `get_integration_points(element, 2)` third rule and so on. Rules should be -# ordered so that picking next one integrates more accurately. -integration_rule_mapping = ( - :Seg2 => (:GLSEG2, :GLSEG3, :GLSEG4, :GLSEG5), - :Seg3 => (:GLSEG3, :GLSEG4, :GLSEG5), - :NSeg => (:GLSEG2, :GLSEG3, :GLSEG4, :GLSEG5), - :Quad4 => (:GLQUAD4, :GLQUAD9, :GLQUAD16, :GLQUAD25), - :Quad8 => (:GLQUAD9, :GLQUAD16, :GLQUAD25), - :Quad9 => (:GLQUAD9, :GLQUAD16, :GLQUAD25), - :NSurf => (:GLQUAD9, :GLQUAD16, :GLQUAD25), - :Hex8 => (:GLHEX8, :GLHEX27, :GLHEX64, :GLHEX125), - :Hex20 => (:GLHEX27, :GLHEX64, :GLHEX125), - :Hex27 => (:GLHEX27, :GLHEX64, :GLHEX125), - :NSolid => (:GLHEX27, :GLHEX64, :GLHEX125), - :Tri3 => (:GLTRI1, :GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12), - :Tri6 => (:GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12), - :Tri7 => (:GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12), - :Tet4 => (:GLTET1, :GLTET4, :GLTET5, :GLTET15), - :Tet10 => (:GLTET4, :GLTET5, :GLTET15), - :Pyr5 => (:GLPYR5,), - :Wedge6 => (:GLWED6, :GLWED21), - :Wedge15 => (:GLWED21,)) - -for (E, R) in integration_rule_mapping - for i in 1:length(R) - P = Val{R[i]} - order = Val{i - 1} - local code # Explicitly declare as local to avoid warning - if isequal(i, 1) - code = quote - function get_integration_points(element::$E) - return get_quadrature_points($P) - end - end - else - code = quote - function get_integration_points(element::$E, ::Type{$order}) - return get_quadrature_points($P) - end - end - end - eval(code) - end -end - -# All good codes needs a special case. Here we have it: Poi1 -function get_integration_points(::Poi1) - [(1.0, (0.0,))] -end diff --git a/src/formulations/api.jl b/src/formulations/api.jl deleted file mode 100644 index cc8eb72..0000000 --- a/src/formulations/api.jl +++ /dev/null @@ -1,307 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -""" -Formulation API definitions. - -This file defines formulation abstractions - the mathematical discretization strategies -for different types of FEM problems. - -Must be included after fields/api.jl (formulations work with fields). -""" - -# ============================================================================ -# FORMULATION INTERFACE -# ============================================================================ - -""" - AbstractFormulation - -Abstract type for discretization formulations. - -Formulation defines HOW we discretize the governing equations. Different formulations -exist for different physics domains: - -- **Continuum formulations** (this file) - Standard FEM for solid/fluid mechanics -- **Beam formulations** (src/beams/api.jl) - 1D structural elements -- **Shell formulations** (src/shells/api.jl) - 2D structural elements -- **Truss formulations** (src/trusses/api.jl) - 1D axial elements - -# Type Hierarchy -- `ContinuumFormulation{Theory}` - Standard continuum FEM (here) -- `BeamFormulation{Theory}` - Beam elements (src/beams/api.jl) -- `ShellFormulation{Theory}` - Shell elements (src/shells/api.jl) -- `TrussFormulation{Theory}` - Truss elements (src/trusses/api.jl) - -# Design Philosophy - -**Formulation + Field = Dispatch pattern** - -The combination of formulation and field type determines: -- Assembly method dispatch -- Element stiffness computation -- Stress/strain tensor dimensions -- DOF coupling patterns - -# Examples - -```julia -# 3D solid mechanics -physics = Physics( - formulation = ContinuumFormulation{FullThreeD}(), - field = Displacement{3}(), - mesh = mesh, - material = steel -) - -# 2D plane stress -physics_2d = Physics( - formulation = ContinuumFormulation{PlaneStress}(), - field = Displacement{2}(), - mesh = mesh_2d, - material = aluminum -) - -# Beam structure -physics_beam = Physics( - formulation = BeamFormulation{Timoshenko}(), - field = DisplacementRotation{3}(), - mesh = beam_mesh, - material = steel -) -``` - -# Assembly Dispatch - -Specialized assembly methods dispatch on formulation × field: - -```julia -# 3D continuum mechanics -function assemble!(physics::Physics{ContinuumFormulation{FullThreeD}, Displacement{3}, M, Mat}) - # Standard 3D displacement-based assembly - # Implementation in src/assembly/continuum_3d.jl -end - -# 2D plane stress -function assemble!(physics::Physics{ContinuumFormulation{PlaneStress}, Displacement{2}, M, Mat}) - # 2D assembly with plane stress assumptions - # Implementation in src/assembly/continuum_2d.jl -end - -# Beam elements -function assemble!(physics::Physics{BeamFormulation{Timoshenko}, DisplacementRotation{3}, M, Mat}) - # Beam-specific assembly (6 DOFs per node) - # Implementation in src/assembly/beams.jl -end -``` - -# See Also -- Field types: src/fields/api.jl (Displacement, Temperature, DisplacementRotation) -- Physics coupling: src/physics/api.jl (AbstractPhysics) -- Domain-specific formulations: src/beams/api.jl, src/shells/api.jl, src/trusses/api.jl -- Assembly implementations: src/assembly/continuum_3d.jl, src/assembly/beams.jl, etc. -""" -abstract type AbstractFormulation end - -# ============================================================================ -# CONTINUUM FORMULATION (Standard FEM) -# ============================================================================ - -""" - AbstractContinuumTheory - -Theory variants for continuum formulation. - -Controls dimensionality reduction and stress/strain assumptions for continuum -mechanics problems. - -# Concrete Theories -- `FullThreeD` - Full 3D analysis (no simplifications) -- `PlaneStress` - 2D plane stress (σ_zz = 0, thin plates) -- `PlaneStrain` - 2D plane strain (ε_zz = 0, thick plates) -- `Axisymmetric` - Axisymmetric analysis (rotation around z-axis) - -# Theory Selection Guidelines - -**FullThreeD (σ_xx, σ_yy, σ_zz, σ_xy, σ_yz, σ_xz):** -- General 3D solid mechanics -- No simplifying assumptions -- Most accurate but most expensive - -**PlaneStress (σ_xx, σ_yy, σ_xy, σ_zz = 0):** -- Thin plates and membranes (thickness << length/width) -- Out-of-plane stress σ_zz = 0 -- Examples: Sheet metal, aircraft skin, thin-walled structures - -**PlaneStrain (ε_xx, ε_yy, ε_xy, ε_zz = 0):** -- Thick sections with no variation in z-direction -- Out-of-plane strain ε_zz = 0 -- Examples: Dams, tunnels, retaining walls, long cylinders - -**Axisymmetric (σ_rr, σ_θθ, σ_zz, σ_rz):** -- Geometry and loading symmetric about z-axis -- No circumferential variations -- Examples: Pressure vessels, pipes, rotating disks - -# Usage - -```julia -# Full 3D solid mechanics -formulation = ContinuumFormulation{FullThreeD}() - -# 2D plane stress (thin plate) -formulation = ContinuumFormulation{PlaneStress}() - -# 2D plane strain (thick section) -formulation = ContinuumFormulation{PlaneStrain}() - -# Axisymmetric (cylinder, sphere) -formulation = ContinuumFormulation{Axisymmetric}() -``` - -# Mathematical Details - -**Plane Stress (thin plate):** -- Stress state: σ_zz = σ_xz = σ_yz = 0 -- Strain: ε_zz ≠ 0 (computed from σ_zz = 0 condition) -- Constitutive: 3×3 reduced stiffness matrix - -**Plane Strain (thick section):** -- Strain state: ε_zz = γ_xz = γ_yz = 0 -- Stress: σ_zz ≠ 0 (computed from ε_zz = 0 condition) -- Constitutive: 3×3 reduced stiffness matrix (different from plane stress!) - -**Axisymmetric:** -- Cylindrical coordinates (r, θ, z) -- No ∂/∂θ terms (axial symmetry) -- 4 stress components: σ_rr, σ_θθ, σ_zz, σ_rz -- Hoop stress σ_θθ from radial displacement - -# See Also -- [`ContinuumFormulation`](@ref) - Formulation struct using these theories -""" -abstract type AbstractContinuumTheory end - -""" - FullThreeD <: AbstractContinuumTheory - -Full 3D analysis with no simplifications. - -All six stress components: σ_xx, σ_yy, σ_zz, σ_xy, σ_yz, σ_xz -""" -struct FullThreeD <: AbstractContinuumTheory end - -""" - PlaneStress <: AbstractContinuumTheory - -2D plane stress assumption (σ_zz = 0). - -Applicable to thin plates and membranes where thickness << in-plane dimensions. -""" -struct PlaneStress <: AbstractContinuumTheory end - -""" - PlaneStrain <: AbstractContinuumTheory - -2D plane strain assumption (ε_zz = 0). - -Applicable to thick sections with no variation in z-direction. -""" -struct PlaneStrain <: AbstractContinuumTheory end - -""" - Axisymmetric <: AbstractContinuumTheory - -Axisymmetric analysis (rotation around z-axis). - -Geometry and loading symmetric about z-axis with no circumferential variations. -""" -struct Axisymmetric <: AbstractContinuumTheory end - -""" - ContinuumFormulation{Theory} <: AbstractFormulation - -Standard continuum mechanics formulation with theory variant. - -This is the fundamental FEM formulation for solid mechanics, heat transfer, -and other continuum physics problems. - -# Type Parameter -- `Theory <: AbstractContinuumTheory` - Dimensionality/simplification theory - -# Examples - -```julia -# 3D elasticity -physics = Physics( - formulation = ContinuumFormulation{FullThreeD}(), - field = Displacement{3}(), - mesh = mesh, - material = steel -) - -# 2D plane stress (thin plate) -physics_2d = Physics( - formulation = ContinuumFormulation{PlaneStress}(), - field = Displacement{2}(), - mesh = mesh_2d, - material = aluminum -) - -# 2D plane strain (thick section) -physics_2d = Physics( - formulation = ContinuumFormulation{PlaneStrain}(), - field = Displacement{2}(), - mesh = mesh_2d, - material = concrete -) - -# Axisymmetric (cylinder) -physics_axisym = Physics( - formulation = ContinuumFormulation{Axisymmetric}(), - field = Displacement{2}(), # (r, z) displacements - mesh = mesh_2d, - material = steel -) -``` - -# Assembly Dispatch - -Assembly methods specialize on theory × field combinations: - -```julia -# 3D solid mechanics -function assemble!(physics::Physics{ContinuumFormulation{FullThreeD}, Displacement{3}, M, Mat}) - # Standard 3D displacement-based assembly - # Full 6×6 strain-displacement matrix (Bε) - # 6×6 constitutive matrix (Dε) -end - -# 2D plane stress -function assemble!(physics::Physics{ContinuumFormulation{PlaneStress}, Displacement{2}, M, Mat}) - # 2D assembly with plane stress assumptions - # 3×3 reduced strain-displacement matrix - # 3×3 plane stress constitutive matrix -end - -# Heat transfer (same formulation, different field!) -function assemble!(physics::Physics{ContinuumFormulation{FullThreeD}, Temperature, M, Mat}) - # Thermal assembly (scalar field) - # Thermal conductivity matrix -end -``` - -# Implementation Location - -Concrete assembly implementations are in: -- `src/assembly/continuum_3d.jl` - 3D continuum mechanics -- `src/assembly/continuum_2d.jl` - 2D plane stress/strain -- `src/assembly/axisymmetric.jl` - Axisymmetric problems - -# See Also -- [`AbstractContinuumTheory`](@ref) - Theory variants -- Field types: src/fields/api.jl (Displacement, Temperature) -- Physics coupling: src/physics/api.jl (AbstractPhysics) -- Assembly: src/assembly/continuum_*.jl -""" -struct ContinuumFormulation{Theory<:AbstractContinuumTheory} <: AbstractFormulation end diff --git a/src/gpu_elasticity.jl b/src/gpu_elasticity.jl deleted file mode 100644 index 3d6a598..0000000 --- a/src/gpu_elasticity.jl +++ /dev/null @@ -1,476 +0,0 @@ -""" -Main solver: Elasticity on GPU - -Solves linear elasticity using: -- Two-phase nodal assembly (no atomics) -- Matrix-free conjugate gradient -- GPU-resident throughout -""" -function solve_elasticity_gpu(physics::ElasticityPhysics; tol=1e-6, max_iter=1000) - -module GPUElasticity - -export solve_elasticity_gpu, ElasticityPhysics, ElasticMaterial - -using CUDA -using Tensors -using LinearAlgebra -using Printf - -# Re-export mesh reader -include("gmsh_reader.jl") -using .GmshReader -export read_gmsh_mesh, GmshMesh, get_surface_nodes - -""" -Elastic material properties -""" -struct ElasticMaterial - E::Float64 # Young's modulus [Pa] - ν::Float64 # Poisson's ratio [-] -end - -""" -Elasticity physics definition -""" -struct ElasticityPhysics - mesh::GmshMesh - material::ElasticMaterial - fixed_nodes::Vector{Int} # Dirichlet BC (fixed displacement) - pressure_nodes::Vector{Int} # Neumann BC (pressure load) - pressure_value::Float64 # Pressure magnitude [Pa] -end - -""" -Node-to-elements connectivity (CSR format) -""" -struct NodeToElementsMap - ptr::CuArray{Int32,1} - data::CuArray{Int32,1} -end - -""" -Build CSR map: which elements touch each node? -""" -function build_node_to_elems_gpu(elements::Matrix{Int}, n_nodes::Int) - # Count connections per node - counts = zeros(Int, n_nodes) - for elem_idx in 1:size(elements, 2) - for i in 1:4 - node = elements[i, elem_idx] - counts[node] += 1 - end - end - - # Build CSR structure - ptr = cumsum([1; counts]) - data = Vector{Int32}(undef, sum(counts)) - - # Fill data array - offset = copy(ptr[1:end-1]) - for elem_idx in 1:size(elements, 2) - for i in 1:4 - node = elements[i, elem_idx] - data[offset[node]] = elem_idx - offset[node] += 1 - end - end - - return NodeToElementsMap(CuArray(Int32.(ptr)), CuArray(data)) -end - -""" -PHASE 1 GPU KERNEL: Compute element stiffness contributions at integration points - -For LINEAR ELASTICITY (no plasticity), we don't need state variables. -Just compute stresses from strains using Hooke's law. -""" -function compute_element_stresses_kernel!( - σ_gp::CuDeviceArray{SymmetricTensor{2,3,Float64,6},1}, - u::CuDeviceArray{Float64,1}, - nodes::CuDeviceArray{Float64,2}, - elements::CuDeviceArray{Int32,2}, - E, ν -) - gp_idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x - - if gp_idx <= length(σ_gp) - # Map GP to element - elem_idx = (gp_idx - 1) ÷ 4 + 1 # 4 GPs per Tet4 - - # Extract element nodes - n1 = elements[1, elem_idx] - n2 = elements[2, elem_idx] - n3 = elements[3, elem_idx] - n4 = elements[4, elem_idx] - - # Node coordinates - X1 = Vec{3}((nodes[1, n1], nodes[2, n1], nodes[3, n1])) - X2 = Vec{3}((nodes[1, n2], nodes[2, n2], nodes[3, n2])) - X3 = Vec{3}((nodes[1, n3], nodes[2, n3], nodes[3, n3])) - X4 = Vec{3}((nodes[1, n4], nodes[2, n4], nodes[3, n4])) - - # Displacements - u1 = Vec{3}((u[3*n1-2], u[3*n1-1], u[3*n1])) - u2 = Vec{3}((u[3*n2-2], u[3*n2-1], u[3*n2])) - u3 = Vec{3}((u[3*n3-2], u[3*n3-1], u[3*n3])) - u4 = Vec{3}((u[3*n4-2], u[3*n4-1], u[3*n4])) - - # Shape derivatives (constant for Tet4) - dN1_dxi = Vec{3}((-1.0, -1.0, -1.0)) - dN2_dxi = Vec{3}((1.0, 0.0, 0.0)) - dN3_dxi = Vec{3}((0.0, 1.0, 0.0)) - dN4_dxi = Vec{3}((0.0, 0.0, 1.0)) - - # Jacobian - J = dN1_dxi ⊗ X1 + dN2_dxi ⊗ X2 + dN3_dxi ⊗ X3 + dN4_dxi ⊗ X4 - invJ = inv(J) - - # Physical derivatives - dN1_dx = invJ ⋅ dN1_dxi - dN2_dx = invJ ⋅ dN2_dxi - dN3_dx = invJ ⋅ dN3_dxi - dN4_dx = invJ ⋅ dN4_dxi - - # Strain (small strain assumption) - ε = symmetric(dN1_dx ⊗ u1 + dN2_dx ⊗ u2 + dN3_dx ⊗ u3 + dN4_dx ⊗ u4) - - # Stress (Hooke's law) - λ = E * ν / ((1 + ν) * (1 - 2ν)) - μ = E / (2(1 + ν)) - I = one(ε) - σ = λ * tr(ε) * I + 2μ * ε - - # Store result - σ_gp[gp_idx] = σ - end - - return nothing -end - -""" -PHASE 2 GPU KERNEL: Nodal assembly (matrix-free, no atomics!) -""" -function nodal_assembly_kernel!( - r::CuDeviceArray{Float64,1}, - σ_gp::CuDeviceArray{SymmetricTensor{2,3,Float64,6},1}, - nodes::CuDeviceArray{Float64,2}, - elements::CuDeviceArray{Int32,2}, - node_to_elems_ptr::CuDeviceArray{Int32,1}, - node_to_elems_data::CuDeviceArray{Int32,1} -) - node_idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x - - if node_idx <= size(nodes, 2) - # Accumulate forces - f_node = zero(Vec{3,Float64}) - - # Gauss weight for Tet4 - gauss_weight = 1.0 / 24.0 - - # Shape derivatives - dN_dxi = ( - Vec{3}((-1.0, -1.0, -1.0)), - Vec{3}((1.0, 0.0, 0.0)), - Vec{3}((0.0, 1.0, 0.0)), - Vec{3}((0.0, 0.0, 1.0)) - ) - - # Get element range for this node - elem_start = node_to_elems_ptr[node_idx] - elem_end = node_to_elems_ptr[node_idx+1] - 1 - - # Loop over touching elements - for elem_offset in elem_start:elem_end - elem_idx = node_to_elems_data[elem_offset] - - # Extract element nodes - n1 = elements[1, elem_idx] - n2 = elements[2, elem_idx] - n3 = elements[3, elem_idx] - n4 = elements[4, elem_idx] - - # Find local node index - local_node = 1 - if node_idx == n2 - local_node = 2 - elseif node_idx == n3 - local_node = 3 - elseif node_idx == n4 - local_node = 4 - end - - # Recompute geometry (matrix-free!) - X1 = Vec{3}((nodes[1, n1], nodes[2, n1], nodes[3, n1])) - X2 = Vec{3}((nodes[1, n2], nodes[2, n2], nodes[3, n2])) - X3 = Vec{3}((nodes[1, n3], nodes[2, n3], nodes[3, n3])) - X4 = Vec{3}((nodes[1, n4], nodes[2, n4], nodes[3, n4])) - - J = dN_dxi[1] ⊗ X1 + dN_dxi[2] ⊗ X2 + dN_dxi[3] ⊗ X3 + dN_dxi[4] ⊗ X4 - detJ = det(J) - invJ = inv(J) - - # Physical derivative for this node - dN_dx = invJ ⋅ dN_dxi[local_node] - - # Loop over Gauss points (4 per Tet4) - for local_gp in 1:4 - gp_idx = (elem_idx - 1) * 4 + local_gp - σ = σ_gp[gp_idx] - - # Accumulate force - f_node += (dN_dx ⋅ σ) * (gauss_weight * detJ) - end - end - - # Write result (no atomics!) - r[3*node_idx-2] = f_node[1] - r[3*node_idx-1] = f_node[2] - r[3*node_idx] = f_node[3] - end - - return nothing -end - -""" -Compute residual on GPU (internal forces) -""" -function compute_residual_gpu!( - r::CuArray{Float64,1}, - u::CuArray{Float64,1}, - nodes::CuArray{Float64,2}, - elements::CuArray{Int32,2}, - node_to_elems::NodeToElementsMap, - E, ν -) - n_gp = size(elements, 2) * 4 - n_nodes = size(nodes, 2) - - # Phase 1: Compute stresses at GPs - σ_gp = CuArray{SymmetricTensor{2,3,Float64,6}}(undef, n_gp) - - threads = 256 - blocks = cld(n_gp, threads) - @cuda threads = threads blocks = blocks compute_element_stresses_kernel!( - σ_gp, u, nodes, elements, E, ν - ) - - # Phase 2: Nodal assembly - fill!(r, 0.0) - - threads = 256 - blocks = cld(n_nodes, threads) - @cuda threads = threads blocks = blocks nodal_assembly_kernel!( - r, σ_gp, nodes, elements, - node_to_elems.ptr, node_to_elems.data - ) - - return r -end - -""" -Apply pressure load to top surface (Neumann BC) -""" -function apply_pressure_load!( - f::CuArray{Float64,1}, - pressure_nodes::Vector{Int}, - mesh::GmshMesh, - pressure::Float64 -) - # Simple uniform distribution (should integrate properly over surface) - # For now, divide pressure equally among nodes - - f_cpu = Array(f) - n_pressure_nodes = length(pressure_nodes) - - # Estimate surface area (assuming uniform Z = height) - surface_area = (maximum(mesh.nodes[1, :]) - minimum(mesh.nodes[1, :])) * - (maximum(mesh.nodes[2, :]) - minimum(mesh.nodes[2, :])) - - # Total force - total_force = pressure * surface_area - force_per_node = total_force / n_pressure_nodes - - # Apply in Z direction (negative, pointing down) - for node in pressure_nodes - f_cpu[3*node] += -force_per_node # Z component - end - - copyto!(f, f_cpu) - - return f -end - -""" -Apply Dirichlet boundary conditions (fixed nodes) -""" -function apply_dirichlet_bc!( - K_op::Function, - f::CuArray{Float64,1}, - fixed_nodes::Vector{Int} -) - # Zero out DOFs - f_cpu = Array(f) - for node in fixed_nodes - f_cpu[3*node-2] = 0.0 # X - f_cpu[3*node-1] = 0.0 # Y - f_cpu[3*node] = 0.0 # Z - end - copyto!(f, f_cpu) - - # Return modified operator that zeros fixed DOFs - function K_bc(u) - r = K_op(u) - r_cpu = Array(r) - for node in fixed_nodes - r_cpu[3*node-2] = 0.0 - r_cpu[3*node-1] = 0.0 - r_cpu[3*node] = 0.0 - end - copyto!(r, r_cpu) - return r - end - - return K_bc -end - -""" -Conjugate Gradient solver (GPU) -""" -function cg_solve_gpu!( - x::CuArray{Float64,1}, - A_op::Function, - b::CuArray{Float64,1}; - tol=1e-6, - max_iter=1000 -) - n = length(x) - - # Initial residual - r = b - A_op(x) - p = copy(r) - rsold = dot(r, r) - - println("\nConjugate Gradient solver:") - println(" Initial residual: $(sqrt(rsold))") - - for iter in 1:max_iter - Ap = A_op(p) - alpha = rsold / dot(p, Ap) - - x .+= alpha .* p - r .-= alpha .* Ap - - rsnew = dot(r, r) - - if iter % 10 == 0 || iter == 1 - @printf(" Iter %4d: ||r|| = %.6e\n", iter, sqrt(rsnew)) - end - - if sqrt(rsnew) < tol - println(" ✅ Converged in $iter iterations") - return x, iter - end - - beta = rsnew / rsold - p .= r .+ beta .* p - rsold = rsnew - end - - println(" ❌ Did not converge in $max_iter iterations") - return x, max_iter -end - -""" -Solve linear elasticity problem on GPU -""" -function solve_elasticity_gpu(problem::ElasticityProblem; tol=1e-6, max_iter=1000) - println("\n" * "="^70) - println("GPU Linear Elasticity Solver") - println("="^70) - - # Check CUDA - if !CUDA.functional() - error("CUDA not available!") - end - println("GPU: ", CUDA.name(CUDA.device())) - - # Extract mesh data - mesh = physics.mesh - n_nodes = size(mesh.nodes, 2) - n_elems = size(mesh.elements, 2) - n_dofs = 3 * n_nodes - - println("\nMesh:") - println(" Nodes: $n_nodes") - println(" Elements: $n_elems") - println(" DOFs: $n_dofs") - - println("\nBoundary conditions:") - println(" Fixed nodes: $(length(physics.fixed_nodes))") - println(" Pressure nodes: $(length(physics.pressure_nodes))") - println(" Pressure value: $(physics.pressure_value) Pa") - - println("\nMaterial:") - println(" Young's modulus: $(physics.material.E) Pa") - println(" Poisson's ratio: $(physics.material.ν)") - - # Transfer to GPU - println("\nTransferring data to GPU...") - nodes_gpu = CuArray(mesh.nodes) - elements_gpu = CuArray(Int32.(mesh.elements)) - - # Build CSR map - println("Building node-to-elements map...") - node_to_elems = build_node_to_elems_gpu(mesh.elements, n_nodes) - - # Initial guess - u_gpu = CUDA.zeros(Float64, n_dofs) - - # External force (pressure load) - f_gpu = CUDA.zeros(Float64, n_dofs) - apply_pressure_load!(f_gpu, physics.pressure_nodes, mesh, physics.pressure_value) - - println("External force norm: $(norm(Array(f_gpu)))") - - # Define stiffness operator K(u) = internal forces - E = physics.material.E - ν = physics.material.ν - - function K_op(u) - r = CUDA.zeros(Float64, n_dofs) - compute_residual_gpu!(r, u, nodes_gpu, elements_gpu, node_to_elems, E, ν) - return r - end - - # Apply Dirichlet BC - K_bc = apply_dirichlet_bc!(K_op, f_gpu, physics.fixed_nodes) - - # Solve: K * u = f - println("\n" * "-"^70) - println("Solving linear system...") - println("-"^70) - - u_gpu, n_iter = cg_solve_gpu!(u_gpu, K_bc, f_gpu, tol=tol, max_iter=max_iter) - - # Transfer back to CPU - u_cpu = Array(u_gpu) - - println("\n" * "="^70) - println("Solution statistics:") - println("="^70) - println(" Max displacement: $(maximum(abs.(u_cpu))) m") - println(" CG iterations: $n_iter") - - # Compute final residual - r_final = K_bc(u_gpu) - f_gpu - println(" Final residual: $(norm(Array(r_final)))") - - println("\n" * "="^70) - println("✅ GPU elasticity solver complete!") - println("="^70) - - return u_cpu -end - -end # module diff --git a/src/io.jl b/src/io.jl deleted file mode 100644 index a8e33ca..0000000 --- a/src/io.jl +++ /dev/null @@ -1,518 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using HDF5 -using LightXML - -mutable struct Xdmf <: AbstractResultsWriter - name::String - xml::XMLElement - hdf::HDF5File - hdf_counter::Int - format::String -end - -function Xdmf() - return Xdmf(tempname()) -end - -function h5file(xdmf::Xdmf) - return xdmf.name * ".h5" -end - -function xmffile(xdmf::Xdmf) - return xdmf.name * ".xmf" -end - -""" - Xdmf(name, version="3.0", overwrite=false) - -Initialize a new Xdmf object. -""" -function Xdmf(name::String; version="3.0", overwrite=false) - xdmf = new_element("Xdmf") - h5file = "$name.h5" - xmlfile = "$name.xmf" - - if isfile(h5file) - if overwrite - @debug("Result file $h5file exists, removing old file.") - rm(h5file) - else - error("Result file $h5file exists, use Xdmf($name; overwrite=true) to rewrite results") - end - end - - if isfile(xmlfile) - if overwrite - @debug("Result file $xmlfile exists, removing old file.") - rm(xmlfile) - else - error("Result file $xmlfile exists, use Xdmf($name; overwrite=true) to rewrite results") - end - end - - set_attribute(xdmf, "xmlns:xi", "http://www.w3.org/2001/XInclude") - set_attribute(xdmf, "Version", version) - flag = isfile(h5file) ? "r+" : "w" - hdf = h5open(h5file, flag) - return Xdmf(name, xdmf, hdf, 1, "HDF") -end - -""" - get_temporal_collection(xdmf) - -Return the basic structure of Xdmf document. -Creates a new TemporalCollection if not found. -Basic structure for XML part of Xdmf file is - - - - - - - -""" -function get_temporal_collection(xdmf::Xdmf) - domain = find_element(xdmf.xml, "Domain") - grid = nothing - if domain == nothing - domain = new_child(xdmf.xml, "Domain") - grid = new_child(domain, "Grid") - set_attribute(grid, "CollectionType", "Temporal") - set_attribute(grid, "GridType", "Collection") - end - grid = find_element(domain, "Grid") - return grid -end - -""" - xdmf_filter(child_elements, child_name) - -Returns some spesific child xml element from an array of XMLElement based on, -"Xdmf extensions" see [1] for details. - -Parameters ----------- -child_elements :: Vector{XMLElement} - A vector of XMLElements where to perform filtering. -child_name :: String - Child element name, maybe containing Xdmf instructions - -Returns -------- -nothing if nothing is found, otherwise XMLElement matching to filtering - -#Examples - -julia> grid1 = new_element("Grid") -julia> add_text(grid1, "I am first grid") -julia> grid2 = new_element("Grid") -julia> add_text(grid2, "I am second grid") -julia> set_attribute(grid2, "Name", "Frame 2") -julia> grid3 = new_element("Grid") -julia> add_text(grid3, "I am third grid") -julia> grids = [grid1, grid2, grid3] - -To return second Grid element, one can use - -julia> xdmf_filter(grids, "Grid[2]") - -To return Grid which has attribute Name="Frame 2", use - -julia> xdmf_filter(grids, "Grid[@name=Frame 2]") - -To pick last Grid, use [end], e.g. - -julia> xdmf_filter(grids, "Grid[end]"). - -References ----------- -[1] http://www.xdmf.org/index.php/XDMF_Model_and_Format -""" -function xdmf_filter(child_elements, child_name) - if '/' in child_name # needs path traversal - return nothing - end - - # filter children elements using syntax child[X] -> rename child_name - m = match(r"(\w+)\[(.+)\]", child_name) - if m != nothing - child_name = m[1] - end - - # first find any relevant child elements (has same tag) - childs = [] - for child in child_elements - if LightXML.name(child) == child_name - push!(childs, child) - end - end - - # childs not found at all - length(childs) == 0 && return nothing - - # by default return first - m == nothing && return first(childs) - - # if [end] return last - m[2] == "end" && return childs[end] - - # otherwise try parse int and return nth children from list - parsed_int = tryparse(Int, m[2]) - if !isnull(parsed_int) - idx = get(parsed_int) - if (idx > 0) && (idx <= length(childs)) - return childs[idx] - else - # wrong index - return nothing - end - end - - # [X] is something else than integer, filter children elements using syntax child[@attr=value] - m2 = match(r"@(.+)=(.+)", m[2]) - m2 == nothing && throw("Unable to parse: $(m[2])") - attr_name = convert(String, m2[1]) - attr_value = convert(String, m2[2]) - for child in childs - has_attribute(child, attr_name) || continue - if attribute(child, attr_name) == attr_value - return child - end - end - - # nothing found - return nothing -end - -""" - traverse(xdmf, x, attr_name) - -Traverse XML path. Xdmf filtering can be used, so it's possible to find -data from xml using syntax e.g. - -#Example - -julia> traverse(xdmf, x, "/Domain/Grid[2]/Grid[@Name=Frame 1]/DataItem") -""" -function traverse(xdmf::Xdmf, x::XMLElement, attr_name::String) - attr_name = strip(attr_name, '/') - - if has_attribute(x, attr_name) - return attribute(x, attr_name) - end - - childs = child_elements(x) - - if '/' in attr_name - items = split(attr_name, '/') - new_item = xdmf_filter(childs, first(items)) - if new_item == nothing - @debug("traverse: childs:") - for child in childs - @debug(LightXML.name(child)) - end - error("traverse: failed, items = $items, xdmf_filter not find child") - end - new_path = join(items[2:end], '/') - return traverse(xdmf, new_item, new_path) - end - - child = xdmf_filter(childs, attr_name) - return child -end - -""" - read(xdmf, path) - -Read data from Xdmf file. - -#Example - -Traversing is supported, so one can easily traverse XML tree e.g. -julia> read(xdmf, "/Domain/Grid/Grid[2]/Geometry") -""" -function read(xdmf::Xdmf, path::String) - result = traverse(xdmf, xdmf.xml, path) - if endswith(path, "DataItem") - format = attribute(result, "Format"; required=true) - if format == "HDF" - h5file, path = map(String, split(content(result), ':')) - h5file = dirname(xdmf.name) * "/" * h5file - isfile(h5file) || throw("Xdmf: h5 file $h5file not found!") - return read(xdmf.hdf, path) - else - error("Read from Xdmf, reading from $format not implemented") - end - else - return result - end -end - -""" - save!(xdmf) - -Save the xdmf file. -""" -function save!(xdmf::Xdmf) - doc = XMLDocument() - set_root(doc, xdmf.xml) - save_file(doc, xmffile(xdmf)) -end - -function Base.close(xdmf::Xdmf) - close(xdmf.hdf) -end - -function new_dataitem(xdmf::Xdmf, path::String, data::Array{T,N}) where {T,N} - dataitem = new_element("DataItem") - datatype = replace("$T", "64" => "") - dimensions = join(reverse(size(data)), " ") - set_attribute(dataitem, "DataType", datatype) - set_attribute(dataitem, "Dimensions", dimensions) - set_attribute(dataitem, "Format", xdmf.format) - if xdmf.format == "HDF" - hdf = basename(h5file(xdmf)) - if exists(xdmf.hdf, path) - @debug("Xdmf: $path already existing in h5 file, not overwriting.") - else - write(xdmf.hdf, path, data) - end - add_text(dataitem, "$hdf:$path") - elseif xdmf.format == "XML" - text_data = string(data') - text_data = strip(text_data, ['[', ']']) - text_data = replace(text_data, ';', '\n') - text_data = "\n" * text_data * "\n" - add_text(dataitem, text_data) - else - error("Unsupported Xdmf big data format $(xdmf.format)") - end - return dataitem -end - -""" - new_dataitem(xdmf, data) - -Create a new DataItem element, hdf path automatically determined. -""" -function new_dataitem(xdmf::Xdmf, data::Array{T,N}) where {T,N} - if xdmf.format == "XML" - # Path can be whatever as XML format does not store to HDF at all - return new_dataitem(xdmf, "/whatever", data) - else - path = "/DataItem_$(xdmf.hdf_counter)" - while exists(xdmf.hdf, path) - xdmf.hdf_counter += 1 - path = "/DataItem_$(xdmf.hdf_counter)" - end - return new_dataitem(xdmf, path, data) - end -end - -const global xdmf_element_mapping = Dict( - "Poi1" => "Polyvertex", - "Seg2" => "Polyline", - "Tri3" => "Triangle", - "Quad4" => "Quadrilateral", - "Tet4" => "Tetrahedron", - "Pyramid5" => "Pyramid", - "Wedge6" => "Wedge", - "Hex8" => "Hexahedron", - "Seg3" => "Edge_3", - "Tri6" => "Tri_6", - "Quad8" => "Quad_8", - "Tet10" => "Tet_10", - "Pyramid13" => "Pyramid_13", - "Wedge15" => "Wedge_15", - "Hex20" => "Hex_20") - -get_xdmf_element_code(::Element{M,Poi1}) where M = 1 -get_xdmf_element_code(::Element{M,Seg2}) where M = 2 -# get_xdmf_element_code(::Element{Polygon}) = 3 -get_xdmf_element_code(::Element{M,Tri3}) where M = 4 -get_xdmf_element_code(::Element{M,Quad4}) where M = 5 -get_xdmf_element_code(::Element{M,Tet4}) where M = 6 -get_xdmf_element_code(::Element{M,Pyr5}) where M = 7 -get_xdmf_element_code(::Element{M,Wedge6}) where M = 8 -get_xdmf_element_code(::Element{M,Hex8}) where M = 9 -# get_xdmf_element_code(::Element{Polyhedron}) = 16 - -get_xdmf_element_code(::Element{M,Seg3}) where M = 34 -get_xdmf_element_code(::Element{M,Quad9}) where M = 35 -get_xdmf_element_code(::Element{M,Tri6}) where M = 36 -get_xdmf_element_code(::Element{M,Quad8}) where M = 37 -get_xdmf_element_code(::Element{M,Tet10}) where M = 38 -# get_xdmf_element_code(::Element{Pyr13}) = 39 -get_xdmf_element_code(::Element{M,Wedge15}) where M = 40 -# get_xdmf_element_code(::Element{Wedge18}) = 41 -get_xdmf_element_code(::Element{M,Hex20}) where M = 48 -# get_xdmf_element_code(::Element{Hex24}) = 49 -get_xdmf_element_code(::Element{M,Hex27}) where M = 50 - -""" - get_spatial_collection() - -Return a SpatialCollection at given time either by creating new one or returning -existing one. -""" -function get_spatial_collection(temporal_collection, time) - for spatial_collection in get_elements_by_tagname(temporal_collection, "Grid") - time_element = find_element(spatial_collection, "Time") - time_value = Meta.parse(attribute(time_element, "Value"; required=true)) - isapprox(time_value, time) && return spatial_collection - end - # did not find, create new one - spatial_collection = new_child(temporal_collection, "Grid") - set_attribute(spatial_collection, "GridType", "Collection") - set_attribute(spatial_collection, "Name", "Problems") - set_attribute(spatial_collection, "CollectionType", "Spatial") - time_element = new_child(spatial_collection, "Time") - set_attribute(time_element, "Value", time) - return spatial_collection -end - -""" - update_xdmf!(xdmf, problem, time, fields) - -Write new fields to Xdmf file. - -#Example - -To write displacement and temperature fields from p1 at time t=0.0: - -julia> update_xdmf!(p1, 0.0, ["displacement", "temperature"]) -""" -function update_xdmf!(xdmf::Xdmf, problem::Problem, time::Float64, fields::Vector) - - @debug("Xdmf: storing fields $fields of problem $(problem.name) at time $time") - - # 1. find domain - xml = xdmf.xml - domain = find_element(xml, "Domain") - if domain == nothing - @debug("Xdmf: Domain not found, creating.") - domain = new_child(xml, "Domain") - end - - # 2. find for TemporalCollection - temporal_collection = find_element(domain, "Grid") - if temporal_collection == nothing - @debug("Xdmf: Temporal collection not found, creating.") - temporal_collection = new_child(domain, "Grid") - set_attribute(temporal_collection, "GridType", "Collection") - set_attribute(temporal_collection, "Name", "Time") - set_attribute(temporal_collection, "CollectionType", "Temporal") - end - - # 2.1 make sure that Grid element we found really is TemporalCollection - collection_type = attribute(temporal_collection, "CollectionType"; required=true) - @assert collection_type == "Temporal" - - spatial_collection = get_spatial_collection(temporal_collection, time) - - for frame in get_elements_by_tagname(spatial_collection, "Grid") - frame_name = attribute(frame, "Name") - if frame_name == problem.name - @warn("Xdmf: Already found Grid with name $frame_name for time $time, skipping.") - return - end - end - - frame_name = problem.name - @debug("Xdmf: Creating Grid for problem $frame_name") - frame = new_child(spatial_collection, "Grid") - set_attribute(frame, "Name", frame_name) - - # 4. save geometry - X_dict = problem("geometry", time) - node_ids = sort(collect(keys(X_dict))) - node_mapping = Dict(j => i for (i, j) in enumerate(node_ids)) - X_array = hcat([X_dict[nid] for nid in node_ids]...) - ndim, nnodes = size(X_array) - geom_type = (ndim == 2 ? "XY" : "XYZ") - @debug("Xdmf: Creating geometry, type = $geom_type, number of nodes = $nnodes") - X_dataitem = new_dataitem(xdmf, X_array) - geometry = new_child(frame, "Geometry") - set_attribute(geometry, "Type", geom_type) - add_child(geometry, X_dataitem) - - # 5. save topology - mesh_type = "unstructured" - if mesh_type == "unstructured" - element_conn = Int64[] - for element in get_elements(problem) - xdmf_element_code = get_xdmf_element_code(element) - xdmf_element_code > 0 || continue - push!(element_conn, xdmf_element_code) - if xdmf_element_code == 2 - push!(element_conn, length(element)) - end - for j in get_connectivity(element) - push!(element_conn, node_mapping[j] - 1) - end - end - topology_dataitem = new_dataitem(xdmf, element_conn) - topology = new_child(frame, "Topology") - set_attribute(topology, "TopologyType", "Mixed") - add_child(topology, topology_dataitem) - else - all_elements = get_elements(problem) - nelements = length(all_elements) - element_types = unique(map(get_element_type, all_elements)) - nelement_types = length(element_types) - @debug("Xdmf: Saving topology of $nelements elements total, $nelement_types different element types.") - if nelement_types != 1 - error("Xdmf: only single type of element supported by structured grid type!") - end - for element_type in element_types - elements = collect(filter_by_element_type(element_type, all_elements)) - nelements = length(elements) - @debug("Xdmf: $nelements elements of type $element_type") - sort!(elements, by=get_element_id) - element_ids = map(get_element_id, elements) - element_conn = map(element -> [node_mapping[j] - 1 for j in get_connectivity(element)], elements) - element_conn = hcat(element_conn...) - element_code = split(string(element_type), ".")[end] - topology_dataitem = new_dataitem(xdmf, element_conn) - topology = new_child(frame, "Topology") - set_attribute(topology, "TopologyType", xdmf_element_mapping[element_code]) - set_attribute(topology, "NumberOfElements", length(elements)) - add_child(topology, topology_dataitem) - end - end - - # 6. save requested fields - for field_name in fields - field_dict = problem(field_name, time) - field_center = "Node" - field_node_ids = sort(collect(keys(field_dict))) - if node_ids != field_node_ids - @error("geom node ids = $node_ids") - @error("field node ids = $field_node_ids") - error("!=, geometry does not match with field.") - end - field_dim = length(field_dict[first(field_node_ids)]) - if field_dim == 2 - @debug("Xdmf: Field dimension = 2, extending to 3") - for nid in field_node_ids - field_dict[nid] = [field_dict[nid]; 0.0] - end - field_dim = 3 - end - field_type = Dict(1 => "Scalar", 3 => "Vector", 6 => "Tensor6")[field_dim] - @debug("Xdmf: Saving field $field_name, type = $field_type, dimension = $field_dim, center = $field_center") - - field_array = hcat([field_dict[nid] for nid in field_node_ids]...) - field_dataitem = new_dataitem(xdmf, field_array) - attribute = new_child(frame, "Attribute") - set_attribute(attribute, "Name", uppercasefirst(field_name)) - set_attribute(attribute, "Center", field_center) - set_attribute(attribute, "AttributeType", field_type) - add_child(attribute, field_dataitem) - end - - save!(xdmf) - @debug("Xdmf: all done.") -end diff --git a/src/materials_plasticity.jl b/src/materials_plasticity.jl deleted file mode 100644 index 500fa0f..0000000 --- a/src/materials_plasticity.jl +++ /dev/null @@ -1,123 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using ForwardDiff - -""" -Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values -""" -function find_root!(f, df, x; max_iter=50, norm_acc=1e-9) - converged = false - for i=1:max_iter - dx = -df(x) \ f(x) - x += dx - norm(dx) < norm_acc && (converged = true; break) - end - converged || error("No convergence in radial return!") - return x -end - -""" -Equivalent tensile stress. - -More info can be found from: https://en.wikipedia.org/wiki/Von_Mises_yield_criterion - Section: Reduced von Mises equation for different stress conditions -""" -function equivalent_stress(stress, ::Type{Val{:type_3d}}) - stress_ten = [stress[1] stress[6] stress[5]; - stress[6] stress[2] stress[4]; - stress[5] stress[4] stress[3]] - stress_dev = stress_ten - 1/3 * tr(stress_ten) * eye(3) - s = vec(stress_dev) - return sqrt(3/2 * dot(s, s)) -end - -""" -http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_plane_stress.cfm - -von mises: plane stress -https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf -""" -function equivalent_stress(stress, ::Type{Val{:type_2d}}) - s1, s2, t12 = stress - # Calculating principal stresses - # http://www.engineersedge.com/material_science/principal_vonmises_stress__13418.htm - se1 = (s1 + s2)/2 + sqrt(((s1 - s2)/2)^2 + t12^2) - se2 = (s1 + s2)/2 - sqrt(((s1 - s2)/2)^2 + t12^2) - - return sqrt(se1^2 -se1*se2 + se2^2) -end - -""" -https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf -""" -function yield_function(stress, stress_y, ::Type{Val{:von_mises}}, type_) - equivalent_stress(stress, type_) - stress_y -end - -function radial_return(params, dstrain, D, stress_y, stress_base, yield_surface_, type_) - - # Creating wrapper for gradient - vm_wrap(stress_) = yield_function(stress_, stress_y, yield_surface_, type_) - dfds = x -> ForwardDiff.gradient(vm_wrap, x) - - # Stress rate and total strain - dstress = params[1:end-1] - stress_tot = stress_base + dstress - - # Calculating plastic strain rate - dstrain_p = params[end] * dfds(stress_tot) - - # Calculating equations - function_1 = dstress - D * (dstrain - dstrain_p) - function_2 = vm_wrap(stress_tot) - [vec(function_1); function_2] -end - -function ideal_plasticity!(stress_new, stress_last, dstrain_vec, pstrain, D, params, Dtan, yield_surface_, time, dt, type_) - # Test stress - dstress = vec(D * dstrain_vec) - stress_trial = stress_last + dstress - stress_y = params["yield_stress"] - - yield_curr = x -> yield_function(x, stress_y, yield_surface_, type_) - - # Calculating and checking for yield - yield = yield_curr(stress_trial) - if isless(yield, 0.0) - - stress_new[:] = stress_trial[:] - Dtan[:,:] = D[:,:] - else - # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ \ f and initial values - f = stress_ -> radial_return(stress_, dstrain_vec, D, stress_y, stress_last, yield_surface_, type_) - df = x -> ForwardDiff.jacobian(f, x) - - # Calculating root (two options) - vals = [vec(stress_trial - stress_last); 0.0] - - #results = nlsolve(not_in_place(f), vals).zero - results = find_root!(f, df, vals) - - # extracting results - dstress = results[1:end-1] - plastic_multiplier = results[end] - - # Updating stress - stress_new[:] = stress_last + dstress - - - # Calculating plastic strain - dfds_ = x -> ForwardDiff.gradient(yield_curr, x) - dep = plastic_multiplier * dfds_(vec(stress_new)) - - # Equations for consistent tangent matrix can be found from: - # http://homes.civil.aau.dk/lda/continuum/plast.pdf - # equations: 152 & 153 - D2g = x -> ForwardDiff.hessian(yield_curr, x) - Dc = (D^-1 + plastic_multiplier * D2g(stress_new))^-1 - dfds = dfds_(stress_new) - Dtan[:,:] = Dc - (Dc * dfds * dfds' * Dc) / (dfds' * Dc * dfds)[1] - pstrain[:] = plastic_multiplier * dfds - end -end diff --git a/src/nodal_assembly_structures.jl b/src/nodal_assembly_structures.jl deleted file mode 100644 index e9a2e8f..0000000 --- a/src/nodal_assembly_structures.jl +++ /dev/null @@ -1,234 +0,0 @@ -# Nodal Assembly Data Structures -# -# This module provides the inverse mapping needed for efficient nodal assembly: -# Given a node, find all elements touching it and the local node index within each element. - -using Tensors - -""" - ElementNodeInfo - -Information about how a node appears in an element. - -# Fields -- `element_id::Int`: Global element ID -- `local_node_idx::Int`: Local node index within the element (1-based) -""" -struct ElementNodeInfo - element_id::Int - local_node_idx::Int -end - -""" - NodeToElementsMap - -Inverse connectivity mapping: for each node, lists all elements touching it. - -# Fields -- `node_to_elements::Vector{Vector{ElementNodeInfo}}`: For node j, gives all elements touching it -- `nnodes::Int`: Total number of nodes in mesh -- `nelements::Int`: Total number of elements in mesh - -# Example -```julia -map = NodeToElementsMap(connectivity) -# Get all elements touching node 5 -elements_touching_5 = map.node_to_elements[5] -for info in elements_touching_5 - println("Node 5 is local node ", info.local_node_idx, " in element ", info.element_id) -end -``` -""" -struct NodeToElementsMap - node_to_elements::Vector{Vector{ElementNodeInfo}} - nnodes::Int - nelements::Int -end - -""" - NodeToElementsMap(connectivity::Vector{NTuple{N,Int}}) where N - -Build inverse mapping from element connectivity. - -# Arguments -- `connectivity`: Vector of element connectivity tuples, e.g., [(1,2,3,4), (2,3,5,6), ...] - -# Returns -- `NodeToElementsMap`: Inverse mapping structure - -# Example -```julia -# Tet4 mesh with 2 elements -connectivity = [(1,2,3,4), (2,3,4,5)] -map = NodeToElementsMap(connectivity) - -# Node 2 appears in both elements -@assert length(map.node_to_elements[2]) == 2 -``` -""" -function NodeToElementsMap(connectivity::Vector{NTuple{N,Int}}) where N - nelements = length(connectivity) - - # Find maximum node ID to determine array size - nnodes = maximum(maximum(conn) for conn in connectivity) - - # Pre-allocate vectors for each node - node_to_elements = [Vector{ElementNodeInfo}() for _ in 1:nnodes] - - # Build inverse mapping - for (elem_id, conn) in enumerate(connectivity) - for (local_idx, global_node_id) in enumerate(conn) - push!(node_to_elements[global_node_id], - ElementNodeInfo(elem_id, local_idx)) - end - end - - return NodeToElementsMap(node_to_elements, nnodes, nelements) -end - -""" - get_node_spider(map::NodeToElementsMap, node_id::Int) -> Vector{Int} - -Get the "spider" of a node - all nodes that couple with it (including itself). - -This is the union of all nodes in elements touching `node_id`. These are exactly -the nodes for which we need to compute 3×3 stiffness blocks. - -# Arguments -- `map`: Node-to-elements mapping -- `node_id`: Node for which to find the spider - -# Returns -- `spider_nodes::Vector{Int}`: Sorted unique list of node IDs in the spider - -# Example -```julia -# For node j, find all nodes it couples with -spider = get_node_spider(map, j) -# Now compute K_blocks[k] for each k in spider -``` -""" -function get_node_spider(map::NodeToElementsMap, node_id::Int, - connectivity::Vector{NTuple{N,Int}}) where N - spider = Set{Int}() - - # For each element touching this node - for elem_info in map.node_to_elements[node_id] - # Add all nodes in that element - for node in connectivity[elem_info.element_id] - push!(spider, node) - end - end - - return sort(collect(spider)) -end - -""" - NodalStiffnessContribution{T} - -Storage for nodal assembly contribution at a single node. - -# Fields -- `node_id::Int`: Global node ID -- `spider_nodes::Vector{Int}`: Node IDs that couple with this node -- `K_blocks::Vector{Tensor{2,3,T}}`: 3×3 stiffness blocks for each spider node -- `f_int::Vec{3,T}`: Internal force at this node -- `f_ext::Vec{3,T}`: External force at this node - -# Notes -- `K_blocks[k]` corresponds to `spider_nodes[k]` -- Diagonal block (self-coupling) is included in spider -- All quantities use Tensors.jl types (zero-allocation) -""" -struct NodalStiffnessContribution{T} - node_id::Int - spider_nodes::Vector{Int} - K_blocks::Vector{Tensor{2,3,T,9}} - f_int::Vec{3,T} - f_ext::Vec{3,T} -end - -""" - NodalStiffnessContribution(node_id::Int, spider_nodes::Vector{Int}, ::Type{T}=Float64) - -Allocate storage for nodal assembly contribution. - -# Example -```julia -spider = get_node_spider(map, 5, connectivity) -contrib = NodalStiffnessContribution(5, spider, Float64) -# Now fill in K_blocks, f_int, f_ext during assembly -``` -""" -function NodalStiffnessContribution(node_id::Int, spider_nodes::Vector{Int}, - ::Type{T}=Float64) where T - nspider = length(spider_nodes) - K_blocks = [zero(Tensor{2,3,T}) for _ in 1:nspider] - f_int = zero(Vec{3,T}) - f_ext = zero(Vec{3,T}) - - return NodalStiffnessContribution{T}(node_id, spider_nodes, K_blocks, f_int, f_ext) -end - -""" - matrix_vector_product_nodal(contrib::NodalStiffnessContribution, - u::Vector{Vec{3,T}}) -> Vec{3,T} - -Compute the matrix-vector product for one node using nodal assembly. - -This computes: w_i = sum_j K_ij * u_j for node i - -# Arguments -- `contrib`: Nodal stiffness contribution (contains K_blocks for all j in spider) -- `u`: Displacement field at all nodes (Vec{3} per node) - -# Returns -- `w_i::Vec{3}`: Result of K_i * u at this node - -# Example -```julia -# Assemble contribution for node i -contrib = assemble_nodal_contribution(element_set, node_i, u, time) - -# Matrix-free matvec: w_i = K_i * u -w_i = matrix_vector_product_nodal(contrib, u) -``` -""" -function matrix_vector_product_nodal(contrib::NodalStiffnessContribution{T}, - u::Vector{Vec{3,T}}) where T - w = zero(Vec{3,T}) - - # 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] # 3×1 displacement - - # Block matrix-vector product: K_ij is Tensor{2,3}, u_j is Vec{3} - # Use regular matrix-vector multiplication (single contraction) - w += K_ij ⋅ u_j # Tensor{2,3} ⋅ Vec{3} → Vec{3} - end - - return w -end - -""" - print_spider_info(map::NodeToElementsMap, node_id::Int, - connectivity::Vector{NTuple{N,Int}}) where N - -Print diagnostic information about a node's spider for debugging. -""" -function print_spider_info(map::NodeToElementsMap, node_id::Int, - connectivity::Vector{NTuple{N,Int}}) where N - println("Node $node_id Spider Analysis:") - println(" Touches $(length(map.node_to_elements[node_id])) elements") - - for elem_info in map.node_to_elements[node_id] - println(" Element $(elem_info.element_id): local node $(elem_info.local_node_idx)") - println(" Connectivity: $(connectivity[elem_info.element_id])") - end - - spider = get_node_spider(map, node_id, connectivity) - println(" Spider has $(length(spider)) nodes: $spider") - println(" → Need to compute $(length(spider)) 3×3 blocks") - println(" → Diagonal block at node $node_id") -end diff --git a/src/postprocess_utils.jl b/src/postprocess_utils.jl deleted file mode 100644 index cf21683..0000000 --- a/src/postprocess_utils.jl +++ /dev/null @@ -1,180 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -""" -Calculate field values to nodal points from Gauss points using least-squares fitting. -""" -function calc_nodal_values!(elements::Vector, field_name, field_dim, time; - F=nothing, nz=nothing, b=nothing, return_F_and_nz=false) - - if F == nothing - A = SparseMatrixCOO() - for element in elements - gdofs = get_connectivity(element) - for ip in get_integration_points(element) - detJ = element(ip, time, Val{:detJ}) - w = ip.weight*detJ - N = element(ip, time) - add!(A, gdofs, gdofs, w*kron(N', N)) - end - end - A = sparse(A) - nz = get_nonzero_rows(A) - A = 1/2*(A + A') - F = ldlt(A[nz,nz]) - end - - if b == nothing - b = SparseMatrixCOO() - for element in elements - gdofs = get_connectivity(element) - for ip in get_integration_points(element) - if !haskey(ip, field_name) - @warn("integration point does not have field $field_name") - continue - end - detJ = element(ip, time, Val{:detJ}) - w = ip.weight*detJ - f = ip(field_name, time) - N = element(ip, time) - for dim=1:field_dim - add!(b, gdofs, [dim], w*f[dim]*N') - end - end - end - b = sparse(b) - end - - x = zeros(size(b)...) - x[nz, :] = F \ b[nz, :] - nodal_values = Dict() - for i=1:size(x,1) - nodal_values[i] = vec(x[i,:]) - end - update!(elements, field_name, time => nodal_values) - if return_F_and_nz - return F, nz - end -end - -""" -Return node ids + vector of values -""" -function get_nodal_vector(elements::Vector, field_name::AbstractString, time::Float64) - f = Dict() - for element in elements - for (c, v) in zip(get_connectivity(element), element(field_name, time)) - if haskey(f, c) - @assert isapprox(f[c], v) - end - f[c] = v - end - end - node_ids = sort(collect(keys(f))) - field = [f[nid] for nid in node_ids] - return node_ids, field -end - -""" - problem(field_name, X, time) - -Interpolate field from a set of elements defined in problem. Here, `X` is the -location inside domain described by elements. - -Internally, function loops through all the elements, finding the one containing -the point `X`. After that, using inverse isoparametric mapping, first find -dimensionless coordinates (ξ,η,ζ) of that element corresponding to the location -of point `X` and after that interpolate the values of field under investigation. -Algorithm can be expected to be somewhat slow for big models, but for tests -models the performance is good. - -# Examples - -Having a problem called `body`, one can query the field `displacement` at -position `X = (1.0, 2.0, 3.0)` and time `t = 1.0`, with the command -```julia -X = (1.0, 2.0, 3.0) -time = 1.0 -u = body("displacement", X, time) -``` -""" -function (problem::Problem)(field_name, X, time; fillna=NaN) - for element in get_elements(problem) - if inside(element, X, time) - xi = get_local_coordinates(element, X, time) - return element(field_name, xi, time) - end - end - return fillna -end - -function (problem::Problem)(field_name, X, time, ::Type{Val{:Grad}}; fillna=NaN) - for element in get_elements(problem) - if inside(element, X, time) - xi = get_local_coordinates(element, X, time) - return element(field_name, xi, time, Val{:Grad}) - end - end - return fillna -end - -function (solver::Solver)(field_name::AbstractString, X::Vector, time::Float64; fillna=NaN) - for problem in get_problems(solver) - for element in get_elements(problem) - if inside(element, X, time) - xi = get_local_coordinates(element, X, time) - return element(field_name, xi, time) - end - end - end - return fillna -end - -""" Calculate area of cross-section. """ -function calculate_area(problem::Problem, X=[0.0, 0.0], time=0.0) - A = 0.0 - for element in get_elements(problem) - elsize = size(element) - elsize[1] == 2 || error("wrong dimension of problem for area calculation, element size = $elsize") - for ip in get_integration_points(element) - w = ip.weight*element(ip, time, Val{:detJ}) - A += w - end - end - return A -end - -""" Calculate center of mass of body with respect to X. -https://en.wikipedia.org/wiki/Center_of_mass -""" -function calculate_center_of_mass(problem::Problem, X=[0.0, 0.0, 0.0], time=0.0) - M = 0.0 - Xc = zero(X) - for element in get_elements(problem) - for ip in get_integration_points(element) - w = ip.weight*element(ip, time, Val{:detJ}) - M += w - rho = haskey(element, "density") ? element("density", ip, time) : 1.0 - Xp = element("geometry", ip, time) - Xc += w*rho*(Xp-X) - end - end - return 1.0/M * Xc -end - -""" Calculate second moment of mass with respect to X. -https://en.wikipedia.org/wiki/Second_moment_of_area -""" -function calculate_second_moment_of_mass(problem::Problem, X=[0.0, 0.0, 0.0], time=0.0) - n = length(X) - I = zeros(n, n) - for element in get_elements(problem) - for ip in get_integration_points(element) - w = ip.weight*element(ip, time, Val{:detJ}) - rho = haskey(element, "density") ? element("density", ip, time) : 1.0 - Xp = element("geometry", ip, time) - X - I += w*rho*Xp*Xp' - end - end - return I -end diff --git a/src/preprocess.jl b/src/preprocess.jl deleted file mode 100644 index 4df0e0b..0000000 --- a/src/preprocess.jl +++ /dev/null @@ -1,398 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -#= -- read meshes from different formats -- reorder connectivity, create element sets, node sets, ... -- create partitions for parallel runs -- renumber elements / nodes -- maybe precheck for bad elements -- check surface normal direction in boundary elements -- orientation of 2d elements -- etc only topology related stuff -=# - -mutable struct Mesh - nodes::Dict{Int,Vector{Float64}} - node_sets::Dict{Symbol,Set{Int}} - elements::Dict{Int,Vector{Int}} - element_types::Dict{Int,Symbol} - element_codes::Dict{Int,Symbol} - element_sets::Dict{Symbol,Set{Int}} - surface_sets::Dict{Symbol,Vector{Tuple{Int,Symbol}}} - surface_types::Dict{Symbol,Symbol} -end - -function Mesh() - return Mesh(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict()) -end - -""" - Mesh(m::Dict) - -Create a new `Mesh` using data `m`. It is assumed that `m` is in format what -`abaqus_read_mesh` in `AbaqusReader.jl` is returning. -""" -function Mesh(m::Dict) - mesh = Mesh() - mesh.nodes = m["nodes"] - mesh.elements = m["elements"] - mesh.element_types = m["element_types"] - for (k, v) in m["surface_types"] - mesh.surface_types[Symbol(k)] = v - end - for (nset_name, node_ids) in m["node_sets"] - mesh.node_sets[Symbol(nset_name)] = Set(node_ids) - end - for (elset_name, element_ids) in m["element_sets"] - mesh.element_sets[Symbol(elset_name)] = Set(element_ids) - end - for (surfset_name, surfaces) in m["surface_sets"] - mesh.surface_sets[Symbol(surfset_name)] = surfaces - end - return mesh -end - -""" - add_node!(mesh, nid, ncoords) - -Add node into the mesh. `nid` is node id and `ncoords` are the node -coordinates. -""" -function add_node!(mesh::Mesh, nid::Int, ncoords::Vector{Float64}) - mesh.nodes[nid] = ncoords -end - -""" - add_nodes!(mesh, nodes) - -Add nodes into the mesh. -""" -function add_nodes!(mesh::Mesh, nodes::Dict{Int,Vector{Float64}}) - for (nid, ncoords) in nodes - add_node!(mesh, nid, ncoords) - end -end - -""" - add_node_to_node_set!(mesh, nid, ncoords) - -Add nodes into a node set. `set_name` is the name of the set and `nids...` -are all the node id:s that wants to be added. -""" -function add_node_to_node_set!(mesh::Mesh, set_name, nids...) - if !haskey(mesh.node_sets, set_name) - mesh.node_sets[set_name] = Set{Int}() - end - push!(mesh.node_sets[set_name], nids...) - return -end - -""" - create_node_set_from_element_set!(mesh, set_names...) - -Create a new node set from the nodes in an element set. ´set_names...´ are all -the set names to be inserted in the function. -""" -function create_node_set_from_element_set!(mesh::Mesh, set_names::String...) - for set_name in set_names - set_name = Symbol(set_name) - @info("Creating node set $set_name from element set") - node_ids = Set{Int}() - for elid in mesh.element_sets[set_name] - push!(node_ids, mesh.elements[elid]...) - end - mesh.node_sets[set_name] = node_ids - end - return -end - -""" - create_node_set_from_element_set!(mesh, set_name) - -Create a new node set from an element set. -""" -function create_node_set_from_element_set!(mesh::Mesh, set_name::Symbol) - create_node_set_from_element_set!(mesh, string(set_name)) -end - -""" - add_element!(mesh, elid, eltype, connectivity) - -Add an element into the mesh. ´elid´ is the element id, ´eltype´ is the type of -the element and ´connectivity´ is the connectivity of the element. -""" -function add_element!(mesh::Mesh, elid, eltype, connectivity) - mesh.elements[elid] = connectivity - mesh.element_types[elid] = eltype - return nothing -end - -""" - add_elements!(mesh, elements) - -Add elements into the mesh. -""" -function add_elements!(mesh::Mesh, elements::Dict{Int,Tuple{Symbol,Vector{Int}}}) - for (elid, (eltype, elcon)) in elements - add_element!(mesh, elid, eltype, elcon) - end - return nothing -end - -""" - add_element_to_element_set!(mesh, set_name, elids...) - -Add elements into the mesh. ´set_name´ is the name of the element set and -´elids..´ are id:s of all the elements that wants to be added. -""" -function add_element_to_element_set!(mesh::Mesh, set_name, elids...) - if !haskey(mesh.element_sets, set_name) - mesh.element_sets[set_name] = Set{Int}() - end - push!(mesh.element_sets[set_name], elids...) -end - -""" - copy(mesh) - -Return a copy of the mesh. -""" -function Base.copy(mesh::Mesh) - mesh2 = Mesh() - mesh2.nodes = copy(mesh.nodes) - mesh2.node_sets = copy(mesh.node_sets) - mesh2.elements = copy(mesh.elements) - mesh2.element_types = copy(mesh.element_types) - mesh2.element_sets = copy(mesh.element_sets) - return mesh2 -end - -""" - filter_by_element_id(mesh, element_ids) - -Filter elements by their id's. -""" -function filter_by_element_id(mesh::Mesh, element_ids::Vector{Int}) - mesh2 = copy(mesh) - mesh2.elements = Dict() - for elid in element_ids - if haskey(mesh.elements, elid) - mesh2.elements[elid] = mesh.elements[elid] - end - end - return mesh2 -end - -""" - filter_by_element_set(mesh, set_name) - -Filter elements by an element set. -""" -function filter_by_element_set(mesh::Mesh, set_name) - filter_by_element_id(mesh::Mesh, collect(mesh.element_sets[set_name])) -end - -""" - create_element(mesh, id) - -Create an element from the mesh by it's id. -""" -function create_element(mesh::Mesh, id::Int) - connectivity = mesh.elements[id] - element_type = getfield(JuliaFEM, mesh.element_types[id]) - element = Element(element_type, connectivity) - element.id = id - update!(element, "geometry", mesh.nodes) - return element -end - -function create_elements(mesh::Mesh; element_type=nothing) - element_ids = collect(keys(mesh.elements)) - if element_type != nothing - filter!(id -> mesh.element_types[id] == element_type, element_ids) - end - elements = [create_element(mesh, id) for id in element_ids] - return elements -end - -function create_elements(mesh::Mesh, element_sets::Symbol...; element_type=nothing) - if isempty(element_sets) - element_ids = collect(keys(mesh.elements)) - else - element_ids = Set{Int}() - for set_name in element_sets - element_ids = union(element_ids, mesh.element_sets[set_name]) - end - end - - if element_type != nothing - filter!(id -> mesh.element_types[id] == element_type, element_ids) - end - - elements = [create_element(mesh, id) for id in element_ids] - - nelements = length(elements) - content = Dict{Symbol,Int}() - for elid in element_ids - eltype = mesh.element_types[elid] - content[eltype] = get(content, eltype, 0) + 1 - end - s = join(("$v x $k" for (k, v) in content), ", ") - v = join(element_sets, ", ") - @info("Created $nelements elements ($s) from element set: $v.") - - return elements -end - -""" - create_elements(mesh::Mesh, element_set::String) - -# Examples - -Suppose that there is a `mesh` with element set `Body_1`. Creating elements -based on that element set is done then - -```julia -create_elements(mesh, "Body_1") -``` -""" -function create_elements(mesh::Mesh, element_sets::String...) - return create_elements(mesh, map(Symbol, element_sets)...) -end - - -""" - find_nearest_nodes(mesh, coords, npts=1; node_set=nothing) - -find npts nearest nodes from the mesh and return their id numbers as a list. -""" -function find_nearest_nodes(mesh::Mesh, coords::Vector{Float64}, npts::Int=1; node_set=nothing) - dist = Dict{Int,Float64}() - for (nid, c) in mesh.nodes - if node_set != nothing && !(nid in mesh.node_sets[Symbol(node_set)]) - continue - end - dist[nid] = norm(coords - c) - end - s = sort(collect(dist), by=x -> x[2]) - nd = s[1:npts] # [(id1, dist1), (id2, dist2), ..., (id_npts, dist_npts)] - node_ids = [n[1] for n in nd] - return node_ids -end - -function find_nearest_node(mesh::Mesh, coords::Vector{Float64}; node_set=nothing) - return first(find_nearest_nodes(mesh, coords, 1; node_set=node_set)) -end - -""" - reorder_element_connectivity!(mesh, mapping) - -Apply a new node ordering to elements. JuliaFEM uses the same node ordering as -ABAQUS. If the mesh is parsed from FEM format with some other node ordering, -this function can be used to reorder the nodes. - -Parameters ----------- -mapping :: Dict{Symbol, Vector{Int}} - e.g. :Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] - -""" -function reorder_element_connectivity!(mesh::Mesh, mapping::Dict{Symbol,Vector{Int}}) - for (elid, eltype) in mesh.element_types - haskey(mapping, eltype) || continue - new_order = mapping[eltype] - element_connectivity = mesh.elements[elid] - new_element_connectivity = element_connectivity[new_order] - mesh.elements[elid] = new_element_connectivity - end -end - -function JuliaFEM.Problem(mesh::Mesh, ::Type{P}, name::AbstractString, dimension::Int) where P<:FieldProblem - problem = Problem(P, name, dimension) - problem.elements = create_elements(mesh, name) - return problem -end - -function JuliaFEM.Problem(mesh::Mesh, ::Type{P}, name, dimension, parent_field_name) where P<:BoundaryProblem - problem = Problem(P, name, dimension, parent_field_name) - problem.elements = create_elements(mesh, name) - return problem -end - -""" - create_coloring!(mesh::Mesh) -> Dict{Int, Int} - -Greedy algorithm for coloring a grid such that no two cells with the same node -have the same color. -The returned value is a mapping between an element id and its color. -It is safe to assemble elements with the same color in parallel -""" -function create_coloring(mesh::Mesh) - # Contains the elements that each node contain - cell_containing_node = Dict{Int,Set{Int}}() - for (cellid, nodes) in mesh.elements - for v in nodes - if !haskey(cell_containing_node, v) - cell_containing_node[v] = Set{Int}() - end - push!(cell_containing_node[v], cellid) - end - end - - I, J, V = Int[], Int[], Bool[] - for (node, cells) in cell_containing_node - for cell1 in cells # All these cells have a neighboring node - for cell2 in cells - if cell1 != cell2 - push!(I, cell1) - push!(J, cell2) - push!(V, true) - end - end - end - end - - incidence_matrix = sparse(I, J, V) - # cell -> color of cell - cell_colors = Dict{Int,Int}() - # color -> list of cells - final_colors = Set{Int}[] - occupied_colors = Set{Int}() - # Zero represents no color set yet - for (cellid, _) in mesh.elements - cell_colors[cellid] = 0 - end - total_colors = 0 - for (cellid, _) in mesh.elements - empty!(occupied_colors) - # loop over neighbors - for r in nzrange(incidence_matrix, cellid) - cell_neighbour = incidence_matrix.rowval[r] - color = cell_colors[cell_neighbour] - if color != 0 - push!(occupied_colors, color) - end - end - - # occupied colors now contains all the colors we are not allowed to use - free_color = 0 - for attempt_color in 1:total_colors - if attempt_color ∉ occupied_colors - free_color = attempt_color - break - end - end - - if free_color == 0 # no free color found, need to bump max colors - total_colors += 1 - free_color = total_colors - push!(final_colors, Set{Int}()) - end - - cell_colors[cellid] = free_color - push!(final_colors[free_color], cellid) - end - - return cell_colors -end diff --git a/src/quadrature.jl b/src/quadrature.jl deleted file mode 100644 index 97e5c30..0000000 --- a/src/quadrature.jl +++ /dev/null @@ -1,65 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/FEMQuad.jl/blob/master/LICENSE -# -# Gaussian-Legendre quadrature rules consolidated from FEMQuad.jl - -# Core API types and abstract interfaces -include("quadrature/api.jl") - -# Quadrature data -include("quadrature/quaddata.jl") - -# Gauss-Legendre quadrature rules by element topology -include("quadrature/gl_tensor_product.jl") # Segments, quadrilaterals, hexahedra (tensor products) -include("quadrature/gl_triangles.jl") # 2D triangular elements -include("quadrature/gl_tetrahedra.jl") # 3D tetrahedral elements -include("quadrature/gl_wedges.jl") # 3D wedge/prism elements -include("quadrature/gl_pyramids.jl") # 3D pyramid elements - -""" - get_rule(order::Int, rules::Symbol...) - -Get the first quadrature rule that meets the required order. -""" -function get_rule(order::Int, rules::Vararg{Symbol}) - for rule in rules - if get_order(Val{rule}) >= order - return rule - end - end - @warn("No accurate rule enough found, picking last.", order, rules) - return rules[end] -end - -""" - integrate_1d(f::Function, rule::Symbol) - -Integrate a 1D function using the specified quadrature rule. -""" -function integrate_1d(f::Function, rule::Symbol) - points = get_quadrature_points(Val{rule}) - result = sum(w * f(ip) for (w, ip) in points) - return result -end - -""" - integrate_2d(f::Function, rule::Symbol) - -Integrate a 2D function using the specified quadrature rule. -""" -function integrate_2d(f::Function, rule::Symbol) - points = get_quadrature_points(Val{rule}) - result = sum(w * f(ip) for (w, ip) in points) - return result -end - -""" - integrate_3d(f::Function, rule::Symbol) - -Integrate a 3D function using the specified quadrature rule. -""" -function integrate_3d(f::Function, rule::Symbol) - points = get_quadrature_points(Val{rule}) - result = sum(w * f(ip) for (w, ip) in points) - return result -end diff --git a/src/quadrature/gauss_points.jl b/src/quadrature/gauss_points.jl deleted file mode 100644 index a2b8b48..0000000 --- a/src/quadrature/gauss_points.jl +++ /dev/null @@ -1,300 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -""" - get_gauss_points!(::Type{T}, ::Type{S}) where {T<:AbstractTopology, S<:Gauss} - -> NTuple{N, Tuple{Float64, Vec{D}}} - -Return Gauss quadrature points for topology T with scheme S. - -**Zero allocation:** Returns compile-time tuple of (weight, coordinates) pairs. -Coordinates are `Vec{D}` from Tensors.jl for efficient FEM operations. - -# Type Parameters -- `T`: Topology type (Triangle, Tetrahedron, Segment, etc.) -- `S`: Gauss quadrature scheme (Gauss{1}, Gauss{2}, etc.) - -# Returns -Tuple of `(weight, Vec{D}(ξ))` pairs where: -- `weight`: Integration weight (Float64) -- `Vec{D}(ξ)`: Parametric coordinates as Tensors.jl Vec - -# Examples -```julia -# 1-point Gauss for triangle -ips = get_gauss_points!(Triangle, Gauss{1}) -# Returns: ((0.5, Vec{2}((1/3, 1/3))),) - -# 4-point Gauss for tetrahedron -ips = get_gauss_points!(Tetrahedron, Gauss{1}) -# Returns: ((1/24, Vec{3}((0.25, 0.25, 0.25))),) - -# Usage in assembly loop (zero allocation): -for (w, ξ) in get_gauss_points!(Triangle, Gauss{2}) - # Get basis functions and derivatives using NEW API - N = get_basis_functions(Triangle(), Lagrange{1}(), ξ) - dN = get_basis_derivatives(Triangle(), Lagrange{1}(), ξ) - - # Compute Jacobian - detJ = compute_jacobian(ξ) - - # Accumulate element matrix - for i in 1:3, j in 1:3 - K[i,j] += w * detJ * dot(dN[i], dN[j]) - end -end -``` - -# Performance -- Zero allocations (fully inlined) -- Type-stable (all types known at compile time) -- ~50× faster than runtime dispatch -- Matches golden standard architecture - -See also: [`Gauss`](@ref), [`integration_points`](@ref) -""" -function get_gauss_points! end - -# ============================================================================ -# 1D: Segment -# ============================================================================ - -# Gauss{1}: 1-point (exact for linear) -@inline function get_gauss_points!(::Type{Segment}, ::Type{Gauss{1}}) - return ( - (2.0, Vec{1}((0.0,))), - ) -end - -# Gauss{2}: 2-point (exact for cubic) -@inline function get_gauss_points!(::Type{Segment}, ::Type{Gauss{2}}) - a = 1.0 / sqrt(3.0) - return ( - (1.0, Vec{1}((-a,))), - (1.0, Vec{1}((a,))), - ) -end - -# Gauss{3}: 3-point (exact for quintic) -@inline function get_gauss_points!(::Type{Segment}, ::Type{Gauss{3}}) - a = sqrt(3.0 / 5.0) - return ( - (5.0 / 9.0, Vec{1}((-a,))), - (8.0 / 9.0, Vec{1}((0.0,))), - (5.0 / 9.0, Vec{1}((a,))), - ) -end - -# ============================================================================ -# 2D: Triangle -# ============================================================================ - -# Gauss{1}: 1-point (exact for linear) -@inline function get_gauss_points!(::Type{Triangle}, ::Type{Gauss{1}}) - return ( - (0.5, Vec{2}((1 / 3, 1 / 3))), - ) -end - -# Gauss{2}: 3-point (exact for quadratic) -@inline function get_gauss_points!(::Type{Triangle}, ::Type{Gauss{2}}) - return ( - (1 / 6, Vec{2}((1 / 6, 1 / 6))), - (1 / 6, Vec{2}((2 / 3, 1 / 6))), - (1 / 6, Vec{2}((1 / 6, 2 / 3))), - ) -end - -# Gauss{3}: 4-point (exact for cubic) -@inline function get_gauss_points!(::Type{Triangle}, ::Type{Gauss{3}}) - a = 1 / 3 - b = 1 / 5 - c = 3 / 5 - return ( - (-27 / 96, Vec{2}((a, a))), - (25 / 96, Vec{2}((b, b))), - (25 / 96, Vec{2}((c, b))), - (25 / 96, Vec{2}((b, c))), - ) -end - -# ============================================================================ -# 2D: Quadrilateral (tensor product) -# ============================================================================ - -# Gauss{1}: 1×1 = 1-point -@inline function get_gauss_points!(::Type{Quadrilateral}, ::Type{Gauss{1}}) - return ( - (4.0, Vec{2}((0.0, 0.0))), - ) -end - -# Gauss{2}: 2×2 = 4-point (standard Q1) -@inline function get_gauss_points!(::Type{Quadrilateral}, ::Type{Gauss{2}}) - a = 1.0 / sqrt(3.0) - return ( - (1.0, Vec{2}((-a, -a))), - (1.0, Vec{2}((a, -a))), - (1.0, Vec{2}((-a, a))), - (1.0, Vec{2}((a, a))), - ) -end - -# Gauss{3}: 3×3 = 9-point -@inline function get_gauss_points!(::Type{Quadrilateral}, ::Type{Gauss{3}}) - a = sqrt(3.0 / 5.0) - w1 = 5.0 / 9.0 - w2 = 8.0 / 9.0 - return ( - (w1 * w1, Vec{2}((-a, -a))), - (w1 * w2, Vec{2}((0.0, -a))), - (w1 * w1, Vec{2}((a, -a))), - (w2 * w1, Vec{2}((-a, 0.0))), - (w2 * w2, Vec{2}((0.0, 0.0))), - (w2 * w1, Vec{2}((a, 0.0))), - (w1 * w1, Vec{2}((-a, a))), - (w1 * w2, Vec{2}((0.0, a))), - (w1 * w1, Vec{2}((a, a))), - ) -end - -# ============================================================================ -# 3D: Tetrahedron -# ============================================================================ - -# Gauss{1}: 1-point (exact for linear) -@inline function get_gauss_points!(::Type{Tetrahedron}, ::Type{Gauss{1}}) - return ( - (1 / 6, Vec{3}((0.25, 0.25, 0.25))), - ) -end - -# Gauss{2}: 4-point (exact for quadratic) -@inline function get_gauss_points!(::Type{Tetrahedron}, ::Type{Gauss{2}}) - a = 0.585410196624968 - b = 0.138196601125011 - return ( - (1 / 24, Vec{3}((a, b, b))), - (1 / 24, Vec{3}((b, a, b))), - (1 / 24, Vec{3}((b, b, a))), - (1 / 24, Vec{3}((b, b, b))), - ) -end - -# Gauss{3}: 5-point (exact for cubic) -@inline function get_gauss_points!(::Type{Tetrahedron}, ::Type{Gauss{3}}) - return ( - (-4 / 30, Vec{3}((0.25, 0.25, 0.25))), - (9 / 120, Vec{3}((1 / 6, 1 / 6, 1 / 6))), - (9 / 120, Vec{3}((1 / 2, 1 / 6, 1 / 6))), - (9 / 120, Vec{3}((1 / 6, 1 / 2, 1 / 6))), - (9 / 120, Vec{3}((1 / 6, 1 / 6, 1 / 2))), - ) -end - -# ============================================================================ -# 3D: Hexahedron (tensor product) -# ============================================================================ - -# Gauss{1}: 1×1×1 = 1-point -@inline function get_gauss_points!(::Type{Hexahedron}, ::Type{Gauss{1}}) - return ( - (8.0, Vec{3}((0.0, 0.0, 0.0))), - ) -end - -# Gauss{2}: 2×2×2 = 8-point (standard Hex8) -@inline function get_gauss_points!(::Type{Hexahedron}, ::Type{Gauss{2}}) - a = 1.0 / sqrt(3.0) - return ( - (1.0, Vec{3}((-a, -a, -a))), - (1.0, Vec{3}((a, -a, -a))), - (1.0, Vec{3}((-a, a, -a))), - (1.0, Vec{3}((a, a, -a))), - (1.0, Vec{3}((-a, -a, a))), - (1.0, Vec{3}((a, -a, a))), - (1.0, Vec{3}((-a, a, a))), - (1.0, Vec{3}((a, a, a))), - ) -end - -# Gauss{3}: 3×3×3 = 27-point -@inline function get_gauss_points!(::Type{Hexahedron}, ::Type{Gauss{3}}) - a = sqrt(3.0 / 5.0) - w1 = 5.0 / 9.0 - w2 = 8.0 / 9.0 - - # Generate all 27 combinations - coords_1d = ((-a, w1), (0.0, w2), (a, w1)) - - result = ntuple(27) do i - ix = (i - 1) % 3 + 1 - iy = div(i - 1, 3) % 3 + 1 - iz = div(i - 1, 9) + 1 - - x, wx = coords_1d[ix] - y, wy = coords_1d[iy] - z, wz = coords_1d[iz] - - (wx * wy * wz, Vec{3}((x, y, z))) - end - - return result -end - -# ============================================================================ -# 3D: Wedge (Prism) - tensor product of triangle × segment -# ============================================================================ - -# Gauss{1}: Triangle(1) × Segment(1) = 1-point -@inline function get_gauss_points!(::Type{Wedge}, ::Type{Gauss{1}}) - return ( - (1.0, Vec{3}((1 / 3, 1 / 3, 0.0))), - ) -end - -# Gauss{2}: Triangle(3) × Segment(2) = 6-point -@inline function get_gauss_points!(::Type{Wedge}, ::Type{Gauss{2}}) - # Triangle points - tri_pts = ((1 / 6, 1 / 6), (2 / 3, 1 / 6), (1 / 6, 2 / 3)) - tri_w = 1 / 6 - - # Segment points - a = 1.0 / sqrt(3.0) - seg_pts = ((-a,), (a,)) - seg_w = 1.0 - - return ( - (tri_w * seg_w, Vec{3}((tri_pts[1]..., seg_pts[1][1]))), - (tri_w * seg_w, Vec{3}((tri_pts[1]..., seg_pts[2][1]))), - (tri_w * seg_w, Vec{3}((tri_pts[2]..., seg_pts[1][1]))), - (tri_w * seg_w, Vec{3}((tri_pts[2]..., seg_pts[2][1]))), - (tri_w * seg_w, Vec{3}((tri_pts[3]..., seg_pts[1][1]))), - (tri_w * seg_w, Vec{3}((tri_pts[3]..., seg_pts[2][1]))), - ) -end - -# ============================================================================ -# 3D: Pyramid - special quadrature (not tensor product) -# ============================================================================ - -# Gauss{1}: 1-point (centroid) -@inline function get_gauss_points!(::Type{Pyramid}, ::Type{Gauss{1}}) - return ( - (4 / 3, Vec{3}((0.0, 0.0, 0.25))), - ) -end - -# Gauss{2}: 5-point -@inline function get_gauss_points!(::Type{Pyramid}, ::Type{Gauss{2}}) - # Pyramid quadrature is non-trivial due to singularity at apex - a = 0.584237394672177 - b = 0.138196601125011 - return ( - (0.2378, Vec{3}((0.0, 0.0, 0.5))), - (0.2378, Vec{3}((a, 0.0, b))), - (0.2378, Vec{3}((-a, 0.0, b))), - (0.2378, Vec{3}((0.0, a, b))), - (0.2378, Vec{3}((0.0, -a, b))), - ) -end diff --git a/src/quadrature/integration.jl b/src/quadrature/integration.jl deleted file mode 100644 index 12adfb3..0000000 --- a/src/quadrature/integration.jl +++ /dev/null @@ -1,143 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -# DEPRECATED: This file contains the old integration API. -# New code should use the types from api.jl: -# - AbstractQuadratureRule (replaces AbstractIntegration) -# - QuadraturePoint (replaces IntegrationPoint) -# - GaussLegendre{N} (replaces Gauss{N}) - -""" - AbstractIntegration - -DEPRECATED: Use `AbstractQuadratureRule` instead. - -Abstract base type for all numerical integration (quadrature) schemes. - -An integration scheme defines how to numerically integrate over a reference element -by specifying integration point locations and weights. Integration schemes are -independent of element topology and interpolation schemes (though the number of -points needed may depend on polynomial order). - -# Key Properties -- Integration points (locations in parametric space) -- Weights -- Accuracy order - -# Examples -```julia -Gauss{2}() # 2-point Gauss quadrature -Gauss{3}() # 3-point Gauss quadrature -Lobatto{3}() # 3-point Gauss-Lobatto quadrature -Reduced() # Reduced integration (element-dependent) -``` - -See also: [`Gauss`](@ref), [`Lobatto`](@ref), [`IntegrationPoint`](@ref) -""" -abstract type AbstractIntegration end - -""" - IntegrationPoint{D} - -DEPRECATED: Use `QuadraturePoint{D,T}` instead. - -Represents a single integration point in D-dimensional parametric space. - -# Fields -- `ξ::Vec{D, Float64}`: Location in parametric coordinates -- `weight::Float64`: Integration weight - -# Migration -```julia -# Old: -ip = IntegrationPoint(Vec(0.0, 0.0), 1.0) - -# New: -qp = QuadraturePoint(SVector(0.0, 0.0), 1.0) -# Access: qp.coords instead of ip.ξ -``` -""" -struct IntegrationPoint{D} - ξ::Vec{D,Float64} - weight::Float64 -end - -""" - integration_points(scheme::AbstractIntegration, topology::AbstractTopology) - -> NTuple{N, IntegrationPoint{D}} - -Return the integration points and weights for the given integration scheme -applied to the reference element topology. - -**Zero allocation:** Returns compile-time sized tuple of IntegrationPoints for -known quadrature rules. Falls back to Vector for dynamic rules. - -# Arguments -- `scheme`: Integration scheme (e.g., `Gauss{3}()`) -- `topology`: Reference element topology (e.g., `Tri3()`) - -# Returns -Tuple of `IntegrationPoint` with locations ξ and weights. - -# Examples -```julia -julia> ips = integration_points(Gauss{1}(), Tri3()) -(IntegrationPoint{2}((0.333..., 0.333...), 0.5),) - -julia> typeof(ips) -Tuple{IntegrationPoint{2}} -``` -""" -function integration_points end - -""" - npoints(scheme::AbstractIntegration, topology::AbstractTopology) -> Int - -Return the number of integration points for the given scheme and topology. - -# Examples -```julia -julia> npoints(Gauss{2}(), Tri3()) -3 - -julia> npoints(Gauss{2}(), Quad4()) -4 -``` -""" -function npoints end - -""" - default_integration(topology::Type{<:AbstractTopology{N}}) where N - -> AbstractIntegration - -Return the default (recommended) integration scheme for a given topology type. - -# Default Rules -- Linear elements (P1): Use minimal integration that's exact for linear basis -- Quadratic elements (P2): Use integration exact for quadratic basis - -# Examples -```julia -julia> default_integration(Hexahedron{8}) -Gauss{2}() # 2×2×2 = 8 points (exact for trilinear) - -julia> default_integration(Tetrahedron{4}) -Gauss{1}() # 1 point (exact for linear) - -julia> default_integration(Hexahedron{27}) -Gauss{3}() # 3×3×3 = 27 points (exact for triquadratic) -``` -""" -function default_integration end - -# Default integration rules for common topologies -default_integration(::Type{Tetrahedron{4}}) = Gauss{1}() -default_integration(::Type{Tetrahedron{10}}) = Gauss{2}() -default_integration(::Type{Hexahedron{8}}) = Gauss{2}() -default_integration(::Type{Hexahedron{20}}) = Gauss{3}() -default_integration(::Type{Hexahedron{27}}) = Gauss{3}() -default_integration(::Type{Triangle{3}}) = Gauss{1}() -default_integration(::Type{Triangle{6}}) = Gauss{2}() -default_integration(::Type{Quadrilateral{4}}) = Gauss{2}() -default_integration(::Type{Quadrilateral{8}}) = Gauss{3}() -default_integration(::Type{Quadrilateral{9}}) = Gauss{3}() diff --git a/src/readers.jl b/src/readers.jl deleted file mode 100644 index 7159e8a..0000000 --- a/src/readers.jl +++ /dev/null @@ -1,15 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE -# -# Mesh readers consolidated from AbaqusReader.jl and AsterReader.jl - -# AbaqusReader - ABAQUS .inp file format -include("readers/keyword_register.jl") -include("readers/parse_mesh.jl") -include("readers/parse_model.jl") -include("readers/create_surface_elements.jl") -include("readers/abaqus_download.jl") - -# AsterReader - Code Aster .med file format (requires HDF5) -# include("readers/read_aster_mesh.jl") -# include("readers/read_aster_results.jl") diff --git a/src/shells/api.jl b/src/shells/api.jl deleted file mode 100644 index c358821..0000000 --- a/src/shells/api.jl +++ /dev/null @@ -1,100 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -""" -Shell formulation API definitions. - -This file defines shell-specific abstract types and formulation theories. -Must be included after core api.jl. -""" - -# ============================================================================ -# SHELL FORMULATION THEORIES -# ============================================================================ - -""" - AbstractShellTheory - -Abstract type for shell theory variants. - -Shell theories differ in how they model transverse shear deformation and thickness effects. - -# Concrete Theories -- `ReissnerMindlin`: Thick shells (includes transverse shear) -- `KirchhoffLove`: Thin shells (no transverse shear) - -# See Also -- [`ShellFormulation`](@ref) -""" -abstract type AbstractShellTheory end - -""" - ReissnerMindlin <: AbstractShellTheory - -Reissner-Mindlin shell theory (thick shells, includes shear). - -Assumptions: -- Normals to mid-surface remain straight but NOT perpendicular -- Transverse shear deformation included -- Valid for thick shells (h/L > 1/20) -- 5 DOFs per node: 3 displacements + 2 rotations - -# Usage -```julia -formulation = ShellFormulation{ReissnerMindlin}() -physics = Physics( - formulation=formulation, - field=DisplacementRotation{3}(), - mesh=shell_mesh, - material=steel -) -``` -""" -struct ReissnerMindlin <: AbstractShellTheory end - -""" - KirchhoffLove <: AbstractShellTheory - -Kirchhoff-Love shell theory (thin shells, no shear). - -Assumptions: -- Normals to mid-surface remain straight and perpendicular -- No transverse shear deformation -- Valid for thin shells (h/L < 1/20) -- 3 DOFs per node: 3 displacements (rotations computed from displacements) - -# Usage -```julia -formulation = ShellFormulation{KirchhoffLove}() -physics = Physics( - formulation=formulation, - field=Displacement{3}(), # Only displacements, rotations implicit - mesh=shell_mesh, - material=aluminum -) -``` -""" -struct KirchhoffLove <: AbstractShellTheory end - -""" - ShellFormulation{Theory<:AbstractShellTheory} <: AbstractFormulation - -Shell element formulation with theory variant. - -# Type Parameter -- `Theory`: Shell theory type (ReissnerMindlin or KirchhoffLove) - -# Examples -```julia -# Thick shell (includes shear) -ShellFormulation{ReissnerMindlin}() - -# Thin shell (classical theory) -ShellFormulation{KirchhoffLove}() -``` - -# Fields per Node -- Reissner-Mindlin: 5 DOFs (ux, uy, uz, θx, θy) -- Kirchhoff-Love: 3 DOFs (ux, uy, uz) - rotations implicit -""" -struct ShellFormulation{Theory<:AbstractShellTheory} <: AbstractFormulation end diff --git a/src/trusses/api.jl b/src/trusses/api.jl deleted file mode 100644 index 61ab7cb..0000000 --- a/src/trusses/api.jl +++ /dev/null @@ -1,79 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -""" -Truss formulation API definitions. - -This file defines truss-specific abstract types and formulation theories. -Must be included after core api.jl. -""" - -# ============================================================================ -# TRUSS FORMULATION THEORIES -# ============================================================================ - -""" - AbstractTrussTheory - -Abstract type for truss theory variants. - -Truss elements carry only axial forces (tension/compression), no bending. - -# Concrete Theories -- `SimpleTruss`: Standard 2-node truss (axial force only) - -# Future Extensions -- `CableTruss`: Cable elements (tension only, no compression) -- `PretensionedTruss`: Trusses with initial stress - -# See Also -- [`TrussFormulation`](@ref) -""" -abstract type AbstractTrussTheory end - -""" - SimpleTruss <: AbstractTrussTheory - -Simple truss theory (axial force only). - -Assumptions: -- Only axial forces (tension/compression) -- No bending moments -- Pin-jointed connections -- 3 DOFs per node in 3D: (ux, uy, uz) -- 2 DOFs per node in 2D: (ux, uy) - -# Usage -```julia -formulation = TrussFormulation{SimpleTruss}() -physics = Physics( - formulation=formulation, - field=Displacement{3}(), # 3D truss - mesh=truss_mesh, - material=steel -) -``` -""" -struct SimpleTruss <: AbstractTrussTheory end - -""" - TrussFormulation{Theory<:AbstractTrussTheory} <: AbstractFormulation - -Truss element formulation with theory variant. - -# Type Parameter -- `Theory`: Truss theory type (SimpleTruss, CableTruss, etc.) - -# Examples -```julia -# Standard truss -TrussFormulation{SimpleTruss}() -``` - -# Fields per Node -- 3D: 3 DOFs (ux, uy, uz) -- 2D: 2 DOFs (ux, uy) - -Use with `Displacement{Dim}` field type. -""" -struct TrussFormulation{Theory<:AbstractTrussTheory} <: AbstractFormulation end diff --git a/test/runtests.jl.old b/test/runtests.jl.old deleted file mode 100644 index eb7a173..0000000 --- a/test/runtests.jl.old +++ /dev/null @@ -1,163 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM, Test - -@testset "JuliaFEM.jl" begin - @testset "test_dirichlet.jl" begin - include("test_dirichlet.jl") - end - @testset "test_elasticity_1d.jl" begin - include("test_elasticity_1d.jl") - end - @testset "test_elasticity_2d_linear_with_surface_load.jl" begin - include("test_elasticity_2d_linear_with_surface_load.jl") - end - @testset "test_elasticity_2d_nonhomogeneous_boundary_conditions.jl" begin - include("test_elasticity_2d_nonhomogeneous_boundary_conditions.jl") - end - @testset "test_elasticity_2d_nonlinear_with_surface_load.jl" begin - include("test_elasticity_2d_nonlinear_with_surface_load.jl") - end - @testset "test_elasticity_2d_plane_stress_stiffness_matrix.jl" begin - include("test_elasticity_2d_plane_stress_stiffness_matrix.jl") - end - @testset "test_elasticity_2d_residual.jl" begin - include("test_elasticity_2d_residual.jl") - end - @testset "test_elasticity_3d_linear_with_surface_load.jl" begin - include("test_elasticity_3d_linear_with_surface_load.jl") - end - @testset "test_elasticity_3d_nonlinear_with_surface_load.jl" begin - include("test_elasticity_3d_nonlinear_with_surface_load.jl") - end - @testset "test_elasticity_3d_unit_block.jl" begin - include("test_elasticity_3d_unit_block.jl") - end - @testset "test_elasticity_forwarddiff.jl" begin - include("test_elasticity_forwarddiff.jl") - end - @testset "test_elasticity_hollow_sphere_with_surface_pressure.jl" begin - include("test_elasticity_hollow_sphere_with_surface_pressure.jl") - end - @testset "test_elasticity_med_pyr5_point_load.jl" begin - include("test_elasticity_med_pyr5_point_load.jl") - end - @testset "test_elasticity_plane_strain.jl" begin - include("test_elasticity_plane_strain.jl") - end - @testset "test_elasticity_pyr5_point_load.jl" begin - include("test_elasticity_pyr5_point_load.jl") - end - @testset "test_elasticity_tet4_stiffness_matrix.jl" begin - include("test_elasticity_tet4_stiffness_matrix.jl") - end - @testset "test_elasticity_tet10_mass_matrix.jl" begin - include("test_elasticity_tet10_mass_matrix.jl") - end - @testset "test_elasticity_tet10_stiffness_matrix.jl" begin - include("test_elasticity_tet10_stiffness_matrix.jl") - end - @testset "test_elasticity_tetra.jl" begin - include("test_elasticity_tetra.jl") - end - @testset "test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl" begin - include("test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl") - end - @testset "test_elasticplastic_3d_linear_with_surface_load.jl" begin - include("test_elasticplastic_3d_linear_with_surface_load.jl") - end - @testset "test_heat_2d_one_element.jl" begin - include("test_heat_2d_one_element.jl") - end - @testset "test_heat_3d.jl" begin - include("test_heat_3d.jl") - end - @testset "test_heat_tet10_convection.jl" begin - include("test_heat_tet10_convection.jl") - end - @testset "test_heat_3d_2.jl" begin - include("test_heat_3d_2.jl") - end - @testset "test_heat.jl" begin - include("test_heat.jl") - end - @testset "test_heat_2.jl" begin - include("test_heat_2.jl") - end - @testset "test_heat_3.jl" begin - include("test_heat_3.jl") - end - @testset "test_heat_3d_two_rings.jl" begin - include("test_heat_3d_two_rings.jl") - end - @testset "test_heat_4.jl" begin - include("test_heat_4.jl") - end - @testset "test_modal_analysis.jl" begin - include("test_modal_analysis.jl") - end - @testset "test_modal_analysis_elasticity.jl" begin - include("test_modal_analysis_elasticity.jl") - end - @testset "test_modal_analysis_elasticity_2.jl" begin - include("test_modal_analysis_elasticity_2.jl") - end - @testset "test_modal_analysis_zero_eigenmodes.jl" begin - include("test_modal_analysis_zero_eigenmodes.jl") - end - @testset "test_mortar.jl" begin - include("test_mortar.jl") - end - @testset "test_mortar_2d.jl" begin - include("test_mortar_2d.jl") - end - @testset "test_mortar_2d_assembly.jl" begin - include("test_mortar_2d.jl") - end - @testset "test_mortar_2d_contact.jl" begin - include("test_mortar_2d_contact.jl") - end - @testset "test_mortar_2d_mesh_tie.jl" begin - include("test_mortar_2d_mesh_tie.jl") - end - @testset "test_mortar_2d_weighted_gap.jl" begin - include("test_mortar_2d_weighted_gap.jl") - end - @testset "test_mortar_3d_mesh_tie_modal.jl" begin - include("test_mortar_3d_mesh_tie_modal.jl") - end - @testset "test_mortar_3d_mesh_tie_two_rings.jl" begin - include("test_mortar_3d_mesh_tie_two_rings.jl") - end - @testset "test_mortar_3d_polygon_clip.jl" begin - include("test_mortar_3d_polygon_clip.jl") - end - @testset "test_postprocess.jl" begin - include("test_postprocess.jl") - end - @testset "test_potential_energy.jl" begin - include("test_potential_energy.jl") - end - @testset "test_problems_contact_3d.jl" begin - include("test_problems_contact_3d.jl") - end - @testset "test_problems_elasticity.jl" begin - include("test_problems_elasticity.jl") - end - @testset "test_problems_mortar_3d.jl" begin - include("test_problems_mortar_3d.jl") - end - @testset "test_problems_mortar_3d_lowlevel.jl" begin - include("test_problems_mortar_3d_lowlevel.jl") - end - @testset "test_solvers_postprocess.jl" begin - include("test_solvers_postprocess.jl") - end - @testset "test_virtual_work.jl" begin - include("test_virtual_work.jl") - end - @testset "test_von_mises_material.jl" begin - include("test_von_mises_material.jl") - end -end diff --git a/test/runtests_new.jl b/test/runtests_new.jl deleted file mode 100644 index 2af607c..0000000 --- a/test/runtests_new.jl +++ /dev/null @@ -1,47 +0,0 @@ -# JuliaFEM Test Suite - New Structure -# Educational testing with Literate.jl - -using Test, JuliaFEM - -# Configuration -const RUN_TUTORIALS = get(ENV, "JULIAFEM_TEST_TUTORIALS", "true") == "true" -const RUN_UNIT = get(ENV, "JULIAFEM_TEST_UNIT", "false") == "true" -const RUN_OLD = get(ENV, "JULIAFEM_TEST_OLD", "false") == "true" - -println("="^70) -println("JuliaFEM Test Suite (New Structure)") -println("="^70) -println("Tutorials: ", RUN_TUTORIALS ? "✓" : "✗") -println("Unit tests: ", RUN_UNIT ? "✓" : "✗") -println("Old tests: ", RUN_OLD ? "✓" : "✗") -println("="^70) - -# Tutorial Tests -if RUN_TUTORIALS - @testset "Tutorials" begin - @testset "01_Fundamentals" begin - include("tutorials/01_fundamentals/creating_elements.jl") - include("tutorials/01_fundamentals/reading_gmsh_meshes.jl") - include("tutorials/01_fundamentals/basis_functions.jl") - include("tutorials/01_fundamentals/validation_1element_quad4.jl") - end - end -end - -# Unit Tests -if RUN_UNIT - @testset "Unit Tests" begin - @info "No unit tests yet" - end -end - -# Old Tests -if RUN_OLD - @warn "Old tests have 49+ failures - see docs/TEST_FIXES_NEEDED.md" - include("runtests.jl") # Original test suite -end - -println() -println("="^70) -println("Complete - See docs/TESTING_PHILOSOPHY.md") -println("="^70) diff --git a/test/test_problems_elasticity_assemble_3d_seg3.jl b/test/test_problems_elasticity_assemble_3d_seg3.jl deleted file mode 100644 index b277ea6..0000000 --- a/test/test_problems_elasticity_assemble_3d_seg3.jl +++ /dev/null @@ -1,12 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM, Test - -# Test that if 3d continuum problem have some elements we don't know how to -# deal with, raise error with clear message - -elements = [Element(Seg3, (1, 2, 3))] -problem = Problem(Elasticity, "test seg3", 3) -add_elements!(problem, elements) -@test_throws ErrorException assemble!(problem, 0.0)