Test coverage for PointEval with nonlinear geometries

[termi-official]

Point evaluation on master just does not work on nonlinear geometries. Here is an attempt to make it work. Fixes #763 but might not completely fix the point eval.

TODO

• docs
• flexibility for inner solver not sure how. Should be a separate PR.
• expose inner solver parameters
• fix point eval for mixed-dimensional and embedded grids

[codecov-commenter]

Codecov Report

Attention: 16 lines in your changes are missing coverage. Please review.

Comparison is base (e806ed9) 93.27% compared to head (8b93b9b) 93.04%.

Files	Patch %	Lines
src/FEValues/GeometryMapping.jl	56.66%	13 Missing ⚠️
src/PointEvalHandler.jl	94.59%	2 Missing ⚠️
src/FEValues/common_values.jl	90.00%	1 Missing ⚠️

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #820      +/-   ##
==========================================
- Coverage   93.27%   93.04%   -0.24%     
==========================================
  Files          36       36              
  Lines        5235     5276      +41     
==========================================
+ Hits         4883     4909      +26     
- Misses        352      367      +15     

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

[termi-official]

This hotfixes dof distribution on embedded elements.

[termi-official]

I moved the computation of J and x without the caches (i.e. GeometryMapping) into GeometryMapping.jl and common_values.jl and changed the signature for the PointEvalHandler to support also mixed-dimensional grids (which are e.g. obtained by mixing shells, beams and solid elements in one grid). These functions are up to further optimization in subsequent PRs (i.e. the stuff which I do in the GPU support PR).

I also fixed an oopsie in the line searcher (which we need to identify the cell boundary for nonlinear cells).

Docs have been added as suggested.

[KnutAM]

Nice work - really cool with the ever-increasing support for embedded elements!

• I would suggest leaving the inclusion of sample_random_point into src to a separate pr for now.

• I like the renaming of x to cell_coordinates, but on the other hand such a renaming is unrelated? If renamed it should probably also be renamed in the reinit! methods? Should this renaming take place in a separate pr then?

I'm unsure regarding iterations for embedded elements. I would have thought we should find the value of $\xi$ which gives zero distance to the sought coordinate, $\boldsymbol{x}_0$, projected onto the element. And then potentially say didn't find if the distance to the element surface is too large.

$$\begin{align} \boldsymbol{x}(\boldsymbol{\xi}) &= \boldsymbol{x}_i N_i(\boldsymbol{\xi})\\\ \boldsymbol{r}(\boldsymbol{\xi}) &= [\boldsymbol{x}(\boldsymbol{\xi}) - \boldsymbol{x}_0] \cdot \frac{\partial\boldsymbol{x}}{\partial \boldsymbol{\xi}}, \quad \text{Potentially normalizing each column in }\frac{\partial\boldsymbol{x}}{\partial \boldsymbol{\xi}} \ \\\ \frac{\partial\boldsymbol{r}}{\partial \boldsymbol{\xi}} &= [\boldsymbol{x}(\boldsymbol{\xi}) - \boldsymbol{x}_0] \cdot \frac{\partial^2\boldsymbol{x}}{\partial \boldsymbol{\xi}^2} + \left[\frac{\partial\boldsymbol{x}}{\partial \boldsymbol{\xi}}\right]^T \cdot \left[ \frac{\partial\boldsymbol{x}}{\partial \boldsymbol{\xi}} \right] \end{align}$$

I suppose the pinv updates are then doing this in some quasi-newton fashion? (But last time I checked, this is quite slow, so probably the above algorithm is faster even with the higher order derivative calculations?)

But since we are using autodiff anyway to compute the gradients of hessians of the base functions, perhaps it would be equally efficient to just define the residual function and just iterate on this? I.e.

