48 #ifndef __INTREPID2_ORIENTATIONTOOLS_DEF_MODIFY_BASIS_HPP__
49 #define __INTREPID2_ORIENTATIONTOOLS_DEF_MODIFY_BASIS_HPP__
52 #if defined (__clang__) && !defined (__INTEL_COMPILER)
53 #pragma clang system_header
61 template<
typename elemOrtValueType,
class ...elemOrtProperties,
62 typename elemNodeValueType,
class ...elemNodeProperties>
65 getOrientation( Kokkos::DynRankView<elemOrtValueType,elemOrtProperties...> elemOrts,
66 const Kokkos::DynRankView<elemNodeValueType,elemNodeProperties...> elemNodes,
67 const shards::CellTopology cellTopo,
70 auto elemOrtsHost = Kokkos::create_mirror_view(elemOrts);
71 auto elemNodesHost = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), elemNodes);
73 const ordinal_type numCells = elemNodes.extent(0);
74 for (
auto cell=0;cell<numCells;++cell) {
75 const auto nodes = Kokkos::subview(elemNodesHost, cell, Kokkos::ALL());
76 elemOrtsHost(cell) = Orientation::getOrientation(cellTopo, nodes, isSide);
79 Kokkos::deep_copy(elemOrts, elemOrtsHost);
82 template<
typename ortViewType,
83 typename OutputViewType,
84 typename inputViewType,
87 typename dataViewType>
90 OutputViewType output;
92 o2tViewType ordinalToTag;
93 t2oViewType tagToOrdinal;
95 const dataViewType matData;
96 const ordinal_type cellDim, numVerts, numEdges, numFaces, numPoints, dimBasis;
97 const bool leftMultiply;
100 const bool transpose;
103 OutputViewType output_,
104 inputViewType input_,
105 o2tViewType ordinalToTag_,
106 t2oViewType tagToOrdinal_,
107 const dataViewType matData_,
108 const ordinal_type cellDim_,
109 const ordinal_type numVerts_,
110 const ordinal_type numEdges_,
111 const ordinal_type numFaces_,
112 const ordinal_type numPoints_,
113 const ordinal_type dimBasis_,
114 const bool leftMultiply_ =
true,
115 const bool transpose_ =
false)
119 ordinalToTag(ordinalToTag_),
120 tagToOrdinal(tagToOrdinal_),
126 numPoints(numPoints_),
128 leftMultiply(leftMultiply_),
129 transpose(transpose_)
132 KOKKOS_INLINE_FUNCTION
133 void operator()(
const ordinal_type cell)
const {
134 typedef typename inputViewType::non_const_value_type input_value_type;
136 auto out = Kokkos::subview(output, cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
137 auto in = (input.rank() == output.rank()) ?
138 Kokkos::subview(input, cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL())
139 : Kokkos::subview(input, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
142 ordinal_type existEdgeDofs = 0;
144 ordinal_type ortEdges[12];
145 orts(cell).getEdgeOrientation(ortEdges, numEdges);
148 for (ordinal_type edgeId=0;edgeId<numEdges;++edgeId) {
149 const ordinal_type ordEdge = (1 < tagToOrdinal.extent(0) ? (
static_cast<size_type
>(edgeId) < tagToOrdinal.extent(1) ? tagToOrdinal(1, edgeId, 0) : -1) : -1);
153 const ordinal_type ndofEdge = ordinalToTag(ordEdge, 3);
154 const auto mat = Kokkos::subview(matData,
155 edgeId, ortEdges[edgeId],
156 Kokkos::ALL(), Kokkos::ALL());
158 for (ordinal_type j=0;j<numPoints;++j)
159 for (ordinal_type i=0;i<ndofEdge;++i) {
160 const ordinal_type ii = tagToOrdinal(1, edgeId, i);
162 for (ordinal_type k=0;k<dimBasis;++k) {
163 input_value_type temp = 0.0;
164 for (ordinal_type l=0;l<ndofEdge;++l) {
165 const ordinal_type ll = tagToOrdinal(1, edgeId, l);
166 auto & input_ = leftMultiply ? in(ll, j, k) : in(j, ll, k);
167 auto & mat_il = transpose ? mat(l,i) : mat(i,l);
168 temp += mat_il*input_;
170 auto & output_ = leftMultiply ? out(ii, j, k) : out(j, ii, k);
180 ordinal_type ortFaces[12];
181 orts(cell).getFaceOrientation(ortFaces, numFaces);
184 for (ordinal_type faceId=0;faceId<numFaces;++faceId) {
185 const ordinal_type ordFace = (2 < tagToOrdinal.extent(0) ? (
static_cast<size_type
>(faceId) < tagToOrdinal.extent(1) ? tagToOrdinal(2, faceId, 0) : -1) : -1);
188 const ordinal_type ndofFace = ordinalToTag(ordFace, 3);
189 const auto mat = Kokkos::subview(matData,
190 numEdges*existEdgeDofs+faceId, ortFaces[faceId],
191 Kokkos::ALL(), Kokkos::ALL());
193 for (ordinal_type j=0;j<numPoints;++j)
194 for (ordinal_type i=0;i<ndofFace;++i) {
195 const ordinal_type ii = tagToOrdinal(2, faceId, i);
197 for (ordinal_type k=0;k<dimBasis;++k) {
198 input_value_type temp = 0.0;
199 for (ordinal_type l=0;l<ndofFace;++l) {
200 const ordinal_type ll = tagToOrdinal(2, faceId, l);
201 auto & input_ = leftMultiply ? in(ll, j, k) : in(j, ll, k);
202 auto & mat_il = transpose ? mat(l,i) : mat(i,l);
203 temp += mat_il*input_;
206 auto & output_ = leftMultiply ? out(ii, j, k) : out(j, ii, k);
215 ordinal_type faceOrt(0), edgeOrt(0);
216 if(cellDim == 2) orts(cell).getFaceOrientation(&faceOrt, 1);
218 const ordinal_type ordFace = (2 < tagToOrdinal.extent(0) ? (
static_cast<size_type
>(0) < tagToOrdinal.extent(1) ? tagToOrdinal(2, 0, 0) : -1) : -1);
221 const ordinal_type ndofFace = ordinalToTag(ordFace, 3);
222 const auto mat = Kokkos::subview(matData,
223 numEdges*existEdgeDofs, faceOrt,
224 Kokkos::ALL(), Kokkos::ALL());
226 for (ordinal_type j=0;j<numPoints;++j)
227 for (ordinal_type i=0;i<ndofFace;++i) {
228 const ordinal_type ii = tagToOrdinal(2, 0, i);
230 for (ordinal_type k=0;k<dimBasis;++k) {
231 input_value_type temp = 0.0;
232 for (ordinal_type l=0;l<ndofFace;++l) {
233 const ordinal_type ll = tagToOrdinal(2, 0, l);
234 auto & input_ = leftMultiply ? in(ll, j, k) : in(j, ll, k);
235 auto & mat_il = transpose ? mat(l,i) : mat(i,l);
236 temp += mat_il*input_;
238 auto & output_ = leftMultiply ? out(ii, j, k) : out(j, ii, k);
245 if(cellDim == 1) orts(cell).getEdgeOrientation(&edgeOrt, 1);
247 const ordinal_type ordEdge = (1 < tagToOrdinal.extent(0) ? (
static_cast<size_type
>(0) < tagToOrdinal.extent(1) ? tagToOrdinal(1, 0, 0) : -1) : -1);
250 const ordinal_type ndofEdge = ordinalToTag(ordEdge, 3);
251 const auto mat = Kokkos::subview(matData,
253 Kokkos::ALL(), Kokkos::ALL());
255 for (ordinal_type j=0;j<numPoints;++j)
256 for (ordinal_type i=0;i<ndofEdge;++i) {
257 const ordinal_type ii = tagToOrdinal(1, 0, i);
259 for (ordinal_type k=0;k<dimBasis;++k) {
260 input_value_type temp = 0.0;
261 for (ordinal_type l=0;l<ndofEdge;++l) {
262 const ordinal_type ll = tagToOrdinal(1, 0, l);
263 auto & input_ = leftMultiply ? in(ll, j, k) : in(j, ll, k);
264 auto & mat_il = transpose ? mat(l,i) : mat(i,l);
265 temp += mat_il*input_;
267 auto & output_ = leftMultiply ? out(ii, j, k) : out(j, ii, k);
276 template<
typename DT>
277 template<
typename outputValueType,
class ...outputProperties,
278 typename inputValueType,
class ...inputProperties,
279 typename OrientationViewType,
284 const Kokkos::DynRankView<inputValueType, inputProperties...> input,
285 const OrientationViewType orts,
286 const BasisType* basis,
287 const bool transpose) {
288 #ifdef HAVE_INTREPID2_DEBUG
290 if (input.rank() == output.rank())
292 for (size_type i=0;i<input.rank();++i)
293 INTREPID2_TEST_FOR_EXCEPTION( input.extent(i) != output.extent(i), std::invalid_argument,
294 ">>> ERROR (OrientationTools::modifyBasisByOrientation): Input and output dimensions do not match.");
296 else if (input.rank() == output.rank() - 1)
298 for (size_type i=0;i<input.rank();++i)
299 INTREPID2_TEST_FOR_EXCEPTION( input.extent(i) != output.extent(i+1), std::invalid_argument,
300 ">>> ERROR (OrientationTools::modifyBasisByOrientation): Input dimensions must match output dimensions exactly, or else match all but the first dimension (in the case that input does not have a 'cell' dimension).");
304 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
305 ">>> ERROR (OrientationTools::modifyBasisByOrientation): input and output ranks must either match, or input rank must be one less than that of output.")
308 INTREPID2_TEST_FOR_EXCEPTION( static_cast<ordinal_type>(output.extent(1)) != basis->getCardinality(), std::invalid_argument,
309 ">>> ERROR (OrientationTools::modifyBasisByOrientation): Field dimension of input/output does not match to basis cardinality.");
313 const shards::CellTopology cellTopo = basis->getBaseCellTopology();
314 const ordinal_type cellDim = cellTopo.getDimension();
317 if(input.rank() == output.rank())
318 Kokkos::deep_copy(output, input);
322 if ((cellDim < 3) || basis->requireOrientation()) {
323 auto ordinalToTag = Kokkos::create_mirror_view_and_copy(
typename DT::memory_space(), basis->getAllDofTags());
324 auto tagToOrdinal = Kokkos::create_mirror_view_and_copy(
typename DT::memory_space(), basis->getAllDofOrdinal());
327 numCells = output.extent(0),
329 numPoints = output.extent(2),
330 dimBasis = output.extent(3);
334 ordinal_type numVerts(0), numEdges(0), numFaces(0);
336 if (basis->requireOrientation()) {
337 numVerts = cellTopo.getVertexCount()*ordinal_type(basis->getDofCount(0, 0) > 0);
338 numEdges = cellTopo.getEdgeCount()*ordinal_type(basis->getDofCount(1, 0) > 0);
339 numFaces = cellTopo.getFaceCount()*ordinal_type(basis->getDofCount(2, 0) > 0);
342 bool leftMultiply =
true;
344 const Kokkos::RangePolicy<typename DT::execution_space> policy(0, numCells);
347 decltype(output),decltype(input),
348 decltype(ordinalToTag),decltype(tagToOrdinal),
349 decltype(matData)> FunctorType;
354 ordinalToTag, tagToOrdinal,
356 cellDim, numVerts, numEdges, numFaces,
357 numPoints, dimBasis, leftMultiply, transpose));
361 template<
typename DT>
362 template<
typename outputValueType,
class ...outputProperties,
363 typename inputValueType,
class ...inputProperties,
364 typename OrientationViewType,
369 const Kokkos::DynRankView<inputValueType, inputProperties...> input,
370 const OrientationViewType orts,
371 const BasisType* basis ) {
372 bool transpose =
true;
373 modifyBasisByOrientation(output, input, orts, basis, transpose);
376 template<
typename DT>
377 template<
typename outputValueType,
class ...outputProperties,
378 typename inputValueType,
class ...inputProperties,
379 typename OrientationViewType,
384 const Kokkos::DynRankView<inputValueType, inputProperties...> input,
385 const OrientationViewType orts,
386 const BasisType* basis,
387 const bool transpose ) {
388 #ifdef HAVE_INTREPID2_DEBUG
390 if (input.rank() == output.rank())
392 for (size_type i=0;i<input.rank();++i)
393 INTREPID2_TEST_FOR_EXCEPTION( input.extent(i) != output.extent(i), std::invalid_argument,
394 ">>> ERROR (OrientationTools::modifyBasisByOrientation): Input and output dimensions do not match.");
396 else if (input.rank() == output.rank() - 1)
398 for (size_type i=0;i<input.rank();++i)
399 INTREPID2_TEST_FOR_EXCEPTION( input.extent(i) != output.extent(i+1), std::invalid_argument,
400 ">>> ERROR (OrientationTools::modifyBasisByOrientation): Input dimensions must match output dimensions exactly, or else match all but the first dimension (in the case that input does not have a 'cell' dimension).");
404 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
405 ">>> ERROR (OrientationTools::modifyBasisByOrientation): input and output ranks must either match, or input rank must be one less than that of output.")
408 INTREPID2_TEST_FOR_EXCEPTION( static_cast<ordinal_type>(output.extent(1)) != basis->getCardinality(), std::invalid_argument,
409 ">>> ERROR (OrientationTools::modifyBasisByOrientation): Field dimension of input/output does not match to basis cardinality.");
413 const shards::CellTopology cellTopo = basis->getBaseCellTopology();
414 const ordinal_type cellDim = cellTopo.getDimension();
417 if(input.rank() == output.rank())
418 Kokkos::deep_copy(output, input);
422 if ((cellDim < 3) || basis->requireOrientation()) {
423 auto ordinalToTag = Kokkos::create_mirror_view_and_copy(
typename DT::memory_space(), basis->getAllDofTags());
424 auto tagToOrdinal = Kokkos::create_mirror_view_and_copy(
typename DT::memory_space(), basis->getAllDofOrdinal());
427 numCells = output.extent(0),
429 numPoints = output.extent(2),
430 dimBasis = output.extent(3);
434 ordinal_type numVerts(0), numEdges(0), numFaces(0);
436 if (basis->requireOrientation()) {
437 numVerts = cellTopo.getVertexCount()*ordinal_type(basis->getDofCount(0, 0) > 0);
438 numEdges = cellTopo.getEdgeCount()*ordinal_type(basis->getDofCount(1, 0) > 0);
439 numFaces = cellTopo.getFaceCount()*ordinal_type(basis->getDofCount(2, 0) > 0);
442 bool leftMultiply =
true;
444 const Kokkos::RangePolicy<typename DT::execution_space> policy(0, numCells);
447 decltype(output),decltype(input),
448 decltype(ordinalToTag),decltype(tagToOrdinal),
449 decltype(matData)> FunctorType;
454 ordinalToTag, tagToOrdinal,
456 cellDim, numVerts, numEdges, numFaces,
457 numPoints, dimBasis, leftMultiply, transpose));
461 template<
typename DT>
462 template<
typename outputValueType,
class ...outputProperties,
463 typename inputValueType,
class ...inputProperties,
464 typename OrientationViewType,
465 typename BasisTypeLeft,
466 typename BasisTypeRight>
470 const Kokkos::DynRankView<inputValueType, inputProperties...> input,
471 const OrientationViewType orts,
472 const BasisTypeLeft* basisLeft,
473 const BasisTypeRight* basisRight)
475 const ordinal_type numCells = output.extent(0);
476 const ordinal_type numFieldsLeft = basisLeft->getCardinality();
477 const ordinal_type numFieldsRight = basisRight->getCardinality();
478 #ifdef HAVE_INTREPID2_DEBUG
480 if (input.rank() == output.rank())
482 for (size_type i=0;i<input.rank();++i)
483 INTREPID2_TEST_FOR_EXCEPTION( input.extent(i) != output.extent(i), std::invalid_argument,
484 ">>> ERROR (OrientationTools::modifyMatrixByOrientation): Input and output dimensions do not match.");
486 else if (input.rank() == output.rank() - 1)
488 for (size_type i=0;i<input.rank();++i)
489 INTREPID2_TEST_FOR_EXCEPTION( input.extent(i) != output.extent(i+1), std::invalid_argument,
490 ">>> ERROR (OrientationTools::modifyMatrixByOrientation): Input dimensions must match output dimensions exactly, or else match all but the first dimension (in the case that input does not have a 'cell' dimension).");
494 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
495 ">>> ERROR (OrientationTools::modifyMatrixByOrientation): input and output ranks must either match, or input rank must be one less than that of output.")
498 INTREPID2_TEST_FOR_EXCEPTION( static_cast<ordinal_type>(output.extent(1)) != numFieldsLeft, std::invalid_argument,
499 ">>> ERROR (OrientationTools::modifyMatrixByOrientation): First field dimension of input/output does not match left basis cardinality.");
500 INTREPID2_TEST_FOR_EXCEPTION( static_cast<ordinal_type>(output.extent(2)) != numFieldsRight, std::invalid_argument,
501 ">>> ERROR (OrientationTools::modifyMatrixByOrientation): Second field dimension of input/output does not match right basis cardinality.");
502 INTREPID2_TEST_FOR_EXCEPTION( static_cast<ordinal_type>(output.extent(3)) != 1, std::invalid_argument,
503 ">>> ERROR (OrientationTools::modifyMatrixByOrientation): Third dimension of output must be 1.");
507 const shards::CellTopology cellTopo = basisLeft->getBaseCellTopology();
508 const ordinal_type cellDim = cellTopo.getDimension();
511 decltype(output) outputLeft(
"temp view - output from left application", numCells, numFieldsLeft, numFieldsRight);
514 if(input.rank() == output.rank())
515 Kokkos::deep_copy(outputLeft, input);
519 if ((cellDim < 3) || basisLeft->requireOrientation()) {
520 bool leftMultiply =
true;
521 auto ordinalToTag = Kokkos::create_mirror_view_and_copy(
typename DT::memory_space(), basisLeft->getAllDofTags());
522 auto tagToOrdinal = Kokkos::create_mirror_view_and_copy(
typename DT::memory_space(), basisLeft->getAllDofOrdinal());
525 numOtherFields = output.extent(2),
526 dimBasis = output.extent(3);
530 ordinal_type numVerts(0), numEdges(0), numFaces(0);
532 if (basisLeft->requireOrientation()) {
533 numVerts = cellTopo.getVertexCount()*ordinal_type(basisLeft->getDofCount(0, 0) > 0);
534 numEdges = cellTopo.getEdgeCount()*ordinal_type(basisLeft->getDofCount(1, 0) > 0);
535 numFaces = cellTopo.getFaceCount()*ordinal_type(basisLeft->getDofCount(2, 0) > 0);
538 const Kokkos::RangePolicy<typename DT::execution_space> policy(0, numCells);
541 decltype(outputLeft),decltype(input),
542 decltype(ordinalToTag),decltype(tagToOrdinal),
543 decltype(matData)> FunctorType;
548 ordinalToTag, tagToOrdinal,
550 cellDim, numVerts, numEdges, numFaces,
551 numOtherFields, dimBasis, leftMultiply));
556 Kokkos::deep_copy(output, outputLeft);
557 if ((cellDim < 3) || basisRight->requireOrientation()) {
558 bool leftMultiply =
false;
559 auto ordinalToTag = Kokkos::create_mirror_view_and_copy(
typename DT::memory_space(), basisRight->getAllDofTags());
560 auto tagToOrdinal = Kokkos::create_mirror_view_and_copy(
typename DT::memory_space(), basisRight->getAllDofOrdinal());
563 numOtherFields = output.extent(1),
564 dimBasis = output.extent(3);
568 ordinal_type numVerts(0), numEdges(0), numFaces(0);
570 if (basisRight->requireOrientation()) {
571 numVerts = cellTopo.getVertexCount()*ordinal_type(basisRight->getDofCount(0, 0) > 0);
572 numEdges = cellTopo.getEdgeCount()*ordinal_type(basisRight->getDofCount(1, 0) > 0);
573 numFaces = cellTopo.getFaceCount()*ordinal_type(basisRight->getDofCount(2, 0) > 0);
576 const Kokkos::RangePolicy<typename DT::execution_space> policy(0, numCells);
579 decltype(output),decltype(outputLeft),
580 decltype(ordinalToTag),decltype(tagToOrdinal),
581 decltype(matData)> FunctorType;
586 ordinalToTag, tagToOrdinal,
588 cellDim, numVerts, numEdges, numFaces,
589 numOtherFields, dimBasis, leftMultiply));
Header file for the Intrepid2::Orientation class.