C++ CellTopology::getDimension方法代码示例

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

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

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

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

  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 BilinearFormUtility::computeOptimalStiffnessMatrix(FieldContainer<double> &stiffness, 
                                                        FieldContainer<double> &optimalTestWeights,
                                                        BilinearFormPtr bilinearForm,
                                                        Teuchos::RCP<DofOrdering> trialOrdering, Teuchos::RCP<DofOrdering> testOrdering,
                                                        shards::CellTopology &cellTopo, FieldContainer<double> &physicalCellNodes,
                                                        FieldContainer<double> &cellSideParities) {
  // lots of code copied and pasted from the very similar computeStiffnessMatrix.  The difference here is that for each optimal test function,
  // we need to ask the bilinear form about each of its components (it's a vector whereas the other guy had just a single basis function
  // for each...), and then apply the appropriate weights....
  // physicalCellNodes: the nodal points for the element(s) with topology cellTopo
  //                 The dimensions are (numCells, numNodesPerElement, spaceDimension)
  // optimalTestWeights dimensions are: (numCells, numTrial, numTest) -- numTrial is the optTest index
  // stiffness dimensions are: (numCells, # trialOrdering Dofs, # trialOrdering Dofs)
  // (while (cell,trial,test) is more natural conceptually, I believe the above ordering makes
  //  more sense given the inversion that we must do to compute the optimal test functions...)
  
  // steps:
  // 0. Set up BasisCache
  // 3. For each (test, trial) combination:
  //   a. Apply the specified operators to the basis in the DofOrdering, at the cubature points
  //   b. Multiply the two bases together, weighted with Jacobian/Piola transform and cubature weights
  //   c. Pass the result to bilinearForm's applyBilinearFormData method
  //   d. Sum up (integrate) and place in stiffness matrix according to DofOrdering indices
  
  // 0. Set up Cubature
  
  unsigned numCells = physicalCellNodes.dimension(0);
  unsigned numNodesPerElem = physicalCellNodes.dimension(1);
  unsigned spaceDim = physicalCellNodes.dimension(2);
  
  // Check that cellTopo and physicalCellNodes agree
  TEUCHOS_TEST_FOR_EXCEPTION( ( numNodesPerElem != cellTopo.getNodeCount() ),
                     std::invalid_argument,
                     "Second dimension of physicalCellNodes and cellTopo.getNodeCount() do not match.");
  TEUCHOS_TEST_FOR_EXCEPTION( ( spaceDim != cellTopo.getDimension() ),
                     std::invalid_argument,
                     "Third dimension of physicalCellNodes and cellTopo.getDimension() do not match.");
  
  int numOptTestFunctions = optimalTestWeights.dimension(1); // should also == numTrialDofs
  
  TEUCHOS_TEST_FOR_EXCEPTION( ( optimalTestWeights.dimension(1) != stiffness.dimension(2) ),
                     std::invalid_argument,
                     "optimalTestWeights.dimension(1) (=" << optimalTestWeights.dimension(1) << ") and stiffness.dimension(2) (=" << stiffness.dimension(2) << ") do not match.");
  TEUCHOS_TEST_FOR_EXCEPTION( ( stiffness.dimension(1) != stiffness.dimension(2) ),
                     std::invalid_argument,
                     "stiffness.dimension(1) (=" << stiffness.dimension(1) << ") and stiffness.dimension(2) (=" << stiffness.dimension(2) << ") do not match.");
  
  // Set up BasisCache
  int cubDegreeTrial = trialOrdering->maxBasisDegree();
  int cubDegreeTest = testOrdering->maxBasisDegree();
  int cubDegree = cubDegreeTrial + cubDegreeTest;
  
  BasisCache basisCache(physicalCellNodes, cellTopo, *trialOrdering, cubDegreeTest, true); // DO create side caches, too
  
  unsigned numSides = CamelliaCellTools::getSideCount(cellTopo);

  vector<int> testIDs = bilinearForm->testIDs();
  vector<int>::iterator testIterator;
  
  vector<int> trialIDs = bilinearForm->trialIDs();
  vector<int>::iterator trialIterator;
  
  BasisPtr trialBasis,testBasis;
  
  stiffness.initialize(0.0);
  
  for (trialIterator = trialIDs.begin(); trialIterator != trialIDs.end(); trialIterator++) {
    int trialID = *trialIterator;
    
    for (int optTestIndex=0; optTestIndex < numOptTestFunctions; optTestIndex++) {
      FieldContainer<double> weights(numCells,testOrdering->totalDofs());
      for (unsigned i=0; i<numCells; i++) {
        for (int j=0; j<testOrdering->totalDofs(); j++) {
          weights(i,j) = optimalTestWeights(i,optTestIndex,j);
        }
      }
      for (testIterator = testIDs.begin(); testIterator != testIDs.end(); testIterator++) {
        int testID = *testIterator;
        
        vector<EOperatorExtended> trialOperators, testOperators;
        bilinearForm->trialTestOperators(trialID, testID, trialOperators, testOperators);
        vector<EOperatorExtended>::iterator trialOpIt, testOpIt;
        testOpIt = testOperators.begin();
        
        int operatorIndex = -1;
        for (trialOpIt = trialOperators.begin(); trialOpIt != trialOperators.end(); trialOpIt++) {
          IntrepidExtendedTypes::EOperatorExtended trialOperator = *trialOpIt;
          IntrepidExtendedTypes::EOperatorExtended testOperator = *testOpIt;
          operatorIndex++;
          
          if (testOperator==OP_TIMES_NORMAL) {
            TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,"OP_TIMES_NORMAL not supported for tests.  Use for trial only");
          }
          
          Teuchos::RCP < const FieldContainer<double> > testValuesTransformed;
          Teuchos::RCP < const FieldContainer<double> > trialValuesTransformed;
          Teuchos::RCP < const FieldContainer<double> > testValuesTransformedWeighted;

          if (! bilinearForm->isFluxOrTrace(trialID)) {
            trialBasis = trialOrdering->getBasis(trialID);
//.........这里部分代码省略.........

      void getLocalToGlobalMap(int (*local2global)[2],
                               int &off_global,
                               const Basis<Scalar,ArrayType> &basis,
                               const int *element) {
        const int local = 0, global = 1;
        const int nbf = basis.getCardinality();
        const shards::CellTopology cell = basis.getBaseCellTopology();
        const int dim = cell.getDimension();

        int cnt = 0, off_element = 0;
        int subcell_verts[4], nids;

        const int nvert = cell.getVertexCount();
        for (int i=0;i<nvert;++i) {
          const int ord_vert = (off_element < nbf ? basis.getDofOrdinal(0, i, 0) : 0);
          const int dof_vert = (off_element < nbf ? basis.getDofTag(ord_vert)[3] : 0);
      
          local2global[cnt][local] = off_element;
          off_element += dof_vert;
          Orientation::getElementNodeMap(subcell_verts, nids,
                                         cell, element,
                                         0, i);
      
          if (!findNode(local2global[cnt][global], subcell_verts, nids, true)) {
            addNode(subcell_verts, nids, off_global);
            local2global[cnt][global] = off_global;
            off_global += dof_vert;
          }
          ++cnt;
        }
        const int nedge = cell.getEdgeCount();
        for (int i=0;i<nedge;++i) {
          const int ord_edge = (off_element < nbf ? basis.getDofOrdinal(1, i, 0) : 0);
          const int dof_edge = (off_element < nbf ? basis.getDofTag(ord_edge)[3] : 0);
      
          local2global[cnt][local] = off_element;
          off_element += dof_edge;
          Orientation::getElementNodeMap(subcell_verts, nids,
                                         cell, element,
                                         1, i);
      
          if (!findNode(local2global[cnt][global], subcell_verts, nids, true)) {
            addNode(subcell_verts, nids, off_global);
            local2global[cnt][global] = off_global;
            off_global += dof_edge;
          }
          ++cnt;
        }
        const int nface = cell.getFaceCount();
        for (int i=0;i<nface;++i) {
          const int ord_face = (off_element < nbf ? basis.getDofOrdinal(2, i, 0) : 0);
          const int dof_face = (off_element < nbf ? basis.getDofTag(ord_face)[3] : 0);
      
          local2global[cnt][local] = off_element;
          off_element += dof_face;
          Orientation::getElementNodeMap(subcell_verts, nids,
                                         cell, element,
                                         2, i);
      
          if (!findNode(local2global[cnt][global], subcell_verts, nids, true)) {
            addNode(subcell_verts, nids, off_global);
            local2global[cnt][global] = off_global;
            off_global += dof_face;
          }
          ++cnt;
        }
        {
          const int i = 0;
          const int ord_intr = (off_element < nbf ? basis.getDofOrdinal(dim, i, 0) : 0);
          const int dof_intr = (off_element < nbf ? basis.getDofTag(ord_intr)[3]   : 0);
      
          local2global[cnt][local] = off_element;
          off_element += dof_intr;
          Orientation::getElementNodeMap(subcell_verts, nids,
                                         cell, element,
                                         dim, i);
      
          if (!findNode(local2global[cnt][global], subcell_verts, nids, true)) {
            addNode(subcell_verts, nids, off_global);
            local2global[cnt][global] = off_global;
            off_global += dof_intr;
          }
          ++cnt;
        }
    
        // add the last offset
        local2global[cnt][local] = off_element;
        local2global[cnt][global] = -1; // invalid values
      }

int NodalFieldPattern::getDimension() const
{
   const shards::CellTopology ct = getCellTopology();
   return ct.getDimension();
}

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 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
  CellTools<SpT>::
  setSubcellParametrization( subcellParamViewType      &subcellParam,
                             const ordinal_type         subcellDim,
                             const shards::CellTopology parentCell ) {
#ifdef HAVE_INTREPID2_DEBUG
    INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(parentCell), std::invalid_argument, 
                                  ">>> ERROR (Intrepid2::CellTools::setSubcellParametrization): the specified cell topology does not have a reference cell.");
#endif
    // subcellParametrization is rank-3 FieldContainer with dimensions (SC, PCD, COEF) where:
    //  - SC    is the subcell count of subcells with the specified dimension in the parent cell
    //  - PCD   is Parent Cell Dimension, which gives the number of coordinate functions in the map
    //          PCD = 2 for standard 2D cells and non-standard 2D cells: shell line and beam
    //          PCD = 3 for standard 3D cells and non-standard 3D cells: shell Tri and Quad
    //  - COEF  is number of coefficients needed to specify a coordinate function:
    //          COEFF = 2 for edge parametrizations
    //          COEFF = 3 for both Quad and Tri face parametrizations. Because all Quad reference faces
    //          are affine, the coefficient of the bilinear term u*v is zero and is not stored, i.e.,
    //          3 coefficients are sufficient to store Quad face parameterization maps.
    //  
    // Edge parametrization maps [-1,1] to edge defined by (v0, v1)
    // Face parametrization maps [-1,1]^2 to quadrilateral face (v0, v1, v2, v3), or
    // standard 2-simplex  {(0,0),(1,0),(0,1)} to traingle face (v0, v1, v2).
    // This defines orientation-preserving parametrizations with respect to reference edge and
    // face orientations induced by their vertex order. 

    // get subcellParametrization dimensions: (sc, pcd, coeff)
    const auto sc    = parentCell.getSubcellCount(subcellDim);
    const auto pcd   = parentCell.getDimension();   
    const auto coeff = (subcellDim == 1) ? 2 : 3;

    INTREPID2_TEST_FOR_EXCEPTION( subcellDim < 1 || subcellDim > static_cast<ordinal_type>(pcd-1), std::invalid_argument, 
                                  ">>> ERROR (Intrepid2::CellTools::setSubcellParametrization): Parametrizations defined in a range between 1 and (dim-1)");


    // create a view
    subcellParam = subcellParamViewType("CellTools::setSubcellParametrization",
                                        sc, pcd, coeff);

    referenceNodeDataViewType 
      v0("CellTools::setSubcellParametrization::v0", Parameters::MaxDimension), 
      v1("CellTools::setSubcellParametrization::v1", Parameters::MaxDimension),
      v2("CellTools::setSubcellParametrization::v1", Parameters::MaxDimension),
      v3("CellTools::setSubcellParametrization::v1", Parameters::MaxDimension);

    if (subcellDim == 1) {
      // Edge parametrizations of 2D and 3D cells (shell lines and beams are 2D cells with edges)
      for (size_type subcellOrd=0;subcellOrd<sc;++subcellOrd) {
        // vertexK[0] = x_k; vertexK[1] = y_k; vertexK[2] = z_k; z_k = 0 for 2D cells
        // Note that ShellLine and Beam are 2D cells!
        const auto v0ord = parentCell.getNodeMap(subcellDim, subcellOrd, 0);
        const auto v1ord = parentCell.getNodeMap(subcellDim, subcellOrd, 1);

        getReferenceVertex(v0, parentCell, v0ord);
        getReferenceVertex(v1, parentCell, v1ord);
        
        // x(t) = (x0 + x1)/2 + t*(x1 - x0)/2 
        subcellParam(subcellOrd, 0, 0) = (v0[0] + v1[0])/2.0;
        subcellParam(subcellOrd, 0, 1) = (v1[0] - v0[0])/2.0;
        
        // y(t) = (y0 + y1)/2 + t*(y1 - y0)/2 
        subcellParam(subcellOrd, 1, 0) = (v0[1] + v1[1])/2.0;
        subcellParam(subcellOrd, 1, 1) = (v1[1] - v0[1])/2.0;
        
        if( pcd == 3 ) {
          // z(t) = (z0 + z1)/2 + t*(z1 - z0)/2 
          subcellParam(subcellOrd, 2, 0) = (v0[2] + v1[2])/2.0;
          subcellParam(subcellOrd, 2, 1) = (v1[2] - v0[2])/2.0;
        }
      }
    }
    else if (subcellDim == 2) {
      // Face parametrizations of 3D cells: (shell Tri and Quad are 3D cells with faces)
      // A 3D cell can have both Tri and Quad faces, but because they are affine images of the
      // parametrization domain, 3 coefficients are enough to store them in both cases.
      for (size_type subcellOrd=0;subcellOrd<sc;++subcellOrd) {
        
        switch (parentCell.getKey(subcellDim,subcellOrd)) {
          
        case shards::Triangle<3>::key:
        case shards::Triangle<4>::key:
        case shards::Triangle<6>::key: {
          const auto v0ord = parentCell.getNodeMap(subcellDim, subcellOrd, 0);
          const auto v1ord = parentCell.getNodeMap(subcellDim, subcellOrd, 1);
          const auto v2ord = parentCell.getNodeMap(subcellDim, subcellOrd, 2);

          getReferenceVertex(v0, parentCell, v0ord);
          getReferenceVertex(v1, parentCell, v1ord);
          getReferenceVertex(v2, parentCell, v2ord);
          
          // x(u,v) = x0 + (x1 - x0)*u + (x2 - x0)*v
          subcellParam(subcellOrd, 0, 0) = v0[0];
          subcellParam(subcellOrd, 0, 1) = v1[0] - v0[0];
          subcellParam(subcellOrd, 0, 2) = v2[0] - v0[0];
          
          // y(u,v) = y0 + (y1 - y0)*u + (y2 - y0)*v
          subcellParam(subcellOrd, 1, 0) = v0[1];
          subcellParam(subcellOrd, 1, 1) = v1[1] - v0[1];
          subcellParam(subcellOrd, 1, 2) = v2[1] - v0[1];
          
//.........这里部分代码省略.........

  void
  CellTools<SpT>::
  getSubcellParametrization(       subcellParamViewType &subcellParam,  
                             const ordinal_type         subcellDim,
                             const shards::CellTopology parentCell ) {
#ifdef HAVE_INTREPID2_DEBUG
    INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(parentCell), std::invalid_argument, 
                                  ">>> ERROR (Intrepid2::CellTools::getSubcellParametrization): the specified cell topology does not have a reference cell.");
#endif

    if (!isSubcellParametrizationSet_)
      setSubcellParametrization();
    
    // Select subcell parametrization according to its parent cell type
    const auto pcd = parentCell.getDimension(); // parent cell dim
    INTREPID2_TEST_FOR_EXCEPTION( subcellDim < 1 || subcellDim > static_cast<ordinal_type>(pcd-1), std::invalid_argument, 
                                  ">>> ERROR (Intrepid2::CellTools::getSubcellParametrization): Parametrizations defined in a range between 1 and (dim-1)");
    
    switch (parentCell.getKey() ) {
    case shards::Tetrahedron<4>::key:
    case shards::Tetrahedron<8>::key:
    case shards::Tetrahedron<10>::key:
    case shards::Tetrahedron<11>::key:       subcellParam = ( subcellDim == 2 ? subcellParamData_.tetFaces : subcellParamData_.tetEdges ); break;
      
    case shards::Hexahedron<8>::key:
    case shards::Hexahedron<20>::key:
    case shards::Hexahedron<27>::key:        subcellParam = ( subcellDim == 2 ? subcellParamData_.hexFaces : subcellParamData_.hexEdges ); break;
      
    case shards::Pyramid<5>::key:
    case shards::Pyramid<13>::key:
    case shards::Pyramid<14>::key:           subcellParam = ( subcellDim == 2 ? subcellParamData_.pyrFaces : subcellParamData_.pyrEdges ); break;
      
    case shards::Wedge<6>::key:
    case shards::Wedge<15>::key:
    case shards::Wedge<18>::key:             subcellParam = ( subcellDim == 2 ? subcellParamData_.wedgeFaces : subcellParamData_.wedgeEdges ); break;

    case shards::Triangle<3>::key:
    case shards::Triangle<4>::key:
    case shards::Triangle<6>::key:           subcellParam = subcellParamData_.triEdges; break;
                  
    case shards::Quadrilateral<4>::key:
    case shards::Quadrilateral<8>::key:
    case shards::Quadrilateral<9>::key:      subcellParam = subcellParamData_.quadEdges; break;

    // case shards::ShellTriangle<3>::key:
    // case shards::ShellTriangle<6>::key:      subcellParam = ( subcellDim == 2 ? subcellParamData_.shellTriFaces : subcellParamData_.shellTriEdges ); break;

    // case shards::ShellQuadrilateral<4>::key:
    // case shards::ShellQuadrilateral<8>::key:
    // case shards::ShellQuadrilateral<9>::key: subcellParam = ( subcellDim == 2 ? subcellParamData_.shellQuadFaces : subcellParamData_.shellQuadEdges ); break;

    case shards::ShellLine<2>::key:
    case shards::ShellLine<3>::key:
    case shards::Beam<2>::key:
    case shards::Beam<3>::key:               subcellParam = subcellParamData_.lineEdges; break;
    default: {
      INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument, 
                                    ">>> ERROR (Intrepid2::CellTools::getSubcellParametrization): invalid cell topology.");
    }
    }
  }