Combining Dirichlet and Neumann conditions#
Author: Jørgen S. Dokken
Let’s return to the Poisson problem from the Fundamentals chapter and see how to extend the mathematics and the implementation to handle Dirichlet condition in combination with a Neumann condition.
The domain is still the unit square, but now we set the Dirichlet condition
is applied to the remaining sides
The PDE problem#
Let
Again, we choose
For the ease of programming, we define
The variational formulation#
The first task is to derive the variational formulatin. This time we cannot omit the boundary term arising from integration by parts, because
and since
by applying the boundary condition on
Expressing this equation in the standard notation
Implementation#
As in the previous example, we define our mesh,function space and bilinear form
from dolfinx import default_scalar_type
from dolfinx.fem import (Constant, Function, functionspace,
assemble_scalar, dirichletbc, form, locate_dofs_geometrical)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square
from dolfinx.plot import vtk_mesh
from mpi4py import MPI
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad
import numpy as np
import pyvista
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
V = functionspace(mesh, ("Lagrange", 1))
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v)) * dx
Now we get to the Neumann and Dirichlet boundary condition. As previously, we use a Python-function to define the boundary where we should have a Dirichlet condition. Then, with this function, we locate degrees of freedom that fulfill this condition.
def u_exact(x):
return 1 + x[0]**2 + 2 * x[1]**2
def boundary_D(x):
return np.logical_or(np.isclose(x[0], 0), np.isclose(x[0], 1))
dofs_D = locate_dofs_geometrical(V, boundary_D)
u_bc = Function(V)
u_bc.interpolate(u_exact)
bc = dirichletbc(u_bc, dofs_D)
The next step is to define the Neumann condition. We first define UFL
s SpatialCoordinate
-function, and then in turn create a boundary integration measure ds
. As the test function g*v*ds
over the entire boundary.
x = SpatialCoordinate(mesh)
g = -4 * x[1]
f = Constant(mesh, default_scalar_type(-6))
L = f * v * dx - g * v * ds
We can now assemble and solve the linear system of equations
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
V2 = functionspace(mesh, ("Lagrange", 2))
uex = Function(V2)
uex.interpolate(u_exact)
error_L2 = assemble_scalar(form((uh - uex)**2 * dx))
error_L2 = np.sqrt(MPI.COMM_WORLD.allreduce(error_L2, op=MPI.SUM))
u_vertex_values = uh.x.array
uex_1 = Function(V)
uex_1.interpolate(uex)
u_ex_vertex_values = uex_1.x.array
error_max = np.max(np.abs(u_vertex_values - u_ex_vertex_values))
error_max = MPI.COMM_WORLD.allreduce(error_max, op=MPI.MAX)
print(f"Error_L2 : {error_L2:.2e}")
print(f"Error_max : {error_max:.2e}")
Error_L2 : 5.27e-03
Error_max : 2.44e-15
Visualization#
To look at the actual solution, run the script as a python script with off_screen=True
or as a Jupyter notebook with off_screen=False
pyvista.start_xvfb()
pyvista_cells, cell_types, geometry = vtk_mesh(V)
grid = pyvista.UnstructuredGrid(pyvista_cells, cell_types, geometry)
grid.point_data["u"] = uh.x.array
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
plotter.add_text("uh", position="upper_edge", font_size=14, color="black")
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
figure = plotter.screenshot("neumann_dirichlet.png")