using Base.Iterators
using Scapin.Brick

function expected_avg_gradient_operator(h::NTuple{d,T}) where {d,T<:Number}
    ℒ = cell_vertices(d)
    D = zeros(T, d, fill(2, d)...)
    for i = 1:d, p in ℒ
        D[i, p] = (-1)^p[i] / 2^(d - 1) / h[i]
    end
    return D
end

function expected_avg_strain_displacement_operator(h::NTuple{d,T}) where {d,T<:Number}
    ℒ = cell_vertices(d)
    D = expected_avg_gradient_operator(h)
    B = zeros(T, d, d, d, fill(2, d)...)
    for i = 1:d, j = 1:d, k = 1:d, p in ℒ
        B[i, j, k, p] = (T(j == k) * D[i, p] + T(i == k) * D[j, p]) / 2
    end
    return B
end

function expected_K_λ(h::NTuple{d,T}) where {d,T<:Number}
    K_λ = zeros(Float64, d, fill(2, d)..., d, fill(2, d)...)
    if d == 2
        value = -1 / 4
        K_λ[1, 1, 1, 2, 1, 2] = value
        K_λ[1, 1, 1, 2, 2, 2] = value
        K_λ[1, 1, 2, 2, 1, 2] = value
        K_λ[1, 1, 2, 2, 2, 2] = value
        K_λ[1, 2, 1, 2, 1, 1] = value
        K_λ[1, 2, 1, 2, 2, 1] = value
        K_λ[1, 2, 2, 2, 1, 1] = value
        K_λ[1, 2, 2, 2, 2, 1] = value
        K_λ[2, 1, 1, 1, 2, 1] = value
        K_λ[2, 1, 1, 1, 2, 2] = value
        K_λ[2, 1, 2, 1, 1, 1] = value
        K_λ[2, 1, 2, 1, 1, 2] = value
        K_λ[2, 2, 1, 1, 2, 1] = value
        K_λ[2, 2, 1, 1, 2, 2] = value
        K_λ[2, 2, 2, 1, 1, 1] = value
        K_λ[2, 2, 2, 1, 1, 2] = value

        value = 1 / 4
        K_λ[1, 1, 1, 2, 1, 1] = value
        K_λ[1, 1, 1, 2, 2, 1] = value
        K_λ[1, 1, 2, 2, 1, 1] = value
        K_λ[1, 1, 2, 2, 2, 1] = value
        K_λ[1, 2, 1, 2, 1, 2] = value
        K_λ[1, 2, 1, 2, 2, 2] = value
        K_λ[1, 2, 2, 2, 1, 2] = value
        K_λ[1, 2, 2, 2, 2, 2] = value
        K_λ[2, 1, 1, 1, 1, 1] = value
        K_λ[2, 1, 1, 1, 1, 2] = value
        K_λ[2, 1, 2, 1, 2, 1] = value
        K_λ[2, 1, 2, 1, 2, 2] = value
        K_λ[2, 2, 1, 1, 1, 1] = value
        K_λ[2, 2, 1, 1, 1, 2] = value
        K_λ[2, 2, 2, 1, 2, 1] = value
        K_λ[2, 2, 2, 1, 2, 2] = value

        value = -h[1] / (3 * h[2])
        K_λ[2, 1, 1, 2, 1, 2] = value
        K_λ[2, 1, 2, 2, 1, 1] = value
        K_λ[2, 2, 1, 2, 2, 2] = value
        K_λ[2, 2, 2, 2, 2, 1] = value

        value = -h[1] / (6 * h[2])
        K_λ[2, 1, 1, 2, 2, 2] = value
        K_λ[2, 1, 2, 2, 2, 1] = value
        K_λ[2, 2, 1, 2, 1, 2] = value
        K_λ[2, 2, 2, 2, 1, 1] = value

        value = h[1] / (6 * h[2])
        K_λ[2, 1, 1, 2, 2, 1] = value
        K_λ[2, 1, 2, 2, 2, 2] = value
        K_λ[2, 2, 1, 2, 1, 1] = value
        K_λ[2, 2, 2, 2, 1, 2] = value

        value = h[1] / (3 * h[2])
        K_λ[2, 1, 1, 2, 1, 1] = value
        K_λ[2, 1, 2, 2, 1, 2] = value
        K_λ[2, 2, 1, 2, 2, 1] = value
        K_λ[2, 2, 2, 2, 2, 2] = value

        value = -h[2] / (3 * h[1])
        K_λ[1, 1, 1, 1, 2, 1] = value
        K_λ[1, 1, 2, 1, 2, 2] = value
        K_λ[1, 2, 1, 1, 1, 1] = value
        K_λ[1, 2, 2, 1, 1, 2] = value

        value = -h[2] / (6 * h[1])
        K_λ[1, 1, 1, 1, 2, 2] = value
        K_λ[1, 1, 2, 1, 2, 1] = value
        K_λ[1, 2, 1, 1, 1, 2] = value
        K_λ[1, 2, 2, 1, 1, 1] = value

        value = h[2] / (6 * h[1])
        K_λ[1, 1, 1, 1, 1, 2] = value
        K_λ[1, 1, 2, 1, 1, 1] = value
        K_λ[1, 2, 1, 1, 2, 2] = value
        K_λ[1, 2, 2, 1, 2, 1] = value

        value = h[2] / (3 * h[1])
        K_λ[1, 1, 1, 1, 1, 1] = value
        K_λ[1, 1, 2, 1, 1, 2] = value
        K_λ[1, 2, 1, 1, 2, 1] = value
        K_λ[1, 2, 2, 1, 2, 2] = value
    elseif d == 3
        value = -h[1] / 12
        K_λ[2, 1, 1, 1, 3, 1, 1, 2] = value
        K_λ[2, 1, 1, 1, 3, 1, 2, 2] = value
        K_λ[2, 1, 1, 2, 3, 1, 1, 2] = value
        K_λ[2, 1, 1, 2, 3, 1, 2, 2] = value
        K_λ[2, 1, 2, 1, 3, 1, 1, 1] = value
        K_λ[2, 1, 2, 1, 3, 1, 2, 1] = value
        K_λ[2, 1, 2, 2, 3, 1, 1, 1] = value
        K_λ[2, 1, 2, 2, 3, 1, 2, 1] = value
        K_λ[2, 2, 1, 1, 3, 2, 1, 2] = value
        K_λ[2, 2, 1, 1, 3, 2, 2, 2] = value
        K_λ[2, 2, 1, 2, 3, 2, 1, 2] = value
        K_λ[2, 2, 1, 2, 3, 2, 2, 2] = value
        K_λ[2, 2, 2, 1, 3, 2, 1, 1] = value
        K_λ[2, 2, 2, 1, 3, 2, 2, 1] = value
        K_λ[2, 2, 2, 2, 3, 2, 1, 1] = value
        K_λ[2, 2, 2, 2, 3, 2, 2, 1] = value
        K_λ[3, 1, 1, 1, 2, 1, 2, 1] = value
        K_λ[3, 1, 1, 1, 2, 1, 2, 2] = value
        K_λ[3, 1, 1, 2, 2, 1, 1, 1] = value
        K_λ[3, 1, 1, 2, 2, 1, 1, 2] = value
        K_λ[3, 1, 2, 1, 2, 1, 2, 1] = value
        K_λ[3, 1, 2, 1, 2, 1, 2, 2] = value
        K_λ[3, 1, 2, 2, 2, 1, 1, 1] = value
        K_λ[3, 1, 2, 2, 2, 1, 1, 2] = value
        K_λ[3, 2, 1, 1, 2, 2, 2, 1] = value
        K_λ[3, 2, 1, 1, 2, 2, 2, 2] = value
        K_λ[3, 2, 1, 2, 2, 2, 1, 1] = value
        K_λ[3, 2, 1, 2, 2, 2, 1, 2] = value
        K_λ[3, 2, 2, 1, 2, 2, 2, 1] = value
        K_λ[3, 2, 2, 1, 2, 2, 2, 2] = value
        K_λ[3, 2, 2, 2, 2, 2, 1, 1] = value
        K_λ[3, 2, 2, 2, 2, 2, 1, 2] = value

        value = -h[1] / 24
        K_λ[2, 1, 1, 1, 3, 2, 1, 2] = value
        K_λ[2, 1, 1, 1, 3, 2, 2, 2] = value
        K_λ[2, 1, 1, 2, 3, 2, 1, 2] = value
        K_λ[2, 1, 1, 2, 3, 2, 2, 2] = value
        K_λ[2, 1, 2, 1, 3, 2, 1, 1] = value
        K_λ[2, 1, 2, 1, 3, 2, 2, 1] = value
        K_λ[2, 1, 2, 2, 3, 2, 1, 1] = value
        K_λ[2, 1, 2, 2, 3, 2, 2, 1] = value
        K_λ[2, 2, 1, 1, 3, 1, 1, 2] = value
        K_λ[2, 2, 1, 1, 3, 1, 2, 2] = value
        K_λ[2, 2, 1, 2, 3, 1, 1, 2] = value
        K_λ[2, 2, 1, 2, 3, 1, 2, 2] = value
        K_λ[2, 2, 2, 1, 3, 1, 1, 1] = value
        K_λ[2, 2, 2, 1, 3, 1, 2, 1] = value
        K_λ[2, 2, 2, 2, 3, 1, 1, 1] = value
        K_λ[2, 2, 2, 2, 3, 1, 2, 1] = value
        K_λ[3, 1, 1, 1, 2, 2, 2, 1] = value
        K_λ[3, 1, 1, 1, 2, 2, 2, 2] = value
        K_λ[3, 1, 1, 2, 2, 2, 1, 1] = value
        K_λ[3, 1, 1, 2, 2, 2, 1, 2] = value
        K_λ[3, 1, 2, 1, 2, 2, 2, 1] = value
        K_λ[3, 1, 2, 1, 2, 2, 2, 2] = value
        K_λ[3, 1, 2, 2, 2, 2, 1, 1] = value
        K_λ[3, 1, 2, 2, 2, 2, 1, 2] = value
        K_λ[3, 2, 1, 1, 2, 1, 2, 1] = value
        K_λ[3, 2, 1, 1, 2, 1, 2, 2] = value
        K_λ[3, 2, 1, 2, 2, 1, 1, 1] = value
        K_λ[3, 2, 1, 2, 2, 1, 1, 2] = value
        K_λ[3, 2, 2, 1, 2, 1, 2, 1] = value
        K_λ[3, 2, 2, 1, 2, 1, 2, 2] = value
        K_λ[3, 2, 2, 2, 2, 1, 1, 1] = value
        K_λ[3, 2, 2, 2, 2, 1, 1, 2] = value

        value = h[1] / 24
        K_λ[2, 1, 1, 1, 3, 2, 1, 1] = value
        K_λ[2, 1, 1, 1, 3, 2, 2, 1] = value
        K_λ[2, 1, 1, 2, 3, 2, 1, 1] = value
        K_λ[2, 1, 1, 2, 3, 2, 2, 1] = value
        K_λ[2, 1, 2, 1, 3, 2, 1, 2] = value
        K_λ[2, 1, 2, 1, 3, 2, 2, 2] = value
        K_λ[2, 1, 2, 2, 3, 2, 1, 2] = value
        K_λ[2, 1, 2, 2, 3, 2, 2, 2] = value
        K_λ[2, 2, 1, 1, 3, 1, 1, 1] = value
        K_λ[2, 2, 1, 1, 3, 1, 2, 1] = value
        K_λ[2, 2, 1, 2, 3, 1, 1, 1] = value
        K_λ[2, 2, 1, 2, 3, 1, 2, 1] = value
        K_λ[2, 2, 2, 1, 3, 1, 1, 2] = value
        K_λ[2, 2, 2, 1, 3, 1, 2, 2] = value
        K_λ[2, 2, 2, 2, 3, 1, 1, 2] = value
        K_λ[2, 2, 2, 2, 3, 1, 2, 2] = value
        K_λ[3, 1, 1, 1, 2, 2, 1, 1] = value
        K_λ[3, 1, 1, 1, 2, 2, 1, 2] = value
        K_λ[3, 1, 1, 2, 2, 2, 2, 1] = value
        K_λ[3, 1, 1, 2, 2, 2, 2, 2] = value
        K_λ[3, 1, 2, 1, 2, 2, 1, 1] = value
        K_λ[3, 1, 2, 1, 2, 2, 1, 2] = value
        K_λ[3, 1, 2, 2, 2, 2, 2, 1] = value
        K_λ[3, 1, 2, 2, 2, 2, 2, 2] = value
        K_λ[3, 2, 1, 1, 2, 1, 1, 1] = value
        K_λ[3, 2, 1, 1, 2, 1, 1, 2] = value
        K_λ[3, 2, 1, 2, 2, 1, 2, 1] = value
        K_λ[3, 2, 1, 2, 2, 1, 2, 2] = value
        K_λ[3, 2, 2, 1, 2, 1, 1, 1] = value
        K_λ[3, 2, 2, 1, 2, 1, 1, 2] = value
        K_λ[3, 2, 2, 2, 2, 1, 2, 1] = value
        K_λ[3, 2, 2, 2, 2, 1, 2, 2] = value

        value = h[1] / 12
        K_λ[2, 1, 1, 1, 3, 1, 1, 1] = value
        K_λ[2, 1, 1, 1, 3, 1, 2, 1] = value
        K_λ[2, 1, 1, 2, 3, 1, 1, 1] = value
        K_λ[2, 1, 1, 2, 3, 1, 2, 1] = value
        K_λ[2, 1, 2, 1, 3, 1, 1, 2] = value
        K_λ[2, 1, 2, 1, 3, 1, 2, 2] = value
        K_λ[2, 1, 2, 2, 3, 1, 1, 2] = value
        K_λ[2, 1, 2, 2, 3, 1, 2, 2] = value
        K_λ[2, 2, 1, 1, 3, 2, 1, 1] = value
        K_λ[2, 2, 1, 1, 3, 2, 2, 1] = value
        K_λ[2, 2, 1, 2, 3, 2, 1, 1] = value
        K_λ[2, 2, 1, 2, 3, 2, 2, 1] = value
        K_λ[2, 2, 2, 1, 3, 2, 1, 2] = value
        K_λ[2, 2, 2, 1, 3, 2, 2, 2] = value
        K_λ[2, 2, 2, 2, 3, 2, 1, 2] = value
        K_λ[2, 2, 2, 2, 3, 2, 2, 2] = value
        K_λ[3, 1, 1, 1, 2, 1, 1, 1] = value
        K_λ[3, 1, 1, 1, 2, 1, 1, 2] = value
        K_λ[3, 1, 1, 2, 2, 1, 2, 1] = value
        K_λ[3, 1, 1, 2, 2, 1, 2, 2] = value
        K_λ[3, 1, 2, 1, 2, 1, 1, 1] = value
        K_λ[3, 1, 2, 1, 2, 1, 1, 2] = value
        K_λ[3, 1, 2, 2, 2, 1, 2, 1] = value
        K_λ[3, 1, 2, 2, 2, 1, 2, 2] = value
        K_λ[3, 2, 1, 1, 2, 2, 1, 1] = value
        K_λ[3, 2, 1, 1, 2, 2, 1, 2] = value
        K_λ[3, 2, 1, 2, 2, 2, 2, 1] = value
        K_λ[3, 2, 1, 2, 2, 2, 2, 2] = value
        K_λ[3, 2, 2, 1, 2, 2, 1, 1] = value
        K_λ[3, 2, 2, 1, 2, 2, 1, 2] = value
        K_λ[3, 2, 2, 2, 2, 2, 2, 1] = value
        K_λ[3, 2, 2, 2, 2, 2, 2, 2] = value

        value = -h[2] / 12
        K_λ[1, 1, 1, 1, 3, 1, 1, 2] = value
        K_λ[1, 1, 1, 1, 3, 2, 1, 2] = value
        K_λ[1, 1, 1, 2, 3, 1, 1, 2] = value
        K_λ[1, 1, 1, 2, 3, 2, 1, 2] = value
        K_λ[1, 1, 2, 1, 3, 1, 2, 2] = value
        K_λ[1, 1, 2, 1, 3, 2, 2, 2] = value
        K_λ[1, 1, 2, 2, 3, 1, 2, 2] = value
        K_λ[1, 1, 2, 2, 3, 2, 2, 2] = value
        K_λ[1, 2, 1, 1, 3, 1, 1, 1] = value
        K_λ[1, 2, 1, 1, 3, 2, 1, 1] = value
        K_λ[1, 2, 1, 2, 3, 1, 1, 1] = value
        K_λ[1, 2, 1, 2, 3, 2, 1, 1] = value
        K_λ[1, 2, 2, 1, 3, 1, 2, 1] = value
        K_λ[1, 2, 2, 1, 3, 2, 2, 1] = value
        K_λ[1, 2, 2, 2, 3, 1, 2, 1] = value
        K_λ[1, 2, 2, 2, 3, 2, 2, 1] = value
        K_λ[3, 1, 1, 1, 1, 2, 1, 1] = value
        K_λ[3, 1, 1, 1, 1, 2, 1, 2] = value
        K_λ[3, 1, 1, 2, 1, 1, 1, 1] = value
        K_λ[3, 1, 1, 2, 1, 1, 1, 2] = value
        K_λ[3, 1, 2, 1, 1, 2, 2, 1] = value
        K_λ[3, 1, 2, 1, 1, 2, 2, 2] = value
        K_λ[3, 1, 2, 2, 1, 1, 2, 1] = value
        K_λ[3, 1, 2, 2, 1, 1, 2, 2] = value
        K_λ[3, 2, 1, 1, 1, 2, 1, 1] = value
        K_λ[3, 2, 1, 1, 1, 2, 1, 2] = value
        K_λ[3, 2, 1, 2, 1, 1, 1, 1] = value
        K_λ[3, 2, 1, 2, 1, 1, 1, 2] = value
        K_λ[3, 2, 2, 1, 1, 2, 2, 1] = value
        K_λ[3, 2, 2, 1, 1, 2, 2, 2] = value
        K_λ[3, 2, 2, 2, 1, 1, 2, 1] = value
        K_λ[3, 2, 2, 2, 1, 1, 2, 2] = value

        value = -h[2] / 24
        K_λ[1, 1, 1, 1, 3, 1, 2, 2] = value
        K_λ[1, 1, 1, 1, 3, 2, 2, 2] = value
        K_λ[1, 1, 1, 2, 3, 1, 2, 2] = value
        K_λ[1, 1, 1, 2, 3, 2, 2, 2] = value
        K_λ[1, 1, 2, 1, 3, 1, 1, 2] = value
        K_λ[1, 1, 2, 1, 3, 2, 1, 2] = value
        K_λ[1, 1, 2, 2, 3, 1, 1, 2] = value
        K_λ[1, 1, 2, 2, 3, 2, 1, 2] = value
        K_λ[1, 2, 1, 1, 3, 1, 2, 1] = value
        K_λ[1, 2, 1, 1, 3, 2, 2, 1] = value
        K_λ[1, 2, 1, 2, 3, 1, 2, 1] = value
        K_λ[1, 2, 1, 2, 3, 2, 2, 1] = value
        K_λ[1, 2, 2, 1, 3, 1, 1, 1] = value
        K_λ[1, 2, 2, 1, 3, 2, 1, 1] = value
        K_λ[1, 2, 2, 2, 3, 1, 1, 1] = value
        K_λ[1, 2, 2, 2, 3, 2, 1, 1] = value
        K_λ[3, 1, 1, 1, 1, 2, 2, 1] = value
        K_λ[3, 1, 1, 1, 1, 2, 2, 2] = value
        K_λ[3, 1, 1, 2, 1, 1, 2, 1] = value
        K_λ[3, 1, 1, 2, 1, 1, 2, 2] = value
        K_λ[3, 1, 2, 1, 1, 2, 1, 1] = value
        K_λ[3, 1, 2, 1, 1, 2, 1, 2] = value
        K_λ[3, 1, 2, 2, 1, 1, 1, 1] = value
        K_λ[3, 1, 2, 2, 1, 1, 1, 2] = value
        K_λ[3, 2, 1, 1, 1, 2, 2, 1] = value
        K_λ[3, 2, 1, 1, 1, 2, 2, 2] = value
        K_λ[3, 2, 1, 2, 1, 1, 2, 1] = value
        K_λ[3, 2, 1, 2, 1, 1, 2, 2] = value
        K_λ[3, 2, 2, 1, 1, 2, 1, 1] = value
        K_λ[3, 2, 2, 1, 1, 2, 1, 2] = value
        K_λ[3, 2, 2, 2, 1, 1, 1, 1] = value
        K_λ[3, 2, 2, 2, 1, 1, 1, 2] = value

        value = h[2] / 24
        K_λ[1, 1, 1, 1, 3, 1, 2, 1] = value
        K_λ[1, 1, 1, 1, 3, 2, 2, 1] = value
        K_λ[1, 1, 1, 2, 3, 1, 2, 1] = value
        K_λ[1, 1, 1, 2, 3, 2, 2, 1] = value
        K_λ[1, 1, 2, 1, 3, 1, 1, 1] = value
        K_λ[1, 1, 2, 1, 3, 2, 1, 1] = value
        K_λ[1, 1, 2, 2, 3, 1, 1, 1] = value
        K_λ[1, 1, 2, 2, 3, 2, 1, 1] = value
        K_λ[1, 2, 1, 1, 3, 1, 2, 2] = value
        K_λ[1, 2, 1, 1, 3, 2, 2, 2] = value
        K_λ[1, 2, 1, 2, 3, 1, 2, 2] = value
        K_λ[1, 2, 1, 2, 3, 2, 2, 2] = value
        K_λ[1, 2, 2, 1, 3, 1, 1, 2] = value
        K_λ[1, 2, 2, 1, 3, 2, 1, 2] = value
        K_λ[1, 2, 2, 2, 3, 1, 1, 2] = value
        K_λ[1, 2, 2, 2, 3, 2, 1, 2] = value
        K_λ[3, 1, 1, 1, 1, 1, 2, 1] = value
        K_λ[3, 1, 1, 1, 1, 1, 2, 2] = value
        K_λ[3, 1, 1, 2, 1, 2, 2, 1] = value
        K_λ[3, 1, 1, 2, 1, 2, 2, 2] = value
        K_λ[3, 1, 2, 1, 1, 1, 1, 1] = value
        K_λ[3, 1, 2, 1, 1, 1, 1, 2] = value
        K_λ[3, 1, 2, 2, 1, 2, 1, 1] = value
        K_λ[3, 1, 2, 2, 1, 2, 1, 2] = value
        K_λ[3, 2, 1, 1, 1, 1, 2, 1] = value
        K_λ[3, 2, 1, 1, 1, 1, 2, 2] = value
        K_λ[3, 2, 1, 2, 1, 2, 2, 1] = value
        K_λ[3, 2, 1, 2, 1, 2, 2, 2] = value
        K_λ[3, 2, 2, 1, 1, 1, 1, 1] = value
        K_λ[3, 2, 2, 1, 1, 1, 1, 2] = value
        K_λ[3, 2, 2, 2, 1, 2, 1, 1] = value
        K_λ[3, 2, 2, 2, 1, 2, 1, 2] = value

        value = h[2] / 12
        K_λ[1, 1, 1, 1, 3, 1, 1, 1] = value
        K_λ[1, 1, 1, 1, 3, 2, 1, 1] = value
        K_λ[1, 1, 1, 2, 3, 1, 1, 1] = value
        K_λ[1, 1, 1, 2, 3, 2, 1, 1] = value
        K_λ[1, 1, 2, 1, 3, 1, 2, 1] = value
        K_λ[1, 1, 2, 1, 3, 2, 2, 1] = value
        K_λ[1, 1, 2, 2, 3, 1, 2, 1] = value
        K_λ[1, 1, 2, 2, 3, 2, 2, 1] = value
        K_λ[1, 2, 1, 1, 3, 1, 1, 2] = value
        K_λ[1, 2, 1, 1, 3, 2, 1, 2] = value
        K_λ[1, 2, 1, 2, 3, 1, 1, 2] = value
        K_λ[1, 2, 1, 2, 3, 2, 1, 2] = value
        K_λ[1, 2, 2, 1, 3, 1, 2, 2] = value
        K_λ[1, 2, 2, 1, 3, 2, 2, 2] = value
        K_λ[1, 2, 2, 2, 3, 1, 2, 2] = value
        K_λ[1, 2, 2, 2, 3, 2, 2, 2] = value
        K_λ[3, 1, 1, 1, 1, 1, 1, 1] = value
        K_λ[3, 1, 1, 1, 1, 1, 1, 2] = value
        K_λ[3, 1, 1, 2, 1, 2, 1, 1] = value
        K_λ[3, 1, 1, 2, 1, 2, 1, 2] = value
        K_λ[3, 1, 2, 1, 1, 1, 2, 1] = value
        K_λ[3, 1, 2, 1, 1, 1, 2, 2] = value
        K_λ[3, 1, 2, 2, 1, 2, 2, 1] = value
        K_λ[3, 1, 2, 2, 1, 2, 2, 2] = value
        K_λ[3, 2, 1, 1, 1, 1, 1, 1] = value
        K_λ[3, 2, 1, 1, 1, 1, 1, 2] = value
        K_λ[3, 2, 1, 2, 1, 2, 1, 1] = value
        K_λ[3, 2, 1, 2, 1, 2, 1, 2] = value
        K_λ[3, 2, 2, 1, 1, 1, 2, 1] = value
        K_λ[3, 2, 2, 1, 1, 1, 2, 2] = value
        K_λ[3, 2, 2, 2, 1, 2, 2, 1] = value
        K_λ[3, 2, 2, 2, 1, 2, 2, 2] = value

        value = -(h[1] * h[2]) / (9 * h[3])
        K_λ[3, 1, 1, 1, 3, 1, 1, 2] = value
        K_λ[3, 1, 1, 2, 3, 1, 1, 1] = value
        K_λ[3, 1, 2, 1, 3, 1, 2, 2] = value
        K_λ[3, 1, 2, 2, 3, 1, 2, 1] = value
        K_λ[3, 2, 1, 1, 3, 2, 1, 2] = value
        K_λ[3, 2, 1, 2, 3, 2, 1, 1] = value
        K_λ[3, 2, 2, 1, 3, 2, 2, 2] = value
        K_λ[3, 2, 2, 2, 3, 2, 2, 1] = value

        value = -(h[1] * h[2]) / (18 * h[3])
        K_λ[3, 1, 1, 1, 3, 1, 2, 2] = value
        K_λ[3, 1, 1, 1, 3, 2, 1, 2] = value
        K_λ[3, 1, 1, 2, 3, 1, 2, 1] = value
        K_λ[3, 1, 1, 2, 3, 2, 1, 1] = value
        K_λ[3, 1, 2, 1, 3, 1, 1, 2] = value
        K_λ[3, 1, 2, 1, 3, 2, 2, 2] = value
        K_λ[3, 1, 2, 2, 3, 1, 1, 1] = value
        K_λ[3, 1, 2, 2, 3, 2, 2, 1] = value
        K_λ[3, 2, 1, 1, 3, 1, 1, 2] = value
        K_λ[3, 2, 1, 1, 3, 2, 2, 2] = value
        K_λ[3, 2, 1, 2, 3, 1, 1, 1] = value
        K_λ[3, 2, 1, 2, 3, 2, 2, 1] = value
        K_λ[3, 2, 2, 1, 3, 1, 2, 2] = value
        K_λ[3, 2, 2, 1, 3, 2, 1, 2] = value
        K_λ[3, 2, 2, 2, 3, 1, 2, 1] = value
        K_λ[3, 2, 2, 2, 3, 2, 1, 1] = value

        value = -(h[1] * h[2]) / (36 * h[3])
        K_λ[3, 1, 1, 1, 3, 2, 2, 2] = value
        K_λ[3, 1, 1, 2, 3, 2, 2, 1] = value
        K_λ[3, 1, 2, 1, 3, 2, 1, 2] = value
        K_λ[3, 1, 2, 2, 3, 2, 1, 1] = value
        K_λ[3, 2, 1, 1, 3, 1, 2, 2] = value
        K_λ[3, 2, 1, 2, 3, 1, 2, 1] = value
        K_λ[3, 2, 2, 1, 3, 1, 1, 2] = value
        K_λ[3, 2, 2, 2, 3, 1, 1, 1] = value

        value = (h[1] * h[2]) / (36 * h[3])
        K_λ[3, 1, 1, 1, 3, 2, 2, 1] = value
        K_λ[3, 1, 1, 2, 3, 2, 2, 2] = value
        K_λ[3, 1, 2, 1, 3, 2, 1, 1] = value
        K_λ[3, 1, 2, 2, 3, 2, 1, 2] = value
        K_λ[3, 2, 1, 1, 3, 1, 2, 1] = value
        K_λ[3, 2, 1, 2, 3, 1, 2, 2] = value
        K_λ[3, 2, 2, 1, 3, 1, 1, 1] = value
        K_λ[3, 2, 2, 2, 3, 1, 1, 2] = value

        value = (h[1] * h[2]) / (18 * h[3])
        K_λ[3, 1, 1, 1, 3, 1, 2, 1] = value
        K_λ[3, 1, 1, 1, 3, 2, 1, 1] = value
        K_λ[3, 1, 1, 2, 3, 1, 2, 2] = value
        K_λ[3, 1, 1, 2, 3, 2, 1, 2] = value
        K_λ[3, 1, 2, 1, 3, 1, 1, 1] = value
        K_λ[3, 1, 2, 1, 3, 2, 2, 1] = value
        K_λ[3, 1, 2, 2, 3, 1, 1, 2] = value
        K_λ[3, 1, 2, 2, 3, 2, 2, 2] = value
        K_λ[3, 2, 1, 1, 3, 1, 1, 1] = value
        K_λ[3, 2, 1, 1, 3, 2, 2, 1] = value
        K_λ[3, 2, 1, 2, 3, 1, 1, 2] = value
        K_λ[3, 2, 1, 2, 3, 2, 2, 2] = value
        K_λ[3, 2, 2, 1, 3, 1, 2, 1] = value
        K_λ[3, 2, 2, 1, 3, 2, 1, 1] = value
        K_λ[3, 2, 2, 2, 3, 1, 2, 2] = value
        K_λ[3, 2, 2, 2, 3, 2, 1, 2] = value

        value = (h[1] * h[2]) / (9 * h[3])
        K_λ[3, 1, 1, 1, 3, 1, 1, 1] = value
        K_λ[3, 1, 1, 2, 3, 1, 1, 2] = value
        K_λ[3, 1, 2, 1, 3, 1, 2, 1] = value
        K_λ[3, 1, 2, 2, 3, 1, 2, 2] = value
        K_λ[3, 2, 1, 1, 3, 2, 1, 1] = value
        K_λ[3, 2, 1, 2, 3, 2, 1, 2] = value
        K_λ[3, 2, 2, 1, 3, 2, 2, 1] = value
        K_λ[3, 2, 2, 2, 3, 2, 2, 2] = value

        value = -h[3] / 12
        K_λ[1, 1, 1, 1, 2, 1, 2, 1] = value
        K_λ[1, 1, 1, 1, 2, 2, 2, 1] = value
        K_λ[1, 1, 1, 2, 2, 1, 2, 2] = value
        K_λ[1, 1, 1, 2, 2, 2, 2, 2] = value
        K_λ[1, 1, 2, 1, 2, 1, 2, 1] = value
        K_λ[1, 1, 2, 1, 2, 2, 2, 1] = value
        K_λ[1, 1, 2, 2, 2, 1, 2, 2] = value
        K_λ[1, 1, 2, 2, 2, 2, 2, 2] = value
        K_λ[1, 2, 1, 1, 2, 1, 1, 1] = value
        K_λ[1, 2, 1, 1, 2, 2, 1, 1] = value
        K_λ[1, 2, 1, 2, 2, 1, 1, 2] = value
        K_λ[1, 2, 1, 2, 2, 2, 1, 2] = value
        K_λ[1, 2, 2, 1, 2, 1, 1, 1] = value
        K_λ[1, 2, 2, 1, 2, 2, 1, 1] = value
        K_λ[1, 2, 2, 2, 2, 1, 1, 2] = value
        K_λ[1, 2, 2, 2, 2, 2, 1, 2] = value
        K_λ[2, 1, 1, 1, 1, 2, 1, 1] = value
        K_λ[2, 1, 1, 1, 1, 2, 2, 1] = value
        K_λ[2, 1, 1, 2, 1, 2, 1, 2] = value
        K_λ[2, 1, 1, 2, 1, 2, 2, 2] = value
        K_λ[2, 1, 2, 1, 1, 1, 1, 1] = value
        K_λ[2, 1, 2, 1, 1, 1, 2, 1] = value
        K_λ[2, 1, 2, 2, 1, 1, 1, 2] = value
        K_λ[2, 1, 2, 2, 1, 1, 2, 2] = value
        K_λ[2, 2, 1, 1, 1, 2, 1, 1] = value
        K_λ[2, 2, 1, 1, 1, 2, 2, 1] = value
        K_λ[2, 2, 1, 2, 1, 2, 1, 2] = value
        K_λ[2, 2, 1, 2, 1, 2, 2, 2] = value
        K_λ[2, 2, 2, 1, 1, 1, 1, 1] = value
        K_λ[2, 2, 2, 1, 1, 1, 2, 1] = value
        K_λ[2, 2, 2, 2, 1, 1, 1, 2] = value
        K_λ[2, 2, 2, 2, 1, 1, 2, 2] = value

        value = -h[3] / 24
        K_λ[1, 1, 1, 1, 2, 1, 2, 2] = value
        K_λ[1, 1, 1, 1, 2, 2, 2, 2] = value
        K_λ[1, 1, 1, 2, 2, 1, 2, 1] = value
        K_λ[1, 1, 1, 2, 2, 2, 2, 1] = value
        K_λ[1, 1, 2, 1, 2, 1, 2, 2] = value
        K_λ[1, 1, 2, 1, 2, 2, 2, 2] = value
        K_λ[1, 1, 2, 2, 2, 1, 2, 1] = value
        K_λ[1, 1, 2, 2, 2, 2, 2, 1] = value
        K_λ[1, 2, 1, 1, 2, 1, 1, 2] = value
        K_λ[1, 2, 1, 1, 2, 2, 1, 2] = value
        K_λ[1, 2, 1, 2, 2, 1, 1, 1] = value
        K_λ[1, 2, 1, 2, 2, 2, 1, 1] = value
        K_λ[1, 2, 2, 1, 2, 1, 1, 2] = value
        K_λ[1, 2, 2, 1, 2, 2, 1, 2] = value
        K_λ[1, 2, 2, 2, 2, 1, 1, 1] = value
        K_λ[1, 2, 2, 2, 2, 2, 1, 1] = value
        K_λ[2, 1, 1, 1, 1, 2, 1, 2] = value
        K_λ[2, 1, 1, 1, 1, 2, 2, 2] = value
        K_λ[2, 1, 1, 2, 1, 2, 1, 1] = value
        K_λ[2, 1, 1, 2, 1, 2, 2, 1] = value
        K_λ[2, 1, 2, 1, 1, 1, 1, 2] = value
        K_λ[2, 1, 2, 1, 1, 1, 2, 2] = value
        K_λ[2, 1, 2, 2, 1, 1, 1, 1] = value
        K_λ[2, 1, 2, 2, 1, 1, 2, 1] = value
        K_λ[2, 2, 1, 1, 1, 2, 1, 2] = value
        K_λ[2, 2, 1, 1, 1, 2, 2, 2] = value
        K_λ[2, 2, 1, 2, 1, 2, 1, 1] = value
        K_λ[2, 2, 1, 2, 1, 2, 2, 1] = value
        K_λ[2, 2, 2, 1, 1, 1, 1, 2] = value
        K_λ[2, 2, 2, 1, 1, 1, 2, 2] = value
        K_λ[2, 2, 2, 2, 1, 1, 1, 1] = value
        K_λ[2, 2, 2, 2, 1, 1, 2, 1] = value

        value = h[3] / 24
        K_λ[1, 1, 1, 1, 2, 1, 1, 2] = value
        K_λ[1, 1, 1, 1, 2, 2, 1, 2] = value
        K_λ[1, 1, 1, 2, 2, 1, 1, 1] = value
        K_λ[1, 1, 1, 2, 2, 2, 1, 1] = value
        K_λ[1, 1, 2, 1, 2, 1, 1, 2] = value
        K_λ[1, 1, 2, 1, 2, 2, 1, 2] = value
        K_λ[1, 1, 2, 2, 2, 1, 1, 1] = value
        K_λ[1, 1, 2, 2, 2, 2, 1, 1] = value
        K_λ[1, 2, 1, 1, 2, 1, 2, 2] = value
        K_λ[1, 2, 1, 1, 2, 2, 2, 2] = value
        K_λ[1, 2, 1, 2, 2, 1, 2, 1] = value
        K_λ[1, 2, 1, 2, 2, 2, 2, 1] = value
        K_λ[1, 2, 2, 1, 2, 1, 2, 2] = value
        K_λ[1, 2, 2, 1, 2, 2, 2, 2] = value
        K_λ[1, 2, 2, 2, 2, 1, 2, 1] = value
        K_λ[1, 2, 2, 2, 2, 2, 2, 1] = value
        K_λ[2, 1, 1, 1, 1, 1, 1, 2] = value
        K_λ[2, 1, 1, 1, 1, 1, 2, 2] = value
        K_λ[2, 1, 1, 2, 1, 1, 1, 1] = value
        K_λ[2, 1, 1, 2, 1, 1, 2, 1] = value
        K_λ[2, 1, 2, 1, 1, 2, 1, 2] = value
        K_λ[2, 1, 2, 1, 1, 2, 2, 2] = value
        K_λ[2, 1, 2, 2, 1, 2, 1, 1] = value
        K_λ[2, 1, 2, 2, 1, 2, 2, 1] = value
        K_λ[2, 2, 1, 1, 1, 1, 1, 2] = value
        K_λ[2, 2, 1, 1, 1, 1, 2, 2] = value
        K_λ[2, 2, 1, 2, 1, 1, 1, 1] = value
        K_λ[2, 2, 1, 2, 1, 1, 2, 1] = value
        K_λ[2, 2, 2, 1, 1, 2, 1, 2] = value
        K_λ[2, 2, 2, 1, 1, 2, 2, 2] = value
        K_λ[2, 2, 2, 2, 1, 2, 1, 1] = value
        K_λ[2, 2, 2, 2, 1, 2, 2, 1] = value

        value = h[3] / 12
        K_λ[1, 1, 1, 1, 2, 1, 1, 1] = value
        K_λ[1, 1, 1, 1, 2, 2, 1, 1] = value
        K_λ[1, 1, 1, 2, 2, 1, 1, 2] = value
        K_λ[1, 1, 1, 2, 2, 2, 1, 2] = value
        K_λ[1, 1, 2, 1, 2, 1, 1, 1] = value
        K_λ[1, 1, 2, 1, 2, 2, 1, 1] = value
        K_λ[1, 1, 2, 2, 2, 1, 1, 2] = value
        K_λ[1, 1, 2, 2, 2, 2, 1, 2] = value
        K_λ[1, 2, 1, 1, 2, 1, 2, 1] = value
        K_λ[1, 2, 1, 1, 2, 2, 2, 1] = value
        K_λ[1, 2, 1, 2, 2, 1, 2, 2] = value
        K_λ[1, 2, 1, 2, 2, 2, 2, 2] = value
        K_λ[1, 2, 2, 1, 2, 1, 2, 1] = value
        K_λ[1, 2, 2, 1, 2, 2, 2, 1] = value
        K_λ[1, 2, 2, 2, 2, 1, 2, 2] = value
        K_λ[1, 2, 2, 2, 2, 2, 2, 2] = value
        K_λ[2, 1, 1, 1, 1, 1, 1, 1] = value
        K_λ[2, 1, 1, 1, 1, 1, 2, 1] = value
        K_λ[2, 1, 1, 2, 1, 1, 1, 2] = value
        K_λ[2, 1, 1, 2, 1, 1, 2, 2] = value
        K_λ[2, 1, 2, 1, 1, 2, 1, 1] = value
        K_λ[2, 1, 2, 1, 1, 2, 2, 1] = value
        K_λ[2, 1, 2, 2, 1, 2, 1, 2] = value
        K_λ[2, 1, 2, 2, 1, 2, 2, 2] = value
        K_λ[2, 2, 1, 1, 1, 1, 1, 1] = value
        K_λ[2, 2, 1, 1, 1, 1, 2, 1] = value
        K_λ[2, 2, 1, 2, 1, 1, 1, 2] = value
        K_λ[2, 2, 1, 2, 1, 1, 2, 2] = value
        K_λ[2, 2, 2, 1, 1, 2, 1, 1] = value
        K_λ[2, 2, 2, 1, 1, 2, 2, 1] = value
        K_λ[2, 2, 2, 2, 1, 2, 1, 2] = value
        K_λ[2, 2, 2, 2, 1, 2, 2, 2] = value

        value = -(h[1] * h[3]) / (9 * h[2])
        K_λ[2, 1, 1, 1, 2, 1, 2, 1] = value
        K_λ[2, 1, 1, 2, 2, 1, 2, 2] = value
        K_λ[2, 1, 2, 1, 2, 1, 1, 1] = value
        K_λ[2, 1, 2, 2, 2, 1, 1, 2] = value
        K_λ[2, 2, 1, 1, 2, 2, 2, 1] = value
        K_λ[2, 2, 1, 2, 2, 2, 2, 2] = value
        K_λ[2, 2, 2, 1, 2, 2, 1, 1] = value
        K_λ[2, 2, 2, 2, 2, 2, 1, 2] = value

        value = -(h[1] * h[3]) / (18 * h[2])
        K_λ[2, 1, 1, 1, 2, 1, 2, 2] = value
        K_λ[2, 1, 1, 1, 2, 2, 2, 1] = value
        K_λ[2, 1, 1, 2, 2, 1, 2, 1] = value
        K_λ[2, 1, 1, 2, 2, 2, 2, 2] = value
        K_λ[2, 1, 2, 1, 2, 1, 1, 2] = value
        K_λ[2, 1, 2, 1, 2, 2, 1, 1] = value
        K_λ[2, 1, 2, 2, 2, 1, 1, 1] = value
        K_λ[2, 1, 2, 2, 2, 2, 1, 2] = value
        K_λ[2, 2, 1, 1, 2, 1, 2, 1] = value
        K_λ[2, 2, 1, 1, 2, 2, 2, 2] = value
        K_λ[2, 2, 1, 2, 2, 1, 2, 2] = value
        K_λ[2, 2, 1, 2, 2, 2, 2, 1] = value
        K_λ[2, 2, 2, 1, 2, 1, 1, 1] = value
        K_λ[2, 2, 2, 1, 2, 2, 1, 2] = value
        K_λ[2, 2, 2, 2, 2, 1, 1, 2] = value
        K_λ[2, 2, 2, 2, 2, 2, 1, 1] = value

        value = -(h[1] * h[3]) / (36 * h[2])
        K_λ[2, 1, 1, 1, 2, 2, 2, 2] = value
        K_λ[2, 1, 1, 2, 2, 2, 2, 1] = value
        K_λ[2, 1, 2, 1, 2, 2, 1, 2] = value
        K_λ[2, 1, 2, 2, 2, 2, 1, 1] = value
        K_λ[2, 2, 1, 1, 2, 1, 2, 2] = value
        K_λ[2, 2, 1, 2, 2, 1, 2, 1] = value
        K_λ[2, 2, 2, 1, 2, 1, 1, 2] = value
        K_λ[2, 2, 2, 2, 2, 1, 1, 1] = value

        value = (h[1] * h[3]) / (36 * h[2])
        K_λ[2, 1, 1, 1, 2, 2, 1, 2] = value
        K_λ[2, 1, 1, 2, 2, 2, 1, 1] = value
        K_λ[2, 1, 2, 1, 2, 2, 2, 2] = value
        K_λ[2, 1, 2, 2, 2, 2, 2, 1] = value
        K_λ[2, 2, 1, 1, 2, 1, 1, 2] = value
        K_λ[2, 2, 1, 2, 2, 1, 1, 1] = value
        K_λ[2, 2, 2, 1, 2, 1, 2, 2] = value
        K_λ[2, 2, 2, 2, 2, 1, 2, 1] = value

        value = (h[1] * h[3]) / (18 * h[2])
        K_λ[2, 1, 1, 1, 2, 1, 1, 2] = value
        K_λ[2, 1, 1, 1, 2, 2, 1, 1] = value
        K_λ[2, 1, 1, 2, 2, 1, 1, 1] = value
        K_λ[2, 1, 1, 2, 2, 2, 1, 2] = value
        K_λ[2, 1, 2, 1, 2, 1, 2, 2] = value
        K_λ[2, 1, 2, 1, 2, 2, 2, 1] = value
        K_λ[2, 1, 2, 2, 2, 1, 2, 1] = value
        K_λ[2, 1, 2, 2, 2, 2, 2, 2] = value
        K_λ[2, 2, 1, 1, 2, 1, 1, 1] = value
        K_λ[2, 2, 1, 1, 2, 2, 1, 2] = value
        K_λ[2, 2, 1, 2, 2, 1, 1, 2] = value
        K_λ[2, 2, 1, 2, 2, 2, 1, 1] = value
        K_λ[2, 2, 2, 1, 2, 1, 2, 1] = value
        K_λ[2, 2, 2, 1, 2, 2, 2, 2] = value
        K_λ[2, 2, 2, 2, 2, 1, 2, 2] = value
        K_λ[2, 2, 2, 2, 2, 2, 2, 1] = value

        value = (h[1] * h[3]) / (9 * h[2])
        K_λ[2, 1, 1, 1, 2, 1, 1, 1] = value
        K_λ[2, 1, 1, 2, 2, 1, 1, 2] = value
        K_λ[2, 1, 2, 1, 2, 1, 2, 1] = value
        K_λ[2, 1, 2, 2, 2, 1, 2, 2] = value
        K_λ[2, 2, 1, 1, 2, 2, 1, 1] = value
        K_λ[2, 2, 1, 2, 2, 2, 1, 2] = value
        K_λ[2, 2, 2, 1, 2, 2, 2, 1] = value
        K_λ[2, 2, 2, 2, 2, 2, 2, 2] = value

        value = -(h[2] * h[3]) / (9 * h[1])
        K_λ[1, 1, 1, 1, 1, 2, 1, 1] = value
        K_λ[1, 1, 1, 2, 1, 2, 1, 2] = value
        K_λ[1, 1, 2, 1, 1, 2, 2, 1] = value
        K_λ[1, 1, 2, 2, 1, 2, 2, 2] = value
        K_λ[1, 2, 1, 1, 1, 1, 1, 1] = value
        K_λ[1, 2, 1, 2, 1, 1, 1, 2] = value
        K_λ[1, 2, 2, 1, 1, 1, 2, 1] = value
        K_λ[1, 2, 2, 2, 1, 1, 2, 2] = value

        value = -(h[2] * h[3]) / (18 * h[1])
        K_λ[1, 1, 1, 1, 1, 2, 1, 2] = value
        K_λ[1, 1, 1, 1, 1, 2, 2, 1] = value
        K_λ[1, 1, 1, 2, 1, 2, 1, 1] = value
        K_λ[1, 1, 1, 2, 1, 2, 2, 2] = value
        K_λ[1, 1, 2, 1, 1, 2, 1, 1] = value
        K_λ[1, 1, 2, 1, 1, 2, 2, 2] = value
        K_λ[1, 1, 2, 2, 1, 2, 1, 2] = value
        K_λ[1, 1, 2, 2, 1, 2, 2, 1] = value
        K_λ[1, 2, 1, 1, 1, 1, 1, 2] = value
        K_λ[1, 2, 1, 1, 1, 1, 2, 1] = value
        K_λ[1, 2, 1, 2, 1, 1, 1, 1] = value
        K_λ[1, 2, 1, 2, 1, 1, 2, 2] = value
        K_λ[1, 2, 2, 1, 1, 1, 1, 1] = value
        K_λ[1, 2, 2, 1, 1, 1, 2, 2] = value
        K_λ[1, 2, 2, 2, 1, 1, 1, 2] = value
        K_λ[1, 2, 2, 2, 1, 1, 2, 1] = value

        value = -(h[2] * h[3]) / (36 * h[1])
        K_λ[1, 1, 1, 1, 1, 2, 2, 2] = value
        K_λ[1, 1, 1, 2, 1, 2, 2, 1] = value
        K_λ[1, 1, 2, 1, 1, 2, 1, 2] = value
        K_λ[1, 1, 2, 2, 1, 2, 1, 1] = value
        K_λ[1, 2, 1, 1, 1, 1, 2, 2] = value
        K_λ[1, 2, 1, 2, 1, 1, 2, 1] = value
        K_λ[1, 2, 2, 1, 1, 1, 1, 2] = value
        K_λ[1, 2, 2, 2, 1, 1, 1, 1] = value

        value = (h[2] * h[3]) / (36 * h[1])
        K_λ[1, 1, 1, 1, 1, 1, 2, 2] = value
        K_λ[1, 1, 1, 2, 1, 1, 2, 1] = value
        K_λ[1, 1, 2, 1, 1, 1, 1, 2] = value
        K_λ[1, 1, 2, 2, 1, 1, 1, 1] = value
        K_λ[1, 2, 1, 1, 1, 2, 2, 2] = value
        K_λ[1, 2, 1, 2, 1, 2, 2, 1] = value
        K_λ[1, 2, 2, 1, 1, 2, 1, 2] = value
        K_λ[1, 2, 2, 2, 1, 2, 1, 1] = value

        value = (h[2] * h[3]) / (18 * h[1])
        K_λ[1, 1, 1, 1, 1, 1, 1, 2] = value
        K_λ[1, 1, 1, 1, 1, 1, 2, 1] = value
        K_λ[1, 1, 1, 2, 1, 1, 1, 1] = value
        K_λ[1, 1, 1, 2, 1, 1, 2, 2] = value
        K_λ[1, 1, 2, 1, 1, 1, 1, 1] = value
        K_λ[1, 1, 2, 1, 1, 1, 2, 2] = value
        K_λ[1, 1, 2, 2, 1, 1, 1, 2] = value
        K_λ[1, 1, 2, 2, 1, 1, 2, 1] = value
        K_λ[1, 2, 1, 1, 1, 2, 1, 2] = value
        K_λ[1, 2, 1, 1, 1, 2, 2, 1] = value
        K_λ[1, 2, 1, 2, 1, 2, 1, 1] = value
        K_λ[1, 2, 1, 2, 1, 2, 2, 2] = value
        K_λ[1, 2, 2, 1, 1, 2, 1, 1] = value
        K_λ[1, 2, 2, 1, 1, 2, 2, 2] = value
        K_λ[1, 2, 2, 2, 1, 2, 1, 2] = value
        K_λ[1, 2, 2, 2, 1, 2, 2, 1] = value

        value = (h[2] * h[3]) / (9 * h[1])
        K_λ[1, 1, 1, 1, 1, 1, 1, 1] = value
        K_λ[1, 1, 1, 2, 1, 1, 1, 2] = value
        K_λ[1, 1, 2, 1, 1, 1, 2, 1] = value
        K_λ[1, 1, 2, 2, 1, 1, 2, 2] = value
        K_λ[1, 2, 1, 1, 1, 2, 1, 1] = value
        K_λ[1, 2, 1, 2, 1, 2, 1, 2] = value
        K_λ[1, 2, 2, 1, 1, 2, 2, 1] = value
        K_λ[1, 2, 2, 2, 1, 2, 2, 2] = value
    else
    end
    return K_λ
