Update membrane code with crosslinks

_config.yml

@@ -42,19 +42,20 @@ sphinx:
     #   navigation_with_keys: false
     codeautolink_concat_default: True
     intersphinx_mapping:
-        basix : ["https://docs.fenicsproject.org/basix/main/python/", null]
-        ffcx : ["https://docs.fenicsproject.org/ffcx/main/", null]
-        ufl : ["https://docs.fenicsproject.org/ufl/main/", null]
-        dolfinx : ["https://docs.fenicsproject.org/dolfinx/main/python", null]
-        petsc4py: ["https://petsc.org/release/petsc4py", null]
-        mpi4py: ["https://mpi4py.readthedocs.io/en/stable", null]
-        numpy: ["https://numpy.org/doc/stable/", null]
-        pyvista: ["https://docs.pyvista.org/", null]
-        packaging: ["https://packaging.pypa.io/en/stable/", null]
+      basix: ["https://docs.fenicsproject.org/basix/main/python/", null]
+      ffcx: ["https://docs.fenicsproject.org/ffcx/main/", null]
+      ufl: ["https://docs.fenicsproject.org/ufl/main/", null]
+      dolfinx: ["https://docs.fenicsproject.org/dolfinx/main/python", null]
+      petsc4py: ["https://petsc.org/release/petsc4py", null]
+      mpi4py: ["https://mpi4py.readthedocs.io/en/stable", null]
+      numpy: ["https://numpy.org/doc/stable/", null]
+      pyvista: ["https://docs.pyvista.org/", null]
+      packaging: ["https://packaging.pypa.io/en/stable/", null]
+      matplotlib: ["https://matplotlib.org/stable/", null]
 
   extra_extensions:
-    - 'sphinx.ext.autodoc'
-    - 'sphinx.ext.intersphinx'
+    - "sphinx.ext.autodoc"
+    - "sphinx.ext.intersphinx"
     - "sphinx_codeautolink"
 
 parse:

chapter1/membrane_code.py

@@ -19,58 +19,70 @@
 # In this section, we will solve the deflection of the membrane problem.
 # After finishing this section, you should be able to:
 # - Create a simple mesh using the GMSH Python API and load it into DOLFINx
-# - Create constant boundary conditions using a geometrical identifier
-# - Use `ufl.SpatialCoordinate` to create a spatially varying function
-# - Interpolate a `ufl.Expression` into an appropriate function space
-# - Evaluate a `dolfinx.fem.Function` at any point $x$
+# - Create constant boundary conditions using a {py:func}`geometrical identifier<dolfinx.fem.locate_dofs_geometrical>`
+# - Use {py:class}`ufl.SpatialCoordinate` to create a spatially varying function
+# - Interpolate a {py:class}`ufl-Expression<ufl.core.expr.Expr>` into an appropriate function space
+# - Evaluate a {py:class}`dolfinx.fem.Function` at any point $x$
 # - Use Paraview to visualize the solution of a PDE
 #
 # ## Creating the mesh
 #
-# To create the computational geometry, we use the Python-API of [GMSH](https://gmsh.info/). We start by importing the gmsh-module and initializing it.
+# To create the computational geometry, we use the Python-API of [GMSH](https://gmsh.info/).
+# We start by importing the gmsh-module and initializing it.
 
 # +
 import gmsh
 
 gmsh.initialize()
 # -
 
-# The next step is to create the membrane and start the computations by the GMSH CAD kernel, to generate the relevant underlying data structures. The first arguments of `addDisk` are the x, y and z coordinate of the center of the circle, while the two last arguments are the x-radius and y-radius.
+# The next step is to create the membrane and start the computations by the GMSH CAD kernel,
+# to generate the relevant underlying data structures.
+# The first arguments of `addDisk` are the x, y and z coordinate of the center of the circle,
+# while the two last arguments are the x-radius and y-radius.
 
 membrane = gmsh.model.occ.addDisk(0, 0, 0, 1, 1)
 gmsh.model.occ.synchronize()
 
-# After that, we make the membrane a physical surface, such that it is recognized by `gmsh` when generating the mesh. As a surface is a two-dimensional entity, we add `2` as the first argument, the entity tag of the membrane as the second argument, and the physical tag as the last argument. In a later demo, we will get into when this tag matters.
+# After that, we make the membrane a physical surface, such that it is recognized by `gmsh` when generating the mesh.
+# As a surface is a two-dimensional entity, we add `2` as the first argument,
+# the entity tag of the membrane as the second argument, and the physical tag as the last argument.
+# In a later demo, we will get into when this tag matters.
 
 gdim = 2
 gmsh.model.addPhysicalGroup(gdim, [membrane], 1)
 
-# Finally, we generate the two-dimensional mesh. We set a uniform mesh size by modifying the GMSH options.
+# Finally, we generate the two-dimensional mesh.
+# We set a uniform mesh size by modifying the GMSH options.
 
 gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.05)
 gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.05)
 gmsh.model.mesh.generate(gdim)
 
 # # Interfacing with GMSH in DOLFINx
-# We will import the GMSH-mesh directly from GMSH into DOLFINx via the `dolfinx.io.gmsh` interface.
+# We will import the GMSH-mesh directly from GMSH into DOLFINx via the {py:mod}`dolfinx.io.gmsh` interface.
 # The {py:mod}`dolfinx.io.gmsh` module contains two functions
 # 1. {py:func}`model_to_mesh<dolfinx.io.gmsh.model_to_mesh>` which takes in a `gmsh.model`
 #   and returns a {py:class}`dolfinx.io.gmsh.MeshData` object.
 # 2. {py:func}`read_from_msh<dolfinx.io.gmsh.read_from_msh>` which takes in a path to a `.msh`-file
 #  and returns a {py:class}`dolfinx.io.gmsh.MeshData` object.
 #
-# The `MeshData` object will contain a `dolfinx.mesh.Mesh`, under the attribute `mesh`.
+# The {py:class}`MeshData` object will contain a {py:class}`dolfinx.mesh.Mesh`,
+# under the attribute {py:attr}`mesh<dolfinx.io.gmsh.MeshData.mesh>`.
 # This mesh will contain all GMSH Physical Groups of the highest topological dimension.
 # ```{note}
 # If you do not use `gmsh.model.addPhysicalGroup` when creating the mesh with GMSH, it can not be read into DOLFINx.
 # ```
-# The `MeshData` object can also contain tags for all other `PhysicalGroups` that has been added to the mesh,
-# that being `vertex_tags`, `edge_tags`, `facet_tags` and `cell_tags`.
+# The {py:class}`MeshData<dolfinx.io.gmsh.MeshData>` object can also contain tags for
+# all other `PhysicalGroups` that has been added to the mesh, that being
+# {py:attr}`cell_tags<dolfinx.io.gmsh.MeshData.cell_tags>`, {py:attr}`facet_tags<dolfinx.io.gmsh.MeshData.facet_tags>`,
+# {py:attr}`ridge_tags<dolfinx.io.gmsh.MeshData.ridge_tags>` and
+# {py:attr}`peak_tags<dolfinx.io.gmsh.MeshData.peak_tags>`.
 # To read either `gmsh.model` or a `.msh`-file, one has to distribute the mesh to all processes used by DOLFINx.
 # As GMSH does not support mesh creation with MPI, we currently have a `gmsh.model.mesh` on each process.
 # To distribute the mesh, we have to specify which process the mesh was created on,
 # and which communicator rank should distribute the mesh.
-# The `model_to_mesh` will then load the mesh on the specified rank,
+# The {py:func}`model_to_mesh<dolfinx.io.gmsh.model_to_mesh>` will then load the mesh on the specified rank,
 # and distribute it to the communicator using a mesh partitioner.
 
 # +
@@ -95,18 +107,25 @@
 # -
 
 # ## Defining a spatially varying load
-# The right hand side pressure function is represented using `ufl.SpatialCoordinate` and two constants, one for $\beta$ and one for $R_0$.
+# The right hand side pressure function is represented using {py:class}`ufl.SpatialCoordinate` and two constants,
+# one for $\beta$ and one for $R_0$.
 
+# +
 import ufl
 from dolfinx import default_scalar_type
 
 x = ufl.SpatialCoordinate(domain)
 beta = fem.Constant(domain, default_scalar_type(12))
 R0 = fem.Constant(domain, default_scalar_type(0.3))
 p = 4 * ufl.exp(-(beta**2) * (x[0] ** 2 + (x[1] - R0) ** 2))
