Merge pull request #93 from nschloe/better-remove-cells

Better remove cells

Makefile

@@ -11,7 +11,7 @@ tag:
 upload: clean
 	@if [ "$(shell git rev-parse --abbrev-ref HEAD)" != "master" ]; then exit 1; fi
 	# https://stackoverflow.com/a/58756491/353337
-	python3 -m pep517.build --source --binary .
+	python3 -m build --sdist --wheel .
 	twine upload dist/*
 
 publish: tag upload

logo/logo.py

@@ -19,7 +19,7 @@ def _main():
     # circumcircle
     circle1 = plt.Circle(
         mesh.cell_circumcenters[0],
-        mesh.circumradius[0],
+        mesh.cell_circumradius[0],
         color=col,
         fill=False,
         linewidth=lw,
@@ -32,7 +32,7 @@ def _main():
         v1 = m1 - mesh.cell_circumcenters[0]
         e1 = (
             mesh.cell_circumcenters[0]
-            + v1 / numpy.linalg.norm(v1) * mesh.circumradius[0]
+            + v1 / numpy.linalg.norm(v1) * mesh.cell_circumradius[0]
         )
         plt.plot(
             [mesh.cell_circumcenters[0, 0], e1[0]],

meshplex/base.py

@@ -5,7 +5,7 @@
 
 
 class _BaseMesh:
-    def __init__(self, points, cells_points):
+    def __init__(self):
         self._edge_lengths = None
 
     def write(self, filename, point_data=None, cell_data=None, field_data=None):
@@ -87,7 +87,7 @@ def get_face_mask(self, subdomain):
 
         if subdomain.is_boundary_only:
             # Filter for boundary
-            is_inside = is_inside & self.is_boundary_facet
+            is_inside = is_inside & self.is_boundary_facet_local
 
         return is_inside
 
@@ -121,4 +121,3 @@ def _mark_vertices(self, subdomain):
                 is_inside = is_inside & self.is_boundary_point
 
         self.subdomains[subdomain] = {"vertices": is_inside}
-        return

meshplex/helpers.py

@@ -18,9 +18,9 @@ def get_signed_simplex_volumes(cells, pts):
 
 
 def grp_start_len(a):
-    """Given a sorted 1D input array `a`, e.g., [0 0, 1, 2, 3, 4, 4, 4], this
-    routine returns the indices where the blocks of equal integers start and
-    how long the blocks are.
+    """Given a sorted 1D input array `a`, e.g., [0 0, 1, 2, 3, 4, 4, 4], this routine
+    returns the indices where the blocks of equal integers start and how long the blocks
+    are.
     """
     # https://stackoverflow.com/a/50394587/353337
     m = numpy.concatenate([[True], a[:-1] != a[1:], [True]])
@@ -46,7 +46,7 @@ def compute_tri_areas(ei_dot_ej):
         + ei_dot_ej[1] * ei_dot_ej[2]
     )
     # vol2 is the squared volume, but can be slightly negative if it comes to round-off
-    # errors. Corrrect those.
+    # errors. Correct those.
     assert numpy.all(vol2 > -1.0e-14)
     vol2[vol2 < 0] = 0.0
     return numpy.sqrt(vol2)

meshplex/mesh_tetra.py

@@ -12,18 +12,18 @@ class MeshTetra(_BaseMesh):
 
     def __init__(self, points, cells):
         """Initialization."""
-        # Sort cells and points, first every row, then the rows themselves. This helps in
-        # many downstream applications, e.g., when constructing linear systems with the
-        # cells/edges. (When converting to CSR format, the I/J entries must be sorted.)
-        # Don't use cells.sort(axis=1) to avoid
+        # Sort cells and points, first every row, then the rows themselves. This helps
+        # in many downstream applications, e.g., when constructing linear systems with
+        # the cells/edges. (When converting to CSR format, the I/J entries must be
+        # sorted.) Don't use cells.sort(axis=1) to avoid
         # ```
         # ValueError: sort array is read-only
         # ```
         cells = numpy.sort(cells, axis=1)
         cells = cells[cells[:, 0].argsort()]
 
         self._points = points
-        super().__init__(points, cells)
+        super().__init__()
 
         # Assert that all vertices are used.
         # If there are vertices which do not appear in the cells list, this
@@ -84,19 +84,19 @@ def __init__(self, points, cells):
         ]
 
         # create ei_dot_ei, ei_dot_ej
-        self.edge_coords = (
+        self.half_edge_coords = (
             self.points[self.idx_hierarchy[1]] - self.points[self.idx_hierarchy[0]]
         )
         self.ei_dot_ei = numpy.einsum(
-            "ijkl, ijkl->ijk", self.edge_coords, self.edge_coords
+            "ijkl, ijkl->ijk", self.half_edge_coords, self.half_edge_coords
         )
         self.ei_dot_ej = numpy.einsum(
             "ijkl, ijkl->ijk",
-            self.edge_coords[[1, 2, 0]],
-            self.edge_coords[[2, 0, 1]]
+            self.half_edge_coords[[1, 2, 0]],
+            self.half_edge_coords[[2, 0, 1]]
             # This is equivalent:
-            # numpy.roll(self.edge_coords, 1, axis=0),
-            # numpy.roll(self.edge_coords, 2, axis=0),
+            # numpy.roll(self.half_edge_coords, 1, axis=0),
+            # numpy.roll(self.half_edge_coords, 2, axis=0),
         )
 
         self.ce_ratios = self._compute_ce_ratios_geometric()
@@ -105,8 +105,8 @@ def __init__(self, points, cells):
         self.is_boundary_point = None
         self._inv_faces = None
         self.edges = None
-        self.is_boundary_facet_individual = None
         self.is_boundary_facet = None
+        self.is_boundary_facet_local = None
         self.faces = None
 
         self._cell_centroids = None
@@ -137,10 +137,7 @@ def mark_boundary(self):
             self.create_cell_face_relationships()
 
         self.is_boundary_point = numpy.zeros(len(self.points), dtype=bool)
-        self.is_boundary_point[
-            self.faces["points"][self.is_boundary_facet_individual]
-        ] = True
-        return
+        self.is_boundary_point[self.faces["points"][self.is_boundary_facet]] = True
 
     def create_cell_face_relationships(self):
         # Reshape into individual faces, and take the first point per edge. (The face is
@@ -162,8 +159,8 @@ def create_cell_face_relationships(self):
         # listed twice.
         assert all(cts < 3)
 
-        self.is_boundary_facet = (cts[inv] == 1).reshape(s[2:])
-        self.is_boundary_facet_individual = cts == 1
+        self.is_boundary_facet_local = (cts[inv] == 1).reshape(s[2:])
+        self.is_boundary_facet = cts == 1
 
         self.faces = {"points": a[idx]}
 
@@ -178,8 +175,6 @@ def create_cell_face_relationships(self):
         # save for create_edge_cells
         self._inv_faces = inv
 
-        return
-
     def create_face_edge_relationships(self):
         a = numpy.vstack(
             [
@@ -203,8 +198,6 @@ def create_face_edge_relationships(self):
         face_edges = inv.reshape([3, num_faces]).T
         self.faces["edges"] = face_edges
 
-        return
-
     def _compute_cell_circumcenters(self):
         """Computes the center of the circumsphere of each cell."""
         # Just like for triangular cells, tetrahedron circumcenters are most easily
@@ -235,7 +228,6 @@ def _compute_cell_circumcenters(self):
         self._circumcenters = numpy.sum(
             alpha[None].T * self.points[self.cells["points"]], axis=1
         )
-        return
 
     # Question:
     # We're looking for an explicit expression for the algebraic c/e ratios. Might it be