Merge pull request #96 from nschloe/more-tests

More tests

meshplex/mesh_tri.py

@@ -121,7 +121,7 @@ def remove_cells(self, remove_array):
             if self._edges_cells is None:
                 self._compute_edges_cells()
 
-            # Set edge to is_boundary_edge_local=True if it was adjacent to a removed
+            # Set edge to is_boundary_edge_local=True if it is adjacent to a removed
             # cell.
             edge_ids = self.cells["edges"][~keep].flatten()
             # only consider interior edges
@@ -221,9 +221,6 @@ def remove_cells(self, remove_array):
             self._boundary_edges = None
             self._interior_edges = None
 
-        if self._is_boundary_point is not None:
-            self._is_boundary_point[self.cells["points"][~keep].flatten()] = True
-
         self.cells["points"] = self.cells["points"][keep]
         self.idx_hierarchy = self.idx_hierarchy[..., keep]
 
@@ -249,7 +246,7 @@ def remove_cells(self, remove_array):
             self._cell_circumcenters = self._cell_circumcenters[keep]
 
         if self._cell_partitions is not None:
-            self._cell_partitions = self._cell_partitions[keep]
+            self._cell_partitions = self._cell_partitions[:, keep]
 
         if self._signed_cell_areas is not None:
             self._signed_cell_areas = self._signed_cell_areas[keep]
@@ -261,6 +258,8 @@ def remove_cells(self, remove_array):
         self._cv_centroids = None
         self._cvc_cell_mask = None
         self._is_point_used = None
+        self._is_interior_point = None
+        self._is_boundary_point = None
 
         return numpy.sum(~keep)
 
@@ -501,8 +500,8 @@ def _compute_edges_cells(self):
         assert numpy.all((count == 1) == self.is_boundary_edge)
         assert numpy.all((count == 2) == self.is_interior_edge)
 
-        idx_start_count_1 = idx_start[count == 1]
-        idx_start_count_2 = idx_start[count == 2]
+        idx_start_count_1 = idx_start[self.is_boundary_edge]
+        idx_start_count_2 = idx_start[self.is_interior_edge]
         res1 = idx_sort[idx_start_count_1]
         res2 = idx_sort[numpy.array([idx_start_count_2, idx_start_count_2 + 1])]
 
@@ -511,7 +510,8 @@ def _compute_edges_cells(self):
 
         # It'd be nicer if we could organize the data differently, e.g., as a structured
         # array or as a dict. Those possibilities are slower, unfortunately, for some
