ordinal_type compareToAnalytic( std::ifstream &inputFile,
                                  const Kokkos::DynRankView<ValueType,testMatProperties...> testMat,
                                  const ValueType reltol,
                                  const ordinal_type iprint,
                                  const TypeOfExactData analyticDataType ) {
    INTREPID2_TEST_FOR_EXCEPTION( testMat.rank() != 2, std::invalid_argument,
                                  ">>> ERROR (compareToAnalytic): testMat must have rank 2");

    Teuchos::RCP<std::ostream> outStream;
    Teuchos::oblackholestream  outNothing;
    if (iprint > 0)
      outStream = Teuchos::rcp(&std::cout, false);
    else
      outStream = Teuchos::rcp(&outNothing, false);

    // Save the format state of the original std::cout.
    Teuchos::oblackholestream oldFormatState;
    oldFormatState.copyfmt(std::cout);

    std::string line;
    ValueType testentry;
    ValueType abstol;
    ValueType absdiff;
    ordinal_type i = 0, j = 0;
    ordinal_type err = 0;

    while (! inputFile.eof() && i < static_cast<ordinal_type>(testMat.dimension(0)) ) {
      std::getline(inputFile,line);
      std::istringstream linestream(line);
      std::string chunk;
      j = 0;
      while( linestream >> chunk ) {
        ordinal_type num1;
        ordinal_type num2;
        std::string::size_type loc = chunk.find( "/", 0);
        if( loc != std::string::npos ) {
          chunk.replace( loc, 1, " ");
          std::istringstream chunkstream(chunk);
          chunkstream >> num1;
          chunkstream >> num2;
          testentry = (ValueType)(num1)/(ValueType)(num2);
          abstol = ( std::fabs(testentry) < reltol ? reltol : std::fabs(reltol*testentry) );
          absdiff = std::fabs(testentry - testMat(i, j));
          if (absdiff > abstol) {
            ++err;
            *outStream << "FAILURE --> ";
          }
          *outStream << "entry[" << i << "," << j << "]:" << "   "
                     << testMat(i, j) << "   " << num1 << "/" << num2 << "   "
                     << absdiff << "   " << "<?" << "   " << abstol << "\n";
        }
        else {
          std::istringstream chunkstream(chunk);
          if (analyticDataType == INTREPID2_UTILS_FRACTION) {
            chunkstream >> num1;
            testentry = (ValueType)(num1);
          }
          else if (analyticDataType == INTREPID2_UTILS_SCALAR)

    void
    Basis_HGRAD_LINE_Cn_FEM::
    getValues(       Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
               const Kokkos::DynRankView<inputPointValueType, inputPointProperties...>  inputPoints,
               const Kokkos::DynRankView<vinvValueType,       vinvProperties...>        vinv,
               const EOperator operatorType ) {
      typedef          Kokkos::DynRankView<outputValueValueType,outputValueProperties...>         outputValueViewType;
      typedef          Kokkos::DynRankView<inputPointValueType, inputPointProperties...>          inputPointViewType;
      typedef          Kokkos::DynRankView<vinvValueType,       vinvProperties...>                vinvViewType;
      typedef typename ExecSpace<typename inputPointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;

      // loopSize corresponds to cardinality
      const auto loopSizeTmp1 = (inputPoints.dimension(0)/numPtsPerEval);
      const auto loopSizeTmp2 = (inputPoints.dimension(0)%numPtsPerEval != 0);
      const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
      Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);

      typedef typename inputPointViewType::value_type inputPointType;

      const ordinal_type cardinality = outputValues.dimension(0);

      auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
      typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
      workViewType  work(Kokkos::view_alloc("Basis_HGRAD_LINE_Cn_FEM::getValues::work", vcprop), cardinality, inputPoints.dimension(0));

      switch (operatorType) {
      case OPERATOR_VALUE: {
        typedef Functor<outputValueViewType,inputPointViewType,vinvViewType,workViewType,
            OPERATOR_VALUE,numPtsPerEval> FunctorType;
        Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work) );
        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: {
        typedef Functor<outputValueViewType,inputPointViewType,vinvViewType,workViewType,
            OPERATOR_Dn,numPtsPerEval> FunctorType;
        Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work,
                                                  getOperatorOrder(operatorType)) );
        break;
      }
      default: {
        INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument,
                                      ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM): Operator type not implemented" );
        //break; commented out because this always throws
      }
      }
    }

  void 
  PointTools::
  getGaussPoints( /**/  Kokkos::DynRankView<pointValueType,pointProperties...> points,
                  const ordinal_type order ) {
#ifdef HAVE_INTREPID2_DEBUG    
    INTREPID2_TEST_FOR_EXCEPTION( points.rank() != 2,
                                  std::invalid_argument ,
                                  ">>> ERROR (PointTools::getGaussPoints): points rank must be 1." );
    INTREPID2_TEST_FOR_EXCEPTION( order < 0,
                                  std::invalid_argument ,
                                  ">>> ERROR (PointTools::getGaussPoints): order must be positive value." );
#endif
    const ordinal_type np = order + 1;
    const double alpha = 0.0, beta = 0.0;
    
    // until view and dynrankview inter-operatible, we use views in a consistent way
    Kokkos::View<pointValueType*,Kokkos::HostSpace> 
      zHost("PointTools::getGaussPoints::z", np), 
      wHost("PointTools::getGaussPoints::w", np);
    
    // sequential means that the code is decorated with KOKKOS_INLINE_FUNCTION 
    // and it does not invoke parallel for inside (cheap operation), which means 
    // that gpu memory is not accessible unless this is called inside of functor.
    Polylib::Serial::Cubature<POLYTYPE_GAUSS>::getValues(zHost, wHost, np, alpha, beta);

    typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
    auto pts = Kokkos::subview( points, range_type(0,np), 0 );
    // should be fixed after view and dynrankview are inter-operatible
    auto z   = Kokkos::DynRankView<pointValueType,Kokkos::HostSpace>(zHost.data(), np);
    Kokkos::deep_copy(pts, z);
  }

void 
Basis_HCURL_TRI_I1_FEM<SpT,OT,PT>::Internal::
getDofCoords( Kokkos::DynRankView<dofCoordValueType,dofCoordProperties...> dofCoords ) const {
#ifdef HAVE_INTREPID2_DEBUG
    // Verify rank of output array.
    INTREPID2_TEST_FOR_EXCEPTION( dofCoords.rank() != 2, std::invalid_argument,
                                  ">>> ERROR: (Intrepid2::Basis_HCURL_TRI_I1_FEM::getDofCoords) rank = 2 required for dofCoords array");
    // Verify 0th dimension of output array.
    INTREPID2_TEST_FOR_EXCEPTION( dofCoords.dimension(0) != obj_->basisCardinality_, std::invalid_argument,
                                  ">>> ERROR: (Intrepid2::Basis_HCURL_TRI_I1_FEM::getDofCoords) mismatch in number of dof and 0th dimension of dofCoords array");
    // Verify 1st dimension of output array.
    INTREPID2_TEST_FOR_EXCEPTION( dofCoords.dimension(1) != obj_->basisCellTopology_.getDimension(), std::invalid_argument,
                                  ">>> ERROR: (Intrepid2::Basis_HCURL_TRI_I1_FEM::getDofCoords) incorrect reference cell (1st) dimension in dofCoords array");
#endif
    Kokkos::deep_copy(dofCoords, obj_->dofCoords_);
}

    void
    Basis_HGRAD_QUAD_C1_FEM::
    getValues(       Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
               const Kokkos::DynRankView<inputPointValueType, inputPointProperties...>  inputPoints,
               const EOperator operatorType )  {
      typedef          Kokkos::DynRankView<outputValueValueType,outputValueProperties...>         outputValueViewType;
      typedef          Kokkos::DynRankView<inputPointValueType, inputPointProperties...>          inputPointViewType;
      typedef typename ExecSpace<typename inputPointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;

      // Number of evaluation points = dim 0 of inputPoints
      const auto loopSize = inputPoints.dimension(0);
      Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);

      switch (operatorType) {

      case OPERATOR_VALUE: {
        typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
        Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
        break;
      }

      case OPERATOR_GRAD:
      case OPERATOR_D1: {
        typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
        Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
        break;
      }
      case OPERATOR_CURL: {
        typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_CURL> FunctorType;
        Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
        break;
      }
      case OPERATOR_DIV: {
        INTREPID2_TEST_FOR_EXCEPTION( operatorType == OPERATOR_DIV, std::invalid_argument,
                                      ">>> ERROR (Basis_HGRAD_QUAD_C1_FEM): DIV is invalid operator for rank-0 (scalar) functions in 2D");
        break;
      }
      case OPERATOR_D2: {
        typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D2> FunctorType;
        Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
        break;
      }
      case OPERATOR_D3:
      case OPERATOR_D4:
      case OPERATOR_D5:
      case OPERATOR_D6:
      case OPERATOR_D7:
      case OPERATOR_D8:
      case OPERATOR_D9:
      case OPERATOR_D10: {
        typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
        Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
        break;
      }
      default: {
        INTREPID2_TEST_FOR_EXCEPTION( !Intrepid2::isValidOperator(operatorType), std::invalid_argument,
                                      ">>> ERROR (Basis_HGRAD_QUAD_C1_FEM): Invalid operator type");
      }
      }
    }

  void
  Basis_HGRAD_LINE_C1_FEM<SpT>::  
  getDofCoords( Kokkos::DynRankView<dofCoordValueType,dofCoordProperties...> dofCoords ) const {

#ifdef HAVE_INTREPID2_DEBUG
    // Verify rank of output array.
    INTREPID2_TEST_FOR_EXCEPTION( dofCoords.rank() != 2, std::invalid_argument,
                                  ">>> ERROR: (Intrepid2::Basis_HGRAD_LINE_C1_FEM::getDofCoords) rank = 2 required for dofCoords array");
    // Verify 0th dimension of output array.
    INTREPID2_TEST_FOR_EXCEPTION( dofCoords.dimension(0) != this->basisCardinality_, std::invalid_argument,
                                  ">>> ERROR: (Intrepid2::Basis_HGRAD_LINE_C1_FEM::getDofCoords) mismatch in number of dof and 0th dimension of dofCoords array");
    // Verify 1st dimension of output array.
    INTREPID2_TEST_FOR_EXCEPTION( dofCoords.dimension(1) != this->basisCellTopology_.getDimension(), std::invalid_argument,
                                  ">>> ERROR: (Intrepid2::Basis_HGRAD_LINE_C1_FEM::getDofCoords) incorrect reference cell (1st) dimension in dofCoords array");
#endif
    dofCoords(0,0) = -1.0;   
    dofCoords(1,0) =  1.0;   
  }

  void
  ArrayTools<SpT>::Internal::
  dotMultiply( /**/  Kokkos::DynRankView<outputValueType,    outputProperties...>      output,
               const Kokkos::DynRankView<leftInputValueType, leftInputProperties...>   leftInput,
               const Kokkos::DynRankView<rightInputValueType,rightInputProperties...>  rightInput, 
               const bool hasField ) {

    typedef Kokkos::DynRankView<outputValueType,    outputProperties...>      outputViewType;
    typedef Kokkos::DynRankView<leftInputValueType, leftInputProperties...>   leftInputViewType;
    typedef Kokkos::DynRankView<rightInputValueType,rightInputProperties...>  rightInputViewType;
    typedef FunctorArrayTools::F_dotMultiply<outputViewType, leftInputViewType, rightInputViewType> FunctorType;
    typedef typename ExecSpace< typename leftInputViewType::execution_space , SpT >::ExecSpaceType ExecSpaceType;

    const size_type loopSize = ( hasField ? output.dimension(0)*output.dimension(1)*output.dimension(2) :
                                 /**/       output.dimension(0)*output.dimension(1) );
    Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
    Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, hasField) );
  }

  KOKKOS_INLINE_FUNCTION
  void
  Basis_HGRAD_TRI_C1_FEM<SpT,OT,PT>::Serial<opType>::
  getValues( /**/  Kokkos::DynRankView<outputValueValueType,outputValueProperties...> output,
             const Kokkos::DynRankView<inputPointValueType, inputPointProperties...>  input ) {
    switch (opType) {
    case OPERATOR_VALUE: {
      const auto x = input(0);
      const auto y = input(1);
      
      // output is a rank-2 array with dimensions (basisCardinality_)
      output(0) = 1.0 - x - y;
      output(1) = x;
      output(2) = y;
      break;
    }
    case OPERATOR_GRAD: {
      // output is a rank-3 array with dimensions (basisCardinality_, spaceDim)
      output(0, 0) = -1.0;
      output(0, 1) = -1.0;
      
      output(1, 0) =  1.0;
      output(1, 1) =  0.0;
      
      output(2, 0) =  0.0;
      output(2, 1) =  1.0;
      break;
    }
    case OPERATOR_CURL: {
      output(0, 0) = -1.0;
      output(0, 1) =  1.0;
      
      output(1, 0) =  0.0;
      output(1, 1) = -1.0;
      
      output(2, 0) =  1.0;
      output(2, 1) =  0.0;
      break;
    }
    case OPERATOR_MAX: {
      const auto jend = output.dimension(1);
      const auto iend = output.dimension(0);

      for (size_type j=0;j<jend;++j)
        for (size_type i=0;i<iend;++i)
          output(i, j) = 0.0;
      break;
    }
    default: {
      INTREPID2_TEST_FOR_ABORT( opType != OPERATOR_VALUE &&
                                opType != OPERATOR_GRAD &&
                                opType != OPERATOR_CURL &&
                                opType != OPERATOR_MAX,
                                ">>> ERROR: (Intrepid2::Basis_HGRAD_TRI_C1_FEM::Serial::getValues) operator is not supported");
    }
    }
  }

    void
    Basis_HGRAD_TRI_Cn_FEM_ORTH::
    getValues( /**/  Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
               const Kokkos::DynRankView<inputPointValueType, inputPointProperties...>  inputPoints,
               const ordinal_type order,
               const EOperator operatorType ) {
      typedef          Kokkos::DynRankView<outputValueValueType,outputValueProperties...>         outputValueViewType;
      typedef          Kokkos::DynRankView<inputPointValueType, inputPointProperties...>          inputPointViewType;
      typedef typename ExecSpace<typename inputPointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;

      // loopSize corresponds to the # of points
      const auto loopSizeTmp1 = (inputPoints.dimension(0)/numPtsPerEval);
      const auto loopSizeTmp2 = (inputPoints.dimension(0)%numPtsPerEval != 0);
      const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
      Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);

      switch (operatorType) {
      case OPERATOR_VALUE: {
        typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE,numPtsPerEval> FunctorType;
        Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints,
                                                  order) );
        break;
      }
      case OPERATOR_GRAD:
      case OPERATOR_D1: {
        typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD,numPtsPerEval> FunctorType;
        Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints,
                                                  order) );
        break;
      }
      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: {
        typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_Dn,numPtsPerEval> FunctorType;
        Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints,
                                                  order, getOperatorOrder(operatorType)) );
        break;
      }
      case OPERATOR_DIV:
      case OPERATOR_CURL: {
        INTREPID2_TEST_FOR_EXCEPTION( operatorType == OPERATOR_DIV ||
                                      operatorType == OPERATOR_CURL, std::invalid_argument,
                                      ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM_ORTH): invalid operator type (div and curl).");
        break;
      }
      default: {
        INTREPID2_TEST_FOR_EXCEPTION( !Intrepid2::isValidOperator(operatorType), std::invalid_argument,
                                      ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM_ORTH): invalid operator type");
      }
      }
    }

  void 
  PointTools::
  getLattice( /**/  Kokkos::DynRankView<pointValueType,pointProperties...> points,
              const shards::CellTopology cell,
              const ordinal_type         order,
              const ordinal_type         offset,
              const EPointType           pointType ) {
#ifdef HAVE_INTREPID2_DEBUG
    INTREPID2_TEST_FOR_EXCEPTION( points.rank() != 2,
                                  std::invalid_argument ,
                                  ">>> ERROR (PointTools::getLattice): points rank must be 2." );
    INTREPID2_TEST_FOR_EXCEPTION( order < 0 || offset < 0,
                                  std::invalid_argument ,
                                  ">>> ERROR (PointTools::getLattice): order and offset must be positive values." );

    const size_type latticeSize = getLatticeSize( cell, order, offset );
    const size_type spaceDim = cell.getDimension();
    
    INTREPID2_TEST_FOR_EXCEPTION( points.dimension(0) != latticeSize ||
                                  points.dimension(1) != spaceDim,
                                  std::invalid_argument ,
                                  ">>> ERROR (PointTools::getLattice): dimension does not match to lattice size." );
#endif

    // const auto latticeSize = getLatticeSize( cell, order, offset );
    // const auto spaceDim = cell.getDimension();
    
    // // the interface assumes that the input array follows the cell definition
    // // so, let's match all dimensions according to the cell specification
    // typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
    // auto pts = Kokkos::subview( points, 
    //                                    range_type(0, latticeSize), 
    //                                    range_type(0, spaceDim) );   
    switch (pointType) {
    case POINTTYPE_EQUISPACED:  getEquispacedLattice( points, cell, order, offset ); break;
    case POINTTYPE_WARPBLEND:   getWarpBlendLattice ( points, cell, order, offset ); break;
    default: {
      INTREPID2_TEST_FOR_EXCEPTION( true ,
                                    std::invalid_argument ,
                                    ">>> ERROR (PointTools::getLattice): invalid EPointType." );
    }
    }
  }

  void
  Basis_HGRAD_LINE_C1_FEM<SpT,OT,PT>::Internal::
  getValues( /**/  Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
             const Kokkos::DynRankView<inputPointValueType, inputPointProperties...>  inputPoints,
             const EOperator operatorType ) const {
#ifdef HAVE_INTREPID2_DEBUG
    Intrepid2::getValues_HGRAD_Args(outputValues,
                                    inputPoints,
                                    operatorType,
                                    obj_->getBaseCellTopology(),
                                    obj_->getCardinality() );
#endif

    typedef          Kokkos::DynRankView<outputValueValueType,outputValueProperties...>         outputValueViewType;
    typedef          Kokkos::DynRankView<inputPointValueType, inputPointProperties...>          inputPointViewType;
    typedef typename ExecSpace<typename inputPointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;

    // Number of evaluation points = dim 0 of inputPoints
    const auto loopSize = inputPoints.dimension(0);
    Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);

    switch (operatorType) {

    case OPERATOR_VALUE: {
      typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
      Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
      break;
    }
    case OPERATOR_GRAD:
    case OPERATOR_DIV:
    case OPERATOR_CURL:
    case OPERATOR_D1: {
      typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
      Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
      break;
    }
    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: {
      typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
      Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
      break;
    }
    default: {
      INTREPID2_TEST_FOR_EXCEPTION( !Intrepid2::isValidOperator(operatorType), std::invalid_argument,
                                    ">>> ERROR (Basis_HGRAD_LINE_C1_FEM): Invalid operator type");
    }
    }
  }

  /** Get the local coordinates for this field. This is independent of element
   * locations.
   *
   * \param[in,out] coords   Coordinates associated with this field type.
   */
  void 
  Intrepid2FieldPattern::
  getInterpolatoryCoordinates(const Kokkos::DynRankView<double,PHX::Device> & cellVertices,
                              Kokkos::DynRankView<double,PHX::Device> & coords) const
  {
    TEUCHOS_ASSERT(cellVertices.rank()==3);

    int numCells = cellVertices.dimension(0);

    // grab the local coordinates
    Kokkos::DynRankView<double,PHX::Device> localCoords;
    getInterpolatoryCoordinates(localCoords);

    // resize the coordinates field container
    coords = Kokkos::DynRankView<double,PHX::Device>("coords",numCells,localCoords.dimension(0),getDimension());

    if(numCells>0) {
      Intrepid2::CellTools<PHX::Device> cellTools;
      cellTools.mapToPhysicalFrame(coords,localCoords,cellVertices,intrepidBasis_->getBaseCellTopology());
    }
  }

 void
 OrientationTools<SpT>::
 getOrientation(/**/  Kokkos::DynRankView<elemOrtValueType,elemOrtProperties...> elemOrts,
                const Kokkos::DynRankView<elemNodeValueType,elemNodeProperties...> elemNodes,
                const shards::CellTopology cellTopo) {
   // small meta data modification and it uses shards; let's do this on host
   typedef typename Kokkos::Impl::is_space<SpT>::host_mirror_space::execution_space host_space_type;
   auto elemOrtsHost = Kokkos::create_mirror_view(typename host_space_type::memory_space(), elemOrts);
   auto elemNodesHost = Kokkos::create_mirror_view(typename host_space_type::memory_space(), elemNodes);
   
   const ordinal_type numCells = elemNodes.dimension(0);
   for (auto cell=0;cell<numCells;++cell) {
     const auto nodes = Kokkos::subview(elemNodesHost, cell, Kokkos::ALL());
     elemOrtsHost(cell) = Orientation::getOrientation(cellTopo, nodes);
   }
   
   Kokkos::deep_copy(elemOrts, elemOrtsHost);
 }

  void
  CellTools<SpT>::
  mapToPhysicalFrame( /**/  Kokkos::DynRankView<physPointValueType,physPointProperties...>     physPoints,
                      const Kokkos::DynRankView<refPointValueType,refPointProperties...>       refPoints,
                      const Kokkos::DynRankView<worksetCellValueType,worksetCellProperties...> worksetCell,
                      const HGradBasisPtrType basis ) {
#ifdef HAVE_INTREPID2_DEBUG
    CellTools_mapToPhysicalFrameArgs( physPoints, refPoints, worksetCell, basis->getBaseCellTopology() );
#endif
    const auto cellTopo = basis->getBaseCellTopology();
    const auto numCells = worksetCell.dimension(0);

    //points can be rank-2 (P,D), or rank-3 (C,P,D)
    const auto refPointRank = refPoints.rank();
    const auto numPoints = (refPointRank == 2 ? refPoints.dimension(0) : refPoints.dimension(1));
    const auto basisCardinality = basis->getCardinality();

    typedef Kokkos::DynRankView<physPointValueType,physPointProperties...> physPointViewType;
    typedef Kokkos::DynRankView<decltype(basis->getDummyOutputValue()),SpT> valViewType;

    valViewType vals;

    switch (refPointRank) {
    case 2: {
      // refPoints is (P,D): single set of ref. points is mapped to one or multiple physical cells
      vals = valViewType("CellTools::mapToPhysicalFrame::vals", basisCardinality, numPoints);
      basis->getValues(vals,
                       refPoints,
                       OPERATOR_VALUE);
      break;
    }
    case 3: {
      // refPoints is (C,P,D): multiple sets of ref. points are mapped to matching number of physical cells.
      vals = valViewType("CellTools::mapToPhysicalFrame::vals", numCells, basisCardinality, numPoints);
      for (size_type cell=0;cell<numCells;++cell)
        basis->getValues(Kokkos::subdynrankview( vals,      cell, Kokkos::ALL(), Kokkos::ALL() ),
                         Kokkos::subdynrankview( refPoints, cell, Kokkos::ALL(), Kokkos::ALL() ),
                         OPERATOR_VALUE);
      break;
    }
    }

    typedef Kokkos::DynRankView<worksetCellValueType,worksetCellProperties...> worksetCellViewType;
    typedef FunctorCellTools::F_mapToPhysicalFrame<physPointViewType,worksetCellViewType,valViewType> FunctorType;
    typedef typename ExecSpace<typename worksetCellViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;

    const auto loopSize = physPoints.dimension(0)*physPoints.dimension(1);
    Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
    Kokkos::parallel_for( policy, FunctorType(physPoints, worksetCell, vals) );
  }

  void
  OrientationTools<SpT>::
  modifyBasisByOrientation(/**/  Kokkos::DynRankView<outputValueType,outputProperties...> output,
                           const Kokkos::DynRankView<inputValueType, inputProperties...>  input,
                           const Kokkos::DynRankView<ortValueType,   ortProperties...>    orts,
                           const BasisPtrType basis ) {
#ifdef HAVE_INTREPID2_DEBUG
    {
      INTREPID2_TEST_FOR_EXCEPTION( input.rank() != output.rank(), std::invalid_argument,
                                    ">>> ERROR (OrientationTools::modifyBasisByOrientation): Input and output rank are not 3.");
      for (ordinal_type i=0;i<input.rank();++i)
        INTREPID2_TEST_FOR_EXCEPTION( input.dimension(i) != output.dimension(i), std::invalid_argument,
                                      ">>> ERROR (OrientationTools::modifyBasisByOrientation): Input and output dimension does not match.");

      INTREPID2_TEST_FOR_EXCEPTION( input.dimension(1) != basis->getCardinality(), std::invalid_argument,
                                    ">>> ERROR (OrientationTools::modifyBasisByOrientation): Field dimension of input/output does not match to basis cardinality.");
      INTREPID2_TEST_FOR_EXCEPTION( input.dimension(3) != basis->getBaseCellTopology().getDimension(), std::invalid_argument,
                                    ">>> ERROR (OrientationTools::modifyBasisByOrientation): Space dimension of input/output does not match to topology dimension.");
    }
#endif
    if (basis->requireOrientation()) {
      auto ordinalToTag = Kokkos::create_mirror_view(typename SpT::memory_space(), basis->getAllDofTags());
      auto tagToOrdinal = Kokkos::create_mirror_view(typename SpT::memory_space(), basis->getAllDofOrdinal());
      
      Kokkos::deep_copy(ordinalToTag, basis->getAllDofTags());
      Kokkos::deep_copy(tagToOrdinal, basis->getAllDofOrdinal());
      
      const ordinal_type 
        numCells  = output.dimension(0),
        //numBasis  = output.dimension(1),
        numPoints = output.dimension(2),
        dimBasis  = output.dimension(3);
      
      const CoeffMatrixDataViewType matData = createCoeffMatrix(basis);
      const shards::CellTopology cellTopo = basis->getBaseCellTopology();
      
      const ordinal_type 
        numVerts = cellTopo.getVertexCount(), 
        numEdges = cellTopo.getEdgeCount(),
        numFaces = cellTopo.getFaceCount();
      
      const ordinal_type intrDim = ( numEdges == 0 ? 1 : 
                                     numFaces == 0 ? 2 : 
                                     /**/            3 );
      
      for (auto cell=0;cell<numCells;++cell) {
        auto out = Kokkos::subview(output, cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
        auto in  = Kokkos::subview(input,  cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
        
        // vertex copy (no orientation)
        for (auto vertId=0;vertId<numVerts;++vertId) {
          const auto i = tagToOrdinal(0, vertId, 0);
          if (i != -1) // if dof does not exist i returns with -1
            for (auto j=0;j<numPoints;++j)
              for (auto k=0;k<dimBasis;++k)
                out(i, j, k) = in(i, j, k);
        }
        
        // interior copy
        {
          const auto ordIntr = tagToOrdinal(intrDim, 0, 0);
          if (ordIntr != -1) {
            const auto ndofIntr = ordinalToTag(ordIntr, 3);
            for (auto i=0;i<ndofIntr;++i) {
              const auto ii = tagToOrdinal(intrDim, 0, i);
              for (auto j=0;j<numPoints;++j)
                for (auto k=0;k<dimBasis;++k)
                  out(ii, j, k) = in(ii, j, k);
            }
          }
        }
        
        // edge transformation
        if (numEdges > 0) {
          ordinal_type ortEdges[12];
          orts(cell).getEdgeOrientation(ortEdges, numEdges);
          
          // apply coeff matrix
          for (auto edgeId=0;edgeId<numEdges;++edgeId) {
            const auto ordEdge = tagToOrdinal(1, edgeId, 0);
            
            if (ordEdge != -1) {
              const auto ndofEdge = ordinalToTag(ordEdge, 3);
              const auto mat = Kokkos::subview(matData, 
                                               edgeId, ortEdges[edgeId], 
                                               Kokkos::ALL(), Kokkos::ALL());
              
              for (auto j=0;j<numPoints;++j) 
                for (auto i=0;i<ndofEdge;++i) {
                  const auto ii = tagToOrdinal(1, edgeId, i);
                  
                  for (auto k=0;k<dimBasis;++k) {
                    double temp = 0.0;
                    for (auto l=0;l<ndofEdge;++l) {
                      const auto ll = tagToOrdinal(1, edgeId, l);
                      temp += mat(i,l)*in(ll, j, k);
                    }
                    out(ii, j, k) = temp;
                  }
                }
//.........这里部分代码省略.........

  void
  CellTools<SpT>::
  mapToReferenceSubcell( /**/  Kokkos::DynRankView<refSubcellPointValueType,refSubcellPointProperties...> refSubcellPoints,
                         const Kokkos::DynRankView<paramPointValueType,paramPointProperties...>           paramPoints,
                         const ordinal_type subcellDim,
                         const ordinal_type subcellOrd,
                         const shards::CellTopology parentCell ) {
#ifdef HAVE_INTREPID2_DEBUG
    INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(parentCell), std::invalid_argument,
                                  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): the specified cell topology does not have a reference cell.");

    INTREPID2_TEST_FOR_EXCEPTION( subcellDim != 1 &&
                                  subcellDim != 2, std::invalid_argument,
                                  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): method defined only for 1 and 2-dimensional subcells.");

    INTREPID2_TEST_FOR_EXCEPTION( subcellOrd <  0 ||
                                  subcellOrd >= parentCell.getSubcellCount(subcellDim), std::invalid_argument,
                                  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): subcell ordinal out of range.");

    // refSubcellPoints is rank-2 (P,D1), D1 = cell dimension
    INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.rank() != 2, std::invalid_argument,
                                  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints must have rank 2.");
    INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.dimension(1) != parentCell.getDimension(), std::invalid_argument,
                                  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints dimension (1) does not match to parent cell dimension.");

    // paramPoints is rank-2 (P,D2) with D2 = subcell dimension
    INTREPID2_TEST_FOR_EXCEPTION( paramPoints.rank() != 2, std::invalid_argument,
                                  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): paramPoints must have rank 2.");
    INTREPID2_TEST_FOR_EXCEPTION( paramPoints.dimension(1) != subcellDim, std::invalid_argument,
                                  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): paramPoints dimension (1) does not match to subcell dimension.");

    // cross check: refSubcellPoints and paramPoints: dimension 0 must match
    INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.dimension(0) < paramPoints.dimension(0), std::invalid_argument,
                                  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints dimension (0) does not match to paramPoints dimension(0).");
#endif


    const auto cellDim = parentCell.getDimension();
    const auto numPts  = paramPoints.dimension(0);

    // Get the subcell map, i.e., the coefficients of the parametrization function for the subcell

    // can i get this map from devices ?
    subcellParamViewType subcellMap;
    getSubcellParametrization( subcellMap,
                               subcellDim,
                               parentCell );

    // subcell parameterization should be small computation (numPts is small) and it should be decorated with
    // kokkos inline... let's not do this yet

    // Apply the parametrization map to every point in parameter domain
    switch (subcellDim) {
    case 2: {
      for (size_type pt=0;pt<numPts;++pt) {
        const auto u = paramPoints(pt, 0);
        const auto v = paramPoints(pt, 1);

        // map_dim(u,v) = c_0(dim) + c_1(dim)*u + c_2(dim)*v because both Quad and Tri ref faces are affine!
        for (size_type i=0;i<cellDim;++i)
          refSubcellPoints(pt, i) = subcellMap(subcellOrd, i, 0) + ( subcellMap(subcellOrd, i, 1)*u +
                                                                     subcellMap(subcellOrd, i, 2)*v );
      }
      break;
    }
    case 1: {
      for (size_type pt=0;pt<numPts;++pt) {
        const auto u = paramPoints(pt, 0);
        for (size_type i=0;i<cellDim;++i)
          refSubcellPoints(pt, i) = subcellMap(subcellOrd, i, 0) + ( subcellMap(subcellOrd, i, 1)*u );
      }
      break;
    }
    default: {
      INTREPID2_TEST_FOR_EXCEPTION( subcellDim != 1 &&
                                    subcellDim != 2, std::invalid_argument,
                                    ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): method defined only for 1 and 2-subcells");
    }
    }
  }

  KOKKOS_INLINE_FUNCTION
  void
  Basis_HGRAD_HEX_C1_FEM<SpT>::Serial<opType>::
  getValues( /**/  Kokkos::DynRankView<outputValueValueType,outputValueProperties...> output,
             const Kokkos::DynRankView<inputPointValueType, inputPointProperties...>  input ) {
    switch (opType) {
    case OPERATOR_VALUE : {
      const auto x = input(0);
      const auto y = input(1);
      const auto z = input(2);

      // output is a rank-2 array with dimensions (basisCardinality_, dim0)
      output(0) = (1.0 - x)*(1.0 - y)*(1.0 - z)/8.0;
      output(1) = (1.0 + x)*(1.0 - y)*(1.0 - z)/8.0;
      output(2) = (1.0 + x)*(1.0 + y)*(1.0 - z)/8.0;
      output(3) = (1.0 - x)*(1.0 + y)*(1.0 - z)/8.0;
          
      output(4) = (1.0 - x)*(1.0 - y)*(1.0 + z)/8.0;
      output(5) = (1.0 + x)*(1.0 - y)*(1.0 + z)/8.0;
      output(6) = (1.0 + x)*(1.0 + y)*(1.0 + z)/8.0;
      output(7) = (1.0 - x)*(1.0 + y)*(1.0 + z)/8.0;        
      break;
    }
    case OPERATOR_GRAD : {
      const auto x = input(0);
      const auto y = input(1);
      const auto z = input(2);

      // output is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
      output(0, 0) = -(1.0 - y)*(1.0 - z)/8.0;
      output(0, 1) = -(1.0 - x)*(1.0 - z)/8.0;
      output(0, 2) = -(1.0 - x)*(1.0 - y)/8.0;
          
      output(1, 0) =  (1.0 - y)*(1.0 - z)/8.0;
      output(1, 1) = -(1.0 + x)*(1.0 - z)/8.0;
      output(1, 2) = -(1.0 + x)*(1.0 - y)/8.0;
          
      output(2, 0) =  (1.0 + y)*(1.0 - z)/8.0;
      output(2, 1) =  (1.0 + x)*(1.0 - z)/8.0;
      output(2, 2) = -(1.0 + x)*(1.0 + y)/8.0;
          
      output(3, 0) = -(1.0 + y)*(1.0 - z)/8.0;
      output(3, 1) =  (1.0 - x)*(1.0 - z)/8.0;
      output(3, 2) = -(1.0 - x)*(1.0 + y)/8.0;
          
      output(4, 0) = -(1.0 - y)*(1.0 + z)/8.0;
      output(4, 1) = -(1.0 - x)*(1.0 + z)/8.0;
      output(4, 2) =  (1.0 - x)*(1.0 - y)/8.0;
          
      output(5, 0) =  (1.0 - y)*(1.0 + z)/8.0;
      output(5, 1) = -(1.0 + x)*(1.0 + z)/8.0;
      output(5, 2) =  (1.0 + x)*(1.0 - y)/8.0;
          
      output(6, 0) =  (1.0 + y)*(1.0 + z)/8.0;
      output(6, 1) =  (1.0 + x)*(1.0 + z)/8.0;
      output(6, 2) =  (1.0 + x)*(1.0 + y)/8.0;
          
      output(7, 0) = -(1.0 + y)*(1.0 + z)/8.0;
      output(7, 1) =  (1.0 - x)*(1.0 + z)/8.0;
      output(7, 2) =  (1.0 - x)*(1.0 + y)/8.0;
      break;
    }
    case OPERATOR_D2 : {
      const auto x = input(0);
      const auto y = input(1);
      const auto z = input(2);

      // output is a rank-3 array with dimensions (basisCardin