Merge pull request trilinos#7201 from mperego/projections

Intrepid2: Enable projections for Continuous Hierarchical Basis
tcad-charon · May 26, 2020 · 84f4da0 · 84f4da0
2 parents f3c5db3 + f37c63b
commit 84f4da0
Show file tree

Hide file tree

Showing 40 changed files with 5,417 additions and 4,287 deletions.
diff --git a/packages/intrepid2/cmake/Dependencies.cmake b/packages/intrepid2/cmake/Dependencies.cmake
@@ -1,5 +1,5 @@
 SET(LIB_REQUIRED_DEP_PACKAGES TeuchosCore TeuchosNumerics Shards KokkosCore KokkosContainers KokkosAlgorithms)
-SET(LIB_OPTIONAL_DEP_PACKAGES Sacado)
+SET(LIB_OPTIONAL_DEP_PACKAGES Sacado KokkosKernels)
 SET(TEST_REQUIRED_DEP_PACKAGES)
 SET(TEST_OPTIONAL_DEP_PACKAGES Sacado)
 SET(LIB_REQUIRED_DEP_TPLS)

diff --git a/packages/intrepid2/cmake/Intrepid2_config.h.in b/packages/intrepid2/cmake/Intrepid2_config.h.in
@@ -7,3 +7,6 @@
 /* Define if want to build with sacado enabled */
 #cmakedefine HAVE_INTREPID2_SACADO
 
+/* Define if want to build with KokkosKernels enabled */
+#cmakedefine HAVE_INTREPID2_KOKKOSKERNELS
+
diff --git a/packages/intrepid2/src/Cell/Intrepid2_CellTools.hpp b/packages/intrepid2/src/Cell/Intrepid2_CellTools.hpp
@@ -154,6 +154,8 @@ namespace Intrepid2 {
       return r_val;
     }
 
+    typedef Kokkos::DynRankView<double,ExecSpaceType> subcellParamViewType;
+
   private:
 
     /** \brief Generates default HGrad basis based on cell topology 
@@ -250,7 +252,6 @@ namespace Intrepid2 {
     /** \struct Intrepid2::CellTools::SubcellParamData
         \brief Parametrization coefficients of edges and faces of reference cells
      */
-    typedef Kokkos::DynRankView<double,ExecSpaceType> subcellParamViewType;
     struct SubcellParamData {
       subcellParamViewType dummy;
       subcellParamViewType lineEdges;  // edge maps for 2d non-standard cells; shell line and beam
@@ -265,6 +266,24 @@ namespace Intrepid2 {
 
     static bool isSubcellParametrizationSet_;
 
+
+    /** \brief  Sets orientation-preserving parametrizations of reference edges and faces of cell
+        topologies with reference cells. Used to populate Intrepid2::CellTools::SubcellParamData.
+
+        See Intrepid2::CellTools::setSubcellParametrization and Section \ref sec_cell_topology_subcell_map
+        more information about parametrization maps.
+
+        \param  subcellParam           [out]  - array with the coefficients of the parametrization map
+        \param  subcellDim             [in]   - dimension of the subcells being parametrized (1 or 2)
+        \param  parentCell             [in]   - topology of the parent cell owning the subcells.
+    */
+    static void
+    setSubcellParametrization(       subcellParamViewType &subcellParam,
+                               const ordinal_type          subcellDim,
+                               const shards::CellTopology  parentCell );
+
+  public:
+
     /** \brief  Defines orientation-preserving parametrizations of reference edges and faces of cell
         topologies with reference cells.
 
@@ -301,23 +320,6 @@ namespace Intrepid2 {
     */
     static void setSubcellParametrization();
 
-    /** \brief  Sets orientation-preserving parametrizations of reference edges and faces of cell
-        topologies with reference cells. Used to populate Intrepid2::CellTools::SubcellParamData.
-
-        See Intrepid2::CellTools::setSubcellParametrization and Section \ref sec_cell_topology_subcell_map
-        more information about parametrization maps.
-
-        \param  subcellParametrization [out]  - array with the coefficients of the parametrization map
-        \param  subcellDim             [in]   - dimension of the subcells being parametrized (1 or 2)
-        \param  parentCell             [in]   - topology of the parent cell owning the subcells.
-    */
-    static void
-    setSubcellParametrization(       subcellParamViewType &subcellParam,
-                               const ordinal_type          subcellDim,
-                               const shards::CellTopology  parentCell );
-
-  public:
-
    /** \brief  Default constructor.
       */
     CellTools() = default;
@@ -1141,6 +1143,28 @@ namespace Intrepid2 {
                            const shards::CellTopology parentCell );
 
 
+    /**
+     \brief  Overload of mapToReferenceSubcell that runs on device.
+
+        \param  refSubcellPoints  [out] - rank-2 (P,D1) array with images of parameter space points
+        \param  paramPoints       [in]  - rank-2 (P,D2) array with points in 1D or 2D parameter domain
+        \param  subcellMap        [in]  - array with the coefficients of the subcell parametrization map
+        \param  subcellDim        [in]  - dimension of the subcell where points are mapped to
+        \param  subcellOrd        [in]  - subcell ordinal
+        \param  parentCellDim     [in]  - dimension of the parent cell.
+     */
+    template<typename refSubcellPointValueType, class ...refSubcellPointProperties,
+             typename paramPointValueType, class ...paramPointProperties>
+    static void
+    KOKKOS_INLINE_FUNCTION
+    mapToReferenceSubcell(       Kokkos::DynRankView<refSubcellPointValueType,refSubcellPointProperties...> refSubcellPoints,
+                           const Kokkos::DynRankView<paramPointValueType,paramPointProperties...>           paramPoints,
+                           const subcellParamViewType subcellMap,
+                           const ordinal_type subcellDim,
+                           const ordinal_type subcellOrd,
+                           const ordinal_type parentCellDim);
+
+
     //============================================================================================//
     //                                                                                            //
     //                      Physical-to-reference frame mapping and its inverse                   //

diff --git a/packages/intrepid2/src/Cell/Intrepid2_CellToolsDefParametrization.hpp b/packages/intrepid2/src/Cell/Intrepid2_CellToolsDefParametrization.hpp
@@ -67,6 +67,9 @@ namespace Intrepid2 {
   void
   CellTools<SpT>::
   setSubcellParametrization() {
+    if(isSubcellParametrizationSet_)
+      return;
+
     { 
       const auto tet = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >());
       setSubcellParametrization( subcellParamData_.tetFaces,   2, tet );
@@ -119,7 +122,7 @@ namespace Intrepid2 {
       subcellParamData_.wedgeFaces = subcellParamViewType();
     });
 
-    isReferenceNodeDataSet_ = true;
+    isSubcellParametrizationSet_= true;
   }
 
   // template<typename SpT>

diff --git a/packages/intrepid2/src/Cell/Intrepid2_CellToolsDefRefToPhys.hpp b/packages/intrepid2/src/Cell/Intrepid2_CellToolsDefRefToPhys.hpp
@@ -223,20 +223,39 @@ namespace Intrepid2 {
                                   ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints dimension (0) does not match to paramPoints dimension(0).");
 #endif
 
+    if(!isSubcellParametrizationSet_)
+      setSubcellParametrization();
 
-    const ordinal_type cellDim = parentCell.getDimension();
-    const ordinal_type numPts  = paramPoints.extent(0);
+    INTREPID2_TEST_FOR_EXCEPTION( subcellDim != 1 &&
+                                  subcellDim != 2, std::invalid_argument,
+                                  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): method defined only for 1 and 2-subcells");
 
-    // Get the subcell map, i.e., the coefficients of the parametrization function for the subcell
 
-    // can i get this map from devices ?
+    // Get the subcell map, i.e., the coefficients of the parametrization function for the subcell
     subcellParamViewType subcellMap;
     getSubcellParametrization( subcellMap,
                                subcellDim,
                                parentCell );
 
-    // subcell parameterization should be small computation (numPts is small) and it should be decorated with
-    // kokkos inline... let's not do this yet
+    // Apply the parametrization map to every point in parameter domain
+    mapToReferenceSubcell( refSubcellPoints, paramPoints, subcellMap, subcellDim, subcellOrd, parentCell.getDimension());
+  }
+
+
+  template<typename SpT>
+  template<typename refSubcellPointValueType, class ...refSubcellPointProperties,
+           typename paramPointValueType, class ...paramPointProperties>
+  void
+  KOKKOS_INLINE_FUNCTION
+  CellTools<SpT>::
+  mapToReferenceSubcell(       Kokkos::DynRankView<refSubcellPointValueType,refSubcellPointProperties...> refSubcellPoints,
+                         const Kokkos::DynRankView<paramPointValueType,paramPointProperties...>           paramPoints,
+                         const subcellParamViewType subcellMap,
+                         const ordinal_type subcellDim,
+                         const ordinal_type subcellOrd,
+                         const ordinal_type parentCellDim ) {
+
+    const ordinal_type numPts  = paramPoints.extent(0);
 
     // Apply the parametrization map to every point in parameter domain
     switch (subcellDim) {
@@ -246,7 +265,7 @@ namespace Intrepid2 {
         const auto v = paramPoints(pt, 1);
 
         // map_dim(u,v) = c_0(dim) + c_1(dim)*u + c_2(dim)*v because both Quad and Tri ref faces are affine!
-        for (ordinal_type i=0;i<cellDim;++i)
+        for (ordinal_type i=0;i<parentCellDim;++i)
           refSubcellPoints(pt, i) = subcellMap(subcellOrd, i, 0) + ( subcellMap(subcellOrd, i, 1)*u +
                                                                      subcellMap(subcellOrd, i, 2)*v );
       }
@@ -255,16 +274,12 @@ namespace Intrepid2 {
     case 1: {
       for (ordinal_type pt=0;pt<numPts;++pt) {
         const auto u = paramPoints(pt, 0);
-        for (ordinal_type i=0;i<cellDim;++i)
+        for (ordinal_type i=0;i<parentCellDim;++i)
           refSubcellPoints(pt, i) = subcellMap(subcellOrd, i, 0) + ( subcellMap(subcellOrd, i, 1)*u );
       }
       break;
     }
-    default: {
-      INTREPID2_TEST_FOR_EXCEPTION( subcellDim != 1 &&
-                                    subcellDim != 2, std::invalid_argument,
-                                    ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): method defined only for 1 and 2-subcells");
-    }
+    default: {}
     }
   }
 }

diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HCURL_HEX.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HCURL_HEX.hpp
@@ -424,7 +424,7 @@ namespace Intrepid2
     /** \brief True if orientation is required
     */
     virtual bool requireOrientation() const {
-      return true;
+      return (this->getDofCount(1,0) > 0); //if it has edge DOFs, than it needs orientations
     }
   };
 } // end namespace Intrepid2

diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HCURL_QUAD.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HCURL_QUAD.hpp
@@ -257,7 +257,7 @@ namespace Intrepid2
     /** \brief True if orientation is required
     */
     virtual bool requireOrientation() const {
-      return true;
+      return (this->getDofCount(1,0) > 0); //if it has edge DOFs, than it needs orientations
     }
 
   };

diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HDIV_HEX.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HDIV_HEX.hpp
@@ -290,7 +290,7 @@ namespace Intrepid2
         op3 = Intrepid2::OPERATOR_VALUE;
 
         // family 3 goes in the z component; 0 in the x and y components
-        auto outputValuesComponent_xy = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),std::make_pair(0,1));
+        auto outputValuesComponent_xy = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),std::make_pair(0,2));
         auto outputValuesComponent_z = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),2);
 
         // 0 in x and y components
@@ -388,7 +388,7 @@ namespace Intrepid2
     /** \brief True if orientation is required
     */
     virtual bool requireOrientation() const {
-      return true;
+      return (this->getDofCount(2,0) > 0); //if it has side DOFs, than it needs orientations
     }
   };
 } // end namespace Intrepid2

diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HDIV_QUAD.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HDIV_QUAD.hpp
@@ -251,7 +251,7 @@ namespace Intrepid2
     /** \brief True if orientation is required
     */
     virtual bool requireOrientation() const {
-      return true;
+      return (this->getDofCount(1,0) > 0); //if it has side DOFs, than it needs orientations
     }
   };
 } // end namespace Intrepid2

diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_HEX.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_HEX.hpp
@@ -113,7 +113,7 @@ namespace Intrepid2
     /** \brief True if orientation is required
     */
     virtual bool requireOrientation() const override {
-      return (this->getDegree() > 2);
+      return (this->getDofCount(1,0) > 1); //if it has more than 1 DOF per edge, than it needs orientations
     }
 
     using Basis<ExecutionSpace,OutputValueType,PointValueType>::getValues;

diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_QUAD.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_QUAD.hpp
@@ -101,7 +101,7 @@ namespace Intrepid2
     /** \brief True if orientation is required
     */
     virtual bool requireOrientation() const override {
-      return (this->getDegree() > 2);
+      return (this->getDofCount(1,0) > 1); //if it has more than 1 DOF per edge, than it needs orientations
     }
 
     using Basis<ExecutionSpace,OutputValueType,PointValueType>::getValues;

diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HVOL_HEX.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HVOL_HEX.hpp
@@ -70,14 +70,18 @@ namespace Intrepid2
   // TODO: make this a subclass of TensorBasis3 instead, following what we've done for H(curl) and H(div)
   {
   public:
+    using ExecutionSpace  = typename HVOL_LINE::ExecutionSpace;
+    using OutputValueType = typename HVOL_LINE::OutputValueType;
+    using PointValueType  = typename HVOL_LINE::PointValueType;
+
     using OutputViewType = typename HVOL_LINE::OutputViewType;
     using PointViewType  = typename HVOL_LINE::PointViewType ;
     using ScalarViewType = typename HVOL_LINE::ScalarViewType;
 
     using LineBasis = HVOL_LINE;
     using QuadBasis = Intrepid2::Basis_Derived_HVOL_QUAD<HVOL_LINE>;
     using TensorBasis = Basis_TensorBasis<QuadBasis,LineBasis>;
-  public:
+
     /** \brief  Constructor.
         \param [in] polyOrder_x - the polynomial order in the x dimension.
         \param [in] polyOrder_y - the polynomial order in the y dimension.
@@ -96,6 +100,22 @@ namespace Intrepid2
      */
     Basis_Derived_HVOL_HEX(int polyOrder) : Basis_Derived_HVOL_HEX(polyOrder, polyOrder, polyOrder) {}
 
+    /** \brief  Returns basis name
+
+        \return the name of the basis
+    */
+    virtual
+    const char*
+    getName() const {
+      return "Intrepid2_DerivedBasis_HVOL_HEX";
+    }
+
+    /** \brief True if orientation is required
+    */
+    virtual bool requireOrientation() const {
+      return false;
+    }
+
     using TensorBasis::getValues;
 
     /** \brief  multi-component getValues() method (required/called by TensorBasis)

diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HVOL_QUAD.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HVOL_QUAD.hpp
@@ -66,13 +66,19 @@ namespace Intrepid2
   :
   public Basis_TensorBasis<HVOL_LINE, HVOL_LINE>
   {
+    using LineBasis = HVOL_LINE;
+    using TensorBasis = Basis_TensorBasis<LineBasis,LineBasis>;
+
+  public:
+
+   using ExecutionSpace  = typename HVOL_LINE::ExecutionSpace;
+   using OutputValueType = typename HVOL_LINE::OutputValueType;
+   using PointValueType  = typename HVOL_LINE::PointValueType;
+
     using OutputViewType = typename HVOL_LINE::OutputViewType;
     using PointViewType  = typename HVOL_LINE::PointViewType ;
     using ScalarViewType = typename HVOL_LINE::ScalarViewType;
 
-    using LineBasis = HVOL_LINE;
-    using TensorBasis = Basis_TensorBasis<LineBasis,LineBasis>;
-  public:
     /** \brief  Constructor.
         \param [in] polyOrder_x - the polynomial order in the x dimension.
         \param [in] polyOrder_y - the polynomial order in the y dimension.
@@ -90,6 +96,22 @@ namespace Intrepid2
      */
     Basis_Derived_HVOL_QUAD(int polyOrder) : Basis_Derived_HVOL_QUAD(polyOrder,polyOrder) {}
 
+    /** \brief  Returns basis name
+
+        \return the name of the basis
+    */
+    virtual
+    const char*
+    getName() const {
+      return "Intrepid2_DerivedBasis_HVOL_QUAD";
+    }
+
+    /** \brief True if orientation is required
+    */
+    virtual bool requireOrientation() const {
+      return false;
+    }
+
     using TensorBasis::getValues;
 
     /** \brief  multi-component getValues() method (required/called by TensorBasis)

diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_LINE.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_LINE.hpp
@@ -338,6 +338,20 @@ namespace Intrepid2
       }
     }
 
+    /** \brief  Returns basis name
+
+     \return the name of the basis
+    */
+   const char* getName() const override {
+     return "Intrepid2_IntegratedLegendreBasis_HGRAD_LINE";
+   }
+
+   /** \brief True if orientation is required
+   */
+   virtual bool requireOrientation() const override {
+     return false;
+   }
+
     // since the getValues() below only overrides the FEM variant, we specify that
     // we use the base class's getValues(), which implements the FVD variant by throwing an exception.
     // (It's an error to use the FVD variant on this basis.)