mesh initialisation with iterators

[Richard Katz]

(using latest release of Mac binary for osx 10.6)

The code below uses a function on the boundary facets to mark certain facets for a Dirichlet boundary condition. This ONLY works when I uncomment the line
mesh.init()
in the code below.

If this line is commented, I get the error:
*** Warning: Found no facets matching domain for boundary condition.

Rich

#################### CODE BELOW #####################

from dolfin import *

# Create mesh
mesh = UnitCircleMesh(20)

# Create function space
degree = 2
V = VectorFunctionSpace(mesh,"Lagrange",degree)
Q = FunctionSpace(mesh,"Lagrange",degree-1)
W = MixedFunctionSpace([V,Q])
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0.0, 0.0))
U = Function(W)

facet_function = FacetFunction("size_t", mesh, 0)
#mesh.init() ## This MUST be uncommented to make the code work
D = mesh.topology().dim()
for facet in facets(mesh):
    if facet.num_global_entities(D) == 1:
        if facet.normal()[1]>0.5: facet_function[facet.index()] = 1

class BCExpression(Expression):
       def __init__(self, mesh):
           self.mesh = mesh
       def eval_cell(self, values, x, ufc_cell):
           cell = Cell(self.mesh, ufc_cell.index)
           n = cell.normal(ufc_cell.local_facet)
           values[0] = n[0]
           values[1] = n[1]
       def value_shape(self):
           return (2,)
bcvel = BCExpression(mesh)

# Dirichlet bc
bcs = DirichletBC(W.sub(0), bcvel, facet_function, 1)

# Procedure B: Set operators and solve
sgu = sym(grad(u))
sgv = sym(grad(v))
a = ( + sgu[i,j]*sgv[i,j]*dx # stress tensor
      - p*v[i].dx(i)*dx # pressure gradient
      + q*u[i].dx(i)*dx ) # velocity divergence
L = f[i]*v[i]*dx # body-force
solve(a==L,U,bcs)