From 13469ef217e31fa7c06d67750a186ce404df841b Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Sun, 1 Nov 2026 23:18:21 +0100 Subject: [PATCH] Fix mixed poisson for nightly pullback + a typo in fundamentals --- chapter1/fundamentals_code.ipynb | 2 +- chapter1/fundamentals_code.py | 2 +- chapter4/mixed_poisson.ipynb | 33 ++++++++++++++++---------------- chapter4/mixed_poisson.py | 30 +++++++++++++++-------------- 4 files changed, 34 insertions(+), 33 deletions(-) diff --git a/chapter1/fundamentals_code.ipynb b/chapter1/fundamentals_code.ipynb index 8f71f9a6..53b7099c 100644 --- a/chapter1/fundamentals_code.ipynb +++ b/chapter1/fundamentals_code.ipynb @@ -441,7 +441,7 @@ "Secondly, we compute the maximum error at any degree of freedom.\n", "As the finite element function $u$ can be expressed as a linear combination of basis functions $\\phi_j$,\n", "spanning the space $V$: $ u = \\sum_{j=1}^N U_j\\phi_j.$\n", - "By writing {py:meth}`problem.solve()`\n", + "By writing {py:meth}`problem.solve()`\n", "we compute all the coefficients $U_1,\\dots, U_N$.\n", "These values are known as the _degrees of freedom_ (dofs).\n", "We can access the degrees of freedom by accessing the underlying vector in `uh`.\n", diff --git a/chapter1/fundamentals_code.py b/chapter1/fundamentals_code.py index d22d1253..f5dca4f8 100644 --- a/chapter1/fundamentals_code.py +++ b/chapter1/fundamentals_code.py @@ -291,7 +291,7 @@ # Secondly, we compute the maximum error at any degree of freedom. # As the finite element function $u$ can be expressed as a linear combination of basis functions $\phi_j$, # spanning the space $V$: $ u = \sum_{j=1}^N U_j\phi_j.$ -# By writing {py:meth}`problem.solve()` +# By writing {py:meth}`problem.solve()` # we compute all the coefficients $U_1,\dots, U_N$. # These values are known as the _degrees of freedom_ (dofs). # We can access the degrees of freedom by accessing the underlying vector in `uh`. diff --git a/chapter4/mixed_poisson.ipynb b/chapter4/mixed_poisson.ipynb index be934fec..cbcf1e8b 100644 --- a/chapter4/mixed_poisson.ipynb +++ b/chapter4/mixed_poisson.ipynb @@ -6,7 +6,7 @@ "metadata": {}, "source": [ "# Mixed Poisson with a Schur complement pre-conditioner\n", - "This example demonstrates how to use PETSc fieldsplits with custom preconditioners in DOLFINx.\n", + "This example demonstrates how to use {py:mod}`PETSc` fieldsplits with custom preconditioners in DOLFINx.\n", "This example is heavily inspired by the [FEniCSx PCTools example](https://rafinex-external-rifle.gitlab.io/fenicsx-pctools/demo/demo_mixed-poisson.html)\n", "which was presented in {cite}`mp-rehor2025pctools`." ] @@ -32,9 +32,6 @@ "As we will see below, the roles of these conditions are swapped compared to the primal formulation of the Poisson equation.\n", "The Dirichlet condition becomes a natural boundary condition, which enters the variational formulation through a boundary integral,\n", "while the Neumann condition becomes an essential boundary condition, which is enforced strongly on the function space for $\\sigma$.\n", - "In particular, we do not impose the Dirichlet condition strongly with a `dolfinx.fem.dirichletbc`, as in the previous Poisson examples.\n", - "This would not be possible either, as we will approximate $u$ in the zeroth order discontinuous Lagrange space,\n", - "whose degrees of freedom are associated with the cells, and not with the boundary facets.\n", "\n", "As in previous examples, we pick a manufactured solution to ensure that we can verify\n", "the correctness of our implementation.\n", @@ -136,7 +133,7 @@ "id": "9", "metadata": {}, "source": [ - "We define a `ufl.MixedFunctionSpace` to automatically handle the block structure of the problem" + "We define a {py:class}`ufl.MixedFunctionSpace` to automatically handle the block structure of the problem" ] }, { @@ -176,7 +173,7 @@ "- We enforce the continuity of $u$ weakly by removing the jump term.\n", "- On $\\Gamma_D$, we insert the Dirichlet condition $u=u_D$.\n", " The Dirichlet condition is thus enforced weakly, as a natural boundary condition.\n", - "- On $\\Gamma_N$, we enforce the Neumann condition $\\sigma\\cdot \\mathbf{n}=g$ strongly, by using a `dolfinx.fem.dirichletbc` for $\\sigma$.\n", + "- On $\\Gamma_N$, we enforce the Neumann condition $\\sigma\\cdot \\mathbf{n}=g$ strongly, by using a {py:func}`dolfinx.fem.dirichletbc` for $\\sigma$.\n", " The test function $\\tau$ then satisfies $\\tau\\cdot\\mathbf{n}=0$ on $\\Gamma_N$, and the integral over $\\Gamma_N$ vanishes.\n", "\n", "Thus we end up with:\n", @@ -252,7 +249,7 @@ "b_1\n", "\\end{pmatrix}\n", "\\end{align}\n", - "We can extract the block structure of the bilinear form using `ufl.extract_blocks`, which returns a nested list of bilinear forms.\n", + "We can extract the block structure of the bilinear form using {py:func}`ufl.extract_blocks`, which returns a nested list of bilinear forms.\n", "You can also build this nested list by hand if you want to, but it is usually more error-prone." ] }, @@ -279,14 +276,14 @@ "metadata": {}, "source": [ "Next, we create the essential boundary condition $\\sigma\\cdot\\mathbf{n}=g$ on $\\Gamma_N$,\n", - "which in DOLFINx is a `dolfinx.fem.dirichletbc` for $\\sigma$.\n", + "which in DOLFINx is a {py:func}`dolfinx.fem.dirichletbc` for $\\sigma$.\n", "The degrees of freedom of the first order BDM space are moments of the normal component $\\sigma\\cdot\\mathbf{n}$ on each facet.\n", "Thus, we can interpolate any vector field whose normal component is equal to $g$ on $\\Gamma_N$.\n", "We use $g\\mathbf{n}$, which only requires knowledge of $g$, and not of the full flux $\\sigma$ on the boundary.\n", "As we are using manufactured solutions for this problem, we could manually derive the explicit expression\n", "for $\\sigma$ on the boundary $\\Gamma_N$.\n", "However, in general this is not possible (especially for curved boundaries), and we have to use a more generic approach.\n", - "For this we will use the `dolfinx.fem.Expression` class to interpolate the expression into the function space $Q$.\n", + "For this we will use the {py:class}`dolfinx.fem.Expression` class to interpolate the expression into the function space $Q$.\n", "This is done by evaluating the expression at the physical interpolation points of the mesh.\n", "A convenience function for this is provided in the `interpolate_facet_expression` function below." ] @@ -333,13 +330,15 @@ " facet_cmap = basix.ufl.element(\n", " \"Lagrange\", facet_type, c_el.degree, shape=(domain.geometry.dim,), dtype=x_type\n", " )\n", + " # Wrap the C++ coordinate element in the Python interface, as the signature of\n", + " # the C++ `pull_back` differs between DOLFINx versions\n", " if np.issubdtype(x_type, np.float32):\n", - " facet_cel = dolfinx.cpp.fem.CoordinateElement_float32(\n", - " facet_cmap.basix_element._e\n", + " facet_cel = dolfinx.fem.CoordinateElement(\n", + " dolfinx.cpp.fem.CoordinateElement_float32(facet_cmap.basix_element._e)\n", " )\n", " elif np.issubdtype(x_type, np.float64):\n", - " facet_cel = dolfinx.cpp.fem.CoordinateElement_float64(\n", - " facet_cmap.basix_element._e\n", + " facet_cel = dolfinx.fem.CoordinateElement(\n", + " dolfinx.cpp.fem.CoordinateElement_float64(facet_cmap.basix_element._e)\n", " )\n", " else:\n", " raise TypeError(\n", @@ -353,7 +352,7 @@ " interpolation_points = Q_el.basix_element.x\n", " for i, points in enumerate(interpolation_points[fdim]):\n", " geom = ref_geom[ref_top[fdim][i]]\n", - " ref_points = facet_cel.pull_back(points, geom, tol=1e-6, maxit=25)\n", + " ref_points = facet_cel.pull_back(points, geom)\n", " # Assert that interpolation points are all equal on all facets\n", " if reference_facet_points is None:\n", " reference_facet_points = ref_points\n", @@ -427,7 +426,7 @@ "metadata": {}, "source": [ "Now that we have created the bilinear and linear form, and the boundary conditions,\n", - "we turn to solving the problem. For this we use the `dolfinx.fem.petsc.LinearProblem` class.\n", + "we turn to solving the problem. For this we use the {py:class}`dolfinx.fem.petsc.LinearProblem` class.\n", "As opposed to the previous examples, we now have an explicit block structure, which we would like to\n", "exploit when solving the problem. However, first we will solve the problem without any preconditioner\n", "to have a baseline performance." @@ -464,7 +463,7 @@ "id": "22", "metadata": {}, "source": [ - "Note that we have specified `kind=\"mpi\"` in the initialization of the `LinearProblem`.\n", + "Note that we have specified `kind=\"mpi\"` in the initialization of the {py:class}`LinearProblem`.\n", "This is to inform DOLFINx that we want to preserve the block structure of the problem when assembling." ] }, @@ -531,7 +530,7 @@ "As mentioned earlier, there are more efficient ways of solving this problem, than using a direct solver.\n", "Especially with the saddle point structure of the problem, we can use a Schur complement preconditioner.\n", "As described in [FEniCSx PCTools: Mixed Poisson](https://rafinex-external-rifle.gitlab.io/fenicsx-pctools/demo/demo_mixed-poisson.html),\n", - "instead of wrapping the matrices in a custom wrapper, we can use `dolfinx.fem.petsc.LinearProblem` to solve the problem." + "instead of wrapping the matrices in a custom wrapper, we can use {py:class}`dolfinx.fem.petsc.LinearProblem` to solve the problem." ] }, { diff --git a/chapter4/mixed_poisson.py b/chapter4/mixed_poisson.py index def0b46e..07439c22 100644 --- a/chapter4/mixed_poisson.py +++ b/chapter4/mixed_poisson.py @@ -11,7 +11,7 @@ # --- # # Mixed Poisson with a Schur complement pre-conditioner -# This example demonstrates how to use PETSc fieldsplits with custom preconditioners in DOLFINx. +# This example demonstrates how to use {py:mod}`PETSc` fieldsplits with custom preconditioners in DOLFINx. # This example is heavily inspired by the [FEniCSx PCTools example](https://rafinex-external-rifle.gitlab.io/fenicsx-pctools/demo/demo_mixed-poisson.html) # which was presented in {cite}`mp-rehor2025pctools`. @@ -82,7 +82,7 @@ def Gamma_N(x): V = dolfinx.fem.functionspace(mesh, ("DG", 0)) Q = dolfinx.fem.functionspace(mesh, ("BDM", 1)) -# We define a `ufl.MixedFunctionSpace` to automatically handle the block structure of the problem +# We define a {py:class}`ufl.MixedFunctionSpace` to automatically handle the block structure of the problem # + import ufl @@ -110,7 +110,7 @@ def Gamma_N(x): # - We enforce the continuity of $u$ weakly by removing the jump term. # - On $\Gamma_D$, we insert the Dirichlet condition $u=u_D$. # The Dirichlet condition is thus enforced weakly, as a natural boundary condition. -# - On $\Gamma_N$, we enforce the Neumann condition $\sigma\cdot \mathbf{n}=g$ strongly, by using a `dolfinx.fem.dirichletbc` for $\sigma$. +# - On $\Gamma_N$, we enforce the Neumann condition $\sigma\cdot \mathbf{n}=g$ strongly, by using a {py:func}`dolfinx.fem.dirichletbc` for $\sigma$. # The test function $\tau$ then satisfies $\tau\cdot\mathbf{n}=0$ on $\Gamma_N$, and the integral over $\Gamma_N$ vanishes. # # Thus we end up with: @@ -163,7 +163,7 @@ def Gamma_N(x): # b_1 # \end{pmatrix} # \end{align} -# We can extract the block structure of the bilinear form using `ufl.extract_blocks`, which returns a nested list of bilinear forms. +# We can extract the block structure of the bilinear form using {py:func}`ufl.extract_blocks`, which returns a nested list of bilinear forms. # You can also build this nested list by hand if you want to, but it is usually more error-prone. # + @@ -178,14 +178,14 @@ def Gamma_N(x): # - # Next, we create the essential boundary condition $\sigma\cdot\mathbf{n}=g$ on $\Gamma_N$, -# which in DOLFINx is a `dolfinx.fem.dirichletbc` for $\sigma$. +# which in DOLFINx is a {py:func}`dolfinx.fem.dirichletbc` for $\sigma$. # The degrees of freedom of the first order BDM space are moments of the normal component $\sigma\cdot\mathbf{n}$ on each facet. # Thus, we can interpolate any vector field whose normal component is equal to $g$ on $\Gamma_N$. # We use $g\mathbf{n}$, which only requires knowledge of $g$, and not of the full flux $\sigma$ on the boundary. # As we are using manufactured solutions for this problem, we could manually derive the explicit expression # for $\sigma$ on the boundary $\Gamma_N$. # However, in general this is not possible (especially for curved boundaries), and we have to use a more generic approach. -# For this we will use the `dolfinx.fem.Expression` class to interpolate the expression into the function space $Q$. +# For this we will use the {py:class}`dolfinx.fem.Expression` class to interpolate the expression into the function space $Q$. # This is done by evaluating the expression at the physical interpolation points of the mesh. # A convenience function for this is provided in the `interpolate_facet_expression` function below. @@ -217,13 +217,15 @@ def interpolate_facet_expression( facet_cmap = basix.ufl.element( "Lagrange", facet_type, c_el.degree, shape=(domain.geometry.dim,), dtype=x_type ) + # Wrap the C++ coordinate element in the Python interface, as the signature of + # the C++ `pull_back` differs between DOLFINx versions if np.issubdtype(x_type, np.float32): - facet_cel = dolfinx.cpp.fem.CoordinateElement_float32( - facet_cmap.basix_element._e + facet_cel = dolfinx.fem.CoordinateElement( + dolfinx.cpp.fem.CoordinateElement_float32(facet_cmap.basix_element._e) ) elif np.issubdtype(x_type, np.float64): - facet_cel = dolfinx.cpp.fem.CoordinateElement_float64( - facet_cmap.basix_element._e + facet_cel = dolfinx.fem.CoordinateElement( + dolfinx.cpp.fem.CoordinateElement_float64(facet_cmap.basix_element._e) ) else: raise TypeError( @@ -237,7 +239,7 @@ def interpolate_facet_expression( interpolation_points = Q_el.basix_element.x for i, points in enumerate(interpolation_points[fdim]): geom = ref_geom[ref_top[fdim][i]] - ref_points = facet_cel.pull_back(points, geom, tol=1e-6, maxit=25) + ref_points = facet_cel.pull_back(points, geom) # Assert that interpolation points are all equal on all facets if reference_facet_points is None: reference_facet_points = ref_points @@ -298,7 +300,7 @@ def interpolate_facet_expression( ) # Now that we have created the bilinear and linear form, and the boundary conditions, -# we turn to solving the problem. For this we use the `dolfinx.fem.petsc.LinearProblem` class. +# we turn to solving the problem. For this we use the {py:class}`dolfinx.fem.petsc.LinearProblem` class. # As opposed to the previous examples, we now have an explicit block structure, which we would like to # exploit when solving the problem. However, first we will solve the problem without any preconditioner # to have a baseline performance. @@ -323,7 +325,7 @@ def interpolate_facet_expression( ) # - -# Note that we have specified `kind="mpi"` in the initialization of the `LinearProblem`. +# Note that we have specified `kind="mpi"` in the initialization of the {py:class}`LinearProblem`. # This is to inform DOLFINx that we want to preserve the block structure of the problem when assembling. import time @@ -352,7 +354,7 @@ def interpolate_facet_expression( # As mentioned earlier, there are more efficient ways of solving this problem, than using a direct solver. # Especially with the saddle point structure of the problem, we can use a Schur complement preconditioner. # As described in [FEniCSx PCTools: Mixed Poisson](https://rafinex-external-rifle.gitlab.io/fenicsx-pctools/demo/demo_mixed-poisson.html), -# instead of wrapping the matrices in a custom wrapper, we can use `dolfinx.fem.petsc.LinearProblem` to solve the problem. +# instead of wrapping the matrices in a custom wrapper, we can use {py:class}`dolfinx.fem.petsc.LinearProblem` to solve the problem. # We start by defining the $S$ matrix in the Schur complement (see the aforementioned link for details on the variational formulation).