Distinguish between interpolatory quadrature and quadrature element groups

meshmode/discretization/poly_element.py

@@ -31,6 +31,24 @@ from meshmode.mesh import SimplexElementGroup as _MeshSimplexElementGroup
 import modepy as mp
 
 __doc__ = """
+
+Group types
+^^^^^^^^^^^
+
+.. autclass:: InterpolatoryQuadratureSimplexElementGroup
+.. autoclass:: QuadratureSimplexElementGroup
+.. autoclass:: PolynomialWarpAndBlendElementGroup
+
+Group factories
+^^^^^^^^^^^^^^^
+
+.. autclass:: InterpolatoryQuadratureSimplexGroupFactory
+.. autoclass:: QuadratureSimplexGroupFactory
+.. autoclass:: PolynomialWarpAndBlendGroupFactory
+
+Discretization class
+^^^^^^^^^^^^^^^^^^^^
+
 .. autoclass:: PolynomialElementDiscretization
 """
 
@@ -80,7 +98,14 @@ class PolynomialSimplexElementGroupBase(ElementGroupBase):
                 self.unit_nodes, meg.unit_nodes)
 
 
-class QuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase):
+class InterpolatoryQuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase):
+    """Elemental discretization supplying a high-order quadrature rule
+    with a number of nodes matching the number of polynomials in :math:`P^k`,
+    hence usable for differentiation and interpolation.
+
+    No interpolation nodes are present on the boundary of the simplex.
+    """
+
     @memoize_method
     def _quadrature_rule(self):
         dims = self.mesh_el_group.dim
@@ -106,7 +131,40 @@ class QuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase):
         return self._quadrature_rule().weights
 
 
+class QuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase):
+    """Elemental discretization supplying a high-order quadrature rule
+    with a number of nodes which does not necessarily match the number of
+    polynomials in :math:`P^k`. This discretization therefore excels at
+    quadarature, but is not necessarily usable for interpolation.
+
+    No interpolation nodes are present on the boundary of the simplex.
+    """
+
+    @memoize_method
+    def _quadrature_rule(self):
+        dims = self.mesh_el_group.dim
+        if dims == 1:
+            return mp.LegendreGaussQuadrature(self.order)
+        else:
+            return mp.XiaoGimbutasSimplexQuadrature(self.order, dims)
+
+    @property
+    @memoize_method
+    def unit_nodes(self):
+        return self._quadrature_rule().nodes
+
+    @property
+    @memoize_method
+    def weights(self):
+        return self._quadrature_rule().weights
+
+
 class PolynomialWarpAndBlendElementGroup(PolynomialSimplexElementGroupBase):
+    """Elemental discretization with a number of nodes matching the number of
+    polynomials in :math:`P^k`, hence usable for differentiation and
+    interpolation. Interpolation nodes are present on the boundary of the
+    simplex.
+    """
     @property
     @memoize_method
     def unit_nodes(self):
@@ -145,6 +203,11 @@ class OrderBasedGroupFactory(ElementGroupFactory):
         return self.group_class(mesh_el_group, self.order, node_nr_base)
 
 
+class InterpolatoryQuadratureSimplexGroupFactory(OrderBasedGroupFactory):
+    mesh_group_class = _MeshSimplexElementGroup
+    group_class = InterpolatoryQuadratureSimplexElementGroup
+
+
 class QuadratureSimplexGroupFactory(OrderBasedGroupFactory):
     mesh_group_class = _MeshSimplexElementGroup
     group_class = QuadratureSimplexElementGroup

test/test_meshmode.py

@@ -45,7 +45,7 @@ def test_boundary_interpolation(ctx_getter):
     from meshmode.mesh.io import generate_gmsh, FileSource
     from meshmode.discretization import Discretization
     from meshmode.discretization.poly_element import \
-            QuadratureSimplexGroupFactory
+            InterpolatoryQuadratureSimplexGroupFactory
     from meshmode.discretization.connection import make_boundary_restriction
 
     from pytools.convergence import EOCRecorder
@@ -64,7 +64,7 @@ def test_boundary_interpolation(ctx_getter):
         print "END GEN"
 
         vol_discr = Discretization(cl_ctx, mesh,
-                QuadratureSimplexGroupFactory(order))
+                InterpolatoryQuadratureSimplexGroupFactory(order))
         print "h=%s -> %d elements" % (
                 h, sum(mgrp.nelements for mgrp in mesh.groups))
 
@@ -72,7 +72,7 @@ def test_boundary_interpolation(ctx_getter):
         f = 0.1*cl.clmath.sin(30*x)
 
         bdry_mesh, bdry_discr, bdry_connection = make_boundary_restriction(
-                queue, vol_discr, QuadratureSimplexGroupFactory(order))
+                queue, vol_discr, InterpolatoryQuadratureSimplexGroupFactory(order))
 
         bdry_x = bdry_discr.nodes()[0].with_queue(queue)
         bdry_f = 0.1*cl.clmath.sin(30*bdry_x)
@@ -233,13 +233,13 @@ def test_sanity_balls(ctx_getter, src_file, dim, mesh_order,
 
         from meshmode.discretization import Discretization
         from meshmode.discretization.poly_element import \
-                QuadratureSimplexGroupFactory
+                InterpolatoryQuadratureSimplexGroupFactory
         vol_discr = Discretization(ctx, mesh,
-                QuadratureSimplexGroupFactory(quad_order))
+                InterpolatoryQuadratureSimplexGroupFactory(quad_order))
 
         from meshmode.discretization.connection import make_boundary_restriction
         bdry_mesh, bdry_discr, bdry_connection = make_boundary_restriction(
-                queue, vol_discr, QuadratureSimplexGroupFactory(quad_order))
+                queue, vol_discr, InterpolatoryQuadratureSimplexGroupFactory(quad_order))
 
         # }}}