+# -
 
 # ## Create a Dirichlet boundary condition using geometrical conditions
-# The next step is to create the homogeneous boundary condition. As opposed to the [first tutorial](./fundamentals_code.ipynb) we will use `dolfinx.fem.locate_dofs_geometrical` to locate the degrees of freedom on the boundary. As we know that our domain is a circle with radius 1, we know that any degree of freedom should be located at a coordinate $(x,y)$ such that $\sqrt{x^2+y^2}=1$.
+# The next step is to create the homogeneous boundary condition.
+# As opposed to the [first tutorial](./fundamentals_code.ipynb) we will use
+# {py:func}`locate_dofs_geometrical<dolfinx.fem.locate_dofs_geometrical>` to locate the degrees of freedom on the boundary.
+# As we know that our domain is a circle with radius 1, we know that any degree of freedom should be
+# located at a coordinate $(x,y)$ such that $\sqrt{x^2+y^2}=1$.
 
 # +
 import numpy as np
@@ -119,7 +138,9 @@ def on_boundary(x):
 boundary_dofs = fem.locate_dofs_geometrical(V, on_boundary)
 # -
 
-# As our Dirichlet condition is homogeneous (`u=0` on the whole boundary), we can initialize the `dolfinx.fem.dirichletbc` with a constant value, the degrees of freedom and the function space to apply the boundary condition on.
+# As our Dirichlet condition is homogeneous (`u=0` on the whole boundary), we can initialize the
+# {py:class}`dolfinx.fem.DirichletBC` with a constant value, the degrees of freedom and the function
+# space to apply the boundary condition on. We use the constructor {py:func}`dolfinx.fem.dirichletbc`.
 
 bc = fem.dirichletbc(default_scalar_type(0), boundary_dofs, V)
 
@@ -139,8 +160,13 @@ def on_boundary(x):
 )
 uh = problem.solve()
 
-# ## Interpolation of a `ufl`-expression
-# As we previously defined the load `p` as a spatially varying function, we would like to interpolate this function into an appropriate function space for visualization. To do this we use the `dolfinx.Expression`. The expression takes in any `ufl`-expression, and a set of points on the reference element. We will use the interpolation points of the space we want to interpolate in to.
+# ## Interpolation of a UFL-expression
+# As we previously defined the load `p` as a spatially varying function,
+# we would like to interpolate this function into an appropriate function space for visualization.
+# To do this we use the class {py:class}`Expression<dolfinx.fem.Expression>`.
+# The expression takes in any UFL-expression, and a set of points on the reference element.
+# We will use the {py:attr}`interpolation points<dolfinx.fem.FiniteElement.interpolation_points>`
+# of the space we want to interpolate in to.
 # We choose a high order function space to represent the function `p`, as it is rapidly varying in space.
 
 Q = fem.functionspace(domain, ("Lagrange", 5))
@@ -154,8 +180,7 @@ def on_boundary(x):
 from dolfinx.plot import vtk_mesh
 import pyvista
 
-
-# Extract topology from mesh and create pyvista mesh
+# Extract topology from mesh and create {py:class}`pyvista.UnstructuredGrid`
 
 topology, cell_types, x = vtk_mesh(V)
 grid = pyvista.UnstructuredGrid(topology, cell_types, x)
@@ -190,7 +215,8 @@ def on_boundary(x):
 
 # ## Making curve plots throughout the domain
 # Another way to compare the deflection and the load is to make a plot along the line $x=0$.
-# This is just a matter of defining a set of points along the $y$-axis and evaluating the finite element functions $u$ and $p$ at these points.
+# This is just a matter of defining a set of points along the $y$-axis and evaluating the
+# finite element functions $u$ and $p$ at these points.
 
 tol = 0.001  # Avoid hitting the outside of the domain
 y = np.linspace(-1 + tol, 1 - tol, 101)
@@ -199,22 +225,40 @@ def on_boundary(x):
 u_values = []
 p_values = []
 
