docs/intrepid2/Intrepid2__HierarchicalBasis__HDIV__TET_8hpp_source.html

// @HEADER

// *****************************************************************************

//                           Intrepid2 Package

//

// Copyright 2007 NTESS and the Intrepid2 contributors.

// SPDX-License-Identifier: BSD-3-Clause

// *****************************************************************************

// @HEADER


#ifndef Intrepid2_HierarchicalBasis_HDIV_TET_h

#define Intrepid2_HierarchicalBasis_HDIV_TET_h


#include <Kokkos_DynRankView.hpp>


#include <Intrepid2_config.h>


#include "Intrepid2_Basis.hpp"

#include "Intrepid2_LegendreBasis_HVOL_TRI.hpp"

#include "Intrepid2_Polynomials.hpp"

#include "Intrepid2_Utils.hpp"


namespace Intrepid2

{

  template<class DeviceType, class OutputScalar, class PointScalar,

           class OutputFieldType, class InputPointsType>


  struct Hierarchical_HDIV_TET_Functor

  {

    using ExecutionSpace     = typename DeviceType::execution_space;

    using ScratchSpace       = typename ExecutionSpace::scratch_memory_space;

    using OutputScratchView  = Kokkos::View<OutputScalar*,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;

    using PointScratchView   = Kokkos::View<PointScalar*, ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;


    using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;

    using TeamMember = typename TeamPolicy::member_type;


    EOperator opType_;


    OutputFieldType  output_;      // F,P

    InputPointsType  inputPoints_; // P,D


    ordinal_type polyOrder_;

    ordinal_type numFields_, numPoints_;


    size_t fad_size_output_;


    static constexpr ordinal_type numVertices = 4;

    static constexpr ordinal_type numFaces = 4;

    static constexpr ordinal_type numVerticesPerFace = 3;

    static constexpr ordinal_type numInteriorFamilies = 3;


    // index into face_vertices with faceOrdinal * numVerticesPerFace + vertexNumber

    const ordinal_type face_vertices[numFaces*numVerticesPerFace] = {0,1,2, // face 0

                                                                     0,1,3, // face 1

                                                                     1,2,3, // face 2

                                                                     0,2,3  // face 3

                                                                    };


    const ordinal_type numFaceFunctionsPerFace_;

    const ordinal_type numFaceFunctions_;

    const ordinal_type numInteriorFunctionsPerFamily_;

    const ordinal_type numInteriorFunctions_;


    // interior basis functions are computed in terms of certain face basis functions.

    const ordinal_type faceOrdinalForInterior_[numInteriorFamilies] = {0,2,3};

    const ordinal_type faceFamilyForInterior_[numInteriorFamilies]  = {0,0,1};

    const ordinal_type interiorCoordinateOrdinal_[numInteriorFamilies] = {3,0,1}; // m, where V^b_{ijk} is computed in terms of [L^{2(i+j)}_k](1-lambda_m, lambda_m)


    //

    const ordinal_type interior_face_family_start_ [numInteriorFamilies] = {0,0,1};

    const ordinal_type interior_face_family_middle_[numInteriorFamilies] = {1,1,2};

    const ordinal_type interior_face_family_end_   [numInteriorFamilies] = {2,2,0};


    KOKKOS_INLINE_FUNCTION

    ordinal_type dofOrdinalForFace(const ordinal_type &faceOrdinal,

                                   const ordinal_type &zeroBasedFaceFamily,

                                   const ordinal_type &i,

                                   const ordinal_type &j) const

    {

      // determine where the functions for this face start

      const ordinal_type faceDofOffset = faceOrdinal * numFaceFunctionsPerFace_;


      // rather than deriving a closed formula in terms of i and j (which is potentially error-prone),

      // we simply step through a for loop much as we do in the basis computations themselves.  (This method

      // is not expected to be called so much as to be worth optimizing.)


      const ordinal_type max_ij_sum = polyOrder_ - 1;


      ordinal_type fieldOrdinal = faceDofOffset + zeroBasedFaceFamily; // families are interleaved on the face.


      for (ordinal_type ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)

      {

        for (ordinal_type ii=0; ii<ij_sum; ii++)

        {

          // j will be ij_sum - i; j >= 1.

          const ordinal_type jj = ij_sum - ii; // jj >= 1

          if ( (ii == i) && (jj == j))

          {

            // have reached the (i,j) we're looking for

            return fieldOrdinal;

          }

          fieldOrdinal++;

        }

      }

      return -1; // error: not found.

    }


    Hierarchical_HDIV_TET_Functor(EOperator opType, OutputFieldType output, InputPointsType inputPoints, int polyOrder)

    : opType_(opType), output_(output), inputPoints_(inputPoints),

      polyOrder_(polyOrder),

      fad_size_output_(getScalarDimensionForView(output)),

      numFaceFunctionsPerFace_(polyOrder * (polyOrder+1)/2), // p*(p+1)/2 functions per face

      numFaceFunctions_(numFaceFunctionsPerFace_*numFaces),  // 4 faces

      numInteriorFunctionsPerFamily_((polyOrder-1)*polyOrder*(polyOrder+1)/6),    // (p+1) choose 3

      numInteriorFunctions_(numInteriorFunctionsPerFamily_ * numInteriorFamilies) // 3 families of interior functions

    {

      numFields_ = output.extent_int(0);

      numPoints_ = output.extent_int(1);


      const ordinal_type expectedCardinality  = numFaceFunctions_ + numInteriorFunctions_;


      // interior family I: computed in terms of face 012 (face ordinal 0), ordinal 0 in face family I.

      // interior family II: computed in terms of face 123 (face ordinal 2), ordinal 2 in face family I.

      // interior family III: computed in terms of face 230 (face ordinal 3), ordinal 3 in face family II.


      INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument, "point counts need to match!");

      INTREPID2_TEST_FOR_EXCEPTION(numFields_ != expectedCardinality, std::invalid_argument, "output field size does not match basis cardinality");

    }


