"IndexError: list index out of range" problem

Question

I'm having difficulty understanding what generates this error message, and how I can avoid it. I'm running a code in FEniCS where I solve consecutive PDEs, and when it comes to attributing different values to the parameter m in different problem subdomains I get this message. What I find unusual is that the subdomains plot looks all fine, but the parameter attribution goes wrong somehow. The problematic piece of code is here:

mesh = Rectangle(0, 0, 1, 0.5, 250, 250)

# define a meshfunction for numbering subdomains
subdomains = MeshFunction("uint", mesh, 2)

#default_value = 0.0
#subdomains.set_all(0)

# define the subdomains
class Biomass(SubDomain):
        def inside(self, x, on_boundary):
                return True if phi1(x) <= 0 else False
class Interface(SubDomain):
        def inside(self, x, on_boundary):
                return True if phi1(x) >= 0 else False
# note: it is essential to use <= and >= in the comparisons

# mark the subdomains with numbers
subdomain0 = Biomass()
subdomain0.mark(subdomains, 0)
subdomain1 = Interface()
subdomain1.mark(subdomains, 1)

# plot the subdomains
plot(subdomains, title="P_t1 subdomains")

# define functionspace for m
V0P1 = FunctionSpace(mesh, "DG", 0)
mP1 = Function(V0P1)

# different parameter values in different subdomains
mP1_values = [1, 0]  # values of m in the two subdomains
for cell_no in range(len(subdomains.array())):
    subdomain_no = subdomains.array()[cell_no]
    mP1.vector()[cell_no] = mP1_values[subdomain_no]

# print values of m
print 'mP1 degree of freedoms:', mP1.vector().array()

(phi1 is the obtained solution to a PDE solved earlier on in the code)

Many thanks,

Comments

1) please provide a minimal runnable code
2) if the code is not runnable you need to point out where the issue is raised, and also a stack trace.

Answer 1

There are cells in your domain that are not marked. I think that:

subdomain0.mark(subdomains, 0)

marks cells with a zero if "inside" returns True for all of the vertices. Hence for:

subdomain0 = Biomass()
subdomain0.mark(subdomains, 0)
subdomain1 = Interface()
subdomain1.mark(subdomains, 1)

then there is an interface between the SubDomain s, consisting of cells which only have some of their vertices in subdomain0 and subdomain1 (but not all of them in only one). Since these are never marked, subdomains contains some uninitialised data.

Comments

Precisely, SubDomain.mark method marks entity iff all its vertices and midpoint lie inside.

---

Many thanks for responding. Problem resolved!

Answer 2

On lines

for cell_no in range(len(subdomains.array())):
    subdomain_no = subdomains.array()[cell_no]
    mP1.vector()[cell_no] = mP1_values[subdomain_no]

you assume that vector DOFs are indexed as cells are which may not be true since FEniCS 1.1.0 (and actually is not, at least in parallel). See http://fenicsproject.org/qa/78/assigning-expansion-coefficients-function-does-parallel how to resolve this.

Comments

Many thanks for responding. Problem resolved!