From ed30a2dcea81105c38877363288482053ad2886a Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Sun, 25 Oct 2015 14:34:43 -0500
Subject: [PATCH] Detect and make use of permutation-matrix resampling matrices
 in connection

---
 meshmode/discretization/connection.py | 127 +++++++++++++++++++++++---
 test/test_meshmode.py                 |  21 +++--
 2 files changed, 130 insertions(+), 18 deletions(-)

diff --git a/meshmode/discretization/connection.py b/meshmode/discretization/connection.py
index 0eb5568a..fec50f0e 100644
--- a/meshmode/discretization/connection.py
+++ b/meshmode/discretization/connection.py
@@ -120,6 +120,11 @@ class DiscretizationConnection(object):
         a list of :class:`MeshConnectionGroup` instances, with
         a one-to-one correspondence to the groups in
         :attr:`to_discr`.
+
+    .. automethod:: __call__
+
+    .. automethod:: full_resample_matrix
+
     """
 
     def __init__(self, from_discr, to_discr, groups):
@@ -129,6 +134,16 @@ class DiscretizationConnection(object):
 
         self.cl_context = from_discr.cl_context
 
+        if from_discr.mesh.vertex_id_dtype != to_discr.mesh.vertex_id_dtype:
+            raise ValueError("from_discr and to_discr must agree on the "
+                    "vertex_id_dtype")
+
+        if from_discr.mesh.element_id_dtype != to_discr.mesh.element_id_dtype:
+            raise ValueError("from_discr and to_discr must agree on the "
+                    "element_id_dtype")
+
+        self.cl_context = from_discr.cl_context
+
         self.from_discr = from_discr
         self.to_discr = to_discr
         self.groups = groups
@@ -139,11 +154,60 @@ class DiscretizationConnection(object):
         ibatch = self.groups[to_group_index].batches[ibatch_index]
         from_grp = self.from_discr.groups[ibatch.from_group_index]
 
-        return mp.resampling_matrix(
+        result = mp.resampling_matrix(
                 mp.simplex_onb(self.from_discr.dim, from_grp.order),
                 ibatch.result_unit_nodes, from_grp.unit_nodes)
 
+        with cl.CommandQueue(self.cl_context) as queue:
+            return cl.array.to_device(queue, result).with_queue(None)
+
+    @memoize_method
+    def _resample_point_pick_indices(self, to_group_index, ibatch_index,
+            tol_multiplier=None):
+        """If :meth:`_resample_matrix` *R* is a row subset of a permutation matrix *P*,
+        return the index subset I so that, loosely, ``x[I] == R @ x``.
+
+        Will return *None* if no such index array exists, or a
+        :class:`pyopencl.array.Array` containing the index subset.
+        """
+
+        with cl.CommandQueue(self.cl_context) as queue:
+            mat = self._resample_matrix(to_group_index, ibatch_index).get(
+                    queue=queue)
+
+        nrows, ncols = mat.shape
+        result = np.zeros(nrows, dtype=self.to_discr.mesh.element_id_dtype)
+
+        if tol_multiplier is None:
+            tol_multiplier = 50
+
+        tol = np.finfo(mat.dtype).eps * tol_multiplier
+
+        for irow in range(nrows):
+            one_indices, = np.where(np.abs(mat[irow] - 1) < tol)
+            zero_indices, = np.where(np.abs(mat[irow]) < tol)
+
+            if len(one_indices) != 1:
+                return None
+            if len(zero_indices) != ncols - 1:
+                return None
+
+            one_index, = one_indices
+            result[irow] = one_index
+
+        with cl.CommandQueue(self.cl_context) as queue:
+            return cl.array.to_device(queue, result).with_queue(None)
+
     def full_resample_matrix(self, queue):
+        """Build a dense matrix representing this discretization connection.
+
+        .. warning::
+
+            On average, this will be exceedingly expensive (:math:`O(N^2)` in
+            the number *N* of discretization points) in terms of memory usage
+            and thus not what you'd typically want.
+        """
+
         @memoize_in(self, "oversample_mat_knl")
         def knl():
             import loopy as lp
@@ -193,8 +257,8 @@ class DiscretizationConnection(object):
         return result
 
     def __call__(self, queue, vec):
-        @memoize_in(self, "oversample_knl")
-        def knl():
+        @memoize_in(self, "resample_by_mat_knl")
+        def mat_knl():
             import loopy as lp
             knl = lp.make_kernel(
                 """{[k,i,j]:
@@ -215,7 +279,33 @@ class DiscretizationConnection(object):
                     lp.ValueArg("nelements_vec", np.int32),
                     "...",
                     ],
-                name="oversample")
+                name="resample_by_mat")
+
+            knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
+            return lp.tag_inames(knl, dict(k="g.0"))
+
+        @memoize_in(self, "resample_by_picking_knl")
+        def pick_knl():
+            import loopy as lp
+            knl = lp.make_kernel(
+                """{[k,i,j]:
+                    0<=k<nelements and
+                    0<=i<n_to_nodes}""",
+                "result[to_element_indices[k], i] \
+                    = vec[from_element_indices[k], pick_list[i]]",
+                [
+                    lp.GlobalArg("result", None,
+                        shape="nelements_result, n_to_nodes",
+                        offset=lp.auto),
+                    lp.GlobalArg("vec", None,
+                        shape="nelements_vec, n_from_nodes",
+                        offset=lp.auto),
+                    lp.ValueArg("nelements_result", np.int32),
+                    lp.ValueArg("nelements_vec", np.int32),
+                    lp.ValueArg("n_from_nodes", np.int32),
+                    "...",
+                    ],
+                name="resample_by_picking")
 
             knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
             return lp.tag_inames(knl, dict(k="g.0"))
