Failing boundary integral

Reported by Gabriel Balaban
##!/usr/bin/env python
from dolfin import *

parameters["form_compiler"]["representation"] = "uflacs"

BOUNDNUM = 1
class BoundarySubDomain(SubDomain):
    def setlocation(self, location):
        self.location = location
        return self

    def inside(self,x, on_boundary):
                return on_boundary and self.location(x)

T = as_vector((1.0, 0.0))
mesh = UnitSquareMesh(2, 2)
facetFunc = FacetFunction("size_t", mesh)

neumannBoundary = BoundarySubDomain().setlocation(lambda x: near(x[0], 0.0))
neumannBoundary.mark(facetFunc, BOUNDNUM)

ds = Measure("exterior_facet", domain_data = facetFunc)
v = TestFunction(VectorFunctionSpace(mesh, "CG", 1))
f = inner(T, v)*ds(BOUNDNUM)
F = assemble(f)
#Only 3 points should be nonzero
print F.array()
assert(len(filter(lambda x: x != 0, F.array())) == 3)
Comments (2)
