diff --git a/include/fe/fe_base.h b/include/fe/fe_base.h index 5afe0a3f62f..72ab71b970f 100644 --- a/include/fe/fe_base.h +++ b/include/fe/fe_base.h @@ -208,6 +208,27 @@ class FEGenericBase : public FEAbstract { libmesh_assert(!calculations_started || calculate_phi); calculate_phi = true; return phi; } + /** + * \returns The dual (biorthogonal) shape function values at the quadrature + * points on the element. + * + * The dual basis is constructed to satisfy + * \f$ \int \Phi_j N_k \, d\gamma = \delta_{jk} \int N_k \, d\gamma \f$ + * against the primal basis \f$ N \f$. That biorthogonality is what makes a + * mortar coupling matrix diagonal, so a Lagrange multiplier can be condensed + * with a diagonal inverse. + * + * \note On QUAD8 and TRI6 these functions are biorthogonal to a locally + * transformed primal basis rather than to \f$ N \f$ itself: \f$ \int N_k \f$ + * is not positive on those faces (exactly 0 at a TRI6 vertex, -1/3 at a QUAD8 + * corner), so the usual construction degenerates. For those two element types + * \f$ \int \Phi_j N_k \f$ is therefore not diagonal; it factors as + * \f$ D = \tilde{D} T^{-1} \f$ with \f$ \tilde{D} \f$ diagonal, leaving + * \f$ D^{-1} = T \tilde{D}^{-1} \f$ sparse and cheap to apply. The functions + * returned here are still expressed in the primal basis, so the type and + * meaning of this return value are unchanged. See Popp, Wohlmuth, Gee and + * Wall, SIAM J. Sci. Comput. 34(4):B421-B446, 2012, Sec. 4.4.1 and Eq. (4.12). + */ const std::vector> & get_dual_phi() const { libmesh_assert(!calculations_started || calculate_dual); @@ -231,6 +252,11 @@ class FEGenericBase : public FEAbstract { libmesh_assert(!calculations_started || calculate_dphi); calculate_dphi = calculate_dphiref = true; return dphi; } + /** + * \returns The dual shape function derivatives at the quadrature points. + * + * \note See get_dual_phi() for the QUAD8/TRI6 caveat. + */ const std::vector> & get_dual_dphi() const { libmesh_assert(!calculations_started || calculate_dphi); calculate_dphi = calculate_dual = calculate_dphiref = true; return dual_dphi; } @@ -241,6 +267,11 @@ class FEGenericBase : public FEAbstract virtual void request_dual_dphi() const override { get_dual_dphi(); } + /** + * \returns The coefficients expressing the dual basis in the primal basis. + * + * \note See get_dual_phi() for the QUAD8/TRI6 caveat. + */ const DenseMatrix & get_dual_coeff() const { return dual_coeff; } @@ -320,6 +351,12 @@ class FEGenericBase : public FEAbstract { libmesh_assert(!calculations_started || calculate_d2phi); calculate_d2phi = calculate_dphiref = true; return d2phi; } + /** + * \returns The dual shape function second derivatives at the quadrature + * points. + * + * \note See get_dual_phi() for the QUAD8/TRI6 caveat. + */ const std::vector> & get_dual_d2phi() const { libmesh_assert(!calculations_started || calculate_d2phi); calculate_d2phi = calculate_dual = calculate_dphiref = true; return dual_d2phi; } diff --git a/src/fe/fe_base.C b/src/fe/fe_base.C index bea9ef46596..92dbd181f6a 100644 --- a/src/fe/fe_base.C +++ b/src/fe/fe_base.C @@ -821,6 +821,29 @@ void FEGenericBase::compute_dual_shape_coeffs (const std::vector & J A(i,j) += JxW[qp]*phi_vals[i][qp]*phi_vals[j][qp]; } + // The integral positivity condition on the dual basis (Popp, Wohlmuth, Gee and Wall, SIAM J. Sci. + // Comput. 34(4):B421-B446, 2012, Eq. (4.2)) fails on QUAD8 and TRI6: D(k,k) = \int N_k is -1/3 at a + // QUAD8 corner and exactly 0 at a TRI6 vertex (their Eqs. (4.3) and (4.4)), the latter leaving those + // dual shape functions identically zero. Biorthogonalize instead against the locally quadratic + // transformed basis Ntilde = T N of their 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)) { diff --git a/tests/fe/dual_shape_verification_test.C b/tests/fe/dual_shape_verification_test.C index ea4d4a123bd..72338dae225 100644 --- a/tests/fe/dual_shape_verification_test.C +++ b/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 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 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 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);