Introduce a convenience method to determine the dimensionality of a CellValues object

[Joroks]

I think it would be good to have convenient acces to the dimensionality of a CellValues object. Especially when someone chooses to use a Scalar CellValues object instead of a Vector CellValues object for the element assembly.

In v0.3.14 I could do something like that:

function assemble_element!(cv::CellScalarValues{dim}, ...) where dim

Which although not documented was quite convenient. However due to the canges made to CellValues this is no longer possible. Now I have to do this:

function assemble_element!(cv::CellValues, ...)
    dim = cv |> Ferrite.function_interpolation |> Ferrite.getrefdim

Maybe it's possible to overload getrefdim() to work for CellValues as well?

Furthermore, the pretty-printing of the CellValues Object already displays this information

julia>     cellValues = CellValues(
               QuadratureRule{RefHexahedron}(:legendre, 2),
               Lagrange{RefHexahedron, 1}()^3
           )
CellValues(vdim=3, rdim=3, and sdim=3): 8 quadrature points
 Function interpolation: Lagrange{RefHexahedron, 1}()^3
 Geometric interpolation: Lagrange{RefHexahedron, 1}()^3

julia>     cellValues = CellValues(
               QuadratureRule{RefHexahedron}(:legendre, 2),
               Lagrange{RefHexahedron, 1}()
           )
CellValues(scalar, rdim=3, and sdim=3): 8 quadrature points
 Function interpolation: Lagrange{RefHexahedron, 1}()
 Geometric interpolation: Lagrange{RefHexahedron, 1}()^3

So it's a bit confusing why you can't easily acces this through the API.

I can create a pull request for this, if you think that this is a good inclusion.

[fredrikekre]

Out of curiousity, what do you use the dimension for?

Maybe it's possible to overload getrefdim() to work for CellValues as well?

Seems like a good idea. We have moved away a bit from extracting things like this from the type parameters (because it becomes a bit brittle if you need to make changes to the parametrization) and instead should prefer extractor functions like getrefdim and such.

I can create a pull request for this, if you think that this is a good inclusion.

That would be great!

[KnutAM]

I think geometric_interpolation(cv) works on master, which could be used to create getrefdim.

[Joroks]

Out of curiousity, what do you use the dimension for?

Originally I used it for this function to make the structure of the resulting element matrix for the FE-LM coupling Matrix more obvious. I also thought this would be slightly faster but after benchmarking it, it turns out, that using vector CellValues is acctually twice as fast.

#using scalar CellValues
function fill_Gcl_elem!(Gcl_elem::Matrix, cv::CellValues, β1::Real, β2::Real)
    fill!(Gcl_elem, 0.0)

    dim = cv |> Ferrite.function_interpolation |> Ferrite.getrefdim
    base_functions = 1:getnbasefunctions(cv)

    for gp in 1:getnquadpoints(cv)
        dΩ = getdetJdV(cv, gp)

        for j in base_functions
            rows = local_dofs(dim, j)

            for k in base_functions
                cols = local_dofs(dim, k)

                valueContribution = shape_value(cv, gp, j) * shape_value(cv, gp, k)
                gradientContribution = shape_gradient(cv, gp, j) ⋅ shape_gradient(cv, gp, k)

                Gcl_elem[rows,cols] += I*(β1*valueContribution + β2*gradientContribution)*dΩ
            end
        end
    end

    return Gcl_elem
end


local_dofs(dim, basefunction) = (1:dim) .+ dim*(basefunction-1)

# using Vector CellValues
function fill_Gcl_elem!(Gcl_elem::Matrix, cv::CellValues, β1::Real, β2::Real)
    fill!(Gcl_elem, 0.0)

    for gp in 1:getnquadpoints(cv)
        for j in 1:getnbasefunctions(cv)
            for k in 1:getnbasefunctions(cv)

                valueContribution = shape_value(cv, gp, j) ⋅ shape_value(cv, gp, k)
                gradientContribution = shape_gradient(cv, gp, j) ⊡ shape_gradient(cv, gp, k)

                Gcl_elem[j,k] += (β1*valueContribution + β2*gradientContribution)*getdetJdV(cv, gp)
            end
        end
    end

    return Gcl_elem
end

I also do something similar for my LM-MD coupling Matrix. But due to how that works, it's more convenient to use the Scaler Values here instead.

[termi-official]

I think it makes sense to provide dispatches for the individual dimensions and document them. It is quite easy to mess up here for some problems by making some implicit assumption about the dimensions being the same.

[fredrikekre]

I also thought this would be slightly faster but after benchmarking it, it turns out, that using vector CellValues is acctually twice as fast.

Yea that doesn't seem so surprising if you dig into the details. This

Gcl_elem[rows,cols] += I*(β1*valueContribution + β2*gradientContribution)*dΩ

is equivalent to this (there is no "update-in-place" operator):

Gcl_elem[rows,cols] = Gcl_elem[rows,cols] + I*(β1*valueContribution + β2*gradientContribution)*dΩ

If we expand a bit more it becomes clear why this is expensive:

tmp1 = Gcl_elem[rows,cols]  # Allocates a new matrix, expensive
tmp2 = I*(...) * dΩ         # UniformScaling with a number, cheap
tmp3 = tmp1 + tmp2          # Allocates a new matrix, expensive
Gcl_elem[rows,cols] = tmp3  # Insert matrix, probably a bit more expensive than the equivalent number of scalar inserts

[Joroks]

Interesting! I wasn't quite sure how julia handles UniformScaling internally.
I've now changed the code to use scalar CellValues but without UniformScaling. It's again about twice as fast as the Vector CellValues solution.

function fill_Gcl_elem!(Gcl_elem::Matrix, cv::CellValues, β1::Real, β2::Real)
    fill!(Gcl_elem, 0)

    dim = cv |> Ferrite.function_interpolation |> Ferrite.getrefdim
    base_functions = 1:getnbasefunctions(cv)

    for gp in 1:getnquadpoints(cv)
        dΩ = getdetJdV(cv, gp)

        for j in base_functions
            rows = local_dofs(dim, j)

            for k in base_functions
                cols = local_dofs(dim, k)

                valueContribution = shape_value(cv, gp, j) * shape_value(cv, gp, k)
                gradientContribution = shape_gradient(cv, gp, j) ⋅ shape_gradient(cv, gp, k)

                value = (β1*valueContribution + β2*gradientContribution)*dΩ

               for (row, col) in zip(rows, cols) 
                    Gcl_elem[row,col] = value
               end
            end
        end
    end

    return Gcl_elem
end

local_dofs(dim, basefunction) = (1:dim) .+ dim*(basefunction-1)

(The matrix assembly is now also almost exactly 100 times faster than the current Matlab Implementation)