Hi, I am trying to use some function 'u_e' as an exact solution to a DG problem over a unit square mesh, that is divided into two subdomains, in these subdomains, the coefficient "aa" takes the values [1, 100] respectively. I used "u_e" to find the formulas for the source function and the boundary values.
I am not receiving any error massages, but the computed Error seems to be a bet big! which is not expected! I couldn't find why the error is big, can any one help me?
Thank you,
here is an example
from dolfin import *
import numpy as np
mesh = UnitSquareMesh(10, 10, 'crossed')
class Omega1(SubDomain):
def inside(self,x, on_boundary):
return (between(x[0], (0.0, 1.0)) and between(x[1], (x[0], 1.-x[0])))
class Omega2(SubDomain):
def inside(self,x, on_boundary):
return (between(x[1], (0.0, 1.0)) and between(x[0], (x[1], 1.-x[1])))
mesh.init()
# Marking the domains
subdomains = CellFunction("size_t",mesh)
subdomain1 = Omega1()
subdomain1.mark(subdomains, 0)
subdomain2 = Omega2()
subdomain2.mark(subdomains, 1)
#values of aa
a1 = 1.0
a2 = 100.0
V0 = FunctionSpace(mesh, 'DG', 0)
aa = Function(V0)
aa_values = [a1, a2] #values of a in each subdomain
#the following is to make sure a takes the right value in each subdomain:
for cell_no in range(len(subdomains.array())):
subdomain_no = subdomains.array()[cell_no]
aa.vector()[cell_no] = aa_values[subdomain_no]
help = np.asarray(subdomains.array(), dtype=np.int32)
aa.vector()[:] = np.choose(help, aa_values)
# plot (subdomains, interactive = True)
# Exact solution and Source
u_e = Expression('(x[1]-x[0])*(1-x[0]-x[1])')
f = Expression('(x[1]-x[0])*(1-x[0]-x[1])')
# boundaries
class DirichletBoundary_top(SubDomain):
def inside(self ,x, on_boundary):
return x[1] > 1.0 - DOLFIN_EPS
class DirichletBoundary_bottom(SubDomain):
def inside(self ,x, on_boundary):
return x[1] < DOLFIN_EPS
class DirichletBoundary_left(SubDomain):
def inside(self, x, ob_boundary):
return x[0] < DOLFIN_EPS
class DirichletBoundary_right(SubDomain):
def inside(self, x, ob_boundary):
return x[0] > 1.0 - DOLFIN_EPS
boundaries = FacetFunction('size_t', mesh)
Gamma_D1 = DirichletBoundary_top()
Gamma_D1.mark(boundaries, 2)
Gamma_D2 = DirichletBoundary_bottom()
Gamma_D2.mark(boundaries, 4)
Gamma_N1 = DirichletBoundary_left()
Gamma_N1.mark(boundaries, 1)
Gamma_N2 = DirichletBoundary_right()
Gamma_N2.mark(boundaries, 3)
# boundary values
u_top = Expression('-x[0]*(1-x[0])')
u_bottom = Expression('-x[0]*(1-x[0])')
u_left = Expression('x[1]*(1-x[1])')
u_right = Expression('x[1]*(1-x[1])')
# solve the problem
V = FunctionSpace(mesh, "DG", 3)
u = TrialFunction(V)
v = TestFunction(V)
dx = Measure("dx")[subdomains]
ds = Measure("ds")[boundaries]
alpha = Constant(40)
h_T = CellSize(mesh)
h_T_avg = (h_T('+') + h_T('-'))/2.0
n = FacetNormal(mesh)
a = dot(a1*grad(u),grad(v))*dx(0) + dot(a2*grad(u),grad(v))*dx(1)\
+ u*v*dx()- dot(avg(aa*grad(u)),jump(v,n))*dS\
- dot(avg(aa*grad(v)),jump(u,n))*dS\
- dot(aa*grad(u),v*n)*ds()\
- dot(avg(aa*grad(v)),jump(u,n))*dS\
- dot(aa*grad(v),u*n)*ds()\
+ alpha('+')/h_T_avg*dot(jump(u,n),jump(v,n))*dS\
+ alpha/h_T*u*v*ds()
L = f*v*dx() - u_left*dot(aa*grad(v),n)*ds(1) - u_top*dot(aa*grad(v),n)*ds(2)\
- u_right*dot(aa*grad(v),n)*ds(3)- u_bottom*dot(aa*grad(v),n)*ds(4) \
+ alpha/h_T*u_left*v*ds(1) + alpha/h_T*u_top*v*ds(2) \
+ alpha/h_T*u_right*v*ds(3) + alpha/h_T*u_bottom*v*ds(4)
bcs = [ DirichletBC(V, u_left, Gamma_N1),DirichletBC(V, u_top, Gamma_D1),
DirichletBC(V, u_right, Gamma_N2), DirichletBC(V, u_bottom, Gamma_D2)]
u = Function(V)
A, b = assemble_system(a, L,bcs, mesh= mesh, cell_domains = subdomains,
exterior_facet_domains = boundaries)
solve(A, u.vector(), b)
# compute the error
u_e_V = interpolate(u_e, V)
error1 = (u - u_e_V)**2*dx
E1 = sqrt(assemble(error1))
# H1 seminorm
Ve = FunctionSpace(mesh, 'Lagrange', 5)
u_Ve = interpolate(u, Ve)
u_e_Ve = interpolate(u_e, Ve)
e_Ve = Function(Ve)
e_Ve.vector()[:] = u_e_Ve.vector().array() - u_Ve.vector().array()
error2 = inner(aa*grad(e_Ve), grad(e_Ve))*dx()
E2 = sqrt(assemble(error2))
print 'E1:', E1
print 'E2:', E2
plot(u_e_V,title="u_e", interactive=True)
plot(u,title="u_h", interactive=True)
