Intrepid2_Basis.hpp

Go to the documentation of this file.

// @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

#ifndef __INTREPID2_BASIS_HPP__
#define __INTREPID2_BASIS_HPP__

#include "Intrepid2_ConfigDefs.hpp"
#include "Intrepid2_Types.hpp"
#include "Intrepid2_Utils.hpp"

#include "Intrepid2_CellTopologyTags.hpp"
#include "Kokkos_Vector.hpp"
#include "Shards_CellTopology.hpp"

#include <vector>

namespace Intrepid2 {

  template<typename ExecSpaceType = void,
           typename outputValueType = double,
           typename pointValueType = double>
  class Basis {
  public:
    using ExecutionSpace  = ExecSpaceType;
    
    using OutputValueType = outputValueType;
    
    using PointValueType  = pointValueType;
    
    using OrdinalViewType = Kokkos::View<ordinal_type,ExecSpaceType>;

    using EBasisViewType = Kokkos::View<EBasis,ExecSpaceType>;

    using ECoordinatesViewType = Kokkos::View<ECoordinates,ExecSpaceType>;

    // ** tag interface
    //  - tag interface is not decorated with Kokkos inline so it should be allocated on hostspace

    using OrdinalTypeArray1DHost = Kokkos::View<ordinal_type*,typename ExecSpaceType::array_layout,Kokkos::HostSpace>;

    using OrdinalTypeArray2DHost = Kokkos::View<ordinal_type**,typename ExecSpaceType::array_layout,Kokkos::HostSpace>;

    using OrdinalTypeArray3DHost = Kokkos::View<ordinal_type***,typename ExecSpaceType::array_layout,Kokkos::HostSpace>;

    using OrdinalTypeArrayStride1DHost = Kokkos::View<ordinal_type*, Kokkos::LayoutStride, Kokkos::HostSpace>;

    using OrdinalTypeArray1D = Kokkos::View<ordinal_type*,ExecSpaceType>;

    using OrdinalTypeArray2D = Kokkos::View<ordinal_type**,ExecSpaceType>;

    using OrdinalTypeArray3D = Kokkos::View<ordinal_type***,ExecSpaceType>;

    using OrdinalTypeArrayStride1D = Kokkos::View<ordinal_type*, Kokkos::LayoutStride, ExecSpaceType>;

    typedef typename ScalarTraits<pointValueType>::scalar_type scalarType;
  protected:

    ordinal_type basisCardinality_;

    ordinal_type basisDegree_;

    shards::CellTopology basisCellTopology_;

    EBasis basisType_;

    ECoordinates basisCoordinates_;
    
    EFunctionSpace functionSpace_ = FUNCTION_SPACE_MAX;
    
    //Kokkos::View<bool,ExecSpaceType> basisTagsAreSet_;

    OrdinalTypeArray2DHost ordinalToTag_;

    OrdinalTypeArray3DHost tagToOrdinal_;

    template<typename OrdinalTypeView3D,
             typename OrdinalTypeView2D,
             typename OrdinalTypeView1D>
    void setOrdinalTagData(       OrdinalTypeView3D &tagToOrdinal,
                                  OrdinalTypeView2D &ordinalToTag,
                            const OrdinalTypeView1D  tags,
                            const ordinal_type       basisCard,
                            const ordinal_type       tagSize,
                            const ordinal_type       posScDim,
                            const ordinal_type       posScOrd,
                            const ordinal_type       posDfOrd ) {
      // Create ordinalToTag
      ordinalToTag = OrdinalTypeView2D("ordinalToTag", basisCard, tagSize);

      // Initialize with -1
      Kokkos::deep_copy( ordinalToTag, -1 );

      // Copy tags
      for (ordinal_type i=0;i<basisCard;++i)
        for (ordinal_type j=0;j<tagSize;++j)
          ordinalToTag(i, j) = tags(i*tagSize + j);

      // Find out dimension of tagToOrdinal
      auto maxScDim = 0;  // first dimension of tagToOrdinal
      for (ordinal_type i=0;i<basisCard;++i)
        if (maxScDim < tags(i*tagSize + posScDim))
          maxScDim = tags(i*tagSize + posScDim);
      ++maxScDim;

      auto maxScOrd = 0; // second dimension of tagToOrdinal
      for (ordinal_type i=0;i<basisCard;++i)
        if (maxScOrd < tags(i*tagSize + posScOrd))
          maxScOrd = tags(i*tagSize + posScOrd);
      ++maxScOrd;

      auto maxDfOrd = 0;  // third dimension of tagToOrdinal
      for (ordinal_type i=0;i<basisCard;++i)
        if (maxDfOrd < tags(i*tagSize + posDfOrd))
          maxDfOrd = tags(i*tagSize + posDfOrd);
      ++maxDfOrd;

      // Create tagToOrdinal
      tagToOrdinal = OrdinalTypeView3D("tagToOrdinal", maxScDim, maxScOrd, maxDfOrd);

      // Initialize with -1
      Kokkos::deep_copy( tagToOrdinal, -1 );

      // Overwrite elements of the array corresponding to tags with local DoF Id's, leave all other = -1
      for (ordinal_type i=0;i<basisCard;++i)
        tagToOrdinal(tags(i*tagSize), tags(i*tagSize+1), tags(i*tagSize+2)) = i;
    }

    // dof coords
    Kokkos::DynRankView<scalarType,ExecSpaceType> dofCoords_;

    // dof coeffs
    Kokkos::DynRankView<scalarType,ExecSpaceType> dofCoeffs_;
    
    OrdinalTypeArray2DHost fieldOrdinalPolynomialDegree_;
  public:

    Basis() = default;
    virtual~Basis() = default;

    // receives input arguments
    OutputValueType getDummyOutputValue() { return outputValueType(); }

    PointValueType getDummyPointValue() { return pointValueType(); }

    using OutputViewType = Kokkos::DynRankView<OutputValueType,Kokkos::LayoutStride,ExecSpaceType>;

    using PointViewType = Kokkos::DynRankView<PointValueType,Kokkos::LayoutStride,ExecSpaceType>;

    using ScalarViewType = Kokkos::DynRankView<scalarType,Kokkos::LayoutStride,ExecSpaceType>;
    
    virtual
    void
    getValues(       OutputViewType /* outputValues */,
               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 overridden accordingly by derived classes.");
    }

    virtual
    void
    getValues(        OutputViewType /* outputValues */,
                const PointViewType  /* inputPoints */,
                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 overridden accordingly by derived classes.");
    }


    virtual
    void
    getDofCoords( ScalarViewType /* dofCoords */ ) const {
      INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
                                    ">>> ERROR (Basis::getDofCoords): this method is not supported or should be overridden accordingly by derived classes.");
    }

