LaPy test cases (Remaining) (#44)

* test_visualization_meshes
* test_TriaMesh_Geodesics

Author: engrosamaali91

data/cubeTetra.ev

@@ -115,15 +115,8 @@
             ],
             "tolerance": 1e-4
         },
-        "test_compute_distance_tet": {"exp_compute_distance": 0.0},
-        "test_triangle_mesh_visualization_using_checksum": {
-            "real_img_checksum": "abaaa757a1e5c1f42a8302fbfe068f49"
-        },
-        "test_TriaMesh_Geodesics": {
-            "real_img_checksum": "8166f546e6630f017d3ec1f655c0ac85"
-        },
-        "test_Laplace_Geodesics": {
-            "real_img_checksum": "25f9b10a5d417de3854cf8b8faf77001"
+        "test_compute_distance_tet": {
+            "exp_compute_distance": 0.0
         },
         "test_Geodesics_format": {
             "expected_matrix_format": "csc",
@@ -150,6 +143,16 @@
                 6.000013,
                 6.000008
             ]
+        },
+        "test_visualization_triangle_mesh": {
+            "expected_elements": 4800,
+            "expected_dof": 2402,
+            "expected_ev": [-4.1549094e-05, 4.169634, 4.170457]
+        },
+        "test_visualization_tetrahedral_mesh": {
+            "expected_elements": 48000,
+            "expected_dof": 9261,
+            "expected_ev": [8.4652565e-05, 9.889787, 9.889887]
         }
     }
 }

data/cubeTria.ev

@@ -72,6 +72,16 @@ def test_evals_evec_dimension(load_tet_mesh, loaded_data):
 
 
 def test_gradients_normalization_and_divergence(load_tet_mesh, loaded_data):
+    """
+    Test the computation of gradients, normalization, and divergence for a TetMesh.
+
+    Parameters:
+        load_tet_mesh (fixture): Fixture to load a TetMesh for testing.
+        loaded_data (dict): Dictionary containing loaded test data.
+
+    Raises:
+        AssertionError: If any test condition is not met.
+    """
     T = load_tet_mesh
     tria = T.boundary_tria()
     bvert = np.unique(tria.t)
@@ -157,7 +167,6 @@ def test_tetMesh_Geodesics_format(load_tet_mesh, loaded_data):
     b0 = -divx
 
     # solve H x = b0
-    # print("Matrix Format now: "+H.getformat())
     if useCholmod:
         print("Solver: cholesky decomp - performance optimal ...")
         chol = cholesky(H)
@@ -174,7 +183,6 @@ def test_tetMesh_Geodesics_format(load_tet_mesh, loaded_data):
     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

lapy/utils/tests/expected_outcomes.json

