Allow evaluating FE functions at arbitrary points #523
Conversation
This is a first cut only, still very inefficient and restrictive.
|
Hi, @eschnett,
I have answered you in #425. Cheers! |
Which function do you need? Perhaps a gridap equivalent is already there.
I would use the caching mechanism we have in gridap. I.e., implement points = ... # Vector{Point{D,Float64}}
uh = ... # CellField
cache = return_cache(uh,testitem(points)) # Generate the cache. i.e., do the setup
for x in points
uh_at_x = evaluate!(cache,uh,x)
end
Yes, the final implementation should only rely on the CellField interface. |
# Conflicts: # src/CellData/CellFields.jl
|
The function I'm using converts a point in a simplex from Cartesian to barycentric coordinates so that I know whether it is inside the simplex or not. This is necessary because the k-d tree search returns the nearest vertex to a give point, and I then have to search all neighbouring simplices. Is there a function for contained-ness? |
|
Perhaps we can consider this package https://github.com/JuliaGeometry/GeometricalPredicates.jl It has no dependencies, so including it is not harmful |
|
That package only supports 2d and 3d simplices. I'll copy the respective code – it's not much, essentially just constructing and inverting a small matrix. |
| @notimplemented """\n | ||
| Evaluation of a CellField at a given Point is not implemented yet. | ||
| cache::CellFieldCache | ||
|
|
There was a problem hiding this comment.
We could probably encapsulate this code in a function and do dispatching with respect to the triangulation type. The code below is for an unstructured (simplicial) mesh. We could also consider the case for a structured Cartesian mesh (which is much easier, no tree involved).
Since this is unrelated to Berstein polynomials I would rather add the code here (I guess it is not much)
No, I am not sure what is the best geometrical predicate for determining whether a a point belongs to an n-simplex. The one you have implemented seems a good option. |
Codecov Report
@@ Coverage Diff @@
## master #523 +/- ##
==========================================
+ Coverage 86.37% 86.53% +0.15%
==========================================
Files 145 145
Lines 11386 11446 +60
==========================================
+ Hits 9835 9905 +70
+ Misses 1551 1541 -10
Continue to review full report at Codecov.
|
|
@santiagobadia Are you volunteering @Balaje to update the PR? Let me know if you need help or get stuck. I'd be happy to do this work myself. |
Evaluate CellField on arbitrary points
|
Hi @eschnett I have already made a PR to your repo with the suggested changes. I see that you have merged the PR. @santiagobadia @ericneiva I have updated |
More tests for evaluate
| @@ -57,6 +61,10 @@ export CellDof | |||
| export get_normal_vector | |||
| export get_cell_measure | |||
|
|
|||
| export make_inverse_table | |||
| export point_to_cell! | |||
There was a problem hiding this comment.
You are exporting point_to_cell! but the function to generate its cache object seems to be private.
Here, there are 2 options:
- Do not export
point_to_cell!and rename it to_point_to_cell!so that it is clear that it's private. - Make public also the mechanism to generate the cache. In this case, follow the API of the abstract type
Mapjust to be consistent with other parts of the code.
struct PointToCellMap{T<:Triangulation} <: Map
trian::T
end
return_cache(k::PointToCellMap,x::Point) = ...
evaluate!(cache,k::PointToCellMap,x::Point) = ...
|
I have created an issue with the last minor comment. I proceed to merge this PR. |
I'm looking for feedback, in particular:
Bernstein.jl, but only call a single function in it. Should I instead import the respective code?kdtree? Can I somehow attach it to the triangulation?I'm aware that the code is still quite restrictive (only handles
SingleFieldFEFunctions with simplicial meshes), and that it uses linear interpolation instead of actually evaluating the basis functions.Will eventually close #425, I hope.