Repair the dual basis on QUAD8 and TRI6 faces

src/fe/fe_base.C

@@ -821,6 +821,28 @@ void FEGenericBase<Real>::compute_dual_shape_coeffs (const std::vector<Real> & J
         A(i,j) += JxW[qp]*phi_vals[i][qp]*phi_vals[j][qp];
     }
 
+  // D(k,k) = \int N_k is not positive on the serendipity quadratic faces: 0 at a TRI6 vertex, which
+  // leaves those dual shape functions identically zero, and -1/3 at a QUAD8 corner. Biorthogonalize
+  // instead against the locally quadratic transformed basis Ntilde = T N of Popp et al., SIAM J. Sci.
+  // Comput. 34(4):B421-B446, 2012, Sec. 4.4.1, in which each vertex absorbs a fraction alpha of its
+  // adjacent mid-edge shapes. That replaces D by T^-1 diag(T d), which for this T is the sparse
+  // update below. alpha = 1/5 is recommended there, and makes the weights strictly positive (QUAD8
+  // 1/5 and 4/5, TRI6 1/15 and 1/10) while preserving the partition of unity. T = I elsewhere.
+  if (_elem && (_elem->type() == TRI6 || _elem->type() == QUAD8) &&
+      get_family() == LAGRANGE && sz == _elem->n_nodes())
+    {
+      const Real alpha = Real(1)/5;
+      // Mid-edge nodes are the trailing indices, and only vertex entries are written, so D(m,m) here
+      // is never a value an earlier iteration modified.
+      for (const auto m : make_range(_elem->n_vertices(), sz))
+        for (const auto v : make_range(_elem->n_second_order_adjacent_vertices(m)))
+          {
+            const auto vertex = _elem->second_order_adjacent_vertex(m, v);
+            D(vertex, vertex) += alpha*D(m,m);
+            D(vertex, m)      -= alpha*D(m,m);
+          }
+    }
+
   // dual_coeff = A^-1*D
   for (const auto j : index_range(phi_vals))
   {

tests/fe/dual_shape_verification_test.C

@@ -1,5 +1,6 @@
 // libmesh includes
 #include "libmesh/libmesh.h"
+#include "libmesh/dense_matrix.h"
 #include "libmesh/edge_edge2.h"
 #include "libmesh/fe.h"
 #include "libmesh/quadrature_gauss.h"
@@ -21,6 +22,8 @@ class DualShapeTest : public CppUnit::TestCase
 public:
   LIBMESH_CPPUNIT_TEST_SUITE( DualShapeTest );
   CPPUNIT_TEST( testEdge2Lagrange );
+  CPPUNIT_TEST( testQuad8Lagrange );
+  CPPUNIT_TEST( testTri6Lagrange );
   CPPUNIT_TEST_SUITE_END();
 
 private:
@@ -71,6 +74,107 @@ public:
         my_tol);
   }
 
+  /**
+   * Check what the transformed dual basis on QUAD8/TRI6 guarantees, none of which the untransformed
+   * construction gives: strictly positive weights, reproduction of constants, and biorthogonality
+   * against Ntilde = T N rather than N. \p vertex_over_mid pins down alpha.
+   */
+  void testTransformedDual (const ElemType elem_type, const Real vertex_over_mid)
+  {
+    Mesh mesh(*TestCommWorld);
+    MeshTools::Generation::build_square(mesh, 1, 1, -1., 1., -1., 1., elem_type);
+
+    auto rng = mesh.active_local_element_ptr_range();
+    if (rng.begin() == rng.end())
+      return;
+    const Elem * elem = *(rng.begin());
+
+    FEType fe_type(SECOND, LAGRANGE);
+    std::unique_ptr<FEBase> fe = FEBase::build(2, fe_type);
+    const auto & JxW = fe->get_JxW();
+    const auto & phi = fe->get_phi();
+    const auto & dual_phi = fe->get_dual_phi();
+    QGauss qrule(2, fe_type.default_quadrature_order());
+    fe->attach_quadrature_rule(&qrule);
+    fe->reinit(elem);
+
+    const unsigned int n = elem->n_nodes();
+    CPPUNIT_ASSERT_EQUAL(std::size_t(n), phi.size());
+    CPPUNIT_ASSERT_EQUAL(std::size_t(n), dual_phi.size());
+
+    // TOLERANCE*TOLERANCE works with double but not float128
+    const Real my_tol = TOLERANCE*std::sqrt(TOLERANCE);
+
+    // dtilde_j = \int dual_phi_j; these are the weights that must be positive.
+    std::vector<Real> dtilde(n, 0);
+    for (const auto j : make_range(n))
+      for (const auto qp : index_range(JxW))
+        dtilde[j] += JxW[qp]*dual_phi[j][qp];
+
+    Real sum_dtilde = 0, volume = 0;
+    for (const auto j : make_range(n))
+      {
+        CPPUNIT_ASSERT(dtilde[j] > 0);
+        sum_dtilde += dtilde[j];
+      }
+    for (const auto qp : index_range(JxW))
+      volume += JxW[qp];
+
+    // Partition of unity is preserved, so constants are reproduced and the weights sum to the volume.
+    LIBMESH_ASSERT_FP_EQUAL(volume, sum_dtilde, my_tol);
+    for (const auto qp : index_range(JxW))
+      {
+        Real sum = 0;
+        for (const auto j : make_range(n))
+          sum += dual_phi[j][qp];
+        LIBMESH_ASSERT_FP_EQUAL(1., sum, my_tol);
+      }
+
+    // Scale-free signature of alpha.
+    LIBMESH_ASSERT_FP_EQUAL(vertex_over_mid,
+                            dtilde[0]/dtilde[elem->n_vertices()], my_tol);
+
+    // T from the element's own adjacency, not an assumed index convention.
+    const Real alpha = Real(1)/5;
+    DenseMatrix<Real> T(n, n);
+    for (const auto i : make_range(n))
+      T(i,i) = 1.;
+    for (const auto m : make_range(elem->n_vertices(), n))
+      {
+        T(m,m) = 1. - 2.*alpha;
+        for (const auto v : make_range(elem->n_second_order_adjacent_vertices(m)))
+          T(elem->second_order_adjacent_vertex(m,v), m) += alpha;
+      }
+
+    for (const auto k : make_range(n))
+      for (const auto j : make_range(n))
+        {
+          Real integral = 0;
+          for (const auto qp : index_range(JxW))
+            {
+              Real ntilde_k = 0;
+              for (const auto l : make_range(n))
+                ntilde_k += T(k,l)*phi[l][qp];
+              integral += JxW[qp]*dual_phi[j][qp]*ntilde_k;
+            }
+          LIBMESH_ASSERT_FP_EQUAL(k == j ? dtilde[j] : 0., integral, my_tol);
+        }
+  }
+
+  void testQuad8Lagrange ()
+  {
+    LOG_UNIT_TEST;
+    // QUAD8 weights 1/5 (corner) and 4/5 (mid-edge) on the reference element
+    testTransformedDual(QUAD8, Real(1)/4);
+  }
+
+  void testTri6Lagrange ()
+  {
+    LOG_UNIT_TEST;
+    // TRI6 weights 1/15 (vertex) and 1/10 (mid-edge) on the reference element
+    testTransformedDual(TRI6, Real(2)/3);
+  }
+
   void setUp()
     {
       FEType fe_type(FIRST, LAGRANGE);