Comment 0 for bug 1206559

Hi,
I modified the

# Heavily borrowed from
# http://fenicsproject.org/pub/book/tutorial/stationary/poisson/mat2_p2D.py
"""
FEniCS tutorial demo: Poisson problem in 2D with 2 materials.
"""

Hi,
I tried to apply the protocol for marking cells using the based on the poisson problem with two materials example (http://fenicsproject.org/pub/book/tutorial/stationary/poisson/mat2_p2D.py) for a 2D mesh I created with gmsh (attached). It appears that the classes for marking the subdomains do not work for my mesh (several elements are missed); unfortunately it is not obvious to me why this is failing. Is there something incorrect about the mesh I am using (it runs fine with other dolfin PDE solves) or the code I have below?

Thanks for your help
pete

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

class Omega0(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[1] <= 0.5 else False

class Omega1(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[1] >= 0.5 else False

if(1):
    nx = 4; ny = 6
    mesh = UnitSquare(nx, ny)
    #mesh = Mesh('nohole.xml')
    mesh = Mesh('hole.xml')

    # Define a MeshFunction over two subdomains
    subdomains = MeshFunction('uint', mesh, 2)

    # Mark subdomains with numbers 0 and 1
    subdomain0 = Omega0()
    subdomain0.mark(subdomains, 0)
    subdomain1 = Omega1()
    subdomain1.mark(subdomains, 1)

    # Discontinuous basis is used for boundary condition
    # LAter, continuous basis used for PDE solution
    V0 = FunctionSpace(mesh, 'DG', 0)
    k = Function(V0)

    print 'mesh:', mesh
    print 'subdomains:', subdomains
    print 'k:', k

    import numpy
    help = numpy.asarray(subdomains.array(), dtype=numpy.int32)
    #print help
    print "bad ", np.where(help>1)
    k.vector()[:] = numpy.choose(help, k_values)