Shifting a dolfin.Function by a vector (with UserExpression)

1

I’m trying to define a new function on a finite element space by shifting another function by a vector:

theta_shifted(x) := theta(x + shift)

In case the point x + shift is outside of the domain, these values should be taken from another function.

My idea was to define an expression and interpolate it in the function space. Here is the relevant piece of the code:

class ShiftedExpression(UserExpression):
    def __init__(
            self,
            func_interpolate,
            func_extrapolate,
            shift,
            bounding_box_tree,
            **kwargs
            ):
        self.func_interpolate = func_interpolate
        self.func_extrapolate = func_extrapolate
        self.shift = shift
        self.bounding_box_tree = bounding_box_tree
        super().__init__(**kwargs)

    def eval(self, values, x):
        def is_inside(point):
            collisions = self.bounding_box_tree.compute_collisions(point)
            return len(collisions) > 0

        point = Point(*(x[i] + self.shift[i] for i in range(3)))
        values[0] = self.func_interpolate(point) if is_inside(point) else self.func_extrapolate(point)

mesh = ...
V = FunctionSpace(mesh, "CG", 1)

theta = Function(V)
theta.vector().set_local(...)

theta_shifted_expression = ShiftedExpression(
    func_interpolate=(lambda point: theta(point)),
    func_extrapolate=(lambda _: temp_amb),
    shift=(d, 0, 0),
    bounding_box_tree=mesh.bounding_box_tree(),
    )

theta_shifted_function = interpolate(theta_shifted_expression, V)

The expression theta_shifted_expression works as expected, but I’m not able to get a function on V out of it:

>>> theta_init_new = dolfin.interpolate(theta_shifted, problem.V)

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/interpolation.py", line 71, in interpolate
    Pv.interpolate(v._cpp_object)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/function/function.py", line 365, in interpolate
    self._cpp_object.interpolate(u)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/function/expression.py", line 53, in wrapped_eval
    self.user_expression.eval(values, x)
  File "run.py", line 163, in eval
    values[0] = self.func_interpolate(point) if is_inside(point) else self.func_extrapolate(point)
  File "run.py", line 178, in <lambda>
    func_interpolate=(lambda point: theta_final(point)),
  File "/usr/local/lib/python3.6/dist-packages/dolfin/function/function.py", line 337, in __call__
    self._cpp_object.eval(values, x)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to evaluate function at point.
*** Reason:  The point is not inside the domain. Consider calling "Function::set_allow_extrapolation(true)" on this Function to allow extrapolation.
*** Where:   This error was encountered inside Function.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

I cannot understand at which line and why I’m getting out of the domain.

2

Please make sure that your minimal code example can reproduce your error message for anyone reading the post. See: Read before posting: How do I get my question answered? - #2 for details.

Please note that the bounds of a boundingboxtree is not the same as the bounds of an actual mesh entitiy.

3

Thank you for the answer! You are right, it’s my fault I didn’t provide a reproducible minimal example. However, my question is rather “what is the right way to define a shifted function?” instead of “how to fix this error?”.

Please note that the bounds of a boundingboxtree is not the same as the bounds of an actual mesh entitiy.

Could you please give me a hint on a more correct way to check if a point lies within the boundaries of the mesh?

4

I managed to fix the error and get the result I expected. The function theta needed the set_allow_extrapolation=True parameter (as the error message suggests), which makes sense with @dokken’s remark concerning the bounding_box_tree.

It is interesting that

theta = Function(V, set_allow_extrapolation=True)

does not work, but

theta = Function(V)
theta.set_allow_extrapolation(True)

works as expected.

5

Per the source code for Function, the only kwarg used by Function.__init__ is name, so it is unsurprising that Function(V, set_allow_extrapolation=True) does not work.

6

Yes, you are right, I just realized it and was about to write this remark. I don’t know from where I got the idea to use the keyword parameter, probably just a confusion. Anyway, thanks!

1 Like
7

As a side-note, you could use compute_entity_collisions as shown in: Bitbucket, which would check for actual collision with the cell

8

Thanks for this hint!