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
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)
a = dot(grad(u), grad(v)) * dx
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