Merge pull request #3430 from jwpeterson/allow_nodal_pyramid_quadrature

Add flag that controls whether nodal pyramid quadrature is allowed

Author: jwpeterson

include/quadrature/quadrature.h

@@ -245,6 +245,19 @@ class QBase : public ReferenceCountedObject<QBase>
    */
   bool allow_rules_with_negative_weights;
 
+  /**
+   * The flag's value defaults to false so that one does not accidentally
+   * use a nodal quadrature rule on Pyramid elements, since evaluating the
+   * inverse element Jacobian (e.g. dphi) is not well-defined at the
+   * Pyramid apex because the element Jacobian is zero there.
+   *
+   * We do not want to completely prevent someone from using a nodal
+   * quadrature rule on Pyramids, however, since there are legitimate
+   * use cases (lumped mass matrix) so the flag can be set to true to
+   * override this behavior.
+   */
+  bool allow_nodal_pyramid_quadrature;
+
 protected:
 
 

src/geom/elem.C

@@ -2424,6 +2424,9 @@ bool Elem::has_invertible_map(Real /*tol*/) const
   FEMap fe_map;
   auto & jac = fe_map.get_jacobian();
 
+  // We have a separate check for is_singular_node() below, so in this
+  // case its "OK" to do nodal quadrature on pyramids.
+  qnodal.allow_nodal_pyramid_quadrature = true;
   qnodal.init(this->type());
   auto & qp = qnodal.get_points();
   libmesh_assert_equal_to(qp.size(), this->n_nodes());

src/quadrature/quadrature.C

@@ -27,6 +27,7 @@ namespace libMesh
 QBase::QBase(unsigned int d,
              Order o) :
   allow_rules_with_negative_weights(true),
+  allow_nodal_pyramid_quadrature(false),
   _dim(d),
   _order(o),
   _type(INVALID_ELEM),

src/quadrature/quadrature_nodal_3D.C

@@ -37,7 +37,9 @@ void QNodal::init_3D(const ElemType, unsigned int)
     case HEX8:
     case PYRAMID5:
       {
+        // Construct QTrap rule that matches our own nodal pyramid quadrature permissions
         QTrap rule(/*dim=*/3, /*ignored*/_order);
+        rule.allow_nodal_pyramid_quadrature = this->allow_nodal_pyramid_quadrature;
         rule.init(_type, /*ignored*/_p_level);
         _points.swap (rule.get_points());
         _weights.swap(rule.get_weights());
@@ -107,7 +109,9 @@ void QNodal::init_3D(const ElemType, unsigned int)
     case PYRAMID13:
     case PYRAMID14:
       {
+        // Construct QSimpson rule that matches our own nodal pyramid quadrature permissions
         QSimpson rule(/*dim=*/3, /*ignored*/_order);
+        rule.allow_nodal_pyramid_quadrature = this->allow_nodal_pyramid_quadrature;
         rule.init(_type, /*ignored*/_p_level);
         _points.swap (rule.get_points());
         _weights.swap(rule.get_weights());

src/quadrature/quadrature_simpson_3D.C

@@ -140,6 +140,11 @@ void QSimpson::init_3D(const ElemType, unsigned int)
     case PYRAMID13:
     case PYRAMID14:
       {
+        libmesh_error_msg_if(!allow_nodal_pyramid_quadrature,
+                             "Nodal quadrature on Pyramid elements is not allowed by default since\n"
+                             "the Jacobian of the inverse element map is not well-defined at the Pyramid apex.\n"
+                             "Set the QBase::allow_nodal_pyramid_quadrature flag to true to ignore skip this check.");
+
         _points.resize(14);
         _weights.resize(14);
 

src/quadrature/quadrature_trap_3D.C

@@ -92,6 +92,11 @@ void QTrap::init_3D(const ElemType, unsigned int)
     case PYRAMID13:
     case PYRAMID14:
       {
+        libmesh_error_msg_if(!allow_nodal_pyramid_quadrature,
+                             "Nodal quadrature on Pyramid elements is not allowed by default since\n"
+                             "the Jacobian of the inverse element map is not well-defined at the Pyramid apex.\n"
+                             "Set the QBase::allow_nodal_pyramid_quadrature flag to true to ignore skip this check.");
+
         _points.resize(5);
         _weights.resize(5);
 

tests/quadrature/quadrature_test.C

@@ -177,6 +177,12 @@ private:
     std::unique_ptr<QBase> qrule =
       QBase::build(qtype, dim, static_cast<Order>(order));
 
+    // We are testing integration on the reference elements here, so
+    // the element map is not relevant, and we can safely use nodal
+    // Pyramid quadrature.
+    if (elem_type == PYRAMID5 || elem_type == PYRAMID13 || elem_type == PYRAMID14 || elem_type == PYRAMID18)
+      qrule->allow_nodal_pyramid_quadrature = true;
+
     qrule->init (elem_type);
 
     const int max_y_order = dim>1 ? 0 : exactorder;