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. How is now included during review.
• expose inner solver parameters
• fix point eval for mixed-dimensional and embedded grids

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 96.00000% with 3 lines in your changes missing coverage. Please review.

Project coverage is 93.67%. Comparing base (1dad1fe) to head (4ee7c78).
Report is 2 commits behind head on master.

Files	Patch %	Lines
src/PointEvalHandler.jl	94.23%	3 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #820      +/-   ##
==========================================
+ Coverage   93.58%   93.67%   +0.08%     
==========================================
  Files          39       39              
  Lines        5895     5990      +95     
==========================================
+ Hits         5517     5611      +94     
- Misses        378      379       +1     

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

[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?

[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)

[KnutAM]

Long time since I looked at this (sorry!)
I think it would be good to move some unrelated changes to a separate PR, that will also keep the git blames easier to follow.

[KnutAM]

A few more comments, but otherwise lgtm!

[fredrikekre]

What is the reason for this change?

[termi-official]

Support for mixed-dimensional meshes. E.g. shells with beams.

[fredrikekre]

Okay, and if there is a field that lives on both domains then?

[termi-official]

I am not sure if I understand the question correctly. The point finder just gives you the reference coordinate and cell index for some spatial coordinate. If you have some spatial coordinate which lives on the boundary of two elements, then it is always up to the prayers of convergence which element is returned for now. Technically we need to distinguish this case from the case when the point is strictly in the interior of some element. But note that resolving this issue is orthogonal to the goal of this PR, which is about improving the handling of nonlinear geometries. Before this PR using either mixed-dimensional or embedded elements resulted in a cryptic error for users and it looks like changing this (and the solver to some pseudo-inverse instead of the actual inverse below) is sufficient to handle these cases. I can separate these 2 lines and the tests into a separate PR if you think it convolutes this PR too much