Make FEFunction objects callable (as functions of coordinates)

[stevengj]

Given an f::FEFunction, it would be nice to be able to evaluate it at a position x::VectorValue just by calling f(x). (Right now this is a MethodError)

(Presumably you have this functionality already, and this is just a one-line method addition?)

[santiagobadia]

Not really, this issue is very related to the Dirac delta one.

A FE function is a cell-wise function. To evaluate f(x) we need first to find which is the cell in the mesh that contains x. (This value can be undefined when is x is on the skeleton mesh and the space is not C0, which is the case of Nedelec, Raviart-Thomas, etc). This is something we do not use because it is expensive and we do not need it in 99% of cases. All the time, our points are already defined in a cell (integration points) and we can skip that part.

But this functionality, i.e., to define a FE function in an arbitrary point, is needed, e.g., when we have a Dirac delta as forcing term.

Once we have defined the algorithm such that given x returns a cell in the mesh such that its closure contains x for the different Gridap meshes, the implementation will be straightforward. That is the missing ingredient. The implementation for Cartesian mesh should be pretty simple if this is what you are using. But I think you were interested on unstructured meshes, from what I read in the Nedelec issue.

[stevengj]

This is something we do not use because it is expensive and we do not need it in 99% of cases. All the time, our points are already defined in a cell (integration points) and we can skip that part.

I agree that you don't need it during assembly etcetera, but the ability to evaluate the solution at any given point is pretty convenient for users doing post-processing and visualization. (e.g. in our case, we wanted to validate our solution by comparing it to a known solution at some set of points.) For example, Fenics has an eval method and also overloads __call__ I think?

(It would also be useful if you needed to interpolate a solution from one mesh to another; e.g. Fenics supports this IIRC.)

To implement it reasonably efficiently (with log(N) complexity) you'd need a tree of some sort, e.g. a kd tree, of course, but an O(N) fallback to start with would be better than nothing. (You could always construct a tree lazily, the first time the user tries to evaluate the function at an arbitrary point.)

(I agree that non-C0 functions are tricky here, but since we are talking about Sobolev spaces here it seems reasonable to return a value that is equal to the finite-element function in the weak sense — e.g. where you just pick one side of the discontinuity at mesh boundaries.)

[santiagobadia]

I agree, I also consider that it is something to be done. But the point to cell map is the part that prevents us to easily provide this functionality now. There is already a k-d-tree searching algorithm implemented in Julia, this package.

[eschnett]

Santiago Badia @santiagobadia said on Gitter: Finite element functions are cell-wise defined, so the first missing ingredient is an efficient searching algorithm that given an x in the physical domain returns the cell it belongs to. With that info, you can just evaluate the (finite element) solution at this point at this cell. There is another minor issue. If you define the finite element space in the physical space, i.e., DomainStyle being PhysicalDomain(), that is all. However, if DomainStyle equal to ReferenceDomain() one needs the inverse geometrical map that takes the physical cell and maps it to the reference one (in which we evaluate shape functions). For structured meshes that is simple and I think that @eneiva has already implemented it in one of his branches.

[ericneiva]

For structured meshes that is simple and I think that @eneiva has already implemented it in one of his branches.

Hi, @eschnett, this functionality is contained in two commits. To bring it into your own devs, you just need to run on your local branch:

git cherry-pick 09804ba2
git cherry-pick b967a908

The first cherry-pick will likely produce a conflict, but it's rather straightforward to solve.

Please, let me know if you have any issues. Regards!

[ericneiva]

@ericneiva I think @eschnett means you

Indeed, @eneiva, thank you and sorry for the inconvenience.

[eschnett]

@ericneiva Thanks!

[fverdugo]

Hi, @eschnett, this functionality is contained in two commits. To bring it into your own devs, you just need to run on your local branch:

git cherry-pick 09804ba
git cherry-pick b967a90
The first cherry-pick will likely produce a conflict, but it's rather straightforward to solve.

These are going to be merged into master via PR #552

[eschnett]

@fverdugo Thanks. I think I have these on my branch https://github.com/eschnett/Gridap.jl/tree/eschnett/evaluate already.