49#ifndef Intrepid2_TensorBasis_h
50#define Intrepid2_TensorBasis_h
52#include <Kokkos_DynRankView.hpp>
54#include <Teuchos_RCP.hpp>
56#include <Intrepid2_config.h>
72 template<ordinal_type spaceDim>
73 KOKKOS_INLINE_FUNCTION
76 template<ordinal_type spaceDim>
77 KOKKOS_INLINE_FUNCTION
78 void getDkEnumerationInverse(Kokkos::Array<int,spaceDim> &entries,
const ordinal_type dkEnum,
const ordinal_type operatorOrder);
81 KOKKOS_INLINE_FUNCTION
82 void getDkEnumerationInverse<1>(Kokkos::Array<int,1> &entries,
const ordinal_type dkEnum,
const ordinal_type operatorOrder)
84 entries[0] = operatorOrder;
88 KOKKOS_INLINE_FUNCTION
89 void getDkEnumerationInverse<2>(Kokkos::Array<int,2> &entries,
const ordinal_type dkEnum,
const ordinal_type operatorOrder)
91 entries[0] = operatorOrder - dkEnum;
96 KOKKOS_INLINE_FUNCTION
97 void getDkEnumerationInverse<3>(Kokkos::Array<int,3> &entries,
const ordinal_type dkEnum,
const ordinal_type operatorOrder)
102 for (ordinal_type yMult=0; yMult<=operatorOrder; yMult++)
104 for (ordinal_type zMult=0; zMult<=operatorOrder-yMult; zMult++)
106 const ordinal_type xMult = operatorOrder-(zMult+yMult);
107 if (dkEnum == getDkEnumeration<3>(xMult,yMult,zMult))
118 template<ordinal_type spaceDim>
119 KOKKOS_INLINE_FUNCTION
120 void getDkEnumerationInverse(Kokkos::Array<int,spaceDim> &entries,
const ordinal_type dkEnum,
const ordinal_type operatorOrder)
127 for (
int k0=0; k0<=operatorOrder; k0++)
130 for (
int d=1; d<spaceDim-1; d++)
135 if (spaceDim > 1) entries[spaceDim-1] = operatorOrder - k0;
136 else if (k0 != operatorOrder)
continue;
137 const ordinal_type dkEnumFor_k0 = getDkEnumeration<spaceDim>(entries);
139 if (dkEnumFor_k0 > dkEnum)
continue;
140 else if (dkEnumFor_k0 == dkEnum)
return;
149 constexpr ordinal_type sizeForSubArray = (spaceDim > 2) ? spaceDim - 1 : 1;
150 Kokkos::Array<int,sizeForSubArray> subEntries;
153 getDkEnumerationInverse<spaceDim-1>(subEntries, dkEnum - dkEnumFor_k0 - 1, operatorOrder - entries[0]);
155 for (
int i=1; i<spaceDim; i++)
157 entries[i] = subEntries[i-1];
166 KOKKOS_INLINE_FUNCTION
167 ordinal_type getDkEnumeration<1>(Kokkos::Array<int,1> &entries)
169 return getDkEnumeration<1>(entries[0]);
172 template<ordinal_type spaceDim>
173 KOKKOS_INLINE_FUNCTION
176 ordinal_type k_minus_k0 = 0;
180 constexpr ordinal_type sizeForSubArray = (spaceDim > 2) ? spaceDim - 1 : 1;
181 Kokkos::Array<int,sizeForSubArray> remainingEntries;
182 for (
int i=1; i<spaceDim; i++)
184 k_minus_k0 += entries[i];
185 remainingEntries[i-1] = entries[i];
194 EOperator opFor_k_minus_k0_minus_1 = (k_minus_k0 > 1) ? EOperator(OPERATOR_D1 + k_minus_k0 - 2) : EOperator(OPERATOR_VALUE);
195 const ordinal_type dkCardinality =
getDkCardinality(opFor_k_minus_k0_minus_1, spaceDim);
196 const ordinal_type dkEnum = dkCardinality + getDkEnumeration<sizeForSubArray>(remainingEntries);
201 template<ordinal_type spaceDim1, ordinal_type spaceDim2>
202 KOKKOS_INLINE_FUNCTION
203 ordinal_type getDkTensorIndex(
const ordinal_type dkEnum1,
const ordinal_type operatorOrder1,
204 const ordinal_type dkEnum2,
const ordinal_type operatorOrder2)
206 Kokkos::Array<int,spaceDim1> entries1;
207 getDkEnumerationInverse<spaceDim1>(entries1, dkEnum1, operatorOrder1);
209 Kokkos::Array<int,spaceDim2> entries2;
210 getDkEnumerationInverse<spaceDim2>(entries2, dkEnum2, operatorOrder2);
212 const int spaceDim = spaceDim1 + spaceDim2;
213 Kokkos::Array<int,spaceDim> entries;
215 for (ordinal_type d=0; d<spaceDim1; d++)
217 entries[d] = entries1[d];
220 for (ordinal_type d=0; d<spaceDim2; d++)
222 entries[d+spaceDim1] = entries2[d];
225 return getDkEnumeration<spaceDim>(entries);
228template<
typename BasisBase>
229class Basis_TensorBasis;
237 std::vector< std::vector<EOperator> > ops;
238 std::vector<double> weights;
239 ordinal_type numBasisComponents_;
240 bool rotateXYNinetyDegrees_ =
false;
242 OperatorTensorDecomposition(
const std::vector<EOperator> &opsBasis1,
const std::vector<EOperator> &opsBasis2,
const std::vector<double> vectorComponentWeights)
244 weights(vectorComponentWeights),
245 numBasisComponents_(2)
247 const ordinal_type size = opsBasis1.size();
248 const ordinal_type opsBasis2Size = opsBasis2.size();
249 const ordinal_type weightsSize = weights.size();
250 INTREPID2_TEST_FOR_EXCEPTION(size != opsBasis2Size, std::invalid_argument,
"opsBasis1.size() != opsBasis2.size()");
251 INTREPID2_TEST_FOR_EXCEPTION(size != weightsSize, std::invalid_argument,
"opsBasis1.size() != weights.size()");
253 for (ordinal_type i=0; i<size; i++)
255 ops.push_back(std::vector<EOperator>{opsBasis1[i],opsBasis2[i]});
259 OperatorTensorDecomposition(
const std::vector< std::vector<EOperator> > &vectorEntryOps,
const std::vector<double> &vectorComponentWeights)
262 weights(vectorComponentWeights)
264 const ordinal_type numVectorComponents = ops.size();
265 const ordinal_type weightsSize = weights.size();
266 INTREPID2_TEST_FOR_EXCEPTION(numVectorComponents != weightsSize, std::invalid_argument,
"opsBasis1.size() != weights.size()");
268 INTREPID2_TEST_FOR_EXCEPTION(numVectorComponents == 0, std::invalid_argument,
"must have at least one entry!");
270 ordinal_type numBases = 0;
271 for (ordinal_type i=0; i<numVectorComponents; i++)
275 numBases = ops[i].size();
277 else if (ops[i].size() != 0)
279 const ordinal_type opsiSize = ops[i].size();
280 INTREPID2_TEST_FOR_EXCEPTION(numBases != opsiSize, std::invalid_argument,
"must have one operator for each basis in each nontrivial entry in vectorEntryOps");
283 INTREPID2_TEST_FOR_EXCEPTION(numBases == 0, std::invalid_argument,
"at least one vectorEntryOps entry must be non-trivial");
284 numBasisComponents_ = numBases;
291 numBasisComponents_(basisOps.size())
296 ops({ std::vector<EOperator>{opBasis1, opBasis2} }),
298 numBasisComponents_(2)
303 ops({ std::vector<EOperator>{opBasis1, opBasis2, opBasis3} }),
305 numBasisComponents_(3)
308 ordinal_type numVectorComponents()
const
313 ordinal_type numBasisComponents()
const
315 return numBasisComponents_;
318 double weight(
const ordinal_type &vectorComponentOrdinal)
const
320 return weights[vectorComponentOrdinal];
323 bool identicallyZeroComponent(
const ordinal_type &vectorComponentOrdinal)
const
327 return ops[vectorComponentOrdinal].size() == 0;
330 EOperator op(
const ordinal_type &vectorComponentOrdinal,
const ordinal_type &basisOrdinal)
const
334 if (identicallyZeroComponent(vectorComponentOrdinal))
342 return ops[vectorComponentOrdinal][basisOrdinal];
347 template<
typename DeviceType,
typename OutputValueType,
class Po
intValueType>
350 const ordinal_type basesSize = bases.size();
351 INTREPID2_TEST_FOR_EXCEPTION(basesSize != numBasisComponents_, std::invalid_argument,
"The number of bases provided must match the number of basis components in this decomposition");
353 ordinal_type numExpandedBasisComponents = 0;
356 std::vector<TensorBasis*> basesAsTensorBasis(numBasisComponents_);
357 for (ordinal_type basisComponentOrdinal=0; basisComponentOrdinal<numBasisComponents_; basisComponentOrdinal++)
359 TensorBasis* basisAsTensorBasis =
dynamic_cast<TensorBasis*
>(bases[basisComponentOrdinal].get());
360 basesAsTensorBasis[basisComponentOrdinal] = basisAsTensorBasis;
361 if (basisAsTensorBasis)
363 numExpandedBasisComponents += basisAsTensorBasis->getTensorBasisComponents().size();
367 numExpandedBasisComponents += 1;
371 std::vector< std::vector<EOperator> > expandedOps;
372 std::vector<double> expandedWeights;
373 const ordinal_type opsSize = ops.size();
374 for (ordinal_type simpleVectorEntryOrdinal=0; simpleVectorEntryOrdinal<opsSize; simpleVectorEntryOrdinal++)
376 if (identicallyZeroComponent(simpleVectorEntryOrdinal))
378 expandedOps.push_back(std::vector<EOperator>{});
379 expandedWeights.push_back(0.0);
383 std::vector< std::vector<EOperator> > expandedBasisOpsForSimpleVectorEntry(1);
386 auto addExpandedOp = [&expandedBasisOpsForSimpleVectorEntry](
const EOperator &op)
388 const ordinal_type size = expandedBasisOpsForSimpleVectorEntry.size();
389 for (ordinal_type i=0; i<size; i++)
391 expandedBasisOpsForSimpleVectorEntry[i].push_back(op);
397 auto vectorizeExpandedOps = [&expandedBasisOpsForSimpleVectorEntry](
const int &numSubVectors)
400 INTREPID2_TEST_FOR_EXCEPTION(expandedBasisOpsForSimpleVectorEntry.size() != 1, std::invalid_argument,
"multiple basis components may not be vector-valued!");
401 for (ordinal_type i=1; i<numSubVectors; i++)
403 expandedBasisOpsForSimpleVectorEntry.push_back(expandedBasisOpsForSimpleVectorEntry[0]);
407 std::vector<EOperator> subVectorOps;
408 std::vector<double> subVectorWeights {weights[simpleVectorEntryOrdinal]};
409 for (ordinal_type basisComponentOrdinal=0; basisComponentOrdinal<numBasisComponents_; basisComponentOrdinal++)
411 const auto &op = ops[simpleVectorEntryOrdinal][basisComponentOrdinal];
413 if (! basesAsTensorBasis[basisComponentOrdinal])
420 if (basisOpDecomposition.numVectorComponents() > 1)
423 INTREPID2_TEST_FOR_EXCEPTION(subVectorWeights.size() > 1, std::invalid_argument,
"Unhandled case: multiple component bases are vector-valued");
425 ordinal_type numSubVectors = basisOpDecomposition.numVectorComponents();
426 vectorizeExpandedOps(numSubVectors);
428 double weightSoFar = subVectorWeights[0];
429 for (ordinal_type subVectorEntryOrdinal=1; subVectorEntryOrdinal<numSubVectors; subVectorEntryOrdinal++)
431 subVectorWeights.push_back(weightSoFar * basisOpDecomposition.weight(subVectorEntryOrdinal));
433 subVectorWeights[0] *= basisOpDecomposition.weight(0);
434 for (ordinal_type subVectorEntryOrdinal=0; subVectorEntryOrdinal<numSubVectors; subVectorEntryOrdinal++)
436 for (ordinal_type subComponentBasis=0; subComponentBasis<basisOpDecomposition.numBasisComponents(); subComponentBasis++)
438 const auto &basisOp = basisOpDecomposition.op(subVectorEntryOrdinal, subComponentBasis);
439 expandedBasisOpsForSimpleVectorEntry[subVectorEntryOrdinal].push_back(basisOp);
445 double componentWeight = basisOpDecomposition.weight(0);
446 const ordinal_type size = subVectorWeights.size();
447 for (ordinal_type i=0; i<size; i++)
449 subVectorWeights[i] *= componentWeight;
451 ordinal_type subVectorEntryOrdinal = 0;
452 const ordinal_type numBasisComponents = basisOpDecomposition.numBasisComponents();
453 for (ordinal_type subComponentBasis=0; subComponentBasis<numBasisComponents; subComponentBasis++)
455 const auto &basisOp = basisOpDecomposition.op(subVectorEntryOrdinal, basisComponentOrdinal);
456 addExpandedOp( basisOp );
463 for (ordinal_type i=0; i<static_cast<ordinal_type>(expandedBasisOpsForSimpleVectorEntry.size()); i++)
465 const ordinal_type size = expandedBasisOpsForSimpleVectorEntry[i].size();
466 INTREPID2_TEST_FOR_EXCEPTION(size != numExpandedBasisComponents, std::logic_error,
"each vector in expandedBasisOpsForSimpleVectorEntry should have as many entries as there are expanded basis components");
469 expandedOps.insert(expandedOps.end(), expandedBasisOpsForSimpleVectorEntry.begin(), expandedBasisOpsForSimpleVectorEntry.end());
470 expandedWeights.insert(expandedWeights.end(), subVectorWeights.begin(), subVectorWeights.end());
473 INTREPID2_TEST_FOR_EXCEPTION(expandedOps.size() != expandedWeights.size(), std::logic_error,
"expandedWeights and expandedOps do not agree on the number of vector components");
476 result.setRotateXYNinetyDegrees(rotateXYNinetyDegrees_);
483 return rotateXYNinetyDegrees_;
486 void setRotateXYNinetyDegrees(
const bool &value)
488 rotateXYNinetyDegrees_ = value;
497 template<
class ExecutionSpace,
class OutputScalar,
class OutputFieldType>
500 using ScratchSpace =
typename ExecutionSpace::scratch_memory_space;
501 using OutputScratchView = Kokkos::View<OutputScalar*,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
503 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
504 using TeamMember =
typename TeamPolicy::member_type;
507 using RankCombinationType =
typename TensorViewIteratorType::RankCombinationType;
509 OutputFieldType output_;
510 OutputFieldType input1_;
511 OutputFieldType input2_;
513 int numFields_, numPoints_;
514 int numFields1_, numPoints1_;
515 int numFields2_, numPoints2_;
519 using RankCombinationViewType =
typename TensorViewIteratorType::RankCombinationViewType;
520 RankCombinationViewType rank_combinations_;
526 TensorViewFunctor(OutputFieldType output, OutputFieldType inputValues1, OutputFieldType inputValues2,
527 bool tensorPoints,
double weight)
528 : output_(output), input1_(inputValues1), input2_(inputValues2), tensorPoints_(tensorPoints), weight_(weight)
530 numFields_ = output.extent_int(0);
531 numPoints_ = output.extent_int(1);
533 numFields1_ = inputValues1.extent_int(0);
534 numPoints1_ = inputValues1.extent_int(1);
536 numFields2_ = inputValues2.extent_int(0);
537 numPoints2_ = inputValues2.extent_int(1);
542 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != numPoints1_, std::invalid_argument,
"incompatible point counts");
543 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != numPoints2_, std::invalid_argument,
"incompatible point counts");
547 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != numPoints1_ * numPoints2_, std::invalid_argument,
"incompatible point counts");
550 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != numFields1_ * numFields2_, std::invalid_argument,
"incompatible field sizes");
552 const ordinal_type max_rank = std::max(inputValues1.rank(),inputValues2.rank());
555 const ordinal_type outputRank = output.rank();
556 INTREPID2_TEST_FOR_EXCEPTION(outputRank > max_rank, std::invalid_argument,
"Unsupported view combination.");
557 rank_combinations_ = RankCombinationViewType(
"Rank_combinations_", max_rank);
558 auto rank_combinations_host = Kokkos::create_mirror_view(rank_combinations_);
560 rank_combinations_host[0] = TensorViewIteratorType::TENSOR_PRODUCT;
561 rank_combinations_host[1] = tensorPoints ? TensorViewIteratorType::TENSOR_PRODUCT : TensorViewIteratorType::DIMENSION_MATCH;
562 for (ordinal_type d=2; d<max_rank; d++)
566 if ((inputValues1.extent_int(d) == inputValues2.extent_int(d)) && (output.extent_int(d) == 1))
568 rank_combinations_host[d] = TensorViewIteratorType::TENSOR_CONTRACTION;
570 else if (((inputValues1.extent_int(d) == output.extent_int(d)) && (inputValues2.extent_int(d) == 1))
571 || ((inputValues2.extent_int(d) == output.extent_int(d)) && (inputValues1.extent_int(d) == 1))
576 rank_combinations_host[d] = TensorViewIteratorType::TENSOR_PRODUCT;
578 else if ((inputValues1.extent_int(d) == inputValues2.extent_int(d)) && (output.extent_int(d) == inputValues1.extent_int(d) * inputValues2.extent_int(d)))
581 rank_combinations_host[d] = TensorViewIteratorType::TENSOR_PRODUCT;
583 else if ((inputValues1.extent_int(d) == inputValues2.extent_int(d)) && (output.extent_int(d) == inputValues1.extent_int(d)))
587 rank_combinations_host[d] = TensorViewIteratorType::DIMENSION_MATCH;
591 std::cout <<
"inputValues1.extent_int(" << d <<
") = " << inputValues1.extent_int(d) << std::endl;
592 std::cout <<
"inputValues2.extent_int(" << d <<
") = " << inputValues2.extent_int(d) << std::endl;
593 std::cout <<
"output.extent_int(" << d <<
") = " << output.extent_int(d) << std::endl;
594 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"unable to find an interpretation for this combination of views");
597 Kokkos::deep_copy(rank_combinations_,rank_combinations_host);
600 KOKKOS_INLINE_FUNCTION
601 void operator()(
const TeamMember & teamMember )
const
603 auto fieldOrdinal1 = teamMember.league_rank();
605 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (
const int& fieldOrdinal2) {
607 const int FIELD_ORDINAL_DIMENSION = 0;
608 it.
setLocation({fieldOrdinal1,0,0,0,0,0,0},{fieldOrdinal2,0,0,0,0,0,0});
609 int next_increment_rank = FIELD_ORDINAL_DIMENSION;
610 OutputScalar accumulator = 0;
617 if ((next_increment_rank < 0) || (rank_combinations_[next_increment_rank] != TensorViewIteratorType::TENSOR_CONTRACTION))
624 }
while (it.
increment() > FIELD_ORDINAL_DIMENSION);
642 template<
typename BasisBaseClass =
void>
645 public BasisBaseClass
648 using BasisBase = BasisBaseClass;
649 using BasisPtr = Teuchos::RCP<BasisBase>;
655 std::vector<BasisPtr> tensorComponents_;
659 int numTensorialExtrusions_;
661 using DeviceType =
typename BasisBase::DeviceType;
662 using ExecutionSpace =
typename BasisBase::ExecutionSpace;
663 using OutputValueType =
typename BasisBase::OutputValueType;
664 using PointValueType =
typename BasisBase::PointValueType;
666 using OrdinalTypeArray1DHost =
typename BasisBase::OrdinalTypeArray1DHost;
667 using OrdinalTypeArray2DHost =
typename BasisBase::OrdinalTypeArray2DHost;
668 using OutputViewType =
typename BasisBase::OutputViewType;
669 using PointViewType =
typename BasisBase::PointViewType;
678 Basis_TensorBasis(BasisPtr basis1, BasisPtr basis2, EFunctionSpace functionSpace = FUNCTION_SPACE_MAX,
679 const bool useShardsCellTopologyAndTags =
false)
681 basis1_(basis1),basis2_(basis2)
683 this->functionSpace_ = functionSpace;
688 auto basis1Components = basis1AsTensor->getTensorBasisComponents();
689 tensorComponents_.insert(tensorComponents_.end(), basis1Components.begin(), basis1Components.end());
693 tensorComponents_.push_back(basis1_);
699 auto basis2Components = basis2AsTensor->getTensorBasisComponents();
700 tensorComponents_.insert(tensorComponents_.end(), basis2Components.begin(), basis2Components.end());
704 tensorComponents_.push_back(basis2_);
707 this->basisCardinality_ = basis1->getCardinality() * basis2->getCardinality();
708 this->basisDegree_ = std::max(basis1->getDegree(), basis2->getDegree());
711 std::ostringstream basisName;
712 basisName << basis1->getName() <<
" x " << basis2->getName();
713 name_ = basisName.str();
717 this->basisCellTopology_ = tensorComponents_[0]->getBaseCellTopology();
718 this->numTensorialExtrusions_ = tensorComponents_.size() - 1;
720 this->basisType_ = basis1_->getBasisType();
721 this->basisCoordinates_ = COORDINATES_CARTESIAN;
723 ordinal_type spaceDim1 = basis1_->getDomainDimension();
724 ordinal_type spaceDim2 = basis2_->getDomainDimension();
726 INTREPID2_TEST_FOR_EXCEPTION(spaceDim2 != 1, std::invalid_argument,
"TensorBasis only supports 1D bases in basis2_ position");
728 if (this->getBasisType() == BASIS_FEM_HIERARCHICAL)
731 int degreeSize = basis1_->getPolynomialDegreeLength() + basis2_->getPolynomialDegreeLength();
732 this->fieldOrdinalPolynomialDegree_ = OrdinalTypeArray2DHost(
"TensorBasis - field ordinal polynomial degree", this->basisCardinality_, degreeSize);
733 this->fieldOrdinalH1PolynomialDegree_ = OrdinalTypeArray2DHost(
"TensorBasis - field ordinal polynomial H^1 degree", this->basisCardinality_, degreeSize);
735 const ordinal_type basis1Cardinality = basis1_->getCardinality();
736 const ordinal_type basis2Cardinality = basis2_->getCardinality();
738 int degreeLengthField1 = basis1_->getPolynomialDegreeLength();
739 int degreeLengthField2 = basis2_->getPolynomialDegreeLength();
741 for (ordinal_type fieldOrdinal1 = 0; fieldOrdinal1 < basis1Cardinality; fieldOrdinal1++)
743 OrdinalTypeArray1DHost degreesField1 = basis1_->getPolynomialDegreeOfField(fieldOrdinal1);
744 OrdinalTypeArray1DHost h1DegreesField1 = basis1_->getH1PolynomialDegreeOfField(fieldOrdinal1);
745 for (ordinal_type fieldOrdinal2 = 0; fieldOrdinal2 < basis2Cardinality; fieldOrdinal2++)
747 OrdinalTypeArray1DHost degreesField2 = basis2_->getPolynomialDegreeOfField(fieldOrdinal2);
748 OrdinalTypeArray1DHost h1DegreesField2 = basis2_->getH1PolynomialDegreeOfField(fieldOrdinal2);
749 const ordinal_type tensorFieldOrdinal = fieldOrdinal2 * basis1Cardinality + fieldOrdinal1;
751 for (
int d3=0; d3<degreeLengthField1; d3++)
753 this->fieldOrdinalPolynomialDegree_ (tensorFieldOrdinal,d3) = degreesField1(d3);
754 this->fieldOrdinalH1PolynomialDegree_(tensorFieldOrdinal,d3) = h1DegreesField1(d3);
756 for (
int d3=0; d3<degreeLengthField2; d3++)
758 this->fieldOrdinalPolynomialDegree_ (tensorFieldOrdinal,d3+degreeLengthField1) = degreesField2(d3);
759 this->fieldOrdinalH1PolynomialDegree_(tensorFieldOrdinal,d3+degreeLengthField1) = h1DegreesField2(d3);
765 if (useShardsCellTopologyAndTags)
767 setShardsTopologyAndTags();
773 const auto & cardinality = this->basisCardinality_;
776 const ordinal_type tagSize = 4;
777 const ordinal_type posScDim = 0;
778 const ordinal_type posScOrd = 1;
779 const ordinal_type posDfOrd = 2;
780 const ordinal_type posDfCnt = 3;
782 OrdinalTypeArray1DHost tagView(
"tag view", cardinality*tagSize);
786 auto basis1Topo = cellTopo->getTensorialComponent();
788 const ordinal_type spaceDim = spaceDim1 + spaceDim2;
789 const ordinal_type sideDim = spaceDim - 1;
791 const OrdinalTypeArray2DHost ordinalToTag1 = basis1_->getAllDofTags();
792 const OrdinalTypeArray2DHost ordinalToTag2 = basis2_->getAllDofTags();
794 for (
int fieldOrdinal1=0; fieldOrdinal1<basis1_->getCardinality(); fieldOrdinal1++)
796 ordinal_type subcellDim1 = ordinalToTag1(fieldOrdinal1,posScDim);
797 ordinal_type subcellOrd1 = ordinalToTag1(fieldOrdinal1,posScOrd);
798 ordinal_type subcellDfCnt1 = ordinalToTag1(fieldOrdinal1,posDfCnt);
799 for (
int fieldOrdinal2=0; fieldOrdinal2<basis2_->getCardinality(); fieldOrdinal2++)
801 ordinal_type subcellDim2 = ordinalToTag2(fieldOrdinal2,posScDim);
802 ordinal_type subcellOrd2 = ordinalToTag2(fieldOrdinal2,posScOrd);
803 ordinal_type subcellDfCnt2 = ordinalToTag2(fieldOrdinal2,posDfCnt);
805 ordinal_type subcellDim = subcellDim1 + subcellDim2;
806 ordinal_type subcellOrd;
807 if (subcellDim2 == 0)
811 ordinal_type sideOrdinal = cellTopo->getTensorialComponentSideOrdinal(subcellOrd2);
813 subcellDim1, subcellOrd1);
818 subcellOrd = cellTopo->getExtrudedSubcellOrdinal(subcellDim1, subcellOrd1);
819 if (subcellOrd == -1)
821 std::cout <<
"ERROR: -1 subcell ordinal.\n";
822 subcellOrd = cellTopo->getExtrudedSubcellOrdinal(subcellDim1, subcellOrd1);
825 ordinal_type tensorFieldOrdinal = fieldOrdinal2 * basis1_->getCardinality() + fieldOrdinal1;
827 ordinal_type dofOffsetOrdinal1 = ordinalToTag1(fieldOrdinal1,posDfOrd);
828 ordinal_type dofOffsetOrdinal2 = ordinalToTag2(fieldOrdinal2,posDfOrd);
829 ordinal_type dofsForSubcell1 = ordinalToTag1(fieldOrdinal1,posDfCnt);
830 ordinal_type dofOffsetOrdinal = dofOffsetOrdinal2 * dofsForSubcell1 + dofOffsetOrdinal1;
831 tagView(tagSize*tensorFieldOrdinal + posScDim) = subcellDim;
832 tagView(tagSize*tensorFieldOrdinal + posScOrd) = subcellOrd;
833 tagView(tagSize*tensorFieldOrdinal + posDfOrd) = dofOffsetOrdinal;
834 tagView(tagSize*tensorFieldOrdinal + posDfCnt) = subcellDfCnt1 * subcellDfCnt2;
840 this->setOrdinalTagData(this->tagToOrdinal_,
843 this->basisCardinality_,
851 void setShardsTopologyAndTags()
853 shards::CellTopology cellTopo1 = basis1_->getBaseCellTopology();
854 shards::CellTopology cellTopo2 = basis2_->getBaseCellTopology();
856 auto cellKey1 = basis1_->getBaseCellTopology().getKey();
857 auto cellKey2 = basis2_->getBaseCellTopology().getKey();
859 const int numTensorialExtrusions = basis1_->getNumTensorialExtrusions() + basis2_->getNumTensorialExtrusions();
860 if ((cellKey1 == shards::Line<2>::key) && (cellKey2 == shards::Line<2>::key) && (numTensorialExtrusions == 0))
862 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
864 else if ( ((cellKey1 == shards::Quadrilateral<4>::key) && (cellKey2 == shards::Line<2>::key))
865 || ((cellKey2 == shards::Quadrilateral<4>::key) && (cellKey1 == shards::Line<2>::key))
866 || ((cellKey1 == shards::Line<2>::key) && (cellKey2 == shards::Line<2>::key) && (numTensorialExtrusions == 1))
869 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
871 else if ((cellKey1 == shards::Triangle<3>::key) && (cellKey2 == shards::Line<2>::key))
873 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Wedge<6> >() );
877 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Cell topology combination not yet supported");
881 numTensorialExtrusions_ = 0;
885 const auto & cardinality = this->basisCardinality_;
888 const ordinal_type tagSize = 4;
889 const ordinal_type posScDim = 0;
890 const ordinal_type posScOrd = 1;
891 const ordinal_type posDfOrd = 2;
893 OrdinalTypeArray1DHost tagView(
"tag view", cardinality*tagSize);
895 shards::CellTopology cellTopo = this->basisCellTopology_;
897 ordinal_type tensorSpaceDim = cellTopo.getDimension();
898 ordinal_type spaceDim1 = cellTopo1.getDimension();
899 ordinal_type spaceDim2 = cellTopo2.getDimension();
901 TensorTopologyMap topoMap(cellTopo1, cellTopo2);
903 for (ordinal_type d=0; d<=tensorSpaceDim; d++)
905 ordinal_type d2_max = std::min(spaceDim2,d);
906 int subcellOffset = 0;
907 for (ordinal_type d2=0; d2<=d2_max; d2++)
909 ordinal_type d1 = d-d2;
910 if (d1 > spaceDim1)
continue;
912 ordinal_type subcellCount2 = cellTopo2.getSubcellCount(d2);
913 ordinal_type subcellCount1 = cellTopo1.getSubcellCount(d1);
914 for (ordinal_type subcellOrdinal2=0; subcellOrdinal2<subcellCount2; subcellOrdinal2++)
916 ordinal_type subcellDofCount2 = basis2_->getDofCount(d2, subcellOrdinal2);
917 for (ordinal_type subcellOrdinal1=0; subcellOrdinal1<subcellCount1; subcellOrdinal1++)
919 ordinal_type subcellDofCount1 = basis1_->getDofCount(d1, subcellOrdinal1);
920 ordinal_type tensorLocalDofCount = subcellDofCount1 * subcellDofCount2;
921 for (ordinal_type localDofID2 = 0; localDofID2<subcellDofCount2; localDofID2++)
923 ordinal_type fieldOrdinal2 = basis2_->getDofOrdinal(d2, subcellOrdinal2, localDofID2);
924 OrdinalTypeArray1DHost degreesField2;
925 if (this->basisType_ == BASIS_FEM_HIERARCHICAL) degreesField2 = basis2_->getPolynomialDegreeOfField(fieldOrdinal2);
926 for (ordinal_type localDofID1 = 0; localDofID1<subcellDofCount1; localDofID1++)
928 ordinal_type fieldOrdinal1 = basis1_->getDofOrdinal(d1, subcellOrdinal1, localDofID1);
929 ordinal_type tensorLocalDofID = localDofID2 * subcellDofCount1 + localDofID1;
930 ordinal_type tensorFieldOrdinal = fieldOrdinal2 * basis1_->getCardinality() + fieldOrdinal1;
931 tagView(tensorFieldOrdinal*tagSize+0) = d;
932 tagView(tensorFieldOrdinal*tagSize+1) = topoMap.getCompositeSubcellOrdinal(d1, subcellOrdinal1, d2, subcellOrdinal2);
933 tagView(tensorFieldOrdinal*tagSize+2) = tensorLocalDofID;
934 tagView(tensorFieldOrdinal*tagSize+3) = tensorLocalDofCount;
939 subcellOffset += subcellCount1 * subcellCount2;
945 this->setOrdinalTagData(this->tagToOrdinal_,
948 this->basisCardinality_,
956 virtual int getNumTensorialExtrusions()
const override
958 return numTensorialExtrusions_;
970 ordinal_type dkEnum2, ordinal_type operatorOrder2)
const
972 ordinal_type spaceDim1 = basis1_->getDomainDimension();
973 ordinal_type spaceDim2 = basis2_->getDomainDimension();
980 INTREPID2_TEST_FOR_EXCEPTION(operatorOrder1 > 0, std::invalid_argument,
"For spaceDim1 = 0, operatorOrder1 must be 0.");
986 case 1:
return getDkTensorIndex<1, 1>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
987 case 2:
return getDkTensorIndex<1, 2>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
988 case 3:
return getDkTensorIndex<1, 3>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
989 case 4:
return getDkTensorIndex<1, 4>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
990 case 5:
return getDkTensorIndex<1, 5>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
991 case 6:
return getDkTensorIndex<1, 6>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
993 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported dimension combination");
998 case 1:
return getDkTensorIndex<2, 1>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
999 case 2:
return getDkTensorIndex<2, 2>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1000 case 3:
return getDkTensorIndex<2, 3>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1001 case 4:
return getDkTensorIndex<2, 4>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1002 case 5:
return getDkTensorIndex<2, 5>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1004 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported dimension combination");
1009 case 1:
return getDkTensorIndex<3, 1>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1010 case 2:
return getDkTensorIndex<3, 2>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1011 case 3:
return getDkTensorIndex<3, 3>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1012 case 4:
return getDkTensorIndex<3, 4>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1014 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported dimension combination");
1019 case 1:
return getDkTensorIndex<4, 1>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1020 case 2:
return getDkTensorIndex<4, 2>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1021 case 3:
return getDkTensorIndex<4, 3>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1023 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported dimension combination");
1028 case 1:
return getDkTensorIndex<5, 1>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1029 case 2:
return getDkTensorIndex<5, 2>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1031 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported dimension combination");
1036 case 1:
return getDkTensorIndex<6, 1>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1038 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported dimension combination");
1041 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported dimension combination");
1051 const int spaceDim = this->getDomainDimension();
1053 const EOperator VALUE = Intrepid2::OPERATOR_VALUE;
1055 std::vector< std::vector<EOperator> > opsVALUE{{VALUE, VALUE}};
1057 std::vector< std::vector<EOperator> > ops(spaceDim);
1059 switch (operatorType)
1067 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported operator type - TensorBasis subclass should override");
1085 ops = std::vector< std::vector<EOperator> >(dkCardinality);
1089 for (
int derivativeCountComp2=0; derivativeCountComp2<=opOrder; derivativeCountComp2++)
1091 int derivativeCountComp1=opOrder-derivativeCountComp2;
1092 EOperator op1 = (derivativeCountComp1 == 0) ? OPERATOR_VALUE : EOperator(OPERATOR_D1 + (derivativeCountComp1 - 1));
1093 EOperator op2 = (derivativeCountComp2 == 0) ? OPERATOR_VALUE : EOperator(OPERATOR_D1 + (derivativeCountComp2 - 1));
1098 for (
int dkEnum1=0; dkEnum1<dkCardinality1; dkEnum1++)
1100 for (
int dkEnum2=0; dkEnum2<dkCardinality2; dkEnum2++)
1102 ordinal_type dkTensorIndex = getDkTensorIndex<1, 1>(dkEnum1, derivativeCountComp1, dkEnum2, derivativeCountComp2);
1103 ops[dkTensorIndex] = std::vector<EOperator>{op1, op2};
1111 std::vector<double> weights(ops.size(), 1.0);
1119 if (((operatorType >= OPERATOR_D1) && (operatorType <= OPERATOR_D10)) || (operatorType == OPERATOR_GRAD))
1126 ordinal_type numBasisComponents = tensorComponents_.size();
1129 const int dkCardinality =
getDkCardinality(operatorType, numBasisComponents);
1131 std::vector< std::vector<EOperator> > ops(dkCardinality);
1133 std::vector<EOperator> prevEntry(numBasisComponents, OPERATOR_VALUE);
1134 prevEntry[0] = operatorType;
1138 for (ordinal_type dkOrdinal=1; dkOrdinal<dkCardinality; dkOrdinal++)
1140 std::vector<EOperator> entry = prevEntry;
1149 ordinal_type cumulativeOpOrder = 0;
1151 for (ordinal_type compOrdinal=0; compOrdinal<numBasisComponents; compOrdinal++)
1154 cumulativeOpOrder += thisOpOrder;
1155 if (cumulativeOpOrder + finalOpOrder == opOrder)
1158 EOperator decrementedOp;
1159 if (thisOpOrder == 1)
1161 decrementedOp = OPERATOR_VALUE;
1165 decrementedOp =
static_cast<EOperator
>(OPERATOR_D1 + ((thisOpOrder - 1) - 1));
1167 entry[compOrdinal] = decrementedOp;
1168 const ordinal_type remainingOpOrder = opOrder - cumulativeOpOrder + 1;
1169 entry[compOrdinal+1] =
static_cast<EOperator
>(OPERATOR_D1 + (remainingOpOrder - 1));
1170 for (ordinal_type i=compOrdinal+2; i<numBasisComponents; i++)
1172 entry[i] = OPERATOR_VALUE;
1177 ops[dkOrdinal] = entry;
1180 std::vector<double> weights(dkCardinality, 1.0);
1187 std::vector<BasisPtr> componentBases {basis1_, basis2_};
1198 const bool operatorIsDk = (operatorType >= OPERATOR_D1) && (operatorType <= OPERATOR_D10);
1199 const bool operatorSupported = (operatorType == OPERATOR_VALUE) || (operatorType == OPERATOR_GRAD) || (operatorType == OPERATOR_CURL) || (operatorType == OPERATOR_DIV) || operatorIsDk;
1200 INTREPID2_TEST_FOR_EXCEPTION(!operatorSupported, std::invalid_argument,
"operator is not supported by allocateBasisValues");
1203 const int spaceDim = this->getDomainDimension();
1204 INTREPID2_TEST_FOR_EXCEPTION(spaceDim != points.
extent_int(1), std::invalid_argument,
"points must be shape (P,D), with D equal to the dimension of the basis domain");
1207 ordinal_type numBasisComponents = tensorComponents_.size();
1213 INTREPID2_TEST_FOR_EXCEPTION(points.
numTensorComponents() != 1, std::invalid_argument,
"If points does not have the same number of tensor components as the basis, then it should have trivial tensor structure.");
1214 const ordinal_type numPoints = points.
extent_int(0);
1215 auto outputView = this->allocateOutputView(numPoints, operatorType);
1222 INTREPID2_TEST_FOR_EXCEPTION(numBasisComponents > points.
numTensorComponents(), std::invalid_argument,
"points must have at least as many tensorial components as basis.");
1226 ordinal_type numVectorComponents = opDecomposition.numVectorComponents();
1227 const bool useVectorData = numVectorComponents > 1;
1229 std::vector<ordinal_type> componentPointCounts(numBasisComponents);
1230 ordinal_type pointComponentNumber = 0;
1231 for (ordinal_type r=0; r<numBasisComponents; r++)
1233 const ordinal_type compSpaceDim = tensorComponents_[r]->getDomainDimension();
1234 ordinal_type dimsSoFar = 0;
1235 ordinal_type numPointsForBasisComponent = 1;
1236 while (dimsSoFar < compSpaceDim)
1238 INTREPID2_TEST_FOR_EXCEPTION(pointComponentNumber >= points.
numTensorComponents(), std::invalid_argument,
"Error in processing points container; perhaps it is mis-sized?");
1240 const int numComponentDims = points.
getTensorComponent(pointComponentNumber).extent_int(1);
1241 numPointsForBasisComponent *= numComponentPoints;
1242 dimsSoFar += numComponentDims;
1243 INTREPID2_TEST_FOR_EXCEPTION(dimsSoFar > points.
numTensorComponents(), std::invalid_argument,
"Error in processing points container; perhaps it is mis-sized?");
1244 pointComponentNumber++;
1246 componentPointCounts[r] = numPointsForBasisComponent;
1251 const int numFamilies = 1;
1254 const int familyOrdinal = 0;
1255 for (ordinal_type vectorComponentOrdinal=0; vectorComponentOrdinal<numVectorComponents; vectorComponentOrdinal++)
1257 if (!opDecomposition.identicallyZeroComponent(vectorComponentOrdinal))
1259 std::vector< Data<OutputValueType,DeviceType> > componentData;
1260 for (ordinal_type r=0; r<numBasisComponents; r++)
1262 const int numComponentPoints = componentPointCounts[r];
1263 const EOperator op = opDecomposition.op(vectorComponentOrdinal, r);
1264 auto componentView = tensorComponents_[r]->allocateOutputView(numComponentPoints, op);
1276 std::vector< Data<OutputValueType,DeviceType> > componentData;
1278 const ordinal_type vectorComponentOrdinal = 0;
1279 for (ordinal_type r=0; r<numBasisComponents; r++)
1281 const int numComponentPoints = componentPointCounts[r];
1282 const EOperator op = opDecomposition.op(vectorComponentOrdinal, r);
1283 auto componentView = tensorComponents_[r]->allocateOutputView(numComponentPoints, op);
1288 const Kokkos::Array<int,7> extents {componentView.extent_int(0), componentView.extent_int(1), 1,1,1,1,1};
1295 std::vector< TensorData<OutputValueType,DeviceType> > tensorDataEntries {tensorData};
1303 using BasisBase::getValues;
1317 PointViewType & inputPoints1, PointViewType & inputPoints2,
bool &tensorDecompositionSucceeded)
const
1319 INTREPID2_TEST_FOR_EXCEPTION(attemptTensorDecomposition, std::invalid_argument,
"tensor decomposition not yet supported");
1333 int spaceDim1 = basis1_->getDomainDimension();
1334 int spaceDim2 = basis2_->getDomainDimension();
1336 int totalSpaceDim = inputPoints.extent_int(1);
1338 TEUCHOS_ASSERT(spaceDim1 + spaceDim2 == totalSpaceDim);
1342 inputPoints1 = Kokkos::subview(inputPoints,Kokkos::ALL(),std::make_pair(0,spaceDim1));
1343 inputPoints2 = Kokkos::subview(inputPoints,Kokkos::ALL(),std::make_pair(spaceDim1,totalSpaceDim));
1349 tensorDecompositionSucceeded =
false;
1360 virtual void getDofCoords(
typename BasisBase::ScalarViewType dofCoords )
const override
1362 int spaceDim1 = basis1_->getBaseCellTopology().getDimension();
1363 int spaceDim2 = basis2_->getBaseCellTopology().getDimension();
1365 using ValueType =
typename BasisBase::ScalarViewType::value_type;
1366 using ResultLayout =
typename DeduceLayout< typename BasisBase::ScalarViewType >::result_layout;
1367 using ViewType = Kokkos::DynRankView<ValueType, ResultLayout, DeviceType >;
1369 const ordinal_type basisCardinality1 = basis1_->getCardinality();
1370 const ordinal_type basisCardinality2 = basis2_->getCardinality();
1372 ViewType dofCoords1(
"dofCoords1",basisCardinality1,spaceDim1);
1373 ViewType dofCoords2(
"dofCoords2",basisCardinality2,spaceDim2);
1375 basis1_->getDofCoords(dofCoords1);
1376 basis2_->getDofCoords(dofCoords2);
1378 Kokkos::RangePolicy<ExecutionSpace> policy(0, basisCardinality2);
1379 Kokkos::parallel_for(policy, KOKKOS_LAMBDA (
const int fieldOrdinal2)
1381 for (
int fieldOrdinal1=0; fieldOrdinal1<basisCardinality1; fieldOrdinal1++)
1383 const ordinal_type fieldOrdinal = fieldOrdinal1 + fieldOrdinal2 * basisCardinality1;
1384 for (
int d1=0; d1<spaceDim1; d1++)
1386 dofCoords(fieldOrdinal,d1) = dofCoords1(fieldOrdinal1,d1);
1388 for (
int d2=0; d2<spaceDim2; d2++)
1390 dofCoords(fieldOrdinal,spaceDim1+d2) = dofCoords2(fieldOrdinal2,d2);
1404 virtual void getDofCoeffs(
typename BasisBase::ScalarViewType dofCoeffs )
const override
1406 using ValueType =
typename BasisBase::ScalarViewType::value_type;
1407 using ResultLayout =
typename DeduceLayout< typename BasisBase::ScalarViewType >::result_layout;
1408 using ViewType = Kokkos::DynRankView<ValueType, ResultLayout, DeviceType >;
1410 ViewType dofCoeffs1(
"dofCoeffs1",basis1_->getCardinality());
1411 ViewType dofCoeffs2(
"dofCoeffs2",basis2_->getCardinality());
1413 basis1_->getDofCoeffs(dofCoeffs1);
1414 basis2_->getDofCoeffs(dofCoeffs2);
1416 const ordinal_type basisCardinality1 = basis1_->getCardinality();
1417 const ordinal_type basisCardinality2 = basis2_->getCardinality();
1419 Kokkos::RangePolicy<ExecutionSpace> policy(0, basisCardinality2);
1420 Kokkos::parallel_for(policy, KOKKOS_LAMBDA (
const int fieldOrdinal2)
1422 for (
int fieldOrdinal1=0; fieldOrdinal1<basisCardinality1; fieldOrdinal1++)
1424 const ordinal_type fieldOrdinal = fieldOrdinal1 + fieldOrdinal2 * basisCardinality1;
1425 dofCoeffs(fieldOrdinal) = dofCoeffs1(fieldOrdinal1);
1426 dofCoeffs(fieldOrdinal) *= dofCoeffs2(fieldOrdinal2);
1438 return name_.c_str();
1441 std::vector<BasisPtr> getTensorBasisComponents()
const
1443 return tensorComponents_;
1461 const EOperator operatorType = OPERATOR_VALUE )
const override
1463 const ordinal_type numTensorComponents = tensorComponents_.size();
1467 INTREPID2_TEST_FOR_EXCEPTION( inputPoints.
numTensorComponents() != 1, std::invalid_argument,
"If inputPoints differs from the tensor basis in component count, then inputPoints must have trivial tensor product structure" );
1468 INTREPID2_TEST_FOR_EXCEPTION( outputValues.
numFamilies() != 1, std::invalid_argument,
"If inputPoints differs from the tensor basis in component count, outputValues must have a single family with trivial tensor product structure" );
1469 INTREPID2_TEST_FOR_EXCEPTION( outputValues.
tensorData().
numTensorComponents() != 1, std::invalid_argument,
"If inputPoints differs from the tensor basis in component count, outputValues must have a single family with trivial tensor product structure" );
1473 this->
getValues(outputView, pointView, operatorType);
1479 const ordinal_type numVectorComponents = operatorDecomposition.numVectorComponents();
1480 const bool useVectorData = numVectorComponents > 1;
1481 const ordinal_type numBasisComponents = operatorDecomposition.numBasisComponents();
1483 for (ordinal_type vectorComponentOrdinal=0; vectorComponentOrdinal<numVectorComponents; vectorComponentOrdinal++)
1485 const double weight = operatorDecomposition.weight(vectorComponentOrdinal);
1486 ordinal_type pointComponentOrdinal = 0;
1487 for (ordinal_type basisOrdinal=0; basisOrdinal<numBasisComponents; basisOrdinal++, pointComponentOrdinal++)
1489 const EOperator op = operatorDecomposition.op(vectorComponentOrdinal, basisOrdinal);
1491 if (op != OPERATOR_MAX)
1493 const int vectorFamily = 0;
1495 INTREPID2_TEST_FOR_EXCEPTION( ! tensorData.getTensorComponent(basisOrdinal).isValid(), std::invalid_argument,
"Invalid output component encountered");
1501 const ordinal_type basisDomainDimension = tensorComponents_[basisOrdinal]->getDomainDimension();
1502 if (pointView.extent_int(1) == basisDomainDimension)
1504 tensorComponents_[basisOrdinal]->getValues(basisValueView, pointView, op);
1511 ordinal_type dimsSoFar = 0;
1512 std::vector< ScalarView<PointValueType,DeviceType> > basisPointComponents;
1513 while (dimsSoFar < basisDomainDimension)
1515 INTREPID2_TEST_FOR_EXCEPTION(pointComponentOrdinal >= inputPoints.
numTensorComponents(), std::invalid_argument,
"Error in processing points container; perhaps it is mis-sized?");
1517 const ordinal_type numComponentDims = pointComponent.extent_int(1);
1518 dimsSoFar += numComponentDims;
1519 INTREPID2_TEST_FOR_EXCEPTION(dimsSoFar > inputPoints.
numTensorComponents(), std::invalid_argument,
"Error in processing points container; perhaps it is mis-sized?");
1520 basisPointComponents.push_back(pointComponent);
1521 if (dimsSoFar < basisDomainDimension)
1524 pointComponentOrdinal++;
1530 bool useVectorData2 = (basisValueView.rank() == 3);
1544 tensorComponents_[basisOrdinal]->getValues(basisValues, basisPoints, op);
1552 const bool spansXY = (vectorComponentOrdinal == 0) && (basisValueView.extent_int(2) == 2);
1556 auto policy = Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<2>>({0,0},{basisValueView.extent_int(0),basisValueView.extent_int(1)});
1557 Kokkos::parallel_for(
"rotateXYNinetyDegrees", policy,
1558 KOKKOS_LAMBDA (
const int &fieldOrdinal,
const int &pointOrdinal) {
1559 const auto f_x = basisValueView(fieldOrdinal,pointOrdinal,0);
1560 const auto &f_y = basisValueView(fieldOrdinal,pointOrdinal,1);
1561 basisValueView(fieldOrdinal,pointOrdinal,0) = -f_y;
1562 basisValueView(fieldOrdinal,pointOrdinal,1) = f_x;
1568 if ((weight != 1.0) && (basisOrdinal == 0))
1570 if (basisValueView.rank() == 2)
1572 auto policy = Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<2>>({0,0},{basisValueView.extent_int(0),basisValueView.extent_int(1)});
1573 Kokkos::parallel_for(
"multiply basisValueView by weight", policy,
1574 KOKKOS_LAMBDA (
const int &fieldOrdinal,
const int &pointOrdinal) {
1575 basisValueView(fieldOrdinal,pointOrdinal) *= weight;
1578 else if (basisValueView.rank() == 3)
1580 auto policy = Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<3>>({0,0,0},{basisValueView.extent_int(0),basisValueView.extent_int(1),basisValueView.extent_int(2)});
1581 Kokkos::parallel_for(
"multiply basisValueView by weight", policy,
1582 KOKKOS_LAMBDA (
const int &fieldOrdinal,
const int &pointOrdinal,
const int &d) {
1583 basisValueView(fieldOrdinal,pointOrdinal,d) *= weight;
1588 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported rank for basisValueView");
1614 void getValues( OutputViewType outputValues,
const PointViewType inputPoints,
1615 const EOperator operatorType = OPERATOR_VALUE )
const override
1618 bool attemptTensorDecomposition =
false;
1619 PointViewType inputPoints1, inputPoints2;
1620 getComponentPoints(inputPoints, attemptTensorDecomposition, inputPoints1, inputPoints2, tensorPoints);
1622 const auto functionSpace = this->getFunctionSpace();
1624 if ((functionSpace == FUNCTION_SPACE_HVOL) || (functionSpace == FUNCTION_SPACE_HGRAD))
1627 switch (operatorType)
1629 case OPERATOR_VALUE:
1645 for (
int derivativeCountComp2=0; derivativeCountComp2<=opOrder; derivativeCountComp2++)
1647 int derivativeCountComp1=opOrder-derivativeCountComp2;
1648 EOperator op1 = (derivativeCountComp1 == 0) ? OPERATOR_VALUE : EOperator(OPERATOR_D1 + (derivativeCountComp1 - 1));
1649 EOperator op2 = (derivativeCountComp2 == 0) ? OPERATOR_VALUE : EOperator(OPERATOR_D1 + (derivativeCountComp2 - 1));
1651 int spaceDim1 = inputPoints1.extent_int(1);
1652 int spaceDim2 = inputPoints2.extent_int(1);
1654 int dkCardinality1 = (op1 != OPERATOR_VALUE) ?
getDkCardinality(op1, spaceDim1) : 1;
1655 int dkCardinality2 = (op2 != OPERATOR_VALUE) ?
getDkCardinality(op2, spaceDim2) : 1;
1657 int basisCardinality1 = basis1_->getCardinality();
1658 int basisCardinality2 = basis2_->getCardinality();
1660 int totalPointCount = tensorPoints ? inputPoints1.extent_int(0) * inputPoints2.extent_int(0) : inputPoints1.extent_int(0);
1662 int pointCount1, pointCount2;
1665 pointCount1 = inputPoints1.extent_int(0);
1666 pointCount2 = inputPoints2.extent_int(0);
1670 pointCount1 = totalPointCount;
1671 pointCount2 = totalPointCount;
1674 OutputViewType outputValues1, outputValues2;
1675 if (op1 == OPERATOR_VALUE)
1678 outputValues1 =
getMatchingViewWithLabel(outputValues,
"output values - basis 1",basisCardinality1,pointCount1,dkCardinality1);
1680 if (op2 == OPERATOR_VALUE)
1683 outputValues2 =
getMatchingViewWithLabel(outputValues,
"output values - basis 2",basisCardinality2,pointCount2,dkCardinality2);
1685 basis1_->getValues(outputValues1,inputPoints1,op1);
1686 basis2_->getValues(outputValues2,inputPoints2,op2);
1688 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputValueType>();
1689 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointValueType>();
1690 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
1692 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(basisCardinality1,Kokkos::AUTO(),vectorSize);
1694 double weight = 1.0;
1697 for (
int dkEnum1=0; dkEnum1<dkCardinality1; dkEnum1++)
1699 auto outputValues1_dkEnum1 = (op1 != OPERATOR_VALUE) ? Kokkos::subview(outputValues1,Kokkos::ALL(),Kokkos::ALL(),dkEnum1)
1700 : Kokkos::subview(outputValues1,Kokkos::ALL(),Kokkos::ALL());
1701 for (
int dkEnum2=0; dkEnum2<dkCardinality2; dkEnum2++)
1703 auto outputValues2_dkEnum2 = (op2 != OPERATOR_VALUE) ? Kokkos::subview(outputValues2,Kokkos::ALL(),Kokkos::ALL(),dkEnum2)
1704 : Kokkos::subview(outputValues2,Kokkos::ALL(),Kokkos::ALL());
1706 ordinal_type dkTensorIndex =
getTensorDkEnumeration(dkEnum1, derivativeCountComp1, dkEnum2, derivativeCountComp2);
1707 auto outputValues_dkTensor = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),dkTensorIndex);
1713 FunctorType functor(outputValues_dkTensor, outputValues1_dkEnum1, outputValues2_dkEnum2, tensorPoints, weight);
1714 Kokkos::parallel_for(
"TensorViewFunctor", policy , functor);
1721 this->
getValues(outputValues, operatorType, inputPoints1, inputPoints2, tensorPoints);
1727 this->
getValues(outputValues, operatorType, inputPoints1, inputPoints2, tensorPoints);
1756 virtual void getValues(OutputViewType outputValues,
const EOperator operatorType,
1757 const PointViewType inputPoints1,
const PointViewType inputPoints2,
1758 bool tensorPoints)
const
1760 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"one-operator, two-inputPoints getValues should be overridden by TensorBasis subclasses");
1787 const PointViewType inputPoints1,
const EOperator operatorType1,
1788 const PointViewType inputPoints2,
const EOperator operatorType2,
1789 bool tensorPoints,
double weight=1.0)
const
1791 int basisCardinality1 = basis1_->getCardinality();
1792 int basisCardinality2 = basis2_->getCardinality();
1794 int totalPointCount = tensorPoints ? inputPoints1.extent_int(0) * inputPoints2.extent_int(0) : inputPoints1.extent_int(0);
1796 int pointCount1, pointCount2;
1799 pointCount1 = inputPoints1.extent_int(0);
1800 pointCount2 = inputPoints2.extent_int(0);
1804 pointCount1 = totalPointCount;
1805 pointCount2 = totalPointCount;
1808 const ordinal_type spaceDim1 = inputPoints1.extent_int(1);
1809 const ordinal_type spaceDim2 = inputPoints2.extent_int(1);
1811 INTREPID2_TEST_FOR_EXCEPTION(!tensorPoints && (totalPointCount != inputPoints2.extent_int(0)),
1812 std::invalid_argument,
"If tensorPoints is false, the point counts must match!");
1814 const ordinal_type opRank1 =
getOperatorRank(basis1_->getFunctionSpace(), operatorType1, spaceDim1);
1815 const ordinal_type opRank2 =
getOperatorRank(basis2_->getFunctionSpace(), operatorType2, spaceDim2);
1817 const ordinal_type outputRank1 = opRank1 +
getFieldRank(basis1_->getFunctionSpace());
1818 const ordinal_type outputRank2 = opRank2 +
getFieldRank(basis2_->getFunctionSpace());
1820 OutputViewType outputValues1, outputValues2;
1821 if (outputRank1 == 0)
1825 else if (outputRank1 == 1)
1827 outputValues1 =
getMatchingViewWithLabel(outputValues,
"output values - basis 1",basisCardinality1,pointCount1,spaceDim1);
1831 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported opRank1");
1834 if (outputRank2 == 0)
1838 else if (outputRank2 == 1)
1840 outputValues2 =
getMatchingViewWithLabel(outputValues,
"output values - basis 2",basisCardinality2,pointCount2,spaceDim2);
1844 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported opRank2");
1847 basis1_->getValues(outputValues1,inputPoints1,operatorType1);
1848 basis2_->getValues(outputValues2,inputPoints2,operatorType2);
1850 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputValueType>();
1851 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointValueType>();
1852 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
1854 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(basisCardinality1,Kokkos::AUTO(),vectorSize);
1858 FunctorType functor(outputValues, outputValues1, outputValues2, tensorPoints, weight);
1859 Kokkos::parallel_for(
"TensorViewFunctor", policy , functor);
1868 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"TensorBasis subclasses must override getHostBasis");
1879 template<
class ExecutionSpace,
class OutputScalar,
class OutputFieldType>
1882 using ScratchSpace =
typename ExecutionSpace::scratch_memory_space;
1883 using OutputScratchView = Kokkos::View<OutputScalar*,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
1885 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
1886 using TeamMember =
typename TeamPolicy::member_type;
1888 OutputFieldType output_;
1889 OutputFieldType input1_;
1890 OutputFieldType input2_;
1891 OutputFieldType input3_;
1893 int numFields_, numPoints_;
1894 int numFields1_, numPoints1_;
1895 int numFields2_, numPoints2_;
1896 int numFields3_, numPoints3_;
1902 TensorBasis3_Functor(OutputFieldType output, OutputFieldType inputValues1, OutputFieldType inputValues2, OutputFieldType inputValues3,
1903 bool tensorPoints,
double weight)
1904 : output_(output), input1_(inputValues1), input2_(inputValues2), input3_(inputValues3), tensorPoints_(tensorPoints), weight_(weight)
1906 numFields_ = output.extent_int(0);
1907 numPoints_ = output.extent_int(1);
1909 numFields1_ = inputValues1.extent_int(0);
1910 numPoints1_ = inputValues1.extent_int(1);
1912 numFields2_ = inputValues2.extent_int(0);
1913 numPoints2_ = inputValues2.extent_int(1);
1915 numFields3_ = inputValues3.extent_int(0);
1916 numPoints3_ = inputValues3.extent_int(1);
1923 INTREPID2_TEST_FOR_EXCEPTION(inputValues1.rank() > 3, std::invalid_argument,
"ranks greater than 3 not yet supported");
1924 INTREPID2_TEST_FOR_EXCEPTION(inputValues2.rank() > 3, std::invalid_argument,
"ranks greater than 3 not yet supported");
1925 INTREPID2_TEST_FOR_EXCEPTION(inputValues3.rank() > 3, std::invalid_argument,
"ranks greater than 3 not yet supported");
1926 INTREPID2_TEST_FOR_EXCEPTION((inputValues1.rank() == 3) && (inputValues2.rank() == 3), std::invalid_argument,
"two vector-valued input ranks not yet supported");
1927 INTREPID2_TEST_FOR_EXCEPTION((inputValues1.rank() == 3) && (inputValues3.rank() == 3), std::invalid_argument,
"two vector-valued input ranks not yet supported");
1928 INTREPID2_TEST_FOR_EXCEPTION((inputValues2.rank() == 3) && (inputValues3.rank() == 3), std::invalid_argument,
"two vector-valued input ranks not yet supported");
1933 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != numPoints1_, std::invalid_argument,
"incompatible point counts");
1934 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != numPoints2_, std::invalid_argument,
"incompatible point counts");
1935 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != numPoints3_, std::invalid_argument,
"incompatible point counts");
1939 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != numPoints1_ * numPoints2_ * numPoints3_, std::invalid_argument,
"incompatible point counts");
1942 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != numFields1_ * numFields2_ * numFields3_, std::invalid_argument,
"incompatible field sizes");
1945 KOKKOS_INLINE_FUNCTION
1946 void operator()(
const TeamMember & teamMember )
const
1948 auto fieldOrdinal1 = teamMember.league_rank();
1952 if ((input1_.rank() == 2) && (input2_.rank() == 2) && (input3_.rank() == 2))
1954 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (
const int& fieldOrdinal2) {
1955 for (
int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
1957 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
1958 for (
int pointOrdinal=0; pointOrdinal<numPoints_; pointOrdinal++)
1960 output_(fieldOrdinal,pointOrdinal) = weight_ * input1_(fieldOrdinal1,pointOrdinal) * input2_(fieldOrdinal2,pointOrdinal) * input3_(fieldOrdinal3,pointOrdinal);
1965 else if (input1_.rank() == 3)
1967 int spaceDim = input1_.extent_int(2);
1968 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (
const int& fieldOrdinal2) {
1969 for (
int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
1971 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
1972 for (
int pointOrdinal=0; pointOrdinal<numPoints_; pointOrdinal++)
1974 for (
int d=0; d<spaceDim; d++)
1976 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal,d) * input2_(fieldOrdinal2,pointOrdinal) * input3_(fieldOrdinal3,pointOrdinal);
1982 else if (input2_.rank() == 3)
1984 int spaceDim = input2_.extent_int(2);
1985 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (
const int& fieldOrdinal2) {
1986 for (
int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
1988 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
1989 for (
int pointOrdinal=0; pointOrdinal<numPoints_; pointOrdinal++)
1991 for (
int d=0; d<spaceDim; d++)
1993 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal) * input2_(fieldOrdinal2,pointOrdinal,d) * input3_(fieldOrdinal3,pointOrdinal);
1999 else if (input3_.rank() == 3)
2001 int spaceDim = input3_.extent_int(2);
2002 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (
const int& fieldOrdinal2) {
2003 for (
int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
2005 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
2006 for (
int pointOrdinal=0; pointOrdinal<numPoints_; pointOrdinal++)
2008 for (
int d=0; d<spaceDim; d++)
2010 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal) * input2_(fieldOrdinal2,pointOrdinal) * input3_(fieldOrdinal3,pointOrdinal,d);
2023 if ((input1_.rank() == 2) && (input2_.rank() == 2) && (input3_.rank() == 2) )
2025 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (
const int& fieldOrdinal2) {
2026 for (
int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
2028 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
2029 for (
int pointOrdinal3=0; pointOrdinal3<numPoints3_; pointOrdinal3++)
2031 for (
int pointOrdinal2=0; pointOrdinal2<numPoints2_; pointOrdinal2++)
2033 for (
int pointOrdinal1=0; pointOrdinal1<numPoints1_; pointOrdinal1++)
2035 int pointOrdinal = (pointOrdinal3 * numPoints2_ + pointOrdinal2) * numPoints1_ + pointOrdinal1;
2036 output_(fieldOrdinal,pointOrdinal) = weight_ * input1_(fieldOrdinal1,pointOrdinal1) * input2_(fieldOrdinal2,pointOrdinal2) * input3_(fieldOrdinal3,pointOrdinal3);
2043 else if (input1_.rank() == 3)
2045 int spaceDim = input1_.extent_int(2);
2046 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (
const int& fieldOrdinal2) {
2047 for (
int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
2049 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
2050 for (
int pointOrdinal3=0; pointOrdinal3<numPoints3_; pointOrdinal3++)
2052 for (
int pointOrdinal2=0; pointOrdinal2<numPoints2_; pointOrdinal2++)
2054 for (
int pointOrdinal1=0; pointOrdinal1<numPoints1_; pointOrdinal1++)
2056 int pointOrdinal = (pointOrdinal3 * numPoints2_ + pointOrdinal2) * numPoints1_ + pointOrdinal1;
2057 for (
int d=0; d<spaceDim; d++)
2059 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal1,d) * input2_(fieldOrdinal2,pointOrdinal2) * input3_(fieldOrdinal3,pointOrdinal3);
2067 else if (input2_.rank() == 3)
2069 int spaceDim = input2_.extent_int(2);
2070 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (
const int& fieldOrdinal2) {
2071 for (
int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
2073 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
2074 for (
int pointOrdinal3=0; pointOrdinal3<numPoints3_; pointOrdinal3++)
2076 for (
int pointOrdinal2=0; pointOrdinal2<numPoints2_; pointOrdinal2++)
2078 for (
int pointOrdinal1=0; pointOrdinal1<numPoints1_; pointOrdinal1++)
2080 int pointOrdinal = (pointOrdinal3 * numPoints2_ + pointOrdinal2) * numPoints1_ + pointOrdinal1;
2081 for (
int d=0; d<spaceDim; d++)
2083 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal1) * input2_(fieldOrdinal2,pointOrdinal2,d) * input3_(fieldOrdinal3,pointOrdinal3);
2091 else if (input3_.rank() == 3)
2093 int spaceDim = input3_.extent_int(2);
2094 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (
const int& fieldOrdinal2) {
2095 for (
int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
2097 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
2098 for (
int pointOrdinal3=0; pointOrdinal3<numPoints3_; pointOrdinal3++)
2100 for (
int pointOrdinal2=0; pointOrdinal2<numPoints2_; pointOrdinal2++)
2102 for (
int pointOrdinal1=0; pointOrdinal1<numPoints1_; pointOrdinal1++)
2104 int pointOrdinal = (pointOrdinal3 * numPoints2_ + pointOrdinal2) * numPoints1_ + pointOrdinal1;
2105 for (
int d=0; d<spaceDim; d++)
2107 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal1) * input2_(fieldOrdinal2,pointOrdinal2) * input3_(fieldOrdinal3,pointOrdinal3,d);
2124 template<
typename BasisBaseClass =
void>
2128 using BasisBase = BasisBaseClass;
2131 using typename BasisBase::OutputViewType;
2132 using typename BasisBase::PointViewType;
2133 using typename BasisBase::ScalarViewType;
2135 using typename BasisBase::OutputValueType;
2136 using typename BasisBase::PointValueType;
2138 using typename BasisBase::ExecutionSpace;
2140 using BasisPtr = Teuchos::RCP<BasisBase>;
2146 Basis_TensorBasis3(BasisPtr basis1, BasisPtr basis2, BasisPtr basis3,
const bool useShardsCellTopologyAndTags =
false)
2148 TensorBasis(Teuchos::rcp(
new TensorBasis(basis1,basis2,FUNCTION_SPACE_MAX,useShardsCellTopologyAndTags)),
2150 FUNCTION_SPACE_MAX,useShardsCellTopologyAndTags),
2165 std::vector<BasisPtr> componentBases {basis1_, basis2_, basis3_};
2195 virtual void getValues(OutputViewType outputValues,
const EOperator operatorType,
2196 const PointViewType inputPoints12,
const PointViewType inputPoints3,
2197 bool tensorPoints)
const override
2201 int spaceDim1 = basis1_->getDomainDimension();
2202 int spaceDim2 = basis2_->getDomainDimension();
2204 int totalSpaceDim12 = inputPoints12.extent_int(1);
2206 TEUCHOS_ASSERT(spaceDim1 + spaceDim2 == totalSpaceDim12);
2210 auto inputPoints1 = Kokkos::subview(inputPoints12,Kokkos::ALL(),std::make_pair(0,spaceDim1));
2211 auto inputPoints2 = Kokkos::subview(inputPoints12,Kokkos::ALL(),std::make_pair(spaceDim1,totalSpaceDim12));
2213 this->
getValues(outputValues, operatorType, inputPoints1, inputPoints2, inputPoints3, tensorPoints);
2220 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"This method does not yet handle tensorPoints=true");
2251 virtual void getValues(OutputViewType outputValues,
const EOperator operatorType,
2252 const PointViewType inputPoints1,
const PointViewType inputPoints2,
const PointViewType inputPoints3,
2253 bool tensorPoints)
const = 0;
2283 const PointViewType inputPoints1,
const EOperator operatorType1,
2284 const PointViewType inputPoints2,
const EOperator operatorType2,
2285 const PointViewType inputPoints3,
const EOperator operatorType3,
2286 bool tensorPoints,
double weight=1.0)
const
2288 int basisCardinality1 = basis1_->getCardinality();
2289 int basisCardinality2 = basis2_->getCardinality();
2290 int basisCardinality3 = basis3_->getCardinality();
2292 int spaceDim1 = inputPoints1.extent_int(1);
2293 int spaceDim2 = inputPoints2.extent_int(1);
2294 int spaceDim3 = inputPoints3.extent_int(1);
2296 int totalPointCount;
2297 int pointCount1, pointCount2, pointCount3;
2300 pointCount1 = inputPoints1.extent_int(0);
2301 pointCount2 = inputPoints2.extent_int(0);
2302 pointCount3 = inputPoints3.extent_int(0);
2303 totalPointCount = pointCount1 * pointCount2 * pointCount3;
2307 totalPointCount = inputPoints1.extent_int(0);
2308 pointCount1 = totalPointCount;
2309 pointCount2 = totalPointCount;
2310 pointCount3 = totalPointCount;
2312 INTREPID2_TEST_FOR_EXCEPTION((totalPointCount != inputPoints2.extent_int(0)) || (totalPointCount != inputPoints3.extent_int(0)),
2313 std::invalid_argument,
"If tensorPoints is false, the point counts must match!");
2332 OutputViewType outputValues1, outputValues2, outputValues3;
2336 if ((spaceDim1 == 1) && (operatorType1 == OPERATOR_VALUE))
2343 outputValues1 =
getMatchingViewWithLabel(outputValues,
"output values - basis 1",basisCardinality1,pointCount1,spaceDim1);
2345 if ((spaceDim2 == 1) && (operatorType2 == OPERATOR_VALUE))
2352 outputValues2 =
getMatchingViewWithLabel(outputValues,
"output values - basis 2",basisCardinality2,pointCount2,spaceDim2);
2354 if ((spaceDim3 == 1) && (operatorType3 == OPERATOR_VALUE))
2361 outputValues3 =
getMatchingViewWithLabel(outputValues,
"output values - basis 3",basisCardinality3,pointCount3,spaceDim3);
2364 basis1_->getValues(outputValues1,inputPoints1,operatorType1);
2365 basis2_->getValues(outputValues2,inputPoints2,operatorType2);
2366 basis3_->getValues(outputValues3,inputPoints3,operatorType3);
2368 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputValueType>();
2369 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointValueType>();
2370 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
2372 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(basisCardinality1,Kokkos::AUTO(),vectorSize);
2375 FunctorType functor(outputValues, outputValues1, outputValues2, outputValues3, tensorPoints, weight);
2376 Kokkos::parallel_for(
"TensorBasis3_Functor", policy , functor);
2385 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"TensorBasis3 subclasses must override getHostBasis");
Header file for the abstract base class Intrepid2::Basis.
BasisPtr< typename Kokkos::HostSpace::device_type, OutputType, PointType > HostBasisPtr
Pointer to a Basis whose device type is on the host (Kokkos::HostSpace::device_type),...
KOKKOS_INLINE_FUNCTION ordinal_type getOperatorRank(const EFunctionSpace spaceType, const EOperator operatorType, const ordinal_type spaceDim)
Returns rank of an operator.
KOKKOS_INLINE_FUNCTION ordinal_type getFieldRank(const EFunctionSpace spaceType)
Returns the rank of fields in a function space of the specified type.
KOKKOS_INLINE_FUNCTION ordinal_type getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
KOKKOS_INLINE_FUNCTION ordinal_type getDkCardinality(const EOperator operatorType, const ordinal_type spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
KOKKOS_INLINE_FUNCTION ordinal_type getDkEnumeration(const ordinal_type xMult, const ordinal_type yMult=-1, const ordinal_type zMult=-1)
Returns the ordinal of a partial derivative of order k based on the multiplicities of the partials dx...
Implements arbitrary-dimensional extrusion of a base shards::CellTopology.
@ GENERAL
arbitrary variation
Implementation of an assert that can safely be called from device code.
Class that defines mappings from component cell topologies to their tensor product topologies.
Implementation of support for traversing component views alongside a view that represents a combinati...
Header function for Intrepid2::Util class and other utility functions.
#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.
The data containers in Intrepid2 that support sum factorization and other reduced-data optimizations ...
const VectorDataType & vectorData() const
VectorData accessor.
KOKKOS_INLINE_FUNCTION int numFamilies() const
For valid vectorData, returns the number of families in vectorData; otherwise, returns number of Tens...
TensorDataType & tensorData()
TensorData accessor for single-family scalar data.
virtual HostBasisPtr< OutputValueType, PointValueType > getHostBasis() const override
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
virtual void getValues(OutputViewType outputValues, const EOperator operatorType, const PointViewType inputPoints12, const PointViewType inputPoints3, bool tensorPoints) const override
Evaluation of a tensor FEM basis on a reference cell.
virtual void getValues(OutputViewType outputValues, const EOperator operatorType, const PointViewType inputPoints1, const PointViewType inputPoints2, const PointViewType inputPoints3, bool tensorPoints) const =0
Evaluation of a tensor FEM basis on a reference cell; subclasses should override this.
void getValues(OutputViewType outputValues, const PointViewType inputPoints1, const EOperator operatorType1, const PointViewType inputPoints2, const EOperator operatorType2, const PointViewType inputPoints3, const EOperator operatorType3, bool tensorPoints, double weight=1.0) const
Evaluation of a tensor FEM basis on a reference cell; subclasses should override this.
virtual OperatorTensorDecomposition getOperatorDecomposition(const EOperator operatorType) const override
Returns a full decomposition of the specified operator. (Full meaning that all TensorBasis components...
Basis defined as the tensor product of two component bases.
virtual OperatorTensorDecomposition getSimpleOperatorDecomposition(const EOperator &operatorType) const
Returns a simple decomposition of the specified operator: what operator(s) should be applied to basis...
virtual OperatorTensorDecomposition getOperatorDecomposition(const EOperator operatorType) const
Returns a full decomposition of the specified operator. (Full meaning that all TensorBasis components...
virtual void getDofCoords(typename BasisBase::ScalarViewType dofCoords) const override
Fills in spatial locations (coordinates) of degrees of freedom (nodes) on the reference cell.
virtual void getDofCoeffs(typename BasisBase::ScalarViewType dofCoeffs) const override
Fills in coefficients of degrees of freedom on the reference cell.
virtual BasisValues< OutputValueType, DeviceType > allocateBasisValues(TensorPoints< PointValueType, DeviceType > points, const EOperator operatorType=OPERATOR_VALUE) const override
Allocate BasisValues container suitable for passing to the getValues() variant that takes a TensorPoi...
virtual void getValues(OutputViewType outputValues, const EOperator operatorType, const PointViewType inputPoints1, const PointViewType inputPoints2, bool tensorPoints) const
Evaluation of a tensor FEM basis on a reference cell; subclasses should override this.
virtual const char * getName() const override
Returns basis name.
virtual HostBasisPtr< OutputValueType, PointValueType > getHostBasis() const override
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
void getValues(OutputViewType outputValues, const PointViewType inputPoints1, const EOperator operatorType1, const PointViewType inputPoints2, const EOperator operatorType2, bool tensorPoints, double weight=1.0) const
Evaluation of a tensor FEM basis on a reference cell.
ordinal_type getTensorDkEnumeration(ordinal_type dkEnum1, ordinal_type operatorOrder1, ordinal_type dkEnum2, ordinal_type operatorOrder2) const
Given "Dk" enumeration indices for the component bases, returns a Dk enumeration index for the compos...
Basis_TensorBasis(BasisPtr basis1, BasisPtr basis2, EFunctionSpace functionSpace=FUNCTION_SPACE_MAX, const bool useShardsCellTopologyAndTags=false)
Constructor.
void getComponentPoints(const PointViewType inputPoints, const bool attemptTensorDecomposition, PointViewType &inputPoints1, PointViewType &inputPoints2, bool &tensorDecompositionSucceeded) const
Method to extract component points from composite points.
void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
virtual void getValues(BasisValues< OutputValueType, DeviceType > outputValues, const TensorPoints< PointValueType, DeviceType > inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell, using point and output value containers that allow pre...
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
static CellTopoPtr cellTopology(const shards::CellTopology &shardsCellTopo, ordinal_type tensorialDegree=0)
static accessor that returns a CellTopoPtr; these are lazily constructed and cached.
static ordinal_type getSubcellOrdinalMap(CellTopoPtr cellTopo, ordinal_type subcdim, ordinal_type subcord, ordinal_type subsubcdim, ordinal_type subsubcord)
Maps the from a subcell within a subcell of the present CellTopology to the subcell in the present Ce...
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.
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 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 ordinal_type numTensorComponents() const
Returns the number of tensorial components.
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.
ordinal_type componentPointCount(const ordinal_type &tensorComponentOrdinal) const
Returns the number of points in the indicated component.
Functor for computing values for the TensorBasis class.
A helper class that allows iteration over three Kokkos Views simultaneously, according to tensor comb...
KOKKOS_INLINE_FUNCTION ScalarType getView1Entry()
KOKKOS_INLINE_FUNCTION int nextIncrementRank()
KOKKOS_INLINE_FUNCTION void setLocation(const Kokkos::Array< int, 7 > location)
KOKKOS_INLINE_FUNCTION int increment()
KOKKOS_INLINE_FUNCTION void set(ScalarType value)
KOKKOS_INLINE_FUNCTION ScalarType getView2Entry()
Reference-space field values for a basis, designed to support typical vector-valued bases.
KOKKOS_INLINE_FUNCTION const TensorData< Scalar, DeviceType > & getComponent(const int &componentOrdinal) const
Single-argument component accessor for the axial-component or the single-family case; in this case,...
For a multi-component tensor basis, specifies the operators to be applied to the components to produc...
bool rotateXYNinetyDegrees() const
If true, this flag indicates that a vector component that spans the first two dimensions should be ro...
OperatorTensorDecomposition expandedDecomposition(std::vector< Teuchos::RCP< Basis< DeviceType, OutputValueType, PointValueType > > > &bases)
takes as argument bases that are components in this decomposition, and decomposes them further if the...
Functor for computing values for the TensorBasis3 class.