diff --git a/include/fe/hdiv_fe_transformation.h b/include/fe/hdiv_fe_transformation.h index 8712250df32..ba06f9bf658 100644 --- a/include/fe/hdiv_fe_transformation.h +++ b/include/fe/hdiv_fe_transformation.h @@ -74,19 +74,21 @@ class HDivFETransformation : public FETransformationBase /** * Evaluates shape function gradients in physical coordinates for - * \f$ H(div) \f$ conforming elements. + * \f$ H(div) \f$ conforming elements, via the contravariant Piola + * map's derivative (the gradient of \f$ \phi = J^{-1} (dx/d\xi) + * \hat{\phi} \f$ with respect to physical coordinates). Requires + * \p LIBMESH_ENABLE_SECOND_DERIVATIVES, since the map's own second + * derivatives are needed to differentiate \f$ J^{-1} \f$ and + * \f$ (dx/d\xi) \f$. */ - virtual void map_dphi(const unsigned int /*dim*/, - const Elem * const /*elem*/, - const std::vector & /*qp*/, - const FEGenericBase & /*fe*/, - std::vector::OutputGradient>> & /*dphi*/, - std::vector> & /*dphidx*/, - std::vector> & /*dphidy*/, - std::vector> & /*dphidz*/) const override - { - libmesh_warning("WARNING: Shape function gradients for HDiv elements are not currently being computed!"); - } + virtual void map_dphi(const unsigned int dim, + const Elem * const elem, + const std::vector & qp, + const FEGenericBase & fe, + std::vector::OutputGradient>> & dphi, + std::vector> & dphidx, + std::vector> & dphidy, + std::vector> & dphidz) const override; #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES /** diff --git a/src/fe/hdiv_fe_transformation.C b/src/fe/hdiv_fe_transformation.C index 889a9938f72..89c232b8081 100644 --- a/src/fe/hdiv_fe_transformation.C +++ b/src/fe/hdiv_fe_transformation.C @@ -19,6 +19,44 @@ #include "libmesh/fe_interface.h" #include "libmesh/int_range.h" +namespace { + +using namespace libMesh; + +// The contravariant Piola map phi = J^{-1} * (dx/dxi) * phihat, shared by +// map_phi and map_dphi (which also needs the physical shape function value +// for the term coming from the derivative of J^{-1}). +template +OutputShape hdiv_piola_map(const RealGradient & dx_dxi, + const RealGradient & dx_deta, + Real J, + const OutputShape & phi_ref) +{ + OutputShape phi_val; + phi_val(0) = (dx_dxi(0)*phi_ref(0) + dx_deta(0)*phi_ref(1))/J; + phi_val(1) = (dx_dxi(1)*phi_ref(0) + dx_deta(1)*phi_ref(1))/J; +#if LIBMESH_DIM > 2 + phi_val(2) = (dx_dxi(2)*phi_ref(0) + dx_deta(2)*phi_ref(1))/J; +#endif + return phi_val; +} + +template +OutputShape hdiv_piola_map(const RealGradient & dx_dxi, + const RealGradient & dx_deta, + const RealGradient & dx_dzeta, + Real J, + const OutputShape & phi_ref) +{ + OutputShape phi_val; + phi_val(0) = (dx_dxi(0)*phi_ref(0) + dx_deta(0)*phi_ref(1) + dx_dzeta(0)*phi_ref(2))/J; + phi_val(1) = (dx_dxi(1)*phi_ref(0) + dx_deta(1)*phi_ref(1) + dx_dzeta(1)*phi_ref(2))/J; + phi_val(2) = (dx_dxi(2)*phi_ref(0) + dx_deta(2)*phi_ref(1) + dx_dzeta(2)*phi_ref(2))/J; + return phi_val; +} + +} // anonymous namespace + namespace libMesh { @@ -34,6 +72,9 @@ template void HDivFETransformation::init_map_dphi(const FEGenericBase & fe) const { fe.get_fe_map().get_dxidx(); +#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES + fe.get_fe_map().get_d2xyzdxi2(); +#endif } @@ -74,24 +115,13 @@ void HDivFETransformation::map_phi(const unsigned int dim, for (auto i : index_range(phi)) for (auto p : index_range(phi[i])) { - Real dx_dxi = dxyz_dxi[p](0); - Real dx_deta = dxyz_deta[p](0); - - Real dy_dxi = dxyz_dxi[p](1); - Real dy_deta = dxyz_deta[p](1); - - Real dz_dxi = dxyz_dxi[p](2); - Real dz_deta = dxyz_deta[p](2); - // Need to temporarily cache reference shape functions // We are computing mapping basis functions, so we explicitly ignore // any non-zero p_level() the Elem might have. OutputShape phi_ref; FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref); - phi[i][p](0) = (dx_dxi*phi_ref(0) + dx_deta*phi_ref(1))/J[p]; - phi[i][p](1) = (dy_dxi*phi_ref(0) + dy_deta*phi_ref(1))/J[p]; - phi[i][p](2) = (dz_dxi*phi_ref(0) + dz_deta*phi_ref(1))/J[p]; + phi[i][p] = hdiv_piola_map(dxyz_dxi[p], dxyz_deta[p], J[p], phi_ref); } break; @@ -108,27 +138,229 @@ void HDivFETransformation::map_phi(const unsigned int dim, for (auto i : index_range(phi)) for (auto p : index_range(phi[i])) { - Real dx_dxi = dxyz_dxi[p](0); - Real dx_deta = dxyz_deta[p](0); - Real dx_dzeta = dxyz_dzeta[p](0); + // Need to temporarily cache reference shape functions + // We are computing mapping basis functions, so we explicitly ignore + // any non-zero p_level() the Elem might have. + OutputShape phi_ref; + FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref); + + phi[i][p] = hdiv_piola_map(dxyz_dxi[p], dxyz_deta[p], dxyz_dzeta[p], J[p], phi_ref); + } - Real dy_dxi = dxyz_dxi[p](1); - Real dy_deta = dxyz_deta[p](1); - Real dy_dzeta = dxyz_dzeta[p](1); + break; + } - Real dz_dxi = dxyz_dxi[p](2); - Real dz_deta = dxyz_deta[p](2); - Real dz_dzeta = dxyz_dzeta[p](2); + default: + libmesh_error_msg("Invalid dim = " << dim); + } // switch(dim) +} +template +void HDivFETransformation::map_dphi(const unsigned int dim, + const Elem * const elem, + const std::vector & qp, + const FEGenericBase & fe, + std::vector::OutputGradient>> & dphi, + std::vector> & dphidx, + std::vector> & dphidy, + std::vector> & dphidz) const +{ +#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES + switch (dim) + { + case 0: + case 1: + libmesh_error_msg("These element transformations only make sense in 2D and 3D."); + + case 2: + { + const std::vector & dxyz_dxi = fe.get_fe_map().get_dxyzdxi(); + const std::vector & dxyz_deta = fe.get_fe_map().get_dxyzdeta(); + + const std::vector & d2xyz_dxi2 = fe.get_fe_map().get_d2xyzdxi2(); + const std::vector & d2xyz_deta2 = fe.get_fe_map().get_d2xyzdeta2(); + const std::vector & d2xyz_dxideta = fe.get_fe_map().get_d2xyzdxideta(); + + const std::vector & J = fe.get_fe_map().get_jacobian(); + + const std::vector & dxi_dx = fe.get_fe_map().get_dxidx(); + const std::vector & dxi_dy = fe.get_fe_map().get_dxidy(); + const std::vector & deta_dx = fe.get_fe_map().get_detadx(); + const std::vector & deta_dy = fe.get_fe_map().get_detady(); +#if LIBMESH_DIM > 2 + const std::vector & dxi_dz = fe.get_fe_map().get_dxidz(); + const std::vector & deta_dz = fe.get_fe_map().get_detadz(); +#endif + + const std::vector> & dphi_dxi = fe.get_dphidxi(); + const std::vector> & dphi_deta = fe.get_dphideta(); + + for (auto i : index_range(dphi)) + for (auto p : index_range(dphi[i])) + { // Need to temporarily cache reference shape functions // We are computing mapping basis functions, so we explicitly ignore // any non-zero p_level() the Elem might have. OutputShape phi_ref; FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref); - phi[i][p](0) = (dx_dxi*phi_ref(0) + dx_deta*phi_ref(1) + dx_dzeta*phi_ref(2))/J[p]; - phi[i][p](1) = (dy_dxi*phi_ref(0) + dy_deta*phi_ref(1) + dy_dzeta*phi_ref(2))/J[p]; - phi[i][p](2) = (dz_dxi*phi_ref(0) + dz_deta*phi_ref(1) + dz_dzeta*phi_ref(2))/J[p]; + const OutputShape phi_val = + hdiv_piola_map(dxyz_dxi[p], dxyz_deta[p], J[p], phi_ref); + + // Jacobi's formula: dJ/dxi_n = J * tr(F^{-1} * dF/dxi_n) + Real dJ_dxi = + J[p] * (dxi_dx[p]*d2xyz_dxi2[p](0) + deta_dx[p]*d2xyz_dxideta[p](0) + + dxi_dy[p]*d2xyz_dxi2[p](1) + deta_dy[p]*d2xyz_dxideta[p](1)); + Real dJ_deta = + J[p] * (dxi_dx[p]*d2xyz_dxideta[p](0) + deta_dx[p]*d2xyz_deta2[p](0) + + dxi_dy[p]*d2xyz_dxideta[p](1) + deta_dy[p]*d2xyz_deta2[p](1)); +#if LIBMESH_DIM > 2 + dJ_dxi += J[p] * (dxi_dz[p]*d2xyz_dxi2[p](2) + deta_dz[p]*d2xyz_dxideta[p](2)); + dJ_deta += J[p] * (dxi_dz[p]*d2xyz_dxideta[p](2) + deta_dz[p]*d2xyz_deta2[p](2)); +#endif + + for (unsigned int k = 0; k < LIBMESH_DIM; ++k) + { + // dphi_k/dx_l = A + B + C, where (n, m summed over reference directions): + // A = -phi_k(x) * SUM_n (dxi_n/dx_l) * (1/J) * (dJ/dxi_n) + // B = (1/J) * SUM_n (dxi_n/dx_l) * SUM_m (d^2 x_k/dxi_m dxi_n) * phihat_m + // C = (1/J) * SUM_n (dxi_n/dx_l) * SUM_m F_{km} * dphihat_m/dxi_n + + // Term A: -phi_k(x) * SUM_n (dxi_n/dx_l)*(1/J)*(dJ/dxi_n) + const Real phik = (k == 0) ? phi_val(0) : (k == 1) ? phi_val(1) : phi_val(2); + + // Term B: (1/J) * SUM_n (dxi_n/dx_l) * SUM_m (d^2 x_k/dxi_m dxi_n) * phihat_m + const Real d2x_k_dxi2 = (k == 0) ? d2xyz_dxi2[p](0) : (k == 1) ? d2xyz_dxi2[p](1) : d2xyz_dxi2[p](2); + const Real d2x_k_deta2 = (k == 0) ? d2xyz_deta2[p](0) : (k == 1) ? d2xyz_deta2[p](1) : d2xyz_deta2[p](2); + const Real d2x_k_dxideta = (k == 0) ? d2xyz_dxideta[p](0) : (k == 1) ? d2xyz_dxideta[p](1) : d2xyz_dxideta[p](2); + + const Real dphik_dxi_B = (d2x_k_dxi2*phi_ref(0) + d2x_k_dxideta*phi_ref(1))/J[p]; + const Real dphik_deta_B = (d2x_k_dxideta*phi_ref(0) + d2x_k_deta2*phi_ref(1))/J[p]; + + // Term C: (1/J) * SUM_n (dxi_n/dx_l) * SUM_m F_{km} * dphihat_m/dxi_n + const Real Fk_xi = (k == 0) ? dxyz_dxi[p](0) : (k == 1) ? dxyz_dxi[p](1) : dxyz_dxi[p](2); + const Real Fk_eta = (k == 0) ? dxyz_deta[p](0) : (k == 1) ? dxyz_deta[p](1) : dxyz_deta[p](2); + + const Real dphik_dxi_C = (Fk_xi*dphi_dxi[i][p](0) + Fk_eta*dphi_dxi[i][p](1))/J[p]; + const Real dphik_deta_C = (Fk_xi*dphi_deta[i][p](0) + Fk_eta*dphi_deta[i][p](1))/J[p]; + + const Real dphik_dxi_total = -phik*dJ_dxi/J[p] + dphik_dxi_B + dphik_dxi_C; + const Real dphik_deta_total = -phik*dJ_deta/J[p] + dphik_deta_B + dphik_deta_C; + + dphidx[i][p](k) = dxi_dx[p]*dphik_dxi_total + deta_dx[p]*dphik_deta_total; + dphidy[i][p](k) = dxi_dy[p]*dphik_dxi_total + deta_dy[p]*dphik_deta_total; +#if LIBMESH_DIM > 2 + dphidz[i][p](k) = dxi_dz[p]*dphik_dxi_total + deta_dz[p]*dphik_deta_total; +#endif + } + + dphi[i][p].slice(0) = dphidx[i][p]; + dphi[i][p].slice(1) = dphidy[i][p]; +#if LIBMESH_DIM > 2 + dphi[i][p].slice(2) = dphidz[i][p]; +#endif + } + + break; + } + case 3: + { + const std::vector & dxyz_dxi = fe.get_fe_map().get_dxyzdxi(); + const std::vector & dxyz_deta = fe.get_fe_map().get_dxyzdeta(); + const std::vector & dxyz_dzeta = fe.get_fe_map().get_dxyzdzeta(); + + const std::vector & d2xyz_dxi2 = fe.get_fe_map().get_d2xyzdxi2(); + const std::vector & d2xyz_deta2 = fe.get_fe_map().get_d2xyzdeta2(); + const std::vector & d2xyz_dzeta2 = fe.get_fe_map().get_d2xyzdzeta2(); + const std::vector & d2xyz_dxideta = fe.get_fe_map().get_d2xyzdxideta(); + const std::vector & d2xyz_dxidzeta = fe.get_fe_map().get_d2xyzdxidzeta(); + const std::vector & d2xyz_detadzeta = fe.get_fe_map().get_d2xyzdetadzeta(); + + const std::vector & J = fe.get_fe_map().get_jacobian(); + + const std::vector & dxi_dx = fe.get_fe_map().get_dxidx(); + const std::vector & dxi_dy = fe.get_fe_map().get_dxidy(); + const std::vector & dxi_dz = fe.get_fe_map().get_dxidz(); + const std::vector & deta_dx = fe.get_fe_map().get_detadx(); + const std::vector & deta_dy = fe.get_fe_map().get_detady(); + const std::vector & deta_dz = fe.get_fe_map().get_detadz(); + const std::vector & dzeta_dx = fe.get_fe_map().get_dzetadx(); + const std::vector & dzeta_dy = fe.get_fe_map().get_dzetady(); + const std::vector & dzeta_dz = fe.get_fe_map().get_dzetadz(); + + const std::vector> & dphi_dxi = fe.get_dphidxi(); + const std::vector> & dphi_deta = fe.get_dphideta(); + const std::vector> & dphi_dzeta = fe.get_dphidzeta(); + + for (auto i : index_range(dphi)) + for (auto p : index_range(dphi[i])) + { + // Need to temporarily cache reference shape functions + // We are computing mapping basis functions, so we explicitly ignore + // any non-zero p_level() the Elem might have. + OutputShape phi_ref; + FEInterface::shape(fe.get_fe_type(), /*extra_order=*/0, elem, i, qp[p], phi_ref); + + const OutputShape phi_val = + hdiv_piola_map(dxyz_dxi[p], dxyz_deta[p], dxyz_dzeta[p], J[p], phi_ref); + + // Jacobi's formula: dJ/dxi_n = J * tr(F^{-1} * dF/dxi_n) + const Real dJ_dxi = + J[p] * (dxi_dx[p]*d2xyz_dxi2[p](0) + deta_dx[p]*d2xyz_dxideta[p](0) + dzeta_dx[p]*d2xyz_dxidzeta[p](0) + + dxi_dy[p]*d2xyz_dxi2[p](1) + deta_dy[p]*d2xyz_dxideta[p](1) + dzeta_dy[p]*d2xyz_dxidzeta[p](1) + + dxi_dz[p]*d2xyz_dxi2[p](2) + deta_dz[p]*d2xyz_dxideta[p](2) + dzeta_dz[p]*d2xyz_dxidzeta[p](2)); + const Real dJ_deta = + J[p] * (dxi_dx[p]*d2xyz_dxideta[p](0) + deta_dx[p]*d2xyz_deta2[p](0) + dzeta_dx[p]*d2xyz_detadzeta[p](0) + + dxi_dy[p]*d2xyz_dxideta[p](1) + deta_dy[p]*d2xyz_deta2[p](1) + dzeta_dy[p]*d2xyz_detadzeta[p](1) + + dxi_dz[p]*d2xyz_dxideta[p](2) + deta_dz[p]*d2xyz_deta2[p](2) + dzeta_dz[p]*d2xyz_detadzeta[p](2)); + const Real dJ_dzeta = + J[p] * (dxi_dx[p]*d2xyz_dxidzeta[p](0) + deta_dx[p]*d2xyz_detadzeta[p](0) + dzeta_dx[p]*d2xyz_dzeta2[p](0) + + dxi_dy[p]*d2xyz_dxidzeta[p](1) + deta_dy[p]*d2xyz_detadzeta[p](1) + dzeta_dy[p]*d2xyz_dzeta2[p](1) + + dxi_dz[p]*d2xyz_dxidzeta[p](2) + deta_dz[p]*d2xyz_detadzeta[p](2) + dzeta_dz[p]*d2xyz_dzeta2[p](2)); + + for (unsigned int k = 0; k < 3; ++k) + { + // dphi_k/dx_l = A + B + C, where (n, m summed over reference directions): + // A = -phi_k(x) * SUM_n (dxi_n/dx_l) * (1/J) * (dJ/dxi_n) + // B = (1/J) * SUM_n (dxi_n/dx_l) * SUM_m (d^2 x_k/dxi_m dxi_n) * phihat_m + // C = (1/J) * SUM_n (dxi_n/dx_l) * SUM_m F_{km} * dphihat_m/dxi_n + + // Term A: -phi_k(x) * SUM_n (dxi_n/dx_l)*(1/J)*(dJ/dxi_n) + const Real phik = (k == 0) ? phi_val(0) : (k == 1) ? phi_val(1) : phi_val(2); + + // Term B: (1/J) * SUM_n (dxi_n/dx_l) * SUM_m (d^2 x_k/dxi_m dxi_n) * phihat_m + const Real d2x_k_dxi2 = (k == 0) ? d2xyz_dxi2[p](0) : (k == 1) ? d2xyz_dxi2[p](1) : d2xyz_dxi2[p](2); + const Real d2x_k_deta2 = (k == 0) ? d2xyz_deta2[p](0) : (k == 1) ? d2xyz_deta2[p](1) : d2xyz_deta2[p](2); + const Real d2x_k_dzeta2 = (k == 0) ? d2xyz_dzeta2[p](0) : (k == 1) ? d2xyz_dzeta2[p](1) : d2xyz_dzeta2[p](2); + const Real d2x_k_dxideta = (k == 0) ? d2xyz_dxideta[p](0) : (k == 1) ? d2xyz_dxideta[p](1) : d2xyz_dxideta[p](2); + const Real d2x_k_dxidzeta = (k == 0) ? d2xyz_dxidzeta[p](0) : (k == 1) ? d2xyz_dxidzeta[p](1) : d2xyz_dxidzeta[p](2); + const Real d2x_k_detadzeta = (k == 0) ? d2xyz_detadzeta[p](0) : (k == 1) ? d2xyz_detadzeta[p](1) : d2xyz_detadzeta[p](2); + + const Real dphik_dxi_B = (d2x_k_dxi2*phi_ref(0) + d2x_k_dxideta*phi_ref(1) + d2x_k_dxidzeta*phi_ref(2))/J[p]; + const Real dphik_deta_B = (d2x_k_dxideta*phi_ref(0) + d2x_k_deta2*phi_ref(1) + d2x_k_detadzeta*phi_ref(2))/J[p]; + const Real dphik_dzeta_B = (d2x_k_dxidzeta*phi_ref(0) + d2x_k_detadzeta*phi_ref(1) + d2x_k_dzeta2*phi_ref(2))/J[p]; + + // Term C: (1/J) * SUM_n (dxi_n/dx_l) * SUM_m F_{km} * dphihat_m/dxi_n + const Real Fk_xi = (k == 0) ? dxyz_dxi[p](0) : (k == 1) ? dxyz_dxi[p](1) : dxyz_dxi[p](2); + const Real Fk_eta = (k == 0) ? dxyz_deta[p](0) : (k == 1) ? dxyz_deta[p](1) : dxyz_deta[p](2); + const Real Fk_zeta = (k == 0) ? dxyz_dzeta[p](0) : (k == 1) ? dxyz_dzeta[p](1) : dxyz_dzeta[p](2); + + const Real dphik_dxi_C = (Fk_xi*dphi_dxi[i][p](0) + Fk_eta*dphi_dxi[i][p](1) + Fk_zeta*dphi_dxi[i][p](2))/J[p]; + const Real dphik_deta_C = (Fk_xi*dphi_deta[i][p](0) + Fk_eta*dphi_deta[i][p](1) + Fk_zeta*dphi_deta[i][p](2))/J[p]; + const Real dphik_dzeta_C = (Fk_xi*dphi_dzeta[i][p](0) + Fk_eta*dphi_dzeta[i][p](1) + Fk_zeta*dphi_dzeta[i][p](2))/J[p]; + + const Real dphik_dxi_total = -phik*dJ_dxi/J[p] + dphik_dxi_B + dphik_dxi_C; + const Real dphik_deta_total = -phik*dJ_deta/J[p] + dphik_deta_B + dphik_deta_C; + const Real dphik_dzeta_total = -phik*dJ_dzeta/J[p] + dphik_dzeta_B + dphik_dzeta_C; + + dphidx[i][p](k) = dxi_dx[p]*dphik_dxi_total + deta_dx[p]*dphik_deta_total + dzeta_dx[p]*dphik_dzeta_total; + dphidy[i][p](k) = dxi_dy[p]*dphik_dxi_total + deta_dy[p]*dphik_deta_total + dzeta_dy[p]*dphik_dzeta_total; + dphidz[i][p](k) = dxi_dz[p]*dphik_dxi_total + deta_dz[p]*dphik_deta_total + dzeta_dz[p]*dphik_dzeta_total; + } + + dphi[i][p].slice(0) = dphidx[i][p]; + dphi[i][p].slice(1) = dphidy[i][p]; + dphi[i][p].slice(2) = dphidz[i][p]; } break; @@ -137,6 +369,11 @@ void HDivFETransformation::map_phi(const unsigned int dim, default: libmesh_error_msg("Invalid dim = " << dim); } // switch(dim) +#else + libmesh_ignore(dim, elem, qp, fe, dphi, dphidx, dphidy, dphidz); + libmesh_error_msg("HDiv shape function gradients require the library to be configured " + "with --enable-second-derivatives (LIBMESH_ENABLE_SECOND_DERIVATIVES)."); +#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES } template @@ -222,6 +459,19 @@ void HDivFETransformation::map_phi(const unsigned int, libmesh_error_msg("HDiv transformations only make sense for vector-valued elements."); } +template<> +void HDivFETransformation::map_dphi(const unsigned int, + const Elem * const, + const std::vector &, + const FEGenericBase &, + std::vector::OutputGradient>> &, + std::vector> &, + std::vector> &, + std::vector> &) const +{ + libmesh_error_msg("HDiv transformations only make sense for vector-valued elements."); +} + template<> void HDivFETransformation::map_div(const unsigned int, const Elem * const, diff --git a/src/mesh/unstructured_mesh.C b/src/mesh/unstructured_mesh.C index fa581d14545..5cff72f1097 100644 --- a/src/mesh/unstructured_mesh.C +++ b/src/mesh/unstructured_mesh.C @@ -119,8 +119,10 @@ void transfer_elem(Elem & lo_elem, const unsigned int hon_end = hi_elem->n_nodes(); libmesh_assert_less (hon_begin, hon_end); +#ifdef LIBMESH_ENABLE_UNIQUE_ID libmesh_assert_less_equal (hon_end-hon_begin, max_new_nodes_per_elem); +#endif for (unsigned int hon=hon_begin; honget_dphi() against a central-difference approximation built from + * map_phi's own output (fe->get_phi()), evaluated at physical points + * perturbed via FEMap::inverse_map. + */ +class HDivGradTest : public CppUnit::TestCase +{ +public: + LIBMESH_CPPUNIT_TEST_SUITE(HDivGradTest); + +#if LIBMESH_DIM > 1 + CPPUNIT_TEST(testQuad8Affine); + CPPUNIT_TEST(testQuad8NonAffine); + CPPUNIT_TEST(testTri6Affine); + CPPUNIT_TEST(testTri6NonAffine); +#endif +#if LIBMESH_DIM > 2 + CPPUNIT_TEST(testHex27Affine); + CPPUNIT_TEST(testHex27NonAffine); + CPPUNIT_TEST(testTet14Affine); + CPPUNIT_TEST(testTet14NonAffine); +#endif + + CPPUNIT_TEST_SUITE_END(); + +public: + void setUp() {} + void tearDown() {} + +#if LIBMESH_DIM > 1 + // RAVIART_THOMAS dof placement piggybacks on the second-order edge/face + // nodes, so the geometric element types must be the second-order ones + // (QUAD8, TRI6, ...) even though the RT approximation order used here + // is FIRST. + void testQuad8Affine() { testGradByFiniteDifference(QUAD8, false); } + void testQuad8NonAffine() { testGradByFiniteDifference(QUAD8, true); } + void testTri6Affine() { testGradByFiniteDifference(TRI6, false); } + void testTri6NonAffine() { testGradByFiniteDifference(TRI6, true); } +#endif + +#if LIBMESH_DIM > 2 + void testHex27Affine() { testGradByFiniteDifference(HEX27, false); } + void testHex27NonAffine() { testGradByFiniteDifference(HEX27, true); } + void testTet14Affine() { testGradByFiniteDifference(TET14, false); } + void testTet14NonAffine() { testGradByFiniteDifference(TET14, true); } +#endif + +private: + // Builds a single-element mesh of the given type, optionally distorting + // one node so the element's map is non-affine, then checks the HDiv + // shape function gradients returned by fe->get_dphi() at every + // quadrature point and coordinate direction. + void testGradByFiniteDifference(const ElemType elem_type, const bool distort) + { + LOG_UNIT_TEST; + + Mesh mesh(*TestCommWorld); + unsigned int dim = 0; + + switch (elem_type) + { + case QUAD8: + case TRI6: + dim = 2; + MeshTools::Generation::build_square + (mesh, 1, 1, 0., 1., 0., 1., elem_type); + break; + case HEX27: + case TET14: + dim = 3; + MeshTools::Generation::build_cube + (mesh, 1, 1, 1, 0., 1., 0., 1., 0., 1., elem_type); + break; + default: + libmesh_error_msg("Unsupported element type " << elem_type); + } + + if (distort) + { + // Move a single vertex node without touching the associated + // edge/face nodes, turning the element's boundary curved and + // giving the map nonzero second derivatives -- this exercises + // the terms in map_dphi (A and B in the derivation) that vanish + // identically on an affine element. + Node & node = mesh.node_ref(0); + for (auto d : make_range(dim)) + node(d) += Real(0.1) * (d + 1); + } + + auto elem_range = mesh.active_local_element_ptr_range(); + if (elem_range.begin() == elem_range.end()) + return; + const Elem * elem = *(elem_range.begin()); + + const FEType fe_type(FIRST, RAVIART_THOMAS); + std::unique_ptr fe(FEVectorBase::build(dim, fe_type)); + + const std::vector> & dphi = fe->get_dphi(); + const std::vector> & phi = fe->get_phi(); + const std::vector & xyz = fe->get_xyz(); + + QGauss qrule(dim, fe_type.default_quadrature_order()); + fe->attach_quadrature_rule(&qrule); + + fe->reinit(elem); + + // Copy out what we need before the perturbed reinit() calls below + // overwrite fe's internal state. + const std::vector> dphi_base = dphi; + const std::vector qpoints = xyz; + + const Real h = 1e-6; + const Real tol = 5e-4; + + for (auto qp : index_range(qpoints)) + { + const Point x0 = qpoints[qp]; + + for (auto l : make_range(dim)) + { + Point x_plus = x0; + x_plus(l) += h; + Point x_minus = x0; + x_minus(l) -= h; + + std::vector pts_plus + (1, FEMap::inverse_map(dim, elem, x_plus)); + fe->reinit(elem, &pts_plus); + const std::vector> phi_plus = phi; + + std::vector pts_minus + (1, FEMap::inverse_map(dim, elem, x_minus)); + fe->reinit(elem, &pts_minus); + const std::vector> phi_minus = phi; + + for (auto i : index_range(dphi_base)) + for (auto k : make_range(LIBMESH_DIM)) + { + const Real fd = + (phi_plus[i][0](k) - phi_minus[i][0](k)) / (2*h); + + LIBMESH_ASSERT_FP_EQUAL(fd, dphi_base[i][qp](k, l), tol); + } + } + } + } +}; + +CPPUNIT_TEST_SUITE_REGISTRATION(HDivGradTest);