Dirichlet condition with normal vector

Question

Hi all,

I need to solve the Laplace/Poisson equation on some custom mesh omega (generated by gmsh) where in some of my DirichletBCs the normal vector is used. See:

V = VectorFunctionSpace(omega, "Lagrange", 1)

w = TrialFunction(V)
v = TestFunction(V)
n = FacetNormal(omega)

#not working:
#n = project(n,V)
#n = interpolate(n,V)

f = Constant((0.0, 0.0))

bc1 = DirichletBC(V, zero_v, boundaries, 2)
bc2 = DirichletBC(V, zero_v, boundaries, 3)
bc3 = DirichletBC(V, zero_v, boundaries, 4)
bc4 = DirichletBC(V, zero_v, boundaries, 5)
bc5 = DirichletBC(V, n, boundaries, 6)
bc6 = DirichletBC(V, n, boundaries, 7)
bc7 = DirichletBC(V, zero_v, boundaries, 8)

bcs = [bc1, bc2, bc3, bc4, bc5, bc6, bc7]

a = inner(grad(w), grad(v))*dx
L = inner(f,v)*dx

w = Function(V)
solve(a == L, w, bcs)

However, using

n = FacetNormal(omega)

doesn't work, as I know by searching around that there is no eval - Routine for the FacetNormal.
These boundary conditions however are criticial to my project, so I need to have atleast some workaround, maybe calculating the FacetNormal manually so that I can use them in my DirichletBC, but I don't have an idea how to do this or adding some evalute-routine manually, but I don't know where to start with the latter (already seen: https://fenicsproject.org/qa/11723/boundary-condition-that-depends-on-facet-normal )

Any ideas?

Comments

Would defining your normal like n = Constant((0,0,1)) suffice? Since you can project that on your functionspace.

---

Hi,

thanks for your answer. Sadly this doesn't work, as I'm getting unexpected results -which isn't too bad, but the resulting function with these "new" BCs doesn't fit my needs.

Answer 1

Hi there,

if you want to impose a Neumann BC on parts 6 and 7 of your boundary in FEniCS, have a look in the relevant part of the tutorial.

Does this help?

Comments

Just to be clear, a Neumann BC is one where one specifies the normal derivative of the solution u on that part of the boundary- in contrast to a Dirichlet BC, where the a function for u itself is prescribed.

Answer 2

Someone helpfully wrote a function here which manually calculates a the normal vector function along the boundary (as pictured in their output image).

You can run the following code for illustration; I've made a couple of modifications so that it works in version 2016.2 (e.g. SphereMesh is no longer part of the library).

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]'), degree=1))
    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 = UnitCubeMesh(10, 10, 10)
bmesh = BoundaryMesh(mesh, 'exterior')
n = get_facet_normal(bmesh)

plot(n)
interactive()

Comments

Hi,

thanks for your reply.

This does indeed work for a 3D mesh, I should have mentioned my mesh is only 2D.

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

My guess is, that this is some sort of typo, because from the dimensions in the vector I think it should rather be

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

as this is not working for my 2D mesh. I have tried to adjust it to get it working with my 2D mesh, but no success so far.

---

Does this work?

from dolfin import *
import numpy as np

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

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

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

    vec1 = vertices[cells[:, 1]] - vertices[cells[:, 0]]
    normals = vec1[:,[1,0]]*np.array([1,-1])
    normals /= np.sqrt((normals**2).sum(axis=1))[:, np.newaxis]

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

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

    for n in (0,1):
        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 = UnitSquareMesh(10, 10)
bmesh = BoundaryMesh(mesh, 'exterior')
n = get_facet_normal(bmesh)

fid = File('normal.pvd')
fid << n

If you plot the vectors as glyphs in Paraview (since FEniCS does not seem to be able to render it properly), you should see this:

---

Hi,

thanks for that approach and sorry for my late reply. It produces the normal vector just fine, but it seems I still can't put it into a DirichletBC for some reason. I even tried it with "method=pointwise" in the DirichletBC, but still no success.