end

function expected_K_μ(h::NTuple{d,T}) where {d,T<:Number}
    K_μ = zeros(Float64, d, fill(2, d)..., d, fill(2, d)...)
    if d == 2
        value = -1 / 4
        K_μ[1, 1, 1, 2, 2, 1] = value
        K_μ[1, 1, 1, 2, 2, 2] = value
        K_μ[1, 1, 2, 2, 1, 1] = value
        K_μ[1, 1, 2, 2, 1, 2] = value
        K_μ[1, 2, 1, 2, 2, 1] = value
        K_μ[1, 2, 1, 2, 2, 2] = value
        K_μ[1, 2, 2, 2, 1, 1] = value
        K_μ[1, 2, 2, 2, 1, 2] = value
        K_μ[2, 1, 1, 1, 1, 2] = value
        K_μ[2, 1, 1, 1, 2, 2] = value
        K_μ[2, 1, 2, 1, 1, 2] = value
        K_μ[2, 1, 2, 1, 2, 2] = value
        K_μ[2, 2, 1, 1, 1, 1] = value
        K_μ[2, 2, 1, 1, 2, 1] = value
        K_μ[2, 2, 2, 1, 1, 1] = value
        K_μ[2, 2, 2, 1, 2, 1] = value

        value = 1 / 4
        K_μ[1, 1, 1, 2, 1, 1] = value
        K_μ[1, 1, 1, 2, 1, 2] = value
        K_μ[1, 1, 2, 2, 2, 1] = value
        K_μ[1, 1, 2, 2, 2, 2] = value
        K_μ[1, 2, 1, 2, 1, 1] = value
        K_μ[1, 2, 1, 2, 1, 2] = value
        K_μ[1, 2, 2, 2, 2, 1] = value
        K_μ[1, 2, 2, 2, 2, 2] = value
        K_μ[2, 1, 1, 1, 1, 1] = value
        K_μ[2, 1, 1, 1, 2, 1] = value
        K_μ[2, 1, 2, 1, 1, 1] = value
        K_μ[2, 1, 2, 1, 2, 1] = value
        K_μ[2, 2, 1, 1, 1, 2] = value
        K_μ[2, 2, 1, 1, 2, 2] = value
        K_μ[2, 2, 2, 1, 1, 2] = value
        K_μ[2, 2, 2, 1, 2, 2] = value

        value = h[1] / (6 * h[2]) - (2 * h[2]) / (3 * h[1])
        K_μ[1, 1, 1, 1, 2, 1] = value
        K_μ[1, 1, 2, 1, 2, 2] = value
        K_μ[1, 2, 1, 1, 1, 1] = value
        K_μ[1, 2, 2, 1, 1, 2] = value

        value = (-h[2] / (3 * h[1])) - h[1] / (6 * h[2])
        K_μ[1, 1, 1, 1, 2, 2] = value
        K_μ[1, 1, 2, 1, 2, 1] = value
        K_μ[1, 2, 1, 1, 1, 2] = value
        K_μ[1, 2, 2, 1, 1, 1] = value

        value = h[1] / (3 * h[2]) - h[2] / (3 * h[1])
        K_μ[2, 1, 1, 2, 2, 1] = value
        K_μ[2, 1, 2, 2, 2, 2] = value
        K_μ[2, 2, 1, 2, 1, 1] = value
        K_μ[2, 2, 2, 2, 1, 2] = value

        value = (-h[2] / (6 * h[1])) - h[1] / (3 * h[2])
        K_μ[2, 1, 1, 2, 2, 2] = value
        K_μ[2, 1, 2, 2, 2, 1] = value
        K_μ[2, 2, 1, 2, 1, 2] = value
        K_μ[2, 2, 2, 2, 1, 1] = value

        value = h[2] / (6 * h[1]) - (2 * h[1]) / (3 * h[2])
        K_μ[2, 1, 1, 2, 1, 2] = value
        K_μ[2, 1, 2, 2, 1, 1] = value
        K_μ[2, 2, 1, 2, 2, 2] = value
        K_μ[2, 2, 2, 2, 2, 1] = value

        value = h[2] / (3 * h[1]) - h[1] / (3 * h[2])
        K_μ[1, 1, 1, 1, 1, 2] = value
        K_μ[1, 1, 2, 1, 1, 1] = value
        K_μ[1, 2, 1, 1, 2, 2] = value
        K_μ[1, 2, 2, 1, 2, 1] = value

        value = h[2] / (3 * h[1]) + (2 * h[1]) / (3 * h[2])
        K_μ[2, 1, 1, 2, 1, 1] = value
        K_μ[2, 1, 2, 2, 1, 2] = value
        K_μ[2, 2, 1, 2, 2, 1] = value
        K_μ[2, 2, 2, 2, 2, 2] = value

        value = (2 * h[2]) / (3 * h[1]) + h[1] / (3 * h[2])
        K_μ[1, 1, 1, 1, 1, 1] = value
        K_μ[1, 1, 2, 1, 1, 2] = value
        K_μ[1, 2, 1, 1, 2, 1] = value
        K_μ[1, 2, 2, 1, 2, 2] = value
    elseif d == 3
        value = -h[1] / 12
        K_μ[2, 1, 1, 1, 3, 1, 2, 1] = value
        K_μ[2, 1, 1, 1, 3, 1, 2, 2] = value
        K_μ[2, 1, 1, 2, 3, 1, 1, 1] = value
        K_μ[2, 1, 1, 2, 3, 1, 1, 2] = value
        K_μ[2, 1, 2, 1, 3, 1, 2, 1] = value
        K_μ[2, 1, 2, 1, 3, 1, 2, 2] = value
        K_μ[2, 1, 2, 2, 3, 1, 1, 1] = value
        K_μ[2, 1, 2, 2, 3, 1, 1, 2] = value
        K_μ[2, 2, 1, 1, 3, 2, 2, 1] = value
        K_μ[2, 2, 1, 1, 3, 2, 2, 2] = value
        K_μ[2, 2, 1, 2, 3, 2, 1, 1] = value
        K_μ[2, 2, 1, 2, 3, 2, 1, 2] = value
        K_μ[2, 2, 2, 1, 3, 2, 2, 1] = value
        K_μ[2, 2, 2, 1, 3, 2, 2, 2] = value
        K_μ[2, 2, 2, 2, 3, 2, 1, 1] = value
        K_μ[2, 2, 2, 2, 3, 2, 1, 2] = value
        K_μ[3, 1, 1, 1, 2, 1, 1, 2] = value
        K_μ[3, 1, 1, 1, 2, 1, 2, 2] = value
        K_μ[3, 1, 1, 2, 2, 1, 1, 2] = value
        K_μ[3, 1, 1, 2, 2, 1, 2, 2] = value
        K_μ[3, 1, 2, 1, 2, 1, 1, 1] = value
        K_μ[3, 1, 2, 1, 2, 1, 2, 1] = value
        K_μ[3, 1, 2, 2, 2, 1, 1, 1] = value
        K_μ[3, 1, 2, 2, 2, 1, 2, 1] = value
        K_μ[3, 2, 1, 1, 2, 2, 1, 2] = value
        K_μ[3, 2, 1, 1, 2, 2, 2, 2] = value
        K_μ[3, 2, 1, 2, 2, 2, 1, 2] = value
        K_μ[3, 2, 1, 2, 2, 2, 2, 2] = value
        K_μ[3, 2, 2, 1, 2, 2, 1, 1] = value
        K_μ[3, 2, 2, 1, 2, 2, 2, 1] = value
        K_μ[3, 2, 2, 2, 2, 2, 1, 1] = value
        K_μ[3, 2, 2, 2, 2, 2, 2, 1] = value

        value = -h[1] / 24
        K_μ[2, 1, 1, 1, 3, 2, 2, 1] = value
        K_μ[2, 1, 1, 1, 3, 2, 2, 2] = value
        K_μ[2, 1, 1, 2, 3, 2, 1, 1] = value
        K_μ[2, 1, 1, 2, 3, 2, 1, 2] = value
        K_μ[2, 1, 2, 1, 3, 2, 2, 1] = value
        K_μ[2, 1, 2, 1, 3, 2, 2, 2] = value
        K_μ[2, 1, 2, 2, 3, 2, 1, 1] = value
        K_μ[2, 1, 2, 2, 3, 2, 1, 2] = value
        K_μ[2, 2, 1, 1, 3, 1, 2, 1] = value
        K_μ[2, 2, 1, 1, 3, 1, 2, 2] = value
        K_μ[2, 2, 1, 2, 3, 1, 1, 1] = value
        K_μ[2, 2, 1, 2, 3, 1, 1, 2] = value
        K_μ[2, 2, 2, 1, 3, 1, 2, 1] = value
        K_μ[2, 2, 2, 1, 3, 1, 2, 2] = value
        K_μ[2, 2, 2, 2, 3, 1, 1, 1] = value
        K_μ[2, 2, 2, 2, 3, 1, 1, 2] = value
        K_μ[3, 1, 1, 1, 2, 2, 1, 2] = value
        K_μ[3, 1, 1, 1, 2, 2, 2, 2] = value
        K_μ[3, 1, 1, 2, 2, 2, 1, 2] = value
        K_μ[3, 1, 1, 2, 2, 2, 2, 2] = value
        K_μ[3, 1, 2, 1, 2, 2, 1, 1] = value
        K_μ[3, 1, 2, 1, 2, 2, 2, 1] = value
        K_μ[3, 1, 2, 2, 2, 2, 1, 1] = value
        K_μ[3, 1, 2, 2, 2, 2, 2, 1] = value
        K_μ[3, 2, 1, 1, 2, 1, 1, 2] = value
        K_μ[3, 2, 1, 1, 2, 1, 2, 2] = value
        K_μ[3, 2, 1, 2, 2, 1, 1, 2] = value
        K_μ[3, 2, 1, 2, 2, 1, 2, 2] = value
        K_μ[3, 2, 2, 1, 2, 1, 1, 1] = value
        K_μ[3, 2, 2, 1, 2, 1, 2, 1] = value
        K_μ[3, 2, 2, 2, 2, 1, 1, 1] = value
        K_μ[3, 2, 2, 2, 2, 1, 2, 1] = value

        value = h[1] / 24
        K_μ[2, 1, 1, 1, 3, 2, 1, 1] = value
        K_μ[2, 1, 1, 1, 3, 2, 1, 2] = value
        K_μ[2, 1, 1, 2, 3, 2, 2, 1] = value
        K_μ[2, 1, 1, 2, 3, 2, 2, 2] = value
        K_μ[2, 1, 2, 1, 3, 2, 1, 1] = value
        K_μ[2, 1, 2, 1, 3, 2, 1, 2] = value
        K_μ[2, 1, 2, 2, 3, 2, 2, 1] = value
        K_μ[2, 1, 2, 2, 3, 2, 2, 2] = value
        K_μ[2, 2, 1, 1, 3, 1, 1, 1] = value
        K_μ[2, 2, 1, 1, 3, 1, 1, 2] = value
        K_μ[2, 2, 1, 2, 3, 1, 2, 1] = value
        K_μ[2, 2, 1, 2, 3, 1, 2, 2] = value
        K_μ[2, 2, 2, 1, 3, 1, 1, 1] = value
        K_μ[2, 2, 2, 1, 3, 1, 1, 2] = value
        K_μ[2, 2, 2, 2, 3, 1, 2, 1] = value
        K_μ[2, 2, 2, 2, 3, 1, 2, 2] = value
        K_μ[3, 1, 1, 1, 2, 2, 1, 1] = value
        K_μ[3, 1, 1, 1, 2, 2, 2, 1] = value
        K_μ[3, 1, 1, 2, 2, 2, 1, 1] = value
        K_μ[3, 1, 1, 2, 2, 2, 2, 1] = value
        K_μ[3, 1, 2, 1, 2, 2, 1, 2] = value
        K_μ[3, 1, 2, 1, 2, 2, 2, 2] = value
        K_μ[3, 1, 2, 2, 2, 2, 1, 2] = value
        K_μ[3, 1, 2, 2, 2, 2, 2, 2] = value
        K_μ[3, 2, 1, 1, 2, 1, 1, 1] = value
        K_μ[3, 2, 1, 1, 2, 1, 2, 1] = value
        K_μ[3, 2, 1, 2, 2, 1, 1, 1] = value
        K_μ[3, 2, 1, 2, 2, 1, 2, 1] = value
        K_μ[3, 2, 2, 1, 2, 1, 1, 2] = value
        K_μ[3, 2, 2, 1, 2, 1, 2, 2] = value
        K_μ[3, 2, 2, 2, 2, 1, 1, 2] = value
        K_μ[3, 2, 2, 2, 2, 1, 2, 2] = value

        value = h[1] / 12
        K_μ[2, 1, 1, 1, 3, 1, 1, 1] = value
        K_μ[2, 1, 1, 1, 3, 1, 1, 2] = value
        K_μ[2, 1, 1, 2, 3, 1, 2, 1] = value
        K_μ[2, 1, 1, 2, 3, 1, 2, 2] = value
        K_μ[2, 1, 2, 1, 3, 1, 1, 1] = value
        K_μ[2, 1, 2, 1, 3, 1, 1, 2] = value
        K_μ[2, 1, 2, 2, 3, 1, 2, 1] = value
        K_μ[2, 1, 2, 2, 3, 1, 2, 2] = value
        K_μ[2, 2, 1, 1, 3, 2, 1, 1] = value
        K_μ[2, 2, 1, 1, 3, 2, 1, 2] = value
        K_μ[2, 2, 1, 2, 3, 2, 2, 1] = value
        K_μ[2, 2, 1, 2, 3, 2, 2, 2] = value
        K_μ[2, 2, 2, 1, 3, 2, 1, 1] = value
        K_μ[2, 2, 2, 1, 3, 2, 1, 2] = value
        K_μ[2, 2, 2, 2, 3, 2, 2, 1] = value
        K_μ[2, 2, 2, 2, 3, 2, 2, 2] = value
        K_μ[3, 1, 1, 1, 2, 1, 1, 1] = value
        K_μ[3, 1, 1, 1, 2, 1, 2, 1] = value
        K_μ[3, 1, 1, 2, 2, 1, 1, 1] = value
        K_μ[3, 1, 1, 2, 2, 1, 2, 1] = value
        K_μ[3, 1, 2, 1, 2, 1, 1, 2] = value
        K_μ[3, 1, 2, 1, 2, 1, 2, 2] = value
        K_μ[3, 1, 2, 2, 2, 1, 1, 2] = value
        K_μ[3, 1, 2, 2, 2, 1, 2, 2] = value
        K_μ[3, 2, 1, 1, 2, 2, 1, 1] = value
        K_μ[3, 2, 1, 1, 2, 2, 2, 1] = value
        K_μ[3, 2, 1, 2, 2, 2, 1, 1] = value
        K_μ[3, 2, 1, 2, 2, 2, 2, 1] = value
        K_μ[3, 2, 2, 1, 2, 2, 1, 2] = value
        K_μ[3, 2, 2, 1, 2, 2, 2, 2] = value
        K_μ[3, 2, 2, 2, 2, 2, 1, 2] = value
        K_μ[3, 2, 2, 2, 2, 2, 2, 2] = value

        value = -h[2] / 12
        K_μ[1, 1, 1, 1, 3, 2, 1, 1] = value
        K_μ[1, 1, 1, 1, 3, 2, 1, 2] = value
        K_μ[1, 1, 1, 2, 3, 1, 1, 1] = value
        K_μ[1, 1, 1, 2, 3, 1, 1, 2] = value
        K_μ[1, 1, 2, 1, 3, 2, 2, 1] = value
        K_μ[1, 1, 2, 1, 3, 2, 2, 2] = value
        K_μ[1, 1, 2, 2, 3, 1, 2, 1] = value
        K_μ[1, 1, 2, 2, 3, 1, 2, 2] = value
        K_μ[1, 2, 1, 1, 3, 2, 1, 1] = value
        K_μ[1, 2, 1, 1, 3, 2, 1, 2] = value
        K_μ[1, 2, 1, 2, 3, 1, 1, 1] = value
        K_μ[1, 2, 1, 2, 3, 1, 1, 2] = value
        K_μ[1, 2, 2, 1, 3, 2, 2, 1] = value
        K_μ[1, 2, 2, 1, 3, 2, 2, 2] = value
        K_μ[1, 2, 2, 2, 3, 1, 2, 1] = value
        K_μ[1, 2, 2, 2, 3, 1, 2, 2] = value
        K_μ[3, 1, 1, 1, 1, 1, 1, 2] = value
        K_μ[3, 1, 1, 1, 1, 2, 1, 2] = value
        K_μ[3, 1, 1, 2, 1, 1, 1, 2] = value
        K_μ[3, 1, 1, 2, 1, 2, 1, 2] = value
        K_μ[3, 1, 2, 1, 1, 1, 2, 2] = value
        K_μ[3, 1, 2, 1, 1, 2, 2, 2] = value
        K_μ[3, 1, 2, 2, 1, 1, 2, 2] = value
        K_μ[3, 1, 2, 2, 1, 2, 2, 2] = value
        K_μ[3, 2, 1, 1, 1, 1, 1, 1] = value
        K_μ[3, 2, 1, 1, 1, 2, 1, 1] = value
        K_μ[3, 2, 1, 2, 1, 1, 1, 1] = value
        K_μ[3, 2, 1, 2, 1, 2, 1, 1] = value
        K_μ[3, 2, 2, 1, 1, 1, 2, 1] = value
        K_μ[3, 2, 2, 1, 1, 2, 2, 1] = value
        K_μ[3, 2, 2, 2, 1, 1, 2, 1] = value
        K_μ[3, 2, 2, 2, 1, 2, 2, 1] = value

        value = -h[2] / 24
        K_μ[1, 1, 1, 1, 3, 2, 2, 1] = value
        K_μ[1, 1, 1, 1, 3, 2, 2, 2] = value
        K_μ[1, 1, 1, 2, 3, 1, 2, 1] = value
        K_μ[1, 1, 1, 2, 3, 1, 2, 2] = value
        K_μ[1, 1, 2, 1, 3, 2, 1, 1] = value
        K_μ[1, 1, 2, 1, 3, 2, 1, 2] = value
        K_μ[1, 1, 2, 2, 3, 1, 1, 1] = value
        K_μ[1, 1, 2, 2, 3, 1, 1, 2] = value
        K_μ[1, 2, 1, 1, 3, 2, 2, 1] = value
        K_μ[1, 2, 1, 1, 3, 2, 2, 2] = value
        K_μ[1, 2, 1, 2, 3, 1, 2, 1] = value
        K_μ[1, 2, 1, 2, 3, 1, 2, 2] = value
        K_μ[1, 2, 2, 1, 3, 2, 1, 1] = value
        K_μ[1, 2, 2, 1, 3, 2, 1, 2] = value
        K_μ[1, 2, 2, 2, 3, 1, 1, 1] = value
        K_μ[1, 2, 2, 2, 3, 1, 1, 2] = value
        K_μ[3, 1, 1, 1, 1, 1, 2, 2] = value
        K_μ[3, 1, 1, 1, 1, 2, 2, 2] = value
        K_μ[3, 1, 1, 2, 1, 1, 2, 2] = value
        K_μ[3, 1, 1, 2, 1, 2, 2, 2] = value
        K_μ[3, 1, 2, 1, 1, 1, 1, 2] = value
        K_μ[3, 1, 2, 1, 1, 2, 1, 2] = value
        K_μ[3, 1, 2, 2, 1, 1, 1, 2] = value
        K_μ[3, 1, 2, 2, 1, 2, 1, 2] = value
        K_μ[3, 2, 1, 1, 1, 1, 2, 1] = value
        K_μ[3, 2, 1, 1, 1, 2, 2, 1] = value
        K_μ[3, 2, 1, 2, 1, 1, 2, 1] = value
        K_μ[3, 2, 1, 2, 1, 2, 2, 1] = value
        K_μ[3, 2, 2, 1, 1, 1, 1, 1] = value
        K_μ[3, 2, 2, 1, 1, 2, 1, 1] = value
        K_μ[3, 2, 2, 2, 1, 1, 1, 1] = value
        K_μ[3, 2, 2, 2, 1, 2, 1, 1] = value

        value = h[2] / 24
        K_μ[1, 1, 1, 1, 3, 1, 2, 1] = value
        K_μ[1, 1, 1, 1, 3, 1, 2, 2] = value
        K_μ[1, 1, 1, 2, 3, 2, 2, 1] = value
        K_μ[1, 1, 1, 2, 3, 2, 2, 2] = value
        K_μ[1, 1, 2, 1, 3, 1, 1, 1] = value
        K_μ[1, 1, 2, 1, 3, 1, 1, 2] = value
        K_μ[1, 1, 2, 2, 3, 2, 1, 1] = value
        K_μ[1, 1, 2, 2, 3, 2, 1, 2] = value
        K_μ[1, 2, 1, 1, 3, 1, 2, 1] = value
        K_μ[1, 2, 1, 1, 3, 1, 2, 2] = value
        K_μ[1, 2, 1, 2, 3, 2, 2, 1] = value
        K_μ[1, 2, 1, 2, 3, 2, 2, 2] = value
        K_μ[1, 2, 2, 1, 3, 1, 1, 1] = value
        K_μ[1, 2, 2, 1, 3, 1, 1, 2] = value
        K_μ[1, 2, 2, 2, 3, 2, 1, 1] = value
        K_μ[1, 2, 2, 2, 3, 2, 1, 2] = value
        K_μ[3, 1, 1, 1, 1, 1, 2, 1] = value
        K_μ[3, 1, 1, 1, 1, 2, 2, 1] = value
        K_μ[3, 1, 1, 2, 1, 1, 2, 1] = value
        K_μ[3, 1, 1, 2, 1, 2, 2, 1] = value
        K_μ[3, 1, 2, 1, 1, 1, 1, 1] = value
        K_μ[3, 1, 2, 1, 1, 2, 1, 1] = value
        K_μ[3, 1, 2, 2, 1, 1, 1, 1] = value
        K_μ[3, 1, 2, 2, 1, 2, 1, 1] = value
        K_μ[3, 2, 1, 1, 1, 1, 2, 2] = value
        K_μ[3, 2, 1, 1, 1, 2, 2, 2] = value
        K_μ[3, 2, 1, 2, 1, 1, 2, 2] = value
        K_μ[3, 2, 1, 2, 1, 2, 2, 2] = value
        K_μ[3, 2, 2, 1, 1, 1, 1, 2] = value
        K_μ[3, 2, 2, 1, 1, 2, 1, 2] = value
        K_μ[3, 2, 2, 2, 1, 1, 1, 2] = value
        K_μ[3, 2, 2, 2, 1, 2, 1, 2] = value

        value = h[2] / 12
        K_μ[1, 1, 1, 1, 3, 1, 1, 1] = value
        K_μ[1, 1, 1, 1, 3, 1, 1, 2] = value
        K_μ[1, 1, 1, 2, 3, 2, 1, 1] = value
        K_μ[1, 1, 1, 2, 3, 2, 1, 2] = value
        K_μ[1, 1, 2, 1, 3, 1, 2, 1] = value
        K_μ[1, 1, 2, 1, 3, 1, 2, 2] = value
        K_μ[1, 1, 2, 2, 3, 2, 2, 1] = value
        K_μ[1, 1, 2, 2, 3, 2, 2, 2] = value
        K_μ[1, 2, 1, 1, 3, 1, 1, 1] = value
        K_μ[1, 2, 1, 1, 3, 1, 1, 2] = value
        K_μ[1, 2, 1, 2, 3, 2, 1, 1] = value
        K_μ[1, 2, 1, 2, 3, 2, 1, 2] = value
        K_μ[1, 2, 2, 1, 3, 1, 2, 1] = value
        K_μ[1, 2, 2, 1, 3, 1, 2, 2] = value
        K_μ[1, 2, 2, 2, 3, 2, 2, 1] = value
        K_μ[1, 2, 2, 2, 3, 2, 2, 2] = value
        K_μ[3, 1, 1, 1, 1, 1, 1, 1] = value
        K_μ[3, 1, 1, 1, 1, 2, 1, 1] = value
        K_μ[3, 1, 1, 2, 1, 1, 1, 1] = value
        K_μ[3, 1, 1, 2, 1, 2, 1, 1] = value
        K_μ[3, 1, 2, 1, 1, 1, 2, 1] = value
        K_μ[3, 1, 2, 1, 1, 2, 2, 1] = value
        K_μ[3, 1, 2, 2, 1, 1, 2, 1] = value
        K_μ[3, 1, 2, 2, 1, 2, 2, 1] = value
        K_μ[3, 2, 1, 1, 1, 1, 1, 2] = value
        K_μ[3, 2, 1, 1, 1, 2, 1, 2] = value
        K_μ[3, 2, 1, 2, 1, 1, 1, 2] = value
        K_μ[3, 2, 1, 2, 1, 2, 1, 2] = value
        K_μ[3, 2, 2, 1, 1, 1, 2, 2] = value
        K_μ[3, 2, 2, 1, 1, 2, 2, 2] = value
        K_μ[3, 2, 2, 2, 1, 1, 2, 2] = value
        K_μ[3, 2, 2, 2, 1, 2, 2, 2] = value

        value = -h[3] / 12
        K_μ[1, 1, 1, 1, 2, 2, 1, 1] = value
        K_μ[1, 1, 1, 1, 2, 2, 2, 1] = value
        K_μ[1, 1, 1, 2, 2, 2, 1, 2] = value
        K_μ[1, 1, 1, 2, 2, 2, 2, 2] = value
        K_μ[1, 1, 2, 1, 2, 1, 1, 1] = value
        K_μ[1, 1, 2, 1, 2, 1, 2, 1] = value
        K_μ[1, 1, 2, 2, 2, 1, 1, 2] = value
        K_μ[1, 1, 2, 2, 2, 1, 2, 2] = value
        K_μ[1, 2, 1, 1, 2, 2, 1, 1] = value
        K_μ[1, 2, 1, 1, 2, 2, 2, 1] = value
        K_μ[1, 2, 1, 2, 2, 2, 1, 2] = value
        K_μ[1, 2, 1, 2, 2, 2, 2, 2] = value
        K_μ[1, 2, 2, 1, 2, 1, 1, 1] = value
        K_μ[1, 2, 2, 1, 2, 1, 2, 1] = value
        K_μ[1, 2, 2, 2, 2, 1, 1, 2] = value
        K_μ[1, 2, 2, 2, 2, 1, 2, 2] = value
        K_μ[2, 1, 1, 1, 1, 1, 2, 1] = value
        K_μ[2, 1, 1, 1, 1, 2, 2, 1] = value
        K_μ[2, 1, 1, 2, 1, 1, 2, 2] = value
        K_μ[2, 1, 1, 2, 1, 2, 2, 2] = value
        K_μ[2, 1, 2, 1, 1, 1, 2, 1] = value
        K_μ[2, 1, 2, 1, 1, 2, 2, 1] = value
        K_μ[2, 1, 2, 2, 1, 1, 2, 2] = value
        K_μ[2, 1, 2, 2, 1, 2, 2, 2] = value
        K_μ[2, 2, 1, 1, 1, 1, 1, 1] = value
        K_μ[2, 2, 1, 1, 1, 2, 1, 1] = value
        K_μ[2, 2, 1, 2, 1, 1, 1, 2] = value
        K_μ[2, 2, 1, 2, 1, 2, 1, 2] = value
        K_μ[2, 2, 2, 1, 1, 1, 1, 1] = value
        K_μ[2, 2, 2, 1, 1, 2, 1, 1] = value
        K_μ[2, 2, 2, 2, 1, 1, 1, 2] = value
        K_μ[2, 2, 2, 2, 1, 2, 1, 2] = value

        value = -h[3] / 24
        K_μ[1, 1, 1, 1, 2, 2, 1, 2] = value
        K_μ[1, 1, 1, 1, 2, 2, 2, 2] = value
        K_μ[1, 1, 1, 2, 2, 2, 1, 1] = value
        K_μ[1, 1, 1, 2, 2, 2, 2, 1] = value
        K_μ[1, 1, 2, 1, 2, 1, 1, 2] = value
        K_μ[1, 1, 2, 1, 2, 1, 2, 2] = value
        K_μ[1, 1, 2, 2, 2, 1, 1, 1] = value
        K_μ[1, 1, 2, 2, 2, 1, 2, 1] = value
        K_μ[1, 2, 1, 1, 2, 2, 1, 2] = value
        K_μ[1, 2, 1, 1, 2, 2, 2, 2] = value
        K_μ[1, 2, 1, 2, 2, 2, 1, 1] = value
        K_μ[1, 2, 1, 2, 2, 2, 2, 1] = value
        K_μ[1, 2, 2, 1, 2, 1, 1, 2] = value
        K_μ[1, 2, 2, 1, 2, 1, 2, 2] = value
        K_μ[1, 2, 2, 2, 2, 1, 1, 1] = value
        K_μ[1, 2, 2, 2, 2, 1, 2, 1] = value
        K_μ[2, 1, 1, 1, 1, 1, 2, 2] = value
        K_μ[2, 1, 1, 1, 1, 2, 2, 2] = value
        K_μ[2, 1, 1, 2, 1, 1, 2, 1] = value
        K_μ[2, 1, 1, 2, 1, 2, 2, 1] = value
        K_μ[2, 1, 2, 1, 1, 1, 2, 2] = value
        K_μ[2, 1, 2, 1, 1, 2, 2, 2] = value
        K_μ[2, 1, 2, 2, 1, 1, 2, 1] = value
        K_μ[2, 1, 2, 2, 1, 2, 2, 1] = value
        K_μ[2, 2, 1, 1, 1, 1, 1, 2] = value
        K_μ[2, 2, 1, 1, 1, 2, 1, 2] = value
        K_μ[2, 2, 1, 2, 1, 1, 1, 1] = value
        K_μ[2, 2, 1, 2, 1, 2, 1, 1] = value
        K_μ[2, 2, 2, 1, 1, 1, 1, 2] = value
        K_μ[2, 2, 2, 1, 1, 2, 1, 2] = value
        K_μ[2, 2, 2, 2, 1, 1, 1, 1] = value
        K_μ[2, 2, 2, 2, 1, 2, 1, 1] = value

        value = h[3] / 24
        K_μ[1, 1, 1, 1, 2, 1, 1, 2] = value
        K_μ[1, 1, 1, 1, 2, 1, 2, 2] = value
        K_μ[1, 1, 1, 2, 2, 1, 1, 1] = value
        K_μ[1, 1, 1, 2, 2, 1, 2, 1] = value
        K_μ[1, 1, 2, 1, 2, 2, 1, 2] = value
        K_μ[1, 1, 2, 1, 2, 2, 2, 2] = value
        K_μ[1, 1, 2, 2, 2, 2, 1, 1] = value
        K_μ[1, 1, 2, 2, 2, 2, 2, 1] = value
        K_μ[1, 2, 1, 1, 2, 1, 1, 2] = value
        K_μ[1, 2, 1, 1, 2, 1, 2, 2] = value
        K_μ[1, 2, 1, 2, 2, 1, 1, 1] = value
        K_μ[1, 2, 1, 2, 2, 1, 2, 1] = value
        K_μ[1, 2, 2, 1, 2, 2, 1, 2] = value
        K_μ[1, 2, 2, 1, 2, 2, 2, 2] = value
        K_μ[1, 2, 2, 2, 2, 2, 1, 1] = value
        K_μ[1, 2, 2, 2, 2, 2, 2, 1] = value
        K_μ[2, 1, 1, 1, 1, 1, 1, 2] = value
        K_μ[2, 1, 1, 1, 1, 2, 1, 2] = value
        K_μ[2, 1, 1, 2, 1, 1, 1, 1] = value
        K_μ[2, 1, 1, 2, 1, 2, 1, 1] = value
        K_μ[2, 1, 2, 1, 1, 1, 1, 2] = value
        K_μ[2, 1, 2, 1, 1, 2, 1, 2] = value
        K_μ[2, 1, 2, 2, 1, 1, 1, 1] = value
        K_μ[2, 1, 2, 2, 1, 2, 1, 1] = value
        K_μ[2, 2, 1, 1, 1, 1, 2, 2] = value
        K_μ[2, 2, 1, 1, 1, 2, 2, 2] = value
        K_μ[2, 2, 1, 2, 1, 1, 2, 1] = value
        K_μ[2, 2, 1, 2, 1, 2, 2, 1] = value
        K_μ[2, 2, 2, 1, 1, 1, 2, 2] = value
        K_μ[2, 2, 2, 1, 1, 2, 2, 2] = value
        K_μ[2, 2, 2, 2, 1, 1, 2, 1] = value
        K_μ[2, 2, 2, 2, 1, 2, 2, 1] = value

        value = h[3] / 12
        K_μ[1, 1, 1, 1, 2, 1, 1, 1] = value
        K_μ[1, 1, 1, 1, 2, 1, 2, 1] = value
        K_μ[1, 1, 1, 2, 2, 1, 1, 2] = value
        K_μ[1, 1, 1, 2, 2, 1, 2, 2] = value
        K_μ[1, 1, 2, 1, 2, 2, 1, 1] = value
        K_μ[1, 1, 2, 1, 2, 2, 2, 1] = value
        K_μ[1, 1, 2, 2, 2, 2, 1, 2] = value
        K_μ[1, 1, 2, 2, 2, 2, 2, 2] = value
        K_μ[1, 2, 1, 1, 2, 1, 1, 1] = value
        K_μ[1, 2, 1, 1, 2, 1, 2, 1] = value
        K_μ[1, 2, 1, 2, 2, 1, 1, 2] = value
        K_μ[1, 2, 1, 2, 2, 1, 2, 2] = value
        K_μ[1, 2, 2, 1, 2, 2, 1, 1] = value
        K_μ[1, 2, 2, 1, 2, 2, 2, 1] = value
        K_μ[1, 2, 2, 2, 2, 2, 1, 2] = value
        K_μ[1, 2, 2, 2, 2, 2, 2, 2] = value
        K_μ[2, 1, 1, 1, 1, 1, 1, 1] = value
        K_μ[2, 1, 1, 1, 1, 2, 1, 1] = value
        K_μ[2, 1, 1, 2, 1, 1, 1, 2] = value
        K_μ[2, 1, 1, 2, 1, 2, 1, 2] = value
        K_μ[2, 1, 2, 1, 1, 1, 1, 1] = value
        K_μ[2, 1, 2, 1, 1, 2, 1, 1] = value
        K_μ[2, 1, 2, 2, 1, 1, 1, 2] = value
        K_μ[2, 1, 2, 2, 1, 2, 1, 2] = value
        K_μ[2, 2, 1, 1, 1, 1, 2, 1] = value
        K_μ[2, 2, 1, 1, 1, 2, 2, 1] = value
        K_μ[2, 2, 1, 2, 1, 1, 2, 2] = value
        K_μ[2, 2, 1, 2, 1, 2, 2, 2] = value
        K_μ[2, 2, 2, 1, 1, 1, 2, 1] = value
        K_μ[2, 2, 2, 1, 1, 2, 2, 1] = value
        K_μ[2, 2, 2, 2, 1, 1, 2, 2] = value
        K_μ[2, 2, 2, 2, 1, 2, 2, 2] = value

        value =
            (-(2 * h[2] * h[3]) / (9 * h[1])) +
            (h[1] * h[3]) / (18 * h[2]) +
            (h[1] * h[2]) / (18 * h[3])
        K_μ[1, 1, 1, 1, 1, 2, 1, 1] = value
        K_μ[1, 1, 1, 2, 1, 2, 1, 2] = value
        K_μ[1, 1, 2, 1, 1, 2, 2, 1] = value
        K_μ[1, 1, 2, 2, 1, 2, 2, 2] = value
        K_μ[1, 2, 1, 1, 1, 1, 1, 1] = value
        K_μ[1, 2, 1, 2, 1, 1, 1, 2] = value
        K_μ[1, 2, 2, 1, 1, 1, 2, 1] = value
        K_μ[1, 2, 2, 2, 1, 1, 2, 2] = value

        value =
            (-(h[2] * h[3]) / (9 * h[1])) - (h[1] * h[3]) / (18 * h[2]) +
            (h[1] * h[2]) / (36 * h[3])
        K_μ[1, 1, 1, 1, 1, 2, 2, 1] = value
        K_μ[1, 1, 1, 2, 1, 2, 2, 2] = value
        K_μ[1, 1, 2, 1, 1, 2, 1, 1] = value
        K_μ[1, 1, 2, 2, 1, 2, 1, 2] = value
        K_μ[1, 2, 1, 1, 1, 1, 2, 1] = value
        K_μ[1, 2, 1, 2, 1, 1, 2, 2] = value
        K_μ[1, 2, 2, 1, 1, 1, 1, 1] = value
        K_μ[1, 2, 2, 2, 1, 1, 1, 2] = value

        value =
            (-(h[2] * h[3]) / (9 * h[1])) + (h[1] * h[3]) / (36 * h[2]) -
            (h[1] * h[2]) / (18 * h[3])
        K_μ[1, 1, 1, 1, 1, 2, 1, 2] = value
        K_μ[1, 1, 1, 2, 1, 2, 1, 1] = value
        K_μ[1, 1, 2, 1, 1, 2, 2, 2] = value
        K_μ[1, 1, 2, 2, 1, 2, 2, 1] = value
        K_μ[1, 2, 1, 1, 1, 1, 1, 2] = value
        K_μ[1, 2, 1, 2, 1, 1, 1, 1] = value
        K_μ[1, 2, 2, 1, 1, 1, 2, 2] = value
        K_μ[1, 2, 2, 2, 1, 1, 2, 1] = value

        value =
            (-(h[2] * h[3]) / (9 * h[1])) +
            (h[1] * h[3]) / (18 * h[2]) +
            (h[1] * h[2]) / (9 * h[3])
        K_μ[3, 1, 1, 1, 3, 2, 1, 1] = value
        K_μ[3, 1, 1, 2, 3, 2, 1, 2] = value
        K_μ[3, 1, 2, 1, 3, 2, 2, 1] = value
        K_μ[3, 1, 2, 2, 3, 2, 2, 2] = value
        K_μ[3, 2, 1, 1, 3, 1, 1, 1] = value
        K_μ[3, 2, 1, 2, 3, 1, 1, 2] = value
        K_μ[3, 2, 2, 1, 3, 1, 2, 1] = value
        K_μ[3, 2, 2, 2, 3, 1, 2, 2] = value

        value =
            (-(h[2] * h[3]) / (9 * h[1])) +
            (h[1] * h[3]) / (9 * h[2]) +
            (h[1] * h[2]) / (18 * h[3])
        K_μ[2, 1, 1, 1, 2, 2, 1, 1] = value
        K_μ[2, 1, 1, 2, 2, 2, 1, 2] = value
        K_μ[2, 1, 2, 1, 2, 2, 2, 1] = value
        K_μ[2, 1, 2, 2, 2, 2, 2, 2] = value
        K_μ[2, 2, 1, 1, 2, 1, 1, 1] = value
        K_μ[2, 2, 1, 2, 2, 1, 1, 2] = value
        K_μ[2, 2, 2, 1, 2, 1, 2, 1] = value
        K_μ[2, 2, 2, 2, 2, 1, 2, 2] = value

        value =
            (-(h[2] * h[3]) / (18 * h[1])) - (h[1] * h[3]) / (9 * h[2]) +
            (h[1] * h[2]) / (36 * h[3])
        K_μ[2, 1, 1, 1, 2, 2, 2, 1] = value
        K_μ[2, 1, 1, 2, 2, 2, 2, 2] = value
        K_μ[2, 1, 2, 1, 2, 2, 1, 1] = value
        K_μ[2, 1, 2, 2, 2, 2, 1, 2] = value
        K_μ[2, 2, 1, 1, 2, 1, 2, 1] = value
        K_μ[2, 2, 1, 2, 2, 1, 2, 2] = value
        K_μ[2, 2, 2, 1, 2, 1, 1, 1] = value
        K_μ[2, 2, 2, 2, 2, 1, 1, 2] = value

        value =
            (-(h[2] * h[3]) / (18 * h[1])) - (h[1] * h[3]) / (18 * h[2]) +
            (h[1] * h[2]) / (18 * h[3])
        K_μ[3, 1, 1, 1, 3, 2, 2, 1] = value
        K_μ[3, 1, 1, 2, 3, 2, 2, 2] = value
        K_μ[3, 1, 2, 1, 3, 2, 1, 1] = value
        K_μ[3, 1, 2, 2, 3, 2, 1, 2] = value
        K_μ[3, 2, 1, 1, 3, 1, 2, 1] = value
        K_μ[3, 2, 1, 2, 3, 1, 2, 2] = value
        K_μ[3, 2, 2, 1, 3, 1, 1, 1] = value
        K_μ[3, 2, 2, 2, 3, 1, 1, 2] = value

        value =
            (-(h[2] * h[3]) / (18 * h[1])) - (h[1] * h[3]) / (36 * h[2]) -
            (h[1] * h[2]) / (36 * h[3])
        K_μ[1, 1, 1, 1, 1, 2, 2, 2] = value
        K_μ[1, 1, 1, 2, 1, 2, 2, 1] = value
        K_μ[1, 1, 2, 1, 1, 2, 1, 2] = value
        K_μ[1, 1, 2, 2, 1, 2, 1, 1] = value
        K_μ[1, 2, 1, 1, 1, 1, 2, 2] = value
        K_μ[1, 2, 1, 2, 1, 1, 2, 1] = value
        K_μ[1, 2, 2, 1, 1, 1, 1, 2] = value
        K_μ[1, 2, 2, 2, 1, 1, 1, 1] = value

        value =
            (-(h[2] * h[3]) / (18 * h[1])) + (h[1] * h[3]) / (36 * h[2]) -
            (h[1] * h[2]) / (9 * h[3])
        K_μ[3, 1, 1, 1, 3, 2, 1, 2] = value
        K_μ[3, 1, 1, 2, 3, 2, 1, 1] = value
        K_μ[3, 1, 2, 1, 3, 2, 2, 2] = value
        K_μ[3, 1, 2, 2, 3, 2, 2, 1] = value
        K_μ[3, 2, 1, 1, 3, 1, 1, 2] = value
        K_μ[3, 2, 1, 2, 3, 1, 1, 1] = value
        K_μ[3, 2, 2, 1, 3, 1, 2, 2] = value
        K_μ[3, 2, 2, 2, 3, 1, 2, 1] = value

        value =
            (-(h[2] * h[3]) / (18 * h[1])) + (h[1] * h[3]) / (18 * h[2]) -
            (h[1] * h[2]) / (18 * h[3])
        K_μ[2, 1, 1, 1, 2, 2, 1, 2] = value
        K_μ[2, 1, 1, 2, 2, 2, 1, 1] = value
        K_μ[2, 1, 2, 1, 2, 2, 2, 2] = value
        K_μ[2, 1, 2, 2, 2, 2, 2, 1] = value
        K_μ[2, 2, 1, 1, 2, 1, 1, 2] = value
        K_μ[2, 2, 1, 2, 2, 1, 1, 1] = value
        K_μ[2, 2, 2, 1, 2, 1, 2, 2] = value
        K_μ[2, 2, 2, 2, 2, 1, 2, 1] = value

        value =
            (-(h[2] * h[3]) / (36 * h[1])) - (h[1] * h[3]) / (18 * h[2]) -
            (h[1] * h[2]) / (36 * h[3])
        K_μ[2, 1, 1, 1, 2, 2, 2, 2] = value
        K_μ[2, 1, 1, 2, 2, 2, 2, 1] = value
        K_μ[2, 1, 2, 1, 2, 2, 1, 2] = value
        K_μ[2, 1, 2, 2, 2, 2, 1, 1] = value
        K_μ[2, 2, 1, 1, 2, 1, 2, 2] = value
        K_μ[2, 2, 1, 2, 2, 1, 2, 1] = value
        K_μ[2, 2, 2, 1, 2, 1, 1, 2] = value
        K_μ[2, 2, 2, 2, 2, 1, 1, 1] = value

        value =
            (-(h[2] * h[3]) / (36 * h[1])) - (h[1] * h[3]) / (36 * h[2]) -
            (h[1] * h[2]) / (18 * h[3])
        K_μ[3, 1, 1, 1, 3, 2, 2, 2] = value
        K_μ[3, 1, 1, 2, 3, 2, 2, 1] = value
        K_μ[3, 1, 2, 1, 3, 2, 1, 2] = value
        K_μ[3, 1, 2, 2, 3, 2, 1, 1] = value
        K_μ[3, 2, 1, 1, 3, 1, 2, 2] = value
        K_μ[3, 2, 1, 2, 3, 1, 2, 1] = value
        K_μ[3, 2, 2, 1, 3, 1, 1, 2] = value
        K_μ[3, 2, 2, 2, 3, 1, 1, 1] = value

        value =
            (h[2] * h[3]) / (36 * h[1]) - (h[1] * h[3]) / (9 * h[2]) -
            (h[1] * h[2]) / (18 * h[3])
        K_μ[2, 1, 1, 1, 2, 1, 2, 2] = value
        K_μ[2, 1, 1, 2, 2, 1, 2, 1] = value
        K_μ[2, 1, 2, 1, 2, 1, 1, 2] = value
        K_μ[2, 1, 2, 2, 2, 1, 1, 1] = value
        K_μ[2, 2, 1, 1, 2, 2, 2, 2] = value
        K_μ[2, 2, 1, 2, 2, 2, 2, 1] = value
        K_μ[2, 2, 2, 1, 2, 2, 1, 2] = value
        K_μ[2, 2, 2, 2, 2, 2, 1, 1] = value

        value =
            (h[2] * h[3]) / (36 * h[1]) - (h[1] * h[3]) / (18 * h[2]) -
            (h[1] * h[2]) / (9 * h[3])
        K_μ[3, 1, 1, 1, 3, 1, 2, 2] = value
        K_μ[3, 1, 1, 2, 3, 1, 2, 1] = value
        K_μ[3, 1, 2, 1, 3, 1, 1, 2] = value
        K_μ[3, 1, 2, 2, 3, 1, 1, 1] = value
        K_μ[3, 2, 1, 1, 3, 2, 2, 2] = value
        K_μ[3, 2, 1, 2, 3, 2, 2, 1] = value
        K_μ[3, 2, 2, 1, 3, 2, 1, 2] = value
        K_μ[3, 2, 2, 2, 3, 2, 1, 1] = value

        value =
            (h[2] * h[3]) / (18 * h[1]) - (2 * h[1] * h[3]) / (9 * h[2]) +
            (h[1] * h[2]) / (18 * h[3])
        K_μ[2, 1, 1, 1, 2, 1, 2, 1] = value
        K_μ[2, 1, 1, 2, 2, 1, 2, 2] = value
        K_μ[2, 1, 2, 1, 2, 1, 1, 1] = value
        K_μ[2, 1, 2, 2, 2, 1, 1, 2] = value
        K_μ[2, 2, 1, 1, 2, 2, 2, 1] = value
        K_μ[2, 2, 1, 2, 2, 2, 2, 2] = value
        K_μ[2, 2, 2, 1, 2, 2, 1, 1] = value
        K_μ[2, 2, 2, 2, 2, 2, 1, 2] = value

        value =
            (h[2] * h[3]) / (18 * h[1]) - (h[1] * h[3]) / (9 * h[2]) +
            (h[1] * h[2]) / (9 * h[3])
        K_μ[3, 1, 1, 1, 3, 1, 2, 1] = value
        K_μ[3, 1, 1, 2, 3, 1, 2, 2] = value
        K_μ[3, 1, 2, 1, 3, 1, 1, 1] = value
        K_μ[3, 1, 2, 2, 3, 1, 1, 2] = value
        K_μ[3, 2, 1, 1, 3, 2, 2, 1] = value
        K_μ[3, 2, 1, 2, 3, 2, 2, 2] = value
        K_μ[3, 2, 2, 1, 3, 2, 1, 1] = value
        K_μ[3, 2, 2, 2, 3, 2, 1, 2] = value

        value =
            (h[2] * h[3]) / (18 * h[1]) - (h[1] * h[3]) / (18 * h[2]) -
            (h[1] * h[2]) / (18 * h[3])
        K_μ[1, 1, 1, 1, 1, 1, 2, 2] = value
        K_μ[1, 1, 1, 2, 1, 1, 2, 1] = value
        K_μ[1, 1, 2, 1, 1, 1, 1, 2] = value
        K_μ[1, 1, 2, 2, 1, 1, 1, 1] = value
        K_μ[1, 2, 1, 1, 1, 2, 2, 2] = value
        K_μ[1, 2, 1, 2, 1, 2, 2, 1] = value
        K_μ[1, 2, 2, 1, 1, 2, 1, 2] = value
        K_μ[1, 2, 2, 2, 1, 2, 1, 1] = value

        value =
            (h[2] * h[3]) / (18 * h[1]) + (h[1] * h[3]) / (18 * h[2]) -
            (2 * h[1] * h[2]) / (9 * h[3])
        K_μ[3, 1, 1, 1, 3, 1, 1, 2] = value
        K_μ[3, 1, 1, 2, 3, 1, 1, 1] = value
        K_μ[3, 1, 2, 1, 3, 1, 2, 2] = value
        K_μ[3, 1, 2, 2, 3, 1, 2, 1] = value
        K_μ[3, 2, 1, 1, 3, 2, 1, 2] = value
        K_μ[3, 2, 1, 2, 3, 2, 1, 1] = value
        K_μ[3, 2, 2, 1, 3, 2, 2, 2] = value
        K_μ[3, 2, 2, 2, 3, 2, 2, 1] = value

        value =
            (h[2] * h[3]) / (18 * h[1]) + (h[1] * h[3]) / (9 * h[2]) -
            (h[1] * h[2]) / (9 * h[3])
        K_μ[2, 1, 1, 1, 2, 1, 1, 2] = value
        K_μ[2, 1, 1, 2, 2, 1, 1, 1] = value
        K_μ[2, 1, 2, 1, 2, 1, 2, 2] = value
        K_μ[2, 1, 2, 2, 2, 1, 2, 1] = value
        K_μ[2, 2, 1, 1, 2, 2, 1, 2] = value
        K_μ[2, 2, 1, 2, 2, 2, 1, 1] = value
        K_μ[2, 2, 2, 1, 2, 2, 2, 2] = value
        K_μ[2, 2, 2, 2, 2, 2, 2, 1] = value

        value =
            (h[2] * h[3]) / (9 * h[1]) - (h[1] * h[3]) / (9 * h[2]) +
            (h[1] * h[2]) / (18 * h[3])
        K_μ[1, 1, 1, 1, 1, 1, 2, 1] = value
        K_μ[1, 1, 1, 2, 1, 1, 2, 2] = value
        K_μ[1, 1, 2, 1, 1, 1, 1, 1] = value
        K_μ[1, 1, 2, 2, 1, 1, 1, 2] = value
        K_μ[1, 2, 1, 1, 1, 2, 2, 1] = value
        K_μ[1, 2, 1, 2, 1, 2, 2, 2] = value
        K_μ[1, 2, 2, 1, 1, 2, 1, 1] = value
        K_μ[1, 2, 2, 2, 1, 2, 1, 2] = value

        value =
            (h[2] * h[3]) / (9 * h[1]) + (h[1] * h[3]) / (18 * h[2]) -
            (h[1] * h[2]) / (9 * h[3])
        K_μ[1, 1, 1, 1, 1, 1, 1, 2] = value
        K_μ[1, 1, 1, 2, 1, 1, 1, 1] = value
        K_μ[1, 1, 2, 1, 1, 1, 2, 2] = value
        K_μ[1, 1, 2, 2, 1, 1, 2, 1] = value
        K_μ[1, 2, 1, 1, 1, 2, 1, 2] = value
        K_μ[1, 2, 1, 2, 1, 2, 1, 1] = value
        K_μ[1, 2, 2, 1, 1, 2, 2, 2] = value
        K_μ[1, 2, 2, 2, 1, 2, 2, 1] = value

        value =
            (h[2] * h[3]) / (9 * h[1]) +
            (h[1] * h[3]) / (9 * h[2]) +
            (2 * h[1] * h[2]) / (9 * h[3])
        K_μ[3, 1, 1, 1, 3, 1, 1, 1] = value
        K_μ[3, 1, 1, 2, 3, 1, 1, 2] = value
        K_μ[3, 1, 2, 1, 3, 1, 2, 1] = value
        K_μ[3, 1, 2, 2, 3, 1, 2, 2] = value
        K_μ[3, 2, 1, 1, 3, 2, 1, 1] = value
        K_μ[3, 2, 1, 2, 3, 2, 1, 2] = value
        K_μ[3, 2, 2, 1, 3, 2, 2, 1] = value
        K_μ[3, 2, 2, 2, 3, 2, 2, 2] = value

        value =
            (h[2] * h[3]) / (9 * h[1]) +
            (2 * h[1] * h[3]) / (9 * h[2]) +
            (h[1] * h[2]) / (9 * h[3])
        K_μ[2, 1, 1, 1, 2, 1, 1, 1] = value
        K_μ[2, 1, 1, 2, 2, 1, 1, 2] = value
        K_μ[2, 1, 2, 1, 2, 1, 2, 1] = value
        K_μ[2, 1, 2, 2, 2, 1, 2, 2] = value
        K_μ[2, 2, 1, 1, 2, 2, 1, 1] = value
        K_μ[2, 2, 1, 2, 2, 2, 1, 2] = value
        K_μ[2, 2, 2, 1, 2, 2, 2, 1] = value
        K_μ[2, 2, 2, 2, 2, 2, 2, 2] = value

        value =
            (2 * h[2] * h[3]) / (9 * h[1]) +
            (h[1] * h[3]) / (9 * h[2]) +
            (h[1] * h[2]) / (9 * h[3])
        K_μ[1, 1, 1, 1, 1, 1, 1, 1] = value
        K_μ[1, 1, 1, 2, 1, 1, 1, 2] = value
        K_μ[1, 1, 2, 1, 1, 1, 2, 1] = value
        K_μ[1, 1, 2, 2, 1, 1, 2, 2] = value
        K_μ[1, 2, 1, 1, 1, 2, 1, 1] = value
        K_μ[1, 2, 1, 2, 1, 2, 1, 2] = value
        K_μ[1, 2, 2, 1, 1, 2, 2, 1] = value
        K_μ[1, 2, 2, 2, 1, 2, 2, 2] = value
    else
    end
    return K_μ
