Hessian in 3D

[bmxam]

MWE failing (watch out for your RAM!) Error:

ERROR: LoadError: MethodError: no method matching _reshape_cell_shape_functions(::StaticArraysCore.SizedMatrix{36, 9, Float64, 2, Matrix{Float64}})

Closest candidates are:
_reshape_cell_shape_functions(::StaticArraysCore.SVector)
@ Bcube /scratchm/mbouyges/.julia/packages/Bcube/Nsl8Q/src/cellfunction/cellfunction.jl:190
_reshape_cell_shape_functions(::StaticArraysCore.SMatrix{Ndof, Ndim}) where {Ndof, Ndim}
@ Bcube /scratchm/mbouyges/.julia/packages/Bcube/Nsl8Q/src/cellfunction/cellfunction.jl:185

using Bcube

degree = 1

mesh = one_cell_mesh(:tetra)

fs = FunctionSpace(:Lagrange, degree)
U = TrialFESpace(fs, mesh)
u = FEFunction(U)

dΩ = Measure(CellDomain(mesh), 2 * degree + 1)

Ugrad = TrialFESpace(fs, mesh; size = 3)
∇ϕ = FEFunction(Ugrad)
projection_l2!(∇ϕ, ∇(u), dΩ)

Uhess = TrialFESpace(fs, mesh; size = 3^2)
∇∇ϕ = FEFunction(Uhess)
projection_l2!(∇∇ϕ, vec ∘ ∇(∇ϕ), dΩ)

[bmxam]

Easier to debug:

module test
using Bcube
using LinearAlgebra

DIM = 3

degree = 1

if DIM == 2
    mesh = one_cell_mesh(:quad)
elseif DIM == 3
    mesh = one_cell_mesh(:tetra)
end

fs = FunctionSpace(:Lagrange, degree)
U = TrialFESpace(fs, mesh)
u = FEFunction(U)

dΩ = Measure(CellDomain(mesh), 2 * degree + 1)

Ugrad = TrialFESpace(fs, mesh; size = DIM)
∇ϕ = FEFunction(Ugrad)
projection_l2!(∇ϕ, ∇(u), dΩ)

Uhess = TrialFESpace(fs, mesh; size = DIM^2)
∇∇ϕ = FEFunction(Uhess)
cInfo = Bcube.CellInfo(mesh, 1)
cPoint = Bcube.CellPoint(fill(0.5, DIM), cInfo, Bcube.ReferenceDomain())

op_cell = Bcube.materialize(vec ∘ ∇(∇ϕ), cInfo)
@show Bcube.materialize(op_cell, cPoint)

op_cell = Bcube.materialize(∇∇ϕ, cInfo)
@show Bcube.materialize(op_cell, cPoint)

l(v) = ∫((vec ∘ ∇(∇ϕ)) ⋅ v)dΩ
b = assemble_linear(l, TestFESpace(Uhess))

end

[bmxam]

problem comes from here : https://github.com/JuliaArrays/StaticArrays.jl/blob/b62e2579f77983f48f7bc52e96d5d31d5d8e42b2/src/linalg.jl#L368

[bmxam]

The original error comes from the fact that Bcube._reshape_cell_shape_functions is expecting a SMatrix and we're giving it a SizedMatrix. This SizedMatrix is actually produced by the kron call in lagrange.jl. This is because StaticArrays has a specific implementation of kron that produces SizedMatrix when the result is larged than a certain limit (see my previous comment).

@ghislainb what do you prefer between:

1. force StaticArrays.kron to give a SMatrix by accessing their internal _kron implementation
2. modify _reshape_cell_shape_functions to accept SizedMatrix

[bmxam]

next problem is the Tuple max number of elements. We have 36 in the present case, leading to endless compil times:

module test
using Bcube
using LinearAlgebra

DIM = 3

degree = 1

if DIM == 2
    mesh = one_cell_mesh(:quad)
elseif DIM == 3
    mesh = one_cell_mesh(:tetra)
end

fs = FunctionSpace(:Lagrange, degree)
U = TrialFESpace(fs, mesh)
u = FEFunction(U)

dΩ = Measure(CellDomain(mesh), 2 * degree + 1)

Ugrad = TrialFESpace(fs, mesh; size = DIM)
∇ϕ = FEFunction(Ugrad)
projection_l2!(∇ϕ, ∇(u), dΩ)

Uhess = TrialFESpace(fs, mesh; size = DIM^2)
∇∇ϕ = FEFunction(Uhess)
cInfo = Bcube.CellInfo(mesh, 1)
cPoint = Bcube.CellPoint(fill(0.5, DIM), cInfo, Bcube.ReferenceDomain())

f(u, v) = u ⋅ v
Vhess = TestFESpace(Uhess)

quadrature = Bcube.get_quadrature(dΩ)
λu, λv = Bcube.blockmap_bilinear_shape_functions(Uhess, Vhess, cInfo)
g1 = Bcube.materialize(f(λu, λv), cInfo)
values = Bcube.integrate_on_ref_element(g1, cInfo, quadrature)
@show values

end

[bmxam]

Origin of the problem : during the assembly, Tuple's of size 36 are created. But Julia Tuple machinery is implemented up to 32 elements, leading to tremendous compilations when the size is larger than this number. An alternative is to perform the L2-projection using a matrix free approach:

function my_l2_projection!(u::Bcube.AbstractSingleFieldFEFunction, f, mesh)
    degree = 2 * Bcube.get_degree(Bcube.get_function_space(Bcube.get_fespace(u))) + 1
    dΩ = Measure(CellDomain(mesh), degree)

    U = get_fespace(u)
    V = TestFESpace(U)
    l(v) = ∫(f ⋅ v)dΩ
    T, = Bcube.get_return_type_and_codim(f, get_domain(dΩ))
    b = assemble_linear(l, V; T = T)

    n = get_ndofs(U)
    function a(x)
        _u = FEFunction(U, x)
        y = assemble_linear(v -> ∫(_u ⋅ v)dΩ, V)
        return y
    end
    A = LinearMap(a, n)

    x = gmres(A, b)
    set_dof_values!(u, x)
end

[bmxam]

Alternative way to proceed (because previous solution with matrix free is very long...). We have to compute the gradient of a gradient. Instead of doing this in one step, we compute the gradient of each direction of the gradient: gradient(gradient_x), gradient(gradient_y), ... And finally we recombine the components in a FEFunction with the following function:

"""
Restricted to case where u is of size `n^2` and a Tuple of size `n` of FEFunction function of size `n`
"""
function combine_components!(u, u_components::Tuple, U, U_components, mesh)
    n_u = length(u_components)
    @assert Bcube.get_ncomponents(u) == n_u^2
    dhl = Bcube._get_dhl(U)
    dhl_components = Bcube._get_dhl(U_components)
    val_u = get_dof_values(u)
    val_u_components = map(get_dof_values, u_components)
    for icell in 1:ncells(mesh)
        for (i, x) in enumerate(val_u_components)
            for j in 1:n_u
                dofs_components = Bcube.get_dof(dhl_components, icell, j)
                dofs = Bcube.get_dof(dhl, icell, (j - 1) * n_u + i)
                val_u[dofs] .= x[dofs_components]
            end
        end
    end
end