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 __CUDA_ARCH__ 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 __CUDA_ARCH__ 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();
218 blockPlusDiagonalLastNonDiagonal = -1;
223 else if (cellGeometryType_ == TENSOR_GRID)
225 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(
true, std::invalid_argument,
"tensor grid support not yet implemented");
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);
268 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(
true, std::invalid_argument,
"support for this CellGeometryType is not yet implemented");
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)
333 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(
true, std::invalid_argument,
"tensor grid support not yet implemented");
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");
388 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(!refData.
isValid(), std::invalid_argument,
"refData should be a valid container for cases with non-affine geometry");
389 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(refData.
rank() != 3, std::invalid_argument,
"refData should have shape (F,P,D)");
390 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(refData.
extent_int(0) != basisForNodes->getCardinality(), std::invalid_argument,
"refData should have shape (F,P,D)");
391 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(refData.
extent_int(1) != points.
extent_int(0), std::invalid_argument,
"refData should have shape (F,P,D)");
392 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(refData.
extent_int(2) != spaceDim, std::invalid_argument,
"refData should have shape (F,P,D)");
401 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(
true, std::invalid_argument,
"support for this CellGeometryType is not yet implemented");
406 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
408 const Kokkos::Array<PointScalar,spaceDim> &domainExtents,
409 const Kokkos::Array<int,spaceDim> &gridCellCounts,
413 nodeOrdering_(nodeOrdering),
414 cellGeometryType_(UNIFORM_GRID),
415 subdivisionStrategy_(subdivisionStrategy),
418 domainExtents_(domainExtents),
419 gridCellCounts_(gridCellCounts)
422 for (
int d=0; d<spaceDim; d++)
424 numCells_ *= gridCellCounts_[d];
428 shards::CellTopology cellTopo;
432 numNodesPerCell_ = 1 << spaceDim;
436 cellTopo = shards::CellTopology(shards::getCellTopologyData<shards::Line<> >());
438 else if (spaceDim == 2)
440 cellTopo = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<> >());
442 else if (spaceDim == 3)
444 cellTopo = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<> >());
454 numNodesPerCell_ = spaceDim + 1;
457 cellTopo = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<> >());
459 else if (spaceDim == 3)
461 cellTopo = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<> >());
470 const int linearPolyOrder = 1;
471 BasisPtr
basisForNodes = getBasis<BasisFamily>(cellTopo, FUNCTION_SPACE_HGRAD, linearPolyOrder);
477 if (cellTopo.getKey() == shards::Quadrilateral<>::key)
481 else if (cellTopo.getKey() == shards::Hexahedron<>::key)
488 HostMemberLookup::constructorCalled(
this, cellTopo,
basisForNodes);
494 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
496 ScalarView<int,DeviceType> cellToNodes,
497 ScalarView<PointScalar,DeviceType> nodes,
498 const bool claimAffine,
501 nodeOrdering_(nodeOrdering),
502 cellGeometryType_(FIRST_ORDER),
503 cellToNodes_(cellToNodes),
506 numCells_ = cellToNodes.extent_int(0);
507 numNodesPerCell_ = cellToNodes.extent_int(1);
508 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(numNodesPerCell_ != cellTopo.getNodeCount(), std::invalid_argument,
"cellToNodes.extent(1) does not match the cell topology node count");
513 const bool simplicialTopo = (cellTopo.getNodeCount() == cellTopo.getDimension() + 1);
514 affine_ = simplicialTopo;
522 const int linearPolyOrder = 1;
523 BasisPtr
basisForNodes = getBasis<BasisFamily>(cellTopo, FUNCTION_SPACE_HGRAD, linearPolyOrder);
529 if (cellTopo.getKey() == shards::Quadrilateral<>::key)
533 else if (cellTopo.getKey() == shards::Hexahedron<>::key)
540 HostMemberLookup::constructorCalled(
this, cellTopo,
basisForNodes);
544 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
546 ScalarView<PointScalar,DeviceType> cellNodes)
548 nodeOrdering_(HYPERCUBE_NODE_ORDER_TENSOR),
549 cellGeometryType_(HIGHER_ORDER),
552 numCells_ = cellNodes.extent_int(0);
553 numNodesPerCell_ = cellNodes.extent_int(1);
556 const bool firstOrderGeometry = (
basisForNodes->getDegree() == 1);
559 shards::CellTopology cellTopo =
basisForNodes->getBaseCellTopology();
561 if (firstOrderGeometry && (cellTopo.getNodeCount() == spaceDim + 1))
570 HostMemberLookup::constructorCalled(
this, cellTopo,
basisForNodes);
573 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
574 KOKKOS_INLINE_FUNCTION
580 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
589 INTREPID2_TEST_FOR_EXCEPTION(cubatureWeights.
rank() != 1, std::invalid_argument,
"cubatureWeights container must have shape (P)");
592 std::vector< Data<PointScalar,DeviceType> > tensorComponents(numTensorComponents);
596 const int cellExtent = jacobianDet.
extent_int(0);
599 Kokkos::Array<int,7> cellExtents{cellExtent,1,1,1,1,1,1};
601 ScalarView<PointScalar,DeviceType> detDataView (
"cell relative volumes", cellDataDim);
604 for (
int cubTensorComponent=0; cubTensorComponent<numTensorComponents-1; cubTensorComponent++)
607 const auto cubatureExtents = cubatureComponent.getExtents();
608 const auto cubatureVariationTypes = cubatureComponent.getVariationTypes();
609 const int numPoints = cubatureComponent.getDataExtent(0);
610 ScalarView<PointScalar,DeviceType> cubatureWeightView (
"cubature component weights", numPoints);
611 const int pointComponentRank = 1;
612 tensorComponents[cubTensorComponent+1] =
Data<PointScalar,DeviceType>(cubatureWeightView,pointComponentRank,cubatureExtents,cubatureVariationTypes);
617 const int cellExtent = jacobianDet.
extent_int(0);
621 const int numPoints = cubatureWeights.
extent_int(0);
622 Kokkos::Array<int,7> extents{cellExtent,numPoints,1,1,1,1,1};
624 ScalarView<PointScalar,DeviceType> cubatureWeightView;
627 cubatureWeightView = ScalarView<PointScalar,DeviceType>(
"cell measure", cellDataDim, numPoints);
631 cubatureWeightView = ScalarView<PointScalar,DeviceType>(
"cell measure", numPoints);
633 const int cellMeasureRank = 2;
636 const bool separateFirstComponent = (numTensorComponents > 1);
640 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
650 "cellMeasure must either have a tensor component count of 1 or a tensor component count that is one higher than that of cubatureWeights");
652 INTREPID2_TEST_FOR_EXCEPTION(cubatureWeights.
rank() != 1, std::invalid_argument,
"cubatureWeights container must have shape (P)");
659 for (
int i=1; i<numTensorDimensions+1; i++)
668 const bool detCellVaries = detVaries[0] !=
CONSTANT;
669 const bool detPointVaries = detVaries[1] !=
CONSTANT;
671 if (detCellVaries && detPointVaries)
675 const int numCells = detData.extent_int(0);
676 const int numPoints = detData.extent_int(1);
677 INTREPID2_TEST_FOR_EXCEPTION(numCells != cellMeasureData.extent_int(0), std::invalid_argument,
"cellMeasureData doesn't match jacobianDet in cell dimension");
678 INTREPID2_TEST_FOR_EXCEPTION(numPoints != cellMeasureData.extent_int(1), std::invalid_argument,
"cellMeasureData doesn't match jacobianDet in point dimension");
683 using ExecutionSpace =
typename DeviceType::execution_space;
684 Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<2>> rangePolicy({0,0},{numCells,numPoints});
685 Kokkos::parallel_for(rangePolicy, cellMeasureFunctor);
687 else if (detCellVaries && !detPointVaries)
691 using ExecutionSpace =
typename DeviceType::execution_space;
692 Kokkos::parallel_for(
693 Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<2>>({0,0},{detData.extent_int(0),cubatureWeights.
extent_int(0)}),
694 KOKKOS_LAMBDA (
int cellOrdinal,
int pointOrdinal) {
695 cellMeasureData(cellOrdinal,pointOrdinal) = detData(cellOrdinal) * cubatureWeights(pointOrdinal);
704 using ExecutionSpace =
typename DeviceType::execution_space;
705 Kokkos::parallel_for(Kokkos::RangePolicy<ExecutionSpace>(0,cellMeasureData.extent_int(0)),
706 KOKKOS_LAMBDA (
const int &pointOrdinal) {
707 cellMeasureData(pointOrdinal) = detData(0) * cubatureWeights(pointOrdinal);
713 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
714 typename CellGeometry<PointScalar,spaceDim,DeviceType>::BasisPtr
718 return HostMemberLookup::getBasis(
this);
721 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
725 return HostMemberLookup::getCellTopology(
this);
728 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
729 KOKKOS_INLINE_FUNCTION
732 if (cellGeometryType_ == UNIFORM_GRID)
734 const int numSubdivisions = numCellsPerGridCell(subdivisionStrategy_);
735 if (numSubdivisions == 1)
747 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
751 if (cellGeometryType_ == UNIFORM_GRID)
756 else if (cellGeometryType_ == TENSOR_GRID)
761 else if ((cellGeometryType_ == FIRST_ORDER) || (cellGeometryType_ == HIGHER_ORDER))
763 const bool simplex = (spaceDim + 1 == cellToNodes_.extent_int(1));
771 auto basisForNodes = this->basisForNodes();
780 auto basisGradients = basisForNodes->allocateBasisValues(points, OPERATOR_GRAD);
781 basisForNodes->getValues(basisGradients, points, OPERATOR_GRAD);
784 int numFields = basisForNodes->getCardinality();
792 auto basisGradientsView =
getMatchingViewWithLabel(firstPointComponentView,
"CellGeometryProvider: temporary basisGradients", numFields, numPoints, spaceDim);
794 using ExecutionSpace =
typename DeviceType::execution_space;
795 auto policy = Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<3>>({0,0,0},{numFields,numPoints,spaceDim});
797 Kokkos::parallel_for(
"copy basis gradients", policy,
798 KOKKOS_LAMBDA (
const int &fieldOrdinal,
const int &pointOrdinal,
const int &d) {
799 basisGradientsView(fieldOrdinal,pointOrdinal,d) = basisGradients(fieldOrdinal,pointOrdinal,d);
801 ExecutionSpace().fence();
811 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(
true, std::invalid_argument,
"support for this CellGeometryType is not yet implemented");
816 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
817 KOKKOS_INLINE_FUNCTION
820 if (nodeOrdering_ == HYPERCUBE_NODE_ORDER_CLASSIC_SHARDS)
826 if ((hypercubeNodeNumber % 4 == 1) || (hypercubeNodeNumber % 4 == 2))
833 if ((hypercubeNodeNumber % 4 == 2) || (hypercubeNodeNumber % 4 == 3))
840 const int nodesForPriorDimensions = 1 << d;
841 if ((hypercubeNodeNumber / nodesForPriorDimensions) % 2 == 1)
847 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
850 using HostExecSpace =
typename Kokkos::Impl::is_space<DeviceType>::host_mirror_space::execution_space ;
852 const bool isGridType = (cellGeometryType_ == TENSOR_GRID) || (cellGeometryType_ == UNIFORM_GRID);
853 const int numOrientations = isGridType ? numCellsPerGridCell(subdivisionStrategy_) : numCells();
855 const int nodesPerCell = numNodesPerCell();
857 ScalarView<Orientation, DeviceType> orientationsView(
"orientations", numOrientations);
858 auto orientationsHost = Kokkos::create_mirror_view(
typename HostExecSpace::memory_space(), orientationsView);
860 ScalarView<PointScalar, HostExecSpace> cellNodesHost(
"cellNodesHost",numOrientations,nodesPerCell);
869 const int numSubdivisions = numCellsPerGridCell(subdivisionStrategy_);
870 const int gridCellOrdinal = 0;
871 auto hostPolicy = Kokkos::MDRangePolicy<HostExecSpace,Kokkos::Rank<2>>({0,0},{numSubdivisions,nodesPerCell});
872 Kokkos::parallel_for(
"fill cellNodesHost", hostPolicy,
873 KOKKOS_LAMBDA (
const int &subdivisionOrdinal,
const int &nodeInCell) {
874 auto node = gridCellNodeForSubdivisionNode(gridCellOrdinal, subdivisionOrdinal, nodeInCell);
875 cellNodesHost(subdivisionOrdinal,nodeInCell) = node;
883 auto cellToNodesHost = Kokkos::create_mirror_view_and_copy(
typename HostExecSpace::memory_space(), cellToNodes_);
887 Kokkos::deep_copy(orientationsView,orientationsHost);
889 const int orientationsRank = 1;
890 const Kokkos::Array<int,7> orientationExtents {
static_cast<int>(numCells_),1,1,1,1,1,1};
895 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
896 KOKKOS_INLINE_FUNCTION
904 return numNodesPerCell_;
916 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
917 template <
typename iType>
918 KOKKOS_INLINE_FUNCTION
919 typename std::enable_if<std::is_integral<iType>::value,
int>::type
922 return static_cast<int>(extent(r));
925 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
926 KOKKOS_INLINE_FUNCTION
930 return nodeOrdering_;
933 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
934 KOKKOS_INLINE_FUNCTION
940 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
941 KOKKOS_INLINE_FUNCTION
944 if (cellGeometryType_ == UNIFORM_GRID)
946 return gridCellCounts_[dim];
948 else if (cellGeometryType_ == TENSOR_GRID)
950 return tensorVertices_.extent_int(dim);
958 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
959 KOKKOS_INLINE_FUNCTION
962 return numNodesPerCell_;
965 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
966 KOKKOS_INLINE_FUNCTION
969 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(!orientations_.isValid(), std::invalid_argument,
"orientations_ not initialized; call initializeOrientations() first");
970 return orientations_(cellNumber);
973 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
976 if (!orientations_.isValid())
978 initializeOrientations();
980 return orientations_;
983 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
984 KOKKOS_INLINE_FUNCTION
987 const int componentNode = hypercubeComponentNodeNumber(localNodeNumber, dim);
988 int cellCountForPriorDimensions = 1;
989 for (
int d=0; d<dim; d++)
991 cellCountForPriorDimensions *= numCellsInDimension(d);
993 const int componentGridCellOrdinal = (gridCellOrdinal / cellCountForPriorDimensions) % numCellsInDimension(dim);
994 const int vertexOrdinal = componentGridCellOrdinal + componentNode;
995 if (cellGeometryType_ == UNIFORM_GRID)
997 return origin_[dim] + (vertexOrdinal * domainExtents_[dim]) / gridCellCounts_[dim];
999 else if (cellGeometryType_ == TENSOR_GRID)
1001 Kokkos::Array<int,spaceDim> pointOrdinalComponents;
1002 for (
int d=0; d<spaceDim; d++)
1004 pointOrdinalComponents[d] = 0;
1006 pointOrdinalComponents[dim] = vertexOrdinal;
1007 return tensorVertices_(pointOrdinalComponents,dim);
1011 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(
true, std::invalid_argument,
"Unsupported geometry type");
1016 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1017 KOKKOS_INLINE_FUNCTION
1023 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1024 KOKKOS_INLINE_FUNCTION
1026 const int &subdivisionNodeNumber)
const 1029 switch (subdivisionStrategy_)
1031 case NO_SUBDIVISION:
1032 return subdivisionNodeNumber;
1033 case TWO_TRIANGLES_RIGHT:
1034 case TWO_TRIANGLES_LEFT:
1035 case FOUR_TRIANGLES:
1037 Kokkos::Array<int,3> nodeLookup;
1038 if (subdivisionStrategy_ == TWO_TRIANGLES_RIGHT)
1040 if (subdivisionOrdinal == 0)
1043 nodeLookup = {0,1,2};
1045 else if (subdivisionOrdinal == 1)
1050 nodeLookup = {2,3,0};
1054 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(
true, std::invalid_argument,
"Unsupported subdivision ordinal");
1057 else if (subdivisionStrategy_ == TWO_TRIANGLES_LEFT)
1059 if (subdivisionOrdinal == 0)
1065 nodeLookup = {3,0,1};
1067 else if (subdivisionOrdinal == 1)
1072 nodeLookup = {2,3,0};
1076 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(
true, std::invalid_argument,
"Unsupported subdivision ordinal");
1088 if (subdivisionNodeNumber == 1)
1095 nodeLookup = {(subdivisionOrdinal + 1) % 4, -1, subdivisionOrdinal};
1098 const int gridCellNodeNumber = nodeLookup[subdivisionNodeNumber];
1099 return gridCellNodeNumber;
1101 case FIVE_TETRAHEDRA:
1102 case SIX_TETRAHEDRA:
1104 Kokkos::Array<int,4> nodeLookup;
1105 if (subdivisionStrategy_ == FIVE_TETRAHEDRA)
1135 switch (subdivisionOrdinal) {
1137 nodeLookup = {1,3,4,6};
1140 nodeLookup = {0,1,3,4};
1143 nodeLookup = {1,2,3,6};
1146 nodeLookup = {1,4,5,6};
1149 nodeLookup = {3,4,6,7};
1152 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(
true, std::invalid_argument,
"invalid subdivisionOrdinal");
1156 else if (subdivisionStrategy_ == SIX_TETRAHEDRA)
1158 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(
true, std::invalid_argument,
"support for SIX_TETRAHEDRA not yet implemented");
1160 const int gridCellNodeNumber = nodeLookup[subdivisionNodeNumber];
1161 return gridCellNodeNumber;
1164 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(
true, std::invalid_argument,
"Subdivision strategy not yet implemented!");
1170 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1171 KOKKOS_INLINE_FUNCTION
1173 const int &subdivisionNodeNumber,
const int &d)
const 1175 int gridCellNode = gridCellNodeForSubdivisionNode(gridCellOrdinal, subdivisionOrdinal, subdivisionNodeNumber);
1177 if (subdivisionStrategy_ == FOUR_TRIANGLES)
1180 if (gridCellNode == 4)
1184 const int gridVertex0 = 0;
1185 const int gridVertex1 = (d == 0) ? 1 : 3;
1186 return 0.5 * (gridCellCoordinate(gridCellOrdinal, gridVertex0, d) + gridCellCoordinate(gridCellOrdinal, gridVertex1, d));
1189 return gridCellCoordinate(gridCellOrdinal, gridCellNode, d);
1192 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1193 KOKKOS_INLINE_FUNCTION
1196 if ((cellGeometryType_ == UNIFORM_GRID) || (cellGeometryType_ == TENSOR_GRID))
1198 const int numSubdivisions = numCellsPerGridCell(subdivisionStrategy_);
1199 if (numSubdivisions == 1)
1202 return gridCellCoordinate(cell, node, dim);
1206 const int subdivisionOrdinal = cell % numSubdivisions;
1207 const int gridCellOrdinal = cell / numSubdivisions;
1208 return subdivisionCoordinate(gridCellOrdinal, subdivisionOrdinal, node, dim);
1213 #ifdef HAVE_INTREPID2_DEBUG 1214 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE((cell < 0), std::invalid_argument,
"cell out of bounds");
1215 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(static_cast<unsigned>(cell) > numCells_, std::invalid_argument,
"cell out of bounds");
1216 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE((node < 0), std::invalid_argument,
"node out of bounds");
1217 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(static_cast<unsigned>(node) > numNodesPerCell_, std::invalid_argument,
"node out of bounds");
1218 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE((dim < 0), std::invalid_argument,
"dim out of bounds" );
1219 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(dim > spaceDim, std::invalid_argument,
"dim out of bounds" );
1221 if (cellToNodes_.is_allocated())
1223 const int nodeNumber = cellToNodes_(cell,node);
1224 return nodes_(nodeNumber,dim);
1228 return nodes_(cell,node,dim);
1233 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1234 KOKKOS_INLINE_FUNCTION
1237 if (cellGeometryType_ == UNIFORM_GRID)
1239 return numCellsPerGridCell(subdivisionStrategy_);
1247 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1250 const int pointsPerCell = points.
extent_int(0);
1251 return allocateJacobianDataPrivate(points,pointsPerCell,startCell,endCell);
1254 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1258 const int pointDimension = (points.rank() == 3) ? 1 : 0;
1259 const int pointsPerCell = points.extent_int(pointDimension);
1261 return allocateJacobianDataPrivate(tensorPoints,pointsPerCell,startCell,endCell);
1264 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1267 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(!affine_, std::invalid_argument,
"this version of allocateJacobianData() is only supported for affine CellGeometry");
1270 return allocateJacobianDataPrivate(emptyPoints,numPoints,startCell,endCell);
1273 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1277 const int pointsPerCell = points.
extent_int(0);
1278 setJacobianDataPrivate(jacobianData,points,pointsPerCell,refData,startCell,endCell);
1281 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1286 const int pointDimension = (points.rank() == 3) ? 1 : 0;
1287 const int pointsPerCell = points.extent_int(pointDimension);
1289 setJacobianDataPrivate(jacobianData,tensorPoints,pointsPerCell,refData,startCell,endCell);
1292 template<
class Po
intScalar,
int spaceDim,
typename DeviceType>
1295 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(!affine_, std::invalid_argument,
"this version of setJacobian() is only supported for affine CellGeometry");
1299 setJacobianDataPrivate(jacobianData,emptyPoints,numPoints,emptyRefData,startCell,endCell);
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Quadrilateral cell...
geometry expressible in terms of vertices of the cell
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 int numCells() const
Returns the number of cells.
KOKKOS_INLINE_FUNCTION int hypercubeComponentNodeNumber(int hypercubeNodeNumber, int d) const
For hypercube vertex number hypercubeNodeNumber, returns the component node number in specified dimen...
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 int numCellsInDimension(const int &dim) const
For uniform grid and tensor grid CellGeometry, returns the number of cells in the specified component...
CellGeometry provides the nodes for a set of cells; has options that support efficient definition of ...
KOKKOS_INLINE_FUNCTION bool affine() const
Returns true if Jacobian is constant within each cell.
Functor for full (C,P) Jacobian determinant container. CUDA compiler issues led us to avoid lambdas f...
View-like interface to tensor points; point components are stored separately; the appropriate coordin...
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
void setJacobianDataPrivate(Data< PointScalar, DeviceType > &jacobianData, const TensorPoints< PointScalar, DeviceType > &points, 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...
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.
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 constexpr bool isValid() const
returns true for containers that have data; false for those that don't (namely, those that have been ...
one of two dimensions in a matrix; bottom-right part of matrix is diagonal
KOKKOS_INLINE_FUNCTION unsigned rank() const
Returns the logical rank of this container. This is always 3.
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...
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.
Wrapper around a Kokkos::View that allows data that is constant or repeating in various logical dimen...
varies according to modulus of the index
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, this is simply the physical coordinate of the vertex. For all geometries, this can be understood as a weight on the corresponding H^1 basis function used in the reference-to-physical map.
KOKKOS_INLINE_FUNCTION const Kokkos::View< DataScalar *, DeviceType > & getUnderlyingView1() const
returns the View that stores the unique data. For rank-1 underlying containers.
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 ...
KOKKOS_INLINE_FUNCTION unsigned rank() const
Returns the logical rank of the Data container.
KOKKOS_INLINE_FUNCTION DataVariationType cellVariationType() const
KOKKOS_INLINE_FUNCTION HypercubeNodeOrdering nodeOrderingForHypercubes() const
Returns the node ordering used for hypercubes.
KOKKOS_INLINE_FUNCTION ordinal_type rank() const
Returns the rank of the container.
Orientation encoding and decoding.
Data< PointScalar, DeviceType > getJacobianRefData(const TensorPoints< PointScalar, DeviceType > &points) const
Computes reference-space data for the specified points, to be used in setJacobian().
KOKKOS_INLINE_FUNCTION int extent_int(const int &r) const
Returns the logical extent in the specified dimension.
Store host-only "members" of CellGeometry using a static map indexed on the CellGeometry pointer...
geometry expressible in terms of a higher-order basis (must be specified)
KOKKOS_INLINE_FUNCTION int uniformJacobianModulus() const
Returns an integer indicating the number of distinct cell types vis-a-vis Jacobians.
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...
void initializeOrientations()
Initialize the internal orientations_ member with the orientations of each member cell...
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 ...
Data< Orientation, DeviceType > getOrientations()
Returns the orientations for all cells. Calls initializeOrientations() if it has not previously been ...
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 const Data< Scalar, DeviceType > & getTensorComponent(const ordinal_type &r) const
Returns the requested tensor component.
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 ...
A family of nodal basis functions which is related to, but not identical with, the Lagrangian basis f...
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.
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.
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 ...
KOKKOS_INLINE_FUNCTION ScalarView< PointScalar, DeviceType > getTensorComponent(const ordinal_type &r) const
Returns the requested tensor component.
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 ...
KOKKOS_INLINE_FUNCTION ordinal_type numTensorComponents() const
Return the number of tensorial components.
BasisPtr basisForNodes() const
H^1 Basis used in the reference-to-physical transformation. Linear for straight-edged geometry; highe...
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 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 const Kokkos::View< DataScalar **, DeviceType > & getUnderlyingView2() const
returns the View that stores the unique data. For rank-2 underlying containers.
KOKKOS_INLINE_FUNCTION ~CellGeometry()
Destructor.
Data< PointScalar, DeviceType > allocateJacobianDataPrivate(const TensorPoints< PointScalar, DeviceType > &points, const int &pointsPerCell, const int startCell, const int endCell) const
Notionally-private method that provides a common interface for multiple public-facing allocateJacobia...
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Hexahedron cell...
DataVariationType
Enumeration to indicate how data varies in a particular dimension of an Intrepid2::Data object...
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 Orientation getOrientation(int &cellNumber) const
Returns the orientation for the specified cell. Requires that initializeOrientations() has been calle...
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...
KOKKOS_INLINE_FUNCTION const Kokkos::Array< DataVariationType, 7 > & getVariationTypes() const
Returns an array with the variation types in each logical dimension.