diff --git a/src/darsia/measure/wasserstein.py b/src/darsia/measure/wasserstein.py
index 5dc39f47..c5167299 100644
--- a/src/darsia/measure/wasserstein.py
+++ b/src/darsia/measure/wasserstein.py
@@ -583,7 +583,7 @@ def transport_density(
             # Apply numerical integration of RT0 extensions into cells.
             # Underlying functional for mixed finite element method (MFEM).
             quad_pts, quad_weights = darsia.quadrature.gauss_reference_cell(
-                self.grid.dim, 3
+                self.grid.dim, "max"
             )
 
         elif self.l1_mode == "constant_subcell_projection":
@@ -1216,6 +1216,9 @@ def _solve(self, flat_mass_diff: np.ndarray) -> tuple[float, np.ndarray, dict]:
             self.darcy_init.copy(), rhs.copy(), solution_i
         )
 
+        # Initialize distance in case below iteration fails
+        new_distance = 0
+
         # Initialize container for storing the convergence history
         convergence_history = {
             "distance": [],
@@ -1493,6 +1496,9 @@ def _solve(self, flat_mass_diff: np.ndarray) -> tuple[float, np.ndarray, dict]:
             self.darcy_init.copy(), rhs.copy(), solution_i
         )
 
+        # Initialize distance in case below iteration fails
+        new_distance = 0
+
         # Initialize container for storing the convergence history
         convergence_history = {
             "distance": [],
diff --git a/src/darsia/utils/quadrature.py b/src/darsia/utils/quadrature.py
index 5dddd57b..bbc674ab 100644
--- a/src/darsia/utils/quadrature.py
+++ b/src/darsia/utils/quadrature.py
@@ -1,9 +1,11 @@
 """Quadrature rules for numerical integration."""
 
+from typing import Union
+
 import numpy as np
 
 
-def gauss(dim: int, order: int) -> tuple[np.ndarray, np.ndarray]:
+def gauss(dim: int, order: Union[int, str]) -> tuple[np.ndarray, np.ndarray]:
     """Return the Gauss points and weights for the given dimension and order.
 
     These are the Gauss points and weights for the reference element [-1, 1]^dim.
@@ -17,6 +19,8 @@ def gauss(dim: int, order: int) -> tuple[np.ndarray, np.ndarray]:
 
     """
     if dim == 1:
+        if order == "max":
+            order = 4
         if order == 0:
             return np.array([0.0]), np.array([2.0])
         elif order == 1:
@@ -73,6 +77,8 @@ def gauss(dim: int, order: int) -> tuple[np.ndarray, np.ndarray]:
                 f"Gauss points of order {order} not implemented for dimension {dim}."
             )
     elif dim == 2:
+        if order == "max":
+            order = 3
         if order == 0:
             return (
                 np.array([[0.0, 0.0]]),
@@ -215,6 +221,8 @@ def gauss(dim: int, order: int) -> tuple[np.ndarray, np.ndarray]:
                 f"Gauss points of order {order} not implemented for dimension {dim}."
             )
     elif dim == 3:
+        if order == "max":
+            order = 2
         if order == 0:
             return (
                 np.array([[0.0, 0.0, 0.0]]),