diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 6d6191aa394544b47f8b994ddc6ea7c8691b88a0..211e764c87b95310d5b68995677f0d958a8c9550 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) @@ -568,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 @@ -720,6 +723,9 @@ class Mesh(Record): # }}} + if vertices is None: + is_conforming = None + if not is_conforming: if nodal_adjacency is None: nodal_adjacency = False @@ -753,7 +759,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 +819,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): @@ -821,6 +829,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 @@ -829,13 +841,12 @@ class Mesh(Record): @property def nodal_adjacency(self): + from meshmode import DataUnavailable 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 +863,12 @@ class Mesh(Record): @property def facial_adjacency_groups(self): + from meshmode import DataUnavailable 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 +914,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 70279bd948163d501e3b5294f6afb2397f57ba9b..f37f66215ded63f530db8d5a2cb465d38771b5db 100644 --- a/meshmode/mesh/refinement/no_adjacency.py +++ b/meshmode/mesh/refinement/no_adjacency.py @@ -146,17 +146,19 @@ class RefinerWithoutAdjacency(object): """ mesh = self._current_mesh - 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( @@ -195,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 # }}} @@ -277,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 794c1c5a2362d1255f5e3d329c48fb766c6229a9..4631c9caa7220ee7a742670d10b0a7a5998003bd 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 915b89e8c0a569c064ba41c801743d223fc5bce3..c73c029264ecf625d2c0b2e684c153befa2221f9 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1139,6 +1139,35 @@ 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) + + # 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 \ + 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: