int 
 Intrepid2FieldPattern::
 getSubcellCount(int dim) const
 {
   const shards::CellTopology ct = intrepidBasis_->getBaseCellTopology();
   return ct.getSubcellCount(dim);
 }

void Intrepid2FieldPattern::buildSubcellClosure(const shards::CellTopology & cellTopo,unsigned dim,unsigned subCell,
                                               std::set<std::pair<unsigned,unsigned> > & closure)
{
   switch(dim) {
   case 0:
      closure.insert(std::make_pair(0,subCell));
      break;
   case 1:
      closure.insert(std::make_pair(0,cellTopo.getNodeMap(dim,subCell,0)));
      closure.insert(std::make_pair(0,cellTopo.getNodeMap(dim,subCell,1)));
      closure.insert(std::make_pair(1,subCell));
      break;
   case 2:
      {
      unsigned cnt = (shards::CellTopology(cellTopo.getCellTopologyData(dim,subCell))).getSubcellCount(dim-1);
      for(unsigned i=0;i<cnt;i++) {
         int edge = mapCellFaceEdge(cellTopo.getCellTopologyData(),subCell,i);
         buildSubcellClosure(cellTopo,dim-1,edge,closure);
      }
      closure.insert(std::make_pair(2,subCell));
      }
      break;
   default:
      // beyond a two dimension surface this thing crashes!
      TEUCHOS_ASSERT(false);
   };
}

    inline
    void 
    OrientationTools::
    mapToModifiedReference(outPointViewType outPoints,
                           const refPointViewType refPoints,
                           const shards::CellTopology cellTopo,
                           const ordinal_type cellOrt) {
#ifdef HAVE_INTREPID2_DEBUG
      {
        const auto cellDim = cellTopo.getDimension();
        INTREPID2_TEST_FOR_EXCEPTION( !( (1 <= cellDim) && (cellDim <= 2 ) ), std::invalid_argument,
                                      ">>> ERROR (Intrepid::OrientationTools::mapToModifiedReference): " \
                                      "Method defined only for 1 and 2-dimensional subcells.");
        
        INTREPID2_TEST_FOR_EXCEPTION( !( outPoints.dimension(0) == refPoints.dimension(0) ), std::invalid_argument,
                                      ">>> ERROR (Intrepid::OrientationTools::mapToModifiedReference): " \
                                      "Size of input and output point arrays does not match each other.");
      }
#endif
      
      // Apply the parametrization map to every point in parameter domain
      const auto numPts = outPoints.dimension(0);
      const auto key = cellTopo.getBaseCellTopologyData()->key;
      switch (key) {
      case shards::Line<>::key : {
        for (auto pt=0;pt<numPts;++pt)
          getModifiedLinePoint(outPoints(pt, 0),
                               refPoints(pt, 0),
                               cellOrt);
        break;
      }
      case shards::Triangle<>::key : {
        for (auto pt=0;pt<numPts;++pt)
          getModifiedTrianglePoint(outPoints(pt, 0), outPoints(pt, 1),
                                   refPoints(pt, 0), refPoints(pt, 1),
                                   cellOrt);
        break;
      }
      case shards::Quadrilateral<>::key : {
        for (auto pt=0;pt<numPts;++pt)
          getModifiedQuadrilateralPoint(outPoints(pt, 0), outPoints(pt, 1),
                                        refPoints(pt, 0), refPoints(pt, 1),
                                        cellOrt);
        break;
      }
      default: {
        INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument,
                                      ">>> ERROR (Intrepid2::OrientationTools::mapToModifiedReference): " \
                                      "Invalid cell topology." );
        break;
      }
      }
    }

  inline
  ordinal_type
  PointTools::
  getLatticeSize( const shards::CellTopology cellType,
                  const ordinal_type order,
                  const ordinal_type offset ) {
#ifdef HAVE_INTREPID2_DEBUG    
    INTREPID2_TEST_FOR_EXCEPTION( order < 0 || offset < 0,
                                  std::invalid_argument ,
                                  ">>> ERROR (PointTools::getLatticeSize): order and offset must be positive values." );
#endif
    ordinal_type r_val = 0;
    switch (cellType.getBaseKey()) {
    case shards::Tetrahedron<>::key: {
      const auto effectiveOrder = order - 4 * offset;
      r_val = (effectiveOrder < 0 ? 0 :(effectiveOrder+1)*(effectiveOrder+2)*(effectiveOrder+3)/6);
      break;
    }
    case shards::Triangle<>::key: {
      const auto effectiveOrder = order - 3 * offset;
      r_val = (effectiveOrder < 0 ? 0 : (effectiveOrder+1)*(effectiveOrder+2)/2);
      break;
    }
    case shards::Line<>::key: {
      const auto effectiveOrder = order - 2 * offset;
      r_val = (effectiveOrder < 0 ? 0 : (effectiveOrder+1));
      break;
    }
    default: {
      INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument ,
                                    ">>> ERROR (Intrepid2::PointTools::getLatticeSize): the specified cell type is not supported." );
    }
    }
    return r_val;
  }

void Intrepid2FieldPattern::getSubcellNodes(const shards::CellTopology & cellTopo,unsigned dim,unsigned subCell,
                                           std::vector<unsigned> & nodes)
{
   if(dim==0) {
      nodes.push_back(subCell);
      return;
   }

   // get all nodes on requested sub cell
   unsigned subCellNodeCount = cellTopo.getNodeCount(dim,subCell);
   for(unsigned node=0;node<subCellNodeCount;++node)
      nodes.push_back(cellTopo.getNodeMap(dim,subCell,node));

   // sort them so they are ordered correctly for "includes" call
   std::sort(nodes.begin(),nodes.end());
}

void Intrepid2FieldPattern::findContainedSubcells(const shards::CellTopology & cellTopo,unsigned dim,
                                                 const std::vector<unsigned> & nodes,
                                                 std::set<std::pair<unsigned,unsigned> > & subCells)
{

   unsigned subCellCount = cellTopo.getSubcellCount(dim); 
   for(unsigned subCellOrd=0;subCellOrd<subCellCount;++subCellOrd) {
      // get all nodes in sub cell
      std::vector<unsigned> subCellNodes;
      getSubcellNodes(cellTopo,dim,subCellOrd,subCellNodes);

      // if subCellNodes \subset nodes => add (dim,subCellOrd) to subCells
      bool isSubset = std::includes(       nodes.begin(),        nodes.end(),
                                    subCellNodes.begin(), subCellNodes.end());
      if(isSubset)
         subCells.insert(std::make_pair(dim,subCellOrd));
       
   }

   // stop recursion base case
   if(dim==0) return;

   // find subcells in next sub dimension
   findContainedSubcells(cellTopo,dim-1,nodes,subCells);
}

