### Expression subclass too slow, can't figure out how to transform to cpp_code

294

views

0

Hi:

I am doing a fluid structure interaction project, where the fluid part is from another fluid package and the solid part is solved by fenics, to get the fluid force acting on solid, I use expression subclass to represent the traction force. Due to the mesh unmatchness between the fluid and solid, I use kdtree to find the traction at x, here is my traction class:

class _tracion(Expression):

def __init__(self, force_coords, force_x, force_y, force_z):

from scipy import spatial

x = force_coords[...,0].flatten()

y = force_coords[...,1].flatten()

z = force_coords[...,2].flatten()

coords=np.column_stack((x,y,z))

self.force_x = force_x.flatten()

self.force_y = force_y.flatten()

self.force_z = force_z.flatten()

self.tree = spatial.KDTree(coords)

def eval(self, value, x):

value[:] = 0

ind = self.tree.query(x)[1]

value[0] = self.force_y[ind] # nearest neighbor

value[1] = self.force_z[ind] # nearest neighbor

value[2] = self.force_x[ind] # nearest neighbor

def value_shape(self):

return (3,)

This class is apparently too slow to evaluate, because of the a callback from C++ to Python and the query of KD tree. If without kd_tree this expression could be easily rewritten to the cpp_code as suggested by document, anyone has some idea of how to make this expression evaluate faster?

Thanks.

Yuxiang

I am doing a fluid structure interaction project, where the fluid part is from another fluid package and the solid part is solved by fenics, to get the fluid force acting on solid, I use expression subclass to represent the traction force. Due to the mesh unmatchness between the fluid and solid, I use kdtree to find the traction at x, here is my traction class:

class _tracion(Expression):

def __init__(self, force_coords, force_x, force_y, force_z):

from scipy import spatial

x = force_coords[...,0].flatten()

y = force_coords[...,1].flatten()

z = force_coords[...,2].flatten()

coords=np.column_stack((x,y,z))

self.force_x = force_x.flatten()

self.force_y = force_y.flatten()

self.force_z = force_z.flatten()

self.tree = spatial.KDTree(coords)

def eval(self, value, x):

value[:] = 0

ind = self.tree.query(x)[1]

value[0] = self.force_y[ind] # nearest neighbor

value[1] = self.force_z[ind] # nearest neighbor

value[2] = self.force_x[ind] # nearest neighbor

def value_shape(self):

return (3,)

This class is apparently too slow to evaluate, because of the a callback from C++ to Python and the query of KD tree. If without kd_tree this expression could be easily rewritten to the cpp_code as suggested by document, anyone has some idea of how to make this expression evaluate faster?

Thanks.

Yuxiang

Community: FEniCS Project

### 1 Answer

1

You can evaluate a Function at an arbitrary point \((x,y,z)\) like

```
from fenics import *
mesh = UnitCubeMesh(2,2,2
V = VectorFunctionSpace(mesh, 'CG', 1)
f = Function(V)
x,y,z = [0.5]*3
print f(x,y,z)
# allow evaluation outside of mesh :
f.set_allow_extrapolation(True)
x,y,z = [1.5]*3
print f(x,y,z)
```

1

Careful with

`f.set_allow_extrapolation(True)`

. Its effect on meshes with non-matching parallel partitionining is generally undesirable.
written
9 months ago by
Jan Blechta

Let me note that evaluating function at point (without cell argument) involves

`BoundingBoxTree`

query. The class can also be used directly if in a need.
written
9 months ago by
Jan Blechta

The binary tree search associated with BBT will be the most efficient way to do this though, correct?

written
9 months ago by
pf4d

Hard to answer. I don't know how

`KDTree`

is implemented.
written
9 months ago by
Jan Blechta

Maybe Yuxiang will perform an analysis.

written
9 months ago by
pf4d

the force and coordinates are from another fluid package not from fenics, so f(x, y, z) won't work in this case. Only solid part are solved by fenics. One has to manually take force out from fluid package and evaluate the x by nearest point(in this class) or fem interpolation

written
9 months ago by
Yuxiang Lin

Any reason not to use FEniCS for your fluid?

written
9 months ago by
pf4d

Well, the fluid is a bit complicated. The fluid is a two phase fluid which involves solving navier-stokes, vof, level-set and auxillary equation capturing the free surface, and include ALE for coupling boundary moving. There is already some codes handling this, therefore I won't start from scratch to build from fenics.

written
9 months ago by
Yuxiang Lin

Well, FEniCS can handle these thing you mention (not sure what vof is though), if you change your mind.

written
9 months ago by
pf4d

Oh, that is volume of fluid, basically for capturing the free surface together with level set equaion

written
9 months ago by
Yuxiang Lin

I wonder do you have a demo case for two phase flow using Fenics?

written
9 months ago by
Yuxiang Lin

There is this:
http://www.karlin.mff.cuni.cz/~hron/fenics-tutorial/multiphase/doc.html
Also the fenics book has a chapter on it I recall...

written
9 months ago by
pf4d

Great, many thanks..

written
9 months ago by
Yuxiang Lin

Please login to add an answer/comment or follow this question.