Refactor towards more knowledge in element groups

doc/mesh.rst

 Mesh management
 ===============
 
-.. automodule:: pytential.mesh
+.. automodule:: meshmode.mesh
 
 .. autoclass:: MeshElementGroup
     :members:
@@ -14,16 +14,11 @@ Mesh management
 Mesh generation
 ---------------
 
-.. automodule:: pytential.mesh.generation
+.. automodule:: meshmode.mesh.generation
 
 Mesh input/output
 -----------------
 
-.. automodule:: pytential.mesh.io
-
-Mesh connections (for interpolation)
-------------------------------------
-
-.. automodule:: pytential.mesh.connection
+.. automodule:: meshmode.mesh.io
 
 .. vim: sw=4

meshmode/discretization/__init__.py

@@ -22,19 +22,96 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-
 import numpy as np
-from pytools import memoize_method
-
+from pytools import memoize_method, memoize_method_nested
+import loopy as lp
+import pyopencl as cl
 
 __doc__ = """
-
+.. autoclass:: ElementGroupBase
+.. autoclass:: ElementGroupFactory
 .. autoclass:: Discretization
 """
 
 
+# {{{ element group base
+
+class ElementGroupBase(object):
+    """Container for the :class:`QBXGrid` data corresponding to
+    one :class:`meshmode.mesh.MeshElementGroup`.
+
+    .. attribute :: mesh_el_group
+    .. attribute :: node_nr_base
+    """
+
+    def __init__(self, mesh_el_group, order, node_nr_base):
+        """
+        :arg mesh_el_group: an instance of
+            :class:`meshmode.mesh.MeshElementGroup`
+        """
+        self.mesh_el_group = mesh_el_group
+        self.order = order
+        self.node_nr_base = node_nr_base
+
+    @property
+    def nelements(self):
+        return self.mesh_el_group.nelements
+
+    @property
+    def nunit_nodes(self):
+        return self.unit_nodes.shape[-1]
+
+    @property
+    def nnodes(self):
+        return self.nunit_nodes * self.nelements
+
+    @property
+    def dim(self):
+        return self.mesh_el_group.dim
+
+    def _nodes(self):
+        # Not cached, because the global nodes array is what counts.
+        # This is just used to build that.
+
+        return np.tensordot(
+                self.mesh_el_group.nodes,
+                self._from_mesh_interp_matrix(),
+                (-1, -1))
+
+    def view(self, global_array):
+        return global_array[
+                ..., self.node_nr_base:self.node_nr_base + self.nnodes] \
+                .reshape(
+                        global_array.shape[:-1]
+                        + (self.nelements, self.nunit_nodes))
+
+# }}}
+
+
+# {{{ group factories
+
+class ElementGroupFactory(object):
+    """
+    .. function:: __call__(mesh_ele_group, node_nr_base)
+    """
+
+
+class OrderBasedGroupFactory(ElementGroupFactory):
+    def __init__(self, order):
+        self.order = order
+
+    def __call__(self, mesh_el_group, node_nr_base):
+        if not isinstance(mesh_el_group, self.mesh_group_class):
+            raise TypeError("only mesh element groups of type '%s' "
+                    "are supported" % self.mesh_group_class.__name__)
+
+        return self.group_class(mesh_el_group, self.order, node_nr_base)
+
+# }}}
+
+
 class Discretization(object):
-    """Abstract interface for discretizations.
+    """An unstructured composite discretization.
 
     .. attribute:: real_dtype
 
@@ -46,7 +123,11 @@ class Discretization(object):
 
     .. attribute:: ambient_dim
 
-    .. method::  empty(dtype, queue=None, extra_dims=None)
+    .. attribute :: nnodes
+
+    .. attribute :: groups
+
+    .. method:: empty(dtype, queue=None, extra_dims=None)
 
     .. method:: nodes()
 
@@ -57,51 +138,124 @@ class Discretization(object):
     .. method:: quad_weights(queue)
 
         shape: ``(nnodes)``
-
-    .. rubric:: Layer potential source discretizations only
-
-    .. method:: preprocess_optemplate(name, expr)
-
-    .. method:: op_group_features(expr)
-
-        Return a characteristic tuple by which operators that can be
-        executed together can be grouped.
-
-        *expr* is a subclass of
-        :class:`pymbolic.primitives.LayerPotentialOperatorBase`.
     """
 
-    @memoize_method
-    def _integral_op(self):
-        from pytential import sym, bind
-        return bind(self, sym.integral(sym.var("integrand")))
-
-    def integral(self, queue, x):
-        return self._integral_op()(queue, integrand=x)
-
-    @memoize_method
-    def _norm_op(self, num_components):
-        from pytential import sym, bind
-        if num_components is not None:
-            from pymbolic.primitives import make_sym_vector
-            v = make_sym_vector("integrand", num_components)
-            integrand = sym.real(np.dot(sym.conj(v), v))
+    def __init__(self, cl_ctx, mesh, group_factory, real_dtype=np.float64):
+        """
+        :arg order: A polynomial-order-like parameter passed unmodified to
+            :attr:`group_class`. See subclasses for more precise definition.
+        """
+
+        self.cl_context = cl_ctx
+
+        self.mesh = mesh
+        self.nnodes = 0
+        self.groups = []
+        for mg in mesh.groups:
+            ng = group_factory(mg, self.nnodes)
+            self.groups.append(ng)
+            self.nnodes += ng.nnodes
+
+        self.real_dtype = np.dtype(real_dtype)
+        self.complex_dtype = {
+                np.float32: np.complex64,
+                np.float64: np.complex128
+                }[self.real_dtype.type]
+
+    @property
+    def dim(self):
+        return self.mesh.dim
+
+    @property
+    def ambient_dim(self):
+        return self.mesh.ambient_dim
+
+    def empty(self, dtype, queue=None, extra_dims=None):
+        if queue is None:
+            first_arg = self.cl_context
         else:
-            integrand = sym.abs(sym.var("integrand"))**2
-
-        return bind(self, sym.integral(integrand))
+            first_arg = queue
+
+        shape = (self.nnodes,)
+        if extra_dims is not None:
+            shape = extra_dims + shape
+
+        return cl.array.empty(first_arg, shape, dtype=dtype)
+
+    def num_reference_derivative(
+            self, queue, ref_axes, vec):
+        @memoize_method_nested
+        def knl():
+            knl = lp.make_kernel(
+                """{[k,i,j]:
+                    0<=k<nelements and
+                    0<=i,j<ndiscr_nodes}""",
+                "result[k,i] = sum(j, diff_mat[i, j] * vec[k, j])",
+                default_offset=lp.auto, name="diff")
+
+            knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
+            return lp.tag_inames(knl, dict(k="g.0"))
+
+        result = self.empty(vec.dtype)
+
+        for grp in self.groups:
+            mat = None
+            for ref_axis in ref_axes:
+                next_mat = grp.diff_matrices()[ref_axis]
+                if mat is None:
+                    mat = next_mat
+                else:
+                    mat = np.dot(next_mat, mat)
+
+            knl()(queue, diff_mat=mat, result=grp.view(result), vec=grp.view(vec))
+
+        return result
+
+    def quad_weights(self, queue):
+        @memoize_method_nested
+        def knl():
+            knl = lp.make_kernel(
+                "{[k,i]: 0<=k<nelements and 0<=i<ndiscr_nodes}",
+                "result[k,i] = weights[i]",
+                name="quad_weights")
+
+            knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
+            return lp.tag_inames(knl, dict(k="g.0"))
+
+        result = self.empty(self.real_dtype)
+        for grp in self.groups:
+            knl()(queue, result=grp.view(result), weights=grp.weights)
+        return result
 
-    def norm(self, queue, x):
-        from pymbolic.geometric_algebra import MultiVector
-        if isinstance(x, MultiVector):
-            x = x.as_vector(np.object)
-
-        num_components = None
-        if isinstance(x, np.ndarray):
-            num_components, = x.shape
+    @memoize_method
+    def nodes(self):
+        @memoize_method_nested
+        def knl():
+            knl = lp.make_kernel(
+                """{[d,k,i,j]:
+                    0<=d<dims and
+                    0<=k<nelements and
+                    0<=i<ndiscr_nodes and
+                    0<=j<nmesh_nodes}""",
+                """
+                    result[d, k, i] = \
+                        sum(j, resampling_mat[i, j] * nodes[d, k, j])
+                    """,
+                name="nodes")
+
+            knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
+            return lp.tag_inames(knl, dict(k="g.0"))
+
+        result = self.empty(self.real_dtype, extra_dims=(self.ambient_dim,))
+
+        with cl.CommandQueue(self.cl_context) as queue:
+            for grp in self.groups:
+                meg = grp.mesh_el_group
+                knl()(queue,
+                        resampling_mat=grp.resampling_matrix(),
+                        result=grp.view(result), nodes=meg.nodes)
+
+        return result
 
-        norm_op = self._norm_op(num_components)
-        from math import sqrt
-        return sqrt(norm_op(queue, integrand=x))
 
 # vim: fdm=marker

meshmode/discretization/poly_element.py

@@ -27,6 +27,7 @@ import numpy as np
 #import numpy.linalg as la
 from pytools import memoize_method, memoize_method_nested
 from meshmode.discretization import Discretization
+from meshmode.mesh import SimplexElementGroup as _MeshSimplexElementGroup
 
 import pyopencl as cl
 
@@ -46,70 +47,44 @@ __doc__ = """
 # loopy sees strides that it doesn't expect and complains.
 
 
-# {{{ element group base
+from meshmode.discretization import ElementGroupBase
 
-class PolynomialElementGroupBase(object):
-    """Container for the :class:`QBXGrid` data corresponding to
-    one :class:`pytential.mesh.MeshElementGroup`.
 
-    .. attribute :: grid
-    .. attribute :: mesh_el_group
-    .. attribute :: node_nr_base
-    """
-
-    def __init__(self, grid, mesh_el_group, order, node_nr_base):
-        """
-        :arg grid: an instance of :class:`QBXGridBase`
-        :arg mesh_el_group: an instance of
-            :class:`pytential.mesh.MeshElementGroup`
-        """
-        self.grid = grid
-        self.mesh_el_group = mesh_el_group
-        self.order = order
-        self.node_nr_base = node_nr_base
-
-    @property
-    def nelements(self):
-        return self.mesh_el_group.nelements
+# {{{ concrete element groups
 
-    @property
-    def nunit_nodes(self):
-        return self.unit_nodes.shape[-1]
-
-    @property
-    def nnodes(self):
-        return self.nunit_nodes * self.nelements
+class PolynomialSimplexElementGroupBase(ElementGroupBase):
+    def basis(self):
+        return mp.simplex_onb(self.dim, self.order)
 
     @memoize_method
-    def _from_mesh_interp_matrix(self):
+    def from_mesh_interp_matrix(self):
         meg = self.mesh_el_group
         return mp.resampling_matrix(
                 mp.simplex_onb(meg.dim, meg.order),
                 self.unit_nodes,
                 meg.unit_nodes)
 
-    def _nodes(self):
-        # Not cached, because the global nodes array is what counts.
-        # This is just used to build that.
-
-        return np.tensordot(
-                self.mesh_el_group.nodes,
-                self._from_mesh_interp_matrix(),
-                (-1, -1))
-
-    def view(self, global_array):
-        return global_array[
-                ..., self.node_nr_base:self.node_nr_base + self.nnodes] \
-                .reshape(
-                        global_array.shape[:-1]
-                        + (self.nelements, self.nunit_nodes))
+    @memoize_method
+    def diff_matrices(self):
+        result = mp.differentiation_matrices(
+                self.basis(),
+                mp.grad_simplex_onb(self.dim, self.order),
+                self.unit_nodes)
 
-# }}}
+        if not isinstance(result, tuple):
+            return (result,)
+        else:
+            return result
 
+    @memoize_method
+    def resampling_matrix(self):
+        meg = self.mesh_el_group
+        return mp.resampling_matrix(
+                mp.simplex_onb(self.dim, meg.order),
+                self.unit_nodes, meg.unit_nodes)
 
-# {{{ concrete element groups
 
-class PolynomialQuadratureElementGroup(PolynomialElementGroupBase):
+class QuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase):
     @memoize_method
     def _quadrature_rule(self):
         dims = self.mesh_el_group.dim
@@ -129,7 +104,7 @@ class PolynomialQuadratureElementGroup(PolynomialElementGroupBase):
         return self._quadrature_rule().weights
 
 
-class PolynomialWarpAndBlendElementGroup(PolynomialElementGroupBase):
+class PolynomialWarpAndBlendElementGroup(PolynomialSimplexElementGroupBase):
     @property
     @memoize_method
     def unit_nodes(self):
@@ -144,165 +119,34 @@ class PolynomialWarpAndBlendElementGroup(PolynomialElementGroupBase):
 # }}}
 
 
-# {{{ discretization
-
-class PolynomialElementDiscretizationBase(Discretization):
-    """An (unstructured) composite polynomial discretization without
-    any specific opinion on how to evaluate layer potentials.
-
-    .. attribute :: mesh
-    .. attribute :: groups
-    .. attribute :: nnodes
-
-    .. autoattribute :: nodes
-    """
+# {{{ group factories
 
-    def __init__(self, cl_ctx, mesh, order, real_dtype=np.float64):
-        """
-        :arg order: A polynomial-order-like parameter passed unmodified to
-            :attr:`group_class`. See subclasses for more precise definition.
-        """
+class ElementGroupFactory(object):
+    pass
 
-        self.cl_context = cl_ctx
 
-        self.mesh = mesh
-        self.nnodes = 0
-        self.groups = []
-        for mg in mesh.groups:
-            ng = self.group_class(self, mg, order, self.nnodes)
-            self.groups.append(ng)
-            self.nnodes += ng.nnodes
-
-        self.real_dtype = np.dtype(real_dtype)
-        self.complex_dtype = {
-                np.float32: np.complex64,
-                np.float64: np.complex128
-                }[self.real_dtype.type]
-
-    @property
-    def dim(self):
-        return self.mesh.dim
-
-    @property
-    def ambient_dim(self):
-        return self.mesh.ambient_dim
-
-    def empty(self, dtype, queue=None, extra_dims=None):
-        if queue is None:
-            first_arg = self.cl_context
-        else:
-            first_arg = queue
-
-        shape = (self.nnodes,)
-        if extra_dims is not None:
-            shape = extra_dims + shape
+class OrderBasedGroupFactory(ElementGroupFactory):
+    def __init__(self, order):
+        self.order = order
 
-        return cl.array.empty(first_arg, shape, dtype=dtype)
+    def __call__(self, mesh_el_group, node_nr_base):
+        if not isinstance(mesh_el_group, self.mesh_group_class):
+            raise TypeError("only mesh element groups of type '%s' "
+                    "are supported" % self.mesh_group_class.__name__)
 
-    @memoize_method
-    def _diff_matrices(self, grp):
-        result = mp.differentiation_matrices(
-                mp.simplex_onb(self.dim, grp.order),
-                mp.grad_simplex_onb(self.dim, grp.order),
-                grp.unit_nodes)
+        return self.group_class(mesh_el_group, self.order, node_nr_base)
 
-        if not isinstance(result, tuple):
-            return (result,)
-        else:
-            return result
 
-    def num_reference_derivative(
-            self, queue, ref_axes, vec):
-        @memoize_method_nested
-        def knl():
-            knl = lp.make_kernel(
-                """{[k,i,j]:
-                    0<=k<nelements and
-                    0<=i,j<ndiscr_nodes}""",
-                "result[k,i] = sum(j, diff_mat[i, j] * vec[k, j])",
-                default_offset=lp.auto, name="diff")
-
-            knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
-            return lp.tag_inames(knl, dict(k="g.0"))
-
-        result = self.empty(vec.dtype)
-
-        for grp in self.groups:
-            mat = None
-            for ref_axis in ref_axes:
-                next_mat = self._diff_matrices(grp)[ref_axis]
-                if mat is None:
-                    mat = next_mat
-                else:
-                    mat = np.dot(next_mat, mat)
-
-            knl()(queue, diff_mat=mat, result=grp.view(result), vec=grp.view(vec))
-
-        return result
-
-    def quad_weights(self, queue):
-        @memoize_method_nested
-        def knl():
-            knl = lp.make_kernel(
-                "{[k,i]: 0<=k<nelements and 0<=i<ndiscr_nodes}",
-                "result[k,i] = weights[i]",
-                name="quad_weights")
-
-            knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
-            return lp.tag_inames(knl, dict(k="g.0"))
-
-        result = self.empty(self.real_dtype)
-        for grp in self.groups:
-            knl()(queue, result=grp.view(result), weights=grp.weights)
-        return result
+class QuadratureSimplexGroupFactory(OrderBasedGroupFactory):
+    mesh_group_class = _MeshSimplexElementGroup
+    group_class = QuadratureSimplexElementGroup
 
-    @memoize_method
-    def _resampling_matrix(self, grp):
-        meg = grp.mesh_el_group
-        return mp.resampling_matrix(
-                mp.simplex_onb(self.dim, meg.order),
-                grp.unit_nodes, meg.unit_nodes)
 
-    @memoize_method
-    def nodes(self):
-        @memoize_method_nested
-        def knl():
-            knl = lp.make_kernel(
-                """{[d,k,i,j]:
-                    0<=d<dims and
-                    0<=k<nelements and
-                    0<=i<ndiscr_nodes and
-                    0<=j<nmesh_nodes}""",
-                """
-                    result[d, k, i] = \
-                        sum(j, resampling_mat[i, j] * nodes[d, k, j])
-                    """,
-                name="nodes")
-
-            knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
-            return lp.tag_inames(knl, dict(k="g.0"))
-
-        result = self.empty(self.real_dtype, extra_dims=(self.ambient_dim,))
-
-        with cl.CommandQueue(self.cl_context) as queue:
-            for grp in self.groups:
-                meg = grp.mesh_el_group
-                knl()(queue,
-                        resampling_mat=self._resampling_matrix(grp),
-                        result=grp.view(result), nodes=meg.nodes)
-
-        return result
+class PolynomialWarpAndBlendGroupFactory(OrderBasedGroupFactory):
+    mesh_group_class = _MeshSimplexElementGroup
+    group_class = PolynomialWarpAndBlendElementGroup
 
 # }}}
 
 
-class PolynomialQuadratureElementDiscretization(
-        PolynomialElementDiscretizationBase):
-    group_class = PolynomialQuadratureElementGroup
-
-
-class PolynomialWarpAndBlendElementDiscretization(
-        PolynomialElementDiscretizationBase):
-    group_class = PolynomialWarpAndBlendElementGroup
-
 # vim: fdm=marker

meshmode/discretization/visualization.py

@@ -205,10 +205,12 @@ class Visualizer(object):
 
 
 def make_visualizer(queue, discr, vis_order):
+    from meshmode.discretization import Discretization
     from meshmode.discretization.poly_element import \