IntrepidSideCell<MDArray>::IntrepidSideCell( 
    const shards::CellTopology& side_topology,
    const unsigned side_id,
    const shards::CellTopology& parent_topology,
    const unsigned degree )
    : Base( side_topology, degree )
    , d_side_id( side_id )
    , d_parent_topology( parent_topology )
{
    // Map the side cubature points to the cell frame.
    MDArray mapped_cub_points( this->d_cubature->getNumPoints(), 
			       parent_topology.getDimension() );
    Intrepid::CellTools<Scalar>::mapToReferenceSubcell(
	mapped_cub_points, this->d_cub_points, side_topology.getDimension(),
	d_side_id, d_parent_topology );
    this->d_cub_points = mapped_cub_points;
}

 Basis_HGRAD_POLY_C1_FEM<Scalar, ArrayScalar>::Basis_HGRAD_POLY_C1_FEM(const shards::CellTopology& cellTopology){
   this -> basisCardinality_  = cellTopology.getNodeCount();
   this -> basisDegree_       = 1;
   this -> basisCellTopology_ = cellTopology;
   this -> basisType_         = BASIS_FEM_DEFAULT;
   this -> basisCoordinates_  = COORDINATES_CARTESIAN;
   this -> basisTagsAreSet_   = false;
 }

   unsigned
   getLocalSubcellIndexFromGlobalNodeList(const ArrayCellGIDs& cellGIDs,
                                          const ArraySideGIDs& subcellGIDs,
                                          const shards::CellTopology& cell,unsigned subcell_dim)
 {
   unsigned local_subcell;
   bool found_local_subcell = false;
   unsigned subcell = 0;
   while ( (subcell < cell.getSubcellCount(subcell_dim)) && (!found_local_subcell) ) {
 
     unsigned num_subcell_nodes =
 	cell.getCellTopologyData()->subcell[subcell_dim][subcell].topology->node_count;
 
     std::list<unsigned> tmp_subcell_gid_list;
     for (unsigned node = 0; node < num_subcell_nodes; ++node)
       tmp_subcell_gid_list.push_back(cellGIDs[cell.getNodeMap(subcell_dim,
                                                            subcell, node)]);
 
     bool subcell_matches = true;
     unsigned node = 0;
     while ( subcell_matches && (node < num_subcell_nodes) ) {
 
       std::list<unsigned>::iterator search =
         std::find(tmp_subcell_gid_list.begin(), tmp_subcell_gid_list.end(),
                   subcellGIDs[node]);
 
       if (search == tmp_subcell_gid_list.end())
         subcell_matches = false;
 
       ++node;
     }
 
     if (subcell_matches) {
       found_local_subcell = true;
       local_subcell = subcell;
     }
 
     ++subcell;
   }
 
   TEUCHOS_TEST_FOR_EXCEPTION(!found_local_subcell, std::runtime_error,
                      "Failed to find subcell!");
 
   return local_subcell;
 }

 void getBoundaryFlags(int *boundary,
                       const shards::CellTopology cell,
                       const int *element) {
   int subcell_verts[4], nids;
   const int dim = cell.getDimension();
   const int nside = cell.getSideCount();
   for (int i=0;i<nside;++i) {
     Orientation::getElementNodeMap(subcell_verts, nids,
                                    cell, element,
                                    dim-1, i);
     
     if (!findBoundary(boundary[i], subcell_verts, nids)) {
       TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error,
                                   ">>> ERROR (Intrepid::HGRAD_TRI_Cn::Test 04): " \
                                   "Side node is not found");
     }
   }
 }

 void PointTools::
 getWarpBlendLattice( /**/  Kokkos::DynRankView<pointValueType,pointProperties...> points,
                      const shards::CellTopology cell,
                      const ordinal_type order,
                      const ordinal_type offset ) {
   switch (cell.getBaseKey()) {
   // case shards::Tetrahedron<>::key: getWarpBlendLatticeTetrahedron( points, order, offset );  break;
   // case shards::Triangle<>::key:     getWarpBlendLatticeTriangle   ( points, order, offset );  break;
   case shards::Line<>::key:         getWarpBlendLatticeLine       ( points, order, offset );  break;
   default: {
     INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument ,
                                   ">>> ERROR (Intrepid2::PointTools::getWarpBlendLattice): the specified cell type is not supported." );
   }
   }
 }

    unsigned 
    getLocalSideIndexFromGlobalNodeList(const ArrayCellGIDs& cellGIDs, 
					const ArraySideGIDs& sideGIDs,
					const shards::CellTopology& cell)
  {
    unsigned cell_dim = cell.getDimension();
    //TEUCHOS_TEST_FOR_EXCEPTION(!cell.getSubcellHomogeneity(cell_dim - 1),
    //	       std::runtime_error, "Sides are not homogeneous!");
    
    unsigned local_side;
    bool found_local_side = false;
    unsigned side = 0;  
    while ( (side < cell.getSideCount()) && (!found_local_side) ) {
      
      const shards::CellTopology 
	side_topo(cell.getCellTopologyData(cell.getDimension()-1, side));
      
      unsigned num_side_nodes = 
	cell.getCellTopologyData()->side[side].topology->node_count;
 

      std::list<unsigned> tmp_side_gid_list;
      for (unsigned node = 0; node < num_side_nodes; ++node)
	tmp_side_gid_list.push_back(cellGIDs[cell.getNodeMap(cell_dim - 1, 
							     side, node)]);
     
      bool side_matches = true;
      unsigned node = 0;
      while ( side_matches && (node < num_side_nodes) ) {

	std::list<unsigned>::iterator search = 
	  std::find(tmp_side_gid_list.begin(), tmp_side_gid_list.end(),
		    sideGIDs[node]);
	
	if (search == tmp_side_gid_list.end())
	  side_matches = false;
	  
	++node;
      }
      
      if (side_matches) {
	found_local_side = true;
	local_side = side;
      }
      
      ++side;
    }
    
    TEUCHOS_TEST_FOR_EXCEPTION(!found_local_side, std::runtime_error,
		       "Failed to find side!");
    
    return local_side;
  }

  Basis_Constant_FEM<SpT,OT,PT>::
  Basis_Constant_FEM(const shards::CellTopology cellTopo) {
    const ordinal_type spaceDim = cellTopo.getDimension();

    this->basisCardinality_  = 1;
    this->basisDegree_       = 0;
    this->basisCellTopology_ = cellTopo;
    this->basisType_         = Intrepid2::BASIS_FEM_DEFAULT;
    this->basisCoordinates_  = Intrepid2::COORDINATES_CARTESIAN;

    // initialize tags
    {
      // Basis-dependent intializations
      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

      // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
      ordinal_type tags[4] = { spaceDim, 0, 0, 1 };

      ordinal_type_array_1d_host tagView(&tags[0], 4);

      this->setOrdinalTagData(this->tagToOrdinal_,
                              this->ordinalToTag_,
                              tagView,
                              this->basisCardinality_,
                              tagSize,
                              posScDim,
                              posScOrd,
                              posDfOrd);
    }

    // dofCoords on host and create its mirror view to device
    Kokkos::DynRankView<typename scalarViewType::value_type,typename SpT::array_layout,Kokkos::HostSpace>
      dofCoords("dofCoordsHost", this->basisCardinality_, spaceDim), cellVerts("cellVerts", spaceDim);

    CellTools<SpT>::getReferenceCellCenter(Kokkos::subview(dofCoords, 0, Kokkos::ALL()),
                                           cellVerts,
                                           cellTopo);

    this->dofCoords_ = Kokkos::create_mirror_view(typename SpT::memory_space(), dofCoords);
    Kokkos::deep_copy(this->dofCoords_, dofCoords);
  }

