Bug when using Expression and subdomains

Hello,

I have been trying to use subdomains and I noticed a very weird behavior. The minimal version of my code is the following:

from dolfin import *
from math import *
import numpy as np
import sys

sys.getrecursionlimit()
sys.setrecursionlimit(5000)
#parameters["form_compiler"]["precision"] = 100

Xmin = 0
Zmin = 0
Xmax = 7.0e-2
Zmax = 2.0e-2

I = 70
K = 20

# Create mesh
mesh = RectangleMesh(Xmin, Zmin, Xmax, Zmax, I, K, 'left')
gdim = mesh.geometry().dim()

# CROSS-SECTIONS
sigma_a_thin = 2e0
sigma_a_thick = 2e3

# VOLUMETRIC SOURCE
f1 = Expression('1.0')
f2 = Expression('2.0')

class Omega0(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[0] <= 3.0e-2 else False

class Omega1(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[0] >= 3.0e-2 else False

subdomains = MeshFunction('size_t', mesh, 2) # 'size_t' for unsigned int (deprecated version: 'uint'), '2' for dimension of the space
# Mark subdomains with numbers 0 and 1
subdomain0 = Omega0()
subdomain0.mark(subdomains, 0)
subdomain1 = Omega1()
subdomain1.mark(subdomains, 1)


# Define function spaces and mixed (product) space
V0 = FunctionSpace(mesh, 'DG', 1)      # to define the material properties using subdomains

Sigma_a  = Function(V0)
sigma_a_values = [sigma_a_thin,sigma_a_thick]

for cell_no in range(len(subdomains.array())):
    dofs_0 = V0.dofmap().cell_dofs(cell_no)
    subdomain_no = subdomains.array()[cell_no]
    print cell_no,subdomain_no, dofs_0
    Sigma_a.vector()[cell_no] = sigma_a_values[subdomain_no]

print "Everything worked fine"

If I run it, it returns the following error:

...

 1317 0 [4447 4448 4446]  
 1318 37937976 [4444 4445 4443] 

 Traceback (most recent call last):
  File "bugReport1.py", line 56, in <module>
    Sigma_a.vector()[cell_no] = sigma_a_values[subdomain_no]
 IndexError: list index out of range

The value of subdomain_no becomes 37937976 for the cell #1318 for an unknown reason. Surprisingly, if I rerun the code, this behavior can happen for a different cell.

However, if I change

f1 = Expression('1.0')
f2 = Expression('2.0')

in:

f1 = Expression('1.0')
f2 = Constant('2.0')

everything works fine, although f1 and f2 are not used in this portion of code!!

I have tried this for both the version 1.4 and 1.4.0+ of dolfin.

I think I am not the only one to have experienced this bug: see http://fenicsproject.org/qa/5230/defining-a-function-on-a-subdomain?show=5230#q5230

What do you think of that?

Thanks! Vincent

Comments (3)