Notebook

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

In [1]:

import dolfinx
import numpy as np
import ufl

from mpi4py import MPI

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

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

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

In [2]:

dofs_L = dolfinx.fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
u_L = dolfinx.Function(V)
u_L.interpolate(lambda x: 1 + 2*x[1]**2)
dolfinx.cpp.la.scatter_forward(u_L.x)
bc_L = dolfinx.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.

In [3]:

dofs_R = dolfinx.fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 1))
u_R = dolfinx.Function(V)
u_R.interpolate(lambda x: 2 + 2*x[1]**2)
dolfinx.cpp.la.scatter_forward(u_R.x)
bc_R = dolfinx.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.

In [4]:

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

V2 = dolfinx.FunctionSpace(mesh, ("CG", 2))
uex = dolfinx.Function(V2)
uex.interpolate(u_exact)
dolfinx.cpp.la.scatter_forward(uex.x)
error_L2 = dolfinx.fem.assemble_scalar((uh - uex)**2 * ufl.dx)
error_L2 = np.sqrt(MPI.COMM_WORLD.allreduce(error_L2, op=MPI.SUM))

u_vertex_values = uh.compute_point_values()
u_ex_vertex_values = uex.compute_point_values()
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 : 6.22e-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.

In [5]:

import pyvista
# Start virtual framebuffer
pyvista.start_xvfb(wait=0.0)

import dolfinx.plot
pyvista_cells, cell_types = dolfinx.plot.create_vtk_topology(mesh, mesh.topology.dim)
grid = pyvista.UnstructuredGrid(pyvista_cells, cell_types, mesh.geometry.x)
point_values = uh.compute_point_values()
if np.iscomplexobj(point_values):
    point_values = point_values.real
grid.point_arrays["u"] = point_values
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()
figure = plotter.screenshot("multiple_dirichlet.png")

In [ ]: