JuliaFEM.jl/test/test_nodal_vs_element_assembly.jl

# Comparison tests: Nodal Assembly vs Traditional Element Assembly
#
# This file implements BOTH assembly methods for the same problem to:
# 1. Validate nodal assembly gives correct results
# 2. Compare performance characteristics
# 3. Demonstrate the architectural differences

using Test
using Tensors
using LinearAlgebra
using SparseArrays

include("../src/nodal_assembly_structures.jl")

# =============================================================================
# Test Problem: Linear Elasticity on Simple Tet4 Mesh
# =============================================================================

"""
Simple 3D linear elasticity material model for testing.
"""
struct TestLinearElastic
    E::Float64   # Young's modulus
    ν::Float64   # Poisson's ratio
end

function compute_lame_parameters(mat::TestLinearElastic)
    λ = mat.E * mat.ν / ((1 + mat.ν) * (1 - 2 * mat.ν))
    μ = mat.E / (2 * (1 + mat.ν))
    return λ, μ
end

"""
Compute 4th-order elasticity tensor using Tensors.jl.
"""
function elasticity_tensor(mat::TestLinearElastic)
    λ, μ = compute_lame_parameters(mat)

    # δ_ij δ_kl + μ(δ_ik δ_jl + δ_il δ_jk)
    δ = one(Tensor{2,3})  # Identity tensor
    I = one(SymmetricTensor{4,3})  # Symmetric 4th-order identity

    # C = λ δ ⊗ δ + 2μ I_sym
    C = λ * δ ⊗ δ + 2μ * I

    return C
end

"""
Compute stress from strain using linear elasticity.
"""
function compute_stress(mat::TestLinearElastic, ε::SymmetricTensor{2,3})
    C = elasticity_tensor(mat)
    return C ⊡ ε  # Double contraction
end

"""
Convert 6-component Voigt vector to symmetric tensor.
"""
function voigt_to_tensor(v::Vector{Float64})
    return SymmetricTensor{2,3}((v[1], v[4], v[6],
        v[4], v[2], v[5],
        v[6], v[5], v[3]))
end

"""
Convert symmetric tensor to 6-component Voigt vector.
"""
function tensor_to_voigt(σ::SymmetricTensor{2,3})
    return [σ[1, 1], σ[2, 2], σ[3, 3], σ[1, 2], σ[2, 3], σ[1, 3]]
end

# =============================================================================
# Traditional Element Assembly
# =============================================================================

"""
Compute B matrix (strain-displacement) for a single node in 3D.

Maps nodal displacements to strain via: ε = B * u

# Returns
- `B_node::Matrix{Float64}`: 6×3 matrix for this node
"""
function compute_B_matrix_node(dN::Vec{3,Float64})
    B = zeros(6, 3)

    # ε_11 = ∂u_x/∂x
    B[1, 1] = dN[1]
    # ε_22 = ∂u_y/∂y
    B[2, 2] = dN[2]
    # ε_33 = ∂u_z/∂z
    B[3, 3] = dN[3]
    # 2ε_12 = ∂u_x/∂y + ∂u_y/∂x
    B[4, 1] = dN[2]
    B[4, 2] = dN[1]
    # 2ε_23 = ∂u_y/∂z + ∂u_z/∂y
    B[5, 2] = dN[3]
    B[5, 3] = dN[2]
    # 2ε_13 = ∂u_x/∂z + ∂u_z/∂x
    B[6, 1] = dN[3]
    B[6, 3] = dN[1]

    return B
end

