Domain decomposition solvers#
This demo illustrates how to:
Assemble an operator in PETSc’s unassembled
MATISformat.Solve it with a Balancing Domain Decomposition by Constraints (BDDC) preconditioner.
Set up a mesh so that each MPI process owns one non-overlapping subdomain.
Attach the rigid body modes that BDDC needs to build an effective coarse space for elasticity.
Solve the Poisson and elasticity problems with PCHPDDM, an algebraic overlapping Schwarz method whose coarse space is computed from local eigenproblems.
Precondition a curl-curl problem in \(H({\rm curl})\) on a cube, giving BDDC the discrete gradient so that the near-kernel of the curl goes into its coarse space.
Download sources
Equations and problem definitions#
Three problems are solved, each by its own function, so that the parts specific to a problem are separated from the domain decomposition machinery they share. The first two share a mesh of the unit square; the third is posed in 3D and has its own.
The first is the Poisson equation with homogeneous Dirichlet conditions on the whole boundary,
with \(f = 2 \pi^{2} \sin(\pi x_{0}) \sin(\pi x_{1})\), for which the exact solution is \(u = \sin(\pi x_{0}) \sin(\pi x_{1})\). The variational problem is: find \(u \in V\) such that
with \(V\) the space of piecewise linear Lagrange functions.
The second is linearised elasticity for a displacement \(u\), clamped on \(x_{0} = 0\) and loaded by its own weight,
with \(\sigma(u) = 2 \mu \, \epsilon(u) + \lambda \, {\rm tr}(\epsilon(u)) I\) and \(\epsilon(u) = (\nabla u + (\nabla u)^{T}) / 2\). Its variational problem is: find \(u \in W\) such that
with \(W\) the space of vector-valued piecewise linear Lagrange functions.
Domain decomposition and the unassembled operator#
BDDC is a non-overlapping domain decomposition method. Each MPI process owns one subdomain, and the preconditioner is built from local solves on those subdomains plus a small coarse problem that couples them. It therefore needs the operator in unassembled form,
where \(A_{i}\) is the operator assembled on subdomain \(i\) alone and
\(R_{i}\) restricts from the global to the local degrees of freedom.
PETSc stores operators this way as MATIS, which DOLFINx creates
when the matrix kind is "is". Nothing is summed across processes:
a matrix entry belonging to a degree of freedom on a subdomain
interface is the sum of the contributions held by each process that
shares it.
Two consequences shape the demo.
The mesh is created with GhostMode.none, so that the degrees of
freedom a process holds are exactly those of the cells it owns. With
ghost cells, a process would also hold degrees of freedom that
receive no contribution from its own cells, leaving empty rows in
\(A_{i}\) and making the subdomain solves singular.
Dirichlet degrees of freedom on a subdomain interface are shared. Wherever an interface meets the constrained boundary, the degrees of freedom on the entities they have in common are held by every process meeting there: the vertices for the \(P_{1}\) problems below, and the edges and faces for the Nedelec one. Each of those processes must place a value on its local diagonal, or \(A_{i}\) is again singular — but those values are summed, so each contributes a share rather than the whole value. The demo reports how many such degrees of freedom there are.
Implementation#
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
from dolfinx import common, fem, mesh
from dolfinx.fem.petsc import (
apply_lifting,
assemble_matrix,
assemble_vector,
discrete_gradient,
set_bc,
)
dtype = PETSc.ScalarType
xdtype = PETSc.RealType
With CG, BDDC solves the subdomain interior problems once up front
and assumes the interior residual then stays zero. In single
precision rounding breaks this and CG breaks down, so both BDDC
solves set pc_bddc_switch_static there, which redoes the interior
correction in every application at the cost of an extra local solve.
Every subdomain and coarse problem below is factorised on a single process, and all of them are symmetric positive definite. PETSc’s own factorisation is always available, but an external package is usually faster, so the first one PETSc has is used instead, paired with Cholesky where it offers one. SuperLU_DIST’s Cholesky is an LU with a symmetric ordering, and sequential SuperLU has no Cholesky at all, so both rank below SuiteSparse, whose Cholesky solver is named CHOLMOD.
factor_pc_type, mat_solver_type = "cholesky", "petsc"
for pc, solver, package in (
("cholesky", "mumps", "mumps"),
("cholesky", "cholmod", "suitesparse"),
("lu", "superlu_dist", "superlu_dist"),
("lu", "superlu", "superlu"),
):
if PETSc.Sys.hasExternalPackage(package):
factor_pc_type, mat_solver_type = pc, solver
break
if MPI.COMM_WORLD.rank == 0:
print(f"Direct solver for the local problems: {mat_solver_type}")
def direct_solver(prefix: str) -> PCOptions:
"""PETSc options pointing a sub-solver at the selected direct solver.
Args:
prefix: Option prefix of the sub-solver to configure.
Returns:
Factorisation type and solver package for ``prefix``.
"""
return {
f"{prefix}_pc_type": factor_pc_type,
f"{prefix}_pc_factor_mat_solver_type": mat_solver_type,
}
Each problem reports how many of its constrained degrees of freedom
lie on a subdomain interface, which
:func:~dolfinx.common.num_sharing_ranks answers directly.
def norm_L2(v, quadrature_degree: int = 6) -> float:
"""L2 norm of a UFL expression over the whole mesh.
The quadrature degree is set rather than estimated, because the
integrands are not all polynomial: the Poisson error measures a
piecewise linear solution against a product of sines, and the
curl-curl norm squares a second-order Nedelec function. Degree 6
agrees with every higher degree tried.
Args:
v: Expression to measure.
quadrature_degree: Degree of the quadrature rule.
Returns:
The L2 norm over the whole mesh.
"""
dx = ufl.dx(metadata={"quadrature_degree": quadrature_degree})
form = fem.form(ufl.inner(v, v) * dx, dtype=dtype)
return np.sqrt(form.mesh.comm.allreduce(fem.assemble_scalar(form), MPI.SUM).real)
def num_shared_dofs(V: fem.FunctionSpace, dofs: np.ndarray) -> int:
"""Count the entries of ``dofs`` held by more than one process.
Args:
V: Space the degrees of freedom belong to.
dofs: Degrees of freedom to test, as local block indices.
Returns:
How many of ``dofs`` this process shares with another.
"""
sharers = common.num_sharing_ranks(V.dofmap.index_map, dofs)
return int(np.count_nonzero(sharers > 1))
All three problems are assembled and solved the same way: build the
operator in the format the preconditioner needs, carry the boundary
conditions over to the right-hand side, and run CG under an options
prefix of the solver’s own. Only prepare differs, and only where a
preconditioner needs more than the operator itself.
def solve_cg(
V: fem.FunctionSpace,
a: ufl.Form,
L: ufl.Form,
bcs: list,
dofs: np.ndarray,
kind: str | None,
pc_options: PCOptions,
name: str,
prepare=None,
) -> tuple[fem.Function, int, int]:
"""Assemble a problem and solve it by CG.
Args:
V: Space the solution lives in.
a: Bilinear form.
L: Linear form.
bcs: Dirichlet boundary conditions.
dofs: Constrained degrees-of-freedom, for the shared count.
kind: PETSc matrix kind, ``"is"`` for ``MATIS`` or ``None``
for the default assembled matrix.
pc_options: PETSc options that select and configure the
preconditioner, without a prefix.
name: Names this solver's options prefix.
prepare: Called with the matrix and the solver once both
exist and before the solve, for whatever the
preconditioner needs beyond the operator.
Returns:
The solution, the number of Krylov iterations, and the number
of constrained degrees-of-freedom shared with another process.
"""
a_form, L_form = fem.form(a, dtype=dtype), fem.form(L, dtype=dtype)
A = assemble_matrix(a_form, bcs=bcs, kind=kind)
A.assemble()
A.setOption(PETSc.Mat.Option.SPD, True) # type: ignore[arg-type]
b = assemble_vector(L_form)
apply_lifting(b, [a_form], bcs=[bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) # type: ignore[arg-type]
set_bc(b, bcs)
# Everything, Krylov method included, goes through the options
# database under a prefix of this solver's own, so that nothing is
# left behind for the next solve to inherit
pc_type = pc_options["pc_type"]
options: PCOptions = {
"ksp_type": "cg",
"ksp_rtol": rtol,
"ksp_max_it": 100,
"ksp_error_if_not_converged": True,
**pc_options,
}
prefix = f"demo_dd_{name}_{pc_type}_{V.dofmap.index_map.size_global}_"
ksp = PETSc.KSP().create(V.mesh.comm) # type: ignore[arg-type]
ksp.setOperators(A)
ksp.setOptionsPrefix(prefix)
opts = PETSc.Options(prefix)
for key, value in options.items():
opts[key] = value # type: ignore[index]
ksp.setFromOptions()
if prepare is not None:
prepare(A, ksp)
uh = fem.Function(V, name="u", dtype=dtype)
ksp.solve(b, uh.x.petsc_vec)
# Cleared after the solve, not before: PCHPDDM reads some of its
# options when it is first applied rather than at setFromOptions
for key in options:
del opts[key] # type: ignore[arg-type]
uh.x.scatter_forward()
its = ksp.getIterationNumber()
num_shared = num_shared_dofs(V, dofs)
ksp.destroy()
A.destroy()
b.destroy()
return uh, its, num_shared
solve_poisson builds the Poisson problem on a given mesh and solves
it with the preconditioner it is given, on a matrix of the kind that
preconditioner needs: MATIS for BDDC, the default assembled matrix
for PCHPDDM.
def solve_poisson(
msh: mesh.Mesh, kind: str | None, pc_options: PCOptions
) -> tuple[fem.Function, int, int]:
"""Solve the Poisson problem on ``msh`` using CG.
Args:
msh: Mesh, which must have been built without ghost cells.
kind: PETSc matrix kind, ``"is"`` for ``MATIS`` or ``None``
for the default assembled matrix.
pc_options: PETSc options that select and configure the
preconditioner, without a prefix.
Returns:
The solution, the number of Krylov iterations, and the number
of constrained degrees of freedom shared with another process.
"""
V = fem.functionspace(msh, ("Lagrange", 1))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
f = 2 * ufl.pi**2 * ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1])
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = ufl.inner(f, v) * ufl.dx
tdim = msh.topology.dim
msh.topology.create_connectivity(tdim - 1, tdim)
dofs = fem.locate_dofs_topological(V, tdim - 1, mesh.exterior_facet_indices(msh.topology))
bcs = [fem.dirichletbc(dtype(0), dofs, V)] # type: ignore[operator]
return solve_cg(V, a, L, bcs, dofs, kind, pc_options, "poisson")
Elasticity needs more than this. The rigid body modes — the displacements that cost no energy — span the kernel of the operator on a subdomain that the Dirichlet condition does not touch, and at any useful number of processes some subdomain is always interior. Poisson has the same difficulty in principle, but its kernel is one dimensional (the constants) and a single primal vertex removes it. In elasticity the kernel is three dimensional in 2D, and the vertices alone leave the rotation, so the subdomain solves stay singular.
Two things fix it. PETSc reads the modes from the matrix as its
near null space and uses them to build the coarse space. And
pc_bddc_use_change_of_basis turns the edge constraints into
explicit primal degrees of freedom, which are then removed from the
subdomain problems; without it those problems keep their rigid body
modes and the factorisation fails.
BDDC factorises three problems of its own: the subdomain interior
(dirichlet) and correction (neumann) problems, and the coarse
problem, which it solves redundantly on every process. Each takes the
direct solver selected above, in place of BDDC’s default LU.
bddc_options: PCOptions = {
"pc_type": "bddc",
"pc_bddc_use_change_of_basis": True,
"pc_bddc_switch_static": single_precision,
**direct_solver("pc_bddc_dirichlet"),
**direct_solver("pc_bddc_neumann"),
**direct_solver("pc_bddc_coarse_redundant"),
}
PCHPDDM, PETSc’s interface to the HPDDM library, takes a different route. It is an overlapping Schwarz method that works on the ordinary assembled matrix: the rows each process owns are extended into overlapping subdomains, and the preconditioner solves on every one of them. Its coarse space comes from local eigenproblems, in the spirit of GenEO: the eigenvectors with the smallest eigenvalues, which include the rigid body modes of a subdomain the Dirichlet condition does not touch, span the coarse space. No near null space or primal constraints are needed.
The options select:
pc_hpddm_harmonic_overlap: build the eigenproblems algebraically, on subdomains extended by one layer. It is required here: with only an assembled matrix, PCHPDDM otherwise has no local operator to build the eigenproblems from and falls back to one-level Schwarz.pc_hpddm_levels_1_eps_threshold_relative: keep the eigenvectors whose eigenvalues fall below a relative threshold, which is what sizes the coarse space. How many to compute is left to PCHPDDM.pc_hpddm_levels_1_st_pc_typeandpc_hpddm_levels_1_eps_pc_type: the factorisations inside the eigensolver.pc_hpddm_levels_1_pc_typeandpc_hpddm_levels_1_pc_asm_overlap: additive Schwarz on subdomains overlapping by two layers, solved by direct factorisation (pc_hpddm_levels_1_sub_pc_type).pc_hpddm_define_subdomainsis off so that the Schwarz method builds these subdomains itself, rather than reusing those of the eigenproblems.pc_hpddm_levels_1_pc_asm_typeandpc_hpddm_coarse_correction: symmetric variants of the Schwarz method and of the coarse correction, which CG requires.the
_pc_factor_mat_solver_typeentries: the direct solver behind each of those factorisations.
PCHPDDM is available only when PETSc is configured with HPDDM and SLEPc.
hpddm_options: PCOptions = {
"pc_type": "hpddm",
"pc_hpddm_harmonic_overlap": 1,
"pc_hpddm_levels_1_eps_threshold_relative": 100,
**direct_solver("pc_hpddm_levels_1_st"),
**direct_solver("pc_hpddm_levels_1_eps"),
"pc_hpddm_define_subdomains": False,
"pc_hpddm_levels_1_pc_type": "asm",
"pc_hpddm_levels_1_pc_asm_overlap": 2,
**direct_solver("pc_hpddm_levels_1_sub"),
"pc_hpddm_levels_1_pc_asm_type": "basic",
"pc_hpddm_coarse_correction": "balanced",
}
def rigid_body_modes(V: fem.FunctionSpace) -> PETSc.NullSpace:
"""Build the rigid body modes of a displacement space.
Args:
V: Vector-valued displacement space.
Returns:
The translations and rotations, which PETSc builds from the
coordinates of the owned degrees of freedom.
"""
gdim = V.mesh.geometry.dim
num_owned = V.dofmap.index_map.size_local
x = V.tabulate_dof_coordinates()[:num_owned, :gdim].copy()
coords = PETSc.Vec().createWithArray(x.ravel(), bsize=gdim, comm=V.mesh.comm) # type: ignore[arg-type]
modes = PETSc.NullSpace().createRigidBody(coords)
coords.destroy()
return modes
solve_elasticity is assembled and solved like the Poisson problem,
with a prepare step that hands BDDC the rigid body modes.
def solve_elasticity(
msh: mesh.Mesh, kind: str | None, pc_options: PCOptions
) -> tuple[fem.Function, int, int]:
"""Solve the elasticity problem on ``msh`` using CG.
Args:
msh: Mesh, which must have been built without ghost cells.
kind: PETSc matrix kind, ``"is"`` for ``MATIS`` or ``None``
for the default assembled matrix.
pc_options: PETSc options that select and configure the
preconditioner, without a prefix.
Returns:
The displacement, the number of Krylov iterations, and the
number of constrained degrees of freedom shared with another
process.
"""
gdim = msh.geometry.dim
V = fem.functionspace(msh, ("Lagrange", 1, (gdim,)))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
# Non-dimensionalised, so that the matrix entries and the diagonal
# written on constrained rows are of the same order
ν = 0.3
μ = 1.0 / (2.0 * (1.0 + ν))
λ = ν / ((1.0 + ν) * (1.0 - 2.0 * ν))
def σ(w):
"""Return an expression for the stress σ given a displacement w."""
return 2.0 * μ * ufl.sym(ufl.grad(w)) + λ * ufl.div(w) * ufl.Identity(gdim)
f = ufl.as_vector((0.0, -1.0e-2))
a = fem.form(ufl.inner(σ(u), ufl.sym(ufl.grad(v))) * ufl.dx, dtype=dtype)
L = fem.form(ufl.inner(f, v) * ufl.dx, dtype=dtype)
# Clamp the x0 = 0 edge
tdim = msh.topology.dim
facets = mesh.locate_entities_boundary(msh, tdim - 1, lambda x: np.isclose(x[0], 0.0))
dofs = fem.locate_dofs_topological(V, tdim - 1, facets)
bc = fem.dirichletbc(np.zeros(gdim, dtype=dtype), dofs, V)
def prepare(A, _ksp) -> None:
"""Give BDDC the rigid body modes.
PCHPDDM recovers them itself, from the eigenproblems it solves
on each subdomain, so they are neither built nor attached for
it. PETSc takes its own reference on attachment, so this side
keeps none.
"""
if pc_options["pc_type"] == "bddc":
modes = rigid_body_modes(V)
A.setNearNullSpace(modes)
modes.destroy()
return solve_cg(V, a, L, [bc], dofs, kind, pc_options, "elasticity", prepare)
The third problem is a curl-curl, or definite Maxwell, operator on a cube,
discretised with second-order Nedelec elements of the first kind, and with the tangential component of \(u\) set to zero on the boundary. It is solved in 3D, where \(\nabla \times\) is a vector. The source is taken as the curl of a smooth field, so it has no component along the gradients.
The curl of a gradient vanishes, so the curl-curl term alone has the whole range of the gradient in its kernel. The mass term is what makes the operator positive definite, and \(\varepsilon\) weights it: the smaller it is, the closer the operator is to the singular one and the harder the problem. It is set well below one here, so the curl-curl term dominates.
That near-kernel is what makes Maxwell problems hard, and PCBDDC has
dedicated support for it. It takes the discrete gradient, whose
range is the kernel of the curl, and puts that kernel in its coarse
space. :func:~dolfinx.fem.petsc.discrete_gradient assembles the
matrix, and the solver below passes it on.
def solve_curl_curl(
msh: mesh.Mesh, kind: str | None, pc_options: PCOptions
) -> tuple[fem.Function, int, int]:
"""Solve the curl-curl problem on ``msh`` using CG.
Args:
msh: Mesh, which must have been built without ghost cells.
kind: PETSc matrix kind, ``"is"`` for ``MATIS`` or ``None``
for the default assembled matrix.
pc_options: PETSc options that select and configure the
preconditioner, without a prefix.
Returns:
The solution, the number of Krylov iterations, and the number
of constrained degrees of freedom shared with another process.
"""
V = fem.functionspace(msh, ("N1curl", 2))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
g = ufl.as_vector((ufl.sin(ufl.pi * x[1]), ufl.sin(ufl.pi * x[2]), ufl.sin(ufl.pi * x[0])))
f = ufl.curl(g)
eps = 1.0e-2
a = (ufl.inner(ufl.curl(u), ufl.curl(v)) + eps * ufl.inner(u, v)) * ufl.dx
L = ufl.inner(f, v) * ufl.dx
# Zero tangential component on the boundary
tdim = msh.topology.dim
msh.topology.create_connectivity(tdim - 1, tdim)
dofs = fem.locate_dofs_topological(V, tdim - 1, mesh.exterior_facet_indices(msh.topology))
bcs = [fem.dirichletbc(fem.Function(V, dtype=dtype), dofs)]
def prepare(_A, ksp) -> None:
"""Give BDDC the discrete gradient.
Its range is the kernel of the curl, which BDDC puts into its
coarse space. PCBDDC analyses the subdomain edges with it, and
`conforming` asserts that each is a simple chain of
degrees-of-freedom between two subdomain corners. That is a
property of the partition rather than of the mesh, and holds
for the partitions this demo is run on. PETSc takes its own
reference, so this side keeps none.
"""
if pc_options["pc_type"] == "bddc":
# The gradient maps the H1 space of the same degree into V
W = fem.functionspace(msh, ("Lagrange", 2))
G = discrete_gradient(W, V)
G.assemble()
# `conforming=False` is the fallback where the partition
# does not give simple edges: PCBDDC otherwise stops with
# "Unexpected SIZE OF EDGE > EXTCOL SECOND PASS"
ksp.getPC().setBDDCDiscreteGradient(G, order=2, conforming=True)
G.destroy()
return solve_cg(V, a, L, bcs, dofs, kind, pc_options, "curl", prepare)
The number of iterations of both methods is close to independent of the mesh size, so refining the mesh does not slow convergence the way it would for a one-level method. Solving on a sequence of meshes shows this. Each mesh is built once and handed to every solver on it. The Poisson and elasticity problems are also solved with PCHPDDM when PETSc provides it and there is more than one process; with one there is no decomposition.
def report(label, n, uh, its, shared, pc_options, metric, value) -> None:
"""Print one solver's result, on rank 0 only."""
V = uh.function_space
num_dofs = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
shared = V.mesh.comm.allreduce(shared, MPI.SUM)
if V.mesh.comm.rank == 0:
pc_name = str(pc_options["pc_type"]).upper()
print(
f"{label}, n = {n:>3d}: {num_dofs:>7d} dofs, {its:>3d} CG iterations "
f"({pc_name}), {metric} = {value:.3e}, shared Dirichlet dofs = {shared}"
)
def has_hpddm() -> bool:
"""Whether a PCHPDDM preconditioner can be created.
PETSc reporting HPDDM among its packages is not enough: PCHPDDM
loads SLEPc from ``$SLEPC_DIR/lib`` when it is first created, and
that fails for a SLEPc built in place, whose library sits under
``$PETSC_ARCH``. Creating one is the only reliable test.
"""
if not PETSc.Sys.hasExternalPackage("hpddm"):
return False
pc = PETSc.PC().create(PETSc.COMM_SELF) # type: ignore[arg-type]
try:
# PETSc reports the failure to load before raising, so an
# error here is printed once and then explained below
pc.setType("hpddm")
return True
except PETSc.Error:
return False
finally:
pc.destroy()
comm = MPI.COMM_WORLD
preconditioners: list[tuple[str | None, PCOptions]] = [("is", bddc_options)]
# PCHPDDM has nothing to decompose on one process: it builds only its
# coarse level, leaving the Schwarz options unused, so it is added only
# in parallel.
if has_hpddm() and comm.size > 1:
preconditioners.append((None, hpddm_options))
elif comm.rank == 0 and comm.size > 1:
print(
"PCHPDDM could not be created, so only BDDC is shown. PCHPDDM loads "
"SLEPc from $SLEPC_DIR/lib when first created; for a SLEPc built in "
"place that is $SLEPC_DIR/$PETSC_ARCH/lib."
)
for n in (32, 64):
# BDDC requires one non-overlapping subdomain per process, so the
# mesh is built without ghost cells
msh = mesh.create_unit_square(
comm, n, n, mesh.CellType.triangle, ghost_mode=mesh.GhostMode.none, dtype=xdtype
)
x = ufl.SpatialCoordinate(msh)
u_exact = ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1])
for kind, pc_options in preconditioners:
uh, its, shared = solve_poisson(msh, kind, pc_options)
report("Poisson ", n, uh, its, shared, pc_options, "L2 error", norm_L2(uh - u_exact))
for kind, pc_options in preconditioners:
uh, its, shared = solve_elasticity(msh, kind, pc_options)
report("Elasticity", n, uh, its, shared, pc_options, "|u|_L2", norm_L2(uh))
# The curl-curl problem is posed in 3D, so it gets meshes of its own,
# and it is solved with BDDC alone.
for n in (8, 12):
msh = mesh.create_unit_cube(
comm, n, n, n, mesh.CellType.tetrahedron, ghost_mode=mesh.GhostMode.none, dtype=xdtype
)
uh, its, shared = solve_curl_curl(msh, "is", bddc_options)
report("Curl-curl ", n, uh, its, shared, bddc_options, "|u|_L2", norm_L2(uh))
With more than one process, shared Dirichlet dofs is non-zero: those
are the constrained degrees of freedom where a subdomain interface
meets the Dirichlet boundary, and they are the reason the boundary
condition value has to be distributed across the processes that share
them rather than written by the owner alone. How many there are
follows from how much of the boundary is constrained and from the
element: the Poisson problem is constrained on the whole boundary and
the elasticity problem on one edge, while the curl-curl problem
constrains the tangential component over the whole boundary of a 3D
mesh, where an interface meets it along a line of edges rather than
at a point.