Dirichlet boundary conditions of the form dot(b, n): Integral of type cell cannot contain a ReferenceNormal and segfault

[Nico Schlömer]

The following code produces a segfault:

from dolfin import *

b = Constant([1.0, 2.0])

mesh = UnitSquareMesh(20, 20)
V = FunctionSpace(mesh, 'CG', 1)

n = FacetNormal(mesh)
bcs = DirichletBC(V, dot(b, n), 'on_boundary')

Right before that:

Integral of type cell cannot contain a ReferenceNormal.

[Johan Hake]

This is a compound issue. It is the projection of a FacetFunction that fails. The boundary form dot(b,n) gets projected using dot(b,n)vdx. Changing dx to ds will pass the compilation and assemble. So we could expect an error message in UFL. Then it would be nifty to be able to project a boundary form to a Function.

The former is a bug either in UFL or in DOLFIN (Martin?), the latter is a feature request.

[Martin Sandve Alnæs]

The former is this one:

https://bitbucket.org/fenics-project/ufl/issue/6/add-error-message-if-facet-normal-is-used

The latter is a feature request for generating expressions from ufl code that can be evaluated in boundary dofs. Something like that is on my horizon but not my highest priority.

[Jan Blechta]

Note that normal can't be evaluated also in cell integrals on boundary mesh (where it could work)

from dolfin import *

mesh = UnitSquareMesh(4, 4)
bmesh = BoundaryMesh(mesh, 'exterior')

n = FacetNormal(mesh)
assemble(inner(n, n)*dx, mesh=bmesh)

error: ‘n1’ was not declared in this scope
error: ‘n0’ was not declared in this scope

[Jan Blechta]

In the meantime the enclosed hack can be used for a calculation of RT1 outter normal which can be used in DirichletBC

from dolfin import *

#mesh = UnitSquareMesh(4, 4) # Produces funny values in corners
mesh = EllipseMesh(Point(1.0, 2.0), [1.0, 2.2], 0.2)

RT1 = FunctionSpace(mesh, 'RT', 1)
n = FacetNormal(mesh)

# Not working
#n = project(n, RT1)
#n = interpolate(n, RT1)

n_trial, n_test = TrialFunction(RT1), TestFunction(RT1)
a = inner(n_trial, n_test)*ds
L = inner(n,       n_test)*ds

A = assemble(a, keep_diagonal=True)
A.ident_zeros() # Regularize the matrix
b = assemble(L)

#parameters['allow_extrapolation'] = True # May be needed for evaluation tests below
n = Function(RT1)
#parameters['allow_extrapolation'] = False
solve(A, n.vector(), b, 'cg', 'amg') 

# Check values on facet midpoints
#for f in facets(mesh):
#    if f.exterior():
#        p = f.midpoint()
#        print f.index(), [v.index() for v in vertices(f)], n(p)

# Check values at vertices
#for f in facets(mesh):
#    if f.exterior():
#        for v in vertices(f):
#            p = f.midpoint()
#            print v.index(), n(p)

# Check correctness of mean value
area = assemble(Constant(1.)*ds, mesh=mesh)
nds  = assemble(inner(n, n)*ds)
print "Average value of RT1 normal on boundary:", nds / area

print "Plotting. Note that RT1 field is first interpolated onto CG1."
plot(n)
interactive()

print "Now you have approximate RT1 outter normal field which can be used in DirichletBC."

[Nikolaus Rath]

Thanks for the workaround Jan! Are you sure that this is working correctly in 3-D as well? When I replace mesh = EllipseMesh(Point(1.0, 2.0), [1.0, 2.2], 0.2) with mesh = SphereMesh(Point(1.0, 2.0, 0.0), 1, 0.1) the normal vectors in the plot look odd...

[Jan Blechta]

@taerath

1. RT1 normal is definitely only an approximation to proper normal field which is piece-wise constant field on facets.

2. Plots of RT fields may be funny as interpolation to plottable CG1 done behind the scene is often oscillatory near boundary to my opinion. See also this question by @dnarnold.

[Douglas Arnold]

I believe there are at least two quite different problems here. One is that it is not possible to compile expressions like

dot(v, n)*ds

where v is an RT on a manifold and n is the facet normal. This means that you cannot solve a mixed Poisson problem with inhomogeneous Dirichlet data on manifold. This is the issue discussed here.

The second issue is not a problem of compilation. The code runs, but gives incorrect results when solving a mixed Poisson problem with inhomogeneous Neumann on a manifold. This issue is discussed here.

[Jan Blechta]

@dnarnold But it is not obvious that any of your issues make @taerath's normals looking funny. Remember, the example does not concern manifold.

[Nikolaus Rath]

@blechta I need the facet normals to calculate the flux of a vector field through a surface. When I'm using the above workaround, the results look very noisy (compared to the results from a different code), so I don't think the problem is just with the plotting. Is there a particular reason to approximate the normals with RT? For testing, I replaced FunctionSpace(mesh, 'RT', 1) with VectorFunctionSpace(mesh, 'CG', 2). Both the plot and the results of my flux calculation look a lot better with that (though still not perfect).

[Jan Blechta]

@taerath No good reason generally. Just it may be in a good duality to what you're multiplying and integrating.

Yes, I'm not fully convinced that these are just post-processing oscillations. Maybe the procedure for calculation of RT1 normal introduces spurious oscillations or there is some bug connected with Doug's examples. But these seem to be manifold exclusive.

[Nikolaus Rath]

For what it's worth, here is a different workaround that I came up with. At least for my specific case, it gives much better results. It probably makes all sorts of implicit assumptions about my mesh and problem though:

from dolfin import *
import numpy as np

def get_facet_normal(bmesh):
    '''Manually calculate FacetNormal function'''

    if not bmesh.type().dim() == 2:
        raise ValueError('Only works for 2-D mesh')

    vertices = bmesh.coordinates()
    cells = bmesh.cells()

    vec1 = vertices[cells[:, 1]] - vertices[cells[:, 0]]
    vec2 = vertices[cells[:, 2]] - vertices[cells[:, 0]]

    normals = np.cross(vec1, vec2)
    normals /= np.sqrt((normals**2).sum(axis=1))[:, np.newaxis]

    # Ensure outward pointing normal
    bmesh.init_cell_orientations(Expression(('x[0]', 'x[1]', 'x[2]')))
    normals[bmesh.cell_orientations() == 1] *= -1

    V = VectorFunctionSpace(bmesh, 'DG', 0)
    norm = Function(V)
    nv = norm.vector()

    for n in (0,1,2):
        dofmap = V.sub(n).dofmap()
        for i in xrange(dofmap.global_dimension()):
            dof_indices = dofmap.cell_dofs(i)
            assert len(dof_indices) == 1
            nv[dof_indices[0]] = normals[i, n]

    return norm


mesh = SphereMesh(Point(1.0, 2.0, 0.0), 1, 0.1)
bmesh = BoundaryMesh(mesh, 'exterior')
n = get_facet_normal(bmesh)

# Check values on facet midpoints
for f in facets(bmesh):
    p = f.midpoint()
    print f.index(), [v.index() for v in vertices(f)], n(p)

# Check correctness of mean value
area = assemble(Constant(1.)*dx, mesh=bmesh)
nds  = assemble(inner(n, n)*dx)
print "Average value of RT1 normal on boundary:", nds / area

plot(n)
interactive()

[Anders Logg (Chalmers)]

Removing milestone: 1.7 (automated comment)

[Nico Schlömer]

Update MWE