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?