end

@testset "Brick module" begin

    @testset "Gauss–Legendre integration" begin
        order = 2
        for d = 1:3
            @testset "Gauss-Legendre integration, $(d)d" begin
                for n in map(Tuple, CartesianIndices(Tuple(fill(1:(2*order-1), d))))
                    h = Tuple([1.1, 1.2, 1.3][1:d])
                    actual = integrate(x -> prod(x .^ n), h)
                    n1 = n .+ 1
                    expected = prod(((h ./ 2) .^ n1 .- (h ./ (-2)) .^ n1) ./ n1)
                    @test actual ≈ expected rtol = 1e-15 atol = 1e-15
                end
            end
        end
    end

    @testset "Node numbering" begin
        for d = 1:3
            h = Tuple([1.1, 1.2, 1.3][1:d])
            ℒ = cell_vertices(d)
            for p ∈ ℒ
                x = h .* (-1) .^ Tuple(p) ./ 2
                N = shape(x, h)
                for q ∈ ℒ
                    @test N[q] == Float64(p == q)
                end
            end
        end
    end

    @testset "Brick element operators" begin
        for d = 2:3
            @testset "Brick element operators, $(d)d" begin
                # Note: reference values where computed with maxima
                # (see `ref/nd_brick_elasticity.mac`).
                h = Tuple([1.1, 1.2, 1.3][1:d])
                C = Hooke{d,Float64}(5.6, 0.3)

                @testset "Average gradient_operator, $(d)d" begin
                    D_act = avg_gradient_operator(h)
                    D_exp = expected_avg_gradient_operator(h)
                    @test D_act ≈ D_exp rtol = 1e-15 atol = 1e-15
                end

                @testset "Average strain-displacement matrix, $(d)d" begin
                    B_act = avg_strain_displacement_operator(h)
                    B_exp = expected_avg_strain_displacement_operator(h)
                    @test B_act ≈ B_exp rtol = 1e-15 atol = 1e-15
                end

                @testset "Stiffness operator, $(d)d" begin
                    K_act = stiffness_operator(h, C)
                    K_λ = expected_K_λ(h)
                    K_μ = expected_K_μ(h)
                    K_exp = C.λ * K_λ + C.μ * K_μ
                    @test K_act ≈ K_exp rtol = 1e-15 atol = 1e-15
                end
            end
        end
    end
end
