Expression depending on surface-normals of a sphere

Question

I want to create a dolfin Expression on the BoundaryMesh of a 3D Volume mesh. The Expression depends on the surface-normals. Based on the comment https://bitbucket.org/fenics-project/dolfin/issue/53#comment-9001361 , how can I use the Function returned by get_facet_normal(bmesh) in the eval_cell method of a dolfin Expression?

from dolfin import *
import numpy as np
import mshr
# for simplicity I use a sphere here, but the code should work with any 3d volume mesh
c = mshr.Sphere(dolfin.Point(0., 0., 0.), 1.)

mesh = mshr.generate_mesh(c, 10)
bmesh = BoundaryMesh(mesh, 'exterior')


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

n = get_facet_normal(bmesh)

class BoundarySource(Expression):
    def __init__(self, mesh):
        self.mesh = mesh

    def eval_cell(self, values, x):
        cur_cell_normal = ?
        a, b, c = 1.0
        values[0] = a*cur_cell_normal[0] + b*cur_cell_normal[1] + c*cur_cell_normal[2]

f = BoundarySource(bmesh)
f(bmesh.coordinates()[0])

Answer 1

Hi Thomas,

since the normal is constant for every cell of your mesh you are correct that you should use eval_cell in the Expression, however, the correct function arguments in this case include a ufc_cell. You can use the ufc_cell to then evaluate the normal of the facet for the particular ufc_cell. Note, that the value for x doesn't matter in this case, since the normal is constant on the facet anyway:

class BoundarySource(Expression):
    def __init__(self, mesh):
        self.mesh = mesh

    def eval_cell(self, values, x, ufc_cell):
        cur_cell_normal = np.empty((3,), dtype=np.float)
        n.eval_cell(cur_cell_normal, np.array([0.0, 0.0, 0.0]), ufc_cell)
        a, b, c = 1.0, 1.0, 1.0
        values[0] = a*cur_cell_normal[0] + b*cur_cell_normal[1] + c*cur_cell_normal[2]

Note, that you cannot evaluate this function by simply providing a coordinate. Instead you should always evaluate the function using eval_cell. This is automatically used when passing the Expression to interpolate or any other dolfin routine.