Hello everyone,
I am playing around with the following problem: Find u satisfying - \Delta u + u = f in \Omega with \nabla u\cdot n = g on \Gamma.
Here \Omega =(0,1)\times (0,1) is the unit square and \Gamma = \partial \Omega is its boundary. The boundary \Gamma is split up in 4 parts (Left, Top, Right, Bottom) and the function g is given on each part.
I am currently aware of 4 different methods to implement this problem (see also the code below):
- Define g with a
UserExpression; - Define g as an
Expressionwith a C++ string; - Use classes to split the boundary and make four expressions for g on each part;
- Implement the boundary parts with
AutoSubDomain.
To test the accuracy of these methods, I prescribe an exact solution u_e and calculate the corresponding functions f_e and g_e such that u_e solves the problem. Then I try to recover u with each method and I compare to the exact solution by means of the L2-norm. I found that Method 1 and 2 yield the same errors. Method 3 and 4 also give the same values.
(A) If the exact function is u(x,y) = x^2y, method 3 (and 4) yield smaller error than method 1 (and 2):
M1-abs 0.006743777927294081
M1-rel 0.026118539443312112
M2-abs 0.006743777927294081
M2-rel 0.026118539443312112
M3-abs 1.563460472663945e-08
M3-rel 6.05525633607536e-08
M4-abs 1.563460472663945e-08
M4-rel 6.05525633607536e-08
(B) In case I use as exact function u(x,y) = \exp\left(-\left(x-\frac{1}{2}\right)^2-\left(y-\frac{1}{2}\right)^2\right) then all methods give the same errors. Note that here g=\nabla u \cdot n is continuous on \Gamma.
M1-abs 1.9954417769876844e-08
M1-rel 2.3321469060676773e-08
M2-abs 1.9954417769876844e-08
M2-rel 2.3321469060676773e-08
M3-abs 1.9954417769876844e-08
M3-rel 2.3321469060676773e-08
M4-abs 1.9954417769876844e-08
M4-rel 2.3321469060676773e-08
Is there a reason why in (A) there is a difference, while in (B) there is none (or is it only the continuity of g that plays a role here)?
Many thanks in advance!
from fenics import *
from dolfin import *
from ufl import nabla_grad
import sympy as sym
###############################################################################
# Problem Statement: Find u such that
#
# -Delta(u) + u = f on Omega
# nabla(u).n = g on Gamma
#
###############################################################################
# Exact solution
x,y = sym.symbols('x[0], x[1]')
u_e_str = 'x**2*y'
#u_e_str = 'sin(x*y)'
#u_e_str = 'exp(-(x-.5)**2-(y-.5)**2)'
#u_e_str = 'exp(-(x-.5)**2 + (y-.5)**2)'
u_e = sym.sympify(u_e_str, locals={'x':x,'y':y})
u_cpp = sym.printing.ccode(u_e)
# Sources and boundary (C++ strings)
f_e = sym.simplify(-(sym.diff(u_e,x,2)+sym.diff(u_e,y,2)) + u_e)
f_cpp = sym.printing.ccode(f_e)
gL_e = sym.simplify(-sym.diff(u_e,x,1))
gL_cpp = sym.printing.ccode(gL_e)
gR_e = sym.simplify(sym.diff(u_e,x,1))
gR_cpp = sym.printing.ccode(gR_e)
gT_e = sym.simplify(sym.diff(u_e,y,1))
gT_cpp = sym.printing.ccode(gT_e)
gB_e = sym.simplify(-sym.diff(u_e,y,1))
gB_cpp = sym.printing.ccode(gB_e)
###############################################################################
# Mesh - Finite Elements - FunctionSpace
mesh = UnitSquareMesh(50,50)
dX = Measure('dx',mesh)
dS = Measure('ds',mesh)
order = 2
el = FiniteElement('CG',mesh.ufl_cell(),order)
V = FunctionSpace(mesh,el)
# Expressions from C++ objects
f = Expression(f_cpp,degree=2)
g_T = Expression(gT_cpp,degree=2)
g_R = Expression(gR_cpp,degree=2)
g_L = Expression(gL_cpp,degree=2)
g_B = Expression(gB_cpp,degree=2)
# Exact solution
exact = interpolate(Expression(u_cpp,degree=2),V)
###############################################################################
# Variational Formulation
u = TrialFunction(V)
phi = TestFunction(V)
a = inner(nabla_grad(u),nabla_grad(phi))*dX + inner(u,phi)*dX
###############################################################################
###############################################################################
# METHOD 1 - Via UserExpression for g
class g_expr(UserExpression):
def __init__(self,**kwargs):
super().__init__(**kwargs)
def eval(self,value,x):
tol = 1E-16
if abs(x[0])<tol: # left
value[0] = eval(gL_cpp)
elif abs(1-x[0])<tol: # right
value[0] = eval(gR_cpp)
elif abs(1-x[1])<tol: # top
value[0] = eval(gT_cpp)
elif abs(x[1])<tol: # bottom
value[0] = eval(gB_cpp)
else:
value[0] = 0.
def value_shape(self):
return()
g = g_expr()
## Variational Formulation
L1 = inner(f,phi)*dx + inner(g,phi)*dS
## Solve
u1 = Function(V)
solve(a == L1, u1)
abs_L2_error1 = errornorm(u1,exact,'L2',degree_rise=3,mesh=mesh)
rel_L2_error1 = abs_L2_error1/norm(exact,'L2',mesh=mesh)
print('M1-abs',abs_L2_error1)
print('M1-rel',rel_L2_error1)
###############################################################################
###############################################################################
# METHOD 2 - Use of C++ in Expression for g
## Set data function g
tol = 1E-16
G = Expression(
'x[0]<tol ?' + gL_cpp + ':' +
'1-x[0]<tol ?' + gR_cpp +':' +
'1-x[1]<tol ?' + gT_cpp + ':' +
'x[1]<tol ?' +gB_cpp + ':' + '.0',
tol=tol,degree=2)
GG = interpolate(G,V)
file = File('G.pvd')
file << GG
## Variational Formulation
L2 = inner(f,phi)*dx + inner(G,phi)*dS
## Solve
u2 = Function(V)
solve(a == L2, u2)
abs_L2_error2 = errornorm(u2,exact,'L2',degree_rise=3,mesh=mesh)
rel_L2_error2 = abs_L2_error2/norm(exact,'L2',mesh=mesh)
print('M2-abs',abs_L2_error2)
print('M2-rel',rel_L2_error2)
###############################################################################
###############################################################################
# METHOD 3 - Boundary via Classes
## define classes + initialize
tol = 1E-16
class Left(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and x[0] < tol
class Top(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and x[1] > 1-tol
class Right(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and x[0] > 1-tol
class Bottom(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and x[1] < tol
BL = Left()
BT = Top()
BR = Right()
BB = Bottom()
## associated marker
marker = MeshFunction('size_t',mesh,mesh.topology().dim()-1)
marker.set_all(0)
BL.mark(marker,0)
BT.mark(marker,1)
BR.mark(marker,2)
BB.mark(marker,3)
## measure
dSc = Measure('ds',domain=mesh,subdomain_data=marker)
## Variatonal formulation
L3 = inner(f,phi)*dX + inner(g_L,phi)*dSc(0) + inner(g_T,phi)*dSc(1) + inner(g_R,phi)*dSc(2) + inner(g_B,phi)*dSc(3)
## Solve
u3 = Function(V)
solve(a == L3, u3)
abs_L2_error3 = errornorm(u3,exact,'L2',degree_rise=3,mesh=mesh)
rel_L2_error3 = abs_L2_error3/norm(exact,'L2',mesh=mesh)
print('M3-abs',abs_L2_error3)
print('M3-rel',rel_L2_error3)
###############################################################################
###############################################################################
# METHOD 4 - AutoSubDomain
## Define boundary
def right(x, on_boundary):
return on_boundary and abs(x[0]-1) < tol
def left(x, on_boundary):
return on_boundary and abs(x[0]) < tol
def top(x, on_boundary):
return on_boundary and abs(x[1]-1) < tol
def bottom(x, on_boundary):
return on_boundary and abs(x[1]) < tol
## Mark the boundary
boundary_marker = MeshFunction('size_t',mesh,mesh.topology().dim()-1)
boundary_marker.set_all(0)
auto_L = AutoSubDomain(left)
auto_R = AutoSubDomain(right)
auto_T = AutoSubDomain(top)
auto_B = AutoSubDomain(bottom)
auto_L.mark(boundary_marker,0)
auto_R.mark(boundary_marker,2)
auto_T.mark(boundary_marker,1)
auto_B.mark(boundary_marker,3)
## Measure
dss = ds(subdomain_data = boundary_marker)
## Variational Formulation
L4 = inner(f,phi)*dx + inner(g_L,phi)*dss(0) + inner(g_T,phi)*dss(1) + inner(g_R,phi)*dss(2) + inner(g_B,phi)*dss(3)
## Solve
u4 = Function(V)
solve(a==L4,u4)
abs_L2_error4 = errornorm(u4,exact,'L2',degree_rise=3,mesh=mesh)
rel_L2_error4 = abs_L2_error4/norm(exact,'L2',mesh=mesh)
print('M4-abs',abs_L2_error4)
print('M4-rel',rel_L2_error4)
###############################################################################