Discontinuous Data in Neumann Boundary Condition

variational formulation
1

Hello,

if I have discontinuous Neumann boundary data, but I guess that FEniCS interpolates these data with a Courant function (that is continuous). This leads to difficulties, for example, the solution to the following problem is u(x,y) = x and should be computed exactly. But the algorithm does not lead to the exact solution. Can I circumvent this issue? Best regards, Johannes

from dolfin import *

mesh = UnitSquareMesh(4,4)
f = Expression(‘x[0]’,degree = 1)
u_exact = Expression(‘x[0]’,degree = 1)
g = Expression(’-(x[0]<=0.000000000001)+(0.99999999999999<=x[0])’,degree = 0)
V = FunctionSpace(mesh,‘P’,1)
u = TrialFunction(V)
v = TestFunction(V)

a = (dot(grad(u),grad(v)) + uv)dx
L = fvdx + gvds
u_h = Function(V)
solve(a == L, u_h)

error = assemble((u_exact-u_h)**2*dx)
print(error)

2

You can directly use spatial coordinates in UFL as follows:

#f = Expression('x[0]',degree = 1)
#u_exact = Expression('x[0]',degree = 1)
#g = Expression('-(x[0]<=0.000000000001)+(0.99999999999999<=x[0])',degree = 0)
x = SpatialCoordinate(mesh)
f = x[0]
u_exact = x[0]
g = conditional(lt(x[0],DOLFIN_EPS),-1.0,
                conditional(gt(x[0],1.0-DOLFIN_EPS),1.0,Constant(0.0)))
3

Hi,

this works fine, thank you.