I have the following script:
from dolfin import *
from mshr import *
nx, ny = 30, 30
plx, prx = 0.2, 0.8
u0, u1 = 0.0, 0.05
mesh = UnitSquareMesh(nx, ny)
#-------------refinement-start
p, q = Point(0.2, 1.), Point(0.8, 1.)
tol = 0.1
tol2= 0.5
# Selecting edges to refine
class Border(SubDomain):
def inside(self, x, on_boundary):
return (near(x[0], p.x(), tol) and near(x[1], p.y(), tol) or
near(x[0], q.x(), tol) and near(x[1], q.y(), tol)) and near(x[1], 1.0, tol2)
Border = Border()
# Number of refinements
nor = 3
for i in range(nor):
edge_markers = EdgeFunction("bool", mesh)
Border.mark(edge_markers, True)
adapt(mesh, edge_markers)
mesh = mesh.child()
#---------------refinement-end
def plate_boundary(x, on_boundary):
return (plx < x[0] < prx) and near(x[1], 1.0)
def empty_boundary(x, on_boundary):
return on_boundary and not plate_boundary(x, on_boundary)
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)
bcs = [DirichletBC(V, Constant(u1), plate_boundary),
DirichletBC(V, Constant(u0), empty_boundary)]
a = inner(grad(u), grad(v)) * dx
L = Constant(0.0) * v * dx
u = Function(V)
problem = LinearVariationalProblem(a, L, u, bcs)
solver = LinearVariationalSolver(problem)
solver.solve()
plot(u)
plot(mesh)
interactive()
I want to calculate errornorm
of numerical solution and solution that I got by the separation of variables. Analytic solution has emissions in little subdomains near the points of break (points p
and q
) caused by Gibbs phenomenon. So I need to cut out these subdomains for calculating an errornorm
.
How can I perform this?