Compatible projections onto the boundary mesh

Question

Hello,

I have been trying to find a way around the following problem:

Given a scalar field $ \omega $ defined in a 2D mesh, I wish to compute the form \( ( \nabla \times \omega ) \times \vec{n} \) on the boundary mesh, for the purpose of solving a 1D ODE further on.

$\nabla\times\omega\times\vec{n}$ along the boundary will be a DG function. I would expect to be able to reconstruct it exactly using DG basis functions along the boundary.

I have attempted to do this using CG and DG elements, as can be observed in the example code I have included below. Difficulties arise when considering the components of the normal vector $\vec{n}$ on the boundary mesh. I have had to define two expressions which correspond with the components of $\vec{n}$, such that I can compute the tangential component of $\nabla \times \omega$.

However, I have found that errors appear on the corners, probably because I do not distinguish between the different normal vectors of each adjacent element to the corner.

I would really appreciate if somebody could suggest a better way of carrying out this computation, as I am out of ideas at the moment. How can the components of the normal vector to the mesh be included in an integration along the boundary mesh?

Thank you very much.

Gonzalo

from dolfin import *

# Discretization parameters
p = 2
N = 20

# Arbitrary scalar field
w_exp = Expression('x[0]*sin(2*pi*x[0]*x[1])', degree = p)

# Generation of mesh
mesh = UnitSquareMesh(N,N,'left')
n = FacetNormal(mesh)
n_perp = as_vector((n[1], -n[0]))

# Definition of function spaces
W = FunctionSpace(mesh,'CG',p)
P = FunctionSpace(mesh,'DG',p)

# Interpolate scalar field
w = interpolate(w_exp, W)

# Integrate tangential component of curl w along the boundary
curlw_tang = dot(curl(w), n_perp)
integral_1 = dolfin.assemble(curlw_tang*ds)
print 'Integral of curlw_tang along boundary :: ' + str(integral_1)

# Initialize boundary mesh
boundary_mesh = BoundaryMesh(mesh, "exterior")   

# Initialize Discontinous Lagrange elements on boundary
P_b = FunctionSpace(boundary_mesh, 'DG', p)

# Project components of curlw onto P
curlw_x = project(curl(w)[0],P)
curlw_y = project(curl(w)[1],P)

# Project components of curlw onto P_b
curlw_x_b = interpolate(curlw_x, P_b)
curlw_y_b = interpolate(curlw_y, P_b)

# Define components of outwards facing normal vector as dolfin expressions
class n1_exp(Expression):
    def __init__(self, degree = 0):
        self.degree = degree

    def eval(self, value, x):
        if x[0] == 0: # Left boundary
            value[0] = -1
        elif x[0] == 1: # Right boundary
            value[0] = 1
        else:
            value[0] = 0        

class n2_exp(Expression):
    def __init__(self, degree = 0):
        self.degree = degree
    def eval(self, value, x):
        if x[1] == 0: # Lower boundary
            value[0] = -1
        elif x[1] == 1: # Upper boundary
            value[0] = 1
        else:
            value[0] = 0


n1_b = n1_exp(degree = 0)
n2_b = n2_exp(degree = 0)

# Integrate tangential component of w_b (defined in boundary mesh) along the boundary
curlw_tang_b = project(curlw_x_b*n2_b - curlw_y_b*n1_b, P_b)
integral_2 = dolfin.assemble(curlw_tang_b*dx)
print 'Integral of curlw_tang_b along boundary :: ' + str(integral_2)