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), or
//                    Mauro Perego  (mperego@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 "Shards_CellTopology.hpp"

namespace Intrepid2 {

  template<typename ExecSpaceType = void,
           typename outputValueType = double,
           typename pointValueType = double>
  class Basis {
  public:

    typedef Kokkos::View<ordinal_type,ExecSpaceType> ordinal_view_type;

    typedef Kokkos::View<EBasis,ExecSpaceType> ebasis_view_type;

    typedef Kokkos::View<ECoordinates,ExecSpaceType> ecoordinates_view_type;

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

    typedef Kokkos::View<ordinal_type*  ,typename ExecSpaceType::array_layout,Kokkos::HostSpace> ordinal_type_array_1d_host;

    typedef Kokkos::View<ordinal_type** ,typename ExecSpaceType::array_layout,Kokkos::HostSpace> ordinal_type_array_2d_host;

    typedef Kokkos::View<ordinal_type***,typename ExecSpaceType::array_layout,Kokkos::HostSpace> ordinal_type_array_3d_host;

    typedef Kokkos::View<ordinal_type*  , Kokkos::LayoutStride, Kokkos::HostSpace> ordinal_type_array_stride_1d_host;

    typedef Kokkos::View<ordinal_type*  ,ExecSpaceType> ordinal_type_array_1d;

    typedef Kokkos::View<ordinal_type** ,ExecSpaceType> ordinal_type_array_2d;

    typedef Kokkos::View<ordinal_type***,ExecSpaceType> ordinal_type_array_3d;

    typedef Kokkos::View<ordinal_type*  , Kokkos::LayoutStride, ExecSpaceType> ordinal_type_array_stride_1d;

    typedef typename ScalarTraits<pointValueType>::scalar_type scalarType;

  protected:

    ordinal_type basisCardinality_;

    ordinal_type basisDegree_;

    shards::CellTopology basisCellTopology_;

    EBasis basisType_;

    ECoordinates basisCoordinates_;

    //Kokkos::View<bool,ExecSpaceType> basisTagsAreSet_;

    ordinal_type_array_2d_host ordinalToTag_;

    ordinal_type_array_3d_host 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_;

  public:

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

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

    pointValueType getDummyPointValue() { return pointValueType(); }

    typedef Kokkos::DynRankView<outputValueType,Kokkos::LayoutStride,ExecSpaceType> outputViewType;

    typedef Kokkos::DynRankView<pointValueType,Kokkos::LayoutStride,ExecSpaceType>  pointViewType;

    typedef Kokkos::DynRankView<scalarType,Kokkos::LayoutStride,ExecSpaceType>      scalarViewType;

    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 over-riden 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 over-riden 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 over-riden 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 over-riden accordingly by derived classes.");
    }


    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_;
    }


    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 ordinal_type_array_3d_host
    getAllDofOrdinal() const {
      return tagToOrdinal_;
    }


    const ordinal_type_array_stride_1d_host
    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 ordinal_type_array_2d_host
    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