doc/misc.rst

@@ -3,28 +3,53 @@
 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
+#.  Install `miniforge <https://github.com/conda-forge/miniforge>`_::
 
-For a more manual installation, `download the source
-<http://pypi.python.org/pypi/meshmode>`_, unpack it, and say::
+        curl -L -O "https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-$(uname)-$(uname -m).sh"
+        # then run
+        bash ./Miniforge3-*.sh
 
-    python setup.py install
+#.  ``export CONDA=/WHERE/YOU/INSTALLED/miniforge3``
 
-You may also clone its git repository::
+    If you accepted the default location, this should work:
 
-    git clone --recursive git://github.com/inducer/meshmode
-    git clone --recursive http://git.tiker.net/trees/meshmode.git
+    ``export CONDA=$HOME/miniforge3``
+
+#.  ``$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 arraycontext meshmode; do python -m pip install --editable "git+https://github.com/inducer/$i.git#egg=$i"; done
+
+.. note::
+
+    In each case, you may leave out the ``--editable`` flag if you would not like
+    a checkout of the source code.
+
+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
 ====================
@@ -69,11 +94,15 @@ OTHER DEALINGS IN THE SOFTWARE.
 Acknowledgments
 ===============
 
-Andreas Klöckner's work on :mod:`meshmode` was supported in part by
+Work on meshmode was supported in part by
+
+* the Department of Energy, National Nuclear Security Administration,
+  under Award Number DE-NA0003963,
+* the US Navy ONR, under grant number N00014-14-1-0117, and
+* the US National Science Foundation under grant numbers DMS-1418961, CCF-1524433,
+  DMS-1654756, SHF-1911019, and OAC-1931577.
 
-* US Navy ONR grant number N00014-14-1-0117
-* the US National Science Foundation under grant numbers DMS-1418961 and CCF-1524433.
+AK also gratefully acknowledges a hardware gift from Nvidia Corporation.
 
-AK also gratefully acknowledges a hardware gift from Nvidia Corporation.  The
-views and opinions expressed herein do not necessarily reflect those of the
+The views and opinions expressed herein do not necessarily reflect those of the
 funding agencies.

doc/upload-docs.sh

 #! /bin/sh
 
