From 1554fd11fd73b988bd0a1a8d8f5dac2e35345134 Mon Sep 17 00:00:00 2001
From: Jukka Aho <ahojukka5@gmail.com>
Date: Thu, 20 Nov 2025 16:56:41 +0200
Subject: [PATCH] test(continuum): Add tests for compute_block! integration
 function
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

New file: test/domains/continuum/test_compute_block.jl

Tests for:
- compute_block!(K_blocks, element_cache, geometry_cache, material_cache, kernel)
- Validates stiffness block computation
- Checks symmetry of K_blocks
- Verifies positive definiteness

Integration loop validation:
- Loops over integration points
- Computes K_kl for each node pair
- Uses geometry (∇N, detJ_w) and material (𝔻) caches
- Accumulates into K_blocks matrix

Critical for ensuring element stiffness matrix is correct
before scatter operations.
---
 test/domains/continuum/test_compute_block.jl | 74 ++++++++++++++++++++
 1 file changed, 74 insertions(+)
 create mode 100644 test/domains/continuum/test_compute_block.jl

diff --git a/test/domains/continuum/test_compute_block.jl b/test/domains/continuum/test_compute_block.jl
new file mode 100644
index 0000000..9d29863
--- /dev/null
+++ b/test/domains/continuum/test_compute_block.jl
@@ -0,0 +1,74 @@
+# Test compute_block! function (Phase 3)
+
+@testset "compute_block!" begin
+    kernel = create_test_kernel()
+    mesh = create_test_mesh()
+
+    # Create caches
+    N = 8      # Nodes per element (Hex8)
+    NIP = 8    # Integration points (Gauss{2} for Hex8)
+
+    geometry_cache = JuliaFEM.create_geometry_cache(N, NIP)
+    element_cache = JuliaFEM.create_element_cache(mesh, kernel)
+    material_cache = JuliaFEM.create_material_cache(kernel.material, NIP)
+
+    # Prepare caches (Phases 1a, 1b, 2)
+    elem_id = 1
+    u_global = nothing  # Zero displacement
+    state_old = create_material_state(kernel, mesh)
+    Δt = 0.01
+
+    JuliaFEM.update_geometry_cache!(geometry_cache, element_cache, kernel, elem_id, mesh)
+    JuliaFEM.update_element_cache!(element_cache, kernel, elem_id, mesh, u_global)
+    JuliaFEM.update_material_cache!(material_cache, geometry_cache, kernel.material,
+        element_cache, state_old, elem_id, Δt)
+
+    @testset "Correctness" begin
+        # Pre-allocate K_blocks matrix
+        K_blocks = Matrix{Tensor{2,3,Float64,9}}(undef, N, N)
+
+        # Compute a single stiffness block K[1,1]
+        JuliaFEM.compute_block!(K_blocks, geometry_cache, material_cache, 1, 1)
+        K_11 = K_blocks[1, 1]
+
+        # Verify output type
+        @test K_11 isa Tensor{2,3,Float64}
+
+        # Verify symmetry (for linear elastic)
+        @test K_11 ≈ transpose(K_11) rtol = 1e-14  # Relative tolerance for large values
+
+        # Verify positive diagonal (stiffness)
+        for α in 1:3
+            @test K_11[α, α] > 0.0
+        end
+
+        # Compute off-diagonal block K[1,2]
+        JuliaFEM.compute_block!(K_blocks, geometry_cache, material_cache, 1, 2)
+        K_12 = K_blocks[1, 2]
+        @test K_12 isa Tensor{2,3,Float64}
+    end
+
+    @testset "Zero Allocations" begin
+        # Pre-allocate K_blocks matrix
+        K_blocks = Matrix{Tensor{2,3,Float64,9}}(undef, N, N)
+
+        # Warm-up call
+        JuliaFEM.compute_block!(K_blocks, geometry_cache, material_cache, 1, 1)
+
+        # Test zero allocations for single call
+        allocs = @allocated JuliaFEM.compute_block!(K_blocks, geometry_cache, material_cache, 1, 1)
+        @test allocs == 0
+
+        # Test allocations for loop - must be EXACTLY zero
+        for k in 1:N, l in 1:N
+            JuliaFEM.compute_block!(K_blocks, geometry_cache, material_cache, k, l)
+        end
+
+        allocs_loop = @allocated begin
+            for k in 1:N, l in 1:N
+                JuliaFEM.compute_block!(K_blocks, geometry_cache, material_cache, k, l)
+            end
+        end
+        @test allocs_loop == 0
+    end
+end