shape function evaluation, higher derivatives

[edljk]

Is there a way to evaluate second order (for instance) derivatives of shape functions at quadrature points?
Thanks

[termi-official]

We use automatic differentiation internally. For scalar valued problems it should be something like hessian(ξ -> value(func_interpol, i, ξ), ξ, :all). For vector valued problems we have the open issue that 3rd order tensors are not yet supported in Tensors.jl. If your problem depends on matrix valued input, then it should work out fine.

[edljk]

Right, this should give the hessian of a shape function at a point in the reference domain. I was thinking to the evaluation (in the grid space) of the hessian of a fem function at quadrature points. Something similar to a function_gradient call.

[edljk]

OK I understand now that is not so simple: the key part of the evaluation seems to be in the reinit! function where the jacobian is used. A second order derivative would require to collect the hessians cv.d2Ndξ2 as suggested but also the hessians of the mappings cv.d2Mdξ2..
Perhaps there is a naive (less efficient) way to do that. I will have a try.

[lijas]

Yeah, it would require cv.d2Ndξ2 like you said, which we dont store in the CellValues. So you would need to implement something yourself. The expression for d^2NdX^2 can be found here: https://scicomp.stackexchange.com/questions/36778/how-can-i-implement-second-order-derivatives-of-shape-functions-of-a-3d-elements

It would be perhaps be nice to have it as an option in the reinit! function, reinit!(cv, coords, UpdateFlags(; second=true)).

[KristofferC]

Note that you need H2 continuity for your function space in order to be allowed to integrate second derivatives over your domain.

[edljk]

Fine. The use of hessian in fem approach is pretty nonstandard. I was essentially thinking to postprocessing purpose for which efficiency is not the most important issue.

[edljk]

Here is a very first tentative tested on a simple interpolation example.

test_hessian(nel = 50, degcurve = 1, degfem = 3)
  2.912276 seconds (21.54 M allocations: 2.121 GiB, 22.54% gc time)
error_eval = norm(uq - qeval, Inf) = 1.6045110484697034e-6
error_eval∂x = norm(∇uq - q∇eval, Inf) = 0.0006354386768165307
error_evalH = norm(Huq - qHeval, Inf) = 0.1489596236421943
errormean_evalH = mean(norm.(Huq - qHeval, Inf)) = 0.029545458888939925

The code is still very slow with too much allocations. Any suggestions are welcome..

using LinearAlgebra, SparseArrays, StaticArrays
using Ferrite, ForwardDiff
"""
    test_hessian(;              
                 nel::Int64 = 20, degfem::Int64 = 3, degcurve::Int64 = 2, 
                 degquad::Int64 = 10)
        
"""
function test_hessian(;                 
                      nel::Int64 = 20, degfem::Int64 = 3, degcurve::Int64 = 2, 
                      degquad::Int64 = 10, dim::Int64 = 2)
    # initial grid, interpolations, etc.
    celltype = degcurve == 1 ? Ferrite.Triangle : Ferrite.QuadraticTriangle
    grid = generate_grid(celltype, (nel, nel))
    geom_interpol = Lagrange{dim, RefTetrahedron, degcurve}() 
    func_interpol = Lagrange{dim, RefTetrahedron, degfem}()
    quad_rule = QuadratureRule{dim, RefTetrahedron}(:legendre, degquad)
    cellvalues = CellScalarValues(quad_rule, func_interpol, geom_interpol)
    dh = DofHandler(grid)
    push!(dh, :u, 1, func_interpol)
    close!(dh)
    # evaluate points associated to degrees of freedom
    pdof = dof_points_ferrite(dim, dh, cellvalues)
    nbdof = size(pdof, 1)
    # interpoalte function at dof points
    finterp(x) = sin(2 * x[1]) * cos(5 * x[2])
    ∇finterp(x) = ForwardDiff.gradient(finterp, x)
    Hfinterp(x) = ForwardDiff.hessian(finterp, x)
    u = [finterp(pdof[k, :]) for k = 1:nbdof]
    # eval interpolation at quadrature points
    qpoints, qeval, q∇eval, qHeval = derivatives_at_quadpoints(dh, quad_rule, 
                                                func_interpol, geom_interpol,
                                                cellvalues, u) 
    # evalutation of f and its derivatives at quadrature points
    uq = finterp.(qpoints)
    ∇uq = ∇finterp.(qpoints)
    Huq = Hfinterp.(qpoints)
    # compare
    @show error_eval = norm(uq - qeval, Inf)
    @show error_eval∇  = norm(∇uq  - q∇eval, Inf)
    @show error_evalH = norm(Huq - qHeval, Inf)
    @show errormean_evalH = mean(norm.(Huq - qHeval, Inf))
    nothing