    virtual
    void
    getDofCoeffs( ScalarViewType /* dofCoeffs */ ) const {
      INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
                                    ">>> ERROR (Basis::getDofCoeffs): this method is not supported or should be overridden accordingly by derived classes.");
    }

    OrdinalTypeArray1DHost getFieldOrdinalsForDegree(OrdinalTypeArray1DHost &degrees) const
    {
      INTREPID2_TEST_FOR_EXCEPTION( basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
                                   ">>> ERROR (Basis::getFieldOrdinalsForDegree): this method is not supported for non-hierarchical bases.");
      int degreeEntryLength     = fieldOrdinalPolynomialDegree_.extent_int(1);
      int requestedDegreeLength = degrees.extent_int(0);
      INTREPID2_TEST_FOR_EXCEPTION(degreeEntryLength != requestedDegreeLength, std::invalid_argument, "length of degrees does not match the entries in fieldOrdinalPolynomialDegree_");
      std::vector<int> fieldOrdinalsVector;
      for (int basisOrdinal=0; basisOrdinal<fieldOrdinalPolynomialDegree_.extent_int(0); basisOrdinal++)
      {
        bool matches = true;
        for (int d=0; d<degreeEntryLength; d++)
        {
          if (fieldOrdinalPolynomialDegree_(basisOrdinal,d) > degrees(d)) matches = false;
        }
        if (matches) fieldOrdinalsVector.push_back(basisOrdinal);
      }
      OrdinalTypeArray1DHost fieldOrdinals("fieldOrdinalsForDegree",fieldOrdinalsVector.size());
      for (unsigned i=0; i<fieldOrdinalsVector.size(); i++)
      {
        fieldOrdinals(i) = fieldOrdinalsVector[i];
      }
      return fieldOrdinals;
    }
    
    std::vector<int> getFieldOrdinalsForDegree(std::vector<int> &degrees) const
    {
      INTREPID2_TEST_FOR_EXCEPTION( basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
                                   ">>> ERROR (Basis::getFieldOrdinalsForDegree): this method is not supported for non-hierarchical bases.");
      OrdinalTypeArray1DHost degreesView("degrees",degrees.size());
      for (unsigned d=0; d<degrees.size(); d++)
      {
        degreesView(d) = degrees[d];
      }
      auto fieldOrdinalsView = getFieldOrdinalsForDegree(degreesView);
      std::vector<int> fieldOrdinalsVector(fieldOrdinalsView.extent_int(0));
      for (int i=0; i<fieldOrdinalsView.extent_int(0); i++)
      {
        fieldOrdinalsVector[i] = fieldOrdinalsView(i);
      }
      return fieldOrdinalsVector;
    }
    
    OrdinalTypeArray1DHost getPolynomialDegreeOfField(int fieldOrdinal) const
    {
      INTREPID2_TEST_FOR_EXCEPTION( basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
                                   ">>> ERROR (Basis::getPolynomialDegreeOfField): this method is not supported for non-hierarchical bases.");
      INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinal < 0, std::invalid_argument, "field ordinal must be non-negative");
      INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinal >= fieldOrdinalPolynomialDegree_.extent_int(0), std::invalid_argument, "field ordinal out of bounds");
      
      int polyDegreeLength = getPolynomialDegreeLength(); 
      OrdinalTypeArray1DHost polyDegree("polynomial degree", polyDegreeLength);
      for (int d=0; d<polyDegreeLength; d++)
      {
        polyDegree(d) = fieldOrdinalPolynomialDegree_(fieldOrdinal,d);
      }
      return polyDegree;
    }
    
    std::vector<int> getPolynomialDegreeOfFieldAsVector(int fieldOrdinal) const
    {
      INTREPID2_TEST_FOR_EXCEPTION( basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
                                   ">>> ERROR (Basis::getPolynomialDegreeOfFieldAsVector): this method is not supported for non-hierarchical bases.");
      auto polynomialDegreeView = getPolynomialDegreeOfField(fieldOrdinal);
      std::vector<int> polynomialDegree(polynomialDegreeView.extent_int(0));
      
      for (unsigned d=0; d<polynomialDegree.size(); d++)
      {
        polynomialDegree[d] = polynomialDegreeView(d);
      }
      return polynomialDegree;
    }
    
    int getPolynomialDegreeLength() const
    {
      INTREPID2_TEST_FOR_EXCEPTION( basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
                                   ">>> ERROR (Basis::getPolynomialDegreeLength): this method is not supported for non-hierarchical bases.");
      return fieldOrdinalPolynomialDegree_.extent_int(1);
    }
    
    virtual
    const char*
    getName() const {
      return "Intrepid2_Basis";
    }

    virtual
    bool
    requireOrientation() const {
      return false;
    }

    ordinal_type
    getCardinality() const {
      return basisCardinality_;
    }


    ordinal_type
    getDegree() const {
      return basisDegree_;
    }
    
    EFunctionSpace
    getFunctionSpace() const {
      return functionSpace_;
    }
    
    shards::CellTopology
    getBaseCellTopology() const {
      return basisCellTopology_;
    }


    EBasis
    getBasisType() const {
      return basisType_;
    }


    ECoordinates
    getCoordinateSystem() const {
      return basisCoordinates_;
    }

    ordinal_type
    getDofCount( const ordinal_type subcDim,
                 const ordinal_type subcOrd ) const {
      if ( subcDim >= 0   &&   subcDim < static_cast<ordinal_type>(tagToOrdinal_.extent(0)) &&
           subcOrd >= 0   &&   subcOrd < static_cast<ordinal_type>(tagToOrdinal_.extent(1)) )
      {
        int firstDofOrdinal = tagToOrdinal_(subcDim, subcOrd, 0); // will be -1 if no dofs for subcell
        if (firstDofOrdinal == -1) return static_cast<ordinal_type>(0);
        // otherwise, lookup its tag and return the dof count stored there
        return static_cast<ordinal_type>(this->getDofTag(firstDofOrdinal)[3]);
      }
      else
      {
        // subcDim and/or subcOrd out of bounds -> no dofs associated with subcell
        return static_cast<ordinal_type>(0);
      }
    }

    ordinal_type
    getDofOrdinal( const ordinal_type subcDim,
                   const ordinal_type subcOrd,
                   const ordinal_type subcDofOrd ) const {
      // this should be abort and able to be called as a device function
#ifdef HAVE_INTREPID2_DEBUG
      INTREPID2_TEST_FOR_EXCEPTION( subcDim < 0 || subcDim >= static_cast<ordinal_type>(tagToOrdinal_.extent(0)), std::out_of_range,
                                    ">>> ERROR (Basis::getDofOrdinal): subcDim is out of range");
      INTREPID2_TEST_FOR_EXCEPTION( subcOrd < 0 || subcOrd >= static_cast<ordinal_type>(tagToOrdinal_.extent(1)), std::out_of_range,
                                    ">>> ERROR (Basis::getDofOrdinal): subcOrd is out of range");
      INTREPID2_TEST_FOR_EXCEPTION( subcDofOrd < 0 || subcDofOrd >= static_cast<ordinal_type>(tagToOrdinal_.extent(2)), std::out_of_range,
                                    ">>> ERROR (Basis::getDofOrdinal): subcDofOrd is out of range");
#endif
      ordinal_type r_val = -1;
      if ( subcDim    < static_cast<ordinal_type>(tagToOrdinal_.extent(0)) &&
           subcOrd    < static_cast<ordinal_type>(tagToOrdinal_.extent(1)) &&
           subcDofOrd < static_cast<ordinal_type>(tagToOrdinal_.extent(2)) )
        r_val = tagToOrdinal_(subcDim, subcOrd, subcDofOrd);
#ifdef HAVE_INTREPID2_DEBUG
      INTREPID2_TEST_FOR_EXCEPTION( r_val == -1, std::runtime_error,
                                    ">>> ERROR (Basis::getDofOrdinal): Invalid DoF tag is found.");
#endif
      return r_val;
    }

    const OrdinalTypeArray3DHost
    getAllDofOrdinal() const {
      return tagToOrdinal_;
    }


    const OrdinalTypeArrayStride1DHost
    getDofTag( const ordinal_type dofOrd ) const {
#ifdef HAVE_INTREPID2_DEBUG
      INTREPID2_TEST_FOR_EXCEPTION( dofOrd < 0 || dofOrd >= static_cast<ordinal_type>(ordinalToTag_.extent(0)), std::out_of_range,
                                    ">>> ERROR (Basis::getDofTag): dofOrd is out of range");
#endif
      return Kokkos::subview(ordinalToTag_, dofOrd, Kokkos::ALL());
    }


    const OrdinalTypeArray2DHost
    getAllDofTags() const {
      return ordinalToTag_;
    }

  }; // class Basis


  //--------------------------------------------------------------------------------------------//
  //                                                                                            //
  //                            Helper functions of the Basis class                             //
  //                                                                                            //
  //--------------------------------------------------------------------------------------------//

  //
  // functions for orders, cardinality, etc.
  //


  KOKKOS_INLINE_FUNCTION
  ordinal_type getFieldRank(const EFunctionSpace spaceType);

  KOKKOS_INLINE_FUNCTION
  ordinal_type getOperatorRank(const EFunctionSpace spaceType,
                               const EOperator      operatorType,
                               const ordinal_type   spaceDim);

  KOKKOS_INLINE_FUNCTION
  ordinal_type getOperatorOrder(const EOperator operatorType);

  template<EOperator operatorType>
  KOKKOS_INLINE_FUNCTION
  constexpr ordinal_type getOperatorOrder();

  template<ordinal_type spaceDim>
  KOKKOS_INLINE_FUNCTION
  ordinal_type getDkEnumeration(const ordinal_type xMult,
                                const ordinal_type yMult = -1,
                                const ordinal_type zMult = -1);  


  template<ordinal_type spaceDim>
  KOKKOS_INLINE_FUNCTION
  ordinal_type getPnEnumeration(const ordinal_type p,
                                const ordinal_type q = 0,
                                const ordinal_type r = 0);



template<typename value_type>
KOKKOS_INLINE_FUNCTION
void getJacobyRecurrenceCoeffs (
          value_type  &an,
          value_type  &bn,
          value_type  &cn,
    const ordinal_type alpha,
    const ordinal_type beta ,
    const ordinal_type n);





  // /** \brief  Returns multiplicities of dx, dy, and dz based on the enumeration of the partial
  //     derivative, its order and the space dimension. Inverse of the getDkEnumeration() method.

  //     \param  partialMult      [out]    - array with the multiplicities f dx, dy and dz
  //     \param  derivativeEnum   [in]     - enumeration of the partial derivative
  //     \param  operatorType     [in]     - k-th partial derivative Dk
  //     \param  spaceDim         [in]     - space dimension
  // */
  // template<typename OrdinalArrayType>
  // KOKKOS_INLINE_FUNCTION
  // void getDkMultiplicities(OrdinalArrayType   partialMult,
  //                          const ordinal_type derivativeEnum,
  //                          const EOperator    operatorType,
  //                          const ordinal_type spaceDim);

  KOKKOS_INLINE_FUNCTION
  ordinal_type getDkCardinality(const EOperator    operatorType,
                                const ordinal_type spaceDim);

  template<EOperator operatorType, ordinal_type spaceDim>
  KOKKOS_INLINE_FUNCTION
  constexpr ordinal_type getDkCardinality();



   template<ordinal_type spaceDim>
   KOKKOS_INLINE_FUNCTION
   ordinal_type getPnCardinality (ordinal_type n);

   template<ordinal_type spaceDim, ordinal_type n>
   KOKKOS_INLINE_FUNCTION
   constexpr ordinal_type getPnCardinality ();



  //
  // Argument check
  //


  template<typename outputValueViewType,
           typename inputPointViewType>
  void getValues_HGRAD_Args( const outputValueViewType   outputValues,
                             const inputPointViewType    inputPoints,
                             const EOperator             operatorType,
                             const shards::CellTopology  cellTopo,
                             const ordinal_type          basisCard );

  template<typename outputValueViewType,
           typename inputPointViewType>
  void getValues_HCURL_Args( const outputValueViewType   outputValues,
                             const inputPointViewType    inputPoints,
                             const EOperator             operatorType,
                             const shards::CellTopology  cellTopo,
                             const ordinal_type          basisCard );

  template<typename outputValueViewType,
           typename inputPointViewType>
  void getValues_HDIV_Args( const  outputValueViewType  outputValues,
                            const inputPointViewType    inputPoints,
                            const EOperator             operatorType,
                            const shards::CellTopology  cellTopo,
                            const ordinal_type          basisCard );

  template<typename outputValueViewType,
           typename inputPointViewType>
  void getValues_HVOL_Args( const outputValueViewType   outputValues,
                             const inputPointViewType    inputPoints,
                             const EOperator             operatorType,
                             const shards::CellTopology  cellTopo,
                             const ordinal_type          basisCard );

}// namespace Intrepid2

#include <Intrepid2_BasisDef.hpp>

//--------------------------------------------------------------------------------------------//
//                                                                                            //
//                            D O C U M E N T A T I O N   P A G E S                           //
//                                                                                            //
//--------------------------------------------------------------------------------------------//
#endif