From 4002c6b22e7cf5cc6ae30b4615f688475e1bcb00 Mon Sep 17 00:00:00 2001 From: Nate Roberts Date: Mon, 1 Mar 2021 15:17:51 -0700 Subject: [PATCH] Implement getHostBasis() in newer Basis subclasses (#8811) Intrepid2: added getHostBasis() implementations for several Basis subclasses, along with support for DeviceType template arguments in the same. Specifically, the following classes now support getHostBasis and DeviceType template arguments: - DerivedBasis_HCURL_HEX - DerivedBasis_HCURL_QUAD - DerivedBasis_HDIV_HEX - DerivedBasis_HDIV_QUAD - DerivedBasis_HGRAD_HEX - DerivedBasis_HGRAD_QUAD - DerivedBasis_HVOL_HEX - DerivedBasis_HVOL_QUAD - HGRAD_LINE_Cn_FEM - HVOL_LINE_Cn_FEM - IntegratedLegendreBasis_HGRAD_LINE - LegendreBasis_HVOL_LINE Other changes: - added typedef Intrepid2::DefaultTestDeviceType to TestUtils; this is UVM-free on CUDA, and Kokkos::DefaultExecutionSpace::device_type otherwise. - Made use of DefaultTestDeviceType as well as getHostBasis in the BasisValues test group (in MonolithicExecutable) - revised template arguments of Basis_DirectSumBasis, Basis_TensorBasis, andBasis_TensorBasis3; these now depend only on the base Basis class - added typedef for BasisBase in several Basis subclasses - added typedef for HostBasis in several Basis subclasses (specifically, those defined on the line which are used in DerivedBasis instantiations) The modified tests within MonolithicExecutable exercise getHostBasis() on all modified bases. These also test the modified bases with the UVM-free DeviceType when run under CUDA. --- .../Discretization/Basis/Intrepid2_Basis.hpp | 15 +- .../Basis/Intrepid2_DerivedBasisFamily.hpp | 2 +- .../Intrepid2_DerivedBasis_HCURL_HEX.hpp | 43 +++--- .../Intrepid2_DerivedBasis_HCURL_QUAD.hpp | 32 +++- .../Basis/Intrepid2_DerivedBasis_HDIV_HEX.hpp | 40 +++-- .../Intrepid2_DerivedBasis_HDIV_QUAD.hpp | 33 +++- .../Intrepid2_DerivedBasis_HGRAD_HEX.hpp | 22 ++- .../Intrepid2_DerivedBasis_HGRAD_QUAD.hpp | 22 ++- .../Basis/Intrepid2_DerivedBasis_HVOL_HEX.hpp | 25 ++- .../Intrepid2_DerivedBasis_HVOL_QUAD.hpp | 30 +++- .../Basis/Intrepid2_DirectSumBasis.hpp | 48 +++--- .../Basis/Intrepid2_HGRAD_LINE_C1_FEM.hpp | 20 +-- .../Basis/Intrepid2_HGRAD_LINE_Cn_FEM.hpp | 32 ++-- .../Basis/Intrepid2_HGRAD_LINE_Cn_FEMDef.hpp | 1 + .../Basis/Intrepid2_HVOL_LINE_Cn_FEM.hpp | 21 +-- .../Intrepid2_HierarchicalBasisFamily.hpp | 46 +++--- ...id2_IntegratedLegendreBasis_HGRAD_LINE.hpp | 43 ++++-- ...pid2_IntegratedLegendreBasis_HGRAD_TET.hpp | 45 ++++-- ...pid2_IntegratedLegendreBasis_HGRAD_TRI.hpp | 43 ++++-- .../Intrepid2_LegendreBasis_HVOL_LINE.hpp | 41 +++-- .../Basis/Intrepid2_TensorBasis.hpp | 145 ++++++++++-------- .../src/Shared/Intrepid2_TestUtils.hpp | 7 + .../MonolithicExecutable/BasisValuesTests.cpp | 113 ++++++++++++-- 23 files changed, 596 insertions(+), 273 deletions(-) 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};