-rsync --verbose --archive --delete _build/html/* doc-upload:doc/meshmode
+rsync --verbose --archive --delete _build/html/ doc-upload:doc/meshmode

examples/from_firedrake.py

+__copyright__ = "Copyright (C) 2020 Benjamin Sepanski"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import pyopencl as cl
+
+from meshmode.array_context import PyOpenCLArrayContext
+
+
+# This example provides a brief template for bringing information in
+# from firedrake and makes some plots to help you better understand
+# what a FromBoundaryFiredrakeConnection does
+def main(visualize=True):
+    # If can't import firedrake, do nothing
+    #
+    # filename MUST include "firedrake" (i.e. match *firedrake*.py) in order
+    # to be run during CI
+    try:
+        import firedrake  # noqa : F401
+    except ImportError:
+        return 0
+
+    from firedrake import (
+        Function,
+        FunctionSpace,
+        SpatialCoordinate,
+        UnitSquareMesh,
+        cos,
+    )
+
+    from meshmode.interop.firedrake import build_connection_from_firedrake
+
+    # Create a firedrake mesh and interpolate cos(x+y) onto it
+    fd_mesh = UnitSquareMesh(10, 10)
+    fd_fspace = FunctionSpace(fd_mesh, "DG", 2)
+    spatial_coord = SpatialCoordinate(fd_mesh)
+    fd_fntn = Function(fd_fspace).interpolate(cos(sum(spatial_coord)))
+
+    # Make connections
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+    actx = PyOpenCLArrayContext(queue)
+
+    fd_connection = build_connection_from_firedrake(actx, fd_fspace)
+    fd_bdy_connection = \
+        build_connection_from_firedrake(actx,
+                                        fd_fspace,
+                                        restrict_to_boundary="on_boundary")
+
+    # Plot the meshmode meshes that the connections connect to
+    import matplotlib.pyplot as plt
+
+    from meshmode.mesh.visualization import draw_2d_mesh
+    fig, (ax1, ax2) = plt.subplots(1, 2)
+    ax1.set_title("FiredrakeConnection")
+    plt.sca(ax1)
+    draw_2d_mesh(fd_connection.discr.mesh,
+                 draw_vertex_numbers=False,
+                 draw_element_numbers=False,
+                 set_bounding_box=True)
+    ax2.set_title("FiredrakeConnection 'on_boundary'")
+    plt.sca(ax2)
+    draw_2d_mesh(fd_bdy_connection.discr.mesh,
+                 draw_vertex_numbers=False,
+                 draw_element_numbers=False,
+                 set_bounding_box=True)
+    plt.show()
+
+    # Plot fd_fntn using unrestricted FiredrakeConnection
+    from meshmode.discretization.visualization import make_visualizer
+    discr = fd_connection.discr
+    vis = make_visualizer(actx, discr, discr.groups[0].order+3)
+    field = fd_connection.from_firedrake(fd_fntn, actx=actx)
+
+    if visualize:
+        fig = plt.figure()
+        ax1 = fig.add_subplot(1, 2, 1, projection="3d")
+        ax1.set_title("cos(x+y) in\nFiredrakeConnection")
+        vis.show_scalar_in_matplotlib_3d(field, do_show=False)
+
+    # Now repeat using FiredrakeConnection restricted to "on_boundary"
+    bdy_discr = fd_bdy_connection.discr
+    bdy_vis = make_visualizer(actx, bdy_discr, bdy_discr.groups[0].order+3)
+    bdy_field = fd_bdy_connection.from_firedrake(fd_fntn, actx=actx)
+
+    if visualize:
+        ax2 = fig.add_subplot(1, 2, 2, projection="3d")
+        plt.sca(ax2)
+        ax2.set_title("cos(x+y) in\nFiredrakeConnection 'on_boundary'")
+        bdy_vis.show_scalar_in_matplotlib_3d(bdy_field, do_show=False)
+
+        import matplotlib.cm as cm
+        fig.colorbar(cm.ScalarMappable(), ax=ax2)
+        plt.show()
+
+
+if __name__ == "__main__":
+    main()
+
+# vim: foldmethod=marker

examples/mesh-to-tikz.py

+from meshmode.mesh.io import FileSource, generate_gmsh
+from meshmode.mesh.visualization import mesh_to_tikz
+
+
+h = 0.3
+order = 1
+
+mesh = generate_gmsh(
+        FileSource("../test/blob-2d.step"), 2, order=order,
+        force_ambient_dim=2,
+        other_options=[
+            "-string", "Mesh.CharacteristicLengthMax = %s;" % h]
+        )
+
+print(mesh_to_tikz(mesh))

examples/moving-geometry.py

+__copyright__ = "Copyright (C) 2021 Alexandru Fikl"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import logging
+
+import numpy as np
+
+import pyopencl as cl
+from pytools import keyed_memoize_in
+from pytools.obj_array import make_obj_array
+
+from meshmode.array_context import PyOpenCLArrayContext
+from meshmode.transform_metadata import FirstAxisIsElementsTag
+
+
+logger = logging.getLogger(__name__)
+logging.basicConfig(level=logging.INFO)
+
+
+def plot_solution(actx, vis, filename, discr, t, x):
+    names_and_fields = []
+
+    try:
+        from pytential import bind, sym
+        kappa = bind(discr, sym.mean_curvature(discr.ambient_dim))(actx)
+        names_and_fields.append(("kappa", kappa))
+    except ImportError:
+        pass
+
+    vis.write_vtk_file(filename, names_and_fields, overwrite=True)
+
+
+def reconstruct_discr_from_nodes(actx, discr, x):
+    @keyed_memoize_in(actx,
+            (reconstruct_discr_from_nodes, "to_mesh_interp_matrix"),
+            lambda grp: grp.discretization_key())
+    def to_mesh_interp_matrix(grp) -> np.ndarray:
+        import modepy as mp
+        mat = mp.resampling_matrix(
+                grp.basis_obj().functions,
+                grp.mesh_el_group.unit_nodes,
+                grp.unit_nodes)
+
+        return actx.freeze(actx.from_numpy(mat))
+
+    def resample_nodes_to_mesh(grp, igrp, iaxis):
+        discr_nodes = x[iaxis][igrp]
+
+        grp_unit_nodes = grp.unit_nodes.reshape(-1)
+        meg_unit_nodes = grp.mesh_el_group.unit_nodes.reshape(-1)
+
+        tol = 10 * np.finfo(grp_unit_nodes.dtype).eps
+        if (grp_unit_nodes.shape == meg_unit_nodes.shape
+                and np.linalg.norm(grp_unit_nodes - meg_unit_nodes) < tol):
+            return discr_nodes
+
+        return actx.einsum("ij,ej->ei",
+                           to_mesh_interp_matrix(grp),
+                           discr_nodes,
+                           tagged=(FirstAxisIsElementsTag(),))
+
+    from dataclasses import replace
+
+    megs = []
+    for igrp, grp in enumerate(discr.groups):
+        nodes = np.stack([
+            actx.to_numpy(resample_nodes_to_mesh(grp, igrp, iaxis))
+            for iaxis in range(discr.ambient_dim)
+            ])
+
+        meg = replace(grp.mesh_el_group, vertex_indices=None, nodes=nodes)
+        megs.append(meg)
+
+    mesh = discr.mesh.copy(groups=megs, vertices=None)
+    return discr.copy(actx, mesh=mesh)
+
+
+def advance(actx, dt, t, x, fn):
+    # NOTE: everybody's favorite three stage SSP RK3 method
+    k1 = x + dt * fn(t, x)
+    k2 = 3.0 / 4.0 * x + 1.0 / 4.0 * (k1 + dt * fn(t + dt, k1))
+    return 1.0 / 3.0 * x + 2.0 / 3.0 * (k2 + dt * fn(t + 0.5 * dt, k2))
+
+
+def run(actx, *,
+        ambient_dim: int = 3,
+        resolution: int | None = None,
+        target_order: int = 4,
+        tmax: float = 1.0,
+        timestep: float = 1.0e-2,
+        group_factory_name: str = "warp_and_blend",
+        visualize: bool = True):
+    if ambient_dim not in (2, 3):
+        raise ValueError(f"unsupported dimension: {ambient_dim}")
+
+    mesh_order = target_order
+    radius = 1.0
+
+    # {{{ geometry
+
+    # {{{ element groups
+
+    import modepy as mp
+
+    import meshmode.discretization.poly_element as poly
+
+    # NOTE: picking the same unit nodes for the mesh and the discr saves
+    # a bit of work when reconstructing after a time step
+
+    if group_factory_name == "warp_and_blend":
+        group_factory_cls: type[poly.HomogeneousOrderBasedGroupFactory] = (
+            poly.PolynomialWarpAndBlend2DRestrictingGroupFactory)
+
+        unit_nodes = mp.warp_and_blend_nodes(ambient_dim - 1, mesh_order)
+    elif group_factory_name == "quadrature":
+        group_factory_cls = poly.InterpolatoryQuadratureSimplexGroupFactory
+
+        if ambient_dim == 2:
+            unit_nodes = mp.LegendreGaussQuadrature(
+                    mesh_order, force_dim_axis=True).nodes
+        else:
+            unit_nodes = mp.VioreanuRokhlinSimplexQuadrature(mesh_order, 2).nodes
+    else:
+        raise ValueError(f"unknown group factory: '{group_factory_name}'")
+
+    # }}}
+
+    # {{{ discretization
+
+    import meshmode.mesh.generation as gen
+    if ambient_dim == 2:
+        nelements = 8192 if resolution is None else resolution
+        mesh = gen.make_curve_mesh(
+                lambda t: radius * gen.ellipse(1.0, t),
+                np.linspace(0.0, 1.0, nelements + 1),
+                order=mesh_order,
+                unit_nodes=unit_nodes)
+    else:
+        nrounds = 4 if resolution is None else resolution
+        mesh = gen.generate_sphere(radius,
+                uniform_refinement_rounds=nrounds,
+                order=mesh_order,
+                unit_nodes=unit_nodes)
+
+    from meshmode.discretization import Discretization
+    discr0 = Discretization(actx, mesh, group_factory_cls(target_order))
+
+    logger.info("ndofs:     %d", discr0.ndofs)
+    logger.info("nelements: %d", discr0.mesh.nelements)
+
+    # }}}
+
+    if visualize:
+        from meshmode.discretization.visualization import make_visualizer
+        vis = make_visualizer(actx, discr0,
+                vis_order=target_order,
+                # NOTE: setting this to True will add some unnecessary
+                # resampling in Discretization.nodes for the vis_discr underneath
+                force_equidistant=False)
+
+    # }}}
+
+    # {{{ ode
+
+    def velocity_field(nodes, alpha=1.0):
+        return make_obj_array([
+            alpha * nodes[0], -alpha * nodes[1], 0.0 * nodes[0]
+            ][:ambient_dim])
+
+    def source(t, x):
+        discr = reconstruct_discr_from_nodes(actx, discr0, x)
+        u = velocity_field(actx.thaw(discr.nodes()))
+
+        # {{{
+
+        # NOTE: these are just here because this was at some point used to
+        # profile some more operators (turned out well!)
+
+        from meshmode.discretization import num_reference_derivative
+        x = actx.thaw(discr.nodes()[0])
+        gradx = sum(
+                num_reference_derivative(discr, (i,), x)
+                for i in range(discr.dim))
+        intx = sum(actx.np.sum(xi * wi)
+                   for xi, wi in zip(x, discr.quad_weights(), strict=True))
+
+        assert gradx is not None
+        assert intx is not None
+
+        # }}}
+
+        return u
+
+    # }}}
+
+    # {{{ evolve
+
+    maxiter = int(tmax // timestep) + 1
+    dt = tmax / maxiter + 1.0e-15
+
+    x = actx.thaw(discr0.nodes())
+    t = 0.0
+
+    if visualize:
+        filename = f"moving-geometry-{0:09d}.vtu"
+        plot_solution(actx, vis, filename, discr0, t, x)
+
+    for n in range(1, maxiter + 1):
+        x = advance(actx, dt, t, x, source)
+        t += dt
+
+        if visualize:
+            discr = reconstruct_discr_from_nodes(actx, discr0, x)
+            vis = make_visualizer(actx, discr, vis_order=target_order)
+            # vis = vis.copy_with_same_connectivity(actx, discr)
+
+            filename = f"moving-geometry-{n:09d}.vtu"
+            plot_solution(actx, vis, filename, discr, t, x)
+
+        logger.info("[%05d/%05d] t = %.5e/%.5e dt = %.5e",
+                n, maxiter, t, tmax, dt)
+
+    # }}}
+
+
+if __name__ == "__main__":
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+    actx = PyOpenCLArrayContext(queue)
+
+    from pytools import ProcessTimer
+    for _ in range(1):
+        with ProcessTimer() as p:
+            run(actx,
+                    ambient_dim=3,
+                    group_factory_name="warp_and_blend",
+                    tmax=1.0,
+                    timestep=1.0e-2,
+                    visualize=False)
+
+        logger.info("elapsed: %.3fs wall %.2fx cpu",
+                p.wall_elapsed, p.process_elapsed / p.wall_elapsed)

examples/multiple-meshes.py

-from __future__ import division
+import sys
 
-import numpy as np  # noqa
+import numpy as np
 
 
 order = 4
 
 
 def main():
-    from meshmode.mesh.generation import (  # noqa
-            make_curve_mesh, starfish)
+    from meshmode.mesh.generation import make_curve_mesh, starfish
     mesh1 = make_curve_mesh(starfish, np.linspace(0, 1, 20), 4)
 
     from meshmode.mesh.processing import affine_map, merge_disjoint_meshes
@@ -20,7 +19,11 @@ def main():
     draw_2d_mesh(mesh, set_bounding_box=True)
 
     import matplotlib.pyplot as pt
-    pt.show()
+    if sys.stdin.isatty():
+        pt.show()
+    else:
+        pt.savefig("plot.pdf")
+
 
 if __name__ == "__main__":
     main()

examples/parallel-vtkhdf.py

+import logging
+
+import numpy as np
+
+from meshmode.mesh import Mesh
+
+
+logger = logging.getLogger(__file__)
+
+
+def make_example_mesh(ambient_dim: int, nelements: int, order: int) -> Mesh:
+    import meshmode.mesh.generation as mgen
+    if ambient_dim == 2:
+        return mgen.make_curve_mesh(
+            mgen.starfish,
+            np.linspace(0.0, 1.0, nelements + 1),
+            order=order)
+    else:
+        return mgen.generate_torus(4.0, 2.0,
+            n_major=nelements, n_minor=nelements // 2,
+            order=order)
+
+
+def main(*, ambient_dim: int) -> None:
+    logging.basicConfig(level=logging.INFO)
+
+    import h5py
+    if "mpio" not in h5py.registered_drivers():
+        logger.info("h5py does not have the 'mpio' driver")
+        return
+
+    from meshmode import _acf
+    actx = _acf()
+
+    from mpi4py import MPI
+    comm = MPI.COMM_WORLD
+
+    from meshmode.distributed import membership_list_to_map
+    from meshmode.mesh.processing import partition_mesh
+
+    order = 5
+    nelements = 64 if ambient_dim == 3 else 256
+
+    logger.info("[%4d] distributing mesh: started", comm.rank)
+
+    if comm.rank == 0:
+        mesh = make_example_mesh(ambient_dim, nelements, order=order)
+        logger.info("[%4d] mesh: nelements %d nvertices %d",
+                    comm.rank, mesh.nelements, mesh.nvertices)
+
+        rng = np.random.default_rng()
+
+        part_id_to_part = partition_mesh(mesh,
+                       membership_list_to_map(
+                           rng.integers(comm.size, size=mesh.nelements)))
+        parts = [part_id_to_part[i] for i in range(comm.size)]
+        local_mesh = comm.scatter(parts)
+    else:
+        local_mesh = comm.scatter(None)
+
+    logger.info("[%4d] distributing mesh: finished", comm.rank)
+
+    from meshmode.discretization import Discretization
+    from meshmode.discretization.poly_element import default_simplex_group_factory
+    discr = Discretization(actx, local_mesh,
+        default_simplex_group_factory(local_mesh.dim, order=order))
+
+    logger.info("[%4d] discretization: finished", comm.rank)
+
+    vector_field = actx.thaw(discr.nodes())
+    scalar_field = actx.np.sin(vector_field[0])
+    part_id = 1.0 + comm.rank + discr.zeros(actx)
+    logger.info("[%4d] fields: finished", comm.rank)
+
+    from meshmode.discretization.visualization import make_visualizer
+    vis = make_visualizer(actx, discr, force_equidistant=False)
+    logger.info("[%4d] make_visualizer: finished", comm.rank)
+
+    filename = f"parallel-vtkhdf-example-{ambient_dim}d.hdf"
+    vis.write_vtkhdf_file(filename, [
+        ("scalar", scalar_field),
+        ("vector", vector_field),
+        ("part_id", part_id)
+        ], comm=comm, overwrite=True, use_high_order=False)
+
+    logger.info("[%4d] write: finished: %s", comm.rank, filename)
+
+
+if __name__ == "__main__":
+    try:
+        main(ambient_dim=2)
+        main(ambient_dim=3)
+    except ImportError as exc:
+        logger.info(exc)

examples/plot-connectivity.py

-from __future__ import division
-
 import numpy as np  # noqa
+
 import pyopencl as cl
 
+from meshmode.array_context import PyOpenCLArrayContext
+
 
 order = 4
 
@@ -10,30 +11,34 @@ order = 4
 def main():
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
+    actx = PyOpenCLArrayContext(queue)
+
+    from meshmode.mesh.generation import (  # noqa: F401
+        generate_icosahedron,
+        generate_sphere,
+        generate_torus,
+    )
 
-    from meshmode.mesh.generation import (  # noqa
-            generate_icosphere, generate_icosahedron,
-            generate_torus)
-    #mesh = generate_icosphere(1, order=order)
+    #mesh = generate_sphere(1, order=order)
     mesh = generate_icosahedron(1, order=order)
     #mesh = generate_torus(3, 1, order=order)
 
     from meshmode.discretization import Discretization
-    from meshmode.discretization.poly_element import \
-            PolynomialWarpAndBlendGroupFactory
+    from meshmode.discretization.poly_element import (
+        PolynomialWarpAndBlend3DRestrictingGroupFactory,
+    )
 
     discr = Discretization(
-            cl_ctx, mesh, PolynomialWarpAndBlendGroupFactory(order))
+            actx, mesh, PolynomialWarpAndBlend3DRestrictingGroupFactory(order))
 
     from meshmode.discretization.visualization import make_visualizer
-    vis = make_visualizer(queue, discr, order)
+    vis = make_visualizer(actx, discr, order)
 
     vis.write_vtk_file("geometry.vtu", [
-        ("f", discr.nodes()[0]),
+        ("f", actx.thaw(discr.nodes()[0])),
         ])
 
-    from meshmode.discretization.visualization import \
-            write_nodal_adjacency_vtk_file
+    from meshmode.discretization.visualization import write_nodal_adjacency_vtk_file
 
     write_nodal_adjacency_vtk_file("adjacency.vtu",
             mesh)

examples/refinement-playground.py

-from __future__ import division, print_function
-
-from six.moves import range
-import numpy as np  # noqa
-import pyopencl as cl
-import random
-import os
-import logging
-order = 1
-
-from meshmode.mesh.refinement.utils import check_nodal_adj_against_geometry
-from meshmode.mesh.refinement import Refiner
-
-
-#construct vertex vertex_index
-def remove_if_exists(name):
-    from errno import ENOENT
-    try:
-        os.remove(name)
-    except OSError as e:
-        if e.errno == ENOENT:
-            pass
-        else:
-            raise
-
-
-def linear_func(vert):
-    csum = 0
-    for i in vert:
-        csum += i
-        #print csum
-    return csum
-
-
-def sine_func(vert):
-    #print vert
-    import math
-    res = 1
-    for i in vert:
-        res *= math.sin(i*2.0*math.pi)
-    #print res
-    return abs(res) * 0.2 + 0.2
-
-
-def quadratic_func(vert):
-    csum = 0
-    for i in vert:
-        csum += i * i
-    return csum * 1.5
-
-
-def get_function_flags(mesh, function):
-    from math import sqrt
-    flags = np.zeros(len(mesh.groups[0].vertex_indices))
-    for grp in mesh.groups:
-        for iel_grp in range(grp.nelements):
-            vertex_indices = grp.vertex_indices[iel_grp]
-            max_edge_len = 0
-            for i in range(len(vertex_indices)):
-                for j in range(i+1, len(vertex_indices)):
-                    edge_len = 0
-                    for k in range(len(mesh.vertices)):
-                        edge_len += (
-                                (mesh.vertices[k, vertex_indices[i]]
-                                    - mesh.vertices[k, vertex_indices[j]])
-                                * (mesh.vertices[k, vertex_indices[i]]
-                                    - mesh.vertices[k, vertex_indices[j]]))
-                    edge_len = sqrt(edge_len)
-                    max_edge_len = max(max_edge_len, edge_len)
-                #print(edge_lens[0], mesh.vertices[0, vertex_indices[i]], mesh.vertices[1, vertex_indices[i]], mesh.vertices[2, vertex_indices[i]])  # noqa
-                centroid = [0] * len(mesh.vertices)
-                for j in range(len(mesh.vertices)):
-                    centroid[j] += mesh.vertices[j, vertex_indices[i]]
-            for i in range(len(mesh.vertices)):
-                centroid[i] /= len(vertex_indices)
-            val = function(centroid)
-            if max_edge_len > val:
-                flags[iel_grp] = True
-    return flags
-
-
-def get_corner_flags(mesh):
-    flags = np.zeros(len(mesh.groups[0].vertex_indices))
-    for grp in mesh.groups:
-        for iel_grp in range(grp.nelements):
-            is_corner_el = False
-            vertex_indices = grp.vertex_indices[iel_grp]
-            for i in range(len(vertex_indices)):
-                cur_vertex_corner = True
-                for j in range(len(mesh.vertices)):
-                    print(iel_grp, i, mesh.vertices[j, vertex_indices[i]])
-                    if mesh.vertices[j, vertex_indices[i]] != 0.0:
-                        cur_vertex_corner = False
-                if cur_vertex_corner:
-                    is_corner_el = True
-                    break
-            if is_corner_el:
-                print(iel_grp)
-                flags[iel_grp] = True
-    return flags
-
-
-def get_random_flags(mesh):
-    flags = np.zeros(len(mesh.groups[0].vertex_indices))
-    for i in range(0, len(flags)):
-            flags[i] = random.randint(0, 1)
-    return flags
-
-
-def refine_and_generate_chart_function(mesh, filename, function):
-    from time import clock
-    cl_ctx = cl.create_some_context()
-    queue = cl.CommandQueue(cl_ctx)
-    print("NELEMENTS: ", mesh.nelements)
-    #print mesh
-    for i in range(len(mesh.groups[0].vertex_indices[0])):
-        for k in range(len(mesh.vertices)):
-            print(mesh.vertices[k, i])
-
-    #check_nodal_adj_against_geometry(mesh);
-    r = Refiner(mesh)
-    #random.seed(0)
-    #times = 3
-    num_elements = []
-    time_t = []
-    #nelements = mesh.nelements
-    while True:
-        print("NELS:", mesh.nelements)
-        #flags = get_corner_flags(mesh)
-        flags = get_function_flags(mesh, function)
-        nels = 0
-        for i in flags:
-            if i:
-                nels += 1
-        if nels == 0:
-            break
-        print("LKJASLFKJALKASF:", nels)
-        num_elements.append(nels)
-        #flags = get_corner_flags(mesh)
-        beg = clock()
-        mesh = r.refine(flags)
-        end = clock()
-        time_taken = end - beg
-        time_t.append(time_taken)
-        #if nelements == mesh.nelements:
-            #break
-        #nelements = mesh.nelements
-        #from meshmode.mesh.visualization import draw_2d_mesh
-        #draw_2d_mesh(mesh, True, True, True, fill=None)
-        #import matplotlib.pyplot as pt
-        #pt.show()
-
-        #poss_flags = np.zeros(len(mesh.groups[0].vertex_indices))
-        #for i in range(0, len(flags)):
-        #    poss_flags[i] = flags[i]
-        #for i in range(len(flags), len(poss_flags)):
-        #    poss_flags[i] = 1
-
-    import matplotlib.pyplot as pt
-    pt.xlabel('Number of elements being refined')
-    pt.ylabel('Time taken')
-    pt.plot(num_elements, time_t, "o")
-    pt.savefig(filename, format='pdf')
-    pt.clf()
-    print('DONE REFINING')
-    '''
-    flags = np.zeros(len(mesh.groups[0].vertex_indices))
-    flags[0] = 1
-    flags[1] = 1
-    mesh = r.refine(flags)
-    flags = np.zeros(len(mesh.groups[0].vertex_indices))
-    flags[0] = 1
-    flags[1] = 1
-    flags[2] = 1
-    mesh = r.refine(flags)
-    '''
-    #check_nodal_adj_against_geometry(mesh)
-    #r.print_rays(70)
-    #r.print_rays(117)
-    #r.print_hanging_elements(10)
-    #r.print_hanging_elements(117)
-    #r.print_hanging_elements(757)
-    #from meshmode.mesh.visualization import draw_2d_mesh
-    #draw_2d_mesh(mesh, False, False, False, fill=None)
-    #import matplotlib.pyplot as pt
-    #pt.show()
-
-    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)
-    remove_if_exists("connectivity2.vtu")
-    remove_if_exists("geometry2.vtu")
-    vis.write_vtk_file("geometry2.vtu", [
-        ("f", discr.nodes()[0]),
-        ])
-
-    from meshmode.discretization.visualization import \
-            write_nodal_adjacency_vtk_file
-
-    write_nodal_adjacency_vtk_file("connectivity2.vtu",
-            mesh)
-
-
-def main2():
-    from meshmode.mesh.generation import (  # noqa
-            generate_icosphere, generate_icosahedron,
-            generate_torus, generate_regular_rect_mesh,
-            generate_box_mesh)
-#    mesh = generate_icosphere(1, order=order)
-    #mesh = generate_icosahedron(1, order=order)
-    #mesh = generate_torus(3, 1, order=order)
-    #mesh = generate_regular_rect_mesh()
-    #mesh =  generate_box_mesh(3*(np.linspace(0, 3, 5),))
-    #mesh =  generate_box_mesh(3*(np.linspace(0, 1, 3),))
-    mesh = generate_box_mesh(3*(np.linspace(0, 1, 5),))
-    refine_and_generate_chart_function(mesh, "plot.pdf", sine_func)
-
-
-def all_refine(num_mesh, depth, fname):
-    from meshmode.mesh.generation import (  # noqa
-            generate_icosphere, generate_icosahedron,
-            generate_torus, generate_regular_rect_mesh,
-            generate_box_mesh)
-    import timeit
-    nelements = []
-    runtimes = []
-    for el_fact in range(2, num_mesh+2):
-        mesh = generate_box_mesh(3*(np.linspace(0, 1, el_fact),))
-        r = Refiner(mesh)
-        for time in range(depth):
-            flags = np.ones(len(mesh.groups[0].vertex_indices))
-            if time < depth-1:
-                mesh = r.refine(flags)
-            else:
-                start = timeit.default_timer()
-                mesh = r.refine(flags)
-                stop = timeit.default_timer()
-                nelements.append(mesh.nelements)
-                runtimes.append(stop-start)
-        check_nodal_adj_against_geometry(mesh)
-    import matplotlib.pyplot as pt
-    pt.plot(nelements, runtimes, "o")
-    pt.savefig(fname)
-    pt.clf()
-    #pt.show()
-
-
-def uniform_refine(num_mesh, fract, depth, fname):
-    from meshmode.mesh.generation import (  # noqa
-            generate_icosphere, generate_icosahedron,
-            generate_torus, generate_regular_rect_mesh,
-            generate_box_mesh)
-    import timeit
-    nelements = []
-    runtimes = []
-    for el_fact in range(2, num_mesh+2):
-        mesh = generate_box_mesh(3*(np.linspace(0, 1, el_fact),))
-        r = Refiner(mesh)
-        all_els = list(range(mesh.nelements))
-        for time in range(depth):
-            print("EL_FACT", el_fact, "TIME", time)
-            flags = np.zeros(mesh.nelements)
-            from random import shuffle, seed
-            seed(1)
-            shuffle(all_els)
-            nels_this_round = 0
-            for i in range(len(all_els)):
-                if i / len(flags) > fract:
-                    break
-                flags[all_els[i]] = 1
-                nels_this_round += 1
-
-            if time < depth-1:
-                mesh = r.refine(flags)
-            else:
-                start = timeit.default_timer()
-                mesh = r.refine(flags)
-                stop = timeit.default_timer()
-                nelements.append(mesh.nelements)
-                runtimes.append(stop-start)
-            all_els = []
-            for i in range(len(flags)):
-                if flags[i]:
-                    all_els.append(i)
-            for i in range(len(flags), mesh.nelements):
-                all_els.append(i)
-            check_nodal_adj_against_geometry(mesh)
-
-    import matplotlib.pyplot as pt
-    pt.plot(nelements, runtimes, "o")
-    pt.savefig(fname)
-    pt.clf()
-
-
-if __name__ == "__main__":
-    logging.basicConfig(level="DEBUG")
-    print("HEREERERE")
-    #all_refine(3, 2, 'all_a.pdf')
-    all_refine(3, 3, 'all_b.pdf')
-    #uniform_refine(3, 0.2, 3, 'uniform_a.pdf')
-    #main2()

examples/simple-dg.py

+__copyright__ = "Copyright (C) 2020 Andreas Kloeckner"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import logging
+from dataclasses import dataclass
+
+import numpy as np
+import numpy.linalg as la  # noqa
+
+import pyopencl as cl
+import pyopencl.array as cla  # noqa
+from arraycontext import (
+    ArrayContainer,
+    dataclass_array_container,
+    map_array_container,
+    with_container_arithmetic,
+)
+from pytools import log_process, memoize_method
+from pytools.obj_array import flat_obj_array, make_obj_array
+
+from meshmode.array_context import PyOpenCLArrayContext, PytatoPyOpenCLArrayContext
+from meshmode.dof_array import DOFArray, flat_norm
+from meshmode.mesh import BTAG_ALL, BTAG_NONE  # noqa
+from meshmode.transform_metadata import FirstAxisIsElementsTag
+
+
+logger = logging.getLogger(__name__)
+
+
+# Features lost vs. https://github.com/inducer/grudge:
+# - dimension independence / differential geometry
+# - overintegration
+# - operator fusion
+# - distributed-memory
+
+
+# {{{ discretization
+
+def parametrization_derivative(actx, discr):
+    thawed_nodes = actx.thaw(discr.nodes())
+
+    from meshmode.discretization import num_reference_derivative
+    result = np.zeros((discr.ambient_dim, discr.dim), dtype=object)
+    for iambient in range(discr.ambient_dim):
+        for idim in range(discr.dim):
+            result[iambient, idim] = num_reference_derivative(discr,
+                    (idim,), thawed_nodes[iambient])
+
+    return result
+
+
+class DGDiscretization:
+    def __init__(self, actx, mesh, order):
+        self.order = order
+
+        from meshmode.discretization import Discretization
+        from meshmode.discretization.poly_element import (
+            PolynomialWarpAndBlend2DRestrictingGroupFactory,
+        )
+        self.group_factory = PolynomialWarpAndBlend2DRestrictingGroupFactory(
+                order=order)
+        self.volume_discr = Discretization(actx, mesh, self.group_factory)
+
+        assert self.volume_discr.dim == 2
+
+    @property
+    def _setup_actx(self):
+        return self.volume_discr._setup_actx
+
+    @property
+    def dim(self):
+        return self.volume_discr.dim
+
+    # {{{ discretizations/connections
+
+    @memoize_method
+    def boundary_connection(self, boundary_tag):
+        from meshmode.discretization.connection import make_face_restriction
+        return make_face_restriction(
+                        self.volume_discr._setup_actx,
+                        self.volume_discr,
+                        self.group_factory,
+                        boundary_tag=boundary_tag)
+
+    @memoize_method
+    def interior_faces_connection(self):
+        from meshmode.discretization.connection import (
+            FACE_RESTR_INTERIOR,
+            make_face_restriction,
+        )
+        return make_face_restriction(
+                        self.volume_discr._setup_actx,
+                        self.volume_discr,
+                        self.group_factory,
+                        FACE_RESTR_INTERIOR,
+                        per_face_groups=False)
+
+    @memoize_method
+    def opposite_face_connection(self):
+        from meshmode.discretization.connection import make_opposite_face_connection
+
+        return make_opposite_face_connection(
+                self._setup_actx, self.interior_faces_connection())
+
+    @memoize_method
+    def all_faces_connection(self):
+        from meshmode.discretization.connection import (
+            FACE_RESTR_ALL,
+            make_face_restriction,
+        )
+        return make_face_restriction(
+                        self.volume_discr._setup_actx,
+                        self.volume_discr,
+                        self.group_factory,
+                        FACE_RESTR_ALL,
+                        per_face_groups=False)
+
+    @memoize_method
+    def get_to_all_face_embedding(self, where):
+        from meshmode.discretization.connection import make_face_to_all_faces_embedding
+
+        faces_conn = self.get_connection("vol", where)
+        return make_face_to_all_faces_embedding(
+                self._setup_actx, faces_conn, self.get_discr("all_faces"))
+
+    def get_connection(self, src, tgt):
+        src_tgt = (src, tgt)
+
+        if src_tgt == ("vol", "int_faces"):
+            return self.interior_faces_connection()
+        elif src_tgt == ("vol", "all_faces"):
+            return self.all_faces_connection()
+        elif src_tgt == ("vol", BTAG_ALL):
+            return self.boundary_connection(tgt)
+        elif src_tgt == ("int_faces", "all_faces"):
+            return self.get_to_all_face_embedding(src)
+        elif src_tgt == (BTAG_ALL, "all_faces"):
+            return self.get_to_all_face_embedding(src)
+        else:
+            raise ValueError(f"locations '{src}'->'{tgt}' not understood")
+
+    def interp(self, src, tgt, vec):
+        return self.get_connection(src, tgt)(vec)
+
+    def get_discr(self, where):
+        if where == "vol":
+            return self.volume_discr
+        elif where == "all_faces":
+            return self.all_faces_connection().to_discr
+        elif where == "int_faces":
+            return self.interior_faces_connection().to_discr
+        elif where == BTAG_ALL:
+            return self.boundary_connection(where).to_discr
+        else:
+            raise ValueError(f"location '{where}' not understood")
+
+    # }}}
+
+    @memoize_method
+    def parametrization_derivative(self):
+        return self._setup_actx.freeze(
+                parametrization_derivative(self._setup_actx, self.volume_discr))
+
+    @memoize_method
+    def vol_jacobian(self):
+        [a, b], [c, d] = self._setup_actx.thaw(self.parametrization_derivative())
+        return self._setup_actx.freeze(a*d - b*c)
+
+    @memoize_method
+    def inverse_parametrization_derivative(self):
+        [a, b], [c, d] = self._setup_actx.thaw(self.parametrization_derivative())
+
+        result = np.zeros((2, 2), dtype=object)
+        det = a*d-b*c
+        result[0, 0] = d/det
+        result[0, 1] = -b/det
+        result[1, 0] = -c/det
+        result[1, 1] = a/det
+
+        return self._setup_actx.freeze(result)
+
+    def zeros(self, actx):
+        return self.volume_discr.zeros(actx)
+
+    def grad(self, vec):
+        ipder = vec.array_context.thaw(self.inverse_parametrization_derivative())
+
+        from meshmode.discretization import num_reference_derivative
+        dref = [
+                num_reference_derivative(self.volume_discr, (idim,), vec)
+                for idim in range(self.volume_discr.dim)]
+
+        return make_obj_array([
+            sum(dref_i*ipder_i
+                for dref_i, ipder_i in zip(dref, ipder[iambient], strict=True))
+            for iambient in range(self.volume_discr.ambient_dim)])
+
+    def div(self, vecs):
+        return sum(
+                self.grad(vec_i)[i] for i, vec_i in enumerate(vecs))
+
+    @memoize_method
+    def normal(self, where):
+        bdry_discr = self.get_discr(where)
+
+        ((a,), (b,)) = parametrization_derivative(self._setup_actx, bdry_discr)
+
+        nrm = 1/(a**2+b**2)**0.5
+        return self._setup_actx.freeze(flat_obj_array(b*nrm, -a*nrm))
+
+    @memoize_method
+    def face_jacobian(self, where):
+        bdry_discr = self.get_discr(where)
+
+        ((a,), (b,)) = parametrization_derivative(self._setup_actx, bdry_discr)
+
+        return self._setup_actx.freeze((a**2 + b**2)**0.5)
+
+    @memoize_method
+    def get_inverse_mass_matrix(self, grp, dtype):
+        import modepy as mp
+        matrix = mp.inverse_mass_matrix(grp.basis_obj(), grp.unit_nodes)
+
+        actx = self._setup_actx
+        return actx.freeze(actx.from_numpy(matrix))
+
+    def inverse_mass(self, vec):
+        if not isinstance(vec, DOFArray):
+            return map_array_container(self.inverse_mass, vec)
+
+        actx = vec.array_context
+        dtype = vec.entry_dtype
+        discr = self.volume_discr
+
+        return DOFArray(
+            actx,
+            data=tuple(
+                actx.einsum(
+                    "ij,ej->ei",
+                    self.get_inverse_mass_matrix(grp, dtype),
+                    vec_i,
+                    arg_names=("mass_inv_mat", "vec"),
+                    tagged=(FirstAxisIsElementsTag(),)
+                ) for grp, vec_i in zip(discr.groups, vec, strict=True)
+            )
+        ) / actx.thaw(self.vol_jacobian())
+
+    @memoize_method
+    def get_local_face_mass_matrix(self, afgrp, volgrp, dtype):
+        nfaces = volgrp.mesh_el_group.nfaces
+        assert afgrp.nelements == nfaces * volgrp.nelements
+
+        matrix = np.empty(
+                (volgrp.nunit_dofs,
+                    nfaces,
+                    afgrp.nunit_dofs),
+                dtype=dtype)
+
+        import modepy as mp
+        shape = mp.Simplex(volgrp.dim)
+
+        for face in mp.faces_for_shape(shape):
+            from modepy import PN, Simplex, basis_for_space, quadrature_for_space
+            face_space = PN(volgrp.dim - 1, volgrp.order)
+            face_shape = Simplex(volgrp.dim-1)
+            face_quad = quadrature_for_space(face_space, face_shape)
+            face_basis = basis_for_space(face_space, face_shape)
+
+            matrix[:, face.face_index, :] = mp.nodal_mass_matrix_for_face(
+                    face, face_quad,
+                    face_basis.functions,
+                    volgrp.basis_obj().functions,
+                    volgrp.unit_nodes, afgrp.unit_nodes)
+
+        actx = self._setup_actx
+        return actx.freeze(actx.from_numpy(matrix))
+
+    def face_mass(self, vec):
+        if not isinstance(vec, DOFArray):
+            return map_array_container(self.face_mass, vec)
+
+        actx = vec.array_context
+        dtype = vec.entry_dtype
+
+        all_faces_conn = self.get_connection("vol", "all_faces")
+        all_faces_discr = all_faces_conn.to_discr
+        vol_discr = all_faces_conn.from_discr
+
+        fj = vec.array_context.thaw(self.face_jacobian("all_faces"))
+        vec = vec*fj
+
+        assert len(all_faces_discr.groups) == len(vol_discr.groups)
+
+        return DOFArray(
+            actx,
+            data=tuple(
+                actx.einsum("ifj,fej->ei",
+                            self.get_local_face_mass_matrix(afgrp, volgrp, dtype),
+                            vec_i.reshape(
+                                volgrp.mesh_el_group.nfaces,
+                                volgrp.nelements,
+                                afgrp.nunit_dofs
+                            ),
+                            tagged=(FirstAxisIsElementsTag(),))
+                for afgrp, volgrp, vec_i in zip(all_faces_discr.groups,
+                                                vol_discr.groups,
+                                                vec, strict=True)
+            )
+        )
+
+# }}}
+
+
+# {{{ trace pair
+
+@with_container_arithmetic(
+        bcast_obj_array=False, eq_comparison=False, rel_comparison=False)
+@dataclass_array_container
+@dataclass(frozen=True)
+class TracePair:
+    where: str
+    interior: ArrayContainer
+    exterior: ArrayContainer
+
+    # NOTE: let the container do the broadcasting + arithmetic
+    __array_ufunc__ = None
+
+    def __getattr__(self, name):
+        return map_array_container(
+                lambda ary: getattr(ary, name),
+                self)
+
+    @property
+    def int(self):
+        return self.interior
+
+    @property
+    def ext(self):
+        return self.exterior
+
+    @property
+    def avg(self):
+        return 0.5*(self.int + self.ext)
+
+
+def interior_trace_pair(discr, vec):
+    i = discr.interp("vol", "int_faces", vec)
+    e = discr.opposite_face_connection()(i)
+    return TracePair("int_faces", interior=i, exterior=e)
+
+# }}}
+
+
+# {{{ wave equation bits
+
+def wave_flux(actx, discr, c, q_tpair):
+    u = q_tpair.u
+    v = q_tpair.v
+
+    normal = actx.thaw(discr.normal(q_tpair.where))
+
+    flux_weak = WaveState(
+            u=np.dot(v.avg, normal),
+            v=normal * u.avg)
+
+    # upwind
+    v_jump = np.dot(normal, v.ext-v.int)
+    flux_weak += WaveState(
+            u=0.5*(u.ext-u.int),
+            v=0.5*normal*v_jump)
+
+    flux_strong = WaveState(
+            u=np.dot(v.int, normal),
+            v=u.int * normal,
+            ) - flux_weak
+
+    return discr.interp(q_tpair.where, "all_faces", c*flux_strong)
+
+
+def wave_operator(actx, discr, c, q):
+    dir_q = discr.interp("vol", BTAG_ALL, q)
+    dir_bc = WaveState(u=-dir_q.u, v=dir_q.v)
+
+    return (
+            WaveState(
+                u=c*discr.div(q.v),
+                v=c*discr.grad(q.u)
+                )
+            -
+            discr.inverse_mass(
+                discr.face_mass(
+                    wave_flux(actx, discr, c=c,
+                        q_tpair=interior_trace_pair(discr, q))
+                    + wave_flux(actx, discr, c=c,
+                        q_tpair=TracePair(BTAG_ALL, dir_q, dir_bc))
+                    ))
+                )
+
+# }}}
+
+
+def rk4_step(y, t, h, f):
+    k1 = f(t, y)
+    k2 = f(t+h/2, y + h/2*k1)
+    k3 = f(t+h/2, y + h/2*k2)
+    k4 = f(t+h, y + h*k3)
+    return y + h/6*(k1 + 2*k2 + 2*k3 + k4)
+
+
+def bump(actx, discr, t=0):
+    source_center = np.array([0.0, 0.05])
+    source_width = 0.05
+    source_omega = 3
+
+    nodes = actx.thaw(discr.volume_discr.nodes())
+    center_dist = flat_obj_array([
+        nodes[0] - source_center[0],
+        nodes[1] - source_center[1],
+        ])
+
+    return (
+        np.cos(source_omega*t)
+        * actx.np.exp(
+            -np.dot(center_dist, center_dist)
+            / source_width**2))
+
+
+@with_container_arithmetic(bcast_obj_array=True,
+                           rel_comparison=True,
+                           _cls_has_array_context_attr=True)
+@dataclass_array_container
+@dataclass(frozen=True)
+class WaveState:
+    u: DOFArray
+    v: np.ndarray  # [object]
+
+    # NOTE: let the container do the broadcasting + arithmetic
+    __array_ufunc__ = None
+
+    def __post_init__(self):
+        assert isinstance(self.v, np.ndarray) and self.v.dtype.char == "O"
+
+    @property
+    def array_context(self):
+        return self.u.array_context
+
+
+@log_process(logger)
+def main(lazy=False):
+    logging.basicConfig(level=logging.INFO)
+
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+    actx_outer = PyOpenCLArrayContext(queue)
+    if lazy:
+        actx_rhs = PytatoPyOpenCLArrayContext(queue)
+    else:
+        actx_rhs = actx_outer
+
+    nel_1d = 16
+    from meshmode.mesh.generation import generate_regular_rect_mesh
+    mesh = generate_regular_rect_mesh(
+            a=(-0.5, -0.5),
+            b=(0.5, 0.5),
+            nelements_per_axis=(nel_1d, nel_1d))
+
+    order = 3
+
+    # no deep meaning here, just a fudge factor
+    dt = 0.7/(nel_1d*order**2)
+
+    logger.info("%d elements", mesh.nelements)
+
+    discr = DGDiscretization(actx_outer, mesh, order=order)
+
+    fields = WaveState(
+            u=bump(actx_outer, discr),
+            v=make_obj_array([discr.zeros(actx_outer) for i in range(discr.dim)]),
+            )
+
+    from meshmode.discretization.visualization import make_visualizer
+    vis = make_visualizer(actx_outer, discr.volume_discr)
+
+    def rhs(t, q):
+        return wave_operator(actx_rhs, discr, c=1, q=q)
+
+    compiled_rhs = actx_rhs.compile(rhs)
+
+    def rhs_wrapper(t, q):
+        r = compiled_rhs(t, actx_rhs.thaw(actx_outer.freeze(q)))
+        return actx_outer.thaw(actx_rhs.freeze(r))
+
+    t = np.float64(0)
+    t_final = 3
+    istep = 0
+    while t < t_final:
+        fields = rk4_step(fields, t, dt, rhs_wrapper)
+
+        if istep % 10 == 0:
+            # FIXME: Maybe an integral function to go with the
+            # DOFArray would be nice?
+            assert len(fields.u) == 1
+            logger.info("[%05d] t %.5e / %.5e norm %.5e",
+                    istep, t, t_final, actx_outer.to_numpy(flat_norm(fields.u, 2)))
+            vis.write_vtk_file("fld-wave-min-%04d.vtu" % istep, [
+                ("q", fields),
+                ])
+
+        t += dt
+        istep += 1
+
+    assert flat_norm(fields.u, 2) < 100
+
+
+if __name__ == "__main__":
+    import argparse
+    parser = argparse.ArgumentParser(description="Wave Equation Solver")
+    parser.add_argument("--lazy", action="store_true",
+                        help="switch to a lazy computation mode")
+    args = parser.parse_args()
+    main(lazy=args.lazy)
+
+# vim: foldmethod=marker

examples/to_firedrake.py

+__copyright__ = "Copyright (C) 2020 Benjamin Sepanski"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+
+import numpy as np
+
+import pyopencl as cl
+
+from meshmode.array_context import PyOpenCLArrayContext
+
+
+# Nb: Some of the initial setup was adapted from meshmode/examplse/simple-dg.py
+#     written by Andreas Klockner:
+#    https://gitlab.tiker.net/inducer/meshmode/-/blob/7826fa5e13854bf1dae425b4226865acc10ee01f/examples/simple-dg.py  # noqa : E501
+def main():
+    # If can't import firedrake, do nothing
+    #
+    # filename MUST include "firedrake" (i.e. match *firedrake*.py) in order
+    # to be run during CI
+    try:
+        import firedrake  # noqa : F401
+    except ImportError:
+        return 0
+
+    # For this example, imagine we wish to solve the Laplace equation
+    # on a meshmode mesh with some given Dirichlet boundary conditions,
+    # and decide to use firedrake.
+    #
+    # To verify this is working, we use a solution to the wave equation
+    # to get our boundary conditions
+
+    # {{{ First we make a discretization in meshmode and get our bcs
+
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+    actx = PyOpenCLArrayContext(queue)
+
+    nel_1d = 16
+    from meshmode.mesh.generation import generate_regular_rect_mesh
+    mesh = generate_regular_rect_mesh(
+            a=(-0.5, -0.5),
+            b=(0.5, 0.5),
+            nelements_per_axis=(nel_1d, nel_1d))
+
+    order = 3
+
+    from meshmode.discretization import Discretization
+    from meshmode.discretization.poly_element import (
+        InterpolatoryQuadratureSimplexGroupFactory,
+    )
+    group_factory = InterpolatoryQuadratureSimplexGroupFactory(order=order)
+    discr = Discretization(actx, mesh, group_factory)
+
+    # Get our solution: we will use
+    # Real(e^z) = Real(e^{x+i y})
+    #           = e^x Real(e^{i y})
+    #           = e^x cos(y)
+    nodes = discr.nodes()
+    for i in range(len(nodes)):
+        nodes[i] = actx.thaw(nodes[i])
+    # First index is dimension
+    candidate_sol = actx.np.exp(nodes[0]) * actx.np.cos(nodes[1])
+
+    # }}}
+
+    # {{{ Now send candidate_sol into firedrake and use it for boundary conditions
+
+    from meshmode.interop.firedrake import build_connection_to_firedrake
+    fd_connection = build_connection_to_firedrake(discr, group_nr=0)
+    # convert candidate_sol to firedrake
+    fd_candidate_sol = fd_connection.from_meshmode(candidate_sol)
+    # get the firedrake function space
+    fd_fspace = fd_connection.firedrake_fspace()
+
+    # set up dirichlet laplace problem in fd and solve
+    from firedrake import (
+        Constant,
+        DirichletBC,
+        Function,
+        FunctionSpace,
+        TestFunction,
+        TrialFunction,
+        dx,
+        grad,
+        inner,
+        project,
+        solve,
+    )
+
+    # because it's easier to write down the variational problem,
+    # we're going to project from our "DG" space
+    # into a continuous one.
+    cfd_fspace = FunctionSpace(fd_fspace.mesh(), "CG", order)
+    u = TrialFunction(cfd_fspace)
+    v = TestFunction(cfd_fspace)
+    sol = Function(cfd_fspace)
+
+    a = inner(grad(u), grad(v)) * dx
+    rhs = inner(Constant(0.0), v) * dx
+    bc_value = project(fd_candidate_sol, cfd_fspace)
+    bc = DirichletBC(cfd_fspace, bc_value, "on_boundary")
+    params = {"ksp_monitor": None}
+    solve(a == rhs, sol, bcs=[bc], solver_parameters=params)
+
+    # project back into our "DG" space
+    sol = project(sol, fd_fspace)
+
+    # }}}
+
+    # {{{ Take the solution from firedrake and compare it to candidate_sol
+
+    true_sol = fd_connection.from_firedrake(sol, actx=actx)
+    # pull back into numpy
+    true_sol = actx.to_numpy(true_sol[0])
+    candidate_sol = actx.to_numpy(candidate_sol[0])
+    print("l^2 difference between candidate solution and firedrake solution=",
+          np.linalg.norm(true_sol - candidate_sol))
+
+    # }}}
+
+
+if __name__ == "__main__":
+    main()
+
+# vim: foldmethod=marker

examples/tp-lagrange-stl.py

+# Requires numpy-stl
+# pip install numpy-stl
+
+import numpy as np
+import numpy.linalg as la
+
+import modepy
+import pyopencl as cl
+
+import meshmode.discretization.connection as conn
+import meshmode.discretization.poly_element as poly
+import meshmode.mesh.generation as mgen
+from meshmode.array_context import PyOpenCLArrayContext
+from meshmode.discretization import Discretization
+from meshmode.mesh import BTAG_ALL, make_mesh
+
+
+def main():
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+    actx = PyOpenCLArrayContext(queue)
+
+    npts1d = 100  # use 300 for actual print
+    rs_coords = np.linspace(-1, 1, npts1d)
+    mesh = mgen.generate_box_mesh((
+                rs_coords,
+                rs_coords,
+                np.linspace(-0.1, 0, 10),
+                ))
+
+    group_factory = poly.PolynomialWarpAndBlend2DRestrictingGroupFactory(1)
+    discr = Discretization(actx, mesh, group_factory)
+
+    frestr = conn.make_face_restriction(actx, discr, group_factory, BTAG_ALL)
+
+    bdry_mesh = frestr.to_discr.mesh
+
+    order = 10
+    hc_shape = modepy.Hypercube(2)
+    space = modepy.QN(2, order)
+    basis = modepy.basis_for_space(space, hc_shape)
+
+    from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes
+    nodes_1d, = legendre_gauss_lobatto_nodes(order, force_dim_axis=True)
+
+    nodes = modepy.tensor_product_nodes(2, nodes_1d)
+    vdm_inv = la.inv(modepy.vandermonde(basis.functions, nodes))
+
+    vertices = bdry_mesh.vertices
+
+    top_surface = np.abs(vertices[2] - 0) < 1e-12
+
+    rs = vertices[:2, top_surface]
+
+    val = 0*rs[0]
+
+    for i, bf in enumerate(basis.functions):
+        val += 0.5*vdm_inv[i, len(nodes_1d)*4 + 3] * bf(rs)
+
+    for i, bf in enumerate(basis.functions):
+        val += 0.5*vdm_inv[i, len(nodes_1d)*6 + 4] * bf(rs)
+
+    for r1d in nodes_1d:
+        r1d_discrete = rs_coords[np.argmin(np.abs(r1d - rs_coords))]
+        assert abs(r1d - r1d_discrete) < 2/npts1d
+
+        dimple = np.zeros_like(val)
+        at_node_r = np.abs(rs[0] - r1d_discrete) < 1e-12
+        dimple[at_node_r] = 1
+
+        at_node_s = np.abs(rs[1] - r1d_discrete) < 1e-12
+        dimple[at_node_s] = 1
+
+        val -= 0.005 * dimple
+
+    vertices[2, top_surface] = 0.1+val
+
+    grp, = bdry_mesh.groups
+    from meshmode.mesh.generation import make_group_from_vertices
+    mod_grp = make_group_from_vertices(vertices, grp.vertex_indices, order=grp.order)
+
+    mod_mesh = make_mesh(
+            vertices=vertices, groups=[mod_grp],
+            is_conforming=bdry_mesh.is_conforming)
+
+    from meshmode.mesh.visualization import write_stl_file
+    write_stl_file(mod_mesh, "tp-lagrange.stl")
+
+
+if __name__ == "__main__":
+    main()
+
+# vim: foldmethod=marker

meshmode/__init__.py

-from __future__ import division
-
 __copyright__ = "Copyright (C) 2014 Andreas Kloeckner"
 
 __license__ = """
@@ -24,14 +22,67 @@ THE SOFTWARE.
 
 
 __doc__ = """
-.. exception:: Error
-.. exception:: DataUnavailable
+.. autoexception:: Error
+.. autoexception:: DataUnavailableError
+
+.. autoexception:: InconsistentMeshError
+.. autoexception:: InconsistentArrayDTypeError
+.. autoexception:: InconsistentVerticesError
+.. autoexception:: InconsistentAdjacencyError
 """
 
+from builtins import FileExistsError  # noqa: F401
+
+from meshmode.mesh.tools import AffineMap  # noqa: F401
+
 
 class Error(RuntimeError):
-    pass
+    """Exception base for :mod:`meshmode` errors."""
+
+
+class DataUnavailableError(Error):
+    """Raised when some data on the mesh or the discretization is not available.
+
+    This error should not be raised when the specific data simply fails to be
+    computed for other reasons.
+    """
+
+
+DataUnavailable = DataUnavailableError
+
+
+class InconsistentMeshError(Error):
+    """Raised when the mesh is inconsistent in some fashion.
+
+    Prefer the more specific exceptions, e.g. :exc:`InconsistentVerticesError`
+    when possible.
+    """
+
+
+class InconsistentArrayDTypeError(InconsistentMeshError):
+    """Raised when a mesh (or group) array does not match the provided
+    :class:`~numpy.dtype`.
+    """
+
+
+class InconsistentVerticesError(InconsistentMeshError):
+    """Raised when an element's local-to-global mapping does not map the unit
+    vertices to the corresponding values in the mesh's *vertices* array.
+    """
+
+
+class InconsistentAdjacencyError(InconsistentMeshError):
+    """Raised when the nodal or the facial adjacency is inconsistent."""
+
+
+def _acf():
+    """A tiny undocumented function to pass to tests that take an ``actx_factory``
+    argument when running them from the command line.
+    """
+    import pyopencl as cl
 
+    from meshmode.array_context import PyOpenCLArrayContext
 
-class DataUnavailable(Error):
-    pass
+    context = cl._csc()
+    queue = cl.CommandQueue(context)
+    return PyOpenCLArrayContext(queue)

meshmode/array_context.py

+"""
+.. autoclass:: PyOpenCLArrayContext
+.. autoclass:: PytatoPyOpenCLArrayContext
+"""
+
+__copyright__ = "Copyright (C) 2020 Andreas Kloeckner"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+from warnings import warn
+
+from arraycontext import (
+    PyOpenCLArrayContext as PyOpenCLArrayContextBase,
+    PytatoPyOpenCLArrayContext as PytatoPyOpenCLArrayContextBase,
+)
+from arraycontext.pytest import (
+    _PytestPyOpenCLArrayContextFactoryWithClass,
+    _PytestPytatoPyOpenCLArrayContextFactory,
+    register_pytest_array_context_factory,
+)
+
+
+def thaw(actx, ary):
+    warn("meshmode.array_context.thaw is deprecated. Use arraycontext.thaw instead. "
+            "WARNING: The argument order is reversed between these two functions. "
+            "meshmode.array_context.thaw will continue to work until 2022.",
+            DeprecationWarning, stacklevel=2)
+
+    return actx.thaw(ary)
+
+
+# {{{ kernel transform function
+
+def _transform_loopy_inner(t_unit):
+    import loopy as lp
+    from arraycontext.transform_metadata import ElementwiseMapKernelTag
+    from pymbolic.primitives import Subscript, Variable
+
+    from meshmode.transform_metadata import FirstAxisIsElementsTag
+
+    default_ep = t_unit.default_entrypoint
+
+    # FIXME: Firedrake branch lacks kernel tags
+    kernel_tags = getattr(default_ep, "tags", ())
+
+    # {{{ FirstAxisIsElementsTag on kernel (compatibility)
+
+    if any(isinstance(tag, FirstAxisIsElementsTag) for tag in kernel_tags):
+        if (len(default_ep.instructions) != 1
+                or not isinstance(
+                    default_ep.instructions[0], lp.Assignment)):
+            raise ValueError("FirstAxisIsElementsTag may only be applied to "
+                    "a kernel if the kernel contains a single assignment.")
+
+        stmt, = default_ep.instructions
+
+        if not isinstance(stmt.assignee, Subscript):
+            raise ValueError("single assignment in FirstAxisIsElementsTag kernel "
+                    "must be a subscript")
+
+        output_name = stmt.assignee.aggregate.name
+        new_args = [
+                arg.tagged(FirstAxisIsElementsTag())
+                if arg.name == output_name else arg
+                for arg in default_ep.args]
+        default_ep = default_ep.copy(args=new_args)
+        t_unit = t_unit.with_kernel(default_ep)
+
+    # }}}
+
+    # {{{ ElementwiseMapKernelTag on kernel
+
+    if any(isinstance(tag, ElementwiseMapKernelTag) for tag in kernel_tags):
+        el_inames = []
+        dof_inames = []
+        for stmt in default_ep.instructions:
+            if isinstance(stmt, lp.MultiAssignmentBase):
+                for assignee in stmt.assignees:
+                    if isinstance(assignee, Variable):
+                        # some scalar assignee kernel => no concurrency in the
+                        # workload => skip
+                        continue
+                    if not isinstance(assignee, Subscript):
+                        raise ValueError("assignees in "
+                                "ElementwiseMapKernelTag-tagged kernels must be "
+                                "subscripts")
+
+                    for i, subscript in enumerate(assignee.index_tuple[:2]):
+                        if (not isinstance(subscript, Variable)
+                                or subscript.name not in default_ep.all_inames()):
+                            raise ValueError("subscripts in "
+                                    "ElementwiseMapKernelTag-tagged kernels must be "
+                                    "inames")
+
+                        if i == 0:
+                            el_inames.append(subscript.name)
+                        elif i == 1:
+                            dof_inames.append(subscript.name)
+
+        return _transform_with_element_and_dof_inames(t_unit, el_inames, dof_inames)
+
+    # }}}
+
+    # {{{ FirstAxisIsElementsTag on output variable
+
+    first_axis_el_args = [arg.name for arg in default_ep.args
+            if any(isinstance(tag, FirstAxisIsElementsTag) for tag in arg.tags)]
+
+    if first_axis_el_args:
+        el_inames = []
+        dof_inames = []
+
+        for stmt in default_ep.instructions:
+            if isinstance(stmt, lp.MultiAssignmentBase):
+                for assignee in stmt.assignees:
+                    if not isinstance(assignee, Subscript):
+                        raise ValueError("assignees in "
+                                "FirstAxisIsElementsTag-tagged kernels must be "
+                                "subscripts")
+
+                    if assignee.aggregate.name not in first_axis_el_args:
+                        continue
+
+                    subscripts = assignee.index_tuple[:2]
+
+                    for i, subscript in enumerate(subscripts):
+                        if (not isinstance(subscript, Variable)
+                                or subscript.name not in default_ep.all_inames()):
+                            raise ValueError("subscripts in "
+                                    "FirstAxisIsElementsTag-tagged kernels must be "
+                                    "inames")
+
+                        if i == 0:
+                            el_inames.append(subscript.name)
+                        elif i == 1:
+                            dof_inames.append(subscript.name)
+        return _transform_with_element_and_dof_inames(t_unit, el_inames, dof_inames)
+
+    # }}}
+
+    # {{{ element/dof iname tag
+
+    from meshmode.transform_metadata import (
+        ConcurrentDOFInameTag,
+        ConcurrentElementInameTag,
+    )
+    el_inames = [iname.name
+            for iname in default_ep.inames.values()
+            if ConcurrentElementInameTag() in iname.tags]
+    dof_inames = [iname.name
+            for iname in default_ep.inames.values()
+            if ConcurrentDOFInameTag() in iname.tags]
+
+    if el_inames:
+        return _transform_with_element_and_dof_inames(t_unit, el_inames, dof_inames)
+
+    # }}}
+
+    # *shrug* no idea how to transform this thing.
+    return None
+
+
+def _transform_with_element_and_dof_inames(t_unit, el_inames, dof_inames):
+    import loopy as lp
+
+    if set(el_inames) & set(dof_inames):
+        raise ValueError("Some inames are marked as both 'element' and 'dof' "
+                "inames. These must be disjoint.")
+
+    # Sorting ensures the same order of transformations is used every
+    # time; avoids accidentally generating cache misses or kernel
+    # hash conflicts.
+
+    for dof_iname in sorted(dof_inames):
+        t_unit = lp.split_iname(t_unit, dof_iname, 32, inner_tag="l.0")
+    for el_iname in sorted(el_inames):
+        t_unit = lp.tag_inames(t_unit, {el_iname: "g.0"})
+    return t_unit
+
+# }}}
+
+
+# {{{ pyopencl array context subclass
+
+class PyOpenCLArrayContext(PyOpenCLArrayContextBase):
+    """Extends :class:`arraycontext.PyOpenCLArrayContext` with knowledge about
+    program transformation for finite element programs.
+
+    See :mod:`meshmode.transform_metadata` for relevant metadata.
+    """
+
+    def transform_loopy_program(self, t_unit):
+        default_ep = t_unit.default_entrypoint
+        options = default_ep.options
+        if not (options.return_dict and options.no_numpy):
+            raise ValueError("Loopy kernel passed to call_loopy must "
+                    "have return_dict and no_numpy options set. "
+                    "Did you use arraycontext.make_loopy_program "
+                    "to create this kernel?")
+
+        transformed_t_unit = _transform_loopy_inner(t_unit)
+
+        if transformed_t_unit is not None:
+            return transformed_t_unit
+
+        warn("meshmode.array_context.PyOpenCLArrayContext."
+                "transform_loopy_program fell back on "
+                "arraycontext.PyOpenCLArrayContext to find a transform for "
+                f"'{default_ep.name}'. "
+                "Please update your program to use metadata from "
+                "meshmode.transform_metadata. "
+                "This code path will stop working in 2022.",
+                DeprecationWarning, stacklevel=3)
+
+        return super().transform_loopy_program(t_unit)
+
+# }}}
+
+
+# {{{ pytato pyopencl array context subclass
+
+class PytatoPyOpenCLArrayContext(PytatoPyOpenCLArrayContextBase):
+    def transform_dag(self, dag):
+        dag = super().transform_dag(dag)
+
+        # {{{ /!\ Remove tags from NamedArrays
+        # See <https://www.github.com/inducer/pytato/issues/195>
+
+        import pytato as pt
+
+        def untag_loopy_call_results(expr):
+            if isinstance(expr, pt.NamedArray):
+                return expr.copy(tags=frozenset(),
+                                 axes=(pt.Axis(frozenset()),)*expr.ndim)
+            else:
+                return expr
+
+        dag = pt.transform.map_and_copy(dag, untag_loopy_call_results)
+
+        # }}}
+
+        return dag
+
+    def transform_loopy_program(self, t_unit):
+        # FIXME: Do not parallelize for now.
+        return t_unit
+
+# }}}
+
+
+# {{{ pytest actx factory
+
+class PytestPyOpenCLArrayContextFactory(
+        _PytestPyOpenCLArrayContextFactoryWithClass):
+    actx_class = PyOpenCLArrayContext
+
+
+class PytestPytatoPyOpenCLArrayContextFactory(
+        _PytestPytatoPyOpenCLArrayContextFactory):
+
+    @property
+    def actx_class(self):
+        return PytatoPyOpenCLArrayContext
+
+
+register_pytest_array_context_factory("meshmode.pyopencl",
+        PytestPyOpenCLArrayContextFactory)
+register_pytest_array_context_factory("meshmode.pytato_cl",
+        PytestPytatoPyOpenCLArrayContextFactory)
+
+# }}}
+
+
+# {{{ handle move deprecation
+
+_actx_names = (
+        "ArrayContext",
+
+        "CommonSubexpressionTag",
+        "FirstAxisIsElementsTag",
+
+        "ArrayContainer",
+        "is_array_container", "is_array_container_type",
+        "serialize_container", "deserialize_container",
+        "get_container_context", "get_container_context_recursively",
+        "with_container_arithmetic",
+        "dataclass_array_container",
+
+        "map_array_container", "multimap_array_container",
+        "rec_map_array_container", "rec_multimap_array_container",
+        "mapped_over_array_containers",
+        "multimapped_over_array_containers",
+        "freeze",
+
+        "make_loopy_program",
+
+        "pytest_generate_tests_for_pyopencl_array_context"
+        )
+
+
+def __getattr__(name):
+    if name not in _actx_names:
+        raise AttributeError(name)
+
+    import arraycontext
+    result = getattr(arraycontext, name)
+
+    warn(f"meshmode.array_context.{name} is deprecated. "
+         f"Use arraycontext.{name} instead. "
+         f"meshmode.array_context.{name} will continue to work until 2022.",
+         DeprecationWarning, stacklevel=2)
+
+    return result
+
+# }}}
+
+
+# vim: foldmethod=marker

meshmode/discretization/connection/chained.py

+__copyright__ = """Copyright (C) 2018 Alexandru Fikl"""
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+from dataclasses import dataclass
+
+import numpy as np
+
+import modepy as mp
+
+from meshmode.discretization.connection.direct import DiscretizationConnection
+
+
+# {{{ chained discretization connection
+
+class ChainedDiscretizationConnection(DiscretizationConnection):
+    """Aggregates multiple :class:`DiscretizationConnection` instances
+    into a single one.
+
+    .. attribute:: connections
+    """
+
+    def __init__(self, connections, from_discr=None):
+        if connections:
+            if from_discr is not None:
+                assert from_discr is connections[0].from_discr
+            else:
+                from_discr = connections[0].from_discr
+            is_surjective = all(cnx.is_surjective for cnx in connections)
+            to_discr = connections[-1].to_discr
+        else:
+            if from_discr is None:
+                raise ValueError("connections may not be empty if from_discr "
+                        "is not specified")
+
+            to_discr = from_discr
+
+            # It's an identity
+            is_surjective = True
+
+        super().__init__(
+                from_discr, to_discr, is_surjective=is_surjective)
+
+        self.connections = connections
+
+    def __call__(self, ary):
+        for cnx in self.connections:
+            ary = cnx(ary)
+
+        return ary
+
+# }}}
+
+
+# {{{ flatten chained connection
+
+@dataclass(frozen=True)
+class _ConnectionBatchData:
+    result_unit_nodes: np.ndarray
+    from_group_index: int
+    to_element_face: int
+
+
+def _iterbatches(groups):
+    for igrp, grp in enumerate(groups):
+        for ibatch, batch in enumerate(grp.batches):
+            yield (igrp, ibatch), (grp, batch)
+
+
+def _build_element_lookup_table(actx, conn):
+    el_table = [np.full(g.nelements, -1, dtype=np.int64)
+                for g in conn.to_discr.groups]
+
+    for (igrp, _), (_, batch) in _iterbatches(conn.groups):
+        el_table[igrp][actx.to_numpy(batch.to_element_indices)] = \
+                actx.to_numpy(batch.from_element_indices)
+
+    return el_table
+
+
+def _build_new_group_table(from_conn, to_conn):
+    def find_batch(nodes, gtb):
+        for igrp, batches in enumerate(gtb):
+            for ibatch, batch in enumerate(batches):
+                if np.allclose(nodes, batch.result_unit_nodes):
+                    return (igrp, ibatch)
+        return (-1, -1)
+
+    nfrom_groups = len(from_conn.groups)
+    n_to_groups = len(to_conn.groups)
+
+    # construct a table from (old groups) -> (new groups)
+    # NOTE: we try to reduce the number of new groups and batches by matching
+    # the `result_unit_nodes` and only adding a new batch if necessary
+    grp_to_grp = {}
+    batch_info = [[] for i in range(nfrom_groups * n_to_groups)]
+    for (igrp, ibatch), (_, fbatch) in _iterbatches(from_conn.groups):
+        for (jgrp, jbatch), (_, tbatch) in _iterbatches(to_conn.groups):
+            # compute result_unit_nodes
+            ffgrp = from_conn.from_discr.groups[fbatch.from_group_index]
+            from_matrix = mp.resampling_matrix(
+                    ffgrp.basis_obj().functions,
+                    fbatch.result_unit_nodes,
+                    ffgrp.unit_nodes)
+            result_unit_nodes = from_matrix.dot(ffgrp.unit_nodes.T)
+
+            tfgrp = to_conn.from_discr.groups[tbatch.from_group_index]
+            to_matrix = mp.resampling_matrix(
+                    tfgrp.basis_obj().functions,
+                    tbatch.result_unit_nodes,
+                    tfgrp.unit_nodes)
+            result_unit_nodes = to_matrix.dot(result_unit_nodes).T
+
+            # find new (group, batch)
+            (igrp_new, ibatch_new) = find_batch(result_unit_nodes, batch_info)
+            if igrp_new < 0:
+                igrp_new = n_to_groups * igrp + jgrp
+                ibatch_new = len(batch_info[igrp_new])
+
+                batch_info[igrp_new].append(_ConnectionBatchData(
+                    from_group_index=fbatch.from_group_index,
+                    result_unit_nodes=result_unit_nodes,
+                    to_element_face=tbatch.to_element_face))
+
+            grp_to_grp[igrp, ibatch, jgrp, jbatch] = (igrp_new, ibatch_new)
+
+    return grp_to_grp, batch_info
+
+
+def _build_batches(actx, from_bins, to_bins, batch):
+    from meshmode.discretization.connection.direct import InterpolationBatch
+
+    def to_device(x):
+        return actx.freeze(actx.from_numpy(np.asarray(x)))
+
+    for ibatch, (from_bin, to_bin) in enumerate(zip(from_bins, to_bins, strict=True)):
+        yield InterpolationBatch(
+                from_group_index=batch[ibatch].from_group_index,
+                from_element_indices=to_device(from_bin),
+                to_element_indices=to_device(to_bin),
+                result_unit_nodes=batch[ibatch].result_unit_nodes,
+                to_element_face=batch[ibatch].to_element_face)
+
+
+def flatten_chained_connection(actx, connection):
+    """Collapse a connection into a direct connection.
+
+    If the given connection is already a
+    :class:`~meshmode.discretization.connection.DirectDiscretizationConnection`
+    nothing is done. However, if the connection is a
+    :class:`~meshmode.discretization.connection.ChainedDiscretizationConnection`,
+    a new direct connection is constructed that transports from
+    :attr:`~meshmode.discretization.connection.DiscretizationConnection.from_discr`
+    to
+    :attr:`~meshmode.discretization.connection.DiscretizationConnection.to_discr`.
+
+    The new direct connection will have a number of groups and batches that
+    is, at worse, the product of all the connections in the chain. For
+    example, if we consider a connection between a discretization and a
+    two-level refinement, both levels will have :math:`n` groups and
+    :math:`m + 1` batches per group, where :math:`m` is the number of
+    subdivisions of an element (exact number depends on implementation
+    details in
+    :func:`~meshmode.discretization.connection.make_refinement_connection`).
+    However, a direct connection from level :math:`0` to level :math:`2`
+    will have at worst :math:`n^2` groups and each group will have
+    :math:`(m + 1)^2` batches.
+
+    .. warning::
+
+        If a large number of connections is chained, the number of groups and
+        batches can become very large.
+
+    :arg actx: An instance of :class:`arraycontext.ArrayContext`.
+    :arg connection: An instance of
+        :class:`~meshmode.discretization.connection.DiscretizationConnection`.
+    :return: An instance of
+        :class:`~meshmode.discretization.connection.DirectDiscretizationConnection`.
+    """
+    from meshmode.discretization.connection import (
+        DirectDiscretizationConnection,
+        DiscretizationConnectionElementGroup,
+        IdentityDiscretizationConnection,
+        make_same_mesh_connection,
+    )
+
+    if not hasattr(connection, "connections"):
+        return connection
+
+    if not connection.connections:
+        assert connection.to_discr is connection.from_discr
+        return make_same_mesh_connection(actx, connection.to_discr,
+                                         connection.from_discr)
+
+    # recursively build direct connections
+    connections = connection.connections
+    direct_connections = []
+    for conn in connections:
+        direct_connections.append(flatten_chained_connection(actx, conn))
+
+    direct_connections = [conn for conn in direct_connections
+            if not isinstance(conn, IdentityDiscretizationConnection)]
+
+    # merge all the direct connections
+    from_conn = direct_connections[0]
+    for to_conn in direct_connections[1:]:
+        el_table = _build_element_lookup_table(actx, from_conn)
+        grp_to_grp, batch_info = _build_new_group_table(from_conn, to_conn)
+
+        # distribute the indices to new groups and batches
+        from_bins = [[np.empty(0, dtype=np.int64) for _ in g] for g in batch_info]
+        to_bins = [[np.empty(0, dtype=np.int64) for _ in g] for g in batch_info]
+
+        for (igrp, ibatch), (_, from_batch) in _iterbatches(from_conn.groups):
+            from_to_element_indices = actx.to_numpy(from_batch.to_element_indices)
+
+            for (jgrp, jbatch), (_, to_batch) in _iterbatches(to_conn.groups):
+                igrp_new, ibatch_new = grp_to_grp[igrp, ibatch, jgrp, jbatch]
+
+                jfrom = actx.to_numpy(to_batch.from_element_indices)
+                jto = actx.to_numpy(to_batch.to_element_indices)
+
+                mask = np.isin(jfrom, from_to_element_indices)
+                from_bins[igrp_new][ibatch_new] = \
+                    np.hstack([from_bins[igrp_new][ibatch_new],
+                               el_table[igrp][jfrom[mask]]])
+                to_bins[igrp_new][ibatch_new] = \
+                    np.hstack([to_bins[igrp_new][ibatch_new],
+                               jto[mask]])
+
+        # build new groups
+        groups = []
+        for igrp, (from_bin, to_bin) in enumerate(zip(from_bins, to_bins, strict=True)):
+            groups.append(DiscretizationConnectionElementGroup(
+                list(_build_batches(actx, from_bin, to_bin,
+                                    batch_info[igrp]))))
+
+        from_conn = DirectDiscretizationConnection(
+            from_discr=from_conn.from_discr,
+            to_discr=to_conn.to_discr,
+            groups=groups,
+            is_surjective=connection.is_surjective)
+
+    return from_conn
+
+# }}}
+
+
+# {{{ build chained resample matrix
+
+def make_full_resample_matrix(actx, connection):
+    """Build a dense matrix representing the discretization connection.
+
+    This is based on
+    :func:`~meshmode.discretization.connection.direct.make_direct_full_resample_matrix`.
+    If a chained connection is given, the matrix is constructed recursively
+    for each connection and multiplied left to right.
+
+    .. warning::
+
+        This method will be very slow, both in terms of speed and memory
+        usage, and should only be used for testing or if absolutely necessary.
+
+    :arg actx: a :class:`arraycontext.ArrayContext`.
+    :arg connection: a
+        :class:`~meshmode.discretization.connection.DiscretizationConnection`.
+    :return: a :class:`pyopencl.array.Array` of shape
+        `(connection.from_discr.ndofs, connection.to_discr.ndofs)`.
+    """
+    from meshmode.discretization.connection.direct import (
+        DirectDiscretizationConnection,
+        make_direct_full_resample_matrix,
+    )
+
+    if isinstance(connection, DirectDiscretizationConnection):
+        return make_direct_full_resample_matrix(actx, connection)
+
+    if not isinstance(connection, ChainedDiscretizationConnection):
+        raise TypeError("only 'ChainedDiscretizationConnection's are supported")
+
+    if not connection.connections:
+        result = np.eye(connection.to_discr.ndofs)
+        return actx.from_numpy(result)
+
+    acc = actx.to_numpy(
+            make_full_resample_matrix(actx, connection.connections[0]))
+    for conn in connection.connections[1:]:
+        resampler = actx.to_numpy(make_full_resample_matrix(actx, conn))
+        acc = resampler @ acc
+
+    return actx.from_numpy(acc)
+
+# }}}
+
+
+# vim: foldmethod=marker