Unable to use conditional module

I/O
1

Dear all,
I’m trying to use the conditional Module but, obviously, I’m doing something wrong.
My problem could be summarised as follow: given a certain condition, I would like to get the value 1 if the condition is respected, 0 otherwise.
I tried to do so with the minimal code posted below:

from dolfin import *
from mshr import *
import numpy as np

mesh = UnitSquareMesh(1, 1)
V = FunctionSpace(mesh, "DG", 1)
rv = Function(V)

# conditon to be checked: rv >= 2
cond = ge(rv-2,DOLFIN_EPS) 

# imposing the values above the cond
rv.vector()[:] = 2.1

# setting only the frist DoF below the cond
rv.vector()[0] = 1.9 
# 1.9 gives 1 in position 0 
# 1.8 gives 0 in positions 0, 1 and 2

# conditional statement
a = conditional(cond,1,0)

# output
np.savetxt("rv.csv", project(rv,V).vector().get_local(), delimiter=",")
np.savetxt("a.csv", project(a,V).vector().get_local(), delimiter=",")

I was expecting to get the value 0 in the first position of a. Instead, in the case (rv.vector()[0] = 1.9 ) I obtain 1 and in the other case (rv.vector()[0] = 1.8 ) the value 0 appears in the first three positions of a.
Any hint?
Thank you very much.

$ fenics -version
2019.1.0

2

This is due to the fact that a projection does not evaluate the function at the degrees of freedom, but at the quadrature points. You observe this when simply assembling your conditional with various quadrature degrees:

from dolfin import *
import numpy as np

mesh = UnitSquareMesh(1, 1)
V = FunctionSpace(mesh, "DG", 1)
rv = Function(V)

# conditon to be checked: rv >= 2
cond = ge(rv-2,DOLFIN_EPS) 

# imposing the values above the cond
rv.vector()[:] = 2.1

# setting only the frist DoF below the cond
rv.vector()[0] = 1.9 
# 1.9 gives 1 in position 0 
# 1.8 gives 0 in positions 0, 1 and 2

# conditional statement
a = conditional(cond,1,0)
for degree in range(15):
    print(degree, assemble(a*dx(domain=mesh, metadata={"quadrature_degree": degree})))

giving:

0 1.0
1 1.0
2 0.8333333333333335
3 0.8333333333333333
4 0.9450241281723379
5 0.9370304097275864
6 0.833333333333335
7 0.8968584753861697
8 0.8882208207693354
9 0.8882208207693354
10 0.8586167080281739
11 0.8586167080281739
12 0.8768368548676453
13 0.8768368548676453
14 0.8791552305348609

If you use the conditional in a variational form, it will be evaluated at the quadrature point, and you will get the exact value at this point.

2 Likes