Normals on the boundary mesh not pointing in the same direction

1

Good day

I have been running 2d DG code that determines the normals on the boundary mesh.
what I have noticed and can not understand is that the normals are pointing in opposite directions. please see the code below and the output.

from dolfin import *
from mshr import *
import numpy as np
import matplotlib.pyplot as pyplot

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 range(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)
channel = Rectangle(Point(-2.0, -2.0),Point(5.0, 2.0))
cylinder = Circle(Point(-0.5, 0.5), 0.5)
domain = channel - cylinder
mesh = generate_mesh(domain,64)
bmesh = BoundaryMesh(mesh, ‘exterior’)
n = get_facet_normal(bmesh)
plot(n)
pyplot.savefig(‘normals.png’)
pyplot.show()