    KOKKOS_INLINE_FUNCTION

    void computeFaceJacobi(OutputScratchView &P_2ip1,

                           const ordinal_type &zeroBasedFaceOrdinal,

                           const ordinal_type &i,

                           const PointScalar* lambda) const

    {

      // index into face_vertices with faceOrdinal * numVerticesPerFace + vertexNumber

      const auto &s0_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 0];

      const auto &s1_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 1];

      const auto &s2_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 2];


      const auto & s0 = lambda[s0_index];

      const auto & s1 = lambda[s1_index];

      const auto & s2 = lambda[s2_index];

      const PointScalar jacobiScaling = s0 + s1 + s2;


      const double alpha = i*2.0 + 1;

      Polynomials::shiftedScaledJacobiValues(P_2ip1, alpha, polyOrder_-1, s2, jacobiScaling);

    }


    KOKKOS_INLINE_FUNCTION


    void computeFaceJacobiForInterior(OutputScratchView &P_2ip1,

                                      const ordinal_type &zeroBasedInteriorFamilyOrdinal,

                                      const ordinal_type &i,

                                      const PointScalar* lambda) const

    {

      const ordinal_type & relatedFaceOrdinal = faceOrdinalForInterior_[zeroBasedInteriorFamilyOrdinal];

      const auto &s0_vertex_number = interior_face_family_start_ [zeroBasedInteriorFamilyOrdinal];

      const auto &s1_vertex_number = interior_face_family_middle_[zeroBasedInteriorFamilyOrdinal];

      const auto &s2_vertex_number = interior_face_family_end_   [zeroBasedInteriorFamilyOrdinal];


      // index into face_vertices with faceOrdinal * numVerticesPerFace + vertexNumber

      const auto &s0_index = face_vertices[relatedFaceOrdinal * numVerticesPerFace + s0_vertex_number];

      const auto &s1_index = face_vertices[relatedFaceOrdinal * numVerticesPerFace + s1_vertex_number];

      const auto &s2_index = face_vertices[relatedFaceOrdinal * numVerticesPerFace + s2_vertex_number];


      const auto & s0 = lambda[s0_index];

      const auto & s1 = lambda[s1_index];

      const auto & s2 = lambda[s2_index];

      const PointScalar jacobiScaling = s0 + s1 + s2;


      const double alpha = i*2.0 + 1;

      Polynomials::shiftedScaledJacobiValues(P_2ip1, alpha, polyOrder_-1, s2, jacobiScaling);

    }


    KOKKOS_INLINE_FUNCTION

    void computeFaceLegendre(OutputScratchView &P,

                             const ordinal_type &zeroBasedFaceOrdinal,

                             const PointScalar* lambda) const

    {

      // index into face_vertices with faceOrdinal * numVerticesPerFace + vertexNumber

      const auto &s0_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 0];

      const auto &s1_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 1];


      const auto & s0 = lambda[s0_index];

      const auto & s1 = lambda[s1_index];

      const PointScalar legendreScaling = s0 + s1;


      Polynomials::shiftedScaledLegendreValues(P, polyOrder_-1, s1, legendreScaling);

    }


    KOKKOS_INLINE_FUNCTION

    void computeFaceLegendreForInterior(OutputScratchView &P,

                                        const ordinal_type &zeroBasedInteriorFamilyOrdinal,

                                        const PointScalar* lambda) const

    {

      const ordinal_type & relatedFaceOrdinal = faceOrdinalForInterior_[zeroBasedInteriorFamilyOrdinal];

      const auto &s0_vertex_number = interior_face_family_start_ [zeroBasedInteriorFamilyOrdinal];

      const auto &s1_vertex_number = interior_face_family_middle_[zeroBasedInteriorFamilyOrdinal];


      // index into face_vertices with faceOrdinal * numVerticesPerFace + vertexNumber

      const auto &s0_index = face_vertices[relatedFaceOrdinal * numVerticesPerFace + s0_vertex_number];

      const auto &s1_index = face_vertices[relatedFaceOrdinal * numVerticesPerFace + s1_vertex_number];


      const auto & s0 = lambda[s0_index];

      const auto & s1 = lambda[s1_index];

      const PointScalar legendreScaling = s0 + s1;


      Polynomials::shiftedScaledLegendreValues(P, polyOrder_-1, s1, legendreScaling);

    }


    KOKKOS_INLINE_FUNCTION

    void computeFaceVectorWeight(OutputScalar &vectorWeight_x,

                                 OutputScalar &vectorWeight_y,

                                 OutputScalar &vectorWeight_z,

                                 const ordinal_type &zeroBasedFaceOrdinal,

                                 const PointScalar* lambda,

                                 const PointScalar* lambda_dx,

                                 const PointScalar* lambda_dy,

                                 const PointScalar* lambda_dz) const

    {

      // compute s0 (grad s1 x grad s2) + s1 (grad s2 x grad s0) + s2 (grad s0 x grad s1)


      const auto &s0_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 0];

      const auto &s1_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 1];

      const auto &s2_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 2];


      const auto & s0    = lambda   [s0_index];

      const auto & s0_dx = lambda_dx[s0_index];

      const auto & s0_dy = lambda_dy[s0_index];

      const auto & s0_dz = lambda_dz[s0_index];


      const auto & s1    = lambda   [s1_index];

      const auto & s1_dx = lambda_dx[s1_index];

      const auto & s1_dy = lambda_dy[s1_index];

      const auto & s1_dz = lambda_dz[s1_index];


      const auto & s2    = lambda   [s2_index];

      const auto & s2_dx = lambda_dx[s2_index];

      const auto & s2_dy = lambda_dy[s2_index];

      const auto & s2_dz = lambda_dz[s2_index];


      vectorWeight_x = s0 * (s1_dy * s2_dz - s1_dz * s2_dy)

                     + s1 * (s2_dy * s0_dz - s2_dz * s0_dy)

                     + s2 * (s0_dy * s1_dz - s0_dz * s1_dy);


      vectorWeight_y = s0 * (s1_dz * s2_dx - s1_dx * s2_dz)

                     + s1 * (s2_dz * s0_dx - s2_dx * s0_dz)

                     + s2 * (s0_dz * s1_dx - s0_dx * s1_dz);


      vectorWeight_z = s0 * (s1_dx * s2_dy - s1_dy * s2_dx)

                     + s1 * (s2_dx * s0_dy - s2_dy * s0_dx)

                     + s2 * (s0_dx * s1_dy - s0_dy * s1_dx);

    }


    // This is the "Ancillary Operator" V^{tri}_{ij} on p. 433 of Fuentes et al.

    KOKKOS_INLINE_FUNCTION

    void faceFunctionValue(OutputScalar &value_x,

                           OutputScalar &value_y,

                           OutputScalar &value_z,

                           const ordinal_type &i, // i >= 0

                           const ordinal_type &j, // j >= 0

                           const OutputScratchView &P,      // container in which shiftedScaledLegendreValues have been computed for the appropriate face

                           const OutputScratchView &P_2ip1, // container in which shiftedScaledJacobiValues have been computed for (2i+1) for the appropriate face

                           const OutputScalar &vectorWeight_x, // x component of s0 (grad s1 x grad s2) + s1 (grad s2 x grad s0) + s2 (grad s0 x grad s1)

                           const OutputScalar &vectorWeight_y, // y component

                           const OutputScalar &vectorWeight_z, // z component

                           const PointScalar* lambda) const

    {

      const auto &P_i      = P(i);

      const auto &P_2ip1_j = P_2ip1(j);


      value_x = P_i * P_2ip1_j * vectorWeight_x;

      value_y = P_i * P_2ip1_j * vectorWeight_y;

      value_z = P_i * P_2ip1_j * vectorWeight_z;

    }


    KOKKOS_INLINE_FUNCTION

    void computeFaceDivWeight(OutputScalar &divWeight,

                              const ordinal_type &zeroBasedFaceOrdinal,

                              const PointScalar* lambda_dx,

                              const PointScalar* lambda_dy,

                              const PointScalar* lambda_dz) const

    {

      // grad s0 \dot (grad s1 x grad s2)


      const auto &s0_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 0];

      const auto &s1_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 1];

      const auto &s2_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 2];


      const auto & s0_dx = lambda_dx[s0_index];

      const auto & s0_dy = lambda_dy[s0_index];

      const auto & s0_dz = lambda_dz[s0_index];


      const auto & s1_dx = lambda_dx[s1_index];

      const auto & s1_dy = lambda_dy[s1_index];

      const auto & s1_dz = lambda_dz[s1_index];


      const auto & s2_dx = lambda_dx[s2_index];

      const auto & s2_dy = lambda_dy[s2_index];

      const auto & s2_dz = lambda_dz[s2_index];


      divWeight = s0_dx * (s1_dy * s2_dz - s1_dz * s2_dy)

                + s0_dy * (s1_dz * s2_dx - s1_dx * s2_dz)

                + s0_dz * (s1_dx * s2_dy - s1_dy * s2_dx);

    }


    KOKKOS_INLINE_FUNCTION

    void computeInteriorIntegratedJacobi(OutputScratchView &L_2ipjp1,

                                         const ordinal_type &i,

                                         const ordinal_type &j,

                                         const ordinal_type &zeroBasedFamilyOrdinal,

                                         const PointScalar* lambda) const

    {

      const auto &lambda_m = lambda[interiorCoordinateOrdinal_[zeroBasedFamilyOrdinal]];


      const double alpha = 2 * (i + j + 1);


      const PointScalar jacobiScaling = 1.0;

      Polynomials::shiftedScaledIntegratedJacobiValues(L_2ipjp1, alpha, polyOrder_-1, lambda_m, jacobiScaling);

    }


    KOKKOS_INLINE_FUNCTION

    void computeInteriorJacobi(OutputScratchView &P_2ipjp1,

                               const ordinal_type &i,

                               const ordinal_type &j,

                               const ordinal_type &zeroBasedFamilyOrdinal,

                               const PointScalar* lambda) const

    {

      const auto &lambda_m = lambda[interiorCoordinateOrdinal_[zeroBasedFamilyOrdinal]];


      const double alpha = 2 * (i + j + 1);


      const PointScalar jacobiScaling = 1.0;

      Polynomials::shiftedScaledJacobiValues(P_2ipjp1, alpha, polyOrder_-1, lambda_m, jacobiScaling);

    }


    // Divergence of the "Ancillary Operator" V^{tri}_{ij} on p. 433 of Fuentes et al.

    KOKKOS_INLINE_FUNCTION

    void faceFunctionDiv(OutputScalar &divValue,

                         const ordinal_type &i, // i >= 0

                         const ordinal_type &j, // j >= 0

                         const OutputScratchView &P,      // container in which shiftedScaledLegendreValues have been computed for the appropriate face

                         const OutputScratchView &P_2ip1, // container in which shiftedScaledJacobiValues have been computed for (2i+1) for the appropriate face

                         const OutputScalar &divWeight,   // grad s0 \dot (grad s1 x grad s2)

                         const PointScalar* lambda) const

    {

      const auto &P_i      = P(i);

      const auto &P_2ip1_j = P_2ip1(j);


      divValue = (i + j + 3.) * P_i * P_2ip1_j * divWeight;

    }


    // grad ([L^{2(i+j+1)}_k](1-lambda_m,lambda_m)), used in divergence of interior basis functions

    KOKKOS_INLINE_FUNCTION

    void gradInteriorIntegratedJacobi(OutputScalar &L_2ipjp1_dx,

                                      OutputScalar &L_2ipjp1_dy,

                                      OutputScalar &L_2ipjp1_dz,

                                      const ordinal_type &zeroBasedFamilyOrdinal,

                                      const ordinal_type &j,

                                      const ordinal_type &k,

                                      const OutputScratchView &P_2ipjp1, // container in which shiftedScaledJacobiValues have been computed for alpha=2(i+j+1), t0=1-lambda_m, t1=lambda_m

                                      const PointScalar* lambda,

                                      const PointScalar* lambda_dx,

                                      const PointScalar* lambda_dy,

                                      const PointScalar* lambda_dz) const

    {

      // grad [L^alpha_k](t0,t1) = [P^alpha_{k-1}](t0,t1) grad(t1) + [R^alpha_{k-1}](t0,t1) grad(t1+t0)

      // here, t0 = 1-lambda_m, t1 = lambda_m ==> t1 + t0 = 1 ==> grad(t1+t0) = 0 ==> the R term vanishes.


      const ordinal_type &m = interiorCoordinateOrdinal_[zeroBasedFamilyOrdinal];


      L_2ipjp1_dx = P_2ipjp1(k-1) * lambda_dx[m];

      L_2ipjp1_dy = P_2ipjp1(k-1) * lambda_dy[m];

      L_2ipjp1_dz = P_2ipjp1(k-1) * lambda_dz[m];

    }


    KOKKOS_INLINE_FUNCTION

    void interiorFunctionDiv(OutputScalar &outputDiv,

                             OutputScalar &L_2ipjp1_k,

                             OutputScalar &faceDiv,

                             OutputScalar &L_2ipjp1_k_dx,

                             OutputScalar &L_2ipjp1_k_dy,

                             OutputScalar &L_2ipjp1_k_dz,

                             OutputScalar &faceValue_x,

                             OutputScalar &faceValue_y,

                             OutputScalar &faceValue_z) const

    {

      outputDiv = L_2ipjp1_k * faceDiv + L_2ipjp1_k_dx * faceValue_x + L_2ipjp1_k_dy * faceValue_y + L_2ipjp1_k_dz * faceValue_z;

    }


    KOKKOS_INLINE_FUNCTION

    void operator()(const TeamMember &teamMember) const {

      auto pointOrdinal = teamMember.league_rank();

      OutputScratchView scratch0, scratch1, scratch2, scratch3;

      if (fad_size_output_ > 0) {

        scratch0 = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);

        scratch1 = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);

        scratch2 = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);

        scratch3 = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);

      }

      else {

        scratch0 = OutputScratchView(teamMember.team_shmem(), polyOrder_);

        scratch1 = OutputScratchView(teamMember.team_shmem(), polyOrder_);

        scratch2 = OutputScratchView(teamMember.team_shmem(), polyOrder_);

        scratch3 = OutputScratchView(teamMember.team_shmem(), polyOrder_);

      }


      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[4]    = {1. - x - y - z, x, y, z};

      const PointScalar lambda_dx[4] = {-1., 1., 0., 0.};

      const PointScalar lambda_dy[4] = {-1., 0., 1., 0.};

      const PointScalar lambda_dz[4] = {-1., 0., 0., 1.};


      switch (opType_)

      {

        case OPERATOR_VALUE:

        {

          // face functions

          {

            // relabel scratch views

            auto &scratchP = scratch0;

            auto &scratchP_2ip1 = scratch1;


            const ordinal_type max_ij_sum = polyOrder_ - 1;


            for (ordinal_type faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)

            {

              OutputScalar divWeight;

              computeFaceDivWeight(divWeight, faceOrdinal, lambda_dx, lambda_dy, lambda_dz);


              OutputScalar vectorWeight_x, vectorWeight_y, vectorWeight_z;

              computeFaceVectorWeight(vectorWeight_x, vectorWeight_y, vectorWeight_z, faceOrdinal, lambda, lambda_dx, lambda_dy, lambda_dz);


              ordinal_type fieldOrdinal = faceOrdinal * numFaceFunctionsPerFace_;

              computeFaceLegendre(scratchP, faceOrdinal, lambda);


              for (int ij_sum=0; ij_sum <= max_ij_sum; ij_sum++)

              {

                for (int i=0; i<=ij_sum; i++)

                {

                  computeFaceJacobi(scratchP_2ip1, faceOrdinal, i, lambda);


                  const int j = ij_sum - i; // j >= 1


                  auto & output_x = output_(fieldOrdinal,pointOrdinal,0);

                  auto & output_y = output_(fieldOrdinal,pointOrdinal,1);

                  auto & output_z = output_(fieldOrdinal,pointOrdinal,2);


                  faceFunctionValue(output_x, output_y, output_z, i, j,

                                    scratchP, scratchP_2ip1, vectorWeight_x,

                                    vectorWeight_y, vectorWeight_z, lambda);


                  fieldOrdinal++;

                } // i

              } // ij_sum

            } // faceOrdinal

          } // face functions block


          // interior functions

          {

            // relabel scratch views

            auto &scratchP = scratch0;

            auto &scratchP_2ip1 = scratch1;

            auto &scratchL_2ipjp1 = scratch2; // L^{2(i+j+1)}, integrated Jacobi


            const ordinal_type min_ijk_sum = 1;

            const ordinal_type max_ijk_sum = polyOrder_-1;

            const ordinal_type min_ij_sum  = 0;

            const ordinal_type min_k       = 1;

            const ordinal_type min_j       = 0;

            const ordinal_type min_i       = 0;


            OutputScalar vectorWeight_x, vectorWeight_y, vectorWeight_z;


            for (int interiorFamilyOrdinal=1; interiorFamilyOrdinal<=numInteriorFamilies; interiorFamilyOrdinal++)

            {

              // following ESEAS, we interleave the interior families.  This groups all the interior dofs of a given degree together.


              ordinal_type interiorFamilyFieldOrdinal =

                  numFaceFunctions_ + interiorFamilyOrdinal - 1;


              const ordinal_type relatedFaceOrdinal = faceOrdinalForInterior_[interiorFamilyOrdinal-1];


              computeFaceLegendreForInterior(scratchP,

                                             interiorFamilyOrdinal - 1, lambda);

              computeFaceVectorWeight(vectorWeight_x, vectorWeight_y, vectorWeight_z, relatedFaceOrdinal, lambda, lambda_dx, lambda_dy, lambda_dz);


              for (int ijk_sum=min_ijk_sum; ijk_sum <= max_ijk_sum; ijk_sum++)

              {

                for (int ij_sum=min_ij_sum; ij_sum<=ijk_sum-min_k; ij_sum++)

                {

                  for (int i=min_i; i<=ij_sum-min_j; i++)

                  {

                    const ordinal_type j = ij_sum-i;

                    const ordinal_type k = ijk_sum - ij_sum;


                    computeFaceJacobiForInterior(

                        scratchP_2ip1, interiorFamilyOrdinal - 1, i, lambda);

                    computeInteriorIntegratedJacobi(scratchL_2ipjp1, i, j,

                                                    interiorFamilyOrdinal - 1,

                                                    lambda);


                    OutputScalar V_x, V_y, V_z;


                    faceFunctionValue(V_x, V_y, V_z, i, j, scratchP,

                                      scratchP_2ip1, vectorWeight_x,

                                      vectorWeight_y, vectorWeight_z, lambda);


                    auto &output_x =

                        output_(interiorFamilyFieldOrdinal, pointOrdinal, 0);

                    auto &output_y =

                        output_(interiorFamilyFieldOrdinal, pointOrdinal, 1);

                    auto &output_z =

                        output_(interiorFamilyFieldOrdinal, pointOrdinal, 2);


                    output_x = V_x * scratchL_2ipjp1(k);

                    output_y = V_y * scratchL_2ipjp1(k);

                    output_z = V_z * scratchL_2ipjp1(k);


                    interiorFamilyFieldOrdinal +=

                        numInteriorFamilies; // increment due to the

                                             // interleaving.

                  }

                }

              }

            }

          } // interior functions block


        } // end OPERATOR_VALUE

          break;

        case OPERATOR_DIV:

        {

          // rename the scratch memory to match our usage here:

          auto &scratchP = scratch0;

          auto &scratchP_2ip1 = scratch1;


          // following ESEAS, we interleave the face families.  This groups all the face dofs of a given degree together.

          ordinal_type fieldOrdinal = 0;

          for (int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)

          {

            const int max_ij_sum = polyOrder_ - 1;

            computeFaceLegendre(scratchP, faceOrdinal, lambda);

            OutputScalar divWeight;

            computeFaceDivWeight(divWeight, faceOrdinal, lambda_dx, lambda_dy, lambda_dz);

            for (int ij_sum=0; ij_sum <= max_ij_sum; ij_sum++)

            {

              for (int i=0; i<=ij_sum; i++)

              {

                const int j = ij_sum - i; // j >= 0


                computeFaceJacobi(scratchP_2ip1, faceOrdinal, i, lambda);

                auto &outputValue = output_(fieldOrdinal,pointOrdinal);

                faceFunctionDiv(outputValue, i, j, scratchP, scratchP_2ip1,

                                divWeight, lambda);


                fieldOrdinal++;

              } // i

            } // ij_sum

          } // faceOrdinal


          // interior functions

          {

            // rename the scratch memory to match our usage here:

            auto &scratchL_2ipjp1 = scratch2;

            auto &scratchP_2ipjp1 = scratch3;


            const int interiorFieldOrdinalOffset = numFaceFunctions_;

            for (int interiorFamilyOrdinal=1; interiorFamilyOrdinal<=numInteriorFamilies; interiorFamilyOrdinal++)

            {

              // following ESEAS, we interleave the interior families.  This groups all the interior dofs of a given degree together.


              const ordinal_type relatedFaceOrdinal = faceOrdinalForInterior_[interiorFamilyOrdinal-1];


              computeFaceLegendreForInterior(scratchP,

                                             interiorFamilyOrdinal - 1, lambda);

              OutputScalar divWeight;

              computeFaceDivWeight(divWeight, relatedFaceOrdinal, lambda_dx, lambda_dy, lambda_dz);


              OutputScalar vectorWeight_x, vectorWeight_y, vectorWeight_z;

              computeFaceVectorWeight(vectorWeight_x, vectorWeight_y, vectorWeight_z, relatedFaceOrdinal, lambda, lambda_dx, lambda_dy, lambda_dz);


              ordinal_type interiorFieldOrdinal = interiorFieldOrdinalOffset + interiorFamilyOrdinal - 1;


              const ordinal_type min_ijk_sum = 1;

              const ordinal_type max_ijk_sum = polyOrder_-1;

              const ordinal_type min_ij_sum  = 0;

              const ordinal_type min_k       = 1;

              const ordinal_type min_j       = 0;

              const ordinal_type min_i       = 0;

              for (int ijk_sum=min_ijk_sum; ijk_sum <= max_ijk_sum; ijk_sum++)

              {

                for (int ij_sum=min_ij_sum; ij_sum<=ijk_sum-min_k; ij_sum++)

                {

                  for (int i=min_i; i<=ij_sum-min_j; i++)

                  {

                    const ordinal_type j = ij_sum-i;

                    const ordinal_type k = ijk_sum - ij_sum;

                    computeFaceJacobiForInterior(

                        scratchP_2ip1, interiorFamilyOrdinal - 1, i, lambda);


                    OutputScalar faceDiv;

                    faceFunctionDiv(faceDiv, i, j, scratchP, scratchP_2ip1,

                                    divWeight, lambda);


                    OutputScalar faceValue_x, faceValue_y, faceValue_z;


                    faceFunctionValue(faceValue_x, faceValue_y, faceValue_z, i,

                                      j, scratchP, scratchP_2ip1,

                                      vectorWeight_x, vectorWeight_y,

                                      vectorWeight_z, lambda);

                    computeInteriorJacobi(scratchP_2ipjp1, i, j,

                                          interiorFamilyOrdinal - 1, lambda);


                    computeInteriorIntegratedJacobi(scratchL_2ipjp1, i, j,

                                                    interiorFamilyOrdinal - 1,

                                                    lambda);


                    OutputScalar L_2ipjp1_k_dx, L_2ipjp1_k_dy, L_2ipjp1_k_dz;

                    gradInteriorIntegratedJacobi(

                        L_2ipjp1_k_dx, L_2ipjp1_k_dy, L_2ipjp1_k_dz,

                        interiorFamilyOrdinal - 1, j, k, scratchP_2ipjp1,

                        lambda, lambda_dx, lambda_dy, lambda_dz);


                    auto &outputDiv = output_(interiorFieldOrdinal, pointOrdinal);

                    interiorFunctionDiv(outputDiv, scratchL_2ipjp1(k), faceDiv,

                                        L_2ipjp1_k_dx, L_2ipjp1_k_dy,

                                        L_2ipjp1_k_dz, faceValue_x, faceValue_y,

                                        faceValue_z);


                    interiorFieldOrdinal += numInteriorFamilies;  // increment due to the interleaving.

                  }

                }

              }

            }

          } // interior functions block

        } // OPERATOR_DIV

          break;

        case OPERATOR_GRAD:

        case OPERATOR_D1:

        case OPERATOR_D2:

        case OPERATOR_D3:

        case OPERATOR_D4:

        case OPERATOR_D5:

        case OPERATOR_D6:

        case OPERATOR_D7:

        case OPERATOR_D8:

        case OPERATOR_D9:

        case OPERATOR_D10:

          INTREPID2_TEST_FOR_ABORT( true,

                                   ">>> ERROR: (Intrepid2::Hierarchical_HDIV_TET_Functor) Unsupported differential operator");

        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;

    }

  };


  template<typename DeviceType,

           typename OutputScalar = double,

           typename PointScalar  = double,

           bool useCGBasis = true> // if useCGBasis is false, all basis functions will be associated with the interior


  class HierarchicalBasis_HDIV_TET

  : public Basis<DeviceType,OutputScalar,PointScalar>

  {

  public:

    using BasisBase = Basis<DeviceType,OutputScalar,PointScalar>;


    using typename BasisBase::OrdinalTypeArray1DHost;

    using typename BasisBase::OrdinalTypeArray2DHost;


    using typename BasisBase::OutputViewType;

    using typename BasisBase::PointViewType;

    using typename BasisBase::ScalarViewType;


    using typename BasisBase::ExecutionSpace;


  protected:

    int polyOrder_; // the maximum order of the polynomial

  public:


    HierarchicalBasis_HDIV_TET(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)

    :

    polyOrder_(polyOrder)

    {

      const shards::CellTopology cellTopo(shards::getCellTopologyData<shards::Tetrahedron<> >());

      this->basisCellTopologyKey_ = shards::Tetrahedron<>::key;

      const int numFaces          = cellTopo.getFaceCount();


      const int numVertexFunctions   = 0;

      const int numEdgeFunctions     = 0;

      const int numFaceFunctions     = numFaces * polyOrder * (polyOrder+1) / 2;  // 4 faces; 2 families, each with p*(p-1)/2 functions per face

      const int numInteriorFunctionsPerFamily = (polyOrder-1)*polyOrder*(polyOrder+1)/6;

      const int numInteriorFunctions = numInteriorFunctionsPerFamily * 3; // 3 families of interior functions

      this->basisCardinality_  = numVertexFunctions + numEdgeFunctions + numFaceFunctions + numInteriorFunctions;

      this->basisDegree_       = polyOrder;


      this->basisType_         = BASIS_FEM_HIERARCHICAL;

      this->basisCoordinates_  = COORDINATES_CARTESIAN;

      this->functionSpace_     = FUNCTION_SPACE_HDIV;


      const int degreeLength = 1;

      this->fieldOrdinalPolynomialDegree_ = OrdinalTypeArray2DHost("Hierarchical H(div) triangle polynomial degree lookup", this->basisCardinality_, degreeLength);


      // **** vertex functions **** //

      // no vertex functions in H(div)


      // **** edge functions **** //

      // no edge functions in H(div)


      // **** face functions **** //

      const int max_ij_sum = polyOrder-1;

      ordinal_type fieldOrdinal = 0;

      for (int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)

      {

        for (int ij_sum=0; ij_sum <= max_ij_sum; ij_sum++)

        {

          for (int i=0; i<=ij_sum; i++)

          {

            this->fieldOrdinalPolynomialDegree_(fieldOrdinal,0) = ij_sum+1;

            fieldOrdinal++;

          }

        }

      }

      INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinal != numEdgeFunctions + numFaceFunctions, std::invalid_argument, "Internal error: basis enumeration is incorrect");


      const int numInteriorFamilies = 3;

      const int interiorFieldOrdinalOffset = fieldOrdinal;

      const ordinal_type min_ijk_sum = 1;

      const ordinal_type max_ijk_sum = polyOrder_-1;

      const ordinal_type min_ij_sum  = 0;

      const ordinal_type min_k       = 1;

      const ordinal_type min_j       = 0;

      const ordinal_type min_i       = 0;

      for (int interiorFamilyOrdinal=1; interiorFamilyOrdinal<=numInteriorFamilies; interiorFamilyOrdinal++)

      {

        // following ESEAS, we interleave the interior families.  This groups all the interior dofs of a given degree together.

        fieldOrdinal = interiorFieldOrdinalOffset + interiorFamilyOrdinal - 1;

        for (int ijk_sum=min_ijk_sum; ijk_sum <= max_ijk_sum; ijk_sum++)

        {

          for (int ij_sum=min_ij_sum; ij_sum<=ijk_sum-min_k; ij_sum++)

          {

            for (int i=min_i; i<=ij_sum-min_j; i++)

            {

              this->fieldOrdinalPolynomialDegree_(fieldOrdinal,0) = ijk_sum+1;

              fieldOrdinal += numInteriorFamilies; // increment due to the interleaving.

            }

          }

        }

        fieldOrdinal = fieldOrdinal - numInteriorFamilies + 1; // due to the interleaving increment, we've gone numInteriorFamilies past the last interior ordinal.  Set fieldOrdinal to be one past.

      }


      INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinal != 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 ordinal_type faceDim = 2, volumeDim = 3;


        if (useCGBasis) {

          {

            int tagNumber = 0;

            const int numFunctionsPerFace = numFaceFunctions / numFaces;

            for (int faceOrdinalESEAS=0; faceOrdinalESEAS<numFaces; faceOrdinalESEAS++)

            {

              int faceOrdinalIntrepid2 = intrepid2FaceOrdinals[faceOrdinalESEAS];

              for (int functionOrdinal=0; functionOrdinal<numFunctionsPerFace; functionOrdinal++)

              {

                tagView(tagNumber*tagSize+0) = faceDim;               // face dimension

                tagView(tagNumber*tagSize+1) = faceOrdinalIntrepid2;  // face id

                tagView(tagNumber*tagSize+2) = functionOrdinal;       // local dof id

                tagView(tagNumber*tagSize+3) = numFunctionsPerFace;   // total number of dofs on this face

                tagNumber++;

              }

            }


            // interior

            for (int functionOrdinal=0; functionOrdinal<numInteriorFunctions; functionOrdinal++)

            {

              tagView(tagNumber*tagSize+0) = volumeDim;            // interior dimension

              tagView(tagNumber*tagSize+1) = 0;                    // volume id

              tagView(tagNumber*tagSize+2) = functionOrdinal;      // local dof id

              tagView(tagNumber*tagSize+3) = numInteriorFunctions; // total number of interior dofs

              tagNumber++;

            }

          }

        }

        else

        {

          // DG basis: all functions are associated with interior

          for (ordinal_type i=0;i<cardinality;++i) {

            tagView(i*tagSize+0) = volumeDim;   // face dimension

            tagView(i*tagSize+1) = 0;           // interior/volume id

            tagView(i*tagSize+2) = i;           // local dof id

            tagView(i*tagSize+3) = cardinality; // total number of dofs on this face

          }

        }


        // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:

        // tags are constructed on host

        this->setOrdinalTagData(this->tagToOrdinal_,

                                this->ordinalToTag_,

                                tagView,

                                this->basisCardinality_,

                                tagSize,

                                posScDim,

                                posScOrd,

                                posDfOrd);

      }

    }


    const char* getName() const override {

      return "Intrepid2_HierarchicalBasis_HDIV_TET";

    }


    virtual bool requireOrientation() const override {

      return (this->getDegree() > 2);

    }


    // 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 BasisBase::getValues;


    virtual void getValues( OutputViewType outputValues, const PointViewType inputPoints,

                           const EOperator operatorType = OPERATOR_VALUE ) const override

    {

      auto numPoints = inputPoints.extent_int(0);


      using FunctorType = Hierarchical_HDIV_TET_Functor<DeviceType, OutputScalar, PointScalar, OutputViewType, PointViewType>;


      FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_);


      const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();

      const int pointVectorSize  = getVectorSizeForHierarchicalParallelism<PointScalar>();

      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<ExecutionSpace>(numPoints,teamSize,vectorSize);

      Kokkos::parallel_for("Hierarchical_HDIV_TET_Functor", policy , functor);

    }


    BasisPtr<DeviceType,OutputScalar,PointScalar>


      getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override

    {

      using HVOL_Tri = LegendreBasis_HVOL_TRI<DeviceType,OutputScalar,PointScalar>;

      if (subCellDim == 2)

      {

        return Teuchos::rcp(new HVOL_Tri(this->basisDegree_-1));

      }

      INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds");

    }


    virtual BasisPtr<typename Kokkos::HostSpace::device_type, OutputScalar, PointScalar>


    getHostBasis() const override {

      using HostDeviceType = typename Kokkos::HostSpace::device_type;

      using HostBasisType  = HierarchicalBasis_HDIV_TET<HostDeviceType, OutputScalar, PointScalar, useCGBasis>;

      return Teuchos::rcp( new HostBasisType(polyOrder_) );

    }


  };


} // end namespace Intrepid2


#endif /* Intrepid2_HierarchicalBasis_HDIV_TET_h */

Intrepid2_Basis.hpp
Header file for the abstract base class Intrepid2::Basis.

Intrepid2::BasisPtr
Teuchos::RCP< Basis< DeviceType, OutputType, PointType > > BasisPtr
Basis Pointer.
Definition Intrepid2_Basis.hpp:43

Intrepid2::device_assert
KOKKOS_INLINE_FUNCTION void device_assert(bool val)
Definition Intrepid2_DeviceAssert.hpp:31

Intrepid2_LegendreBasis_HVOL_TRI.hpp
H(vol) basis on the triangle based on integrated Legendre polynomials.

Intrepid2_Polynomials.hpp
Free functions, callable from device code, that implement various polynomials useful in basis definit...

Intrepid2_Utils.hpp
Header function for Intrepid2::Util class and other utility functions.

Intrepid2::getScalarDimensionForView
KOKKOS_INLINE_FUNCTION constexpr unsigned getScalarDimensionForView(const ViewType &view)
Returns the size of the Scalar dimension for the View. This is 0 for non-AD types....
Definition Intrepid2_Utils.hpp:805

Intrepid2::Basis
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Definition Intrepid2_Basis.hpp:89

Intrepid2::Basis::basisCoordinates_
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
Definition Intrepid2_Basis.hpp:183

Intrepid2::Basis::PointViewType
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Definition Intrepid2_Basis.hpp:338

Intrepid2::Basis::getValues
virtual KOKKOS_INLINE_FUNCTION void getValues(OutputViewType, const PointViewType, const EOperator, const typename Kokkos::TeamPolicy< ExecutionSpace >::member_type &teamMember, const typename ExecutionSpace::scratch_memory_space &scratchStorage, const ordinal_type subcellDim=-1, const ordinal_type subcellOrdinal=-1) const
Team-level evaluation of basis functions on a reference cell.
Definition Intrepid2_Basis.hpp:426

Intrepid2::Basis::OutputViewType
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
Definition Intrepid2_Basis.hpp:334

Intrepid2::Basis::basisType_
EBasis basisType_
Type of the basis.
Definition Intrepid2_Basis.hpp:179

Intrepid2::Basis::basisDegree_
ordinal_type basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
Definition Intrepid2_Basis.hpp:166

Intrepid2::Basis::getDegree
ordinal_type getDegree() const
Returns the degree of the basis.
Definition Intrepid2_Basis.hpp:815

Intrepid2::Basis::setOrdinalTagData
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)
Fills ordinalToTag_ and tagToOrdinal_ by basis-specific tag data.
Definition Intrepid2_Basis.hpp:233