@@ -228,17 +318,28 @@ class DiscretizationConnection(object):
         if vec.shape != (self.from_discr.nnodes,):
             raise ValueError("invalid shape of incoming resampling data")
 
-        for i_grp, (tgrp, sgrp, cgrp) in enumerate(
+        for i_tgrp, (tgrp, sgrp, cgrp) in enumerate(
                 zip(self.to_discr.groups, self.from_discr.groups, self.groups)):
             for i_batch, batch in enumerate(cgrp.batches):
                 if not len(batch.from_element_indices):
                     continue
 
-                knl()(queue,
-                        resample_mat=self._resample_matrix(i_grp, i_batch),
-                        result=tgrp.view(result), vec=sgrp.view(vec),
-                        from_element_indices=batch.from_element_indices,
-                        to_element_indices=batch.to_element_indices)
+                point_pick_indices = self._resample_point_pick_indices(
+                        i_tgrp, i_batch)
+
+                if point_pick_indices is None:
+                    mat_knl()(queue,
+                            resample_mat=self._resample_matrix(i_tgrp, i_batch),
+                            result=tgrp.view(result), vec=sgrp.view(vec),
+                            from_element_indices=batch.from_element_indices,
+                            to_element_indices=batch.to_element_indices)
+
+                else:
+                    pick_knl()(queue,
+                            pick_list=point_pick_indices,
+                            result=tgrp.view(result), vec=sgrp.view(vec),
+                            from_element_indices=batch.from_element_indices,
+                            to_element_indices=batch.to_element_indices)
 
         return result
 
@@ -318,12 +419,12 @@ def _build_boundary_connection(queue, vol_discr, bdry_discr, connection_data):
             vol_discr, bdry_discr, connection_groups)
 
 
+# {{{ pull together boundary vertices
+
 def _get_face_vertices(mesh, boundary_tag):
     # a set of volume vertex numbers
     bdry_vertex_vol_nrs = set()
 
-    # {{{ pull together boundary vertices
-
     if boundary_tag is not None:
         # {{{ boundary faces
 
@@ -357,6 +458,8 @@ def _get_face_vertices(mesh, boundary_tag):
 
         return np.arange(mesh.nvertices, dtype=np.intp)
 
+# }}}
+
 
 def make_face_restriction(queue, discr, group_factory, boundary_tag):
     """Create a mesh, a discretization and a connection to restrict
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 55602e0f..86009d7a 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -33,6 +33,11 @@ from pyopencl.tools import (  # noqa
         pytest_generate_tests_for_pyopencl
         as pytest_generate_tests)
 
+from meshmode.discretization.poly_element import (
+        InterpolatoryQuadratureSimplexGroupFactory,
+        PolynomialWarpAndBlendGroupFactory
+        )
+
 import pytest
 
 import logging
@@ -62,15 +67,17 @@ def test_circle_mesh(do_plot=False):
         pt.show()
 
 
-def test_boundary_interpolation(ctx_getter):
+@pytest.mark.parametrize("group_factory", [
+    InterpolatoryQuadratureSimplexGroupFactory,
+    PolynomialWarpAndBlendGroupFactory
+    ])
+def test_boundary_interpolation(ctx_getter, group_factory):
     cl_ctx = ctx_getter()
     queue = cl.CommandQueue(cl_ctx)
 
     from meshmode.mesh.io import generate_gmsh, FileSource
     from meshmode.mesh import BTAG_ALL
     from meshmode.discretization import Discretization
-    from meshmode.discretization.poly_element import \
-            InterpolatoryQuadratureSimplexGroupFactory
     from meshmode.discretization.connection import make_face_restriction
 
     from pytools.convergence import EOCRecorder
@@ -89,7 +96,7 @@ def test_boundary_interpolation(ctx_getter):
         print("END GEN")
 
         vol_discr = Discretization(cl_ctx, mesh,
-                InterpolatoryQuadratureSimplexGroupFactory(order))
+                group_factory(order))
         print("h=%s -> %d elements" % (
                 h, sum(mgrp.nelements for mgrp in mesh.groups)))
 
@@ -97,7 +104,7 @@ def test_boundary_interpolation(ctx_getter):
         f = 0.1*cl.clmath.sin(30*x)
 
         bdry_mesh, bdry_discr, bdry_connection = make_face_restriction(
-                queue, vol_discr, InterpolatoryQuadratureSimplexGroupFactory(order),
+                queue, vol_discr, group_factory(order),
                 BTAG_ALL)
 
         bdry_x = bdry_discr.nodes()[0].with_queue(queue)
@@ -113,7 +120,9 @@ def test_boundary_interpolation(ctx_getter):
         eoc_rec.add_data_point(h, err)
 
     print(eoc_rec)
-    assert eoc_rec.order_estimate() >= order-0.5
+    assert (
+            eoc_rec.order_estimate() >= order-0.5
+            or eoc_rec.max_error() < 1e-14)
 
 
 def test_element_orientation():
-- 
GitLab