50#ifndef Intrepid2_CellGeometryDef_h
51#define Intrepid2_CellGeometryDef_h
60 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
63 using BasisPtr = Teuchos::RCP<Intrepid2::Basis<DeviceType,PointScalar,PointScalar> >;
67 static std::map<const CellGeometryType *, shards::CellTopology> cellTopology_;
68 static std::map<const CellGeometryType *, BasisPtr> basisForNodes_;
71 static void constructorCalled(
const CellGeometryType *cellGeometry,
const shards::CellTopology &cellTopo, BasisPtr basisForNodes)
73 cellTopology_[cellGeometry] = cellTopo;
74 basisForNodes_[cellGeometry] = basisForNodes;
79 cellTopology_.erase(cellGeometry);
80 basisForNodes_.erase(cellGeometry);
85 return basisForNodes_[cellGeometry];
88 static const shards::CellTopology & getCellTopology(
const CellGeometryType *cellGeometry)
90 return cellTopology_[cellGeometry];
101 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
104 Kokkos::View<PointScalar**, DeviceType> cellMeasures_;
105 Kokkos::View<PointScalar**, DeviceType> detData_;
111 cellMeasures_(cellMeasures),
113 cubatureWeights_(cubatureWeights)
116 KOKKOS_INLINE_FUNCTION
void
117 operator () (
const ordinal_type cellOrdinal,
const ordinal_type pointOrdinal)
const
119 cellMeasures_(cellOrdinal,pointOrdinal) = detData_(cellOrdinal,pointOrdinal) * cubatureWeights_(pointOrdinal);
124 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
125 KOKKOS_INLINE_FUNCTION
128 nodeOrdering_(cellGeometry.nodeOrdering_),
129 cellGeometryType_(cellGeometry.cellGeometryType_),
130 subdivisionStrategy_(cellGeometry.subdivisionStrategy_),
131 affine_(cellGeometry.affine_),
132 orientations_(cellGeometry.orientations_),
133 origin_(cellGeometry.origin_),
134 domainExtents_(cellGeometry.domainExtents_),
135 gridCellCounts_(cellGeometry.gridCellCounts_),
136 tensorVertices_(cellGeometry.tensorVertices_),
137 cellToNodes_(cellGeometry.cellToNodes_),
138 nodes_(cellGeometry.nodes_),
139 numCells_(cellGeometry.numCells_),
140 numNodesPerCell_(cellGeometry.numNodesPerCell_)
143#ifndef INTREPID2_COMPILE_DEVICE_CODE
144 shards::CellTopology cellTopo = cellGeometry.
cellTopology();
147 HostMemberLookup::constructorCalled(
this, cellTopo,
basisForNodes);
151 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
152 KOKKOS_INLINE_FUNCTION
156#ifndef INTREPID2_COMPILE_DEVICE_CODE
158 HostMemberLookup::destructorCalled(
this);
162 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
163 KOKKOS_INLINE_FUNCTION
166 switch (subdivisionStrategy) {
169 case TWO_TRIANGLES_LEFT:
170 case TWO_TRIANGLES_RIGHT:
174 case FIVE_TETRAHEDRA:
182 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
186 ScalarView<PointScalar,DeviceType> data;
188 const int CELL_DIM = 0;
189 const int POINT_DIM = 1;
190 const int D1_DIM = 2;
191 const int D2_DIM = 3;
193 const int numCellsWorkset = (endCell == -1) ? (numCells_ - startCell) : (endCell - startCell);
195 Kokkos::Array<int,7> extents { numCellsWorkset, pointsPerCell, spaceDim, spaceDim, 1, 1, 1 };
198 int blockPlusDiagonalLastNonDiagonal = -1;
200 if (cellGeometryType_ == UNIFORM_GRID)
202 if (uniformJacobianModulus() != 1)
204 variationType[CELL_DIM] =
MODULAR;
205 variationType[POINT_DIM] =
CONSTANT;
206 variationType[D1_DIM] =
GENERAL;
207 variationType[D2_DIM] =
GENERAL;
209 int cellTypeModulus = uniformJacobianModulus();
211 data =
getMatchingViewWithLabel(pointComponentView,
"CellGeometryProvider: Jacobian data", cellTypeModulus, spaceDim, spaceDim);
218 blockPlusDiagonalLastNonDiagonal = -1;
223 else if (cellGeometryType_ == TENSOR_GRID)
227 else if (cellGeometryType_ == FIRST_ORDER)
229 const bool simplex = (spaceDim + 1 == cellToNodes_.extent_int(1));
232 variationType[CELL_DIM] =
GENERAL;
233 variationType[POINT_DIM] =
CONSTANT;
234 variationType[D1_DIM] =
GENERAL;
235 variationType[D2_DIM] =
GENERAL;
241 variationType[CELL_DIM] =
GENERAL;
242 variationType[D1_DIM] =
GENERAL;
243 variationType[D2_DIM] =
GENERAL;
247 variationType[POINT_DIM] =
CONSTANT;
252 variationType[POINT_DIM] =
GENERAL;
253 data =
getMatchingViewWithLabel(data,
"CellGeometryProvider: Jacobian data", numCellsWorkset, pointsPerCell, spaceDim, spaceDim);
257 else if (cellGeometryType_ == HIGHER_ORDER)
260 variationType[CELL_DIM] =
GENERAL;
261 variationType[POINT_DIM] =
GENERAL;
262 variationType[D1_DIM] =
GENERAL;
263 variationType[D2_DIM] =
GENERAL;
264 data =
getMatchingViewWithLabel(data,
"CellGeometryProvider: Jacobian data", numCellsWorkset, pointsPerCell, spaceDim, spaceDim);
275 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
278 const int startCell,
const int endCell)
const
280 const int numCellsWorkset = (endCell == -1) ? (numCells_ - startCell) : (endCell - startCell);
282 if (cellGeometryType_ == UNIFORM_GRID)
284 if (uniformJacobianModulus() != 1)
286 int cellTypeModulus = uniformJacobianModulus();
289 auto dataHost = Kokkos::create_mirror_view(dataView3);
291 const int startCellType = startCell % cellTypeModulus;
292 const int endCellType = (numCellsWorkset >= cellTypeModulus) ? startCellType + cellTypeModulus : startCellType + numCellsWorkset;
293 const int gridCellOrdinal = 0;
294 for (
int cellType=startCellType; cellType<endCellType; cellType++)
296 const int subdivisionOrdinal = cellType % cellTypeModulus;
297 const int nodeZero = 0;
299 for (
int i=0; i<spaceDim; i++)
301 for (
int j=0; j<spaceDim; j++)
303 const int node = j+1;
305 const auto J_ij = subdivisionCoordinate(gridCellOrdinal, subdivisionOrdinal, node, i) - subdivisionCoordinate(gridCellOrdinal, subdivisionOrdinal, nodeZero, i);
306 dataHost(cellType,i,j) = J_ij;
311 Kokkos::deep_copy(dataView3,dataHost);
317 const auto domainExtents = domainExtents_;
318 const auto gridCellCounts = gridCellCounts_;
320 using ExecutionSpace =
typename DeviceType::execution_space;
321 auto policy = Kokkos::RangePolicy<>(ExecutionSpace(),0,spaceDim);
322 Kokkos::parallel_for(
"fill jacobian", policy, KOKKOS_LAMBDA(
const int d1)
325 const double REF_SPACE_EXTENT = 2.0;
326 dataView1(d1) = (domainExtents[d1] / REF_SPACE_EXTENT) / gridCellCounts[d1];
328 ExecutionSpace().fence();
331 else if (cellGeometryType_ == TENSOR_GRID)
335 else if ((cellGeometryType_ == FIRST_ORDER) || (cellGeometryType_ == HIGHER_ORDER))
337 const bool simplex = (spaceDim + 1 == cellToNodes_.extent_int(1));
343 auto cellToNodes = cellToNodes_;
346 using ExecutionSpace =
typename DeviceType::execution_space;
347 auto policy = Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<3>>({startCell,0,0},{numCellsWorkset,spaceDim,spaceDim});
349 Kokkos::parallel_for(
"compute first-order simplex Jacobians", policy,
350 KOKKOS_LAMBDA (
const int &cellOrdinal,
const int &d1,
const int &d2) {
351 const int nodeZero = 0;
352 const int node = d2+1;
353 const auto & nodeCoord = nodes(cellToNodes(cellOrdinal,node), d1);
354 const auto & nodeZeroCoord = nodes(cellToNodes(cellOrdinal,nodeZero), d1);
355 const PointScalar J_ij = nodeCoord - nodeZeroCoord;
356 dataView3(cellOrdinal,d1,d2) = (spaceDim != 1) ? J_ij : J_ij * 0.5;
362 auto basisForNodes = this->basisForNodes();
370 const int onePoint = 1;
372 auto tempData =
getMatchingViewWithLabel(dataView3,
"CellGeometryProvider: temporary Jacobian data", numCellsWorkset, onePoint, spaceDim, spaceDim);
374 Kokkos::deep_copy(testPointView, 0.0);
378 auto tempDataSubview = Kokkos::subview(tempData, Kokkos::ALL(), 0, Kokkos::ALL(), Kokkos::ALL());
379 Kokkos::deep_copy(dataView3, tempDataSubview);
384 TEUCHOS_TEST_FOR_EXCEPTION(basisForNodes == Teuchos::null, std::invalid_argument,
"basisForNodes must not be null");
385 TEUCHOS_TEST_FOR_EXCEPTION(dataView.size() == 0, std::invalid_argument,
"underlying view is not valid");
390 if (refData.
rank() == 3)
416 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
418 const Kokkos::Array<PointScalar,spaceDim> &domainExtents,
419 const Kokkos::Array<int,spaceDim> &gridCellCounts,
423 nodeOrdering_(nodeOrdering),
424 cellGeometryType_(UNIFORM_GRID),
425 subdivisionStrategy_(subdivisionStrategy),
428 domainExtents_(domainExtents),
429 gridCellCounts_(gridCellCounts)
432 for (
int d=0; d<spaceDim; d++)
434 numCells_ *= gridCellCounts_[d];
438 shards::CellTopology cellTopo;
442 numNodesPerCell_ = 1 << spaceDim;
446 cellTopo = shards::CellTopology(shards::getCellTopologyData<shards::Line<> >());
448 else if (spaceDim == 2)
450 cellTopo = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<> >());
452 else if (spaceDim == 3)
454 cellTopo = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<> >());
464 numNodesPerCell_ = spaceDim + 1;
467 cellTopo = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<> >());
469 else if (spaceDim == 3)
471 cellTopo = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<> >());
480 const int linearPolyOrder = 1;
481 BasisPtr
basisForNodes = getBasis<BasisFamily>(cellTopo, FUNCTION_SPACE_HGRAD, linearPolyOrder);
487 if (cellTopo.getKey() == shards::Quadrilateral<>::key)
491 else if (cellTopo.getKey() == shards::Hexahedron<>::key)
498 HostMemberLookup::constructorCalled(
this, cellTopo,
basisForNodes);
504 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
506 ScalarView<int,DeviceType> cellToNodes,
507 ScalarView<PointScalar,DeviceType> nodes,
508 const bool claimAffine,
511 nodeOrdering_(nodeOrdering),
512 cellGeometryType_(FIRST_ORDER),
513 cellToNodes_(cellToNodes),
516 if(cellToNodes.is_allocated())
518 numCells_ = cellToNodes.extent_int(0);
519 numNodesPerCell_ = cellToNodes.extent_int(1);
524 numCells_ = nodes.extent_int(0);
525 numNodesPerCell_ = nodes.extent_int(1);
533 const bool simplicialTopo = (cellTopo.getNodeCount() == cellTopo.getDimension() + 1);
534 affine_ = simplicialTopo;
542 const int linearPolyOrder = 1;
543 BasisPtr
basisForNodes = getBasis<BasisFamily>(cellTopo, FUNCTION_SPACE_HGRAD, linearPolyOrder);
549 if (cellTopo.getKey() == shards::Quadrilateral<>::key)
553 else if (cellTopo.getKey() == shards::Hexahedron<>::key)
560 HostMemberLookup::constructorCalled(
this, cellTopo,
basisForNodes);
564 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
566 ScalarView<PointScalar,DeviceType> cellNodes)
568 nodeOrdering_(HYPERCUBE_NODE_ORDER_TENSOR),
569 cellGeometryType_(HIGHER_ORDER),
572 numCells_ = cellNodes.extent_int(0);
573 numNodesPerCell_ = cellNodes.extent_int(1);
576 const bool firstOrderGeometry = (
basisForNodes->getDegree() == 1);
579 shards::CellTopology cellTopo =
basisForNodes->getBaseCellTopology();
581 if (firstOrderGeometry && (cellTopo.getNodeCount() == spaceDim + 1))
590 HostMemberLookup::constructorCalled(
this, cellTopo,
basisForNodes);
593 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
594 KOKKOS_INLINE_FUNCTION
600 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
609 INTREPID2_TEST_FOR_EXCEPTION(cubatureWeights.
rank() != 1, std::invalid_argument,
"cubatureWeights container must have shape (P)");
612 std::vector< Data<PointScalar,DeviceType> > tensorComponents(numTensorComponents);
616 const int cellExtent = jacobianDet.
extent_int(0);
619 Kokkos::Array<int,7> cellExtents{cellExtent,1,1,1,1,1,1};
621 ScalarView<PointScalar,DeviceType> detDataView (
"cell relative volumes", cellDataDim);
624 for (
int cubTensorComponent=0; cubTensorComponent<numTensorComponents-1; cubTensorComponent++)
627 const auto cubatureExtents = cubatureComponent.
getExtents();
628 const auto cubatureVariationTypes = cubatureComponent.getVariationTypes();
629 const int numPoints = cubatureComponent.getDataExtent(0);
630 ScalarView<PointScalar,DeviceType> cubatureWeightView (
"cubature component weights", numPoints);
631 const int pointComponentRank = 1;
632 tensorComponents[cubTensorComponent+1] =
Data<PointScalar,DeviceType>(cubatureWeightView,pointComponentRank,cubatureExtents,cubatureVariationTypes);
637 const int cellExtent = jacobianDet.
extent_int(0);
641 const int numPoints = cubatureWeights.
extent_int(0);
642 Kokkos::Array<int,7> extents{cellExtent,numPoints,1,1,1,1,1};
644 ScalarView<PointScalar,DeviceType> cubatureWeightView;
647 cubatureWeightView = ScalarView<PointScalar,DeviceType>(
"cell measure", cellDataDim, numPoints);
651 cubatureWeightView = ScalarView<PointScalar,DeviceType>(
"cell measure", numPoints);
653 const int cellMeasureRank = 2;
656 const bool separateFirstComponent = (numTensorComponents > 1);
660 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
670 "cellMeasure must either have a tensor component count of 1 or a tensor component count that is one higher than that of cubatureWeights");
672 INTREPID2_TEST_FOR_EXCEPTION(cubatureWeights.
rank() != 1, std::invalid_argument,
"cubatureWeights container must have shape (P)");
679 for (
int i=1; i<numTensorDimensions+1; i++)
688 const bool detCellVaries = detVaries[0] !=
CONSTANT;
689 const bool detPointVaries = detVaries[1] !=
CONSTANT;
691 if (detCellVaries && detPointVaries)
695 const int numCells = detData.extent_int(0);
696 const int numPoints = detData.extent_int(1);
697 INTREPID2_TEST_FOR_EXCEPTION(numCells != cellMeasureData.extent_int(0), std::invalid_argument,
"cellMeasureData doesn't match jacobianDet in cell dimension");
698 INTREPID2_TEST_FOR_EXCEPTION(numPoints != cellMeasureData.extent_int(1), std::invalid_argument,
"cellMeasureData doesn't match jacobianDet in point dimension");
703 using ExecutionSpace =
typename DeviceType::execution_space;
704 Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<2>> rangePolicy({0,0},{numCells,numPoints});
705 Kokkos::parallel_for(rangePolicy, cellMeasureFunctor);
707 else if (detCellVaries && !detPointVaries)
711 using ExecutionSpace =
typename DeviceType::execution_space;
712 Kokkos::parallel_for(
713 Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<2>>({0,0},{detData.extent_int(0),cubatureWeights.
extent_int(0)}),
714 KOKKOS_LAMBDA (
int cellOrdinal,
int pointOrdinal) {
715 cellMeasureData(cellOrdinal,pointOrdinal) = detData(cellOrdinal) * cubatureWeights(pointOrdinal);
724 using ExecutionSpace =
typename DeviceType::execution_space;
725 Kokkos::parallel_for(Kokkos::RangePolicy<ExecutionSpace>(0,cellMeasureData.extent_int(0)),
726 KOKKOS_LAMBDA (
const int &pointOrdinal) {
727 cellMeasureData(pointOrdinal) = detData(0) * cubatureWeights(pointOrdinal);
733 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
734 typename CellGeometry<PointScalar,spaceDim,DeviceType>::BasisPtr
738 return HostMemberLookup::getBasis(
this);
741 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
745 return HostMemberLookup::getCellTopology(
this);
748 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
749 KOKKOS_INLINE_FUNCTION
752 if (cellGeometryType_ == UNIFORM_GRID)
754 const int numSubdivisions = numCellsPerGridCell(subdivisionStrategy_);
755 if (numSubdivisions == 1)
767 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
771 if (cellGeometryType_ == UNIFORM_GRID)
776 else if (cellGeometryType_ == TENSOR_GRID)
781 else if ((cellGeometryType_ == FIRST_ORDER) || (cellGeometryType_ == HIGHER_ORDER))
783 const bool simplex = (spaceDim + 1 == cellToNodes_.extent_int(1));
791 auto basisForNodes = this->basisForNodes();
802 if (points.rank() == 2)
805 return getJacobianRefData(tensorPoints);
809 const int numCells = points.extent_int(0);
810 const int numPoints = points.extent_int(1);
811 const int numFields = basisForNodes->getCardinality();
813 auto cellBasisGradientsView =
getMatchingViewWithLabel(points,
"CellGeometryProvider: cellBasisGradients", numCells, numFields, numPoints, spaceDim);
814 auto basisGradientsView =
getMatchingViewWithLabel(points,
"CellGeometryProvider: basisGradients", numFields, numPoints, spaceDim);
816 for (
int cellOrdinal=0; cellOrdinal<numCells; cellOrdinal++)
818 auto refPointsForCell = Kokkos::subview(points, cellOrdinal, Kokkos::ALL(), Kokkos::ALL());
819 basisForNodes->getValues(basisGradientsView, refPointsForCell, OPERATOR_GRAD);
826 using ExecutionSpace =
typename DeviceType::execution_space;
827 auto policy = Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<3>>({0,0,0},{numFields,numPoints,spaceDim});
829 Kokkos::parallel_for(
"copy basis gradients", policy,
830 KOKKOS_LAMBDA (
const int &fieldOrdinal,
const int &pointOrdinal,
const int &d) {
831 cellBasisGradientsView(cellOrdinal,fieldOrdinal,pointOrdinal,d) = basisGradientsView(fieldOrdinal,pointOrdinal,d);
833 ExecutionSpace().fence();
851 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
855 if (cellGeometryType_ == UNIFORM_GRID)
860 else if (cellGeometryType_ == TENSOR_GRID)
865 else if ((cellGeometryType_ == FIRST_ORDER) || (cellGeometryType_ == HIGHER_ORDER))
867 const bool simplex = (spaceDim + 1 == cellToNodes_.extent_int(1));
875 auto basisForNodes = this->basisForNodes();
884 auto basisGradients = basisForNodes->allocateBasisValues(points, OPERATOR_GRAD);
885 basisForNodes->getValues(basisGradients, points, OPERATOR_GRAD);
888 int numFields = basisForNodes->getCardinality();
896 auto basisGradientsView =
getMatchingViewWithLabel(firstPointComponentView,
"CellGeometryProvider: temporary basisGradients", numFields, numPoints, spaceDim);
898 using ExecutionSpace =
typename DeviceType::execution_space;
899 auto policy = Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<3>>({0,0,0},{numFields,numPoints,spaceDim});
901 Kokkos::parallel_for(
"copy basis gradients", policy,
902 KOKKOS_LAMBDA (
const int &fieldOrdinal,
const int &pointOrdinal,
const int &d) {
903 basisGradientsView(fieldOrdinal,pointOrdinal,d) = basisGradients(fieldOrdinal,pointOrdinal,d);
905 ExecutionSpace().fence();
920 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
921 KOKKOS_INLINE_FUNCTION
924 if (nodeOrdering_ == HYPERCUBE_NODE_ORDER_CLASSIC_SHARDS)
930 if ((hypercubeNodeNumber % 4 == 1) || (hypercubeNodeNumber % 4 == 2))
937 if ((hypercubeNodeNumber % 4 == 2) || (hypercubeNodeNumber % 4 == 3))
944 const int nodesForPriorDimensions = 1 << d;
945 if ((hypercubeNodeNumber / nodesForPriorDimensions) % 2 == 1)
951 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
954 using HostExecSpace = Kokkos::DefaultHostExecutionSpace;
956 const bool isGridType = (cellGeometryType_ == TENSOR_GRID) || (cellGeometryType_ == UNIFORM_GRID);
957 const int numOrientations = isGridType ? numCellsPerGridCell(subdivisionStrategy_) : numCells();
959 const int nodesPerCell = numNodesPerCell();
961 ScalarView<Orientation, DeviceType> orientationsView(
"orientations", numOrientations);
962 auto orientationsHost = Kokkos::create_mirror_view(
typename HostExecSpace::memory_space(), orientationsView);
964 ScalarView<PointScalar, HostExecSpace> cellNodesHost(
"cellNodesHost",numOrientations,nodesPerCell);
972 const int numSubdivisions = numCellsPerGridCell(subdivisionStrategy_);
974#if defined(INTREPID2_COMPILE_DEVICE_CODE)
977 const int gridCellOrdinal = 0;
978 auto hostPolicy = Kokkos::MDRangePolicy<HostExecSpace,Kokkos::Rank<2>>({0,0},{numSubdivisions,nodesPerCell});
979 Kokkos::parallel_for(
"fill cellNodesHost", hostPolicy,
980 [
this,gridCellOrdinal,cellNodesHost] (
const int &subdivisionOrdinal,
const int &nodeInCell) {
981 auto node = this->gridCellNodeForSubdivisionNode(gridCellOrdinal, subdivisionOrdinal, nodeInCell);
982 cellNodesHost(subdivisionOrdinal,nodeInCell) = node;
990 auto cellToNodesHost = Kokkos::create_mirror_view_and_copy(
typename HostExecSpace::memory_space(), cellToNodes_);
994 Kokkos::deep_copy(orientationsView,orientationsHost);
996 const int orientationsRank = 1;
997 const Kokkos::Array<int,7> orientationExtents {
static_cast<int>(numCells_),1,1,1,1,1,1};
1002 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1003 KOKKOS_INLINE_FUNCTION
1011 return numNodesPerCell_;
1023 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1024 template <
typename iType>
1025 KOKKOS_INLINE_FUNCTION
1026 typename std::enable_if<std::is_integral<iType>::value,
int>::type
1029 return static_cast<int>(extent(r));
1032 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1033 KOKKOS_INLINE_FUNCTION
1037 return nodeOrdering_;
1040 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1041 KOKKOS_INLINE_FUNCTION
1047 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1048 KOKKOS_INLINE_FUNCTION
1051 if (cellGeometryType_ == UNIFORM_GRID)
1053 return gridCellCounts_[dim];
1055 else if (cellGeometryType_ == TENSOR_GRID)
1057 return tensorVertices_.extent_int(dim);
1065 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1066 KOKKOS_INLINE_FUNCTION
1069 return numNodesPerCell_;
1072 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1073 KOKKOS_INLINE_FUNCTION
1077 return orientations_(cellNumber);
1080 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1083 if (!orientations_.isValid())
1085 initializeOrientations();
1087 return orientations_;
1090 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1091 KOKKOS_INLINE_FUNCTION
1094 const int componentNode = hypercubeComponentNodeNumber(localNodeNumber, dim);
1095 int cellCountForPriorDimensions = 1;
1096 for (
int d=0; d<dim; d++)
1098 cellCountForPriorDimensions *= numCellsInDimension(d);
1100 const int componentGridCellOrdinal = (gridCellOrdinal / cellCountForPriorDimensions) % numCellsInDimension(dim);
1101 const int vertexOrdinal = componentGridCellOrdinal + componentNode;
1102 if (cellGeometryType_ == UNIFORM_GRID)
1104 return origin_[dim] + (vertexOrdinal * domainExtents_[dim]) / gridCellCounts_[dim];
1106 else if (cellGeometryType_ == TENSOR_GRID)
1108 Kokkos::Array<int,spaceDim> pointOrdinalComponents;
1109 for (
int d=0; d<spaceDim; d++)
1111 pointOrdinalComponents[d] = 0;
1113 pointOrdinalComponents[dim] = vertexOrdinal;
1114 return tensorVertices_(pointOrdinalComponents,dim);
1123 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1124 KOKKOS_INLINE_FUNCTION
1130 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1131 KOKKOS_INLINE_FUNCTION
1133 const int &subdivisionNodeNumber)
const
1136 switch (subdivisionStrategy_)
1138 case NO_SUBDIVISION:
1139 return subdivisionNodeNumber;
1140 case TWO_TRIANGLES_RIGHT:
1141 case TWO_TRIANGLES_LEFT:
1142 case FOUR_TRIANGLES:
1144 Kokkos::Array<int,3> nodeLookup;
1145 if (subdivisionStrategy_ == TWO_TRIANGLES_RIGHT)
1147 if (subdivisionOrdinal == 0)
1150 nodeLookup = {0,1,2};
1152 else if (subdivisionOrdinal == 1)
1157 nodeLookup = {2,3,0};
1164 else if (subdivisionStrategy_ == TWO_TRIANGLES_LEFT)
1166 if (subdivisionOrdinal == 0)
1172 nodeLookup = {3,0,1};
1174 else if (subdivisionOrdinal == 1)
1179 nodeLookup = {2,3,0};
1195 if (subdivisionNodeNumber == 1)
1202 nodeLookup = {(subdivisionOrdinal + 1) % 4, -1, subdivisionOrdinal};
1205 const int gridCellNodeNumber = nodeLookup[subdivisionNodeNumber];
1206 return gridCellNodeNumber;
1208 case FIVE_TETRAHEDRA:
1209 case SIX_TETRAHEDRA:
1211 Kokkos::Array<int,4> nodeLookup;
1212 if (subdivisionStrategy_ == FIVE_TETRAHEDRA)
1242 switch (subdivisionOrdinal) {
1244 nodeLookup = {1,3,4,6};
1247 nodeLookup = {0,1,3,4};
1250 nodeLookup = {1,2,3,6};
1253 nodeLookup = {1,4,5,6};
1256 nodeLookup = {3,4,6,7};
1263 else if (subdivisionStrategy_ == SIX_TETRAHEDRA)
1267 const int gridCellNodeNumber = nodeLookup[subdivisionNodeNumber];
1268 return gridCellNodeNumber;
1277 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1278 KOKKOS_INLINE_FUNCTION
1280 const int &subdivisionNodeNumber,
const int &d)
const
1282 int gridCellNode = gridCellNodeForSubdivisionNode(gridCellOrdinal, subdivisionOrdinal, subdivisionNodeNumber);
1284 if (subdivisionStrategy_ == FOUR_TRIANGLES)
1287 if (gridCellNode == 4)
1291 const int gridVertex0 = 0;
1292 const int gridVertex1 = (d == 0) ? 1 : 3;
1293 return 0.5 * (gridCellCoordinate(gridCellOrdinal, gridVertex0, d) + gridCellCoordinate(gridCellOrdinal, gridVertex1, d));
1296 return gridCellCoordinate(gridCellOrdinal, gridCellNode, d);
1299 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1300 KOKKOS_INLINE_FUNCTION
1303 if ((cellGeometryType_ == UNIFORM_GRID) || (cellGeometryType_ == TENSOR_GRID))
1305 const int numSubdivisions = numCellsPerGridCell(subdivisionStrategy_);
1306 if (numSubdivisions == 1)
1309 return gridCellCoordinate(cell, node, dim);
1313 const int subdivisionOrdinal = cell % numSubdivisions;
1314 const int gridCellOrdinal = cell / numSubdivisions;
1315 return subdivisionCoordinate(gridCellOrdinal, subdivisionOrdinal, node, dim);
1320#ifdef HAVE_INTREPID2_DEBUG
1328 if (cellToNodes_.is_allocated())
1330 const int nodeNumber = cellToNodes_(cell,node);
1331 return nodes_(nodeNumber,dim);
1335 return nodes_(cell,node,dim);
1340 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1341 KOKKOS_INLINE_FUNCTION
1344 if (cellGeometryType_ == UNIFORM_GRID)
1346 return numCellsPerGridCell(subdivisionStrategy_);
1354 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1357 const int pointsPerCell = points.
extent_int(0);
1358 return allocateJacobianDataPrivate(points.
getTensorComponent(0),pointsPerCell,startCell,endCell);
1361 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1365 const int pointDimension = (points.rank() == 3) ? 1 : 0;
1366 const int pointsPerCell = points.extent_int(pointDimension);
1367 return allocateJacobianDataPrivate(points,pointsPerCell,startCell,endCell);
1370 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1375 ScalarView<PointScalar,DeviceType> emptyPoints;
1376 return allocateJacobianDataPrivate(emptyPoints,numPoints,startCell,endCell);
1379 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1383 const int pointsPerCell = points.
extent_int(0);
1384 setJacobianDataPrivate(jacobianData,pointsPerCell,refData,startCell,endCell);
1387 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1392 const int pointDimension = (points.rank() == 3) ? 1 : 0;
1393 const int pointsPerCell = points.extent_int(pointDimension);
1394 setJacobianDataPrivate(jacobianData,pointsPerCell,refData,startCell,endCell);
1397 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1403 setJacobianDataPrivate(jacobianData,numPoints,emptyRefData,startCell,endCell);
@ GENERAL
arbitrary variation
@ BLOCK_PLUS_DIAGONAL
one of two dimensions in a matrix; bottom-right part of matrix is diagonal
@ MODULAR
varies according to modulus of the index
#define INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(test, x, msg)
Kokkos::DynRankView< typename ViewType::value_type, typename DeduceLayout< ViewType >::result_layout, typename ViewType::device_type > getMatchingViewWithLabel(const ViewType &view, const std::string &label, DimArgs... dims)
Creates and returns a view that matches the provided view in Kokkos Layout.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Hexahedron cell.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Quadrilateral cell.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
CellGeometry provides the nodes for a set of cells; has options that support efficient definition of ...
void computeCellMeasure(TensorData< PointScalar, DeviceType > &cellMeasure, const Data< PointScalar, DeviceType > &jacobianDet, const TensorData< PointScalar, DeviceType > &cubatureWeights) const
Compute cell measures that correspond to provided Jacobian determinants and.
void setJacobianDataPrivate(Data< PointScalar, DeviceType > &jacobianData, const int &pointsPerCell, const Data< PointScalar, DeviceType > &refData, const int startCell, const int endCell) const
Notionally-private method that provides a common interface for multiple public-facing setJacobianData...
BasisPtr basisForNodes() const
H^1 Basis used in the reference-to-physical transformation. Linear for straight-edged geometry; highe...
KOKKOS_INLINE_FUNCTION int numCellsPerGridCell(SubdivisionStrategy subdivisionStrategy) const
Helper method that returns the number of cells into which each grid cell will be subdivided based on ...
KOKKOS_INLINE_FUNCTION size_t extent(const int &r) const
Returns the logical extent of the container in the specified dimension; the shape of CellGeometry is ...
void setJacobian(Data< PointScalar, DeviceType > &jacobianData, const TensorPoints< PointScalar, DeviceType > &points, const Data< PointScalar, DeviceType > &refData, const int startCell=0, const int endCell=-1) const
Compute Jacobian values for the reference-to-physical transformation, and place them in the provided ...
TensorData< PointScalar, DeviceType > allocateCellMeasure(const Data< PointScalar, DeviceType > &jacobianDet, const TensorData< PointScalar, DeviceType > &cubatureWeights) const
Allocate a TensorData object appropriate for passing to computeCellMeasure().
KOKKOS_INLINE_FUNCTION int numCells() const
Returns the number of cells.
KOKKOS_INLINE_FUNCTION int numNodesPerCell() const
Returns the number of nodes per cell; may be more than the number of vertices in the corresponding Ce...
Data< PointScalar, DeviceType > getJacobianRefData(const ScalarView< PointScalar, DeviceType > &points) const
Computes reference-space data for the specified points, to be used in setJacobian().
KOKKOS_INLINE_FUNCTION int numCellsInDimension(const int &dim) const
For uniform grid and tensor grid CellGeometry, returns the number of cells in the specified component...
KOKKOS_INLINE_FUNCTION Orientation getOrientation(int &cellNumber) const
Returns the orientation for the specified cell. Requires that initializeOrientations() has been calle...
KOKKOS_INLINE_FUNCTION PointScalar gridCellCoordinate(const int &gridCellOrdinal, const int &localNodeNumber, const int &dim) const
returns coordinate in dimension dim of the indicated node in the indicated grid cell
CellGeometry(const Kokkos::Array< PointScalar, spaceDim > &origin, const Kokkos::Array< PointScalar, spaceDim > &domainExtents, const Kokkos::Array< int, spaceDim > &gridCellCounts, SubdivisionStrategy subdivisionStrategy=NO_SUBDIVISION, HypercubeNodeOrdering nodeOrdering=HYPERCUBE_NODE_ORDER_TENSOR)
Uniform grid constructor, with optional subdivision into simplices.
void initializeOrientations()
Initialize the internal orientations_ member with the orientations of each member cell....
Data< PointScalar, DeviceType > allocateJacobianDataPrivate(const ScalarView< PointScalar, DeviceType > &pointComponentView, const int &pointsPerCell, const int startCell, const int endCell) const
Notionally-private method that provides a common interface for multiple public-facing allocateJacobia...
KOKKOS_INLINE_FUNCTION DataVariationType cellVariationType() const
KOKKOS_INLINE_FUNCTION int hypercubeComponentNodeNumber(int hypercubeNodeNumber, int d) const
For hypercube vertex number hypercubeNodeNumber, returns the component node number in specified dimen...
KOKKOS_INLINE_FUNCTION PointScalar subdivisionCoordinate(const int &gridCellOrdinal, const int &subdivisionOrdinal, const int &subdivisionNodeNumber, const int &d) const
returns coordinate in dimension d for the indicated subdivision of the indicated grid cell
@ NO_SUBDIVISION
no subdivision
KOKKOS_INLINE_FUNCTION unsigned rank() const
Returns the logical rank of this container. This is always 3.
@ FIRST_ORDER
geometry expressible in terms of vertices of the cell
@ HIGHER_ORDER
geometry expressible in terms of a higher-order basis (must be specified)
KOKKOS_INLINE_FUNCTION bool affine() const
Returns true if Jacobian is constant within each cell.
Data< PointScalar, DeviceType > allocateJacobianData(const TensorPoints< PointScalar, DeviceType > &points, const int startCell=0, const int endCell=-1) const
Allocate a container into which Jacobians of the reference-to-physical mapping can be placed.
@ HYPERCUBE_NODE_ORDER_CLASSIC_SHARDS
classic shards ordering
KOKKOS_INLINE_FUNCTION ~CellGeometry()
Destructor.
KOKKOS_INLINE_FUNCTION HypercubeNodeOrdering nodeOrderingForHypercubes() const
Returns the node ordering used for hypercubes.
Data< Orientation, DeviceType > getOrientations()
Returns the orientations for all cells. Calls initializeOrientations() if it has not previously been ...
KOKKOS_INLINE_FUNCTION int uniformJacobianModulus() const
Returns an integer indicating the number of distinct cell types vis-a-vis Jacobians.
KOKKOS_INLINE_FUNCTION int gridCellNodeForSubdivisionNode(const int &gridCellOrdinal, const int &subdivisionOrdinal, const int &subdivisionNodeNumber) const
returns coordinate in dimension d for the indicated subdivision of the indicated grid cell
const shards::CellTopology & cellTopology() const
The shards CellTopology for each cell within the CellGeometry object. Note that this is always a lowe...
KOKKOS_INLINE_FUNCTION PointScalar operator()(const int &cell, const int &node, const int &dim) const
Return the coordinate (weight) of the specified node. For straight-edged geometry,...
KOKKOS_INLINE_FUNCTION std::enable_if< std::is_integral< iType >::value, int >::type extent_int(const iType &r) const
Returns the logical extent of the container in the specified dimension as an int; the shape of CellGe...
Wrapper around a Kokkos::View that allows data that is constant or repeating in various logical dimen...
KOKKOS_INLINE_FUNCTION enable_if_t< rank==1, const Kokkos::View< typename RankExpander< DataScalar, rank >::value_type, DeviceType > & > getUnderlyingView() const
Returns the underlying view. Throws an exception if the underlying view is not rank 1.
KOKKOS_INLINE_FUNCTION int extent_int(const int &r) const
Returns the logical extent in the specified dimension.
KOKKOS_INLINE_FUNCTION const Kokkos::View< DataScalar **, DeviceType > & getUnderlyingView2() const
returns the View that stores the unique data. For rank-2 underlying containers.
KOKKOS_INLINE_FUNCTION constexpr bool isValid() const
returns true for containers that have data; false for those that don't (namely, those that have been ...
KOKKOS_INLINE_FUNCTION const Kokkos::View< DataScalar ***, DeviceType > & getUnderlyingView3() const
returns the View that stores the unique data. For rank-3 underlying containers.
KOKKOS_INLINE_FUNCTION const Kokkos::Array< DataVariationType, 7 > & getVariationTypes() const
Returns an array with the variation types in each logical dimension.
KOKKOS_INLINE_FUNCTION const Kokkos::View< DataScalar *, DeviceType > & getUnderlyingView1() const
returns the View that stores the unique data. For rank-1 underlying containers.
KOKKOS_INLINE_FUNCTION int getDataExtent(const ordinal_type &d) const
returns the true extent of the data corresponding to the logical dimension provided; if the data does...
KOKKOS_INLINE_FUNCTION Kokkos::Array< int, 7 > getExtents() const
Returns an array containing the logical extents in each dimension.
KOKKOS_INLINE_FUNCTION unsigned rank() const
Returns the logical rank of the Data container.
A family of basis functions, constructed from H(vol) and H(grad) bases on the line.
Store host-only "members" of CellGeometry using a static map indexed on the CellGeometry pointer....
Functor for full (C,P) Jacobian determinant container. CUDA compiler issues led us to avoid lambdas f...
Orientation encoding and decoding.
View-like interface to tensor data; tensor components are stored separately and multiplied together a...
KOKKOS_INLINE_FUNCTION const Data< Scalar, DeviceType > & getTensorComponent(const ordinal_type &r) const
Returns the requested tensor component.
KOKKOS_INLINE_FUNCTION ordinal_type rank() const
Returns the rank of the container.
KOKKOS_INLINE_FUNCTION std::enable_if< std::is_integral< iType >::value, ordinal_type >::type extent_int(const iType &d) const
Returns the logical extent in the requested dimension.
KOKKOS_INLINE_FUNCTION ordinal_type numTensorComponents() const
Return the number of tensorial components.
View-like interface to tensor points; point components are stored separately; the appropriate coordin...
KOKKOS_INLINE_FUNCTION ScalarView< PointScalar, DeviceType > getTensorComponent(const ordinal_type &r) const
Returns the requested tensor component.
KOKKOS_INLINE_FUNCTION std::enable_if< std::is_integral< iType >::value, int >::type extent_int(const iType &r) const
Returns the logical extent in the requested dimension.