function residual(ξ::Vec{dim}, x0::Vec{dim}, ip, cell_coordinates::AbstractVector{<:Vec{dim}) where dim
    return spatial_coordinate(ip, ξ, cell_coordinates) - x0
end
function residual(ξ::Vec{dim}, x0::Vec{dim}, ip, cell_coordinates::AbstractVector{<:Vec{dim}) where dim
    dxdξ, x = gradient(ξ_ -> spatial_coordinate(ip, ξ_, cell_coordinates), :all)
    # do scaling if desired or needed?
    return (x - x0) \cdot dxdξ
end

(And even if not used this modified residual, I suppose the first definition of residual can be used to remove the need for the calculate_mapping_and_spatial_coordinate functions?

[KnutAM]

Change to cell_coordinates here too then?

[termi-official]

Here or separate PR?

[termi-official]

Thanks for taking the time @KnutAM and @fredrikekre !

* I would suggest leaving the inclusion of `sample_random_point` into `src` to a separate pr for now.

* I like the renaming of `x` to `cell_coordinates`, but on the other hand such a renaming is unrelated? If renamed it should probably also be renamed in the `reinit!` methods? Should this renaming take place in a separate pr then?

I am not a big fan of making many small 1-5 line PRs (i.e. add docs, move function with same name from one file to another, ..) due to the associated overhead and dependencies, but I can make many tiny PRs if it is preferred.

I'm unsure regarding iterations for embedded elements. I would have thought we should find the value of ξ which gives zero distance to the sought coordinate, x0, projected onto the element. And then potentially say didn't find if the distance to the element surface is too large.

I suppose the pinv updates are then doing this in some quasi-newton fashion? (But last time I checked, this is quite slow, so probably the above algorithm is faster even with the higher order derivative calculations?)

My initial idea was also to have better modularity for the inner solver, but I could not find a good way how to do it. The rationale behind the current Newton solver for embedded (nonlinear) stuff is that it does not converge to a root if the coordinate is not in the element. You can interpret the solution of the linear system in a least square sense, as pinv returns by default the Moore-Penrose pseudoinverse. The problem with the projection is that you need additional infrastructure to handle the projections properly. Both approaches are valid tho.

But since we are using autodiff anyway to compute the gradients of hessians of the base functions, perhaps it would be equally efficient to just define the residual function and just iterate on this? I.e.

function residual(ξ::Vec{dim}, x0::Vec{dim}, ip, cell_coordinates::AbstractVector{<:Vec{dim}) where dim
    return spatial_coordinate(ip, ξ, cell_coordinates) - x0
end
function residual(ξ::Vec{dim}, x0::Vec{dim}, ip, cell_coordinates::AbstractVector{<:Vec{dim}) where dim
    dxdξ, x = gradient(ξ_ -> spatial_coordinate(ip, ξ_, cell_coordinates), :all)
    # do scaling if desired or needed?
    return (x - x0) \cdot dxdξ
end

(And even if not used this modified residual, I suppose the first definition of residual can be used to remove the need for the calculate_mapping_and_spatial_coordinate functions?

I am not sure if I follow the idea here. So you are proposing to use a gradient descend over a Newton with line search as default? I think this will end up with the same issues as master for points near (or even on) element boundaries.

[KnutAM]

I am not a big fan of making many small 1-5 line PRs (i.e. add docs, move function with same name from one file to another, ..) due to the associated overhead and dependencies, but I can make many tiny PRs if it is preferred.

I rather prefer the opposite, when one PR does one thing (of course - introducing a new feature should also include updated docs and test) and refactoring is left as a standalone pr as this makes the review easier. Furthermore, some changes such as the change from x to cell_coordinates in the function signature are easy and non-controversial and I feel confident review + merge on those. But algorithmic changes / new features etc. might require more discussions.

The reason for suggesting to wait with the inclusion of the testing functions was simply because I think a utils package or including these functions in core is worthwhile discussing, and would be bad to have that discussion stopping this pr.

I am not sure if I follow the idea here. So you are proposing to use a gradient descend over a Newton with line search as default? I think this will end up with the same issues as master for points near (or even on) element boundaries.

No, rather the opposite - the current implementation is a quasi-newton for embedded elements due to the pseudo-inverse. With that code the regular elements would have no change in algorithm, but with easy extension to the embedded elements (but there, it is not guaranteed that the problem is convex / has a solution, not sure if the pseudo inverse approach can circumvent this). I'm struggling with makie-issues, so to take a quick break I'll add a PR to this branch showing what I meant (and check if it actually works :)) (Edit: #876)