From ffe42be74e789991a2573221ea4745f9d3474b51 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Tue, 8 Jul 2014 08:31:09 -0500
Subject: [PATCH] Mesh connectivity representation, ico{sahedron,sphere}
generation
---
examples/plot-connectivity.py | 41 ++++++++
meshmode/discretization/visualization.py | 61 ++++++++++++
meshmode/mesh/__init__.py | 121 ++++++++++++++++++++++-
meshmode/mesh/generation.py | 92 +++++++++--------
meshmode/mesh/io.py | 2 +-
5 files changed, 270 insertions(+), 47 deletions(-)
create mode 100644 examples/plot-connectivity.py
diff --git a/examples/plot-connectivity.py b/examples/plot-connectivity.py
new file mode 100644
index 0000000..a888636
--- /dev/null
+++ b/examples/plot-connectivity.py
@@ -0,0 +1,41 @@
+from __future__ import division
+
+import numpy as np # noqa
+import pyopencl as cl
+
+
+order = 3
+
+
+def main():
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ from meshmode.mesh.generation import ( # noqa
+ generate_icosphere, generate_icosahedron)
+ #mesh = generate_icosphere(1, order=order)
+ mesh = generate_icosahedron(1, order=order)
+
+ from meshmode.discretization import Discretization
+ from meshmode.discretization.poly_element import \
+ PolynomialWarpAndBlendGroupFactory
+
+ discr = Discretization(
+ cl_ctx, mesh, PolynomialWarpAndBlendGroupFactory(order))
+
+ from meshmode.discretization.visualization import make_visualizer
+ vis = make_visualizer(queue, discr, order)
+
+ vis.write_vtk_file("geometry.vtu", [
+ ("f", discr.nodes()[0]),
+ ])
+
+ from meshmode.discretization.visualization import \
+ write_mesh_connectivity_vtk_file
+
+ write_mesh_connectivity_vtk_file("connectivity.vtu",
+ mesh)
+
+
+if __name__ == "__main__":
+ main()
diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py
index 815df95..a2b251e 100644
--- a/meshmode/discretization/visualization.py
+++ b/meshmode/discretization/visualization.py
@@ -31,9 +31,13 @@ __doc__ = """
.. autofunction:: make_visualizer
.. autoclass:: Visualizer
+
+.. autofunction:: write_mesh_connectivity_vtk_file
"""
+# {{{ visualizer
+
def separate_by_real_and_imag(data, real_only):
for name, field in data:
from pytools.obj_array import log_shape
@@ -66,6 +70,11 @@ def separate_by_real_and_imag(data, real_only):
class Visualizer(object):
+ """
+ .. automethod:: show_scalar_in_mayavi
+ .. automethod:: write_vtk_file
+ """
+
def __init__(self, discr, vis_discr, connection):
self.discr = discr
self.vis_discr = vis_discr
@@ -233,4 +242,56 @@ def make_visualizer(queue, discr, vis_order):
return Visualizer(discr, vis_discr, cnx)
+# }}}
+
+
+# {{{ connectivity
+
+def write_mesh_connectivity_vtk_file(file_name, mesh, compressor=None,):
+ from pyvisfile.vtk import (
+ UnstructuredGrid, DataArray,
+ AppendedDataXMLGenerator,
+ VTK_LINE,
+ VF_LIST_OF_COMPONENTS)
+
+ centroids = np.empty(
+ (mesh.ambient_dim, mesh.nelements),
+ dtype=mesh.vertices.dtype)
+
+ for grp in mesh.groups:
+ iel_base = grp.element_nr_base
+ centroids[:, iel_base:iel_base+grp.nelements] = (
+ np.sum(mesh.vertices[:, grp.vertex_indices], axis=-1)
+ / grp.vertex_indices.shape[-1])
+
+ cnx = mesh.element_connectivity
+
+ nconnections = len(cnx.neighbors)
+ connections = np.empty((nconnections, 2), dtype=np.int32)
+
+ nb_starts = cnx.neighbors_starts
+ for iel in xrange(mesh.nelements):
+ connections[nb_starts[iel]:nb_starts[iel+1], 0] = iel
+
+ connections[:, 1] = cnx.neighbors
+
+ grid = UnstructuredGrid(
+ (mesh.nelements,
+ DataArray("points",
+ centroids,
+ vector_format=VF_LIST_OF_COMPONENTS)),
+ cells=connections.reshape(-1),
+ cell_types=np.asarray([VTK_LINE] * nconnections,
+ dtype=np.uint8))
+
+ from os.path import exists
+ if exists(file_name):
+ raise RuntimeError("output file '%s' already exists"
+ % file_name)
+
+ with open(file_name, "w") as outf:
+ AppendedDataXMLGenerator(compressor)(grid).write(outf)
+
+# }}}
+
# vim: foldmethod=marker
diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index 7d19b1a..bef25c2 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -219,6 +219,24 @@ class SimplexElementGroup(MeshElementGroup):
# {{{ mesh
+class ElementConnectivity(Record):
+ """
+ .. attribute:: neighbors_starts
+
+ ``element_id_t [nelments+1]``
+
+ Use together with :attr:`neighbors`. ``neighbors_starts[iel]`` and
+ ``neighbors_starts[iel+1]`` together indicate a ranges of element indices
+ :attr:`neighbors` which are adjacent to *iel*.
+
+ .. attribute:: neighbors
+
+ ``element_id_t []``
+
+ See :attr:`neighbors_starts`.
+ """
+
+
class Mesh(Record):
"""
.. attribute:: vertices
@@ -229,12 +247,38 @@ class Mesh(Record):
.. attribute:: groups
A list of :class:`MeshElementGroup` instances.
+
+ .. attribute:: element_connectivity
+
+ An instance of :class:`ElementConnectivity`.
+
+ Referencing this attribute may raise
+ :exc:`meshmode.ConnectivityUnavailable`.
+
+ .. attribute:: vertex_id_dtype
+
+ .. attribute:: element_id_dtype
"""
- def __init__(self, vertices, groups, skip_tests=False):
+ def __init__(self, vertices, groups, skip_tests=False,
+ element_connectivity=False,
+ vertex_id_dtype=np.int32,
+ element_id_dtype=np.int32):
"""
+ The following are keyword-only:
+
:arg skip_tests: Skip mesh tests, in case you want to load a broken
mesh anyhow and then fix it inside of this data structure.
+ :arg element_connectivity: One of three options:
+ *None*, in which case this information
+ will be deduced from vertex adjacency. *False*, in which case
+ this information will be marked unavailable (such as if there are
+ hanging nodes in the geometry, so that vertex adjacency does not convey
+ the full picture), and references to
+ :attr:`element_neighbors_starts` and :attr:`element_neighbors`
+ will result in exceptions. Lastly, a tuple
+ *(element_neighbors_starts, element_neighbors)*, representing the
+ correspondingly-named attributes.
"""
el_nr = 0
node_nr = 0
@@ -246,11 +290,28 @@ class Mesh(Record):
el_nr += ng.nelements
node_nr += ng.nnodes
- Record.__init__(self, vertices=vertices, groups=new_groups)
+ if element_connectivity is not False and element_connectivity is not None:
+ nb_starts, nbs = element_connectivity
+ element_connectivity = ElementConnectivity(
+ element_neighbors_starts=nb_starts,
+ element_neighbors=nbs)
+
+ del nb_starts
+ del nbs
+
+ Record.__init__(
+ self, vertices=vertices, groups=new_groups,
+ _element_connectivity=element_connectivity,
+ vertex_id_dtype=np.dtype(vertex_id_dtype),
+ element_id_dtype=np.dtype(element_id_dtype),
+ )
if not skip_tests:
assert _test_node_vertex_consistency(self)
+ for g in self.groups:
+ assert g.vertex_indices.dtype == self.vertex_id_dtype
+
from meshmode.mesh.processing import \
test_volume_mesh_element_orientations
@@ -272,9 +333,21 @@ class Mesh(Record):
def nelements(self):
return sum(grp.nelements for grp in self.groups)
- # Design experience: Try not to add too many global data structures
- # to the mesh. Let the element groups be responsible for that.
- # The interpolatory discretization has these big global things.
+ @property
+ def element_connectivity(self):
+ if self._element_connectivity is False:
+ from meshmode import ConnectivityUnavailable
+ raise ConnectivityUnavailable()
+ elif self._element_connectivity is None:
+ self._element_connectivity = _compute_connectivity_from_vertices(self)
+
+ return self._element_connectivity
+
+
+ # Design experience: Try not to add too many global data structures to the
+ # mesh. Let the element groups be responsible for that at the mesh level.
+ #
+ # There are more big, global structures on the discretization level.
# }}}
@@ -328,4 +401,42 @@ def _test_node_vertex_consistency(mesh):
# }}}
+# {{{ vertex-based connectivity
+
+def _compute_connectivity_from_vertices(mesh):
+ # FIXME Native code would make this faster
+
+ _, nvertices = mesh.vertices.shape
+ vertex_to_element = [[] for i in xrange(nvertices)]
+
+ for grp in mesh.groups:
+ iel_base = grp.element_nr_base
+ for iel_grp in xrange(grp.nelements):
+ for ivertex in grp.vertex_indices[iel_grp]:
+ vertex_to_element[ivertex].append(iel_base + iel_grp)
+
+ element_to_element = [set() for i in xrange(mesh.nelements)]
+ for grp in mesh.groups:
+ iel_base = grp.element_nr_base
+ for iel_grp in xrange(grp.nelements):
+ for ivertex in grp.vertex_indices[iel_grp]:
+ element_to_element[iel_base + iel_grp].update(
+ vertex_to_element[ivertex])
+
+ lengths = [len(el_list) for el_list in element_to_element]
+ neighbors_starts = np.cumsum(
+ np.array([0] + lengths, dtype=mesh.element_id_dtype))
+ from pytools import flatten
+ neighbors = np.array(
+ list(flatten(element_to_element)),
+ dtype=mesh.element_id_dtype)
+
+ assert neighbors_starts[-1] == len(neighbors)
+
+ return ElementConnectivity(
+ neighbors_starts=neighbors_starts,
+ neighbors=neighbors)
+
+# }}}
+
# vim: foldmethod=marker
diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index cbd4363..5df058f 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -138,6 +138,34 @@ def qbx_peanut(t):
# }}}
+def make_group_from_vertices(vertices, vertex_indices, order):
+ el_vertices = vertices[:, vertex_indices]
+
+ el_origins = el_vertices[:, :, 0][:, :, np.newaxis]
+ # ambient_dim, nelements, nspan_vectors
+ spanning_vectors = (
+ el_vertices[:, :, 1:] - el_origins)
+
+ nspan_vectors = spanning_vectors.shape[-1]
+ dim = nspan_vectors
+
+ # dim, nunit_nodes
+ unit_nodes = mp.warp_and_blend_nodes(dim, order)
+ unit_nodes_01 = 0.5 + 0.5*unit_nodes
+
+ nodes = np.einsum(
+ "si,des->dei",
+ unit_nodes_01, spanning_vectors) + el_origins
+
+ # make contiguous
+ nodes = nodes.copy()
+
+ from meshmode.mesh import SimplexElementGroup
+ return SimplexElementGroup(
+ order, vertex_indices, nodes,
+ unit_nodes=unit_nodes)
+
+
def make_curve_mesh(curve_f, element_boundaries, order):
"""
:arg curve_f: A callable representing a parametrization for a curve,
@@ -175,23 +203,23 @@ def make_curve_mesh(curve_f, element_boundaries, order):
nodes=nodes,
unit_nodes=unodes)
- return Mesh(vertices=vertices, groups=[egroup])
+ return Mesh(vertices=vertices, groups=[egroup],
+ element_connectivity=None)
-def generate_icosahedron(r):
+def generate_icosahedron(r, order):
# http://en.wikipedia.org/w/index.php?title=Icosahedron&oldid=387737307
phi = (1+5**(1/2))/2
from pytools import flatten
- pts = np.array(sorted(flatten([
+ vertices = np.array(sorted(flatten([
(0, pm1*1, pm2*phi),
(pm1*1, pm2*phi, 0),
(pm1*phi, 0, pm2*1)]
for pm1 in [-1, 1]
- for pm2 in [-1, 1])))
+ for pm2 in [-1, 1]))).T.copy()
- from meshmode.mesh import Mesh
top_ring = [11, 7, 1, 2, 8]
bottom_ring = [10, 9, 3, 0, 4]
bottom_point = 6
@@ -205,55 +233,34 @@ def generate_icosahedron(r):
tris.append([bottom_ring[i], bottom_ring[(i+1) % l], top_ring[i]])
tris.append([top_ring[i], bottom_ring[(i+1) % l], top_ring[(i+1) % l]])
- pts *= r/la.norm(pts[0])
-
- return Mesh(pts, tris)
+ vertices *= r/la.norm(vertices[:, 0])
+ vertex_indices = np.array(tris, dtype=np.int32)
-def generate_icosphere(r, refinements, order=1):
- mesh = generate_icosahedron(r)
+ grp = make_group_from_vertices(vertices, vertex_indices, order)
from meshmode.mesh import Mesh
- from meshmode.mesh import bisect_geometry
- for i in range(refinements):
- mesh, _ = bisect_geometry(mesh)
-
- pts = mesh.vertex_coordinates
- pt_norms = np.sum(pts**2, axis=-1)**0.5
- pts = pts * r / pt_norms[:, np.newaxis]
-
- mesh = Mesh(pts, mesh.elements)
-
- if order > 1:
- from meshmode.mesh.trianglemaps import interv_pts_on_unit
- eq_nodes = interv_pts_on_unit(order, 0, 1)
-
- # indices: triangle #, node #, xyz axis
- mapped_eq_nodes = np.empty((len(mesh), len(eq_nodes), 3))
+ return Mesh(vertices, [grp],
+ element_connectivity=None)
- for iel, el_nodes in enumerate(mesh.elements):
- n1, n2, n3 = mesh.vertex_coordinates[el_nodes]
- a = np.vstack([n2-n1, n3-n1]).T
- b = n1
- mapped_eq_nodes[iel] = np.dot(a, eq_nodes.T).T + b
+def generate_icosphere(r, order):
+ mesh = generate_icosahedron(r, order)
- mn = mapped_eq_nodes
- mapped_eq_nodes = \
- r * mn / np.sqrt(np.sum(mn**2, axis=-1))[:, :, np.newaxis]
+ grp, = mesh.groups
- reshaped_nodes = mapped_eq_nodes.reshape(-1, 3)
+ grp = grp.copy(
+ nodes=grp.nodes * r / np.sqrt(np.sum(grp.nodes**2, axis=0)))
- mesh = Mesh(
- vertex_coordinates=reshaped_nodes,
- elements=np.arange(len(reshaped_nodes)).reshape(len(mesh), -1),
- order=order)
-
- return mesh
+ from meshmode.mesh import Mesh
+ return Mesh(mesh.vertices, [grp],
+ element_connectivity=None)
def generate_torus_and_cycle_vertices(r_outer, r_inner,
n_outer=20, n_inner=10, order=1):
+ raise NotImplementedError("not yet updated for meshmode")
+
a = r_outer
b = r_inner
u, v = np.mgrid[0:2*np.pi:2*np.pi/n_outer, 0:2*np.pi:2*np.pi/n_inner]
@@ -333,6 +340,9 @@ def generate_torus_and_cycle_vertices(r_outer, r_inner,
def generate_torus(r_outer, r_inner, n_outer=20, n_inner=10, order=1):
+ # FIXME
+ raise NotImplementedError("not yet updated for meshmode")
+
geo, a_cycle, b_cycle = generate_torus_and_cycle_vertices(
r_outer, r_inner, n_outer, n_inner, order)
return geo
diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py
index bdfbf28..dc540a8 100644
--- a/meshmode/mesh/io.py
+++ b/meshmode/mesh/io.py
@@ -167,7 +167,7 @@ class GmshMeshReceiver(GmshMeshReceiverBase):
groups.append(group)
- return Mesh(vertices, groups)
+ return Mesh(vertices, groups, element_connectivity=None)
# }}}
--
GitLab