void mapToPhysicalFrame(ArrayPhysPoint      &        physPoints,
                                           const ArrayRefPoint &        refPoints,
                                           const ArrayCell     &        cellWorkset,
                                           const shards::CellTopology & cellTopo,
                                           const int &                  whichCell)
{
  int spaceDim  = (int)cellTopo.getDimension();
  int numCells  = cellWorkset.dimension(0);
  //points can be rank-2 (P,D), or rank-3 (C,P,D)
  int numPoints = (refPoints.rank() == 2) ? refPoints.dimension(0) : refPoints.dimension(1);

  // Initialize physPoints
  for(int i = 0; i < physPoints.dimentions(0); i++)
    for(int j = 0; j < physPoints.dimentions(1); j++)  
      for(int k = 0; k < physPoints.dimentions(2); k++)
         physPoints(i,j,k) = 0.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." );
    }
    }
  }

  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 mapToPhysicalFrame(ArrayPhysPoint      &        physPoints,
                                           const ArrayRefPoint &        refPoints,
                                           const ArrayCell     &        cellWorkset,
                                           const shards::CellTopology & cellTopo,
                                           const int &                  whichCell)
{
  int spaceDim  = (int)cellTopo.getDimension();
  int numCells  = cellWorkset.dimension(0);
  //points can be rank-2 (P,D), or rank-3 (C,P,D)
  int numPoints = (refPoints.rank() == 2) ? refPoints.dimension(0) : refPoints.dimension(1);

/*  // Mapping is computed using an appropriate H(grad) basis function: define RCP to the base class
  Teuchos::RCP<Basis<Scalar, FieldContainer<Scalar> > > HGRAD_Basis;

  // Choose the H(grad) basis depending on the cell topology. \todo define maps for shells and beams
  switch( cellTopo.getKey() ){

    // Standard Base topologies (number of cellWorkset = number of vertices)
      case shards::Line<2>::key:
      HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_LINE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
      break;

    case shards::Triangle<3>::key:
      HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C1_FEM<Scalar, FieldContainer<Scalar> >() );
      break;

    case shards::Quadrilateral<4>::key:
      HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C1_FEM<Scalar, FieldContainer<Scalar> >() );
      break;

    case shards::Tetrahedron<4>::key:
      HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C1_FEM<Scalar, FieldContainer<Scalar> >() );
      break;

    case shards::Hexahedron<8>::key:
      HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C1_FEM<Scalar, FieldContainer<Scalar> >() );
      break;

    case shards::Wedge<6>::key:
      HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
      break;

    case shards::Pyramid<5>::key:
          HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_C1_FEM<Scalar, FieldContainer<Scalar> >() );
          break;

    // Standard Extended topologies
    case shards::Triangle<6>::key:
      HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C2_FEM<Scalar, FieldContainer<Scalar> >() );
      break;

    case shards::Quadrilateral<9>::key:
      HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C2_FEM<Scalar, FieldContainer<Scalar> >() );
      break;

    case shards::Tetrahedron<10>::key:
      HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C2_FEM<Scalar, FieldContainer<Scalar> >() );
      break;

    case shards::Tetrahedron<11>::key:
      HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_COMP12_FEM<Scalar, FieldContainer<Scalar> >() );
      break;

    case shards::Hexahedron<27>::key:
      HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C2_FEM<Scalar, FieldContainer<Scalar> >() );
      break;

    case shards::Wedge<18>::key:
      HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C2_FEM<Scalar, FieldContainer<Scalar> >() );
      break;

    // These extended topologies are not used for mapping purposes
    case shards::Quadrilateral<8>::key:
    case shards::Hexahedron<20>::key:
    case shards::Wedge<15>::key:
      TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
                          ">>> ERROR (Intrepid::CellTools::mapToPhysicalFrame): Cell topology not supported. ");
      break;

    // Base and Extended Line, Beam and Shell topologies 
     case shards::Line<3>::key:
    case shards::Beam<2>::key:
    case shards::Beam<3>::key:
    case shards::ShellLine<2>::key:
    case shards::ShellLine<3>::key:
    case shards::ShellTriangle<3>::key:
    case shards::ShellTriangle<6>::key:
    case shards::ShellQuadrilateral<4>::key:
    case shards::ShellQuadrilateral<8>::key:
    case shards::ShellQuadrilateral<9>::key:
      TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
                          ">>> ERROR (Intrepid::CellTools::mapToPhysicalFrame): Cell topology not supported. ");
      break;
    default:
      TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
                          ">>> ERROR (Intrepid::CellTools::mapToPhysicalFrame): Cell topology not supported.");
  }// switch  

  // Temp (F,P) array for the values of nodal basis functions at the reference points
  int basisCardinality = HGRAD_Basis -> getCardinality();
//.........这里部分代码省略.........

void testSubcellParametrizations(int&                               errorFlag,
                                 const shards::CellTopology&        parentCell,
                                 const FieldContainer<double>&      subcParamVert_A,
                                 const FieldContainer<double>&      subcParamVert_B,
                                 const int                          subcDim,
                                 const Teuchos::RCP<std::ostream>&  outStream){
  
  // Get cell dimension and subcell count
  int cellDim      = parentCell.getDimension();
  int subcCount    = parentCell.getSubcellCount(subcDim);
  
  
  // Loop over subcells of the specified dimension
  for(int subcOrd = 0; subcOrd < subcCount; subcOrd++){
    int subcVertexCount = parentCell.getVertexCount(subcDim, subcOrd);
    
    
    // Storage for correct reference subcell vertices and for the images of the parametrization domain points
    FieldContainer<double> refSubcellVertices(subcVertexCount, cellDim);
    FieldContainer<double> mappedParamVertices(subcVertexCount, cellDim);
    
    
    // Retrieve correct reference subcell vertices
    CellTools<double>::getReferenceSubcellVertices(refSubcellVertices, subcDim, subcOrd, parentCell);
    
    
    // Map vertices of the parametrization domain to 1 or 2-subcell with ordinal subcOrd
    // For edges parametrization domain is always 1-cube passed as "subcParamVert_A"
    if(subcDim == 1) {
      CellTools<double>::mapToReferenceSubcell(mappedParamVertices,
                                               subcParamVert_A,
                                               subcDim,
                                               subcOrd,
                                               parentCell);
    }
    // For faces need to treat Triangle and Quadrilateral faces separately
    else if(subcDim == 2) {
      
      // domain "subcParamVert_A" is the standard 2-simplex  
      if(subcVertexCount == 3){
        CellTools<double>::mapToReferenceSubcell(mappedParamVertices,
                                                 subcParamVert_A,
                                                 subcDim,
                                                 subcOrd,
                                                 parentCell);
      }
      // Domain "subcParamVert_B" is the standard 2-cube
      else if(subcVertexCount == 4){
        CellTools<double>::mapToReferenceSubcell(mappedParamVertices,
                                                 subcParamVert_B,
                                                 subcDim,
                                                 subcOrd,
                                                 parentCell);
      }
    }
    
    // Compare the images of the parametrization domain vertices with the true vertices.
    for(int subcVertOrd = 0; subcVertOrd < subcVertexCount; subcVertOrd++){
      for(int dim = 0; dim <  cellDim; dim++){
        
        if(mappedParamVertices(subcVertOrd, dim) != refSubcellVertices(subcVertOrd, dim) ) {
          errorFlag++; 
          *outStream 
            << std::setw(70) << "^^^^----FAILURE!" << "\n"
            << " Cell Topology = " << parentCell.getName() << "\n"
            << " Parametrization of subcell " << subcOrd << " which is "
            << parentCell.getName(subcDim,subcOrd) << " failed for vertex " << subcVertOrd << ":\n"
            << " parametrization map fails to map correctly coordinate " << dim << " of that vertex\n\n";
          
        }//if
      }// for dim 
    }// for subcVertOrd      
  }// for subcOrd
  
}

  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;
                  }
                }
//.........这里部分代码省略.........