Throw when asking for dof indices of vectorized interpolations

[KnutAM]

Master

ip = Lagrange{RefQuadrilateral,2}()
Ferrite.facedof_indices(ip)   # ((1, 2, 3, 4, 5, 6, 7, 8, 9),)
Ferrite.facedof_indices(ip^2) # ((),)

The last call errors with this PR.
This works for <entity>dof_indices, dirichlet_<entity>dof_indices, and <entity>dof_interior_indices.

Example for when this behavior gives a problem,

julia> Ferrite.BCValues(ip, ip, FacetIndex)
Ferrite.BCValues{Float64}([1.0 0.0 0.0; 0.0 1.0 0.0; … ; 0.0 0.0 0.0; 0.0 0.0 0.0;;; 0.0 0.0 0.0; 1.0 0.0 0.0; … ; 0.0 0.0 0.0; 0.0 0.0 0.0;;; 0.0 0.0 0.0; 0.0 0.0 0.0; … ; 0.0 0.0 0.0; 0.0 0.0 0.0;;; 0.0 1.0 0.0; 0.0 0.0 0.0; … ; 0.0 0.0 1.0; 0.0 0.0 0.0], [3, 3, 3, 3], 0)

julia> Ferrite.BCValues(ip^2, ip^2, FacetIndex)
Ferrite.BCValues{Float64}(Array{Float64, 3}(undef, 18, 0, 4), [0, 0, 0, 0], 0)

The last call errors with this PR.

Of course, it could be discussed if these calls should work and return the correct indices as a "true" vector interpolation should (that does make sense to me), but there are some performance pitfalls that could occur if that is passed around, so for now I suggest erroring, and then we can add support for this in the future if desired.

[termi-official]

I think this goes into the wrong direction. Shoudn't we rather fix the dof indices for vectorized interpolations than making them error?

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 92.94%. Comparing base (8c897da) to head (5ef58c9).

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1010      +/-   ##
==========================================
- Coverage   92.94%   92.94%   -0.01%     
==========================================
  Files          38       38              
  Lines        5757     5769      +12     
==========================================
+ Hits         5351     5362      +11     
- Misses        406      407       +1     

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

[KnutAM]

Currently, we special-case a lot of code by using the base interpolation. So I agree that we should support this somehow. But since we currently give the wrong result, I think this solution improves the situation. When including #798 and supporting more features for this (such as L2 projection and vtk export) it will be easier to see how, where and if this makes sense.

Right now, this actually caught the error that the dof-distribution check was skipped for vectorized interpolations here.

While mathematically nice to use VectorizedInterpolations just as a standard vector interpolation, my guess is that there could be several places where this is vdim times more costly (at least, cf. e.g. L2Projection).

For future reference though, here is some code I had to do rather generalize the results.
But I don't think it makes sense to change to code-base this much now right before 1.0 (since this should be possible to do in a no-breaking way later).

# Don't know any function for `jointuple`, ought to exist. 
jointuple(t1::Tuple, t2::Tuple, args::Tuple...) = jointuple((t1..., t2...), args...)
jointuple(t::Tuple) = t

function dofinds(dof_ind_fun::F, ip::VectorizedInterpolation{dim}) where {F, dim}
    return map(inds -> jointuple(map(i -> ntuple(j -> (i-1)*dim + j, dim), inds)...), dof_ind_fun(ip.ip))
end
Ferrite.edgedof_indices(ip::VectorizedInterpolation) = dofinds(Ferrite.edgedof_indices, ip)

[termi-official]

Currently, we special-case a lot of code by using the base interpolation. So I agree that we should support this somehow. But since we currently give the wrong result, I think this solution improves the situation. When including #798 and supporting more features for this (such as L2 projection and vtk export) it will be easier to see how, where and if this makes sense.

[...]
While mathematically nice to use VectorizedInterpolations just as a standard vector interpolation, my guess is that there could be several places where this is vdim times more costly (at least, cf. e.g. L2Projection).

We can still support specializations for these cases. I am more concerned about the API doing the right thing if you want to check e.g. consistency or you want to try out something without worring about the correct way to vectorization (/tensorization) manually.

For future reference though, here is some code I had to do rather generalize the results. But I don't think it makes sense to change to code-base this much now right before 1.0 (since this should be possible to do in a no-breaking way later).

# Don't know any function for `jointuple`, ought to exist. 
jointuple(t1::Tuple, t2::Tuple, args::Tuple...) = jointuple((t1..., t2...), args...)
jointuple(t::Tuple) = t

function dofinds(dof_ind_fun::F, ip::VectorizedInterpolation{dim}) where {F, dim}
    return map(inds -> jointuple(map(i -> ntuple(j -> (i-1)*dim + j, dim), inds)...), dof_ind_fun(ip.ip))
end
Ferrite.edgedof_indices(ip::VectorizedInterpolation) = dofinds(Ferrite.edgedof_indices, ip)

Yea I also ran into the issue on how to represent the vdim properly when designing this, but I stopped when I realized that access to the interpolation is not type stable in the dof handler. Will try to approach this in the future again (or Abd will pick up #829 again).

Should be non-breaking for most users tho.

[KnutAM]

I guess we we both agree that in the future we would like to allow this!
Are you ok with including this fix for now to prevent silent errors?

We could open a separate issue "treat VectorizedInterpolation as true vector interpolation" which would be cool if that doesn't degrade performance and we could get rid of some of the special-casing.

[termi-official]

I guess we we both agree that in the future we would like to allow this! Are you ok with including this fix for now to prevent silent errors?

Yes, good. :)

We could open a separate issue "treat VectorizedInterpolation as true vector interpolation" which would be cool if that doesn't degrade performance and we could get rid of some of the special-casing.

Agreed that we should add an issue (independent of this PR). We will see a performance hit for some operations if we do not use the property, because we increase the number of arithmetic operations. However, I think this should be fine in cases which are not performance-critical to keep the code easier.

[termi-official]

Example for when this behavior gives a problem,

julia> Ferrite.BCValues(ip, ip, FacetIndex)
Ferrite.BCValues{Float64}([1.0 0.0 0.0; 0.0 1.0 0.0; … ; 0.0 0.0 0.0; 0.0 0.0 0.0;;; 0.0 0.0 0.0; 1.0 0.0 0.0; … ; 0.0 0.0 0.0; 0.0 0.0 0.0;;; 0.0 0.0 0.0; 0.0 0.0 0.0; … ; 0.0 0.0 0.0; 0.0 0.0 0.0;;; 0.0 1.0 0.0; 0.0 0.0 0.0; … ; 0.0 0.0 1.0; 0.0 0.0 0.0], [3, 3, 3, 3], 0)

julia> Ferrite.BCValues(ip^2, ip^2, FacetIndex)
Ferrite.BCValues{Float64}(Array{Float64, 3}(undef, 18, 0, 4), [0, 0, 0, 0], 0)

The last call errors with this PR.

Should the last one error? If not, could you add the fix in this PR? If I see it correctly then it should be a straightforward, as there seems to be some dispatch missing in BCValues to de-vectorize the interpolation correctly.

[KnutAM]

Should the last one error?

BCValues are internal, and we de-vectorize here, and the devectorized ip is used for multiple places.

And in the future, if we want to go down the path of treating a vectorized ip as a true vector ip, then BCValues should probably duplicate the reference coordinates etc. vdim times since that is what a vector interpolation would have, so I'm not sure if it is correct to simply devectorize.