end
#-------------------------------------------------------------------------------
value(interpol, j, ξ) = Ferrite.value(interpol, j, Tensors.Vec(ξ...))
function derivatives_at_quadpoints(dh, quad_rule, 
                                   func_interpol::Interpolation,
                                   geom_interpol::Interpolation,
                                   cellvalues::CellValues{dim}, u::Array{T}) where {dim, T}
    n_basefuncs = getnbasefunctions(cellvalues)
    n_geom_basefuncs = getnbasefunctions(geom_interpol)
    global_dofs = zeros(Int, ndofs_per_cell(dh))
    # allocate
    n_qpoints = length(getweights(quad_rule))
    N = fill(zero(T), n_basefuncs, n_qpoints)
    dNdξ = fill(zeros(T, dim), n_basefuncs, n_qpoints)
    d2Ndξ2 = fill(zeros(T, dim, dim), n_basefuncs, n_qpoints)
    M = fill(zero(T), n_geom_basefuncs, n_qpoints)
    dMdξ = fill(zeros(T, dim), n_geom_basefuncs, n_qpoints)
    dM2dξ2 = fill(zeros(T, dim, dim), n_geom_basefuncs, n_qpoints)
    # initialization 
    for (qp, ξq) in enumerate(quad_rule.points)
        for basefunc in 1:getnbasefunctions(func_interpol)
            N[basefunc, qp] = value(func_interpol, basefunc, ξq)
            dNdξ[basefunc, qp] = gradient(ξ -> value(func_interpol, basefunc, 
                                                     ξ), ξq)
            d2Ndξ2[basefunc, qp] = hessian(ξ -> value(func_interpol, basefunc,
                                                      ξ), ξq)
        end
        for basefunc in 1:n_geom_basefuncs
            M[basefunc, qp] = value(geom_interpol, basefunc, ξq) 
            dMdξ[basefunc, qp] = gradient(ξ -> value(geom_interpol, basefunc, 
                                                     ξ), ξq)
            dM2dξ2[basefunc, qp] = hessian(ξ -> value(geom_interpol, basefunc, 
                                                      ξ), ξq)
        end
    end
    d2M = [zeros(T, dim, dim) for _ = 1:dim]
    d2Minvdx2 = [zeros(T, dim, dim) for _ = 1:dim]
    ceval = 0
    nbeval = length(dh.grid.cells) * n_qpoints
    qeval, q∇eval = fill(T(0), nbeval), [zeros(T, dim) for _ = 1:nbeval]
    qHeval = [zeros(T, dim, dim)  for _ = 1:nbeval]
    qpoints = [zeros(T, dim) for _ = 1:nbeval]
    vtmp = zeros(T, dim)
    fecv_J = SizedMatrix{dim, dim}(zeros(T, dim, dim))
    Mtmp = SizedMatrix{dim, dim}(zeros(T, dim, dim))
    nbquadpoints, nbfaces = 0, 0
    #@time begin
    @inbounds for cell in CellIterator(dh)
        reinit!(cellvalues, cell)
        celldofs!(global_dofs, cell)
        for qp = 1:getnquadpoints(cellvalues)
            # global numbering
            ceval += 1
            # compute spatial coordinates of quadpoint
            for (i, x) in enumerate(getcoordinates(cell))
                #qpoints[ceval] += M[i, qp] * x
                copy!(vtmp, x)
                vtmp .*= M[i, qp]
                qpoints[ceval] .+= vtmp
            end
            # compute evaluation
            for i in 1:getnbasefunctions(func_interpol)
                qeval[ceval] += N[i, qp] * u[global_dofs[i]]
            end
            # compute gradient
            fill!(fecv_J, 0)
            for (i, x) in enumerate(getcoordinates(cell))
                #fecv_J += x * dMdξ[i, qp]'
                for j = 1:dim, k = 1:dim
                    fecv_J[j, k] += x[j] * dMdξ[i, qp][k]
                end
            end
            dMinvdx = inv(fecv_J)'
            for i in 1:n_basefuncs
                # q∇eval[ceval] += u[global_dofs[i]] *  dMinvdξ * dNdξ[i, qp]
                mul!(vtmp, dMinvdx, dNdξ[i, qp])
                vtmp .*= u[global_dofs[i]]
                q∇eval[ceval] .+= vtmp
            end
            # compute hessian
            for l = 1:dim
                fill!(d2M[l], 0.)
                for (i, x) in enumerate(getcoordinates(cell))
                    for j = 1:dim, k = 1:dim
                        d2M[l][j, k] += x[k] * dM2dξ2[i, qp][j, l]
                    end                 
                end
            end
            for l = 1:dim
                fill!(d2Minvdx2[l], 0) 
                mul!(Mtmp, d2M[l], dMinvdx)
                mul!(d2Minvdx2[l], dMinvdx, Mtmp)   
                d2Minvdx2[l] .*= - 1.             
            end
            for i in 1:n_basefuncs
                for j = 1:dim, k = 1:dim
                    for l = 1:dim
                        for m = 1:dim
                            qHeval[ceval][j, k] += u[global_dofs[i]] * 
                                                   dNdξ[i, qp][m] * 
                                                   d2Minvdx2[l][k, m] * 
                                                   dMinvdx[j, l]
                            qHeval[ceval][j, k] += u[global_dofs[i]] * 
                                                   d2Ndξ2[i, qp][l, m] *
                                                   dMinvdx[j, l] * 
                                                   dMinvdx[k, m]
                        end
                    end 
                end
            end        
        end 
    end
    #end
    return qpoints, qeval, q∇eval, qHeval
