Dear All,
I use FEniCS to implement a wide stencil scheme and have the following bottleneck in my computation: Given a point v somewhere in the computational domain, I need to find the element containing it and the local (barycentric) coordinates of v in that element. I have written the following code, which works but is very slowly:
def interpolate3(self, v):
# given point coordinates return the interpolation weights
cell = -1 # set to (non-negative) cell number when found
x = [0.0, 0.0, 0.0]
# Bounding box tree to determine relevant cells
for c in self.bbt.compute_collisions(Point(v[0], v[1], v[2])):
# The barycentric coordinates are the interpolation weights
localvertices = self.coord[self.c2v[c]] # vertex coordinates of cell
lA = localvertices[1:4] - localvertices[0] # translate vertices 'towards origin'
lb = v - localvertices[0] # translate first stencil point
x = np.linalg.solve(lA.T, lb) # x in local coordinates
if x.min() >= 0.0 - self.TOL and x.sum() <= 1.0 + self.TOL: # check if in cell
cell = c
break
return [cell, x]
In words, I search with the bounding box tree which elements c may contain v and then solve a linear system to check whether (a) that is really the case and (b) find the local coordinates. However, for a mesh with a 10^6 DoFs, say, I might call this function 10^8 times, leading to a run time in days on my iMac.
Is there a faster way to do this?
