C++ shards::CellTopology类代码示例

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