Dear all,
I want to solve the Poisson equation on a line d^2 u / dx^2 = f in a line x \in [0,1], where f = - (2 \pi)^2 \cos(2 \pi x). by imposing
- a periodic boundary condition (BC) u(0) = u(1)
- a Dirichlet condition in the middle u(1/2) = -1.
I know the exact solution u_{\rm ex} = \cos (2 \pi x).
I try to solve this in two ways (Options 2 and 3 in the minimal working example code attached). I use a periodic function space and impose:
- a Dirichlet condition in the middle vertex, ds_{\rm m} (option 2)
- a penalty term \alpha/h (u - u_{\rm ex}) \nu_u ds_{\rm m} (option 3)
The solver does not converge in neither of these Options. On the other hand, it converges if I impose, on top of the periodic space, a Dirichlet BC u = 1 on the left vertex (Option 1).
Here is the minimal working example solve.py, which you can run with python3 solve.py. Please note that this code needs the .h5 mesh files, which you can find here. I apologize for the Dropbox link, but I could not find a way to generate a line mesh with an inner vertex on the fly.
from fenics import *
import numpy as np
import ufl
import sys
mesh_path = 'mesh_files/'
mesh = Mesh()
with HDF5File(mesh.mpi_comm(), mesh_path + "line_mesh.h5", "r") as infile:
infile.read(mesh, 'mesh', False)
r_mesh = mesh.hmin()
# read lines and vertices
cf = MeshFunction("size_t", mesh, mesh.topology().dim())
with HDF5File(mesh.mpi_comm(), mesh_path + "line_mesh.h5", "r") as infile:
infile.read(cf, 'cf')
vf = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
with HDF5File(mesh.mpi_comm(), mesh_path + "vertex_mesh.h5", "r") as infile:
infile.read(vf, 'vf')
class PeriodicBoundary(SubDomain):
# Identify the "target domain": the left edge of the line
def inside(self, x, on_boundary):
return near(x[0], 0)
# Map the right edge of the line into the "target domain"
def map(self, x, y):
if near(x[0], 1):
y[0] = 0
else:
# Required: set unmapped points to identity
y[0] = x[0]
periodic_boundary = PeriodicBoundary()
# radius of the smallest cell in the mesh
r_mesh = mesh.hmin()
alpha = 1e2
vertex_l_id = 2
vertex_r_id = 3
vertex_m_id = 4
# measures for the line, for left and right points and for middle point
dx = Measure("dx", domain=mesh, subdomain_data=cf, subdomain_id=1)
ds_l = Measure("ds", domain=mesh, subdomain_data=vf, subdomain_id=vertex_l_id)
ds_r = Measure("ds", domain=mesh, subdomain_data=vf, subdomain_id=vertex_r_id)
ds_m = Measure("dS", domain=mesh, subdomain_data=vf, subdomain_id=vertex_m_id)
ds_lr = Measure("ds", domain=mesh)
Q = FunctionSpace(mesh, 'P', 4, constrained_domain=periodic_boundary)
u = Function(Q)
nu_u = TestFunction(Q)
f = Function(Q)
J_u = TrialFunction(Q)
u_exact = Function(Q)
# exact solution
class u_exact_expression(UserExpression):
def eval(self, values, x):
# test case 1
values[0] = np.cos(2 * np.pi * x[0])
def value_shape(self):
return (1,)
class f_expression(UserExpression):
def eval(self, values, x):
# test case 1
values[0] = -( 2 * np.pi )**2 * np.cos(2 * np.pi * x[0])
def value_shape(self):
return (1,)
u_exact.interpolate(u_exact_expression(element=Q.ufl_element()))
f.interpolate(f_expression(element=Q.ufl_element()))
# Dirichlet boundary conditions on left and right edge
bc_l = DirichletBC(Q, u_exact, vf, vertex_l_id)
bc_r = DirichletBC(Q, u_exact, vf, vertex_r_id )
#Dirichlet condition to enforce u = u _exact at the middle (m) point
bc_m = DirichletBC(Q, u_exact, vf, vertex_m_id )
i = ufl.indices(1)
facet_normal = FacetNormal(mesh)
# variational problem
F_u = (u.dx(i) * nu_u.dx(i) + f * nu_u) * dx \
- facet_normal[i] * (u.dx(i)) * nu_u * ds_lr
F_p = alpha/r_mesh * ( (u('+')+u('-'))/2.0 - (u_exact('+')+u_exact('-'))/2.0 ) * (nu_u('+') + nu_u('-'))/2.0 * ds_m
'''
# Option 1:
# This works
bcs = [bc_l]
F = F_u
'''
# Option 2:
# This does not work
bcs = [bc_m]
F = F_u
'''
# Option 3:
# This does not work
bcs = []
F = F_u + F_p
'''
J = derivative(F, u, J_u)
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()
print(f'\nerror = {assemble((u - u_exact)**2 * dx)/assemble(Constant(1)* dx)}')
Do you know what is going wrong here ?
I guess you will tell me that I should drop the term - facet_normal[i] * (u.dx(i)) * nu_u * ds_lr. I do know that if I fix the value of a field on a boundary, then the test function vanishes. But here I am not fixing the value of u: I am saying that u(0) = u(1). Is there any reason to think that even in this case that term should be dropped ? Why ?
Thanks