Intrepid2
Intrepid2_CellGeometryDef.hpp
1// @HEADER
2// ************************************************************************
3//
4// Intrepid2 Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Kyungjoo Kim (kyukim@sandia.gov),
38// Mauro Perego (mperego@sandia.gov), or
39// Nate Roberts (nvrober@sandia.gov),
40//
41// ************************************************************************
42// @HEADER
43
50#ifndef Intrepid2_CellGeometryDef_h
51#define Intrepid2_CellGeometryDef_h
52
53namespace Intrepid2
54{
55
56 namespace Impl
57 {
60 template<class PointScalar, int spaceDim, typename DeviceType>
62 {
63 using BasisPtr = Teuchos::RCP<Intrepid2::Basis<DeviceType,PointScalar,PointScalar> >;
65 public:
66 // conceptually, these should be private members, but for the definition of these, we need them to be externally accessible.
67 static std::map<const CellGeometryType *, shards::CellTopology> cellTopology_;
68 static std::map<const CellGeometryType *, BasisPtr> basisForNodes_;
69
70 public:
71 static void constructorCalled(const CellGeometryType *cellGeometry, const shards::CellTopology &cellTopo, BasisPtr basisForNodes)
72 {
73 cellTopology_[cellGeometry] = cellTopo;
74 basisForNodes_[cellGeometry] = basisForNodes;
75 }
76
77 static void destructorCalled(const CellGeometryType *cellGeometry)
78 {
79 cellTopology_.erase(cellGeometry);
80 basisForNodes_.erase(cellGeometry);
81 }
82
83 static BasisPtr getBasis(const CellGeometryType *cellGeometry)
84 {
85 return basisForNodes_[cellGeometry];
86 }
87
88 static const shards::CellTopology & getCellTopology(const CellGeometryType *cellGeometry)
89 {
90 return cellTopology_[cellGeometry];
91 }
92 };
93
94 // member lookup map definitions for CellGeometryHostMembers:
95 template< class PointScalar, int spaceDim, typename DeviceType > typename std::map<const CellGeometry<PointScalar,spaceDim,DeviceType> *, shards::CellTopology> CellGeometryHostMembers< PointScalar,spaceDim,DeviceType>::cellTopology_;
96
97 template< class PointScalar, int spaceDim, typename DeviceType > typename std::map<const CellGeometry<PointScalar,spaceDim,DeviceType> *, Teuchos::RCP<Intrepid2::Basis<DeviceType,PointScalar,PointScalar> >> CellGeometryHostMembers< PointScalar,spaceDim,DeviceType>::basisForNodes_;
98
101 template<class PointScalar, int spaceDim, typename DeviceType>
103 {
104 Kokkos::View<PointScalar**, DeviceType> cellMeasures_; // (C,P)
105 Kokkos::View<PointScalar**, DeviceType> detData_; // (C,P)
106 TensorData<PointScalar,DeviceType> cubatureWeights_; // (P)
107 public:
108 CellMeasureFunctor(Kokkos::View<PointScalar**, DeviceType> cellMeasures,
109 Kokkos::View<PointScalar**, DeviceType> detData, TensorData<PointScalar,DeviceType> cubatureWeights)
110 :
111 cellMeasures_(cellMeasures),
112 detData_(detData),
113 cubatureWeights_(cubatureWeights)
114 {}
115
116 KOKKOS_INLINE_FUNCTION void
117 operator () (const ordinal_type cellOrdinal, const ordinal_type pointOrdinal) const
118 {
119 cellMeasures_(cellOrdinal,pointOrdinal) = detData_(cellOrdinal,pointOrdinal) * cubatureWeights_(pointOrdinal);
120 }
121 };
122 }
123
124 template<class PointScalar, int spaceDim, typename DeviceType>
125 KOKKOS_INLINE_FUNCTION
127:
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_)
141 {
142 // host-only registration with HostMemberLookup:
143#ifndef INTREPID2_COMPILE_DEVICE_CODE
144 shards::CellTopology cellTopo = cellGeometry.cellTopology();
145 BasisPtr basisForNodes = cellGeometry.basisForNodes();
147 HostMemberLookup::constructorCalled(this, cellTopo, basisForNodes);
148#endif
149 }
150
151 template<class PointScalar, int spaceDim, typename DeviceType>
152 KOKKOS_INLINE_FUNCTION
154 {
155 // host-only deregistration with HostMemberLookup:
156#ifndef INTREPID2_COMPILE_DEVICE_CODE
158 HostMemberLookup::destructorCalled(this);
159#endif
160 }
161
162 template<class PointScalar, int spaceDim, typename DeviceType>
163 KOKKOS_INLINE_FUNCTION
165 {
166 switch (subdivisionStrategy) {
167 case NO_SUBDIVISION:
168 return 1;
169 case TWO_TRIANGLES_LEFT:
170 case TWO_TRIANGLES_RIGHT:
171 return 2;
172 case FOUR_TRIANGLES:
173 return 4;
174 case FIVE_TETRAHEDRA:
175 return 5;
176 case SIX_TETRAHEDRA:
177 return 6;
178 }
179 return -1;
180 }
181
182 template<class PointScalar, int spaceDim, typename DeviceType>
184 CellGeometry<PointScalar,spaceDim,DeviceType>::allocateJacobianDataPrivate(const ScalarView<PointScalar,DeviceType> &pointComponentView, const int &pointsPerCell, const int startCell, const int endCell) const
185 {
186 ScalarView<PointScalar,DeviceType> data;
187 const int rank = 4; // C,P,D,D
188 const int CELL_DIM = 0;
189 const int POINT_DIM = 1;
190 const int D1_DIM = 2;
191 const int D2_DIM = 3;
192
193 const int numCellsWorkset = (endCell == -1) ? (numCells_ - startCell) : (endCell - startCell);
194
195 Kokkos::Array<int,7> extents { numCellsWorkset, pointsPerCell, spaceDim, spaceDim, 1, 1, 1 };
196 Kokkos::Array<DataVariationType,7> variationType { CONSTANT, CONSTANT, CONSTANT, CONSTANT, CONSTANT, CONSTANT, CONSTANT };
197
198 int blockPlusDiagonalLastNonDiagonal = -1;
199
200 if (cellGeometryType_ == UNIFORM_GRID)
201 {
202 if (uniformJacobianModulus() != 1)
203 {
204 variationType[CELL_DIM] = MODULAR;
205 variationType[POINT_DIM] = CONSTANT;
206 variationType[D1_DIM] = GENERAL;
207 variationType[D2_DIM] = GENERAL;
208
209 int cellTypeModulus = uniformJacobianModulus();
210
211 data = getMatchingViewWithLabel(pointComponentView, "CellGeometryProvider: Jacobian data", cellTypeModulus, spaceDim, spaceDim);
212 }
213 else
214 {
215 // diagonal Jacobian
216 variationType[D1_DIM] = BLOCK_PLUS_DIAGONAL;
217 variationType[D2_DIM] = BLOCK_PLUS_DIAGONAL;
218 blockPlusDiagonalLastNonDiagonal = -1;
219
220 data = getMatchingViewWithLabel(pointComponentView, "CellGeometryProvider: Jacobian data", spaceDim);
221 }
222 }
223 else if (cellGeometryType_ == TENSOR_GRID)
224 {
225 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "tensor grid support not yet implemented");
226 }
227 else if (cellGeometryType_ == FIRST_ORDER)
228 {
229 const bool simplex = (spaceDim + 1 == cellToNodes_.extent_int(1));
230 if (simplex)
231 {
232 variationType[CELL_DIM] = GENERAL;
233 variationType[POINT_DIM] = CONSTANT; // affine: no point variation
234 variationType[D1_DIM] = GENERAL;
235 variationType[D2_DIM] = GENERAL;
236
237 data = getMatchingViewWithLabel(data, "CellGeometryProvider: Jacobian data", numCells_, spaceDim, spaceDim);
238 }
239 else
240 {
241 variationType[CELL_DIM] = GENERAL;
242 variationType[D1_DIM] = GENERAL;
243 variationType[D2_DIM] = GENERAL;
244 if (affine_)
245 {
246 // no point variation
247 variationType[POINT_DIM] = CONSTANT;
248 data = getMatchingViewWithLabel(data, "CellGeometryProvider: Jacobian data", numCellsWorkset, spaceDim, spaceDim);
249 }
250 else
251 {
252 variationType[POINT_DIM] = GENERAL;
253 data = getMatchingViewWithLabel(data, "CellGeometryProvider: Jacobian data", numCellsWorkset, pointsPerCell, spaceDim, spaceDim);
254 }
255 }
256 }
257 else if (cellGeometryType_ == HIGHER_ORDER)
258 {
259 // most general case: varies in all 4 dimensions
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);
265 }
266 else
267 {
268 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "support for this CellGeometryType is not yet implemented");
269 }
270
271 Data<PointScalar,DeviceType> jacobianData(data,rank,extents,variationType,blockPlusDiagonalLastNonDiagonal);
272 return jacobianData;
273 }
274
275 template<class PointScalar, int spaceDim, typename DeviceType>
277 const int &pointsPerCell, const Data<PointScalar,DeviceType> &refData,
278 const int startCell, const int endCell) const
279 {
280 const int numCellsWorkset = (endCell == -1) ? (numCells_ - startCell) : (endCell - startCell);
281
282 if (cellGeometryType_ == UNIFORM_GRID)
283 {
284 if (uniformJacobianModulus() != 1)
285 {
286 int cellTypeModulus = uniformJacobianModulus();
287
288 auto dataView3 = jacobianData.getUnderlyingView3(); // (cellTypeModulus, spaceDim, spaceDim) allocated in allocateJacobianDataPrivate()
289 auto dataHost = Kokkos::create_mirror_view(dataView3);
290
291 const int startCellType = startCell % cellTypeModulus;
292 const int endCellType = (numCellsWorkset >= cellTypeModulus) ? startCellType + cellTypeModulus : startCellType + numCellsWorkset;
293 const int gridCellOrdinal = 0; // sample cell
294 for (int cellType=startCellType; cellType<endCellType; cellType++)
295 {
296 const int subdivisionOrdinal = cellType % cellTypeModulus;
297 const int nodeZero = 0;
298 // simplex Jacobian formula is J_00 = x1 - x0, J_01 = x2 - x0, etc.
299 for (int i=0; i<spaceDim; i++)
300 {
301 for (int j=0; j<spaceDim; j++)
302 {
303 const int node = j+1; // this is the only node other than the 0 node that has non-zero derivative in the j direction -- and this has unit derivative
304 // nodeZero has derivative -1 in every dimension.
305 const auto J_ij = subdivisionCoordinate(gridCellOrdinal, subdivisionOrdinal, node, i) - subdivisionCoordinate(gridCellOrdinal, subdivisionOrdinal, nodeZero, i);
306 dataHost(cellType,i,j) = J_ij;
307 }
308 }
309 }
310
311 Kokkos::deep_copy(dataView3,dataHost);
312 }
313 else
314 {
315 // diagonal Jacobian
316 auto dataView1 = jacobianData.getUnderlyingView1(); // (spaceDim) allocated in allocateJacobianDataPrivate()
317 const auto domainExtents = domainExtents_;
318 const auto gridCellCounts = gridCellCounts_;
319
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)
323 {
324 // diagonal jacobian
325 const double REF_SPACE_EXTENT = 2.0;
326 dataView1(d1) = (domainExtents[d1] / REF_SPACE_EXTENT) / gridCellCounts[d1];
327 });
328 ExecutionSpace().fence();
329 }
330 }
331 else if (cellGeometryType_ == TENSOR_GRID)
332 {
333 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "tensor grid support not yet implemented");
334 }
335 else if ((cellGeometryType_ == FIRST_ORDER) || (cellGeometryType_ == HIGHER_ORDER))
336 {
337 const bool simplex = (spaceDim + 1 == cellToNodes_.extent_int(1));
338 if (simplex)
339 {
340 auto dataView3 = jacobianData.getUnderlyingView3(); // (numCells_, spaceDim, spaceDim) allocated in allocateJacobianDataPrivate()
341
342 // get local (shallow) copies to avoid implicit references to this
343 auto cellToNodes = cellToNodes_;
344 auto nodes = nodes_;
345
346 using ExecutionSpace = typename DeviceType::execution_space;
347 auto policy = Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<3>>({startCell,0,0},{numCellsWorkset,spaceDim,spaceDim});
348
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; // nodeZero has derivative -1 in every dimension.
352 const int node = d2+1; // this is the only node other than the 0 node that has non-zero derivative in the d2 direction -- and this has unit derivative (except in 1D, where each derivative is ±0.5)
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;
357 });
358 }
359 else
360 {
362 auto basisForNodes = this->basisForNodes();
363
364 if (affine_)
365 {
366 // no point variation
367 auto dataView3 = jacobianData.getUnderlyingView3(); // (numCellsWorkset, spaceDim, spaceDim) allocated in allocateJacobianDataPrivate()
368
369 // TODO: find an allocation-free way to do this… (consider modifying CellTools::setJacobian() to support affine case.)
370 const int onePoint = 1;
371 auto testPointView = getMatchingViewWithLabel(dataView3, "CellGeometryProvider: test point", onePoint, spaceDim);
372 auto tempData = getMatchingViewWithLabel(dataView3, "CellGeometryProvider: temporary Jacobian data", numCellsWorkset, onePoint, spaceDim, spaceDim);
373
374 Kokkos::deep_copy(testPointView, 0.0);
375
376 CellTools::setJacobian(tempData, testPointView, *this, basisForNodes, startCell, endCell);
377
378 auto tempDataSubview = Kokkos::subview(tempData, Kokkos::ALL(), 0, Kokkos::ALL(), Kokkos::ALL());
379 Kokkos::deep_copy(dataView3, tempDataSubview);
380 }
381 else
382 {
383 auto dataView = jacobianData.getUnderlyingView(); // (numCellsWorkset, pointsPerCell, spaceDim, spaceDim) allocated in allocateJacobianDataPrivate()
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");
386
387 // refData should contain the basis gradients; shape is (F,P,D) or (C,F,P,D)
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) && (refData.rank() != 4), std::invalid_argument, "refData should have shape (F,P,D) or (C,F,P,D)");
390 if (refData.rank() == 3)
391 {
392 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(refData.extent_int(0) != basisForNodes->getCardinality(), std::invalid_argument, "refData should have shape (F,P,D) or (C,F,P,D)");
393 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(refData.extent_int(1) != pointsPerCell, std::invalid_argument, "refData should have shape (F,P,D) or (C,F,P,D)");
394 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(refData.extent_int(2) != spaceDim, std::invalid_argument, "refData should have shape (F,P,D) or (C,F,P,D)");
395 }
396 else
397 {
398 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(refData.extent_int(0) != numCellsWorkset, std::invalid_argument, "refData should have shape (F,P,D) or (C,F,P,D)");
399 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(refData.extent_int(1) != basisForNodes->getCardinality(), std::invalid_argument, "refData should have shape (F,P,D) or (C,F,P,D)");
400 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(refData.extent_int(2) != pointsPerCell, std::invalid_argument, "refData should have shape (F,P,D) or (C,F,P,D)");
401 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(refData.extent_int(3) != spaceDim, std::invalid_argument, "refData should have shape (F,P,D) or (C,F,P,D)");
402 }
403
404 CellTools::setJacobian(dataView, *this, refData, startCell, endCell);
405 }
406 }
407 }
408 else
409 {
410 // TODO: handle the other cases
411 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "support for this CellGeometryType is not yet implemented");
412 }
413 }
414
415 // Uniform grid constructor, with optional subdivision into simplices
416 template<class PointScalar, int spaceDim, typename DeviceType>
417 CellGeometry<PointScalar,spaceDim,DeviceType>::CellGeometry(const Kokkos::Array<PointScalar,spaceDim> &origin,
418 const Kokkos::Array<PointScalar,spaceDim> &domainExtents,
419 const Kokkos::Array<int,spaceDim> &gridCellCounts,
420 SubdivisionStrategy subdivisionStrategy,
421 HypercubeNodeOrdering nodeOrdering)
422 :
423 nodeOrdering_(nodeOrdering),
424 cellGeometryType_(UNIFORM_GRID),
425 subdivisionStrategy_(subdivisionStrategy),
426 affine_(true),
427 origin_(origin),
428 domainExtents_(domainExtents),
429 gridCellCounts_(gridCellCounts)
430 {
431 numCells_ = 1;
432 for (int d=0; d<spaceDim; d++)
433 {
434 numCells_ *= gridCellCounts_[d];
435 }
436 numCells_ *= numCellsPerGridCell(subdivisionStrategy_);
437
438 shards::CellTopology cellTopo; // will register with HostMemberLookup below
439 if (subdivisionStrategy_ == NO_SUBDIVISION)
440 {
441 // hypercube
442 numNodesPerCell_ = 1 << spaceDim; // 2^D vertices in a D-dimensional hypercube
443
444 if (spaceDim == 1)
445 {
446 cellTopo = shards::CellTopology(shards::getCellTopologyData<shards::Line<> >());
447 }
448 else if (spaceDim == 2)
449 {
450 cellTopo = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<> >());
451 }
452 else if (spaceDim == 3)
453 {
454 cellTopo = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<> >());
455 }
456 else
457 {
458 // TODO: Once shards supports higher-dimensional hypercubes, initialize cellTopo accordingly
459 }
460 }
461 else
462 {
463 // simplex
464 numNodesPerCell_ = spaceDim + 1; // D+1 vertices in a D-dimensional simplex
465 if (spaceDim == 2)
466 {
467 cellTopo = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<> >());
468 }
469 else if (spaceDim == 3)
470 {
471 cellTopo = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<> >());
472 }
473 else
474 {
475 // TODO: Once shards supports higher-dimensional simplices, initialize cellTopo_ accordingly
476 }
477 }
478
480 const int linearPolyOrder = 1;
481 BasisPtr basisForNodes = getBasis<BasisFamily>(cellTopo, FUNCTION_SPACE_HGRAD, linearPolyOrder);
482
483 if (nodeOrdering_ == HYPERCUBE_NODE_ORDER_CLASSIC_SHARDS)
484 {
485 // override basisForNodes for quad, hexahedron. Apparently the lowest-order bases below are *not* in the same order as their
486 // arbitrary-polynomial-order counterparts; the latter do not match the order of the shards::CellTopology nodes.
487 if (cellTopo.getKey() == shards::Quadrilateral<>::key)
488 {
490 }
491 else if (cellTopo.getKey() == shards::Hexahedron<>::key)
492 {
494 }
495 }
496
498 HostMemberLookup::constructorCalled(this, cellTopo, basisForNodes);
499 }
500
501 // Node-based constructor for straight-edged cell geometry.
502 // If claimAffine is true, we assume (without checking) that the mapping from reference space is affine.
503 // (If claimAffine is false, we check whether the topology is simplicial; if so, we conclude that the mapping must be affine.)
504 template<class PointScalar, int spaceDim, typename DeviceType>
506 ScalarView<int,DeviceType> cellToNodes,
507 ScalarView<PointScalar,DeviceType> nodes,
508 const bool claimAffine,
509 const HypercubeNodeOrdering nodeOrdering)
510 :
511 nodeOrdering_(nodeOrdering),
512 cellGeometryType_(FIRST_ORDER),
513 cellToNodes_(cellToNodes),
514 nodes_(nodes)
515 {
516 if(cellToNodes.is_allocated())
517 {
518 numCells_ = cellToNodes.extent_int(0);
519 numNodesPerCell_ = cellToNodes.extent_int(1);
520 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(numNodesPerCell_ != cellTopo.getNodeCount(), std::invalid_argument, "cellToNodes.extent(1) does not match the cell topology node count");
521 }
522 else
523 {
524 numCells_ = nodes.extent_int(0);
525 numNodesPerCell_ = nodes.extent_int(1);
526 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(numNodesPerCell_ != cellTopo.getNodeCount(), std::invalid_argument, "nodes.extent(1) does not match the cell topology node count");
527 }
528
529
530 if (!claimAffine)
531 {
532 // if cellTopo is simplicial, then since the geometry is straight-edged, it is also affine
533 const bool simplicialTopo = (cellTopo.getNodeCount() == cellTopo.getDimension() + 1);
534 affine_ = simplicialTopo;
535 }
536 else
537 {
538 affine_ = true;
539 }
540
542 const int linearPolyOrder = 1;
543 BasisPtr basisForNodes = getBasis<BasisFamily>(cellTopo, FUNCTION_SPACE_HGRAD, linearPolyOrder);
544
545 if (nodeOrdering_ == HYPERCUBE_NODE_ORDER_CLASSIC_SHARDS)
546 {
547 // override basisForNodes for quad, hexahedron. Apparently the lowest-order bases below are *not* in the same order as their
548 // arbitrary-polynomial-order counterparts; the latter do not match the order of the shards::CellTopology nodes.
549 if (cellTopo.getKey() == shards::Quadrilateral<>::key)
550 {
552 }
553 else if (cellTopo.getKey() == shards::Hexahedron<>::key)
554 {
556 }
557 }
558
560 HostMemberLookup::constructorCalled(this, cellTopo, basisForNodes);
561 }
562
563 // Constructor for higher-order geometry
564 template<class PointScalar, int spaceDim, typename DeviceType>
566 ScalarView<PointScalar,DeviceType> cellNodes)
567 :
568 nodeOrdering_(HYPERCUBE_NODE_ORDER_TENSOR),
569 cellGeometryType_(HIGHER_ORDER),
570 nodes_(cellNodes)
571 {
572 numCells_ = cellNodes.extent_int(0);
573 numNodesPerCell_ = cellNodes.extent_int(1);
574
575 // if basis degree is 1, mark as first-order geometry
576 const bool firstOrderGeometry = (basisForNodes->getDegree() == 1);
577 cellGeometryType_ = firstOrderGeometry ? FIRST_ORDER : HIGHER_ORDER;
578
579 shards::CellTopology cellTopo = basisForNodes->getBaseCellTopology();
580
581 if (firstOrderGeometry && (cellTopo.getNodeCount() == spaceDim + 1)) // lowest-order and simplicial
582 {
583 affine_ = true;
584 }
585 else
586 {
587 affine_ = false;
588 }
590 HostMemberLookup::constructorCalled(this, cellTopo, basisForNodes);
591 }
592
593 template<class PointScalar, int spaceDim, typename DeviceType>
594 KOKKOS_INLINE_FUNCTION
596 {
597 return affine_;
598 }
599
600 template<class PointScalar, int spaceDim, typename DeviceType>
603 const TensorData<PointScalar,DeviceType> & cubatureWeights ) const
604 {
605 // Output possibilities for a cubatureWeights with N components:
606 // 1. For AFFINE elements (jacobianDet cell-wise constant), returns a container with N+1 tensorial components; the first component corresponds to cells
607 // 2. Otherwise, returns a container with 1 tensorial component
608
609 INTREPID2_TEST_FOR_EXCEPTION(cubatureWeights.rank() != 1, std::invalid_argument, "cubatureWeights container must have shape (P)");
610
611 const int numTensorComponents = affine_ ? cubatureWeights.numTensorComponents() + 1 : 1;
612 std::vector< Data<PointScalar,DeviceType> > tensorComponents(numTensorComponents);
613
614 if (affine_)
615 {
616 const int cellExtent = jacobianDet.extent_int(0);
617 Kokkos::Array<DataVariationType,7> cellVariationTypes {jacobianDet.getVariationTypes()[0], CONSTANT, CONSTANT, CONSTANT, CONSTANT, CONSTANT, CONSTANT};
618 const int cellDataDim = jacobianDet.getDataExtent(0);
619 Kokkos::Array<int,7> cellExtents{cellExtent,1,1,1,1,1,1};
620
621 ScalarView<PointScalar,DeviceType> detDataView ("cell relative volumes", cellDataDim);
622 tensorComponents[0] = Data<PointScalar,DeviceType>(detDataView,1,cellExtents,cellVariationTypes);
623
624 for (int cubTensorComponent=0; cubTensorComponent<numTensorComponents-1; cubTensorComponent++)
625 {
626 auto cubatureComponent = cubatureWeights.getTensorComponent(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);
633 }
634 }
635 else
636 {
637 const int cellExtent = jacobianDet.extent_int(0);
638 Kokkos::Array<DataVariationType,7> variationTypes {jacobianDet.getVariationTypes()[0], GENERAL, CONSTANT, CONSTANT, CONSTANT, CONSTANT, CONSTANT};
639 const int cellDataDim = jacobianDet.getDataExtent(0);
640
641 const int numPoints = cubatureWeights.extent_int(0);
642 Kokkos::Array<int,7> extents{cellExtent,numPoints,1,1,1,1,1};
643
644 ScalarView<PointScalar,DeviceType> cubatureWeightView;
645 if (variationTypes[0] != CONSTANT)
646 {
647 cubatureWeightView = ScalarView<PointScalar,DeviceType>("cell measure", cellDataDim, numPoints);
648 }
649 else
650 {
651 cubatureWeightView = ScalarView<PointScalar,DeviceType>("cell measure", numPoints);
652 }
653 const int cellMeasureRank = 2;
654 tensorComponents[0] = Data<PointScalar,DeviceType>(cubatureWeightView,cellMeasureRank,extents,variationTypes);
655 }
656 const bool separateFirstComponent = (numTensorComponents > 1);
657 return TensorData<PointScalar,DeviceType>(tensorComponents, separateFirstComponent);
658 }
659
660 template<class PointScalar, int spaceDim, typename DeviceType>
662 const Data<PointScalar,DeviceType> & jacobianDet,
663 const TensorData<PointScalar,DeviceType> & cubatureWeights ) const
664 {
665 // Output possibilities for a cubatureWeights with N components:
666 // 1. For AFFINE elements (jacobianDet constant on each cell), returns a container with N+1 tensorial components; the first component corresponds to cells
667 // 2. Otherwise, returns a container with 1 tensorial component
668
669 INTREPID2_TEST_FOR_EXCEPTION((cellMeasure.numTensorComponents() != cubatureWeights.numTensorComponents() + 1) && (cellMeasure.numTensorComponents() != 1), std::invalid_argument,
670 "cellMeasure must either have a tensor component count of 1 or a tensor component count that is one higher than that of cubatureWeights");
671
672 INTREPID2_TEST_FOR_EXCEPTION(cubatureWeights.rank() != 1, std::invalid_argument, "cubatureWeights container must have shape (P)");
673
674 if (cellMeasure.numTensorComponents() == cubatureWeights.numTensorComponents() + 1)
675 {
676 // affine case; the first component should contain the cell volume divided by ref cell volume; this should be stored in jacobianDet
677 Kokkos::deep_copy(cellMeasure.getTensorComponent(0).getUnderlyingView1(), jacobianDet.getUnderlyingView1()); // copy point-invariant data from jacobianDet to the first tensor component of cell measure container
678 const int numTensorDimensions = cubatureWeights.numTensorComponents();
679 for (int i=1; i<numTensorDimensions+1; i++)
680 {
681 Kokkos::deep_copy(cellMeasure.getTensorComponent(i).getUnderlyingView1(), cubatureWeights.getTensorComponent(i-1).getUnderlyingView1());
682 }
683 }
684 else
685 {
686 auto detVaries = jacobianDet.getVariationTypes();
687
688 const bool detCellVaries = detVaries[0] != CONSTANT;
689 const bool detPointVaries = detVaries[1] != CONSTANT;
690
691 if (detCellVaries && detPointVaries)
692 {
693 auto cellMeasureData = cellMeasure.getTensorComponent(0).getUnderlyingView2();
694 auto detData = jacobianDet.getUnderlyingView2();
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");
699
700 // We implement this case as a functor (rather than a lambda) to work around an apparent CUDA 10.1.243 compiler bug
701 Impl::CellMeasureFunctor<PointScalar,spaceDim,DeviceType> cellMeasureFunctor(cellMeasureData, detData, cubatureWeights);
702
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);
706 }
707 else if (detCellVaries && !detPointVaries)
708 {
709 auto cellMeasureData = cellMeasure.getTensorComponent(0).getUnderlyingView2();
710 auto detData = jacobianDet.getUnderlyingView1();
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);
716 });
717 }
718 else
719 {
720 // constant jacobian det case
721 // cell measure data has shape (P)
722 auto cellMeasureData = cellMeasure.getTensorComponent(0).getUnderlyingView1();
723 auto detData = jacobianDet.getUnderlyingView1();
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);
728 });
729 }
730 }
731 }
732
733 template<class PointScalar, int spaceDim, typename DeviceType>
734 typename CellGeometry<PointScalar,spaceDim,DeviceType>::BasisPtr
736 {
738 return HostMemberLookup::getBasis(this);
739 }
740
741 template<class PointScalar, int spaceDim, typename DeviceType>
743 {
745 return HostMemberLookup::getCellTopology(this);
746 }
747
748 template<class PointScalar, int spaceDim, typename DeviceType>
749 KOKKOS_INLINE_FUNCTION
751 {
752 if (cellGeometryType_ == UNIFORM_GRID)
753 {
754 const int numSubdivisions = numCellsPerGridCell(subdivisionStrategy_);
755 if (numSubdivisions == 1)
756 {
757 return CONSTANT;
758 }
759 else
760 {
761 return MODULAR;
762 }
763 }
764 else return GENERAL;
765 }
766
767 template<class PointScalar, int spaceDim, typename DeviceType>
769 {
770 Data<PointScalar,DeviceType> emptyRefData;
771 if (cellGeometryType_ == UNIFORM_GRID)
772 {
773 // no need for basis computations
774 return emptyRefData;
775 }
776 else if (cellGeometryType_ == TENSOR_GRID)
777 {
778 // no need for basis values
779 return emptyRefData;
780 }
781 else if ((cellGeometryType_ == FIRST_ORDER) || (cellGeometryType_ == HIGHER_ORDER))
782 {
783 const bool simplex = (spaceDim + 1 == cellToNodes_.extent_int(1));
784 if (simplex)
785 {
786 // no need for precomputed basis values
787 return emptyRefData;
788 }
789 else
790 {
791 auto basisForNodes = this->basisForNodes();
792
793 if (affine_)
794 {
795 // no need for precomputed basis values
796 return emptyRefData;
797 }
798 else
799 {
800 // 2 use cases: (P,D) and (C,P,D).
801 // if (P,D), call the TensorPoints variant
802 if (points.rank() == 2)
803 {
804 TensorPoints<PointScalar,DeviceType> tensorPoints(points);
805 return getJacobianRefData(tensorPoints);
806 }
807 else
808 {
809 const int numCells = points.extent_int(0);
810 const int numPoints = points.extent_int(1);
811 const int numFields = basisForNodes->getCardinality();
812
813 auto cellBasisGradientsView = getMatchingViewWithLabel(points, "CellGeometryProvider: cellBasisGradients", numCells, numFields, numPoints, spaceDim);
814 auto basisGradientsView = getMatchingViewWithLabel(points, "CellGeometryProvider: basisGradients", numFields, numPoints, spaceDim);
815
816 for (int cellOrdinal=0; cellOrdinal<numCells; cellOrdinal++)
817 {
818 auto refPointsForCell = Kokkos::subview(points, cellOrdinal, Kokkos::ALL(), Kokkos::ALL());
819 basisForNodes->getValues(basisGradientsView, refPointsForCell, OPERATOR_GRAD);
820
821 // At some (likely relatively small) memory cost, we copy the BasisGradients into an explicit (F,P,D) container.
822 // Given that we expect to reuse this for a non-trivial number of cell in the common use case, the extra memory
823 // cost is likely worth the increased flop count, etc. (One might want to revisit this in cases of high spaceDim
824 // and/or very high polynomial order.)
825
826 using ExecutionSpace = typename DeviceType::execution_space;
827 auto policy = Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<3>>({0,0,0},{numFields,numPoints,spaceDim});
828
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);
832 });
833 ExecutionSpace().fence();
834 }
835 Data<PointScalar,DeviceType> basisRefData(cellBasisGradientsView);
836 return basisRefData;
837 }
838 }
839 }
840 }
841 else
842 {
843 // TODO: handle the other cases
844 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "support for this CellGeometryType is not yet implemented");
845 }
846 return emptyRefData;
847
848
849 }
850
851 template<class PointScalar, int spaceDim, typename DeviceType>
853 {
854 Data<PointScalar,DeviceType> emptyRefData;
855 if (cellGeometryType_ == UNIFORM_GRID)
856 {
857 // no need for basis computations
858 return emptyRefData;
859 }
860 else if (cellGeometryType_ == TENSOR_GRID)
861 {
862 // no need for basis values
863 return emptyRefData;
864 }
865 else if ((cellGeometryType_ == FIRST_ORDER) || (cellGeometryType_ == HIGHER_ORDER))
866 {
867 const bool simplex = (spaceDim + 1 == cellToNodes_.extent_int(1));
868 if (simplex)
869 {
870 // no need for precomputed basis values
871 return emptyRefData;
872 }
873 else
874 {
875 auto basisForNodes = this->basisForNodes();
876
877 if (affine_)
878 {
879 // no need for precomputed basis values
880 return emptyRefData;
881 }
882 else
883 {
884 auto basisGradients = basisForNodes->allocateBasisValues(points, OPERATOR_GRAD);
885 basisForNodes->getValues(basisGradients, points, OPERATOR_GRAD);
886
887 int numPoints = points.extent_int(0);
888 int numFields = basisForNodes->getCardinality();
889
890 // At some (likely relatively small) memory cost, we copy the BasisGradients into an explicit (F,P,D) container.
891 // Given that we expect to reuse this for a non-trivial number of cell in the common use case, the extra memory
892 // cost is likely worth the increased flop count, etc. (One might want to revisit this in cases of high spaceDim
893 // and/or very high polynomial order.)
894
895 auto firstPointComponentView = points.getTensorComponent(0); // (P,D0)
896 auto basisGradientsView = getMatchingViewWithLabel(firstPointComponentView, "CellGeometryProvider: temporary basisGradients", numFields, numPoints, spaceDim);
897
898 using ExecutionSpace = typename DeviceType::execution_space;
899 auto policy = Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<3>>({0,0,0},{numFields,numPoints,spaceDim});
900
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);
904 });
905 ExecutionSpace().fence();
906
907 Data<PointScalar,DeviceType> basisRefData(basisGradientsView);
908 return basisRefData;
909 }
910 }
911 }
912 else
913 {
914 // TODO: handle the other cases
915 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "support for this CellGeometryType is not yet implemented");
916 }
917 return emptyRefData;
918 }
919
920 template<class PointScalar, int spaceDim, typename DeviceType>
921 KOKKOS_INLINE_FUNCTION
923 {
924 if (nodeOrdering_ == HYPERCUBE_NODE_ORDER_CLASSIC_SHARDS)
925 {
926 // note that Shards numbers nodes for quad counter-clockwise
927 // cube is tensor-product of the (x,y) quad with a z line segment
928 if (d==0)
929 {
930 if ((hypercubeNodeNumber % 4 == 1) || (hypercubeNodeNumber % 4 == 2))
931 return 1;
932 else
933 return 0;
934 }
935 else if (d==1)
936 {
937 if ((hypercubeNodeNumber % 4 == 2) || (hypercubeNodeNumber % 4 == 3))
938 return 1;
939 else
940 return 0;
941 }
942 }
943 // tensor formula coincides with shards formula for d ≥ 2
944 const int nodesForPriorDimensions = 1 << d;
945 if ((hypercubeNodeNumber / nodesForPriorDimensions) % 2 == 1)
946 return 1;
947 else
948 return 0;
949 }
950
951 template<class PointScalar, int spaceDim, typename DeviceType>
953 {
954 using HostExecSpace = Kokkos::DefaultHostExecutionSpace;
955
956 const bool isGridType = (cellGeometryType_ == TENSOR_GRID) || (cellGeometryType_ == UNIFORM_GRID);
957 const int numOrientations = isGridType ? numCellsPerGridCell(subdivisionStrategy_) : numCells();
958
959 const int nodesPerCell = numNodesPerCell();
960
961 ScalarView<Orientation, DeviceType> orientationsView("orientations", numOrientations);
962 auto orientationsHost = Kokkos::create_mirror_view(typename HostExecSpace::memory_space(), orientationsView);
963
964 ScalarView<PointScalar, HostExecSpace> cellNodesHost("cellNodesHost",numOrientations,nodesPerCell); // (C,N) -- where C = numOrientations
965
966 DataVariationType cellVariationType;
967
968 if (isGridType)
969 {
970 // then there are as many distinct orientations possible as there are there are cells per grid cell
971 // fill cellNodesHost with sample nodes from grid cell 0
972 const int numSubdivisions = numCellsPerGridCell(subdivisionStrategy_); // can be up to 6
973
974#if defined(INTREPID2_COMPILE_DEVICE_CODE)
976#else
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;
983 });
984#endif
985 cellVariationType = (numSubdivisions == 1) ? CONSTANT : MODULAR;
986 }
987 else
988 {
989 cellVariationType = GENERAL;
990 auto cellToNodesHost = Kokkos::create_mirror_view_and_copy(typename HostExecSpace::memory_space(), cellToNodes_);
991 }
992
993 OrientationTools<HostExecSpace>::getOrientation(orientationsHost,cellNodesHost,this->cellTopology());
994 Kokkos::deep_copy(orientationsView,orientationsHost);
995
996 const int orientationsRank = 1; // shape (C)
997 const Kokkos::Array<int,7> orientationExtents {static_cast<int>(numCells_),1,1,1,1,1,1};
998 const Kokkos::Array<DataVariationType,7> orientationVariationTypes { cellVariationType, CONSTANT, CONSTANT, CONSTANT, CONSTANT, CONSTANT, CONSTANT};
999 orientations_ = Data<Orientation,DeviceType>(orientationsView, orientationsRank, orientationExtents, orientationVariationTypes);
1000 }
1001
1002 template<class PointScalar, int spaceDim, typename DeviceType>
1003 KOKKOS_INLINE_FUNCTION
1005 if (r == 0)
1006 {
1007 return numCells_;
1008 }
1009 else if (r == 1)
1010 {
1011 return numNodesPerCell_;
1012 }
1013 else if (r == 2)
1014 {
1015 return spaceDim;
1016 }
1017 else
1018 {
1019 return 1;
1020 }
1021 }
1022
1023 template<class PointScalar, int spaceDim, typename DeviceType>
1024 template <typename iType>
1025 KOKKOS_INLINE_FUNCTION
1026 typename std::enable_if<std::is_integral<iType>::value, int>::type
1028 {
1029 return static_cast<int>(extent(r));
1030 }
1031
1032 template<class PointScalar, int spaceDim, typename DeviceType>
1033 KOKKOS_INLINE_FUNCTION
1039
1040 template<class PointScalar, int spaceDim, typename DeviceType>
1041 KOKKOS_INLINE_FUNCTION
1043 {
1044 return numCells_;
1045 }
1046
1047 template<class PointScalar, int spaceDim, typename DeviceType>
1048 KOKKOS_INLINE_FUNCTION
1050 {
1051 if (cellGeometryType_ == UNIFORM_GRID)
1052 {
1053 return gridCellCounts_[dim];
1054 }
1055 else if (cellGeometryType_ == TENSOR_GRID)
1056 {
1057 return tensorVertices_.extent_int(dim);
1058 }
1059 else
1060 {
1061 return -1; // not valid for this cell geometry type
1062 }
1063 }
1064
1065 template<class PointScalar, int spaceDim, typename DeviceType>
1066 KOKKOS_INLINE_FUNCTION
1068 {
1069 return numNodesPerCell_;
1070 }
1071
1072 template<class PointScalar, int spaceDim, typename DeviceType>
1073 KOKKOS_INLINE_FUNCTION
1075 {
1076 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(!orientations_.isValid(), std::invalid_argument, "orientations_ not initialized; call initializeOrientations() first");
1077 return orientations_(cellNumber);
1078 }
1079
1080 template<class PointScalar, int spaceDim, typename DeviceType>
1082 {
1083 if (!orientations_.isValid())
1084 {
1085 initializeOrientations();
1086 }
1087 return orientations_;
1088 }
1089
1090 template<class PointScalar, int spaceDim, typename DeviceType>
1091 KOKKOS_INLINE_FUNCTION
1092 PointScalar CellGeometry<PointScalar,spaceDim,DeviceType>::gridCellCoordinate(const int &gridCellOrdinal, const int &localNodeNumber, const int &dim) const
1093 {
1094 const int componentNode = hypercubeComponentNodeNumber(localNodeNumber, dim);
1095 int cellCountForPriorDimensions = 1;
1096 for (int d=0; d<dim; d++)
1097 {
1098 cellCountForPriorDimensions *= numCellsInDimension(d);
1099 }
1100 const int componentGridCellOrdinal = (gridCellOrdinal / cellCountForPriorDimensions) % numCellsInDimension(dim);
1101 const int vertexOrdinal = componentGridCellOrdinal + componentNode;
1102 if (cellGeometryType_ == UNIFORM_GRID)
1103 {
1104 return origin_[dim] + (vertexOrdinal * domainExtents_[dim]) / gridCellCounts_[dim];
1105 }
1106 else if (cellGeometryType_ == TENSOR_GRID)
1107 {
1108 Kokkos::Array<int,spaceDim> pointOrdinalComponents;
1109 for (int d=0; d<spaceDim; d++)
1110 {
1111 pointOrdinalComponents[d] = 0;
1112 }
1113 pointOrdinalComponents[dim] = vertexOrdinal;
1114 return tensorVertices_(pointOrdinalComponents,dim);
1115 }
1116 else
1117 {
1118 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "Unsupported geometry type");
1119 return 0; // unreachable; here to avoid compiler warnings
1120 }
1121 }
1122
1123 template<class PointScalar, int spaceDim, typename DeviceType>
1124 KOKKOS_INLINE_FUNCTION
1126 {
1127 return 3; // (C,N,D)
1128 }
1129
1130 template<class PointScalar, int spaceDim, typename DeviceType>
1131 KOKKOS_INLINE_FUNCTION
1132 int CellGeometry<PointScalar,spaceDim,DeviceType>::gridCellNodeForSubdivisionNode(const int &gridCellOrdinal, const int &subdivisionOrdinal,
1133 const int &subdivisionNodeNumber) const
1134 {
1135 // TODO: do something to reuse the nodeLookup containers
1136 switch (subdivisionStrategy_)
1137 {
1138 case NO_SUBDIVISION:
1139 return subdivisionNodeNumber;
1140 case TWO_TRIANGLES_RIGHT:
1141 case TWO_TRIANGLES_LEFT:
1142 case FOUR_TRIANGLES:
1143 {
1144 Kokkos::Array<int,3> nodeLookup;
1145 if (subdivisionStrategy_ == TWO_TRIANGLES_RIGHT)
1146 {
1147 if (subdivisionOrdinal == 0)
1148 {
1149 // bottom-right cell: node numbers coincide with quad node numbers
1150 nodeLookup = {0,1,2};
1151 }
1152 else if (subdivisionOrdinal == 1)
1153 {
1154 // node 0 --> node 2
1155 // node 1 --> node 3
1156 // node 2 --> node 0
1157 nodeLookup = {2,3,0};
1158 }
1159 else
1160 {
1161 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "Unsupported subdivision ordinal");
1162 }
1163 }
1164 else if (subdivisionStrategy_ == TWO_TRIANGLES_LEFT)
1165 {
1166 if (subdivisionOrdinal == 0)
1167 {
1168 // bottom-left cell:
1169 // node 0 --> node 3
1170 // node 1 --> node 0
1171 // node 2 --> node 1
1172 nodeLookup = {3,0,1};
1173 }
1174 else if (subdivisionOrdinal == 1)
1175 {
1176 // node 0 --> node 2
1177 // node 1 --> node 3
1178 // node 2 --> node 0
1179 nodeLookup = {2,3,0};
1180 }
1181 else
1182 {
1183 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "Unsupported subdivision ordinal");
1184 }
1185 }
1186 else // FOUR_TRIANGLES
1187 {
1188 // counter-clockwise, bottom triangle first
1189 // bottom triangle goes:
1190 // 0 --> 1
1191 // 1 --> center
1192 // 2 --> 0
1193 // and similarly for the other triangles, proceeding counter-clockwise
1194 // node 1 always being the center, we special-case that
1195 if (subdivisionNodeNumber == 1)
1196 {
1197 // center coordinate: we call this node 4 in the quadrilateral
1198 return 4;
1199 }
1200 else
1201 {
1202 nodeLookup = {(subdivisionOrdinal + 1) % 4, -1, subdivisionOrdinal};
1203 }
1204 }
1205 const int gridCellNodeNumber = nodeLookup[subdivisionNodeNumber];
1206 return gridCellNodeNumber;
1207 }
1208 case FIVE_TETRAHEDRA:
1209 case SIX_TETRAHEDRA:
1210 {
1211 Kokkos::Array<int,4> nodeLookup;
1212 if (subdivisionStrategy_ == FIVE_TETRAHEDRA)
1213 {
1214 /*
1215 // to discretize a unit cube into 5 tetrahedra, we can take the four vertices
1216 // (1,1,1)
1217 // (0,0,1)
1218 // (0,1,0)
1219 // (1,0,0)
1220 // as an interior tetrahedron. Call this cell 0. The remaining 4 cells can be determined
1221 // by selecting three of the above points (there are exactly 4 such combinations) and then selecting
1222 // from the remaining four vertices of the cube the one nearest the plane defined by those three points.
1223 // The remaining four vertices are:
1224 // (0,0,0)
1225 // (1,1,0)
1226 // (1,0,1)
1227 // (0,1,1)
1228 // For each of these four, we pick one, and then take the three nearest vertices from cell 0 to form a new tetrahedron.
1229 // We enumerate as follows:
1230 // cell 0: (1,1,1), (0,0,1), (0,1,0), (1,0,0)
1231 // cell 1: (0,0,0), (1,0,0), (0,1,0), (0,0,1)
1232 // cell 2: (1,1,0), (1,1,1), (0,1,0), (1,0,0)
1233 // cell 3: (1,0,1), (1,1,1), (0,0,1), (1,0,0)
1234 // cell 4: (0,1,1), (1,1,1), (0,0,1), (0,1,0)
1235 */
1236 // tetrahedra are as follows:
1237 // 0: {1,3,4,6}
1238 // 1: {0,1,3,4}
1239 // 2: {1,2,3,6}
1240 // 3: {1,4,5,6}
1241 // 4: {3,4,6,7}
1242 switch (subdivisionOrdinal) {
1243 case 0:
1244 nodeLookup = {1,3,4,6};
1245 break;
1246 case 1:
1247 nodeLookup = {0,1,3,4};
1248 break;
1249 case 2:
1250 nodeLookup = {1,2,3,6};
1251 break;
1252 case 3:
1253 nodeLookup = {1,4,5,6};
1254 break;
1255 case 4:
1256 nodeLookup = {3,4,6,7};
1257 break;
1258 default:
1259 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "invalid subdivisionOrdinal");
1260 break;
1261 }
1262 }
1263 else if (subdivisionStrategy_ == SIX_TETRAHEDRA)
1264 {
1265 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "support for SIX_TETRAHEDRA not yet implemented");
1266 }
1267 const int gridCellNodeNumber = nodeLookup[subdivisionNodeNumber];
1268 return gridCellNodeNumber;
1269 }
1270 default:
1271 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "Subdivision strategy not yet implemented!");
1272 // some compilers complain about missing return
1273 return 0; // statement should be unreachable...
1274 }
1275 }
1276
1277 template<class PointScalar, int spaceDim, typename DeviceType>
1278 KOKKOS_INLINE_FUNCTION
1279 PointScalar CellGeometry<PointScalar,spaceDim,DeviceType>::subdivisionCoordinate(const int &gridCellOrdinal, const int &subdivisionOrdinal,
1280 const int &subdivisionNodeNumber, const int &d) const
1281 {
1282 int gridCellNode = gridCellNodeForSubdivisionNode(gridCellOrdinal, subdivisionOrdinal, subdivisionNodeNumber);
1283
1284 if (subdivisionStrategy_ == FOUR_TRIANGLES)
1285 {
1286 // this is the one case in which the gridCellNode may not actually be a node in the grid cell
1287 if (gridCellNode == 4) // center vertex
1288 {
1289 // d == 0 means quad vertices 0 and 1 suffice;
1290 // d == 1 means quad vertices 0 and 3 suffice
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));
1294 }
1295 }
1296 return gridCellCoordinate(gridCellOrdinal, gridCellNode, d);
1297 }
1298
1299 template<class PointScalar, int spaceDim, typename DeviceType>
1300 KOKKOS_INLINE_FUNCTION
1301 PointScalar
1302 CellGeometry<PointScalar,spaceDim,DeviceType>::operator()(const int& cell, const int& node, const int& dim) const {
1303 if ((cellGeometryType_ == UNIFORM_GRID) || (cellGeometryType_ == TENSOR_GRID))
1304 {
1305 const int numSubdivisions = numCellsPerGridCell(subdivisionStrategy_);
1306 if (numSubdivisions == 1)
1307 {
1308 // hypercube
1309 return gridCellCoordinate(cell, node, dim);
1310 }
1311 else
1312 {
1313 const int subdivisionOrdinal = cell % numSubdivisions;
1314 const int gridCellOrdinal = cell / numSubdivisions;
1315 return subdivisionCoordinate(gridCellOrdinal, subdivisionOrdinal, node, dim);
1316 }
1317 }
1318 else
1319 {
1320#ifdef HAVE_INTREPID2_DEBUG
1321 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE((cell < 0), std::invalid_argument, "cell out of bounds");
1322 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(static_cast<unsigned>(cell) > numCells_, std::invalid_argument, "cell out of bounds");
1323 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE((node < 0), std::invalid_argument, "node out of bounds");
1324 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(static_cast<unsigned>(node) > numNodesPerCell_, std::invalid_argument, "node out of bounds");
1325 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE((dim < 0), std::invalid_argument, "dim out of bounds" );
1326 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(dim > spaceDim, std::invalid_argument, "dim out of bounds" );
1327#endif
1328 if (cellToNodes_.is_allocated())
1329 {
1330 const int nodeNumber = cellToNodes_(cell,node);
1331 return nodes_(nodeNumber,dim);
1332 }
1333 else
1334 {
1335 return nodes_(cell,node,dim);
1336 }
1337 }
1338 }
1339
1340 template<class PointScalar, int spaceDim, typename DeviceType>
1341 KOKKOS_INLINE_FUNCTION
1343 {
1344 if (cellGeometryType_ == UNIFORM_GRID)
1345 {
1346 return numCellsPerGridCell(subdivisionStrategy_);
1347 }
1348 else
1349 {
1350 return numCells_;
1351 }
1352 }
1353
1354 template<class PointScalar, int spaceDim, typename DeviceType>
1356 {
1357 const int pointsPerCell = points.extent_int(0);
1358 return allocateJacobianDataPrivate(points.getTensorComponent(0),pointsPerCell,startCell,endCell);
1359 }
1360
1361 template<class PointScalar, int spaceDim, typename DeviceType>
1362 Data<PointScalar,DeviceType> CellGeometry<PointScalar,spaceDim,DeviceType>::allocateJacobianData(const ScalarView<PointScalar,DeviceType> &points, const int startCell, const int endCell) const
1363 {
1364 // if points is rank 3, it has shape (C,P,D). If it's rank 2, (P,D).
1365 const int pointDimension = (points.rank() == 3) ? 1 : 0;
1366 const int pointsPerCell = points.extent_int(pointDimension);
1367 return allocateJacobianDataPrivate(points,pointsPerCell,startCell,endCell);
1368 }
1369
1370 template<class PointScalar, int spaceDim, typename DeviceType>
1371 Data<PointScalar,DeviceType> CellGeometry<PointScalar,spaceDim,DeviceType>::allocateJacobianData(const int &numPoints, const int startCell, const int endCell) const
1372 {
1373 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(!affine_, std::invalid_argument, "this version of allocateJacobianData() is only supported for affine CellGeometry");
1374
1375 ScalarView<PointScalar,DeviceType> emptyPoints;
1376 return allocateJacobianDataPrivate(emptyPoints,numPoints,startCell,endCell);
1377 }
1378
1379 template<class PointScalar, int spaceDim, typename DeviceType>
1381 const Data<PointScalar,DeviceType> &refData, const int startCell, const int endCell) const
1382 {
1383 const int pointsPerCell = points.extent_int(0);
1384 setJacobianDataPrivate(jacobianData,pointsPerCell,refData,startCell,endCell);
1385 }
1386
1387 template<class PointScalar, int spaceDim, typename DeviceType>
1388 void CellGeometry<PointScalar,spaceDim,DeviceType>::setJacobian(Data<PointScalar,DeviceType> &jacobianData, const ScalarView<PointScalar,DeviceType> &points,
1389 const Data<PointScalar,DeviceType> &refData, const int startCell, const int endCell) const
1390 {
1391 // if points is rank 3, it has shape (C,P,D). If it's rank 2, (P,D).
1392 const int pointDimension = (points.rank() == 3) ? 1 : 0;
1393 const int pointsPerCell = points.extent_int(pointDimension);
1394 setJacobianDataPrivate(jacobianData,pointsPerCell,refData,startCell,endCell);
1395 }
1396
1397 template<class PointScalar, int spaceDim, typename DeviceType>
1398 void CellGeometry<PointScalar,spaceDim,DeviceType>::setJacobian(Data<PointScalar,DeviceType> &jacobianData, const int &numPoints, const int startCell, const int endCell) const
1399 {
1400 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(!affine_, std::invalid_argument, "this version of setJacobian() is only supported for affine CellGeometry");
1401
1402 Data<PointScalar,DeviceType> emptyRefData;
1403 setJacobianDataPrivate(jacobianData,numPoints,emptyRefData,startCell,endCell);
1404 }
1405} // namespace Intrepid2
1406
1407#endif /* Intrepid2_CellGeometryDef_h */
@ CONSTANT
does not vary
@ 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
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...
A stateless class for operations on cell data. Provides methods for:
static void setJacobian(Kokkos::DynRankView< jacobianValueType, jacobianProperties... > jacobian, const Kokkos::DynRankView< pointValueType, pointProperties... > points, const WorksetType worksetCell, const Teuchos::RCP< HGradBasisType > basis, const int startCell=0, const int endCell=-1)
Computes the Jacobian matrix DF of the reference-to-physical frame map F.
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...
static void getOrientation(Kokkos::DynRankView< elemOrtValueType, elemOrtProperties... > elemOrts, const Kokkos::DynRankView< elemNodeValueType, elemNodeProperties... > elemNodes, const shards::CellTopology cellTopo, bool isSide=false)
Compute orientations of cells in a workset.
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.