Update all notebooks to render correctly on web page and in notebooks.

python/demo/demo_axis/demo_axis.py

@@ -99,10 +99,12 @@
 # cylindrical harmonics:
 #
 # $$
+# \begin{align}
 # \mathbf{E}_s(\rho, z, \phi) &= \sum_m\mathbf{E}^{(m)}_s(\rho, z)e^{-jm\phi} \\
 # \mathbf{E}_b(\rho, z, \phi) &= \sum_m\mathbf{E}^{(m)}_b(\rho, z)e^{-jm\phi} \\
 # \bar{\mathbf{v}}(\rho, z, \phi) &=
 # \sum_m\bar{\mathbf{v}}^{(m)}(\rho, z)e^{+jm\phi}
+# \end{align}
 # $$
 #
 # The curl operator $\nabla\times$ in cylindrical coordinates becomes:
@@ -230,11 +232,13 @@ def background_field_p(theta: float, n_bkg: float, k0: float, m: int, x):
 # PML:
 #
 # $$
+# \begin{align}
 # &\rho^{\prime} = \rho\left[1 +j \alpha/k_0 \left(\frac{r
 # - r_{dom}}{r~r_{pml}}\right)\right] \\
 # &z^{\prime} = z\left[1 +j \alpha/k_0 \left(\frac{r
 # - r_{dom}}{r~r_{pml}}\right)\right] \\
 # &\phi^{\prime} = \phi
+# \end{align}
 # $$
 #
 # with $\alpha$ tuning the absorption inside the PML, and $r =
@@ -265,10 +269,12 @@ def background_field_p(theta: float, n_bkg: float, k0: float, m: int, x):
 # ${\boldsymbol{\mu}_{pml}}$:
 #
 # $$
+# \begin{align}
 # & {\boldsymbol{\varepsilon}_{pml}} =
 # A^{-1} \mathbf{A} {\boldsymbol{\varepsilon}_b}\mathbf{A}^{T}\\
 # & {\boldsymbol{\mu}_{pml}} =
 # A^{-1} \mathbf{A} {\boldsymbol{\mu}_b}\mathbf{A}^{T}
+# \end{align}
 # $$
 #
 # For doing these calculations, we define the `pml_coordinate` and
@@ -414,17 +420,21 @@ def create_eps_mu(pml, rho, eps_bkg, mu_bkg):
 # to few harmonic numbers, e.g., $m = -1, 0, 1$. Besides, we have that:
 #
 # $$
+# \begin{align}
 # &J_{-m}=(-1)^m J_m \\
 # &J_{-m}^{\prime}=(-1)^m J_m^{\prime} \\
 # &j^{-m}=(-1)^m j^m
+# \end{align}
 # $$
 #
 # and therefore:
 #
 # $$
+# \begin{align}
 # &E_{b, \rho}^{(m)}=E_{b, \rho}^{(-m)} \\
 # &E_{b, \phi}^{(m)}=-E_{b, \phi}^{(-m)} \\
 # &E_{b, z}^{(m)}=E_{b, z}^{(-m)}
+# \end{align}
 # $$
 #
 # In light of this, we can solve the problem for $m\geq 0$.
@@ -479,9 +489,11 @@ def create_eps_mu(pml, rho, eps_bkg, mu_bkg):
 # following way:
 #
 # $$
+# \begin{align}
 # &E_{s, \rho}^{(m)}(\phi)=E_{s, \rho}^{(m)}(e^{-jm\phi}+e^{jm\phi}) \\
 # &E_{s, \phi}^{(m)}(\phi)=E_{s, \phi}^{(m)}(e^{-jm\phi}-e^{jm\phi}) \\
 # &E_{s, z}^{(m)}(\phi)=E_{s, z}^{(m)}(e^{-jm\phi}+e^{jm\phi})
+# \end{align}
 # $$
 #
 # For this reason, we also add a `phase` constant for the above phase

python/demo/demo_biharmonic.py

@@ -65,8 +65,10 @@
 # and considering the boundary conditions
 #
 # $$
