Setting multiple Dirichlet condition

In the previous section, we used a single function to \(u_d\) to setting Dirichlet conditions on two parts of the boundary. However, it is often more practical to use multiple functins, one for each subdomain of the boundary. We consider a similar example to the previous example and redefine it consist of two Dirichlet boundary conditions

\[ -\nabla^2 u =f \quad \text{in } \Omega, \]

\[ u=u_L \quad \text{on } \Lambda_D^L \]

\[ u=u_R \quad \text{on } \Lambda_D^R \]

\[ -\frac{\partial u}{\partial n} = g \quad \text{on } \Lambda_N. \]

Here, \(\Lambda_D^L\) is the left boundary \(x=0\), while \(\Lambda_D^R\) is the right boundary \(x=1\). We note that \(u_L(y)=1+2y^2\), \(u_R(y)=2+2y^2\) and \(g(y)=-4y\) using the same analytical example as in the previous section.

We start by defining the mesh, function space and variational formulation as in the previous exercise

import numpy as np
import pyvista

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 mpi4py import MPI
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, dx, ds, grad
from petsc4py.PETSc import ScalarType

def u_exact(x):
    return 1 + x[0]**2 + 2*x[1]**2

mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
V = FunctionSpace(mesh, ("CG", 1))
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
x = SpatialCoordinate(mesh)
g = - 4 * x[1]
f = Constant(mesh, ScalarType(-6))
L = f * v * dx - g * v * ds

We next mark the two boundaries separately, starting with the left boundary

dofs_L = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
u_L = Function(V)
u_L.interpolate(lambda x: 1 + 2*x[1]**2)
bc_L = dirichletbc(u_L, dofs_L)

Note that we have used lambda-functions to compactly define the functions returning the subdomain evaluation and function evaluation. We can use a similar procedure for the right boundary condition, and gather both boundary conditions in a vector bcs.

dofs_R = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 1))
u_R = Function(V)
u_R.interpolate(lambda x: 2 + 2*x[1]**2)
bc_R = dirichletbc(u_R, dofs_R)
bcs = [bc_R, bc_L]

We are now ready to again solve the problem, and check the \(L^2\) and max error at the mesh vertices.

problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

V2 = FunctionSpace(mesh, ("CG", 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 visualize the solution, run the script with in a Jupyter notebook with off_screen=False or as a python script with off_screen=True.

pyvista.start_xvfb()
from dolfinx.plot import create_vtk_mesh
pyvista_cells, cell_types, geometry = create_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("multiple_dirichlet.png")