Unable to use conditional module

I/O
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