diff --git a/hog/forms_facets.py b/hog/forms_facets.py
index daff842ae003bd86ca62323cbeb3d6c10665fab6..47292dad5766e463b81536529de5304aaa0c67de 100644
--- a/hog/forms_facets.py
+++ b/hog/forms_facets.py
@@ -179,7 +179,9 @@ def stokes_p0_stabilization(
     r"""
     Interface integrals for the stabilization term
 
-        \gamma |\meas{e}| ( [[p_h]]_e, [[q_h]]_e )_e
+    .. math::
+
+        \gamma |\text{meas}(e)| ( [[p_h]]_e, [[q_h]]_e )_e
 
     for the P1-P0 Stokes discretization as found in
 
@@ -189,13 +191,18 @@ def stokes_p0_stabilization(
 
     Interface integrals are performed from both sides. interface_type has to be
 
-        'inner' coupling of element unknowns with themselves
+    'inner' coupling of element unknowns with themselves
+
+     .. math::
+
+         \int_e \gamma \cdot \text{meas}(e) p q
 
-            \int_e \gamma \meas{e} p q
+    'outer' coupling of element unknowns with those opposite of the interface
 
-        'outer' coupling of element unknowns with those opposite of the interface
+     .. math::
+
+         - \int_e \gamma \cdot \text{meas}(e) p q
 
-            - \int_e \gamma \meas{e} p q
     """
 
     if not isinstance(blending, IdentityMap):
@@ -306,24 +313,29 @@ def diffusion_sip_facet(
     All "interface types" are handled by this function. The type has to be specified by the interface_type argument:
     Possible types:
 
-        'inner': coupling of element unknowns with themselves (both, v and w restricted to E1)
+    'inner': coupling of element unknowns with themselves (both, v and w restricted to E1)
+
+    .. math::
 
-            - 0.5 * \int_e J_v^{-\top} \nabla v \cdot n w
-            - 0.5 * \int_e J_w^{-\top} \nabla w \cdot n v
-            + \frac{\sigma}{\meas(e)^\beta} \int_e v w
+        - 0.5 \cdot \int_e J_v^{-\top} \nabla v \cdot n w
+        - 0.5 \cdot \int_e J_w^{-\top} \nabla w \cdot n v
+        + \frac{\sigma}{\text{meas}(e)^\beta} \int_e v w
 
-        'outer': coupling of element unknowns with those opposite of the interface (v restricted to E1, w restricted to
-                 E2)
+    'outer': coupling of element unknowns with those opposite of the interface (v restricted to E1, w restricted to E2)
 
-            + 0.5 * \int_e J_v^{-\top} \nabla v \cdot n w
-            - 0.5 * \int_e J_w^{-\top} \nabla w \cdot n v
-            - \frac{\sigma}{\meas(e)^\beta} \int_e v w
+    .. math::
 
-        'dirichlet': interface is at Dirichlet boundary (both, v and w restricted to E1)
+        + 0.5 \cdot \int_e J_v^{-\top} \nabla v \cdot n w
+        - 0.5 \cdot \int_e J_w^{-\top} \nabla w \cdot n v
+        - \frac{\sigma}{\text{meas}(e)^\beta} \int_e v w
 
-            - \int_e J_v^{-\top} \nabla v \cdot n w
-            - \int_e J_w^{-\top} \nabla w \cdot n v
-            + 2 \frac{\sigma}{\meas(e)^\beta} \int_e v w
+    'dirichlet': interface is at Dirichlet boundary (both, v and w restricted to E1)
+
+    .. math::
+
+        - \int_e J_v^{-\top} \nabla v \cdot n w
+        - \int_e J_w^{-\top} \nabla w \cdot n v
+        + 2 \frac{\sigma}{\text{meas}(e)^\beta} \int_e v w
 
     Implementation:
 
@@ -336,12 +348,13 @@ def diffusion_sip_facet(
     The integrals are evaluated over e in reference space of e (S_e).
 
     v, w, \nabla v and \nabla w are evaluated on the _element_ reference space. This involves trafos
-    T_e^{E1}: S_e \rightarrow S_E1, and T_e^{E2}: S_e \rightarrow S_E2. So v above is actually v( T_e^{E1}( x ) ).
+    :math:`T_e^{E1}: S_e \rightarrow S_{E1}`, and :math:`T_e^{E2}: S_e \rightarrow S_{E2}`. So v above is actually
+    :math:`v( T_e^{E1}( x ) )`.
 
     The Jacobians are formed from the transformation from S_E1 to S and S_E2 to S.
 
     Parameters: the first set of affine vertices belongs to the interior element, the second set to the neighboring
-    element. An additional set of affine vertices corresponds to the interface vertices. Overall, e.g. in 2D we have
+    element. An additional set of affine vertices corresponds to the interface vertices. Overall, e.g. in 2D we have::
 
         p_0, p_1, p_2 for the inner triangular element
         p_3, p_4, p_5 for the outer triangular element
@@ -350,7 +363,7 @@ def diffusion_sip_facet(
         p_9           outer triangle vertex that is not on the interface
         p_10          for the normal from E1 to E2
 
-    In 3D we get
+    In 3D we get::
 
         p_0, p_1, p_2, p_3 for the inner tet element
         p_4, p_5, p_6, p_7 for the outer tet element
@@ -512,7 +525,9 @@ def diffusion_sip_rhs_dirichlet(
 
     This evaluates
 
-        - \int_e ( \nabla v \cdot n + \frac{\sigma}{\meas(e)^\beta} v ) g_D
+    .. math::
+
+        - \int_e \left( \nabla v \cdot n + \frac{\sigma}{\text{meas}(e)^\beta} v \right) g_D
 
     where g_D is an external function that evaluates to the Dirichlet boundary condition.
 
@@ -532,14 +547,14 @@ def diffusion_sip_rhs_dirichlet(
     The Jacobians are formed from the transformation from S_E1 to S and S_E2 to S.
 
     Parameters: the first set of affine vertices belongs to the interior element, the second set to the neighboring
-    element. An additional set of affine vertices corresponds to the interface vertices. Overall, e.g. in 2D we have
+    element. An additional set of affine vertices corresponds to the interface vertices. Overall, e.g. in 2D we have::
 
         p_0, p_1, p_2 for the inner triangular element
         p_6, p_7      for the interface
         p_8           inner triangle vertex that is not on the interface
         p_10          for the normal from E1 to E2
 
-    In 3D we get
+    In 3D we get::
 
         p_0, p_1, p_2, p_3 for the inner tet element
         p_8, p_9, p_10     for the interface
diff --git a/hog/function_space.py b/hog/function_space.py
index af4088369061277242670f6742b05fd230a9cdc1..1b576c275b454e12d2d561b9f70c04ddb0f660a4 100644
--- a/hog/function_space.py
+++ b/hog/function_space.py
@@ -80,14 +80,14 @@ class FunctionSpace(ABC):
     ) -> List[sp.MatrixBase]:
         """Returns a list containing the gradients of the shape functions on the element.
 
-        Particularly, for each (vector- or scalar-valued) shape function N = (N_1, ..., N_n) returns the _transposed_
-        Jacobian
+        Particularly, for each (vector- or scalar-valued) shape function N = (N_1, ..., N_n) returns the *transposed*
+        Jacobian::
 
-                          ⎡∂N₁/∂x₁  ···  ∂Nₙ/∂x₁⎤
-                          ⎢   ·             ·   ⎥
-        grad(N) = J(N)ᵀ = ⎢   ·             ·   ⎥
-                          ⎢   ·             ·   ⎥
-                          ⎣∂N₁/∂xₖ  ···  ∂Nₙ/∂xₖ ⎦
+                            ⎡∂N₁/∂x₁  ···  ∂Nₙ/∂x₁⎤
+                            ⎢   ·             ·   ⎥
+          grad(N) = J(N)ᵀ = ⎢   ·             ·   ⎥
+                            ⎢   ·             ·   ⎥
+                            ⎣∂N₁/∂xₖ  ···  ∂Nₙ/∂xₖ⎦
 
         i.e., the returned gradient is a column vector for scalar shape functions.
 
@@ -364,7 +364,7 @@ class TensorialVectorFunctionSpace(FunctionSpace):
     shape functions of the given scalar function space.
 
     For instance, a tensorial function space of dim d = 3 that is based on a P1 (Lagrangian) scalar function space with
-    n shape functions N_1, ..., N_n has the shape functions
+    n shape functions N_1, ..., N_n has the shape functions::
 
       (N_1, 0, 0),
       ...
@@ -377,7 +377,7 @@ class TensorialVectorFunctionSpace(FunctionSpace):
       (0, 0, N_n).
 
     This class also enables specifying only a single component to get the shape functions, e.g. only for component 1
-    (starting to count components at 0)
+    (starting to count components at 0)::
 
       (0, N_1, 0),
       ...