fix: Bug with elastic wave propagator at higher orders

src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp

@@ -328,6 +328,51 @@ class Qk_Hexahedron_Lagrange_GaussLobatto final : public FiniteElementBase
                            real64 const (&X)[numNodes][3],
                            StackVariables const & stack,
                            real64 ( &gradN )[numNodes][3] );
+  /**
+   * @brief Calculate the shape functions derivatives wrt the physical
+   *   coordinates.
+   * @param q Index of the quadrature point.
+   * @param X Array containing the coordinates of the mesh corners.
+   * @param gradN Array to contain the shape function derivatives for all
+   *   support points at the coordinates of the quadrature point @p q.
+   * @return The determinant of the parent/physical transformation matrix.
+   */
+  GEOS_HOST_DEVICE
+  GEOS_FORCE_INLINE
+  static real64 calcGradNWithCorners( localIndex const q,
+                                      real64 const (&X)[8][3],
+                                      real64 ( &gradN )[numNodes][3] );
+  /**
+   * @brief Calculate the shape functions derivatives wrt the physical
+   *   coordinates at a single point.
+   * @param[in] coords The parent coordinates at which to evaluate the shape function value
+   * @param[in] X Array containing the coordinates of the mesh corners.
+   * @param[out] gradN Array to contain the shape function derivatives for all
+   *   support points at the coordinates of the quadrature point @p q.
+   * @return The determinant of the parent/physical transformation matrix.
+   */
+  GEOS_HOST_DEVICE
+  GEOS_FORCE_INLINE
+  static real64 calcGradNWithCorners( real64 const (&coords)[3],
+                                      real64 const (&X)[8][3],
+                                      real64 ( &gradN )[numNodes][3] );
+
+  /**
+   * @brief Calculate the shape functions derivatives wrt the physical
+   *   coordinates.
+   * @param q Index of the quadrature point.
+   * @param X Array containing the coordinates of the mesh corners.
+   * @param stack Variables allocated on the stack as filled by @ref setupStack.
+   * @param gradN Array to contain the shape function derivatives for all
+   *   support points at the coordinates of the quadrature point @p q.
+   * @return The determinant of the parent/physical transformation matrix.
+   */
+  GEOS_HOST_DEVICE
+  GEOS_FORCE_INLINE
+  static real64 calcGradNWithCorners( localIndex const q,
+                                      real64 const (&X)[8][3],
+                                      StackVariables const & stack,
+                                      real64 ( &gradN )[numNodes][3] );
 
   /**
    * @brief Calculate the integration weights for a quadrature point.
@@ -494,6 +539,22 @@ class Qk_Hexahedron_Lagrange_GaussLobatto final : public FiniteElementBase
   static void jacobianTransformation( real64 const (&coords)[3],
                                       real64 const (&X)[numNodes][3],
                                       real64 ( &J )[3][3] );
+
+  /**
+   * @brief Calculates the isoparametric "Jacobian" transformation
+   *   matrix/mapping from the parent space to the physical space at a single point.
+   *   Assumes that the coordinate of high-order nodes are given by trilinear
+   *   interpolation of the mesh corners.
+   * @param coords The parent coordinates at which to evaluate the shape function value
+   * @param X Array containing the coordinates of the mesh corners.
+   * @param J Array to store the Jacobian transformation.
+   */
+  GEOS_HOST_DEVICE
+  GEOS_FORCE_INLINE
+  static void jacobianTransformationWithCorners( real64 const (&coords)[3],
+                                                 real64 const (&X)[8][3],
+                                                 real64 ( &J )[3][3] );
+
   /**
    * @brief performs a trilinear interpolation to determine the real-world coordinates of a
    *   vertex
@@ -923,6 +984,57 @@ calcGradN( localIndex const q,
   return calcGradN( q, X, gradN );
 }
 
+template< typename GL_BASIS >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+real64
+Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::calcGradNWithCorners( localIndex const q,
+                                                                       real64 const (&X)[8][3],
+                                                                       real64 (& gradN)[numNodes][3] )
+{
+  int qa, qb, qc;
+  GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
+
+  real64 J[3][3] = {{0}};
+
+  jacobianTransformation( qa, qb, qc, X, J );
+
+  real64 const detJ = LvArray::tensorOps::invert< 3 >( J );
+
+  applyTransformationToParentGradients( q, J, gradN );
+
+  return detJ;
+}
+//*************************************************************************************************
+template< typename GL_BASIS >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+real64
+Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::calcGradNWithCorners( real64 const (&coords)[3],
+                                                                       real64 const (&X)[8][3],
+                                                                       real64 (& gradN)[numNodes][3] )
+{
+  real64 J[3][3] = {{0}};
+
+  jacobianTransformationWithCorners( coords, X, J );
+
+  real64 const detJ = LvArray::tensorOps::invert< 3 >( J );
+
+  applyTransformationToParentGradients( coords, J, gradN );
+
+  return detJ;
+}
+template< typename GL_BASIS >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+real64 Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
+calcGradNWithCorners( localIndex const q,
+                      real64 const (&X)[8][3],
+                      StackVariables const & GEOS_UNUSED_PARAM( stack ),
+                      real64 ( & gradN )[numNodes][3] )
+{
+  return calcGradN( q, X, gradN );
+}
 //*************************************************************************************************
 #if __GNUC__
 #pragma GCC diagnostic push
@@ -983,6 +1095,37 @@ jacobianTransformation( real64 const (&coords)[3],
   }, X, J );
 }
 
+template< typename GL_BASIS >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+void
+Qk_Hexahedron_Lagrange_GaussLobatto< GL_BASIS >::
+jacobianTransformationWithCorners( real64 const (&coords)[3],
+                                   real64 const (&X)[8][3],
+                                   real64 ( & J )[3][3] )
+{
+  supportLoop( coords, [] GEOS_HOST_DEVICE ( real64 const (&dNdXi)[3],
+                                             int const nodeIndex,
+                                             real64 const (&X)[8][3],
+                                             real64 (& J)[3][3] )
+  {
+    int qa, qb, qc;
+    GL_BASIS::TensorProduct3D::multiIndex( nodeIndex, qa, qb, qc );
+    real64 Xnode[3];
+    real64 alpha = ( GL_BASIS::parentSupportCoord( qa ) + 1.0 ) / 2.0;
+    real64 beta = ( GL_BASIS::parentSupportCoord( qb ) + 1.0 ) / 2.0;
+    real64 gamma = ( GL_BASIS::parentSupportCoord( qc ) + 1.0 ) / 2.0;
+    trilinearInterp( alpha, beta, gamma, X, Xnode );
+    for( int i = 0; i < 3; ++i )
+    {
+      for( int j = 0; j < 3; ++j )
+      {
+        J[i][j] = J[i][j] + dNdXi[ j ] * Xnode[i];
+      }
+    }
+  }, X, J );
+}
+
 template< typename GL_BASIS >
 GEOS_HOST_DEVICE
 GEOS_FORCE_INLINE

src/coreComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp

@@ -460,13 +460,13 @@ struct PreComputeSourcesAndReceivers
                                      sourceCoordinates[isrc][1],
                                      sourceCoordinates[isrc][2] };
 
-          real64 xLocal[numNodesPerElem][3];
+          real64 xLocal[8][3];
 
-          for( localIndex a=0; a< numNodesPerElem; ++a )
+          for( localIndex a = 0; a < 8; ++a )
           {
-            for( localIndex i=0; i<3; ++i )
+            for( localIndex i = 0; i < 3; ++i )
             {
-              xLocal[a][i] = baseNodeCoords( elemsToNodes( k, a ), i );
+              xLocal[a][i] = baseNodeCoords( baseElemsToNodes( k, a ), i );
             }
           }
 
@@ -496,7 +496,7 @@ struct PreComputeSourcesAndReceivers
             real64 N[numNodesPerElem];
             real64 gradN[numNodesPerElem][3];
             FE_TYPE::calcN( coordsOnRefElem, N );
-            FE_TYPE::calcGradN( coordsOnRefElem, xLocal, gradN );
+            FE_TYPE::calcGradNWithCorners( coordsOnRefElem, xLocal, gradN );
             R2SymTensor moment = sourceMoment;
             for( localIndex q=0; q< numNodesPerElem; ++q )
             {
@@ -593,13 +593,13 @@ struct PreComputeSourcesAndReceivers
           if( sampleFound && elemGhostRank[k] < 0 )
           {
             real64 coordsOnRefElem[3]{};
-            real64 xLocal[numNodesPerElem][3];
+            real64 xLocal[8][3];
 
-            for( localIndex a=0; a< numNodesPerElem; ++a )
+            for( localIndex a = 0; a < 8; ++a )
             {
-              for( localIndex i=0; i<3; ++i )
+              for( localIndex i=0; i < 3; ++i )
               {
-                xLocal[a][i] = baseNodeCoords( elemsToNodes( k, a ), i );
+                xLocal[a][i] = baseNodeCoords( baseElemsToNodes( k, a ), i );
               }
             }
 
@@ -610,7 +610,7 @@ struct PreComputeSourcesAndReceivers
             real64 N[numNodesPerElem];
             real64 gradN[numNodesPerElem][3];
             FE_TYPE::calcN( coordsOnRefElem, N );
-            FE_TYPE::calcGradN( coordsOnRefElem, xLocal, gradN );
+            FE_TYPE::calcGradNWithCorners( coordsOnRefElem, xLocal, gradN );
             for( localIndex a = 0; a < numNodesPerElem; ++a )
             {
               receiverNodeIds[ircv][iSample * numNodesPerElem + a] = elemsToNodes( k,