I am trying to implement periodic boundaries on curved boundaries. As a test I modified the periodic Poisson example implementing the constraints from the periodic boundaries as null-space. The problem: I obtain nearly the same solution as in the example but the values u(0,y)=u(1,y) at the boundaries differ slightly from the ones obtained in the demo! They should be exactly the same! I have verified by using another FEM software that the solution of the demo is correct (unless that software has the same problem ;-)
I have checked the null-space. As far as I can see it is correct. So what is the problem?
from dolfin import *
import numpy as N
parameters['allow_extrapolation'] = True
class Constraint:
"""
Constraint implements a tie between the values at two points p1 and p2.
Example:
Create a tie between the values at (0.0, 0.5) and (1.0, 0.5)
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'CG', 2)
constraint = Constraint(V)
tie = constraint.vector(Point(0.0, 0.5), Point(1.0, 0.5))
The constraint equation is given by tie.inner(u.vector()) == 0,
i.e. all ties span the nullspace of the linear equation.
"""
def __init__(self, V):
self.V = V
self.mesh = V.mesh()
self.dofmap = V.dofmap()
self.finite_element = V.element()
def evaluate_basis(self, p):
bbt = mesh.bounding_box_tree()
id = bbt.compute_first_entity_collision(p)
if id >= mesh.num_cells():
id = bbt.compute_closest_entity(p)[0]
c = Cell(self.mesh, id)
vc = c.get_vertex_coordinates()
dofs = self.dofmap.cell_dofs(id)
no_basis_fns = self.finite_element.space_dimension()
value_dimension = self.finite_element.value_dimension(0)
basis = N.zeros((no_basis_fns, value_dimension), dtype=N.float64)
coords = N.zeros(3, dtype=N.float64)
coords[0], coords[1], coords[2] = p.x(), p.y(), p.z()
self.finite_element.evaluate_basis_all(basis, coords, vc, 0)
u = Function(self.V)
v = u.vector()
# fixme: implement mixed spaces
for k in range(value_dimension):
for j in range(no_basis_fns):
l = no_basis_fns*(k-1)+j
v[dofs[l]] = basis[j][k]
return v
def vector(self, p1, p2):
v = (self.evaluate_basis(p1)-self.evaluate_basis(p2))
v /= v.norm('l2')
return v
class PeriodicBoundaryCondition(SubDomain):
def __init__(self):
SubDomain.__init__(self)
pass
def inside(self, x, on_boundary):
return near(x[0], 0.0) and on_boundary
def map(self, x, y):
y[0] = x[0]-1.0
y[1] = x[1]
def nullspace(self, V):
r = Constraint(V)
dofmap = V.dofmap()
coords = dofmap.tabulate_all_coordinates(V.mesh())
s = []
x = [0, 0]
for x[0], x[1] in zip(coords[0::2], coords[1::2]):
y = [0, 0]
self.map(x, y)
if self.inside(y, True):
A = Point(x[0], x[1])
B = Point(y[0], y[1])
s.append(r.vector(A, B))
return VectorSpaceBasis(s)
class DirichletBoundary(SubDomain):
def inside(self, x, on_boundary):
return bool((near(x[1], 0.0) or near(x[1], 1.0)) and on_boundary)
class Source(Expression):
def eval(self, values, x):
dx = x[0] - 0.5
dy = x[1] - 0.5
values[0] = x[0]*sin(5.0*pi*x[1]) + 1.0*exp(-(dx*dx + dy*dy)/0.02)
mesh = UnitSquareMesh(32, 32)
mesh.init()
V = FunctionSpace(mesh, 'CG', 1)
pbc = PeriodicBoundaryCondition()
nullspace = pbc.nullspace(V)
f = Source()
u, v = TrialFunction(V), TestFunction(V)
a = inner(grad(u), grad(v))*dx
L = f*v*dx
# Create Dirichlet boundary condition
zero = Constant(0.0)
db = DirichletBoundary()
dbc = DirichletBC(V, zero, db)
# Compute solution
u = Function(V)
U = u.vector()
A, b = assemble_system(a, L, [dbc])
prm = parameters['krylov_solver']
prm['relative_tolerance'] = 1.0e-16
prm['absolute_tolerance'] = 1.0e-16
nullspace.orthogonalize(b)
solver = KrylovSolver()
solver.set_nullspace(nullspace)
solver.solve(A, U, b)
print "dimension of nullspace", nullspace.dim()
# Save solution to file
file = File("periodic.pvd")
file << u
# Plot solution
plot(u, interactive=True)
The following plot shows the discrepancy:
