Intrepid2: add hierarchical bases for H(div) pyramids (#12286)

@trilinos/intrepid2 

## Motivation
This PR adds an arbitrary-order, hierarchical basis for H(div) on the pyramid.  This follows #12079, which added hierarchical bases for H(grad) on the pyramid, and #12136, which added bases for H(vol).

This PR also adds support in CellGeometry for a pyramidal subdivision of regular hexahedral meshes, and support for a rudimentary node numbering within CellGeometry.  This is in support of new projection convergence tests, so far implemented for H(grad) and H(vol).  H(div) support for such projections actually requires H(curl) bases to be implemented as well; those will be coming in a subsequent PR.

Additionally, this PR adds orientation tests against H(div) on the pyramid.

## Testing
This PR has orientation and basis cardinality tests for the new H(div) basis.  The new CellGeometry features have some limited tests against them as well.  It also adds projection convergence tests for H(grad) and H(vol) on the pyramid.

Additionally, I have done offline comparison with the ESEAS implementation of this basis, testing up to 9th order, with excellent agreement.  (There were actually a couple of bugs I found in ESEAS's implementation, which I plan to submit fixes for there shortly.)  I hope to include these tests with Intrepid2 soon; this has been prevented previously by the license for ESEAS, but they are changing the license to one that will allow inclusion.

packages/intrepid2/src/Cell/Intrepid2_CellGeometry.hpp

@@ -98,6 +98,7 @@ namespace Intrepid2
       , FOUR_TRIANGLES      /// square --> four triangles, with a new vertex at center
       , FIVE_TETRAHEDRA     /// cube   --> five tetrahedra
       , SIX_TETRAHEDRA      /// cube   --> six tetrahedra
+      , SIX_PYRAMIDS        /// cube   --> six pyramids
     };
 
     /*! Distinguish between "classic" (counterclockwise in 2D, counterclockwise bottom followed by counterclockwise top in 3D) Shards ordering and a more natural tensor ordering.  The tensor ordering is such that the lowest-order bit in the composite vertex number corresponds to the x dimension, the second-lowest-order bit corresponds to the y dimension, etc.  (There is not yet a non-"classic" Shards ordering, but we hope to add one.) */
@@ -222,6 +223,10 @@ namespace Intrepid2
     //! The shards CellTopology for each cell within the CellGeometry object.  Note that this is always a lowest-order CellTopology, even for higher-order curvilinear geometry.  Higher-order geometry for CellGeometry is expressed in terms of the basis returned by basisForNodes().
     const shards::CellTopology & cellTopology() const;
 
+    //! Returns a global node number corresponding to the specified (cell, node) combination.  If uniform grid (possibly subdivided), this number is computed dynamically; for more general meshes, this simply returns the result of a lookup in the cellToNodes container provided at construction.
+    KOKKOS_INLINE_FUNCTION
+    ordinal_type cellToNode(ordinal_type cellIndex, ordinal_type cellNodeNumber) const;
+
     //! Indicates the type of geometric variation from one cell to the next.  If all cells are known to have the same geometry, returns CONSTANT.
     //! Returns CONSTANT for uniform grids with no subdivisions, MODULAR for uniform grids with subdivisions, GENERAL for all others.
     KOKKOS_INLINE_FUNCTION

packages/intrepid2/src/Cell/Intrepid2_CellGeometryDef.hpp

@@ -173,6 +173,7 @@ namespace Intrepid2
         return 4;
       case FIVE_TETRAHEDRA:
         return 5;
+      case SIX_PYRAMIDS:
       case SIX_TETRAHEDRA:
         return 6;
     }
@@ -745,6 +746,69 @@ namespace Intrepid2
     return HostMemberLookup::getCellTopology(this);
   }
 
+  template<class PointScalar, int spaceDim, typename DeviceType>
+  KOKKOS_INLINE_FUNCTION
+  ordinal_type CellGeometry<PointScalar,spaceDim,DeviceType>::cellToNode(ordinal_type cell, ordinal_type node) const
+  {
+    if (cellToNodes_.is_allocated())
+    {
+      return cellToNodes_(cell,node);
+    }
+    else if (cellGeometryType_ == UNIFORM_GRID)
+    {
+      const ordinal_type numSubdivisions    = numCellsPerGridCell(subdivisionStrategy_);
+      const ordinal_type gridCell           = cell / numSubdivisions;
+      const ordinal_type subdivisionOrdinal = cell % numSubdivisions;
+
+      Kokkos::Array<ordinal_type,spaceDim> gridCellCoords;
+      ordinal_type gridCellDivided = gridCell;
+      ordinal_type numGridNodes    = 1;
+      ordinal_type numGridCells    = 1;
+      for (int d=0; d<spaceDim; d++)
+      {
+        gridCellCoords[d] = gridCellDivided % gridCellCounts_[d];
+        gridCellDivided   = gridCellDivided / gridCellCounts_[d];
+        numGridCells *= gridCellCounts_[d];
+        numGridNodes *= (gridCellCounts_[d] + 1);
+      }
+
+      const ordinal_type gridCellNode = gridCellNodeForSubdivisionNode(gridCell, subdivisionOrdinal, node);
+
+      // most subdivision strategies don't add internal node(s), but a couple do.  If so, the gridCellNode
+      // will be equal to the number of vertices in the grid cell.
+      const ordinal_type numInteriorNodes = ((subdivisionStrategy_ == FOUR_TRIANGLES) || (subdivisionStrategy_ == SIX_PYRAMIDS)) ? 1 : 0;
+
+      // the global nodes list the grid-cell-interior nodes first.
+      const ordinal_type gridNodeNumberingOffset = numInteriorNodes * numGridCells;
+
+      const ordinal_type numNodesPerGridCell = 1 << spaceDim;
+      if (gridCellNode >= numNodesPerGridCell)
+      {
+        const ordinal_type interiorGridNodeOffset = gridCellNode - numNodesPerGridCell;
+        return numNodesPerGridCell * gridCell + interiorGridNodeOffset;
+      }
+      else
+      {
+        // use gridCellCoords, plus hypercubeComponentNodeNumber(gridCellNode, d), to set up a Cartesian vertex numbering of the grid as a whole.  Then, add gridNodeNumberingOffset to account for interior nodes.
+        // let x be the fastest-moving index
+        ordinal_type d_index_stride = 1; // number of node indices we move by moving 1 node in dimension d
+        ordinal_type cartesianNodeNumber = 0;
+        for (int d=0; d<spaceDim; d++)
+        {
+          const ordinal_type d_index = gridCellCoords[d] + hypercubeComponentNodeNumber(gridCellNode,d);
+          cartesianNodeNumber += d_index * d_index_stride;
+          d_index_stride *= (gridCellCounts_[d] + 1);
+        }
+        return cartesianNodeNumber + gridNodeNumberingOffset;
+      }
+    }
+    else
+    {
+      // cellToNode() not supported
+      return -1;
+    }
+  }
+
   template<class PointScalar, int spaceDim, typename DeviceType>
   KOKKOS_INLINE_FUNCTION
   DataVariationType CellGeometry<PointScalar,spaceDim,DeviceType>::cellVariationType() const
@@ -1266,6 +1330,65 @@ namespace Intrepid2
         const int gridCellNodeNumber = nodeLookup[subdivisionNodeNumber];
         return gridCellNodeNumber;
       }
+      case SIX_PYRAMIDS:
+      {
+        Kokkos::Array<int,5> nodeLookup; // nodeLookup[4] = 8; // center
+        // we order the subcell divisions as bottom, top, left, right, front, back
+         if (nodeOrdering_ == HYPERCUBE_NODE_ORDER_CLASSIC_SHARDS)
+         {
+          switch (subdivisionOrdinal)
+          {
+            case 0: // bottom (min z)
+              nodeLookup = {0,1,2,3,8};
+              break;
+            case 1: // top (max z)
+              nodeLookup = {4,5,6,7,8};
+              break;
+            case 2: // left (min y)
+              nodeLookup = {0,1,5,4,8};
+              break;
+            case 3: // right (max y)
+              nodeLookup = {3,2,6,7,8};
+              break;
+            case 4: // front (max x)
+              nodeLookup = {1,2,6,5,8};
+              break;
+            case 5: // back (min x)
+              nodeLookup = {0,3,7,4,8};
+              break;
+            default:
+              INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "Invalid subdivisionOrdinal!");
+          }
+        }
+        else // tensor ordering
+        {
+          switch (subdivisionOrdinal)
+          {
+            case 0: // bottom (min z)
+              nodeLookup = {0,1,3,2,8};
+              break;
+            case 1: // top (max z)
+              nodeLookup = {4,5,7,6,8};
+              break;
+            case 2: // left (min y)
+              nodeLookup = {0,1,5,4,8};
+              break;
+            case 3: // right (max y)
+              nodeLookup = {2,3,7,6,8};
+              break;
+            case 4: // front (max x)
+              nodeLookup = {1,3,7,5,8};
+              break;
+            case 5: // back (min x)
+              nodeLookup = {0,2,6,4,8};
+              break;
+            default:
+              INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "Invalid subdivisionOrdinal!");
+        }
+        }
+        const int gridCellNodeNumber = nodeLookup[subdivisionNodeNumber];
+        return gridCellNodeNumber;
+      }
       default:
         INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "Subdivision strategy not yet implemented!");
         // some compilers complain about missing return
@@ -1282,7 +1405,7 @@ namespace Intrepid2
 
     if (subdivisionStrategy_ == FOUR_TRIANGLES)
     {
-      // this is the one case in which the gridCellNode may not actually be a node in the grid cell
+      // this is a case in which the gridCellNode may not actually be a node in the grid cell
       if (gridCellNode == 4) // center vertex
       {
         // d == 0 means quad vertices 0 and 1 suffice;
@@ -1292,6 +1415,18 @@ namespace Intrepid2
         return 0.5 * (gridCellCoordinate(gridCellOrdinal, gridVertex0, d) + gridCellCoordinate(gridCellOrdinal, gridVertex1, d));
       }
     }
+    else if (subdivisionStrategy_ == SIX_PYRAMIDS)
+    {
+      // this is a case in which the gridCellNode may not actually be a node in the grid cell
+      if (gridCellNode == 8) // center vertex
+      {
+        // can compute the center vertex coord in dim d by averaging diagonally opposite vertices'
+        // coords in dimension d.
+        const int gridVertex0 = 0; // 0 is diagonally opposite to 6 or 7, depending on nodeOrdering_
+        const int gridVertex1 = (nodeOrdering_ == HYPERCUBE_NODE_ORDER_CLASSIC_SHARDS) ? 6 : 7;
+        return 0.5 * (gridCellCoordinate(gridCellOrdinal, gridVertex0, d) + gridCellCoordinate(gridCellOrdinal, gridVertex1, d));
+      }
+    }
     return gridCellCoordinate(gridCellOrdinal, gridCellNode, d);
   }
 
@@ -1343,7 +1478,7 @@ namespace Intrepid2
     const int numCellsWorkset = (endCell == -1) ? (numCells_ - startCell) : (endCell - startCell);
 
     using ExecutionSpace = typename DeviceType::execution_space;
-    auto policy = Kokkos::RangePolicy<>(ExecutionSpace(),0,numCellsWorkset);
+    auto policy = Kokkos::RangePolicy<ExecutionSpace>(ExecutionSpace(),0,numCellsWorkset);
     Kokkos::parallel_for("copy orientations", policy, KOKKOS_LAMBDA(const int cellOrdinal)
     {
       orientationsView(cellOrdinal) = orientationsData(cellOrdinal);

packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasisFamily.hpp

@@ -409,7 +409,7 @@ namespace Intrepid2
     {
       case FUNCTION_SPACE_HVOL:  return rcp(new typename BasisFamily::HVOL_PYR (polyOrder, pointType));
 //      case FUNCTION_SPACE_HCURL: return rcp(new typename BasisFamily::HCURL_PYR(polyOrder, pointType));
-//      case FUNCTION_SPACE_HDIV:  return rcp(new typename BasisFamily::HDIV_PYR (polyOrder, pointType));
+      case FUNCTION_SPACE_HDIV:  return rcp(new typename BasisFamily::HDIV_PYR (polyOrder, pointType));
       case FUNCTION_SPACE_HGRAD: return rcp(new typename BasisFamily::HGRAD_PYR(polyOrder, pointType));
       default:
         INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported function space");

packages/intrepid2/src/Discretization/Basis/Intrepid2_HierarchicalBasisFamily.hpp

@@ -54,6 +54,7 @@
 #include "Intrepid2_HierarchicalBasis_HCURL_TET.hpp"
 #include "Intrepid2_HierarchicalBasis_HDIV_TRI.hpp"
 #include "Intrepid2_HierarchicalBasis_HDIV_TET.hpp"
+#include "Intrepid2_HierarchicalBasis_HDIV_PYR.hpp"
 #include "Intrepid2_IntegratedLegendreBasis_HGRAD_LINE.hpp"
 #include "Intrepid2_IntegratedLegendreBasis_HGRAD_TRI.hpp"
 #include "Intrepid2_IntegratedLegendreBasis_HGRAD_TET.hpp"
@@ -115,7 +116,7 @@ namespace Intrepid2 {
     // we will fill these in as we implement them
     using HGRAD = IntegratedLegendreBasis_HGRAD_PYR<DeviceType,OutputScalar,PointScalar,defineVertexFunctions>;
     using HCURL = void;
-    using HDIV  = void;
+    using HDIV  = HierarchicalBasis_HDIV_PYR<DeviceType,OutputScalar,PointScalar>;
     using HVOL  = LegendreBasis_HVOL_PYR<DeviceType,OutputScalar,PointScalar>;
   };