# Setting multiple Dirichlet condition#

In the previous section, we used a single function for $$u_D$$ to set Dirichlet conditions on two parts of the boundary. However, it is often more practical to use multiple functions, one for each subdomain of the boundary. We consider a similar example to the previous example and redefine it to 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

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, dx, ds, grad

import numpy as np
import pyvista

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, ("Lagrange", 1))
u = TrialFunction(V)
v = TestFunction(V)
x = SpatialCoordinate(mesh)
g = - 4 * x[1]
f = Constant(mesh, default_scalar_type(-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, ("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 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()
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()