Plot FacetNormal

Question

Hi. I'm solving a PDE with a Neumann boundary value that is a scalar product of a function with the normal to the boundary facets. Because I get faulty results, I'm trying to debug my code. My first step is to plot FacetNormal(mesh), but I'm not succeeding.

from dolfin import *

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

f = Expression(('x[0]','x[1]'))
n = FacetNormal(mesh)

# works
File('f.pvd') << project(f, VectorFunctionSpace(bmesh, 'Lagrange', 1))
# error in instant compilation, see below
File('n.pvd') << project(n, VectorFunctionSpace(bmesh, 'Lagrange', 1))

[...]
‘n1’ was not declared in this scope
A[j] += (FE0_C1[0][j]n1 + FE0_C0[0][j]n0)W1det;
[...]
‘n0’ was not declared in this scope
A[j] += (FE0_C1[0][j]n1 + FE0_C0[0][j]n0)W1det;

I expect that I'm missing something obvious. Thanks for your help!

Answer 1

The FacetNormal "lives" only on the facets, so it cannot be evaluated inside the cells which is needed for a projection to normal Lagrange space (although I am not 100% sure how boundary meshes work in this regard)

You could, if you wanted, project to a Discontinuous Lagrange Trace space, that "lives" only on facets, but I do not think dolfin knows how to plot these, so you may not get you far, but you could try.

The FacetNormal has length 1 and is outwards pointing normal to each external facet. For internal facets it is multiple-valued and you need to specify which direction you want. You can potentially use this information to plot it by use of matplotlib. The mesh cells can be queried for their facet normals and facet midpoints for plotting vectors, Something like this code could be adapted:
https://bitbucket.org/trlandet/ocellaris/src/08c0249e8dd827e6dab763d74cf74817b7e77900/ocellaris/postprocess/plot/plot_CR1_2D_vector.py?at=master&fileviewer=file-view-default

Comments

By the way, I do not think the project function supports projection on facets only, so you would need to write your own. I have done so a few times since I have been working with DGT (Discontinuous Galerkin Trace) functions. My (very simplified) pure Python implementation of FEniCS has such a projection implemented which you can see here:
https://bitbucket.org/trlandet/phonyx/src/ad81391ce5de8a9e231eb0d77eb5bd10e66a7e75/phonyx/public_functions.py?at=master&fileviewer=file-view-default
This code should work with with some small changes also in dolfin, at least if you restrict it to external facets only.

Another possibility is that the FacetNormal machinery just does not handle boundary meshes very well. During quadrature the assembler tells the function (n in this case) which coordinate position to calculate and whether or not this is on a facet, and in that case the ID number of the facet. In this case it is not on a facet (in the boundary mesh), but is on a facet (in the mesh used to create the facet normal). You can look at the dolfin code to verify, but I suspect this may not work as intended.