diff --git a/src/assemblers/element_based_csc.jl b/src/assemblers/element_based_csc.jl deleted file mode 100644 index 19f2ae8..0000000 --- a/src/assemblers/element_based_csc.jl +++ /dev/null @@ -1,375 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -""" -CSC (Compressed Sparse Column) assembly with pre-built structure. - -Optimized element-by-element assembly using pre-allocated CSC sparse matrix. -Uses two-pointer merge algorithm to insert element contributions directly -into CSC arrays without intermediate COO format. - -**Performance**: 4.1x faster than COO, 16.6x less memory. -**Best for**: Production code, nonlinear problems (repeated assembly). - -**Algorithm inspiration**: Ferrite.jl CSC assembly, adapted for JuliaFEM. -""" - -using SparseArrays - -""" - assemble!( - cache::CSCCache, - assembler::CSCAssembler, - kernel::AbstractKernel, - mesh::AbstractMesh - ) -> Nothing - -Assemble global system using pre-built CSC structure **in-place, zero allocations**. - -# Algorithm - -1. Reset cache (zero matrix values, keep structure) -2. Loop over elements: - a. Compute element stiffness in-place: `compute_element_stiffness!(cache.element_cache, ...)` - b. Get DOF mapping: `get_dof_mapping!(cache.element_cache.dofs, ...)` - c. Merge Ke into CSC matrix: two-pointer algorithm, O(nnz) per element - d. Scatter fe to global force vector: `f[dofs] += fe` - -# Two-Pointer Merge - -For each column j in element matrix: -1. Find range in CSC column: `K.colptr[j]:K.colptr[j+1]-1` -2. Find range in element column: rows matching DOF indices -3. Merge element values into CSC values using two pointers - -**Key insight**: CSC structure pre-built, so we know where every (i,j) entry is. -Just need to find correct position and add value. - -# Arguments -- `cache`: Pre-allocated CSC cache (with pre-built sparsity pattern) -- `assembler`: CSC assembler -- `kernel`: Domain kernel (continuum, plate, beam, etc.) -- `mesh`: Finite element mesh - -# Zero-Allocation Guarantee - -**Absolute zero allocations** during assembly. All arrays pre-allocated. -CSC structure built once in `CSCCache(mesh, kernel)`, reused forever. - -# Example - -```julia -# Setup (one-time cost: build sparsity pattern) -mesh = create_cantilever_mesh(10, 2, 2) -kernel = ContinuumKernel(formulation, material, field) -assembler = CSCAssembler() -cache = CSCCache(mesh, kernel) # Pre-builds K structure - -# Nonlinear loop (zero allocations every iteration) -for iteration in 1:max_iter - assemble!(cache, assembler, kernel, mesh) # Zero allocations! - K, f = extract_system(cache) # Just returns references - u = K \\ f - # ... update ... -end -``` - -# Performance - -For 2500 Tet4 elements: -- Time: 2.36 ms (4.1x faster than COO) -- Memory: 506 KB (16.6x less than COO) -- Speedup: 4.1x -""" -function assemble!( - cache::CSCCache, - assembler::CSCAssembler, - kernel::AbstractKernel, - mesh::AbstractMesh -) - # Reset cache for new assembly - reset!(cache) - - nelems = nelements(mesh) - element_cache = cache.element_cache - K = cache.K - - # Loop over elements - for elem_id in 1:nelems - # Compute element stiffness in-place (zero allocations) - compute_element_stiffness!(element_cache, kernel, elem_id, mesh) - - # Get element nodes and DOF count - nodes = mesh.connectivity[elem_id] - nnodes_elem = length(nodes) - ndofs_per_node = dofs_per_node(kernel) - ndofs_elem = nnodes_elem * ndofs_per_node - - # Get DOF mapping in-place (zero allocations) - dofs = @view element_cache.dofs[1:ndofs_elem] - get_dof_mapping!(dofs, kernel, elem_id, mesh) - - # Merge element stiffness into CSC matrix (zero allocations) - Ke = @view element_cache.Ke[1:ndofs_elem, 1:ndofs_elem] - merge_into_csc!(K, Ke, dofs) - - # Scatter element force to global force vector - fe = @view element_cache.fe[1:ndofs_elem] - for (i_local, i_global) in enumerate(dofs) - cache.f[i_global] += fe[i_local] - end - end - - return nothing -end - -""" - merge_into_csc!(K::SparseMatrixCSC, Ke::AbstractMatrix, dofs::AbstractVector{Int}) - -Merge element stiffness matrix into CSC sparse matrix **in-place, zero allocations**. - -Uses two-pointer algorithm to find correct positions in CSC arrays and -accumulate element contributions. - -# Algorithm - -For each column `j_local` in element matrix `Ke`: -1. Get global column index: `j_global = dofs[j_local]` -2. Find CSC column range: `range = K.colptr[j_global]:K.colptr[j_global+1]-1` -3. For each row `i_local` in element column: - a. Get global row index: `i_global = dofs[i_local]` - b. Find position in CSC column where `K.rowval[pos] == i_global` (binary search) - c. Add element contribution: `K.nzval[pos] += Ke[i_local, j_local]` - -# Optimization - -CSC columns are sorted by row index, so we can use: -- Binary search to find position: O(log n) per entry -- Linear scan with two pointers: O(n) per column (faster in practice) - -Current implementation: Linear scan (simpler, still very fast). - -# Arguments -- `K`: CSC sparse matrix (modified in-place) -- `Ke`: Element stiffness matrix [ndofs_elem × ndofs_elem] -- `dofs`: Global DOF indices [ndofs_elem] - -# Zero-Allocation Guarantee - -No allocations - modifies `K.nzval` in-place. Structure unchanged. -""" -function merge_into_csc!( - K::SparseMatrixCSC{Float64,Int}, - Ke::AbstractMatrix, - dofs::AbstractVector{Int} -) - ndofs_elem = length(dofs) - rowval = K.rowval - nzval = K.nzval - colptr = K.colptr - - # Loop over columns of element matrix - for j_local in 1:ndofs_elem - j_global = dofs[j_local] - - # Get range of CSC column j_global - col_start = colptr[j_global] - col_end = colptr[j_global + 1] - 1 - - # Loop over rows of element column - for i_local in 1:ndofs_elem - i_global = dofs[i_local] - value = Ke[i_local, j_local] - - # Find position in CSC column where rowval[pos] == i_global - pos = find_row_in_column(rowval, col_start, col_end, i_global) - - if pos != -1 - # Found: add contribution - nzval[pos] += value - else - # Not found: this should never happen if sparsity pattern is correct - error("CSC assembly error: entry ($i_global, $j_global) not in sparsity pattern. " * - "This indicates a bug in sparsity pattern construction.") - end - end - end - - return nothing -end - -""" - find_row_in_column(rowval::Vector{Int}, start::Int, stop::Int, row::Int) -> Int - -Find position where `rowval[pos] == row` in range `[start, stop]`. - -Uses linear scan (CSC columns are sorted, could use binary search, but -linear is fast enough for typical element sizes). - -# Arguments -- `rowval`: CSC row indices array -- `start`: Start of column range -- `stop`: End of column range -- `row`: Target row index - -# Returns -- Position where `rowval[pos] == row`, or `-1` if not found - -# Complexity - -O(n) where n = column length. For typical elements (4-27 nodes), this is -very fast (< 100 comparisons). -""" -function find_row_in_column( - rowval::Vector{Int}, - start::Int, - stop::Int, - row::Int -) - for pos in start:stop - if rowval[pos] == row - return pos - end - end - return -1 # Not found -end - -""" - find_row_in_column_binary(rowval::Vector{Int}, start::Int, stop::Int, row::Int) -> Int - -Find position using binary search (alternative to linear scan). - -CSC columns are sorted by row index, so binary search is O(log n). -In practice, linear scan is often faster for small n (< 100 entries per column). - -# Arguments -- `rowval`: CSC row indices array (sorted) -- `start`: Start of column range -- `stop`: End of column range -- `row`: Target row index - -# Returns -- Position where `rowval[pos] == row`, or `-1` if not found - -# Complexity - -O(log n) where n = column length. -""" -function find_row_in_column_binary( - rowval::Vector{Int}, - start::Int, - stop::Int, - row::Int -) - left = start - right = stop - - while left <= right - mid = (left + right) ÷ 2 - mid_row = rowval[mid] - - if mid_row == row - return mid - elseif mid_row < row - left = mid + 1 - else - right = mid - 1 - end - end - - return -1 # Not found -end - -# ============================================================================ -# HELPER FUNCTIONS -# ============================================================================ - -""" - create_cache(assembler::CSCAssembler, mesh::AbstractMesh, kernel::AbstractKernel) -> CSCCache - -Create pre-allocated cache for CSC assembly. - -**One-time cost**: Builds sparsity pattern by analyzing mesh connectivity. -Structure is cached and reused for all subsequent assemblies. - -# Arguments -- `assembler`: CSC assembler -- `mesh`: Finite element mesh -- `kernel`: Domain kernel - -# Returns -- Pre-allocated CSC cache with pre-built sparse matrix structure - -# Example - -```julia -cache = create_cache(CSCAssembler(), mesh, kernel) # Build structure once - -# Reuse in nonlinear loop (zero allocations) -for iter in 1:max_iter - assemble!(cache, CSCAssembler(), kernel, mesh) - K, f = extract_system(cache) - # ... -end -``` -""" -function create_cache(assembler::CSCAssembler, mesh::AbstractMesh, kernel::AbstractKernel) - return CSCCache(mesh, kernel) -end - -""" - verify_sparsity_pattern(K::SparseMatrixCSC, mesh::AbstractMesh, kernel::AbstractKernel) -> Bool - -Verify that CSC sparsity pattern contains all element connections. - -Debug function to check that `build_sparsity_pattern` is correct. - -# Arguments -- `K`: CSC sparse matrix -- `mesh`: Finite element mesh -- `kernel`: Domain kernel - -# Returns -- `true` if all element DOF pairs are in K structure - -# Throws -- `ErrorException` if missing entries found -""" -function verify_sparsity_pattern( - K::SparseMatrixCSC, - mesh::AbstractMesh, - kernel::AbstractKernel -) - nelems = nelements(mesh) - ndofs_per_node = dofs_per_node(kernel) - dof_buffer = Int[] - - for elem_id in 1:nelems - nodes = mesh.connectivity[elem_id] - nnodes_elem = length(nodes) - ndofs_elem = nnodes_elem * ndofs_per_node - - # Get DOF mapping - resize!(dof_buffer, ndofs_elem) - get_dof_mapping!(dof_buffer, kernel, elem_id, mesh) - - # Check all (i,j) pairs exist in K - for j_local in 1:ndofs_elem - j_global = dof_buffer[j_local] - col_start = K.colptr[j_global] - col_end = K.colptr[j_global + 1] - 1 - - for i_local in 1:ndofs_elem - i_global = dof_buffer[i_local] - - pos = find_row_in_column(K.rowval, col_start, col_end, i_global) - if pos == -1 - error("Sparsity pattern missing entry ($i_global, $j_global) " * - "from element $elem_id") - end - end - end - end - - return true -end