[BREAKING] Embedded elements

[termi-official]

Cleaner integration of embedded elements. Decouples the reference from the spatial dimension. Closes #619 .

cc @lijas may need your help regarding the shell example.

TODO

1. Make things functional
2. Fix performance regression
3. Merge individual reinit! functions into a single one to enforce DRY.
4. Migrate point values to the new interface After going through the related code I think this is a PR by itself.
5. Improved doc strings
6. Update shell example Not sure how.

[codecov-commenter]

Codecov Report

Patch coverage: 99.13% and project coverage change: +0.11 🎉

Comparison is base (9b3b827) 92.27% compared to head (ff81ee5) 92.38%.

❗ Current head ff81ee5 differs from pull request most recent head 137cd6c. Consider uploading reports for the commit 137cd6c to get more accurate results

❗ Your organization is not using the GitHub App Integration. As a result you may experience degraded service beginning May 15th. Please install the Github App Integration for your organization. Read more.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #651      +/-   ##
==========================================
+ Coverage   92.27%   92.38%   +0.11%     
==========================================
  Files          30       30              
  Lines        4439     4520      +81     
==========================================
+ Hits         4096     4176      +80     
- Misses        343      344       +1     

Impacted Files	Coverage Δ
src/Ferrite.jl	100.00% <ø> (ø)
src/interpolations.jl	97.69% <94.11%> (-0.14%)	⬇️
src/FEValues/cell_values.jl	100.00% <100.00%> (ø)
src/FEValues/face_values.jl	100.00% <100.00%> (ø)
src/PointEval/point_values.jl	97.36% <100.00%> (ø)
src/Quadrature/quadrature.jl	78.57% <100.00%> (+0.31%)	⬆️

☔ View full report in Codecov by Sentry.
📢 Do you have feedback about the report comment? Let us know in this issue.

[fredrikekre]

Perhaps we can change the storage in CellValues to SVectors and SMatrix unconditionally, and then cast to Tensors.Vec and Tensors.Tensor when possible? Since the storage of StaticArrays and Tensors is equivalent this should be a no-op. Something like this:

using Tensors
using Tensors.StaticArrays

struct CellVectorValues{rdim,sdim,vdim,T,M1,M2}
    N::Matrix{SVector{vdim,T}}
    dNdx::Matrix{SMatrix{vdim,sdim,T,M1}}
    dNdξ::Matrix{SMatrix{vdim,rdim,T,M2}}
end

function shape_value(cv::CellVectorValues, qp::Int, i::Int)
    return tensor_cast(cv.N[qp, i])
end
function shape_derivative(cv::CellVectorValues, qp::Int, i::Int)
    return tensor_cast(cv.dNdx[qp, i])
end

# Cast to Tensors.(Vec|Tensor) when possible
for dim in 1:3 @eval begin
    @inline function tensor_cast(v::SVector{$dim,T}) where T
        return Vec{$dim,T}(Tuple(v))
    end
    @inline function tensor_cast(v::SMatrix{$dim,$dim,T,M}) where {T,M}
        return Tensor{2,$dim,T,M}(Tuple(v))
    end
end end

@inline tensor_cast(x) = x # Fallback

function make_values(rdim,sdim,vdim)
    n_qp = 2
    n_shapefuncs = 2
    N = rand(SVector{vdim,Float64}, n_qp, n_shapefuncs)
    dNdx = rand(SMatrix{vdim,sdim,Float64}, n_qp, n_shapefuncs)
    dNdξ = rand(SMatrix{vdim,rdim,Float64}, n_qp, n_shapefuncs)
    return CellVectorValues{rdim,sdim,vdim,Float64,vdim*sdim,vdim*rdim}(N, dNdx, dNdξ)
end

cv_standard = make_values(3, 3, 3)
cv_embedded = make_values(2, 3, 4)

shape_value(cv_standard, 1, 1)        # -> Vec{3}
shape_derivative(cv_standard, 1, 1)   # -> Tensor{2, 3}

shape_value(cv_embedded, 1, 1)        # -> SVector{4}
shape_derivative(cv_embedded, 1, 1)   # -> SMatrix{4, 3}

[termi-official]

