diff --git a/lapy/utils/tests/expected_outcomes.json b/lapy/utils/tests/expected_outcomes.json
index 81de57e..03a6868 100644
--- a/lapy/utils/tests/expected_outcomes.json
+++ b/lapy/utils/tests/expected_outcomes.json
@@ -122,7 +122,7 @@
         "test_TriaMesh_Geodesics": {
             "real_img_checksum": "8166f546e6630f017d3ec1f655c0ac85"
         },
-        "test_Laplase_Geodesics": {
+        "test_Laplace_Geodesics": {
             "real_img_checksum": "25f9b10a5d417de3854cf8b8faf77001"
         },
         "test_Geodesics_format": {
diff --git a/lapy/utils/tests/test_TetMesh_Geodesics.py b/lapy/utils/tests/test_TetMesh_Geodesics.py
index 136c6dc..4fa0b9b 100644
--- a/lapy/utils/tests/test_TetMesh_Geodesics.py
+++ b/lapy/utils/tests/test_TetMesh_Geodesics.py
@@ -70,15 +70,14 @@ def test_evals_evec_dimension(load_tet_mesh, loaded_data):
 
 # Geodesics
 
-T = TetMesh.read_vtk("data/cubeTetra.vtk")
-tria = T.boundary_tria()
-bvert = np.unique(tria.t)
-u = diffusion(T, bvert, m=1)
-
 
 def test_gradients_normalization_and_divergence(load_tet_mesh, loaded_data):
-    # Compute gradients
     T = load_tet_mesh
+    tria = T.boundary_tria()
+    bvert = np.unique(tria.t)
+    u = diffusion(T, bvert, m=1)
+
+    # Compute gradients
     tfunc = compute_gradient(T, u)
 
     # Define the expected shape of tfunc (gradient)
@@ -119,56 +118,7 @@ def test_gradients_normalization_and_divergence(load_tet_mesh, loaded_data):
     assert divx.shape == expected_divx_shape
 
 
-# get gradients
-tfunc = compute_gradient(T, u)
-# flip and normalize
-X = -tfunc / np.sqrt((tfunc**2).sum(1))[:, np.newaxis]
-X = np.nan_to_num(X)
-# compute divergence
-divx = compute_divergence(T, X)
-
-# compute distance
-useCholmod = True
-try:
-    from sksparse.cholmod import cholesky
-except ImportError:
-    useCholmod = False
-
-
-fem = Solver(T, lump=True)
-A, B = fem.stiffness, fem.mass  # computed above when creating Solver
-
-H = A
-b0 = -divx
-
-# solve H x = b0
-# print("Matrix Format now: "+H.getformat())
-if useCholmod:
-    print("Solver: cholesky decomp - performance optimal ...")
-    chol = cholesky(H)
-    x = chol(b0)
-else:
-    print("Solver: spsolve (LU decomp) - performance not optimal ...")
-    lu = splu(H)
-    x = lu.solve(b0)
-
-x = x - np.min(x)
-
-
-# get heat diffusion
-gu = compute_geodesic_f(T, u)
-
-v1func = T.v[:, 0] * T.v[:, 0] + T.v[:, 1] * T.v[:, 1] + T.v[:, 2] * T.v[:, 2]
-grad = compute_gradient(T, v1func)
-glength = np.sqrt(np.sum(grad * grad, axis=1))
-# fcols=glength
-A, B = fem.stiffness, fem.mass
-Bi = B.copy()
-Bi.data **= -1
-divx2 = Bi * divx
-
-
-def test_tetMesh_Geodesics_format(loaded_data):
+def test_tetMesh_Geodesics_format(load_tet_mesh, loaded_data):
     """
     Test if matrix format, solver settings, max distance,
     and computed values match the expected outcomes.
@@ -179,6 +129,58 @@ def test_tetMesh_Geodesics_format(loaded_data):
     Raises:
     - AssertionError: If any test condition is not met.
     """
+
+    T = load_tet_mesh
+    tria = T.boundary_tria()
+    bvert = np.unique(tria.t)
+    u = diffusion(T, bvert, m=1)
+
+    # get gradients
+    tfunc = compute_gradient(T, u)
+    # flip and normalize
+    X = -tfunc / np.sqrt((tfunc**2).sum(1))[:, np.newaxis]
+    X = np.nan_to_num(X)
+    # compute divergence
+    divx = compute_divergence(T, X)
+
+    # compute distance
+    useCholmod = True
+    try:
+        from sksparse.cholmod import cholesky
+    except ImportError:
+        useCholmod = False
+
+    fem = Solver(T, lump=True)
+    A, B = fem.stiffness, fem.mass  # computed above when creating Solver
+
+    H = A
+    b0 = -divx
+
+    # solve H x = b0
+    # print("Matrix Format now: "+H.getformat())
+    if useCholmod:
+        print("Solver: cholesky decomp - performance optimal ...")
+        chol = cholesky(H)
+        x = chol(b0)
+    else:
+        print("Solver: spsolve (LU decomp) - performance not optimal ...")
+        lu = splu(H)
+        x = lu.solve(b0)
+
+    x = x - np.min(x)
+
+    # get heat diffusion
+    gu = compute_geodesic_f(T, u)
+
+    v1func = T.v[:, 0] * T.v[:, 0] + T.v[:, 1] * T.v[:, 1] + T.v[:, 2] * T.v[:, 2]
+    grad = compute_gradient(T, v1func)
+    glength = np.sqrt(np.sum(grad * grad, axis=1))
+    # fcols=glength
+    A, B = fem.stiffness, fem.mass
+    Bi = B.copy()
+    Bi.data **= -1
+    divx2 = Bi * divx
+
     expected_matrix_format = loaded_data["expected_outcomes"]["test_TetMesh_Geodesics"][
         "expected_matrix_format"
     ]
diff --git a/lapy/utils/tests/test_shape_DNA.py b/lapy/utils/tests/test_shape_DNA.py
index a70cbc6..fbf98dc 100644
--- a/lapy/utils/tests/test_shape_DNA.py
+++ b/lapy/utils/tests/test_shape_DNA.py
@@ -10,10 +10,6 @@
 tria = TriaMesh.read_vtk("data/cubeTria.vtk")
 tet = TetMesh.read_vtk("data/cubeTetra.vtk")
 
-ev = compute_shapedna(tria, k=3)
-ev["Eigenvectors"]
-ev["Eigenvalues"]
-
 
 @pytest.fixture
 def loaded_data():
@@ -36,6 +32,8 @@ def test_compute_shapedna(loaded_data):
         AssertionError: If computed eigenvalues don't match expected values within tolerance.
         AssertionError: If eigenvalues' dtype isn't float32.
     """
+    ev = compute_shapedna(tria, k=3)
+
     expected_Eigenvalues = np.array(
         loaded_data["expected_outcomes"]["test_compute_shapedna"][
             "expected_eigenvalues"
@@ -57,6 +55,8 @@ def test_normalize_ev_geometry(loaded_data):
         AssertionError: If normalized eigenvalues don't match expected values within tolerance.
         AssertionError: If normalized eigenvalues' dtype isn't float32.
     """
+    ev = compute_shapedna(tria, k=3)
+
     expected_normalized_values = np.array(
         loaded_data["expected_outcomes"]["test_normalize_ev_geometry"][
             "expected_normalized_values"
@@ -83,6 +83,8 @@ def test_reweight_ev(loaded_data):
         AssertionError: If reweighted eigenvalues don't match expected values within tolerance.
         AssertionError: If reweighted eigenvalues' dtype isn't float32.
     """
+    ev = compute_shapedna(tria, k=3)
+
     expected_reweighted_values = np.array(
         loaded_data["expected_outcomes"]["test_reweight_ev"][
             "expected_reweighted_values"
@@ -106,6 +108,8 @@ def test_compute_distance(loaded_data):
     Raises:
         AssertionError: If computed distance doesn't match the expected value.
     """
+    ev = compute_shapedna(tria, k=3)
+
     expected_compute_distance = loaded_data["expected_outcomes"][
         "test_compute_distance"
     ]["expected_compute_distance"]
@@ -114,11 +118,7 @@ def test_compute_distance(loaded_data):
     assert computed_distance == expected_compute_distance
 
 
-# Repeat testing steps for a tetrahedral mesh
-# compute eigenvalues and eigenvectors for tet mesh
-evTet = compute_shapedna(tet, k=3)
-evTet["Eigenvectors"]
-evTet["Eigenvalues"]
+# Repeating test steps for a tetrahedral mesh
 
 
 def test_compute_shapedna_tet(loaded_data):
@@ -132,6 +132,8 @@ def test_compute_shapedna_tet(loaded_data):
         AssertionError: If computed eigenvalues don't match expected values within tolerance.
         AssertionError: If eigenvalues' dtype isn't float32.
     """
+    evTet = compute_shapedna(tet, k=3)
+
     expected_eigen_values = np.array(
         loaded_data["expected_outcomes"]["test_compute_shapedna_tet"][
             "expected_eigen_values"
@@ -156,6 +158,8 @@ def test_normalize_ev_geometry_tet(loaded_data):
         AssertionError: If normalized eigenvalues don't match expected values within tolerance.
         AssertionError: If normalized eigenvalues' dtype isn't float32.
     """
+    evTet = compute_shapedna(tet, k=3)
+
     expected_normalized_values = np.array(
         loaded_data["expected_outcomes"]["test_normalize_ev_geometry_tet"][
             "expected_normalized_values"
@@ -183,6 +187,8 @@ def test_reweight_ev_tet(loaded_data):
     Raises:
         AssertionError: If reweighted eigenvalues don't match expected values within tolerance.
     """
+    evTet = compute_shapedna(tet, k=3)
+
     expected_reweighted_values = np.array(
         loaded_data["expected_outcomes"]["test_reweight_ev_tet"][
             "expected_reweighted_values"
@@ -206,6 +212,8 @@ def test_compute_distance_tet(loaded_data):
     Raises:
         AssertionError: If computed distance doesn't match the expected value.
     """
+    evTet = compute_shapedna(tet, k=3)
+
     # compute distance for tria eigenvalues (trivial case)
     computed_distance = compute_distance(evTet["Eigenvalues"], evTet["Eigenvalues"])
     expected_compute_distance = loaded_data["expected_outcomes"][
diff --git a/lapy/utils/tests/test_tet_mesh.py b/lapy/utils/tests/test_tet_mesh.py
index 6afe6f4..167193d 100644
--- a/lapy/utils/tests/test_tet_mesh.py
+++ b/lapy/utils/tests/test_tet_mesh.py
@@ -89,19 +89,6 @@ def test_is_oriented(tet_mesh_fixture):
     ), f"Expected is_oriented result {expected_result}, but got {result}"
 
 
-# def test_avg_edge_length(tet_mesh_fixture):
-#     """
-#     Testing computation of average edge length
-#     """
-#     expected_result = 1.0543647924813107
-#     mesh = tet_mesh_fixture
-#     result = mesh.avg_edge_length()
-
-#     assert (
-#         result == expected_result
-#     ), f"Expected average edge length {expected_result}, but got {result}"
-
-
 def test_boundary_trai(tet_mesh_fixture):
     """
     Test computation of boundary triangles from tet mesh.