Real function spaces#
Author: Jørgen S. Dokken
License: MIT
In this example we will show how to use the “real” function space to solve a singular Poisson problem.
Mathematical formulation#
The problem at hand is: Find \(u \in H^1(\Omega)\) such that
\[\begin{split}
\begin{align}
-\Delta u &= f \quad \text{in } \Omega, \\
\frac{\partial u}{\partial n} &= g \quad \text{on } \partial \Omega, \\
\int_\Omega u &= h.
\end{align}
\end{split}\]
Lagrange multiplier#
We start by considering the equivalent optimization problem: Find \(u \in H^1(\Omega)\) such that
\[
\begin{align}
\mathop{\mathrm{arg\,min}}_{u \in H^1(\Omega)}~J(u) =
\mathop{\mathrm{arg\,min}}_{u \in H^1(\Omega)}~\frac{1}{2}\int_\Omega \vert \nabla u \cdot \nabla u \vert \mathrm{d}x
- \int_\Omega f u \mathrm{d}x - \int_{\partial \Omega} g u \mathrm{d}s,
\end{align}
\]
such that
\[
\begin{align}
\int_\Omega u = h.
\end{align}
\]
We introduce a Lagrange multiplier \(\lambda\) to enforce the constraint:
\[
\begin{align}
\mathop{\mathrm{arg\,min}}_{u \in H^1(\Omega), \lambda\in \mathbb{R}}~\mathcal{L}(u, \lambda) =
\mathop{\mathrm{arg\,min}}_{u \in H^1(\Omega), \lambda\in \mathbb{R}} J(u) + \lambda (\int_\Omega u \mathrm{d}x-h).
\end{align}
\]
We then compute the optimality conditions for the problem above
\[\begin{split}
\begin{align}
\frac{\partial \mathcal{L}}{\partial u}[\delta u] &= \int_\Omega \nabla u \cdot \nabla \delta u \mathrm{d}x + \lambda\int \delta u \mathrm{d}x - \int_\Omega f \delta u ~\mathrm{d}x - \int_{\partial \Omega} g \delta u~\mathrm{d}s = 0, \\
\frac{\partial \mathcal{L}}{\partial \lambda}[\delta \lambda] &=\delta \lambda (\int_\Omega u \mathrm{d}x -h)= 0.
\end{align}
\end{split}\]
We write the weak formulation:
\[\begin{split}
\begin{align}
\int_\Omega \nabla u \cdot \nabla \delta u~\mathrm{d}x + \int_\Omega \lambda \delta u~\mathrm{d}x &= \int_\Omega f \delta u~\mathrm{d}x + \int_{\partial \Omega} g v \mathrm{d}s\\
\int_\Omega u \delta \lambda \mathrm{d}x &= h \delta \lambda .
\end{align}
\end{split}\]
where we have moved \(\delta\lambda\) into the integral as it is a spatial constant.
Implementation#
We start by import the necessary modules
Clickable functions/classes
Note that for the modules imported in this example, you can click on the function/class name to be redirected to the corresponding documentation page.
from packaging.version import Version
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx.fem.petsc
import numpy as np
from scifem import create_real_functionspace, assemble_scalar
from scifem.petsc import apply_lifting_and_set_bc
import ufl
import pyvista
We start by creating the domain using dolfinx and derive the source terms
\(f\), \(g\) and \(h\) from our manufactured solution using ufl.
For this example we will use the following exact solution
\[
\begin{align}
u_{exact}(x, y) = 0.3y^2 + \sin(2\pi x).
\end{align}
\]
M = 20
mesh = dolfinx.mesh.create_unit_square(
MPI.COMM_WORLD, M, M, dolfinx.mesh.CellType.triangle, dtype=np.float64
)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
def u_exact(x):
return 0.3 * x[1] ** 2 + ufl.sin(2 * ufl.pi * x[0])
x = ufl.SpatialCoordinate(mesh)
n = ufl.FacetNormal(mesh)
g = ufl.dot(ufl.grad(u_exact(x)), n)
f = -ufl.div(ufl.grad(u_exact(x)))
h = assemble_scalar(u_exact(x) * ufl.dx)
Creating the real function space and mixed space#
We then create the Lagrange multiplier space using scifem.create_real_functionspace().
This creates a dolfinx.fem.FunctionSpace with a single degree of freedom (constant) over the whole domain.
/tmp/ipykernel_1724/957554804.py:1: DeprecationWarning: `create_real_functionspace()` is deprecated and will be removed in a future version.Please use `basix.ufl.create_real_element()` instead.
R = create_real_functionspace(mesh)
Next, we can create a mixed-function space for our problem
Note on DOLFINx versions
The API for creating blocked problems in DOLFINx has vastly improved over the last few versions.
It is recommended to use DOLFINx v0.9.0 or later, where one can use ufl.MixedFunctionSpace
to create a mixed-space of M number of spaces. One can then use ufl.TrialFunctions()
and ufl.TestFunctions() as one would do with a basix.ufl.mixed_element().
if Version(dolfinx.__version__) == Version("0.8.0"):
u = ufl.TrialFunction(V)
lmbda = ufl.TrialFunction(R)
du = ufl.TestFunction(V)
dl = ufl.TestFunction(R)
elif Version(dolfinx.__version__) >= Version("0.9.0.0"):
W = ufl.MixedFunctionSpace(V, R)
u, lmbda = ufl.TrialFunctions(W)
du, dl = ufl.TestFunctions(W)
else:
raise RuntimeError("Unsupported version of dolfinx")
Defining and assembling the variational problem#
We can now define the variational problem
zero = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.0))
a00 = ufl.inner(ufl.grad(u), ufl.grad(du)) * ufl.dx
a01 = ufl.inner(lmbda, du) * ufl.dx
a10 = ufl.inner(u, dl) * ufl.dx
L0 = ufl.inner(f, du) * ufl.dx + ufl.inner(g, du) * ufl.ds
L1 = ufl.inner(zero, dl) * ufl.dx
a = [[a00, a01], [a10, None]]
L = [L0, L1]
a_compiled = dolfinx.fem.form(a)
L_compiled = dolfinx.fem.form(L)
Note that we have defined the variational form in a block form, and
that we have not included \(h\) in the variational form. We will enforce this
once we have assembled the right hand side vector.
We can now assemble the matrix and vector usig dolfinx.fem.petsc.assemble_matrix()
and dolfinx.fem.petsc.assemble_vector().
if Version(dolfinx.__version__) < Version("0.10.0"):
A = dolfinx.fem.petsc.assemble_matrix_block(a_compiled)
else:
A = dolfinx.fem.petsc.assemble_matrix(a_compiled)
A.assemble()
In DOLFINx>=v0.10.0, the assemble_vector function for blocked spaces has been rewritten to reflect how
it works for standard assembly and nest assembly. This means that lifting is applied manually.
In this case, with no Dirichlet BC, we could skip those steps.
However, for clarity we include them here.
bcs = []
if Version(dolfinx.__version__) < Version("0.10.0"):
b = dolfinx.fem.petsc.assemble_vector_block(L_compiled, a_compiled, bcs=bcs)
else:
b = dolfinx.fem.petsc.assemble_vector(L_compiled, kind="mpi")
apply_lifting_and_set_bc(b, a_compiled, bcs=bcs)
Next, we modify the second part of the block to contain h
We start by enforcing the multiplier constraint \(h\) by modifying the right hand side vector.
On the main branch, this is greatly simplified
uh = dolfinx.fem.Function(V, name="u")
if Version(dolfinx.__version__) >= Version("0.10.0"):
# We start by inserting the value in the real space
rh = dolfinx.fem.Function(R)
rh.x.array[0] = h
# Next we need to add this value to the existing right hand side vector.
# Therefore we create assign 0s to the primal space
b_real_space = b.duplicate()
uh.x.array[:] = 0
# Transfer the data to `b_real_space`
dolfinx.fem.petsc.assign([uh, rh], b_real_space)
# And accumulate the values in the right hand side vector
b.axpy(1, b_real_space)
# We destroy the temporary work vector after usage
b_real_space.destroy()
else:
from dolfinx.cpp.la.petsc import scatter_local_vectors, get_local_vectors
if Version(dolfinx.__version__) < Version("0.9.0"):
maps = [(V.dofmap.index_map, V.dofmap.index_map_bs), (R.dofmap.index_map, R.dofmap.index_map_bs)]
else:
maps = [(Wi.dofmap.index_map, Wi.dofmap.index_map_bs) for Wi in W.ufl_sub_spaces()]
b_local = get_local_vectors(b, maps)
b_local[1][:] = h
scatter_local_vectors(
b,
b_local,
maps,
)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
Solving the linear system#
We can now solve the linear system using petsc4py.
ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
pc = ksp.getPC()
pc.setType("lu")
pc.setFactorSolverType("mumps")
if Version(dolfinx.__version__) >= Version("0.10.0"):
xh = b.duplicate()
else:
xh = dolfinx.fem.petsc.create_vector_block(L_compiled)
ksp.solve(b, xh)
xh.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
Finally, we extract the solution u from the blocked system and compute the error
uh = dolfinx.fem.Function(V, name="u")
if Version(dolfinx.__version__) >= Version("0.10.0"):
dolfinx.fem.petsc.assign(xh, [uh, rh])
else:
x_local = get_local_vectors(xh, maps)
uh.x.array[: len(x_local[0])] = x_local[0]
uh.x.scatter_forward()
We destroy all PETSc objects
b.destroy()
xh.destroy()
A.destroy()
ksp.destroy()
<petsc4py.PETSc.KSP at 0x7f9c4aee2200>
Post-processing#
Finally, we compare our approximate solution with the exact solution
by computing the \(L^2(\Omega)\) error.
We use the convenience function scifem.assemble_scalar() to compute the error.
diff = uh - u_exact(x)
error = ufl.inner(diff, diff) * ufl.dx
print(f"L2 error: {np.sqrt(assemble_scalar(error)):.2e}")
L2 error: 6.73e-03
We can now plot the solution
vtk_mesh = dolfinx.plot.vtk_mesh(V)
grid = pyvista.UnstructuredGrid(*vtk_mesh)
grid.point_data["u"] = uh.x.array.real
warped = grid.warp_by_scalar("u", factor=1)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, style="wireframe")
plotter.add_mesh(warped)
if not pyvista.OFF_SCREEN:
plotter.show()
2026-08-05 06:37:36.332 ( 1.938s) [ 7F9CB5568140]vtkXOpenGLRenderWindow.:1460 WARN| bad X server connection. DISPLAY=:99.0
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_urldispatcher.py:546, in StaticResource.__init__(self, prefix, directory, name, expect_handler, chunk_size, show_index, follow_symlinks, append_version)
545 try:
--> 546 directory = Path(directory).expanduser().resolve(strict=True)
547 except FileNotFoundError as error:
File /usr/lib/python3.12/pathlib.py:1242, in Path.resolve(self, strict)
1241 try:
-> 1242 s = self._flavour.realpath(self, strict=strict)
1243 except OSError as e:
File <frozen posixpath>:432, in realpath(filename, strict)
File <frozen posixpath>:477, in _joinrealpath(path, rest, strict, seen)
FileNotFoundError: [Errno 2] No such file or directory: '/dolfinx-env/lib/python3.12/site-packages/trame_client/module/vue3-www'
The above exception was the direct cause of the following exception:
ValueError Traceback (most recent call last)
Cell In[14], line 11
7 plotter = pyvista.Plotter()
8 plotter.add_mesh(grid, style="wireframe")
9 plotter.add_mesh(warped)
10 if not pyvista.OFF_SCREEN:
---> 11 plotter.show()
File /dolfinx-env/lib/python3.12/site-packages/pyvista/_deprecate_positional_args.py:243, in _deprecate_positional_args.<locals>._inner_deprecate_positional_args.<locals>.inner_f(*args, **kwargs)
239 warn_external(msg, PyVistaDeprecationWarning)
241 warn_positional_args()
--> 243 return f(*args, **kwargs)
File /dolfinx-env/lib/python3.12/site-packages/pyvista/plotting/plotter.py:8231, in Plotter.show(self, title, window_size, interactive, auto_close, interactive_update, full_screen, screenshot, return_img, cpos, jupyter_backend, return_viewer, return_cpos, before_close_callback, store_image_depth, **kwargs)
8228 jupyter_backend = self._theme.jupyter_backend
8230 if jupyter_backend is None or jupyter_backend.lower() != 'none':
-> 8231 jupyter_disp = handle_plotter(self, backend=jupyter_backend, **jupyter_kwargs)
8233 self.render()
8235 # initial double render needed for certain passes when offscreen
File /dolfinx-env/lib/python3.12/site-packages/pyvista/jupyter/notebook.py:74, in handle_plotter(plotter, backend, screenshot, **kwargs)
71 if backend in ['server', 'client', 'trame', 'html']:
72 from pyvista.trame.jupyter import show_trame # noqa: PLC0415
---> 74 return show_trame(plotter, mode=backend, **kwargs)
76 except ImportError as e:
77 # Trame was explicitly requested but not available
78 _ensure_entry_points()
File /dolfinx-env/lib/python3.12/site-packages/pyvista/trame/jupyter.py:394, in show_trame(plotter, mode, name, server_proxy_enabled, server_proxy_prefix, jupyter_extension_enabled, collapse_menu, add_menu, add_menu_items, default_server_rendering, handler, **kwargs)
391 kwargs.setdefault('height', dh)
393 if mode == 'html':
--> 394 return EmbeddableWidget(plotter, **kwargs)
396 if jupyter_extension_enabled is None:
397 jupyter_extension_enabled = pv.global_theme.trame.jupyter_extension_enabled
File /dolfinx-env/lib/python3.12/site-packages/pyvista/trame/jupyter.py:137, in EmbeddableWidget.__init__(self, plotter, width, height, **kwargs)
135 msg = 'Please install `ipywidgets`.'
136 raise ImportError(msg)
--> 137 scene = plotter.export_html(filename=None)
138 src = scene.getvalue().replace('"', '"')
139 # eventually we could maybe expose this, but for now make sure we're at least
140 # consistent with matplotlib's color (light gray)
File /dolfinx-env/lib/python3.12/site-packages/pyvista/plotting/plotter.py:784, in BasePlotter.export_html(self, filename)
781 msg = 'Please install trame dependencies: pip install "pyvista[jupyter]"'
782 raise ImportError(msg)
--> 784 data = self.export_vtksz(filename=None)
785 buffer = io.StringIO()
786 write_html(data, buffer)
File /dolfinx-env/lib/python3.12/site-packages/pyvista/plotting/plotter.py:837, in BasePlotter.export_vtksz(self, filename, format)
835 server = get_server(pv.global_theme.trame.jupyter_server_name)
836 if not server.running:
--> 837 elegantly_launch(pv.global_theme.trame.jupyter_server_name)
839 view = PyVistaLocalView(self, trame_server=server)
841 content = view.export(format=format)
File /dolfinx-env/lib/python3.12/site-packages/pyvista/trame/jupyter.py:487, in elegantly_launch(*args, **kwargs)
484 # Basically monkey patches asyncio to support this
485 nest_asyncio2.apply()
--> 487 return asyncio.run(launch_it())
File /dolfinx-env/lib/python3.12/site-packages/nest_asyncio2.py:115, in run(main, debug, loop_factory)
113 task = asyncio.ensure_future(main, loop=loop)
114 try:
--> 115 return loop.run_until_complete(task)
116 finally:
117 if not task.done():
File /dolfinx-env/lib/python3.12/site-packages/nest_asyncio2.py:230, in _patch_loop.<locals>.run_until_complete(self, future)
227 if self._ready:
228 self._run_once()
--> 230 return f.result()
File /usr/lib/python3.12/asyncio/futures.py:203, in Future.result(self)
201 self.__log_traceback = False
202 if self._exception is not None:
--> 203 raise self._exception.with_traceback(self._exception_tb)
204 return self._result
File /usr/lib/python3.12/asyncio/tasks.py:314, in Task.__step_run_and_handle_result(***failed resolving arguments***)
310 try:
311 if exc is None:
312 # We use the `send` method directly, because coroutines
313 # don't have `__iter__` and `__next__` methods.
--> 314 result = coro.send(None)
315 else:
316 result = coro.throw(exc)
File /dolfinx-env/lib/python3.12/site-packages/pyvista/trame/jupyter.py:482, in elegantly_launch.<locals>.launch_it()
481 async def launch_it():
--> 482 await launch_server(*args, **kwargs).ready
File /dolfinx-env/lib/python3.12/site-packages/pyvista/trame/jupyter.py:215, in launch_server(server, port, host, wslink_backend, **kwargs)
213 if server._running_stage == 0:
214 server.controller.on_server_ready.add(on_ready)
--> 215 server.start(
216 exec_mode='task',
217 host=host,
218 port=port,
219 open_browser=False,
220 show_connection_info=False,
221 disable_logging=True,
222 timeout=0,
223 backend=wslink_backend,
224 )
225 # else, server is already running or launching
226 return server
File /dolfinx-env/lib/python3.12/site-packages/trame_server/core.py:725, in Server.start(self, port, thread, open_browser, show_connection_info, disable_logging, backend, follow_symlinks, exec_mode, timeout, host, **kwargs)
722 CoreServer.configure(options)
724 self._running_stage = 1
--> 725 task = CoreServer.server_start(
726 options,
727 **{ # Do a proper merging/override
728 **kwargs,
729 "disableLogging": disable_logging,
730 "backend": backend,
731 "exec_mode": exec_mode,
732 },
733 )
735 # Manage exit life cycle unless coroutine
736 if exec_mode == "main":
File /dolfinx-env/lib/python3.12/site-packages/trame_server/protocol.py:52, in CoreServer.server_start(options, disableLogging, backend, exec_mode, **kwargs)
43 @staticmethod
44 def server_start(
45 options,
(...) 50 ):
51 # NOTE: **kwargs to wslink's start_webserver are currently unused
---> 52 return server.start_webserver(
53 options=options,
54 protocol=CoreServer,
55 disableLogging=disableLogging,
56 backend=backend,
57 exec_mode=exec_mode,
58 **kwargs,
59 )
File /dolfinx-env/lib/python3.12/site-packages/wslink/server.py:257, in start_webserver(options, protocol, disableLogging, backend, exec_mode, **_)
254 server_config["handle_signals"] = not options.nosignalhandlers
256 # Create the webserver and start it
--> 257 ws_server = create_webserver(server_config, backend=backend)
259 # Once we have python 3.7 minimum, we can start the server with asyncio.run()
260 # asyncio.run(ws_server.start())
261
262 # Until then, we can start the server this way
263 try:
File /dolfinx-env/lib/python3.12/site-packages/wslink/server.py:168, in create_webserver(server_config, backend)
167 def create_webserver(server_config, backend="aiohttp"):
--> 168 return backends.create_webserver(server_config, backend=backend)
File /dolfinx-env/lib/python3.12/site-packages/wslink/backends/__init__.py:5, in create_webserver(server_config, backend)
2 if backend == "aiohttp":
3 from .aiohttp import create_webserver # noqa: PLC0415
----> 5 return create_webserver(server_config)
7 if backend == "generic":
8 from .generic import create_webserver # noqa: PLC0415
File /dolfinx-env/lib/python3.12/site-packages/wslink/backends/aiohttp/__init__.py:218, in create_webserver(server_config)
215 return ReverseWebAppServer(server_config)
217 # Normal web server
--> 218 return WebAppServer(server_config)
File /dolfinx-env/lib/python3.12/site-packages/wslink/backends/aiohttp/__init__.py:118, in WebAppServer.__init__(self, server_config)
116 # Resolve / => index.html
117 self.app.router.add_route("GET", "/", _root_handler)
--> 118 self.app.add_routes(routes)
120 self.app[STATE_KEY] = {}
File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_app.py:379, in Application.add_routes(self, routes)
378 def add_routes(self, routes: Iterable[AbstractRouteDef]) -> list[AbstractRoute]:
--> 379 return self.router.add_routes(routes)
File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_urldispatcher.py:1259, in UrlDispatcher.add_routes(self, routes)
1257 registered_routes = []
1258 for route_def in routes:
-> 1259 registered_routes.extend(route_def.register(self))
1260 return registered_routes
File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_routedef.py:87, in StaticDef.register(self, router)
86 def register(self, router: UrlDispatcher) -> list[AbstractRoute]:
---> 87 resource = router.add_static(self.prefix, self.path, **self.kwargs)
88 routes = resource.get_info().get("routes", {})
89 return list(routes.values())
File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_urldispatcher.py:1183, in UrlDispatcher.add_static(self, prefix, path, name, expect_handler, chunk_size, show_index, follow_symlinks, append_version)
1181 if prefix.endswith("/"):
1182 prefix = prefix[:-1]
-> 1183 resource = StaticResource(
1184 prefix,
1185 path,
1186 name=name,
1187 expect_handler=expect_handler,
1188 chunk_size=chunk_size,
1189 show_index=show_index,
1190 follow_symlinks=follow_symlinks,
1191 append_version=append_version,
1192 )
1193 self.register_resource(resource)
1194 return resource
File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_urldispatcher.py:548, in StaticResource.__init__(self, prefix, directory, name, expect_handler, chunk_size, show_index, follow_symlinks, append_version)
546 directory = Path(directory).expanduser().resolve(strict=True)
547 except FileNotFoundError as error:
--> 548 raise ValueError(f"'{directory}' does not exist") from error
549 if not directory.is_dir():
550 raise ValueError(f"'{directory}' is not a directory")
ValueError: '/dolfinx-env/lib/python3.12/site-packages/trame_client/module/vue3-www' does not exist