From 9aa1776bd5a029bfcc3536a994aabd1385fcdc25 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 1 Apr 2020 15:35:06 -0500 Subject: [PATCH] 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