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: