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

11 months ago by

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()
        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?

Community: FEniCS Project

1 Answer

11 months ago by
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 :

x,y,z = [1.5]*3
print f(x,y,z)​
Careful with f.set_allow_extrapolation(True). Its effect on meshes with non-matching  parallel partitionining is generally undesirable.
written 11 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 11 months ago by Jan Blechta  
The binary tree search associated with BBT will be the most efficient way to do this though, correct?
written 11 months ago by pf4d  
Hard to answer. I don't know how KDTree is implemented.
written 11 months ago by Jan Blechta  
Maybe Yuxiang will perform an analysis.
written 11 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 11 months ago by Yuxiang Lin  
Any reason not to use FEniCS for your fluid?
written 11 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 11 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 11 months ago by pf4d  
Oh, that is volume of fluid, basically for capturing the free surface together with level set equaion
written 11 months ago by Yuxiang Lin  
I wonder do you have a demo case for two phase flow using Fenics?
written 11 months ago by Yuxiang Lin  
There is this: Also the fenics book has a chapter on it I recall...
written 11 months ago by pf4d  
Great, many thanks..
written 11 months ago by Yuxiang Lin  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »