From ef5553b464c84bc7aa5636e86aeffb7b96be6aff Mon Sep 17 00:00:00 2001 From: "Nathan V. Roberts" Date: Mon, 13 Apr 2020 14:42:40 -0600 Subject: [PATCH] Intrepid2: add support for hierarchical H(grad) bases on triangle and tetrahedron. --- .../Discretization/Basis/Intrepid2_Basis.hpp | 8 +- .../Basis/Intrepid2_DerivedBasisFamily.hpp | 70 +- .../Intrepid2_DerivedBasis_HGRAD_HEX.hpp | 4 +- .../Intrepid2_DerivedBasis_HGRAD_QUAD.hpp | 4 +- .../Basis/Intrepid2_HCURL_HEX_In_FEM.hpp | 4 +- .../Basis/Intrepid2_HDIV_HEX_In_FEM.hpp | 4 +- .../Basis/Intrepid2_HDIV_TET_I1_FEMDef.hpp | 2 +- .../Basis/Intrepid2_HGRAD_WEDGE_C2_FEMDef.hpp | 2 +- .../Basis/Intrepid2_HVOL_C0_FEM.hpp | 2 +- .../Intrepid2_HierarchicalBasisFamily.hpp | 40 +- ...pid2_IntegratedLegendreBasis_HGRAD_TET.hpp | 856 ++++++++++++++++++ ...pid2_IntegratedLegendreBasis_HGRAD_TRI.hpp | 568 ++++++++++++ .../Basis/Intrepid2_NodalBasisFamily.hpp | 27 +- .../Basis/Intrepid2_TensorBasis.hpp | 2 +- .../src/Shared/Intrepid2_Polynomials.hpp | 63 +- .../src/Shared/Intrepid2_TestUtils.hpp | 133 ++- .../Basis/BasisEquivalenceTests.cpp | 230 ++++- .../AnalyticPolynomialsMatchTests.cpp | 441 ++++++++- .../SubBasisInclusionTests.cpp | 108 ++- .../IntegratedJacobiTests.cpp | 81 ++ .../LegendreJacobiPolynomials/JacobiTests.cpp | 49 + 21 files changed, 2579 insertions(+), 119 deletions(-) create mode 100644 packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_TET.hpp create mode 100644 packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_TRI.hpp diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_Basis.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_Basis.hpp index 8976949b5e5d..c1eee3c53913 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_Basis.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_Basis.hpp @@ -359,7 +359,7 @@ namespace Intrepid2 { const PointViewType /* inputPoints */, const EOperator /* operatorType */ = OPERATOR_VALUE ) const { INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error, - ">>> ERROR (Basis::getValues): this method (FEM) is not supported or should be over-riden accordingly by derived classes."); + ">>> ERROR (Basis::getValues): this method (FEM) is not supported or should be overridden accordingly by derived classes."); } /** \brief Evaluation of an FVD basis evaluation on a physical cell. @@ -388,7 +388,7 @@ namespace Intrepid2 { const PointViewType /* cellVertices */, const EOperator /* operatorType */ = OPERATOR_VALUE ) const { INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error, - ">>> ERROR (Basis::getValues): this method (FVM) is not supported or should be over-riden accordingly by derived classes."); + ">>> ERROR (Basis::getValues): this method (FVM) is not supported or should be overridden accordingly by derived classes."); } @@ -399,7 +399,7 @@ namespace Intrepid2 { void getDofCoords( ScalarViewType /* dofCoords */ ) const { INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error, - ">>> ERROR (Basis::getDofCoords): this method is not supported or should be over-riden accordingly by derived classes."); + ">>> ERROR (Basis::getDofCoords): this method is not supported or should be overridden accordingly by derived classes."); } /** \brief Coefficients for computing degrees of freedom for Lagrangian basis @@ -414,7 +414,7 @@ namespace Intrepid2 { void getDofCoeffs( ScalarViewType /* dofCoeffs */ ) const { INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error, - ">>> ERROR (Basis::getDofCoeffs): this method is not supported or should be over-riden accordingly by derived classes."); + ">>> ERROR (Basis::getDofCoeffs): this method is not supported or should be overridden accordingly by derived classes."); } /** \brief For hierarchical bases, returns the field ordinals that have at most the specified degree in each dimension. diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasisFamily.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasisFamily.hpp index 886943cb7e12..3bb2c0d7018e 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasisFamily.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasisFamily.hpp @@ -63,12 +63,22 @@ namespace Intrepid2 { + //! \brief EmptyBasisFamily allows us to set a default void family for a given topology + class EmptyBasisFamily + { + public: + using HGRAD = void; + using HCURL = void; + using HDIV = void; + using HVOL = void; + }; + /** \class Intrepid2::DerivedBasisFamily \brief A family of basis functions, constructed from H(vol) and H(grad) bases on the line. At present, only hypercube topologies (line, quadrilateral, hexahedron) are supported, but other topologies will be supported in the future. */ - template + template class DerivedBasisFamily { public: @@ -89,11 +99,17 @@ namespace Intrepid2 using HDIV_QUAD = Basis_Derived_HDIV_QUAD ; using HVOL_QUAD = Basis_Derived_HVOL_QUAD ; - // hexahedral bases + // hexahedron bases using HGRAD_HEX = Basis_Derived_HGRAD_HEX; using HCURL_HEX = Basis_Derived_HCURL_HEX; using HDIV_HEX = Basis_Derived_HDIV_HEX ; using HVOL_HEX = Basis_Derived_HVOL_HEX ; + + // triangle bases + using HGRAD_TRI = typename TriangleBasisFamily::HGRAD; + + // tetrahedron bases + using HGRAD_TET = typename TetrahedronBasisFamily::HGRAD; }; /** \brief Factory method for line bases in the given family. @@ -152,12 +168,12 @@ namespace Intrepid2 } } - /** \brief Factory method for isotropic hexahedral bases in the given family. + /** \brief Factory method for isotropic bases on the hexahedron in the given family. \param [in] fs - the function space for the basis. \param [in] polyOrder - the polynomial order of the basis. */ template - static typename BasisFamily::BasisPtr getHexahedralBasis(Intrepid2::EFunctionSpace fs, int polyOrder) + static typename BasisFamily::BasisPtr getHexahedronBasis(Intrepid2::EFunctionSpace fs, int polyOrder) { using Teuchos::rcp; switch (fs) @@ -171,14 +187,14 @@ namespace Intrepid2 } } - /** \brief Factory method for potentially anisotropic hexahedral bases in the given family. + /** \brief Factory method for potentially anisotropic hexahedron bases in the given family. \param [in] fs - the function space for the basis. \param [in] polyOrder_x - the polynomial order of the basis in the x dimension. \param [in] polyOrder_y - the polynomial order of the basis in the y dimension. \param [in] polyOrder_z - the polynomial order of the basis in the z dimension. */ template - static typename BasisFamily::BasisPtr getHexahedralBasis(Intrepid2::EFunctionSpace fs, int polyOrder_x, int polyOrder_y, int polyOrder_z) + static typename BasisFamily::BasisPtr getHexahedronBasis(Intrepid2::EFunctionSpace fs, int polyOrder_x, int polyOrder_y, int polyOrder_z) { using Teuchos::rcp; switch (fs) @@ -192,6 +208,44 @@ namespace Intrepid2 } } + /** \brief Factory method for isotropic tetrahedron bases in the given family. + \param [in] fs - the function space for the basis. + \param [in] polyOrder - the polynomial order of the basis. + */ + template + static typename BasisFamily::BasisPtr getTetrahedronBasis(Intrepid2::EFunctionSpace fs, int polyOrder) + { + using Teuchos::rcp; + switch (fs) + { +// case FUNCTION_SPACE_HVOL: return rcp(new typename BasisFamily::HVOL_TET (polyOrder)); +// case FUNCTION_SPACE_HCURL: return rcp(new typename BasisFamily::HCURL_TET(polyOrder)); +// case FUNCTION_SPACE_HDIV: return rcp(new typename BasisFamily::HDIV_TET (polyOrder)); + case FUNCTION_SPACE_HGRAD: return rcp(new typename BasisFamily::HGRAD_TET(polyOrder)); + default: + INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported function space"); + } + } + + /** \brief Factory method for isotropic triangle bases in the given family. + \param [in] fs - the function space for the basis. + \param [in] polyOrder - the polynomial order of the basis. + */ + template + static typename BasisFamily::BasisPtr getTriangleBasis(Intrepid2::EFunctionSpace fs, int polyOrder) + { + using Teuchos::rcp; + switch (fs) + { +// case FUNCTION_SPACE_HVOL: return rcp(new typename BasisFamily::HVOL_TRI (polyOrder)); +// case FUNCTION_SPACE_HCURL: return rcp(new typename BasisFamily::HCURL_TRI(polyOrder)); +// case FUNCTION_SPACE_HDIV: return rcp(new typename BasisFamily::HDIV_TRI (polyOrder)); + case FUNCTION_SPACE_HGRAD: return rcp(new typename BasisFamily::HGRAD_TRI(polyOrder)); + default: + INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported function space"); + } + } + /** \brief Factory method for isotropic bases in the given family on the specified cell topology. \param [in] cellTopo - the cell topology on which the basis is defined. \param [in] fs - the function space for the basis. @@ -208,7 +262,9 @@ namespace Intrepid2 { case shards::Line<>::key: return getLineBasis(fs,polyOrder); case shards::Quadrilateral<>::key: return getQuadrilateralBasis(fs,polyOrder); - case shards::Hexahedron<>::key: return getHexahedralBasis(fs,polyOrder); + case shards::Triangle<>::key: return getTriangleBasis(fs,polyOrder); + case shards::Hexahedron<>::key: return getHexahedronBasis(fs,polyOrder); + case shards::Tetrahedron<>::key: return getTetrahedronBasis(fs,polyOrder); default: INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported cell topology"); } diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_HEX.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_HEX.hpp index fb6dfe7785ba..b8027656fa83 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_HEX.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_HEX.hpp @@ -106,13 +106,13 @@ namespace Intrepid2 */ virtual const char* - getName() const { + getName() const override { return "Intrepid2_DerivedBasis_HGRAD_HEX"; } /** \brief True if orientation is required */ - virtual bool requireOrientation() const { + virtual bool requireOrientation() const override { return (this->getDegree() > 2); } diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_QUAD.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_QUAD.hpp index c0a91bb74942..c79f66bcb88e 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_QUAD.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_QUAD.hpp @@ -94,13 +94,13 @@ namespace Intrepid2 */ virtual const char* - getName() const { + getName() const override { return "Intrepid2_DerivedBasis_HGRAD_QUAD"; } /** \brief True if orientation is required */ - virtual bool requireOrientation() const { + virtual bool requireOrientation() const override { return (this->getDegree() > 2); } diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HCURL_HEX_In_FEM.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HCURL_HEX_In_FEM.hpp index e80c48c0cdb2..da6d60806884 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HCURL_HEX_In_FEM.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HCURL_HEX_In_FEM.hpp @@ -169,9 +169,9 @@ namespace Intrepid2 { } /** \class Intrepid2::Basis_HCURL_HEX_In_FEM - \brief Implementation of the default H(curl)-compatible FEM basis on Hexahedral cell + \brief Implementation of the default H(curl)-compatible FEM basis on Hexahedron cell - Implements Nedelec basis of degree n on the reference Hexahedral cell. The basis has + Implements Nedelec basis of degree n on the reference Hexahedron cell. The basis has cardinality 3n (n+1)^2 and spans a INCOMPLETE polynomial space. */ diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HDIV_HEX_In_FEM.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HDIV_HEX_In_FEM.hpp index 33fc5f60703b..cb385cc8e730 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HDIV_HEX_In_FEM.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HDIV_HEX_In_FEM.hpp @@ -158,9 +158,9 @@ namespace Intrepid2 { } /** \class Intrepid2::Basis_HDIV_HEX_In_FEM - \brief Implementation of the default H(div)-compatible FEM basis on Hexahedral cell + \brief Implementation of the default H(div)-compatible FEM basis on Hexahedron cell - Implements Raviart-Thomas basis of degree n on the reference Hexahedral cell. The basis has + Implements Raviart-Thomas basis of degree n on the reference Hexahedron cell. The basis has cardinality 3(n+1)n^2 and spans a INCOMPLETE polynomial space. */ diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HDIV_TET_I1_FEMDef.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HDIV_TET_I1_FEMDef.hpp index daa2de822f1e..0bde33544797 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HDIV_TET_I1_FEMDef.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HDIV_TET_I1_FEMDef.hpp @@ -241,7 +241,7 @@ namespace Intrepid2 { dofCoords("dofCoordsHost", this->basisCardinality_,this->basisCellTopology_.getDimension()); dofCoords(0,0) = 1.0/3.0; dofCoords(0,1) = 0.0; dofCoords(0,2) = 1.0/3.0; - dofCoords(1,0) = 1.0/3.0; dofCoords(1,1) = 1.0/3.0; dofCoords(1,2) = 1.0/3.0;; + dofCoords(1,0) = 1.0/3.0; dofCoords(1,1) = 1.0/3.0; dofCoords(1,2) = 1.0/3.0; dofCoords(2,0) = 0.0; dofCoords(2,1) = 1.0/3.0; dofCoords(2,2) = 1.0/3.0; dofCoords(3,0) = 1.0/3.0; dofCoords(3,1) = 1.0/3.0; dofCoords(3,2) = 0.0; 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 0a6f7cd9ea16..cd2af7448192 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 @@ -204,7 +204,7 @@ namespace Intrepid2 { output.access(4, 0) = 2.*z*(1. + z); output.access(4, 1) = 0.; - output.access(4, 2) = ((-1. + 4.*x)*(1. + 2.*z))/2.;; + output.access(4, 2) = ((-1. + 4.*x)*(1. + 2.*z))/2.; output.access(4, 3) = 0.; output.access(4, 4) = 0.; output.access(4, 5) = x*(-1. + 2.*x); diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HVOL_C0_FEM.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HVOL_C0_FEM.hpp index f505349d78af..783969e5dbfe 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_HVOL_C0_FEM.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_HVOL_C0_FEM.hpp @@ -124,7 +124,7 @@ namespace Intrepid2 { } /** \class Intrepid2::Basis_HVOL_C0_FEM - \brief Implementation of the default HVOL-compatible FEM contstant basis on triangle, quadrilateral, hexahedral and tetrahedral cells. + \brief Implementation of the default HVOL-compatible FEM contstant basis on triangle, quadrilateral, hexahedron and tetrahedron cells. */ template + class HierarchicalTriangleBasisFamily + { + public: + // we will fill these in as we implement them + using HGRAD = IntegratedLegendreBasis_HGRAD_TRI; + using HCURL = void; + using HDIV = void; + using HVOL = void; + }; + + template + class HierarchicalTetrahedronBasisFamily + { + public: + // we will fill these in as we implement them + using HGRAD = IntegratedLegendreBasis_HGRAD_TET; + using HCURL = void; + using HDIV = void; + using HVOL = void; + }; + /** \class Intrepid2::HierarchicalBasisFamily \brief A family of hierarchical basis functions, constructed in a way that follows work by Fuentes et al. @@ -82,7 +112,10 @@ namespace Intrepid2 { typename OutputScalar = double, typename PointScalar = double> using HierarchicalBasisFamily = DerivedBasisFamily< IntegratedLegendreBasis_HGRAD_LINE, - LegendreBasis_HVOL_LINE >; + LegendreBasis_HVOL_LINE, + HierarchicalTriangleBasisFamily, + HierarchicalTetrahedronBasisFamily + >; /** \class Intrepid2::HierarchicalBasisFamily \brief A family of hierarchical basis functions, constructed in a way that follows work by Fuentes et al., suitable for use in DG contexts. @@ -96,7 +129,10 @@ namespace Intrepid2 { typename OutputScalar = double, typename PointScalar = double> using DGHierarchicalBasisFamily = DerivedBasisFamily< IntegratedLegendreBasis_HGRAD_LINE, - LegendreBasis_HVOL_LINE >; + LegendreBasis_HVOL_LINE, + HierarchicalTriangleBasisFamily, + HierarchicalTetrahedronBasisFamily + >; } diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_TET.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_TET.hpp new file mode 100644 index 000000000000..255c587ada40 --- /dev/null +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_TET.hpp @@ -0,0 +1,856 @@ +// @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), +// Mauro Perego (mperego@sandia.gov), or +// Nate Roberts (nvrober@sandia.gov) +// +// ************************************************************************ +// @HEADER + +/** \file Intrepid2_IntegratedLegendreBasis_HGRAD_TET.hpp + \brief H(grad) basis on the tetrahedon based on integrated Legendre polynomials. + \author Created by N.V. Roberts. + */ + +#ifndef Intrepid2_IntegratedLegendreBasis_HGRAD_TET_h +#define Intrepid2_IntegratedLegendreBasis_HGRAD_TET_h + +#include +#include + +#include + +#include "Intrepid2_Polynomials.hpp" +#include "Intrepid2_Utils.hpp" + +namespace Intrepid2 +{ + /** \class Intrepid2::Hierarchical_HGRAD_TET_Functor + \brief Functor for computing values for the IntegratedLegendreBasis_HGRAD_TET class. + + This functor is not intended for use outside of IntegratedLegendreBasis_HGRAD_TET. + */ + template + struct Hierarchical_HGRAD_TET_Functor + { + using ScratchSpace = Kokkos::DefaultExecutionSpace::scratch_memory_space; + using OutputScratchView = Kokkos::View>; + using OutputScratchView2D = Kokkos::View>; + using PointScratchView = Kokkos::View>; + + using TeamPolicy = Kokkos::TeamPolicy<>; + using TeamMember = TeamPolicy::member_type; + + EOperator opType_; + + OutputFieldType output_; // F,P + InputPointsType inputPoints_; // P,D + + int polyOrder_; + bool defineVertexFunctions_; + int numFields_, numPoints_; + + size_t fad_size_output_; + + static const int numVertices = 4; + static const int numEdges = 6; + // the following ordering of the edges matches that used by ESEAS + const int edge_start_[numEdges] = {0,1,0,0,1,2}; // edge i is from edge_start_[i] to edge_end_[i] + const int edge_end_[numEdges] = {1,2,2,3,3,3}; // edge i is from edge_start_[i] to edge_end_[i] + + static const int numFaces = 4; + const int face_vertex_0[numFaces] = {0,0,1,0}; // faces are abc where 0 ≤ a < b < c ≤ 3 + const int face_vertex_1[numFaces] = {1,1,2,2}; + const int face_vertex_2[numFaces] = {2,3,3,3}; + + // this allows us to look up the edge ordinal of the first edge of a face + // this is useful because face functions are defined using edge basis functions of the first edge of the face + const int face_ordinal_of_first_edge[numFaces] = {0,0,1,2}; + + Hierarchical_HGRAD_TET_Functor(EOperator opType, OutputFieldType output, InputPointsType inputPoints, + int polyOrder, bool defineVertexFunctions) + : opType_(opType), output_(output), inputPoints_(inputPoints), + polyOrder_(polyOrder), defineVertexFunctions_(defineVertexFunctions), + fad_size_output_(getScalarDimensionForView(output)) + { + numFields_ = output.extent_int(0); + numPoints_ = output.extent_int(1); + INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument, "point counts need to match!"); + INTREPID2_TEST_FOR_EXCEPTION(numFields_ != (polyOrder_+1)*(polyOrder_+2)*(polyOrder_+3)/6, std::invalid_argument, "output field size does not match basis cardinality"); + } + + KOKKOS_INLINE_FUNCTION + void operator()( const TeamMember & teamMember ) const + { + const int numFaceBasisFunctionsPerFace = (polyOrder_-2) * (polyOrder_-1) / 2; + const int numInteriorBasisFunctions = (polyOrder_-3) * (polyOrder_-2) * (polyOrder_-1) / 6; + + auto pointOrdinal = teamMember.league_rank(); + OutputScratchView legendre_values1_at_point, legendre_values2_at_point; + OutputScratchView2D jacobi_values1_at_point, jacobi_values2_at_point, jacobi_values3_at_point; + const int numAlphaValues = (polyOrder_-1 > 1) ? (polyOrder_-1) : 1; // make numAlphaValues at least 1 so we can avoid zero-extent allocations… + if (fad_size_output_ > 0) { + legendre_values1_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_); + legendre_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_); + jacobi_values1_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1, fad_size_output_); + jacobi_values2_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1, fad_size_output_); + jacobi_values3_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1, fad_size_output_); + } + else { + legendre_values1_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1); + legendre_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1); + jacobi_values1_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1); + jacobi_values2_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1); + jacobi_values3_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1); + } + + const auto & x = inputPoints_(pointOrdinal,0); + const auto & y = inputPoints_(pointOrdinal,1); + const auto & z = inputPoints_(pointOrdinal,2); + + // write as barycentric coordinates: + const PointScalar lambda[numVertices] = {1. - x - y - z, x, y, z}; + const PointScalar lambda_dx[numVertices] = {-1., 1., 0., 0.}; + const PointScalar lambda_dy[numVertices] = {-1., 0., 1., 0.}; + const PointScalar lambda_dz[numVertices] = {-1., 0., 0., 1.}; + + const int num1DEdgeFunctions = polyOrder_ - 1; + + switch (opType_) + { + case OPERATOR_VALUE: + { + // vertex functions come first, according to vertex ordering: (0,0,0), (1,0,0), (0,1,0), (0,0,1) + for (int vertexOrdinal=0; vertexOrdinal>> ERROR: (Intrepid2::Basis_HGRAD_TET_Cn_FEM_ORTH::OrthPolynomialTri) Computing of second and higher-order derivatives is not currently supported"); + default: + // unsupported operator type + device_assert(false); + } + } + + // Provide the shared memory capacity. + // This function takes the team_size as an argument, + // which allows team_size-dependent allocations. + size_t team_shmem_size (int team_size) const + { + // we will use shared memory to create a fast buffer for basis computations + // for the (integrated) Legendre computations, we just need p+1 values stored + // for the (integrated) Jacobi computations, though, we want (p+1)*(# alpha values) + // alpha is either 2i or 2(i+j), where i=2,…,p or i+j=3,…,p. So there are at most (p-1) alpha values needed. + // We can have up to 3 of the (integrated) Jacobi values needed at once. + const int numAlphaValues = std::max(polyOrder_-1, 1); // make it at least 1 so we can avoid zero-extent ranks… + size_t shmem_size = 0; + if (fad_size_output_ > 0) + { + // Legendre: + shmem_size += 2 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_); + // Jacobi: + shmem_size += 3 * OutputScratchView2D::shmem_size(numAlphaValues, polyOrder_ + 1, fad_size_output_); + } + else + { + // Legendre: + shmem_size += 2 * OutputScratchView::shmem_size(polyOrder_ + 1); + // Jacobi: + shmem_size += 3 * OutputScratchView2D::shmem_size(numAlphaValues, polyOrder_ + 1); + } + + return shmem_size; + } + }; + + /** \class Intrepid2::IntegratedLegendreBasis_HGRAD_TET + \brief Basis defining integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line. + + This is used in the construction of hierarchical bases on higher-dimensional topologies. For + mathematical details of the construction, see: + + Federico Fuentes, Brendan Keith, Leszek Demkowicz, Sriram Nagaraj. + "Orientation embedded high order shape functions for the exact sequence elements of all shapes." + Computers & Mathematics with Applications, Volume 70, Issue 4, 2015, Pages 353-458, ISSN 0898-1221. + https://doi.org/10.1016/j.camwa.2015.04.027. + + Note that the template argument defineVertexFunctions controls whether this basis is defined in a + strictly hierarchical way. If defineVertexFunctions is false, then the first basis function is the + constant 1.0, and this basis will be suitable for discontinuous discretizations. If defineVertexFunctions + is true, then the first basis function will instead be 1.0-x-y, and the basis will be suitable for + continuous discretizations. + */ + template // if defineVertexFunctions is true, first four basis functions are 1-x-y-z, x, y, and z. Otherwise, they are 1, x, y, and z. + class IntegratedLegendreBasis_HGRAD_TET + : public Basis + { + public: + using OrdinalTypeArray1DHost = typename Basis::OrdinalTypeArray1DHost; + using OrdinalTypeArray2DHost = typename Basis::OrdinalTypeArray2DHost; + + using OutputViewType = typename Basis::OutputViewType; + using PointViewType = typename Basis::PointViewType; + using ScalarViewType = typename Basis::ScalarViewType; + protected: + int polyOrder_; // the maximum order of the polynomial + bool defineVertexFunctions_; // if true, first four basis functions are 1-x-y-z, x, y, and z. Otherwise, they are 1, x, y, and z. + public: + /** \brief Constructor. + \param [in] polyOrder - the polynomial order of the basis. + + The basis will have polyOrder + 1 members. + + If defineVertexFunctions is false, then all basis functions are identified with the interior of the line element, and the first four basis functions are 1, x, y, and z. + + If defineVertexFunctions is true, then the first two basis functions are 1-x-y-z, x, y, and z, and these are identified with the left and right vertices of the cell. + + */ + IntegratedLegendreBasis_HGRAD_TET(int polyOrder) + : + polyOrder_(polyOrder) + { + this->basisCardinality_ = ((polyOrder+1) * (polyOrder+2) * (polyOrder+3)) / 6; + this->basisDegree_ = polyOrder; + this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData >() ); + this->basisType_ = BASIS_FEM_HIERARCHICAL; + this->basisCoordinates_ = COORDINATES_CARTESIAN; + this->functionSpace_ = FUNCTION_SPACE_HGRAD; + + const int degreeLength = 1; + this->fieldOrdinalPolynomialDegree_ = OrdinalTypeArray2DHost("Integrated Legendre H(grad) tetrahedron polynomial degree lookup", this->basisCardinality_, degreeLength); + + int fieldOrdinalOffset = 0; + // **** vertex functions **** // + const int numVertices = this->basisCellTopology_.getVertexCount(); + const int numFunctionsPerVertex = 1; + const int numVertexFunctions = numVertices * numFunctionsPerVertex; + for (int i=0; ifieldOrdinalPolynomialDegree_(i,0) = 1; + } + if (!defineVertexFunctions) + { + this->fieldOrdinalPolynomialDegree_(0,0) = 0; + } + fieldOrdinalOffset += numVertexFunctions; + + // **** edge functions **** // + const int numFunctionsPerEdge = polyOrder - 1; // bubble functions: all but the vertices + const int numEdges = this->basisCellTopology_.getEdgeCount(); + for (int edgeOrdinal=0; edgeOrdinalfieldOrdinalPolynomialDegree_(i+fieldOrdinalOffset,0) = i+2; // vertex functions are 1st order; edge functions start at order 2 + } + fieldOrdinalOffset += numFunctionsPerEdge; + } + + // **** face functions **** // + const int numFunctionsPerFace = ((polyOrder-1)*(polyOrder-2))/2; + const int numFaces = 4; + for (int faceOrdinal=0; faceOrdinalfieldOrdinalPolynomialDegree_(fieldOrdinalOffset,0) = totalPolyOrder; + fieldOrdinalOffset++; + } + } + } + + // **** interior functions **** // + const int numFunctionsPerVolume = ((polyOrder-1)*(polyOrder-2)*(polyOrder-3))/6; + const int numVolumes = 1; // interior + for (int volumeOrdinal=0; volumeOrdinalfieldOrdinalPolynomialDegree_(fieldOrdinalOffset,0) = totalPolyOrder; + fieldOrdinalOffset++; + } + } + } + + INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != this->basisCardinality_, std::invalid_argument, "Internal error: basis enumeration is incorrect"); + + // initialize tags + { + // ESEAS numbers tetrahedron faces differently from Intrepid2 + // ESEAS: 012, 013, 123, 023 + // Intrepid2: 013, 123, 032, 021 + const int intrepid2FaceOrdinals[4] {3,0,1,2}; // index is the ESEAS face ordinal; value is the intrepid2 ordinal + + const auto & cardinality = this->basisCardinality_; + + // Basis-dependent initializations + const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag + const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim + const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal + const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell + + OrdinalTypeArray1DHost tagView("tag view", cardinality*tagSize); + const int vertexDim = 0, edgeDim = 1, faceDim = 2, volumeDim = 3; + + if (defineVertexFunctions) { + { + int tagNumber = 0; + for (int vertexOrdinal=0; vertexOrdinalbasisCardinality_, std::invalid_argument, "Internal error: basis tag enumeration is incorrect"); + } + } else { + for (ordinal_type i=0;isetOrdinalTagData(this->tagToOrdinal_, + this->ordinalToTag_, + tagView, + this->basisCardinality_, + tagSize, + posScDim, + posScOrd, + posDfOrd); + } + } + + /** \brief Returns basis name + + \return the name of the basis + */ + const char* getName() const override { + return "Intrepid2_IntegratedLegendreBasis_HGRAD_TET"; + } + + // since the getValues() below only overrides the FEM variant, we specify that + // we use the base class's getValues(), which implements the FVD variant by throwing an exception. + // (It's an error to use the FVD variant on this basis.) + using Basis::getValues; + + /** \brief Evaluation of a FEM basis on a reference cell. + + Returns values of operatorType acting on FEM basis functions for a set of + points in the reference cell for which the basis is defined. + + \param outputValues [out] - variable rank array with the basis values + \param inputPoints [in] - rank-2 array (P,D) with the evaluation points + \param operatorType [in] - the operator acting on the basis functions + + \remark For rank and dimension specifications of the output array see Section + \ref basis_md_array_sec. Dimensions of ArrayScalar arguments are checked + at runtime if HAVE_INTREPID2_DEBUG is defined. + + \remark A FEM basis spans a COMPLETE or INCOMPLETE polynomial space on the reference cell + which is a smooth function space. Thus, all operator types that are meaningful for the + approximated function space are admissible. When the order of the operator exceeds the + degree of the basis, the output array is filled with the appropriate number of zeros. + */ + virtual void getValues( OutputViewType outputValues, const PointViewType inputPoints, + const EOperator operatorType = OPERATOR_VALUE ) const override + { + auto numPoints = inputPoints.extent_int(0); + + using FunctorType = Hierarchical_HGRAD_TET_Functor; + + FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_, defineVertexFunctions); + + const int outputVectorSize = getVectorSizeForHierarchicalParallelism(); + const int pointVectorSize = getVectorSizeForHierarchicalParallelism(); + const int vectorSize = std::max(outputVectorSize,pointVectorSize); + const int teamSize = 1; // because of the way the basis functions are computed, we don't have a second level of parallelism... + + auto policy = Kokkos::TeamPolicy(numPoints,teamSize,vectorSize); + Kokkos::parallel_for( policy , functor, "Hierarchical_HGRAD_TET_Functor"); + } + }; +} // end namespace Intrepid2 + +#endif /* Intrepid2_IntegratedLegendreBasis_HGRAD_TET_h */ diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_TRI.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_TRI.hpp new file mode 100644 index 000000000000..a6923ce624d7 --- /dev/null +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_TRI.hpp @@ -0,0 +1,568 @@ +// @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), +// Mauro Perego (mperego@sandia.gov), or +// Nate Roberts (nvrober@sandia.gov) +// +// ************************************************************************ +// @HEADER + +/** \file Intrepid2_IntegratedLegendreBasis_HGRAD_TRI.hpp + \brief H(grad) basis on the triangle based on integrated Legendre polynomials. + \author Created by N.V. Roberts. + */ + +#ifndef Intrepid2_IntegratedLegendreBasis_HGRAD_TRI_h +#define Intrepid2_IntegratedLegendreBasis_HGRAD_TRI_h + +#include +#include + +#include + +#include "Intrepid2_Polynomials.hpp" +#include "Intrepid2_Utils.hpp" + +namespace Intrepid2 +{ + /** \class Intrepid2::Hierarchical_HGRAD_TRI_Functor + \brief Functor for computing values for the IntegratedLegendreBasis_HGRAD_TRI class. + + This functor is not intended for use outside of IntegratedLegendreBasis_HGRAD_TRI. + */ + template + struct Hierarchical_HGRAD_TRI_Functor + { + using ScratchSpace = Kokkos::DefaultExecutionSpace::scratch_memory_space; + using OutputScratchView = Kokkos::View>; + using PointScratchView = Kokkos::View>; + + using TeamPolicy = Kokkos::TeamPolicy<>; + using TeamMember = TeamPolicy::member_type; + + EOperator opType_; + + OutputFieldType output_; // F,P + InputPointsType inputPoints_; // P,D + + int polyOrder_; + bool defineVertexFunctions_; + int numFields_, numPoints_; + + size_t fad_size_output_; + + static const int numVertices = 3; + static const int numEdges = 3; + const int edge_start_[numEdges] = {0,1,0}; // edge i is from edge_start_[i] to edge_end_[i] + const int edge_end_[numEdges] = {1,2,2}; // edge i is from edge_start_[i] to edge_end_[i] + + Hierarchical_HGRAD_TRI_Functor(EOperator opType, OutputFieldType output, InputPointsType inputPoints, + int polyOrder, bool defineVertexFunctions) + : opType_(opType), output_(output), inputPoints_(inputPoints), + polyOrder_(polyOrder), defineVertexFunctions_(defineVertexFunctions), + fad_size_output_(getScalarDimensionForView(output)) + { + numFields_ = output.extent_int(0); + numPoints_ = output.extent_int(1); + INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument, "point counts need to match!"); + INTREPID2_TEST_FOR_EXCEPTION(numFields_ != (polyOrder_+1)*(polyOrder_+2)/2, std::invalid_argument, "output field size does not match basis cardinality"); + } + + KOKKOS_INLINE_FUNCTION + void operator()( const TeamMember & teamMember ) const + { + auto pointOrdinal = teamMember.league_rank(); + OutputScratchView edge_field_values_at_point, jacobi_values_at_point, other_values_at_point, other_values2_at_point; + if (fad_size_output_ > 0) { + edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_); + jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_); + other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_); + other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_); + } + else { + edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1); + jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1); + other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1); + other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1); + } + + const auto & x = inputPoints_(pointOrdinal,0); + const auto & y = inputPoints_(pointOrdinal,1); + + // write as barycentric coordinates: + const PointScalar lambda[3] = {1. - x - y, x, y}; + const PointScalar lambda_dx[3] = {-1., 1., 0.}; + const PointScalar lambda_dy[3] = {-1., 0., 1.}; + + const int num1DEdgeFunctions = polyOrder_ - 1; + + switch (opType_) + { + case OPERATOR_VALUE: + { + // vertex functions come first, according to vertex ordering: (0,0), (1,0), (0,1) + for (int vertexOrdinal=0; vertexOrdinal>> ERROR: (Intrepid2::Basis_HGRAD_TRI_Cn_FEM_ORTH::OrthPolynomialTri) Computing of second and higher-order derivatives is not currently supported"); + default: + // unsupported operator type + device_assert(false); + } + } + + // Provide the shared memory capacity. + // This function takes the team_size as an argument, + // which allows team_size-dependent allocations. + size_t team_shmem_size (int team_size) const + { + // we will use shared memory to create a fast buffer for basis computations + size_t shmem_size = 0; + if (fad_size_output_ > 0) + shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_); + else + shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1); + + return shmem_size; + } + }; + + /** \class Intrepid2::IntegratedLegendreBasis_HGRAD_TRI + \brief Basis defining integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line. + + This is used in the construction of hierarchical bases on higher-dimensional topologies. For + mathematical details of the construction, see: + + Federico Fuentes, Brendan Keith, Leszek Demkowicz, Sriram Nagaraj. + "Orientation embedded high order shape functions for the exact sequence elements of all shapes." + Computers & Mathematics with Applications, Volume 70, Issue 4, 2015, Pages 353-458, ISSN 0898-1221. + https://doi.org/10.1016/j.camwa.2015.04.027. + + Note that the template argument defineVertexFunctions controls whether this basis is defined in a + strictly hierarchical way. If defineVertexFunctions is false, then the first basis function is the + constant 1.0, and this basis will be suitable for discontinuous discretizations. If defineVertexFunctions + is true, then the first basis function will instead be 1.0-x-y, and the basis will be suitable for + continuous discretizations. + */ + template // if defineVertexFunctions is true, first three basis functions are 1-x-y, x, and y. Otherwise, they are 1, x, and y. + class IntegratedLegendreBasis_HGRAD_TRI + : public Basis + { + public: + using OrdinalTypeArray1DHost = typename Basis::OrdinalTypeArray1DHost; + using OrdinalTypeArray2DHost = typename Basis::OrdinalTypeArray2DHost; + + using OutputViewType = typename Basis::OutputViewType; + using PointViewType = typename Basis::PointViewType; + using ScalarViewType = typename Basis::ScalarViewType; + protected: + int polyOrder_; // the maximum order of the polynomial + bool defineVertexFunctions_; // if true, first and second basis functions are x and 1-x. Otherwise, they are 1 and x. + public: + /** \brief Constructor. + \param [in] polyOrder - the polynomial order of the basis. + + The basis will have polyOrder + 1 members. + + If defineVertexFunctions is false, then all basis functions are identified with the interior of the line element, and the first two basis functions are 1 and x. + + If defineVertexFunctions is true, then the first two basis functions are 1-x and x, and these are identified with the left and right vertices of the cell. + + */ + IntegratedLegendreBasis_HGRAD_TRI(int polyOrder) + : + polyOrder_(polyOrder) + { + this->basisCardinality_ = ((polyOrder+2) * (polyOrder+1)) / 2; + this->basisDegree_ = polyOrder; + this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData >() ); + this->basisType_ = BASIS_FEM_HIERARCHICAL; + this->basisCoordinates_ = COORDINATES_CARTESIAN; + this->functionSpace_ = FUNCTION_SPACE_HGRAD; + + const int degreeLength = 1; + this->fieldOrdinalPolynomialDegree_ = OrdinalTypeArray2DHost("Integrated Legendre H(grad) triangle polynomial degree lookup", this->basisCardinality_, degreeLength); + + int fieldOrdinalOffset = 0; + // **** vertex functions **** // + const int numVertices = this->basisCellTopology_.getVertexCount(); + const int numFunctionsPerVertex = 1; + const int numVertexFunctions = numVertices * numFunctionsPerVertex; + for (int i=0; ifieldOrdinalPolynomialDegree_(i,0) = 1; + } + if (!defineVertexFunctions) + { + this->fieldOrdinalPolynomialDegree_(0,0) = 0; + } + fieldOrdinalOffset += numVertexFunctions; + + // **** edge functions **** // + const int numFunctionsPerEdge = polyOrder - 1; // bubble functions: all but the vertices + const int numEdges = this->basisCellTopology_.getEdgeCount(); + for (int edgeOrdinal=0; edgeOrdinalfieldOrdinalPolynomialDegree_(i+fieldOrdinalOffset,0) = i+2; // vertex functions are 1st order; edge functions start at order 2 + } + fieldOrdinalOffset += numFunctionsPerEdge; + } + + // **** face functions **** // + const int max_ij_sum = polyOrder; + for (int i=2; ifieldOrdinalPolynomialDegree_(fieldOrdinalOffset,0) = i+j; + fieldOrdinalOffset++; + } + } + const int numFaces = 1; + const int numFunctionsPerFace = ((polyOrder-1)*(polyOrder-2))/2; + INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != this->basisCardinality_, std::invalid_argument, "Internal error: basis enumeration is incorrect"); + + // initialize tags + { + const auto & cardinality = this->basisCardinality_; + + // Basis-dependent initializations + const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag + const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim + const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal + const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell + + OrdinalTypeArray1DHost tagView("tag view", cardinality*tagSize); + const int vertexDim = 0, edgeDim = 1, faceDim = 2; + + if (defineVertexFunctions) { + { + int tagNumber = 0; + for (int vertexOrdinal=0; vertexOrdinalsetOrdinalTagData(this->tagToOrdinal_, + this->ordinalToTag_, + tagView, + this->basisCardinality_, + tagSize, + posScDim, + posScOrd, + posDfOrd); + } + } + + /** \brief Returns basis name + + \return the name of the basis + */ + const char* getName() const override { + return "Intrepid2_IntegratedLegendreBasis_HGRAD_TRI"; + } + + // since the getValues() below only overrides the FEM variant, we specify that + // we use the base class's getValues(), which implements the FVD variant by throwing an exception. + // (It's an error to use the FVD variant on this basis.) + using Basis::getValues; + + /** \brief Evaluation of a FEM basis on a reference cell. + + Returns values of operatorType acting on FEM basis functions for a set of + points in the reference cell for which the basis is defined. + + \param outputValues [out] - variable rank array with the basis values + \param inputPoints [in] - rank-2 array (P,D) with the evaluation points + \param operatorType [in] - the operator acting on the basis functions + + \remark For rank and dimension specifications of the output array see Section + \ref basis_md_array_sec. Dimensions of ArrayScalar arguments are checked + at runtime if HAVE_INTREPID2_DEBUG is defined. + + \remark A FEM basis spans a COMPLETE or INCOMPLETE polynomial space on the reference cell + which is a smooth function space. Thus, all operator types that are meaningful for the + approximated function space are admissible. When the order of the operator exceeds the + degree of the basis, the output array is filled with the appropriate number of zeros. + */ + virtual void getValues( OutputViewType outputValues, const PointViewType inputPoints, + const EOperator operatorType = OPERATOR_VALUE ) const override + { + auto numPoints = inputPoints.extent_int(0); + + using FunctorType = Hierarchical_HGRAD_TRI_Functor; + + FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_, defineVertexFunctions); + + const int outputVectorSize = getVectorSizeForHierarchicalParallelism(); + const int pointVectorSize = getVectorSizeForHierarchicalParallelism(); + const int vectorSize = std::max(outputVectorSize,pointVectorSize); + const int teamSize = 1; // because of the way the basis functions are computed, we don't have a second level of parallelism... + + auto policy = Kokkos::TeamPolicy(numPoints,teamSize,vectorSize); + Kokkos::parallel_for( policy , functor, "Hierarchical_HGRAD_TRI_Functor"); + } + }; +} // end namespace Intrepid2 + +#endif /* Intrepid2_IntegratedLegendreBasis_HGRAD_TRI_h */ diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_NodalBasisFamily.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_NodalBasisFamily.hpp index ee80aede0f89..d12ce335641a 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_NodalBasisFamily.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_NodalBasisFamily.hpp @@ -59,11 +59,21 @@ #include #include +#include +#include +#include +#include + #include #include #include #include +#include +#include +#include +#include + namespace Intrepid2 { // the following defines a family of nodal basis functions, derived from a the standard high-order Intrepid2 bases on the line // note that because the standard H(curl) basis uses a lower-order H(grad) basis in place of the H(vol) that is arguably more natural, @@ -112,11 +122,24 @@ namespace Intrepid2 { using HDIV_QUAD = Basis_HDIV_QUAD_In_FEM; using HVOL_QUAD = Basis_HVOL_QUAD_Cn_FEM; - // hexahedral bases + // triangle bases + using HGRAD_TRI = Basis_HGRAD_TRI_Cn_FEM; + using HCURL_TRI = Basis_HCURL_TRI_In_FEM; + using HDIV_TRI = Basis_HDIV_TRI_In_FEM; + using HVOL_TRI = Basis_HVOL_TRI_Cn_FEM; + + // hexahedron bases using HGRAD_HEX = Basis_HGRAD_HEX_Cn_FEM; using HCURL_HEX = Basis_HCURL_HEX_In_FEM; using HDIV_HEX = Basis_HDIV_HEX_In_FEM; using HVOL_HEX = Basis_HVOL_HEX_Cn_FEM; + + // tetrahedron bases + using HGRAD_TET = Basis_HGRAD_TET_Cn_FEM; + using HCURL_TET = Basis_HCURL_TET_In_FEM; + using HDIV_TET = Basis_HDIV_TET_In_FEM; + using HVOL_TET = Basis_HVOL_TET_Cn_FEM; + }; /** \class Intrepid2::DerivedBasisFamilyModified @@ -149,7 +172,7 @@ namespace Intrepid2 { using HDIV_QUAD = Basis_Derived_HDIV_QUAD ; using HVOL_QUAD = Basis_Derived_HVOL_QUAD ; - // hexahedral bases + // hexahedron bases using HGRAD_HEX = Basis_Derived_HGRAD_HEX; using HCURL_HEX = Basis_Derived_HCURL_HEX; using HDIV_HEX = Basis_Derived_HDIV_HEX ; diff --git a/packages/intrepid2/src/Discretization/Basis/Intrepid2_TensorBasis.hpp b/packages/intrepid2/src/Discretization/Basis/Intrepid2_TensorBasis.hpp index 60b832c480a3..f19c03ceacc8 100644 --- a/packages/intrepid2/src/Discretization/Basis/Intrepid2_TensorBasis.hpp +++ b/packages/intrepid2/src/Discretization/Basis/Intrepid2_TensorBasis.hpp @@ -752,7 +752,7 @@ namespace Intrepid2 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "one-operator, two-inputPoints getValues should be overridden by TensorBasis subclasses"); } - /** \brief Evaluation of a tensor FEM basis on a reference cell; subclasses should override this. + /** \brief Evaluation of a tensor FEM basis on a reference cell. Returns values of operatorType acting on FEM basis functions for a set of points in the reference cell for which the basis is defined. diff --git a/packages/intrepid2/src/Shared/Intrepid2_Polynomials.hpp b/packages/intrepid2/src/Shared/Intrepid2_Polynomials.hpp index 676f2560ada1..c751c3b9575a 100644 --- a/packages/intrepid2/src/Shared/Intrepid2_Polynomials.hpp +++ b/packages/intrepid2/src/Shared/Intrepid2_Polynomials.hpp @@ -181,11 +181,12 @@ namespace Intrepid2 Shifted, scaled Legendre polynomials are given by P_i(x;t) = P_i(x/t) * t^i = ~P_i(2x-t;t). */ - template - KOKKOS_INLINE_FUNCTION void shiftedScaledLegendreValues(OutputValueViewType outputValues, Intrepid2::ordinal_type n, ScalarType x, double t) + template + KOKKOS_INLINE_FUNCTION void shiftedScaledLegendreValues(OutputValueViewType outputValues, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t) { - ScalarType two_x_minus_t = 2. * x - t; - ScalarType t_squared = t * t; + using OutputScalar = typename OutputValueViewType::value_type; + OutputScalar two_x_minus_t = 2. * x - t; + OutputScalar t_squared = t * t; if (n >= 0) outputValues(0) = 1.0; if (n >= 1) outputValues(1) = two_x_minus_t; for (int i=2; i<=n; i++) @@ -199,7 +200,7 @@ namespace Intrepid2 \param [out] outputValues - the view into which to place the output values (must have at least n+1 entries) \param [in] n - the maximum polynomial order of integrated Legendre polynomials to compute \param [in] x - point at which to evaluate the polynomials - \param [in] t - scaling parameter + \param [in] t - scaling parameter; may be of type double or may match the type of x These are defined for x in [0,1]. See equation (2.18) in Fuentes et al. (We additionally define L_0 = 1.) Shifted, scaled Legendre polynomials are given by @@ -208,8 +209,8 @@ namespace Intrepid2 The formula in Fuentes et al. is defined in terms of P_i and P_{i-2}. We offer two versions of this computation, one which can reuse an existing P_i computation (in the form of a shiftedScaledLegendreValues input container), and one which reuses space in outputValues. */ - template - KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedLegendreValues(OutputValueViewType outputValues, Intrepid2::ordinal_type n, ScalarType x, double t) + template + KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedLegendreValues(OutputValueViewType outputValues, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t) { // reduced memory version: compute P_i in outputValues shiftedScaledLegendreValues(outputValues,n,x,t); @@ -249,9 +250,9 @@ namespace Intrepid2 The formula in Fuentes et al. is defined in terms of P_i and P_{i-2}. We offer two versions of this computation, one which can reuse an existing P_i computation (in the form of a shiftedScaledLegendreValues input container), and one which reuses space in outputValues. */ - template + template KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedLegendreValues(OutputValueViewType outputValues, const OutputValueViewType shiftedScaledLegendreValues, - Intrepid2::ordinal_type n, ScalarType x, double t) + Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t) { // reduced flops version: rely on previously computed P_i if (n >= 0) outputValues(0) = 1.0; @@ -303,8 +304,8 @@ namespace Intrepid2 \param [in] t - scaling parameter These are defined for x in [0,1]. The x derivative of integrated Legendre is just Legendre; the only distinction is in the index -- outputValues indices are shifted by 1 relative to shiftedScaledLegendreValues, above. */ - template - KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedLegendreValues_dx(OutputValueViewType outputValues, Intrepid2::ordinal_type n, ScalarType x, double t) + template + KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedLegendreValues_dx(OutputValueViewType outputValues, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t) { if (n >= 0) outputValues(0) = 0.0; if (n >= 1) outputValues(1) = 1.0; @@ -325,8 +326,8 @@ namespace Intrepid2 This implementation uses less memory than the one below, but depending on the application may introduce some extra computation, in the form of a call to shiftedScaledLegendreValues(). */ - template - KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedLegendreValues_dt(OutputValueViewType outputValues, Intrepid2::ordinal_type n, ScalarType x, double t) + template + KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedLegendreValues_dt(OutputValueViewType outputValues, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t) { // memory-conserving version -- place the Legendre values in the final output container shiftedScaledLegendreValues(outputValues, n, x, t); @@ -358,9 +359,9 @@ namespace Intrepid2 This implementation uses more memory than the one above, but depending on the application may save some computation, in that it can reuse previously computed shiftedScaledLegendreValues. */ - template + template KOKKOS_INLINE_FUNCTION void shiftedScaledIntegratedLegendreValues_dt(OutputValueViewType outputValues, const OutputValueViewType shiftedScaledLegendreValues, - Intrepid2::ordinal_type n, ScalarType x, double t) + Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t) { // reduced flops version: rely on previously computed P_i if (n >= 0) outputValues(0) = 0.0; @@ -387,11 +388,11 @@ namespace Intrepid2 When alpha = 0, Jacobi coincides with Legendre. */ - template - KOKKOS_INLINE_FUNCTION void shiftedScaledJacobiValues(OutputValueViewType outputValues, double alpha, Intrepid2::ordinal_type n, ScalarType x, double t) + template + KOKKOS_INLINE_FUNCTION void shiftedScaledJacobiValues(OutputValueViewType outputValues, double alpha, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t) { ScalarType two_x_minus_t = 2. * x - t; - double alpha_squared_t = alpha * alpha * t; + ScalarTypeForScaling alpha_squared_t = alpha * alpha * t; if (n >= 0) outputValues(0) = 1.0; if (n >= 1) outputValues(1) = two_x_minus_t + alpha * x; @@ -406,7 +407,7 @@ namespace Intrepid2 double c_i = (2. * i + alpha) * (2. * i + alpha - 2.); double d_i = 2. * (i + alpha - 1.) * (i - 1.) * (2. * i + alpha); - outputValues(i) = (b_i / a_i) * (c_i * two_x_minus_t + alpha_squared_t) * P_i_minus_one - (d_i / a_i) * P_i_minus_two; + outputValues(i) = (b_i / a_i) * (c_i * two_x_minus_t + alpha_squared_t) * P_i_minus_one - (d_i / a_i) * t * t * P_i_minus_two; } } @@ -427,9 +428,9 @@ namespace Intrepid2 Compared with the integratedJacobiValues() below, this version uses more memory, but may require fewer floating point computations by reusing the values in jacobiValues. */ - template + template KOKKOS_INLINE_FUNCTION void integratedJacobiValues(OutputValueViewType outputValues, const OutputValueViewType jacobiValues, - double alpha, Intrepid2::ordinal_type n, ScalarType x, double t) + double alpha, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t) { // reduced flops version: rely on previously computed P_i if (n >= 0) outputValues(0) = 1.0; @@ -446,7 +447,7 @@ namespace Intrepid2 double b_i = alpha / ((2. * i + alpha - 2.) * (2. * i + alpha )); double c_i = (i - 1.) / ((2. * i + alpha - 2.) * (2. * i + alpha - 1.)); - outputValues(i) = a_i * P_i + b_i * P_i_minus_1 + c_i * t_squared * P_i_minus_2; + outputValues(i) = a_i * P_i + b_i * t * P_i_minus_1 - c_i * t_squared * P_i_minus_2; } } @@ -467,9 +468,9 @@ namespace Intrepid2 Compared with the integratedJacobiValues() above, this version uses less memory, but may require more floating point computations. */ - template + template KOKKOS_INLINE_FUNCTION void integratedJacobiValues(OutputValueViewType outputValues, - double alpha, Intrepid2::ordinal_type n, ScalarType x, double t) + double alpha, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t) { // memory-conserving version -- place the Jacobi values in the final output container shiftedScaledJacobiValues(outputValues, alpha, n, x, t); @@ -489,7 +490,7 @@ namespace Intrepid2 double b_i = alpha / ((2. * i + alpha - 2.) * (2. * i + alpha )); double c_i = (i - 1.) / ((2. * i + alpha - 2.) * (2. * i + alpha - 1.)); - ScalarType L_i = a_i * P_i + b_i * P_i_minus_1 + c_i * t_squared * P_i_minus_2; + ScalarType L_i = a_i * P_i + b_i * t * P_i_minus_1 - c_i * t_squared * P_i_minus_2; P_i_minus_2 = P_i_minus_1; P_i_minus_1 = P_i; @@ -507,9 +508,9 @@ namespace Intrepid2 \param [in] t - scaling parameter These are defined for x in [0,1]. The x derivative of integrated Jacobi is just Jacobi; the only distinction is in the index -- outputValues indices are shifted by 1 relative to shiftedScaledJacobiValues, above. */ - template + template KOKKOS_INLINE_FUNCTION void integratedJacobiValues_dx(OutputValueViewType outputValues, - double alpha, Intrepid2::ordinal_type n, ScalarType x, double t) + double alpha, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t) { // rather than repeating the somewhat involved implementation of jacobiValues here, // call with (n-1), and then move values accordingly @@ -537,9 +538,9 @@ namespace Intrepid2 This implementation uses more memory than the one above, but depending on the application may save some computation, in that it can reuse previously computed jacobiValues. */ - template + template KOKKOS_INLINE_FUNCTION void integratedJacobiValues_dt(OutputValueViewType outputValues, const OutputValueViewType jacobiValues, - double alpha, Intrepid2::ordinal_type n, ScalarType x, double t) + double alpha, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t) { // reduced flops version: rely on previously computed P_i if (n >= 0) outputValues(0) = 0.0; @@ -562,9 +563,9 @@ namespace Intrepid2 This implementation requires less memory than the one above, but depending on the application may require some extra computation. */ - template + template KOKKOS_INLINE_FUNCTION void integratedJacobiValues_dt(OutputValueViewType outputValues, - double alpha, Intrepid2::ordinal_type n, ScalarType x, double t) + double alpha, Intrepid2::ordinal_type n, ScalarType x, ScalarTypeForScaling t) { // memory-conserving version -- place the Jacobi values in the final output container shiftedScaledJacobiValues(outputValues, alpha, n, x, t); diff --git a/packages/intrepid2/src/Shared/Intrepid2_TestUtils.hpp b/packages/intrepid2/src/Shared/Intrepid2_TestUtils.hpp index c9d4e43d2ba2..c575b1c70b3c 100644 --- a/packages/intrepid2/src/Shared/Intrepid2_TestUtils.hpp +++ b/packages/intrepid2/src/Shared/Intrepid2_TestUtils.hpp @@ -126,18 +126,28 @@ inline Teuchos::RCP< Intrepid2::Basis::key) + if (cellTopo.getBaseKey() == shards::Line<>::key) { basis = getLineBasis(fs,polyOrder_x); } - else if (cellTopo.getBaseKey() == shards::Quadrilateral<4>::key) + else if (cellTopo.getBaseKey() == shards::Quadrilateral<>::key) { + INTREPID2_TEST_FOR_EXCEPTION(polyOrder_y < 0, std::invalid_argument, "polyOrder_y must be specified"); basis = getQuadrilateralBasis(fs,polyOrder_x,polyOrder_y); } + else if (cellTopo.getBaseKey() == shards::Triangle<>::key) + { + basis = getTriangleBasis(fs,polyOrder_x); + } else if (cellTopo.getBaseKey() == shards::Hexahedron<>::key) { + INTREPID2_TEST_FOR_EXCEPTION(polyOrder_y < 0, std::invalid_argument, "polyOrder_y must be specified"); INTREPID2_TEST_FOR_EXCEPTION(polyOrder_z < 0, std::invalid_argument, "polyOrder_z must be specified"); - basis = getHexahedralBasis(fs,polyOrder_x,polyOrder_y,polyOrder_z); + basis = getHexahedronBasis(fs,polyOrder_x,polyOrder_y,polyOrder_z); + } + else if (cellTopo.getBaseKey() == shards::Tetrahedron<>::key) + { + basis = getTetrahedronBasis(fs, polyOrder_x); } else { @@ -156,8 +166,10 @@ inline Teuchos::RCP< Intrepid2::Basis; using LineBasisVol = Intrepid2::LegendreBasis_HVOL_LINE< ExecSpace, Scalar, Scalar>; + using TriangleBasisFamily = Intrepid2::HierarchicalTriangleBasisFamily; + using TetrahedronBasisFamily = Intrepid2::HierarchicalTetrahedronBasisFamily; - using BasisFamily = DerivedBasisFamily; + using BasisFamily = DerivedBasisFamily; return getBasisUsingFamily(cellTopo, fs, polyOrder_x, polyOrder_y, polyOrder_z); } @@ -177,6 +189,9 @@ inline ViewType getView(const std::string &label, DimArgs... dims) } } +/** \brief Returns a View containing equispaced points on the line. + \param [in] numPointsBase - the number of points that will be defined along each edge. + */ template inline ViewType lineInputPointsView(int numPoints) { @@ -184,11 +199,16 @@ inline ViewType lineInputPointsView(int numPoints) Kokkos::parallel_for(numPoints, KOKKOS_LAMBDA(const int i) { double x_zero_one = i * 1.0 / (numPoints-1); // the x value on [0,1] - inputPoints(i,0) = fromZeroOne(x_zero_one); // standard Intrepid2 element + inputPoints(i,0) = PointValueType(fromZeroOne(x_zero_one)); // standard Intrepid2 element }); return inputPoints; } +/** \brief Returns a View containing equispaced points on the hexahedron. + \param [in] numPointsBase - the number of points that will be defined along each edge. + + The total number of points defined will be a cubic number; if n=numPointsBase, then the point count is n^3. + */ template inline ViewType hexInputPointsView(int numPoints_1D) { @@ -203,15 +223,20 @@ inline ViewType hexInputPointsView(int numPoints_1D) { double z_zero_one = k * 1.0 / (numPoints_1D-1); // the z value on [0,1] int pointOrdinal = (i*numPoints_1D+j)*numPoints_1D+k; - inputPoints(pointOrdinal,0) = fromZeroOne(x_zero_one); // standard Intrepid2 element - inputPoints(pointOrdinal,1) = fromZeroOne(y_zero_one); // standard Intrepid2 element - inputPoints(pointOrdinal,2) = fromZeroOne(z_zero_one); // standard Intrepid2 element + inputPoints(pointOrdinal,0) = PointValueType(fromZeroOne(x_zero_one)); // standard Intrepid2 element + inputPoints(pointOrdinal,1) = PointValueType(fromZeroOne(y_zero_one)); // standard Intrepid2 element + inputPoints(pointOrdinal,2) = PointValueType(fromZeroOne(z_zero_one)); // standard Intrepid2 element } } }); return inputPoints; } +/** \brief Returns a View containing equispaced points on the quadrilateral. + \param [in] numPointsBase - the number of points that will be defined along each edge. + + The total number of points defined will be a square number; if n=numPointsBase, then the point count is n^2. + */ template inline ViewType quadInputPointsView(int numPoints_1D) { @@ -224,8 +249,82 @@ inline ViewType quadInputPointsView(int numPoints_1D) { double y_zero_one = j * 1.0 / (numPoints_1D-1); // the y value on [0,1] int pointOrdinal = i*numPoints_1D+j; - inputPoints(pointOrdinal,0) = fromZeroOne(x_zero_one); // standard Intrepid2 element - inputPoints(pointOrdinal,1) = fromZeroOne(y_zero_one); // standard Intrepid2 element + inputPoints(pointOrdinal,0) = PointValueType(fromZeroOne(x_zero_one)); // standard Intrepid2 element + inputPoints(pointOrdinal,1) = PointValueType(fromZeroOne(y_zero_one)); // standard Intrepid2 element + } + }); + return inputPoints; +} + +/** \brief Returns a View containing regularly-spaced points on the tetrahedron. + \param [in] numPointsBase - the number of points that will be defined along each edge. + + The total number of points defined will be a tetrahedral number; if n=numPointsBase, then the point count is the nth tetrahedral number, given by n*(n+1)*(n+2)/6. + */ +template +inline ViewType tetInputPointsView(int numPointsBase) +{ + const int numPoints = numPointsBase*(numPointsBase+1)*(numPointsBase+2)/6; + ViewType inputPoints = getView("tetrahedron input points",numPoints,3); + Kokkos::parallel_for(numPointsBase, KOKKOS_LAMBDA(const int d0) // d0 generalizes row + { + // the following might be the formula for the nested for loop below, but we need to check this + // for now, we comment this out and do the clearer thing that requires more computation +// const int n = numPointsBase-d0; +// const int pointOrdinalOffset = n*(n+1)*(n+2)/6; + int pointOrdinalOffset = 0; + for (int i=0; i +inline ViewType triInputPointsView(int numPointsBase) +{ + const int numPoints = numPointsBase*(numPointsBase+1)/2; + ViewType inputPoints = getView("triangle input points",numPoints,2); + Kokkos::parallel_for(numPointsBase, KOKKOS_LAMBDA(const int row) + { + int rowPointOrdinalOffset = 0; + for (int i=0; i quadInputPointsView(int numPoints_1D) template inline ViewType getInputPointsView(shards::CellTopology &cellTopo, int numPoints_1D) { - if (cellTopo.getBaseKey() == shards::Line<2>::key) + if (cellTopo.getBaseKey() == shards::Line<>::key) { return lineInputPointsView(numPoints_1D); } - else if (cellTopo.getBaseKey() == shards::Quadrilateral<4>::key) + else if (cellTopo.getBaseKey() == shards::Quadrilateral<>::key) { return quadInputPointsView(numPoints_1D); } - else if (cellTopo.getBaseKey() == shards::Hexahedron<8>::key) + else if (cellTopo.getBaseKey() == shards::Hexahedron<>::key) { return hexInputPointsView(numPoints_1D); } + else if (cellTopo.getBaseKey() == shards::Triangle<>::key) + { + return triInputPointsView(numPoints_1D); + } + else if (cellTopo.getBaseKey() == shards::Tetrahedron<>::key) + { + return tetInputPointsView(numPoints_1D); + } else { INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported Cell Topology"); diff --git a/packages/intrepid2/unit-test/Discretization/Basis/BasisEquivalenceTests.cpp b/packages/intrepid2/unit-test/Discretization/Basis/BasisEquivalenceTests.cpp index a62c2a6ff7c3..5d1959fe702a 100644 --- a/packages/intrepid2/unit-test/Discretization/Basis/BasisEquivalenceTests.cpp +++ b/packages/intrepid2/unit-test/Discretization/Basis/BasisEquivalenceTests.cpp @@ -227,7 +227,7 @@ namespace const double relTol, const double absTol, Teuchos::FancyOStream &out, bool &success) { // first, check that the bases agree on cardinality - TEST_ASSERT(basis1.getCardinality() == basis2.getCardinality()); + TEST_EQUALITY(basis1.getCardinality(), basis2.getCardinality()); const int basisCardinality = basis1.getCardinality(); @@ -246,6 +246,18 @@ namespace ViewType weights = ViewType("quadrature weights 1D ref cell", numRefPoints); quadrature->getCubature(points, weights); + out << "Points being tested:\n"; + for (int pointOrdinal=0; pointOrdinal interiorDim) + { + basis1FirstInterior = basis1.getAllDofOrdinal()(interiorDim, interiorSubcellOrdinal, firstDofOrdinalForSubcell); + } + if (basis2.getAllDofOrdinal().extent_int(0) > interiorDim) + { + basis2FirstInterior = basis2.getAllDofOrdinal()(interiorDim, interiorSubcellOrdinal, firstDofOrdinalForSubcell); + } + + // if there are no interior dofs, we'll get a -1 value + if (basis1FirstInterior != -1) + { + // at index 3, dof tag stores the total dof count associated with the subcell (here, the interior) + basis1InteriorCardinality = basis1.getDofTag(basis1FirstInterior)(3); + } + if (basis2FirstInterior != -1) + { + // at index 3, dof tag stores the total dof count associated with the subcell (here, the interior) + basis2InteriorCardinality = basis2.getDofTag(basis2FirstInterior)(3); + } + const bool basis1IsDG = (basis1InteriorCardinality == basisCardinality); + const bool basis2IsDG = (basis2InteriorCardinality == basisCardinality); + const bool neitherBasisIsDG = !basis1IsDG && !basis2IsDG; + if (neitherBasisIsDG) + { + auto basis1CoefficientsHost = getHostCopy(basis1Coefficients); + // if neither basis is DG, then we can expect them to have matching counts on basis members + // associated with a given subcell. We can also expect the representation of of members of basis1 + // on a given subcell in terms of basis2 to involve at least some members of basis2 associated with the same + // subcell. (We can check this by examining the weights in basis1Coefficients.) + for (int subcellDim=0; subcellDim<=interiorDim; subcellDim++) + { + out << "checking subcells of dimension " << subcellDim << std::endl; + // if there are no dofs for this subcell dimension and greater, then getAllDofOrdinal() won't have + // an entry for subcellDim. The following guards against that: + if (basis1.getAllDofOrdinal().extent_int(0) <= subcellDim) + { + // basis1 has no dofs for this subcell dim. Check that basis2 also does not: + TEST_ASSERT(basis2.getAllDofOrdinal().extent_int(0) <= subcellDim); + break; // we've already checked all subcell dimensions that have any dofs associated with them + } + + const int subcellCount = cellTopo.getSubcellCount(subcellDim); + for (int subcellOrdinal=0; subcellOrdinal basis1SubcellDofOrdinals(basis1SubcellCardinality); + std::vector basis2SubcellDofOrdinals(basis2SubcellCardinality); + for (int subcellDofOrdinal=0; subcellDofOrdinal absTol); + if (nonzeroCoefficient) + { + out << "basis coefficient " << basisCoefficient << " is nonzero (absTol = " << absTol << ").\n"; + hasNonzeroCoefficient = true; + } + } + TEST_ASSERT(hasNonzeroCoefficient); + } + } + } + } + // compute the values for basis2 using basis1Coefficients, basis1Values, and confirm that these agree with basisValues2 ViewType basis2ValuesFromBasis1 = getOutputView(functionSpace, OPERATOR_VALUE, basisCardinality, numRefPoints, spaceDim); deviceGeneralizedMatrixMultiply(basis1Coefficients, true, basis1Values, basis2ValuesFromBasis1); // true: transpose @@ -477,4 +594,115 @@ namespace testBasisEquivalence(cnBasis, c2Basis, opsToTest, relTol, absTol, out, success); } + TEUCHOS_UNIT_TEST( BasisEquivalence, TetrahedronNodalCnVersusNodalC2_HGRAD ) + { + using CnBasis = Intrepid2::Basis_HGRAD_TET_Cn_FEM; + using C2Basis = Intrepid2::Basis_HGRAD_TET_C2_FEM; + + // OPERATOR_D2 and above are not supported by either the nodal or the hierarchical basis at present... + std::vector opsToTest {OPERATOR_GRAD, OPERATOR_D1}; + + // these tolerances are selected such that we have a little leeway for architectural differences + // (It is true, though, that we incur a fair amount of floating point error for higher order bases in higher dimensions) + const double relTol=1e-12; // 2e-14 is sharp on development setup; relaxing for potential architectural differences + const double absTol=1e-13; // 3e-15 is sharp on development setup; relaxing for potential architectural differences + + CnBasis cnBasis(2); + C2Basis c2Basis; + testBasisEquivalence(cnBasis, c2Basis, opsToTest, relTol, absTol, out, success); + } + + TEUCHOS_UNIT_TEST( BasisEquivalence, TetrahedronHierarchicalDGVersusHierarchicalCG_HGRAD ) + { + using CGBasis = HierarchicalBasisFamily<>::HGRAD_TET; + using DGBasis = DGHierarchicalBasisFamily<>::HGRAD_TET; + + // these tolerances are selected such that we have a little leeway for architectural differences + // (It is true, though, that we incur a fair amount of floating point error for higher order bases in higher dimensions) + const double relTol=1e-6; // 3e-8 is sharp on development setup for polyOrder=2; relaxing for potential architectural differences + const double absTol=1e-9; // 5e-11 is sharp on development setup for polyOrder=2; relaxing for potential architectural differences + + std::vector opsToTest {OPERATOR_GRAD, OPERATOR_D1}; + for (int polyOrder=1; polyOrder<7; polyOrder++) + { + CGBasis cgBasis(polyOrder); + DGBasis dgBasis(polyOrder); + testBasisEquivalence(cgBasis, dgBasis, opsToTest, relTol, absTol, out, success); + } + } + + TEUCHOS_UNIT_TEST( BasisEquivalence, TetrahedronNodalVersusHierarchicalCG_HGRAD ) + { + using HierarchicalBasis = HierarchicalBasisFamily<>::HGRAD_TET; + using NodalBasis = NodalBasisFamily<>::HGRAD_TET; + + // OPERATOR_D2 and above are not supported by either the nodal or the hierarchical basis at present... + std::vector opsToTest {OPERATOR_GRAD, OPERATOR_D1}; + + // these tolerances are selected such that we have a little leeway for architectural differences + // (It is true, though, that we incur a fair amount of floating point error for higher order bases in higher dimensions) + const double relTol=1e-6; // 3e-08 is sharp on development setup for polyOrder=6; relaxing for potential architectural differences + const double absTol=1e-10; // 3e-12 is sharp on development setup for polyOrder=6; relaxing for potential architectural differences + + for (int polyOrder=1; polyOrder<7; polyOrder++) + { + HierarchicalBasis hierarchicalBasis(polyOrder); + NodalBasis nodalBasis(polyOrder); + +// auto cellTopo = nodalBasis.getBaseCellTopology(); +// const int faceDim = 2; +// for (int intrepid2FaceOrdinal=0; intrepid2FaceOrdinal<4; intrepid2FaceOrdinal++) +// { +// int vertex0 = cellTopo.getNodeMap(faceDim, intrepid2FaceOrdinal, 0); +// int vertex1 = cellTopo.getNodeMap(faceDim, intrepid2FaceOrdinal, 1); +// int vertex2 = cellTopo.getNodeMap(faceDim, intrepid2FaceOrdinal, 2); +// std::cout << "face " << intrepid2FaceOrdinal << ": " << vertex0 << vertex1 << vertex2 << std::endl; +// } + + testBasisEquivalence(nodalBasis, hierarchicalBasis, opsToTest, relTol, absTol, out, success); + } + } + + TEUCHOS_UNIT_TEST( BasisEquivalence, TriangleNodalVersusHierarchicalCG_HGRAD ) + { + using HierarchicalBasis = HierarchicalBasisFamily<>::HGRAD_TRI; + using NodalBasis = NodalBasisFamily<>::HGRAD_TRI; + + // OPERATOR_D2 and above are not supported by either the nodal or the hierarchical basis at present... + std::vector opsToTest {OPERATOR_GRAD, OPERATOR_D1}; + + // these tolerances are selected such that we have a little leeway for architectural differences + // (It is true, though, that we incur a fair amount of floating point error for higher order bases in higher dimensions) + const double relTol=1e-11; // 7e-13 is sharp on development setup for polyOrder=4; relaxing for potential architectural differences + const double absTol=1e-12; // 2e-14 is sharp on development setup for polyOrder=4; relaxing for potential architectural differences + + for (int polyOrder=1; polyOrder<5; polyOrder++) + { + HierarchicalBasis hierarchicalBasis(polyOrder); + NodalBasis nodalBasis(polyOrder); + testBasisEquivalence(nodalBasis, hierarchicalBasis, opsToTest, relTol, absTol, out, success); + } + } + + TEUCHOS_UNIT_TEST( BasisEquivalence, TriangleHierarchicalCGVersusHierarchicalDG_HGRAD ) + { + using CGBasis = HierarchicalBasisFamily<>::HGRAD_TRI; + using DGBasis = DGHierarchicalBasisFamily<>::HGRAD_TRI; + + // OPERATOR_D2 and above are not supported by either the nodal or the hierarchical basis at present... + std::vector opsToTest {OPERATOR_GRAD, OPERATOR_D1}; + + // these tolerances are selected such that we have a little leeway for architectural differences + // (It is true, though, that we incur a fair amount of floating point error for higher order bases in higher dimensions) + const double relTol=1e-11; // 7e-13 is sharp on development setup for polyOrder=4; relaxing for potential architectural differences + const double absTol=1e-12; // 2e-14 is sharp on development setup for polyOrder=4; relaxing for potential architectural differences + + for (int polyOrder=1; polyOrder<5; polyOrder++) + { + CGBasis cgBasis(polyOrder); + DGBasis dgBasis(polyOrder); + testBasisEquivalence(cgBasis, dgBasis, opsToTest, relTol, absTol, out, success); + } + } + } // namespace diff --git a/packages/intrepid2/unit-test/Discretization/Basis/HierarchicalBases/AnalyticPolynomialsMatchTests.cpp b/packages/intrepid2/unit-test/Discretization/Basis/HierarchicalBases/AnalyticPolynomialsMatchTests.cpp index 528e3a3f5f46..128c969d20ab 100644 --- a/packages/intrepid2/unit-test/Discretization/Basis/HierarchicalBases/AnalyticPolynomialsMatchTests.cpp +++ b/packages/intrepid2/unit-test/Discretization/Basis/HierarchicalBases/AnalyticPolynomialsMatchTests.cpp @@ -52,6 +52,7 @@ #include "Intrepid2_HierarchicalBasisFamily.hpp" #include "Intrepid2_NodalBasisFamily.hpp" +#include "Intrepid2_Polynomials.hpp" #include "Intrepid2_Types.hpp" #include "Intrepid2_TestUtils.hpp" @@ -62,6 +63,174 @@ namespace { using namespace Intrepid2; + template + Scalar integratedJacobi(Scalar x, Scalar t, double alpha, const int n, const int derivativeOrder = 0) + { + // ideally, we would have closed-form expressions somewhere in here, as we do for integrated Legendre below + // for now, though, this is just a thin wrapper around the call to Intrepid2::Polynomials::integratedJacobiValues() + auto values = getView("integrated Jacobi values", n+1); + Polynomials::integratedJacobiValues(values, alpha, n, x, t); + return values(n); + } + + template + Scalar integratedLegendreAnalytic(Scalar x, const int n, bool useMinusOneToOne, const int derivativeOrder = 0) + { + // formulas below are for x in [-1,1]; if we are using [0,1], need to remap appropriately: + const double derivativeScaling = useMinusOneToOne ? 1.0 : pow(2.0, derivativeOrder); + if (!useMinusOneToOne) + { + x = 2.0 * x - 1.0; + } + Scalar value; + switch (derivativeOrder) + { + case 0: + switch (n) + { + case 0: + value = (1.0-x)/2.0; // left vertex function (node at -1) + break; + case 1: + value = (1.0+x)/2.0; // right vertex function (node at 1) + break; + case 2: + value = (x*x-1.0)/4.0; // L_2 : (x^2 - 1) / 4 + break; + case 3: + value = (x*x*x-x)/4.0; // L_3 : (x^3 - x) / 4 + break; + case 4: + value = (5.*x*x*x*x - 6.*x*x+1.)/16.0; // L_4 : (5x^4-6x^2+1) / 16 + break; + default: + TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "unsupported n"); + } + break; + case 1: + switch (n) + { + case 0: + value = -1.0/2.0; // left vertex function (node at -1) + break; + case 1: + value = 1.0/2.0; // right vertex function (node at 1) + break; + case 2: + value = x/2.0; // L_2 : (x^2 - 1) / 4 + break; + case 3: + value = (3.0*x*x-1.0)/4.0; // L_3 : (x^3 - x) / 4 + break; + case 4: + value = (5.*x*x*x - 3.*x)/4.0; // L_4 : (5x^4-6x^2+1) / 16 + break; + default: + TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "unsupported n"); + } + break; + case 2: + switch (n) + { + case 0: + value = 0.0; // left vertex function (node at -1) + break; + case 1: + value = 0.0; // right vertex function (node at 1) + break; + case 2: + value = 1.0/2.0; // L_2 : (x^2 - 1) / 4 + break; + case 3: + value = 3.0*x/2.0; // L_3 : (x^3 - x) / 4 + break; + case 4: + value = (15.*x*x - 3.)/4.; // L_4 : (5x^4-6x^2+1) / 16 + break; + default: + TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "unsupported n"); + } + break; + case 3: + switch (n) + { + case 0: + value = 0.0; // left vertex function (node at -1) + break; + case 1: + value = 0.0; // right vertex function (node at 1) + break; + case 2: + value = 0.0; // L_2 : (x^2 - 1) / 4 + break; + case 3: + value = 3.0/2.0; // L_3 : (x^3 - x) / 4 + break; + case 4: + value = (15.*x)/2.; // L_4 : (5x^4-6x^2+1) / 16 + break; + default: + TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "unsupported n"); + } + break; + case 4: + switch (n) + { + case 0: + value = 0.0; // left vertex function (node at -1) + break; + case 1: + value = 0.0; // right vertex function (node at 1) + break; + case 2: + value = 0.0; // L_2 : (x^2 - 1) / 4 + break; + case 3: + value = 0.0; // L_3 : (x^3 - x) / 4 + break; + case 4: + value = 15./2.; // L_4 : (5x^4-6x^2+1) / 16 + break; + default: + TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "unsupported n"); + } + break; + default: + value = 0.0; + } + return value * derivativeScaling; + } + + template + Scalar shiftedScaledIntegratedLegendreAnalytic(Scalar x, Scalar t, const int n) + { + bool useMinusOneToOne = false; // our domain is [0,1] + const int derivativeOrder = 0; + double tol = 1e-15; + using std::abs; + if ((n == 0) || (n == 1)) + { + // linear in x: scaling by t cancels + return integratedLegendreAnalytic(x, n, useMinusOneToOne, derivativeOrder); + } + else if (abs(t) < tol) + { + // define a scaling by 0 to be 0 (does this match what Fuentes et al do??) + // even if this isn't perfectly general, it works for our tests (at least for OPERATOR_VALUE): + // we only ever will end up with t=0 for edge functions on the triangle, evaluated at a vertex. These vanish. + return 0.0; + } + else + { + Scalar tPower = 1.0; + for (int i=0; i @@ -92,64 +261,36 @@ namespace int pointPassed = true; PointScalar x = inputPointsViewHost(pointOrdinal,0); + const bool useMinusOneToOne = true; // Intrepid2's reference element + int derivativeOrder; + switch (op) { case Intrepid2::OPERATOR_VALUE: - expectedValuesViewHost(0) = (1.0-x)/2.0; // left vertex function (node at -1) - expectedValuesViewHost(1) = (1.0+x)/2.0; // right vertex function (node at 1) - expectedValuesViewHost(2) = (x*x-1.0)/4.0; // L_2 : (x^2 - 1) / 4 - expectedValuesViewHost(3) = (x*x*x-x)/4.0; // L_3 : (x^3 - x) / 4 - expectedValuesViewHost(4) = (5.*x*x*x*x - 6.*x*x+1.)/16.0; // L_4 : (5x^4-6x^2+1) / 16 + derivativeOrder = 0; break; case Intrepid2::OPERATOR_GRAD: - case Intrepid2::OPERATOR_D1: - // first derivatives of the above: - expectedValuesViewHost(0) = -1.0/2.0; // left vertex function (node at -1) - expectedValuesViewHost(1) = 1.0/2.0; // right vertex function (node at 1) - expectedValuesViewHost(2) = x/2.0; // L_2 : (x^2 - 1) / 4 - expectedValuesViewHost(3) = (3.0*x*x-1.0)/4.0; // L_3 : (x^3 - x) / 4 - expectedValuesViewHost(4) = (5.*x*x*x - 3.*x)/4.0; // L_4 : (5x^4-6x^2+1) / 16 + derivativeOrder = 1; break; + case Intrepid2::OPERATOR_D1: case Intrepid2::OPERATOR_D2: - // second derivatives: - expectedValuesViewHost(0) = 0.0; // left vertex function (node at -1) - expectedValuesViewHost(1) = 0.0; // right vertex function (node at 1) - expectedValuesViewHost(2) = 1.0/2.0; // L_2 : (x^2 - 1) / 4 - expectedValuesViewHost(3) = 3.0*x/2.0; // L_3 : (x^3 - x) / 4 - expectedValuesViewHost(4) = (15.*x*x - 3.)/4.; // L_4 : (5x^4-6x^2+1) / 16 - break; case Intrepid2::OPERATOR_D3: - // third derivatives: - expectedValuesViewHost(0) = 0.0; // left vertex function (node at -1) - expectedValuesViewHost(1) = 0.0; // right vertex function (node at 1) - expectedValuesViewHost(2) = 0.0; // L_2 : (x^2 - 1) / 4 - expectedValuesViewHost(3) = 3.0/2.0; // L_3 : (x^3 - x) / 4 - expectedValuesViewHost(4) = (15.*x)/2.; // L_4 : (5x^4-6x^2+1) / 16 - break; case Intrepid2::OPERATOR_D4: - // fourth derivatives: - expectedValuesViewHost(0) = 0.0; // left vertex function (node at -1) - expectedValuesViewHost(1) = 0.0; // right vertex function (node at 1) - expectedValuesViewHost(2) = 0.0; // L_2 : (x^2 - 1) / 4 - expectedValuesViewHost(3) = 0.0; // L_3 : (x^3 - x) / 4 - expectedValuesViewHost(4) = 15./2.; // L_4 : (5x^4-6x^2+1) / 16 - break; case Intrepid2::OPERATOR_D5: case Intrepid2::OPERATOR_D6: case Intrepid2::OPERATOR_D7: case Intrepid2::OPERATOR_D8: case Intrepid2::OPERATOR_D9: case Intrepid2::OPERATOR_D10: - // fourth derivatives: - expectedValuesViewHost(0) = 0.0; // left vertex function (node at -1) - expectedValuesViewHost(1) = 0.0; // right vertex function (node at 1) - expectedValuesViewHost(2) = 0.0; // L_2 : (x^2 - 1) / 4 - expectedValuesViewHost(3) = 0.0; // L_3 : (x^3 - x) / 4 - expectedValuesViewHost(4) = 0.0; // L_4 : (5x^4-6x^2+1) / 16 + derivativeOrder = op - OPERATOR_D1 + 1; break; default: INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported operator"); } + for (int n=0; ngetCardinality(); fieldOrdinal++) { @@ -339,6 +480,208 @@ namespace } } + template + void testHierarchicalHGRAD_TRIANGLE_MatchesAnalyticValues(Intrepid2::EOperator op, const double tol, Teuchos::FancyOStream &out, bool &success) + { + using namespace Intrepid2; + using BasisFamily = HierarchicalBasisFamily; + + const int spaceDim = 2; + const int polyOrder = 4; + auto hgradBasis = getTriangleBasis(FUNCTION_SPACE_HGRAD, polyOrder); + + const int numVertices = 3; + const int numFunctionsPerVertex = 1; + const int numVertexFunctionsExpected = numVertices * numFunctionsPerVertex; + const int numEdges = 3; + const int num1DEdgeFunctions = (polyOrder + 1) - 2; // line basis cardinality, minus two vertex functions + const int numEdgeFunctionsExpected = num1DEdgeFunctions * numEdges; + const int numInteriorFunctionsExpected = (num1DEdgeFunctions-1)*num1DEdgeFunctions/2; // triangular sum + const int expectedCardinality = numVertexFunctionsExpected + numEdgeFunctionsExpected + numInteriorFunctionsExpected; + + TEST_EQUALITY(expectedCardinality, hgradBasis->getCardinality()); + // could also test vertex/edge/face function count individually (worth doing, I think) + + int numPoints_1D = 5; + shards::CellTopology triangleTopo = shards::CellTopology(shards::getCellTopologyData >() ); + auto inputPointsView = getInputPointsView(triangleTopo, numPoints_1D); + + const int numPoints = inputPointsView.extent_int(0); + + auto hgradOutputView = getOutputView(FUNCTION_SPACE_HGRAD, op, hgradBasis->getCardinality(), numPoints, spaceDim); + hgradBasis->getValues(hgradOutputView, inputPointsView, op); + + auto hgradOutputViewHost = getHostCopy(hgradOutputView); + auto inputPointsViewHost = getHostCopy(inputPointsView); + + ViewType expectedValuesView = getOutputView(FUNCTION_SPACE_HGRAD, op, hgradBasis->getCardinality(), numPoints, spaceDim); + auto expectedValuesViewHost = getHostCopy(expectedValuesView); + + auto legendreValuesAtPoint = getView("Legendre values temporary storage", polyOrder+1); + auto legendreValuesAtPointHost = getHostCopy(legendreValuesAtPoint); + + for (int pointOrdinal=0; pointOrdinalgetCardinality(); fieldOrdinal++) + { + for (int d=0; d void testDerivativesMatch(int polyOrder, int derivativeOrder, const double tol, Teuchos::FancyOStream &out, bool &success) @@ -460,8 +803,8 @@ namespace } else if (cellTopo.getKey() == shards::Hexahedron<>::key) { - derivedBasis = getHexahedralBasis (fs, polyOrder); // derived basis supports both isotropic and anisotropic polyOrder - standardBasis = getHexahedralBasis(fs, polyOrder); // isotropic + derivedBasis = getHexahedronBasis (fs, polyOrder); // derived basis supports both isotropic and anisotropic polyOrder + standardBasis = getHexahedronBasis(fs, polyOrder); // isotropic } int standardCardinality = standardBasis->getCardinality(); @@ -694,6 +1037,19 @@ namespace } } + TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL( AnalyticPolynomialsMatch, Hierarchical_HGRAD_TRI, OutputScalar, PointScalar ) + { + const double tol = TEST_TOLERANCE_TIGHT; + using ExecSpace = Kokkos::DefaultExecutionSpace; + +// std::vector operators = {{OPERATOR_VALUE, OPERATOR_GRAD, OPERATOR_D1, OPERATOR_D2, OPERATOR_D3, OPERATOR_D4, OPERATOR_D5, OPERATOR_D6, OPERATOR_D7, OPERATOR_D8, OPERATOR_D9, OPERATOR_D10}}; + std::vector operators = {OPERATOR_VALUE}; + for (auto op : operators) + { + testHierarchicalHGRAD_TRIANGLE_MatchesAnalyticValues(op, tol, out, success); + } + } + TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL( AnalyticPolynomialsMatch, Hierarchical_LineBasisDerivativesAgree, OutputScalar, PointScalar ) { const int maxPolyOrder = 10; @@ -741,7 +1097,7 @@ namespace runNodalBasisComparisonTests(polyOrder_2D, quadTopo, {FUNCTION_SPACE_HGRAD}, {OPERATOR_VALUE,OPERATOR_GRAD}, tol, out, success); runNodalBasisComparisonTests(polyOrder_2D, quadTopo, {FUNCTION_SPACE_HGRAD}, operators_dk, tol, out, success); - out << "Running 3D nodal hexahedral basis comparison tests…\n"; + out << "Running 3D nodal hexahedron basis comparison tests…\n"; runNodalBasisComparisonTests(polyOrder_3D, hexTopo, {FUNCTION_SPACE_HVOL}, {OPERATOR_VALUE}, tol, out, success); runNodalBasisComparisonTests(polyOrder_3D, hexTopo, {FUNCTION_SPACE_HGRAD}, {OPERATOR_VALUE,OPERATOR_GRAD}, tol, out, success); runNodalBasisComparisonTests(polyOrder_3D, hexTopo, {FUNCTION_SPACE_HGRAD}, operators_dk, tol, out, success); @@ -749,5 +1105,6 @@ namespace INTREPID2_OUTPUTSCALAR_POINTSCALAR_TEST_INSTANT( AnalyticPolynomialsMatch, Hierarchical_HGRAD_LINE ) INTREPID2_OUTPUTSCALAR_POINTSCALAR_TEST_INSTANT( AnalyticPolynomialsMatch, Hierarchical_LineBasisDerivativesAgree ) + INTREPID2_OUTPUTSCALAR_POINTSCALAR_TEST_INSTANT( AnalyticPolynomialsMatch, Hierarchical_HGRAD_TRI ) INTREPID2_OUTPUTSCALAR_POINTSCALAR_TEST_INSTANT( AnalyticPolynomialsMatch, HierarchicalNodalComparisons ) } // namespace diff --git a/packages/intrepid2/unit-test/Discretization/Basis/HierarchicalBases/SubBasisInclusionTests.cpp b/packages/intrepid2/unit-test/Discretization/Basis/HierarchicalBases/SubBasisInclusionTests.cpp index b830a60ef97e..eacd9dc80d48 100644 --- a/packages/intrepid2/unit-test/Discretization/Basis/HierarchicalBases/SubBasisInclusionTests.cpp +++ b/packages/intrepid2/unit-test/Discretization/Basis/HierarchicalBases/SubBasisInclusionTests.cpp @@ -90,13 +90,31 @@ namespace int spaceDim = cellTopo.getDimension(); int minDegree = (fs == Intrepid2::FUNCTION_SPACE_HVOL) ? 0 : 1; - auto subBasisDegreeTestCases = getBasisTestCasesUpToDegree(spaceDim, minDegree, polyOrder_x, polyOrder_y, polyOrder_z); + bool isLine = cellTopo.getKey() == shards::Line<>::key; + bool isQuad = cellTopo.getKey() == shards::Quadrilateral<>::key; + bool isHex = cellTopo.getKey() == shards::Hexahedron<>::key; + bool isTri = cellTopo.getKey() == shards::Triangle<>::key; + bool isTet = cellTopo.getKey() == shards::Tetrahedron<>::key; + bool isWedge = cellTopo.getKey() == shards::Wedge<>::key; + int polyOrderDim = -1; // the number of dimensions of p-anisotropy allowed + if (isLine || isQuad || isHex) + { + polyOrderDim = spaceDim; + } + else if (isTri || isTet) + { + polyOrderDim = 1; + } + else if (isWedge) + { + polyOrderDim = 2; // line x tri + } + auto subBasisDegreeTestCases = getBasisTestCasesUpToDegree(polyOrderDim, minDegree, polyOrder_x, polyOrder_y, polyOrder_z); - std::vector degrees(spaceDim); + std::vector degrees(polyOrderDim); degrees[0] = polyOrder_x; - if (spaceDim > 1) degrees[1] = polyOrder_y; - if (spaceDim > 2) degrees[2] = polyOrder_z; - // // TODO: consider what to do when non-hypercubes are being tested + if (polyOrderDim > 1) degrees[1] = polyOrder_y; + if (polyOrderDim > 2) degrees[2] = polyOrder_z; int numPoints_1D = 5; auto inputPoints = getInputPointsView(cellTopo, numPoints_1D); @@ -106,9 +124,23 @@ namespace auto outputValues = getOutputView(fs, op, basis->getCardinality(), numPoints, spaceDim); basis->getValues(outputValues, inputPoints, op); + out << "Testing sub-basis inclusion in degree "; + for (int i=0; igetName() << std::endl; for (auto testCase : subBasisDegreeTestCases) { + out << "testing sub-basis of degree "; + for (int i=0; i >() ); + + // so far, only HGRAD implemented for hierarchical tetrahedron. Once we have full exact sequence, we can + // switch to calling runSubBasisTests(tetTopo, out, success). + std::vector functionSpaces = {FUNCTION_SPACE_HGRAD}; + auto cellTopo = tetTopo; + + const int maxDegree = 5; + const double tol = TEST_TOLERANCE_TIGHT; + + for (auto fs : functionSpaces) + { + std::vector continuousBasisValues; + if (fs != FUNCTION_SPACE_HVOL) + { + continuousBasisValues = {true,false}; + } + else + { + continuousBasisValues = {true}; // false case not supported by the dof tag stuff that testSubBasis() does + } + for (auto continuousBasis : continuousBasisValues) // corresponds to "defineVertexFunctions" in basis definitions + { + for (int degree=1; degree<=maxDegree; degree++) + { + testSubBasis(cellTopo, fs, tol, out, success, degree, -1, -1,continuousBasis); + } + } + } + } + + TEUCHOS_UNIT_TEST( SubBasisInclusion, Triangle ) + { + shards::CellTopology triTopo = shards::CellTopology(shards::getCellTopologyData >() ); + + // so far, only HGRAD implemented for hierarchical triangle. Once we have full exact sequence, we can + // switch to calling runSubBasisTests(triTopo, out, success). + std::vector functionSpaces = {FUNCTION_SPACE_HGRAD}; + auto cellTopo = triTopo; + + const int maxDegree = 5; + const double tol = TEST_TOLERANCE_TIGHT; + + for (auto fs : functionSpaces) + { + std::vector continuousBasisValues; + if (fs != FUNCTION_SPACE_HVOL) + { + continuousBasisValues = {true,false}; + } + else + { + continuousBasisValues = {true}; // false case not supported by the dof tag stuff that testSubBasis() does + } + for (auto continuousBasis : continuousBasisValues) // corresponds to "defineVertexFunctions" in line basis definitions + { + for (int degree=1; degree<=maxDegree; degree++) + { + testSubBasis(cellTopo, fs, tol, out, success, degree, -1, -1,continuousBasis); + } + } + } + } + } // namespace diff --git a/packages/intrepid2/unit-test/Shared/Polylib/LegendreJacobiPolynomials/IntegratedJacobiTests.cpp b/packages/intrepid2/unit-test/Shared/Polylib/LegendreJacobiPolynomials/IntegratedJacobiTests.cpp index b269ff875a5d..d761b1d44447 100644 --- a/packages/intrepid2/unit-test/Shared/Polylib/LegendreJacobiPolynomials/IntegratedJacobiTests.cpp +++ b/packages/intrepid2/unit-test/Shared/Polylib/LegendreJacobiPolynomials/IntegratedJacobiTests.cpp @@ -64,6 +64,73 @@ namespace std::vector t_values = {{0.0, 0.2, 0.4, 0.6, 0.8, 1.0}}; std::vector x_values = {{-1.0,-0.5,-1.0/3.0,0.0,1.0/3.0,0.50,1.0}}; + void testIntegratedJacobiIsZeroAtZero(const int polyOrder, const double tol, Teuchos::FancyOStream &out, bool &success) + { + // for all alpha, t, integrated jacobi should evaluate to 0 at 0. + Kokkos::View integratedJacobiView("integrated jacobi values",polyOrder+1); + + using Intrepid2::Polynomials::integratedJacobiValues; + using Intrepid2::Polynomials::shiftedScaledJacobiValues; + + const double x = 0.0; + + for (auto alpha : alpha_values) + { + for (auto t : t_values) + { + // wrap invocation in parallel_for just to ensure execution on device (for CUDA) + Kokkos::parallel_for(1, KOKKOS_LAMBDA(const int dummy_index) + { + integratedJacobiValues(integratedJacobiView, alpha, polyOrder, x, t); + }); + + for (int i=1; i<=polyOrder; i++) + { + if ( abs(integratedJacobiView(i)) > tol) + { + success = false; + out << "for alpha = " << alpha << ", t = " << t << ", integrated Jacobi for polyOrder " << i; + out << " at x=0 is not zero (it is " << integratedJacobiView(i) << ")\n"; + } + } + } + } + } + + void testIntegratedJacobiAnalyticAlphaZeroJ2(const double tol, Teuchos::FancyOStream &out, bool &success) + { + // for alpha=0, integrated jacobi j=2 should be x^2 - x*t + const int polyOrder = 2; + Kokkos::View integratedJacobiView("integrated jacobi values",polyOrder+1); + + using Intrepid2::Polynomials::integratedJacobiValues; + + const double alpha = 0.0; + + for (auto x : x_values) + { + for (auto t : t_values) + { + double expected_value = x * x - x * t; + // wrap invocation in parallel_for just to ensure execution on device (for CUDA) + Kokkos::parallel_for(1, KOKKOS_LAMBDA(const int dummy_index) + { + integratedJacobiValues(integratedJacobiView, alpha, polyOrder, x, t); + }); + + const int i = 2; + double diff = integratedJacobiView(i) - expected_value; + if ( abs(diff) > tol) + { + success = false; + out << "for alpha = " << alpha << ", t = " << t << ", integrated Jacobi for polyOrder " << i; + out << " at x=" << x << " is not x^2 - x * t (" << expected_value << "); instead, it is " << integratedJacobiView(i); + out << ", a difference of " << abs(diff) << ")\n"; + } + } + } + } + void testIntegratedJacobiTwoPathsMatch(const int polyOrder, const double tol, Teuchos::FancyOStream &out, bool &success) { Kokkos::View jacobiView("jacobi values",polyOrder+1); @@ -170,4 +237,18 @@ namespace testIntegratedJacobi_dtTwoPathsMatch(polyOrderMax, tol, out, success); } + TEUCHOS_UNIT_TEST( IntegratedJacobi, ZeroAtZero) + { + const int polyOrderMax = 10; + const double tol = TEST_TOLERANCE_TIGHT; + testIntegratedJacobiIsZeroAtZero(polyOrderMax, tol, out, success); + } + + TEUCHOS_UNIT_TEST( IntegratedJacobi, AnalyticAlphaZeroJ2) + { + // analytically, we can determine that with alpha=0, the second integrated Jacobi polynomial + // should be x^2 - x*t + const double tol = TEST_TOLERANCE_TIGHT; + testIntegratedJacobiAnalyticAlphaZeroJ2(tol, out, success); + } } // namespace diff --git a/packages/intrepid2/unit-test/Shared/Polylib/LegendreJacobiPolynomials/JacobiTests.cpp b/packages/intrepid2/unit-test/Shared/Polylib/LegendreJacobiPolynomials/JacobiTests.cpp index 10b2266821f2..66412588f5e2 100644 --- a/packages/intrepid2/unit-test/Shared/Polylib/LegendreJacobiPolynomials/JacobiTests.cpp +++ b/packages/intrepid2/unit-test/Shared/Polylib/LegendreJacobiPolynomials/JacobiTests.cpp @@ -60,6 +60,7 @@ namespace using namespace Intrepid2; std::vector alpha_values = {{0.0, 0.2, 0.4, 0.6, 0.8, 1.0}}; std::vector x_values = {{-1.0,-0.5,-1.0/3.0,0.0,1.0/3.0,0.50,1.0}}; + std::vector t_values = {{0.0, 0.2, 0.4, 0.6, 0.8, 1.0}}; void testAgreementWithPolylib(const int polyOrder, const double tol, Teuchos::FancyOStream &out, bool &success) { @@ -121,6 +122,54 @@ namespace } } + void testAgreementWithAnalyticLegendreForP2(const double tol, Teuchos::FancyOStream &out, bool &success) + { + const double alpha = 0; + const int polyOrder = 2; + Kokkos::View jacobiValues("jacobi values from Intrepid2::Polynomials",polyOrder+1); + + for (auto t : t_values) + { + for (auto x_onMinusOneToOne : x_values) + { + double x = (x_onMinusOneToOne + 1.0) / 2.0; // shiftedScaledJacobiValues computes on a [0,1] domain, compared with a [-1,1] domain + + using Intrepid2::Polynomials::shiftedScaledJacobiValues; + + // wrap invocation in parallel_for just to ensure execution on device (for CUDA) + Kokkos::parallel_for(1, KOKKOS_LAMBDA(const int dummyIndex) + { + shiftedScaledJacobiValues(jacobiValues, alpha, polyOrder, x, t); + }); + + auto jacobiValuesHost = getHostCopy(jacobiValues); + + const double expectedValue = 6.0 * x * x - 6.0 * x * t + t*t; + + const int i = 2; + bool valuesAreBothSmall = valuesAreSmall(jacobiValuesHost(i), expectedValue, tol); + + if (!valuesAreBothSmall) + { + if (! approximatelyEqual(jacobiValuesHost(i), expectedValue, tol) ) + { + out << "for polyOrder " << i << ", alpha = " << alpha << ", x = " << x << ", t = " << t << ": "; + out << jacobiValuesHost(i) << " != " << expectedValue; + out << " (diff = " << abs(jacobiValuesHost(i) - expectedValue); + out << "; tol = " << tol << ")\n"; + success = false; + } + } + } + } + } + + TEUCHOS_UNIT_TEST( Jacobi, AgreesWithLegendreForP2 ) + { + const double tol = TEST_TOLERANCE_TIGHT; + testAgreementWithAnalyticLegendreForP2(tol, out, success); + } + TEUCHOS_UNIT_TEST( Jacobi, AgreesWithPolylib ) { const int polyOrderMax = 1;