-            PolynomialWarpAndBlendElementDiscretization
-    vis_discr = PolynomialWarpAndBlendElementDiscretization(
-            discr.cl_context, discr.mesh, order=vis_order,
+            PolynomialWarpAndBlendGroupFactory
+    vis_discr = Discretization(
+            discr.cl_context, discr.mesh,
+            PolynomialWarpAndBlendGroupFactory(vis_order),
             real_dtype=discr.real_dtype)
     from meshmode.discretization.connection import \
             make_same_mesh_connection

meshmode/mesh/__init__.py

@@ -70,31 +70,12 @@ class MeshElementGroup(Record):
             The nodes are assumed to be mapped versions of *unit_nodes*.
         :arg unit_nodes: ``[dim, nunit_nodes]``
             The unit nodes of which *nodes* is a mapped
-            version. If unspecified, the nodes from
-            :func:`modepy.warp_and_blend_nodes` for *dim*
-            are assumed. These must be in unit coordinates
-            as defined in :mod:`modepy.nodes`.
-        :arg dim: only used if *unit_nodes* is None, to get
-            the default unit nodes.
+            version.
 
         Do not supply *element_nr_base* and *node_nr_base*, they will be
         automatically assigned.
         """
 
-        if unit_nodes is None:
-            if dim is None:
-                raise TypeError("'dim' must be passed "
-                        "if 'unit_nodes' is not passed")
-
-            unit_nodes = mp.warp_and_blend_nodes(dim, order)
-
-        dims = unit_nodes.shape[0]
-
-        if vertex_indices.shape[-1] != dims+1:
-            raise ValueError("vertex_indices has wrong number of vertices per "
-                    "element. expected: %d, got: %d" % (dims+1,
-                        vertex_indices.shape[-1]))
-
         Record.__init__(self,
             order=order,
             vertex_indices=vertex_indices,
@@ -127,6 +108,45 @@ class MeshElementGroup(Record):
     def nunit_nodes(self):
         return self.unit_nodes.shape[-1]
 
+
+class SimplexElementGroup(MeshElementGroup):
+    def __init__(self, order, vertex_indices, nodes,
+            element_nr_base=None, node_nr_base=None,
+            unit_nodes=None, dim=None):
+        """
+        :arg order: the mamximum total degree used for interpolation.
+        :arg nodes: ``[ambient_dim, nelements, nunit_nodes]``
+            The nodes are assumed to be mapped versions of *unit_nodes*.
+        :arg unit_nodes: ``[dim, nunit_nodes]``
+            The unit nodes of which *nodes* is a mapped
+            version. If unspecified, the nodes from
+            :func:`modepy.warp_and_blend_nodes` for *dim*
+            are assumed. These must be in unit coordinates
+            as defined in :mod:`modepy.nodes`.
+        :arg dim: only used if *unit_nodes* is None, to get
+            the default unit nodes.
+
+        Do not supply *element_nr_base* and *node_nr_base*, they will be
+        automatically assigned.
+        """
+
+        if unit_nodes is None:
+            if dim is None:
+                raise TypeError("'dim' must be passed "
+                        "if 'unit_nodes' is not passed")
+
+            unit_nodes = mp.warp_and_blend_nodes(dim, order)
+
+        dims = unit_nodes.shape[0]
+
+        if vertex_indices.shape[-1] != dims+1:
+            raise ValueError("vertex_indices has wrong number of vertices per "
+                    "element. expected: %d, got: %d" % (dims+1,
+                        vertex_indices.shape[-1]))
+
+        MeshElementGroup.__init__(self, order, vertex_indices, nodes,
+                element_nr_base, node_nr_base, unit_nodes, dim)
+
 # }}}
 
 

meshmode/mesh/generation.py

@@ -165,8 +165,8 @@ def make_curve_mesh(curve_f, element_boundaries, order):
     t = el_starts[:, np.newaxis] + el_lengths[:, np.newaxis]*nodes_01
     nodes = curve_f(t.ravel()).reshape(vertices.shape[0], nelements, -1)
 
-    from meshmode.mesh import Mesh, MeshElementGroup
-    egroup = MeshElementGroup(
+    from meshmode.mesh import Mesh, SimplexElementGroup
+    egroup = SimplexElementGroup(
             order,
             vertex_indices=np.vstack([
                 np.arange(nelements),

meshmode/mesh/io.py

@@ -122,7 +122,7 @@ class GmshMeshReceiver(GmshMeshReceiverBase):
 
         # }}}
 
-        from meshmode.mesh import MeshElementGroup, Mesh
+        from meshmode.mesh import SimplexElementGroup, Mesh
 
         for group_el_type, ngroup_elements in el_type_hist.iteritems():
             if group_el_type.dimensions != mesh_bulk_dim:
@@ -150,7 +150,7 @@ class GmshMeshReceiver(GmshMeshReceiverBase):
             unit_nodes = (np.array(group_el_type.lexicographic_node_tuples(),
                     dtype=np.float64).T/group_el_type.order)*2 - 1
 
-            groups.append(MeshElementGroup(
+            groups.append(SimplexElementGroup(
                 group_el_type.order,
                 vertex_indices,
                 nodes,