This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Expression depending on surface-normals of a sphere

+6 votes

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])
asked Nov 24, 2014 by thisch FEniCS Novice (580 points)
edited Nov 25, 2014 by thisch

1 Answer

+3 votes
 
Best answer

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.

answered Nov 25, 2014 by cevito FEniCS User (5,550 points)
selected Nov 25, 2014 by thisch
...