diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HCURL_WEDGE.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HCURL_WEDGE.hpp index ec5f5874bb72..e65e483dc63f 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HCURL_WEDGE.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HCURL_WEDGE.hpp @@ -229,8 +229,8 @@ namespace Intrepid2 Note that getDofCoeffs() is supported only for Lagrangian bases. */ virtual void getDofCoeffs( ScalarViewType dofCoeffs ) const override { - auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),0); - auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),1); + auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),std::make_pair(0,2)); + auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),2); this->TensorBasis::getDofCoeffs(dofCoeffs1); Kokkos::deep_copy(dofCoeffs2,0.0); } @@ -376,8 +376,8 @@ namespace Intrepid2 Note that getDofCoeffs() is supported only for Lagrangian bases. */ virtual void getDofCoeffs( ScalarViewType dofCoeffs ) const override { - auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),0); - auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),1); + auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),std::make_pair(0,2)); + auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),2); Kokkos::deep_copy(dofCoeffs1,0.0); this->TensorBasis::getDofCoeffs(dofCoeffs2); } @@ -448,6 +448,46 @@ namespace Intrepid2 return name_.c_str(); } + /** \brief returns the basis associated to a subCell. + + The bases of the subCell are the restriction to the subCell + of the bases of the parent cell. + TODO: test this method when different orders are used in different directions + \param [in] subCellDim - dimension of subCell + \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 + */ + Teuchos::RCP + getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override + { + using LineBasis = HVOL_LINE; + using TriCurlBasis = HCURL_TRI; + using QuadCurlBasis = Basis_Derived_HCURL_QUAD; + if(subCellDim == 1) { + if(subCellOrd < 6) //edges associated to basal and upper triagular faces + return Teuchos::rcp( new LineBasis(order_xy_-1, pointType_) ); + else + return Teuchos::rcp( new LineBasis(order_z_-1, pointType_) ); + } + else if(subCellDim == 2) { + switch(subCellOrd) { + case 0: + return Teuchos::rcp( new QuadCurlBasis(order_xy_, order_z_, pointType_) ); + case 1: + return Teuchos::rcp( new QuadCurlBasis(order_xy_, order_z_, pointType_) ); + case 2: + return Teuchos::rcp( new QuadCurlBasis(order_z_, order_xy_, pointType_) ); + case 3: + return Teuchos::rcp( new TriCurlBasis(order_xy_, pointType_) ); + case 4: + return Teuchos::rcp( new TriCurlBasis(order_xy_, pointType_) ); + default: + INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"subCellOrd is out of bounds"); + } + } + INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"subCellDim is 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. diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HDIV_WEDGE.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HDIV_WEDGE.hpp index abbb94bae10c..4e5765582bd8 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HDIV_WEDGE.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HDIV_WEDGE.hpp @@ -197,10 +197,10 @@ namespace Intrepid2 Note that getDofCoeffs() is supported only for Lagrangian bases. */ virtual void getDofCoeffs( ScalarViewType dofCoeffs ) const override { - auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),0); - auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),1); - Kokkos::deep_copy(dofCoeffs1,0.0); - this->TensorBasis::getDofCoeffs(dofCoeffs2); + auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(), std::make_pair(0,2)); + auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(), 2); + this->TensorBasis::getDofCoeffs(dofCoeffs1); + Kokkos::deep_copy(dofCoeffs2,0.0); } }; @@ -329,10 +329,10 @@ namespace Intrepid2 Note that getDofCoeffs() is supported only for Lagrangian bases. */ virtual void getDofCoeffs( ScalarViewType dofCoeffs ) const override { - auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),0); - auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),1); - this->TensorBasis::getDofCoeffs(dofCoeffs1); - Kokkos::deep_copy(dofCoeffs2,0.0); + auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(), std::make_pair(0,2)); + auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(), 2); + Kokkos::deep_copy(dofCoeffs1,0.0); + this->TensorBasis::getDofCoeffs(dofCoeffs2); } }; @@ -402,6 +402,41 @@ namespace Intrepid2 return name_.c_str(); } + + /** \brief returns the basis associated to a subCell. + + The bases of the subCell are the restriction to the subCell + of the bases of the parent cell. + TODO: test this method when different orders are used in different directions + \param [in] subCellDim - dimension of subCell + \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 + */ + Teuchos::RCP + getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override + { + using QuadBasis = Basis_Derived_HVOL_QUAD; + using TriBasis = HVOL_TRI; + + if(subCellDim == 2) { + switch(subCellOrd) { + case 0: + return Teuchos::rcp( new QuadBasis(order_xy_-1, order_z_-1, pointType_) ); + case 1: + return Teuchos::rcp( new QuadBasis(order_xy_-1, order_z_-1, pointType_) ); + case 2: + return Teuchos::rcp( new QuadBasis(order_z_-1, order_xy_-1, pointType_) ); + case 3: + return Teuchos::rcp( new TriBasis(order_xy_-1, pointType_) ); + case 4: + return Teuchos::rcp( new TriBasis(order_xy_-1, pointType_) ); + default: + INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"subCellOrd is out of bounds"); + } + } + INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"subCellDim is 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. diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_WEDGE.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_WEDGE.hpp index 47ae49a7aaeb..c2556a7e1f9b 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_WEDGE.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_WEDGE.hpp @@ -52,6 +52,7 @@ #define Intrepid2_DerivedBasis_HGRAD_WEDGE_h #include "Intrepid2_TensorBasis.hpp" +#include "Intrepid2_DerivedBasis_HGRAD_QUAD.hpp" namespace Intrepid2 { @@ -137,27 +138,29 @@ namespace Intrepid2 Teuchos::RCP getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override { - INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Method not yet implemented"); - if(subCellDim == 2) { + if(subCellDim == 1) { + if(subCellOrd < 6) //edges associated to basal and upper triagular faces + return Teuchos::rcp( new LineBasis(order_xy_, pointType_) ); + else + return Teuchos::rcp( new LineBasis(order_z_, pointType_) ); + } + else if(subCellDim == 2) { switch(subCellOrd) { - // TODO: check this subcell numbering -- which are the tris, and which are the quads? - // (one way: use the shards CellTopology to look this up…) - case 0: - case 1: - return Teuchos::rcp( new TriBasis(order_xy_, pointType_) ); - case 2: - case 3: - case 4: - // TODO: check what order the poly order arguments should go in... - // (one way: use the shards CellTopology to look this up…) - // TODO: define QuadBasis somehow, somewhere. (Probably use DerivedBasis_HGRAD_QUAD.) -// return Teuchos::rcp( new QuadBasis(order_xy_, order_z_, pointType_) ); - return Teuchos::null; + case 0: + return Teuchos::rcp( new Intrepid2::Basis_Derived_HGRAD_QUAD(order_xy_, order_z_, pointType_) ); + case 1: + return Teuchos::rcp( new Intrepid2::Basis_Derived_HGRAD_QUAD(order_xy_, order_z_, pointType_) ); + case 2: + return Teuchos::rcp( new Intrepid2::Basis_Derived_HGRAD_QUAD(order_z_, order_xy_, pointType_) ); + case 3: + return Teuchos::rcp( new TriBasis(order_xy_, pointType_) ); + case 4: + return Teuchos::rcp( new TriBasis(order_xy_, pointType_) ); + default: + INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"subCellOrd is out of bounds"); } - return Teuchos::null; - } - - INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds"); + } + INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"subCellDim is 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. diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_HEX_C2_FEM.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_HEX_C2_FEM.hpp index 5074dab58dd3..188078e3fbaa 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_HEX_C2_FEM.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_HEX_C2_FEM.hpp @@ -54,12 +54,20 @@ namespace Intrepid2 { - /** \class Intrepid2::Basis_HGRAD_HEX_C2_FEM + /** \class Intrepid2::Basis_HGRAD_HEX_DEG2_FEM \brief Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell - Implements Lagrangian basis of degree 2 on the reference Hexahedron cell. The basis has - cardinality 27 and spans a COMPLETE tri-quadratic polynomial space. Basis functions are dual - to a unisolvent set of degrees-of-freedom (DoF) defined and enumerated as follows: + Implements Lagrangian basis of degree 2 on the reference Hexahedron cell. + + When the serendipity template argument is false, the basis has + cardinality 27 and spans a COMPLETE tri-quadratic polynomial space. + Note, Basis_HGRAD_HEX_C2_FEM = Basis_HGRAD_HEX_DEG2_FEM + + When the serendipity template argument is true, the basis has + cardinality 20 and spans an INCOMPLETE tri-quadratic polynomial space (Dofs are associated only to vertices and edges). + Note, Basis_HGRAD_HEX_I2_Serendipity_FEM = Basis_HGRAD_HEX_DEG2_FEM + + Basis functions are dual to a unisolvent set of degrees-of-freedom (DoF) defined and enumerated as follows: \verbatim ================================================================================================= @@ -139,13 +147,14 @@ namespace Intrepid2 { namespace Impl { /** - \brief See Intrepid2::Basis_HGRAD_HEX_C2_FEM + \brief See Intrepid2::Basis_HGRAD_HEX_DEG2_FEM */ - class Basis_HGRAD_HEX_C2_FEM { + template + class Basis_HGRAD_HEX_DEG2_FEM { public: - typedef struct Hexahedron<27> cell_topology_type; + typedef struct Hexahedron cell_topology_type; /** - \brief See Intrepid2::Basis_HGRAD_HEX_C2_FEM + \brief See Intrepid2::Basis_HGRAD_HEX_DEG2_FEM */ template struct Serial { @@ -167,7 +176,7 @@ namespace Intrepid2 { const EOperator operatorType); /** - \brief See Intrepid2::Basis_HGRAD_HEX_C2_FEM + \brief See Intrepid2::Basis_HGRAD_HEX_DEG2_FEM */ template>> ERROR: (Intrepid2::Basis_HGRAD_HEX_C2_FEM::Serial::getValues) operator is not supported"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::Serial::getValues) operator is not supported"); } } } @@ -217,10 +226,11 @@ namespace Intrepid2 { }; } - template - class Basis_HGRAD_HEX_C2_FEM : public Basis { + template + class Basis_HGRAD_HEX_DEG2_FEM : public Basis { public: using OrdinalTypeArray1DHost = typename Basis::OrdinalTypeArray1DHost; using OrdinalTypeArray2DHost = typename Basis::OrdinalTypeArray2DHost; @@ -228,7 +238,7 @@ namespace Intrepid2 { /** \brief Constructor. */ - Basis_HGRAD_HEX_C2_FEM(); + Basis_HGRAD_HEX_DEG2_FEM(); using OutputViewType = typename Basis::OutputViewType; using PointViewType = typename Basis::PointViewType; @@ -249,10 +259,16 @@ namespace Intrepid2 { this->getBaseCellTopology(), this->getCardinality() ); #endif - Impl::Basis_HGRAD_HEX_C2_FEM:: - getValues( outputValues, - inputPoints, - operatorType ); + if constexpr (serendipity) + Impl::Basis_HGRAD_HEX_DEG2_FEM:: + getValues( outputValues, + inputPoints, + operatorType ); + else + Impl::Basis_HGRAD_HEX_DEG2_FEM:: + getValues( outputValues, + inputPoints, + operatorType ); } virtual @@ -261,13 +277,13 @@ namespace Intrepid2 { #ifdef HAVE_INTREPID2_DEBUG // Verify rank of output array. INTREPID2_TEST_FOR_EXCEPTION( dofCoords.rank() != 2, std::invalid_argument, - ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_C2_FEM::getDofCoords) rank = 2 required for dofCoords array"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::getDofCoords) rank = 2 required for dofCoords array"); // Verify 0th dimension of output array. INTREPID2_TEST_FOR_EXCEPTION( static_cast(dofCoords.extent(0)) != this->getCardinality(), std::invalid_argument, - ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_C2_FEM::getDofCoords) mismatch in number of dof and 0th dimension of dofCoords array"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::getDofCoords) mismatch in number of dof and 0th dimension of dofCoords array"); // Verify 1st dimension of output array. INTREPID2_TEST_FOR_EXCEPTION( dofCoords.extent(1) != this->getBaseCellTopology().getDimension(), std::invalid_argument, - ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_C2_FEM::getDofCoords) incorrect reference cell (1st) dimension in dofCoords array"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::getDofCoords) incorrect reference cell (1st) dimension in dofCoords array"); #endif Kokkos::deep_copy(dofCoords, this->dofCoords_); } @@ -278,10 +294,10 @@ namespace Intrepid2 { #ifdef HAVE_INTREPID2_DEBUG // Verify rank of output array. INTREPID2_TEST_FOR_EXCEPTION( dofCoeffs.rank() != 1, std::invalid_argument, - ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_C2_FEM::getdofCoeffs) rank = 1 required for dofCoeffs array"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::getdofCoeffs) rank = 1 required for dofCoeffs array"); // Verify 0th dimension of output array. INTREPID2_TEST_FOR_EXCEPTION( static_cast(dofCoeffs.extent(0)) != this->getCardinality(), std::invalid_argument, - ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_C2_FEM::getdofCoeffs) mismatch in number of dof and 0th dimension of dofCoeffs array"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::getdofCoeffs) mismatch in number of dof and 0th dimension of dofCoeffs array"); #endif Kokkos::deep_copy(dofCoeffs, 1.0); } @@ -289,7 +305,10 @@ namespace Intrepid2 { virtual const char* getName() const override { - return "Intrepid2_HGRAD_HEX_C2_FEM"; + if constexpr (serendipity) + return "Intrepid2_HGRAD_HEX_I2_Serendipity_FEM"; + else + return "Intrepid2_HGRAD_HEX_C2_FEM"; } /** \brief returns the basis associated to a subCell. @@ -307,16 +326,23 @@ namespace Intrepid2 { Basis_HGRAD_LINE_C2_FEM()); } else if(subCellDim == 2) { return Teuchos::rcp(new - Basis_HGRAD_QUAD_C2_FEM()); + Basis_HGRAD_QUAD_DEG2_FEM()); } INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds"); } BasisPtr getHostBasis() const override{ - return Teuchos::rcp(new Basis_HGRAD_HEX_C2_FEM()); + return Teuchos::rcp(new Basis_HGRAD_HEX_DEG2_FEM()); } }; + + template + using Basis_HGRAD_HEX_C2_FEM = Basis_HGRAD_HEX_DEG2_FEM; + + template + using Basis_HGRAD_HEX_I2_Serendipity_FEM = Basis_HGRAD_HEX_DEG2_FEM; + }// namespace Intrepid2 #include "Intrepid2_HGRAD_HEX_C2_FEMDef.hpp" diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_HEX_C2_FEMDef.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_HEX_C2_FEMDef.hpp index 21ba58bd1eaa..67bed49b2fcc 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_HEX_C2_FEMDef.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_HEX_C2_FEMDef.hpp @@ -54,12 +54,13 @@ namespace Intrepid2 { // ------------------------------------------------------------------------------------- namespace Impl { + template template template KOKKOS_INLINE_FUNCTION void - Basis_HGRAD_HEX_C2_FEM::Serial:: + Basis_HGRAD_HEX_DEG2_FEM::Serial:: getValues( OutputViewType output, const inputViewType input ) { switch (opType) { @@ -89,13 +90,15 @@ namespace Intrepid2 { output.access(17) = 0.25*x*(1.+ x)*(1. - y)*(1. + y)*z*(1.+ z); output.access(18) = 0.25*(1. - x)*(1. + x)*y*(1.+ y)*z*(1.+ z); output.access(19) = 0.25*(-1. + x)*x*(1. - y)*(1. + y)*z*(1.+ z); - output.access(20) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z); - output.access(21) = 0.5*(1. - x)*(1. + x)*(1. - y)*(1. + y)*(-1. + z)*z; - output.access(22) = 0.5*(1. - x)*(1. + x)*(1. - y)*(1. + y)*z*(1.+ z); - output.access(23) = 0.5*(-1. + x)*x*(1. - y)*(1. + y)*(1. - z)*(1. + z); - output.access(24) = 0.5*x*(1.+ x)*(1. - y)*(1. + y)*(1. - z)*(1. + z); - output.access(25) = 0.5*(1. - x)*(1. + x)*(-1. + y)*y*(1. - z)*(1. + z); - output.access(26) = 0.5*(1. - x)*(1. + x)*y*(1.+ y)*(1. - z)*(1. + z); + if constexpr (!serendipity) { + output.access(20) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z); + output.access(21) = 0.5*(1. - x)*(1. + x)*(1. - y)*(1. + y)*(-1. + z)*z; + output.access(22) = 0.5*(1. - x)*(1. + x)*(1. - y)*(1. + y)*z*(1.+ z); + output.access(23) = 0.5*(-1. + x)*x*(1. - y)*(1. + y)*(1. - z)*(1. + z); + output.access(24) = 0.5*x*(1.+ x)*(1. - y)*(1. + y)*(1. - z)*(1. + z); + output.access(25) = 0.5*(1. - x)*(1. + x)*(-1. + y)*y*(1. - z)*(1. + z); + output.access(26) = 0.5*(1. - x)*(1. + x)*y*(1.+ y)*(1. - z)*(1. + z); + } break; } case OPERATOR_GRAD : @@ -185,33 +188,35 @@ namespace Intrepid2 { output.access(19, 1) = (-1. + x)*x*(-0.5*y)*z*(1. + z); output.access(19, 2) = (-1. + x)*x*(1. - y)*(1. + y)*(0.25 + 0.5*z); - output.access(20, 0) = -2.*x*(1. - y)*(1. + y)*(1. - z)*(1. + z); - output.access(20, 1) = (1. - x)*(1. + x)*(-2.*y)*(1. - z)*(1. + z); - output.access(20, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(-2.*z); + if constexpr (!serendipity) { + output.access(20, 0) = -2.*x*(1. - y)*(1. + y)*(1. - z)*(1. + z); + output.access(20, 1) = (1. - x)*(1. + x)*(-2.*y)*(1. - z)*(1. + z); + output.access(20, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(-2.*z); - output.access(21, 0) = -x*(1. - y)*(1. + y)*(-1. + z)*z; - output.access(21, 1) = (1. - x)*(1. + x)*(-y)*(-1. + z)*z; - output.access(21, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(-0.5 + z); + output.access(21, 0) = -x*(1. - y)*(1. + y)*(-1. + z)*z; + output.access(21, 1) = (1. - x)*(1. + x)*(-y)*(-1. + z)*z; + output.access(21, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(-0.5 + z); - output.access(22, 0) = -x*(1. - y)*(1. + y)*z*(1. + z); - output.access(22, 1) = (1. - x)*(1. + x)*(-y)*z*(1. + z); - output.access(22, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(0.5 + z); + output.access(22, 0) = -x*(1. - y)*(1. + y)*z*(1. + z); + output.access(22, 1) = (1. - x)*(1. + x)*(-y)*z*(1. + z); + output.access(22, 2) = (1. - x)*(1. + x)*(1. - y)*(1. + y)*(0.5 + z); - output.access(23, 0) = (-0.5 + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z); - output.access(23, 1) = (-1. + x)*x*(-y)*(1. - z)*(1. + z); - output.access(23, 2) = (-1. + x)*x*(1. - y)*(1. + y)*(-z); + output.access(23, 0) = (-0.5 + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z); + output.access(23, 1) = (-1. + x)*x*(-y)*(1. - z)*(1. + z); + output.access(23, 2) = (-1. + x)*x*(1. - y)*(1. + y)*(-z); - output.access(24, 0) = (0.5 + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z); - output.access(24, 1) = x*(1. + x)*(-y)*(1. - z)*(1. + z); - output.access(24, 2) = x*(1. + x)*(1. - y)*(1. + y)*(-z); + output.access(24, 0) = (0.5 + x)*(1. - y)*(1. + y)*(1. - z)*(1. + z); + output.access(24, 1) = x*(1. + x)*(-y)*(1. - z)*(1. + z); + output.access(24, 2) = x*(1. + x)*(1. - y)*(1. + y)*(-z); - output.access(25, 0) = -x*(-1. + y)*y*(1. - z)*(1. + z); - output.access(25, 1) = (1. - x)*(1. + x)*(-0.5 + y)*(1. - z)*(1. + z); - output.access(25, 2) = (1. - x)*(1. + x)*(-1. + y)*y*(-z); + output.access(25, 0) = -x*(-1. + y)*y*(1. - z)*(1. + z); + output.access(25, 1) = (1. - x)*(1. + x)*(-0.5 + y)*(1. - z)*(1. + z); + output.access(25, 2) = (1. - x)*(1. + x)*(-1. + y)*y*(-z); - output.access(26, 0) = -x*y*(1. + y)*(1. - z)*(1. + z); - output.access(26, 1) = (1. - x)*(1. + x)*(0.5 + y)*(1. - z)*(1. + z); - output.access(26, 2) = (1. - x)*(1. + x)*y*(1. + y)*(-z); + output.access(26, 0) = -x*y*(1. + y)*(1. - z)*(1. + z); + output.access(26, 1) = (1. - x)*(1. + x)*(0.5 + y)*(1. - z)*(1. + z); + output.access(26, 2) = (1. - x)*(1. + x)*y*(1. + y)*(-z); + } break; } case OPERATOR_D2 : { @@ -360,54 +365,56 @@ namespace Intrepid2 { output.access(19, 4) = (x*x)*(y*(-0.5 - z) ) + x*(y*(0.5 + z)); output.access(19, 5) = 0.5*(-1. + x)*x*(1. - y)*(1. + y); - output.access(20, 0) = -2.*(1. - y)*(1. + y)*(1. - z)*(1. + z); - output.access(20, 1) = -4.*x*y*(-1. + z*z); - output.access(20, 2) = x*((y*y)*(-4.*z) + 4.*z); - output.access(20, 3) = -2.*(1. - x)*(1. + x)*(1. - z)*(1. + z); - output.access(20, 4) = (x*x)*(y*(-4.*z) ) + y*(4.*z); - output.access(20, 5) = -2.*(1. - x)*(1. + x)*(1. - y)*(1. + y); - - output.access(21, 0) = -(1. - y)*(1. + y)*(-1. + z)*z; - output.access(21, 1) = (x*(-2.*y) )*z + (0. + x*(2.*y))*(z*z); - output.access(21, 2) = x*(1. - 2.*z + (y*y)*(-1. + 2.*z)); - output.access(21, 3) = -(1. - x)*(1. + x)*(-1. + z)*z; - output.access(21, 4) = y*(1. - 2.*z) + (x*x)*(y*(-1. + 2.*z)); - output.access(21, 5) = (1. - x)*(1. + x)*(1. - y)*(1. + y); - - output.access(22, 0) = -(1. - y)*(1. + y)*z*(1 + z); - output.access(22, 1) = (0. + x*(2.*y))*z + (0. + x*(2.*y))*(z*z); - output.access(22, 2) = x*(-1. - 2.*z + (y*y)*(1. + 2.*z)); - output.access(22, 3) = -(1. - x)*(1. + x)*z*(1 + z); - output.access(22, 4) = y*(-1. - 2.*z) + (x*x)*(y*(1. + 2.*z)); - output.access(22, 5) = (1. - x)*(1. + x)*(1. - y)*(1. + y); - - output.access(23, 0) = (1. - y)*(1. + y)*(1. - z)*(1. + z); - output.access(23, 1) = (-1. + 2.*x)*y*(-1. + z*z); - output.access(23, 2) = (-1. + 2.*x)*(-1. + y*y)*z; - output.access(23, 3) =-(-1. + x)*x*(1. - z)*(1. + z); - output.access(23, 4) = 2.*(-1. + x)*x*y*z; - output.access(23, 5) =-(-1. + x)*x*(1. - y)*(1. + y); - - output.access(24, 0) = (1. - y)*(1. + y)*(1. - z)*(1. + z); - output.access(24, 1) = (1. + 2.*x)*y*(-1. + z*z); - output.access(24, 2) = (1. + 2.*x)*(-1. + y*y)*z; - output.access(24, 3) = x*(1. + x)*(-1. + z)*(1. + z); - output.access(24, 4) = 2.*x*(1. + x)*y*z; - output.access(24, 5) = x*(1. + x)*(-1. + y)*(1. + y); - - output.access(25, 0) = -(-1. + y)*y*(1. - z)*(1. + z); - output.access(25, 1) = x*(-1. + 2.*y)*(-1. + z*z); - output.access(25, 2) = 2.*x*(-1. + y)*y*z; - output.access(25, 3) = (1. - x)*(1. + x)*(1. - z)*(1. + z); - output.access(25, 4) = (-1. + x*x)*(-1. + 2.*y)*z; - output.access(25, 5) =-(1. - x)*(1. + x)*(-1. + y)*y; - - output.access(26, 0) = y*(1. + y)*(-1. + z)*(1. + z); - output.access(26, 1) = x*(1. + 2.*y)*(-1. + z*z); - output.access(26, 2) = 2.*x*y*(1. + y)*z; - output.access(26, 3) = (-1. + x)*(1. + x)*(-1. + z)*(1. + z); - output.access(26, 4) = (-1. + x*x)*(1. + 2.*y)*z; - output.access(26, 5) = (-1. + x)*(1. + x)*y*(1. + y); + if constexpr (!serendipity) { + output.access(20, 0) = -2.*(1. - y)*(1. + y)*(1. - z)*(1. + z); + output.access(20, 1) = -4.*x*y*(-1. + z*z); + output.access(20, 2) = x*((y*y)*(-4.*z) + 4.*z); + output.access(20, 3) = -2.*(1. - x)*(1. + x)*(1. - z)*(1. + z); + output.access(20, 4) = (x*x)*(y*(-4.*z) ) + y*(4.*z); + output.access(20, 5) = -2.*(1. - x)*(1. + x)*(1. - y)*(1. + y); + + output.access(21, 0) = -(1. - y)*(1. + y)*(-1. + z)*z; + output.access(21, 1) = (x*(-2.*y) )*z + (0. + x*(2.*y))*(z*z); + output.access(21, 2) = x*(1. - 2.*z + (y*y)*(-1. + 2.*z)); + output.access(21, 3) = -(1. - x)*(1. + x)*(-1. + z)*z; + output.access(21, 4) = y*(1. - 2.*z) + (x*x)*(y*(-1. + 2.*z)); + output.access(21, 5) = (1. - x)*(1. + x)*(1. - y)*(1. + y); + + output.access(22, 0) = -(1. - y)*(1. + y)*z*(1 + z); + output.access(22, 1) = (0. + x*(2.*y))*z + (0. + x*(2.*y))*(z*z); + output.access(22, 2) = x*(-1. - 2.*z + (y*y)*(1. + 2.*z)); + output.access(22, 3) = -(1. - x)*(1. + x)*z*(1 + z); + output.access(22, 4) = y*(-1. - 2.*z) + (x*x)*(y*(1. + 2.*z)); + output.access(22, 5) = (1. - x)*(1. + x)*(1. - y)*(1. + y); + + output.access(23, 0) = (1. - y)*(1. + y)*(1. - z)*(1. + z); + output.access(23, 1) = (-1. + 2.*x)*y*(-1. + z*z); + output.access(23, 2) = (-1. + 2.*x)*(-1. + y*y)*z; + output.access(23, 3) =-(-1. + x)*x*(1. - z)*(1. + z); + output.access(23, 4) = 2.*(-1. + x)*x*y*z; + output.access(23, 5) =-(-1. + x)*x*(1. - y)*(1. + y); + + output.access(24, 0) = (1. - y)*(1. + y)*(1. - z)*(1. + z); + output.access(24, 1) = (1. + 2.*x)*y*(-1. + z*z); + output.access(24, 2) = (1. + 2.*x)*(-1. + y*y)*z; + output.access(24, 3) = x*(1. + x)*(-1. + z)*(1. + z); + output.access(24, 4) = 2.*x*(1. + x)*y*z; + output.access(24, 5) = x*(1. + x)*(-1. + y)*(1. + y); + + output.access(25, 0) = -(-1. + y)*y*(1. - z)*(1. + z); + output.access(25, 1) = x*(-1. + 2.*y)*(-1. + z*z); + output.access(25, 2) = 2.*x*(-1. + y)*y*z; + output.access(25, 3) = (1. - x)*(1. + x)*(1. - z)*(1. + z); + output.access(25, 4) = (-1. + x*x)*(-1. + 2.*y)*z; + output.access(25, 5) =-(1. - x)*(1. + x)*(-1. + y)*y; + + output.access(26, 0) = y*(1. + y)*(-1. + z)*(1. + z); + output.access(26, 1) = x*(1. + 2.*y)*(-1. + z*z); + output.access(26, 2) = 2.*x*y*(1. + y)*z; + output.access(26, 3) = (-1. + x)*(1. + x)*(-1. + z)*(1. + z); + output.access(26, 4) = (-1. + x*x)*(1. + 2.*y)*z; + output.access(26, 5) = (-1. + x)*(1. + x)*y*(1. + y); + } break; } case OPERATOR_D3 : { @@ -635,82 +642,84 @@ namespace Intrepid2 { output.access(19, 8) = -((-1.+ x)*x*y); output.access(19, 9) = 0.; - output.access(20, 0) = 0.; - output.access(20, 1) = -4*y*(-1.+ (z*z)); - output.access(20, 2) = -4*(-1.+ (y*y))*z; - output.access(20, 3) = -4*x*(-1.+ (z*z)); - output.access(20, 4) = -8*x*y*z; - output.access(20, 5) = -4*x*(-1.+ (y*y)); - output.access(20, 6) = 0.; - output.access(20, 7) = -4*(-1.+ (x*x))*z; - output.access(20, 8) = -4*(-1.+ (x*x))*y; - output.access(20, 9) = 0.; - - output.access(21, 0) = 0.; - output.access(21, 1) = 2.*y*(-1.+ z)*z; - output.access(21, 2) = (-1.+ (y*y))*(-1.+ 2.*z); - output.access(21, 3) = 2.*x*(-1.+ z)*z; - output.access(21, 4) = 2.*x*y*(-1.+ 2.*z); - output.access(21, 5) = 2.*x*(-1.+ (y*y)); - output.access(21, 6) = 0.; - output.access(21, 7) = (-1.+ (x*x))*(-1.+ 2.*z); - output.access(21, 8) = 2.*(-1.+ (x*x))*y; - output.access(21, 9) = 0.; - - output.access(22, 0) = 0.; - output.access(22, 1) = 2.*y*z*(1.+ z); - output.access(22, 2) = (-1.+ (y*y))*(1.+ 2.*z); - output.access(22, 3) = 2.*x*z*(1.+ z); - output.access(22, 4) = 2.*x*y*(1.+ 2.*z); - output.access(22, 5) = 2.*x*(-1.+ (y*y)); - output.access(22, 6) = 0.; - output.access(22, 7) = (-1.+ (x*x))*(1.+ 2.*z); - output.access(22, 8) = 2.*(-1.+ (x*x))*y; - output.access(22, 9) = 0.; - - output.access(23, 0) = 0.; - output.access(23, 1) = 2.*y*(-1.+ (z*z)); - output.access(23, 2) = 2.*(-1.+ (y*y))*z; - output.access(23, 3) = (-1.+ 2.*x)*(-1.+ (z*z)); - output.access(23, 4) = 2.*(-1.+ 2.*x)*y*z; - output.access(23, 5) = (-1.+ 2.*x)*(-1.+ (y*y)); - output.access(23, 6) = 0.; - output.access(23, 7) = 2.*(-1.+ x)*x*z; - output.access(23, 8) = 2.*(-1.+ x)*x*y; - output.access(23, 9) = 0.; - - output.access(24, 0) = 0.; - output.access(24, 1) = 2.*y*(-1.+ (z*z)); - output.access(24, 2) = 2.*(-1.+ (y*y))*z; - output.access(24, 3) = (1.+ 2.*x)*(-1.+ (z*z)); - output.access(24, 4) = 2.*(1.+ 2.*x)*y*z; - output.access(24, 5) = (1.+ 2.*x)*(-1.+ (y*y)); - output.access(24, 6) = 0.; - output.access(24, 7) = 2.*x*(1.+ x)*z; - output.access(24, 8) = 2.*x*(1.+ x)*y; - output.access(24, 9) = 0.; - - output.access(25, 0) = 0.; - output.access(25, 1) = (-1.+ 2.*y)*(-1.+ (z*z)); - output.access(25, 2) = 2.*(-1.+ y)*y*z; - output.access(25, 3) = 2.*x*(-1.+ (z*z)); - output.access(25, 4) = 2.*x*(-1.+ 2.*y)*z; - output.access(25, 5) = 2.*x*(-1.+ y)*y; - output.access(25, 6) = 0.; - output.access(25, 7) = 2.*(-1.+ (x*x))*z; - output.access(25, 8) = (-1.+ (x*x))*(-1.+ 2.*y); - output.access(25, 9) = 0.; - - output.access(26, 0) = 0.; - output.access(26, 1) = (1.+ 2.*y)*(-1.+ (z*z)); - output.access(26, 2) = 2.*y*(1.+ y)*z; - output.access(26, 3) = 2.*x*(-1.+ (z*z)); - output.access(26, 4) = 2.*x*(1.+ 2.*y)*z; - output.access(26, 5) = 2.*x*y*(1.+ y); - output.access(26, 6) = 0.; - output.access(26, 7) = 2.*(-1.+ (x*x))*z; - output.access(26, 8) = (-1.+ (x*x))*(1.+ 2.*y); - output.access(26, 9) = 0.; + if constexpr(!serendipity) { + output.access(20, 0) = 0.; + output.access(20, 1) = -4*y*(-1.+ (z*z)); + output.access(20, 2) = -4*(-1.+ (y*y))*z; + output.access(20, 3) = -4*x*(-1.+ (z*z)); + output.access(20, 4) = -8*x*y*z; + output.access(20, 5) = -4*x*(-1.+ (y*y)); + output.access(20, 6) = 0.; + output.access(20, 7) = -4*(-1.+ (x*x))*z; + output.access(20, 8) = -4*(-1.+ (x*x))*y; + output.access(20, 9) = 0.; + + output.access(21, 0) = 0.; + output.access(21, 1) = 2.*y*(-1.+ z)*z; + output.access(21, 2) = (-1.+ (y*y))*(-1.+ 2.*z); + output.access(21, 3) = 2.*x*(-1.+ z)*z; + output.access(21, 4) = 2.*x*y*(-1.+ 2.*z); + output.access(21, 5) = 2.*x*(-1.+ (y*y)); + output.access(21, 6) = 0.; + output.access(21, 7) = (-1.+ (x*x))*(-1.+ 2.*z); + output.access(21, 8) = 2.*(-1.+ (x*x))*y; + output.access(21, 9) = 0.; + + output.access(22, 0) = 0.; + output.access(22, 1) = 2.*y*z*(1.+ z); + output.access(22, 2) = (-1.+ (y*y))*(1.+ 2.*z); + output.access(22, 3) = 2.*x*z*(1.+ z); + output.access(22, 4) = 2.*x*y*(1.+ 2.*z); + output.access(22, 5) = 2.*x*(-1.+ (y*y)); + output.access(22, 6) = 0.; + output.access(22, 7) = (-1.+ (x*x))*(1.+ 2.*z); + output.access(22, 8) = 2.*(-1.+ (x*x))*y; + output.access(22, 9) = 0.; + + output.access(23, 0) = 0.; + output.access(23, 1) = 2.*y*(-1.+ (z*z)); + output.access(23, 2) = 2.*(-1.+ (y*y))*z; + output.access(23, 3) = (-1.+ 2.*x)*(-1.+ (z*z)); + output.access(23, 4) = 2.*(-1.+ 2.*x)*y*z; + output.access(23, 5) = (-1.+ 2.*x)*(-1.+ (y*y)); + output.access(23, 6) = 0.; + output.access(23, 7) = 2.*(-1.+ x)*x*z; + output.access(23, 8) = 2.*(-1.+ x)*x*y; + output.access(23, 9) = 0.; + + output.access(24, 0) = 0.; + output.access(24, 1) = 2.*y*(-1.+ (z*z)); + output.access(24, 2) = 2.*(-1.+ (y*y))*z; + output.access(24, 3) = (1.+ 2.*x)*(-1.+ (z*z)); + output.access(24, 4) = 2.*(1.+ 2.*x)*y*z; + output.access(24, 5) = (1.+ 2.*x)*(-1.+ (y*y)); + output.access(24, 6) = 0.; + output.access(24, 7) = 2.*x*(1.+ x)*z; + output.access(24, 8) = 2.*x*(1.+ x)*y; + output.access(24, 9) = 0.; + + output.access(25, 0) = 0.; + output.access(25, 1) = (-1.+ 2.*y)*(-1.+ (z*z)); + output.access(25, 2) = 2.*(-1.+ y)*y*z; + output.access(25, 3) = 2.*x*(-1.+ (z*z)); + output.access(25, 4) = 2.*x*(-1.+ 2.*y)*z; + output.access(25, 5) = 2.*x*(-1.+ y)*y; + output.access(25, 6) = 0.; + output.access(25, 7) = 2.*(-1.+ (x*x))*z; + output.access(25, 8) = (-1.+ (x*x))*(-1.+ 2.*y); + output.access(25, 9) = 0.; + + output.access(26, 0) = 0.; + output.access(26, 1) = (1.+ 2.*y)*(-1.+ (z*z)); + output.access(26, 2) = 2.*y*(1.+ y)*z; + output.access(26, 3) = 2.*x*(-1.+ (z*z)); + output.access(26, 4) = 2.*x*(1.+ 2.*y)*z; + output.access(26, 5) = 2.*x*y*(1.+ y); + output.access(26, 6) = 0.; + output.access(26, 7) = 2.*(-1.+ (x*x))*z; + output.access(26, 8) = (-1.+ (x*x))*(1.+ 2.*y); + output.access(26, 9) = 0.; + } break; } case OPERATOR_D4 : { @@ -866,55 +875,56 @@ namespace Intrepid2 { output.access(19, 8) = y - 2.*x*y; output.access(19, 12)= -((-1.+ x)*x); - output.access(20, 3) = 4. - 4.*z*z; - output.access(20, 4) = -8.*y*z; - output.access(20, 5) = 4. - 4.*y*y; - output.access(20, 7) = -8.*x*z; - output.access(20, 8) = -8.*x*y; - output.access(20, 12)= 4. - 4.*x*x; - - output.access(21, 3) = 2.*(-1.+ z)*z; - output.access(21, 4) = 2.*y*(-1.+ 2.*z); - output.access(21, 5) = 2.*(-1.+ y*y); - output.access(21, 7) = 2.*x*(-1.+ 2.*z); - output.access(21, 8) = 4.*x*y; - output.access(21, 12)= 2.*(-1.+ x*x); - - output.access(22, 3) = 2.*z*(1. + z); - output.access(22, 4) = 2.*(y + 2.*y*z); - output.access(22, 5) = 2.*(-1.+ y*y); - output.access(22, 7) = 2.*(x + 2.*x*z); - output.access(22, 8) = 4.*x*y; - output.access(22, 12)= 2.*(-1.+ x*x); - - output.access(23, 3) = 2.*(-1.+ z*z); - output.access(23, 4) = 4.*y*z; - output.access(23, 5) = 2.*(-1.+ y*y); - output.access(23, 7) = 2.*(-1.+ 2.*x)*z; - output.access(23, 8) = 2.*(-1.+ 2.*x)*y; - output.access(23, 12)= 2.*(-1.+ x)*x; - - output.access(24, 3) = 2.*(-1.+ z*z); - output.access(24, 4) = 4.*y*z; - output.access(24, 5) = 2.*(-1.+ y*y); - output.access(24, 7) = 2.*(z + 2.*x*z); - output.access(24, 8) = 2.*(y + 2.*x*y); - output.access(24, 12)= 2.*x*(1. + x); - - output.access(25, 3) = 2.*(-1.+ z*z); - output.access(25, 4) = 2.*(-1.+ 2.*y)*z; - output.access(25, 5) = 2.*(-1.+ y)*y; - output.access(25, 7) = 4.*x*z; - output.access(25, 8) = 2.*x*(-1.+ 2.*y); - output.access(25, 12)= 2.*(-1.+ x*x); - - output.access(26, 3) = 2.*(-1.+ z*z); - output.access(26, 4) = 2.*(z + 2.*y*z); - output.access(26, 5) = 2.*y*(1. + y); - output.access(26, 7) = 4.*x*z; - output.access(26, 8) = 2.*(x + 2.*x*y); - output.access(26, 12)= 2.*(-1.+ x*x); - + if constexpr (!serendipity) { + output.access(20, 3) = 4. - 4.*z*z; + output.access(20, 4) = -8.*y*z; + output.access(20, 5) = 4. - 4.*y*y; + output.access(20, 7) = -8.*x*z; + output.access(20, 8) = -8.*x*y; + output.access(20, 12)= 4. - 4.*x*x; + + output.access(21, 3) = 2.*(-1.+ z)*z; + output.access(21, 4) = 2.*y*(-1.+ 2.*z); + output.access(21, 5) = 2.*(-1.+ y*y); + output.access(21, 7) = 2.*x*(-1.+ 2.*z); + output.access(21, 8) = 4.*x*y; + output.access(21, 12)= 2.*(-1.+ x*x); + + output.access(22, 3) = 2.*z*(1. + z); + output.access(22, 4) = 2.*(y + 2.*y*z); + output.access(22, 5) = 2.*(-1.+ y*y); + output.access(22, 7) = 2.*(x + 2.*x*z); + output.access(22, 8) = 4.*x*y; + output.access(22, 12)= 2.*(-1.+ x*x); + + output.access(23, 3) = 2.*(-1.+ z*z); + output.access(23, 4) = 4.*y*z; + output.access(23, 5) = 2.*(-1.+ y*y); + output.access(23, 7) = 2.*(-1.+ 2.*x)*z; + output.access(23, 8) = 2.*(-1.+ 2.*x)*y; + output.access(23, 12)= 2.*(-1.+ x)*x; + + output.access(24, 3) = 2.*(-1.+ z*z); + output.access(24, 4) = 4.*y*z; + output.access(24, 5) = 2.*(-1.+ y*y); + output.access(24, 7) = 2.*(z + 2.*x*z); + output.access(24, 8) = 2.*(y + 2.*x*y); + output.access(24, 12)= 2.*x*(1. + x); + + output.access(25, 3) = 2.*(-1.+ z*z); + output.access(25, 4) = 2.*(-1.+ 2.*y)*z; + output.access(25, 5) = 2.*(-1.+ y)*y; + output.access(25, 7) = 4.*x*z; + output.access(25, 8) = 2.*x*(-1.+ 2.*y); + output.access(25, 12)= 2.*(-1.+ x*x); + + output.access(26, 3) = 2.*(-1.+ z*z); + output.access(26, 4) = 2.*(z + 2.*y*z); + output.access(26, 5) = 2.*y*(1. + y); + output.access(26, 7) = 4.*x*z; + output.access(26, 8) = 2.*(x + 2.*x*y); + output.access(26, 12)= 2.*(-1.+ x*x); + } break; } case OPERATOR_MAX : { @@ -935,17 +945,18 @@ namespace Intrepid2 { opType != OPERATOR_D3 && opType != OPERATOR_D4 && opType != OPERATOR_MAX, - ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_C2_FEM::Serial::getValues) operator is not supported"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::Serial::getValues) operator is not supported"); } } } + template template void - Basis_HGRAD_HEX_C2_FEM:: + Basis_HGRAD_HEX_DEG2_FEM:: getValues( Kokkos::DynRankView outputValues, const Kokkos::DynRankView inputPoints, const EOperator operatorType ) { @@ -972,12 +983,12 @@ namespace Intrepid2 { } case OPERATOR_CURL: INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument, - ">>> ERROR (Basis_HGRAD_HEX_C2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D"); + ">>> ERROR (Basis_HGRAD_HEX_DEG2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D"); break; case OPERATOR_DIV: INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument, - ">>> ERROR (Basis_HGRAD_HEX_C2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D"); + ">>> ERROR (Basis_HGRAD_HEX_DEG2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D"); break; case OPERATOR_D2: { @@ -998,7 +1009,7 @@ namespace Intrepid2 { case OPERATOR_D5: case OPERATOR_D6: { INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument, - ">>> ERROR (Basis_HGRAD_HEX_C2_FEM): operator not supported"); + ">>> ERROR (Basis_HGRAD_HEX_DEG2_FEM): operator not supported"); // break; commented out becuase this always throws } case OPERATOR_D7: @@ -1011,7 +1022,7 @@ namespace Intrepid2 { } default: { INTREPID2_TEST_FOR_EXCEPTION( !( Intrepid2::isValidOperator(operatorType) ), std::invalid_argument, - ">>> ERROR (Basis_HGRAD_HEX_C2_FEM): Invalid operator type"); + ">>> ERROR (Basis_HGRAD_HEX_DEG2_FEM): Invalid operator type"); } } } @@ -1020,10 +1031,10 @@ namespace Intrepid2 { // ------------------------------------------------------------------------------------- - template - Basis_HGRAD_HEX_C2_FEM:: - Basis_HGRAD_HEX_C2_FEM() { - this -> basisCardinality_ = 27; + template + Basis_HGRAD_HEX_DEG2_FEM:: + Basis_HGRAD_HEX_DEG2_FEM() { + this -> basisCardinality_ = serendipity ? 20 : 27; this -> basisDegree_ = 2; this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData >() ); this -> basisType_ = BASIS_FEM_DEFAULT; @@ -1059,6 +1070,8 @@ namespace Intrepid2 { 1, 5, 0, 1, // Node 17 -> edge 5 1, 6, 0, 1, // Node 18 -> edge 6 1, 7, 0, 1, // Node 19 -> edge 7 + + // following entries not used for serendipity elements 3, 0, 0, 1, // Node 20 -> Hexahedron 2, 4, 0, 1, // Node 21 -> face 4 2, 5, 0, 1, // Node 22 -> face 5 @@ -1069,7 +1082,7 @@ namespace Intrepid2 { }; // host tags - OrdinalTypeArray1DHost tagView(&tags[0], 108); + OrdinalTypeArray1DHost tagView(&tags[0], serendipity ? 80 : 108); // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays: this->setOrdinalTagData(this -> tagToOrdinal_, @@ -1107,14 +1120,16 @@ namespace Intrepid2 { dofCoords(18,0) = 0.0; dofCoords(18,1) = 1.0; dofCoords(18,2) = 1.0; dofCoords(19,0) = -1.0; dofCoords(19,1) = 0.0; dofCoords(19,2) = 1.0; - dofCoords(20,0) = 0.0; dofCoords(20,1) = 0.0; dofCoords(20,2) = 0.0; + if constexpr (!serendipity) { + dofCoords(20,0) = 0.0; dofCoords(20,1) = 0.0; dofCoords(20,2) = 0.0; - dofCoords(21,0) = 0.0; dofCoords(21,1) = 0.0; dofCoords(21,2) = -1.0; - dofCoords(22,0) = 0.0; dofCoords(22,1) = 0.0; dofCoords(22,2) = 1.0; - dofCoords(23,0) = -1.0; dofCoords(23,1) = 0.0; dofCoords(23,2) = 0.0; - dofCoords(24,0) = 1.0; dofCoords(24,1) = 0.0; dofCoords(24,2) = 0.0; - dofCoords(25,0) = 0.0; dofCoords(25,1) = -1.0; dofCoords(25,2) = 0.0; - dofCoords(26,0) = 0.0; dofCoords(26,1) = 1.0; dofCoords(26,2) = 0.0; + dofCoords(21,0) = 0.0; dofCoords(21,1) = 0.0; dofCoords(21,2) = -1.0; + dofCoords(22,0) = 0.0; dofCoords(22,1) = 0.0; dofCoords(22,2) = 1.0; + dofCoords(23,0) = -1.0; dofCoords(23,1) = 0.0; dofCoords(23,2) = 0.0; + dofCoords(24,0) = 1.0; dofCoords(24,1) = 0.0; dofCoords(24,2) = 0.0; + dofCoords(25,0) = 0.0; dofCoords(25,1) = -1.0; dofCoords(25,2) = 0.0; + dofCoords(26,0) = 0.0; dofCoords(26,1) = 1.0; dofCoords(26,2) = 0.0; + } this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords); Kokkos::deep_copy(this->dofCoords_, dofCoords); diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_QUAD_C2_FEM.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_QUAD_C2_FEM.hpp index 456dc961361a..8daeece82ab2 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_QUAD_C2_FEM.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_QUAD_C2_FEM.hpp @@ -54,11 +54,19 @@ namespace Intrepid2 { - /** \class Intrepid2::Basis_HGRAD_QUAD_C2_FEM + /** \class Intrepid2::Basis_HGRAD_QUAD_DEG2_FEM \brief Implementation of the default H(grad)-compatible FEM basis of degree 2 on Quadrilateral cell - Implements Lagrangian basis of degree 2 on the reference Quadrilateral cell. The basis has - cardinality 9 and spans a COMPLETE bi-quadratic polynomial space. Basis functions are dual + Implements Lagrangian basis of degree 2 on the reference Quadrilateral cell. + When the serendipity template argument is false, the basis has + cardinality 9 and spans a COMPLETE bi-quadratic polynomial space. + Note, Basis_HGRAD_QUAD_C2_FEM = Basis_HGRAD_QUAD_DEG2_FEM + + When the serendipity template argument is true, the basis has + cardinality 8 and spans an INCOMPLETE bi-quadratic polynomial space (Dofs are associated only to vertices and edges). + Note, Basis_HGRAD_QUAD_I2_Serendipity_FEM = Basis_HGRAD_QUAD_C2_FEM + + Basis functions are dual to a unisolvent set of degrees-of-freedom (DoF) defined and enumerated as follows: \verbatim @@ -94,13 +102,14 @@ namespace Intrepid2 { namespace Impl { /** - \brief See Intrepid2::Basis_HGRAD_QUAD_C2_FEM + \brief See Intrepid2::Basis_HGRAD_QUAD_DEG2_FEM */ - class Basis_HGRAD_QUAD_C2_FEM { + template + class Basis_HGRAD_QUAD_DEG2_FEM { public: - typedef struct Quadrilateral<9> cell_topology_type; + typedef struct Quadrilateral cell_topology_type; /** - \brief See Intrepid2::Basis_HGRAD_QUAD_C2_FEM + \brief See Intrepid2::Basis_HGRAD_QUAD_DEG2_FEM */ template struct Serial { @@ -122,7 +131,7 @@ namespace Intrepid2 { const EOperator operatorType); /** - \brief See Intrepid2::Basis_HGRAD_QUAD_C2_FEM + \brief See Intrepid2::Basis_HGRAD_QUAD_DEG2_FEM */ template>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_I1_FEM::Serial::getValues) operator is not supported"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_DEG2_FEM::Serial::getValues) operator is not supported"); } } } @@ -174,10 +183,11 @@ namespace Intrepid2 { }; } - template - class Basis_HGRAD_QUAD_C2_FEM : public Basis { + template + class Basis_HGRAD_QUAD_DEG2_FEM : public Basis { public: using OrdinalTypeArray1DHost = typename Basis::OrdinalTypeArray1DHost; using OrdinalTypeArray2DHost = typename Basis::OrdinalTypeArray2DHost; @@ -185,7 +195,7 @@ namespace Intrepid2 { /** \brief Constructor. */ - Basis_HGRAD_QUAD_C2_FEM(); + Basis_HGRAD_QUAD_DEG2_FEM(); using OutputViewType = typename Basis::OutputViewType; using PointViewType = typename Basis::PointViewType; @@ -206,10 +216,16 @@ namespace Intrepid2 { this->getBaseCellTopology(), this->getCardinality() ); #endif - Impl::Basis_HGRAD_QUAD_C2_FEM:: - getValues( outputValues, - inputPoints, - operatorType ); + if constexpr (serendipity) + Impl::Basis_HGRAD_QUAD_DEG2_FEM:: + getValues( outputValues, + inputPoints, + operatorType ); + else + Impl::Basis_HGRAD_QUAD_DEG2_FEM:: + getValues( outputValues, + inputPoints, + operatorType ); } virtual @@ -218,13 +234,13 @@ namespace Intrepid2 { #ifdef HAVE_INTREPID2_DEBUG // Verify rank of output array. INTREPID2_TEST_FOR_EXCEPTION( dofCoords.rank() != 2, std::invalid_argument, - ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_C2_FEM::getDofCoords) rank = 2 required for dofCoords array"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_DEG2_FEM::getDofCoords) rank = 2 required for dofCoords array"); // Verify 0th dimension of output array. INTREPID2_TEST_FOR_EXCEPTION( static_cast(dofCoords.extent(0)) != this->getCardinality(), std::invalid_argument, - ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_C2_FEM::getDofCoords) mismatch in number of dof and 0th dimension of dofCoords array"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_DEG2_FEM::getDofCoords) mismatch in number of dof and 0th dimension of dofCoords array"); // Verify 1st dimension of output array. INTREPID2_TEST_FOR_EXCEPTION( dofCoords.extent(1) != this->getBaseCellTopology().getDimension(), std::invalid_argument, - ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_C2_FEM::getDofCoords) incorrect reference cell (1st) dimension in dofCoords array"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_DEG2_FEM::getDofCoords) incorrect reference cell (1st) dimension in dofCoords array"); #endif Kokkos::deep_copy(dofCoords, this->dofCoords_); } @@ -235,10 +251,10 @@ namespace Intrepid2 { #ifdef HAVE_INTREPID2_DEBUG // Verify rank of output array. INTREPID2_TEST_FOR_EXCEPTION( dofCoeffs.rank() != 1, std::invalid_argument, - ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_C2_FEM::getdofCoeffs) rank = 1 required for dofCoeffs array"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_DEG2_FEM::getdofCoeffs) rank = 1 required for dofCoeffs array"); // Verify 0th dimension of output array. INTREPID2_TEST_FOR_EXCEPTION( static_cast(dofCoeffs.extent(0)) != this->getCardinality(), std::invalid_argument, - ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_C2_FEM::getdofCoeffs) mismatch in number of dof and 0th dimension of dofCoeffs array"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_DEG2_FEM::getdofCoeffs) mismatch in number of dof and 0th dimension of dofCoeffs array"); #endif Kokkos::deep_copy(dofCoeffs, 1.0); } @@ -246,7 +262,10 @@ namespace Intrepid2 { virtual const char* getName() const override { - return "Intrepid2_HGRAD_QUAD_C2_FEM"; + if constexpr (serendipity) + return "Intrepid2_HGRAD_QUAD_I2_Serendipity_FEM"; + else + return "Intrepid2_HGRAD_QUAD_C2_FEM"; } /** \brief returns the basis associated to a subCell. @@ -268,9 +287,16 @@ namespace Intrepid2 { BasisPtr getHostBasis() const override{ - return Teuchos::rcp(new Basis_HGRAD_QUAD_C2_FEM()); + return Teuchos::rcp(new Basis_HGRAD_QUAD_DEG2_FEM()); } }; + + template + using Basis_HGRAD_QUAD_C2_FEM = Basis_HGRAD_QUAD_DEG2_FEM; + + template + using Basis_HGRAD_QUAD_I2_Serendipity_FEM = Basis_HGRAD_QUAD_DEG2_FEM; + }// namespace Intrepid2 #include "Intrepid2_HGRAD_QUAD_C2_FEMDef.hpp" diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_QUAD_C2_FEMDef.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_QUAD_C2_FEMDef.hpp index 9b07b309be08..7fca6fce62e9 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_QUAD_C2_FEMDef.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_QUAD_C2_FEMDef.hpp @@ -54,13 +54,14 @@ namespace Intrepid2 { // ------------------------------------------------------------------------------------- namespace Impl { - + + template template template KOKKOS_INLINE_FUNCTION void - Basis_HGRAD_QUAD_C2_FEM::Serial:: + Basis_HGRAD_QUAD_DEG2_FEM::Serial:: getValues( OutputViewType output, const inputViewType input ) { switch (opType) { @@ -78,8 +79,11 @@ namespace Intrepid2 { output.access(5) = x*(x + 1.0)*(1.0 - y)*(1.0 + y)/2.0; output.access(6) = (1.0 - x)*(1.0 + x)*y*(y + 1.0)/2.0; output.access(7) = x*(x - 1.0)*(1.0 - y)*(1.0 + y)/2.0; + // quad bubble basis function - output.access(8) = (1.0 - x)*(1.0 + x)*(1.0 - y)*(1.0 + y); + if constexpr (!serendipity) { + output.access(8) = (1.0 - x)*(1.0 + x)*(1.0 - y)*(1.0 + y); + } break; } case OPERATOR_D1 : @@ -112,8 +116,10 @@ namespace Intrepid2 { output.access(7, 0) = 0.5*(1.0 - y)*(1.0+ y)*(-1.0 + 2.0*x); output.access(7, 1) = (1.0 - x)*x*y; - output.access(8, 0) =-2.0*(1.0 - y)*(1.0 + y)*x; - output.access(8, 1) =-2.0*(1.0 - x)*(1.0 + x)*y; + if constexpr (!serendipity) { + output.access(8, 0) =-2.0*(1.0 - y)*(1.0 + y)*x; + output.access(8, 1) =-2.0*(1.0 - x)*(1.0 + x)*y; + } break; } case OPERATOR_CURL : { @@ -145,9 +151,11 @@ namespace Intrepid2 { output.access(7, 1) =-0.5*(1.0 - y)*(1.0 + y)*(-1.0 + 2.0*x); output.access(7, 0) = (1.0 - x)*x*y; - - output.access(8, 1) = 2.0*(1.0 - y)*(1.0 + y)*x; - output.access(8, 0) =-2.0*(1.0 - x)*(1.0 + x)*y; + + if constexpr (!serendipity) { + output.access(8, 1) = 2.0*(1.0 - y)*(1.0 + y)*x; + output.access(8, 0) =-2.0*(1.0 - x)*(1.0 + x)*y; + } break; } case OPERATOR_D2 : { @@ -186,9 +194,11 @@ namespace Intrepid2 { output.access(7, 1) = x*(0. - 2.*y) + 1.*y; output.access(7, 2) = (1.0 - x)*x; - output.access(8, 0) =-2.0 + 2.0*y*y; - output.access(8, 1) = 4*x*y; - output.access(8, 2) =-2.0 + 2.0*x*x; + if constexpr (!serendipity) { + output.access(8, 0) =-2.0 + 2.0*y*y; + output.access(8, 1) = 4*x*y; + output.access(8, 2) =-2.0 + 2.0*x*x; + } break; } case OPERATOR_D3 : { @@ -234,10 +244,12 @@ namespace Intrepid2 { output.access(7, 2) = 1.0 - 2.0*x; output.access(7, 3) = 0.0; - output.access(8, 0) = 0.0; - output.access(8, 1) = 4.0*y; - output.access(8, 2) = 4.0*x; - output.access(8, 3) = 0.0; + if constexpr (!serendipity) { + output.access(8, 0) = 0.0; + output.access(8, 1) = 4.0*y; + output.access(8, 2) = 4.0*x; + output.access(8, 3) = 0.0; + } break; } case OPERATOR_D4 : { @@ -289,11 +301,13 @@ namespace Intrepid2 { output.access(7, 3) = 0.0; output.access(7, 4) = 0.0; - output.access(8, 0) = 0.0; - output.access(8, 1) = 0.0; - output.access(8, 2) = 4.0; - output.access(8, 3) = 0.0; - output.access(8, 4) = 0.0; + if constexpr (!serendipity) { + output.access(8, 0) = 0.0; + output.access(8, 1) = 0.0; + output.access(8, 2) = 4.0; + output.access(8, 3) = 0.0; + output.access(8, 4) = 0.0; + } break; } case OPERATOR_MAX : { @@ -319,12 +333,13 @@ namespace Intrepid2 { } } } - + + template template void - Basis_HGRAD_QUAD_C2_FEM:: + Basis_HGRAD_QUAD_DEG2_FEM:: getValues( Kokkos::DynRankView outputValues, const Kokkos::DynRankView inputPoints, const EOperator operatorType ) { @@ -393,17 +408,14 @@ namespace Intrepid2 { } } - - - } // ------------------------------------------------------------------------------------- - template - Basis_HGRAD_QUAD_C2_FEM:: - Basis_HGRAD_QUAD_C2_FEM() { - this->basisCardinality_ = 9; + template + Basis_HGRAD_QUAD_DEG2_FEM:: + Basis_HGRAD_QUAD_DEG2_FEM() { + this->basisCardinality_ = serendipity ? 8 : 9; this->basisDegree_ = 2; this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData >() ); this->basisType_ = BASIS_FEM_DEFAULT; @@ -427,11 +439,11 @@ namespace Intrepid2 { 1, 1, 0, 1, 1, 2, 0, 1, 1, 3, 0, 1, - // quad center + // quad center, not used for serendipity elements 2, 0, 0, 1}; //host view - OrdinalTypeArray1DHost tagView(&tags[0],36); + OrdinalTypeArray1DHost tagView(&tags[0], serendipity ? 32 : 36); // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays: this->setOrdinalTagData(this->tagToOrdinal_, @@ -458,7 +470,9 @@ namespace Intrepid2 { dofCoords(6,0) = 0.0; dofCoords(6,1) = 1.0; dofCoords(7,0) = -1.0; dofCoords(7,1) = 0.0; - dofCoords(8,0) = 0.0; dofCoords(8,1) = 0.0; + if constexpr (!serendipity) { + dofCoords(8,0) = 0.0; dofCoords(8,1) = 0.0; + } this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords); Kokkos::deep_copy(this->dofCoords_, dofCoords); diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_WEDGE_C1_FEM.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_WEDGE_C1_FEM.hpp index fa13dc0df708..021619f02206 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_WEDGE_C1_FEM.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_WEDGE_C1_FEM.hpp @@ -50,6 +50,8 @@ #define __INTREPID2_HGRAD_WEDGE_C1_FEM_HPP__ #include "Intrepid2_Basis.hpp" +#include "Intrepid2_HGRAD_QUAD_C1_FEM.hpp" +#include "Intrepid2_HGRAD_TRI_C1_FEM.hpp" namespace Intrepid2 { @@ -234,6 +236,28 @@ namespace Intrepid2 { return "Intrepid2_HGRAD_WEDGE_C1_FEM"; } + /** \brief returns the basis associated to a subCell. + + The bases of the subCell are the restriction to the subCell + of the bases of the parent cell. + \param [in] subCellDim - dimension of subCell + \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{ + if(subCellDim == 1) { + return Teuchos::rcp(new + Basis_HGRAD_LINE_C1_FEM()); + } else if(subCellDim == 2) { + if(subCellOrd < 3) //lateral faces + return Teuchos::rcp(new Basis_HGRAD_QUAD_C1_FEM()); + else + return Teuchos::rcp(new Basis_HGRAD_TRI_C1_FEM()); + } + INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds"); + } + BasisPtr getHostBasis() const override{ return Teuchos::rcp(new Basis_HGRAD_WEDGE_C1_FEM()); diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_WEDGE_C2_FEM.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_WEDGE_C2_FEM.hpp index e5c66768ede5..c2978d1c257a 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_WEDGE_C2_FEM.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_WEDGE_C2_FEM.hpp @@ -50,15 +50,25 @@ #define __INTREPID2_HGRAD_WEDGE_C2_FEM_HPP__ #include "Intrepid2_Basis.hpp" +#include "Intrepid2_HGRAD_QUAD_C2_FEM.hpp" +#include "Intrepid2_HGRAD_TRI_C2_FEM.hpp" namespace Intrepid2 { - /** \class Intrepid2::Basis_HGRAD_WEDGE_C2_FEM + /** \class Intrepid2::Basis_HGRAD_WEDGE_DEG2_FEM \brief Implementation of the default H(grad)-compatible FEM basis of degree 2 on Wedge cell - Implements Lagrangian basis of degree 2 on the reference Wedge cell. The basis has - cardinality 18 and spans a COMPLETE bi-quadratic polynomial space. Basis functions are dual - to a unisolvent set of degrees-of-freedom (DoF) defined and enumerated as follows: + Implements Lagrangian basis of degree 2 on the reference Wedge cell. + + When the serendipity template argument is false, the basis has + cardinality 18 and spans a COMPLETE bi-quadratic polynomial space. + Note, Basis_HGRAD_WEDGE_C2_FEM = Basis_HGRAD_WEDGE_DEG2_FEM + + When the serendipity template argument is true, the basis has + cardinality 15 and spans an ICOMPLETE bi-quadratic polynomial space. + Note, Basis_HGRAD_WEDGE_I2_Serendipity_FEM = Basis_HGRAD_WEDGE_DEG2_FEM + + Basis functions are dual to a unisolvent set of degrees-of-freedom (DoF) defined and enumerated as follows: \verbatim ================================================================================================= @@ -112,21 +122,22 @@ namespace Intrepid2 { order in this topology does not follow the natural oder of k-subcells where the nodes are located, except for nodes 0 to 5 which coincide with the vertices of the base Wedge<6> topology. As a result, L_0 to L_5 are associated with nodes 0 to 5, but - L_6 to L_14 are not associated with edges 0 to 9 in that order. The last three nodes - are located on 2-subcells (faces) and follow their order. Thus, L_15, L_16 and L17 are - associated with faces 0, 1 and 2 in that order. + L_6 to L_14 are not associated with edges 0 to 9 in that order. The last three nodes, + for the non sereendipity basis, are located on 2-subcells (faces) and follow their order. + Thus, L_15, L_16 and L17 are associated with faces 0, 1 and 2 in that order. */ namespace Impl { /** - \brief See Intrepid2::Basis_HGRAD_WEDGE_C2_FEM + \brief See Intrepid2::Basis_HGRAD_WEDGE_DEG2_FEM */ - class Basis_HGRAD_WEDGE_C2_FEM { + template + class Basis_HGRAD_WEDGE_DEG2_FEM { public: - typedef struct Wedge<18> cell_topology_type; + typedef struct Wedge cell_topology_type; /** - \brief See Intrepid2::Basis_HGRAD_WEDGE_C2_FEM + \brief See Intrepid2::Basis_HGRAD_WEDGE_DEG2_FEM */ template struct Serial { @@ -148,7 +159,7 @@ namespace Intrepid2 { const EOperator operatorType); /** - \brief See Intrepid2::Basis_HGRAD_WEDGE_C2_FEM + \brief See Intrepid2::Basis_HGRAD_WEDGE_DEG2_FEM */ template>> ERROR: (Intrepid2::Basis_HGRAD_WEDGE_C2_FEM::Serial::getValues) operator is not supported"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_WEDGE_DEG2_FEM::Serial::getValues) operator is not supported"); } } } @@ -193,17 +204,18 @@ namespace Intrepid2 { }; } - template - class Basis_HGRAD_WEDGE_C2_FEM : public Basis { + class Basis_HGRAD_WEDGE_DEG2_FEM : public Basis { public: using OrdinalTypeArray1DHost = typename Basis::OrdinalTypeArray1DHost; using OrdinalTypeArray2DHost = typename Basis::OrdinalTypeArray2DHost; using OrdinalTypeArray3DHost = typename Basis::OrdinalTypeArray3DHost; /** \brief Constructor. */ - Basis_HGRAD_WEDGE_C2_FEM(); + Basis_HGRAD_WEDGE_DEG2_FEM(); using OutputViewType = typename Basis::OutputViewType; using PointViewType = typename Basis::PointViewType; @@ -224,10 +236,16 @@ namespace Intrepid2 { this->getBaseCellTopology(), this->getCardinality() ); #endif - Impl::Basis_HGRAD_WEDGE_C2_FEM:: - getValues( outputValues, - inputPoints, - operatorType ); + if constexpr (serendipity) + Impl::Basis_HGRAD_WEDGE_DEG2_FEM:: + getValues( outputValues, + inputPoints, + operatorType ); + else + Impl::Basis_HGRAD_WEDGE_DEG2_FEM:: + getValues( outputValues, + inputPoints, + operatorType );; } virtual @@ -236,13 +254,13 @@ namespace Intrepid2 { #ifdef HAVE_INTREPID2_DEBUG // Verify rank of output array. INTREPID2_TEST_FOR_EXCEPTION( dofCoords.rank() != 2, std::invalid_argument, - ">>> ERROR: (Intrepid2::Basis_HGRAD_WEDGE_C2_FEM::getDofCoords) rank = 2 required for dofCoords array"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_WEDGE_DEG2_FEM::getDofCoords) rank = 2 required for dofCoords array"); // Verify 0th dimension of output array. INTREPID2_TEST_FOR_EXCEPTION( static_cast(dofCoords.extent(0)) != this->getCardinality(), std::invalid_argument, - ">>> ERROR: (Intrepid2::Basis_HGRAD_WEDGE_C2_FEM::getDofCoords) mismatch in number of dof and 0th dimension of dofCoords array"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_WEDGE_DEG2_FEM::getDofCoords) mismatch in number of dof and 0th dimension of dofCoords array"); // Verify 1st dimension of output array. INTREPID2_TEST_FOR_EXCEPTION( dofCoords.extent(1) != this->getBaseCellTopology().getDimension(), std::invalid_argument, - ">>> ERROR: (Intrepid2::Basis_HGRAD_WEDGE_C2_FEM::getDofCoords) incorrect reference cell (1st) dimension in dofCoords array"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_WEDGE_DEG2_FEM::getDofCoords) incorrect reference cell (1st) dimension in dofCoords array"); #endif Kokkos::deep_copy(dofCoords, this->dofCoords_); } @@ -253,10 +271,10 @@ namespace Intrepid2 { #ifdef HAVE_INTREPID2_DEBUG // Verify rank of output array. INTREPID2_TEST_FOR_EXCEPTION( dofCoeffs.rank() != 1, std::invalid_argument, - ">>> ERROR: (Intrepid2::Basis_HGRAD_WEDGE_C2_FEM::getdofCoeffs) rank = 1 required for dofCoeffs array"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_WEDGE_DEG2_FEM::getdofCoeffs) rank = 1 required for dofCoeffs array"); // Verify 0th dimension of output array. INTREPID2_TEST_FOR_EXCEPTION( static_cast(dofCoeffs.extent(0)) != this->getCardinality(), std::invalid_argument, - ">>> ERROR: (Intrepid2::Basis_HGRAD_WEDGE_C2_FEM::getdofCoeffs) mismatch in number of dof and 0th dimension of dofCoeffs array"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_WEDGE_DEG2_FEM::getdofCoeffs) mismatch in number of dof and 0th dimension of dofCoeffs array"); #endif Kokkos::deep_copy(dofCoeffs, 1.0); } @@ -264,15 +282,46 @@ namespace Intrepid2 { virtual const char* getName() const override { - return "Intrepid2_HGRAD_WEDGE_C2_FEM"; + if constexpr (serendipity) + return "Intrepid2_HGRAD_WEDGE_I2_Serendipity_FEM"; + else + return "Intrepid2_HGRAD_WEDGE_C2_FEM"; + } + + /** \brief returns the basis associated to a subCell. + + The bases of the subCell are the restriction to the subCell + of the bases of the parent cell. + \param [in] subCellDim - dimension of subCell + \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{ + if(subCellDim == 1) { + return Teuchos::rcp(new + Basis_HGRAD_LINE_C2_FEM()); + } else if(subCellDim == 2) { + if(subCellOrd < 3) //lateral faces + return Teuchos::rcp(new Basis_HGRAD_QUAD_DEG2_FEM()); + else + return Teuchos::rcp(new Basis_HGRAD_TRI_C2_FEM()); + } + INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds"); } BasisPtr getHostBasis() const override{ - return Teuchos::rcp(new Basis_HGRAD_WEDGE_C2_FEM()); + return Teuchos::rcp(new Basis_HGRAD_WEDGE_DEG2_FEM()); } }; + template + using Basis_HGRAD_WEDGE_C2_FEM = Basis_HGRAD_WEDGE_DEG2_FEM; + + template + using Basis_HGRAD_WEDGE_I2_Serendipity_FEM = Basis_HGRAD_WEDGE_DEG2_FEM; + }// namespace Intrepid2 #include "Intrepid2_HGRAD_WEDGE_C2_FEMDef.hpp" diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_WEDGE_C2_FEMDef.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_WEDGE_C2_FEMDef.hpp index 27869c9ec7c4..1181da1ad5b2 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_WEDGE_C2_FEMDef.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_WEDGE_C2_FEMDef.hpp @@ -53,13 +53,14 @@ namespace Intrepid2 { // ------------------------------------------------------------------------------------- namespace Impl { - + + template template template KOKKOS_INLINE_FUNCTION void - Basis_HGRAD_WEDGE_C2_FEM::Serial:: + Basis_HGRAD_WEDGE_DEG2_FEM::Serial:: getValues( OutputViewType output, const inputViewType input ) { switch (opType) { @@ -85,9 +86,11 @@ namespace Intrepid2 { output.access(12) = -2.*x*(-1. + x + y)*z*(1. + z); output.access(13) = 2.*x*y*z*(1. + z); output.access(14) = -2.*y*(-1. + x + y)*z*(1. + z); - output.access(15) = 4.*x*(-1. + x + y)*(-1. + z)*(1. + z); - output.access(16) = -4.*x*y*(-1. + z)*(1. + z); - output.access(17) = 4.*y*(-1. + x + y)*(-1. + z)*(1. + z); + if constexpr (!serendipity) { + output.access(15) = 4.*x*(-1. + x + y)*(-1. + z)*(1. + z); + output.access(16) = -4.*x*y*(-1. + z)*(1. + z); + output.access(17) = 4.*y*(-1. + x + y)*(-1. + z)*(1. + z); + } break; } case OPERATOR_GRAD: { @@ -156,17 +159,19 @@ namespace Intrepid2 { output.access(14, 1) = -2*(-1 + x + 2*y)*z*(1 + z); output.access(14, 2) = -2*y*(-1 + x + y)*(1 + 2*z); - output.access(15, 0) = 4*(-1 + 2*x + y)*(-1 + z*z); - output.access(15, 1) = 4*x*(-1 + z)*(1 + z); - output.access(15, 2) = 8*x*(-1 + x + y)*z; + if constexpr (!serendipity) { + output.access(15, 0) = 4*(-1 + 2*x + y)*(-1 + z*z); + output.access(15, 1) = 4*x*(-1 + z)*(1 + z); + output.access(15, 2) = 8*x*(-1 + x + y)*z; - output.access(16, 0) = -4*y*(-1 + z)*(1 + z); - output.access(16, 1) = -4*x*(-1 + z)*(1 + z); - output.access(16, 2) = -8*x*y*z; + output.access(16, 0) = -4*y*(-1 + z)*(1 + z); + output.access(16, 1) = -4*x*(-1 + z)*(1 + z); + output.access(16, 2) = -8*x*y*z; - output.access(17, 0) = 4*y*(-1 + z)*(1 + z); - output.access(17, 1) = 4*(-1 + x + 2*y)*(-1 + z*z); - output.access(17, 2) = 8*y*(-1 + x + y)*z; + output.access(17, 0) = 4*y*(-1 + z)*(1 + z); + output.access(17, 1) = 4*(-1 + x + 2*y)*(-1 + z*z); + output.access(17, 2) = 8*y*(-1 + x + y)*z; + } break; } case OPERATOR_D2: { @@ -279,27 +284,29 @@ namespace Intrepid2 { output.access(14, 4) = -2.*(-1. + x + 2.*y)*(1. + 2.*z); output.access(14, 5) = -4.*y*(-1. + x + y); - output.access(15, 0) = 8.*(-1. + z*z); - output.access(15, 1) = 4.*(-1. + z*z); - output.access(15, 2) = 8.*(-1. + 2.*x + y)*z; - output.access(15, 3) = 0.; - output.access(15, 4) = 8.*x*z; - output.access(15, 5) = 8.*x*(-1. + x + y); - - output.access(16, 0) = 0.; - output.access(16, 1) = 4. - 4.*z*z; - output.access(16, 2) = -8.*y*z; - output.access(16, 3) = 0.; - output.access(16, 4) = -8.*x*z; - output.access(16, 5) = -8.*x*y; - - - output.access(17, 0) = 0.; - output.access(17, 1) = 4.*(-1. + z*z); - output.access(17, 2) = 8.*y*z; - output.access(17, 3) = 8.*(-1. + z*z); - output.access(17, 4) = 8.*(-1. + x + 2.*y)*z; - output.access(17, 5) = 8.*y*(-1. + x + y); + if constexpr (!serendipity) { + output.access(15, 0) = 8.*(-1. + z*z); + output.access(15, 1) = 4.*(-1. + z*z); + output.access(15, 2) = 8.*(-1. + 2.*x + y)*z; + output.access(15, 3) = 0.; + output.access(15, 4) = 8.*x*z; + output.access(15, 5) = 8.*x*(-1. + x + y); + + output.access(16, 0) = 0.; + output.access(16, 1) = 4. - 4.*z*z; + output.access(16, 2) = -8.*y*z; + output.access(16, 3) = 0.; + output.access(16, 4) = -8.*x*z; + output.access(16, 5) = -8.*x*y; + + + output.access(17, 0) = 0.; + output.access(17, 1) = 4.*(-1. + z*z); + output.access(17, 2) = 8.*y*z; + output.access(17, 3) = 8.*(-1. + z*z); + output.access(17, 4) = 8.*(-1. + x + 2.*y)*z; + output.access(17, 5) = 8.*y*(-1. + x + y); + } break; } case OPERATOR_D3: { @@ -472,38 +479,40 @@ namespace Intrepid2 { output.access(14, 8) = -4.*(-1 + x + 2*y); output.access(14, 9) = 0.; - output.access(15, 0) = 0.; - output.access(15, 1) = 0.; - output.access(15, 2) = 16.*z; - output.access(15, 3) = 0.; - output.access(15, 4) = 8.*z; - output.access(15, 5) = 8.*(-1 + 2*x + y); - output.access(15, 6) = 0.; - output.access(15, 7) = 0.; - output.access(15, 8) = 8.*x; - output.access(15, 9) = 0.; - - output.access(16, 0) = 0.; - output.access(16, 1) = 0.; - output.access(16, 2) = 0.; - output.access(16, 3) = 0.; - output.access(16, 4) = -8.*z; - output.access(16, 5) = -8.*y; - output.access(16, 6) = 0.; - output.access(16, 7) = 0.; - output.access(16, 8) = -8.*x; - output.access(16, 9) = 0.; - - output.access(17, 0) = 0.; - output.access(17, 1) = 0.; - output.access(17, 2) = 0.; - output.access(17, 3) = 0.; - output.access(17, 4) = 8.*z; - output.access(17, 5) = 8.*y; - output.access(17, 6) = 0.; - output.access(17, 7) = 16.*z; - output.access(17, 8) = 8.*(-1 + x + 2*y); - output.access(17, 9) = 0.; + if constexpr (!serendipity) { + output.access(15, 0) = 0.; + output.access(15, 1) = 0.; + output.access(15, 2) = 16.*z; + output.access(15, 3) = 0.; + output.access(15, 4) = 8.*z; + output.access(15, 5) = 8.*(-1 + 2*x + y); + output.access(15, 6) = 0.; + output.access(15, 7) = 0.; + output.access(15, 8) = 8.*x; + output.access(15, 9) = 0.; + + output.access(16, 0) = 0.; + output.access(16, 1) = 0.; + output.access(16, 2) = 0.; + output.access(16, 3) = 0.; + output.access(16, 4) = -8.*z; + output.access(16, 5) = -8.*y; + output.access(16, 6) = 0.; + output.access(16, 7) = 0.; + output.access(16, 8) = -8.*x; + output.access(16, 9) = 0.; + + output.access(17, 0) = 0.; + output.access(17, 1) = 0.; + output.access(17, 2) = 0.; + output.access(17, 3) = 0.; + output.access(17, 4) = 8.*z; + output.access(17, 5) = 8.*y; + output.access(17, 6) = 0.; + output.access(17, 7) = 16.*z; + output.access(17, 8) = 8.*(-1 + x + 2*y); + output.access(17, 9) = 0.; + } break; } case OPERATOR_D4: { @@ -554,14 +563,16 @@ namespace Intrepid2 { output.access(14, 8) =-4; output.access(14,12) =-8.; - output.access(15, 5) =16.; - output.access(15, 8) = 8.; + if constexpr (!serendipity) { + output.access(15, 5) =16.; + output.access(15, 8) = 8.; - output.access(16, 8) =-8.; + output.access(16, 8) =-8.; - output.access(17, 8) = 8.; - output.access(17,12) =16.; + output.access(17, 8) = 8.; + output.access(17,12) =16.; + } break; } case OPERATOR_MAX : { @@ -581,17 +592,17 @@ namespace Intrepid2 { opType != OPERATOR_D3 && opType != OPERATOR_D4 && opType != OPERATOR_MAX, - ">>> ERROR: (Intrepid2::Basis_HGRAD_WEDGE_C2_FEM::Serial::getValues) operator is not supported"); + ">>> ERROR: (Intrepid2::Basis_HGRAD_WEDGE_DEG2_FEM::Serial::getValues) operator is not supported"); } } } - - template + template void - Basis_HGRAD_WEDGE_C2_FEM:: + Basis_HGRAD_WEDGE_DEG2_FEM:: getValues( Kokkos::DynRankView outputValues, const Kokkos::DynRankView inputPoints, const EOperator operatorType ) { @@ -602,6 +613,7 @@ namespace Intrepid2 { // Number of evaluation points = dim 0 of inputPoints const auto loopSize = inputPoints.extent(0); Kokkos::RangePolicy > policy(0, loopSize); + switch (operatorType) { case OPERATOR_VALUE: { @@ -617,12 +629,12 @@ namespace Intrepid2 { } case OPERATOR_CURL: { INTREPID2_TEST_FOR_EXCEPTION( operatorType == OPERATOR_CURL, std::invalid_argument, - ">>> ERROR (Basis_HGRAD_WEDGE_C2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D"); + ">>> ERROR (Basis_HGRAD_WEDGE_DEG2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D"); break; } case OPERATOR_DIV: { INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument, - ">>> ERROR (Basis_HGRAD_WEDGE_C2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D"); + ">>> ERROR (Basis_HGRAD_WEDGE_DEG2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D"); break; } case OPERATOR_D2: { @@ -652,18 +664,19 @@ namespace Intrepid2 { } default: { INTREPID2_TEST_FOR_EXCEPTION( !( Intrepid2::isValidOperator(operatorType) ), std::invalid_argument, - ">>> ERROR (Basis_HGRAD_WEDGE_C2_FEM): Invalid operator type"); + ">>> ERROR (Basis_HGRAD_WEDGE_DEG2_FEM): Invalid operator type"); } } } } + // ------------------------------------------------------------------------------------- - template - Basis_HGRAD_WEDGE_C2_FEM:: - Basis_HGRAD_WEDGE_C2_FEM() { - this->basisCardinality_ = 18; + template + Basis_HGRAD_WEDGE_DEG2_FEM:: + Basis_HGRAD_WEDGE_DEG2_FEM() { + this->basisCardinality_ = serendipity ? 15 : 18; this->basisDegree_ = 2; this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData >() ); this->basisType_ = BASIS_FEM_DEFAULT; @@ -694,17 +707,16 @@ namespace Intrepid2 { 1, 3, 0, 1, 1, 4, 0, 1, 1, 5, 0, 1, + // following entries not used for serendipity elements 2, 0, 0, 1, 2, 1, 0, 1, 2, 2, 0, 1 }; // host tags - OrdinalTypeArray1DHost tagView(&tags[0], 72); + OrdinalTypeArray1DHost tagView(&tags[0], serendipity ? 60 : 72); // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays: - //OrdinalTypeArray2DHost ordinalToTag; - //OrdinalTypeArray3DHost tagToOrdinal; this->setOrdinalTagData(this->tagToOrdinal_, this->ordinalToTag_, tagView, @@ -736,9 +748,12 @@ namespace Intrepid2 { dofCoords(12,0)= 0.5; dofCoords(12,1)= 0.0; dofCoords(12,2)= 1.0; dofCoords(13,0)= 0.5; dofCoords(13,1)= 0.5; dofCoords(13,2)= 1.0; dofCoords(14,0)= 0.0; dofCoords(14,1)= 0.5; dofCoords(14,2)= 1.0; - dofCoords(15,0)= 0.5; dofCoords(15,1)= 0.0; dofCoords(15,2)= 0.0; - dofCoords(16,0)= 0.5; dofCoords(16,1)= 0.5; dofCoords(16,2)= 0.0; - dofCoords(17,0)= 0.0; dofCoords(17,1)= 0.5; dofCoords(17,2)= 0.0; + + if constexpr (!serendipity) { + dofCoords(15,0)= 0.5; dofCoords(15,1)= 0.0; dofCoords(15,2)= 0.0; + dofCoords(16,0)= 0.5; dofCoords(16,1)= 0.5; dofCoords(16,2)= 0.0; + dofCoords(17,0)= 0.0; dofCoords(17,1)= 0.5; dofCoords(17,2)= 0.0; + } this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords); Kokkos::deep_copy(this->dofCoords_, dofCoords); diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_TensorBasis.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_TensorBasis.hpp index df81a5ed1975..0430c1cbe493 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_TensorBasis.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_TensorBasis.hpp @@ -1397,7 +1397,9 @@ struct OperatorTensorDecomposition /** \brief Fills in coefficients of degrees of freedom on the reference cell \param [out] dofCoeffs - the container into which to place the degrees of freedom. - dofCoeffs is a rank 1 with dimension equal to the cardinality of the basis. + dofCoeffs has the rank 1 or 2 and the first dimension equals the cardinality of the basis. + dofCoeffs has rank 2 if either basis1 or basis2 is a vector basis, in which case the second dimension equals the number of vector components + We do not support the case where both basis1 and basis2 are vector bases Note that getDofCoeffs() is not supported by all bases; in particular, hierarchical bases do not generally support this. */ @@ -1407,23 +1409,35 @@ struct OperatorTensorDecomposition using ResultLayout = typename DeduceLayout< typename BasisBase::ScalarViewType >::result_layout; using ViewType = Kokkos::DynRankView; - ViewType dofCoeffs1("dofCoeffs1",basis1_->getCardinality()); - ViewType dofCoeffs2("dofCoeffs2",basis2_->getCardinality()); + const ordinal_type basisCardinality1 = basis1_->getCardinality(); + const ordinal_type basisCardinality2 = basis2_->getCardinality(); + + bool isVectorBasis1 = getFieldRank(basis1_->getFunctionSpace()) == 1; + bool isVectorBasis2 = getFieldRank(basis2_->getFunctionSpace()) == 1; + + INTREPID2_TEST_FOR_EXCEPTION(isVectorBasis1 && isVectorBasis2, std::invalid_argument, "the case in which basis1 and basis2 are vector bases is not supported"); + int basisDim1 = isVectorBasis1 ? basis1_->getBaseCellTopology().getDimension() : 1; + int basisDim2 = isVectorBasis2 ? basis2_->getBaseCellTopology().getDimension() : 1; + + auto dofCoeffs1 = isVectorBasis1 ? ViewType("dofCoeffs1",basis1_->getCardinality(), basisDim1) : ViewType("dofCoeffs1",basis1_->getCardinality()); + auto dofCoeffs2 = isVectorBasis2 ? ViewType("dofCoeffs2",basis2_->getCardinality(), basisDim2) : ViewType("dofCoeffs2",basis2_->getCardinality()); + basis1_->getDofCoeffs(dofCoeffs1); basis2_->getDofCoeffs(dofCoeffs2); - const ordinal_type basisCardinality1 = basis1_->getCardinality(); - const ordinal_type basisCardinality2 = basis2_->getCardinality(); - Kokkos::RangePolicy policy(0, basisCardinality2); Kokkos::parallel_for(policy, KOKKOS_LAMBDA (const int fieldOrdinal2) { for (int fieldOrdinal1=0; fieldOrdinal1::result_layout; + using ViewType = Kokkos::DynRankView; + + const ordinal_type basisCardinality1 = basis1_->getCardinality(); + const ordinal_type basisCardinality2 = basis2_->getCardinality(); + const ordinal_type basisCardinality3 = basis3_->getCardinality(); + + auto dofCoeffs1 = ViewType("dofCoeffs1",basisCardinality1); + auto dofCoeffs2 = ViewType("dofCoeffs2",basisCardinality2); + auto dofCoeffs3 = ViewType("dofCoeffs3",basisCardinality3); + + basis1_->getDofCoeffs(dofCoeffs1); + basis2_->getDofCoeffs(dofCoeffs2); + basis3_->getDofCoeffs(dofCoeffs3); + + Kokkos::RangePolicy policy(0, basisCardinality3); + Kokkos::parallel_for(policy, KOKKOS_LAMBDA (const int fieldOrdinal3) + { + for (int fieldOrdinal2=0; fieldOrdinal2::createHDivProjectionStruct(const Ba hcurlBasis = new Basis_HCURL_HEX_In_FEM(cellBasis->getDegree()); else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) hcurlBasis = new Basis_HCURL_TET_In_FEM(cellBasis->getDegree()); + else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) + hcurlBasis = new typename DerivedNodalBasisFamily::HCURL_WEDGE(cellBasis->getDegree()); else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) hcurlBasis = new Basis_HGRAD_QUAD_Cn_FEM(cellBasis->getDegree()); else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) diff --git a/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHCURL.hpp b/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHCURL.hpp index 571bb8806296..e21b7d8cf6bb 100644 --- a/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHCURL.hpp +++ b/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHCURL.hpp @@ -650,7 +650,13 @@ ProjectionTools::getHCurlBasisCoeffs(Kokkos::DynRankView(cellBasis->getDegree(),POINTTYPE_WARPBLEND); else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) hgradBasis = new Basis_HGRAD_TRI_Cn_FEM(cellBasis->getDegree(),POINTTYPE_WARPBLEND); - else { + else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) { + if(iface < 3) + hgradBasis = new Basis_HGRAD_QUAD_Cn_FEM(cellBasis->getDegree(),POINTTYPE_WARPBLEND); + else + hgradBasis = new Basis_HGRAD_TRI_Cn_FEM(cellBasis->getDegree(),POINTTYPE_WARPBLEND); + } + else { std::stringstream ss; ss << ">>> ERROR (Intrepid2::ProjectionTools::getHCurlBasisCoeffs): " << "Method not implemented for basis " << name; @@ -760,6 +766,8 @@ ProjectionTools::getHCurlBasisCoeffs(Kokkos::DynRankView(cellBasis->getDegree()); else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) hgradBasis = new Basis_HGRAD_TET_Cn_FEM(cellBasis->getDegree(),POINTTYPE_WARPBLEND); + else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) + hgradBasis = new typename DerivedNodalBasisFamily::HGRAD_WEDGE(cellBasis->getDegree(),POINTTYPE_WARPBLEND); else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) hgradBasis = new Basis_HGRAD_TRI_Cn_FEM(cellBasis->getDegree(),POINTTYPE_WARPBLEND); else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) diff --git a/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHDIV.hpp b/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHDIV.hpp index e640799e186e..04607a847d84 100644 --- a/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHDIV.hpp +++ b/packages/intrepid2/src/Projection/Intrepid2_ProjectionToolsDefHDIV.hpp @@ -462,6 +462,8 @@ ProjectionTools::getHDivBasisCoeffs(Kokkos::DynRankView(cellBasis->getDegree()); else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) hcurlBasis = new Basis_HCURL_TET_In_FEM(cellBasis->getDegree()); + else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) + hcurlBasis = new typename DerivedNodalBasisFamily::HCURL_WEDGE(cellBasis->getDegree()); else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) hcurlBasis = new Basis_HGRAD_QUAD_Cn_FEM(cellBasis->getDegree()); else if(cellTopo.getKey() == shards::getCellTopologyData >()->key) diff --git a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/CMakeLists.txt b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/CMakeLists.txt index c55ac04519c7..08f7f992944c 100644 --- a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/CMakeLists.txt +++ b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/CMakeLists.txt @@ -18,7 +18,10 @@ LIST(LENGTH Intrepid2_TEST_ETI_VALUETYPE_NAME ETI_VALUETYPE_COUNT) MATH(EXPR ETI_VALUETYPE_COUNT "${ETI_VALUETYPE_COUNT}-1") # Host test -SET(Intrepid2_TEST_ETI_FILE "test_01") +SET(Intrepid2_TEST_ETI_FILE "") +LIST(APPEND Intrepid2_TEST_ETI_FILE + "test_01" + "test_01_Serendipity") SET(Intrepid2_TEST_ETI_DEVICE_NAME "") SET(Intrepid2_TEST_ETI_DEVICE "") diff --git a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/eti/test_01_ETI.in b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/eti/test_01_ETI.in index 56b887196100..0c4ee19515e4 100644 --- a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/eti/test_01_ETI.in +++ b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/eti/test_01_ETI.in @@ -41,7 +41,7 @@ // @HEADER /** \file test_01.cpp - \brief Unit tests for the Intrepid2::Basis_HGRAD_LINE_C2_FEM class. + \brief Unit tests for the Intrepid2::Basis_HGRAD_HEX_C2_FEM class. \author Created by P. Bochev, D. Ridzal, K. Peterson. */ @@ -54,7 +54,8 @@ int main(int argc, char *argv[]) { const bool verbose = (argc-1) > 0; Kokkos::initialize(); - const int r_val = Intrepid2::Test::HGRAD_HEX_C2_FEM_Test01<@ETI_VALUETYPE@, @ETI_DEVICE@>(verbose); + constexpr bool serendipity = false; + const int r_val = Intrepid2::Test::HGRAD_HEX_DEG2_FEM_Test01(verbose); Kokkos::finalize(); return r_val; diff --git a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/eti/test_02_ETI.in b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/eti/test_01_Serendipity_ETI.in similarity index 88% rename from packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/eti/test_02_ETI.in rename to packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/eti/test_01_Serendipity_ETI.in index e7b6476dcf2c..a0846e23f9e1 100644 --- a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/eti/test_02_ETI.in +++ b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/eti/test_01_Serendipity_ETI.in @@ -40,21 +40,22 @@ // ************************************************************************ // @HEADER -/** \file test_02.cpp - \brief Unit tests for the Intrepid2::Basis_HGRAD_LINE_C2_FEM class. +/** \file test_01.cpp + \brief Unit tests for the Intrepid2::Basis_HGRAD_HEX_I2_Serendipiy_FEM class. \author Created by P. Bochev, D. Ridzal, K. Peterson. */ #include "Kokkos_Core.hpp" -#include "test_02.hpp" +#include "test_01.hpp" int main(int argc, char *argv[]) { const bool verbose = (argc-1) > 0; Kokkos::initialize(); - const int r_val = Intrepid2::Test::HGRAD_HEX_C2_FEM_Test02<@ETI_VALUETYPE@, @ETI_DEVICE@>(verbose); + constexpr bool serendipity = true; + const int r_val = Intrepid2::Test::HGRAD_HEX_DEG2_FEM_Test01(verbose); Kokkos::finalize(); return r_val; diff --git a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/test_01.hpp b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/test_01.hpp index a6166ba996ac..7db161882d0e 100644 --- a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/test_01.hpp +++ b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_HEX_C2_FEM/test_01.hpp @@ -41,7 +41,7 @@ // @HEADER /** \file test_01.hpp - \brief Unit tests for the Intrepid2::HGRAD_HEX_C2_FEM class. + \brief Unit tests for the Intrepid2::HGRAD_HEX_DEG2_FEM class. \author Created by P. Bochev, D. Ridzal, K. Peterson and Kyungjoo Kim */ #include "Intrepid2_config.h" @@ -72,8 +72,8 @@ namespace Intrepid2 { } - template - int HGRAD_HEX_C2_FEM_Test01(const bool verbose) { + template + int HGRAD_HEX_DEG2_FEM_Test01(const bool verbose) { Teuchos::RCP outStream; Teuchos::oblackholestream bhs; // outputs nothing @@ -96,22 +96,21 @@ namespace Intrepid2 { *outStream << "\n" << "===============================================================================\n" - << "| |\n" - << "| Unit Test (Basis_HGRAD_HEX_C2_FEM) |\n" + << "| |\n"; + + if constexpr (serendipity) + *outStream + << "| Unit Test (Basis_HGRAD_HEX_I2_Serendipity FEM) |\n"; + else + *outStream + << "| Unit Test (Basis_HGRAD_HEX_C2_FEM) |\n"; + *outStream << "| |\n" << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" << "| 2) Basis values for VALUE, GRAD, and Dk operators |\n" << "| |\n" - << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" - << "| Denis Ridzal (dridzal@sandia.gov), |\n" - << "| Kara Peterson (kjpeter@sandia.gov). |\n" - << "| |\n" - << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" - << "| Trilinos website: http://trilinos.sandia.gov |\n" - << "| |\n" << "===============================================================================\n"; - typedef Kokkos::DynRankView DynRankView; typedef Kokkos::DynRankView DynRankViewHost; #define ConstructWithLabel(obj, ...) obj(#obj, __VA_ARGS__) @@ -122,9 +121,11 @@ namespace Intrepid2 { // for virtual function, value and point types are declared in the class typedef ValueType outputValueType; typedef ValueType pointValueType; - Basis_HGRAD_HEX_C2_FEM hexBasis; - //typedef typename decltype(hexBasis)::OutputViewType OutputViewType; - //typedef typename decltype(hexBasis)::PointViewType PointViewType; + BasisPtr hexBasis; + if constexpr (serendipity) + hexBasis = Teuchos::rcp(new Basis_HGRAD_HEX_I2_Serendipity_FEM()); + else + hexBasis = Teuchos::rcp(new Basis_HGRAD_HEX_C2_FEM()); *outStream << "\n" @@ -140,9 +141,9 @@ namespace Intrepid2 { DynRankView ConstructWithLabel( hexNodes, 27, 3); // Generic array for the output values; needs to be properly resized depending on the operator type - const ordinal_type numFields = hexBasis.getCardinality(); + const ordinal_type numFields = hexBasis->getCardinality(); const ordinal_type numPoints = hexNodes.extent(0); - const ordinal_type spaceDim = hexBasis.getBaseCellTopology().getDimension(); + const ordinal_type spaceDim = hexBasis->getBaseCellTopology().getDimension(); const ordinal_type D2Cardin = getDkCardinality(OPERATOR_D2, spaceDim); const ordinal_type workSize = numFields*numPoints*D2Cardin; @@ -153,77 +154,77 @@ namespace Intrepid2 { // exception #1: CURL cannot be applied to scalar functions in 3D // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary) DynRankView tmpvals = DynRankView(work.data(), numFields, numPoints, 4); - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getValues(tmpvals, hexNodes, OPERATOR_CURL) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getValues(tmpvals, hexNodes, OPERATOR_CURL) ); } { // exception #2: DIV cannot be applied to scalar functions in 3D // resize vals to rank-2 container with dimensions (num. basis functions, num. points) DynRankView tmpvals = DynRankView(work.data(), numFields, numPoints); - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getValues(tmpvals, hexNodes, OPERATOR_DIV) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getValues(tmpvals, hexNodes, OPERATOR_DIV) ); } // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and // getDofTag() to access invalid array elements thereby causing bounds check exception { // exception #3 - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getDofOrdinal(3,10,0) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getDofOrdinal(3,10,0) ); // exception #4 - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getDofOrdinal(1,2,1) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getDofOrdinal(1,2,1) ); // exception #5 - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getDofOrdinal(0,4,1) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getDofOrdinal(0,4,1) ); // exception #6 - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getDofTag(numFields) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getDofTag(numFields) ); // exception #7 - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getDofTag(-1) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getDofTag(-1) ); } // Exceptions 8-18 test exception handling with incorrectly dimensioned input/output arrays { // exception #8: input points array must be of rank-2 DynRankView ConstructWithLabel(badPoints1, 4, 5, 3); - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getValues(vals, badPoints1, OPERATOR_VALUE) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getValues(vals, badPoints1, OPERATOR_VALUE) ); } { // exception #9 dimension 1 in the input point array must equal space dimension of the cell DynRankView ConstructWithLabel(badPoints2, 4, spaceDim - 1); - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getValues(vals, badPoints2, OPERATOR_VALUE) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getValues(vals, badPoints2, OPERATOR_VALUE) ); } { // exception #10 output values must be of rank-2 for OPERATOR_VALUE DynRankView ConstructWithLabel(badVals1, 4, 3, 1); - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getValues(badVals1, hexNodes, OPERATOR_VALUE) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getValues(badVals1, hexNodes, OPERATOR_VALUE) ); } { // exception #11 output values must be of rank-3 for OPERATOR_GRAD DynRankView ConstructWithLabel(badVals2, 4, 3); - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getValues(badVals2, hexNodes, OPERATOR_GRAD) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getValues(badVals2, hexNodes, OPERATOR_GRAD) ); // exception #12 output values must be of rank-3 for OPERATOR_D1 - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getValues(badVals2, hexNodes, OPERATOR_D1) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getValues(badVals2, hexNodes, OPERATOR_D1) ); // exception #13 output values must be of rank-3 for OPERATOR_D2 - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getValues(badVals2, hexNodes, OPERATOR_D2) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getValues(badVals2, hexNodes, OPERATOR_D2) ); } { // exception #14 incorrect 0th dimension of output array (must equal number of basis functions) DynRankView ConstructWithLabel(badVals3, numFields + 1, numPoints); - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getValues(badVals3, hexNodes, OPERATOR_VALUE) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getValues(badVals3, hexNodes, OPERATOR_VALUE) ); } { // exception #15 incorrect 1st dimension of output array (must equal number of points) DynRankView ConstructWithLabel(badVals4, numFields, numPoints + 1); - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getValues(badVals4, hexNodes, OPERATOR_VALUE) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getValues(badVals4, hexNodes, OPERATOR_VALUE) ); } { // exception #16: incorrect 2nd dimension of output array (must equal the space dimension) DynRankView ConstructWithLabel(badVals5, numFields, numPoints, spaceDim - 1); - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getValues(badVals5, hexNodes, OPERATOR_GRAD) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getValues(badVals5, hexNodes, OPERATOR_GRAD) ); } { // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 3D) DynRankView ConstructWithLabel(badVals6, numFields, numPoints, 40); - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getValues(badVals6, hexNodes, OPERATOR_D2) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getValues(badVals6, hexNodes, OPERATOR_D2) ); } { // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 3D) DynRankView ConstructWithLabel(badVals7, numFields, numPoints, 50); - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getValues(badVals7, hexNodes, OPERATOR_D3) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getValues(badVals7, hexNodes, OPERATOR_D3) ); } #endif // Check if number of thrown exceptions matches the one we expect @@ -248,15 +249,15 @@ namespace Intrepid2 { try{ - const ordinal_type numFields = hexBasis.getCardinality(); - const auto allTags = hexBasis.getAllDofTags(); + const ordinal_type numFields = hexBasis->getCardinality(); + const auto allTags = hexBasis->getAllDofTags(); // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again const ordinal_type dofTagSize = allTags.extent(0); for (ordinal_type i = 0; i < dofTagSize; ++i) { - const auto bfOrd = hexBasis.getDofOrdinal(allTags(i,0), allTags(i,1), allTags(i,2)); + const auto bfOrd = hexBasis->getDofOrdinal(allTags(i,0), allTags(i,1), allTags(i,2)); - const auto myTag = hexBasis.getDofTag(bfOrd); + const auto myTag = hexBasis->getDofTag(bfOrd); if( !( (myTag(0) == allTags(i,0)) && (myTag(1) == allTags(i,1)) && (myTag(2) == allTags(i,2)) && @@ -278,8 +279,8 @@ namespace Intrepid2 { // Now do the same but loop over basis functions for( ordinal_type bfOrd = 0; bfOrd < numFields; ++bfOrd) { - const auto myTag = hexBasis.getDofTag(bfOrd); - const auto myBfOrd = hexBasis.getDofOrdinal(myTag(0), myTag(1), myTag(2)); + const auto myTag = hexBasis->getDofTag(bfOrd); + const auto myBfOrd = hexBasis->getDofOrdinal(myTag(0), myTag(1), myTag(2)); if( bfOrd != myBfOrd) { errorFlag++; *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; @@ -309,7 +310,8 @@ namespace Intrepid2 { outStream -> precision(20); try{ - // VALUE: Each row gives the 8 correct basis set values at an evaluation point + const int numC2Fields = 27; + // VALUE: Each row gives the 27 correct basis set values at an evaluation point const ValueType basisValues[] = { 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \ 0, 0, 0, 0, 0, 0, 1.000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \ @@ -475,9 +477,9 @@ namespace Intrepid2 { Kokkos::deep_copy(hexNodes, hexNodesHost); // Dimensions for the output arrays: - const ordinal_type numFields = hexBasis.getCardinality(); + const ordinal_type numFields = hexBasis->getCardinality(); const ordinal_type numPoints = hexNodes.extent(0); - const ordinal_type spaceDim = hexBasis.getBaseCellTopology().getDimension(); + const ordinal_type spaceDim = hexBasis->getBaseCellTopology().getDimension(); const ordinal_type D2Cardin = getDkCardinality(OPERATOR_D2, spaceDim); const ordinal_type D3Cardin = getDkCardinality(OPERATOR_D3, spaceDim); const ordinal_type D4Cardin = getDkCardinality(OPERATOR_D4, spaceDim); @@ -486,12 +488,12 @@ namespace Intrepid2 { // Generic array for values, grads, curls, etc. that will be properly sized before each call DynRankView ConstructWithLabel(vals, numFields, numPoints); // Check VALUE of basis functions: resize vals to rank-2 container: - hexBasis.getValues(vals, hexNodes, OPERATOR_VALUE); + hexBasis->getValues(vals, hexNodes, OPERATOR_VALUE); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i = 0; i < numFields; ++i) { for (ordinal_type j = 0; j < numPoints; ++j) { - const ordinal_type l = i + j * numFields; + const ordinal_type l = i + j * numC2Fields; if (std::abs(vals_host(i,j) - basisValues[l]) > tol ) { errorFlag++; *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; @@ -509,7 +511,7 @@ namespace Intrepid2 { { DynRankView ConstructWithLabel(vals, numFields, numPoints, spaceDim); // Check GRAD of basis function: resize vals to rank-3 container - hexBasis.getValues(vals, hexNodes, OPERATOR_GRAD); + hexBasis->getValues(vals, hexNodes, OPERATOR_GRAD); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i = 0; i < numFields; ++i) { @@ -533,7 +535,7 @@ namespace Intrepid2 { } // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD) - hexBasis.getValues(vals, hexNodes, OPERATOR_D1); + hexBasis->getValues(vals, hexNodes, OPERATOR_D1); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i = 0; i < numFields; ++i) { for (ordinal_type j = 0; j < numPoints; ++j) { @@ -559,7 +561,7 @@ namespace Intrepid2 { { // Check D2 of basis function DynRankView ConstructWithLabel(vals, numFields, numPoints, D2Cardin); - hexBasis.getValues(vals, hexNodes, OPERATOR_D2); + hexBasis->getValues(vals, hexNodes, OPERATOR_D2); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i = 0; i < numFields; ++i) { @@ -586,7 +588,7 @@ namespace Intrepid2 { { // Check D3 of basis function DynRankView ConstructWithLabel(vals, numFields, numPoints, D3Cardin); - hexBasis.getValues(vals, hexNodes, OPERATOR_D3); + hexBasis->getValues(vals, hexNodes, OPERATOR_D3); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i = 0; i < numFields; ++i) { @@ -613,7 +615,7 @@ namespace Intrepid2 { { // Check D4 of basis function DynRankView ConstructWithLabel(vals, numFields, numPoints, D4Cardin); - hexBasis.getValues(vals, hexNodes, OPERATOR_D4); + hexBasis->getValues(vals, hexNodes, OPERATOR_D4); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i = 0; i < numFields; i++) { @@ -650,7 +652,7 @@ namespace Intrepid2 { // The last dimension is the number of kth derivatives and needs to be resized for every Dk const ordinal_type DkCardin = getDkCardinality(op, spaceDim); DynRankView ConstructWithLabel(vals, numFields, numPoints, DkCardin); - hexBasis.getValues(vals, hexNodes, op); + hexBasis->getValues(vals, hexNodes, op); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); @@ -682,23 +684,23 @@ namespace Intrepid2 { << "===============================================================================\n"; try{ - const ordinal_type numFields = hexBasis.getCardinality(); - const ordinal_type spaceDim = hexBasis.getBaseCellTopology().getDimension(); + const ordinal_type numFields = hexBasis->getCardinality(); + const ordinal_type spaceDim = hexBasis->getBaseCellTopology().getDimension(); // Check exceptions. ordinal_type nthrow = 0, ncatch = 0; #ifdef HAVE_INTREPID2_DEBUG { DynRankView ConstructWithLabel(badVals, 1, 2, 3); - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getDofCoords(badVals) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getDofCoords(badVals) ); } { DynRankView ConstructWithLabel(badVals, 3, 2); - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getDofCoords(badVals) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getDofCoords(badVals) ); } { DynRankView ConstructWithLabel(badVals, 27, 2); - INTREPID2_TEST_ERROR_EXPECTED( hexBasis.getDofCoords(badVals) ); + INTREPID2_TEST_ERROR_EXPECTED( hexBasis->getDofCoords(badVals) ); } #endif // Check if number of thrown exceptions matches the one we expect @@ -712,8 +714,8 @@ namespace Intrepid2 { DynRankView ConstructWithLabel(cvals, numFields, spaceDim); // Check mathematical correctness. - hexBasis.getDofCoords(cvals); - hexBasis.getValues(bvals, cvals, OPERATOR_VALUE); + hexBasis->getDofCoords(cvals); + hexBasis->getValues(bvals, cvals, OPERATOR_VALUE); auto cvals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), cvals); Kokkos::deep_copy(cvals_host, cvals); auto bvals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), bvals); @@ -747,7 +749,7 @@ namespace Intrepid2 { << "===============================================================================\n"; try { - const EFunctionSpace fs = hexBasis.getFunctionSpace(); + const EFunctionSpace fs = hexBasis->getFunctionSpace(); if (fs != FUNCTION_SPACE_HGRAD) { diff --git a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_QUAD_C2_FEM/CMakeLists.txt b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_QUAD_C2_FEM/CMakeLists.txt index 3df9f211f62d..44c6416915d0 100644 --- a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_QUAD_C2_FEM/CMakeLists.txt +++ b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_QUAD_C2_FEM/CMakeLists.txt @@ -18,7 +18,10 @@ LIST(LENGTH Intrepid2_TEST_ETI_VALUETYPE_NAME ETI_VALUETYPE_COUNT) MATH(EXPR ETI_VALUETYPE_COUNT "${ETI_VALUETYPE_COUNT}-1") # Host test -SET(Intrepid2_TEST_ETI_FILE "test_01") +SET(Intrepid2_TEST_ETI_FILE "") +LIST(APPEND Intrepid2_TEST_ETI_FILE + "test_01" + "test_01_Serendipity") SET(Intrepid2_TEST_ETI_DEVICE_NAME "") SET(Intrepid2_TEST_ETI_DEVICE "") diff --git a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_QUAD_C2_FEM/eti/test_01_ETI.in b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_QUAD_C2_FEM/eti/test_01_ETI.in index d834366421d6..4cd74022533d 100644 --- a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_QUAD_C2_FEM/eti/test_01_ETI.in +++ b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_QUAD_C2_FEM/eti/test_01_ETI.in @@ -41,7 +41,7 @@ // @HEADER /** \file test_01.hpp - \brief Unit tests for the Intrepid2::Basis_HGRAD_LINE_C2_FEM class. + \brief Unit tests for the Intrepid2::Basis_HGRAD_QUAD_C2_FEM class. \author Created by P. Bochev, D. Ridzal, K. Peterson. */ @@ -54,7 +54,8 @@ int main(int argc, char *argv[]) { const bool verbose = (argc-1) > 0; Kokkos::initialize(); - const int r_val = Intrepid2::Test::HGRAD_QUAD_C2_FEM_Test01<@ETI_VALUETYPE@, @ETI_DEVICE@>(verbose); + constexpr bool serendipity = false; + const int r_val = Intrepid2::Test::HGRAD_QUAD_DEG2_FEM_Test01(verbose); Kokkos::finalize(); return r_val; diff --git a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_QUAD_C2_FEM/eti/test_01_Serendipity_ETI.in b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_QUAD_C2_FEM/eti/test_01_Serendipity_ETI.in new file mode 100644 index 000000000000..c022768204b5 --- /dev/null +++ b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_QUAD_C2_FEM/eti/test_01_Serendipity_ETI.in @@ -0,0 +1,63 @@ +// @HEADER +// ************************************************************************ +// +// Intrepid2 Package +// Copyright (2007) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Kyungjoo Kim (kyukim@sandia.gov), or +// Mauro Perego (mperego@sandia.gov) +// +// ************************************************************************ +// @HEADER + +/** \file test_01.hpp + \brief Unit tests for the Intrepid2::Basis_HGRAD_QUAD_I2_Serendipiy_FEM class. + \author Created by P. Bochev, D. Ridzal, K. Peterson. +*/ + +#include "Kokkos_Core.hpp" + +#include "test_01.hpp" + +int main(int argc, char *argv[]) { + + const bool verbose = (argc-1) > 0; + Kokkos::initialize(); + + constexpr bool serendipity = true; + const int r_val = Intrepid2::Test::HGRAD_QUAD_DEG2_FEM_Test01(verbose); + + Kokkos::finalize(); + return r_val; +} + diff --git a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_QUAD_C2_FEM/test_01.hpp b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_QUAD_C2_FEM/test_01.hpp index 7b36ba195404..7e2a6c210acd 100644 --- a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_QUAD_C2_FEM/test_01.hpp +++ b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_QUAD_C2_FEM/test_01.hpp @@ -41,7 +41,7 @@ // @HEADER /** \file test_01.hpp - \brief Unit tests for the Intrepid2::HGRAD_QUAD_C2_FEM class. + \brief Unit tests for the Intrepid2::HGRAD_QUAD_DEG2_FEM class. \author Created by P. Bochev, D. Ridzal, K. Peterson, and Kyungjoo Kim. */ @@ -72,8 +72,8 @@ namespace Intrepid2 { *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \ } - template - int HGRAD_QUAD_C2_FEM_Test01(const bool verbose) { + template + int HGRAD_QUAD_DEG2_FEM_Test01(const bool verbose) { Teuchos::RCP outStream; Teuchos::oblackholestream bhs; // outputs nothing @@ -95,19 +95,18 @@ namespace Intrepid2 { *outStream << "===============================================================================\n" - << "| |\n" - << "| Unit Test (Basis_HGRAD_QUAD_C2_FEM) |\n" + << "| |\n"; + if constexpr (serendipity) + *outStream + << "| Unit Test (Basis_HGRAD_QUAD_I2_Serendipity FEM) |\n"; + else + *outStream + << "| Unit Test (Basis_HGRAD_QUAD_C2_FEM) |\n"; + *outStream << "| |\n" << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" << "| 2) Basis values for VALUE, GRAD, CURL, and Dk operators |\n" - << "| |\n" - << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" - << "| Denis Ridzal (dridzal@sandia.gov), |\n" - << "| Kara Peterson (kjpeter@sandia.gov). |\n" - << "| |\n" - << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" - << "| Trilinos website: http://trilinos.sandia.gov |\n" - << "| |\n" + << "| |\n" << "| |\n" << "===============================================================================\n"; typedef Kokkos::DynRankView DynRankView; @@ -120,9 +119,11 @@ namespace Intrepid2 { // for virtual function, value and point types are declared in the class typedef ValueType outputValueType; typedef ValueType pointValueType; - Basis_HGRAD_QUAD_C2_FEM quadBasis; - //typedef typename decltype(quadBasis)::OutputViewType OutputViewType; - //typedef typename decltype(quadBasis)::PointViewType PointViewType; + BasisPtr quadBasis; + if constexpr (serendipity) + quadBasis = Teuchos::rcp(new Basis_HGRAD_QUAD_I2_Serendipity_FEM()); + else + quadBasis = Teuchos::rcp(new Basis_HGRAD_QUAD_C2_FEM()); *outStream << "\n" @@ -135,9 +136,9 @@ namespace Intrepid2 { #ifdef HAVE_INTREPID2_DEBUG DynRankView ConstructWithLabel(quadNodes, 10, 2); - const ordinal_type numFields = quadBasis.getCardinality(); + const ordinal_type numFields = quadBasis->getCardinality(); const ordinal_type numPoints = quadNodes.extent(0); - const ordinal_type spaceDim = quadBasis.getBaseCellTopology().getDimension(); + const ordinal_type spaceDim = quadBasis->getBaseCellTopology().getDimension(); const ordinal_type D2Cardin = getDkCardinality(OPERATOR_D2, spaceDim); const ordinal_type workSize = numFields*numFields*D2Cardin; @@ -148,68 +149,68 @@ namespace Intrepid2 { { // exception #1: DIV cannot be applied to scalar functions - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getValues(vals, quadNodes, OPERATOR_DIV) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getValues(vals, quadNodes, OPERATOR_DIV) ); } // Exceptions 2-6: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and // getDofTag() to access invalid array elements thereby causing bounds check exception { // exception #2 - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getDofOrdinal(3,0,0) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getDofOrdinal(3,0,0) ); // exception #3 - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getDofOrdinal(1,1,1) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getDofOrdinal(1,1,1) ); // exception #4 - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getDofOrdinal(0,4,0) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getDofOrdinal(0,4,0) ); // exception #5 - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getDofTag(numFields) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getDofTag(numFields) ); // exception #6 - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getDofTag(-1) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getDofTag(-1) ); } // Exceptions 7- test exception handling with incorrectly dimensioned input/output arrays { // exception #7: input points array must be of rank-2 DynRankView ConstructWithLabel(badPoints1, 4, 5, 3); - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getValues(vals, badPoints1, OPERATOR_VALUE) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getValues(vals, badPoints1, OPERATOR_VALUE) ); } { // exception #8 dimension 1 in the input point array must equal space dimension of the cell DynRankView ConstructWithLabel(badPoints2, 4, spaceDim + 1); - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getValues(vals, badPoints2, OPERATOR_VALUE) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getValues(vals, badPoints2, OPERATOR_VALUE) ); } { // exception #9 output values must be of rank-2 for OPERATOR_VALUE DynRankView ConstructWithLabel(badVals1, 4, 3, 1); - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getValues(badVals1, quadNodes, OPERATOR_VALUE) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getValues(badVals1, quadNodes, OPERATOR_VALUE) ); } { // exception #10 output values must be of rank-3 for OPERATOR_GRAD DynRankView ConstructWithLabel(badVals2, 4, 3); - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getValues(badVals2, quadNodes, OPERATOR_GRAD) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getValues(badVals2, quadNodes, OPERATOR_GRAD) ); // exception #11 output values must be of rank-3 for OPERATOR_CURL - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getValues(badVals2, quadNodes, OPERATOR_CURL) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getValues(badVals2, quadNodes, OPERATOR_CURL) ); // exception #12 output values must be of rank-3 for OPERATOR_D2 - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getValues(badVals2, quadNodes, OPERATOR_D2) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getValues(badVals2, quadNodes, OPERATOR_D2) ); } { // exception #13 incorrect 0th dimension of output array (must equal number of basis functions) DynRankView ConstructWithLabel(badVals3, numFields + 1, numPoints); - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getValues(badVals3, quadNodes, OPERATOR_VALUE) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getValues(badVals3, quadNodes, OPERATOR_VALUE) ); } { // exception #14 incorrect 1st dimension of output array (must equal number of points in quadNodes) DynRankView ConstructWithLabel(badVals4, numFields, numPoints + 1); - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getValues(badVals4, quadNodes, OPERATOR_VALUE) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getValues(badVals4, quadNodes, OPERATOR_VALUE) ); } { // exception #15: incorrect 2nd dimension of output array (must equal the space dimension) DynRankView ConstructWithLabel(badVals5, numFields, numPoints, spaceDim + 1); - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getValues(badVals5, quadNodes, OPERATOR_GRAD) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getValues(badVals5, quadNodes, OPERATOR_GRAD) ); } { // exception #16: incorrect 2nd dimension of output array (must equal D2 cardinality in 2D) DynRankView ConstructWithLabel(badVals6, numFields, numPoints, 40); - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getValues(badVals6, quadNodes, OPERATOR_D2) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getValues(badVals6, quadNodes, OPERATOR_D2) ); // exception #17: incorrect 2nd dimension of output array (must equal D3 cardinality in 2D) - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getValues(badVals6, quadNodes, OPERATOR_D3) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getValues(badVals6, quadNodes, OPERATOR_D3) ); } #endif // Check if number of thrown exceptions matches the one we expect @@ -232,15 +233,15 @@ namespace Intrepid2 { << "===============================================================================\n"; try{ - const ordinal_type numFields = quadBasis.getCardinality(); - const auto allTags = quadBasis.getAllDofTags(); + const ordinal_type numFields = quadBasis->getCardinality(); + const auto allTags = quadBasis->getAllDofTags(); // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again const ordinal_type dofTagSize = allTags.extent(0); for (ordinal_type i = 0; i < dofTagSize; ++i) { - const auto bfOrd = quadBasis.getDofOrdinal(allTags(i,0), allTags(i,1), allTags(i,2)); + const auto bfOrd = quadBasis->getDofOrdinal(allTags(i,0), allTags(i,1), allTags(i,2)); - const auto myTag = quadBasis.getDofTag(bfOrd); + const auto myTag = quadBasis->getDofTag(bfOrd); if( !( (myTag(0) == allTags(i,0)) && (myTag(1) == allTags(i,1)) && (myTag(2) == allTags(i,2)) && @@ -262,8 +263,8 @@ namespace Intrepid2 { // Now do the same but loop over basis functions for( ordinal_type bfOrd = 0; bfOrd < numFields; bfOrd++) { - const auto myTag = quadBasis.getDofTag(bfOrd); - const auto myBfOrd = quadBasis.getDofOrdinal(myTag(0), myTag(1), myTag(2)); + const auto myTag = quadBasis->getDofTag(bfOrd); + const auto myBfOrd = quadBasis->getDofOrdinal(myTag(0), myTag(1), myTag(2)); if( bfOrd != myBfOrd) { errorFlag++; *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; @@ -517,9 +518,9 @@ namespace Intrepid2 { Kokkos::deep_copy(quadNodes, quadNodesHost); // Dimensions for the output arrays: - const ordinal_type numFields = quadBasis.getCardinality(); + const ordinal_type numFields = quadBasis->getCardinality(); const ordinal_type numPoints = quadNodes.extent(0); - const ordinal_type spaceDim = quadBasis.getBaseCellTopology().getDimension(); + const ordinal_type spaceDim = quadBasis->getBaseCellTopology().getDimension(); const ordinal_type D2Cardin = getDkCardinality(OPERATOR_D2, spaceDim); const ordinal_type D3Cardin = getDkCardinality(OPERATOR_D3, spaceDim); const ordinal_type D4Cardin = getDkCardinality(OPERATOR_D4, spaceDim); @@ -527,7 +528,7 @@ namespace Intrepid2 { { // Check VALUE of basis functions: resize vals to rank-2 container: DynRankView ConstructWithLabel(vals, numFields, numPoints); - quadBasis.getValues(vals, quadNodes, OPERATOR_VALUE); + quadBasis->getValues(vals, quadNodes, OPERATOR_VALUE); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i = 0; i < numFields; ++i) { @@ -551,7 +552,7 @@ namespace Intrepid2 { { // Check GRAD of basis function: resize vals to rank-3 container DynRankView ConstructWithLabel(vals, numFields, numPoints, spaceDim); - quadBasis.getValues(vals, quadNodes, OPERATOR_GRAD); + quadBasis->getValues(vals, quadNodes, OPERATOR_GRAD); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i = 0; i < numFields; ++i) { @@ -577,7 +578,7 @@ namespace Intrepid2 { { // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD) DynRankView ConstructWithLabel(vals, numFields, numPoints, spaceDim); - quadBasis.getValues(vals, quadNodes, OPERATOR_D1); + quadBasis->getValues(vals, quadNodes, OPERATOR_D1); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i = 0; i < numFields; ++i) { @@ -603,7 +604,7 @@ namespace Intrepid2 { { // Check CURL of basis function: resize vals just for illustration! DynRankView ConstructWithLabel(vals, numFields, numPoints, spaceDim); - quadBasis.getValues(vals, quadNodes, OPERATOR_CURL); + quadBasis->getValues(vals, quadNodes, OPERATOR_CURL); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i = 0; i < numFields; ++i) { @@ -638,7 +639,7 @@ namespace Intrepid2 { { // Check D2 of basis function DynRankView ConstructWithLabel(vals, numFields, numPoints, D2Cardin); - quadBasis.getValues(vals, quadNodes, OPERATOR_D2); + quadBasis->getValues(vals, quadNodes, OPERATOR_D2); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i = 0; i < numFields; ++i) { @@ -664,7 +665,7 @@ namespace Intrepid2 { { // Check D3 of basis function DynRankView ConstructWithLabel(vals, numFields, numPoints, D3Cardin); - quadBasis.getValues(vals, quadNodes, OPERATOR_D3); + quadBasis->getValues(vals, quadNodes, OPERATOR_D3); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); @@ -691,7 +692,7 @@ namespace Intrepid2 { { // Check D4 of basis function DynRankView ConstructWithLabel(vals, numFields, numPoints, D4Cardin); - quadBasis.getValues(vals, quadNodes, OPERATOR_D4); + quadBasis->getValues(vals, quadNodes, OPERATOR_D4); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); @@ -730,7 +731,7 @@ namespace Intrepid2 { // The last dimension is the number of kth derivatives and needs to be resized for every Dk const ordinal_type DkCardin = getDkCardinality(op, spaceDim); DynRankView ConstructWithLabel(vals, numFields, numPoints, DkCardin); - quadBasis.getValues(vals, quadNodes, op); + quadBasis->getValues(vals, quadNodes, op); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); @@ -762,10 +763,9 @@ namespace Intrepid2 { << "===============================================================================\n"; try{ - Basis_HGRAD_QUAD_C2_FEM quadB; - const ordinal_type numFields = quadB.getCardinality(); + const ordinal_type numFields = quadBasis->getCardinality(); const ordinal_type spaceDim = - quadB.getBaseCellTopology().getDimension(); + quadBasis->getBaseCellTopology().getDimension(); // Check exceptions. ordinal_type nthrow = 0, ncatch = 0; @@ -774,28 +774,31 @@ namespace Intrepid2 { #ifdef HAVE_INTREPID2_DEBUG { DynRankView ConstructWithLabel(badVals, 1, 2, 3); - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getDofCoords(badVals) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getDofCoords(badVals) ); } + + if constexpr(!serendipity) { DynRankView ConstructWithLabel(badVals, 8, 2); - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getDofCoords(badVals) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getDofCoords(badVals) ); } + { DynRankView ConstructWithLabel(badVals, 9, 1); - INTREPID2_TEST_ERROR_EXPECTED( quadBasis.getDofCoords(badVals) ); + INTREPID2_TEST_ERROR_EXPECTED( quadBasis->getDofCoords(badVals) ); } #endif if (nthrow != ncatch) { errorFlag++; *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; - *outStream << "# of catch ("<< ncatch << ") is different from # of throw (" << ncatch << ")\n"; + *outStream << "# of catch ("<< ncatch << ") is different from # of throw (" << nthrow << ")\n"; } DynRankView ConstructWithLabel(bvals_dev, numFields, numFields); DynRankView ConstructWithLabel(cvals_dev, numFields, spaceDim); - quadB.getDofCoords(cvals_dev); - quadB.getValues(bvals_dev, cvals_dev, OPERATOR_VALUE); + quadBasis->getDofCoords(cvals_dev); + quadBasis->getValues(bvals_dev, cvals_dev, OPERATOR_VALUE); auto bvals = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), bvals_dev); Kokkos::deep_copy(bvals, bvals_dev); @@ -832,7 +835,7 @@ namespace Intrepid2 { << "===============================================================================\n"; try { - const EFunctionSpace fs = quadBasis.getFunctionSpace(); + const EFunctionSpace fs = quadBasis->getFunctionSpace(); if (fs != FUNCTION_SPACE_HGRAD) { diff --git a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_WEDGE_C2_FEM/CMakeLists.txt b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_WEDGE_C2_FEM/CMakeLists.txt index 9de6c70a099c..45b891e7b3d8 100644 --- a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_WEDGE_C2_FEM/CMakeLists.txt +++ b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_WEDGE_C2_FEM/CMakeLists.txt @@ -18,7 +18,10 @@ LIST(LENGTH Intrepid2_TEST_ETI_VALUETYPE_NAME ETI_VALUETYPE_COUNT) MATH(EXPR ETI_VALUETYPE_COUNT "${ETI_VALUETYPE_COUNT}-1") # Host test -SET(Intrepid2_TEST_ETI_FILE "test_01") +SET(Intrepid2_TEST_ETI_FILE "") +LIST(APPEND Intrepid2_TEST_ETI_FILE + "test_01" + "test_01_Serendipity") SET(Intrepid2_TEST_ETI_DEVICE_NAME "") SET(Intrepid2_TEST_ETI_DEVICE "") diff --git a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_WEDGE_C2_FEM/eti/test_01_ETI.in b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_WEDGE_C2_FEM/eti/test_01_ETI.in index 0350be055ffe..d3a241ef8349 100644 --- a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_WEDGE_C2_FEM/eti/test_01_ETI.in +++ b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_WEDGE_C2_FEM/eti/test_01_ETI.in @@ -41,7 +41,7 @@ // @HEADER /** \file test_01.cpp - \brief Unit tests for the Intrepid2::Basis_HGRAD_WEDGE_C1_FEM class. + \brief Unit tests for the Intrepid2::Basis_HGRAD_WEDGE_C2_FEM class. \author Created by Kyungjoo Kim */ @@ -54,7 +54,8 @@ int main(int argc, char *argv[]) { const bool verbose = (argc-1) > 0; Kokkos::initialize(); - const int r_val = Intrepid2::Test::HGRAD_WEDGE_C2_FEM_Test01<@ETI_VALUETYPE@, @ETI_DEVICE@>(verbose); + constexpr bool serendipity = false; + const int r_val = Intrepid2::Test::HGRAD_WEDGE_DEG2_FEM_Test01(verbose); Kokkos::finalize(); return r_val; diff --git a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_WEDGE_C2_FEM/eti/test_01_Serendipity_ETI.in b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_WEDGE_C2_FEM/eti/test_01_Serendipity_ETI.in new file mode 100644 index 000000000000..8e7c0562bba0 --- /dev/null +++ b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_WEDGE_C2_FEM/eti/test_01_Serendipity_ETI.in @@ -0,0 +1,63 @@ +// @HEADER +// ************************************************************************ +// +// Intrepid2 Package +// Copyright (2007) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Kyungjoo Kim (kyukim@sandia.gov), or +// Mauro Perego (mperego@sandia.gov) +// +// ************************************************************************ +// @HEADER + +/** \file test_01.cpp + \brief Unit tests for the Intrepid2::Basis_HGRAD_WEDGE_I2_Serendipity_FEM classes. + \author Created by Kyungjoo Kim +*/ + +#include "Kokkos_Core.hpp" + +#include "test_01.hpp" + +int main(int argc, char *argv[]) { + + const bool verbose = (argc-1) > 0; + Kokkos::initialize(); + + constexpr bool serendipity = true; + const int r_val = Intrepid2::Test::HGRAD_WEDGE_DEG2_FEM_Test01(verbose); + + Kokkos::finalize(); + return r_val; +} + diff --git a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_WEDGE_C2_FEM/test_01.hpp b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_WEDGE_C2_FEM/test_01.hpp index 66d8be79c231..8657fb61b52a 100644 --- a/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_WEDGE_C2_FEM/test_01.hpp +++ b/packages/intrepid2/unit-test/Discretization/Basis/HGRAD_WEDGE_C2_FEM/test_01.hpp @@ -41,7 +41,7 @@ // @HEADER /** \file test_01.cpp - \brief Unit tests for the Intrepid2::HGRAD_WEDGE_C2_FEM class. + \brief Unit tests for the Intrepid2::HGRAD_WEDGE_DEG2_FEM class. \author Created by P. Bochev, D. Ridzal, K. Peterson and Kyungjoo Kim. */ @@ -71,8 +71,8 @@ namespace Intrepid2 { *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \ } - template - int HGRAD_WEDGE_C2_FEM_Test01(const bool verbose) { + template + int HGRAD_WEDGE_DEG2_FEM_Test01(const bool verbose) { Teuchos::RCP outStream; Teuchos::oblackholestream bhs; // outputs nothing @@ -94,20 +94,19 @@ namespace Intrepid2 { *outStream << "===============================================================================\n" - << "| |\n" - << "| Unit Test (Basis_HGRAD_WEDGE_C2_FEM) |\n" + << "| |\n"; + + if constexpr (serendipity) + *outStream + << "| Unit Test (Basis_HGRAD_WEDGE_I2_Serendipity FEM) |\n"; + else + *outStream + << "| Unit Test (Basis_HGRAD_WEDGE_C2_FEM) |\n"; + *outStream << "| |\n" << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" << "| 2) Basis values for VALUE, GRAD, and Dk operators |\n" << "| |\n" - << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" - << "| Denis Ridzal (dridzal@sandia.gov), |\n" - << "| Kara Peterson (kjpeter@sandia.gov), |\n" - << "| Kyungjoo Kim (kyukim@sandia.gov). |\n" - << "| |\n" - << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" - << "| Trilinos website: http://trilinos.sandia.gov |\n" - << "| |\n" << "===============================================================================\n"; typedef Kokkos::DynRankView DynRankView; @@ -120,7 +119,11 @@ namespace Intrepid2 { // for virtual function, value and point types are declared in the class typedef ValueType outputValueType; typedef ValueType pointValueType; - Basis_HGRAD_WEDGE_C2_FEM wedgeBasis; + BasisPtr wedgeBasis; + if constexpr (serendipity) + wedgeBasis = Teuchos::rcp(new Basis_HGRAD_WEDGE_I2_Serendipity_FEM()); + else + wedgeBasis = Teuchos::rcp(new Basis_HGRAD_WEDGE_C2_FEM()); *outStream << "\n" @@ -137,83 +140,83 @@ namespace Intrepid2 { DynRankView ConstructWithLabel(wedgeNodes, 18, 3); // Generic array for the output values; needs to be properly resized depending on the operator type - const ordinal_type numFields = wedgeBasis.getCardinality(); + const ordinal_type numFields = wedgeBasis->getCardinality(); const ordinal_type numPoints = wedgeNodes.extent(0); - const ordinal_type spaceDim = wedgeBasis.getBaseCellTopology().getDimension(); + const ordinal_type spaceDim = wedgeBasis->getBaseCellTopology().getDimension(); DynRankView vals("vals", numFields, numPoints); DynRankView vals_vec("vals", numFields, numPoints, spaceDim); { // exception #1: CURL cannot be applied to scalar functions - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getValues(vals_vec, wedgeNodes, OPERATOR_DIV) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getValues(vals_vec, wedgeNodes, OPERATOR_DIV) ); // exception #2: DIV cannot be applied to scalar functions - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getValues(vals_vec, wedgeNodes, OPERATOR_DIV) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getValues(vals_vec, wedgeNodes, OPERATOR_DIV) ); } // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and // getDofTag() to access invalid array elements thereby causing bounds check exception { // exception #3 - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getDofOrdinal(3,0,0) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getDofOrdinal(3,0,0) ); // exception #4 - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getDofOrdinal(1,1,1) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getDofOrdinal(1,1,1) ); // exception #5 - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getDofOrdinal(0,9,0) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getDofOrdinal(0,9,0) ); // exception #6 - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getDofTag(numFields) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getDofTag(numFields) ); // exception #7 - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getDofTag(-1) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getDofTag(-1) ); } // Exceptions 8-18 test exception handling with incorrectly dimensioned input/output arrays { // exception #8: input points array must be of rank-2 DynRankView ConstructWithLabel( badPoints1, 4, 5, 3); - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getValues(vals, badPoints1, OPERATOR_VALUE) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getValues(vals, badPoints1, OPERATOR_VALUE) ); } { // exception #9 dimension 1 in the input point array must equal space dimension of the cell DynRankView ConstructWithLabel( badPoints2, 4, spaceDim + 1); - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getValues(vals, badPoints2, OPERATOR_VALUE) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getValues(vals, badPoints2, OPERATOR_VALUE) ); } { // exception #10 output values must be of rank-2 for OPERATOR_VALUE DynRankView ConstructWithLabel( badVals1, 4, 3, 1); - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getValues(badVals1, wedgeNodes, OPERATOR_VALUE) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getValues(badVals1, wedgeNodes, OPERATOR_VALUE) ); } { // exception #11 output values must be of rank-3 for OPERATOR_GRAD DynRankView ConstructWithLabel( badVals2, 4, 3); - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_GRAD) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getValues(badVals2, wedgeNodes, OPERATOR_GRAD) ); // exception #12 output values must be of rank-3 for OPERATOR_D1 - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D1) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getValues(badVals2, wedgeNodes, OPERATOR_D1) ); // exception #13 output values must be of rank-3 for OPERATOR_D2 - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D2) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getValues(badVals2, wedgeNodes, OPERATOR_D2) ); } { // exception #14 incorrect 0th dimension of output array (must equal number of basis functions) DynRankView ConstructWithLabel( badVals3, numFields + 1, numPoints); - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getValues(badVals3, wedgeNodes, OPERATOR_VALUE) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getValues(badVals3, wedgeNodes, OPERATOR_VALUE) ); } { // exception #15 incorrect 1st dimension of output array (must equal number of points) DynRankView ConstructWithLabel( badVals4, numFields, numPoints + 1); - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getValues(badVals4, wedgeNodes, OPERATOR_VALUE) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getValues(badVals4, wedgeNodes, OPERATOR_VALUE) ); } { // exception #16: incorrect 2nd dimension of output array (must equal the space dimension) DynRankView ConstructWithLabel( badVals5, numFields, numPoints, spaceDim - 1); - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getValues(badVals5, wedgeNodes, OPERATOR_GRAD) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getValues(badVals5, wedgeNodes, OPERATOR_GRAD) ); } { // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 3D) DynRankView ConstructWithLabel( badVals6, numFields, numPoints, 40); - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D2) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getValues(badVals6, wedgeNodes, OPERATOR_D2) ); // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 3D) - INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D3) ); + INTREPID2_TEST_ERROR_EXPECTED( wedgeBasis->getValues(badVals6, wedgeNodes, OPERATOR_D3) ); } #endif if (nthrow != ncatch) { @@ -235,15 +238,15 @@ namespace Intrepid2 { << "===============================================================================\n"; try { - const ordinal_type numFields = wedgeBasis.getCardinality(); - const auto allTags = wedgeBasis.getAllDofTags(); + const ordinal_type numFields = wedgeBasis->getCardinality(); + const auto allTags = wedgeBasis->getAllDofTags(); // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again const ordinal_type dofTagSize = allTags.extent(0); for (ordinal_type i=0;igetDofOrdinal(allTags(i,0), allTags(i,1), allTags(i,2)); - const auto myTag = wedgeBasis.getDofTag(bfOrd); + const auto myTag = wedgeBasis->getDofTag(bfOrd); if( !( (myTag(0) == allTags(i,0)) && (myTag(1) == allTags(i,1)) && (myTag(2) == allTags(i,2)) && @@ -265,9 +268,9 @@ namespace Intrepid2 { // Now do the same but loop over basis functions for (ordinal_type bfOrd=0;bfOrdgetDofTag(bfOrd); - const auto myBfOrd = wedgeBasis.getDofOrdinal(myTag(0), myTag(1), myTag(2)); + const auto myBfOrd = wedgeBasis->getDofOrdinal(myTag(0), myTag(1), myTag(2)); if( bfOrd != myBfOrd) { errorFlag++; *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; @@ -296,6 +299,7 @@ namespace Intrepid2 { outStream -> precision(20); // VALUE: correct basis function values in (F,P) format + const int numC2Fields = 18; const ValueType basisValues[] = { 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, \ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, \ @@ -401,6 +405,7 @@ namespace Intrepid2 { wedgeNodesHost(12,0)= 0.5; wedgeNodesHost(12,1)= 0.0; wedgeNodesHost(12,2)= 1.0; wedgeNodesHost(13,0)= 0.5; wedgeNodesHost(13,1)= 0.5; wedgeNodesHost(13,2)= 1.0; wedgeNodesHost(14,0)= 0.0; wedgeNodesHost(14,1)= 0.5; wedgeNodesHost(14,2)= 1.0; + wedgeNodesHost(15,0)= 0.5; wedgeNodesHost(15,1)= 0.0; wedgeNodesHost(15,2)= 0.0; wedgeNodesHost(16,0)= 0.5; wedgeNodesHost(16,1)= 0.5; wedgeNodesHost(16,2)= 0.0; wedgeNodesHost(17,0)= 0.0; wedgeNodesHost(17,1)= 0.5; wedgeNodesHost(17,2)= 0.0; @@ -409,9 +414,9 @@ namespace Intrepid2 { Kokkos::deep_copy(wedgeNodes, wedgeNodesHost); // Dimensions for the output arrays: - const ordinal_type numFields = wedgeBasis.getCardinality(); + const ordinal_type numFields = wedgeBasis->getCardinality(); const ordinal_type numPoints = wedgeNodes.extent(0); - const ordinal_type spaceDim = wedgeBasis.getBaseCellTopology().getDimension(); + const ordinal_type spaceDim = wedgeBasis->getBaseCellTopology().getDimension(); const ordinal_type D2Cardin = getDkCardinality(OPERATOR_D2, spaceDim); const ordinal_type D3Cardin = getDkCardinality(OPERATOR_D3, spaceDim); const ordinal_type D4Cardin = getDkCardinality(OPERATOR_D4, spaceDim); @@ -419,12 +424,12 @@ namespace Intrepid2 { // Check VALUE of basis functions: resize vals to rank-2 container: { DynRankView vals = DynRankView("vals", numFields, numPoints); - wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_VALUE); + wedgeBasis->getValues(vals, wedgeNodes, OPERATOR_VALUE); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i = 0; i < numFields; ++i) { for (ordinal_type j = 0; j < numPoints; ++j) { - const ordinal_type l = i + j * numFields; + const ordinal_type l = i + j * numC2Fields; if (std::abs(vals_host(i,j) - basisValues[l]) > tol) { errorFlag++; *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; @@ -443,7 +448,7 @@ namespace Intrepid2 { // Check GRAD of basis function: resize vals to rank-3 container { DynRankView vals = DynRankView("vals", numFields, numPoints, spaceDim); - wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_GRAD); + wedgeBasis->getValues(vals, wedgeNodes, OPERATOR_GRAD); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i = 0; i < numFields; ++i) { @@ -470,7 +475,7 @@ namespace Intrepid2 { // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD) { DynRankView vals = DynRankView("vals", numFields, numPoints, spaceDim); - wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D1); + wedgeBasis->getValues(vals, wedgeNodes, OPERATOR_D1); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i = 0; i < numFields; ++i) { @@ -498,7 +503,7 @@ namespace Intrepid2 { // Check D2 of basis function { DynRankView vals = DynRankView("vals", numFields, numPoints, D2Cardin); - wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D2); + wedgeBasis->getValues(vals, wedgeNodes, OPERATOR_D2); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i = 0; i < numFields; ++i) { @@ -525,7 +530,7 @@ namespace Intrepid2 { // Check D3 of basis function { DynRankView vals = DynRankView("vals", numFields, numPoints, D3Cardin); - wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D3); + wedgeBasis->getValues(vals, wedgeNodes, OPERATOR_D3); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i = 0; i < numFields; ++i) { @@ -552,7 +557,7 @@ namespace Intrepid2 { // Check D4 of basis function { DynRankView vals = DynRankView("vals", numFields, numPoints, D4Cardin); - wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D4); + wedgeBasis->getValues(vals, wedgeNodes, OPERATOR_D4); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i = 0; i < numFields; ++i) { @@ -590,7 +595,7 @@ namespace Intrepid2 { const ordinal_type DkCardin = getDkCardinality(op, spaceDim); DynRankView vals("vals", numFields, numPoints, DkCardin); - wedgeBasis.getValues(vals, wedgeNodes, op); + wedgeBasis->getValues(vals, wedgeNodes, op); auto vals_host = Kokkos::create_mirror_view(typename HostSpaceType::memory_space(), vals); Kokkos::deep_copy(vals_host, vals); for (ordinal_type i1=0;i1getFunctionSpace(); if (fs != FUNCTION_SPACE_HGRAD) { diff --git a/packages/intrepid2/unit-test/Projection/CMakeLists.txt b/packages/intrepid2/unit-test/Projection/CMakeLists.txt index 9e3779c12705..5b833fe258e5 100644 --- a/packages/intrepid2/unit-test/Projection/CMakeLists.txt +++ b/packages/intrepid2/unit-test/Projection/CMakeLists.txt @@ -26,7 +26,8 @@ LIST(APPEND Intrepid2_TEST_ETI_FILE "test_interpolation_projection_HEX" "test_interpolation_projection_QUAD" "test_interpolation_projection_TET" - "test_interpolation_projection_TRI") + "test_interpolation_projection_TRI" + "test_interpolation_projection_WEDGE") SET(Intrepid2_TEST_ETI_DEVICE_NAME "") SET(Intrepid2_TEST_ETI_DEVICE "") diff --git a/packages/intrepid2/unit-test/Projection/eti/test_interpolation_projection_HEX_ETI.in b/packages/intrepid2/unit-test/Projection/eti/test_interpolation_projection_HEX_ETI.in index 85c16b9569ca..be0252323a6e 100644 --- a/packages/intrepid2/unit-test/Projection/eti/test_interpolation_projection_HEX_ETI.in +++ b/packages/intrepid2/unit-test/Projection/eti/test_interpolation_projection_HEX_ETI.in @@ -41,7 +41,7 @@ // @HEADER /** \file test_01.cpp - \brief Test for checking orientation tools for tetrahedral elements. + \brief Test for checking orientation tools for hexahedral elements. \author Created by Mauro Perego */ diff --git a/packages/intrepid2/unit-test/Projection/eti/test_interpolation_projection_WEDGE_ETI.in b/packages/intrepid2/unit-test/Projection/eti/test_interpolation_projection_WEDGE_ETI.in new file mode 100644 index 000000000000..3b4b35dbe593 --- /dev/null +++ b/packages/intrepid2/unit-test/Projection/eti/test_interpolation_projection_WEDGE_ETI.in @@ -0,0 +1,62 @@ +// @HEADER +// ************************************************************************ +// +// Intrepid2 Package +// Copyright (2007) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Kyungjoo Kim (kyukim@sandia.gov), or +// Mauro Perego (mperego@sandia.gov) +// +// ************************************************************************ +// @HEADER + +/** \file test_01.cpp + \brief Test for checking orientation tools for wedge elements. + \author Created by Mauro Perego +*/ + +#include "test_interpolation_projection_WEDGE.hpp" +#include "Kokkos_Core.hpp" + + +int main(int argc, char *argv[]) { + + const bool verbose = (argc-1) > 0; + Kokkos::initialize(); + + const int r_val = Intrepid2::Test::InterpolationProjectionWedge<@ETI_VALUETYPE@,@ETI_DEVICE@>(verbose); + + Kokkos::finalize(); + return r_val; +} + diff --git a/packages/intrepid2/unit-test/Projection/test_interpolation_projection_HEX.hpp b/packages/intrepid2/unit-test/Projection/test_interpolation_projection_HEX.hpp index 4cf8c270510b..accc84b2462d 100644 --- a/packages/intrepid2/unit-test/Projection/test_interpolation_projection_HEX.hpp +++ b/packages/intrepid2/unit-test/Projection/test_interpolation_projection_HEX.hpp @@ -45,10 +45,10 @@ \brief Test interpolation and projection capabilities for Hexaedral elements The test considers two hexahedra in the physical space sharing a common face. - In order to test significant configurations, we consider 6 mappings of the reference hexahedron + In order to test significant configurations, we consider 6 rotations of the reference hexahedron to the first (physical) hexahedron, so that the common face is mapped from all the 6 faces of the reference hexahedron. - Then, for each of the mappings, the global ids of the vertices of the common face are permuted. + Then, for each of the mappings, the global ids of the vertices of the common face are permuted (8 permutations). The test considers HGRAD, HCURL, HDIV and HVOL, of different degree, and for each of them checks that the Lagrangian interpolation, the interpolation-based projection, and the L2 projection, reproduce the @@ -151,7 +151,7 @@ int InterpolationProjectionHex(const bool verbose) { pickTest = true; /* initialize random seed: */ std::srand (std::time(NULL)); - int configuration = std::rand() % 24; + int configuration = std::rand() % 48; elemPermutation = configuration % 6; sharedSidePermutation = configuration/6; *outStream << "Randomly picked configuration (cellTopo premutation, shared face permutation): (" << elemPermutation << ", " < >()); ordinal_type numNodesPerElem = cellTopo.getNodeCount(); + using faceShape = shards::Quadrilateral<4>; + const CellTopologyData * const faceTopoData = shards::getCellTopologyData(); + *outStream << "===============================================================================\n" << "| |\n" @@ -260,14 +263,17 @@ int InterpolationProjectionHex(const bool verbose) { //reordering of nodes to explore different orientations - ordinal_type reorder[numTotalVertexes] = {0,1,2,3,4,5,6,7,8,9,10,11}; - - int sharedSideCount = 0; - do { - if((sharedSideCount++ != sharedSidePermutation) && pickTest) + for(int sharedSideCount = 0; sharedSideCountpermutation[sharedSideCount].node[i]]; + } + for(ordinal_type i=0;ipermutation[sharedSideCount].node[i]]; + } + for(ordinal_type i=0;ipermutation[sharedSideCount].node[i]]; + } + for(ordinal_type i=0;i +#include +#include +#include + +namespace Intrepid2 { + +namespace Test { + +#define INTREPID2_TEST_ERROR_EXPECTED( S ) \ + try { \ + ++nthrow; \ + S ; \ + } \ + catch (std::exception &err) { \ + ++ncatch; \ + *outStream << "Expected Error ----------------------------------------------------------------\n"; \ + *outStream << err.what() << '\n'; \ + *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \ + } + +template +int InterpolationProjectionWedge(const bool verbose) { + + using ExecSpaceType = typename DeviceType::execution_space; + using MemSpaceType = typename DeviceType::memory_space; + + using DynRankView = Kokkos::DynRankView; + +#define ConstructWithLabel(obj, ...) obj(#obj, __VA_ARGS__) + + Teuchos::RCP outStream; + Teuchos::oblackholestream bhs; // outputs nothing + + if (verbose) + outStream = Teuchos::rcp(&std::cout, false); + else + outStream = Teuchos::rcp(&bhs, false); + + Teuchos::oblackholestream oldFormatState; + oldFormatState.copyfmt(std::cout); + + using HostSpaceType = Kokkos::DefaultHostExecutionSpace; + + using DynRankViewIntHost = Kokkos::DynRankView; + + *outStream << "DeviceSpace:: "; ExecSpaceType().print_configuration(*outStream, false); + *outStream << "HostSpace:: "; HostSpaceType().print_configuration(*outStream, false); + *outStream << "\n"; + + int errorFlag = 0; + const ValueType tol = tolerence(); + + bool pickTest = false; + int elemPermutation=0, sharedSidePermutation=0; + +#ifdef RANDOMLY_PICK_ELEM_PERMUTATION + pickTest = true; + /* initialize random seed: */ + std::srand (std::time(NULL)); + int configuration = std::rand() % 24; + elemPermutation = configuration % 3; + sharedSidePermutation = configuration/3; + *outStream << "Randomly picked configuration (cellTopo premutation, shared face permutation): (" << elemPermutation << ", " < edgeType; + typedef std::array baseFaceType; + typedef std::array latFaceType; + typedef CellTools ct; + typedef OrientationTools ots; + typedef Experimental::ProjectionTools pts; + typedef FunctionSpaceTools fst; + typedef Experimental::LagrangianInterpolation li; + using basisType = Basis; + + constexpr ordinal_type dim = 3; + constexpr ordinal_type numCells = 2; + constexpr ordinal_type numElemVertexes = 6; + constexpr ordinal_type numTotalVertexes = 8; + + ValueType vertices_orig[numTotalVertexes][dim] = {{0,0,-1},{0,1,-1},{1,0,-1},{0,0,1},{0,1,1},{1,0,1},{1,1,-1},{1,1,1}}; + ordinal_type cells_orig[numCells][numElemVertexes] = {{0,1,2,3,4,5},{1,6,2,4,7,5}}; //common face is {1,2,4,5} + latFaceType common_face = {{1,2,5,4}}; + const baseFaceType bottomFace = {{0,2,1}}; + const baseFaceType topFace = {{3,4,5}}; + //faceType latFace1 = {{0,1,4,3}}; + //faceType latFace3 = {{0,3,5,2}}; + ordinal_type cells_rotated[numCells][numElemVertexes]; + baseFaceType bottomFaceOriented, topFaceOriented;//, latFace1Oriented, latFace3Oriented; + + std::set common_edges; + common_edges.insert(edgeType({{1,2}})); common_edges.insert(edgeType({{4,5}})); common_edges.insert(edgeType({{1,4}})); common_edges.insert(edgeType({{2,5}})); + const ordinal_type max_degree = 4; + + //using CG_NBasis = NodalBasisFamily; + using CG_DNBasis = DerivedNodalBasisFamily; + std::vector basis_set; + + shards::CellTopology cellTopo(shards::getCellTopologyData >()); + using faceShape = shards::Quadrilateral<4>; + const CellTopologyData * const faceTopoData = shards::getCellTopologyData(); + + ordinal_type numNodesPerElem = cellTopo.getNodeCount(); + + *outStream + << "===============================================================================\n" + << "| |\n" + << "| Test 1 (Orientation - HGRAD) |\n" + << "| |\n" + << "===============================================================================\n"; + + + try { + + //reordering of nodes to explore different orientations + + + + for(int sharedSideCount = 0; sharedSideCountpermutation[sharedSideCount].node[i]]; + } + + for(ordinal_type i=0;i elemOrts("elemOrts", numCells); + ots::getOrientation(elemOrts, elemNodes, cellTopo); + + for (ordinal_type degree=1; degree <= max_degree; degree++) { + basis_set.clear(); + if(degree==1) + basis_set.push_back(new Basis_HGRAD_WEDGE_C1_FEM()); + if(degree==2) + basis_set.push_back(new Basis_HGRAD_WEDGE_C2_FEM()); + //basis_set.push_back(new typename CG_NBasis::HGRAD_WEDGE(degree,POINTTYPE_EQUISPACED)); + basis_set.push_back(new typename CG_DNBasis::HGRAD_WEDGE(degree,POINTTYPE_WARPBLEND)); + + for (auto basisPtr:basis_set) { + + auto name = basisPtr->getName(); + *outStream << " " << name << ": " << degree << std::endl; + ordinal_type basisCardinality = basisPtr->getCardinality(); + + //compute DofCoords Oriented + DynRankView ConstructWithLabel(dofCoordsOriented, numCells, basisCardinality, dim); + DynRankView ConstructWithLabel(dofCoeffsPhys, numCells, basisCardinality); + DynRankView ConstructWithLabel(physDofCoords, numCells, basisCardinality, dim); + DynRankView ConstructWithLabel(funAtDofCoords, numCells, basisCardinality); + DynRankView ConstructWithLabel(basisCoeffsLI, numCells, basisCardinality); + + //compute Lagrangian Interpolation of fun + { + li::getDofCoordsAndCoeffs(dofCoordsOriented, dofCoeffsPhys, basisPtr, elemOrts); + + //Compute physical Dof Coordinates + DynRankView ConstructWithLabel(linearBasisValuesAtDofCoord, numCells, numNodesPerElem); + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &i) { + Fun fun; + auto basisValuesAtEvalDofCoord = Kokkos::subview(linearBasisValuesAtDofCoord,i,Kokkos::ALL()); + for(ordinal_type j=0; j::getValues(basisValuesAtEvalDofCoord, evalPoint); + for(ordinal_type k=0; kgetValues(outView, inView); + + // modify basis values to account for orientations + ots::modifyBasisByOrientation(basisValuesAtDofCoordsOriented, + basisValuesAtDofCoords, + elemOrts, + basisPtr); + + // transform basis values + fst::HGRADtransformVALUE(transformedBasisValuesAtDofCoordsOriented, + basisValuesAtDofCoordsOriented); + + + auto hostBasisValues = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), transformedBasisValuesAtDofCoordsOriented); + auto hostDofCoeffsPhys = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), dofCoeffsPhys); + ExecSpaceType().fence(); + for(ordinal_type k=0; k 100*tol ) { + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << " Basis function " << k << " of cell " << i << " does not have unit value at its node (" << dofValue <<")\n"; + } + if ( k!=j && std::abs( dofValue ) > tol ) { + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << " Basis function " << k << " of cell " << i << " does not vanish at node " << j << "(" << dofValue <<")\n"; + } + } + } + } + } + + //check that fun values are consistent on common face dofs + auto hostBasisCoeffsLI = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), basisCoeffsLI); + { + bool areDifferent(false); + auto numFaceDOFs = basisPtr->getDofCount(2,0); + for(ordinal_type j=0;jgetDofOrdinal(2,faceIndex[0],j)) + - hostBasisCoeffsLI(1,basisPtr->getDofOrdinal(2,faceIndex[1],j))) > 10*tol; + } + + if(areDifferent) { + auto hostPhysDofCoords = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), physDofCoords); + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "Function DOFs on common face computed using Wedge 0 basis functions are not consistent with those computed using Wedge 1\n"; + *outStream << "Function DOFs for Wedge 0 are:"; + for(ordinal_type j=0;jgetDofOrdinal(2,faceIndex[0],j)) << " | (" << hostPhysDofCoords(0,basisPtr->getDofOrdinal(2,faceIndex[0],j),0) << "," << hostPhysDofCoords(0,basisPtr->getDofOrdinal(2,faceIndex[0],j),1) << ", " << hostPhysDofCoords(0,basisPtr->getDofOrdinal(2,faceIndex[0],j),2) << ") ||"; + *outStream << "\nFunction DOFs for Wedge 1 are:"; + for(ordinal_type j=0;jgetDofOrdinal(2,faceIndex[1],j))<< " | (" << hostPhysDofCoords(1,basisPtr->getDofOrdinal(2,faceIndex[1],j),0) << "," << hostPhysDofCoords(1,basisPtr->getDofOrdinal(2,faceIndex[1],j),1) << ", " << hostPhysDofCoords(1,basisPtr->getDofOrdinal(2,faceIndex[1],j),2) << ") ||"; + *outStream << std::endl; + } + } + + //check that fun values are consistent on common edges dofs + { + bool areDifferent(false); + auto numEdgeDOFs = basisPtr->getDofCount(1,0); + for(std::size_t iEdge=0;iEdgegetDofOrdinal(1,edgeIndexes[0][iEdge],j)) + - hostBasisCoeffsLI(1,basisPtr->getDofOrdinal(1,edgeIndexes[1][iEdge],j))) > 10*tol; + } + if(areDifferent) + { + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "Function DOFs on common edge " << iEdge << " computed using Wedge 0 basis functions are not consistent with those computed using Wedge 1\n"; + *outStream << "Function DOFs for Wedge 0 are:"; + for(ordinal_type j=0;jgetDofOrdinal(1,edgeIndexes[0][iEdge],j)); + *outStream << "\nFunction DOFs for Wedge 1 are:"; + for(ordinal_type j=0;jgetDofOrdinal(1,edgeIndexes[1][iEdge],j)); + *outStream << std::endl; + } + } + } + + //check that fun values at reference points coincide with those computed using basis functions + DynRankView ConstructWithLabel(basisValuesAtDofCoordsOriented, numCells, basisCardinality, basisCardinality); + DynRankView ConstructWithLabel(transformedBasisValuesAtDofCoordsOriented, numCells, basisCardinality, basisCardinality); + DynRankView basisValuesAtDofCoordsCells("inValues", numCells, basisCardinality, basisCardinality); + + for (ordinal_type ic = 0; ic < numCells; ++ic) + basisPtr->getValues(Kokkos::subview(basisValuesAtDofCoordsCells, ic, Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(dofCoordsOriented, ic, Kokkos::ALL(), Kokkos::ALL())); + + // modify basis values to account for orientations + ots::modifyBasisByOrientation(basisValuesAtDofCoordsOriented, + basisValuesAtDofCoordsCells, + elemOrts, + basisPtr); + + // transform basis (pullback) + fst::HGRADtransformVALUE(transformedBasisValuesAtDofCoordsOriented, + basisValuesAtDofCoordsOriented); + + DynRankView ConstructWithLabel(funAtDofCoordsOriented, numCells, basisCardinality); + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &i) { + for(ordinal_type j=0; j100*tol) { + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "Function values at reference points differ from those computed using basis functions of Wedge " << i << "\n"; + *outStream << "Function values at reference points are:\n"; + for(ordinal_type j=0; jgetDegree()),targetDerivCubDegree(basisPtr->getDegree()); + + Experimental::ProjectionStruct projStruct; + projStruct.createHGradProjectionStruct(basisPtr, targetCubDegree, targetDerivCubDegree); + + ordinal_type numPoints = projStruct.getNumTargetEvalPoints(), numGradPoints = projStruct.getNumTargetDerivEvalPoints(); + DynRankView ConstructWithLabel(evaluationPoints, numCells, numPoints, dim); + DynRankView ConstructWithLabel(evaluationGradPoints, numCells, numGradPoints, dim); + + + pts::getHGradEvaluationPoints(evaluationPoints, + evaluationGradPoints, + elemOrts, + basisPtr, + &projStruct); + + + DynRankView ConstructWithLabel(targetAtEvalPoints, numCells, numPoints); + DynRankView ConstructWithLabel(targetGradAtEvalPoints, numCells, numGradPoints, dim); + + DynRankView ConstructWithLabel(hgradBasisAtEvaluationPoints, numCells, basisCardinality , numPoints); + DynRankView ConstructWithLabel(hgradBasisAtEvaluationPointsNonOriented, numCells, basisCardinality , numPoints); + for(int ic=0; icgetValues(Kokkos::subview(hgradBasisAtEvaluationPointsNonOriented, ic, Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(evaluationPoints, ic, Kokkos::ALL(), Kokkos::ALL()), OPERATOR_VALUE); + ots::modifyBasisByOrientation(hgradBasisAtEvaluationPoints, + hgradBasisAtEvaluationPointsNonOriented, + elemOrts, + basisPtr); + + DynRankView ConstructWithLabel(gradOfHGradBasisAtEvaluationPoints, numCells, basisCardinality , numGradPoints, dim); + if(numGradPoints>0) { + DynRankView ConstructWithLabel(gradOfHGradBasisAtEvaluationPointsNonOriented, numCells, basisCardinality , numGradPoints, dim); + for(int ic=0; icgetValues(Kokkos::subview(gradOfHGradBasisAtEvaluationPointsNonOriented, ic, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(evaluationGradPoints, ic, Kokkos::ALL(), Kokkos::ALL()), OPERATOR_GRAD); + ots::modifyBasisByOrientation(gradOfHGradBasisAtEvaluationPoints, + gradOfHGradBasisAtEvaluationPointsNonOriented, + elemOrts, + basisPtr); + } + + + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &ic) { + for(int i=0;i pow(7, degree-1)*tol) { //heuristic relation on how round-off error depends on degree + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "HGRAD_C" << degree << ": The weights recovered with the optimization are different than the one used for generating the functon."<< + "\nThe max The infinite norm of the difference between the weights is: " << diffErr << std::endl; + } + } + + //compute L2 projection of the Lagrangian interpolation + DynRankView ConstructWithLabel(basisCoeffsL2, numCells, basisCardinality); + { + ordinal_type targetCubDegree(basisPtr->getDegree()); + + Experimental::ProjectionStruct projStruct; + projStruct.createL2ProjectionStruct(basisPtr, targetCubDegree); + + ordinal_type numPoints = projStruct.getNumTargetEvalPoints(); + DynRankView ConstructWithLabel(evaluationPoints, numCells, numPoints, dim); + + pts::getL2EvaluationPoints(evaluationPoints, + elemOrts, + basisPtr, + &projStruct); + + + DynRankView ConstructWithLabel(targetAtEvalPoints, numCells, numPoints); + DynRankView ConstructWithLabel(hgradBasisAtEvaluationPoints, numCells, basisCardinality , numPoints); + DynRankView ConstructWithLabel(hgradBasisAtEvaluationPointsNonOriented, numCells, basisCardinality , numPoints); + for(int ic=0; icgetValues(Kokkos::subview(hgradBasisAtEvaluationPointsNonOriented, ic, Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(evaluationPoints, ic, Kokkos::ALL(), Kokkos::ALL()), OPERATOR_VALUE); + ots::modifyBasisByOrientation(hgradBasisAtEvaluationPoints, + hgradBasisAtEvaluationPointsNonOriented, + elemOrts, + basisPtr); + + + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &ic) { + for(int i=0;i pow(7, degree-1)*tol) { //heuristic relation on how round-off error depends on degree + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "HGRAD_C" << degree << ": The weights recovered with the L2 optimization are different than the one used for generating the functon."<< + "\nThe max The infinite norm of the difference between the weights is: " << diffErr << std::endl; + } + } + + //compute DG L2 projection of the Lagrangian interpolation + DynRankView ConstructWithLabel(basisCoeffsL2DG, numCells, basisCardinality); + { + ordinal_type targetCubDegree(basisPtr->getDegree()); + + Experimental::ProjectionStruct projStruct; + projStruct.createL2DGProjectionStruct(basisPtr, targetCubDegree); + + ordinal_type numPoints = projStruct.getNumTargetEvalPoints(); + DynRankView ConstructWithLabel(evaluationPoints, numCells, numPoints, dim); + + pts::getL2DGEvaluationPoints(evaluationPoints, + basisPtr, + &projStruct); + + + DynRankView ConstructWithLabel(targetAtEvalPoints, numCells, numPoints); + DynRankView ConstructWithLabel(hgradBasisAtEvaluationPoints, numCells, basisCardinality , numPoints); + DynRankView ConstructWithLabel(hgradBasisAtEvaluationPointsNonOriented, numCells, basisCardinality , numPoints); + for(int ic=0; icgetValues(Kokkos::subview(hgradBasisAtEvaluationPointsNonOriented, ic, Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(evaluationPoints, ic, Kokkos::ALL(), Kokkos::ALL()), OPERATOR_VALUE); + ots::modifyBasisByOrientation(hgradBasisAtEvaluationPoints, + hgradBasisAtEvaluationPointsNonOriented, + elemOrts, + basisPtr); + + + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &ic) { + + for(int i=0;i 1e4*tol) { //heuristic relation on how round-off error depends on degree + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "HGRAD_C" << degree << ": The weights recovered with the L2DG optimization are different than the one used for generating the functon."<< + "\nThe max The infinite norm of the difference between the weights is: " << diffErr << std::endl; + } + } + delete basisPtr; + } + } + } + } + } catch (std::exception &err) { + std::cout << " Exeption\n"; + *outStream << err.what() << "\n\n"; + errorFlag = -1000; + } + + + + + *outStream + << "===============================================================================\n" + << "| |\n" + << "| Test 2 (Orientation - HCURL) |\n" + << "| |\n" + << "===============================================================================\n"; + + + try { + + for(int sharedSideCount = 0; sharedSideCountpermutation[sharedSideCount].node[i]]; + } + + for(ordinal_type i=0;i elemOrts("elemOrts", numCells); + ots::getOrientation(elemOrts, elemNodes, cellTopo); + + for (ordinal_type degree=1; degree <= max_degree; degree++) { + + basis_set.clear(); + if(degree==1) + basis_set.push_back(new Basis_HCURL_WEDGE_I1_FEM()); + //basis_set.push_back(new typename CG_NBasis::HCURL_WEDGE(degree,POINTTYPE_WARPBLEND)); + basis_set.push_back(new typename CG_DNBasis::HCURL_WEDGE(degree,POINTTYPE_EQUISPACED)); + + for (auto basisPtr:basis_set) { + + auto name = basisPtr->getName(); + *outStream << " " << name << ": " << degree << std::endl; + + ordinal_type basisCardinality = basisPtr->getCardinality(); + + //compute DofCoords Oriented + DynRankView ConstructWithLabel(dofCoordsOriented, numCells, basisCardinality, dim); + DynRankView ConstructWithLabel(physDofCoords, numCells, basisCardinality, dim); + DynRankView ConstructWithLabel(funAtDofCoords, numCells, basisCardinality, dim); + DynRankView ConstructWithLabel(dofCoeffs, numCells, basisCardinality, dim); + DynRankView ConstructWithLabel(basisCoeffsLI, numCells, basisCardinality); + + //compute Lagrangian Interpolation of fun + { + li::getDofCoordsAndCoeffs(dofCoordsOriented, dofCoeffs, basisPtr, elemOrts); + + //Compute physical Dof Coordinates + + DynRankView ConstructWithLabel(jacobian, numCells, basisCardinality, dim, dim); + ct::setJacobian(jacobian, dofCoordsOriented, physVertexes, cellTopo); + + DynRankView ConstructWithLabel(linearBasisValuesAtDofCoord, numCells, numNodesPerElem); + DynRankView ConstructWithLabel(fwdFunAtDofCoords, numCells, basisCardinality, dim); + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &i) { + FunCurl fun; + auto basisValuesAtEvalDofCoord = Kokkos::subview(linearBasisValuesAtDofCoord,i,Kokkos::ALL()); + for(ordinal_type j=0; j::getValues(basisValuesAtEvalDofCoord, evalPoint); + for(ordinal_type k=0; kgetValues(outView, inView); + + // modify basis values to account for orientations + ots::modifyBasisByOrientation(basisValuesAtDofCoordsOriented, + basisValuesAtDofCoords, + elemOrts, + basisPtr); + + auto hostBasisValues = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), basisValuesAtDofCoordsOriented); + auto hostDofCoeffs = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), dofCoeffs); + for(ordinal_type k=0; k 100*tol ) { + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << " Basis function " << k << " of cell " << i << " does not have unit value at its node (" << dofValue <<")\n"; + } + if ( k!=j && std::abs( dofValue ) > tol ) { + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << " Basis function " << k << " of cell " << i << " does not vanish at node " << j << "(" << dofValue <<")\n"; + } + } + } + } + } + + //check that fun values are consistent on common edges dofs + auto hostBasisCoeffsLI = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), basisCoeffsLI); + { + bool areDifferent(false); + auto numEdgeDOFs = basisPtr->getDofCount(1,0); + for(std::size_t iEdge=0;iEdgegetDofOrdinal(1,edgeIndexes[0][iEdge],j)) + - hostBasisCoeffsLI(1,basisPtr->getDofOrdinal(1,edgeIndexes[1][iEdge],j))) > 10*tol; + } + if(areDifferent) { + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "Function DOFs on common edge " << iEdge << " computed using Wedge 0 basis functions are not consistent with those computed using Wedge 1\n"; + *outStream << "Function DOFs for Wedge 0 are:"; + for(ordinal_type j=0;jgetDofOrdinal(1,edgeIndexes[0][iEdge],j)); + *outStream << "\nFunction DOFs for Wedge 1 are:"; + for(ordinal_type j=0;jgetDofOrdinal(1,edgeIndexes[1][iEdge],j)); + *outStream << std::endl; + } + } + } + + //check that fun values are consistent on common face dofs + { + bool areDifferent(false); + auto numFaceDOFs = basisPtr->getDofCount(2,0); + for(ordinal_type j=0;jgetDofOrdinal(2,faceIndex[0],j)) + - hostBasisCoeffsLI(1,basisPtr->getDofOrdinal(2,faceIndex[1],j))) > 10*tol; + } + + if(areDifferent) { + auto hostPhysDofCoords = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), physDofCoords); + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "Function DOFs on common face computed using Wedge 0 basis functions are not consistent with those computed using Wedge 1\n"; + *outStream << "Function DOFs for Wedge 0 are:"; + for(ordinal_type j=0;jgetDofOrdinal(2,faceIndex[0],j)) << " | (" << hostPhysDofCoords(0,basisPtr->getDofOrdinal(2,faceIndex[0],j),0) << "," << hostPhysDofCoords(0,basisPtr->getDofOrdinal(2,faceIndex[0],j),1) << ", " << hostPhysDofCoords(0,basisPtr->getDofOrdinal(2,faceIndex[0],j),2) << ") ||"; + *outStream << "\nFunction DOFs for Wedge 1 are:"; + for(ordinal_type j=0;jgetDofOrdinal(2,faceIndex[1],j))<< " | (" << hostPhysDofCoords(1,basisPtr->getDofOrdinal(2,faceIndex[1],j),0) << "," << hostPhysDofCoords(1,basisPtr->getDofOrdinal(2,faceIndex[1],j),1) << ", " << hostPhysDofCoords(1,basisPtr->getDofOrdinal(2,faceIndex[1],j),2) << ") ||"; + *outStream << std::endl; + } + } + + //check that fun values at reference points coincide with those computed using basis functions + DynRankView ConstructWithLabel(basisValuesAtDofCoordsOriented, numCells, basisCardinality, basisCardinality, dim); + DynRankView ConstructWithLabel(transformedBasisValuesAtDofCoordsOriented, numCells, basisCardinality, basisCardinality, dim); + DynRankView basisValuesAtDofCoordsCells("inValues", numCells, basisCardinality, basisCardinality, dim); + + for (ordinal_type ic = 0; ic < numCells; ++ic) + basisPtr->getValues(Kokkos::subview(basisValuesAtDofCoordsCells, ic, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(dofCoordsOriented, ic, Kokkos::ALL(), Kokkos::ALL())); + + // modify basis values to account for orientations + ots::modifyBasisByOrientation(basisValuesAtDofCoordsOriented, + basisValuesAtDofCoordsCells, + elemOrts, + basisPtr); + + // transform basis values + DynRankView ConstructWithLabel(jacobianAtDofCoords, numCells, basisCardinality, dim, dim); + DynRankView ConstructWithLabel(jacobianAtDofCoords_inv, numCells, basisCardinality, dim, dim); + ct::setJacobian(jacobianAtDofCoords, dofCoordsOriented, physVertexes, cellTopo); + ct::setJacobianInv (jacobianAtDofCoords_inv, jacobianAtDofCoords); + fst::HCURLtransformVALUE(transformedBasisValuesAtDofCoordsOriented, + jacobianAtDofCoords_inv, + basisValuesAtDofCoordsOriented); + DynRankView ConstructWithLabel(funAtDofCoordsOriented, numCells, basisCardinality, dim); + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &i) { + for(ordinal_type j=0; j100*tol) { + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "Function values at reference points differ from those computed using basis functions of Wedge " << i << "\n"; + *outStream << "Function values at reference points are:\n"; + for(ordinal_type j=0; jgetDegree()),targetDerivCubDegree(basisPtr->getDegree()-1); + + Experimental::ProjectionStruct projStruct; + projStruct.createHCurlProjectionStruct(basisPtr, targetCubDegree, targetDerivCubDegree); + + ordinal_type numPoints = projStruct.getNumTargetEvalPoints(), numCurlPoints = projStruct.getNumTargetDerivEvalPoints(); + DynRankView ConstructWithLabel(evaluationPoints, numCells, numPoints, dim); + DynRankView ConstructWithLabel(evaluationCurlPoints, numCells, numCurlPoints, dim); + + pts::getHCurlEvaluationPoints(evaluationPoints, + evaluationCurlPoints, + elemOrts, + basisPtr, + &projStruct); + + + DynRankView ConstructWithLabel(targetAtEvalPoints, numCells, numPoints, dim); + DynRankView ConstructWithLabel(targetCurlAtEvalPoints, numCells, numCurlPoints, dim); + + DynRankView ConstructWithLabel(hcurlBasisAtEvaluationPoints, numCells, basisCardinality , numPoints, dim); + DynRankView ConstructWithLabel(hcurlBasisAtEvaluationPointsNonOriented, numCells, basisCardinality , numPoints, dim); + for(int ic=0; icgetValues(Kokkos::subview(hcurlBasisAtEvaluationPointsNonOriented, ic, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(evaluationPoints, ic, Kokkos::ALL(), Kokkos::ALL()), OPERATOR_VALUE); + ots::modifyBasisByOrientation(hcurlBasisAtEvaluationPoints, + hcurlBasisAtEvaluationPointsNonOriented, + elemOrts, + basisPtr); + + DynRankView ConstructWithLabel(curlOfHCurlBasisAtEvaluationPoints, numCells, basisCardinality , numCurlPoints, dim); + if(numCurlPoints>0) { + DynRankView ConstructWithLabel(curlOfHCurlBasisAtEvaluationPointsNonOriented, numCells, basisCardinality , numCurlPoints, dim); + for(int ic=0; icgetValues(Kokkos::subview(curlOfHCurlBasisAtEvaluationPointsNonOriented, ic, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(evaluationCurlPoints, ic, Kokkos::ALL(), Kokkos::ALL()), OPERATOR_CURL); + ots::modifyBasisByOrientation(curlOfHCurlBasisAtEvaluationPoints, + curlOfHCurlBasisAtEvaluationPointsNonOriented, + elemOrts, + basisPtr); + } + + + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &ic) { + for(int i=0;i pow(10, degree-1)*tol) { //heuristic relation on how round-off error depends on degree + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "HCURL_I" << degree << ": The weights recovered with the optimization are different than the one used for generating the functon."<< + "\nThe max The infinite norm of the difference between the weights is: " << diffErr << " " <getDegree()); + + Experimental::ProjectionStruct projStruct; + projStruct.createL2ProjectionStruct(basisPtr, targetCubDegree); + + ordinal_type numPoints = projStruct.getNumTargetEvalPoints(); + DynRankView ConstructWithLabel(evaluationPoints, numCells, numPoints, dim); + + pts::getL2EvaluationPoints(evaluationPoints, + elemOrts, + basisPtr, + &projStruct); + + + DynRankView ConstructWithLabel(targetAtEvalPoints, numCells, numPoints, dim); + + DynRankView ConstructWithLabel(hcurlBasisAtEvaluationPoints, numCells, basisCardinality , numPoints, dim); + DynRankView ConstructWithLabel(hcurlBasisAtEvaluationPointsNonOriented, numCells, basisCardinality , numPoints, dim); + for(int ic=0; icgetValues(Kokkos::subview(hcurlBasisAtEvaluationPointsNonOriented, ic, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(evaluationPoints, ic, Kokkos::ALL(), Kokkos::ALL()), OPERATOR_VALUE); + ots::modifyBasisByOrientation(hcurlBasisAtEvaluationPoints, + hcurlBasisAtEvaluationPointsNonOriented, + elemOrts, + basisPtr); + + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &ic) { + for(int i=0;i pow(7, degree-1)*tol) { //heuristic relation on how round-off error depends on degree + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "HCURL_I" << degree << ": The weights recovered with the L2 optimization are different than the one used for generating the functon."<< + "\nThe max The infinite norm of the difference between the weights is: " << diffErr << std::endl; + } + } + + //compute L2DG projection of the Lagrangian interpolation + DynRankView ConstructWithLabel(basisCoeffsL2DG, numCells, basisCardinality); + { + ordinal_type targetCubDegree(basisPtr->getDegree()); + + Experimental::ProjectionStruct projStruct; + projStruct.createL2DGProjectionStruct(basisPtr, targetCubDegree); + + ordinal_type numPoints = projStruct.getNumTargetEvalPoints(); + DynRankView ConstructWithLabel(evaluationPoints, numCells, numPoints, dim); + + pts::getL2DGEvaluationPoints(evaluationPoints, + basisPtr, + &projStruct); + + + DynRankView ConstructWithLabel(targetAtEvalPoints, numCells, numPoints, dim); + + DynRankView ConstructWithLabel(hcurlBasisAtEvaluationPoints, numCells, basisCardinality , numPoints, dim); + DynRankView ConstructWithLabel(hcurlBasisAtEvaluationPointsNonOriented, numCells, basisCardinality , numPoints, dim); + for(int ic=0; icgetValues(Kokkos::subview(hcurlBasisAtEvaluationPointsNonOriented, ic, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(evaluationPoints, ic, Kokkos::ALL(), Kokkos::ALL()), OPERATOR_VALUE); + ots::modifyBasisByOrientation(hcurlBasisAtEvaluationPoints, + hcurlBasisAtEvaluationPointsNonOriented, + elemOrts, + basisPtr); + + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &ic) { + for(int i=0;i 1e4*tol) { + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "HCURL_I" << degree << ": The weights recovered with the L2DG optimization are different than the one used for generating the functon."<< + "\nThe max The infinite norm of the difference between the weights is: " << diffErr << std::endl; + } + } + + delete basisPtr; + } + } + } + } //reorder vertices of common face + + } catch (std::exception &err) { + std::cout << " Exeption\n"; + *outStream << err.what() << "\n\n"; + errorFlag = -1000; + } + + + + *outStream + << "===============================================================================\n" + << "| |\n" + << "| Test 3 (Orientation - HDIV) |\n" + << "| |\n" + << "===============================================================================\n"; + + + try { + + for(int sharedSideCount = 0; sharedSideCountpermutation[sharedSideCount].node[i]]; + } + + for(ordinal_type i=0;i elemOrts("elemOrts", numCells); + ots::getOrientation(elemOrts, elemNodes, cellTopo); + + for (ordinal_type degree=1; degree <= max_degree; degree++) { + + basis_set.clear(); + if(degree==1) + basis_set.push_back(new Basis_HDIV_WEDGE_I1_FEM()); + //basis_set.push_back(new typename CG_NBasis::HDIV_WEDGE(degree,POINTTYPE_EQUISPACED)); + basis_set.push_back(new typename CG_DNBasis::HDIV_WEDGE(degree,POINTTYPE_WARPBLEND)); + + for (auto basisPtr:basis_set) { + + auto name = basisPtr->getName(); + *outStream << " " << name << ": "<< degree << std::endl; + ordinal_type basisCardinality = basisPtr->getCardinality(); + + //compute DofCoords Oriented + + DynRankView ConstructWithLabel(dofCoordsOriented, numCells, basisCardinality, dim); + DynRankView ConstructWithLabel(dofCoeffs, numCells, basisCardinality, dim); + DynRankView ConstructWithLabel(physDofCoords, numCells, basisCardinality, dim); + DynRankView ConstructWithLabel(funAtDofCoords, numCells, basisCardinality, dim); + DynRankView ConstructWithLabel(basisCoeffsLI, numCells, basisCardinality); + + //compute Lagrangian Interpolation of fun + { + + li::getDofCoordsAndCoeffs(dofCoordsOriented, dofCoeffs, basisPtr, elemOrts); + //need to transform dofCoeff to physical space (they transform as normals) + DynRankView ConstructWithLabel(jacobian, numCells, basisCardinality, dim, dim); + DynRankView ConstructWithLabel(jacobian_inv, numCells, basisCardinality, dim, dim); + DynRankView ConstructWithLabel(jacobian_det, numCells, basisCardinality); + ct::setJacobian(jacobian, dofCoordsOriented, physVertexes, cellTopo); + ct::setJacobianInv (jacobian_inv, jacobian); + ct::setJacobianDet (jacobian_det, jacobian); + + //Compute physical Dof Coordinates + DynRankView ConstructWithLabel(linearBasisValuesAtDofCoord, numCells, numNodesPerElem); + DynRankView ConstructWithLabel(fwdFunAtDofCoords, numCells, basisCardinality, dim); + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &i) { + FunDiv fun; + auto basisValuesAtEvalDofCoord = Kokkos::subview(linearBasisValuesAtDofCoord,i,Kokkos::ALL()); + for(ordinal_type j=0; j::getValues(basisValuesAtEvalDofCoord, evalPoint); + for(ordinal_type k=0; kgetValues(outView, inView); + + // modify basis values to account for orientations + ots::modifyBasisByOrientation(basisValuesAtDofCoordsOriented, + basisValuesAtDofCoords, + elemOrts, + basisPtr); + + DynRankView ConstructWithLabel(jacobian, numCells, basisCardinality, dim, dim); + DynRankView ConstructWithLabel(jacobian_det, numCells, basisCardinality); + ct::setJacobian(jacobian, dofCoordsOriented, physVertexes, cellTopo); + ct::setJacobianDet (jacobian_det, jacobian); + + auto hostBasisValues = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), basisValuesAtDofCoordsOriented); + auto hostDofCoeffs = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), dofCoeffs); + for(ordinal_type k=0; k 100*tol ) { + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << " Basis function " << k << " of cell " << i << " does not have unit value at its node (" << dofValue <<")\n"; + } + if ( k!=j && std::abs( dofValue ) > tol ) { + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << " Basis function " << k << " of cell " << i << " does not vanish at node " << j << "(" << dofValue <<")\n"; + } + } + } + } + } + + //check that fun values are consistent on common face dofs + auto hostBasisCoeffsLI = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), basisCoeffsLI); + { + bool areDifferent(false); + auto numFaceDOFs = basisPtr->getDofCount(2,0); + for(ordinal_type j=0;jgetDofOrdinal(2,faceIndex[0],j)) + - hostBasisCoeffsLI(1,basisPtr->getDofOrdinal(2,faceIndex[1],j))) > 10*tol; + } + + if(areDifferent) { + errorFlag++; + auto hostPhysDofCoords = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), physDofCoords); + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "Function DOFs on common face computed using Wedge 0 basis functions are not consistent with those computed using Wedge 1\n"; + *outStream << "Function DOFs for Wedge 0 are:"; + for(ordinal_type j=0;jgetDofOrdinal(2,faceIndex[0],j)) << " | (" << hostPhysDofCoords(0,basisPtr->getDofOrdinal(2,faceIndex[0],j),0) << "," << hostPhysDofCoords(0,basisPtr->getDofOrdinal(2,faceIndex[0],j),1) << ", " << hostPhysDofCoords(0,basisPtr->getDofOrdinal(2,faceIndex[0],j),2) << ") ||"; + *outStream << "\nFunction DOFs for Wedge 1 are:"; + for(ordinal_type j=0;jgetDofOrdinal(2,faceIndex[1],j))<< " | (" << hostPhysDofCoords(1,basisPtr->getDofOrdinal(2,faceIndex[1],j),0) << "," << hostPhysDofCoords(1,basisPtr->getDofOrdinal(2,faceIndex[1],j),1) << ", " << hostPhysDofCoords(1,basisPtr->getDofOrdinal(2,faceIndex[1],j),2) << ") ||"; + *outStream << std::endl; + } + } + + //check that fun values at reference points coincide with those computed using basis functions + DynRankView ConstructWithLabel(basisValuesAtDofCoordsOriented, numCells, basisCardinality, basisCardinality, dim); + DynRankView ConstructWithLabel(transformedBasisValuesAtDofCoordsOriented, numCells, basisCardinality, basisCardinality, dim); + DynRankView basisValuesAtDofCoordsCells("inValues", numCells, basisCardinality, basisCardinality, dim); + + for (ordinal_type ic = 0; ic < numCells; ++ic) + basisPtr->getValues(Kokkos::subview(basisValuesAtDofCoordsCells, ic, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(dofCoordsOriented, ic, Kokkos::ALL(), Kokkos::ALL())); + + // modify basis values to account for orientations + ots::modifyBasisByOrientation(basisValuesAtDofCoordsOriented, + basisValuesAtDofCoordsCells, + elemOrts, + basisPtr); + + // transform basis values + DynRankView ConstructWithLabel(jacobianAtDofCoords, numCells, basisCardinality, dim, dim); + DynRankView ConstructWithLabel(jacobianAtDofCoords_det, numCells, basisCardinality); + ct::setJacobian(jacobianAtDofCoords, dofCoordsOriented, physVertexes, cellTopo); + ct::setJacobianDet (jacobianAtDofCoords_det, jacobianAtDofCoords); + fst::HDIVtransformVALUE(transformedBasisValuesAtDofCoordsOriented, + jacobianAtDofCoords, + jacobianAtDofCoords_det, + basisValuesAtDofCoordsOriented); + DynRankView ConstructWithLabel(funAtDofCoordsOriented, numCells, basisCardinality, dim); + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &i) { + for(ordinal_type j=0; j100*tol) { + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "Function values at reference points differ from those computed using basis functions of Wedge " << i << "\n"; + *outStream << "Function values at reference points are:\n"; + for(ordinal_type j=0; jgetDegree()),targetDerivCubDegree(basisPtr->getDegree()-1); + + Experimental::ProjectionStruct projStruct; + projStruct.createHDivProjectionStruct(basisPtr, targetCubDegree, targetDerivCubDegree); + + ordinal_type numPoints = projStruct.getNumTargetEvalPoints(), numDivPoints = projStruct.getNumTargetDerivEvalPoints(); + + DynRankView ConstructWithLabel(evaluationPoints, numCells, numPoints, dim); + DynRankView ConstructWithLabel(evaluationDivPoints, numCells, numDivPoints, dim); + + pts::getHDivEvaluationPoints(evaluationPoints, + evaluationDivPoints, + elemOrts, + basisPtr, + &projStruct); + + DynRankView ConstructWithLabel(targetAtEvalPoints, numCells, numPoints, dim); + DynRankView ConstructWithLabel(targetDivAtEvalPoints, numCells, numDivPoints); + + DynRankView ConstructWithLabel(hdivBasisAtEvaluationPoints, numCells, basisCardinality , numPoints, dim); + DynRankView ConstructWithLabel(hdivBasisAtEvaluationPointsNonOriented, numCells, basisCardinality , numPoints, dim); + for(ordinal_type ic=0; icgetValues(Kokkos::subview(hdivBasisAtEvaluationPointsNonOriented, ic, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(evaluationPoints, ic, Kokkos::ALL(), Kokkos::ALL()), OPERATOR_VALUE); + ots::modifyBasisByOrientation(hdivBasisAtEvaluationPoints, + hdivBasisAtEvaluationPointsNonOriented, + elemOrts, + basisPtr); + + DynRankView ConstructWithLabel(divOfHDivBasisAtEvaluationPoints, numCells, basisCardinality , numDivPoints); + if(numDivPoints>0) { + DynRankView ConstructWithLabel(divOfHDivBasisAtEvaluationPointsNonOriented, numCells, basisCardinality , numDivPoints); + for(ordinal_type ic=0; icgetValues(Kokkos::subview(divOfHDivBasisAtEvaluationPointsNonOriented, ic, Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(evaluationDivPoints, ic, Kokkos::ALL(), Kokkos::ALL()), OPERATOR_DIV); + ots::modifyBasisByOrientation(divOfHDivBasisAtEvaluationPoints, + divOfHDivBasisAtEvaluationPointsNonOriented, + elemOrts, + basisPtr); + } + + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &ic) { + for(int i=0;i pow(8, degree)*tol) { //heuristic relation on how round-off error depends on degree + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "HDIV_I" << degree << ": The weights recovered with the optimization are different than the one used for generating the functon."<< + "\nThe max The infinite norm of the difference between the weights is: " << diffErr << " (tol: " << pow(8, degree)*tol << ")" << std::endl; + } + } + + //compute L2 projection interpolation of the Lagrangian interpolation + DynRankView ConstructWithLabel(basisCoeffsL2, numCells, basisCardinality); + { + ordinal_type targetCubDegree(basisPtr->getDegree()); + + Experimental::ProjectionStruct projStruct; + projStruct.createL2ProjectionStruct(basisPtr, targetCubDegree); + + ordinal_type numPoints = projStruct.getNumTargetEvalPoints(); + + DynRankView ConstructWithLabel(evaluationPoints, numCells, numPoints, dim); + + pts::getL2EvaluationPoints(evaluationPoints, + elemOrts, + basisPtr, + &projStruct); + + DynRankView ConstructWithLabel(targetAtEvalPoints, numCells, numPoints, dim); + + DynRankView ConstructWithLabel(hdivBasisAtEvaluationPoints, numCells, basisCardinality , numPoints, dim); + DynRankView ConstructWithLabel(hdivBasisAtEvaluationPointsNonOriented, numCells, basisCardinality , numPoints, dim); + for(ordinal_type ic=0; icgetValues(Kokkos::subview(hdivBasisAtEvaluationPointsNonOriented, ic, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(evaluationPoints, ic, Kokkos::ALL(), Kokkos::ALL()), OPERATOR_VALUE); + ots::modifyBasisByOrientation(hdivBasisAtEvaluationPoints, + hdivBasisAtEvaluationPointsNonOriented, + elemOrts, + basisPtr); + + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &ic) { + for(int i=0;i pow(7, degree-1)*tol) { //heuristic relation on how round-off error depends on degree + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "HDIV_I" << degree << ": The weights recovered with the L2 optimization are different than the one used for generating the function."<< + "\nThe max The infinite norm of the difference between the weights is: " << diffErr << std::endl; + } + } + + //compute DG L2 projection interpolation of the Lagrangian interpolation + DynRankView ConstructWithLabel(basisCoeffsL2DG, numCells, basisCardinality); + { + ordinal_type targetCubDegree(basisPtr->getDegree()); + + Experimental::ProjectionStruct projStruct; + projStruct.createL2DGProjectionStruct(basisPtr, targetCubDegree); + + ordinal_type numPoints = projStruct.getNumTargetEvalPoints(); + + DynRankView ConstructWithLabel(evaluationPoints, numCells, numPoints, dim); + + pts::getL2DGEvaluationPoints(evaluationPoints, + basisPtr, + &projStruct); + + DynRankView ConstructWithLabel(targetAtEvalPoints, numCells, numPoints, dim); + + DynRankView ConstructWithLabel(hdivBasisAtEvaluationPoints, numCells, basisCardinality , numPoints, dim); + DynRankView ConstructWithLabel(hdivBasisAtEvaluationPointsNonOriented, numCells, basisCardinality , numPoints, dim); + for(ordinal_type ic=0; icgetValues(Kokkos::subview(hdivBasisAtEvaluationPointsNonOriented, ic, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(evaluationPoints, ic, Kokkos::ALL(), Kokkos::ALL()), OPERATOR_VALUE); + ots::modifyBasisByOrientation(hdivBasisAtEvaluationPoints, + hdivBasisAtEvaluationPointsNonOriented, + elemOrts, + basisPtr); + + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &ic) { + for(int i=0;i pow(7, degree-1)*tol) { //heuristic relation on how round-off error depends on degree + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "HDIV_I" << degree << ": The weights recovered with the L2DG optimization are different than the one used for generating the function."<< + "\nThe max The infinite norm of the difference between the weights is: " << diffErr << std::endl; + } + } + + delete basisPtr; + } + } + } + } //reorder vertices of common face + + } catch (std::exception &err) { + std::cout << " Exeption\n"; + *outStream << err.what() << "\n\n"; + errorFlag = -1000; + } + + *outStream + << "===============================================================================\n" + << "| |\n" + << "| Test 4 (Orientation - HVOL) |\n" + << "| |\n" + << "===============================================================================\n"; + + + try { + + + ValueType vertices[numTotalVertexes][dim]; + ordinal_type cells[numCells][numElemVertexes]; + + for(ordinal_type i=0; i elemOrts("elemOrts", numCells); + ots::getOrientation(elemOrts, elemNodes, cellTopo); + + for (ordinal_type degree=1; degree <= max_degree; degree++) { + + basis_set.clear(); + if(degree==1) + basis_set.push_back(new Basis_HVOL_C0_FEM(cellTopo)); + //basis_set.push_back(new typename CG_NBasis::HVOL_WEDGE(degree,POINTTYPE_EQUISPACED)); + basis_set.push_back(new typename CG_DNBasis::HVOL_WEDGE(degree,POINTTYPE_WARPBLEND)); + + for (auto basisPtr:basis_set) { + + auto name = basisPtr->getName(); + *outStream << " " << name << ": "<< degree << std::endl; + + ordinal_type basisCardinality = basisPtr->getCardinality(); + + //compute DofCoords Oriented + DynRankView ConstructWithLabel(dofCoordsOriented, numCells, basisCardinality, dim); + DynRankView ConstructWithLabel(dofCoeffsPhys, numCells, basisCardinality); + DynRankView ConstructWithLabel(physDofCoords, numCells, basisCardinality, dim); + DynRankView ConstructWithLabel(funAtDofCoords, numCells, basisCardinality); + DynRankView ConstructWithLabel(basisCoeffsLI, numCells, basisCardinality); + + //compute Lagrangian Interpolation of fun + { + li::getDofCoordsAndCoeffs(dofCoordsOriented, dofCoeffsPhys, basisPtr, elemOrts); + + //need to transform dofCoeff to physical space (they transform as normals) + DynRankView ConstructWithLabel(jacobian, numCells, basisCardinality, dim, dim); + DynRankView ConstructWithLabel(jacobian_det, numCells, basisCardinality); + ct::setJacobian(jacobian, dofCoordsOriented, physVertexes, cellTopo); + ct::setJacobianDet (jacobian_det, jacobian); + + //Compute physical Dof Coordinates + DynRankView ConstructWithLabel(linearBasisValuesAtDofCoord, numCells, numNodesPerElem); + DynRankView ConstructWithLabel(fwdFunAtDofCoords, numCells, basisCardinality); + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &i) { + Fun fun; + auto basisValuesAtEvalDofCoord = Kokkos::subview(linearBasisValuesAtDofCoord,i,Kokkos::ALL()); + for(ordinal_type j=0; j::getValues(basisValuesAtEvalDofCoord, evalPoint); + for(ordinal_type k=0; kgetValues(outView, inView); + + // modify basis values to account for orientations + ots::modifyBasisByOrientation(basisValuesAtDofCoordsOriented, + basisValuesAtDofCoords, + elemOrts, + basisPtr); + + auto hostBasisValues = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), basisValuesAtDofCoordsOriented); + auto hostDofCoeffsPhys = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), dofCoeffsPhys); + for(ordinal_type k=0; k 100*tol ) { + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << " Basis function " << k << " of cell " << i << " does not have unit value at its node (" << dofValue <<")\n"; + } + if ( k!=j && std::abs( dofValue ) > tol ) { + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << " Basis function " << k << " of cell " << i << " does not vanish at node " << j << "(" << dofValue <<")\n"; + } + } + } + } + } + + //check that fun values at reference points coincide with those computed using basis functions + DynRankView ConstructWithLabel(basisValuesAtDofCoordsOriented, numCells, basisCardinality, basisCardinality); + DynRankView ConstructWithLabel(transformedBasisValuesAtDofCoordsOriented, numCells, basisCardinality, basisCardinality); + DynRankView basisValuesAtDofCoordsCells("inValues", numCells, basisCardinality, basisCardinality); + + for (ordinal_type ic = 0; ic < numCells; ++ic) + basisPtr->getValues(Kokkos::subview(basisValuesAtDofCoordsCells, ic, Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(dofCoordsOriented, ic, Kokkos::ALL(), Kokkos::ALL())); + + // modify basis values to account for orientations + ots::modifyBasisByOrientation(basisValuesAtDofCoordsOriented, + basisValuesAtDofCoordsCells, + elemOrts, + basisPtr); + + // transform basis values + DynRankView ConstructWithLabel(jacobian, numCells, basisCardinality, dim, dim); + DynRankView ConstructWithLabel(jacobian_det, numCells, basisCardinality); + ct::setJacobian(jacobian, dofCoordsOriented, physVertexes, cellTopo); + ct::setJacobianDet (jacobian_det, jacobian); + fst::HVOLtransformVALUE(transformedBasisValuesAtDofCoordsOriented, + jacobian_det, + basisValuesAtDofCoordsOriented); + + DynRankView ConstructWithLabel(funAtDofCoordsOriented, numCells, basisCardinality); + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &i) { + for(ordinal_type j=0; j100*tol) { + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "Function values at reference points differ from those computed using basis functions of Wedge " << i << "\n"; + *outStream << "Function values at reference points are:\n"; + for(ordinal_type j=0; jgetDegree()); + + Experimental::ProjectionStruct projStruct; + projStruct.createHVolProjectionStruct(basisPtr, targetCubDegree); + + ordinal_type numPoints = projStruct.getNumTargetEvalPoints(); + DynRankView ConstructWithLabel(evaluationPoints, numCells, numPoints, dim); + + + pts::getHVolEvaluationPoints(evaluationPoints, + elemOrts, + basisPtr, + &projStruct); + + + DynRankView ConstructWithLabel(targetAtEvalPoints, numCells, numPoints); + + DynRankView ConstructWithLabel(hvolBasisAtEvaluationPoints, numCells, basisCardinality , numPoints); + DynRankView ConstructWithLabel(hvolBasisAtEvaluationPointsNonOriented, numCells, basisCardinality , numPoints); + for(int ic=0; icgetValues(Kokkos::subview(hvolBasisAtEvaluationPointsNonOriented, ic, Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(evaluationPoints, ic, Kokkos::ALL(), Kokkos::ALL()), OPERATOR_VALUE); + ots::modifyBasisByOrientation(hvolBasisAtEvaluationPoints, + hvolBasisAtEvaluationPointsNonOriented, + elemOrts, + basisPtr); + + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &ic) { + for(int i=0;i pow(20, degree)*tol) { //heuristic relation on how round-off error depends on degree + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "HVOL_C" << degree << ": The weights recovered with the optimization are different than the one used for generating the functon."<< + "\nThe max The infinite norm of the difference between the weights is: " << diffErr << std::endl; + } + } + + //compute L2 projection of the Lagrangian interpolation + DynRankView ConstructWithLabel(basisCoeffsL2, numCells, basisCardinality); + { + ordinal_type targetCubDegree(basisPtr->getDegree()); + + Experimental::ProjectionStruct projStruct; + projStruct.createL2ProjectionStruct(basisPtr, targetCubDegree); + + ordinal_type numPoints = projStruct.getNumTargetEvalPoints(); + DynRankView ConstructWithLabel(evaluationPoints, numCells, numPoints, dim); + + + pts::getL2EvaluationPoints(evaluationPoints, + elemOrts, + basisPtr, + &projStruct); + + + DynRankView ConstructWithLabel(targetAtEvalPoints, numCells, numPoints); + + DynRankView ConstructWithLabel(hvolBasisAtEvaluationPoints, numCells, basisCardinality , numPoints); + DynRankView ConstructWithLabel(hvolBasisAtEvaluationPointsNonOriented, numCells, basisCardinality , numPoints); + for(int ic=0; icgetValues(Kokkos::subview(hvolBasisAtEvaluationPointsNonOriented, ic, Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(evaluationPoints, ic, Kokkos::ALL(), Kokkos::ALL()), OPERATOR_VALUE); + ots::modifyBasisByOrientation(hvolBasisAtEvaluationPoints, + hvolBasisAtEvaluationPointsNonOriented, + elemOrts, + basisPtr); + + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &ic) { + for(int i=0;i pow(20, degree)*tol) { //heuristic relation on how round-off error depends on degree + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "HVOL_C" << degree << ": The weights recovered with the L2 optimization are different than the one used for generating the functon."<< + "\nThe max The infinite norm of the difference between the weights is: " << diffErr << std::endl; + } + } + + //compute DG L2 projection of the Lagrangian interpolation + DynRankView ConstructWithLabel(basisCoeffsL2DG, numCells, basisCardinality); + { + ordinal_type targetCubDegree(basisPtr->getDegree()); + + Experimental::ProjectionStruct projStruct; + projStruct.createL2DGProjectionStruct(basisPtr, targetCubDegree); + + ordinal_type numPoints = projStruct.getNumTargetEvalPoints(); + DynRankView ConstructWithLabel(evaluationPoints, numCells, numPoints, dim); + + + pts::getL2DGEvaluationPoints(evaluationPoints, + basisPtr, + &projStruct); + + + DynRankView ConstructWithLabel(targetAtEvalPoints, numCells, numPoints); + + DynRankView ConstructWithLabel(hvolBasisAtEvaluationPoints, numCells, basisCardinality , numPoints); + DynRankView ConstructWithLabel(hvolBasisAtEvaluationPointsNonOriented, numCells, basisCardinality , numPoints); + for(int ic=0; icgetValues(Kokkos::subview(hvolBasisAtEvaluationPointsNonOriented, ic, Kokkos::ALL(), Kokkos::ALL()), Kokkos::subview(evaluationPoints, ic, Kokkos::ALL(), Kokkos::ALL()), OPERATOR_VALUE); + ots::modifyBasisByOrientation(hvolBasisAtEvaluationPoints, + hvolBasisAtEvaluationPointsNonOriented, + elemOrts, + basisPtr); + + Kokkos::parallel_for(Kokkos::RangePolicy(0,numCells), + KOKKOS_LAMBDA (const int &ic) { + for(int i=0;i pow(20, degree)*tol) { //heuristic relation on how round-off error depends on degree + errorFlag++; + *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; + *outStream << "HVOL_C" << degree << ": The weights recovered with the L2DG optimization are different than the one used for generating the functon."<< + "\nThe max The infinite norm of the difference between the weights is: " << diffErr << std::endl; + } + } + + delete basisPtr; + } + } + } catch (std::exception &err) { + std::cout << " Exeption\n"; + *outStream << err.what() << "\n\n"; + errorFlag = -1000; + } + + + if (errorFlag != 0) + std::cout << "End Result: TEST FAILED = " << errorFlag << "\n"; + else + std::cout << "End Result: TEST PASSED\n"; + + // reset format state of std::cout + std::cout.copyfmt(oldFormatState); + + return errorFlag; +} +} +} + diff --git a/packages/sacado/src/KokkosExp_View_Fad.hpp b/packages/sacado/src/KokkosExp_View_Fad.hpp index 390581d07b26..b7d34f467876 100644 --- a/packages/sacado/src/KokkosExp_View_Fad.hpp +++ b/packages/sacado/src/KokkosExp_View_Fad.hpp @@ -1521,10 +1521,11 @@ class ViewMapping< Traits , /* View internal mapping */ m_impl_handle = handle_type( reinterpret_cast< pointer_type >( record->data() ) ); if ( ctor_prop::initialize ) { + auto space = ((ViewCtorProp const &) prop).value; // Assume destruction is only required when construction is requested. // The ViewValueFunctor has both value construction and destruction operators. if (execution_space_specified) - record->m_destroy = functor_type( ( (ViewCtorProp const &) prop).value + record->m_destroy = functor_type( space , (fad_value_type *) m_impl_handle , m_array_offset.span() , record->get_label() @@ -1537,6 +1538,7 @@ class ViewMapping< Traits , /* View internal mapping */ // Construct values record->m_destroy.construct_shared_allocation(); + space.fence(); } } diff --git a/packages/stokhos/src/sacado/kokkos/pce/KokkosExp_View_UQ_PCE_Contiguous.hpp b/packages/stokhos/src/sacado/kokkos/pce/KokkosExp_View_UQ_PCE_Contiguous.hpp index 1905d4dbb9f2..d32555a6c5da 100644 --- a/packages/stokhos/src/sacado/kokkos/pce/KokkosExp_View_UQ_PCE_Contiguous.hpp +++ b/packages/stokhos/src/sacado/kokkos/pce/KokkosExp_View_UQ_PCE_Contiguous.hpp @@ -1397,6 +1397,7 @@ class ViewMapping< Traits , /* View internal mapping */ // Only set the the pointer and initialize if the allocation is non-zero. // May be zero if one of the dimensions is zero. if ( alloc_size ) { + auto space = ((ViewCtorProp const &) prop).value; m_impl_handle.set( reinterpret_cast< pointer_type >( record->data() ), m_impl_offset.span(), m_sacado_size ); @@ -1404,7 +1405,7 @@ class ViewMapping< Traits , /* View internal mapping */ // Assume destruction is only required when construction is requested. // The ViewValueFunctor has both value construction and destruction operators. record->m_destroy = m_impl_handle.create_functor( - ( (ViewCtorProp const &) prop).value + space , ctor_prop::initialize , m_impl_offset.span() , m_sacado_size @@ -1412,6 +1413,7 @@ class ViewMapping< Traits , /* View internal mapping */ // Construct values record->m_destroy.construct_shared_allocation(); + space.fence(); } return record ; diff --git a/packages/stokhos/src/sacado/kokkos/vector/KokkosExp_View_MP_Vector_Contiguous.hpp b/packages/stokhos/src/sacado/kokkos/vector/KokkosExp_View_MP_Vector_Contiguous.hpp index 56767153e0cb..8b2f0e68e25d 100644 --- a/packages/stokhos/src/sacado/kokkos/vector/KokkosExp_View_MP_Vector_Contiguous.hpp +++ b/packages/stokhos/src/sacado/kokkos/vector/KokkosExp_View_MP_Vector_Contiguous.hpp @@ -1155,19 +1155,21 @@ class ViewMapping< Traits , /* View internal mapping */ // May be zero if one of the dimensions is zero. if ( alloc_size ) { + auto space = ((ViewCtorProp const &) prop).value; m_impl_handle.set( reinterpret_cast< pointer_type >( record->data() ), m_impl_offset.span(), m_sacado_size.value ); // Assume destruction is only required when construction is requested. // The ViewValueFunctor has both value construction and destruction operators. record->m_destroy = m_impl_handle.create_functor( - ( (ViewCtorProp const &) prop).value + space , ctor_prop::initialize , m_impl_offset.span() , m_sacado_size.value ); // Construct values record->m_destroy.construct_shared_allocation(); + space.fence(); } return record ;