C++ DynRankView::rank方法代码示例

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

  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
  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 
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_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 
  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." );
    }
    }
  }

  /** 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>::
  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");
    }
    }
  }

  void
  ArrayTools<ExecSpaceType>::
  dotMultiplyDataData( /**/  Kokkos::DynRankView<outputDataValueType,    outputDataProperties...>     outputData,
                       const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...>  inputDataLeft,
                       const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight ) {

#ifdef HAVE_INTREPID2_DEBUG
    {
      if (inputDataRight.rank() >= inputDataLeft.rank()) {
        INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() < 2 || inputDataLeft.rank() > 4, std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataData): Left data input container must have rank 2, 3 or 4.");
        INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() != inputDataLeft.rank(), std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataData): The rank of the right data input container must equal the rank of the left data input container.");
        INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 2, std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataData): Data output container must have rank 2.");
        INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.dimension(1) != inputDataRight.dimension(1) &&
                                        inputDataLeft.dimension(1) != 1, std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataData): First dimensions of the left and right data input containers (number of integration points) must agree or first left data dimension must be 1!");
        for (size_type i=0;i<inputDataLeft.rank();++i) {
          if (i != 1) {
            INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.dimension(i) != inputDataRight.dimension(i), std::invalid_argument,
                                            ">>> ERROR (ArrayTools::dotMultiplyDataData): inputDataLeft dimension (i) does not match to the dimension (i) of inputDataRight");
          }
        }
        for (size_type i=0;i<outputData.rank();++i) {
          INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.dimension(i) != outputData.dimension(i), std::invalid_argument,
                                          ">>> ERROR (ArrayTools::dotMultiplyDataData): inputDataRight dimension (i) does not match to the dimension (i) of outputData");
        }
      } else {
        INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() < 2 || inputDataLeft.rank() > 4, std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataData): Left data input container must have rank 2, 3 or 4.");
        INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() != (inputDataLeft.rank()-1), std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataData): Right data input container must have rank one less than the rank of left data input container.");
        INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 2, std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataData): Data output container must have rank 2.");
        INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.dimension(1) != inputDataRight.dimension(0) &&
                                        inputDataLeft.dimension(1) != 1, std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataData): Zeroth dimension of the right data input container and first dimension of left data input container (number of integration points) must agree or first left data dimension must be 1!");
        INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.dimension(0) != outputData.dimension(1), std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataData): Zeroth dimension of the right data input container and first dimension of output data container (number of integration points) must agree!");
        INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.dimension(0) != outputData.dimension(0), std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataData): Zeroth dimensions of the left data input and data output containers (number of integration domains) must agree!");
        for (size_type i=1;i<inputDataRight.rank();++i) {
          INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.dimension(i+1) != inputDataRight.dimension(i), std::invalid_argument,
                                          ">>> ERROR (ArrayTools::dotMultiplyDataData): inputDataLeft dimension (i+1) does not match to the dimension (i) of inputDataRight");
        }
      }
    }
#endif

    ArrayTools<ExecSpaceType>::Internal::dotMultiply( outputData,
                                                      inputDataLeft,
                                                      inputDataRight,
                                                      false );
  }

  void
  ArrayTools<ExecSpaceType>::
  dotMultiplyDataField( /**/  Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
                        const Kokkos::DynRankView<inputDataValueType,  inputDataProperties...>   inputData,
                        const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields ) {

#ifdef HAVE_INTREPID2_DEBUG
    {
      if (inputFields.rank() > inputData.rank()) {
        INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() < 2 || inputData.rank() > 4, std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataField): Input data container must have rank 2, 3 or 4.");
        INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() != (inputData.rank()+1), std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataField): Input fields container must have rank one larger than the rank of the input data container.");
        INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 3, std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataField): Output fields container must have rank 3.");
        INTREPID2_TEST_FOR_EXCEPTION( inputFields.dimension(0) != inputData.dimension(0), std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataField): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!");
        INTREPID2_TEST_FOR_EXCEPTION( inputData.dimension(1) != inputFields.dimension(2) &&
                                        inputData.dimension(1) != 1, std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataField): Second dimension of the fields input container and first dimension of data input container (number of integration points) must agree or first data dimension must be 1!");
        for (size_type i=2;i<inputData.rank();++i) {
          INTREPID2_TEST_FOR_EXCEPTION( inputData.dimension(i) != inputFields.dimension(i+1), std::invalid_argument,
                                          ">>> ERROR (ArrayTools::dotMultiplyDataField): inputData dimension (i) does not match to the dimension (i+1) of inputFields");
        }
        for (size_t i=0;i<outputFields.rank();++i) {
          INTREPID2_TEST_FOR_EXCEPTION( inputFields.dimension(i) != outputFields.dimension(i), std::invalid_argument,
                                          ">>> ERROR (ArrayTools::dotMultiplyDataField): inputFields dimension (i) does not match to the dimension (i+1) of outputFields");
        }
      } else {
        INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() < 2 || inputData.rank() > 4, std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataField): Input data container must have rank 2, 3 or 4.");
        INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() != inputData.rank(), std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataField): The rank of fields input container must equal the rank of data input container.");
        INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 3, std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataField): Output fields container must have rank 3.");
        INTREPID2_TEST_FOR_EXCEPTION( inputData.dimension(1) != inputFields.dimension(1) &&
                                        inputData.dimension(1) != 1, std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataField): First dimensions of the fields and data input containers (number of integration points) must agree or first data dimension must be 1!");
        INTREPID2_TEST_FOR_EXCEPTION( inputFields.dimension(0) != outputFields.dimension(1), std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataField): Zeroth dimension of the fields input container and first dimension of the fields output container (number of fields) must agree!");
        INTREPID2_TEST_FOR_EXCEPTION( inputFields.dimension(1) != outputFields.dimension(2), std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataField): First dimension of the fields input container and second dimension of the fields output container (number of integration points) must agree!");
        INTREPID2_TEST_FOR_EXCEPTION( outputFields.dimension(0) != inputData.dimension(0), std::invalid_argument,
                                        ">>> ERROR (ArrayTools::dotMultiplyDataField): Zeroth dimensions of the fields output and data input containers (number of integration domains) must agree!");
        for (size_type i=2;i<inputData.rank();++i) {
          INTREPID2_TEST_FOR_EXCEPTION( inputData.dimension(i) != inputFields.dimension(i), std::invalid_argument,
                                          ">>> ERROR (ArrayTools::dotMultiplyDataField): inputData dimension (i) does not match to the dimension (i) of inputFields");
        }
      }
    }
#endif

    ArrayTools<ExecSpaceType>::Internal::dotMultiply( outputFields,
                                                      inputData,
                                                      inputFields,
                                                      true );
  }

  void
  ArrayTools<ExecSpaceType>::
  scalarMultiplyDataData( /**/  Kokkos::DynRankView<outputDataValueType,    outputDataProperties...>     outputData,
                          const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...>  inputDataLeft,
                          const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight,
                          const bool reciprocal ) {

#ifdef HAVE_INTREPID2_DEBUG
    {
      INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() != 2, std::invalid_argument,
                                ">>> ERROR (ArrayTools::scalarMultiplyDataData): Left input data container must have rank 2.");

      if (outputData.rank() <= inputDataRight.rank()) {
        INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 2 || inputDataRight.rank() > 4, std::invalid_argument,
                                  ">>> ERROR (ArrayTools::scalarMultiplyDataData): Right input data container must have rank 2, 3, or 4.");
        INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != inputDataRight.rank(), std::invalid_argument,
                                  ">>> ERROR (ArrayTools::scalarMultiplyDataData): Right input and output data containers must have the same rank.");
        INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.dimension(0) != inputDataLeft.dimension(0), std::invalid_argument,
                                  ">>> ERROR (ArrayTools::scalarMultiplyDataData): Zeroth dimensions (number of integration domains) of the left and right data input containers must agree!");
        INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.dimension(1) != inputDataRight.dimension(1) &&
                                  inputDataLeft.dimension(1) != 1, std::invalid_argument,
                                  ">>> ERROR (ArrayTools::scalarMultiplyDataData): First dimensions of the left and right data input containers (number of integration points) must agree, or first dimension of the left data input container must be 1!");
        for (size_type i=0;i<inputDataRight.rank();++i) {
          INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.dimension(i) != outputData.dimension(i), std::invalid_argument,
                                    ">>> ERROR (ArrayTools::scalarMultiplyDataData): inputDataRight dimension (i) does not match to the dimension (i) of outputData");
        }
      } else {
        INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 1 || inputDataRight.rank() > 3, std::invalid_argument,
                                  ">>> ERROR (ArrayTools::scalarMultiplyDataData): Right input data container must have rank 1, 2, or 3.");
        INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != (inputDataRight.rank()+1), std::invalid_argument,
                                  ">>> ERROR (ArrayTools::scalarMultiplyDataData): The rank of the right input data container must be one less than the rank of the output data container.");
        INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.dimension(1) != inputDataRight.dimension(0)  &&
                                  inputDataLeft.dimension(1) != 1, std::invalid_argument,
                                  ">>> ERROR (ArrayTools::scalarMultiplyDataData): Zeroth dimension of the right input data container and first dimension of the left data input container (number of integration points) must agree or first dimension of the left data input container must be 1!");
        INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.dimension(0) != outputData.dimension(0), std::invalid_argument,
                                  ">>> ERROR (ArrayTools::scalarMultiplyDataData): Zeroth dimensions of data output and left data input containers (number of integration domains) must agree!");
        for (size_type i=0;i<inputDataRight.rank();++i) {
          INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.dimension(i) != outputData.dimension(i+1), std::invalid_argument,
                                    ">>> ERROR (ArrayTools::scalarMultiplyDataData): inputDataRight dimension (i) does not match to the dimension (i+1) of outputData");
        }
      }
    }
#endif
    // body
    ArrayTools<ExecSpaceType>::Internal::scalarMultiply( outputData,
                                                         inputDataLeft,
                                                         inputDataRight,
                                                         false,
                                                         reciprocal );
  }