Reconstructing a continuous function from finite element method. Is there a faster algorithm for doing so?

Question

Lets say I've decomposed a continuous function $y(x)$ over some domain $L_x$ using known finite element method with local basis $Q_i(x)$. Suppose $L$ is divided into $M$ "elements". If I want to know the function $y(x)$ at a point $x=p$ (where $p$ is not at a nodal value) then I will need to first know which element $p$ is in. I can then interpolate between the element nodes and find the value $y(p)$. I can do this for all $p$ in $L$ and get back my continuous function.

To do this, one could search through every element to find the where $p$ lies. This is fine in simple cases, but what about more complex cases when there are 1000's of elements and we are working in 3D? Is there an efficient algorithm that I'm missing that could compute this?

Answer 1

Don't do it.

On general meshes, finding an arbitrary point $p$ requires $O(N)$ operations where $N$ is the number of cells. You can get away with significantly less, say $O(\log N)$, if you have a mesh that is constructed hierachically from a coarse mesh by mesh refinement, but unless you use a completely structured grid, finding an arbitrary point inside a mesh is simply a very expensive operation.

Secondly, once you have the cell $K$ in which $p$ lies, you need to find the location $\hat p$ on the reference cell $\hat K$ that corresponds to $p$. For this, you have to invert the mapping $\phi_K: \hat K \mapsto K$, and this (in general) requires a nonlinear iteration -- which is also expensive.

As a consequence, most algorithms try to avoid this at any cost. For example, when you do integration of functions, you loop over all cells, and on each cell over the quadrature points (which are defined on the reference cell, and so only need to be mapped forward with $\phi_K$). This is a far faster approach that looping over a bunch of quadrature points defined in real space (as opposed to reference space relative to each cell).

Ultimately, the point is that you need to re-think the overall approach to doing things so that you always start with a loop over all cells, if at all possible.

Answer 2

You could use a domain decomposition approach; divide the domain into coarse subdomains recursively and search through each one (like a binary tree search in 1D).

If you know more about the point $p$, you can probably make it more efficient. For example, if the point $p$ does not change, lookup information can be precomputed and stored during the setup stage. If it changes in time but depends on properties of the solution, it can often be determined by transforming to the reference element and doing the search in parallel.

Answer 3

Similar to comments from Jesse Chan, you could consider breaking up the domain into rectangular (prisms in 3-D) shapes. Then, one option could be to pre-compute a spatial hash table with the finite elements by using a 2-D (or 3-D) coordinate to generate the key. You could then create the table by hashing the vertices of each element.

To efficiently find what element to do computation in, you could hash the coordinates of your arbitrary point $p$ and then loop through the elements in the resulting hash table bucket until you find the element you need to do local computation in. If you build the hash function well, you could end up with very quick look ups.