+# \begin{align}
 # u &= 0 \quad {\rm on} \ \partial\Omega, \\
 # \nabla^{2} u &= 0 \quad {\rm on} \ \partial\Omega,
+# \end{align}
 # $$
 #
 # a weak formulation of the biharmonic problem reads: find $u \in V$ such that

python/demo/demo_cahn-hilliard.py

@@ -35,13 +35,15 @@
 # derivatives.  The equation reads:
 #
 # $$
+# \begin{align}
 # \frac{\partial c}{\partial t} -
 #   \nabla \cdot M \left(\nabla\left(\frac{d f}{dc}
 #   - \lambda \nabla^{2}c\right)\right) &= 0 \quad {\rm in} \ \Omega, \\
 # M\left(\nabla\left(\frac{d f}{d c} -
 #   \lambda \nabla^{2}c\right)\right) \cdot n
 #   &= 0 \quad {\rm on} \ \partial\Omega, \\
 # M \lambda \nabla c \cdot n &= 0 \quad {\rm on} \ \partial\Omega.
+# \end{align}
 # $$
 #
 # where $c$ is the unknown field, the function $f$ is usually non-convex
@@ -57,21 +59,25 @@
 # as two coupled second-order equations:
 #
 # $$
+# \begin{align}
 # \frac{\partial c}{\partial t} - \nabla \cdot M \nabla\mu
 #     &= 0 \quad {\rm in} \ \Omega, \\
 # \mu -  \frac{d f}{d c} + \lambda \nabla^{2}c &= 0 \quad {\rm in} \ \Omega.
+# \end{align}
 # $$
 #
 # The unknown fields are now $c$ and $\mu$. The weak (variational) form
 # of the problem reads: find $(c, \mu) \in V \times V$ such that
 #
 # $$
+# \begin{align}
 # \int_{\Omega} \frac{\partial c}{\partial t} q \, {\rm d} x +
 #     \int_{\Omega} M \nabla\mu \cdot \nabla q \, {\rm d} x
 #     &= 0 \quad \forall \ q \in V,  \\
 # \int_{\Omega} \mu v \, {\rm d} x - \int_{\Omega} \frac{d f}{d c} v \, {\rm d} x
 #   - \int_{\Omega} \lambda \nabla c \cdot \nabla v \, {\rm d} x
 #    &= 0 \quad \forall \ v \in V.
+# \end{align}
 # $$
 #
 # ### Time discretisation
@@ -81,12 +87,14 @@
 # equation:
 #
 # $$
+# \begin{align}
 # \int_{\Omega} \frac{c_{n+1} - c_{n}}{dt} q \, {\rm d} x
 # + \int_{\Omega} M \nabla \mu_{n+\theta} \cdot \nabla q \, {\rm d} x
 #        &= 0 \quad \forall \ q \in V  \\
 # \int_{\Omega} \mu_{n+1} v  \, {\rm d} x - \int_{\Omega} \frac{d f_{n+1}}{d c} v  \, {\rm d} x
 # - \int_{\Omega} \lambda \nabla c_{n+1} \cdot \nabla v \, {\rm d} x
 #        &= 0 \quad \forall \ v \in V
+# \end{align}
 # $$
 #
 # where $dt = t_{n+1} - t_{n}$ and $\mu_{n+\theta} = (1-\theta) \mu_{n} + \theta \mu_{n+1}$.

python/demo/demo_mixed-poisson.py

@@ -16,8 +16,13 @@
 # * Use mixed and non-continuous finite element spaces.
 # * Set essential boundary conditions for subspaces and $H(\mathrm{div})$ spaces.
 #
-# {download}`Python script <./demo_mixed-poisson.py>`.\
-# {download}`Jupyter notebook <./demo_mixed-poisson.ipynb>`.
+
+# ```{admonition} Download sources
+# :class: download
+#
+# * {download}`Python script <./demo_mixed-poisson.py>`
+# * {download}`Jupyter notebook <./demo_mixed-poisson.ipynb>`
+# ```
 #
 # ## Equation and problem definition
 #
