Creating an Expression based on a conditional

General
1

Dear FeniCs community,

I’m currently solving a two-phased Fluid Flow problem, using a level-set formulation. This means that I have a function \phi(x,y) that takes the value of -1 or of 1, depending on where on the domain we are. (On one portion of the domain phi = 1, on the other one, phi = -1).

In order to solve my problem, I need to apply an external force, however, said force must only be applied on the parts of the domain where \phi < 0. I’m trying to formulate the force as a FeniCs expression using a conditional, however, I haven’t been able to. My attempt looks like this at the moment.

from dolfin import *
# Define parameters
R = 0.1 # Radius of the circle
delta = 0.05 # Width of the interface

# Define initial condition class
class InitialConditions(UserExpression):
    def __init__(self, R, delta, **kwargs):
        self.R = R
        self.delta = delta
        super().__init__(**kwargs)

    def eval(self, values, x):
        d = sqrt((x[0] - 0.5)**2 + (x[1] - 0.8)**2) - self.R
        values[0] = tanh(d / self.delta)

    def value_shape(self):
        return ()

initial_condition = InitialConditions(R, delta, degree=1)
W = FunctionSpace(mesh, 'P', 2)
phi = TrialFunction(W)
phi_n = interpolate(initial_condition, W)  # Initial condition
w = TestFunction(W)

f  = Constant((0,-1))  #Fuerza constante, la idea es ir cambiando este parametro a ver que ocurre
phi_values = phi_n.vector().get_local()
force_multiplier = conditional(phi < 0, 1.0, 0.0)  # 1.0 where phi < 0, 0.0 otherwise
f_modified = Expression(('f_multiplier * f_x', 'f_multiplier * f_y'), degree=1, f_multiplier=force_multiplier, f_x=f[0], f_y=f[1])

However, I get the following error:
RuntimeError: Unable to compile C++ code with dijitso

Any pointers on how I can either fix this problem, or implement what I need?

2

You’re mixing UFL symbolic algebra with C/C++ code in the expression. Consider using a ternary operator in the Expression’s C/C++ code, or define your own complete C++ dolfin::Expression class yourself with appropriate implementation of conditional statements.