Matrix-valued interpolations

[termi-official]

Foundation for #398 . Also, we basically emulate this in FerriteViz.jl when visualizing stress and in the AMR error estimation when computing the stress error. However, I would like to keep this API internal for now.

TODO

• Tests

[codecov]

Codecov Report

Attention: Patch coverage is 42.85714% with 36 lines in your changes are missing coverage. Please review.

Project coverage is 93.17%. Comparing base (5e894e4) to head (6312f56).
Report is 2 commits behind head on master.

Files	Patch %	Lines
src/interpolations.jl	42.85%	32 Missing ⚠️
src/FEValues/FunctionValues.jl	42.85%	4 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #958      +/-   ##
==========================================
- Coverage   93.71%   93.17%   -0.54%     
==========================================
  Files          36       36              
  Lines        5313     5366      +53     
==========================================
+ Hits         4979     5000      +21     
- Misses        334      366      +32     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[KnutAM]

Very cool to be able to work directly with tensorial fields!
However, would it make sense to rather talk about tensors instead of just a matrix, we could do e.g.

abstract type TensorInterpolation{TB, refshape, iporder} <: Interpolation{refshape, iporder, Nothing} end

struct TensorizedInterpolation{TB, refshape, iporder, SI <: ScalarInterpolation{refshape, iporder}} <: TensorInterpolation{TB, refshape, iporder}
    ip::SI
    function TensorizedInterpolation{TB}(ip::ScalarInterpolation{refshape, iporder}) where {TB<:AbstractTensor{<:Any, <:Any, <:Any}, refshape, iporder}
        return new{TB, refshape, iporder, typeof(ip)}(ip)    
    end
end

where TB is always Tensors.get_base(TT::Type{<:AbstractTensor}).

This would allow

julia> ip = Lagrange{RefTriangle, 1}();

julia> TensorizedInterpolation{Tensor{2,3}}(ip)
TensorizedInterpolation{Tensor{2, 3}, RefTriangle, 1, Lagrange{RefTriangle, 1, Nothing}}(Lagrange{RefTriangle, 1}())

julia> TensorizedInterpolation{SymmetricTensor{2,2}}(ip)
TensorizedInterpolation{SymmetricTensor{2, 2}, RefTriangle, 1, Lagrange{RefTriangle, 1, Nothing}}(Lagrange{RefTriangle, 1}())

This allows re-using e.g. Tensors.n_components and even allows us to initialize the correct type of the shape value with e.g. TB{Float64}

This will also play nicely with mixed tensors, i.e.

julia> struct MixedTensor{order, dims, T, M} <: AbstractTensor{order, dims, T}
           data::NTuple{M, T}
       end

julia> TensorizedInterpolation{MixedTensor{2,(2,3)}}(ip)
TensorizedInterpolation{MixedTensor{2, (2, 3)}, RefTriangle, 1, Lagrange{RefTriangle, 1, Nothing}}(Lagrange{RefTriangle, 1}())

[termi-official]

Very cool to be able to work directly with tensorial fields! However, would it make sense to rather talk about tensors instead of just a matrix, we could do e.g.

abstract type TensorInterpolation{TB, refshape, iporder} <: Interpolation{refshape, iporder, Nothing} end

struct TensorizedInterpolation{TB, refshape, iporder, SI <: ScalarInterpolation{refshape, iporder}} <: TensorInterpolation{TB, refshape, iporder}
    ip::SI
    function TensorizedInterpolation{TB}(ip::ScalarInterpolation{refshape, iporder}) where {TB<:AbstractTensor{<:Any, <:Any, <:Any}, refshape, iporder}
        return new{TB, refshape, iporder, typeof(ip)}(ip)    
    end
end

where TB is always Tensors.get_base(TT::Type{<:AbstractTensor}).

This would allow

julia> ip = Lagrange{RefTriangle, 1}();

julia> TensorizedInterpolation{Tensor{2,3}}(ip)
TensorizedInterpolation{Tensor{2, 3}, RefTriangle, 1, Lagrange{RefTriangle, 1, Nothing}}(Lagrange{RefTriangle, 1}())

julia> TensorizedInterpolation{SymmetricTensor{2,2}}(ip)
TensorizedInterpolation{SymmetricTensor{2, 2}, RefTriangle, 1, Lagrange{RefTriangle, 1, Nothing}}(Lagrange{RefTriangle, 1}())

This allows re-using e.g. Tensors.n_components and even allows us to initialize the correct type of the shape value with e.g. TB{Float64}

This will also play nicely with mixed tensors, i.e.

julia> struct MixedTensor{order, dims, T, M} <: AbstractTensor{order, dims, T}
           data::NTuple{M, T}
       end

julia> TensorizedInterpolation{MixedTensor{2,(2,3)}}(ip)
TensorizedInterpolation{MixedTensor{2, (2, 3)}, RefTriangle, 1, Lagrange{RefTriangle, 1, Nothing}}(Lagrange{RefTriangle, 1}())

Indeed. This is one of the reasons why I would like to keep the API internal for now.I first need to wrap my head around a generic indexing logic and I need to think of a way to properly generalize operations like curl and divergence to be useful by default and to be customizeable if the user works with a slightly different definition.

[KnutAM]

Hm, ok. I'm not sure if I fully understand the purpose? If this is only internal because we expect it to change, then it seems like it could be tested outside as well? I don't see any code that isn't just overloading with new types?

Or is the purpose to move functionality from FerriteViz to Ferrite?

[KnutAM]

I don't think this will do the right thing unfortunately.
But for these matrix interpolations, it probably is better to calculate the gradients wrt. the scalar base-interpolation, and then just insert the values in the right places (and zeros elsewhere).

[termi-official]

Do you have any suggestion how we can cover this case in testing (in addition to the vectorial case)?

[KnutAM]

I think @lijas solution here was quite nice, just having the option of calling that function explicitly as setup here.

Maybe combined with calling something like

for Is in 1:getnbasefunctions(ipm.ip)
    g, v = shape_gradient_and_value(ipm.ip, ξ, Is)
    Im = (Is - 1)*getnbasefunctions(ipm.ip)
    for i in 1:vdim1, j in 1:vdim2
        Im += 1
        G, V = shape_gradient_and_value(ipm.ip, ξ, Im)
        @test G[i,j,:] \approx g
        @test V[i,j] \approx v
    end
end

[termi-official]

If this is only internal because we expect it to change, then it seems like it could be tested outside as well?

What do you mean by outside? If you propose to solve a matrix-valued problem, then I do not have the time currently to implement it.

I don't see any code that isn't just overloading with new types?

Or is the purpose to move functionality from FerriteViz to Ferrite?

We use something very similar to this PR in FerriteViz and this is beneficial to simplify ZZ error estimators. In the future I want to also highlight how we can solve matrix-valued problems (e.g. DG for elasticity) with the API. However, Fredrik repeatedly told me to make smaller PRs to make them more reviewable, which I can fully understand. So here is it. Does this explain it?

[KnutAM]

What do you mean by outside?

I meant that one can use this even without merging this into Ferrite (as opposed to changes to fields of some structs for example). That means that when trying this out (since the way how to implement this is not set as you say), we don't need to merge it in as an experimental feature. But with the code in this PR as a basis, that is great to be able to figure out what the best way is!
And based on this PR, we can then try wrapping our heads around the challenges you mentioned above.