@@ -27,8 +32,10 @@
 # then read
 #
 # $$
+# \begin{align}
 #   \sigma - \nabla u &= 0 \quad {\rm in} \ \Omega, \\
 #   \nabla \cdot \sigma &= - f \quad {\rm in} \ \Omega,
+# \end{align}
 # $$
 # with boundary conditions
 #
@@ -55,11 +62,13 @@
 # where the variational forms $a$ and $L$ are defined as
 #
 # $$
+# \begin{align}
 #   a((\sigma, u), (\tau, v)) &=
 #     \int_{\Omega} \sigma \cdot \tau + \nabla \cdot \tau \ u
 #   + \nabla \cdot \sigma \ v \ {\rm d} x, \\
 #   L((\tau, v)) &= - \int_{\Omega} f v \ {\rm d} x
 #   + \int_{\Gamma_D} u_0 \tau \cdot n  \ {\rm d} s,
+# \end{align}
 # $$
 # and $\Sigma_g = \{ \tau \in H({\rm div})$ such that $\tau \cdot n|_{\Gamma_N}
 # = g \}$ and $V = L^2(\Omega)$.

python/demo/demo_navier-stokes.py

@@ -28,8 +28,10 @@
 # $(0, \infty)$, given by
 #
 # $$
-# \partial_t u - \nu \Delta u + (u \cdot \nabla)u + \nabla p &= f \textnormal{ in } \Omega_t, \\
-# \nabla \cdot u &= 0 \textnormal{ in } \Omega_t,
+# \begin{align}
+# \partial_t u - \nu \Delta u + (u \cdot \nabla)u + \nabla p &= f \text{ in } \Omega_t, \\
+# \nabla \cdot u &= 0 \text{ in } \Omega_t,
+# \end{align}
 # $$
 #
 # where $u: \Omega_t \to \mathbb{R}^d$ is the velocity field,
@@ -40,13 +42,13 @@
 # The problem is supplemented with the initial condition
 #
 # $$
-#     u(x, 0) = u_0(x) \textnormal{ in } \Omega
+#     u(x, 0) = u_0(x) \text{ in } \Omega
 # $$
 #
 # and boundary condition
 #
 # $$
-#     u = u_D \textnormal{ on } \partial \Omega \times (0, \infty),
+#     u = u_D \text{ on } \partial \Omega \times (0, \infty),
 # $$
 #
 # where $u_0: \Omega \to \mathbb{R}^d$ is a prescribed initial velocity field
@@ -64,11 +66,13 @@
 # We begin by introducing the function spaces
 #
 # $$
-#     V_h^g &:= \left\{v \in H(\textnormal{div}; \Omega);
+# \begin{align}
+#     V_h^g &:= \left\{v \in H(\text{div}; \Omega);
 #     v|_K \in V_h(K) \; \forall K \in \mathcal{T}, v \cdot n = g \cdot n
-#     \textnormal{ on } \partial \Omega \right\} \\
+#     \text{ on } \partial \Omega \right\} \\
 #     Q_h &:= \left\{q \in L^2_0(\Omega);
 #     q|_K \in Q_h(K) \; \forall K \in \mathcal{T} \right\}.
+# \end{align}
 # $$
 #
 # The local spaces $V_h(K)$ and $Q_h(K)$ should satisfy
@@ -80,14 +84,14 @@
 # affine simplex cells include
 #
 # $$
-#     V_h(K) := \mathbb{RT}_k(K) \textnormal{ and }
+#     V_h(K) := \mathbb{RT}_k(K) \text{ and }
 #     Q_h(K) := \mathbb{P}_k(K),
 # $$
 #
 # or
 #
 # $$
-#     V_h(K) := \mathbb{BDM}_k(K) \textnormal{ and }
+#     V_h(K) := \mathbb{BDM}_k(K) \text{ and }
 #     Q_h(K) := \mathbb{P}_{k-1}(K).
 # $$
 #
