From 9aa1776bd5a029bfcc3536a994aabd1385fcdc25 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 1 Apr 2020 15:35:06 -0500 Subject: [PATCH 1/3] make Mesh work without explicit vertices --- meshmode/mesh/__init__.py | 56 +++++++++++++++--------- meshmode/mesh/refinement/no_adjacency.py | 2 + test/test_meshmode.py | 25 +++++++++++ 3 files changed, 62 insertions(+), 21 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 6d6191aa..fd458044 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -213,7 +213,7 @@ class MeshElementGroup(Record): @property def nelements(self): - return self.vertex_indices.shape[0] + return self.nodes.shape[1] @property def nnodes(self): @@ -279,9 +279,6 @@ class SimplexElementGroup(MeshElementGroup): automatically assigned. """ - if not issubclass(vertex_indices.dtype.type, np.integer): - raise TypeError("vertex_indices must be integral") - if unit_nodes is None: if dim is None: raise TypeError("'dim' must be passed " @@ -294,10 +291,14 @@ class SimplexElementGroup(MeshElementGroup): dims = unit_nodes.shape[0] - if vertex_indices.shape[-1] != dims+1: - raise ValueError("vertex_indices has wrong number of vertices per " - "element. expected: %d, got: %d" % (dims+1, - vertex_indices.shape[-1])) + if vertex_indices is not None: + if not issubclass(vertex_indices.dtype.type, np.integer): + raise TypeError("vertex_indices must be integral") + + if vertex_indices.shape[-1] != dims+1: + raise ValueError("vertex_indices has wrong number of vertices per " + "element. expected: %d, got: %d" % (dims+1, + vertex_indices.shape[-1])) super(SimplexElementGroup, self).__init__(order, vertex_indices, nodes, element_nr_base, node_nr_base, unit_nodes, dim) @@ -349,15 +350,16 @@ class TensorProductElementGroup(MeshElementGroup): automatically assigned. """ - if not issubclass(vertex_indices.dtype.type, np.integer): - raise TypeError("vertex_indices must be integral") - dims = unit_nodes.shape[0] - if vertex_indices.shape[-1] != 2**dims: - raise ValueError("vertex_indices has wrong number of vertices per " - "element. expected: %d, got: %d" % (2**dims, - vertex_indices.shape[-1])) + if vertex_indices is not None: + if not issubclass(vertex_indices.dtype.type, np.integer): + raise TypeError("vertex_indices must be integral") + + if vertex_indices.shape[-1] != 2**dims: + raise ValueError("vertex_indices has wrong number of vertices per " + "element. expected: %d, got: %d" % (2**dims, + vertex_indices.shape[-1])) super(TensorProductElementGroup, self).__init__(order, vertex_indices, nodes, element_nr_base, node_nr_base, unit_nodes) @@ -720,6 +722,9 @@ class Mesh(Record): # }}} + if vertices is None: + is_conforming = False + if not is_conforming: if nodal_adjacency is None: nodal_adjacency = False @@ -753,7 +758,8 @@ class Mesh(Record): self, node_vertex_consistency_tolerance) for g in self.groups: - assert g.vertex_indices.dtype == self.vertex_id_dtype + if g.vertex_indices is not None: + assert g.vertex_indices.dtype == self.vertex_id_dtype if nodal_adjacency: assert nodal_adjacency.neighbors_starts.shape == (self.nelements+1,) @@ -812,7 +818,8 @@ class Mesh(Record): @property def ambient_dim(self): - return self.vertices.shape[0] + from pytools import single_valued + return single_valued(grp.nodes.shape[0] for grp in self.groups) @property def dim(self): @@ -829,13 +836,15 @@ class Mesh(Record): @property def nodal_adjacency(self): + from meshmode import DataUnavailable + if self.vertices is None: + raise DataUnavailable("vertices") + if self._nodal_adjacency is False: - from meshmode import DataUnavailable raise DataUnavailable("nodal_adjacency") elif self._nodal_adjacency is None: if not self.is_conforming: - from meshmode import DataUnavailable raise DataUnavailable("nodal_adjacency can only " "be computed for known-conforming meshes") @@ -852,13 +861,15 @@ class Mesh(Record): @property def facial_adjacency_groups(self): + from meshmode import DataUnavailable + if self.vertices is None: + raise DataUnavailable("vertices") + if self._facial_adjacency_groups is False: - from meshmode import DataUnavailable raise DataUnavailable("facial_adjacency_groups") elif self._facial_adjacency_groups is None: if not self.is_conforming: - from meshmode import DataUnavailable raise DataUnavailable("facial_adjacency_groups can only " "be computed for known-conforming meshes") @@ -904,6 +915,9 @@ class Mesh(Record): # {{{ node-vertex consistency test def _test_node_vertex_consistency_simplex(mesh, mgrp, tol): + if mesh.vertices is None: + return True + if mgrp.nelements == 0: return True diff --git a/meshmode/mesh/refinement/no_adjacency.py b/meshmode/mesh/refinement/no_adjacency.py index 70279bd9..e2576534 100644 --- a/meshmode/mesh/refinement/no_adjacency.py +++ b/meshmode/mesh/refinement/no_adjacency.py @@ -146,6 +146,8 @@ class RefinerWithoutAdjacency(object): """ mesh = self._current_mesh + if mesh.vertices is None: + raise RuntimeError("Meshes without vertices cannot be refined.") refine_flags = np.asarray(refine_flags, dtype=np.bool) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 915b89e8..ef3f6d4c 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1139,6 +1139,31 @@ def test_affine_map(): assert la.norm(x-m_inv(m(x))) < 1e-10 +def test_mesh_without_vertices(ctx_factory): + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + # create a mesh + from meshmode.mesh.generation import generate_icosphere + mesh = generate_icosphere(r=1.0, order=4) + + # create one without the vertices + from meshmode.mesh import Mesh + grp, = mesh.groups + groups = [grp.copy(nodes=grp.nodes, vertex_indices=None) for grp in mesh.groups] + mesh = Mesh(None, groups, is_conforming=False) + + # make sure the world doesn't end + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory as GroupFactory + discr = Discretization(ctx, mesh, GroupFactory(4)) + discr.nodes().with_queue(queue) + + from meshmode.discretization.visualization import make_visualizer + make_visualizer(queue, discr, 4) + + if __name__ == "__main__": import sys if len(sys.argv) > 1: -- GitLab From c63990be0228da4da7533a67a2acad5393f61b1a Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 1 Apr 2020 20:11:27 -0500 Subject: [PATCH 2/3] update RefinerWithoutAdjacency to work without vertices --- meshmode/mesh/__init__.py | 12 ++- meshmode/mesh/refinement/no_adjacency.py | 95 +++++++++++++----------- meshmode/mesh/refinement/resampler.py | 6 +- test/test_meshmode.py | 4 + 4 files changed, 64 insertions(+), 53 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index fd458044..f3d33c37 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -723,7 +723,7 @@ class Mesh(Record): # }}} if vertices is None: - is_conforming = False + is_conforming = None if not is_conforming: if nodal_adjacency is None: @@ -828,6 +828,10 @@ class Mesh(Record): @property def nvertices(self): + if self.vertices is None: + from meshmode import DataUnavailable + raise DataUnavailable("vertices") + return self.vertices.shape[-1] @property @@ -837,9 +841,6 @@ class Mesh(Record): @property def nodal_adjacency(self): from meshmode import DataUnavailable - if self.vertices is None: - raise DataUnavailable("vertices") - if self._nodal_adjacency is False: raise DataUnavailable("nodal_adjacency") @@ -862,9 +863,6 @@ class Mesh(Record): @property def facial_adjacency_groups(self): from meshmode import DataUnavailable - if self.vertices is None: - raise DataUnavailable("vertices") - if self._facial_adjacency_groups is False: raise DataUnavailable("facial_adjacency_groups") diff --git a/meshmode/mesh/refinement/no_adjacency.py b/meshmode/mesh/refinement/no_adjacency.py index e2576534..f37f6621 100644 --- a/meshmode/mesh/refinement/no_adjacency.py +++ b/meshmode/mesh/refinement/no_adjacency.py @@ -146,19 +146,19 @@ class RefinerWithoutAdjacency(object): """ mesh = self._current_mesh - if mesh.vertices is None: - raise RuntimeError("Meshes without vertices cannot be refined.") - refine_flags = np.asarray(refine_flags, dtype=np.bool) if len(refine_flags) != mesh.nelements: raise ValueError("length of refine_flags does not match " "element count of last generated mesh") + perform_vertex_updates = mesh.vertices is not None + new_el_groups = [] group_refinement_records = [] additional_vertices = [] - inew_vertex = mesh.nvertices + if perform_vertex_updates: + inew_vertex = mesh.nvertices for igrp, group in enumerate(mesh.groups): bisection_info = self._get_bisection_tesselation_info( @@ -197,55 +197,59 @@ class RefinerWithoutAdjacency(object): # {{{ get new vertices together - midpoints = bisection_info.resampler.get_midpoints( - group, bisection_info, refining_el_old_indices) + if perform_vertex_updates: + midpoints = bisection_info.resampler.get_midpoints( + group, bisection_info, refining_el_old_indices) - new_vertex_indices = np.empty( - (new_nelements, group.vertex_indices.shape[1]), - dtype=mesh.vertex_id_dtype) - new_vertex_indices.fill(-17) + new_vertex_indices = np.empty( + (new_nelements, group.vertex_indices.shape[1]), + dtype=mesh.vertex_id_dtype) + new_vertex_indices.fill(-17) - # copy over unchanged vertices - new_vertex_indices[unrefined_el_new_indices] = \ - group.vertex_indices[~grp_flags] + # copy over unchanged vertices + new_vertex_indices[unrefined_el_new_indices] = \ + group.vertex_indices[~grp_flags] - for old_iel in refining_el_old_indices: - new_iel_base = child_el_indices[old_iel] + for old_iel in refining_el_old_indices: + new_iel_base = child_el_indices[old_iel] - refining_vertices = np.empty(len(bisection_info.ref_vertices), + refining_vertices = np.empty(len(bisection_info.ref_vertices), dtype=mesh.vertex_id_dtype) - refining_vertices.fill(-17) + refining_vertices.fill(-17) - # carry over old vertices - refining_vertices[bisection_info.orig_vertex_indices] = \ - group.vertex_indices[old_iel] + # carry over old vertices + refining_vertices[bisection_info.orig_vertex_indices] = \ + group.vertex_indices[old_iel] - for imidpoint, (iref_midpoint, (v1, v2)) in enumerate(zip( - bisection_info.midpoint_indices, - bisection_info.midpoint_vertex_pairs)): + for imidpoint, (iref_midpoint, (v1, v2)) in enumerate(zip( + bisection_info.midpoint_indices, + bisection_info.midpoint_vertex_pairs)): - global_v1 = group.vertex_indices[old_iel, v1] - global_v2 = group.vertex_indices[old_iel, v2] + global_v1 = group.vertex_indices[old_iel, v1] + global_v2 = group.vertex_indices[old_iel, v2] - if global_v1 > global_v2: - global_v1, global_v2 = global_v2, global_v1 + if global_v1 > global_v2: + global_v1, global_v2 = global_v2, global_v1 - try: - global_midpoint = self.global_vertex_pair_to_midpoint[ - global_v1, global_v2] - except KeyError: - global_midpoint = inew_vertex - additional_vertices.append(midpoints[old_iel][:, imidpoint]) - inew_vertex += 1 + try: + global_midpoint = self.global_vertex_pair_to_midpoint[ + global_v1, global_v2] + except KeyError: + global_midpoint = inew_vertex + additional_vertices.append( + midpoints[old_iel][:, imidpoint]) + inew_vertex += 1 - refining_vertices[iref_midpoint] = global_midpoint + refining_vertices[iref_midpoint] = global_midpoint - assert (refining_vertices >= 0).all() + assert (refining_vertices >= 0).all() - new_vertex_indices[new_iel_base:new_iel_base+nchildren] = \ - refining_vertices[bisection_info.children] + new_vertex_indices[new_iel_base:new_iel_base+nchildren] = \ + refining_vertices[bisection_info.children] - assert (new_vertex_indices >= 0).all() + assert (new_vertex_indices >= 0).all() + else: + new_vertex_indices = None # }}} @@ -279,11 +283,14 @@ class RefinerWithoutAdjacency(object): nodes=new_nodes, unit_nodes=group.unit_nodes)) - new_vertices = np.empty( - (mesh.ambient_dim, mesh.nvertices + len(additional_vertices)), - mesh.vertices.dtype) - new_vertices[:, :mesh.nvertices] = mesh.vertices - new_vertices[:, mesh.nvertices:] = np.array(additional_vertices).T + if perform_vertex_updates: + new_vertices = np.empty( + (mesh.ambient_dim, mesh.nvertices + len(additional_vertices)), + mesh.vertices.dtype) + new_vertices[:, :mesh.nvertices] = mesh.vertices + new_vertices[:, mesh.nvertices:] = np.array(additional_vertices).T + else: + new_vertices = None from meshmode.mesh import Mesh new_mesh = Mesh(new_vertices, new_el_groups, is_conforming=( diff --git a/meshmode/mesh/refinement/resampler.py b/meshmode/mesh/refinement/resampler.py index 794c1c5a..4631c9ca 100644 --- a/meshmode/mesh/refinement/resampler.py +++ b/meshmode/mesh/refinement/resampler.py @@ -74,7 +74,8 @@ class SimplexResampler(object): follows their ordering in the tesselation (see also :meth:`SimplexResampler.get_vertex_pair_to_midpoint_order`) """ - assert len(group.vertex_indices[0]) == group.dim + 1 + if group.vertex_indices is not None: + assert len(group.vertex_indices[0]) == group.dim + 1 # Get midpoints, converted to unit coordinates. midpoints = -1 + np.array([vertex for vertex in @@ -106,7 +107,8 @@ class SimplexResampler(object): The ordering of the child nodes follows the ordering of ``tesselation.children.`` """ - assert len(group.vertex_indices[0]) == group.dim + 1 + if group.vertex_indices is not None: + assert len(group.vertex_indices[0]) == group.dim + 1 from meshmode.mesh.refinement.utils import map_unit_nodes_to_children diff --git a/test/test_meshmode.py b/test/test_meshmode.py index ef3f6d4c..c73c0292 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1153,6 +1153,10 @@ def test_mesh_without_vertices(ctx_factory): groups = [grp.copy(nodes=grp.nodes, vertex_indices=None) for grp in mesh.groups] mesh = Mesh(None, groups, is_conforming=False) + # try refining it + from meshmode.mesh.refinement import refine_uniformly + mesh = refine_uniformly(mesh, 1) + # make sure the world doesn't end from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ -- GitLab From c50dac6e98d03b50d1cebc5888b3a1d4f0432f7f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Thu, 16 Apr 2020 07:44:59 +0200 Subject: [PATCH 3/3] Apply suggestion to meshmode/mesh/__init__.py --- meshmode/mesh/__init__.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index f3d33c37..211e764c 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -570,8 +570,9 @@ class Mesh(Record): """ .. attribute:: vertices - An array of vertex coordinates with shape - *(ambient_dim, nvertices)* + *None* or an array of vertex coordinates with shape + *(ambient_dim, nvertices)*. If *None*, vertices are not + known for this mesh. .. attribute:: groups -- GitLab