Intrepid2::Basis::ScalarViewType
Kokkos::DynRankView< scalarType, Kokkos::LayoutStride, DeviceType > ScalarViewType
View type for scalars.
Definition Intrepid2_Basis.hpp:342

Intrepid2::Basis::ordinalToTag_
OrdinalTypeArray2DHost ordinalToTag_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
Definition Intrepid2_Basis.hpp:204

Intrepid2::Basis::OrdinalTypeArray2DHost
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
Definition Intrepid2_Basis.hpp:129

Intrepid2::Basis::basisCardinality_
ordinal_type basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom.
Definition Intrepid2_Basis.hpp:162

Intrepid2::Basis::tagToOrdinal_
OrdinalTypeArray3DHost tagToOrdinal_
DoF tag to ordinal lookup table.
Definition Intrepid2_Basis.hpp:217

Intrepid2::Basis::OrdinalTypeArray1DHost
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Definition Intrepid2_Basis.hpp:125

Intrepid2::Basis::basisCellTopologyKey_
unsigned basisCellTopologyKey_
Identifier of the base topology of the cells for which the basis is defined. See the Shards package f...
Definition Intrepid2_Basis.hpp:175

Intrepid2::Basis::ExecutionSpace
typename DeviceType::execution_space ExecutionSpace
(Kokkos) Execution space for basis.
Definition Intrepid2_Basis.hpp:97

Intrepid2::Basis::fieldOrdinalPolynomialDegree_
OrdinalTypeArray2DHost fieldOrdinalPolynomialDegree_
Polynomial degree for each degree of freedom. Only defined for hierarchical bases right now....
Definition Intrepid2_Basis.hpp:304

Intrepid2::Basis::functionSpace_
EFunctionSpace functionSpace_
The function space in which the basis is defined.
Definition Intrepid2_Basis.hpp:187

Intrepid2::HierarchicalBasis_HDIV_TET
For mathematical details of the construction, see:
Definition Intrepid2_HierarchicalBasis_HDIV_TET.hpp:701

Intrepid2::HierarchicalBasis_HDIV_TET::PointViewType
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Definition Intrepid2_Basis.hpp:338

Intrepid2::HierarchicalBasis_HDIV_TET::OutputViewType
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
Definition Intrepid2_Basis.hpp:334

Intrepid2::HierarchicalBasis_HDIV_TET::getName
const char * getName() const override
Returns basis name.
Definition Intrepid2_HierarchicalBasis_HDIV_TET.hpp:868

Intrepid2::HierarchicalBasis_HDIV_TET::getValues
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
Definition Intrepid2_HierarchicalBasis_HDIV_TET.hpp:901

Intrepid2::HierarchicalBasis_HDIV_TET::getHostBasis
virtual BasisPtr< typename Kokkos::HostSpace::device_type, OutputScalar, PointScalar > getHostBasis() const override
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
Definition Intrepid2_HierarchicalBasis_HDIV_TET.hpp:943

Intrepid2::HierarchicalBasis_HDIV_TET::getSubCellRefBasis
BasisPtr< DeviceType, OutputScalar, PointScalar > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
returns the basis associated to a subCell.
Definition Intrepid2_HierarchicalBasis_HDIV_TET.hpp:928

Intrepid2::HierarchicalBasis_HDIV_TET::OrdinalTypeArray2DHost
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
Definition Intrepid2_Basis.hpp:129

Intrepid2::HierarchicalBasis_HDIV_TET::requireOrientation
virtual bool requireOrientation() const override
True if orientation is required.
Definition Intrepid2_HierarchicalBasis_HDIV_TET.hpp:874

Intrepid2::HierarchicalBasis_HDIV_TET::HierarchicalBasis_HDIV_TET
HierarchicalBasis_HDIV_TET(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
Definition Intrepid2_HierarchicalBasis_HDIV_TET.hpp:721

Intrepid2::HierarchicalBasis_HDIV_TET::OrdinalTypeArray1DHost
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Definition Intrepid2_Basis.hpp:125

Intrepid2::LegendreBasis_HVOL_TRI
Basis defining Legendre basis on the line, a polynomial subspace of H(vol) on the line: extension to ...
Definition Intrepid2_LegendreBasis_HVOL_TRI.hpp:155

Intrepid2::Hierarchical_HDIV_TET_Functor
Functor for computing values for the HierarchicalBasis_HDIV_TET class.
Definition Intrepid2_HierarchicalBasis_HDIV_TET.hpp:37

Intrepid2::Hierarchical_HDIV_TET_Functor::computeFaceJacobiForInterior
KOKKOS_INLINE_FUNCTION void computeFaceJacobiForInterior(OutputScratchView &P_2ip1, const ordinal_type &zeroBasedInteriorFamilyOrdinal, const ordinal_type &i, const PointScalar *lambda) const
The face functions we compute for interior blending can have a different orientation than the ones us...
Definition Intrepid2_HierarchicalBasis_HDIV_TET.hpp:161