From aeadbb2fdd6a4777467b14645ef744380e4f0868 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Thu, 20 Nov 2025 16:56:39 +0200 Subject: [PATCH] feat(assemblers): Add general scatter_blocks_to_triplets! MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit New file: src/assemblers/scatter_blocks_to_triplets.jl (74 lines) Features: - scatter_blocks_to_triplets!(I, J, V, counter, capacity, K_blocks, dofs, N) - General scatter for unsymmetric matrices - Stores all N×N blocks as triplets - Returns updated counter Implementation: - Double loop over node pairs (N × N) - Triple loop over DOF pairs (3 × 3 per block) - Direct triplet storage: I[idx], J[idx], V[idx] - Capacity checking with bounds validation Not currently used (symmetric version preferred), but available for unsymmetric problems like convection or non-symmetric contact. --- src/assemblers/scatter_blocks_to_triplets.jl | 78 ++++++++++++++++++++ 1 file changed, 78 insertions(+) create mode 100644 src/assemblers/scatter_blocks_to_triplets.jl diff --git a/src/assemblers/scatter_blocks_to_triplets.jl b/src/assemblers/scatter_blocks_to_triplets.jl new file mode 100644 index 0000000..ebcc805 --- /dev/null +++ b/src/assemblers/scatter_blocks_to_triplets.jl @@ -0,0 +1,78 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +""" + scatter_blocks_to_triplets!( + cache::COOCache, + K_blocks::AbstractMatrix{<:Tensor{2,3}}, + dofs::AbstractVector{Int}, + N::Int + ) + +Scatter blocked tensor stiffness matrix directly to triplet arrays **in-place**. + +Appends all (i, j, value) triplets from blocked tensor matrix to global triplet arrays. +This avoids the intermediate conversion to Float64 matrix. + +# Arguments +- `cache`: COO cache with triplet arrays +- `K_blocks`: Element stiffness as [N×N] matrix of 3×3 tensor blocks +- `dofs`: Global DOF indices [3*N] +- `N`: Number of nodes in element + +# Zero-Allocation Guarantee + +Writes directly to pre-allocated triplet arrays. No intermediate matrices created. + +# Algorithm + +```julia +for k in 1:N, l in 1:N + block = K_blocks[k, l] + k_offset = 3(k - 1) + l_offset = 3(l - 1) + for α in 1:3, β in 1:3 + counter += 1 + i_global = dofs[k_offset + α] + j_global = dofs[l_offset + β] + I[counter] = i_global + J[counter] = j_global + V[counter] = block[α, β] + end +end +``` +""" +@inline function scatter_blocks_to_triplets!( + cache::COOCache{EC,MC}, + K_blocks::Matrix{Tensor{2,3,Float64,9}}, + dofs::Vector{Int}, + N::Int +) where {EC,MC} + counter = cache.counter[] + ndofs_elem = 3 * N + + # Check capacity + new_triplets = ndofs_elem * ndofs_elem + if counter + new_triplets > cache.capacity + error("COO cache overflow: need $(counter + new_triplets) triplets, " * + "capacity is $(cache.capacity). Increase cache size.") + end + + # Scatter blocks directly to triplets + @inbounds for k in 1:N, l in 1:N + block = K_blocks[k, l] + k_offset = 3(k - 1) + l_offset = 3(l - 1) + for α in 1:3, β in 1:3 + counter += 1 + i_global = dofs[k_offset+α] + j_global = dofs[l_offset+β] + cache.I[counter] = i_global + cache.J[counter] = j_global + cache.V[counter] = block[α, β] + end + end + + cache.counter[] = counter + return nothing +end