From 305c25fbfe68a2a3e2217ae017c3eb2e53f30d0f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 5 Nov 2020 14:11:56 +0100 Subject: [PATCH 1/7] node_coords -> points --- README.md | 4 +- experimental/ce_smoothing.py | 12 ++--- meshplex/base.py | 14 ++--- meshplex/mesh_line.py | 13 ++--- meshplex/mesh_tetra.py | 45 ++++++++-------- meshplex/mesh_tri.py | 100 +++++++++++++++++------------------ test/test_curl.py | 4 +- test/test_gradient.py | 6 +-- test/test_mesh_tetra.py | 16 +++--- test/test_mesh_tri.py | 10 ++-- 10 files changed, 106 insertions(+), 118 deletions(-) diff --git a/README.md b/README.md index 187e7a1..02fe704 100644 --- a/README.md +++ b/README.md @@ -97,7 +97,7 @@ import numpy import meshplex # Generate tetrahedron -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)], @@ -108,7 +108,7 @@ node_coords = numpy.array( 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( 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..64a6ae9 100644 --- a/meshplex/base.py +++ b/meshplex/base.py @@ -156,17 +156,17 @@ def compute_triangle_circumcenters(X, ei_dot_ei, ei_dot_ej): class _base_mesh: def __init__(self, nodes, cells_nodes): - self.node_coords = nodes + self.points = nodes 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: cell_type = "triangle" @@ -262,9 +262,9 @@ 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 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..627f851 100644 --- a/meshplex/mesh_tetra.py +++ b/meshplex/mesh_tetra.py @@ -9,7 +9,7 @@ class MeshTetra(_base_mesh): """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 # many downstream applications, e.g., when constructing linear systems with the @@ -21,7 +21,7 @@ def __init__(self, node_coords, cells): cells = numpy.sort(cells, axis=1) cells = cells[cells[:, 0].argsort()] - super().__init__(node_coords, cells) + super().__init__(points, cells) # Assert that all vertices are used. # If there are vertices which do not appear in the cells list, this @@ -31,7 +31,7 @@ def __init__(self, node_coords, cells): # nodes = nodes[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( numpy.sum(~is_used) @@ -83,8 +83,7 @@ def __init__(self, node_coords, cells): # 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 @@ -112,7 +111,7 @@ def __init__(self, node_coords, cells): return def __repr__(self): - num_nodes = len(self.node_coords) + num_nodes = len(self.points) num_cells = len(self.cells["nodes"]) string = f"" return string @@ -125,7 +124,7 @@ 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 = numpy.zeros(len(self.points), dtype=bool) self.is_boundary_node[ self.faces["nodes"][self.is_boundary_facet_individual] ] = True @@ -222,7 +221,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["nodes"]], axis=1 ) return @@ -235,8 +234,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: @@ -386,7 +385,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["nodes"]], axis=1) / 4.0 ) return self._cell_centroids @@ -408,9 +407,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["nodes"]]) @property def cell_inradius(self): @@ -421,7 +418,7 @@ 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 + 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 # @@ -538,7 +535,7 @@ def control_volumes(self): self._control_volumes = numpy.bincount( self.cells["nodes"].reshape(-1), vals.reshape(-1), - minlength=len(self.node_coords), + minlength=len(self.points), ) return self._control_volumes @@ -580,7 +577,7 @@ def plot(self): if self._circumcenters is None: self._compute_cell_circumcenters() - X = self.node_coords + X = self.points for cell_id in range(len(self.cells["nodes"])): cc = self._circumcenters[cell_id] # @@ -651,16 +648,16 @@ 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["nodes"][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["nodes"][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] # @@ -749,7 +746,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["nodes"][cell_id][ij]] renderer.AddActor(get_line_actor(x0, x1, line_width)) renderer.SetBackground(1.0, 1.0, 1.0) @@ -789,7 +786,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.node_face_cells[..., [cell_id]]] face_ccs = compute_triangle_circumcenters( x, self.ei_dot_ei, self.ei_dot_ej )[:, 0, :] @@ -798,14 +795,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.node_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..2ce8e4c 100644 --- a/meshplex/mesh_tri.py +++ b/meshplex/mesh_tri.py @@ -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. @@ -120,25 +119,25 @@ 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["nodes"][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_nodes = len(self.points) num_cells = len(self.cells["nodes"]) string = f"" return string # def update_node_coordinates(self, X): - # assert X.shape == self.node_coords.shape - # self.node_coords = 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 + # assert X.shape == self.points[self.is_interior_node].shape + # self.points[self.is_interior_node] = X # self.update_values() # return @@ -148,7 +147,7 @@ def euler_characteristic(self): if "edges" not in self.cells: self.create_edges() return ( - self.node_coords.shape[0] + self.points.shape[0] - self.edges["nodes"].shape[0] + self.cells["nodes"].shape[0] ) @@ -168,10 +167,9 @@ 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.points[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]] + self.points[self.idx_hierarchy[1]] - self.points[self.idx_hierarchy[0]] ) if self.ei_dot_ej is not None: @@ -207,7 +205,7 @@ 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) @@ -219,7 +217,7 @@ def remove_cells(self, remove_array): assert len(keep) == len(self.cells["nodes"]), "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] @@ -338,7 +336,7 @@ def get_control_volumes(self, cell_mask=None): self._control_volumes = numpy.bincount( self.cells["nodes"][~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 @@ -384,7 +382,7 @@ def get_control_volume_centroids(self, cell_mask=None): [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 +410,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["nodes"]].T # self._signed_cell_areas = ( +p[0][2] * (p[1][0] - p[1][1]) @@ -435,7 +433,7 @@ 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 = numpy.zeros(len(self.points), dtype=bool) self._is_boundary_node[self.idx_hierarchy[..., self.is_boundary_edge]] = True self._is_interior_node = self.node_is_used & ~self.is_boundary_node @@ -578,7 +576,7 @@ def cell_circumcenters(self): if self._cell_circumcenters is None: node_cells = self.cells["nodes"].T self._cell_circumcenters = compute_triangle_circumcenters( - self.node_coords[node_cells], self.ei_dot_ei, self.ei_dot_ej + self.points[node_cells], self.ei_dot_ei, self.ei_dot_ej ) return self._cell_circumcenters @@ -587,7 +585,7 @@ 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["nodes"]], axis=1) / 3.0 ) return self._cell_centroids @@ -602,7 +600,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["nodes"]]) @property def cell_inradius(self): @@ -668,7 +666,7 @@ def _compute_integral_x(self): node_edges = self.idx_hierarchy - corner = self.node_coords[node_edges] + corner = self.points[node_edges] edge_midpoints = 0.5 * (corner[0] + corner[1]) cc = self.cell_circumcenters @@ -713,18 +711,18 @@ 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['nodes']] # 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) + # num_nodes = len(points2d) # assert len(u) == num_nodes # # gradient = numpy.zeros((num_nodes, 2), dtype=u.dtype) @@ -744,8 +742,8 @@ def _compute_integral_x(self): # if len(self.edges['cells'][edge_gid]) == 1: # # Boundary edge. # edge_midpoint = 0.5 * ( - # node_coords2d[node0] + - # node_coords2d[node1] + # points2d[node0] + + # points2d[node1] # ) # cell0 = self.edges['cells'][edge_gid][0] # coedge_midpoint = 0.5 * ( @@ -770,8 +768,8 @@ def _compute_integral_x(self): # self.control_volumes[self.edges['nodes'][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[node0]) * coeffs[0] + # r1 = (coedge_midpoint - points2d[node1]) * coeffs[1] # # diff = u[node1] - u[node0] # @@ -779,7 +777,7 @@ def _compute_integral_x(self): # gradient[node1] -= r1 * diff # # # Store the boundary correction matrices. - # edge_coords = node_coords2d[node1] - node_coords2d[node0] + # edge_coords = points2d[node1] - points2d[node0] # if node0 in boundary_matrices: # boundary_matrices[node0] += numpy.outer(r0, edge_coords) # if node1 in boundary_matrices: @@ -881,10 +879,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 +895,14 @@ def plot( ax.set_xlim(xmin, xmax) ax.set_ylim(ymin, ymax) - # for k, x in enumerate(self.node_coords): + # for k, x in enumerate(self.points): # if self.is_boundary_node[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): + for i, x in enumerate(self.points): plt.text( x[0], x[1], @@ -929,8 +927,8 @@ 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.points[:, 0], + self.points[:, 1], self.cells["nodes"], self.q_radius_ratio, shading="flat", @@ -945,7 +943,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["nodes"]][:, :, :2] if nondelaunay_edge_color is None: line_segments0 = LineCollection(e, color=mesh_color) @@ -973,8 +971,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["nodes"][:, 0]] + + self.points[self.edges["nodes"][:, 1]] ) # Plot connection of the circumcenter to the midpoint of all three @@ -995,9 +993,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["nodes"][self.is_boundary_edge_gid]][:, :, :2] line_segments1 = LineCollection(e, color=boundary_edge_color) ax.add_collection(line_segments1) @@ -1054,13 +1050,13 @@ def plot_vertex(self, node_id, show_ce_ratio=True): edge_gids = numpy.where((self.edges["nodes"] == node_id).any(axis=1))[0] # ... and plot them for node_ids in self.edges["nodes"][edge_gids]: - x = self.node_coords[node_ids] + x = self.points[node_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["nodes"]] self.cell_circumcenters = self.compute_triangle_circumcenters( X, self.ei_dot_ei, self.ei_dot_ej ) @@ -1074,7 +1070,7 @@ def plot_vertex(self, node_id, show_ce_ratio=True): continue node_ids = self.edges["nodes"][edge_gid] edge_midpoint = 0.5 * ( - self.node_coords[node_ids[0]] + self.node_coords[node_ids[1]] + self.points[node_ids[0]] + self.points[node_ids[1]] ) p = numpy.stack( [self.cell_circumcenters[cell_id], edge_midpoint], axis=1 @@ -1083,7 +1079,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[node_id], ] ) ax.fill(q[0], q[1], color="0.5") @@ -1278,8 +1274,8 @@ def _update_cell_values(self, cell_ids, interior_edge_ids): # 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/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_tetra.py b/test/test_mesh_tetra.py index e030430..ba6dc03 100644 --- a/test/test_mesh_tetra.py +++ b/test/test_mesh_tetra.py @@ -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']) # # -# 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["nodes"]) 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["nodes"], 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..04356c1 100644 --- a/test/test_mesh_tri.py +++ b/test/test_mesh_tri.py @@ -506,7 +506,7 @@ 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) @@ -522,7 +522,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 @@ -550,7 +550,7 @@ def test_flip_delaunay(): 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_coords = mesh.points.copy() # Assert that some key values are updated properly mesh2 = meshplex.MeshTri(new_coords, new_cells) @@ -600,7 +600,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 @@ -802,7 +802,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) From bc74f373d4b224aa9e7235bb1e0c27f25175b037 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 5 Nov 2020 16:07:19 +0100 Subject: [PATCH 2/7] _base_mesh -> _BaseMesh --- meshplex/base.py | 156 +---------------------------------------- meshplex/helpers.py | 150 +++++++++++++++++++++++++++++++++++++++ meshplex/mesh_tetra.py | 4 +- meshplex/mesh_tri.py | 4 +- 4 files changed, 157 insertions(+), 157 deletions(-) diff --git a/meshplex/base.py b/meshplex/base.py index 64a6ae9..8b72a5b 100644 --- a/meshplex/base.py +++ b/meshplex/base.py @@ -1,162 +1,12 @@ 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.points = nodes +class _BaseMesh: + def __init__(self, points, cells_nodes): + self.points = points self._edge_lengths = None def write(self, filename, point_data=None, cell_data=None, field_data=None): 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_tetra.py b/meshplex/mesh_tetra.py index 627f851..c47555b 100644 --- a/meshplex/mesh_tetra.py +++ b/meshplex/mesh_tetra.py @@ -1,12 +1,12 @@ import numpy -from .base import _base_mesh, compute_tri_areas, compute_triangle_circumcenters +from .base import _BaseMesh, 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, points, cells): diff --git a/meshplex/mesh_tri.py b/meshplex/mesh_tri.py index 2ce8e4c..1756ef4 100644 --- a/meshplex/mesh_tri.py +++ b/meshplex/mesh_tri.py @@ -4,7 +4,7 @@ import numpy from .base import ( - _base_mesh, + _BaseMesh, compute_ce_ratios, compute_tri_areas, compute_triangle_circumcenters, @@ -14,7 +14,7 @@ __all__ = ["MeshTri"] -class MeshTri(_base_mesh): +class MeshTri(_BaseMesh): """Class for handling triangular meshes.""" def __init__(self, nodes, cells, sort_cells=False): From fbe0804014e1712bcd7c9c7f98b3082367744f02 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 5 Nov 2020 17:24:55 +0100 Subject: [PATCH 3/7] node -> points, fix point setter in MeshTri --- meshplex/base.py | 15 +- meshplex/mesh_tetra.py | 108 +++++++------ meshplex/mesh_tri.py | 325 +++++++++++++++++++++------------------- test/test_mesh_io.py | 10 +- test/test_mesh_tetra.py | 16 +- test/test_mesh_tri.py | 40 ++--- test/test_subdomain.py | 4 +- 7 files changed, 276 insertions(+), 242 deletions(-) diff --git a/meshplex/base.py b/meshplex/base.py index 8b72a5b..a57e885 100644 --- a/meshplex/base.py +++ b/meshplex/base.py @@ -5,8 +5,7 @@ class _BaseMesh: - def __init__(self, points, cells_nodes): - self.points = points + def __init__(self, points, cells_points): self._edge_lengths = None def write(self, filename, point_data=None, cell_data=None, field_data=None): @@ -18,18 +17,18 @@ def write(self, filename, point_data=None, cell_data=None, field_data=None): else: 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, @@ -59,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))) @@ -80,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) @@ -119,7 +118,7 @@ def _mark_vertices(self, subdomain): 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/mesh_tetra.py b/meshplex/mesh_tetra.py index c47555b..e311349 100644 --- a/meshplex/mesh_tetra.py +++ b/meshplex/mesh_tetra.py @@ -1,6 +1,7 @@ import numpy -from .base import _BaseMesh, compute_tri_areas, compute_triangle_circumcenters +from .base import _BaseMesh +from .helpers import compute_tri_areas, compute_triangle_circumcenters __all__ = ["MeshTetra"] @@ -11,7 +12,7 @@ class MeshTetra(_BaseMesh): 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,6 +22,7 @@ def __init__(self, points, cells): cells = numpy.sort(cells, axis=1) cells = cells[cells[:, 0].argsort()] + self._points = points super().__init__(points, cells) # Assert that all vertices are used. @@ -28,32 +30,32 @@ def __init__(self, points, cells): # ``` # 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(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,11 +76,11 @@ def __init__(self, points, 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 @@ -100,7 +102,7 @@ def __init__(self, points, 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 @@ -111,11 +113,21 @@ def __init__(self, points, cells): return def __repr__(self): - num_nodes = len(self.points) - 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 @@ -124,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.points), 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 @@ -153,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 @@ -171,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]], ] ) @@ -182,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 @@ -221,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.points[self.cells["nodes"]], axis=1 + alpha[None].T * self.points[self.cells["points"]], axis=1 ) return @@ -293,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 @@ -385,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.points[self.cells["nodes"]], axis=1) / 4.0 + numpy.sum(self.points[self.cells["points"]], axis=1) / 4.0 ) return self._cell_centroids @@ -407,7 +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.points[self.cells["nodes"]]) + return numpy.einsum("ij,jik->jk", face_areas, self.points[self.cells["points"]]) @property def cell_inradius(self): @@ -417,7 +429,7 @@ def cell_inradius(self): @property def cell_circumradius(self): - # # Just take the distance of the circumcenter to one of the nodes for now. + # # 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 @@ -515,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( [ @@ -533,7 +545,7 @@ 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.points), ) @@ -578,10 +590,10 @@ def plot(self): self._compute_cell_circumcenters() X = self.points - for cell_id in range(len(self.cells["nodes"])): + 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( @@ -648,11 +660,11 @@ def plot_edge(self, edge_id): ) col = "k" for adj_edge_id in adj_edge_ids: - x = self.points[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.points[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 @@ -661,7 +673,7 @@ def plot_edge(self, edge_id): 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( @@ -746,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.points[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) @@ -786,7 +798,7 @@ def get_sphere_actor(x0, r, rgba): ) if face_circumcenter_rgba is not None: - x = self.points[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, :] @@ -795,7 +807,7 @@ def get_sphere_actor(x0, r, rgba): if control_volume_boundaries_rgba: cell_cc = self.cell_circumcenters[cell_id] - x = self.points[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, :] diff --git a/meshplex/mesh_tri.py b/meshplex/mesh_tri.py index 1756ef4..892b1f3 100644 --- a/meshplex/mesh_tri.py +++ b/meshplex/mesh_tri.py @@ -3,13 +3,14 @@ import numpy -from .base import ( - _BaseMesh, +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"] @@ -17,40 +18,39 @@ 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)]] @@ -109,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) @@ -119,28 +120,55 @@ 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.points + # 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.points) - 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): + # 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.points[self.is_interior_node].shape - # self.points[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 @@ -148,8 +176,8 @@ def euler_characteristic(self): self.create_edges() return ( self.points.shape[0] - - self.edges["nodes"].shape[0] - + self.cells["nodes"].shape[0] + - self.edges["points"].shape[0] + + self.cells["points"].shape[0] ) @property @@ -164,27 +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.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]] - ) + """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 @@ -197,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 @@ -208,19 +232,19 @@ def remove_cells(self, remove_array): 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 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: @@ -251,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 ) @@ -270,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] @@ -294,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), @@ -334,7 +358,7 @@ 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.points), ) @@ -375,9 +399,9 @@ 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]] ) @@ -416,7 +440,7 @@ def signed_cell_areas(self): 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.points[self.cells["nodes"]].T + p = self.points[self.cells["points"]].T # self._signed_cell_areas = ( +p[0][2] * (p[1][0] - p[1][1]) @@ -433,27 +457,27 @@ def mark_boundary(self): assert self.is_boundary_edge is not None - self._is_boundary_node = numpy.zeros(len(self.points), 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): @@ -463,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. @@ -479,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 @@ -508,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) @@ -552,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 @@ -560,7 +583,6 @@ def face_partitions(self): @property def cell_partitions(self): - """""" if self._cell_partitions is None: # Compute the control volumes. Note that # @@ -572,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.points[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 @@ -585,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.points[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 @@ -600,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.points[self.cells["nodes"]]) + return numpy.einsum("ij,jik->jk", abc, self.points[self.cells["points"]]) @property def cell_inradius(self): @@ -656,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 @@ -664,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.points[node_edges] + corner = self.points[point_edges] edge_midpoints = 0.5 * (corner[0] + corner[1]) cc = self.cell_circumcenters @@ -674,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 @@ -682,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] @@ -702,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, @@ -711,7 +732,7 @@ def _compute_integral_x(self): # https://doi.org/10.1002/nme.2187. # ''' # if self.cell_circumcenters is None: - # X = self.points[self.cells['nodes']] + # X = self.points[self.cells['points']] # self.cell_circumcenters = self.compute_triangle_circumcenters(X) # # if 'cells' not in self.edges: @@ -722,28 +743,28 @@ def _compute_integral_x(self): # points2d = self.points[:, :2] # cell_circumcenters2d = self.cell_circumcenters[:, :2] # - # num_nodes = len(points2d) - # 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 * ( - # points2d[node0] + - # points2d[node1] + # points2d[point0] + + # points2d[point1] # ) # cell0 = self.edges['cells'][edge_gid][0] # coedge_midpoint = 0.5 * ( @@ -765,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 - points2d[node0]) * coeffs[0] - # r1 = (coedge_midpoint - points2d[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 = points2d[node1] - points2d[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(): @@ -863,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, @@ -896,12 +917,12 @@ def plot( ax.set_ylim(ymin, ymax) # for k, x in enumerate(self.points): - # if self.is_boundary_node[k]: + # 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: + if show_point_numbers: for i, x in enumerate(self.points): plt.text( x[0], @@ -929,7 +950,7 @@ def plot( plt.tripcolor( self.points[:, 0], self.points[:, 1], - self.cells["nodes"], + self.cells["points"], self.q_radius_ratio, shading="flat", cmap=cmap, @@ -943,7 +964,7 @@ def plot( self.create_edges() # Get edges, cut off z-component. - e = self.points[self.edges["nodes"]][:, :, :2] + e = self.points[self.edges["points"]][:, :, :2] if nondelaunay_edge_color is None: line_segments0 = LineCollection(e, color=mesh_color) @@ -953,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 @@ -971,8 +992,8 @@ def plot( cc = self.cell_circumcenters edge_midpoints = 0.5 * ( - self.points[self.edges["nodes"][:, 0]] - + self.points[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 @@ -993,7 +1014,7 @@ def plot( ax.add_collection(line_segments) if boundary_edge_color: - e = self.points[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) @@ -1027,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. @@ -1047,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.points[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.points[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.points[node_ids[0]] + self.points[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 @@ -1079,7 +1100,7 @@ def plot_vertex(self, node_id, show_ce_ratio=True): [ self.cell_circumcenters[cell_id], edge_midpoint, - self.points[node_id], + self.points[point_id], ] ) ax.fill(q[0], q[1], color="0.5") @@ -1176,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() @@ -1269,7 +1290,7 @@ 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 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 ba6dc03..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, @@ -415,7 +415,7 @@ def test_tetrahedron(): # print(mesh.cells.keys()) # 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.points - mesh2.points # max_diff_coords = max(diff_coords.flatten()) @@ -445,7 +445,7 @@ def test_tetrahedron(): def test_toy_geometric(): mesh = meshplex.read(this_dir / "meshes" / "toy.vtk") - mesh = meshplex.MeshTetra(mesh.points, 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.points) + 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) diff --git a/test/test_mesh_tri.py b/test/test_mesh_tri.py index 04356c1..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 @@ -512,7 +514,7 @@ def test_signed_area2(): 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) @@ -549,7 +551,7 @@ 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_cells = mesh.cells["points"].copy() new_coords = mesh.points.copy() # Assert that some key values are updated properly @@ -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]]) @@ -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(): @@ -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_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(): From 431a6730a00133c031eeb17fd688b6ff2f0520e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 5 Nov 2020 17:29:28 +0100 Subject: [PATCH 4/7] test readme with exdown --- Makefile | 1 + README.md | 29 ++++++++++++++++------------- test/test_readme.py | 19 +++++++++++++++++++ tox.ini | 1 + 4 files changed, 37 insertions(+), 13 deletions(-) create mode 100644 test/test_readme.py 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 02fe704..c90cb0c 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,7 +86,7 @@ mesh.show( # boundary_edge_color=None, # comesh_color=(0.8, 0.8, 0.8), show_axes=False, - ) +) ``` #### Tetrahedra @@ -97,14 +97,17 @@ import numpy import meshplex # Generate tetrahedron -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) +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 @@ -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/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/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 From e7bce2a88221d4835ff850f82daf2978f7dc9fe1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 5 Nov 2020 17:42:41 +0100 Subject: [PATCH 5/7] add changelog --- CHANGELOG.md | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 CHANGELOG.md 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 From 7b75df7cdc3a7a8586940523b9876d236bd1139b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 5 Nov 2020 17:44:09 +0100 Subject: [PATCH 6/7] version bump --- setup.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From d1a50047daab4dbae95ac9a91bb19a38782cf1e1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nico=20Schl=C3=B6mer?= Date: Thu, 5 Nov 2020 17:49:39 +0100 Subject: [PATCH 7/7] skip readme test for need of vtk --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index c90cb0c..cec21ec 100644 --- a/README.md +++ b/README.md @@ -91,7 +91,7 @@ mesh.show( #### Tetrahedra - + ```python import numpy import meshplex