diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_Basis.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_Basis.hpp index 0f9b5dff3c07..055c2e7db630 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_Basis.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_Basis.hpp @@ -75,6 +75,11 @@ class Basis; template using BasisPtr = Teuchos::RCP >; +/** \brief Pointer to a Basis whose device type is on the host (Kokkos::HostSpace::device_type), allowing host access to input and output views, and ensuring host execution of basis evaluation. + */ +template +using HostBasisPtr = BasisPtr; + /** \class Intrepid2::Basis \brief An abstract base class that defines interface for concrete basis implementations for Finite Element (FEM) and Finite Volume/Finite Difference (FVD) discrete spaces. @@ -898,11 +903,11 @@ using BasisPtr = Teuchos::RCP >; ">>> ERROR (Basis::getSubCellRefBasis): this method is not supported or should be overridden accordingly by derived classes."); } - /** \brief creates and returns a basis object allocated on host. - - \return pointer to a basis allocated on host. - */ - virtual BasisPtr + /** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this. + + \return Pointer to the new Basis object. + */ + virtual HostBasisPtr getHostBasis() const { INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error, ">>> ERROR (Basis::getHostBasis): this method is not supported or should be overridden accordingly by derived classes."); diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasisFamily.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasisFamily.hpp index b18e1910d1bc..a1a300af7a10 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasisFamily.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasisFamily.hpp @@ -86,7 +86,7 @@ namespace Intrepid2 using OutputValueType = typename LineBasisHGRAD::OutputValueType; using PointValueType = typename LineBasisHGRAD::PointValueType; - using Basis = ::Intrepid2::Basis; + using Basis = typename LineBasisHGRAD::BasisBase; using BasisPtr = Teuchos::RCP; // line bases diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HCURL_HEX.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HCURL_HEX.hpp index b4005fec9784..91699f6ef26a 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HCURL_HEX.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HCURL_HEX.hpp @@ -67,9 +67,7 @@ namespace Intrepid2 { template class Basis_Derived_HCURL_Family1_HEX - : public Basis_TensorBasis3 + : public Basis_TensorBasis3 { public: using OutputViewType = typename HGRAD_LINE::OutputViewType; @@ -79,7 +77,7 @@ namespace Intrepid2 using LineGradBasis = HGRAD_LINE; using LineVolBasis = HVOL_LINE; - using TensorBasis3 = Basis_TensorBasis3; + using TensorBasis3 = Basis_TensorBasis3; public: /** \brief Constructor. \param [in] polyOrder_x - the polynomial order in the x dimension. @@ -224,9 +222,7 @@ namespace Intrepid2 template class Basis_Derived_HCURL_Family2_HEX - : public Basis_TensorBasis3 + : public Basis_TensorBasis3 { public: using OutputViewType = typename HGRAD_LINE::OutputViewType; @@ -236,7 +232,7 @@ namespace Intrepid2 using LineGradBasis = HGRAD_LINE; using LineVolBasis = HVOL_LINE; - using TensorBasis3 = Basis_TensorBasis3; + using TensorBasis3 = Basis_TensorBasis3; public: /** \brief Constructor. \param [in] polyOrder_x - the polynomial order in the x dimension. @@ -388,9 +384,7 @@ namespace Intrepid2 template class Basis_Derived_HCURL_Family3_HEX - : public Basis_TensorBasis3 + : public Basis_TensorBasis3 { using OutputViewType = typename HGRAD_LINE::OutputViewType; using PointViewType = typename HGRAD_LINE::PointViewType ; @@ -399,7 +393,7 @@ namespace Intrepid2 using LineGradBasis = HGRAD_LINE; using LineVolBasis = HVOL_LINE; - using TensorBasis3 = Basis_TensorBasis3; + using TensorBasis3 = Basis_TensorBasis3; public: /** \brief Constructor. \param [in] polyOrder_x - the polynomial order in the x dimension. @@ -540,11 +534,11 @@ namespace Intrepid2 template class Basis_Derived_HCURL_Family1_Family2_HEX - : public Basis_DirectSumBasis + : public Basis_DirectSumBasis { using Family1 = Basis_Derived_HCURL_Family1_HEX; using Family2 = Basis_Derived_HCURL_Family2_HEX; - using DirectSumBasis = Basis_DirectSumBasis ; + using DirectSumBasis = Basis_DirectSumBasis ; public: /** \brief Constructor. \param [in] polyOrder_x - the polynomial order in the x dimension. @@ -560,11 +554,11 @@ namespace Intrepid2 template class Basis_Derived_HCURL_HEX - : public Basis_DirectSumBasis + : public Basis_DirectSumBasis { using Family12 = Basis_Derived_HCURL_Family1_Family2_HEX; using Family3 = Basis_Derived_HCURL_Family3_HEX ; - using DirectSumBasis = Basis_DirectSumBasis ; + using DirectSumBasis = Basis_DirectSumBasis ; std::string name_; ordinal_type order_x_; @@ -577,6 +571,8 @@ namespace Intrepid2 using ExecutionSpace = typename HGRAD_LINE::ExecutionSpace; using OutputValueType = typename HGRAD_LINE::OutputValueType; using PointValueType = typename HGRAD_LINE::PointValueType; + + using BasisBase = typename HGRAD_LINE::BasisBase; /** \brief Constructor. \param [in] polyOrder_x - the polynomial order in the x dimension. @@ -631,7 +627,7 @@ namespace Intrepid2 \param [in] subCellOrd - position of the subCell among of the subCells having the same dimension \return pointer to the subCell basis of dimension subCellDim and position subCellOrd */ - BasisPtr + Teuchos::RCP getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{ using LineBasis = HVOL_LINE; @@ -674,7 +670,18 @@ namespace Intrepid2 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds"); } - + /** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this. + + \return Pointer to the new Basis object. + */ + virtual HostBasisPtr + getHostBasis() const override { + using HostBasis = Basis_Derived_HCURL_HEX; + + auto hostBasis = Teuchos::rcp(new HostBasis(order_x_, order_y_, order_z_, pointType_)); + + return hostBasis; + } }; } // 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 index 29a24435bc8c..f1a833f72719 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HCURL_QUAD.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HCURL_QUAD.hpp @@ -66,7 +66,7 @@ namespace Intrepid2 { template class Basis_Derived_HCURL_Family1_QUAD - : public Basis_TensorBasis + : public Basis_TensorBasis { public: using ExecutionSpace = typename HGRAD_LINE::ExecutionSpace; @@ -77,10 +77,12 @@ namespace Intrepid2 using PointViewType = typename HGRAD_LINE::PointViewType ; using ScalarViewType = typename HGRAD_LINE::ScalarViewType; + using BasisBase = typename HGRAD_LINE::BasisBase; + using LineGradBasis = HGRAD_LINE; using LineHVolBasis = HVOL_LINE; - using TensorBasis = Basis_TensorBasis; + using TensorBasis = Basis_TensorBasis; public: /** \brief Constructor. \param [in] polyOrder_x - the polynomial order in the x dimension. @@ -191,7 +193,7 @@ namespace Intrepid2 template class Basis_Derived_HCURL_Family2_QUAD - : public Basis_TensorBasis + : public Basis_TensorBasis { public: @@ -206,7 +208,9 @@ namespace Intrepid2 using LineGradBasis = HGRAD_LINE; using LineHVolBasis = HVOL_LINE; - using TensorBasis = Basis_TensorBasis; + using BasisBase = typename HGRAD_LINE::BasisBase; + + using TensorBasis = Basis_TensorBasis; /** \brief Constructor. \param [in] polyOrder_x - the polynomial order in the x dimension. @@ -317,11 +321,13 @@ namespace Intrepid2 template class Basis_Derived_HCURL_QUAD - : public Basis_DirectSumBasis + : public Basis_DirectSumBasis { using Family1 = Basis_Derived_HCURL_Family1_QUAD; using Family2 = Basis_Derived_HCURL_Family2_QUAD; - using DirectSumBasis = Basis_DirectSumBasis ; + using DirectSumBasis = Basis_DirectSumBasis ; + + using BasisBase = typename HGRAD_LINE::BasisBase; protected: std::string name_; @@ -388,7 +394,7 @@ namespace Intrepid2 \param [in] subCellOrd - position of the subCell among of the subCells having the same dimension \return pointer to the subCell basis of dimension subCellDim and position subCellOrd */ - BasisPtr + Teuchos::RCP getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{ if(subCellDim == 1) { switch(subCellOrd) { @@ -404,6 +410,18 @@ namespace Intrepid2 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds"); } + /** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this. + + \return Pointer to the new Basis object. + */ + virtual HostBasisPtr + getHostBasis() const override { + using HostBasis = Basis_Derived_HCURL_QUAD; + + auto hostBasis = Teuchos::rcp(new HostBasis(order_x_, order_y_, pointType_)); + + return hostBasis; + } }; } // end namespace Intrepid2 diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HDIV_HEX.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HDIV_HEX.hpp index c47195092fa8..ba52070895a6 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HDIV_HEX.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HDIV_HEX.hpp @@ -71,7 +71,7 @@ namespace Intrepid2 template class Basis_Derived_HDIV_Family1_HEX : - public Basis_TensorBasis3 + public Basis_TensorBasis3 { public: using OutputViewType = typename HGRAD_LINE::OutputViewType; @@ -81,7 +81,7 @@ namespace Intrepid2 using LineGradBasis = HGRAD_LINE; using LineHVolBasis = HVOL_LINE; - using TensorBasis3 = Basis_TensorBasis3; + using TensorBasis3 = Basis_TensorBasis3; public: /** \brief Constructor. \param [in] polyOrder_x - the polynomial order in the x dimension. @@ -205,7 +205,7 @@ namespace Intrepid2 template class Basis_Derived_HDIV_Family2_HEX : - public Basis_TensorBasis3 + public Basis_TensorBasis3 { public: using OutputViewType = typename HGRAD_LINE::OutputViewType; @@ -215,7 +215,7 @@ namespace Intrepid2 using LineGradBasis = HGRAD_LINE; using LineHVolBasis = HVOL_LINE; - using TensorBasis3 = Basis_TensorBasis3; + using TensorBasis3 = Basis_TensorBasis3; public: /** \brief Constructor. \param [in] polyOrder_x - the polynomial order in the x dimension. @@ -344,7 +344,7 @@ namespace Intrepid2 template class Basis_Derived_HDIV_Family3_HEX - : public Basis_TensorBasis3 + : public Basis_TensorBasis3 { public: using OutputViewType = typename HGRAD_LINE::OutputViewType; @@ -354,7 +354,7 @@ namespace Intrepid2 using LineGradBasis = HGRAD_LINE; using LineHVolBasis = HVOL_LINE; - using TensorBasis3 = Basis_TensorBasis3; + using TensorBasis3 = Basis_TensorBasis3; public: /** \brief Constructor. \param [in] polyOrder_x - the polynomial order in the x dimension. @@ -483,11 +483,11 @@ namespace Intrepid2 // which is to say that we go 3,1,2. template class Basis_Derived_HDIV_Family3_Family1_HEX - : public Basis_DirectSumBasis + : public Basis_DirectSumBasis { using Family3 = Basis_Derived_HDIV_Family3_HEX; using Family1 = Basis_Derived_HDIV_Family1_HEX; - using DirectSumBasis = Basis_DirectSumBasis; + using DirectSumBasis = Basis_DirectSumBasis; public: /** \brief Constructor. \param [in] polyOrder_x - the polynomial order in the x dimension. @@ -505,11 +505,11 @@ namespace Intrepid2 template class Basis_Derived_HDIV_HEX - : public Basis_DirectSumBasis + : public Basis_DirectSumBasis { using Family31 = Basis_Derived_HDIV_Family3_Family1_HEX; using Family2 = Basis_Derived_HDIV_Family2_HEX ; - using DirectSumBasis = Basis_DirectSumBasis; + using DirectSumBasis = Basis_DirectSumBasis; std::string name_; ordinal_type order_x_; @@ -521,6 +521,8 @@ namespace Intrepid2 using ExecutionSpace = typename HGRAD_LINE::ExecutionSpace; using OutputValueType = typename HGRAD_LINE::OutputValueType; using PointValueType = typename HGRAD_LINE::PointValueType; + + using BasisBase = typename HGRAD_LINE::BasisBase; /** \brief Constructor. \param [in] polyOrder_x - the polynomial order in the x dimension. @@ -575,8 +577,8 @@ namespace Intrepid2 \param [in] subCellOrd - position of the subCell among of the subCells having the same dimension \return pointer to the subCell basis of dimension subCellDim and position subCellOrd */ - BasisPtr - getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{ + Teuchos::RCP + getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override { using QuadBasis = Basis_Derived_HVOL_QUAD; @@ -599,7 +601,19 @@ namespace Intrepid2 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds"); } - + + /** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this. + + \return Pointer to the new Basis object. + */ + virtual HostBasisPtr + getHostBasis() const override { + using HostBasis = Basis_Derived_HDIV_HEX; + + auto hostBasis = Teuchos::rcp(new HostBasis(order_x_, order_y_, order_z_, pointType_)); + + return hostBasis; + } }; } // 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 index 16c313c6f605..c3e1179916d5 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HDIV_QUAD.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HDIV_QUAD.hpp @@ -66,7 +66,7 @@ namespace Intrepid2 { template class Basis_Derived_HDIV_Family1_QUAD - : public Basis_TensorBasis + : public Basis_TensorBasis { public: using ExecutionSpace = typename HGRAD_LINE::ExecutionSpace; @@ -80,7 +80,9 @@ namespace Intrepid2 using LineGradBasis = HGRAD_LINE; using LineHVolBasis = HVOL_LINE; - using TensorBasis = Basis_TensorBasis; + using BasisBase = typename HGRAD_LINE::BasisBase; + + using TensorBasis = Basis_TensorBasis; public: /** \brief Constructor. \param [in] polyOrder_x - the polynomial order in the x dimension. @@ -204,7 +206,7 @@ namespace Intrepid2 template class Basis_Derived_HDIV_Family2_QUAD - : public Basis_TensorBasis + : public Basis_TensorBasis { using ExecutionSpace = typename HGRAD_LINE::ExecutionSpace; using OutputValueType = typename HGRAD_LINE::OutputValueType; @@ -217,7 +219,9 @@ namespace Intrepid2 using LineGradBasis = HGRAD_LINE; using LineHVolBasis = HVOL_LINE; - using TensorBasis = Basis_TensorBasis; + using BasisBase = typename HGRAD_LINE::BasisBase; + + using TensorBasis = Basis_TensorBasis; public: /** \brief Constructor. \param [in] polyOrder_x - the polynomial order in the x dimension. @@ -331,16 +335,18 @@ namespace Intrepid2 template class Basis_Derived_HDIV_QUAD - : public Basis_DirectSumBasis + : public Basis_DirectSumBasis { using Family1 = Basis_Derived_HDIV_Family1_QUAD; using Family2 = Basis_Derived_HDIV_Family2_QUAD; - using DirectSumBasis = Basis_DirectSumBasis ; + using DirectSumBasis = Basis_DirectSumBasis ; public: using ExecutionSpace = typename HGRAD_LINE::ExecutionSpace; using OutputValueType = typename HGRAD_LINE::OutputValueType; using PointValueType = typename HGRAD_LINE::PointValueType; + + using BasisBase = typename HGRAD_LINE::BasisBase; protected: std::string name_; @@ -401,7 +407,7 @@ namespace Intrepid2 \param [in] subCellOrd - position of the subCell among of the subCells having the same dimension \return pointer to the subCell basis of dimension subCellDim and position subCellOrd */ - BasisPtr + Teuchos::RCP getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{ if(subCellDim == 1) { switch(subCellOrd) { @@ -416,6 +422,19 @@ namespace Intrepid2 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds"); } + + /** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this. + + \return Pointer to the new Basis object. + */ + virtual HostBasisPtr + getHostBasis() const override { + using HostBasis = Basis_Derived_HDIV_QUAD; + + auto hostBasis = Teuchos::rcp(new HostBasis(order_x_, order_y_, pointType_)); + + return hostBasis; + } }; } // 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 index 9633fb7d2aa2..d1027bf3c8a1 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_HEX.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_HEX.hpp @@ -66,7 +66,7 @@ namespace Intrepid2 // TODO: make this a subclass of TensorBasis3 instead, following what we've done for H(curl) and H(div) template class Basis_Derived_HGRAD_HEX - : public Basis_TensorBasis + : public Basis_TensorBasis { public: using ExecutionSpace = typename HGRAD_LINE::ExecutionSpace; @@ -79,7 +79,8 @@ namespace Intrepid2 using LineBasis = HGRAD_LINE; using QuadBasis = Intrepid2::Basis_Derived_HGRAD_QUAD; - using TensorBasis = Basis_TensorBasis; + using BasisBase = typename HGRAD_LINE::BasisBase; + using TensorBasis = Basis_TensorBasis; std::string name_; ordinal_type order_x_; @@ -155,7 +156,7 @@ namespace Intrepid2 } } - using Basis::getValues; + using BasisBase::getValues; /** \brief multi-component getValues() method (required/called by TensorBasis) \param [out] outputValues - the view into which to place the output values @@ -231,7 +232,7 @@ namespace Intrepid2 \param [in] subCellOrd - position of the subCell among of the subCells having the same dimension \return pointer to the subCell basis of dimension subCellDim and position subCellOrd */ - BasisPtr + Teuchos::RCP getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{ if(subCellDim == 1) { switch(subCellOrd) { @@ -270,6 +271,19 @@ namespace Intrepid2 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds"); } + + /** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this. + + \return Pointer to the new Basis object. + */ + virtual HostBasisPtr + getHostBasis() const override { + using HostBasis = Basis_Derived_HGRAD_HEX; + + auto hostBasis = Teuchos::rcp(new HostBasis(order_x_, order_y_, order_z_, pointType_)); + + return hostBasis; + } }; } // end namespace Intrepid2 diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_QUAD.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_QUAD.hpp index 14a4806d2295..35ce945d0db1 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_QUAD.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_QUAD.hpp @@ -57,7 +57,7 @@ namespace Intrepid2 { template class Basis_Derived_HGRAD_QUAD - : public Basis_TensorBasis + : public Basis_TensorBasis { protected: std::string name_; @@ -72,8 +72,9 @@ namespace Intrepid2 using PointViewType = typename HGRAD_LINE::PointViewType ; using ScalarViewType = typename HGRAD_LINE::ScalarViewType; + using BasisBase = typename HGRAD_LINE::BasisBase; using LineBasis = HGRAD_LINE; - using TensorBasis = Basis_TensorBasis; + using TensorBasis = Basis_TensorBasis; /** \brief Constructor. \param [in] polyOrder_x - the polynomial order in the x dimension. @@ -139,7 +140,7 @@ namespace Intrepid2 } } - using Basis::getValues; + using BasisBase::getValues; /** \brief multi-component getValues() method (required/called by TensorBasis) \param [out] outputValues - the view into which to place the output values @@ -216,7 +217,7 @@ namespace Intrepid2 \param [in] subCellOrd - position of the subCell among of the subCells having the same dimension \return pointer to the subCell basis of dimension subCellDim and position subCellOrd */ - BasisPtr + Teuchos::RCP getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{ if(subCellDim == 1) { switch(subCellOrd) { @@ -231,6 +232,19 @@ namespace Intrepid2 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds"); } + + /** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this. + + \return Pointer to the new Basis object. + */ + virtual HostBasisPtr + getHostBasis() const override { + using HostBasis = Basis_Derived_HGRAD_QUAD; + + auto hostBasis = Teuchos::rcp(new HostBasis(order_x_, order_y_, pointType_)); + + return hostBasis; + } }; } // end namespace Intrepid2 diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HVOL_HEX.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HVOL_HEX.hpp index 80bcad07eaf8..1d908985db41 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HVOL_HEX.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HVOL_HEX.hpp @@ -69,11 +69,12 @@ namespace Intrepid2 template class Basis_Derived_HVOL_HEX : - public Basis_TensorBasis + public Basis_TensorBasis // TODO: make this a subclass of TensorBasis3 instead, following what we've done for H(curl) and H(div) { std::string name_; - + ordinal_type polyOrder_x_, polyOrder_y_, polyOrder_z_; + EPointType pointType_; public: using ExecutionSpace = typename HVOL_LINE::ExecutionSpace; using OutputValueType = typename HVOL_LINE::OutputValueType; @@ -83,9 +84,11 @@ namespace Intrepid2 using PointViewType = typename HVOL_LINE::PointViewType ; using ScalarViewType = typename HVOL_LINE::ScalarViewType; + using BasisBase = typename HVOL_LINE::BasisBase; + using LineBasis = HVOL_LINE; using QuadBasis = Intrepid2::Basis_Derived_HVOL_QUAD; - using TensorBasis = Basis_TensorBasis; + using TensorBasis = Basis_TensorBasis; /** \brief Constructor. \param [in] polyOrder_x - the polynomial order in the x dimension. @@ -96,7 +99,11 @@ namespace Intrepid2 Basis_Derived_HVOL_HEX(int polyOrder_x, int polyOrder_y, int polyOrder_z, const EPointType pointType=POINTTYPE_DEFAULT) : TensorBasis(Teuchos::rcp(new QuadBasis(polyOrder_x,polyOrder_y,pointType)), - Teuchos::rcp(new LineBasis(polyOrder_z,pointType))) + Teuchos::rcp(new LineBasis(polyOrder_z,pointType))), + polyOrder_x_(polyOrder_x), + polyOrder_y_(polyOrder_y), + polyOrder_z_(polyOrder_z), + pointType_(pointType) { this->functionSpace_ = FUNCTION_SPACE_HVOL; @@ -169,6 +176,16 @@ namespace Intrepid2 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"operator not yet supported"); } } + + /** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this. + + \return Pointer to the new Basis object. + */ + virtual BasisPtr + getHostBasis() const override { + using HostBasisType = Basis_Derived_HVOL_HEX; + return Teuchos::rcp( new HostBasisType(polyOrder_x_, polyOrder_y_, polyOrder_z_, pointType_) ); + } }; } // end namespace Intrepid2 diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HVOL_QUAD.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HVOL_QUAD.hpp index dbb741f8e144..061ec8ecd18c 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HVOL_QUAD.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HVOL_QUAD.hpp @@ -64,22 +64,25 @@ namespace Intrepid2 template class Basis_Derived_HVOL_QUAD : - public Basis_TensorBasis + public Basis_TensorBasis { protected: std::string name_; using LineBasis = HVOL_LINE; - using TensorBasis = Basis_TensorBasis; + using TensorBasis = Basis_TensorBasis; + + ordinal_type polyOrder_x_, polyOrder_y_; + EPointType pointType_; public: - - using ExecutionSpace = typename HVOL_LINE::ExecutionSpace; - using OutputValueType = typename HVOL_LINE::OutputValueType; - using PointValueType = typename HVOL_LINE::PointValueType; + 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 BasisBase = typename HVOL_LINE::BasisBase; /** \brief Constructor. \param [in] polyOrder_x - the polynomial order in the x dimension. @@ -89,7 +92,10 @@ namespace Intrepid2 Basis_Derived_HVOL_QUAD(int polyOrder_x, int polyOrder_y, const EPointType pointType=POINTTYPE_DEFAULT) : TensorBasis(Teuchos::rcp( new LineBasis(polyOrder_x, pointType)), - Teuchos::rcp( new LineBasis(polyOrder_y, pointType))) + Teuchos::rcp( new LineBasis(polyOrder_y, pointType))), + polyOrder_x_(polyOrder_x), + polyOrder_y_(polyOrder_y), + pointType_(pointType) { this->functionSpace_ = FUNCTION_SPACE_HVOL; @@ -162,6 +168,16 @@ namespace Intrepid2 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"operator not yet supported"); } } + + /** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this. + + \return Pointer to the new Basis object. + */ + virtual BasisPtr + getHostBasis() const override { + using HostBasisType = Basis_Derived_HVOL_QUAD; + return Teuchos::rcp( new HostBasisType(polyOrder_x_, polyOrder_y_, pointType_) ); + } }; } // end namespace Intrepid2 diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DirectSumBasis.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DirectSumBasis.hpp index 445295a56ae7..73d0683ec916 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DirectSumBasis.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DirectSumBasis.hpp @@ -64,24 +64,22 @@ namespace Intrepid2 The two bases must agree in their BasisType (the return value of getBasisType()). */ - template - class Basis_DirectSumBasis : public Basis + template + class Basis_DirectSumBasis : public BasisBaseClass { public: - using BasisSuper = ::Intrepid2::Basis; - using BasisPtr = Teuchos::RCP; + using BasisBase = BasisBaseClass; + using BasisPtr = Teuchos::RCP; - using ExecutionSpace = typename BasisSuper::ExecutionSpace; - using OutputValueType = typename BasisSuper::OutputValueType; - using PointValueType = typename BasisSuper::PointValueType; + using ExecutionSpace = typename BasisBase::ExecutionSpace; + using OutputValueType = typename BasisBase::OutputValueType; + using PointValueType = typename BasisBase::PointValueType; - using OrdinalTypeArray1DHost = typename BasisSuper::OrdinalTypeArray1DHost; - using OrdinalTypeArray2DHost = typename BasisSuper::OrdinalTypeArray2DHost; - using OutputViewType = typename BasisSuper::OutputViewType; - using PointViewType = typename BasisSuper::PointViewType; - using ScalarViewType = typename BasisSuper::ScalarViewType; + using OrdinalTypeArray1DHost = typename BasisBase::OrdinalTypeArray1DHost; + using OrdinalTypeArray2DHost = typename BasisBase::OrdinalTypeArray2DHost; + using OutputViewType = typename BasisBase::OutputViewType; + using PointViewType = typename BasisBase::PointViewType; + using ScalarViewType = typename BasisBase::ScalarViewType; protected: BasisPtr basis1_; BasisPtr basis2_; @@ -208,10 +206,10 @@ namespace Intrepid2 The default implementation employs a trivial tensor-product structure, for compatibility across all bases. Subclasses that have tensor-product structure should override. Note that only the basic exact-sequence operators are supported at the moment: VALUE, GRAD, DIV, CURL. */ - virtual BasisValues allocateBasisValues( TensorPoints points, const EOperator operatorType = OPERATOR_VALUE) const override + virtual BasisValues allocateBasisValues( TensorPoints points, const EOperator operatorType = OPERATOR_VALUE) const override { - BasisValues basisValues1 = basis1_->allocateBasisValues(points, operatorType); - BasisValues basisValues2 = basis2_->allocateBasisValues(points, operatorType); + BasisValues basisValues1 = basis1_->allocateBasisValues(points, operatorType); + BasisValues basisValues2 = basis2_->allocateBasisValues(points, operatorType); const int numScalarFamilies1 = basisValues1.numTensorDataFamilies(); if (numScalarFamilies1 > 0) @@ -219,7 +217,7 @@ namespace Intrepid2 // then both basis1 and basis2 should be scalar-valued; check that for basis2: const int numScalarFamilies2 = basisValues2.numTensorDataFamilies(); INTREPID2_TEST_FOR_EXCEPTION(basisValues2.numTensorDataFamilies() <=0, std::invalid_argument, "When basis1 has scalar value, basis2 must also"); - std::vector< TensorData > scalarFamilies(numScalarFamilies1 + numScalarFamilies2); + std::vector< TensorData > scalarFamilies(numScalarFamilies1 + numScalarFamilies2); for (int i=0; i(scalarFamilies); + return BasisValues(scalarFamilies); } else { @@ -245,7 +243,7 @@ namespace Intrepid2 const int numFamilies2 = vectorData2.numFamilies(); const int numFamilies = numFamilies1 + numFamilies2; - std::vector< std::vector > > vectorComponents(numFamilies, std::vector >(numComponents)); + std::vector< std::vector > > vectorComponents(numFamilies, std::vector >(numComponents)); for (int i=0; i vectorData(vectorComponents); - return BasisValues(vectorData); + VectorData vectorData(vectorComponents); + return BasisValues(vectorData); } } @@ -323,7 +321,7 @@ namespace Intrepid2 // since the getValues() below only overrides the FEM variants, 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.) - using BasisSuper::getValues; + using BasisBase::getValues; /** \brief Evaluation of a FEM basis on a reference cell, using point and output value containers that allow preservation of tensor-product structure. @@ -338,8 +336,8 @@ namespace Intrepid2 */ virtual void - getValues( BasisValues outputValues, - const TensorPoints inputPoints, + getValues( BasisValues outputValues, + const TensorPoints inputPoints, const EOperator operatorType = OPERATOR_VALUE ) const override { const int fieldStartOrdinal1 = 0; diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_LINE_C1_FEM.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_LINE_C1_FEM.hpp index 10c3b2768403..127115863cf5 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_LINE_C1_FEM.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_LINE_C1_FEM.hpp @@ -155,19 +155,21 @@ namespace Intrepid2 { class Basis_HGRAD_LINE_C1_FEM : public Basis { public: - using OrdinalTypeArray1DHost = typename Basis::OrdinalTypeArray1DHost; - using OrdinalTypeArray2DHost = typename Basis::OrdinalTypeArray2DHost; - using OrdinalTypeArray3DHost = typename Basis::OrdinalTypeArray3DHost; - + using BasisBase = Basis; + + using OrdinalTypeArray1DHost = typename BasisBase::OrdinalTypeArray1DHost; + using OrdinalTypeArray2DHost = typename BasisBase::OrdinalTypeArray2DHost; + using OrdinalTypeArray3DHost = typename BasisBase::OrdinalTypeArray3DHost; + + using OutputViewType = typename BasisBase::OutputViewType; + using PointViewType = typename BasisBase::PointViewType ; + using ScalarViewType = typename BasisBase::ScalarViewType; + /** \brief Constructor. */ Basis_HGRAD_LINE_C1_FEM(); - using OutputViewType = typename Basis::OutputViewType; - using PointViewType = typename Basis::PointViewType; - using ScalarViewType = typename Basis::ScalarViewType; - - using Basis::getValues; + using BasisBase::getValues; virtual void diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_LINE_Cn_FEM.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_LINE_Cn_FEM.hpp index d477c9f65f70..560b8d9f5bbe 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_LINE_Cn_FEM.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_LINE_Cn_FEM.hpp @@ -175,27 +175,31 @@ namespace Intrepid2 { class Basis_HGRAD_LINE_Cn_FEM : public Basis { public: - using OrdinalTypeArray1DHost = typename Basis::OrdinalTypeArray1DHost; - using OrdinalTypeArray2DHost = typename Basis::OrdinalTypeArray2DHost; - using OrdinalTypeArray3DHost = typename Basis::OrdinalTypeArray3DHost; - - using OutputViewType = typename Basis::OutputViewType; - using PointViewType = typename Basis::PointViewType; - using ScalarViewType = typename Basis::ScalarViewType; + using BasisBase = Basis; + + using HostBasis = Basis_HGRAD_LINE_Cn_FEM; + + using OrdinalTypeArray1DHost = typename BasisBase::OrdinalTypeArray1DHost; + using OrdinalTypeArray2DHost = typename BasisBase::OrdinalTypeArray2DHost; + using OrdinalTypeArray3DHost = typename BasisBase::OrdinalTypeArray3DHost; + + using OutputViewType = typename BasisBase::OutputViewType; + using PointViewType = typename BasisBase::PointViewType ; + using ScalarViewType = typename BasisBase::ScalarViewType; private: /** \brief inverse of Generalized Vandermonde matrix, whose columns store the expansion coefficients of the nodal basis in terms of phis_ */ Kokkos::DynRankView vinv_; - + EPointType pointType_; public: /** \brief Constructor. */ Basis_HGRAD_LINE_Cn_FEM(const ordinal_type order, const EPointType pointType = POINTTYPE_EQUISPACED); - using Basis::getValues; + using BasisBase::getValues; virtual void @@ -270,6 +274,16 @@ namespace Intrepid2 { return getPnCardinality<1>(this->basisDegree_); } + /** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this. + + \return Pointer to the new Basis object. + */ + virtual HostBasisPtr + getHostBasis() const override { + auto hostBasis = Teuchos::rcp(new HostBasis(this->basisDegree_, pointType_)); + + return hostBasis; + } }; }// namespace Intrepid2 diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_LINE_Cn_FEMDef.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_LINE_Cn_FEMDef.hpp index 7f47943ae98a..23223b796084 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_LINE_Cn_FEMDef.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_LINE_Cn_FEMDef.hpp @@ -196,6 +196,7 @@ namespace Intrepid2 { Basis_HGRAD_LINE_Cn_FEM:: Basis_HGRAD_LINE_Cn_FEM( const ordinal_type order, const EPointType pointType ) { + this->pointType_ = pointType; this->basisCardinality_ = order+1; this->basisDegree_ = order; this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData >() ); diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HVOL_LINE_Cn_FEM.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HVOL_LINE_Cn_FEM.hpp index 7672df7c89a5..d64c4c82d11b 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HVOL_LINE_Cn_FEM.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HVOL_LINE_Cn_FEM.hpp @@ -177,20 +177,23 @@ namespace Intrepid2 { class Basis_HVOL_LINE_Cn_FEM : public Basis { public: - using OrdinalTypeArray1DHost = typename Basis::OrdinalTypeArray1DHost; - using OrdinalTypeArray2DHost = typename Basis::OrdinalTypeArray2DHost; - using OrdinalTypeArray3DHost = typename Basis::OrdinalTypeArray3DHost; + using BasisBase = Basis; + using HostBasis = Basis_HVOL_LINE_Cn_FEM; + + using OrdinalTypeArray1DHost = typename BasisBase::OrdinalTypeArray1DHost; + using OrdinalTypeArray2DHost = typename BasisBase::OrdinalTypeArray2DHost; + using OrdinalTypeArray3DHost = typename BasisBase::OrdinalTypeArray3DHost; + + using OutputViewType = typename BasisBase::OutputViewType; + using PointViewType = typename BasisBase::PointViewType ; + using ScalarViewType = typename BasisBase::ScalarViewType; /** \brief Constructor. */ Basis_HVOL_LINE_Cn_FEM(const ordinal_type order, - const EPointType pointType = POINTTYPE_EQUISPACED); - - using OutputViewType = typename Basis::OutputViewType; - using PointViewType = typename Basis::PointViewType; - using ScalarViewType = typename Basis::ScalarViewType; + const EPointType pointType = POINTTYPE_EQUISPACED); - using Basis::getValues; + using BasisBase::getValues; virtual void diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HierarchicalBasisFamily.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HierarchicalBasisFamily.hpp index a113df79f031..2b08a170d494 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HierarchicalBasisFamily.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HierarchicalBasisFamily.hpp @@ -69,7 +69,7 @@ namespace Intrepid2 { // the following defines a family of hierarchical basis functions that matches the unpermuted ESEAS basis functions // each basis member is associated with appropriate subcell topologies, making this suitable for continuous Galerkin finite elements. - template @@ -77,24 +77,24 @@ namespace Intrepid2 { { public: // we will fill these in as we implement them - using HGRAD = IntegratedLegendreBasis_HGRAD_TRI; - using HCURL = dummyBasis; - using HDIV = dummyBasis; - using HVOL = dummyBasis; + using HGRAD = IntegratedLegendreBasis_HGRAD_TRI; + using HCURL = dummyBasis; + using HDIV = dummyBasis; + using HVOL = dummyBasis; }; - template + template class HierarchicalTetrahedronBasisFamily { public: // we will fill these in as we implement them - using HGRAD = IntegratedLegendreBasis_HGRAD_TET; - using HCURL = dummyBasis; - using HDIV = dummyBasis; - using HVOL = dummyBasis; + using HGRAD = IntegratedLegendreBasis_HGRAD_TET; + using HCURL = dummyBasis; + using HDIV = dummyBasis; + using HVOL = dummyBasis; }; /** \class Intrepid2::HierarchicalBasisFamily @@ -117,13 +117,13 @@ namespace Intrepid2 { We have offline tests that verify agreement between our implementation and ESEAS. We hope to add these to the Trilinos continuous integration tests in the future. */ - template - using HierarchicalBasisFamily = DerivedBasisFamily< IntegratedLegendreBasis_HGRAD_LINE, - LegendreBasis_HVOL_LINE, - HierarchicalTriangleBasisFamily, - HierarchicalTetrahedronBasisFamily + using HierarchicalBasisFamily = DerivedBasisFamily< IntegratedLegendreBasis_HGRAD_LINE, + LegendreBasis_HVOL_LINE, + HierarchicalTriangleBasisFamily, + HierarchicalTetrahedronBasisFamily >; /** \class Intrepid2::HierarchicalBasisFamily @@ -134,13 +134,13 @@ namespace Intrepid2 { The suitability of this family for DG contexts is primarily due to the fact that the H(grad) basis has a constant member. Note also that in this family, all members are associated with the cell interior; there are no basis functions associated with subcell topologies. */ - template - using DGHierarchicalBasisFamily = DerivedBasisFamily< IntegratedLegendreBasis_HGRAD_LINE, - LegendreBasis_HVOL_LINE, - HierarchicalTriangleBasisFamily, - HierarchicalTetrahedronBasisFamily + using DGHierarchicalBasisFamily = DerivedBasisFamily< IntegratedLegendreBasis_HGRAD_LINE, + LegendreBasis_HVOL_LINE, + HierarchicalTriangleBasisFamily, + HierarchicalTetrahedronBasisFamily >; } diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_LINE.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_LINE.hpp index e44646020cc0..70e950dead58 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_LINE.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_LINE.hpp @@ -64,10 +64,11 @@ namespace Intrepid2 This functor is not intended for use outside of IntegratedLegendreBasis_HGRAD_LINE. */ - template struct Hierarchical_HGRAD_LINE_Functor { + using ExecutionSpace = typename DeviceType::execution_space; using ScratchSpace = typename ExecutionSpace::scratch_memory_space; using OutputScratchView = Kokkos::View>; using PointScratchView = Kokkos::View>; @@ -227,24 +228,28 @@ namespace Intrepid2 is true, then the first basis function will instead be 1.0-x, and the basis will be suitable for continuous discretizations. */ - template // if useMinusOneToOneReferenceElement is true, basis is define on [-1,1]. Otherwise, [0,1]. class IntegratedLegendreBasis_HGRAD_LINE - : public Basis + : public Basis { public: - using OrdinalTypeArray1DHost = typename Basis::OrdinalTypeArray1DHost; - using OrdinalTypeArray2DHost = typename Basis::OrdinalTypeArray2DHost; + using BasisBase = Basis; + using HostBasis = IntegratedLegendreBasis_HGRAD_LINE; - using OutputViewType = typename Basis::OutputViewType; - using PointViewType = typename Basis::PointViewType ; - using ScalarViewType = typename Basis::ScalarViewType; + using OrdinalTypeArray1DHost = typename BasisBase::OrdinalTypeArray1DHost; + using OrdinalTypeArray2DHost = typename BasisBase::OrdinalTypeArray2DHost; + + using OutputViewType = typename BasisBase::OutputViewType; + using PointViewType = typename BasisBase::PointViewType ; + using ScalarViewType = typename BasisBase::ScalarViewType; protected: int polyOrder_; // the maximum order of the polynomial bool defineVertexFunctions_; // if true, first and second basis functions are x and 1-x. Otherwise, they are 1 and x. + EPointType pointType_; public: /** \brief Constructor. \param [in] polyOrder - the polynomial order of the basis. @@ -258,7 +263,8 @@ namespace Intrepid2 */ IntegratedLegendreBasis_HGRAD_LINE(int polyOrder, EPointType pointType=POINTTYPE_DEFAULT) : - polyOrder_(polyOrder) + polyOrder_(polyOrder), + pointType_(pointType) { INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,"PointType not supported"); @@ -357,7 +363,7 @@ namespace Intrepid2 // 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.) - using Basis::getValues; + using BasisBase::getValues; /** \brief Evaluation of a FEM basis on a reference cell. @@ -382,7 +388,7 @@ namespace Intrepid2 { auto numPoints = inputPoints.extent_int(0); - using FunctorType = Hierarchical_HGRAD_LINE_Functor; + using FunctorType = Hierarchical_HGRAD_LINE_Functor; FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_, defineVertexFunctions); @@ -391,8 +397,21 @@ namespace Intrepid2 const int vectorSize = std::max(outputVectorSize,pointVectorSize); const int teamSize = 1; // because of the way the basis functions are computed, we don't have a second level of parallelism... + using ExecutionSpace = typename BasisBase::ExecutionSpace; + auto policy = Kokkos::TeamPolicy(numPoints,teamSize,vectorSize); - Kokkos::parallel_for( policy , functor, "Hierarchical_HGRAD_LINE_Functor"); + Kokkos::parallel_for( policy, functor, "Hierarchical_HGRAD_LINE_Functor"); + } + + /** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this. + + \return Pointer to the new Basis object. + */ + virtual BasisPtr + getHostBasis() const override { + using HostDeviceType = typename Kokkos::HostSpace::device_type; + using HostBasisType = IntegratedLegendreBasis_HGRAD_LINE; + return Teuchos::rcp( new HostBasisType(polyOrder_, pointType_) ); } }; } // end namespace Intrepid2 diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_TET.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_TET.hpp index 80081d4fa0fa..d3218729afc5 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_TET.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_TET.hpp @@ -64,10 +64,11 @@ namespace Intrepid2 This functor is not intended for use outside of IntegratedLegendreBasis_HGRAD_TET. */ - template struct Hierarchical_HGRAD_TET_Functor { + using ExecutionSpace = typename DeviceType::execution_space; using ScratchSpace = typename ExecutionSpace::scratch_memory_space; using OutputScratchView = Kokkos::View>; using OutputScratchView2D = Kokkos::View>; @@ -602,22 +603,25 @@ namespace Intrepid2 is true, then the first basis function will instead be 1.0-x-y, and the basis will be suitable for continuous discretizations. */ - template // if defineVertexFunctions is true, first four basis functions are 1-x-y-z, x, y, and z. Otherwise, they are 1, x, y, and z. class IntegratedLegendreBasis_HGRAD_TET - : public Basis + : public Basis { public: - using OrdinalTypeArray1DHost = typename Basis::OrdinalTypeArray1DHost; - using OrdinalTypeArray2DHost = typename Basis::OrdinalTypeArray2DHost; + using BasisBase = Basis; - using OutputViewType = typename Basis::OutputViewType; - using PointViewType = typename Basis::PointViewType; - using ScalarViewType = typename Basis::ScalarViewType; + using OrdinalTypeArray1DHost = typename BasisBase::OrdinalTypeArray1DHost; + using OrdinalTypeArray2DHost = typename BasisBase::OrdinalTypeArray2DHost; + + using OutputViewType = typename BasisBase::OutputViewType; + using PointViewType = typename BasisBase::PointViewType ; + using ScalarViewType = typename BasisBase::ScalarViewType; protected: int polyOrder_; // the maximum order of the polynomial + EPointType pointType_; public: /** \brief Constructor. \param [in] polyOrder - the polynomial order of the basis. @@ -631,7 +635,8 @@ namespace Intrepid2 */ IntegratedLegendreBasis_HGRAD_TET(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT) : - polyOrder_(polyOrder) + polyOrder_(polyOrder), + pointType_(pointType) { INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,"PointType not supported"); this->basisCardinality_ = ((polyOrder+1) * (polyOrder+2) * (polyOrder+3)) / 6; @@ -819,7 +824,7 @@ namespace Intrepid2 // 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.) - using Basis::getValues; + using BasisBase::getValues; /** \brief Evaluation of a FEM basis on a reference cell. @@ -844,7 +849,7 @@ namespace Intrepid2 { auto numPoints = inputPoints.extent_int(0); - using FunctorType = Hierarchical_HGRAD_TET_Functor; + using FunctorType = Hierarchical_HGRAD_TET_Functor; FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_, defineVertexFunctions); @@ -853,6 +858,8 @@ namespace Intrepid2 const int vectorSize = std::max(outputVectorSize,pointVectorSize); const int teamSize = 1; // because of the way the basis functions are computed, we don't have a second level of parallelism... + using ExecutionSpace = typename BasisBase::ExecutionSpace; + auto policy = Kokkos::TeamPolicy(numPoints,teamSize,vectorSize); Kokkos::parallel_for( policy , functor, "Hierarchical_HGRAD_TET_Functor"); } @@ -865,20 +872,30 @@ namespace Intrepid2 \param [in] subCellOrd - position of the subCell among of the subCells having the same dimension \return pointer to the subCell basis of dimension subCellDim and position subCellOrd */ - BasisPtr + BasisPtr getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{ if(subCellDim == 1) { return Teuchos::rcp(new - IntegratedLegendreBasis_HGRAD_LINE + IntegratedLegendreBasis_HGRAD_LINE (this->basisDegree_)); } else if(subCellDim == 2) { return Teuchos::rcp(new - IntegratedLegendreBasis_HGRAD_TRI + IntegratedLegendreBasis_HGRAD_TRI (this->basisDegree_)); } INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds"); } + /** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this. + + \return Pointer to the new Basis object. + */ + virtual BasisPtr + getHostBasis() const override { + using HostDeviceType = typename Kokkos::HostSpace::device_type; + using HostBasisType = IntegratedLegendreBasis_HGRAD_TET; + return Teuchos::rcp( new HostBasisType(polyOrder_, pointType_) ); + } }; } // end namespace Intrepid2 diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_TRI.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_TRI.hpp index 3d2c2e0bc310..e007c45b0fff 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_TRI.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_TRI.hpp @@ -64,10 +64,11 @@ namespace Intrepid2 This functor is not intended for use outside of IntegratedLegendreBasis_HGRAD_TRI. */ - template struct Hierarchical_HGRAD_TRI_Functor { + using ExecutionSpace = typename DeviceType::execution_space; using ScratchSpace = typename ExecutionSpace::scratch_memory_space; using OutputScratchView = Kokkos::View>; using PointScratchView = Kokkos::View>; @@ -357,22 +358,25 @@ namespace Intrepid2 is true, then the first basis function will instead be 1.0-x-y, and the basis will be suitable for continuous discretizations. */ - template // if defineVertexFunctions is true, first three basis functions are 1-x-y, x, and y. Otherwise, they are 1, x, and y. class IntegratedLegendreBasis_HGRAD_TRI - : public Basis + : public Basis { public: - using OrdinalTypeArray1DHost = typename Basis::OrdinalTypeArray1DHost; - using OrdinalTypeArray2DHost = typename Basis::OrdinalTypeArray2DHost; + using BasisBase = Basis; - using OutputViewType = typename Basis::OutputViewType; - using PointViewType = typename Basis::PointViewType; - using ScalarViewType = typename Basis::ScalarViewType; + using OrdinalTypeArray1DHost = typename BasisBase::OrdinalTypeArray1DHost; + using OrdinalTypeArray2DHost = typename BasisBase::OrdinalTypeArray2DHost; + + using OutputViewType = typename BasisBase::OutputViewType; + using PointViewType = typename BasisBase::PointViewType ; + using ScalarViewType = typename BasisBase::ScalarViewType; protected: int polyOrder_; // the maximum order of the polynomial + EPointType pointType_; public: /** \brief Constructor. \param [in] polyOrder - the polynomial order of the basis. @@ -386,7 +390,8 @@ namespace Intrepid2 */ IntegratedLegendreBasis_HGRAD_TRI(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT) : - polyOrder_(polyOrder) + polyOrder_(polyOrder), + pointType_(pointType) { INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,"PointType not supported"); @@ -532,7 +537,7 @@ namespace Intrepid2 // 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.) - using Basis::getValues; + using BasisBase::getValues; /** \brief Evaluation of a FEM basis on a reference cell. @@ -557,7 +562,7 @@ namespace Intrepid2 { auto numPoints = inputPoints.extent_int(0); - using FunctorType = Hierarchical_HGRAD_TRI_Functor; + using FunctorType = Hierarchical_HGRAD_TRI_Functor; FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_, defineVertexFunctions); @@ -566,6 +571,8 @@ namespace Intrepid2 const int vectorSize = std::max(outputVectorSize,pointVectorSize); const int teamSize = 1; // because of the way the basis functions are computed, we don't have a second level of parallelism... + using ExecutionSpace = typename BasisBase::ExecutionSpace; + auto policy = Kokkos::TeamPolicy(numPoints,teamSize,vectorSize); Kokkos::parallel_for( policy , functor, "Hierarchical_HGRAD_TRI_Functor"); } @@ -578,16 +585,26 @@ namespace Intrepid2 \param [in] subCellOrd - position of the subCell among of the subCells having the same dimension \return pointer to the subCell basis of dimension subCellDim and position subCellOrd */ - BasisPtr + BasisPtr getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{ if(subCellDim == 1) { return Teuchos::rcp(new - IntegratedLegendreBasis_HGRAD_LINE + IntegratedLegendreBasis_HGRAD_LINE (this->basisDegree_)); } INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds"); } + /** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this. + + \return Pointer to the new Basis object. + */ + virtual BasisPtr + getHostBasis() const override { + using HostDeviceType = typename Kokkos::HostSpace::device_type; + using HostBasisType = IntegratedLegendreBasis_HGRAD_TRI; + return Teuchos::rcp( new HostBasisType(polyOrder_, pointType_) ); + } }; } // end namespace Intrepid2 diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_LegendreBasis_HVOL_LINE.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_LegendreBasis_HVOL_LINE.hpp index fcdd21456033..fcaf602a9dd4 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_LegendreBasis_HVOL_LINE.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_LegendreBasis_HVOL_LINE.hpp @@ -69,10 +69,11 @@ namespace Intrepid2 This functor is not intended for use outside of LegendreBasis_HVOL_LINE. */ - template struct Hierarchical_HVOL_LINE_Functor { + using ExecutionSpace = typename DeviceType::execution_space; using ScratchSpace = typename ExecutionSpace::scratch_memory_space; using OutputScratchView = Kokkos::View>; using PointScratchView = Kokkos::View>; @@ -177,21 +178,25 @@ namespace Intrepid2 Computers & Mathematics with Applications, Volume 70, Issue 4, 2015, Pages 353-458, ISSN 0898-1221. https://doi.org/10.1016/j.camwa.2015.04.027. */ - template class LegendreBasis_HVOL_LINE - : public Basis + : public Basis { public: - using OrdinalTypeArray1DHost = typename Basis::OrdinalTypeArray1DHost; - using OrdinalTypeArray2DHost = typename Basis::OrdinalTypeArray2DHost; + using BasisBase = Basis; + using HostBasis = LegendreBasis_HVOL_LINE; - typedef typename Basis::OutputViewType OutputViewType; - typedef typename Basis::PointViewType PointViewType; - typedef typename Basis::ScalarViewType ScalarViewType; + using OrdinalTypeArray1DHost = typename BasisBase::OrdinalTypeArray1DHost; + using OrdinalTypeArray2DHost = typename BasisBase::OrdinalTypeArray2DHost; + + using OutputViewType = typename BasisBase::OutputViewType; + using PointViewType = typename BasisBase::PointViewType ; + using ScalarViewType = typename BasisBase::ScalarViewType; protected: int polyOrder_; // the maximum order of the polynomial + EPointType pointType_; public: /** \brief Constructor. \param [in] polyOrder - the polynomial order of the basis. @@ -203,7 +208,8 @@ namespace Intrepid2 */ LegendreBasis_HVOL_LINE(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT) : - polyOrder_(polyOrder) + polyOrder_(polyOrder), + pointType_(pointType) { INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,"PointType not supported"); this->basisCardinality_ = polyOrder+1; @@ -257,7 +263,7 @@ namespace Intrepid2 // 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.) - using Basis::getValues; + using BasisBase::getValues; /** \brief Returns basis name @@ -292,7 +298,7 @@ namespace Intrepid2 { auto numPoints = inputPoints.extent_int(0); - using FunctorType = Hierarchical_HVOL_LINE_Functor; + using FunctorType = Hierarchical_HVOL_LINE_Functor; FunctorType functor(outputValues, inputPoints, polyOrder_, operatorType); @@ -301,9 +307,22 @@ namespace Intrepid2 const int vectorSize = std::max(outputVectorSize,pointVectorSize); const int teamSize = 1; // because of the way the basis functions are computed, we don't have a second level of parallelism... + using ExecutionSpace = typename BasisBase::ExecutionSpace; + auto policy = Kokkos::TeamPolicy(numPoints,teamSize,vectorSize); Kokkos::parallel_for( policy , functor, "Hierarchical_HVOL_LINE_Functor"); } + + /** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this. + + \return Pointer to the new Basis object. + */ + virtual BasisPtr + getHostBasis() const override { + using HostDeviceType = typename Kokkos::HostSpace::device_type; + using HostBasisType = LegendreBasis_HVOL_LINE; + return Teuchos::rcp( new HostBasisType(polyOrder_, pointType_) ); + } }; } // end namespace Intrepid2 diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_TensorBasis.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_TensorBasis.hpp index 1a8e3c6e8c59..3a48e0c6aca4 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_TensorBasis.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_TensorBasis.hpp @@ -156,7 +156,7 @@ namespace Intrepid2 return getDkEnumeration(entries); } -template +template class Basis_TensorBasis; /** \struct Intrepid2::OperatorTensorDecomposition @@ -274,14 +274,15 @@ struct OperatorTensorDecomposition } //! takes as argument bases that are components in this decomposition, and decomposes them further if they are tensor bases. Returns a fully expanded decomposition. - template - OperatorTensorDecomposition expandedDecomposition(std::vector< Teuchos::RCP > > &bases) + template + OperatorTensorDecomposition expandedDecomposition(std::vector< Teuchos::RCP > > &bases) { const ordinal_type basesSize = bases.size(); INTREPID2_TEST_FOR_EXCEPTION(basesSize != numBasisComponents_, std::invalid_argument, "The number of bases provided must match the number of basis components in this decomposition"); ordinal_type numExpandedBasisComponents = 0; - using TensorBasis = Basis_TensorBasis; + using BasisBase = Basis; + using TensorBasis = Basis_TensorBasis; std::vector basesAsTensorBasis(numBasisComponents_); for (ordinal_type basisComponentOrdinal=0; basisComponentOrdinal + template class Basis_TensorBasis : - public Basis + public BasisBaseClass { public: - using BasisSuper = ::Intrepid2::Basis; - using BasisPtr = Teuchos::RCP; + using BasisBase = BasisBaseClass; + using BasisPtr = Teuchos::RCP; protected: BasisPtr basis1_; @@ -572,24 +571,26 @@ struct OperatorTensorDecomposition std::string name_; // name of the basis public: - using ExecutionSpace = typename BasisSuper::ExecutionSpace; - using OutputValueType = typename BasisSuper::OutputValueType; - using PointValueType = typename BasisSuper::PointValueType; + using ExecutionSpace = typename BasisBase::ExecutionSpace; + using OutputValueType = typename BasisBase::OutputValueType; + using PointValueType = typename BasisBase::PointValueType; - using OrdinalTypeArray1DHost = typename BasisSuper::OrdinalTypeArray1DHost; - using OrdinalTypeArray2DHost = typename BasisSuper::OrdinalTypeArray2DHost; - using OutputViewType = typename BasisSuper::OutputViewType; - using PointViewType = typename BasisSuper::PointViewType; - using TensorBasis = Basis_TensorBasis; + using OrdinalTypeArray1DHost = typename BasisBase::OrdinalTypeArray1DHost; + using OrdinalTypeArray2DHost = typename BasisBase::OrdinalTypeArray2DHost; + using OutputViewType = typename BasisBase::OutputViewType; + using PointViewType = typename BasisBase::PointViewType; + using TensorBasis = Basis_TensorBasis; public: /** \brief Constructor. \param [in] basis1 - the first component basis \param [in] basis2 - the second component basis */ - Basis_TensorBasis(BasisPtr basis1, BasisPtr basis2) + Basis_TensorBasis(BasisPtr basis1, BasisPtr basis2, EFunctionSpace functionSpace = FUNCTION_SPACE_MAX) : basis1_(basis1),basis2_(basis2) { + this->functionSpace_ = functionSpace; + Basis_TensorBasis* basis1AsTensor = dynamic_cast(basis1_.get()); if (basis1AsTensor) { @@ -763,7 +764,7 @@ struct OperatorTensorDecomposition Note that only the basic exact-sequence operators are supported at the moment: VALUE, GRAD, DIV, CURL. */ - virtual BasisValues allocateBasisValues( TensorPoints points, const EOperator operatorType = OPERATOR_VALUE) const override + virtual BasisValues allocateBasisValues( TensorPoints points, const EOperator operatorType = OPERATOR_VALUE) const override { const bool operatorSupported = (operatorType == OPERATOR_VALUE) || (operatorType == OPERATOR_GRAD) || (operatorType == OPERATOR_CURL) || (operatorType == OPERATOR_DIV); INTREPID2_TEST_FOR_EXCEPTION(!operatorSupported, std::invalid_argument, "operator is not supported by allocateBasisValues"); @@ -777,31 +778,31 @@ struct OperatorTensorDecomposition if (useVectorData) { const int numFamilies = 1; - std::vector< std::vector > > vectorComponents(numFamilies, std::vector >(numVectorComponents)); + std::vector< std::vector > > vectorComponents(numFamilies, std::vector >(numVectorComponents)); const int familyOrdinal = 0; for (ordinal_type vectorComponentOrdinal=0; vectorComponentOrdinal > componentData; + std::vector< Data > componentData; for (ordinal_type r=0; rallocateOutputView(numComponentPoints, op); - componentData.push_back(Data(componentView)); + componentData.push_back(Data(componentView)); } - vectorComponents[familyOrdinal][vectorComponentOrdinal] = TensorData(componentData); + vectorComponents[familyOrdinal][vectorComponentOrdinal] = TensorData(componentData); } } - VectorData vectorData(vectorComponents); - return BasisValues(vectorData); + VectorData vectorData(vectorComponents); + return BasisValues(vectorData); } else { // TensorData: single tensor product - std::vector< Data > componentData; + std::vector< Data > componentData; const ordinal_type vectorComponentOrdinal = 0; for (ordinal_type r=0; r extents {componentView.extent_int(0), componentView.extent_int(1), 1,1,1,1,1}; Kokkos::Array variationType {GENERAL, GENERAL, CONSTANT, CONSTANT, CONSTANT, CONSTANT, CONSTANT }; - componentData.push_back(Data(componentView, rank, extents, variationType)); + componentData.push_back(Data(componentView, rank, extents, variationType)); } - TensorData tensorData(componentData); + TensorData tensorData(componentData); - std::vector< TensorData > tensorDataEntries {tensorData}; - return BasisValues(tensorDataEntries); + std::vector< TensorData > tensorDataEntries {tensorData}; + return BasisValues(tensorDataEntries); } } // 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.) - using BasisSuper::getValues; + using BasisBase::getValues; /** \brief Method to extract component points from composite points. \param [in] inputPoints - points defined on the composite cell topology @@ -885,14 +886,14 @@ struct OperatorTensorDecomposition Note that getDofCoords() is not supported by all bases; in particular, hierarchical bases do not generally support this. */ - virtual void getDofCoords( typename BasisSuper::ScalarViewType dofCoords ) const override + virtual void getDofCoords( typename BasisBase::ScalarViewType dofCoords ) const override { int spaceDim1 = basis1_->getBaseCellTopology().getDimension(); int spaceDim2 = basis2_->getBaseCellTopology().getDimension(); - using ValueType = typename BasisSuper::ScalarViewType::value_type; - using ResultLayout = typename DeduceLayout< typename BasisSuper::ScalarViewType >::result_layout; - using DeviceType = typename BasisSuper::ScalarViewType::device_type; + using ValueType = typename BasisBase::ScalarViewType::value_type; + using ResultLayout = typename DeduceLayout< typename BasisBase::ScalarViewType >::result_layout; + using DeviceType = typename BasisBase::ScalarViewType::device_type; using ViewType = Kokkos::DynRankView; const ordinal_type basisCardinality1 = basis1_->getCardinality(); @@ -930,11 +931,11 @@ struct OperatorTensorDecomposition Note that getDofCoeffs() is not supported by all bases; in particular, hierarchical bases do not generally support this. */ - virtual void getDofCoeffs( typename BasisSuper::ScalarViewType dofCoeffs ) const override + virtual void getDofCoeffs( typename BasisBase::ScalarViewType dofCoeffs ) const override { - using ValueType = typename BasisSuper::ScalarViewType::value_type; - using ResultLayout = typename DeduceLayout< typename BasisSuper::ScalarViewType >::result_layout; - using DeviceType = typename BasisSuper::ScalarViewType::device_type; + using ValueType = typename BasisBase::ScalarViewType::value_type; + using ResultLayout = typename DeduceLayout< typename BasisBase::ScalarViewType >::result_layout; + using DeviceType = typename BasisBase::ScalarViewType::device_type; using ViewType = Kokkos::DynRankView; ViewType dofCoeffs1("dofCoeffs1",basis1_->getCardinality()); @@ -1039,8 +1040,8 @@ struct OperatorTensorDecomposition */ virtual void - getValues( BasisValues outputValues, - const TensorPoints inputPoints, + getValues( BasisValues outputValues, + const TensorPoints inputPoints, const EOperator operatorType = OPERATOR_VALUE ) const override { OperatorTensorDecomposition operatorDecomposition = getOperatorDecomposition(operatorType); @@ -1065,7 +1066,7 @@ struct OperatorTensorDecomposition PointViewType pointView = inputPoints.getTensorComponent(basisOrdinal); // Data stores things in fixed-rank Kokkos::View, but Basis requires DynRankView. We allocate a temporary DynRankView, then copy back to Data. - const Data & outputData = tensorData.getTensorComponent(basisOrdinal); + const Data & outputData = tensorData.getTensorComponent(basisOrdinal); auto basisValueView = outputData.getUnderlyingView(); tensorComponents_[basisOrdinal]->getValues(basisValueView, pointView, op); @@ -1076,7 +1077,7 @@ struct OperatorTensorDecomposition { if (basisValueView.rank() == 2) { - auto policy = Kokkos::MDRangePolicy>({0,0},{basisValueView.extent_int(0),basisValueView.extent_int(1)}); + auto policy = Kokkos::MDRangePolicy>({0,0},{basisValueView.extent_int(0),basisValueView.extent_int(1)}); Kokkos::parallel_for("multiply basisValueView by weight", policy, KOKKOS_LAMBDA (const int &fieldOrdinal, const int &pointOrdinal) { basisValueView(fieldOrdinal,pointOrdinal) *= weight; @@ -1084,7 +1085,7 @@ struct OperatorTensorDecomposition } else if (basisValueView.rank() == 3) { - auto policy = Kokkos::MDRangePolicy>({0,0,0},{basisValueView.extent_int(0),basisValueView.extent_int(1),basisValueView.extent_int(2)}); + auto policy = Kokkos::MDRangePolicy>({0,0,0},{basisValueView.extent_int(0),basisValueView.extent_int(1),basisValueView.extent_int(2)}); Kokkos::parallel_for("multiply basisValueView by weight", policy, KOKKOS_LAMBDA (const int &fieldOrdinal, const int &pointOrdinal, const int &d) { basisValueView(fieldOrdinal,pointOrdinal,d) *= weight; @@ -1349,7 +1350,16 @@ struct OperatorTensorDecomposition FunctorType functor(outputValues, outputValues1, outputValues2, tensorPoints, weight); Kokkos::parallel_for( policy , functor, "TensorViewFunctor"); } - }; + + /** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this. + + \return Pointer to the new Basis object. + */ + virtual HostBasisPtr + getHostBasis() const override { + TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "TensorBasis subclasses must override getHostBasis"); + } + }; // Basis_TensorBasis /** \struct Intrepid2::TensorBasis3_Functor \brief Functor for computing values for the TensorBasis3 class. @@ -1600,28 +1610,24 @@ struct OperatorTensorDecomposition } } } - }; + }; // TensorBasis3_Functor - template + template class Basis_TensorBasis3 - : public Basis_TensorBasis + : public Basis_TensorBasis { - using TensorBasis = Basis_TensorBasis; - + using BasisBase = BasisBaseClass; + using TensorBasis = Basis_TensorBasis; public: - using OutputViewType = typename TensorBasis::OutputViewType; - using PointViewType = typename TensorBasis::PointViewType; - using ScalarViewType = typename TensorBasis::ScalarViewType; + using OutputViewType = typename BasisBase::OutputViewType; + using PointViewType = typename BasisBase::PointViewType; + using ScalarViewType = typename BasisBase::ScalarViewType; - using OutputValueType = typename TensorBasis::OutputValueType; - using PointValueType = typename TensorBasis::PointValueType; + using OutputValueType = typename BasisBase::OutputValueType; + using PointValueType = typename BasisBase::PointValueType; - using BasisSuper = ::Intrepid2::Basis; - using BasisPtr = Teuchos::RCP; + using BasisPtr = Teuchos::RCP; protected: BasisPtr basis1_; BasisPtr basis2_; @@ -1847,16 +1853,27 @@ struct OperatorTensorDecomposition basis2_->getValues(outputValues2,inputPoints2,operatorType2); basis3_->getValues(outputValues3,inputPoints3,operatorType3); - const int outputVectorSize = getVectorSizeForHierarchicalParallelism(); - const int pointVectorSize = getVectorSizeForHierarchicalParallelism(); + const int outputVectorSize = getVectorSizeForHierarchicalParallelism(); + const int pointVectorSize = getVectorSizeForHierarchicalParallelism(); const int vectorSize = std::max(outputVectorSize,pointVectorSize); + using ExecutionSpace = typename BasisBase::ExecutionSpace; + auto policy = Kokkos::TeamPolicy(basisCardinality1,Kokkos::AUTO(),vectorSize); - using FunctorType = TensorBasis3_Functor; + using FunctorType = TensorBasis3_Functor; FunctorType functor(outputValues, outputValues1, outputValues2, outputValues3, tensorPoints, weight); Kokkos::parallel_for( policy , functor, "TensorBasis3_Functor"); } + + /** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this. + + \return Pointer to the new Basis object. + */ + virtual HostBasisPtr + getHostBasis() const override { + TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "TensorBasis3 subclasses must override getHostBasis"); + } }; } // end namespace Intrepid2 diff --git a/packages/intrepid2/src/Shared/Intrepid2_TestUtils.hpp b/packages/intrepid2/src/Shared/Intrepid2_TestUtils.hpp index 8eef69416b63..3e282a43c56d 100644 --- a/packages/intrepid2/src/Shared/Intrepid2_TestUtils.hpp +++ b/packages/intrepid2/src/Shared/Intrepid2_TestUtils.hpp @@ -67,6 +67,13 @@ namespace Intrepid2 //! Maximum number of derivatives to track for Fad types in tests. constexpr int MAX_FAD_DERIVATIVES_FOR_TESTS = 3; + //! Default Kokkos::Device to use for tests; depends on platform +#ifdef KOKKOS_ENABLE_CUDA + using DefaultTestDeviceType = Kokkos::Device; +#else + using DefaultTestDeviceType = typename Kokkos::DefaultExecutionSpace::device_type; +#endif + //! Use Teuchos small number determination on host; pass this to Intrepid2::relErr() on device. template const typename Teuchos::ScalarTraits::magnitudeType diff --git a/packages/intrepid2/unit-test/MonolithicExecutable/BasisValuesTests.cpp b/packages/intrepid2/unit-test/MonolithicExecutable/BasisValuesTests.cpp index 73b70b564626..997bf2267a7e 100644 --- a/packages/intrepid2/unit-test/MonolithicExecutable/BasisValuesTests.cpp +++ b/packages/intrepid2/unit-test/MonolithicExecutable/BasisValuesTests.cpp @@ -93,6 +93,13 @@ namespace TensorPoints tensorPoints; TensorData tensorWeights; + using HostPointViewType = Kokkos::DynRankView; + using HostWeightViewType = Kokkos::DynRankView; + auto hostQuadrature = cub_factory.create(cellTopoKey, quadratureDegree); + HostPointViewType hostPoints("quadrature points ref cell - host", numRefPoints, spaceDim); + HostWeightViewType hostWeights("quadrature weights ref cell - host", numRefPoints); + hostQuadrature->getCubature(hostPoints, hostWeights); + using CubatureTensorType = CubatureTensor; CubatureTensorType* tensorQuadrature = dynamic_cast(quadrature.get()); @@ -127,14 +134,22 @@ namespace testFloatingEquality2(points,tensorPoints, relTol, absTol, out, success, "points", "tensorPoints"); + auto hostBasisPtr = basis.getHostBasis(); + for (const auto &op : opsToTest) { auto basisValuesView = basis.allocateOutputView(numRefPoints, op); auto basisValues = basis.allocateBasisValues(tensorPoints, op); + auto hostBasisView = hostBasisPtr->allocateOutputView(numRefPoints, op); + basis.getValues(basisValuesView, points, op); basis.getValues(basisValues, tensorPoints, op); + // copy basisValuesView to host for hostBasis comparison + auto basisValuesViewHostMirror = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), basisValuesView); + hostBasisPtr->getValues(hostBasisView, hostPoints, op); + out << "Comparing getValues() results for " << EOperatorToString(op) << std::endl; TEST_EQUALITY(basisValues.rank(), basisValuesView.rank()); @@ -152,6 +167,16 @@ namespace printFunctor2(basisValues, out, "basisValues"); printFunctor2(basisValuesView, out, "basisValuesView"); } + + localSuccess = true; + testFloatingEquality2(basisValuesView, basisValues, relTol, absTol, out, localSuccess, "DynRankView path - device", "DynRankView path - host"); + success = success && localSuccess; + + if (!localSuccess) + { + printFunctor2(hostBasisView, out, "hostBasisView"); + printFunctor2(basisValuesViewHostMirror, out, "basisValuesViewHostMirror"); + } } else if (basisValuesView.rank() == 3) { @@ -164,6 +189,16 @@ namespace printFunctor3(basisValues, out, "basisValues"); printFunctor3(basisValuesView, out, "basisValuesView"); } + + localSuccess = true; + testFloatingEquality3(basisValuesView, basisValues, relTol, absTol, out, localSuccess, "DynRankView path - device", "DynRankView path - host"); + success = success && localSuccess; + + if (!localSuccess) + { + printFunctor3(hostBasisView, out, "hostBasisView"); + printFunctor3(basisValuesViewHostMirror, out, "basisValuesViewHostMirror"); + } } else { @@ -201,7 +236,8 @@ namespace TEUCHOS_UNIT_TEST( BasisValues, HierarchicalHGRAD_LINE ) { - using Basis = HierarchicalBasisFamily<>::HGRAD_LINE; + using DeviceType = Intrepid2::DefaultTestDeviceType; + using Basis = HierarchicalBasisFamily::HGRAD_LINE; // for now, the BasisValues path only supports the standard exact-sequence operators std::vector opsToTest {OPERATOR_VALUE, OPERATOR_GRAD}; @@ -218,7 +254,8 @@ namespace TEUCHOS_UNIT_TEST( BasisValues, HierarchicalHGRAD_QUAD ) { - using Basis = HierarchicalBasisFamily<>::HGRAD_QUAD; + using DeviceType = Intrepid2::DefaultTestDeviceType; + using Basis = HierarchicalBasisFamily::HGRAD_QUAD; // for now, the BasisValues path only supports the standard exact-sequence operators std::vector opsToTest {OPERATOR_VALUE, OPERATOR_GRAD}; @@ -235,7 +272,8 @@ namespace TEUCHOS_UNIT_TEST( BasisValues, HierarchicalHGRAD_TRI ) { - using Basis = HierarchicalBasisFamily<>::HGRAD_TRI; + using DeviceType = Intrepid2::DefaultTestDeviceType; + using Basis = HierarchicalBasisFamily::HGRAD_TRI; // for now, the BasisValues path only supports the standard exact-sequence operators std::vector opsToTest {OPERATOR_VALUE, OPERATOR_GRAD}; @@ -252,7 +290,8 @@ namespace TEUCHOS_UNIT_TEST( BasisValues, HierarchicalHGRAD_HEX ) { - using Basis = HierarchicalBasisFamily<>::HGRAD_HEX; + using DeviceType = Intrepid2::DefaultTestDeviceType; + using Basis = HierarchicalBasisFamily::HGRAD_HEX; // for now, the BasisValues path only supports the standard exact-sequence operators std::vector opsToTest {OPERATOR_VALUE, OPERATOR_GRAD}; @@ -269,7 +308,8 @@ namespace TEUCHOS_UNIT_TEST( BasisValues, HierarchicalHGRAD_TET ) { - using Basis = HierarchicalBasisFamily<>::HGRAD_TET; + using DeviceType = Intrepid2::DefaultTestDeviceType; + using Basis = HierarchicalBasisFamily::HGRAD_TET; // for now, the BasisValues path only supports the standard exact-sequence operators std::vector opsToTest {OPERATOR_VALUE, OPERATOR_GRAD}; @@ -286,7 +326,8 @@ namespace TEUCHOS_UNIT_TEST( BasisValues, HierarchicalHDIV_QUAD ) { - using Basis = HierarchicalBasisFamily<>::HDIV_QUAD; + using DeviceType = Intrepid2::DefaultTestDeviceType; + using Basis = HierarchicalBasisFamily::HDIV_QUAD; // for now, the BasisValues path only supports the standard exact-sequence operators std::vector opsToTest {OPERATOR_VALUE, OPERATOR_DIV}; @@ -303,7 +344,8 @@ namespace TEUCHOS_UNIT_TEST( BasisValues, NodalHDIV_TRI ) { - using Basis = NodalBasisFamily<>::HDIV_TRI; // Hierarchical basis family does not yet support HDIV_TRI + using DeviceType = Intrepid2::DefaultTestDeviceType; + using Basis = NodalBasisFamily::HDIV_TRI; // Hierarchical basis family does not yet support HDIV_TRI // for now, the BasisValues path only supports the standard exact-sequence operators std::vector opsToTest {OPERATOR_VALUE, OPERATOR_DIV}; @@ -320,7 +362,8 @@ namespace TEUCHOS_UNIT_TEST( BasisValues, HierarchicalHDIV_HEX ) { - using Basis = HierarchicalBasisFamily<>::HDIV_HEX; + using DeviceType = Intrepid2::DefaultTestDeviceType; + using Basis = HierarchicalBasisFamily::HDIV_HEX; // for now, the BasisValues path only supports the standard exact-sequence operators std::vector opsToTest {OPERATOR_VALUE, OPERATOR_DIV}; @@ -337,7 +380,8 @@ namespace TEUCHOS_UNIT_TEST( BasisValues, NodalHDIV_TET ) { - using Basis = NodalBasisFamily<>::HDIV_TET; // Hierarchical basis family does not yet support HDIV_TET + using DeviceType = Intrepid2::DefaultTestDeviceType; + using Basis = NodalBasisFamily::HDIV_TET; // Hierarchical basis family does not yet support HDIV_TET // for now, the BasisValues path only supports the standard exact-sequence operators std::vector opsToTest {OPERATOR_VALUE, OPERATOR_DIV}; @@ -354,7 +398,8 @@ namespace TEUCHOS_UNIT_TEST( BasisValues, HierarchicalHCURL_QUAD ) { - using Basis = HierarchicalBasisFamily<>::HCURL_QUAD; + using DeviceType = Intrepid2::DefaultTestDeviceType; + using Basis = HierarchicalBasisFamily::HCURL_QUAD; // for now, the BasisValues path only supports the standard exact-sequence operators std::vector opsToTest {OPERATOR_VALUE, OPERATOR_CURL}; @@ -371,7 +416,8 @@ namespace TEUCHOS_UNIT_TEST( BasisValues, NodalHCURL_TRI ) { - using Basis = NodalBasisFamily<>::HCURL_TRI; // Hierarchical basis family does not yet support HCURL_TRI + using DeviceType = Intrepid2::DefaultTestDeviceType; + using Basis = NodalBasisFamily::HCURL_TRI; // Hierarchical basis family does not yet support HCURL_TRI // for now, the BasisValues path only supports the standard exact-sequence operators std::vector opsToTest {OPERATOR_VALUE, OPERATOR_CURL}; @@ -388,7 +434,8 @@ namespace TEUCHOS_UNIT_TEST( BasisValues, HierarchicalHCURL_HEX ) { - using Basis = HierarchicalBasisFamily<>::HCURL_HEX; + using DeviceType = Intrepid2::DefaultTestDeviceType; + using Basis = HierarchicalBasisFamily::HCURL_HEX; // for now, the BasisValues path only supports the standard exact-sequence operators std::vector opsToTest {OPERATOR_VALUE, OPERATOR_CURL}; @@ -403,9 +450,46 @@ namespace } } + TEUCHOS_UNIT_TEST( BasisValues, HierarchicalHVOL_QUAD ) + { + using DeviceType = Intrepid2::DefaultTestDeviceType; + using Basis = HierarchicalBasisFamily::HVOL_QUAD; + + // for now, the BasisValues path only supports the standard exact-sequence operators + std::vector opsToTest {OPERATOR_VALUE}; + + const double relTol=1e-13; + const double absTol=1e-13; + + for (int polyOrder=1; polyOrder<5; polyOrder++) + { + Basis basis(polyOrder); + testGetValuesEquality(basis, opsToTest, relTol, absTol, out, success); + } + } + + TEUCHOS_UNIT_TEST( BasisValues, HierarchicalHVOL_HEX ) + { + using DeviceType = Intrepid2::DefaultTestDeviceType; + using Basis = HierarchicalBasisFamily::HVOL_HEX; + + // for now, the BasisValues path only supports the standard exact-sequence operators + std::vector opsToTest {OPERATOR_VALUE}; + + const double relTol=1e-13; + const double absTol=1e-13; + + for (int polyOrder=1; polyOrder<5; polyOrder++) + { + Basis basis(polyOrder); + testGetValuesEquality(basis, opsToTest, relTol, absTol, out, success); + } + } + TEUCHOS_UNIT_TEST( BasisValues, NodalHCURL_TET ) { - using Basis = NodalBasisFamily<>::HCURL_TET; // Hierarchical basis family does not yet support HCURL_TET + using DeviceType = Intrepid2::DefaultTestDeviceType; + using Basis = NodalBasisFamily::HCURL_TET; // Hierarchical basis family does not yet support HCURL_TET // for now, the BasisValues path only supports the standard exact-sequence operators std::vector opsToTest {OPERATOR_VALUE, OPERATOR_CURL}; @@ -422,7 +506,8 @@ namespace TEUCHOS_UNIT_TEST( BasisValues, NodalHGRAD_LINE ) { - using Basis = NodalBasisFamily<>::HGRAD_LINE; + using DeviceType = Intrepid2::DefaultTestDeviceType; + using Basis = NodalBasisFamily::HGRAD_LINE; // for now, the BasisValues path only supports the standard exact-sequence operators std::vector opsToTest {OPERATOR_VALUE, OPERATOR_GRAD};