end
#-------------------------------------------------------------------------------
function doassemble_f_ferrite(cellvalues::CellScalarValues{dim}, 
                              dh::DofHandler;
                              ρ::Union{Function, Nothing} = nothing) where {dim}
    f = zeros(ndofs(dh))
    n_basefuncs = getnbasefunctions(cellvalues)
    fe = zeros(n_basefuncs)
    cell_dofs = zeros(Int, n_basefuncs)
    cellnum, valf = 1, 1.
    @inbounds for cell in CellIterator(dh)
        celldofs!(cell_dofs, dh, cellnum)
        cellnum += 1
        fill!(fe, 0)
        reinit!(cellvalues, cell)
        for q_point in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            if ρ != nothing 
                pq = spatial_coordinate(cellvalues, q_point, 
                getcoordinates(cell))
                valf = ρ(pq)
            end
            for i in 1:n_basefuncs
                v  = shape_value(cellvalues, q_point, i)      
                fe[i] += v * valf * dΩ
            end
        end
        # assemble cell contribution
        for (num, dof) in enumerate(cell_dofs)
            f[dof] += fe[num]
        end
    end
    return f
end
#-------------------------------------------------------------------------------
function doassemble_M_ferrite!(M::SparseMatrixCSC, cellvalues::CellScalarValues,
                               dh::DofHandler) 
    n_basefuncs = getnbasefunctions(cellvalues)
    Me = zeros(n_basefuncs, n_basefuncs)
    assembler = start_assemble(M)
    valρ = 1.
    @inbounds for cell in CellIterator(dh)
        fill!(Me, 0)
        reinit!(cellvalues, cell)
        for q_point in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            for i in 1:n_basefuncs
                v  = shape_value(cellvalues, q_point, i) 
                for j in 1:n_basefuncs
                    u = shape_value(cellvalues, q_point, j)                   
                    Me[i, j] += (v * u) * valρ * dΩ
                end
            end
        end
        assemble!(assembler, celldofs(cell), Me)
    end
    return M
end
#-------------------------------------------------------------------------------
function dof_points_ferrite(dim, dh, cellvalues)
    M = create_sparsity_pattern(dh)
    chM = cholesky(doassemble_M_ferrite!(M, cellvalues, dh))
    fx = x -> x[1]
    fy = x -> x[2]
    fxdof = doassemble_f_ferrite(cellvalues, dh, ρ = fx)
    fydof = doassemble_f_ferrite(cellvalues, dh, ρ = fy)
    xdof = chM \ fxdof 
    ydof = chM \ fydof 
    if dim == 2
        return hcat(xdof, ydof)
    elseif dim == 3
        fz = x -> x[3]
        fzdof = doassemble_f_ferrite(cellvalues, dh, ρ = fz)
        zdof = chM \ fzdof 
        return hcat(xdof, ydof, zdof)
    end
    nothing
end

[edljk]

The second order contribution of M was not correct in main function, sorry. Previous tentative only worked with first order boundary.

The problem should be fixed now (code above updated) and the code a little faster

[termi-official]

Do you mind sharing the functional or problem formulation you want to evaluate?

[edljk]

As I said, this is for post processing purpose.
The example I am working on is very specialized: I have 3 eigenvectors (computed by fem) which give me a parametrization of a minimal surface in the ball. I am interested in higher derivatives to evaluate the curvature of that surface.