@@ -127,27 +131,31 @@
 # $(u_h, p_h) \in V_h^{u_D} \times Q_h$ such that
 #
 # $$
-#     \int_\Omega \partial_t u_h \cdot v + a_h(u_h, v_h) + c_h(u_h; u_h, v_h)
-#     + b_h(v_h, p_h) &= \int_\Omega f \cdot v_h + L_{a_h}(v_h) + L_{c_h}(v_h)
-#      \quad \forall v_h \in V_h^0, \\
-#     b_h(u_h, q_h) &= 0 \quad \forall q_h \in Q_h,
+# \begin{align}
+#   \int_\Omega \partial_t u_h \cdot v + a_h(u_h, v_h) + c_h(u_h; u_h, v_h)
+#   + b_h(v_h, p_h) &= \int_\Omega f \cdot v_h + L_{a_h}(v_h) + L_{c_h}(v_h)
+#   \quad \forall v_h \in V_h^0, \\
+#   b_h(u_h, q_h) &= 0 \quad \forall q_h \in Q_h,
+# \end{align}
 # $$
 #
 # where
 # $\renewcommand{\sumK}[0]{\sum_{K \in \mathcal{T}_h}}$
 # $\renewcommand{\sumF}[0]{\sum_{F \in \mathcal{F}_h}}$
 #
 # $$
-#     a_h(u, v) &= Re^{-1} \left(\sumK \int_K \nabla u : \nabla v
-#     - \sumF \int_F \avg{\nabla u} : \jump{v}
-#     - \sumF \int_F \avg{\nabla v} : \jump{u} \\
-#     + \sumF \int_F \frac{\alpha}{h_K} \jump{u} : \jump{v}\right), \\
-#     c_h(w; u, v) &= - \sumK \int_K u \cdot \nabla \cdot (v \otimes w)
-#     + \sumK \int_{\partial_K} w \cdot n \hat{u}^{w} \cdot v, \\
-#     L_{a_h}(v_h) &= Re^{-1} \left(- \int_{\partial \Omega} u_D \otimes n :
-#     \nabla_h v_h + \frac{\alpha}{h} u_D \otimes n : v_h \otimes n \right), \\
-#     L_{c_h}(v_h) &= - \int_{\partial \Omega} u_D \cdot n \hat{u}_D \cdot v_h, \\
-#     b_h(v, q) &= - \int_K \nabla \cdot v q.
+# \begin{align}
+#   a_h(u, v) &= Re^{-1} \left(\sumK \int_K \nabla u : \nabla v
+#   - \sumF \int_F \avg{\nabla u} : \jump{v}
+#   - \sumF \int_F \avg{\nabla v} : \jump{u} \\
+#   + \sumF \int_F \frac{\alpha}{h_K} \jump{u} : \jump{v}\right), \\
+#   c_h(w; u, v) &= - \sumK \int_K u \cdot \nabla \cdot (v \otimes w)
+#   + \sumK \int_{\partial_K} w \cdot n \hat{u}^{w} \cdot v, \\
+#   L_{a_h}(v_h) &= Re^{-1} \left(- \int_{\partial \Omega} u_D \otimes n :
+#   \nabla_h v_h + \frac{\alpha}{h} u_D \otimes n : v_h \otimes n \right), \\
+#   L_{c_h}(v_h) &= - \int_{\partial \Omega} u_D \cdot n \hat{u}_D \cdot v_h, \\
+#   b_h(v, q) &= - \int_K \nabla \cdot v q.
+# \end{align}
 # $$
 #
 #
@@ -220,7 +228,7 @@ def boundary_marker(x):
 k = 1  # Polynomial degree
 
 # Next, we create a mesh and the required functions spaces over it.
-# Since the velocity uses an $H(\textnormal{div})$-conforming function
+# Since the velocity uses an $H(\text{div})$-conforming function
 # space, we also create a vector valued discontinuous Lagrange space to
 # interpolate into for artifact free visualisation.
 

python/demo/demo_pml/demo_pml.py

@@ -208,11 +208,13 @@ def pml_coordinates(x: ufl.indexed.Indexed, alpha: float, k0: complex, l_dom: fl
 # specified here below:
 #
 # $$
-# \text{PML}_\text{corners} \rightarrow \mathbf{r}^\prime & = (x^\prime, y^\prime) \\
+# \begin{align}
+# \text{PML}_\text{corners} \rightarrow \mathbf{r}^\prime &= (x^\prime, y^\prime) \\
 # \text{PML}_\text{rectangles along x} \rightarrow
-#                                       \mathbf{r}^\prime & = (x^\prime, y) \\
+#                                       \mathbf{r}^\prime &= (x^\prime, y) \\
 # \text{PML}_\text{rectangles along y} \rightarrow
-#                                       \mathbf{r}^\prime & = (x, y^\prime).
+#                                       \mathbf{r}^\prime &= (x, y^\prime).
+# \end{align}
 # $$
 #
 # Now we define some other problem specific parameters:
@@ -335,10 +337,12 @@ def pml_coordinates(x: ufl.indexed.Indexed, alpha: float, k0: complex, l_dom: fl
 # https://www.tandfonline.com/doi/abs/10.1080/09500349608232782):
 #
 # $$
+# \begin{align}
 # {\boldsymbol{\varepsilon}_{pml}} &=
 # A^{-1} \mathbf{A} {\boldsymbol{\varepsilon}_b}\mathbf{A}^{T},\\
 # {\boldsymbol{\mu}_{pml}} &=
 # A^{-1} \mathbf{A} {\boldsymbol{\mu}_b}\mathbf{A}^{T},
+# \end{align}
 # $$
 #
 # with $A^{-1}=\operatorname{det}(\mathbf{J})$.
@@ -388,7 +392,7 @@ def create_eps_mu(pml: ufl.tensors.ListTensor,
 # while in the rest of the domain is:
 #
 # $$
-# & \int_{\Omega_m\cup\Omega_b}-(\nabla \times \mathbf{E}_s)
+# \int_{\Omega_m\cup\Omega_b}-(\nabla \times \mathbf{E}_s)
 # \cdot (\nabla \times \bar{\mathbf{v}})+\varepsilon_{r} k_{0}^{2}
 # \mathbf{E}_s \cdot \bar{\mathbf{v}}+k_{0}^{2}\left(\varepsilon_{r}
 # -\varepsilon_b\right)\mathbf{E}_b \cdot \bar{\mathbf{v}}~\mathrm{d}x.

python/demo/demo_poisson.py

@@ -23,9 +23,11 @@
 # particular boundary conditions reads:
 #
 # $$
-# - \nabla^{2} u &= f \quad {\rm in} \ \Omega, \\
-# u &= 0 \quad {\rm on} \ \Gamma_{D}, \\
-# \nabla u \cdot n &= g \quad {\rm on} \ \Gamma_{N}. \\
+# \begin{align}
+#   - \nabla^{2} u &= f \quad {\rm in} \ \Omega, \\
+#   u &= 0 \quad {\rm on} \ \Gamma_{D}, \\
+#   \nabla u \cdot n &= g \quad {\rm on} \ \Gamma_{N}. \\
+# \end{align}
 # $$
 #
 # where $f$ and $g$ are input data and $n$ denotes the outward directed
@@ -39,8 +41,10 @@
 # where $V$ is a suitable function space and
 #
 # $$
-# a(u, v) &:= \int_{\Omega} \nabla u \cdot \nabla v \, {\rm d} x, \\
-# L(v)    &:= \int_{\Omega} f v \, {\rm d} x + \int_{\Gamma_{N}} g v \, {\rm d} s.
+# \begin{align}
+#   a(u, v) &:= \int_{\Omega} \nabla u \cdot \nabla v \, {\rm d} x, \\
+#   L(v)    &:= \int_{\Omega} f v \, {\rm d} x + \int_{\Gamma_{N}} g v \, {\rm d} s.
+# \end{align}
 # $$
 #
 # The expression $a(u, v)$ is the bilinear form and $L(v)$