Okay, I think I have broken the judgement of the benchmarks when building the makefile in #388 . Will fix this. I follow up with a quick manual benchmark.

using Ferrite
using BenchmarkTools

function setup_dhclosed(grid, T)
    dh = T(grid)
    add!(dh, :v, 2, Lagrange{2,RefCube,2}()) 
    add!(dh, :s, 1, Lagrange{2,RefCube,1}())
    close!(dh)
    return dh
end

dim = 2
ip = Lagrange{dim, RefCube, 1}()
qr = QuadratureRule{dim, RefCube}(2)
cellvalues = CellScalarValues(qr, ip);
grid = generate_grid(Quadrilateral, (1000, 1000));
dh = setup_dhclosed(grid, DofHandler);

@btime CellScalarValues($qr, $ip)

cell = first(CellIterator(dh))
@btime reinit!(cellvalues, cell)

Master

• ctor: 608.819 ns (6 allocations: 1.45 KiB)
• reinit!: 104.871 ns (0 allocations: 0 bytes)

This PR

• ctor: 3.104 μs (18 allocations: 2.11 KiB)
• reinit!: 106.253 ns (0 allocations: 0 bytes)

I think the regression on the ctor is manageable.

[KristofferC]

Would be good to document in what way this is breaking and how to update (e.g. in the changelog).

[termi-official]

I will add more docs once the design settles. We probably want to fix the vectorization of scalar interpolations first, so I can remove the remaining hacks from this PR.

[termi-official]

Converted back to draft as this one has some weird performance regression which I do not understand (docs build times went drastically up). I hope that this is fully resolved now (at least locally).

[termi-official]

Is it intended that all data types for N,M, the derivatives and the quadrature should be of type T?

[termi-official]

Maybe we want to fix #676 first such that we can eliminate the introduced extra parameters sdim and vdim in the CellScalarValues and CellVectorValues constructor.

[fredrikekre]

Maybe we want to fix #676 first such that we can eliminate the introduced extra parameters sdim and vdim in the CellScalarValues and CellVectorValues constructor.

That would let you query vdim from the interpolation, but not sdim, unless the way to opt-in to having sdim > rdim is to vectorize the geometrical interpolation.

For example, 4D function in on a 2D (reference) element embedded in 3D would be something like:

interpolation_fun = VectorizedInterpolation{vdim=4, Lagrange{refdim=2,RefCube,1}()}()
interpolation_geo = VectorizedInterpolation{sdim=3, Lagrange{refdim=2,RefCube,1}()}()

which would give you all the required dimensions. The geometric interpolation can be unpacked internally, because we don't need it to be vectorized, it would just be a way to pass the dimension.

[termi-official]

So I have updated the parts which I am confident about. I was not able to fully update the shell example, because I do not fully understand all assumptions taken here. Embedded problems work now. I can also provide an example if requested.

[fredrikekre]

This is not a very satisfying way of specifying the space dimension, but I am not sure what would be better. Perhaps requiring passing the vectorized geometric interpolation would be an idea?

qr = QuadratureRule{RefTriangle}(1)
ip = Lagrange{RefTriange, 1}()
ipg = Lagrange{RefTriangle, 1}() ^ sdim

CellValues(qr, ip, ipg)

and extracting sdim from ipg while unwrapping it to use the scalar versions, still, internally. Thoughts?

[termi-official]

That is what I wanted to do in the first place, but I also wanted to avoid having these weird inconsistencies that we pass the scalar geometric interpolation at some points and the vectorized at other points.

[fredrikekre]

Looks like construction is still broken for these two cases with vectorized interpolations:

qr = QuadratureRule{2,RefQuadrilateral}(1)
ip = Lagrange{RefQuadrilateral,1}()

CellValues(qr, ip^1)
CellValues(3, qr, ip^1)

[fredrikekre]

I pushed a proposed fix in 291fafb but didn't add any tests just yet.

[fredrikekre]

Very cool! You can now have e.g. a 2D vector field on a 1D line embedded in 3D space :)

[KristofferC]

You can now have e.g. a 2D vector field on a 1D line embedded in 3D space :)

Examples needed! :D

[termi-official]

You can now have e.g. a 2D vector field on a 1D line embedded in 3D space :)

Examples needed! :D

As you wish #715 :P