-# As a finite element function is the linear combination of all degrees of freedom, $u_h(x)=\sum_{i=1}^N c_i \phi_i(x)$ where $c_i$ are the coefficients of $u_h$ and $\phi_i$ is the $i$-th basis function, we can compute the exact solution at any point in $\Omega$.
-# However, as a mesh consists of a large set of degrees of freedom (i.e. $N$ is large), we want to reduce the number of evaluations of the basis function $\phi_i(x)$. We do this by identifying which cell of the mesh $x$ is in.
-# This is efficiently done by creating a bounding box tree of the cells of the mesh, allowing a quick recursive search through the mesh entities.
+# As a finite element function is the linear combination of all degrees of freedom,
+# $u_h(x)=\sum_{i=1}^N c_i \phi_i(x)$ where $c_i$ are the coefficients of $u_h$ and $\phi_i$
+# is the $i$-th basis function, we can compute the exact solution at any point in $\Omega$.
+# However, as a mesh consists of a large set of degrees of freedom (i.e. $N$ is large),
+# we want to reduce the number of evaluations of the basis function $\phi_i(x)$.
+# We do this by identifying which cell of the mesh $x$ is in.
+# This is efficiently done by creating a {py:class}`bounding box tree<dolfinx.geometry.BoundingBoxTree`
+# of the cells of the mesh,
+# allowing a quick recursive search through the mesh entities.
 
 # +
 from dolfinx import geometry
 
 bb_tree = geometry.bb_tree(domain, domain.topology.dim)
 # -
 
-# Now we can compute which cells the bounding box tree collides with using `dolfinx.geometry.compute_collisions_points`. This function returns a list of cells whose bounding box collide for each input point. As different points might have different number of cells, the data is stored in `dolfinx.cpp.graph.AdjacencyList_int32`, where one can access the cells for the `i`th point by calling `links(i)`.
-# However, as the bounding box of a cell spans more of $\mathbb{R}^n$ than the actual cell, we check that the actual cell collides with the input point
-# using `dolfinx.geometry.select_colliding_cells`, which measures the exact distance between the point and the cell (approximated as a convex hull for higher order geometries).
-# This function also returns an adjacency-list, as the point might align with a facet, edge or vertex that is shared between multiple cells in the mesh.
+# Now we can compute which cells the bounding box tree collides with using
+# {py:func}`dolfinx.geometry.compute_collisions_points`.
+# This function returns a list of cells whose bounding box collide for each input point.
+# As different points might have different number of cells, the data is stored in
+# {py:class}`dolfinx.graph.AdjacencyList`, where one can access the cells for the
+# `i`th point by calling `links(i)<dolfinx.graph.AdjacencyList.links>`.
+# However, as the bounding box of a cell spans more of $\mathbb{R}^n$ than the actual cell,
+# we check that the actual cell collides with the input point using
+# {py:func}`dolfinx.geometry.compute_colliding_cells`,
+# which measures the exact distance between the point and the cell
+# (approximated as a convex hull for higher order geometries).
+# This function also returns an adjacency-list, as the point might align with a facet,
+# edge or vertex that is shared between multiple cells in the mesh.
 #
-# Finally, we would like the code below to run in parallel, when the mesh is distributed over multiple processors. In that case, it is not guaranteed that every point in `points` is on each processor. Therefore we create a subset `points_on_proc` only containing the points found on the current processor.
+# Finally, we would like the code below to run in parallel,
+# when the mesh is distributed over multiple processors.
+# In that case, it is not guaranteed that every point in `points` is on each processor.
+# Therefore we create a subset `points_on_proc` only containing the points found on the current processor.
 
 cells = []
 points_on_proc = []
@@ -227,14 +271,16 @@ def on_boundary(x):
         points_on_proc.append(point)
         cells.append(colliding_cells.links(i)[0])
 
-# We now have a list of points on the processor, on in which cell each point belongs. We can then call `uh.eval` and `pressure.eval` to obtain the set of values for all the points.
-#
+# We now have a list of points on the processor, on in which cell each point belongs.
+# We can then call {py:meth}`uh.eval<dolfinx.fem.Function.eval>` and
+# {py:meth}`pressure.eval<dolfinx.fem.Function.eval>` to obtain the set of values for all the points.
 
 points_on_proc = np.array(points_on_proc, dtype=np.float64)
 u_values = uh.eval(points_on_proc, cells)
 p_values = pressure.eval(points_on_proc, cells)
 
