diff --git a/doc/images/facial-adjacency-group.png b/doc/images/facial-adjacency-group.png new file mode 100644 index 0000000000000000000000000000000000000000..a384da6432ecf50073ee30fb02430f7055a2cd6b Binary files /dev/null and b/doc/images/facial-adjacency-group.png differ diff --git a/doc/images/facial-adjacency-group.xoj b/doc/images/facial-adjacency-group.xoj new file mode 100644 index 0000000000000000000000000000000000000000..6b8318eb268469ef96cd3e218f91f0bff5150c00 Binary files /dev/null and b/doc/images/facial-adjacency-group.xoj differ diff --git a/doc/images/nodes-vertices.png b/doc/images/nodes-vertices.png new file mode 100644 index 0000000000000000000000000000000000000000..8ee30bd8e29ff657eb8f2a6c2193956c40f9969b Binary files /dev/null and b/doc/images/nodes-vertices.png differ diff --git a/doc/images/nodes-vertices.xoj b/doc/images/nodes-vertices.xoj new file mode 100644 index 0000000000000000000000000000000000000000..fc0bb2950f4eafb2f1adc226b0b23c2f52ba7fc9 Binary files /dev/null and b/doc/images/nodes-vertices.xoj differ diff --git a/doc/mesh.rst b/doc/mesh.rst index 28a332cf5186b7b95e91baff1d336d545ee47f88..2eb79cbe14869b9624cfd655f65130c3bc3ac451 100644 --- a/doc/mesh.rst +++ b/doc/mesh.rst @@ -8,6 +8,40 @@ Mesh management .. automodule:: meshmode.mesh +Design of the Data Structure +---------------------------- + +Why does a :class:`Mesh` need to be broken into :class:`MeshElementGroup` instances? +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Elements can be of different types (e.g. triangle, quadrilateral, +tetrahedron, what have you). In addition, elements may vary in the +polynomial degree used to represent them (see also below). + +All these bits of information could in principle be stored by element, +but having large, internally homogeneous groups is a good thing from an +efficiency standpoint. (So that you can, e.g., launch one GPU kernel to +deal with all order-3 triangles, instead of maybe having to dispatch +based on type and size inside the kernel) + +What is the difference between 'vertices' and 'nodes'? +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Nodes exist mainly to represent the (potentially non-affine) deformation of +each element, by a one-to-one correspondence with +:attr:`MeshElementGroup.unit_nodes`. They are unique to each element. Vertices +on the other hand exist to clarify whether or not a point shared by two +elements is actually identical (or just happens to be "close"). This is done by +assigning (single, globally shared) vertex numbers and having elements refer to +them. + +Consider the following picture: + +.. image:: images/nodes-vertices.png + +Mesh Data Structure +------------------- + .. autoclass:: MeshElementGroup :members: :undoc-members: diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py index 45fdc22e4b64f93dfd8c04f34d0e5d29cbecf7a6..6da392cae8b58022cf3efd72b3e710fd980b57f5 100644 --- a/meshmode/discretization/__init__.py +++ b/meshmode/discretization/__init__.py @@ -37,6 +37,10 @@ __doc__ = """ # {{{ element group base +class NoninterpolatoryElementGroupError(TypeError): + pass + + class ElementGroupBase(object): """Container for the :class:`Discretization` data corresponding to one :class:`meshmode.mesh.MeshElementGroup`. diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index ba15835834a1dd7afbc53bc8c35ff258d7778e42..bdb0f5250256ad6e8a67a10d9b735c3ff085da4b 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -56,6 +56,7 @@ __all__ = [ __doc__ = """ .. autoclass:: DiscretizationConnection +.. autoclass:: DirectDiscretizationConnection .. autofunction:: make_same_mesh_connection @@ -155,12 +156,12 @@ class DiscretizationConnectionElementGroup(object): # }}} -# {{{ connection class +# {{{ connection classes class DiscretizationConnection(object): - """Data supporting an interpolation-like operation that takes in data on - one discretization and returns it on another. Implemented applications - include: + """Abstract interface for transporting a DOF vector from one + :class:`meshmode.discretization.Discretization` to another. + Possible applications include: * upsampling/downsampling on the same mesh * restricition to the boundary @@ -171,24 +172,14 @@ class DiscretizationConnection(object): .. attribute:: to_discr - .. attribute:: groups - - a list of :class:`DiscretizationConnectionElementGroup` - instances, with a one-to-one correspondence to the groups in - :attr:`to_discr`. - .. attribute:: is_surjective A :class:`bool` indicating whether every output degree of freedom is set by the connection. .. automethod:: __call__ - - .. automethod:: full_resample_matrix - """ - - def __init__(self, from_discr, to_discr, groups, is_surjective): + def __init__(self, from_discr, to_discr, is_surjective): if from_discr.cl_context != to_discr.cl_context: raise ValueError("from_discr and to_discr must live in the " "same OpenCL context") @@ -203,20 +194,86 @@ class DiscretizationConnection(object): raise ValueError("from_discr and to_discr must agree on the " "element_id_dtype") - self.cl_context = from_discr.cl_context - self.from_discr = from_discr self.to_discr = to_discr - self.groups = groups self.is_surjective = is_surjective + def __call__(self, queue, vec): + raise NotImplementedError() + + +class ChainedDiscretizationConnection(DiscretizationConnection): + """Aggregates multiple :class:`DiscretizationConnection` instances + into a single one. + + .. attribute:: connections + """ + + def __init__(self, connections): + if not connections: + raise ValueError("connections may not be empty") + + super(DirectDiscretizationConnection, self).__init__( + connections[0].from_discr, + connections[-1].to_discr, + is_surjective=all( + cnx.is_surjective for cnx in connections)) + + self.connections = connections + + def __call__(self, queue, vec): + for cnx in self.connections: + vec = cnx(queue, vec) + + return vec + + +class DirectDiscretizationConnection(DiscretizationConnection): + """A concrete :class:`DiscretizationConnection` supported by interpolation + data. + + .. attribute:: from_discr + + .. attribute:: to_discr + + .. attribute:: groups + + a list of :class:`DiscretizationConnectionElementGroup` + instances, with a one-to-one correspondence to the groups in + :attr:`to_discr`. + + .. attribute:: is_surjective + + A :class:`bool` indicating whether every output degree + of freedom is set by the connection. + + .. automethod:: __call__ + + .. automethod:: full_resample_matrix + + """ + + def __init__(self, from_discr, to_discr, groups, is_surjective): + super(DirectDiscretizationConnection, self).__init__( + from_discr, to_discr, is_surjective) + + self.groups = groups + @memoize_method def _resample_matrix(self, to_group_index, ibatch_index): import modepy as mp ibatch = self.groups[to_group_index].batches[ibatch_index] from_grp = self.from_discr.groups[ibatch.from_group_index] + if len(from_grp.basis()) != from_grp.unit_nodes.shape[1]: + from meshmode.discretization import NoninterpolatoryElementGroupError + raise NoninterpolatoryElementGroupError( + "%s does not support interpolation because it is not " + "unisolvent (its unit node count does not match its " + "number of basis functions). Using connections requires " + "the ability to interpolate." % type(from_grp).__name__) + result = mp.resampling_matrix( from_grp.basis(), ibatch.result_unit_nodes, from_grp.unit_nodes) diff --git a/meshmode/discretization/connection/face.py b/meshmode/discretization/connection/face.py index 124b4252fbba1f62fe1a89da509f062c2e6e8f4d..88d91d9e615a6102e9e891c76f5948186a7b5f6b 100644 --- a/meshmode/discretization/connection/face.py +++ b/meshmode/discretization/connection/face.py @@ -62,7 +62,7 @@ def _build_boundary_connection(queue, vol_discr, bdry_discr, connection_data, per_face_groups): from meshmode.discretization.connection import ( InterpolationBatch, DiscretizationConnectionElementGroup, - DiscretizationConnection) + DirectDiscretizationConnection) ibdry_grp = 0 batches = [] @@ -105,7 +105,7 @@ def _build_boundary_connection(queue, vol_discr, bdry_discr, connection_data, assert ibdry_grp == len(bdry_discr.groups) - return DiscretizationConnection( + return DirectDiscretizationConnection( vol_discr, bdry_discr, connection_groups, is_surjective=True) @@ -180,10 +180,11 @@ def make_face_restriction(discr, group_factory, boundary_tag, each other one-to-one, and an interpolation batch is created per face. - :return: a :class:`meshmode.discretization.connection.DiscretizationConnection` + :return: a + :class:`meshmode.discretization.connection.DirectDiscretizationConnection` representing the new connection. The new boundary discretization can be obtained from the - :attr:`meshmode.discretization.connection.DiscretizationConnection.to_discr` + :attr:`meshmode.discretization.connection.DirectDiscretizationConnection.to_discr` attribute of the return value, and the corresponding new boundary mesh from that. @@ -415,7 +416,7 @@ def make_face_to_all_faces_embedding(faces_connection, all_faces_discr): "same number of groups") from meshmode.discretization.connection import ( - DiscretizationConnection, + DirectDiscretizationConnection, DiscretizationConnectionElementGroup, InterpolationBatch) @@ -477,7 +478,7 @@ def make_face_to_all_faces_embedding(faces_connection, all_faces_discr): i_faces_grp += 1 - return DiscretizationConnection( + return DirectDiscretizationConnection( faces_connection.to_discr, all_faces_discr, groups, diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py index 70728522b2c06f119a5b20598616e2c464eea249..6ce70b2a11b148ce6a1e2cca6f7588ac2549aada 100644 --- a/meshmode/discretization/connection/opposite_face.py +++ b/meshmode/discretization/connection/opposite_face.py @@ -334,7 +334,7 @@ def _make_el_lookup_table(queue, connection, igrp): def make_opposite_face_connection(volume_to_bdry_conn): """Given a boundary restriction connection *volume_to_bdry_conn*, - return a :class:`DiscretizationConnection` that performs data + return a :class:`DirectDiscretizationConnection` that performs data exchange across opposite faces. """ @@ -381,8 +381,8 @@ def make_opposite_face_connection(volume_to_bdry_conn): vbc_tgt_grp_face_batch, src_grp_el_lookup)) from meshmode.discretization.connection import ( - DiscretizationConnection, DiscretizationConnectionElementGroup) - return DiscretizationConnection( + DirectDiscretizationConnection, DiscretizationConnectionElementGroup) + return DirectDiscretizationConnection( from_discr=bdry_discr, to_discr=bdry_discr, groups=[ diff --git a/meshmode/discretization/connection/refinement.py b/meshmode/discretization/connection/refinement.py index ff009fce6efe7833e7a8bbd43c6d33745e722704..f487b2bf7281fa3e793fbd924e2dd65ab608ab79 100644 --- a/meshmode/discretization/connection/refinement.py +++ b/meshmode/discretization/connection/refinement.py @@ -127,7 +127,8 @@ def make_refinement_connection(refiner, coarse_discr, group_factory): discretizing the fine mesh. """ from meshmode.discretization.connection import ( - DiscretizationConnectionElementGroup, DiscretizationConnection) + DiscretizationConnectionElementGroup, + DirectDiscretizationConnection) coarse_mesh = refiner.get_previous_mesh() fine_mesh = refiner.last_mesh @@ -157,7 +158,7 @@ def make_refinement_connection(refiner, coarse_discr, group_factory): logger.info("building refinement connection: done") - return DiscretizationConnection( + return DirectDiscretizationConnection( from_discr=coarse_discr, to_discr=fine_discr, groups=groups, diff --git a/meshmode/discretization/connection/same_mesh.py b/meshmode/discretization/connection/same_mesh.py index 3b73e34fefd5115e01b97728ab414cc8c77fd178..a470e9298238dbca6874176045ee44425806413c 100644 --- a/meshmode/discretization/connection/same_mesh.py +++ b/meshmode/discretization/connection/same_mesh.py @@ -32,7 +32,7 @@ import pyopencl.array # noqa def make_same_mesh_connection(to_discr, from_discr): from meshmode.discretization.connection import ( InterpolationBatch, DiscretizationConnectionElementGroup, - DiscretizationConnection) + DirectDiscretizationConnection) if from_discr.mesh is not to_discr.mesh: raise ValueError("from_discr and to_discr must be based on " @@ -56,7 +56,7 @@ def make_same_mesh_connection(to_discr, from_discr): groups.append( DiscretizationConnectionElementGroup([ibatch])) - return DiscretizationConnection( + return DirectDiscretizationConnection( from_discr, to_discr, groups, is_surjective=True) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 54aaabef97ac5f9bfc97b642ebe1a476608e0603..5f636e78165f46b32ecb9c41e2883927ece81f82 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -26,7 +26,9 @@ THE SOFTWARE. import numpy as np #import numpy.linalg as la from pytools import memoize_method -from meshmode.mesh import SimplexElementGroup as _MeshSimplexElementGroup +from meshmode.mesh import ( + SimplexElementGroup as _MeshSimplexElementGroup, + TensorProductElementGroup as _MeshTensorProductElementGroup) import modepy as mp @@ -38,7 +40,8 @@ Group types .. autoclass:: InterpolatoryQuadratureSimplexElementGroup .. autoclass:: QuadratureSimplexElementGroup .. autoclass:: PolynomialWarpAndBlendElementGroup -.. autoclass:: PolynomialEquidistantElementGroup +.. autoclass:: PolynomialEquidistantSimplexElementGroup +.. autoclass:: LegendreGaussLobattoTensorProductElementGroup Group factories ^^^^^^^^^^^^^^^ @@ -46,7 +49,8 @@ Group factories .. autoclass:: InterpolatoryQuadratureSimplexGroupFactory .. autoclass:: QuadratureSimplexGroupFactory .. autoclass:: PolynomialWarpAndBlendGroupFactory -.. autoclass:: PolynomialEquidistantGroupFactory +.. autoclass:: PolynomialEquidistantSimplexGroupFactory +.. autoclass:: LegendreGaussLobattoTensorProductGroupFactory """ # FIXME Most of the loopy kernels will break as soon as we start using multiple @@ -61,25 +65,28 @@ Group factories from meshmode.discretization import ElementGroupBase -# {{{ concrete element groups - -class PolynomialSimplexElementGroupBase(ElementGroupBase): - def basis(self): - return mp.simplex_best_available_basis(self.dim, self.order) - - def grad_basis(self): - return mp.grad_simplex_best_available_basis(self.dim, self.order) +# {{{ base class for poynomial elements +class PolynomialElementGroupBase(ElementGroupBase): @memoize_method - def from_mesh_interp_matrix(self): - meg = self.mesh_el_group - return mp.resampling_matrix( - mp.simplex_best_available_basis(meg.dim, meg.order), - self.unit_nodes, - meg.unit_nodes) + def mass_matrix(self): + assert self.is_orthogonal_basis() + + import modepy as mp + return mp.mass_matrix( + self.basis(), + self.unit_nodes) @memoize_method def diff_matrices(self): + if len(self.basis()) != self.unit_nodes.shape[1]: + from meshmode.discretization import NoninterpolatoryElementGroupError + raise NoninterpolatoryElementGroupError( + "%s does not support interpolation because it is not " + "unisolvent (its unit node count does not match its " + "number of basis functions). Differentiation requires " + "the ability to interpolate." % type(self).__name__) + result = mp.differentiation_matrices( self.basis(), self.grad_basis(), @@ -90,6 +97,35 @@ class PolynomialSimplexElementGroupBase(ElementGroupBase): else: return result +# }}} + + +# {{{ concrete element groups for simplices + +class PolynomialSimplexElementGroupBase(PolynomialElementGroupBase): + def is_orthogonal_basis(self): + return self.dim <= 3 + + def basis(self): + if self.dim <= 3: + return mp.simplex_onb(self.dim, self.order) + else: + return mp.simplex_monomial_basis(self.dim, self.order) + + def grad_basis(self): + if self.dim <= 3: + return mp.grad_simplex_onb(self.dim, self.order) + else: + return mp.grad_simplex_monomial_basis(self.dim, self.order) + + @memoize_method + def from_mesh_interp_matrix(self): + meg = self.mesh_el_group + return mp.resampling_matrix( + mp.simplex_best_available_basis(meg.dim, meg.order), + self.unit_nodes, + meg.unit_nodes) + class InterpolatoryQuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase): """Elemental discretization supplying a high-order quadrature rule @@ -163,7 +199,7 @@ class _MassMatrixQuadratureElementGroup(PolynomialSimplexElementGroupBase): @memoize_method def weights(self): return np.dot( - mp.mass_matrix(self.basis(), self.unit_nodes), + self.mass_matrix(), np.ones(len(self.basis()))) @@ -184,7 +220,7 @@ class PolynomialWarpAndBlendElementGroup(_MassMatrixQuadratureElementGroup): return result -class PolynomialEquidistantElementGroup(_MassMatrixQuadratureElementGroup): +class PolynomialEquidistantSimplexElementGroup(_MassMatrixQuadratureElementGroup): """Elemental discretization with a number of nodes matching the number of polynomials in :math:`P^k`, hence usable for differentiation and interpolation. Interpolation nodes are present on the boundary of the @@ -205,13 +241,54 @@ class PolynomialEquidistantElementGroup(_MassMatrixQuadratureElementGroup): # }}} +# {{{ concrete element groups for tensor product elements + +class LegendreGaussLobattoTensorProductElementGroup(PolynomialElementGroupBase): + def is_orthogonal_basis(self): + return True + + def basis(self): + from modepy.modes import tensor_product_basis, jacobi + from functools import partial + return tensor_product_basis( + self.dim, tuple( + partial(jacobi, 0, 0, i) + for i in range(self.order+1))) + + def grad_basis(self): + raise NotImplementedError() + + @property + @memoize_method + def unit_nodes(self): + from modepy.nodes import tensor_product_nodes + from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes + return tensor_product_nodes( + self.dim, legendre_gauss_lobatto_nodes(self.order)) + + @memoize_method + def from_mesh_interp_matrix(self): + from modepy.modes import tensor_product_basis, jacobi + from functools import partial + meg = self.mesh_el_group + + basis = tensor_product_basis( + self.dim, tuple( + partial(jacobi, 0, 0, i) + for i in range(meg.order+1))) + + return mp.resampling_matrix(basis, self.unit_nodes, meg.unit_nodes) + +# }}} + + # {{{ group factories class ElementGroupFactory(object): pass -class OrderBasedGroupFactory(ElementGroupFactory): +class HomogeneousOrderBasedGroupFactory(ElementGroupFactory): def __init__(self, order): self.order = order @@ -223,28 +300,58 @@ class OrderBasedGroupFactory(ElementGroupFactory): return self.group_class(mesh_el_group, self.order, node_nr_base) -class InterpolatoryQuadratureSimplexGroupFactory(OrderBasedGroupFactory): +class OrderAndTypeBasedGroupFactory(ElementGroupFactory): + def __init__(self, order, simplex_group_class, tensor_product_group_class): + self.order = order + self.simplex_group_class = simplex_group_class + self.tensor_product_group_class = tensor_product_group_class + + def __call__(self, mesh_el_group, node_nr_base): + if isinstance(mesh_el_group, _MeshSimplexElementGroup): + group_class = self.simplex_group_class + elif isinstance(mesh_el_group, _MeshTensorProductElementGroup): + group_class = self.tensor_product_group_class + else: + raise TypeError("only mesh element groups of type '%s' and '%s' " + "are supported" % ( + _MeshSimplexElementGroup.__name__, + _MeshTensorProductElementGroup.__name__)) + + return group_class(mesh_el_group, self.order, node_nr_base) + + +# {{{ group factories for simplices + +class InterpolatoryQuadratureSimplexGroupFactory(HomogeneousOrderBasedGroupFactory): mesh_group_class = _MeshSimplexElementGroup group_class = InterpolatoryQuadratureSimplexElementGroup -class QuadratureSimplexGroupFactory(OrderBasedGroupFactory): +class QuadratureSimplexGroupFactory(HomogeneousOrderBasedGroupFactory): mesh_group_class = _MeshSimplexElementGroup group_class = QuadratureSimplexElementGroup -class PolynomialWarpAndBlendGroupFactory(OrderBasedGroupFactory): +class PolynomialWarpAndBlendGroupFactory(HomogeneousOrderBasedGroupFactory): mesh_group_class = _MeshSimplexElementGroup group_class = PolynomialWarpAndBlendElementGroup -class PolynomialEquidistantGroupFactory(OrderBasedGroupFactory): +class PolynomialEquidistantSimplexGroupFactory(HomogeneousOrderBasedGroupFactory): """ .. versionadded:: 2016.1 """ mesh_group_class = _MeshSimplexElementGroup - group_class = PolynomialEquidistantElementGroup + group_class = PolynomialEquidistantSimplexElementGroup + +# }}} + + +class LegendreGaussLobattoTensorProductGroupFactory( + HomogeneousOrderBasedGroupFactory): + mesh_group_class = _MeshTensorProductElementGroup + group_class = LegendreGaussLobattoTensorProductElementGroup # }}} diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py index d68ec0c58078e68417721c8205fc3f24021d0bad..1021d8f4f7556e5306fdd74703dcb94d531b4c5c 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -1,6 +1,4 @@ -from __future__ import division -from __future__ import absolute_import -from six.moves import range +from __future__ import division, absolute_import __copyright__ = "Copyright (C) 2014 Andreas Kloeckner" @@ -24,8 +22,9 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +from six.moves import range import numpy as np -from pytools import memoize_method +from pytools import memoize_method, Record import pyopencl as cl __doc__ = """ @@ -71,17 +70,50 @@ def separate_by_real_and_imag(data, real_only): yield (name, field) +class _VisConnectivityGroup(Record): + """ + .. attribute:: vis_connectivity + + an array of shape ``(group.nelements,nsubelements,primitive_element_size)`` + + .. attribute:: vtk_cell_type + + .. attribute:: subelement_nr_base + """ + + @property + def nsubelements(self): + return self.nelements * self.nsubelements_per_element + + @property + def nelements(self): + return self.vis_connectivity.shape[0] + + @property + def nsubelements_per_element(self): + return self.vis_connectivity.shape[1] + + @property + def primitive_element_size(self): + return self.vis_connectivity.shape[2] + + class Visualizer(object): """ .. automethod:: show_scalar_in_mayavi .. automethod:: write_vtk_file """ - def __init__(self, connection): + def __init__(self, connection, element_shrink_factor=None): self.connection = connection self.discr = connection.from_discr self.vis_discr = connection.to_discr + if element_shrink_factor is None: + element_shrink_factor = 1 + + self.element_shrink_factor = element_shrink_factor + def _resample_and_get(self, queue, vec): from pytools.obj_array import with_object_array_or_scalar @@ -90,42 +122,101 @@ class Visualizer(object): return with_object_array_or_scalar(resample_and_get_one, vec) + # {{{ vis sub-element connectivity + @memoize_method def _vis_connectivity(self): """ - :return: an array of shape - ``(vis_discr.nelements,nsubelements,primitive_element_size)`` + :return: a list of :class:`_VisConnectivityGroup` instances. """ # Assume that we're using modepy's default node ordering. - from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ - as gnitstam, single_valued - vis_order = single_valued( - group.order for group in self.vis_discr.groups) - node_tuples = list(gnitstam(vis_order, self.vis_discr.dim)) + from pytools import ( + generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam, + generate_nonnegative_integer_tuples_below as gnitb) + from meshmode.mesh import TensorProductElementGroup, SimplexElementGroup + + result = [] - from modepy.tools import submesh - el_connectivity = np.array( - submesh(node_tuples), - dtype=np.intp) + from pyvisfile.vtk import ( + VTK_LINE, VTK_TRIANGLE, VTK_TETRA, + VTK_QUAD, VTK_HEXAHEDRON) - nelements = sum(group.nelements for group in self.vis_discr.groups) - vis_connectivity = np.empty( - (nelements,) + el_connectivity.shape, dtype=np.intp) + subel_nr_base = 0 - el_nr_base = 0 for group in self.vis_discr.groups: + if isinstance(group.mesh_el_group, SimplexElementGroup): + node_tuples = list(gnitstam(group.order, group.dim)) + + from modepy.tools import submesh + el_connectivity = np.array( + submesh(node_tuples), + dtype=np.intp) + vtk_cell_type = { + 1: VTK_LINE, + 2: VTK_TRIANGLE, + 3: VTK_TETRA, + }[group.dim] + + elif isinstance(group.mesh_el_group, TensorProductElementGroup): + node_tuples = list(gnitb(group.order+1, group.dim)) + node_tuple_to_index = dict( + (nt, i) for i, nt in enumerate(node_tuples)) + + def add_tuple(a, b): + return tuple(ai+bi for ai, bi in zip(a, b)) + + el_offsets = { + 1: [(0,), (1,)], + 2: [(0, 0), (1, 0), (1, 1), (0, 1)], + 3: [ + (0, 0, 0), + (1, 0, 0), + (1, 1, 0), + (0, 1, 0), + (0, 0, 1), + (1, 0, 1), + (1, 1, 1), + (0, 1, 1), + ] + }[group.dim] + + el_connectivity = np.array([ + [ + node_tuple_to_index[add_tuple(origin, offset)] + for offset in el_offsets] + for origin in gnitb(group.order, group.dim)]) + + vtk_cell_type = { + 1: VTK_LINE, + 2: VTK_QUAD, + 3: VTK_HEXAHEDRON, + }[group.dim] + + else: + raise NotImplementedError("visualization for element groups " + "of type '%s'" % type(group.mesh_el_group).__name__) + assert len(node_tuples) == group.nunit_nodes - vis_connectivity[el_nr_base:el_nr_base+group.nelements] = ( - np.arange( - el_nr_base*group.nunit_nodes, - (el_nr_base+group.nelements)*group.nunit_nodes, - group.nunit_nodes + vis_connectivity = ( + group.node_nr_base + np.arange( + 0, group.nelements*group.nunit_nodes, group.nunit_nodes )[:, np.newaxis, np.newaxis] - + el_connectivity) - el_nr_base += group.nelements + + el_connectivity).astype(np.intp) + + vgrp = _VisConnectivityGroup( + vis_connectivity=vis_connectivity, + vtk_cell_type=vtk_cell_type, + subelement_nr_base=subel_nr_base) + result.append(vgrp) + + subel_nr_base += vgrp.nsubelements + + return result + + # }}} - return vis_connectivity + # {{{ mayavi def show_scalar_in_mayavi(self, field, **kwargs): import mayavi.mlab as mlab @@ -140,7 +231,7 @@ class Visualizer(object): assert nodes.shape[0] == self.vis_discr.ambient_dim #mlab.points3d(nodes[0], nodes[1], 0*nodes[0]) - vis_connectivity = self._vis_connectivity() + vis_connectivity, = self._vis_connectivity() if self.vis_discr.dim == 1: nodes = list(nodes) @@ -153,6 +244,7 @@ class Visualizer(object): # http://docs.enthought.com/mayavi/mayavi/auto/example_plotting_many_lines.html # noqa src = mlab.pipeline.scalar_scatter(*args) + src.mlab_source.dataset.lines = vis_connectivity.reshape(-1, 2) lines = mlab.pipeline.stripper(src) mlab.pipeline.surface(lines, **kwargs) @@ -175,21 +267,17 @@ class Visualizer(object): if do_show: mlab.show() + # }}} + + # {{{ vtk + def write_vtk_file(self, file_name, names_and_fields, compressor=None, real_only=False): from pyvisfile.vtk import ( UnstructuredGrid, DataArray, AppendedDataXMLGenerator, - VTK_LINE, VTK_TRIANGLE, VTK_TETRA, VF_LIST_OF_COMPONENTS) - el_types = { - 1: VTK_LINE, - 2: VTK_TRIANGLE, - 3: VTK_TETRA, - } - - el_type = el_types[self.vis_discr.dim] with cl.CommandQueue(self.vis_discr.cl_context) as queue: nodes = self.vis_discr.nodes().with_queue(queue).get() @@ -198,20 +286,41 @@ class Visualizer(object): (name, self._resample_and_get(queue, fld)) for name, fld in names_and_fields] - connectivity = self._vis_connectivity() + vc_groups = self._vis_connectivity() + + # {{{ create cell_types + + nsubelements = sum(vgrp.nsubelements for vgrp in vc_groups) + cell_types = np.empty(nsubelements, dtype=np.uint8) + cell_types.fill(255) + for vgrp in vc_groups: + cell_types[ + vgrp.subelement_nr_base: + vgrp.subelement_nr_base + vgrp.nsubelements] = \ + vgrp.vtk_cell_type + assert (cell_types < 255).all() + + # }}} - nprimitive_elements = ( - connectivity.shape[0] - * connectivity.shape[1]) + nodes = nodes + if self.element_shrink_factor != 1: + for vgrp in self.vis_discr.groups: + nodes_view = vgrp.view(nodes) + el_centers = np.mean(nodes_view, axis=-1) + nodes_view[:] = ( + (self.element_shrink_factor * nodes_view) + + (1-self.element_shrink_factor) + * el_centers[:, :, np.newaxis]) grid = UnstructuredGrid( (self.vis_discr.nnodes, DataArray("points", nodes.reshape(self.vis_discr.ambient_dim, -1), vector_format=VF_LIST_OF_COMPONENTS)), - cells=connectivity.reshape(-1), - cell_types=np.asarray([el_type] * nprimitive_elements, - dtype=np.uint8)) + cells=np.hstack([ + vgrp.vis_connectivity.reshape(-1) + for vgrp in vc_groups]), + cell_types=cell_types) # for name, field in separate_by_real_and_imag(cell_data, real_only): # grid.add_celldata(DataArray(name, field, @@ -230,18 +339,26 @@ class Visualizer(object): AppendedDataXMLGenerator(compressor)(grid).write(outf) -def make_visualizer(queue, discr, vis_order): +def make_visualizer(queue, discr, vis_order, element_shrink_factor=None): from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import \ - PolynomialWarpAndBlendGroupFactory + from meshmode.discretization.poly_element import ( + PolynomialWarpAndBlendElementGroup, + LegendreGaussLobattoTensorProductElementGroup, + OrderAndTypeBasedGroupFactory) vis_discr = Discretization( discr.cl_context, discr.mesh, - PolynomialWarpAndBlendGroupFactory(vis_order), + OrderAndTypeBasedGroupFactory( + vis_order, + simplex_group_class=PolynomialWarpAndBlendElementGroup, + tensor_product_group_class=( + LegendreGaussLobattoTensorProductElementGroup)), real_dtype=discr.real_dtype) from meshmode.discretization.connection import \ make_same_mesh_connection - return Visualizer(make_same_mesh_connection(vis_discr, discr)) + return Visualizer( + make_same_mesh_connection(vis_discr, discr), + element_shrink_factor=element_shrink_factor) # }}} diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 3f02e7f7162245cdc9b3911aa0a5ae0c91da987e..d3c9c8cfd11e4040ec7d507354776b9682f0fcae 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -1,6 +1,4 @@ -from __future__ import division -from __future__ import absolute_import -from six.moves import range +from __future__ import division, absolute_import __copyright__ = "Copyright (C) 2010,2012,2013 Andreas Kloeckner, Michael Tom" @@ -24,7 +22,9 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +from six.moves import range import six + import numpy as np import modepy as mp import numpy.linalg as la @@ -175,7 +175,9 @@ class MeshElementGroup(Record): def join_mesh(self, element_nr_base, node_nr_base): if self.element_nr_base is not None: raise RuntimeError("this element group has already joined a mesh, " - "cannot join another") + "cannot join another (The element group's element_nr_base " + "is already assigned, and that typically happens when a " + "group joins a Mesh instance.)") return self.copy( element_nr_base=element_nr_base, @@ -326,7 +328,7 @@ class TensorProductElementGroup(MeshElementGroup): if vertex_indices.shape[-1] != 2**dims: raise ValueError("vertex_indices has wrong number of vertices per " - "element. expected: %d, got: %d" % (dims+1, + "element. expected: %d, got: %d" % (2**dims, vertex_indices.shape[-1])) super(TensorProductElementGroup, self).__init__(order, vertex_indices, nodes, @@ -387,6 +389,13 @@ class FacialAdjacencyGroup(Record): :class:`MeshElementGroup`, i.e. information about elements that share (part of) a face. + .. image:: images/facial-adjacency-group.png + + Represents (for example) *one* of the (colored) interfaces between + :class:`MeshElementGroup` instances, or an interface between + :class:`MeshElementGroup` and a boundary. (Note that element groups are not + necessarily contiguous like the figure may suggest.) + .. attribute:: igroup .. attribute:: ineighbor_group @@ -404,7 +413,7 @@ class FacialAdjacencyGroup(Record): ``face_id_t [nfagrp_elements]``. ``element_faces[i]`` indicate what face index of the opposite element indicated in - ``neighbors[iel_grp][iface]`` touches face number *iface* of element + ``neighbors[iface]`` touches face number *iface* of element number *iel_grp* in this element group. .. attribute:: neighbors @@ -481,6 +490,25 @@ class Mesh(Record): Referencing this attribute may raise :exc:`meshmode.DataUnavailable`. + .. image:: images/facial-adjacency-group.png + + For example for the mesh in the figure, the following data structure + would be present:: + + [ + {...}, # connectivity for group 0 + {...}, # connectivity for group 1 + {...}, # connectivity for group 2 + { # connectivity for group 3 + 1: FacialAdjacencyGroup(...) # towards group 1, green + 2: FacialAdjacencyGroup(...) # towards group 2, pink + None: FacialAdjacencyGroup(...) # towards the boundary, orange + } + ] + + (Note that element groups are not necessarily contiguous like the figure + may suggest.) + .. attribute:: boundary_tags A tuple of boundary tag identifiers. :class:`BTAG_ALL` and @@ -609,19 +637,21 @@ class Mesh(Record): assert len(facial_adjacency_groups) == len(self.groups) for fagrp_map in facial_adjacency_groups: for fagrp in six.itervalues(fagrp_map): - grp = self.groups[fagrp.igroup] + nfagrp_elements, = fagrp.elements.shape + + assert fagrp.element_faces.dtype == self.face_id_dtype + assert fagrp.element_faces.shape == (nfagrp_elements,) - fvi = grp.face_vertex_indices() assert fagrp.neighbors.dtype == self.element_id_dtype - assert fagrp.neighbors.shape == ( - grp.nelements, len(fvi)) - assert fagrp.neighbors.dtype == self.face_id_dtype - assert fagrp.neighbor_faces.shape == ( - grp.nelements, len(fvi)) - - is_bdry = fagrp.neighbors < 0 - assert ((1 << btag_to_index[BTAG_REALLY_ALL]) - & fagrp.neighbors[is_bdry]).all(), \ + assert fagrp.neighbors.shape == (nfagrp_elements,) + + assert fagrp.neighbor_faces.dtype == self.face_id_dtype + assert fagrp.neighbor_faces.shape == (nfagrp_elements,) + + if fagrp.ineighbor_group is None: + is_bdry = fagrp.neighbors < 0 + assert ((1 << btag_to_index[BTAG_REALLY_ALL]) + & -fagrp.neighbors[is_bdry]).all(), \ "boundary faces without BTAG_REALLY_ALL found" from meshmode.mesh.processing import \ @@ -632,6 +662,24 @@ class Mesh(Record): assert test_volume_mesh_element_orientations(self), \ "negatively oriented elements found" + def get_copy_kwargs(self, **kwargs): + def set_if_not_present(name, from_name=None): + if from_name is None: + from_name = name + if name not in kwargs: + kwargs[name] = getattr(self, from_name) + + set_if_not_present("vertices") + if "groups" not in kwargs: + kwargs["groups"] = [group.copy() for group in self.groups] + set_if_not_present("boundary_tags") + set_if_not_present("nodal_adjacency", "_nodal_adjacency") + set_if_not_present("facial_adjacency_groups", "_facial_adjacency_groups") + set_if_not_present("vertex_id_dtype") + set_if_not_present("element_id_dtype") + + return kwargs + @property def ambient_dim(self): return self.vertices.shape[0] @@ -822,6 +870,9 @@ def _compute_facial_adjacency_from_vertices(mesh): face_map.setdefault( frozenset(fvi), []).append((igrp, iel_grp, fid)) + del igrp + del grp + # maps tuples (igrp, ineighbor_group) to number of elements group_count = {} for face_tuples in six.itervalues(face_map): @@ -835,6 +886,9 @@ def _compute_facial_adjacency_from_vertices(mesh): else: raise RuntimeError("unexpected number of adjacent faces") + del face_tuples + del igrp + # {{{ build facial_adjacency_groups data structure, still empty facial_adjacency_groups = [] @@ -849,6 +903,11 @@ def _compute_facial_adjacency_from_vertices(mesh): neighbors = np.empty(bdry_count, dtype=mesh.element_id_dtype) neighbor_faces = np.zeros(bdry_count, dtype=mesh.face_id_dtype) + # Ensure uninitialized entries get noticed + elements.fill(-1) + element_faces.fill(-1) + neighbor_faces.fill(-1) + neighbors.fill(-( mesh.boundary_tag_bit(BTAG_ALL) | mesh.boundary_tag_bit(BTAG_REALLY_ALL))) @@ -869,6 +928,12 @@ def _compute_facial_adjacency_from_vertices(mesh): neighbors = np.empty(nb_count, dtype=mesh.element_id_dtype) neighbor_faces = np.empty(nb_count, dtype=mesh.face_id_dtype) + # Ensure uninitialized entries get noticed + elements.fill(-1) + element_faces.fill(-1) + neighbors.fill(-1) + neighbor_faces.fill(-1) + grp_map[ineighbor_group] = FacialAdjacencyGroup( igroup=igroup, ineighbor_group=ineighbor_group, @@ -877,20 +942,24 @@ def _compute_facial_adjacency_from_vertices(mesh): neighbors=neighbors, neighbor_faces=neighbor_faces) + del igroup + del ineighbor_group + del grp_map + # }}} # maps tuples (igrp, ineighbor_group) to number of elements filled in group fill_count = {} for face_tuples in six.itervalues(face_map): if len(face_tuples) == 2: - for (igroup, iel, iface), (inb_group, inb_el, inb_face) in [ + for (igroup, iel, iface), (ineighbor_group, inb_el, inb_face) in [ (face_tuples[0], face_tuples[1]), (face_tuples[1], face_tuples[0]), ]: - idx = fill_count.get((igrp, inb_grp), 0) - fill_count[igrp, inb_grp] = idx + 1 + idx = fill_count.get((igroup, ineighbor_group), 0) + fill_count[igroup, ineighbor_group] = idx + 1 - fagrp = facial_adjacency_groups[igroup][inb_grp] + fagrp = facial_adjacency_groups[igroup][ineighbor_group] fagrp.elements[idx] = iel fagrp.element_faces[idx] = iface fagrp.neighbors[idx] = inb_el @@ -899,8 +968,8 @@ def _compute_facial_adjacency_from_vertices(mesh): elif len(face_tuples) == 1: (igroup, iel, iface), = face_tuples - idx = fill_count.get((igrp, None), 0) - fill_count[igrp, None] = idx + 1 + idx = fill_count.get((igroup, None), 0) + fill_count[igroup, None] = idx + 1 fagrp = facial_adjacency_groups[igroup][None] fagrp.elements[idx] = iel diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index ff5d8581d5d74291e48ecafe5936548b290e0960..3cf18d71eeca11b0d6279895ee2863b2e3433cc7 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -237,34 +237,88 @@ def make_curve_mesh(curve_f, element_boundaries, order, # {{{ make_group_from_vertices -def make_group_from_vertices(vertices, vertex_indices, order): +def make_group_from_vertices(vertices, vertex_indices, order, + group_factory=None): + # shape: (dim, nelements, nvertices) el_vertices = vertices[:, vertex_indices] - el_origins = el_vertices[:, :, 0][:, :, np.newaxis] - # ambient_dim, nelements, nspan_vectors - spanning_vectors = ( - el_vertices[:, :, 1:] - el_origins) + from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup - nspan_vectors = spanning_vectors.shape[-1] - dim = nspan_vectors + if group_factory is None: + group_factory = SimplexElementGroup - # dim, nunit_nodes - if dim <= 3: - unit_nodes = mp.warp_and_blend_nodes(dim, order) - else: - unit_nodes = mp.equidistant_nodes(dim, order) + if issubclass(group_factory, SimplexElementGroup): + 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 + if dim <= 3: + unit_nodes = mp.warp_and_blend_nodes(dim, order) + else: + unit_nodes = mp.equidistant_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 + + elif issubclass(group_factory, TensorProductElementGroup): + nelements, nvertices = vertex_indices.shape + + dim = 0 + while True: + if nvertices == 2**dim: + break + if nvertices < 2**dim: + raise ValueError("invalid number of vertices for tensor-product " + "elements, must be power of two") + dim += 1 + + from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes + from modepy.nodes import tensor_product_nodes + unit_nodes = tensor_product_nodes(dim, legendre_gauss_lobatto_nodes(order)) + # shape: (dim, nnodes) + unit_nodes_01 = 0.5 + 0.5*unit_nodes + + _, nnodes = unit_nodes.shape + + from pytools import generate_nonnegative_integer_tuples_below as gnitb + id_tuples = list(gnitb(2, dim)) + assert len(id_tuples) == nvertices - unit_nodes_01 = 0.5 + 0.5*unit_nodes + vdm = np.empty((nvertices, nvertices)) + for i, vertex_tuple in enumerate(id_tuples): + for j, func_tuple in enumerate(id_tuples): + vertex_ref = np.array(vertex_tuple, dtype=np.float64) + vdm[i, j] = np.prod(vertex_ref**func_tuple) - nodes = np.einsum( - "si,des->dei", - unit_nodes_01, spanning_vectors) + el_origins + # shape: (dim, nelements, nvertices) + coeffs = np.empty((dim, nelements, nvertices)) + for d in range(dim): + coeffs[d] = la.solve(vdm, el_vertices[d].T).T + + vdm_nodes = np.zeros((nnodes, nvertices)) + for j, func_tuple in enumerate(id_tuples): + vdm_nodes[:, j] = np.prod( + unit_nodes_01 ** np.array(func_tuple).reshape(-1, 1), + axis=0) + + nodes = np.einsum("ij,dej->dei", vdm_nodes, coeffs) + + else: + raise ValueError("unsupported value for 'group_factory': %s" + % group_factory) # make contiguous nodes = nodes.copy() - from meshmode.mesh import SimplexElementGroup - return SimplexElementGroup( + return group_factory( order, vertex_indices, nodes, unit_nodes=unit_nodes) @@ -409,12 +463,19 @@ def generate_torus(r_outer, r_inner, n_outer=20, n_inner=10, order=1): # {{{ generate_box_mesh -def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64): +def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64, + group_factory=None): """Create a semi-structured mesh. :param axis_coords: a tuple with a number of entries corresponding to the number of dimensions, with each entry a numpy array specifying the coordinates to be used along that axis. + :param group_factory: One of :class:`meshmode.mesh.SimplexElementGroup` + or :class:`meshmode.mesh.TensorProductElementGroup`. + + .. versionchanged:: 2017.1 + + *group_factory* parameter added. """ for iaxis, axc in enumerate(axis_coords): @@ -438,6 +499,18 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64): vertices = vertices.reshape(dim, -1) + from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup + if group_factory is None: + group_factory = SimplexElementGroup + + if issubclass(group_factory, SimplexElementGroup): + is_tp = False + elif issubclass(group_factory, TensorProductElementGroup): + is_tp = True + else: + raise ValueError("unsupported value for 'group_factory': %s" + % group_factory) + el_vertices = [] if dim == 1: @@ -462,8 +535,11 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64): c = vertex_indices[i, j+1] d = vertex_indices[i+1, j+1] - el_vertices.append((a, b, c)) - el_vertices.append((d, c, b)) + if is_tp: + el_vertices.append((a, b, c, d)) + else: + el_vertices.append((a, b, c)) + el_vertices.append((d, c, b)) elif dim == 3: for i in range(shape[0]-1): @@ -480,13 +556,19 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64): a110 = vertex_indices[i+1, j+1, k] a111 = vertex_indices[i+1, j+1, k+1] - el_vertices.append((a000, a100, a010, a001)) - el_vertices.append((a101, a100, a001, a010)) - el_vertices.append((a101, a011, a010, a001)) + if is_tp: + el_vertices.append( + (a000, a001, a010, a011, + a100, a101, a110, a111)) + + else: + el_vertices.append((a000, a100, a010, a001)) + el_vertices.append((a101, a100, a001, a010)) + el_vertices.append((a101, a011, a010, a001)) - el_vertices.append((a100, a010, a101, a110)) - el_vertices.append((a011, a010, a110, a101)) - el_vertices.append((a011, a111, a101, a110)) + el_vertices.append((a100, a010, a101, a110)) + el_vertices.append((a011, a010, a110, a101)) + el_vertices.append((a011, a111, a101, a110)) else: raise NotImplementedError("box meshes of dimension %d" @@ -495,7 +577,8 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64): el_vertices = np.array(el_vertices, dtype=np.int32) grp = make_group_from_vertices( - vertices.reshape(dim, -1), el_vertices, order) + vertices.reshape(dim, -1), el_vertices, order, + group_factory=group_factory) from meshmode.mesh import Mesh return Mesh(vertices, [grp], diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index f3cdaea853d0f9441d7d463f2155ba8d95004d1e..f4d7acc88bdfd9e526f4f495f0c061f470ff047c 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -177,9 +177,20 @@ class GmshMeshReceiver(GmshMeshReceiverBase): np.ones(ngroup_elements, np.bool)) elif isinstance(group_el_type, GmshTensorProductElementBase): + gmsh_vertex_tuples = type(group_el_type)(order=1).gmsh_node_tuples() + gmsh_vertex_tuples_loc_dict = dict( + (gvt, i) + for i, gvt in enumerate(gmsh_vertex_tuples)) + + from pytools import ( + generate_nonnegative_integer_tuples_below as gnitb) + vertex_shuffle = np.array([ + gmsh_vertex_tuples_loc_dict[vt] + for vt in gnitb(2, group_el_type.dimensions)]) + group = TensorProductElementGroup( group_el_type.order, - vertex_indices, + vertex_indices[:, vertex_shuffle], nodes, unit_nodes=unit_nodes ) diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index ff74ac81abdc6e414bd5296b27e8545a047acd25..31abb9beb27d68090ffad264c901b01c69729b63 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -31,6 +31,7 @@ import modepy as mp __doc__ = """ +.. autofunction:: partition_mesh .. autofunction:: find_volume_mesh_element_orientations .. autofunction:: perform_flips .. autofunction:: find_bounding_box @@ -40,6 +41,113 @@ __doc__ = """ """ +# {{{ partition_mesh + +def partition_mesh(mesh, part_per_element, part_nr): + """ + :arg mesh: A :class:`meshmode.mesh.Mesh` to be partitioned. + :arg part_per_element: A :class:`numpy.ndarray` containing one + integer per element of *mesh* indicating which part of the + partitioned mesh the element is to become a part of. + :arg part_nr: The part number of the mesh to return. + + :returns: A tuple ``(part_mesh, part_to_global)``, where *part_mesh* + is a :class:`meshmode.mesh.Mesh` that is a partition of mesh, and + *part_to_global* is a :class:`numpy.ndarray` mapping element + numbers on *part_mesh* to ones in *mesh*. + + .. versionadded:: 2017.1 + + .. warning:: Interface is not final. Connectivity between elements + across groups needs to be added. + """ + assert len(part_per_element) == mesh.nelements, ( + "part_per_element must have shape (mesh.nelements,)") + + # Contains the indices of the elements requested. + queried_elems = np.where(np.array(part_per_element) == part_nr)[0] + + num_groups = len(mesh.groups) + new_indices = [] + new_nodes = [] + + # The set of vertex indices we need. + # NOTE: There are two methods for producing required_indices. + # Optimizations may come from further exploring these options. + #index_set = np.array([], dtype=int) + index_sets = np.array([], dtype=set) + + skip_groups = [] + num_prev_elems = 0 + start_idx = 0 + for group_nr in range(num_groups): + mesh_group = mesh.groups[group_nr] + + # Find the index of first element in the next group + end_idx = len(queried_elems) + for idx in range(start_idx, len(queried_elems)): + if queried_elems[idx] - num_prev_elems >= mesh_group.nelements: + end_idx = idx + break + + if start_idx == end_idx: + skip_groups.append(group_nr) + new_indices.append(np.array([])) + new_nodes.append(np.array([])) + num_prev_elems += mesh_group.nelements + continue + + elems = queried_elems[start_idx:end_idx] - num_prev_elems + new_indices.append(mesh_group.vertex_indices[elems]) + + new_nodes.append( + np.zeros( + (mesh.ambient_dim, end_idx - start_idx, mesh_group.nunit_nodes))) + for i in range(mesh.ambient_dim): + for j in range(start_idx, end_idx): + elems = queried_elems[j] - num_prev_elems + new_idx = j - start_idx + new_nodes[group_nr][i, new_idx, :] = mesh_group.nodes[i, elems, :] + + #index_set = np.append(index_set, new_indices[group_nr].ravel()) + index_sets = np.append(index_sets, set(new_indices[group_nr].ravel())) + + num_prev_elems += mesh_group.nelements + start_idx = end_idx + + # A sorted np.array of vertex indices we need (without duplicates). + #required_indices = np.unique(np.sort(index_set)) + required_indices = np.array(list(set.union(*index_sets))) + + new_vertices = np.zeros((mesh.ambient_dim, len(required_indices))) + for dim in range(mesh.ambient_dim): + new_vertices[dim] = mesh.vertices[dim][required_indices] + + # Our indices need to be in range [0, len(mesh.nelements)]. + for group_nr in range(num_groups): + if group_nr not in skip_groups: + for i in range(len(new_indices[group_nr])): + for j in range(len(new_indices[group_nr][0])): + original_index = new_indices[group_nr][i, j] + new_indices[group_nr][i, j] = np.where( + required_indices == original_index)[0] + + new_mesh_groups = [] + for group_nr in range(num_groups): + if group_nr not in skip_groups: + mesh_group = mesh.groups[group_nr] + new_mesh_groups.append( + type(mesh_group)(mesh_group.order, new_indices[group_nr], + new_nodes[group_nr], unit_nodes=mesh_group.unit_nodes)) + + from meshmode.mesh import Mesh + part_mesh = Mesh(new_vertices, new_mesh_groups) + + return part_mesh, queried_elems + +# }}} + + # {{{ orientations def find_volume_mesh_element_group_orientation(vertices, grp): @@ -262,15 +370,19 @@ def merge_disjoint_meshes(meshes, skip_tests=False, single_group=False): # {{{ assemble new groups list + nodal_adjacency = None + facial_adjacency_groups = None + if single_group: grp_cls = None order = None unit_nodes = None - nodal_adjacency = None for mesh in meshes: if mesh._nodal_adjacency is not None: nodal_adjacency = False + if mesh._facial_adjacency_groups is not None: + facial_adjacency_groups = False for group in mesh.groups: if grp_cls is None: @@ -296,11 +408,12 @@ def merge_disjoint_meshes(meshes, skip_tests=False, single_group=False): else: new_groups = [] - nodal_adjacency = None for mesh, vert_base in zip(meshes, vert_bases): if mesh._nodal_adjacency is not None: nodal_adjacency = False + if mesh._facial_adjacency_groups is not None: + facial_adjacency_groups = False for group in mesh.groups: new_vertex_indices = group.vertex_indices + vert_base @@ -311,7 +424,8 @@ def merge_disjoint_meshes(meshes, skip_tests=False, single_group=False): from meshmode.mesh import Mesh return Mesh(vertices, new_groups, skip_tests=skip_tests, - nodal_adjacency=nodal_adjacency) + nodal_adjacency=nodal_adjacency, + facial_adjacency_groups=facial_adjacency_groups) # }}} diff --git a/meshmode/mesh/visualization.py b/meshmode/mesh/visualization.py index 601086c1570e1b2a87f57b68361010774a78aeb7..52711a0a8b3654a0ea89291c47a7c14ab19d8ed0 100644 --- a/meshmode/mesh/visualization.py +++ b/meshmode/mesh/visualization.py @@ -41,6 +41,11 @@ def draw_2d_mesh(mesh, draw_vertex_numbers=True, draw_element_numbers=True, for iel, el in enumerate(grp.vertex_indices): elverts = mesh.vertices[:, el] + from meshmode.mesh import TensorProductElementGroup + if isinstance(grp, TensorProductElementGroup) and grp.dim == 2: + elverts = elverts[:, + np.array([0, 1, 3, 2])] + pathdata = [ (Path.MOVETO, (elverts[0, 0], elverts[1, 0])), ] diff --git a/requirements.txt b/requirements.txt index 6ed102e6d1309b0e65dcfec56e15aac08448a7d4..f5bb34647a3c1fe4e8989462447a78416879ba4d 100644 --- a/requirements.txt +++ b/requirements.txt @@ -9,3 +9,6 @@ git+https://github.com/inducer/loopy.git git+https://github.com/inducer/boxtree.git git+https://github.com/inducer/sumpy.git git+https://github.com/inducer/pytential.git + +# requires pymetis for tests for partition_mesh +git+https://github.com/inducer/pymetis.git diff --git a/setup.py b/setup.py index 06d1d63ced0aad542ba53440b716073828ab8eb5..0b15c94320e2073c270f68b247422c0cbbb1b844 100644 --- a/setup.py +++ b/setup.py @@ -42,7 +42,6 @@ def main(): packages=find_packages(), install_requires=[ "numpy", - "llist", "modepy", "meshpy>=2014.1", "six", diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 6195e78706e124c550877ebff5ffb8fc67f81786..413a107d1a3d33b52b1a0fd2eea2047822db54ed 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -36,7 +36,7 @@ from pyopencl.tools import ( # noqa from meshmode.discretization.poly_element import ( InterpolatoryQuadratureSimplexGroupFactory, PolynomialWarpAndBlendGroupFactory, - PolynomialEquidistantGroupFactory, + PolynomialEquidistantSimplexGroupFactory, ) from meshmode.mesh import BTAG_ALL from meshmode.discretization.connection import \ @@ -48,6 +48,55 @@ import logging logger = logging.getLogger(__name__) +# {{{ partition_mesh + +def test_partition_torus_mesh(): + from meshmode.mesh.generation import generate_torus + my_mesh = generate_torus(2, 1, n_outer=2, n_inner=2) + + part_per_element = np.array([0, 1, 2, 1, 1, 2, 1, 0]) + + from meshmode.mesh.processing import partition_mesh + (part_mesh0, _) = partition_mesh(my_mesh, part_per_element, 0) + (part_mesh1, _) = partition_mesh(my_mesh, part_per_element, 1) + (part_mesh2, _) = partition_mesh(my_mesh, part_per_element, 2) + + assert part_mesh0.nelements == 2 + assert part_mesh1.nelements == 4 + assert part_mesh2.nelements == 2 + + +def test_partition_boxes_mesh(): + n = 5 + num_parts = 7 + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh1 = generate_regular_rect_mesh(a=(0, 0, 0), b=(1, 1, 1), n=(n, n, n)) + mesh2 = generate_regular_rect_mesh(a=(2, 2, 2), b=(3, 3, 3), n=(n, n, n)) + + from meshmode.mesh.processing import merge_disjoint_meshes + mesh = merge_disjoint_meshes([mesh1, mesh2]) + + adjacency_list = np.zeros((mesh.nelements,), dtype=set) + for elem in range(mesh.nelements): + adjacency_list[elem] = set() + starts = mesh.nodal_adjacency.neighbors_starts + for n in range(starts[elem], starts[elem + 1]): + adjacency_list[elem].add(mesh.nodal_adjacency.neighbors[n]) + + from pymetis import part_graph + (_, p) = part_graph(num_parts, adjacency=adjacency_list) + part_per_element = np.array(p) + + from meshmode.mesh.processing import partition_mesh + new_meshes = [ + partition_mesh(mesh, part_per_element, i)[0] for i in range(num_parts)] + + assert mesh.nelements == np.sum( + [new_meshes[i].nelements for i in range(num_parts)]) + +# }}} + + # {{{ circle mesh def test_circle_mesh(do_plot=False): @@ -405,32 +454,52 @@ def test_element_orientation(): def test_merge_and_map(ctx_getter, visualize=False): from meshmode.mesh.io import generate_gmsh, FileSource + from meshmode.mesh.generation import generate_box_mesh + from meshmode.mesh import TensorProductElementGroup + from meshmode.discretization.poly_element import ( + PolynomialWarpAndBlendGroupFactory, + LegendreGaussLobattoTensorProductGroupFactory) mesh_order = 3 - mesh = generate_gmsh( - FileSource("blob-2d.step"), 2, order=mesh_order, - force_ambient_dim=2, - other_options=["-string", "Mesh.CharacteristicLengthMax = 0.02;"] - ) + if 1: + mesh = generate_gmsh( + FileSource("blob-2d.step"), 2, order=mesh_order, + force_ambient_dim=2, + other_options=["-string", "Mesh.CharacteristicLengthMax = 0.02;"] + ) + + discr_grp_factory = PolynomialWarpAndBlendGroupFactory(3) + else: + mesh = generate_box_mesh( + ( + np.linspace(0, 1, 4), + np.linspace(0, 1, 4), + np.linspace(0, 1, 4), + ), + 10, group_factory=TensorProductElementGroup) + + discr_grp_factory = LegendreGaussLobattoTensorProductGroupFactory(3) from meshmode.mesh.processing import merge_disjoint_meshes, affine_map - mesh2 = affine_map(mesh, A=np.eye(2), b=np.array([5, 0])) + mesh2 = affine_map(mesh, + A=np.eye(mesh.ambient_dim), + b=np.array([5, 0, 0])[:mesh.ambient_dim]) mesh3 = merge_disjoint_meshes((mesh2, mesh)) + mesh3.facial_adjacency_groups + + mesh3.copy() if visualize: from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import \ - PolynomialWarpAndBlendGroupFactory cl_ctx = ctx_getter() queue = cl.CommandQueue(cl_ctx) - discr = Discretization(cl_ctx, mesh3, - PolynomialWarpAndBlendGroupFactory(3)) + discr = Discretization(cl_ctx, mesh3, discr_grp_factory) from meshmode.discretization.visualization import make_visualizer - vis = make_visualizer(queue, discr, 1) + vis = make_visualizer(queue, discr, 3, element_shrink_factor=0.8) vis.write_vtk_file("merged.vtu", []) # }}} @@ -505,7 +574,7 @@ def test_sanity_single_element(ctx_getter, dim, order, visualize=False): # }}} from pytential import bind, sym - bdry_normals = bind(bdry_discr, sym.normal())(queue).as_vector(dtype=object) + bdry_normals = bind(bdry_discr, sym.normal(dim))(queue).as_vector(dtype=object) if visualize: bdry_vis.write_vtk_file("boundary.vtu", [ @@ -514,9 +583,9 @@ def test_sanity_single_element(ctx_getter, dim, order, visualize=False): from pytential import bind, sym normal_outward_check = bind(bdry_discr, - sym.normal() + sym.normal(dim) | - (sym.Nodes() + 0.5*sym.ones_vec(dim)), + (sym.nodes(dim) + 0.5*sym.ones_vec(dim)), )(queue).as_scalar() > 0 assert normal_outward_check.get().all(), normal_outward_check.get() @@ -546,9 +615,9 @@ def test_sanity_qhull_nd(ctx_getter, dim, order): from meshmode.discretization import Discretization low_discr = Discretization(ctx, mesh, - PolynomialEquidistantGroupFactory(order)) + PolynomialEquidistantSimplexGroupFactory(order)) high_discr = Discretization(ctx, mesh, - PolynomialEquidistantGroupFactory(order+1)) + PolynomialEquidistantSimplexGroupFactory(order+1)) from meshmode.discretization.connection import make_same_mesh_connection cnx = make_same_mesh_connection(high_discr, low_discr) @@ -671,7 +740,7 @@ def test_sanity_balls(ctx_getter, src_file, dim, mesh_order, # {{{ check normals point outward normal_outward_check = bind(bdry_discr, - sym.normal() | sym.Nodes(), + sym.normal(mesh.ambient_dim) | sym.nodes(mesh.ambient_dim), )(queue).as_scalar() > 0 assert normal_outward_check.get().all(), normal_outward_check.get() @@ -727,6 +796,12 @@ def test_box_mesh(ctx_getter, visualize=False): # }}} +def test_mesh_copy(): + from meshmode.mesh.generation import generate_box_mesh + mesh = generate_box_mesh(3*(np.linspace(0, 1, 5),)) + mesh.copy() + + # {{{ as_python stringification def test_as_python(): @@ -857,6 +932,50 @@ def no_test_quad_mesh_3d(): # }}} +def test_quad_single_element(): + from meshmode.mesh.generation import make_group_from_vertices + from meshmode.mesh import Mesh, TensorProductElementGroup + + vertices = np.array([ + [0.91, 1.10], + [2.64, 1.27], + [0.97, 2.56], + [3.00, 3.41], + ]).T + mg = make_group_from_vertices( + vertices, + np.array([[0, 1, 2, 3]], dtype=np.int32), + 30, group_factory=TensorProductElementGroup) + + Mesh(vertices, [mg], nodal_adjacency=None, facial_adjacency_groups=None) + if 0: + import matplotlib.pyplot as plt + plt.plot( + mg.nodes[0].reshape(-1), + mg.nodes[1].reshape(-1), "o") + plt.show() + + +def test_quad_multi_element(): + from meshmode.mesh.generation import generate_box_mesh + from meshmode.mesh import TensorProductElementGroup + mesh = generate_box_mesh( + ( + np.linspace(3, 8, 4), + np.linspace(3, 8, 4), + np.linspace(3, 8, 4), + ), + 10, group_factory=TensorProductElementGroup) + + if 0: + import matplotlib.pyplot as plt + mg = mesh.groups[0] + plt.plot( + mg.nodes[0].reshape(-1), + mg.nodes[1].reshape(-1), "o") + plt.show() + + if __name__ == "__main__": import sys if len(sys.argv) > 1: diff --git a/test/test_refinement.py b/test/test_refinement.py index 8098766d446874b08e676d304d5ac73b88e1163a..558b14088dc012a001cd736d5bc2a0fd854862ed 100644 --- a/test/test_refinement.py +++ b/test/test_refinement.py @@ -38,7 +38,7 @@ from meshmode.mesh.refinement import Refiner from meshmode.discretization.poly_element import ( InterpolatoryQuadratureSimplexGroupFactory, PolynomialWarpAndBlendGroupFactory, - PolynomialEquidistantGroupFactory, + PolynomialEquidistantSimplexGroupFactory, ) import logging @@ -144,7 +144,7 @@ def test_refinement(case_name, mesh_gen, flag_gen, num_generations): @pytest.mark.parametrize("group_factory", [ InterpolatoryQuadratureSimplexGroupFactory, PolynomialWarpAndBlendGroupFactory, - PolynomialEquidistantGroupFactory + PolynomialEquidistantSimplexGroupFactory ]) @pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [ ("circle", 1, [20, 30, 40]),