Merge pull request #122 from nschloe/integrate-line-mesh

Integrate line mesh

meshplex/_helpers.py

@@ -192,3 +192,14 @@ def compute_triangle_circumcenters(X, cell_partitions):
 
     a = X * alpha[..., None]
     return a[0] + a[1] + a[2]
+
+
+def _dot(a, n):
+    """Dot product, preserve the leading n dimensions."""
+    # einsum is faster if the tail survives, e.g., ijk,ijk->jk.
+    # <https://gist.github.com/nschloe/8bc015cc1a9e5c56374945ddd711df7b>
+    # TODO reorganize the data?
+    assert n <= len(a.shape)
+    # Would use -1 as second argument, but <https://github.com/numpy/numpy/issues/18519>
+    b = a.reshape(*a.shape[:n], np.prod(a.shape[n:]).astype(int))
+    return np.einsum("...i,...i->...", b, b)

meshplex/_mesh.py

@@ -6,6 +6,7 @@
 
 from ._exceptions import MeshplexError
 from ._helpers import (
+    _dot,
     compute_tri_areas,
     compute_triangle_circumcenters,
     grp_start_len,
@@ -17,6 +18,9 @@
 
 class Mesh:
     def __init__(self, points, cells, sort_cells=False):
+        points = np.asarray(points)
+        cells = np.asarray(cells)
+
         if sort_cells:
             # Sort cells, first every row, then the rows themselves. This helps in many
             # downstream applications, e.g., when constructing linear systems with the
@@ -28,8 +32,6 @@ def __init__(self, points, cells, sort_cells=False):
             cells = np.sort(cells, axis=1)
             cells = cells[cells[:, 0].argsort()]
 
-        points = np.asarray(points)
-        cells = np.asarray(cells)
         # assert len(points.shape) <= 2, f"Illegal point coordinates shape {points.shape}"
         assert len(cells.shape) == 2, f"Illegal cells shape {cells.shape}"
         self.n = cells.shape[1]
@@ -139,6 +141,7 @@ def _reset_point_data(self):
         self._heights = None
         self._ce_ratios = None
         self._cell_partitions = None
+        self._control_volumes = None
 
         # only used for tetra
         self._zeta = None
@@ -182,14 +185,7 @@ def half_edge_coords(self):
     @property
     def ei_dot_ei(self):
         if self._ei_dot_ei is None:
-            # einsum is faster if the tail survives, e.g., ijk,ijk->jk.
-            # <https://gist.github.com/nschloe/8bc015cc1a9e5c56374945ddd711df7b>
-            # TODO reorganize the data?
-            # TODO if the dimensionality of the points is 0 (e.g., MeshLine), then this
-            # is wrong
-            self._ei_dot_ei = np.einsum(
-                "...k, ...k->...", self.half_edge_coords, self.half_edge_coords
-            )
+            self._ei_dot_ei = _dot(self.half_edge_coords, self.n - 1)
         return self._ei_dot_ei
 
     @property
@@ -381,7 +377,6 @@ def _mark_vertices(self, subdomain):
 
             if subdomain.is_boundary_only:
                 # Filter boundary
-                self.mark_boundary()
                 is_inside = is_inside & self.is_boundary_point
 
         self.subdomains[subdomain] = {"vertices": is_inside}
@@ -568,13 +563,6 @@ def is_interior_point(self):
             self._is_interior_point = self.is_point_used & ~self.is_boundary_point
         return self._is_interior_point
 
-    def mark_boundary(self):
-        warnings.warn(
-            "mark_boundary() does nothing. "
-            "Boundary entities are computed on the fly.",
-            DeprecationWarning,
-        )
-
     @property
     def facets_cells(self):
         if self._facets_cells is None:
@@ -986,3 +974,62 @@ def remove_duplicate_cells(self):
             remove |= rem
 
         return self.remove_cells(remove)
+
+    def get_control_volumes(self, cell_mask=None):
+        """The control volumes around each vertex. Optionally disregard the
+        contributions from particular cells. This is useful, for example, for
+        temporarily disregarding flat cells on the boundary when performing Lloyd mesh
+        optimization.
+        """
+        assert self.n == 3
+        if cell_mask is None:
+            cell_mask = np.zeros(self.cell_partitions.shape[1], dtype=bool)
+
+        if self._control_volumes is None or np.any(cell_mask != self._cv_cell_mask):
+            # Summing up the arrays first makes the work on bincount a bit lighter.
+            v = self.cell_partitions[:, ~cell_mask]
+            vals = np.array([v[1] + v[2], v[2] + v[0], v[0] + v[1]])
+            # sum all the vals into self._control_volumes at ids
+            self.cells["points"][~cell_mask].T.reshape(-1)
+            self._control_volumes = np.bincount(
+                self.cells["points"][~cell_mask].T.reshape(-1),
+                weights=vals.reshape(-1),
+                minlength=len(self.points),
+            )
+            self._cv_cell_mask = cell_mask
+        return self._control_volumes
+
+    @property
+    def control_volumes(self):
+        """The control volumes around each vertex."""
+        if self._control_volumes is None:
+            if self.n == 2:
+                self._control_volumes = np.zeros(len(self.points), dtype=float)
+                for k, node_ids in enumerate(self.cells["points"]):
+                    self._control_volumes[node_ids] += 0.5 * self.cell_volumes[k]
+            elif self.n == 3:
+                return self.get_control_volumes()
+            else:
+                assert self.n == 4
+                #   1/3. * (0.5 * edge_length) * covolume
+                # = 1/6 * edge_length**2 * ce_ratio_edge_ratio
+                v = self.ei_dot_ei * self.ce_ratios / 6
+                # Explicitly sum up contributions per cell first. Makes np.add.at
+                # faster.  For every point k (range(4)), check for which edges k appears
+                # in local_idx, and sum() up the v's from there.
+                vals = np.array(
+                    [
+                        v[0, 2] + v[1, 1] + v[2, 3] + v[0, 1] + v[1, 3] + v[2, 2],
+                        v[0, 3] + v[1, 2] + v[2, 0] + v[0, 2] + v[1, 0] + v[2, 3],
+                        v[0, 0] + v[1, 3] + v[2, 1] + v[0, 3] + v[1, 1] + v[2, 0],
+                        v[0, 1] + v[1, 0] + v[2, 2] + v[0, 0] + v[1, 2] + v[2, 1],
+                    ]
+                ).T
+                #
+                self._control_volumes = np.bincount(
+                    self.cells["points"].reshape(-1),
+                    vals.reshape(-1),
+                    minlength=len(self.points),
+                )
+
+        return self._control_volumes

meshplex/_mesh_line.py

@@ -1,40 +1,11 @@
-import numpy as np
+from ._mesh import Mesh
 
 
-class MeshLine:
+class MeshLine(Mesh):
     """Class for handling line segment "meshes"."""
 
     def __init__(self, points, cells):
-        self.points = np.asarray(points)
-
-        num_cells = len(cells)
-        self.cells = np.empty(num_cells, dtype=np.dtype([("nodes", (int, 2))]))
-        self.cells["nodes"] = cells
-
-        self.create_cell_volumes()
-        self.create_control_volumes()
-
-    def create_cell_volumes(self):
-        """Computes the volumes of the "cells" in the mesh."""
-        self.cell_volumes = np.array(
-            [
-                abs(self.points[cell["nodes"]][1] - self.points[cell["nodes"]][0])
-                for cell in self.cells
-            ]
-        )
-
-    def create_control_volumes(self):
-        """Compute the control volumes of all nodes in the mesh."""
-        self.control_volumes = np.zeros(len(self.points), dtype=float)
-        for k, cell in enumerate(self.cells):
-            node_ids = cell["nodes"]
-            self.control_volumes[node_ids] += 0.5 * self.cell_volumes[k]
-
-        # Sanity checks.
-        sum_cv = sum(self.control_volumes)
-        sum_cells = sum(self.cell_volumes)
-        alpha = sum_cv - sum_cells
-        assert abs(alpha) < 1.0e-6
+        super().__init__(points, cells)
 
     def show_vertex_function(self, u):
         import matplotlib.pyplot as plt

meshplex/_mesh_tetra.py

@@ -12,7 +12,6 @@ class MeshTetra(Mesh):
     def __init__(self, points, cells, sort_cells=False):
         super().__init__(points, cells, sort_cells=sort_cells)
 
-        self._control_volumes = None
         self.subdomains = {}
 
         self.ce_ratios = self._compute_ce_ratios_geometric()
@@ -246,33 +245,6 @@ def q_vol_rms_edgelength3(self):
         alpha = np.sqrt(2) / 12  # normalization factor
         return self.cell_volumes / rms ** 3 / alpha
 
-    @property
-    def control_volumes(self):
-        """Compute the control volumes of all points in the mesh."""
-        if self._control_volumes is None:
-            #   1/3. * (0.5 * edge_length) * covolume
-            # = 1/6 * edge_length**2 * ce_ratio_edge_ratio
-            v = self.ei_dot_ei * self.ce_ratios / 6
-            # Explicitly sum up contributions per cell first. Makes np.add.at faster.
-            # For every point k (range(4)), check for which edges k appears in local_idx,
-            # and sum() up the v's from there.
-            vals = np.array(
-                [
-                    v[0, 2] + v[1, 1] + v[2, 3] + v[0, 1] + v[1, 3] + v[2, 2],
-                    v[0, 3] + v[1, 2] + v[2, 0] + v[0, 2] + v[1, 0] + v[2, 3],
-                    v[0, 0] + v[1, 3] + v[2, 1] + v[0, 3] + v[1, 1] + v[2, 0],
-                    v[0, 1] + v[1, 0] + v[2, 2] + v[0, 0] + v[1, 2] + v[2, 1],
-                ]
-            ).T
-            #
-            self._control_volumes = np.bincount(
-                self.cells["points"].reshape(-1),
-                vals.reshape(-1),
-                minlength=len(self.points),
-            )
-
-        return self._control_volumes
-
     def num_delaunay_violations(self):
         # Delaunay violations are present exactly on the interior faces where the sum of
         # the signed distances between face circumcenter and tetrahedron circumcenter is

meshplex/_mesh_tri.py

@@ -30,7 +30,6 @@ def _reset_point_data(self):
         """Reset all data that changes when point coordinates changes."""
         super()._reset_point_data()
         self._interior_ce_ratios = None
-        self._control_volumes = None
         self._cv_centroids = None
         self._cvc_cell_mask = None
 
@@ -87,34 +86,6 @@ def ce_ratios_per_interior_facet(self):
 
         return self._interior_ce_ratios
 
-    def get_control_volumes(self, cell_mask=None):
-        """The control volumes around each vertex. Optionally disregard the
-        contributions from particular cells. This is useful, for example, for
-        temporarily disregarding flat cells on the boundary when performing Lloyd mesh
-        optimization.
-        """
-        if cell_mask is None:
-            cell_mask = np.zeros(self.cell_partitions.shape[1], dtype=bool)
-
-        if self._control_volumes is None or np.any(cell_mask != self._cv_cell_mask):
-            # Summing up the arrays first makes the work on bincount a bit lighter.
-            v = self.cell_partitions[:, ~cell_mask]
-            vals = np.array([v[1] + v[2], v[2] + v[0], v[0] + v[1]])
-            # sum all the vals into self._control_volumes at ids
-            self.cells["points"][~cell_mask].T.reshape(-1)
-            self._control_volumes = np.bincount(
-                self.cells["points"][~cell_mask].T.reshape(-1),
-                weights=vals.reshape(-1),
-                minlength=len(self.points),
-            )
-            self._cv_cell_mask = cell_mask
-        return self._control_volumes
-
-    @property
-    def control_volumes(self):
-        """The control volumes around each vertex."""
-        return self.get_control_volumes()
-
     def get_control_volume_centroids(self, cell_mask=None):
         """
         The centroid of any volume V is given by

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = meshplex
-version = 0.15.11
+version = 0.15.12
 author = Nico Schlömer
 author_email = nico.schloemer@gmail.com
 description = Fast tools for simplex meshes

test/test_mesh_tetra.py

@@ -99,8 +99,6 @@ def test_regular_tet0(a):
     # cell quality
     assert is_near_equal(mesh.q_radius_ratio, [1.0], tol)
 
-    mesh.mark_boundary()
-
     assert is_near_equal(mesh.cell_barycenters, mesh.cell_centroids, tol)
     assert is_near_equal(mesh.cell_barycenters, [[0.0, 0.0, a * np.sqrt(6) / 12]], tol)