-# As we now have an array of coordinates and two arrays of function values, we can use `matplotlib` to plot them
+# As we now have an array of coordinates and two arrays of function values,
+# we can use {py:mod}`matplotlib<matplotlib.pyplot>` to plot them
 
 # +
 import matplotlib.pyplot as plt
@@ -253,7 +299,7 @@ def on_boundary(x):
 plt.legend()
 # -
 
-# If executed in parallel as a python file, we save a plot per processor
+# If executed in parallel as a Python file, we save a plot per processor
 
 plt.savefig(f"membrane_rank{MPI.COMM_WORLD.rank:d}.png")
 

chapter2/amr.ipynb

@@ -138,7 +138,6 @@
    },
    "outputs": [],
    "source": [
-    "\n",
     "grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh))\n",
     "grid.cell_data[\"ct\"] = ct.values\n",
     "\n",

chapter2/amr.py

@@ -72,7 +72,6 @@
 # We use pyvista to visualize the mesh.
 
 # + tags=["hide-input"]
-
 grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh))
 grid.cell_data["ct"] = ct.values
 

chapter2/helmholtz_code.ipynb

@@ -95,7 +95,7 @@
     "    assert mesh_data.facet_tags is not None\n",
     "    facet_tags = mesh_data.facet_tags\n",
     "else:\n",
-    "    domain, _, facet_tags = mesh_data\n"
+    "    domain, _, facet_tags = mesh_data"
    ]
   },
   {
@@ -113,15 +113,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "\n",
     "V = fem.functionspace(domain, (\"Lagrange\", 1))\n",
     "\n",
     "# Discrete frequency range\n",
     "freq = np.arange(10, 1000, 5)  # Hz\n",
     "\n",
     "# Air parameters\n",
     "rho0 = 1.225  # kg/m^3\n",
-    "c = 340  # m/s\n"
+    "c = 340  # m/s"
    ]
   },
   {

chapter2/helmholtz_code.py

@@ -89,13 +89,11 @@
     facet_tags = mesh_data.facet_tags
 else:
     domain, _, facet_tags = mesh_data
-
 # -
 
 # We define the function space for our unknown $p$ and define the range of frequencies we want to solve the Helmholtz equation for.
 
 # +
-
 V = fem.functionspace(domain, ("Lagrange", 1))
 
 # Discrete frequency range
@@ -104,7 +102,6 @@
 # Air parameters
 rho0 = 1.225  # kg/m^3
 c = 340  # m/s
-
 # -
 
 # ## Boundary conditions

chapter2/nonlinpoisson_code.ipynb

@@ -28,9 +28,8 @@
     "import numpy\n",
     "\n",
     "from mpi4py import MPI\n",
-    "from petsc4py import PETSc\n",
     "\n",
-    "from dolfinx import mesh, fem, log\n",
+    "from dolfinx import mesh, fem\n",
     "from dolfinx.fem.petsc import NonlinearProblem\n",
     "\n",
     "\n",

chapter2/nonlinpoisson_code.py

@@ -31,9 +31,8 @@
 import numpy
 
 from mpi4py import MPI
-from petsc4py import PETSc
 
-from dolfinx import mesh, fem, log
+from dolfinx import mesh, fem
 from dolfinx.fem.petsc import NonlinearProblem
 
 

chapter2/ns_code1.ipynb

@@ -501,7 +501,9 @@
    "cell_type": "code",
    "execution_count": null,
    "id": "26",
-   "metadata": {},
+   "metadata": {
+    "lines_to_end_of_cell_marker": 2
+   },
    "outputs": [],
    "source": [
     "from pathlib import Path\n",
@@ -516,7 +518,9 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "lines_to_next_cell": 2
+   },
    "source": [
     "We also interpolate the analytical solution into our function-space and create a variational formulation for the $L^2$-error.\n"
    ]
@@ -527,8 +531,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "\n",
-    "\n",
     "def u_exact(x):\n",
     "    values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)\n",
     "    values[0] = 4 * x[1] * (1.0 - x[1])\n",