-        # operations. <https://github.com/numpy/numpy/issues/17850>
+        # operations in remove_cells() (and perhaps elsewhere).
+        # <https://github.com/numpy/numpy/issues/17850>
         self._edges_cells = {
             # rows:
             #  0: edge id
@@ -532,14 +532,16 @@ def _compute_edges_cells(self):
     @property
     def edges_cells_idx(self):
         if self._edges_cells_idx is None:
+            if self._edges_cells is None:
+                self._compute_edges_cells()
             assert self.is_boundary_edge is not None
             # For each edge, store the index into the respective edge array.
             num_edges = len(self.edges["points"])
             self._edges_cells_idx = numpy.empty(num_edges, dtype=int)
-            num_b_edges = numpy.sum(self.is_boundary_edge)
-            num_i_edges = numpy.sum(self.is_interior_edge)
-            self._edges_cells_idx[self.is_boundary_edge] = numpy.arange(num_b_edges)
-            self._edges_cells_idx[self.is_interior_edge] = numpy.arange(num_i_edges)
+            num_b = numpy.sum(self.is_boundary_edge)
+            num_i = numpy.sum(self.is_interior_edge)
+            self._edges_cells_idx[self.edges_cells["boundary"][0]] = numpy.arange(num_b)
+            self._edges_cells_idx[self.edges_cells["interior"][0]] = numpy.arange(num_i)
         return self._edges_cells_idx
 
     @property
@@ -548,7 +550,7 @@ def cell_partitions(self):
             # Compute the control volumes. Note that
             #
             #   0.5 * (0.5 * edge_length) * covolume
-            # = 0.25 * edge_length**2 * ce_ratio_edge_ratio
+            # = 0.25 * edge_length ** 2 * ce_ratio_edge_ratio
             #
             self._cell_partitions = self.ei_dot_ei * self.ce_ratios / 4
         return self._cell_partitions

test/mesh_tri/test_remove_cells.py

@@ -8,16 +8,87 @@
 this_dir = pathlib.Path(__file__).resolve().parent
 
 
+def assert_mesh_consistency(mesh):
+    assert numpy.all(numpy.logical_xor(mesh.is_boundary_edge, mesh.is_interior_edge))
+
+    bpts = numpy.array(
+        [
+            mesh.is_boundary_point,
+            mesh.is_interior_point,
+            ~mesh.is_point_used,
+        ]
+    )
+    assert numpy.all(numpy.sum(bpts, axis=0) == 1)
+
+    # consistency check for edges_cells
+    assert numpy.all(mesh.is_boundary_edge[mesh.edges_cells["boundary"][0]])
+    assert not numpy.any(mesh.is_boundary_edge[mesh.edges_cells["interior"][0]])
+
+    for edge_id, cell_id, local_edge_id in mesh.edges_cells["boundary"].T:
+        assert edge_id == mesh.cells["edges"][cell_id][local_edge_id]
+
+    for edge_id, cell_id0, cell_id1, local_edge_id0, local_edge_id1 in mesh.edges_cells[
+        "interior"
+    ].T:
+        assert edge_id == mesh.cells["edges"][cell_id0][local_edge_id0]
+        assert edge_id == mesh.cells["edges"][cell_id1][local_edge_id1]
+
+    # check consistency of edges_cells_idx with edges_cells
+    for edge_id, idx in enumerate(mesh.edges_cells_idx):
+        if mesh.is_boundary_edge[edge_id]:
+            assert edge_id == mesh.edges_cells["boundary"][0, idx]
+        else:
+            assert edge_id == mesh.edges_cells["interior"][0, idx]
+
+    # TODO add more consistency checks
+
+
 def compute_all_entities(mesh):
     mesh.is_boundary_point
+    mesh.is_interior_point
     mesh.is_boundary_edge_local
     mesh.is_boundary_edge
     mesh.is_boundary_cell
     mesh.cell_volumes
     mesh.ce_ratios
     mesh.signed_cell_areas
     mesh.cell_centroids
+    mesh.control_volumes
     mesh.create_edges()
+    mesh.edges_cells
+    mesh.edges_cells_idx
+    mesh.boundary_edges
+    mesh.interior_edges
+    mesh.cell_circumcenters
+    mesh.ce_ratios_per_interior_edge
+    mesh.control_volume_centroids
+
+    assert mesh._cv_cell_mask is not None
+    assert mesh.edges is not None
+    assert mesh.subdomains is not {}
+    assert mesh._is_interior_point is not None
+    assert mesh._is_boundary_point is not None
+    assert mesh._is_boundary_edge_local is not None
+    assert mesh._is_boundary_edge is not None
+    assert mesh._is_boundary_cell is not None
+    assert mesh._edges_cells is not None
+    assert mesh._edges_cells_idx is not None
+    assert mesh._boundary_edges is not None
+    assert mesh._interior_edges is not None
+    assert mesh._is_point_used is not None
+    assert mesh._half_edge_coords is not None
+    assert mesh._ei_dot_ei is not None
+    assert mesh._ei_dot_ej is not None
+    assert mesh._cell_volumes is not None
+    assert mesh._ce_ratios is not None
+    assert mesh._cell_circumcenters is not None
+    assert mesh._interior_ce_ratios is not None
+    assert mesh._control_volumes is not None
+    assert mesh._cell_partitions is not None
+    assert mesh._cv_centroids is not None
+    assert mesh._cvc_cell_mask is not None
+    assert mesh._signed_cell_areas is not None
+    assert mesh._cell_centroids is not None
 
 
 def assert_mesh_equality(mesh0, mesh1):
@@ -27,17 +98,62 @@ def assert_mesh_equality(mesh0, mesh1):
     assert numpy.all(numpy.abs(mesh0.points - mesh1.points) < 1.0e-14)
     assert numpy.all(mesh0.edges["points"] == mesh1.edges["points"])
 
+    # # Assume that in mesh1, the rows are ordered such that the edge indices [0] are in
+    # # order. The mesh0 array isn't order in many cases, e.g., if cells were removed and
+    # # rows appened the boundary array. Hence, justsort the array in mesh0 and compare.
+    # k = numpy.argsort(mesh0.edges_cells["boundary"][0])
+    # assert numpy.all(mesh0.edges_cells["boundary"][:, k] == mesh1.edges_cells["boundary"])
+
+    # # The interior edges_cells are already ordered, even after remove_cells(). (As
+    # # opposed to boundary edges, there can be no new interior edges, just some are
+    # # removed.)
+    # print(mesh0.edges_cells["interior"].T)
+    # print()
+    # print(mesh1.edges_cells["interior"].T)
+    # assert numpy.all(mesh0.edges_cells["interior"] == mesh1.edges_cells["interior"])
+    # exit(1)
+    # # print()
+    # # print(mesh0.edges_cells["interior"])
+    # # print()
+    # # print(mesh1.edges_cells["interior"])
+    #
+    # # These should not be equal; see reordering above
+    # assert numpy.all(mesh0.edges_cells_idx == mesh1.edges_cells_idx)
+
+    assert numpy.all(mesh0.boundary_edges == mesh1.boundary_edges)
+    assert numpy.all(mesh0.interior_edges == mesh1.interior_edges)
+
+    assert numpy.all(mesh0.is_point_used == mesh1.is_point_used)
     assert numpy.all(mesh0.is_boundary_point == mesh1.is_boundary_point)
+    assert numpy.all(mesh0.is_interior_point == mesh1.is_interior_point)
     assert numpy.all(mesh0.is_boundary_edge_local == mesh1.is_boundary_edge_local)
     assert numpy.all(mesh0.is_boundary_edge == mesh1.is_boundary_edge)
     assert numpy.all(mesh0.is_boundary_cell == mesh1.is_boundary_cell)
 
+    assert numpy.all(numpy.abs(mesh0.ei_dot_ei - mesh1.ei_dot_ei) < 1.0e-14)
+    assert numpy.all(numpy.abs(mesh0.ei_dot_ej - mesh1.ei_dot_ej) < 1.0e-14)
     assert numpy.all(numpy.abs(mesh0.cell_volumes - mesh1.cell_volumes) < 1.0e-14)
     assert numpy.all(numpy.abs(mesh0.ce_ratios - mesh1.ce_ratios) < 1.0e-14)
+    assert numpy.all(
+        numpy.abs(mesh0.ce_ratios_per_interior_edge - mesh1.ce_ratios_per_interior_edge)
+        < 1.0e-14
+    )
     assert numpy.all(
         numpy.abs(mesh0.signed_cell_areas - mesh1.signed_cell_areas) < 1.0e-14
     )
     assert numpy.all(numpy.abs(mesh0.cell_centroids - mesh1.cell_centroids) < 1.0e-14)
+    assert numpy.all(
+        numpy.abs(mesh0.cell_circumcenters - mesh1.cell_circumcenters) < 1.0e-14
+    )
+    assert numpy.all(numpy.abs(mesh0.control_volumes - mesh1.control_volumes) < 1.0e-14)
+
+    ipu = mesh0.is_point_used
+    assert numpy.all(
+        numpy.abs(
+            mesh0.control_volume_centroids[ipu] - mesh1.control_volume_centroids[ipu]
+        )
+        < 1.0e-14
+    )
 
 
 def get_mesh0():
@@ -116,6 +232,7 @@ def test_remove_cells(remove_idx, expected_num_cells, expected_num_edges):
     mesh.remove_cells(remove_idx)
     assert len(mesh.cells["points"]) == expected_num_cells
     assert len(mesh.edges["points"]) == expected_num_edges
+    assert_mesh_consistency(mesh)
 
 
 def test_remove_cells_boundary():
@@ -146,6 +263,8 @@ def test_remove_cells_boundary():
     # now lets remove some cells
     mesh.remove_cells([0])
 
+    assert_mesh_consistency(mesh)
+
     assert numpy.all(mesh.is_boundary_point)
     assert numpy.all(mesh.is_boundary_edge_local[0] == [False, False, False])
     assert numpy.all(mesh.is_boundary_edge_local[1] == [True, False, True])
@@ -163,14 +282,16 @@ def test_remove_all():
     assert numpy.all(mesh.is_point_used)
 
     mesh.remove_cells([0])
+    assert_mesh_consistency(mesh)
     assert not numpy.any(mesh.is_point_used)
 
 
 @pytest.mark.parametrize(
     "mesh0, remove_cells",
     [
-        # (get_mesh0(), [0]),
-        # (get_mesh1(), [0, 1]),
+        # TODO enable all
+        (get_mesh0(), [0]),
+        (get_mesh1(), [0, 1]),
         (get_mesh2(), [0, 3, 57, 59, 60, 61, 100]),
     ],
 )
@@ -180,6 +301,9 @@ def test_reference(mesh0, remove_cells):
     compute_all_entities(mesh0)
     # now remove some cells
     mesh0.remove_cells(remove_cells)
+
+    assert_mesh_consistency(mesh0)
+
     # recreate the reduced mesh from scratch
     mesh1 = meshplex.MeshTri(mesh0.points, mesh0.cells["points"])
     mesh1.create_edges()