52#ifndef Intrepid2_CellTopology_h
53#define Intrepid2_CellTopology_h
55#include <Shards_CellTopology.hpp>
64 using CellTopoPtr = Teuchos::RCP<CellTopology>;
65 using CellTopologyKey = std::pair<ordinal_type,ordinal_type>;
67 shards::CellTopology shardsBaseTopology_;
68 ordinal_type tensorialDegree_;
72 std::vector< std::vector<CellTopoPtr> > subcells_;
81 CellTopology(
const shards::CellTopology &baseTopo, ordinal_type tensorialDegree)
83 shardsBaseTopology_(baseTopo),
84 tensorialDegree_(tensorialDegree)
87 if (tensorialDegree_ == 0)
89 name_ = baseTopo.getName();
93 std::ostringstream nameStream;
94 nameStream << baseTopo.getName();
95 for (
int tensorialOrdinal = 0; tensorialOrdinal < tensorialDegree; tensorialOrdinal++)
97 nameStream <<
" x Line_2";
99 name_ = nameStream.str();
102 int baseDim = baseTopo.getDimension();
103 vector<ordinal_type> subcellCounts = vector<ordinal_type>(baseDim + tensorialDegree_ + 1);
104 subcells_ = vector< vector< CellTopoPtr > >(baseDim + tensorialDegree_ + 1);
106 if (tensorialDegree_==0)
108 for (
int d=0; d<=baseDim; d++)
110 subcellCounts[d] =
static_cast<ordinal_type
>(baseTopo.getSubcellCount(d));
116 subcellCounts[0] = 2 * tensorComponentTopo->getSubcellCount(0);
117 for (
int d=1; d < baseDim+tensorialDegree_; d++)
119 subcellCounts[d] = 2 * tensorComponentTopo->getSubcellCount(d) + tensorComponentTopo->getSubcellCount(d-1);
121 subcellCounts[baseDim + tensorialDegree_] = 1;
123 for (
int d=0; d<baseDim+tensorialDegree_; d++)
125 subcells_[d] = vector< CellTopoPtr >(subcellCounts[d]);
126 int subcellCount = subcells_[d].size();
127 for (
int scord=0; scord<subcellCount; scord++)
132 subcells_[baseDim+tensorialDegree_] = vector<CellTopoPtr>(1);
133 subcells_[baseDim+tensorialDegree_][0] = Teuchos::rcp(
this,
false);
139 return shardsBaseTopology_;
145 return tensorialDegree_;
151 return shardsBaseTopology_.getDimension() + tensorialDegree_;
154 static ordinal_type
getNodeCount(
const shards::CellTopology &shardsTopo)
156 if (shardsTopo.getDimension()==0)
return 1;
157 return shardsTopo.getNodeCount();
163 ordinal_type two_pow = 1 << tensorialDegree_;
195 int sideDim = spaceDim - 1;
202 std::pair<ordinal_type,ordinal_type>
getKey()
const
204 return std::make_pair(
static_cast<ordinal_type
>(shardsBaseTopology_.getKey()), tensorialDegree_);
212 const ordinal_type subcell_ord )
const
214 return subcells_[subcell_dim][subcell_ord]->getNodeCount();
229 const ordinal_type subcell_ord )
const
231 return subcells_[subcell_dim][subcell_ord]->getVertexCount();
239 const ordinal_type subcell_ord )
const
241 return subcells_[subcell_dim][subcell_ord]->getEdgeCount();
252 const ordinal_type subcell_ord_in_component_topo )
const
257 if (tensorialDegree_==0)
259 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"getExtrudedSubcellOrdinal() is not valid for un-extruded topologies");
263 ordinal_type componentSubcellCount =
getTensorialComponent()->getSubcellCount(subcell_dim_in_component_topo + 1);
264 return subcell_ord_in_component_topo + componentSubcellCount * 2;
273 const ordinal_type subcell_ord )
const
275 return subcells_[subcell_dim][subcell_ord]->getSideCount();
284 if (subcell_dim >= ordinal_type(subcells_.size()))
return 0;
285 else return subcells_[subcell_dim].size();
294 if (ordinal_type(tensorComponentNodes.size()) != tensorialDegree_ + 1)
296 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"tensorComponentNodes.size() != _tensorialDegree + 1");
308 ordinal_type node = 0;
309 CellTopoPtr line = CellTopology::line();
310 std::vector<CellTopoPtr> componentTopos(tensorialDegree_ + 1, line);
312 for (
int i=tensorComponentNodes.size()-1; i >= 0; i--)
314 ordinal_type componentNode = tensorComponentNodes[i];
315 node *= componentTopos[i]->getNodeCount();
316 node += componentNode;
328 const ordinal_type subcell_ord ,
329 const ordinal_type subcell_node_ord )
const
334 if (subcell_ord != 0)
336 TEUCHOS_TEST_FOR_EXCEPTION(subcell_ord != 0, std::invalid_argument,
"subcell ordinal out of bounds");
338 return subcell_node_ord;
340 else if (subcell_dim==0)
343 if (subcell_node_ord != 0)
345 TEUCHOS_TEST_FOR_EXCEPTION(subcell_node_ord != 0, std::invalid_argument,
"subcell node ordinal out of bounds");
349 if (tensorialDegree_==0)
351 return shardsBaseTopology_.getNodeMap(subcell_dim, subcell_ord, subcell_node_ord);
356 ordinal_type componentSubcellCount = tensorComponentTopo->getSubcellCount(subcell_dim);
357 if (subcell_ord < componentSubcellCount * 2)
359 ordinal_type subcell_ord_comp = subcell_ord % componentSubcellCount;
360 ordinal_type compOrdinal = subcell_ord / componentSubcellCount;
361 ordinal_type mappedNodeInsideComponentTopology = tensorComponentTopo->getNodeMap(subcell_dim, subcell_ord_comp, subcell_node_ord);
362 return mappedNodeInsideComponentTopology + compOrdinal * tensorComponentTopo->getNodeCount();
367 ordinal_type subcell_ord_comp = subcell_ord - componentSubcellCount * 2;
368 ordinal_type subcell_dim_comp = subcell_dim - 1;
369 CellTopoPtr subcellTensorComponent = tensorComponentTopo->getSubcell(subcell_dim_comp, subcell_ord_comp);
371 ordinal_type scCompOrdinal = subcell_node_ord / subcellTensorComponent->getNodeCount();
373 ordinal_type scCompNodeOrdinal = subcell_node_ord % subcellTensorComponent->getNodeCount();
374 ordinal_type mappedNodeInsideComponentTopology = tensorComponentTopo->getNodeMap(subcell_dim_comp, subcell_ord_comp, scCompNodeOrdinal);
375 return mappedNodeInsideComponentTopology + scCompOrdinal * tensorComponentTopo->getNodeCount();
388 const ordinal_type node_ord )
const;
395 const ordinal_type node_ord )
const;
399 CellTopoPtr
getSide( ordinal_type sideOrdinal )
const;
410 CellTopoPtr
getSubcell( ordinal_type scdim, ordinal_type scord )
const
412 if (tensorialDegree_==0)
414 return cellTopology(shardsBaseTopology_.getCellTopologyData(scdim, scord), 0);
419 ordinal_type componentSubcellCount = tensorComponentTopo->getSubcellCount(scdim);
420 if (scord < componentSubcellCount * 2)
422 scord = scord % componentSubcellCount;
423 return tensorComponentTopo->getSubcell(scdim, scord);
426 scord = scord - componentSubcellCount * 2;
428 CellTopoPtr subcellTensorComponent = tensorComponentTopo->getSubcell(scdim, scord);
429 return cellTopology(subcellTensorComponent->getBaseTopology(), subcellTensorComponent->getTensorialDegree() + 1);
437 static ordinal_type
getSubcellOrdinalMap(CellTopoPtr cellTopo, ordinal_type subcdim, ordinal_type subcord, ordinal_type subsubcdim, ordinal_type subsubcord)
443 typedef ordinal_type SubcellOrdinal;
444 typedef ordinal_type SubcellDimension;
445 typedef ordinal_type SubSubcellOrdinal;
446 typedef ordinal_type SubSubcellDimension;
447 typedef ordinal_type SubSubcellOrdinalInCellTopo;
448 typedef std::pair< SubcellDimension, SubcellOrdinal > SubcellIdentifier;
449 typedef std::pair< SubSubcellDimension, SubSubcellOrdinal > SubSubcellIdentifier;
450 typedef std::map< SubcellIdentifier, std::map< SubSubcellIdentifier, SubSubcellOrdinalInCellTopo > > OrdinalMap;
451 static std::map< CellTopologyKey, OrdinalMap > ordinalMaps;
453 if (subsubcdim==subcdim)
461 std::cout <<
"request for subsubcell of the same dimension as subcell, but with subsubcord > 0.\n";
462 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"request for subsubcell of the same dimension as subcell, but with subsubcord > 0.");
466 if (subcdim==cellTopo->getDimension())
474 std::cout <<
"request for subcell of the same dimension as cell, but with subsubcord > 0.\n";
475 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"request for subcell of the same dimension as cell, but with subsubcord > 0.");
479 CellTopologyKey key = cellTopo->getKey();
480 if (ordinalMaps.find(key) == ordinalMaps.end())
483 OrdinalMap ordinalMap;
484 ordinal_type sideDim = cellTopo->getDimension() - 1;
485 typedef ordinal_type NodeOrdinal;
486 std::map< std::set<NodeOrdinal>, SubcellIdentifier > subcellMap;
488 for (ordinal_type d=1; d<=sideDim; d++)
490 ordinal_type subcellCount = cellTopo->getSubcellCount(d);
491 for (ordinal_type subcellOrdinal=0; subcellOrdinal<subcellCount; subcellOrdinal++)
493 std::set<NodeOrdinal> nodes;
494 ordinal_type nodeCount = cellTopo->getNodeCount(d, subcellOrdinal);
495 for (NodeOrdinal subcNode=0; subcNode<nodeCount; subcNode++)
497 nodes.insert(cellTopo->getNodeMap(d, subcellOrdinal, subcNode));
499 SubcellIdentifier subcell = std::make_pair(d, subcellOrdinal);
500 subcellMap[nodes] = subcell;
502 CellTopoPtr subcellTopo = cellTopo->getSubcell(d, subcellOrdinal);
504 for (ordinal_type subsubcellDim=0; subsubcellDim<d; subsubcellDim++)
506 ordinal_type subsubcellCount = subcellTopo->getSubcellCount(subsubcellDim);
507 for (ordinal_type subsubcellOrdinal=0; subsubcellOrdinal<subsubcellCount; subsubcellOrdinal++)
509 SubSubcellIdentifier subsubcell = std::make_pair(subsubcellDim,subsubcellOrdinal);
510 if (subsubcellDim==0)
512 ordinalMap[subcell][subsubcell] = cellTopo->getNodeMap(subcell.first, subcell.second, subsubcellOrdinal);
515 ordinal_type nodeCount_inner = subcellTopo->getNodeCount(subsubcellDim, subsubcellOrdinal);
516 std::set<NodeOrdinal> subcellNodes;
517 for (NodeOrdinal subsubcNode=0; subsubcNode<nodeCount_inner; subsubcNode++)
519 NodeOrdinal subcNode = subcellTopo->getNodeMap(subsubcellDim, subsubcellOrdinal, subsubcNode);
520 NodeOrdinal node = cellTopo->getNodeMap(d, subcellOrdinal, subcNode);
521 subcellNodes.insert(node);
524 SubcellIdentifier subsubcellInCellTopo = subcellMap[subcellNodes];
525 ordinalMap[ subcell ][ subsubcell ] = subsubcellInCellTopo.second;
532 ordinalMaps[key] = ordinalMap;
534 SubcellIdentifier subcell = std::make_pair(subcdim, subcord);
535 SubSubcellIdentifier subsubcell = std::make_pair(subsubcdim, subsubcord);
536 if (ordinalMaps[key][subcell].find(subsubcell) != ordinalMaps[key][subcell].end())
538 return ordinalMaps[key][subcell][subsubcell];
542 std::cout <<
"For topology " << cellTopo->getName() <<
" and subcell " << subcord <<
" of dim " << subcdim;
543 std::cout <<
", subsubcell " << subsubcord <<
" of dim " << subsubcdim <<
" not found.\n";
544 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"subsubcell not found");
554 if (tensorialDegree_ > 0)
556 return cellTopology(shardsBaseTopology_, tensorialDegree_ - 1);
560 return Teuchos::null;
570 return extrusionNodeOrdinal;
579 if (tensorialDegree_ == 0)
return false;
581 return (sideOrdinal > 1) && (sideOrdinal < sideCount);
588 static CellTopoPtr
cellTopology(
const shards::CellTopology &shardsCellTopo, ordinal_type tensorialDegree = 0)
590 ordinal_type shardsKey =
static_cast<ordinal_type
>(shardsCellTopo.getBaseKey());
591 std::pair<ordinal_type,ordinal_type> key = std::make_pair(shardsKey, tensorialDegree);
593 static std::map< CellTopologyKey, CellTopoPtr > tensorizedShardsTopologies;
595 if (tensorizedShardsTopologies.find(key) == tensorizedShardsTopologies.end())
597 tensorizedShardsTopologies[key] = Teuchos::rcp(
new CellTopology(shardsCellTopo, tensorialDegree));
599 return tensorizedShardsTopologies[key];
602 static CellTopoPtr point()
604 return cellTopology(shards::getCellTopologyData<shards::Node >());
607 static CellTopoPtr line()
609 return cellTopology(shards::getCellTopologyData<shards::Line<> >());
612 static CellTopoPtr quad()
614 return cellTopology(shards::getCellTopologyData<shards::Quadrilateral<> >());
617 static CellTopoPtr hexahedron()
619 return cellTopology(shards::getCellTopologyData<shards::Hexahedron<> >());
622 static CellTopoPtr triangle()
624 return cellTopology(shards::getCellTopologyData<shards::Triangle<> >());
627 static CellTopoPtr tetrahedron()
629 return cellTopology(shards::getCellTopologyData<shards::Tetrahedron<> >());
632 static CellTopoPtr wedge()
634 return cellTopology(shards::getCellTopologyData<shards::Wedge<> >());
Implements arbitrary-dimensional extrusion of a base shards::CellTopology.
CellTopoPtr getSide(ordinal_type sideOrdinal) const
Returns a CellTopoPtr for the specified side.
ordinal_type getExtrudedSubcellOrdinal(const ordinal_type subcell_dim_in_component_topo, const ordinal_type subcell_ord_in_component_topo) const
Mapping from the tensorial component CellTopology's subcell ordinal to the corresponding subcell ordi...
ordinal_type getVertexCount() const
Vertex count of this cell topology.
ordinal_type getNodePermutation(const ordinal_type permutation_ord, const ordinal_type node_ord) const
Permutation of a cell's node ordinals.
static CellTopoPtr cellTopology(const shards::CellTopology &shardsCellTopo, ordinal_type tensorialDegree=0)
static accessor that returns a CellTopoPtr; these are lazily constructed and cached.
ordinal_type getTensorialDegree() const
The number of times we have taken a tensor product between a line topology and the shards topology to...
ordinal_type getNodeCount() const
Node count of this cell topology.
CellTopology(const shards::CellTopology &baseTopo, ordinal_type tensorialDegree)
Constructor.
std::string getName() const
Human-readable name of the CellTopology.
ordinal_type getSubcellCount(const ordinal_type subcell_dim) const
Subcell count of subcells of the given dimension.
ordinal_type getTensorialComponentSideOrdinal(ordinal_type extrusionNodeOrdinal)
Returns the side corresponding to the provided node in the final extrusion dimension.
ordinal_type getNodeMap(const ordinal_type subcell_dim, const ordinal_type subcell_ord, const ordinal_type subcell_node_ord) const
Mapping from a subcell's node ordinal to a node ordinal of this parent cell topology.
ordinal_type getFaceCount() const
Face (dimension 2) subcell count of this cell topology.
ordinal_type getVertexCount(const ordinal_type subcell_dim, const ordinal_type subcell_ord) const
Vertex count of a subcell of the given dimension and ordinal.
std::pair< ordinal_type, ordinal_type > getKey() const
Key that's unique for standard shards topologies and any tensorial degree.
ordinal_type getSideCount() const
Side (dimension N-1) subcell count of this cell topology.
ordinal_type getSideCount(const ordinal_type subcell_dim, const ordinal_type subcell_ord) const
Side count of a subcell of the given dimension and ordinal.
CellTopoPtr getSubcell(ordinal_type scdim, ordinal_type scord) const
Get the subcell of dimension scdim with ordinal scord.
bool sideIsExtrudedInFinalDimension(ordinal_type sideOrdinal) const
Returns true if the specified side has extension in the final tensorial dimension....
ordinal_type getNodePermutationCount() const
Number of node permutations defined for this cell.
const shards::CellTopology & getBaseTopology() const
Returns the underlying shards CellTopology.
static ordinal_type getSubcellOrdinalMap(CellTopoPtr cellTopo, ordinal_type subcdim, ordinal_type subcord, ordinal_type subsubcdim, ordinal_type subsubcord)
Maps the from a subcell within a subcell of the present CellTopology to the subcell in the present Ce...
ordinal_type getEdgeCount(const ordinal_type subcell_dim, const ordinal_type subcell_ord) const
Edge count of a subcell of the given dimension and ordinal.
ordinal_type getNodeCount(const ordinal_type subcell_dim, const ordinal_type subcell_ord) const
Node count of a subcell of the given dimension and ordinal.
ordinal_type getNodePermutationInverse(const ordinal_type permutation_ord, const ordinal_type node_ord) const
Inverse permutation of a cell's node ordinals.
ordinal_type getNodeFromTensorialComponentNodes(const std::vector< ordinal_type > &tensorComponentNodes) const
Mapping from the tensorial component node ordinals to the node ordinal of this tensor cell topology.
CellTopoPtr getTensorialComponent() const
For cell topologies of positive tensorial degree, returns the cell topology of tensorial degree one l...
ordinal_type getEdgeCount() const
Edge (dimension 1) subcell count of this cell topology.
ordinal_type getDimension() const
Dimension of this tensor topology.