"""
Assemble element stiffness matrix using traditional element assembly.

K_e = ∫ B^T C B dV

# Arguments
- `X`: Element node coordinates [4×3 for Tet4]
- `material`: Material model
- `gauss_weight`: Integration weight (1/6 for 1-point Tet4)

# Returns
- `K_e::Matrix{Float64}`: 12×12 element stiffness matrix
"""
function assemble_element_stiffness_traditional(
    X::Matrix{Float64},  # 4×3 (4 nodes, 3 coords)
    material::TestLinearElastic,
    gauss_weight::Float64=1.0 / 6.0
)
    nnodes = 4
    ndofs = 12  # 4 nodes × 3 DOF

    # Compute Jacobian and shape function derivatives
    # For Tet4 at centroid (constant derivatives)
    J = zeros(3, 3)
    for i in 1:3
        J[i, :] = X[i+1, :] - X[1, :]
    end

    detJ = det(J)
    invJ = inv(J)

    # Shape function derivatives in reference coordinates
    dN_ref = [
        -1.0 -1.0 -1.0;  # Node 1
        1.0 0.0 0.0;  # Node 2
        0.0 1.0 0.0;  # Node 3
        0.0 0.0 1.0   # Node 4
    ]

    # Transform to physical coordinates: dN = dN_ref * inv(J)
    dN_physical = dN_ref * invJ  # 4×3

    # Build full B matrix (6×12)
    B = zeros(6, ndofs)
    for i in 1:nnodes
        dN_i = Vec{3}((dN_physical[i, 1], dN_physical[i, 2], dN_physical[i, 3]))
        B_i = compute_B_matrix_node(dN_i)
        B[:, 3*(i-1)+1:3*i] = B_i
    end

    # Material stiffness in Voigt notation
    C_tensor = elasticity_tensor(material)
    C = zeros(6, 6)
    for i in 1:3, j in 1:3, k in 1:3, l in 1:3
        # Map to Voigt indices
        voigt_ij = (i == j) ? i : (i + j == 3) ? 4 : (i + j == 4) ? 6 : 5
        voigt_kl = (k == l) ? k : (k + l == 3) ? 4 : (k + l == 4) ? 6 : 5
        C[voigt_ij, voigt_kl] = C_tensor[i, j, k, l]
    end

    # Element stiffness: K_e = w * |J| * B^T * C * B
    w = gauss_weight * abs(detJ)
    K_e = w * (B' * C * B)

    return K_e
end

"""
Assemble global stiffness matrix using traditional element assembly.

# Arguments
- `connectivity`: Element connectivity [(node1, node2, node3, node4), ...]
- `coordinates`: Nodal coordinates [nnodes×3]
- `material`: Material model

# Returns
- `K_global::SparseMatrixCSC`: Global stiffness matrix (nnodes*3 × nnodes*3)
"""
function assemble_global_traditional(
    connectivity::Vector{NTuple{4,Int}},
    coordinates::Matrix{Float64},  # nnodes×3
    material::TestLinearElastic
)
    nnodes = size(coordinates, 1)
    ndof_global = 3 * nnodes

    # Build sparse matrix using COO format
    I_rows = Int[]
    J_cols = Int[]
    values = Float64[]

    # Loop over elements
    for (elem_id, conn) in enumerate(connectivity)
        # Extract element coordinates
        X_elem = coordinates[collect(conn), :]  # 4×3

        # Compute element stiffness
        K_e = assemble_element_stiffness_traditional(X_elem, material)

        # Scatter to global (gather DOF indices)
        gdofs = zeros(Int, 12)
        for (local_i, global_node) in enumerate(conn)
            gdofs[3*(local_i-1)+1:3*local_i] = 3 * (global_node - 1) .+ (1:3)
        end

        # Add to COO lists
        for i in 1:12, j in 1:12
            if abs(K_e[i, j]) > 1e-14  # Skip near-zeros
                push!(I_rows, gdofs[i])
                push!(J_cols, gdofs[j])
                push!(values, K_e[i, j])
            end
        end
    end

    # Assemble sparse matrix (sums duplicate entries automatically)
    K_global = sparse(I_rows, J_cols, values, ndof_global, ndof_global)

    return K_global
end

# =============================================================================
# Nodal Assembly Implementation
# =============================================================================

"""
Compute 3×3 stiffness block between node i and node j in an element.

Simplified approach: compute element K and extract the block.

# Arguments
- `dN_i, dN_j`: Shape function gradients at nodes i and j
- `C_tensor`: 4th-order elasticity tensor
- `weight`: Integration weight × |J|

# Returns
- `K_ij::Tensor{2,3}`: 3×3 stiffness block
"""
function compute_stiffness_block(
    dN_i::Vec{3,Float64},
    dN_j::Vec{3,Float64},
    C_tensor::Union{SymmetricTensor{4,3,Float64},Tensor{4,3,Float64}},
    weight::Float64
)
    # Build B matrices for both nodes (6×3 in Voigt notation)
    B_i = compute_B_matrix_node(dN_i)
    B_j = compute_B_matrix_node(dN_j)

    # Convert C_tensor to 6×6 matrix
    C = zeros(6, 6)
    for i in 1:3, j in 1:3, k in 1:3, l in 1:3
        # Map to Voigt indices (engineering notation)
        voigt_ij = (i == j) ? i : (i + j == 3) ? 4 : (i + j == 4) ? 6 : 5
        voigt_kl = (k == l) ? k : (k + l == 3) ? 4 : (k + l == 4) ? 6 : 5
        C[voigt_ij, voigt_kl] = C_tensor[i, j, k, l]
    end

    # Compute block: K_ij = w * B_i^T * C * B_j
    K_block_mat = weight * (B_i' * C * B_j)  # 3×3 matrix

    # Convert to Tensor{2,3}
    K_ij = Tensor{2,3}((K_block_mat[1, 1], K_block_mat[1, 2], K_block_mat[1, 3],
        K_block_mat[2, 1], K_block_mat[2, 2], K_block_mat[2, 3],
        K_block_mat[3, 1], K_block_mat[3, 2], K_block_mat[3, 3]))

    return K_ij
end

"""
Assemble nodal contribution for a single node using nodal assembly.

# Arguments
- `node_id`: Global node ID
- `map`: Node-to-elements mapping
- `connectivity`: Element connectivity
- `coordinates`: Nodal coordinates
- `material`: Material model

# Returns
- `contrib::NodalStiffnessContribution`: Contains K_blocks for spider nodes
"""
function assemble_nodal_contribution(
    node_id::Int,
    map::NodeToElementsMap,
    connectivity::Vector{NTuple{4,Int}},
    coordinates::Matrix{Float64},
    material::TestLinearElastic
)
    # Get spider nodes
    spider = get_node_spider(map, node_id, connectivity)

    # Allocate storage
    contrib = NodalStiffnessContribution(node_id, spider, Float64)

    # Map spider nodes to indices for fast lookup
    spider_map = Dict(node => idx for (idx, node) in enumerate(spider))

    C_tensor = elasticity_tensor(material)

    # Loop over elements touching this node
    for elem_info in map.node_to_elements[node_id]
        elem_id = elem_info.element_id
        local_i = elem_info.local_node_idx

        conn = connectivity[elem_id]
        X_elem = coordinates[collect(conn), :]

        # Compute Jacobian
        J = zeros(3, 3)
        for i in 1:3
            J[i, :] = X_elem[i+1, :] - X_elem[1, :]
        end
        detJ = det(J)
        invJ = inv(J)

        # Shape function derivatives
        dN_ref = [
            -1.0 -1.0 -1.0;
            1.0 0.0 0.0;
            0.0 1.0 0.0;
            0.0 0.0 1.0
        ]
        dN_physical = dN_ref * invJ

        weight = (1.0 / 6.0) * abs(detJ)

        # Get gradient for our node
        dN_i = Vec{3}((dN_physical[local_i, 1],
            dN_physical[local_i, 2],
            dN_physical[local_i, 3]))

        # Loop over all nodes in this element
        for (local_j, global_j) in enumerate(conn)
            # Get gradient for node j
            dN_j = Vec{3}((dN_physical[local_j, 1],
                dN_physical[local_j, 2],
                dN_physical[local_j, 3]))

            # Compute 3×3 block K_ij
            K_block = compute_stiffness_block(dN_i, dN_j, C_tensor, weight)

            # Add to appropriate spider location
            spider_idx = spider_map[global_j]
            contrib.K_blocks[spider_idx] += K_block
        end
    end

    return contrib
end

"""
Assemble full matrix-vector product using nodal assembly.

# Returns
- `w::Vector{Float64}`: Result of K*u (nnodes*3)
"""
function matvec_nodal_assembly(
    map::NodeToElementsMap,
    connectivity::Vector{NTuple{4,Int}},
    coordinates::Matrix{Float64},
    material::TestLinearElastic,
    u::Vector{Float64}
)
    nnodes = size(coordinates, 1)
    w = zeros(3 * nnodes)

    # Convert u to Vec{3} per node
    u_nodal = [Vec{3}((u[3*(i-1)+1], u[3*(i-1)+2], u[3*(i-1)+3]))
               for i in 1:nnodes]

    # Loop over nodes
    for node_i in 1:nnodes
        # Assemble contribution for this node
        contrib = assemble_nodal_contribution(node_i, map, connectivity,
            coordinates, material)

        # Matrix-vector product for this node
        w_i = matrix_vector_product_nodal(contrib, u_nodal)

        # Store result
        w[3*(node_i-1)+1:3*node_i] = [w_i[1], w_i[2], w_i[3]]
    end

    return w
end

# =============================================================================
# Unit Tests
# =============================================================================

@testset "Nodal vs Element Assembly Comparison" begin

    @testset "Single Tet4 Element" begin
        # Simple single element test
        connectivity = [(1, 2, 3, 4)]

        # Unit tetrahedron
        coordinates = [
            0.0 0.0 0.0;
            1.0 0.0 0.0;
            0.0 1.0 0.0;
            0.0 0.0 1.0
        ]

        material = TestLinearElastic(200e3, 0.3)  # Steel-like

        # Traditional assembly
        K_traditional = assemble_global_traditional(connectivity, coordinates, material)

        # Nodal assembly
        map = NodeToElementsMap(connectivity)
        u_test = randn(12)  # Random displacement

        w_traditional = K_traditional * u_test
        w_nodal = matvec_nodal_assembly(map, connectivity, coordinates, material, u_test)

        # Should match exactly
        @test w_nodal ≈ w_traditional rtol = 1e-10

        println("\nSingle Tet4: ✓ Nodal and traditional assembly match")
    end

    @testset "Two Tet4 Elements Sharing Face" begin
        # Two tetrahedra sharing nodes 2,3,4
        connectivity = [
            (1, 2, 3, 4),
            (2, 3, 4, 5)
        ]

        coordinates = [
            0.0 0.0 0.0;   # Node 1
            1.0 0.0 0.0;   # Node 2
            0.0 1.0 0.0;   # Node 3
            0.0 0.0 1.0;   # Node 4
            1.0 1.0 1.0    # Node 5
        ]

        material = TestLinearElastic(200e3, 0.3)

        # Traditional assembly
        K_traditional = assemble_global_traditional(connectivity, coordinates, material)

        # Nodal assembly
        map = NodeToElementsMap(connectivity)
        u_test = randn(15)  # 5 nodes × 3 DOF

        w_traditional = K_traditional * u_test
        w_nodal = matvec_nodal_assembly(map, connectivity, coordinates, material, u_test)

        @test w_nodal ≈ w_traditional rtol = 1e-10

        println("Two Tet4s: ✓ Nodal and traditional assembly match")
    end

    @testset "Symmetry Check" begin
        connectivity = [(1, 2, 3, 4)]
        coordinates = [
            0.0 0.0 0.0;
            1.0 0.0 0.0;
            0.0 1.0 0.0;
            0.0 0.0 1.0
        ]
        material = TestLinearElastic(200e3, 0.3)

        K = assemble_global_traditional(connectivity, coordinates, material)

        # Stiffness should be symmetric
        @test norm(K - K') / norm(K) < 1e-12

        println("Symmetry: ✓ Stiffness matrix is symmetric")
    end

    @testset "Spider Pattern Verification" begin
        connectivity = [
            (1, 2, 3, 4),
            (2, 3, 4, 5)
        ]

        map = NodeToElementsMap(connectivity)

        # Node 1: corner, only in element 1
        spider_1 = get_node_spider(map, 1, connectivity)
        @test spider_1 == [1, 2, 3, 4]

        # Node 2: interior, in both elements
        spider_2 = get_node_spider(map, 2, connectivity)
        @test spider_2 == [1, 2, 3, 4, 5]

        println("Spider: ✓ Correct coupling pattern detected")
    end

    @testset "Zero Displacement → Zero Forces" begin
        connectivity = [(1, 2, 3, 4)]
        coordinates = [
            0.0 0.0 0.0;
            1.0 0.0 0.0;
            0.0 1.0 0.0;
            0.0 0.0 1.0
        ]
        material = TestLinearElastic(200e3, 0.3)
        map = NodeToElementsMap(connectivity)

        u_zero = zeros(12)
        w = matvec_nodal_assembly(map, connectivity, coordinates, material, u_zero)

        @test norm(w) < 1e-14

        println("Zero test: ✓ Zero displacement → zero forces")
    end

    @testset "Rigid Body Motion → Zero Forces" begin
        connectivity = [(1, 2, 3, 4)]
        coordinates = [
            0.0 0.0 0.0;
            1.0 0.0 0.0;
            0.0 1.0 0.0;
            0.0 0.0 1.0
        ]
        material = TestLinearElastic(200e3, 0.3)

        K = assemble_global_traditional(connectivity, coordinates, material)

        # Translation in x
        u_trans = repeat([1.0, 0.0, 0.0], 4)
        f = K * u_trans
        @test norm(f) < 1e-10  # Should be ~zero (within numerical error)

        println("Rigid body: ✓ Translation produces near-zero forces")
    end
end

println("\n" * "="^70)
println("SUMMARY: Nodal Assembly Implementation Validated ✓")
println("="^70)
println("All tests pass - nodal assembly matches traditional assembly exactly!")
println("Next: Performance benchmarking and GPU implementation")