diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
index 0000000..9505daf
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,10 @@
+# Changelog
+All notable changes to this project will be documented in this file.
+
+The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
+and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
+
+## [0.14.0] - 2020-11-05
+### Changed
+- `node_coords` is now `points`
+- `mesh_tri`: fixed inconsistent state after setting the points
diff --git a/Makefile b/Makefile
index 6f86db9..98acfb8 100644
--- a/Makefile
+++ b/Makefile
@@ -23,6 +23,7 @@ clean:
format:
isort .
black .
+ blacken-docs README.md
lint:
black --check .
diff --git a/README.md b/README.md
index 187e7a1..cec21ec 100644
--- a/README.md
+++ b/README.md
@@ -24,7 +24,7 @@ meshplex is used in [optimesh](https://github.com/nschloe/optimesh) and
### Quickstart
-```python,test
+```python
import numpy
import meshplex
@@ -46,7 +46,7 @@ print(mesh.cell_incenters)
# circumradius, inradius, cell quality, angles
print(mesh.cell_circumradius)
print(mesh.cell_inradius)
-print(mesh.cell_quality) # d * inradius / circumradius (min 0, max 1)
+print(mesh.q_radius_ratio) # d * inradius / circumradius (min 0, max 1)
print(mesh.angles)
# control volumes, centroids
@@ -73,7 +73,7 @@ classes and functions, see [readthedocs](https://meshplex.readthedocs.io/).
#### Triangles
-
+
```python
import meshplex
@@ -86,29 +86,32 @@ mesh.show(
# boundary_edge_color=None,
# comesh_color=(0.8, 0.8, 0.8),
show_axes=False,
- )
+)
```
#### Tetrahedra
-
+
```python
import numpy
import meshplex
# Generate tetrahedron
-node_coords = numpy.array(
- [
- [1.0, 0.0, -1.0 / numpy.sqrt(8)],
- [-0.5, +numpy.sqrt(3.0) / 2.0, -1.0 / numpy.sqrt(8)],
- [-0.5, -numpy.sqrt(3.0) / 2.0, -1.0 / numpy.sqrt(8)],
- [0.0, 0.0, numpy.sqrt(2.0) - 1.0 / numpy.sqrt(8)],
- ]
-) / numpy.sqrt(3.0)
+points = (
+ numpy.array(
+ [
+ [1.0, 0.0, -1.0 / numpy.sqrt(8)],
+ [-0.5, +numpy.sqrt(3.0) / 2.0, -1.0 / numpy.sqrt(8)],
+ [-0.5, -numpy.sqrt(3.0) / 2.0, -1.0 / numpy.sqrt(8)],
+ [0.0, 0.0, numpy.sqrt(2.0) - 1.0 / numpy.sqrt(8)],
+ ]
+ )
+ / numpy.sqrt(3.0)
+)
cells = [[0, 1, 2, 3]]
# Create mesh object
-mesh = meshplex.MeshTetra(node_coords, cells)
+mesh = meshplex.MeshTetra(points, cells)
# Plot cell 0 with control volume boundaries
mesh.show_cell(
@@ -120,7 +123,7 @@ mesh.show_cell(
# insphere_rgba=(1, 0, 1, 1.0),
# face_circumcenter_rgba=(0, 0, 1, 1.0),
control_volume_boundaries_rgba=(1.0, 0.0, 0.0, 1.0),
- line_width=3.0,
+ line_width=3.0,
)
```
diff --git a/experimental/ce_smoothing.py b/experimental/ce_smoothing.py
index edde5d2..ac27714 100644
--- a/experimental/ce_smoothing.py
+++ b/experimental/ce_smoothing.py
@@ -56,7 +56,7 @@ def _grad(mesh):
grad_stack = numpy.array([grad_x0, grad_x1, grad_x2])
# add up all the contributions
- grad = numpy.zeros(mesh.node_coords.shape)
+ grad = numpy.zeros(mesh.points.shape)
numpy.add.at(grad, mesh.cells["nodes"].T, grad_stack)
return grad
@@ -67,9 +67,9 @@ def smooth(mesh, t=1.0e-3, num_iters=10):
for k in range(num_iters):
# mesh = mesh_tri.flip_until_delaunay(mesh)
- x = mesh.node_coords.copy()
+ x = mesh.points.copy()
x -= t * _grad(mesh)
- x[boundary_verts] = mesh.node_coords[boundary_verts]
+ x[boundary_verts] = mesh.points[boundary_verts]
mesh = meshplex.MeshTri(x, mesh.cells["nodes"])
mesh.write("smoo%04d.vtu" % k)
print(_objective(mesh))
@@ -79,11 +79,11 @@ def smooth(mesh, t=1.0e-3, num_iters=10):
def read(filename):
mesh = meshio.read(filename)
- # x = mesh.node_coords.copy()
+ # x = mesh.points.copy()
# x[:, :2] += 1.0e-1 * (
- # numpy.random.rand(len(mesh.node_coords), 2) - 0.5
+ # numpy.random.rand(len(mesh.points), 2) - 0.5
# )
- # x[boundary_verts] = mesh.node_coords[boundary_verts]
+ # x[boundary_verts] = mesh.points[boundary_verts]
# only include nodes which are part of a cell
uvertices, uidx = numpy.unique(mesh.get_cells_type("triangle"), return_inverse=True)
diff --git a/meshplex/base.py b/meshplex/base.py
index 4fc2846..a57e885 100644
--- a/meshplex/base.py
+++ b/meshplex/base.py
@@ -1,185 +1,34 @@
import meshio
import numpy
-from .exceptions import MeshplexError
-
__all__ = []
-def compute_tri_areas(ei_dot_ej):
- vol2 = 0.25 * (
- ei_dot_ej[2] * ei_dot_ej[0]
- + ei_dot_ej[0] * ei_dot_ej[1]
- + 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.
- assert numpy.all(vol2 > -1.0e-14)
- vol2[vol2 < 0] = 0.0
- return numpy.sqrt(vol2)
-
-
-def compute_ce_ratios(ei_dot_ej, tri_areas):
- """Given triangles (specified by their mutual edge projections and the area), this
- routine will return ratio of the signed distances of the triangle circumcenters to
- the edge midpoints and the edge lengths.
- """
- # The input argument are the dot products
- #
- #
- #
- #
- #
- # of the edges
- #
- # e0: x1->x2,
- # e1: x2->x0,
- # e2: x0->x1.
- #
- # Note that edge e_i is opposite of node i and the edges add up to 0.
- # (Those quantities can be shared between numerous methods, so share them.)
- #
-
- # There are multiple ways of deriving a closed form for the
- # covolume-edgelength ratios.
- #
- # * The covolume-edge ratios for the edges of each cell is the solution
- # of the equation system
- #
- # |simplex| ||u||^2 = \sum_i \alpha_i ,
- #
- # where alpha_i are the covolume contributions for the edges. This
- # equation system holds for all vectors u in the plane spanned by the
- # edges, particularly by the edges themselves.
- #
- # For triangles, the exact solution of the system is
- #
- # x_1 = / * |simplex|;
- #
- # see .
- #
- # * In trilinear coordinates
- # , the
- # circumcenter is
- #
- # cos(alpha0) : cos(alpha1) : cos(alpha2)
- #
- # where the alpha_i are the angles opposite of the respective edge.
- # With
- #
- # cos(alpha0) = / ||e1|| / ||e2||
- #
- # and the conversion formula to Cartesian coordinates, ones gets the
- # expression
- #
- # ce0 = * 0.5 / sqrt(alpha)
- #
- # with
- #
- # alpha = *
- # + *
- # + * .
- #
- # Note that some simplifications are achieved by virtue of
- #
- # e1 + e2 + e3 = 0.
- #
- if not numpy.all(tri_areas > 0.0):
- raise MeshplexError("Degenerate cell.")
-
- return -ei_dot_ej * 0.25 / tri_areas[None]
-
-
-def compute_triangle_circumcenters(X, ei_dot_ei, ei_dot_ej):
- """Computes the circumcenters of all given triangles."""
- # The input argument are the dot products
- #
- #
- #
- #
- #
- # of the edges
- #
- # e0: x1->x2,
- # e1: x2->x0,
- # e2: x0->x1.
- #
- # Note that edge e_i is opposite of node i and the edges add up to 0.
-
- # The trilinear coordinates of the circumcenter are
- #
- # cos(alpha0) : cos(alpha1) : cos(alpha2)
- #
- # where alpha_k is the angle at point k, opposite of edge k. The Cartesian
- # coordinates are (see
- # )
- #
- # C = sum_i ||e_i|| * cos(alpha_i)/beta * P_i
- #
- # with
- #
- # beta = sum ||e_i||*cos(alpha_i)
- #
- # Incidentally, the cosines are
- #
- # cos(alpha0) = / ||e1|| / ||e2||,
- #
- # so in total
- #
- # C = / sum_i ( ) P0
- # + ... P1
- # + ... P2.
- #
- # An even nicer formula is given on
- # : The
- # barycentric coordinates of the circumcenter are
- #
- # a^2 (b^2 + c^2 - a^2) : b^2 (c^2 + a^2 - b^2) : c^2 (a^2 + b^2 - c^2).
- #
- # This is only using the squared edge lengths, too!
- #
- alpha = ei_dot_ei * ei_dot_ej
- alpha_sum = alpha[0] + alpha[1] + alpha[2]
- beta = alpha / alpha_sum[None]
- a = X * beta[..., None]
- cc = a[0] + a[1] + a[2]
-
- # alpha = numpy.array([
- # ei_dot_ei[0] * (ei_dot_ei[1] + ei_dot_ei[2] - ei_dot_ei[0]),
- # ei_dot_ei[1] * (ei_dot_ei[2] + ei_dot_ei[0] - ei_dot_ei[1]),
- # ei_dot_ei[2] * (ei_dot_ei[0] + ei_dot_ei[1] - ei_dot_ei[2]),
- # ])
- # alpha /= numpy.sum(alpha, axis=0)
- # cc = (X[0].T * alpha[0] + X[1].T * alpha[1] + X[2].T * alpha[2]).T
- return cc
-
-
-class _base_mesh:
- def __init__(self, nodes, cells_nodes):
- self.node_coords = nodes
+class _BaseMesh:
+ def __init__(self, points, cells_points):
self._edge_lengths = None
def write(self, filename, point_data=None, cell_data=None, field_data=None):
- if self.node_coords.shape[1] == 2:
- n = len(self.node_coords)
+ if self.points.shape[1] == 2:
+ n = len(self.points)
a = numpy.ascontiguousarray(
- numpy.column_stack([self.node_coords, numpy.zeros(n)])
+ numpy.column_stack([self.points, numpy.zeros(n)])
)
else:
- a = self.node_coords
+ a = self.points
- if self.cells["nodes"].shape[1] == 3:
+ if self.cells["points"].shape[1] == 3:
cell_type = "triangle"
else:
assert (
- self.cells["nodes"].shape[1] == 4
+ self.cells["points"].shape[1] == 4
), "Only triangles/tetrahedra supported"
cell_type = "tetra"
meshio.write_points_cells(
filename,
a,
- {cell_type: self.cells["nodes"]},
+ {cell_type: self.cells["points"]},
point_data=point_data,
cell_data=cell_data,
field_data=field_data,
@@ -209,7 +58,7 @@ def get_edge_mask(self, subdomain=None):
self._mark_vertices(subdomain)
# A face is inside if all its edges are in.
- # An edge is inside if all its nodes are in.
+ # An edge is inside if all its points are in.
is_in = self.subdomains[subdomain]["vertices"][self.idx_hierarchy]
# Take `all()` over the first index
is_inside = numpy.all(is_in, axis=tuple(range(1)))
@@ -230,7 +79,7 @@ def get_face_mask(self, subdomain):
self._mark_vertices(subdomain)
# A face is inside if all its edges are in.
- # An edge is inside if all its nodes are in.
+ # An edge is inside if all its points are in.
is_in = self.subdomains[subdomain]["vertices"][self.idx_hierarchy]
# Take `all()` over all axes except the last two (face_ids, cell_ids).
n = len(is_in.shape)
@@ -262,14 +111,14 @@ def get_cell_mask(self, subdomain=None):
def _mark_vertices(self, subdomain):
"""Mark faces/edges which are fully in subdomain."""
if subdomain is None:
- is_inside = numpy.ones(len(self.node_coords), dtype=bool)
+ is_inside = numpy.ones(len(self.points), dtype=bool)
else:
- is_inside = subdomain.is_inside(self.node_coords.T).T
+ is_inside = subdomain.is_inside(self.points.T).T
if subdomain.is_boundary_only:
# Filter boundary
self.mark_boundary()
- is_inside = is_inside & self.is_boundary_node
+ is_inside = is_inside & self.is_boundary_point
self.subdomains[subdomain] = {"vertices": is_inside}
return
diff --git a/meshplex/helpers.py b/meshplex/helpers.py
index 3d5d0db..1828c0f 100644
--- a/meshplex/helpers.py
+++ b/meshplex/helpers.py
@@ -2,6 +2,8 @@
import numpy
+from .exceptions import MeshplexError
+
def get_signed_simplex_volumes(cells, pts):
"""Signed volume of a simplex in nD. Note that signing only makes sense for
@@ -35,3 +37,151 @@ def unique_rows(a):
a_unique, inv, cts = numpy.unique(b, return_inverse=True, return_counts=True)
a_unique = a_unique.view(a.dtype).reshape(-1, a.shape[1])
return a_unique, inv, cts
+
+
+def compute_tri_areas(ei_dot_ej):
+ vol2 = 0.25 * (
+ ei_dot_ej[2] * ei_dot_ej[0]
+ + ei_dot_ej[0] * ei_dot_ej[1]
+ + 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.
+ assert numpy.all(vol2 > -1.0e-14)
+ vol2[vol2 < 0] = 0.0
+ return numpy.sqrt(vol2)
+
+
+def compute_ce_ratios(ei_dot_ej, tri_areas):
+ """Given triangles (specified by their mutual edge projections and the area), this
+ routine will return ratio of the signed distances of the triangle circumcenters to
+ the edge midpoints and the edge lengths.
+ """
+ # The input argument are the dot products
+ #
+ #
+ #
+ #
+ #
+ # of the edges
+ #
+ # e0: x1->x2,
+ # e1: x2->x0,
+ # e2: x0->x1.
+ #
+ # Note that edge e_i is opposite of node i and the edges add up to 0.
+ # (Those quantities can be shared between numerous methods, so share them.)
+ #
+
+ # There are multiple ways of deriving a closed form for the
+ # covolume-edgelength ratios.
+ #
+ # * The covolume-edge ratios for the edges of each cell is the solution
+ # of the equation system
+ #
+ # |simplex| ||u||^2 = \sum_i \alpha_i ,
+ #
+ # where alpha_i are the covolume contributions for the edges. This
+ # equation system holds for all vectors u in the plane spanned by the
+ # edges, particularly by the edges themselves.
+ #
+ # For triangles, the exact solution of the system is
+ #
+ # x_1 = / * |simplex|;
+ #
+ # see .
+ #
+ # * In trilinear coordinates
+ # , the
+ # circumcenter is
+ #
+ # cos(alpha0) : cos(alpha1) : cos(alpha2)
+ #
+ # where the alpha_i are the angles opposite of the respective edge.
+ # With
+ #
+ # cos(alpha0) = / ||e1|| / ||e2||
+ #
+ # and the conversion formula to Cartesian coordinates, ones gets the
+ # expression
+ #
+ # ce0 = * 0.5 / sqrt(alpha)
+ #
+ # with
+ #
+ # alpha = *
+ # + *
+ # + * .
+ #
+ # Note that some simplifications are achieved by virtue of
+ #
+ # e1 + e2 + e3 = 0.
+ #
+ if not numpy.all(tri_areas > 0.0):
+ raise MeshplexError("Degenerate cell.")
+
+ return -ei_dot_ej * 0.25 / tri_areas[None]
+
+
+def compute_triangle_circumcenters(X, ei_dot_ei, ei_dot_ej):
+ """Computes the circumcenters of all given triangles."""
+ # The input argument are the dot products
+ #
+ #
+ #
+ #
+ #
+ # of the edges
+ #
+ # e0: x1->x2,
+ # e1: x2->x0,
+ # e2: x0->x1.
+ #
+ # Note that edge e_i is opposite of node i and the edges add up to 0.
+
+ # The trilinear coordinates of the circumcenter are
+ #
+ # cos(alpha0) : cos(alpha1) : cos(alpha2)
+ #
+ # where alpha_k is the angle at point k, opposite of edge k. The Cartesian
+ # coordinates are (see
+ # )
+ #
+ # C = sum_i ||e_i|| * cos(alpha_i)/beta * P_i
+ #
+ # with
+ #
+ # beta = sum ||e_i||*cos(alpha_i)
+ #
+ # Incidentally, the cosines are
+ #
+ # cos(alpha0) = / ||e1|| / ||e2||,
+ #
+ # so in total
+ #
+ # C = / sum_i ( ) P0
+ # + ... P1
+ # + ... P2.
+ #
+ # An even nicer formula is given on
+ # : The
+ # barycentric coordinates of the circumcenter are
+ #
+ # a^2 (b^2 + c^2 - a^2) : b^2 (c^2 + a^2 - b^2) : c^2 (a^2 + b^2 - c^2).
+ #
+ # This is only using the squared edge lengths, too!
+ #
+ alpha = ei_dot_ei * ei_dot_ej
+ alpha_sum = alpha[0] + alpha[1] + alpha[2]
+ beta = alpha / alpha_sum[None]
+ a = X * beta[..., None]
+ cc = a[0] + a[1] + a[2]
+
+ # alpha = numpy.array([
+ # ei_dot_ei[0] * (ei_dot_ei[1] + ei_dot_ei[2] - ei_dot_ei[0]),
+ # ei_dot_ei[1] * (ei_dot_ei[2] + ei_dot_ei[0] - ei_dot_ei[1]),
+ # ei_dot_ei[2] * (ei_dot_ei[0] + ei_dot_ei[1] - ei_dot_ei[2]),
+ # ])
+ # alpha /= numpy.sum(alpha, axis=0)
+ # cc = (X[0].T * alpha[0] + X[1].T * alpha[1] + X[2].T * alpha[2]).T
+ return cc
diff --git a/meshplex/mesh_line.py b/meshplex/mesh_line.py
index b264eea..26a1567 100644
--- a/meshplex/mesh_line.py
+++ b/meshplex/mesh_line.py
@@ -4,8 +4,8 @@
class MeshLine:
"""Class for handling line segment "meshes"."""
- def __init__(self, node_coords, cells):
- self.node_coords = node_coords
+ def __init__(self, points, cells):
+ self.points = points
num_cells = len(cells)
self.cells = numpy.empty(num_cells, dtype=numpy.dtype([("nodes", (int, 2))]))
@@ -19,10 +19,7 @@ def create_cell_volumes(self):
"""Computes the volumes of the "cells" in the mesh."""
self.cell_volumes = numpy.array(
[
- abs(
- self.node_coords[cell["nodes"]][1]
- - self.node_coords[cell["nodes"]][0]
- )
+ abs(self.points[cell["nodes"]][1] - self.points[cell["nodes"]][0])
for cell in self.cells
]
)
@@ -30,7 +27,7 @@ def create_cell_volumes(self):
def create_control_volumes(self):
"""Compute the control volumes of all nodes in the mesh."""
- self.control_volumes = numpy.zeros(len(self.node_coords), dtype=float)
+ self.control_volumes = numpy.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]
@@ -46,5 +43,5 @@ def show_vertex_function(self, u):
""""""
from matplotlib import pyplot as plt
- plt.plot(self.node_coords, u)
+ plt.plot(self.points, u)
return
diff --git a/meshplex/mesh_tetra.py b/meshplex/mesh_tetra.py
index f598b10..e311349 100644
--- a/meshplex/mesh_tetra.py
+++ b/meshplex/mesh_tetra.py
@@ -1,17 +1,18 @@
import numpy
-from .base import _base_mesh, compute_tri_areas, compute_triangle_circumcenters
+from .base import _BaseMesh
+from .helpers import compute_tri_areas, compute_triangle_circumcenters
__all__ = ["MeshTetra"]
# pylint: disable=too-many-instance-attributes
-class MeshTetra(_base_mesh):
+class MeshTetra(_BaseMesh):
"""Class for handling tetrahedral meshes."""
- def __init__(self, node_coords, cells):
+ def __init__(self, points, cells):
"""Initialization."""
- # Sort cells and nodes, first every row, then the rows themselves. This helps in
+ # 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
@@ -21,39 +22,40 @@ def __init__(self, node_coords, cells):
cells = numpy.sort(cells, axis=1)
cells = cells[cells[:, 0].argsort()]
- super().__init__(node_coords, cells)
+ self._points = points
+ super().__init__(points, cells)
# Assert that all vertices are used.
# If there are vertices which do not appear in the cells list, this
# ```
# uvertices, uidx = numpy.unique(cells, return_inverse=True)
# cells = uidx.reshape(cells.shape)
- # nodes = nodes[uvertices]
+ # points = points[uvertices]
# ```
# helps.
- is_used = numpy.zeros(len(node_coords), dtype=bool)
+ is_used = numpy.zeros(len(points), dtype=bool)
is_used[cells] = True
- assert numpy.all(is_used), "There are {} dangling nodes in the mesh".format(
+ assert numpy.all(is_used), "There are {} dangling points in the mesh".format(
numpy.sum(~is_used)
)
- self.cells = {"nodes": cells}
+ self.cells = {"points": cells}
self._control_volumes = None
self._circumcenters = None
self.subdomains = {}
- # Arrange the node_face_cells such that node k is opposite of face k in each
+ # Arrange the point_face_cells such that point k is opposite of face k in each
# cell.
- nds = self.cells["nodes"].T
+ nds = self.cells["points"].T
idx = numpy.array([[1, 2, 3], [2, 3, 0], [3, 0, 1], [0, 1, 2]]).T
- self.node_face_cells = nds[idx]
+ self.point_face_cells = nds[idx]
- # Arrange the idx_hierarchy (node->edge->face->cells) such that
+ # Arrange the idx_hierarchy (point->edge->face->cells) such that
#
- # * node k is opposite of edge k in each face,
+ # * point k is opposite of edge k in each face,
# * duplicate edges are in the same spot of the each of the faces,
- # * all nodes are in domino order ([1, 2], [2, 3], [3, 1]),
+ # * all points are in domino order ([1, 2], [2, 3], [3, 1]),
# * the same edges are easy to find:
# - edge 0: face+1, edge 2
# - edge 1: face+2, edge 1
@@ -74,17 +76,16 @@ def __init__(self, node_coords, cells):
self.idx_hierarchy = nds[self.local_idx]
# The inverted local index.
- # This array specifies for each of the three nodes which edge endpoints
+ # This array specifies for each of the three points which edge endpoints
# correspond to it.
self.local_idx_inv = [
- [tuple(i) for i in zip(*numpy.where(self.local_idx == node_idx))]
- for node_idx in range(4)
+ [tuple(i) for i in zip(*numpy.where(self.local_idx == point_idx))]
+ for point_idx in range(4)
]
# create ei_dot_ei, ei_dot_ej
self.edge_coords = (
- self.node_coords[self.idx_hierarchy[1]]
- - self.node_coords[self.idx_hierarchy[0]]
+ 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
@@ -101,7 +102,7 @@ def __init__(self, node_coords, cells):
self.ce_ratios = self._compute_ce_ratios_geometric()
# self.ce_ratios = self._compute_ce_ratios_algebraic()
- self.is_boundary_node = None
+ self.is_boundary_point = None
self._inv_faces = None
self.edges = None
self.is_boundary_facet_individual = None
@@ -112,11 +113,21 @@ def __init__(self, node_coords, cells):
return
def __repr__(self):
- num_nodes = len(self.node_coords)
- num_cells = len(self.cells["nodes"])
- string = f""
+ num_points = len(self.points)
+ num_cells = len(self.cells["points"])
+ string = f""
return string
+ @property
+ def points(self):
+ return self._points
+
+ @points.setter
+ def set_points(self, new_points):
+ assert new_points.shape == self._points.shape
+ # reset all computed values
+ self.ei_do_ei = None
+
def get_ce_ratios(self):
"""Covolume-edge length ratios."""
return self.ce_ratios
@@ -125,14 +136,14 @@ def mark_boundary(self):
if "faces" not in self.cells:
self.create_cell_face_relationships()
- self.is_boundary_node = numpy.zeros(len(self.node_coords), dtype=bool)
- self.is_boundary_node[
- self.faces["nodes"][self.is_boundary_facet_individual]
+ 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
def create_cell_face_relationships(self):
- # Reshape into individual faces, and take the first node per edge. (The face is
+ # Reshape into individual faces, and take the first point per edge. (The face is
# fully characterized by it.) Sort the columns to make it possible for
# `unique()` to identify individual faces.
s = self.idx_hierarchy.shape
@@ -154,15 +165,15 @@ def create_cell_face_relationships(self):
self.is_boundary_facet = (cts[inv] == 1).reshape(s[2:])
self.is_boundary_facet_individual = cts == 1
- self.faces = {"nodes": a[idx]}
+ self.faces = {"points": a[idx]}
# cell->faces relationship
- num_cells = len(self.cells["nodes"])
+ num_cells = len(self.cells["points"])
cells_faces = inv.reshape([4, num_cells]).T
self.cells["faces"] = cells_faces
- # Store the opposing nodes too
- self.cells["opposing vertex"] = self.cells["nodes"]
+ # Store the opposing points too
+ self.cells["opposing vertex"] = self.cells["points"]
# save for create_edge_cells
self._inv_faces = inv
@@ -172,9 +183,9 @@ def create_cell_face_relationships(self):
def create_face_edge_relationships(self):
a = numpy.vstack(
[
- self.faces["nodes"][:, [1, 2]],
- self.faces["nodes"][:, [2, 0]],
- self.faces["nodes"][:, [0, 1]],
+ self.faces["points"][:, [1, 2]],
+ self.faces["points"][:, [2, 0]],
+ self.faces["points"][:, [0, 1]],
]
)
@@ -183,12 +194,12 @@ def create_face_edge_relationships(self):
numpy.dtype((numpy.void, a.dtype.itemsize * a.shape[1]))
)
_, idx, inv = numpy.unique(b, return_index=True, return_inverse=True)
- edge_nodes = a[idx]
+ edge_points = a[idx]
- self.edges = {"nodes": edge_nodes}
+ self.edges = {"points": edge_points}
# face->edge relationship
- num_faces = len(self.faces["nodes"])
+ num_faces = len(self.faces["points"])
face_edges = inv.reshape([3, num_faces]).T
self.faces["edges"] = face_edges
@@ -222,7 +233,7 @@ def _compute_cell_circumcenters(self):
alpha = self._zeta / numpy.sum(self._zeta, axis=0)
self._circumcenters = numpy.sum(
- alpha[None].T * self.node_coords[self.cells["nodes"]], axis=1
+ alpha[None].T * self.points[self.cells["points"]], axis=1
)
return
@@ -235,8 +246,8 @@ def _compute_cell_circumcenters(self):
# def _compute_ce_ratios_algebraic(self):
# # Precompute edges.
# half_edges = (
- # self.node_coords[self.idx_hierarchy[1]]
- # - self.node_coords[self.idx_hierarchy[0]]
+ # self.points[self.idx_hierarchy[1]]
+ # - self.points[self.idx_hierarchy[0]]
# )
# # Build the equation system:
@@ -294,7 +305,7 @@ def _compute_ce_ratios_geometric(self):
# + self.ei_dot_ej[0] * self.ei_dot_ej[1] * self.ei_dot_ej[2]
# ).
#
- # for the face [1, 2, 3] (with edges [3, 4, 5]), where nodes and edges are
+ # for the face [1, 2, 3] (with edges [3, 4, 5]), where points and edges are
# ordered like
#
# 3
@@ -386,7 +397,7 @@ def cell_centroids(self):
"""The centroids (barycenters, midpoints of the circumcircles) of all tetrahedra."""
if self._cell_centroids is None:
self._cell_centroids = (
- numpy.sum(self.node_coords[self.cells["nodes"]], axis=1) / 4.0
+ numpy.sum(self.points[self.cells["points"]], axis=1) / 4.0
)
return self._cell_centroids
@@ -408,9 +419,7 @@ def cell_incenters(self):
face_areas = compute_tri_areas(self.ei_dot_ej)
# abc = numpy.sqrt(self.ei_dot_ei)
face_areas /= numpy.sum(face_areas, axis=0)
- return numpy.einsum(
- "ij,jik->jk", face_areas, self.node_coords[self.cells["nodes"]]
- )
+ return numpy.einsum("ij,jik->jk", face_areas, self.points[self.cells["points"]])
@property
def cell_inradius(self):
@@ -420,8 +429,8 @@ def cell_inradius(self):
@property
def cell_circumradius(self):
- # # Just take the distance of the circumcenter to one of the nodes for now.
- dist = self.node_coords[self.idx_hierarchy[0, 0, 0]] - self.cell_circumcenters
+ # # Just take the distance of the circumcenter to one of the points for now.
+ dist = self.points[self.idx_hierarchy[0, 0, 0]] - self.cell_circumcenters
circumradius = numpy.sqrt(numpy.einsum("ij,ij->i", dist, dist))
# https://en.wikipedia.org/wiki/Tetrahedron#Circumradius
#
@@ -518,13 +527,13 @@ def q_vol_rms_edgelength3(self):
@property
def control_volumes(self):
- """Compute the control volumes of all nodes in the mesh."""
+ """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.0
# Explicitly sum up contributions per cell first. Makes numpy.add.at faster.
- # For every node k (range(4)), check for which edges k appears in local_idx,
+ # For every point k (range(4)), check for which edges k appears in local_idx,
# and sum() up the v's from there.
vals = numpy.array(
[
@@ -536,9 +545,9 @@ def control_volumes(self):
).T
#
self._control_volumes = numpy.bincount(
- self.cells["nodes"].reshape(-1),
+ self.cells["points"].reshape(-1),
vals.reshape(-1),
- minlength=len(self.node_coords),
+ minlength=len(self.points),
)
return self._control_volumes
@@ -580,11 +589,11 @@ def plot(self):
if self._circumcenters is None:
self._compute_cell_circumcenters()
- X = self.node_coords
- for cell_id in range(len(self.cells["nodes"])):
+ X = self.points
+ for cell_id in range(len(self.cells["points"])):
cc = self._circumcenters[cell_id]
#
- x = X[self.node_face_cells[..., [cell_id]]]
+ x = X[self.point_face_cells[..., [cell_id]]]
face_ccs = compute_triangle_circumcenters(x, self.ei_dot_ei, self.ei_dot_ej)
# draw the face circumcenters
ax.plot(
@@ -651,20 +660,20 @@ def plot_edge(self, edge_id):
)
col = "k"
for adj_edge_id in adj_edge_ids:
- x = self.node_coords[self.edges["nodes"][adj_edge_id]]
+ x = self.points[self.edges["points"][adj_edge_id]]
ax.plot(x[:, 0], x[:, 1], x[:, 2], col)
# make clear which is edge_id
- x = self.node_coords[self.edges["nodes"][edge_id]]
+ x = self.points[self.edges["points"][edge_id]]
ax.plot(x[:, 0], x[:, 1], x[:, 2], color=col, linewidth=3.0)
# connect the face circumcenters with the corresponding cell
# circumcenters
- X = self.node_coords
+ X = self.points
for cell_id in adj_cell_ids:
cc = self._circumcenters[cell_id]
#
- x = X[self.node_face_cells[..., [cell_id]]]
+ x = X[self.point_face_cells[..., [cell_id]]]
face_ccs = compute_triangle_circumcenters(x, self.ei_dot_ei, self.ei_dot_ej)
# draw the face circumcenters
ax.plot(
@@ -749,7 +758,7 @@ def get_sphere_actor(x0, r, rgba):
render_window_interactor.SetRenderWindow(render_window)
for ij in [[0, 1], [0, 2], [0, 3], [1, 2], [1, 3], [2, 3]]:
- x0, x1 = self.node_coords[self.cells["nodes"][cell_id][ij]]
+ x0, x1 = self.points[self.cells["points"][cell_id][ij]]
renderer.AddActor(get_line_actor(x0, x1, line_width))
renderer.SetBackground(1.0, 1.0, 1.0)
@@ -789,7 +798,7 @@ def get_sphere_actor(x0, r, rgba):
)
if face_circumcenter_rgba is not None:
- x = self.node_coords[self.node_face_cells[..., [cell_id]]]
+ x = self.points[self.point_face_cells[..., [cell_id]]]
face_ccs = compute_triangle_circumcenters(
x, self.ei_dot_ei, self.ei_dot_ej
)[:, 0, :]
@@ -798,14 +807,14 @@ def get_sphere_actor(x0, r, rgba):
if control_volume_boundaries_rgba:
cell_cc = self.cell_circumcenters[cell_id]
- x = self.node_coords[self.node_face_cells[..., [cell_id]]]
+ x = self.points[self.point_face_cells[..., [cell_id]]]
face_ccs = compute_triangle_circumcenters(
x, self.ei_dot_ei, self.ei_dot_ej
)[:, 0, :]
for face, face_cc in zip(range(4), face_ccs):
for edge in range(3):
k0, k1 = self.idx_hierarchy[:, edge, face, cell_id]
- edge_midpoint = 0.5 * (self.node_coords[k0] + self.node_coords[k1])
+ edge_midpoint = 0.5 * (self.points[k0] + self.points[k1])
points = vtk.vtkPoints()
points.InsertNextPoint(*edge_midpoint)
diff --git a/meshplex/mesh_tri.py b/meshplex/mesh_tri.py
index c5f8fe7..892b1f3 100644
--- a/meshplex/mesh_tri.py
+++ b/meshplex/mesh_tri.py
@@ -3,54 +3,54 @@
import numpy
-from .base import (
- _base_mesh,
+from .base import _BaseMesh
+from .helpers import (
compute_ce_ratios,
compute_tri_areas,
compute_triangle_circumcenters,
+ grp_start_len,
+ unique_rows,
)
-from .helpers import grp_start_len, unique_rows
__all__ = ["MeshTri"]
-class MeshTri(_base_mesh):
+class MeshTri(_BaseMesh):
"""Class for handling triangular meshes."""
- def __init__(self, nodes, cells, sort_cells=False):
+ def __init__(self, points, cells, sort_cells=False):
"""Initialization."""
if sort_cells:
- # Sort cells and nodes, 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, 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()]
- assert len(nodes.shape) == 2, "Illegal node coordinates shape {}".format(
- nodes.shape
- )
+ assert len(points.shape) == 2, f"Illegal point coordinates shape {points.shape}"
assert (
len(cells.shape) == 2 and cells.shape[1] == 3
), f"Illegal cells shape {cells.shape}"
- super().__init__(nodes, cells)
+ self._points = points
+ super().__init__(points, cells)
# Assert that all vertices are used.
# If there are vertices which do not appear in the cells list, this
# ```
# uvertices, uidx = numpy.unique(cells, return_inverse=True)
# cells = uidx.reshape(cells.shape)
- # nodes = nodes[uvertices]
+ # points = points[uvertices]
# ```
# helps.
- self.node_is_used = numpy.zeros(len(nodes), dtype=bool)
- self.node_is_used[cells] = True
+ self.point_is_used = numpy.zeros(len(points), dtype=bool)
+ self.point_is_used[cells] = True
- self.cells = {"nodes": cells}
+ self.cells = {"points": cells}
self._interior_ce_ratios = None
self._control_volumes = None
@@ -62,8 +62,8 @@ def __init__(self, nodes, cells, sort_cells=False):
self._cell_circumcenters = None
self._signed_cell_areas = None
self.subdomains = {}
- self._is_interior_node = None
- self._is_boundary_node = None
+ self._is_interior_point = None
+ self._is_boundary_point = None
self.is_boundary_edge = None
self._is_boundary_facet = None
self._interior_edge_lengths = None
@@ -74,19 +74,19 @@ def __init__(self, nodes, cells, sort_cells=False):
self._cell_centroids = None
# compute data
- # Create the idx_hierarchy (nodes->edges->cells), i.e., the value of
- # `self.idx_hierarchy[0, 2, 27]` is the index of the node of cell 27, edge 2,
- # node 0. The shape of `self.idx_hierarchy` is `(2, 3, n)`, where `n` is the
+ # Create the idx_hierarchy (points->edges->cells), i.e., the value of
+ # `self.idx_hierarchy[0, 2, 27]` is the index of the point of cell 27, edge 2,
+ # point 0. The shape of `self.idx_hierarchy` is `(2, 3, n)`, where `n` is the
# number of cells. Make sure that the k-th edge is opposite of the k-th point in
# the triangle.
self.local_idx = numpy.array([[1, 2], [2, 0], [0, 1]]).T
- # Map idx back to the nodes. This is useful if quantities which are in idx shape
- # need to be added up into nodes (e.g., equation system rhs).
- nds = self.cells["nodes"].T
+ # Map idx back to the points. This is useful if quantities which are in idx shape
+ # need to be added up into points (e.g., equation system rhs).
+ nds = self.cells["points"].T
self.idx_hierarchy = nds[self.local_idx]
# The inverted local index.
- # This array specifies for each of the three nodes which edge endpoints
+ # This array specifies for each of the three points which edge endpoints
# correspond to it. For the above local_idx, this should give
#
# [[(1, 1), (0, 2)], [(0, 0), (1, 2)], [(1, 0), (0, 1)]]
@@ -97,8 +97,7 @@ def __init__(self, nodes, cells, sort_cells=False):
# Create the corresponding edge coordinates.
self.half_edge_coords = (
- self.node_coords[self.idx_hierarchy[1]]
- - self.node_coords[self.idx_hierarchy[0]]
+ self.points[self.idx_hierarchy[1]] - self.points[self.idx_hierarchy[0]]
)
# einsum is faster if the tail survives, e.g., ijk,ijk->jk.
@@ -110,8 +109,9 @@ def __init__(self, nodes, cells, sort_cells=False):
self.half_edge_coords[[2, 0, 1]],
)
- e = self.half_edge_coords
- self.ei_dot_ei = numpy.einsum("ijk, ijk->ij", e, e)
+ self.ei_dot_ei = numpy.einsum(
+ "ijk, ijk->ij", self.half_edge_coords, self.half_edge_coords
+ )
self.cell_volumes = compute_tri_areas(self.ei_dot_ej)
@@ -120,37 +120,64 @@ def __init__(self, nodes, cells, sort_cells=False):
# flat_local_edge, self.flat_cells = numpy.where(is_flat_halfedge)
# self.is_flat_cell = numpy.any(is_flat_halfedge, axis=0)
# self.fcc = FlatCellCorrector(
- # self.cells["nodes"][self.fcc_cells], flat_local_edge, self.node_coords
+ # self.cells["points"][self.fcc_cells], flat_local_edge, self.points
# )
# self._ce_ratios[:, self.fcc_cells] = self.fcc.ce_ratios.T
def __repr__(self):
- num_nodes = len(self.node_coords)
- num_cells = len(self.cells["nodes"])
- string = f""
+ num_points = len(self.points)
+ num_cells = len(self.cells["points"])
+ string = f""
return string
- # def update_node_coordinates(self, X):
- # assert X.shape == self.node_coords.shape
- # self.node_coords = X
+ # def update_point_coordinates(self, X):
+ # assert X.shape == self.points.shape
+ # self.points = X
# self._update_values()
# return
- # def update_interior_node_coordinates(self, X):
- # assert X.shape == self.node_coords[self.is_interior_node].shape
- # self.node_coords[self.is_interior_node] = X
+ # def update_interior_point_coordinates(self, X):
+ # assert X.shape == self.points[self.is_interior_point].shape
+ # self.points[self.is_interior_point] = X
# self.update_values()
# return
+ @property
+ def points(self):
+ return self._points
+
+ @points.setter
+ def points(self, new_points):
+ new_points = numpy.asarray(new_points)
+ assert new_points.shape == self._points.shape
+ self._points = new_points
+ # reset all computed values
+ self.ei_dot_ei = None
+ self.half_edge_coords = None
+ self.ei_dot_ej = None
+ self.ei_dot_ei = None
+ self.cell_volumes = None
+ self._ce_ratios = None
+ self._interior_edge_lengths = None
+ self._cell_circumcenters = None
+ self._interior_ce_ratios = None
+ self._control_volumes = None
+ self._cell_partitions = None
+ self._cv_centroids = None
+ self._cvc_cell_mask = None
+ self._surface_areas = None
+ self._signed_cell_areas = None
+ self._cell_centroids = None
+
@property
def euler_characteristic(self):
# number of vertices - number of edges + number of faces
if "edges" not in self.cells:
self.create_edges()
return (
- self.node_coords.shape[0]
- - self.edges["nodes"].shape[0]
- + self.cells["nodes"].shape[0]
+ self.points.shape[0]
+ - self.edges["points"].shape[0]
+ + self.cells["points"].shape[0]
)
@property
@@ -165,28 +192,24 @@ def ce_ratios(self):
return self._ce_ratios
def update_values(self):
- """Update all computes entities around the mesh."""
- if self.half_edge_coords is not None:
- # Constructing the temporary arrays
- # self.node_coords[self.idx_hierarchy] can take quite a while here.
- self.half_edge_coords = (
- self.node_coords[self.idx_hierarchy[1]]
- - self.node_coords[self.idx_hierarchy[0]]
- )
+ """Update all computed values of the mesh (edge lengths etc.)."""
+ # Constructing the temporary arrays
+ # self.points[self.idx_hierarchy] can take quite a while here.
+ self.half_edge_coords = (
+ self.points[self.idx_hierarchy[1]] - self.points[self.idx_hierarchy[0]]
+ )
- if self.ei_dot_ej is not None:
- self.ei_dot_ej = numpy.einsum(
- "ijk, ijk->ij",
- self.half_edge_coords[[1, 2, 0]],
- self.half_edge_coords[[2, 0, 1]],
- )
+ self.ei_dot_ej = numpy.einsum(
+ "ijk, ijk->ij",
+ self.half_edge_coords[[1, 2, 0]],
+ self.half_edge_coords[[2, 0, 1]],
+ )
- if self.ei_dot_ei is not None:
- e = self.half_edge_coords
- self.ei_dot_ei = numpy.einsum("ijk, ijk->ij", e, e)
+ self.ei_dot_ei = numpy.einsum(
+ "ijk, ijk->ij", self.half_edge_coords, self.half_edge_coords
+ )
- if self.cell_volumes is not None or self.ce_ratios is not None:
- self.cell_volumes = compute_tri_areas(self.ei_dot_ej)
+ self.cell_volumes = compute_tri_areas(self.ei_dot_ej)
self._ce_ratios = None
self._interior_edge_lengths = None
@@ -199,7 +222,6 @@ def update_values(self):
self._surface_areas = None
self._signed_cell_areas = None
self._cell_centroids = None
- return
def remove_cells(self, remove_array):
"""Remove cells and take care of all the dependent data structures. The input
@@ -207,22 +229,22 @@ def remove_cells(self, remove_array):
"""
remove_array = numpy.asarray(remove_array)
if len(remove_array) == 0:
- return
+ return 0
if remove_array.dtype == int:
- keep = numpy.ones(len(self.cells["nodes"]), dtype=bool)
+ keep = numpy.ones(len(self.cells["points"]), dtype=bool)
keep[remove_array] = False
else:
assert remove_array.dtype == bool
keep = ~remove_array
- assert len(keep) == len(self.cells["nodes"]), "Wrong length of index array."
+ assert len(keep) == len(self.cells["points"]), "Wrong length of index array."
if numpy.all(keep):
- return
+ return 0
self.cell_volumes = self.cell_volumes[keep]
- self.cells["nodes"] = self.cells["nodes"][keep]
+ self.cells["points"] = self.cells["points"][keep]
self.idx_hierarchy = self.idx_hierarchy[..., keep]
if self._ce_ratios is not None:
@@ -253,11 +275,11 @@ def remove_cells(self, remove_array):
self._cvc_cell_mask = None
self._surface_areas = None
self._signed_cell_areas = None
- self._is_boundary_node = None
+ self._is_boundary_point = None
# handle edges; this is a bit messy
if "edges" in self.cells:
- num_edges_old = len(self.edges["nodes"])
+ num_edges_old = len(self.edges["points"])
adjacent_edges, counts = numpy.unique(
numpy.concatenate(self.cells["edges"][~keep]), return_counts=True
)
@@ -272,9 +294,9 @@ def remove_cells(self, remove_array):
# Now actually remove the edges. This includes a reindexing.
remaining_edges = numpy.ones(len(self.is_boundary_edge_gid), dtype=bool)
remaining_edges[adjacent_edges[is_edge_removed]] = False
- # make sure there is only edges["nodes"], not edges["cells"] etc.
+ # make sure there is only edges["points"], not edges["cells"] etc.
assert len(self.edges) == 1
- self.edges["nodes"] = self.edges["nodes"][remaining_edges]
+ self.edges["points"] = self.edges["points"][remaining_edges]
self.is_boundary_edge_gid = self.is_boundary_edge_gid[remaining_edges]
self.cells["edges"] = self.cells["edges"][keep]
@@ -296,7 +318,7 @@ def ce_ratios_per_interior_edge(self):
if "edges" not in self.cells:
self.create_edges()
- n = self.edges["nodes"].shape[0]
+ n = self.edges["points"].shape[0]
ce_ratios = numpy.bincount(
self.cells["edges"].reshape(-1),
self.ce_ratios.T.reshape(-1),
@@ -336,9 +358,9 @@ def get_control_volumes(self, cell_mask=None):
vals = numpy.array([v[1] + v[2], v[2] + v[0], v[0] + v[1]])
# sum all the vals into self._control_volumes at ids
self._control_volumes = numpy.bincount(
- self.cells["nodes"][~cell_mask].T.reshape(-1),
+ self.cells["points"][~cell_mask].T.reshape(-1),
weights=vals.reshape(-1),
- minlength=len(self.node_coords),
+ minlength=len(self.points),
)
self._cv_cell_mask = cell_mask
return self._control_volumes
@@ -377,14 +399,14 @@ def get_control_volume_centroids(self, cell_mask=None):
_, v = self._compute_integral_x()
v = v[:, :, ~cell_mask, :]
- # Again, make use of the fact that edge k is opposite of node k in every
+ # Again, make use of the fact that edge k is opposite of point k in every
# cell. Adding the arrays first makes the work for bincount lighter.
- ids = self.cells["nodes"][~cell_mask].T
+ ids = self.cells["points"][~cell_mask].T
vals = numpy.array(
[v[1, 1] + v[0, 2], v[1, 2] + v[0, 0], v[1, 0] + v[0, 1]]
)
# add it all up
- n = len(self.node_coords)
+ n = len(self.points)
self._cv_centroids = numpy.array(
[
numpy.bincount(
@@ -412,13 +434,13 @@ def signed_cell_areas(self):
"""Signed area of a triangle in 2D."""
# http://mathworld.wolfram.com/TriangleArea.html
assert (
- self.node_coords.shape[1] == 2
+ self.points.shape[1] == 2
), "Signed areas only make sense for triangles in 2D."
if self._signed_cell_areas is None:
# One could make p contiguous by adding a copy(), but that's not
# really worth it here.
- p = self.node_coords[self.cells["nodes"]].T
+ p = self.points[self.cells["points"]].T
#
self._signed_cell_areas = (
+p[0][2] * (p[1][0] - p[1][1])
@@ -435,27 +457,27 @@ def mark_boundary(self):
assert self.is_boundary_edge is not None
- self._is_boundary_node = numpy.zeros(len(self.node_coords), dtype=bool)
- self._is_boundary_node[self.idx_hierarchy[..., self.is_boundary_edge]] = True
+ self._is_boundary_point = numpy.zeros(len(self.points), dtype=bool)
+ self._is_boundary_point[self.idx_hierarchy[..., self.is_boundary_edge]] = True
- self._is_interior_node = self.node_is_used & ~self.is_boundary_node
+ self._is_interior_point = self.point_is_used & ~self.is_boundary_point
self._is_boundary_facet = self.is_boundary_edge
return
@property
- def is_boundary_node(self):
+ def is_boundary_point(self):
""""""
- if self._is_boundary_node is None:
+ if self._is_boundary_point is None:
self.mark_boundary()
- return self._is_boundary_node
+ return self._is_boundary_point
@property
- def is_interior_node(self):
+ def is_interior_point(self):
""""""
- if self._is_interior_node is None:
+ if self._is_interior_point is None:
self.mark_boundary()
- return self._is_interior_node
+ return self._is_interior_point
@property
def is_boundary_facet(self):
@@ -465,7 +487,7 @@ def is_boundary_facet(self):
return self._is_boundary_facet
def create_edges(self):
- """Set up edge->node and edge->cell relations."""
+ """Set up edge->point and edge->cell relations."""
# Reshape into individual edges.
# Sort the columns to make it possible for `unique()` to identify
# individual edges.
@@ -481,7 +503,7 @@ def create_edges(self):
self.is_boundary_edge_gid = cts == 1
- self.edges = {"nodes": a_unique}
+ self.edges = {"points": a_unique}
# cell->edges relationship
self.cells["edges"] = inv.reshape(3, -1).T
@@ -510,7 +532,7 @@ def _compute_edges_cells(self):
if self.edges is None:
self.create_edges()
- num_edges = len(self.edges["nodes"])
+ num_edges = len(self.edges["points"])
# count = numpy.bincount(self.cells["edges"].reshape(-1), minlength=num_edges)
@@ -554,7 +576,6 @@ def edge_gid_to_edge_list(self):
@property
def face_partitions(self):
- """"""
# face = edge for triangles.
# The partition is simply along the center of the edge.
edge_lengths = self.edge_lengths
@@ -562,7 +583,6 @@ def face_partitions(self):
@property
def cell_partitions(self):
- """"""
if self._cell_partitions is None:
# Compute the control volumes. Note that
#
@@ -574,11 +594,10 @@ def cell_partitions(self):
@property
def cell_circumcenters(self):
- """"""
if self._cell_circumcenters is None:
- node_cells = self.cells["nodes"].T
+ point_cells = self.cells["points"].T
self._cell_circumcenters = compute_triangle_circumcenters(
- self.node_coords[node_cells], self.ei_dot_ei, self.ei_dot_ej
+ self.points[point_cells], self.ei_dot_ei, self.ei_dot_ej
)
return self._cell_circumcenters
@@ -587,13 +606,13 @@ def cell_centroids(self):
"""The centroids (barycenters, midpoints of the circumcircles) of all triangles."""
if self._cell_centroids is None:
self._cell_centroids = (
- numpy.sum(self.node_coords[self.cells["nodes"]], axis=1) / 3.0
+ numpy.sum(self.points[self.cells["points"]], axis=1) / 3.0
)
return self._cell_centroids
@property
def cell_barycenters(self):
- """See cell_centroids()."""
+ """See cell_centroids."""
return self.cell_centroids
@property
@@ -602,7 +621,7 @@ def cell_incenters(self):
# https://en.wikipedia.org/wiki/Incenter#Barycentric_coordinates
abc = numpy.sqrt(self.ei_dot_ei)
abc /= numpy.sum(abc, axis=0)
- return numpy.einsum("ij,jik->jk", abc, self.node_coords[self.cells["nodes"]])
+ return numpy.einsum("ij,jik->jk", abc, self.points[self.cells["points"]])
@property
def cell_inradius(self):
@@ -658,7 +677,7 @@ def _compute_integral_x(self):
#
# \\int_V x,
#
- # over all atomic "triangles", i.e., areas cornered by a node, an edge midpoint,
+ # over all atomic "triangles", i.e., areas cornered by a point, an edge midpoint,
# and a circumcenter.
# The integral of any linear function over a triangle is the average of the
@@ -666,9 +685,9 @@ def _compute_integral_x(self):
# triangle.
right_triangle_vols = self.cell_partitions
- node_edges = self.idx_hierarchy
+ point_edges = self.idx_hierarchy
- corner = self.node_coords[node_edges]
+ corner = self.points[point_edges]
edge_midpoints = 0.5 * (corner[0] + corner[1])
cc = self.cell_circumcenters
@@ -676,7 +695,7 @@ def _compute_integral_x(self):
contribs = right_triangle_vols[None, :, :, None] * average
- return node_edges, contribs
+ return point_edges, contribs
# def _compute_surface_areas(self, cell_ids):
# # For each edge, one half of the the edge goes to each of the end points. Used
@@ -684,9 +703,9 @@ def _compute_integral_x(self):
# # conditions if in the interior.
# #
# # Each of the three edges may contribute to the surface areas of all three
- # # vertices. Here, only the two adjacent nodes receive a contribution, but other
- # # approaches, may contribute to all three nodes.
- # cn = self.cells["nodes"][cell_ids]
+ # # vertices. Here, only the two adjacent points receive a contribution, but other
+ # # approaches, may contribute to all three points.
+ # cn = self.cells["points"][cell_ids]
# ids = numpy.stack([cn, cn, cn], axis=1)
# half_el = 0.5 * self.edge_lengths[..., cell_ids]
@@ -704,7 +723,7 @@ def _compute_integral_x(self):
# def compute_gradient(self, u):
# '''Computes an approximation to the gradient :math:`\\nabla u` of a
- # given scalar valued function :math:`u`, defined in the node points.
+ # given scalar valued function :math:`u`, defined in the point points.
# This is taken from :cite:`NME2187`,
#
# Discrete gradient method in solid mechanics,
@@ -713,39 +732,39 @@ def _compute_integral_x(self):
# https://doi.org/10.1002/nme.2187.
# '''
# if self.cell_circumcenters is None:
- # X = self.node_coords[self.cells['nodes']]
+ # X = self.points[self.cells['points']]
# self.cell_circumcenters = self.compute_triangle_circumcenters(X)
#
# if 'cells' not in self.edges:
# self.edges['cells'] = self.compute_edge_cells()
#
# # This only works for flat meshes.
- # assert (abs(self.node_coords[:, 2]) < 1.0e-10).all()
- # node_coords2d = self.node_coords[:, :2]
+ # assert (abs(self.points[:, 2]) < 1.0e-10).all()
+ # points2d = self.points[:, :2]
# cell_circumcenters2d = self.cell_circumcenters[:, :2]
#
- # num_nodes = len(node_coords2d)
- # assert len(u) == num_nodes
+ # num_points = len(points2d)
+ # assert len(u) == num_points
#
- # gradient = numpy.zeros((num_nodes, 2), dtype=u.dtype)
+ # gradient = numpy.zeros((num_points, 2), dtype=u.dtype)
#
- # # Create an empty 2x2 matrix for the boundary nodes to hold the
+ # # Create an empty 2x2 matrix for the boundary points to hold the
# # edge correction ((17) in [1]).
# boundary_matrices = {}
- # for node in self.get_vertices('boundary'):
- # boundary_matrices[node] = numpy.zeros((2, 2))
+ # for point in self.get_vertices('boundary'):
+ # boundary_matrices[point] = numpy.zeros((2, 2))
#
# for edge_gid, edge in enumerate(self.edges['cells']):
# # Compute edge length.
- # node0 = self.edges['nodes'][edge_gid][0]
- # node1 = self.edges['nodes'][edge_gid][1]
+ # point0 = self.edges['points'][edge_gid][0]
+ # point1 = self.edges['points'][edge_gid][1]
#
# # Compute coedge length.
# if len(self.edges['cells'][edge_gid]) == 1:
# # Boundary edge.
# edge_midpoint = 0.5 * (
- # node_coords2d[node0] +
- # node_coords2d[node1]
+ # points2d[point0] +
+ # points2d[point1]
# )
# cell0 = self.edges['cells'][edge_gid][0]
# coedge_midpoint = 0.5 * (
@@ -767,23 +786,23 @@ def _compute_integral_x(self):
#
# # Compute the coefficient r for both contributions
# coeffs = self.ce_ratios[edge_gid] / \
- # self.control_volumes[self.edges['nodes'][edge_gid]]
+ # self.control_volumes[self.edges['points'][edge_gid]]
#
# # Compute R*_{IJ} ((11) in [1]).
- # r0 = (coedge_midpoint - node_coords2d[node0]) * coeffs[0]
- # r1 = (coedge_midpoint - node_coords2d[node1]) * coeffs[1]
+ # r0 = (coedge_midpoint - points2d[point0]) * coeffs[0]
+ # r1 = (coedge_midpoint - points2d[point1]) * coeffs[1]
#
- # diff = u[node1] - u[node0]
+ # diff = u[point1] - u[point0]
#
- # gradient[node0] += r0 * diff
- # gradient[node1] -= r1 * diff
+ # gradient[point0] += r0 * diff
+ # gradient[point1] -= r1 * diff
#
# # Store the boundary correction matrices.
- # edge_coords = node_coords2d[node1] - node_coords2d[node0]
- # if node0 in boundary_matrices:
- # boundary_matrices[node0] += numpy.outer(r0, edge_coords)
- # if node1 in boundary_matrices:
- # boundary_matrices[node1] += numpy.outer(r1, -edge_coords)
+ # edge_coords = points2d[point1] - points2d[point0]
+ # if point0 in boundary_matrices:
+ # boundary_matrices[point0] += numpy.outer(r0, edge_coords)
+ # if point1 in boundary_matrices:
+ # boundary_matrices[point1] += numpy.outer(r1, -edge_coords)
#
# # Apply corrections to the gradients on the boundary.
# for k, value in boundary_matrices.items():
@@ -865,7 +884,7 @@ def plot(
comesh_color=(0.8, 0.8, 0.8),
show_axes=True,
cell_quality_coloring=None,
- show_node_numbers=False,
+ show_point_numbers=False,
show_cell_numbers=False,
cell_mask=None,
show_edge_numbers=False,
@@ -881,10 +900,10 @@ def plot(
if not show_axes:
ax.set_axis_off()
- xmin = numpy.amin(self.node_coords[:, 0])
- xmax = numpy.amax(self.node_coords[:, 0])
- ymin = numpy.amin(self.node_coords[:, 1])
- ymax = numpy.amax(self.node_coords[:, 1])
+ xmin = numpy.amin(self.points[:, 0])
+ xmax = numpy.amax(self.points[:, 0])
+ ymin = numpy.amin(self.points[:, 1])
+ ymax = numpy.amax(self.points[:, 1])
width = xmax - xmin
xmin -= 0.1 * width
@@ -897,14 +916,14 @@ def plot(
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
- # for k, x in enumerate(self.node_coords):
- # if self.is_boundary_node[k]:
+ # for k, x in enumerate(self.points):
+ # if self.is_boundary_point[k]:
# plt.plot(x[0], x[1], "g.")
# else:
# plt.plot(x[0], x[1], "r.")
- if show_node_numbers:
- for i, x in enumerate(self.node_coords):
+ if show_point_numbers:
+ for i, x in enumerate(self.points):
plt.text(
x[0],
x[1],
@@ -929,9 +948,9 @@ def plot(
if cell_quality_coloring:
cmap, cmin, cmax, show_colorbar = cell_quality_coloring
plt.tripcolor(
- self.node_coords[:, 0],
- self.node_coords[:, 1],
- self.cells["nodes"],
+ self.points[:, 0],
+ self.points[:, 1],
+ self.cells["points"],
self.q_radius_ratio,
shading="flat",
cmap=cmap,
@@ -945,7 +964,7 @@ def plot(
self.create_edges()
# Get edges, cut off z-component.
- e = self.node_coords[self.edges["nodes"]][:, :, :2]
+ e = self.points[self.edges["points"]][:, :, :2]
if nondelaunay_edge_color is None:
line_segments0 = LineCollection(e, color=mesh_color)
@@ -955,7 +974,7 @@ def plot(
ce_ratios = self.ce_ratios_per_interior_edge
pos = ce_ratios >= 0
- is_pos = numpy.zeros(len(self.edges["nodes"]), dtype=bool)
+ is_pos = numpy.zeros(len(self.edges["points"]), dtype=bool)
is_pos[self._edge_to_edge_gid[2][pos]] = True
# Mark Delaunay-conforming boundary edges
@@ -973,8 +992,8 @@ def plot(
cc = self.cell_circumcenters
edge_midpoints = 0.5 * (
- self.node_coords[self.edges["nodes"][:, 0]]
- + self.node_coords[self.edges["nodes"][:, 1]]
+ self.points[self.edges["points"][:, 0]]
+ + self.points[self.edges["points"][:, 1]]
)
# Plot connection of the circumcenter to the midpoint of all three
@@ -995,9 +1014,7 @@ def plot(
ax.add_collection(line_segments)
if boundary_edge_color:
- e = self.node_coords[self.edges["nodes"][self.is_boundary_edge_gid]][
- :, :, :2
- ]
+ e = self.points[self.edges["points"][self.is_boundary_edge_gid]][:, :, :2]
line_segments1 = LineCollection(e, color=boundary_edge_color)
ax.add_collection(line_segments1)
@@ -1031,13 +1048,13 @@ def show_vertex(self, *args, **kwargs):
plt.close()
return
- def plot_vertex(self, node_id, show_ce_ratio=True):
- """Plot the vicinity of a node and its covolume/edgelength ratio.
+ def plot_vertex(self, point_id, show_ce_ratio=True):
+ """Plot the vicinity of a point and its covolume/edgelength ratio.
- :param node_id: Node ID of the node to be shown.
- :type node_id: int
+ :param point_id: Node ID of the point to be shown.
+ :type point_id: int
- :param show_ce_ratio: If true, shows the ce_ratio of the node, too.
+ :param show_ce_ratio: If true, shows the ce_ratio of the point, too.
:type show_ce_ratio: bool, optional
"""
# Importing matplotlib takes a while, so don't do that at the header.
@@ -1051,30 +1068,30 @@ def plot_vertex(self, node_id, show_ce_ratio=True):
self.create_edges()
# Find the edges that contain the vertex
- edge_gids = numpy.where((self.edges["nodes"] == node_id).any(axis=1))[0]
+ edge_gids = numpy.where((self.edges["points"] == point_id).any(axis=1))[0]
# ... and plot them
- for node_ids in self.edges["nodes"][edge_gids]:
- x = self.node_coords[node_ids]
+ for point_ids in self.edges["points"][edge_gids]:
+ x = self.points[point_ids]
ax.plot(x[:, 0], x[:, 1], "k")
# Highlight ce_ratios.
if show_ce_ratio:
if self.cell_circumcenters is None:
- X = self.node_coords[self.cells["nodes"]]
+ X = self.points[self.cells["points"]]
self.cell_circumcenters = self.compute_triangle_circumcenters(
X, self.ei_dot_ei, self.ei_dot_ej
)
# Find the cells that contain the vertex
- cell_ids = numpy.where((self.cells["nodes"] == node_id).any(axis=1))[0]
+ cell_ids = numpy.where((self.cells["points"] == point_id).any(axis=1))[0]
for cell_id in cell_ids:
for edge_gid in self.cells["edges"][cell_id]:
- if node_id not in self.edges["nodes"][edge_gid]:
+ if point_id not in self.edges["points"][edge_gid]:
continue
- node_ids = self.edges["nodes"][edge_gid]
+ point_ids = self.edges["points"][edge_gid]
edge_midpoint = 0.5 * (
- self.node_coords[node_ids[0]] + self.node_coords[node_ids[1]]
+ self.points[point_ids[0]] + self.points[point_ids[1]]
)
p = numpy.stack(
[self.cell_circumcenters[cell_id], edge_midpoint], axis=1
@@ -1083,7 +1100,7 @@ def plot_vertex(self, node_id, show_ce_ratio=True):
[
self.cell_circumcenters[cell_id],
edge_midpoint,
- self.node_coords[node_id],
+ self.points[point_id],
]
)
ax.fill(q[0], q[1], color="0.5")
@@ -1180,27 +1197,27 @@ def flip_interior_edges(self, is_flip_interior_edge):
#
verts = numpy.array(
[
- self.cells["nodes"][adj_cells[0], lids[0]],
- self.cells["nodes"][adj_cells[1], lids[1]],
- self.cells["nodes"][adj_cells[0], (lids[0] + 1) % 3],
- self.cells["nodes"][adj_cells[0], (lids[0] + 2) % 3],
+ self.cells["points"][adj_cells[0], lids[0]],
+ self.cells["points"][adj_cells[1], lids[1]],
+ self.cells["points"][adj_cells[0], (lids[0] + 1) % 3],
+ self.cells["points"][adj_cells[0], (lids[0] + 2) % 3],
]
)
- self.edges["nodes"][edge_gids] = numpy.sort(verts[[0, 1]].T, axis=1)
+ self.edges["points"][edge_gids] = numpy.sort(verts[[0, 1]].T, axis=1)
# No need to touch self.is_boundary_edge,
# self.is_boundary_edge_gid; we're only flipping interior edges.
- # Do the neighboring cells have equal orientation (both node sets
+ # Do the neighboring cells have equal orientation (both point sets
# clockwise/counterclockwise?
equal_orientation = (
- self.cells["nodes"][adj_cells[0], (lids[0] + 1) % 3]
- == self.cells["nodes"][adj_cells[1], (lids[1] + 2) % 3]
+ self.cells["points"][adj_cells[0], (lids[0] + 1) % 3]
+ == self.cells["points"][adj_cells[1], (lids[1] + 2) % 3]
)
# Set new cells
- self.cells["nodes"][adj_cells[0]] = verts[[0, 1, 2]].T
- self.cells["nodes"][adj_cells[1]] = verts[[0, 1, 3]].T
+ self.cells["points"][adj_cells[0]] = verts[[0, 1, 2]].T
+ self.cells["points"][adj_cells[1]] = verts[[0, 1, 3]].T
# Set up new cells->edges relationships.
previous_edges = self.cells["edges"][adj_cells].copy()
@@ -1273,13 +1290,13 @@ def flip_interior_edges(self, is_flip_interior_edge):
def _update_cell_values(self, cell_ids, interior_edge_ids):
"""Updates all sorts of cell information for the given cell IDs."""
# update idx_hierarchy
- nds = self.cells["nodes"][cell_ids].T
+ nds = self.cells["points"][cell_ids].T
self.idx_hierarchy[..., cell_ids] = nds[self.local_idx]
# update self.half_edge_coords
self.half_edge_coords[:, cell_ids, :] = numpy.moveaxis(
- self.node_coords[self.idx_hierarchy[1, ..., cell_ids]]
- - self.node_coords[self.idx_hierarchy[0, ..., cell_ids]],
+ self.points[self.idx_hierarchy[1, ..., cell_ids]]
+ - self.points[self.idx_hierarchy[0, ..., cell_ids]],
0,
1,
)
diff --git a/setup.cfg b/setup.cfg
index 0aff65f..b51fba7 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,6 +1,6 @@
[metadata]
name = meshplex
-version = 0.13.6
+version = 0.14.0
author = Nico Schlömer
author_email = nico.schloemer@gmail.com
description = Fast tools for simplex meshes
diff --git a/test/test_curl.py b/test/test_curl.py
index 34eb769..04d620b 100644
--- a/test/test_curl.py
+++ b/test/test_curl.py
@@ -10,9 +10,7 @@
def _run(mesh):
# Create circular vector field 0.5 * (y, -x, 0)
# which has curl (0, 0, 1).
- A = numpy.array(
- [[-0.5 * coord[1], 0.5 * coord[0], 0.0] for coord in mesh.node_coords]
- )
+ A = numpy.array([[-0.5 * coord[1], 0.5 * coord[0], 0.0] for coord in mesh.points])
# Compute the curl numerically.
B = mesh.compute_curl(A)
diff --git a/test/test_gradient.py b/test/test_gradient.py
index fda2e42..626bd61 100644
--- a/test/test_gradient.py
+++ b/test/test_gradient.py
@@ -11,13 +11,13 @@
# return
#
# def _run_test(self, mesh):
-# num_nodes = len(mesh.node_coords)
+# num_nodes = len(mesh.points)
# # Create function 2*x + 3*y.
# a_x = 7.0
# a_y = 3.0
# a0 = 1.0
-# u = a_x * mesh.node_coords[:, 0] + \
-# a_y * mesh.node_coords[:, 1] + \
+# u = a_x * mesh.points[:, 0] + \
+# a_y * mesh.points[:, 1] + \
# a0 * numpy.ones(num_nodes)
# # Get the gradient analytically.
# sol = numpy.empty((num_nodes, 2))
diff --git a/test/test_mesh_io.py b/test/test_mesh_io.py
index 5413d6b..d2b8bab 100644
--- a/test/test_mesh_io.py
+++ b/test/test_mesh_io.py
@@ -20,8 +20,8 @@ def test_io_2d():
mesh2 = meshplex.read(fname)
- for k in range(len(mesh.cells["nodes"])):
- assert tuple(mesh.cells["nodes"][k]) == tuple(mesh2.cells["nodes"][k])
+ for k in range(len(mesh.cells["points"])):
+ assert tuple(mesh.cells["points"][k]) == tuple(mesh2.cells["points"][k])
# def test_io_3d(self):
@@ -42,8 +42,8 @@ def test_io_2d():
# mesh2, _, _, _ = meshplex.read('test.vtu')
-# for k in range(len(mesh.cells['nodes'])):
+# for k in range(len(mesh.cells['points'])):
# self.assertEqual(
-# tuple(mesh.cells['nodes'][k]),
-# tuple(mesh2.cells['nodes'][k])
+# tuple(mesh.cells['points'][k]),
+# tuple(mesh2.cells['points'][k])
# )
diff --git a/test/test_mesh_tetra.py b/test/test_mesh_tetra.py
index e030430..b33d822 100644
--- a/test/test_mesh_tetra.py
+++ b/test/test_mesh_tetra.py
@@ -27,7 +27,7 @@ def test_regular_tet0(a):
cells = numpy.array([[0, 1, 2, 3]])
mesh = meshplex.MeshTetra(points, cells.copy())
- assert numpy.all(mesh.cells["nodes"] == cells)
+ assert numpy.all(mesh.cells["points"] == cells)
mesh.show()
mesh.show_edge(0)
@@ -163,7 +163,7 @@ def test_unit_tetrahedron_geometric(a):
mesh = meshplex.MeshTetra(points, cells)
- assert all((mesh.cells["nodes"] == cells).flat)
+ assert all((mesh.cells["points"] == cells).flat)
assert near_equal(mesh.cell_circumcenters, [a / 2.0, a / 2.0, a / 2.0], tol)
@@ -218,7 +218,7 @@ def test_regular_tet1_geometric_order():
mesh = meshplex.MeshTetra(points, cells)
- assert all((mesh.cells["nodes"] == [0, 1, 2, 3]).flat)
+ assert all((mesh.cells["points"] == [0, 1, 2, 3]).flat)
assert near_equal(mesh.cell_circumcenters, [a / 2.0, a / 2.0, a / 2.0], tol)
@@ -364,7 +364,7 @@ def test_cubesmall():
def test_arrow3d():
- nodes = numpy.array(
+ points = numpy.array(
[
[+0.0, +0.0, +0.0],
[+2.0, -1.0, +0.0],
@@ -375,7 +375,7 @@ def test_arrow3d():
]
)
cellsNodes = numpy.array([[1, 2, 4, 5], [2, 3, 4, 5], [0, 1, 4, 5], [0, 3, 4, 5]])
- mesh = meshplex.MeshTetra(nodes, cellsNodes)
+ mesh = meshplex.MeshTetra(points, cellsNodes)
run(
mesh,
@@ -413,11 +413,11 @@ def test_tetrahedron():
# # the order of 1e-16. The ce_ratios differ by up to 1e-7.
# if False:
# print(mesh.cells.keys())
-# pts = mesh.node_coords.copy()
+# pts = mesh.points.copy()
# pts += 1.0e-16 * numpy.random.rand(pts.shape[0], pts.shape[1])
-# mesh2 = meshplex.MeshTetra(pts, mesh.cells['nodes'])
+# mesh2 = meshplex.MeshTetra(pts, mesh.cells['points'])
# #
-# diff_coords = mesh.node_coords - mesh2.node_coords
+# diff_coords = mesh.points - mesh2.points
# max_diff_coords = max(diff_coords.flatten())
# print('||coords_1 - coords_2||_inf = %e' % max_diff_coords)
# diff_ce_ratios = mesh.ce_ratios_per_edge - mesh2.ce_ratios
@@ -445,7 +445,7 @@ def test_tetrahedron():
def test_toy_geometric():
mesh = meshplex.read(this_dir / "meshes" / "toy.vtk")
- mesh = meshplex.MeshTetra(mesh.node_coords, mesh.cells["nodes"])
+ mesh = meshplex.MeshTetra(mesh.points, mesh.cells["points"])
run(
mesh,
@@ -466,7 +466,7 @@ def test_toy_geometric():
def test_signed_volume():
mesh = meshplex.read(this_dir / "meshes" / "toy.vtk")
- vols = meshplex.get_signed_simplex_volumes(mesh.cells["nodes"], mesh.node_coords)
+ vols = meshplex.get_signed_simplex_volumes(mesh.cells["points"], mesh.points)
assert numpy.all(abs(abs(vols) - mesh.cell_volumes) < 1.0e-12 * mesh.cell_volumes)
@@ -476,7 +476,7 @@ def test_show_cell():
# filename = download_mesh("toy.vtk", "f48abda972822bab224b91a74d695573")
# mesh = meshplex.read(filename)
- node_coords = (
+ points = (
numpy.array(
[
[1.0, 0.0, -1.0 / numpy.sqrt(8)],
@@ -488,7 +488,7 @@ def test_show_cell():
/ numpy.sqrt(3.0)
)
- # node_coords = numpy.array(
+ # points = numpy.array(
# [
# [1.0, 0.0, -1.0 / numpy.sqrt(8)],
# [-0.5, +numpy.sqrt(3.0) / 2.0, -1.0 / numpy.sqrt(8)],
@@ -496,12 +496,12 @@ def test_show_cell():
# [0.0, 0.5, 1.0],
# ]
# ) / numpy.sqrt(3.0)
- # node_coords = numpy.array(
+ # points = numpy.array(
# [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
# )
cells = [[0, 1, 2, 3]]
- mesh = meshplex.MeshTetra(node_coords, cells)
+ mesh = meshplex.MeshTetra(points, cells)
mesh.show_cell(
0,
barycenter_rgba=(1, 0, 0, 1.0),
diff --git a/test/test_mesh_tri.py b/test/test_mesh_tri.py
index 9bb6c71..636bff5 100644
--- a/test/test_mesh_tri.py
+++ b/test/test_mesh_tri.py
@@ -95,17 +95,19 @@ def test_regular_tri_additional_points():
mesh.mark_boundary()
- assert numpy.array_equal(mesh.node_is_used, [False, True, True, True, False])
- assert numpy.array_equal(mesh.is_boundary_node, [False, True, True, True, False])
- assert numpy.array_equal(mesh.is_interior_node, [False, False, False, False, False])
+ assert numpy.array_equal(mesh.point_is_used, [False, True, True, True, False])
+ assert numpy.array_equal(mesh.is_boundary_point, [False, True, True, True, False])
+ assert numpy.array_equal(
+ mesh.is_interior_point, [False, False, False, False, False]
+ )
tol = 1.0e-14
- assert numpy.array_equal(mesh.cells["nodes"], [[1, 2, 3]])
+ assert numpy.array_equal(mesh.cells["points"], [[1, 2, 3]])
mesh.create_edges()
assert numpy.array_equal(mesh.cells["edges"], [[2, 1, 0]])
- assert numpy.array_equal(mesh.edges["nodes"], [[1, 2], [1, 3], [2, 3]])
+ assert numpy.array_equal(mesh.edges["points"], [[1, 2], [1, 3], [2, 3]])
# ce_ratios
assert near_equal(mesh.ce_ratios.T, [0.0, 0.5, 0.5], tol)
@@ -138,7 +140,7 @@ def test_regular_tri_order():
cells = numpy.array([[0, 1, 2]])
mesh = meshplex.MeshTri(points, cells)
- assert all((mesh.cells["nodes"] == [0, 1, 2]).flat)
+ assert all((mesh.cells["points"] == [0, 1, 2]).flat)
tol = 1.0e-14
@@ -506,13 +508,13 @@ def test_signed_area2():
ref = 0.5
assert abs(mesh.signed_cell_areas[0] - ref) < 1.0e-10 * abs(ref)
- mesh.node_coords = numpy.array([[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
+ mesh.points = numpy.array([[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
mesh.update_values()
ref = -0.5
assert abs(mesh.signed_cell_areas[0] - ref) < 1.0e-10 * abs(ref)
-def test_update_node_coordinates():
+def test_update_point_coordinates():
mesh = meshio.read(this_dir / "meshes" / "pacman.vtk")
assert numpy.all(numpy.abs(mesh.points[:, 2]) < 1.0e-15)
@@ -522,7 +524,7 @@ def test_update_node_coordinates():
X2 = mesh.points + 1.0e-2 * numpy.random.rand(*mesh.points.shape)
mesh2 = meshplex.MeshTri(X2, mesh.get_cells_type("triangle"))
- mesh1.node_coords = X2
+ mesh1.points = X2
mesh1.update_values()
tol = 1.0e-12
@@ -549,8 +551,8 @@ def test_flip_delaunay():
num_adj_cells, edge_id = mesh._edge_gid_to_edge_list[edge_gid]
assert cell_gid in mesh._edges_cells[num_adj_cells][edge_id]
- new_cells = mesh.cells["nodes"].copy()
- new_coords = mesh.node_coords.copy()
+ new_cells = mesh.cells["points"].copy()
+ new_coords = mesh.points.copy()
# Assert that some key values are updated properly
mesh2 = meshplex.MeshTri(new_coords, new_cells)
@@ -570,13 +572,13 @@ def test_flip_delaunay_near_boundary():
mesh.create_edges()
assert mesh.num_delaunay_violations() == 1
- assert numpy.array_equal(mesh.cells["nodes"], [[0, 1, 2], [0, 2, 3]])
+ assert numpy.array_equal(mesh.cells["points"], [[0, 1, 2], [0, 2, 3]])
assert numpy.array_equal(mesh.cells["edges"], [[3, 1, 0], [4, 2, 1]])
mesh.flip_until_delaunay()
assert mesh.num_delaunay_violations() == 0
- assert numpy.array_equal(mesh.cells["nodes"], [[1, 3, 2], [1, 3, 0]])
+ assert numpy.array_equal(mesh.cells["points"], [[1, 3, 2], [1, 3, 0]])
assert numpy.array_equal(mesh.cells["edges"], [[4, 3, 1], [2, 0, 1]])
@@ -600,7 +602,7 @@ def test_flip_same_edge_twice():
new_points = numpy.array(
[[0.0, +0.0, 0.0], [0.1, -0.5, 0.0], [0.2, +0.0, 0.0], [0.1, +0.5, 0.0]]
)
- mesh.node_coords = new_points
+ mesh.points = new_points
mesh.update_values()
assert mesh.num_delaunay_violations() == 1
@@ -624,12 +626,12 @@ def test_flip_two_edges():
assert mesh.num_delaunay_violations() == 0
assert numpy.array_equal(
- mesh.cells["nodes"], [[5, 2, 3], [0, 2, 1], [5, 2, 0], [3, 4, 5]]
+ mesh.cells["points"], [[5, 2, 3], [0, 2, 1], [5, 2, 0], [3, 4, 5]]
)
def test_flip_delaunay_near_boundary_preserve_boundary_count():
- # This test is to make sure meshplex preserves the boundary node count.
+ # This test is to make sure meshplex preserves the boundary point count.
points = numpy.array(
[
[+0.0, +0.0, 0.0],
@@ -647,13 +649,13 @@ def test_flip_delaunay_near_boundary_preserve_boundary_count():
assert mesh.num_delaunay_violations() == 1
mesh.mark_boundary()
- is_boundary_node_ref = [False, True, True, True, True, True]
- assert numpy.array_equal(mesh.is_boundary_node, is_boundary_node_ref)
+ is_boundary_point_ref = [False, True, True, True, True, True]
+ assert numpy.array_equal(mesh.is_boundary_point, is_boundary_point_ref)
mesh.flip_until_delaunay()
mesh.mark_boundary()
- assert numpy.array_equal(mesh.is_boundary_node, is_boundary_node_ref)
+ assert numpy.array_equal(mesh.is_boundary_point, is_boundary_point_ref)
def test_inradius():
@@ -802,7 +804,7 @@ def test_flat_boundary():
cv = numpy.zeros(X.shape[0])
for edges, ce_ratios in zip(mesh.idx_hierarchy.T, mesh.ce_ratios.T):
for i, ce in zip(edges, ce_ratios):
- ei = mesh.node_coords[i[1]] - mesh.node_coords[i[0]]
+ ei = mesh.points[i[1]] - mesh.points[i[0]]
cv[i] += 0.25 * ce * numpy.dot(ei, ei)
assert numpy.all(numpy.abs(cv - mesh.control_volumes) < 1.0e-12 * cv)
@@ -863,12 +865,12 @@ def _get_test_mesh():
)
def test_remove_cells(remove_idx, expected_num_cells, expected_num_edges):
mesh = _get_test_mesh()
- assert len(mesh.cells["nodes"]) == 6
- assert len(mesh.edges["nodes"]) == 13
+ assert len(mesh.cells["points"]) == 6
+ assert len(mesh.edges["points"]) == 13
# remove a corner cell
mesh.remove_cells(remove_idx)
- assert len(mesh.cells["nodes"]) == expected_num_cells
- assert len(mesh.edges["nodes"]) == expected_num_edges
+ assert len(mesh.cells["points"]) == expected_num_cells
+ assert len(mesh.edges["points"]) == expected_num_edges
if __name__ == "__main__":
diff --git a/test/test_readme.py b/test/test_readme.py
new file mode 100644
index 0000000..fef03e8
--- /dev/null
+++ b/test/test_readme.py
@@ -0,0 +1,19 @@
+import pathlib
+
+import exdown
+import pytest
+
+this_dir = pathlib.Path(__file__).resolve().parent
+
+
+@pytest.mark.parametrize(
+ "string,lineno",
+ exdown.extract(this_dir.parent / "README.md", syntax_filter="python"),
+)
+def test_readme(string, lineno):
+ try:
+ # https://stackoverflow.com/a/62851176/353337
+ exec(string, {"__MODULE__": "__main__"})
+ except Exception:
+ print(f"README.md (line {lineno}):\n```\n{string}```")
+ raise
diff --git a/test/test_subdomain.py b/test/test_subdomain.py
index ccd9905..e7c3f59 100644
--- a/test/test_subdomain.py
+++ b/test/test_subdomain.py
@@ -9,8 +9,8 @@ def test_get_edges():
mesh = meshplex.read(this_dir / "meshes" / "pacman.vtk")
mesh.create_edges()
edge_mask = mesh.get_edge_mask()
- edge_nodes = mesh.edges["nodes"][edge_mask]
- assert len(edge_nodes) == 1276
+ edge_points = mesh.edges["points"][edge_mask]
+ assert len(edge_points) == 1276
def test_mark_subdomain2d():
diff --git a/tox.ini b/tox.ini
index afbadad..0e88a1c 100644
--- a/tox.ini
+++ b/tox.ini
@@ -4,6 +4,7 @@ isolated_build = True
[testenv]
deps =
+ exdown
meshzoo
pytest
pytest-cov