@@ -0,0 +1,192 @@
+import json
+
+import numpy as np
+import pytest
+from scipy.sparse.linalg import splu
+
+from ...diffgeo import compute_divergence, compute_geodesic_f, compute_gradient
+from ...heat import diffusion
+from ...solver import Solver
+from ...tria_mesh import TriaMesh
+
+
+@pytest.fixture
+def loaded_data():
+    """
+    Load and provide the expected outcomes data from a JSON file.
+
+    Returns:
+        dict: Dictionary containing the expected outcomes data.
+    """
+    with open("lapy/utils/tests/expected_outcomes.json", "r") as f:
+        expected_outcomes = json.load(f)
+    return expected_outcomes
+
+
+# Fixture to load the TetMesh
+@pytest.fixture
+def load_square_mesh():
+    T = TriaMesh.read_off("data/square-mesh.off")
+    return T
+
+
+def test_tria_qualities(load_square_mesh):
+    """
+    Test triangle mesh quality computation.
+    """
+    T = load_square_mesh
+    computed_q = T.tria_qualities()
+    expected_q_length = 768
+    assert len(computed_q) == expected_q_length
+
+
+# Laplace
+def test_Laplace_Geodesics(load_square_mesh):
+    """
+    Test Laplace solver for geodesics on a mesh.
+    """
+
+    T = load_square_mesh
+
+    # compute first eigenfunction
+    fem = Solver(T, lump=True)
+    eval, evec = fem.eigs()
+    # vfunc = evec[:, 1]
+
+    # Get A,B (lumped), and inverse of B (as it is diagonal due to lumping)
+    A, B = fem.stiffness, fem.mass
+    Bi = B.copy()
+    Bi.data **= -1
+
+    assert B.sum() == 1.0
+    assert Bi is not B
+    # Convert A to a dense NumPy array
+    A_dense = A.toarray()
+
+    # Assert that A is symmetric
+    assert (A_dense == A_dense.T).all()
+
+    expected_eval_length = 10
+    assert len(eval) == expected_eval_length
+
+
+# Geodesics
+def test_Laplace_Geodesics_with_Gradient_Divergence(load_square_mesh):
+    """
+    Test Laplace geodesics using gradient and divergence.
+    """
+    T = load_square_mesh
+
+    # Load eigenfunction
+    fem = Solver(T, lump=True)
+    eval, evec = fem.eigs()
+    vfunc = evec[:, 1]
+
+    # Compute Laplacian using -div(grad(f))
+    grad = compute_gradient(T, vfunc)
+    divx = -compute_divergence(T, grad)
+
+    # Get the lumped mass matrix B
+    fem = Solver(T, lump=True)
+    B = fem.mass
+    Bi = B.copy()
+    Bi.data **= -1
+
+    # Apply Laplacian operator and then the inverse of B
+    Laplacian_result = -divx  # The Laplacian result
+
+    # Apply the inverse of B to recover vfunc
+    recovered_vfunc = Bi.dot(Laplacian_result)
+
+    # Check if the original vfunc and the recovered vfunc length are equal
+    assert len(recovered_vfunc) == len(vfunc)
+
+    expected_eval_length = 10
+    assert len(eval) == expected_eval_length
+
+
+def test_heat_diffusion_shape(load_square_mesh):
+    """
+    Test the shape of the heat diffusion result on a square mesh.
+
+    Parameters:
+        load_square_mesh: Fixture providing a loaded square mesh.
+
+    This test function computes the heat diffusion and verifies that the shape
+    of the result matches the expected shape.
+
+    Returns:
+        None
+    """
+    T = load_square_mesh
+    bvert = T.boundary_loops()
+    u = diffusion(T, bvert, m=1)
+    expected_shape = (len(T.v),)
+    assert u.shape == expected_shape
+
+
+def test_Geodesics_format(loaded_data, load_square_mesh):
+    """
+    Test geodesics format and accuracy.
+    """
+    T = load_square_mesh
+    bvert = T.boundary_loops()
+    u = diffusion(T, bvert, m=1)
+    # compute gradient of heat diffusion
+    tfunc = compute_gradient(T, u)
+
+    # normalize gradient
+    X = -tfunc / np.sqrt((tfunc**2).sum(1))[:, np.newaxis]
+    X = np.nan_to_num(X)
+    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
+    H = -A
+    b0 = divx
+
+    # solve H x = b0
+    # we don't need the B matrix here, as divx is the integrated divergence
+    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)
+
+    # remove shift
+    x = x - min(x)
+
+    Bi = B.copy()
+    vf = fem.poisson(-Bi * divx)
+    vf = vf - min(vf)
+    gf = compute_geodesic_f(T, u)
+    expected_matrix_format = loaded_data["expected_outcomes"]["test_Geodesics_format"][
+        "expected_matrix_format"
+    ]
+    assert H.getformat() == expected_matrix_format
+    assert not useCholmod, "Solver: cholesky decomp - performance optimal ..."
+    expected_max_x = loaded_data["expected_outcomes"]["test_Geodesics_format"][
+        "max_distance"
+    ]
+    expected_sqrt_2_div_2 = loaded_data["expected_outcomes"]["test_Geodesics_format"][
+        "expected_sqrt_2_div_2"
+    ]
+    assert np.isclose(max(x), expected_max_x)
+    computed_sqrt_2_div_2 = np.sqrt(2) / 2
+    assert np.isclose(computed_sqrt_2_div_2, expected_sqrt_2_div_2)
+    expected_max_abs_diff = loaded_data["expected_outcomes"]["test_Geodesics_format"][
+        "expected_max_abs_diff"
+    ]
+    computed_max_abs_diff = max(abs(gf - x))
+    assert np.allclose(computed_max_abs_diff, expected_max_abs_diff)

lapy/utils/tests/test_TetMesh_Geodesics.py

@@ -89,21 +89,19 @@ def test_is_oriented(tet_mesh_fixture):
     ), f"Expected is_oriented result {expected_result}, but got {result}"
 
 
-def test_boundary_trai(tet_mesh_fixture):
+def test_boundary_tria(tet_mesh_fixture):
     """
     Test computation of boundary triangles from tet mesh.
 
-    - `BT.t` represents the array of boundary triangles.
-    - `.shape[0]` counts the number of boundary triangles.
+    `BT.t` represents the array of boundary triangles.
+    `.shape[0]` counts the number of boundary triangles.
     """
     mesh = tet_mesh_fixture
     boundary_tria_mesh = mesh.boundary_tria()
 
     expected_num_traingles = 12
     assert boundary_tria_mesh.t.shape[0] == expected_num_traingles
 
-    # print(f"Found {boundary_tria_mesh.t.shape[0]} triangles on boundary.")
-
     # Check if the boundary triangle mesh is not oriented (this should fail)
     result = boundary_tria_mesh.is_oriented()
     expected_result = False

lapy/utils/tests/test_TriaMesh_Geodesics.py

@@ -297,7 +297,7 @@ def test_volume_(tria_mesh_fixture):
 
 def test_vertex_normals(tria_mesh_fixture, loaded_data):
     """
-    Testing the computation of vertex normals for oriented mesh
+    Testing the computation of vertex normals for oriented mesh.
     """
 
     # Calling tria_mesh_fixture.orient_() will modify the tria_mesh_fixture in-place and
@@ -348,7 +348,7 @@ def test_boundary_mesh(tria_mesh_fixture):
 
 def test_refine_and_boundary_loops(tria_mesh_fixture, loaded_data):
     """
-    Testing boundary loops after refining the mesh
+    Testing boundary loops after refining the mesh.
     """
     # Create a boundary mesh by dropping triangles
     tria_mesh_fixture.orient_()