From fcffecccfe2a6b5acd032717f30deec086995262 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Mon, 5 Oct 2015 10:00:12 -0500
Subject: [PATCH] Rename ElementConnecitvity -> NodalAdjacency, in preparation
for FacialAdjacency
---
examples/plot-connectivity.py | 4 +-
meshmode/__init__.py | 4 +-
meshmode/discretization/visualization.py | 14 ++--
meshmode/mesh/__init__.py | 72 ++++++++---------
meshmode/mesh/generation.py | 10 +--
meshmode/mesh/io.py | 2 +-
meshmode/mesh/processing.py | 2 +-
meshmode/mesh/refinement.py | 99 ++++++++++++++++++++++++
meshmode/mesh/visualization.py | 6 +-
test/test_meshmode.py | 6 +-
10 files changed, 160 insertions(+), 59 deletions(-)
diff --git a/examples/plot-connectivity.py b/examples/plot-connectivity.py
index 5b1a0c1..068983c 100644
--- a/examples/plot-connectivity.py
+++ b/examples/plot-connectivity.py
@@ -33,9 +33,9 @@ def main():
])
from meshmode.discretization.visualization import \
- write_mesh_connectivity_vtk_file
+ write_nodal_adjacency_vtk_file
- write_mesh_connectivity_vtk_file("connectivity.vtu",
+ write_nodal_adjacency_vtk_file("adjacency.vtu",
mesh)
diff --git a/meshmode/__init__.py b/meshmode/__init__.py
index e1a4c5e..1c38d87 100644
--- a/meshmode/__init__.py
+++ b/meshmode/__init__.py
@@ -25,7 +25,7 @@ THE SOFTWARE.
__doc__ = """
.. exception:: Error
-.. exception:: ConnectivityUnavailable
+.. exception:: DataUnavailable
"""
@@ -33,5 +33,5 @@ class Error(RuntimeError):
pass
-class ConnectivityUnavailable(Error):
+class DataUnavailable(Error):
pass
diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py
index 099527c..289d3a6 100644
--- a/meshmode/discretization/visualization.py
+++ b/meshmode/discretization/visualization.py
@@ -34,7 +34,7 @@ __doc__ = """
.. autoclass:: Visualizer
-.. autofunction:: write_mesh_connectivity_vtk_file
+.. autofunction:: write_nodal_adjacency_vtk_file
"""
@@ -265,9 +265,9 @@ def draw_curve(discr):
# }}}
-# {{{ connectivity
+# {{{ adjacency
-def write_mesh_connectivity_vtk_file(file_name, mesh, compressor=None,):
+def write_nodal_adjacency_vtk_file(file_name, mesh, compressor=None,):
from pyvisfile.vtk import (
UnstructuredGrid, DataArray,
AppendedDataXMLGenerator,
@@ -284,16 +284,16 @@ def write_mesh_connectivity_vtk_file(file_name, mesh, compressor=None,):
np.sum(mesh.vertices[:, grp.vertex_indices], axis=-1)
/ grp.vertex_indices.shape[-1])
- cnx = mesh.element_connectivity
+ adj = mesh.nodal_adjacency
- nconnections = len(cnx.neighbors)
+ nconnections = len(adj.neighbors)
connections = np.empty((nconnections, 2), dtype=np.int32)
- nb_starts = cnx.neighbors_starts
+ nb_starts = adj.neighbors_starts
for iel in range(mesh.nelements):
connections[nb_starts[iel]:nb_starts[iel+1], 0] = iel
- connections[:, 1] = cnx.neighbors
+ connections[:, 1] = adj.neighbors
grid = UnstructuredGrid(
(mesh.nelements,
diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index 566f211..4bf68cc 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -39,7 +39,7 @@ __doc__ = """
:members:
:undoc-members:
-.. autoclass:: ElementConnectivity
+.. autoclass:: NodalAdjacency
.. autofunction:: as_python
@@ -239,8 +239,10 @@ class SimplexElementGroup(MeshElementGroup):
# {{{ mesh
-class ElementConnectivity(Record):
- """
+class NodalAdjacency(Record):
+ """Describes element adjacency information, i.e. information about
+ elements that touch in at least one point.
+
.. attribute:: neighbors_starts
``element_id_t [nelments+1]``
@@ -281,12 +283,12 @@ class Mesh(Record):
A list of :class:`MeshElementGroup` instances.
- .. attribute:: element_connectivity
+ .. attribute:: nodal_adjacency
- An instance of :class:`ElementConnectivity`.
+ An instance of :class:`NodalAdjacency`.
Referencing this attribute may raise
- :exc:`meshmode.ConnectivityUnavailable`.
+ :exc:`meshmode.DataUnavailable`.
.. attribute:: vertex_id_dtype
@@ -297,7 +299,7 @@ class Mesh(Record):
"""
def __init__(self, vertices, groups, skip_tests=False,
- element_connectivity=False,
+ nodal_adjacency=False,
vertex_id_dtype=np.int32,
element_id_dtype=np.int32):
"""
@@ -305,7 +307,7 @@ class Mesh(Record):
: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:
+ :arg nodal_adjacency: 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
@@ -326,9 +328,9 @@ class Mesh(Record):
el_nr += ng.nelements
node_nr += ng.nnodes
- if element_connectivity is not False and element_connectivity is not None:
- nb_starts, nbs = element_connectivity
- element_connectivity = ElementConnectivity(
+ if nodal_adjacency is not False and nodal_adjacency is not None:
+ nb_starts, nbs = nodal_adjacency
+ nodal_adjacency = NodalAdjacency(
neighbors_starts=nb_starts,
neighbors=nbs)
@@ -337,7 +339,7 @@ class Mesh(Record):
Record.__init__(
self, vertices=vertices, groups=new_groups,
- _element_connectivity=element_connectivity,
+ _nodal_adjacency=nodal_adjacency,
vertex_id_dtype=np.dtype(vertex_id_dtype),
element_id_dtype=np.dtype(element_id_dtype),
)
@@ -369,25 +371,25 @@ class Mesh(Record):
return sum(grp.nelements for grp in self.groups)
@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)
+ def nodal_adjacency(self):
+ if self._nodal_adjacency is False:
+ from meshmode import DataUnavailable
+ raise DataUnavailable("nodal_adjacency")
+ elif self._nodal_adjacency is None:
+ self._nodal_adjacency = _compute_nodal_adjacency_from_vertices(self)
- return self._element_connectivity
+ return self._nodal_adjacency
- def connectivity_init_arg(self):
- """Returns an 'elment_connectivity' argument that can be
+ def nodal_adjacency_init_arg(self):
+ """Returns an 'nodal_adjacency' argument that can be
passed to a Mesh constructor.
"""
- if isinstance(self._element_connectivity, ElementConnectivity):
- return (self._element_connectivity.neighbors_starts,
- self._element_connectivity.neighbors)
+ if isinstance(self._nodal_adjacency, NodalAdjacency):
+ return (self._nodal_adjacency.neighbors_starts,
+ self._nodal_adjacency.neighbors)
else:
- return self._element_connectivity
+ return self._nodal_adjacency
def __eq__(self, other):
return (
@@ -396,8 +398,8 @@ class Mesh(Record):
and self.groups == other.groups
and self.vertex_id_dtype == other.vertex_id_dtype
and self.element_id_dtype == other.element_id_dtype
- and (self._element_connectivity
- == other._element_connectivity))
+ and (self._nodal_adjacency
+ == other._nodal_adjacency))
def __ne__(self, other):
return not self.__eq__(other)
@@ -459,9 +461,9 @@ def _test_node_vertex_consistency(mesh):
# }}}
-# {{{ vertex-based connectivity
+# {{{ vertex-based adjacency
-def _compute_connectivity_from_vertices(mesh):
+def _compute_nodal_adjacency_from_vertices(mesh):
# FIXME Native code would make this faster
_, nvertices = mesh.vertices.shape
@@ -494,7 +496,7 @@ def _compute_connectivity_from_vertices(mesh):
assert neighbors_starts[-1] == len(neighbors)
- return ElementConnectivity(
+ return NodalAdjacency(
neighbors_starts=neighbors_starts,
neighbors=neighbors)
@@ -545,17 +547,17 @@ def as_python(mesh, function_name="make_mesh"):
cg(" vertex_id_dtype=np.%s," % mesh.vertex_id_dtype.name)
cg(" element_id_dtype=np.%s," % mesh.element_id_dtype.name)
- if isinstance(mesh._element_connectivity, ElementConnectivity):
+ if isinstance(mesh._nodal_adjacency, NodalAdjacency):
el_con_str = "(%s, %s)" % (
_numpy_array_as_python(
- mesh._element_connectivity.neighbors_starts),
+ mesh._nodal_adjacency.neighbors_starts),
_numpy_array_as_python(
- mesh._element_connectivity.neighbors),
+ mesh._nodal_adjacency.neighbors),
)
else:
- el_con_str = repr(mesh._element_connectivity)
+ el_con_str = repr(mesh._nodal_adjacency)
- cg(" element_connectivity=%s)" % el_con_str)
+ cg(" nodal_adjacency=%s)" % el_con_str)
return cg.get()
diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index 9a72c53..0f130a5 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -213,7 +213,7 @@ def make_curve_mesh(curve_f, element_boundaries, order):
unit_nodes=unodes)
return Mesh(vertices=vertices, groups=[egroup],
- element_connectivity=None)
+ nodal_adjacency=None)
# }}}
@@ -286,7 +286,7 @@ def generate_icosahedron(r, order):
from meshmode.mesh import Mesh
return Mesh(vertices, [grp],
- element_connectivity=None)
+ nodal_adjacency=None)
# }}}
@@ -303,7 +303,7 @@ def generate_icosphere(r, order):
from meshmode.mesh import Mesh
return Mesh(mesh.vertices, [grp],
- element_connectivity=None)
+ nodal_adjacency=None)
# }}}
@@ -365,7 +365,7 @@ def generate_torus_and_cycle_vertices(r_outer, r_inner,
nodes[2] = b*np.sin(minor_theta)
from meshmode.mesh import Mesh
- return (Mesh(vertices, [grp.copy(nodes=nodes)], element_connectivity=None),
+ return (Mesh(vertices, [grp.copy(nodes=nodes)], nodal_adjacency=None),
[idx(i, 0) for i in range(n_outer)],
[idx(0, j) for j in range(n_inner)])
@@ -470,7 +470,7 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64):
from meshmode.mesh import Mesh
return Mesh(vertices, [grp],
- element_connectivity=None)
+ nodal_adjacency=None)
# }}}
diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py
index 6c1078e..24bcae5 100644
--- a/meshmode/mesh/io.py
+++ b/meshmode/mesh/io.py
@@ -171,7 +171,7 @@ class GmshMeshReceiver(GmshMeshReceiverBase):
groups.append(group)
- return Mesh(vertices, groups, element_connectivity=None)
+ return Mesh(vertices, groups, nodal_adjacency=None)
# }}}
diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index dfd4bd2..2b98957 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -298,7 +298,7 @@ def affine_map(mesh, A=None, b=None): # noqa
from meshmode.mesh import Mesh
return Mesh(vertices, new_groups, skip_tests=True,
- element_connectivity=mesh.connectivity_init_arg())
+ nodal_adjacency=mesh.nodal_adjacency_init_arg())
# }}}
diff --git a/meshmode/mesh/refinement.py b/meshmode/mesh/refinement.py
index 806fe5a..efe92c7 100644
--- a/meshmode/mesh/refinement.py
+++ b/meshmode/mesh/refinement.py
@@ -25,6 +25,105 @@ THE SOFTWARE.
import numpy as np
+# {{{ simple, no-connectivity refinement
+
+def _refine_simplex_el_group(grp, refine_flags, vertex_nr_base):
+ from pytools import (
+ generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam)
+ from modepy.tools import submesh
+
+ node_tuples = gnitstam(2, grp.dim)
+ template = submesh(node_tuples)
+
+ nrefined_elements = np.sum(refine_flags != 0)
+
+ old_vertex_tuples = (
+ [(0,) * grp.dim]
+ + [(0,)*i + (2,) + (0,)*(grp.dim-1-i) for i in range(grp.dim)])
+ new_vertex_tuples = [
+ vertex_tuple
+ for vertex_tuple in node_tuples
+ if vertex_tuple not in old_vertex_tuples]
+
+ new_unit_nodes = np.array(new_vertex_tuples, dtype=grp.nodes.dtype) - 1
+
+ nnew_vertices = nrefined_elements * len(new_vertex_tuples)
+ ambient_dim = grp.nodes.shape[0]
+ new_vertices = np.empty((nnew_vertices, ambient_dim))
+
+ nnew_elements = (len(template) - 1) * nrefined_elements
+ new_vertex_indices = np.empty(
+ (nnew_elements, grp.dim+1),
+ dtype=grp.vertex_indices.dtype)
+
+ new_nodes = np.empty(
+ (ambient_dim, nnew_elements, grp.nunit_nodes),
+ dtype=grp.nodes.dtype)
+
+ inew_vertex = vertex_nr_base
+ inew_el = 0
+ for iel in xrange(grp.nelements):
+ if not refine_flags[iel]:
+ new_nodes[:, inew_el, :] = grp.nodes[:, iel, :]
+
+ # old vertices always keep their numbers
+ new_vertex_indices[inew_el, :] = grp.vertex_indices[iel, :]
+ continue
+
+ el_vertex_indices = grp.vertex_indices[iel, :]
+
+
+
+
+
+
+
+
+
+
+
+
+ from meshmode.mesh import SimplexElementGroup
+ new_grp = SimplexElementGroup(grp.order, new_vertex_indices,
+ new_nodes, unit_nodes=grp.unit_nodes)
+
+ return new_vertices, new_grp
+
+
+def refine_without_connectivity(mesh, refine_flags):
+ vertex_chunks = [mesh.vertices]
+ nvertices = mesh.vertices.shape[1]
+
+ groups = []
+
+ from meshmode.mesh import SimplexElementGroup, Mesh
+ for grp in mesh.groups:
+ if isinstance(grp, SimplexElementGroup):
+ enb = grp.element_nr_base
+ added_vertices, new_grp = _refine_simplex_el_group(
+ grp, refine_flags[enb:enb+grp.nelements], nvertices)
+
+ vertex_chunks.append(added_vertices)
+ groups.append(new_grp)
+ else:
+ raise NotImplementedError("refining element group of type %s"
+ % type(grp).__name__)
+
+ vertices = np.hstack(vertex_chunks)
+
+ if all(grp.dim == 1 for grp in groups):
+ # For meshes that are exclusively 1D, deducing connectivity from vertex
+ # indices is still OK. For everyone else, not really.
+
+ connectivity = None
+ else:
+ connectivity = False
+
+ return Mesh(vertices, groups, element_connectivity=connectivity)
+
+# }}}
+
+
class _SplitFaceRecord(object):
"""
.. attribute:: neighboring_elements
diff --git a/meshmode/mesh/visualization.py b/meshmode/mesh/visualization.py
index 41f0917..6a2f921 100644
--- a/meshmode/mesh/visualization.py
+++ b/meshmode/mesh/visualization.py
@@ -30,7 +30,7 @@ import numpy as np
# {{{ draw_2d_mesh
def draw_2d_mesh(mesh, draw_vertex_numbers=True, draw_element_numbers=True,
- draw_connectivity=False, set_bounding_box=False, **kwargs):
+ draw_nodal_adjacency=False, set_bounding_box=False, **kwargs):
assert mesh.ambient_dim == 2
import matplotlib.pyplot as pt
@@ -74,7 +74,7 @@ def draw_2d_mesh(mesh, draw_vertex_numbers=True, draw_element_numbers=True,
ha="center", va="center", color="blue",
bbox=dict(facecolor='white', alpha=0.5, lw=0))
- if draw_connectivity:
+ if draw_nodal_adjacency:
def global_iel_to_group_and_iel(global_iel):
for igrp, grp in enumerate(mesh.groups):
if global_iel < grp.nelements:
@@ -83,7 +83,7 @@ def draw_2d_mesh(mesh, draw_vertex_numbers=True, draw_element_numbers=True,
raise ValueError("invalid element nr")
- cnx = mesh.element_connectivity
+ cnx = mesh.nodal_adjacency
nb_starts = cnx.neighbors_starts
for iel_g in range(mesh.nelements):
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 7c6e325..cfed94b 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -56,7 +56,7 @@ def test_circle_mesh(do_plot=False):
if do_plot:
from meshmode.mesh.visualization import draw_2d_mesh
- draw_2d_mesh(mesh, fill=None, draw_connectivity=True,
+ draw_2d_mesh(mesh, fill=None, draw_nodal_adjacency=True,
set_bounding_box=True)
import matplotlib.pyplot as pt
pt.show()
@@ -377,7 +377,7 @@ def test_rect_mesh(do_plot=False):
if do_plot:
from meshmode.mesh.visualization import draw_2d_mesh
- draw_2d_mesh(mesh, fill=None, draw_connectivity=True)
+ draw_2d_mesh(mesh, fill=None, draw_nodal_adjacency=True)
import matplotlib.pyplot as pt
pt.show()
@@ -405,7 +405,7 @@ def test_as_python():
from meshmode.mesh.generation import make_curve_mesh, cloverleaf
mesh = make_curve_mesh(cloverleaf, np.linspace(0, 1, 100), order=3)
- mesh.element_connectivity
+ mesh.nodal_adjacency
from meshmode.mesh import as_python
code = as_python(mesh)
--
GitLab