Calculation of the functions containing the normal-vector at the boundary

Question

ScalarSpace = FunctionSpace(mesh,'CG',1)
VectorSpace = VectorFunctionSpace(mesh,'CG',1)

N = FacetNormal(mesh)

v_exp = Expression(('0.0','x[0]'),degree=2, domain=mesh)
v = interpolate(v_exp, VectorSpace)


SomeFunction_exp = Expression(('0.3'),degree=1, domain=mesh)
SomeFunction = project(SomeFunction_exp, ScalarSpace)

First question:

I cannot use project(1/dot(v,N),ScalarSpace) for calculating 1/dot(v,N). Is there a possibility to do this without using assemble like

assemble(1/dot(v,N)*TestFunction(ScalarSpace)*ds)

and then using the inverse massmatrix to get an approximation of 1/dot(v,N) at the boundary?

Second Question:

Using

assemble(SomeFunction/dot(v,N)*TestFunction(ScalarSpace)*ds)        ( Eq.1)

in the upper example i get reasonable results. But if I do something like this

assemble(SomeFunction*SomeFunction/dot(v,N)*TestFunction(ScalarSpace)*ds)        (Eq. 2)

This calculates very high numbers in some nodes. I dont understand it. 1/dot(v,N) can equal zero or near zero, but i handle that case. And even then its not clear, why (Eq. 1) is working, while (Eq. 2) isn't.

Answer 1

Hi, one way is to represent the normal part of v in DG0 space. Consider

from dolfin import *
from mshr import *
import sys

# mesh = UnitSquareMesh(*[int(sys.argv[1])]*2)
mesh = generate_mesh(Circle(Point(0, 0), 1), int(sys.argv[1]))

V = VectorFunctionSpace(mesh, 'CG', 1)
v = interpolate(Expression(('x[0]', 'x[1]'), degree=1), V)

Q = FunctionSpace(mesh, 'DG', 0)
q = TestFunction(Q)
n = FacetNormal(mesh)
# Projection step accounting for the mass matrix inverse
# TrialFunction(Q)*q*ds -> FacetArea(mesh)
h = FacetArea(mesh)
v_dot_n = assemble(dot(v, n)*q/h*ds)
# As function
v_dot_n = Function(Q, v_dot_n)

# Check: \int_{\Omega} div(v)*dx = \int_{\partial\Omega} v.n*ds
value0 = assemble(div(v)*dx)
value = assemble(v_dot_n*ds)

info('Error %g' % abs(value0-value)) 

Comments

Thank you! That helped me a lot.