Intrepid2_Cubature.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//                           Intrepid2 Package
//
// Copyright 2007 NTESS and the Intrepid2 contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef __INTREPID2_CUBATURE_HPP__
#define __INTREPID2_CUBATURE_HPP__

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

namespace Intrepid2 {

  /* \struct Intrepid2::CubatureTemplate
      \brief  Template for the cubature rules used by Intrepid. Cubature template consists of
              cubature points and cubature weights. Intrepid provides a collection of cubature
              templates for most standard cell topologies. The templates are defined in reference
              coordinates using a standard reference cell for each canonical cell type. Cubature
              points are always specified by a triple of (X,Y,Z) coordinates even if the cell
              dimension is less than 3. The unused dimensions should be padded by zeroes.

              For example, a set of Gauss rules on [-1,1] looks as the following array of CubatureTemplate structs:

  \verbatim
  cubature_rule[4] =
  {                                         // Collection of Gauss rules on [-1,1]
    {
      1,                                      ----> number of points in the rule
      {{0.0,0.0,0.0}},                        ----> X,Y,Z coordinates of the cubature points
      {0.5}                                   ----> the cubature weight
    },
    {
      2,
      {{-sqrt(1.0/3.0),0.0,0.0},
       {+sqrt(1.0/3.0),0.0,0.0}},
      {1.0,1.0}
    },
    {
      3,
      {{-sqrt(3.0/5.0),0.0,0.0},
       {0.0,0.0,0.0},
       {+sqrt(3.0/5.0),0.0,0.0}},
      {5.0/9.0, 8.0/9.0,5.0/9.0}
    },
    {
      4,
      {{-sqrt((3.0+4.0*sqrt(0.3))/7.0),0.0,0.0},
       {-sqrt((3.0-4.0*sqrt(0.3))/7.0),0.0,0.0},
       {+sqrt((3.0-4.0*sqrt(0.3))/7.0),0.0,0.0},
       {+sqrt((3.0+4.0*sqrt(0.3))/7.0),0.0,0.0}},
      //
      {0.5-sqrt(10.0/3.0)/12.0,
       0.5+sqrt(10.0/3.0)/12.0,
       0.5+sqrt(10.0/3.0)/12.0,
       0.5-sqrt(10.0/3.0)/12.0}
    }
  };  // end Gauss
  \endverbatim

      Also see data member documentation.
  */


  template<typename DeviceType = void,
           typename pointValueType = double,
           typename weightValueType = double>
  class Cubature {
  public:
    using ExecSpaceType = typename DeviceType::execution_space;
    using PointViewType             = Kokkos::DynRankView<pointValueType,Kokkos::LayoutStride,DeviceType>;
    using weightViewType            = Kokkos::DynRankView<weightValueType,Kokkos::LayoutStride,DeviceType>;

    using PointViewTypeAllocatable  = Kokkos::DynRankView<pointValueType,DeviceType>;  // uses default layout; allows us to allocate (in contrast to LayoutStride)
    using WeightViewTypeAllocatable = Kokkos::DynRankView<weightValueType,DeviceType>; // uses default layout; allows us to allocate (in contrast to LayoutStride)
    using TensorPointDataType       = TensorPoints<pointValueType,DeviceType>;
    using TensorWeightDataType      = TensorData<weightValueType,DeviceType>;
    
    virtual TensorPointDataType allocateCubaturePoints() const
    {
      // default implementation has trivial tensor structure
      PointViewTypeAllocatable allPointData("cubature points",this->getNumPoints(),this->getDimension());
      Kokkos::Array< PointViewTypeAllocatable, 1> cubaturePointComponents;
      cubaturePointComponents[0] = allPointData;
      return TensorPointDataType(cubaturePointComponents);
    }
    
    virtual TensorWeightDataType allocateCubatureWeights() const
    {
      // default implementation has trivial tensor structure
      using WeightDataType = Data<weightValueType,DeviceType>;
      WeightViewTypeAllocatable allWeightData("cubature weights",this->getNumPoints());
      Kokkos::Array< WeightDataType, 1> cubatureWeightComponents;
      cubatureWeightComponents[0] = WeightDataType(allWeightData);
      return TensorWeightDataType(cubatureWeightComponents);
    }
    
    virtual
    void
    getCubature( PointViewType  /* cubPoints */,
                 weightViewType /* cubWeights */ ) const {
      INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
                                    ">>> ERROR (Cubature::getCubature): this method should be overridden by derived classes.");
    }

    virtual
    void
    getCubature( PointViewType  /* cubPoints */,
                 weightViewType /* cubWeights */,
                 PointViewType  /* cellVertices */) const {
      INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
                                    ">>> ERROR (Cubature::getCubature): this method should be overridden by derived classes.");
    }
    
    virtual
    void
    getCubature( const TensorPointDataType & tensorCubPoints,
                       const TensorWeightDataType & tensorCubWeights) const {
      // default implementation has trivial tensor product structure.
      
      INTREPID2_TEST_FOR_EXCEPTION(1 != tensorCubPoints.numTensorComponents(), std::logic_error, "default implementation of getCubature only supports trivial tensor structure -- numTensorComponents() must be 1");
      INTREPID2_TEST_FOR_EXCEPTION(1 != tensorCubWeights.numTensorComponents(), std::logic_error, "default implementation of getCubature only supports trivial tensor structure -- numTensorComponents() must be 1");
      
      auto underlyingPointView  = tensorCubPoints.getTensorComponent(0);
      
      this->getCubature(underlyingPointView, tensorCubWeights.getTensorComponent(0).getUnderlyingView());
    }

    virtual
    ordinal_type
    getNumPoints() const {
      INTREPID2_TEST_FOR_WARNING( true, 
                                  ">>> ERROR (Cubature::getNumPoints): this method should be overridden by derived classes.");
      return 0;
    }


    virtual
    ordinal_type
    getDimension() const {
      INTREPID2_TEST_FOR_WARNING( true, 
                                  ">>> ERROR (Cubature::getDimension): this method should be overridden by derived classes.");
      return 0;
    }

    virtual
    ordinal_type
    getAccuracy() const {
      INTREPID2_TEST_FOR_WARNING( true, 
                                  ">>> ERROR (Cubature::getDimension): this method should be overridden by derived classes.");
      return 0;
    }

    virtual
    const char*
    getName() const  {
      return "Cubature";
    }

    Cubature() = default;
    virtual ~Cubature() {}

  };

}// end namespace Intrepid2


#endif