diff --git a/doc/misc.rst b/doc/misc.rst
index 488262acbc29882a93c25d4ba137d77ac37e8660..b3fa9ffae0c88eef24b42f1109e463b0ab2021ea 100644
--- a/doc/misc.rst
+++ b/doc/misc.rst
@@ -3,28 +3,44 @@
Installation
============
-This command should install :mod:`meshmode`::
+This set of instructions is intended for 64-bit Linux and MacOS computers.
- pip install meshmode
+#. Make sure your system has the basics to build software.
-(Note the extra "."!)
+ On Debian derivatives (Ubuntu and many more),
+ installing ``build-essential`` should do the trick.
-You may need to run this with :command:`sudo`.
-If you don't already have `pip <https://pypi.python.org/pypi/pip>`_,
-run this beforehand::
+ Everywhere else, just making sure you have the ``g++`` package should be
+ enough.
- curl -O https://raw.github.com/pypa/pip/master/contrib/get-pip.py
- python get-pip.py
+#. Installing `miniforge for Python 3 on your respective system <https://github.com/conda-forge/miniforge>`_.
-For a more manual installation, `download the source
-<http://pypi.python.org/pypi/meshmode>`_, unpack it, and say::
+#. ``export CONDA=/WHERE/YOU/INSTALLED/miniforge3``
- python setup.py install
+ If you accepted the default location, this should work:
-You may also clone its git repository::
+ ``export CONDA=$HOME/miniforge3``
- git clone --recursive git://github.com/inducer/meshmode
- git clone --recursive http://git.tiker.net/trees/meshmode.git
+#. ``$CONDA/bin/conda create -n dgfem``
+
+#. ``source $CONDA/bin/activate dgfem``
+
+#. ``conda config --add channels conda-forge``
+
+#. ``conda install git pip pocl islpy pyopencl``
+
+#. Type the following command::
+
+ hash -r; for i in pymbolic cgen genpy modepy pyvisfile loopy meshmode; do python -m pip install git+https://github.com/inducer/$i.git; done
+
+Next time you want to use :mod:`meshmode`, just run the following command::
+
+ source /WHERE/YOU/INSTALLED/miniforge3/bin/activate dgfem
+
+You may also like to add this to a startup file (like :file:`$HOME/.bashrc`) or create an alias for it.
+
+After this, you should be able to run the `tests <https://github.com/inducer/meshmode/tree/master/test>`_
+or `examples <https://github.com/inducer/meshmode/tree/master/examples>`_.
User-visible Changes
====================
diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py
index 97570b338f9938165004a46c59fb3ec1570e0cbe..e084525ba9a6eb5001548f432c25a0a8f861173e 100644
--- a/meshmode/discretization/connection/opposite_face.py
+++ b/meshmode/discretization/connection/opposite_face.py
@@ -342,22 +342,42 @@ def make_opposite_face_connection(volume_to_bdry_conn):
# {{{ index wrangling
- # Assert that the adjacency group and the restriction
- # interpolation batch and the adjacency group have the same
- # element ordering.
+ # The elements in the adjacency group will be a subset of
+ # the elements in the restriction interpolation batch:
+ # Imagine an inter-group boundary. The volume-to-boundary
+ # connection will include all faces as targets, whereas
+ # there will be separate adjacency groups for intra- and
+ # inter-group connections.
adj_tgt_flags = adj.element_faces == i_face_tgt
+ adj_els = adj.elements[adj_tgt_flags]
+ if adj_els.size == 0:
+ # NOTE: this case can happen for inter-group boundaries
+ # when all elements are adjacent on the same face
+ # index, so all other ones will be empty
+ continue
+
+ vbc_els = vbc_tgt_grp_face_batch.from_element_indices.get(queue)
+
+ if len(adj_els) == len(vbc_els):
+ # Same length: assert (below) that the two use the same
+ # ordering.
+ vbc_used_els = slice(None)
+
+ else:
+ # Genuine subset: figure out an index mapping.
+ vbc_els_sort_idx = np.argsort(vbc_els)
+ vbc_used_els = vbc_els_sort_idx[np.searchsorted(
+ vbc_els, adj_els, sorter=vbc_els_sort_idx
+ )]
- assert (np.array_equal(
- adj.elements[adj_tgt_flags],
- vbc_tgt_grp_face_batch.from_element_indices
- .get(queue=queue)))
+ assert np.array_equal(vbc_els[vbc_used_els], adj_els)
# find to_element_indices
tgt_bdry_element_indices = (
vbc_tgt_grp_face_batch.to_element_indices
- .get(queue=queue))
+ .get(queue=queue)[vbc_used_els])
# find from_element_indices
diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 97b9bbda9c99fd0dc1f537e6cbbbaf043bc6140c..8cb8c13f3b16d4c0d3fbc9862b45ac944e226082 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -37,6 +37,7 @@ __doc__ = """
.. autofunction:: perform_flips
.. autofunction:: find_bounding_box
.. autofunction:: merge_disjoint_meshes
+.. autofunction:: split_mesh_groups
.. autofunction:: map_mesh
.. autofunction:: affine_map
"""
@@ -569,6 +570,64 @@ def merge_disjoint_meshes(meshes, skip_tests=False, single_group=False):
# }}}
+# {{{ split meshes
+
+def split_mesh_groups(mesh, element_flags, return_subgroup_mapping=False):
+ """Split all the groups in *mesh* according to the values of
+ *element_flags*. The element flags are expected to be integers
+ defining, for each group, how the elements are to be split into
+ subgroups. For example, a single-group mesh with flags::
+
+ element_flags = [0, 0, 0, 42, 42, 42, 0, 0, 0, 41, 41, 41]
+
+ will create three subgroups. The integer flags need not be increasing
+ or contiguous and can repeat across different groups (i.e. they are
+ group-local).
+
+ :arg element_flags: a :class:`numpy.ndarray` with
+ :attr:`~meshmode.mesh.Mesh.nelements` entries
+ indicating how the elements in a group are to be split.
+
+ :returns: a :class:`~meshmode.mesh.Mesh` where each group has been split
+ according to flags in *element_flags*. If *return_subgroup_mapping*
+ is *True*, it also returns a mapping of
+ ``(group_index, subgroup) -> new_group_index``.
+
+ """
+ assert element_flags.shape == (mesh.nelements,)
+
+ new_groups = []
+ subgroup_to_group_map = {}
+
+ for igrp, grp in enumerate(mesh.groups):
+ grp_flags = element_flags[
+ grp.element_nr_base:grp.element_nr_base + grp.nelements]
+ unique_grp_flags = np.unique(grp_flags)
+
+ for flag in unique_grp_flags:
+ subgroup_to_group_map[igrp, flag] = len(new_groups)
+
+ # NOTE: making copies to maintain contiguity of the arrays
+ mask = grp_flags == flag
+ new_groups.append(grp.copy(
+ vertex_indices=grp.vertex_indices[mask, :].copy(),
+ nodes=grp.nodes[:, mask, :].copy()
+ ))
+
+ from meshmode.mesh import Mesh
+ mesh = Mesh(
+ vertices=mesh.vertices,
+ groups=new_groups,
+ is_conforming=mesh.is_conforming)
+
+ if return_subgroup_mapping:
+ return mesh, subgroup_to_group_map
+ else:
+ return mesh
+
+# }}}
+
+
# {{{ map
def map_mesh(mesh, f): # noqa
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index ded91fe297a14c1d2241fb67486fbc0843fa237a..047362d20c1b297e22eed758fbe633bee01a82cc 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -1348,6 +1348,89 @@ def test_is_affine_group_check(mesh_name):
assert all(grp.is_affine for grp in mesh.groups) == is_affine
+@pytest.mark.parametrize("ambient_dim", [1, 2, 3])
+def test_mesh_multiple_groups(ctx_factory, ambient_dim, visualize=False):
+ ctx = ctx_factory()
+ queue = cl.CommandQueue(ctx)
+
+ order = 4
+
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+ mesh = generate_regular_rect_mesh(
+ a=(-0.5,)*ambient_dim, b=(0.5,)*ambient_dim,
+ n=(8,)*ambient_dim, order=order)
+ assert len(mesh.groups) == 1
+
+ from meshmode.mesh.processing import split_mesh_groups
+ element_flags = np.any(
+ mesh.vertices[0, mesh.groups[0].vertex_indices] < 0.0,
+ axis=1).astype(np.int)
+ mesh = split_mesh_groups(mesh, element_flags)
+
+ assert len(mesh.groups) == 2
+ assert mesh.facial_adjacency_groups
+ assert mesh.nodal_adjacency
+
+ if visualize and ambient_dim == 2:
+ from meshmode.mesh.visualization import draw_2d_mesh
+ draw_2d_mesh(mesh,
+ draw_vertex_numbers=False,
+ draw_element_numbers=True,
+ draw_face_numbers=False,
+ set_bounding_box=True)
+
+ import matplotlib.pyplot as plt
+ plt.savefig("test_mesh_multiple_groups_2d_elements.png", dpi=300)
+
+ from meshmode.discretization import Discretization
+ from meshmode.discretization.poly_element import \
+ PolynomialWarpAndBlendGroupFactory as GroupFactory
+ discr = Discretization(ctx, mesh, GroupFactory(order))
+
+ if visualize:
+ group_id = discr.empty(queue, dtype=np.int)
+ for igrp, grp in enumerate(discr.groups):
+ group_id_view = grp.view(group_id)
+ group_id_view.fill(igrp)
+
+ from meshmode.discretization.visualization import make_visualizer
+ vis = make_visualizer(queue, discr, vis_order=order)
+ vis.write_vtk_file("test_mesh_multiple_groups.vtu", [
+ ("group_id", group_id)
+ ], overwrite=True)
+
+ # check face restrictions
+ from meshmode.discretization.connection import (
+ make_face_restriction,
+ make_face_to_all_faces_embedding,
+ make_opposite_face_connection,
+ check_connection)
+ for boundary_tag in [BTAG_ALL, FACE_RESTR_INTERIOR, FACE_RESTR_ALL]:
+ conn = make_face_restriction(discr, GroupFactory(order),
+ boundary_tag=boundary_tag,
+ per_face_groups=False)
+ check_connection(conn)
+
+ bdry_f = conn.to_discr.empty(queue)
+ bdry_f.fill(1.0)
+
+ if boundary_tag == FACE_RESTR_INTERIOR:
+ opposite = make_opposite_face_connection(conn)
+ check_connection(opposite)
+
+ op_bdry_f = opposite(queue, bdry_f)
+ error = abs(cl.array.sum(bdry_f - op_bdry_f).get(queue))
+ assert error < 1.0e-11, error
+
+ if boundary_tag == FACE_RESTR_ALL:
+ embedding = make_face_to_all_faces_embedding(conn, conn.to_discr)
+ check_connection(embedding)
+
+ em_bdry_f = embedding(queue, bdry_f)
+ error = abs(cl.array.sum(bdry_f - em_bdry_f).get(queue))
+ assert error < 1.0e-11, error
+
+
if __name__ == "__main__":
import sys
if len(sys.argv) > 1: