Merge pull request #103 from nschloe/fix-incenters

fix incenter computation changing the edge lengths

meshplex/helpers.py

@@ -114,7 +114,7 @@ def compute_ce_ratios(ei_dot_ej, tri_areas):
     #
     # Note that some simplifications are achieved by virtue of
     #
-    #   e1 + e2 + e3 = 0.
+    #   e0 + e1 + e2 = 0.
     #
     if numpy.any(tri_areas <= 0.0):
         raise MeshplexError("Degenerate cells.")
@@ -184,5 +184,6 @@ def compute_triangle_circumcenters(X, ei_dot_ei, ei_dot_ej):
     #     ei_dot_ei[2] * (ei_dot_ei[0] + ei_dot_ei[1] - ei_dot_ei[2]),
     # ])
     # alpha /= numpy.sum(alpha, axis=0)
+    # This last sum is the squared cell volume.
     # cc = (X[0].T * alpha[0] + X[1].T * alpha[1] + X[2].T * alpha[2]).T
     return cc

meshplex/mesh_tri.py

@@ -592,24 +592,21 @@ def cell_circumcenters(self):
     def cell_incenters(self):
         """Get the midpoints of the incircles."""
         # https://en.wikipedia.org/wiki/Incenter#Barycentric_coordinates
-        abc = self.edge_lengths
-        abc /= numpy.sum(abc, axis=0)
+        abc = self.edge_lengths / numpy.sum(self.edge_lengths, axis=0)
         return numpy.einsum("ij,jik->jk", abc, self.points[self.cells["points"]])
 
     @property
     def cell_inradius(self):
         """Get the inradii of all cells"""
         # See <http://mathworld.wolfram.com/Incircle.html>.
-        abc = self.edge_lengths
-        return 2 * self.cell_volumes / numpy.sum(abc, axis=0)
+        return 2 * self.cell_volumes / numpy.sum(self.edge_lengths, axis=0)
 
     @property
     def cell_circumradius(self):
         """Get the circumradii of all cells"""
         # See <http://mathworld.wolfram.com/Circumradius.html> and
         # <https://en.wikipedia.org/wiki/Cayley%E2%80%93Menger_determinant#Finding_the_circumradius_of_a_simplex>.
-        abc = self.edge_lengths
-        return numpy.prod(abc, axis=0) / 4 / self.cell_volumes
+        return numpy.prod(self.edge_lengths, axis=0) / 4 / self.cell_volumes
 
     @property
     def cell_quality(self):