In this tutorial we solve the problem
$$\begin{cases} -\Delta u = f, & \text{in } \Omega,\\ u = g, & \text{on } \Gamma = \partial\Omega, \end{cases}$$
where $\Omega$ is a ball in 2D.
We compare the following two cases: the corresponding weak formulation is
$$ \text{find } u \in V_g \text{ s.t. } \int_\Omega \nabla u \cdot \nabla v = \int_\Omega f v, \quad \forall v \in V_0\\ $$ where $$ V_g = \{v \in H^1(\Omega): v|_\Gamma = g\},\\ V_0 = \{v \in H^1(\Omega): v|_\Gamma = 0\}.\\ $$
This example is a prototypical case of problems containing subdomain/boundary restricted variables (the Lagrange multiplier, in this case).
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import ufl
import viskex
import multiphenicsx.fem
import multiphenicsx.fem.petsc
r = 3
mesh_size = 1. / 4.
gmsh.initialize()
gmsh.model.add("mesh")
p0 = gmsh.model.geo.addPoint(0.0, 0.0, 0.0, mesh_size)
p1 = gmsh.model.geo.addPoint(0.0, +r, 0.0, mesh_size)
p2 = gmsh.model.geo.addPoint(0.0, -r, 0.0, mesh_size)
c0 = gmsh.model.geo.addCircleArc(p1, p0, p2)
c1 = gmsh.model.geo.addCircleArc(p2, p0, p1)
boundary = gmsh.model.geo.addCurveLoop([c0, c1])
domain = gmsh.model.geo.addPlaneSurface([boundary])
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [c0, c1], 1)
gmsh.model.addPhysicalGroup(2, [boundary], 0)
gmsh.model.mesh.generate(2)
Info : Meshing 1D... Info : [ 0%] Meshing curve 1 (Circle) Info : [ 50%] Meshing curve 2 (Circle) Info : Done meshing 1D (Wall 0.000317883s, CPU 0.000548s) Info : Meshing 2D... Info : Meshing surface 1 (Plane, Frontal-Delaunay) Info : Done meshing 2D (Wall 0.0319924s, CPU 0.012821s) Info : 589 nodes 1177 elements
mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(
gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
gmsh.finalize()
# Create connectivities required by the rest of the code
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
facets_Gamma = boundaries.indices[boundaries.values == 1]
viskex.dolfinx.plot_mesh(mesh)
viskex.dolfinx.plot_mesh_tags(mesh, boundaries, "boundaries")
# Define a function space
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))
M = V.clone()
# Define restrictions
dofs_V = np.arange(0, V.dofmap.index_map.size_local + V.dofmap.index_map.num_ghosts)
dofs_M_Gamma = dolfinx.fem.locate_dofs_topological(M, boundaries.dim, facets_Gamma)
restriction_V = multiphenicsx.fem.DofMapRestriction(V.dofmap, dofs_V)
restriction_M_Gamma = multiphenicsx.fem.DofMapRestriction(M.dofmap, dofs_M_Gamma)
restriction = [restriction_V, restriction_M_Gamma]
# Define trial and test functions
(u, l) = (ufl.TrialFunction(V), ufl.TrialFunction(M))
(v, m) = (ufl.TestFunction(V), ufl.TestFunction(M))
# Define problem block forms
g = dolfinx.fem.Function(V)
g.interpolate(lambda x: np.sin(3 * x[0] + 1) * np.sin(3 * x[1] + 1))
a = [[ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx, ufl.inner(l, v) * ufl.ds],
[ufl.inner(u, m) * ufl.ds, None]]
f = [ufl.inner(1, v) * ufl.dx, ufl.inner(g, m) * ufl.ds]
a_cpp = dolfinx.fem.form(a)
f_cpp = dolfinx.fem.form(f)
# Assemble the block linear system
A = multiphenicsx.fem.petsc.assemble_matrix_block(a_cpp, bcs=[], restriction=(restriction, restriction))
A.assemble()
F = multiphenicsx.fem.petsc.assemble_vector_block(f_cpp, a_cpp, bcs=[], restriction=restriction)
# Solve
ul = multiphenicsx.fem.petsc.create_vector_block(f_cpp, restriction=restriction)
ksp = petsc4py.PETSc.KSP()
ksp.create(mesh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
ksp.setFromOptions()
ksp.solve(F, ul)
ul.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
ksp.destroy()
<petsc4py.PETSc.KSP at 0x7fd69c62e5c0>
# Split the block solution in components
(u, l) = (dolfinx.fem.Function(V), dolfinx.fem.Function(M))
with multiphenicsx.fem.petsc.BlockVecSubVectorWrapper(ul, [V.dofmap, M.dofmap], restriction) as ul_wrapper:
for ul_wrapper_local, component in zip(ul_wrapper, (u, l)):
with component.x.petsc_vec.localForm() as component_local:
component_local[:] = ul_wrapper_local
ul.destroy()
<petsc4py.PETSc.Vec at 0x7fd69c62fe70>
viskex.dolfinx.plot_scalar_field(u, "u")
viskex.dolfinx.plot_scalar_field(l, "l")
# Define Dirichlet BC object on Gamma
dofs_V_Gamma = dolfinx.fem.locate_dofs_topological(V, boundaries.dim, facets_Gamma)
bc_ex = dolfinx.fem.dirichletbc(g, dofs_V_Gamma)
# Solve
u_ex = dolfinx.fem.Function(V)
problem_ex = dolfinx.fem.petsc.LinearProblem(
a[0][0], f[0], bcs=[bc_ex], u=u_ex,
petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
problem_ex.solve()
u_ex.x.petsc_vec.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
viskex.dolfinx.plot_scalar_field(u_ex, "u")
u_ex_norm = np.sqrt(mesh.comm.allreduce(
dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(ufl.grad(u_ex), ufl.grad(u_ex)) * ufl.dx)),
op=mpi4py.MPI.SUM))
err_norm = np.sqrt(mesh.comm.allreduce(
dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(ufl.grad(u_ex - u), ufl.grad(u_ex - u)) * ufl.dx)),
op=mpi4py.MPI.SUM))
print("Relative error is equal to", err_norm / u_ex_norm)
assert np.isclose(err_norm / u_ex_norm, 0., atol=1.e-10)
Relative error is equal to 5.75871463647993e-15