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
74 ordinal_type getDkEnumeration(
const Kokkos::Array<int,spaceDim> &entries);
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>(
const Kokkos::Array<int,1> &entries)
169 return getDkEnumeration<1>(entries[0]);
172 template<ordinal_type spaceDim>
173 KOKKOS_INLINE_FUNCTION
174 ordinal_type getDkEnumeration(
const Kokkos::Array<int,spaceDim> &entries)
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);
228 template<
typename BasisBase>
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 for (ordinal_type d2=0; d2<=d2_max; d2++)
908 ordinal_type d1 = d-d2;
909 if (d1 > spaceDim1)
continue;
911 ordinal_type subcellCount2 = cellTopo2.getSubcellCount(d2);
912 ordinal_type subcellCount1 = cellTopo1.getSubcellCount(d1);
913 for (ordinal_type subcellOrdinal2=0; subcellOrdinal2<subcellCount2; subcellOrdinal2++)
915 ordinal_type subcellDofCount2 = basis2_->getDofCount(d2, subcellOrdinal2);
916 for (ordinal_type subcellOrdinal1=0; subcellOrdinal1<subcellCount1; subcellOrdinal1++)
918 ordinal_type subcellDofCount1 = basis1_->getDofCount(d1, subcellOrdinal1);
919 ordinal_type tensorLocalDofCount = subcellDofCount1 * subcellDofCount2;
920 for (ordinal_type localDofID2 = 0; localDofID2<subcellDofCount2; localDofID2++)
922 ordinal_type fieldOrdinal2 = basis2_->getDofOrdinal(d2, subcellOrdinal2, localDofID2);
923 OrdinalTypeArray1DHost degreesField2;
924 if (this->basisType_ == BASIS_FEM_HIERARCHICAL) degreesField2 = basis2_->getPolynomialDegreeOfField(fieldOrdinal2);
925 for (ordinal_type localDofID1 = 0; localDofID1<subcellDofCount1; localDofID1++)
927 ordinal_type fieldOrdinal1 = basis1_->getDofOrdinal(d1, subcellOrdinal1, localDofID1);
928 ordinal_type tensorLocalDofID = localDofID2 * subcellDofCount1 + localDofID1;
929 ordinal_type tensorFieldOrdinal = fieldOrdinal2 * basis1_->getCardinality() + fieldOrdinal1;
930 tagView(tensorFieldOrdinal*tagSize+0) = d;
931 tagView(tensorFieldOrdinal*tagSize+1) = topoMap.getCompositeSubcellOrdinal(d1, subcellOrdinal1, d2, subcellOrdinal2);
932 tagView(tensorFieldOrdinal*tagSize+2) = tensorLocalDofID;
933 tagView(tensorFieldOrdinal*tagSize+3) = tensorLocalDofCount;
943 this->setOrdinalTagData(this->tagToOrdinal_,
946 this->basisCardinality_,
954 virtual int getNumTensorialExtrusions()
const override
956 return numTensorialExtrusions_;
968 ordinal_type dkEnum2, ordinal_type operatorOrder2)
const
970 ordinal_type spaceDim1 = basis1_->getDomainDimension();
971 ordinal_type spaceDim2 = basis2_->getDomainDimension();
978 INTREPID2_TEST_FOR_EXCEPTION(operatorOrder1 > 0, std::invalid_argument,
"For spaceDim1 = 0, operatorOrder1 must be 0.");
984 case 1:
return getDkTensorIndex<1, 1>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
985 case 2:
return getDkTensorIndex<1, 2>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
986 case 3:
return getDkTensorIndex<1, 3>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
987 case 4:
return getDkTensorIndex<1, 4>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
988 case 5:
return getDkTensorIndex<1, 5>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
989 case 6:
return getDkTensorIndex<1, 6>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
991 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported dimension combination");
996 case 1:
return getDkTensorIndex<2, 1>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
997 case 2:
return getDkTensorIndex<2, 2>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
998 case 3:
return getDkTensorIndex<2, 3>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
999 case 4:
return getDkTensorIndex<2, 4>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1000 case 5:
return getDkTensorIndex<2, 5>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1002 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported dimension combination");
1007 case 1:
return getDkTensorIndex<3, 1>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1008 case 2:
return getDkTensorIndex<3, 2>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1009 case 3:
return getDkTensorIndex<3, 3>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1010 case 4:
return getDkTensorIndex<3, 4>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1012 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported dimension combination");
1017 case 1:
return getDkTensorIndex<4, 1>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1018 case 2:
return getDkTensorIndex<4, 2>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1019 case 3:
return getDkTensorIndex<4, 3>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1021 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported dimension combination");
1026 case 1:
return getDkTensorIndex<5, 1>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1027 case 2:
return getDkTensorIndex<5, 2>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1029 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported dimension combination");
1034 case 1:
return getDkTensorIndex<6, 1>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1036 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported dimension combination");
1039 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported dimension combination");
1049 const int spaceDim = this->getDomainDimension();
1051 const EOperator VALUE = Intrepid2::OPERATOR_VALUE;
1053 std::vector< std::vector<EOperator> > opsVALUE{{VALUE, VALUE}};
1055 std::vector< std::vector<EOperator> > ops(spaceDim);
1057 switch (operatorType)
1065 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported operator type - TensorBasis subclass should override");
1080 auto opOrder = getOperatorOrder(operatorType);
1081 const int dkCardinality = getDkCardinality(operatorType, 2);
1083 ops = std::vector< std::vector<EOperator> >(dkCardinality);
1087 for (
int derivativeCountComp2=0; derivativeCountComp2<=opOrder; derivativeCountComp2++)
1089 int derivativeCountComp1=opOrder-derivativeCountComp2;
1090 EOperator op1 = (derivativeCountComp1 == 0) ? OPERATOR_VALUE : EOperator(OPERATOR_D1 + (derivativeCountComp1 - 1));
1091 EOperator op2 = (derivativeCountComp2 == 0) ? OPERATOR_VALUE : EOperator(OPERATOR_D1 + (derivativeCountComp2 - 1));
1093 int dkCardinality1 = getDkCardinality(op1, 1);
1094 int dkCardinality2 = getDkCardinality(op2, 1);
1096 for (
int dkEnum1=0; dkEnum1<dkCardinality1; dkEnum1++)
1098 for (
int dkEnum2=0; dkEnum2<dkCardinality2; dkEnum2++)
1100 ordinal_type dkTensorIndex = getDkTensorIndex<1, 1>(dkEnum1, derivativeCountComp1, dkEnum2, derivativeCountComp2);
1101 ops[dkTensorIndex] = std::vector<EOperator>{op1, op2};
1109 std::vector<double> weights(ops.size(), 1.0);
1117 if (((operatorType >= OPERATOR_D1) && (operatorType <= OPERATOR_D10)) || (operatorType == OPERATOR_GRAD))
1124 ordinal_type numBasisComponents = tensorComponents_.size();
1126 auto opOrder = getOperatorOrder(operatorType);
1127 const int dkCardinality = getDkCardinality(operatorType, numBasisComponents);
1129 std::vector< std::vector<EOperator> > ops(dkCardinality);
1131 std::vector<EOperator> prevEntry(numBasisComponents, OPERATOR_VALUE);
1132 prevEntry[0] = operatorType;
1136 for (ordinal_type dkOrdinal=1; dkOrdinal<dkCardinality; dkOrdinal++)
1138 std::vector<EOperator> entry = prevEntry;
1147 ordinal_type cumulativeOpOrder = 0;
1148 ordinal_type finalOpOrder = getOperatorOrder(entry[numBasisComponents-1]);
1149 for (ordinal_type compOrdinal=0; compOrdinal<numBasisComponents; compOrdinal++)
1151 const ordinal_type thisOpOrder = getOperatorOrder(entry[compOrdinal]);
1152 cumulativeOpOrder += thisOpOrder;
1153 if (cumulativeOpOrder + finalOpOrder == opOrder)
1156 EOperator decrementedOp;
1157 if (thisOpOrder == 1)
1159 decrementedOp = OPERATOR_VALUE;
1163 decrementedOp =
static_cast<EOperator
>(OPERATOR_D1 + ((thisOpOrder - 1) - 1));
1165 entry[compOrdinal] = decrementedOp;
1166 const ordinal_type remainingOpOrder = opOrder - cumulativeOpOrder + 1;
1167 entry[compOrdinal+1] =
static_cast<EOperator
>(OPERATOR_D1 + (remainingOpOrder - 1));
1168 for (ordinal_type i=compOrdinal+2; i<numBasisComponents; i++)
1170 entry[i] = OPERATOR_VALUE;
1175 ops[dkOrdinal] = entry;
1178 std::vector<double> weights(dkCardinality, 1.0);
1185 std::vector<BasisPtr> componentBases {basis1_, basis2_};
1196 const bool operatorIsDk = (operatorType >= OPERATOR_D1) && (operatorType <= OPERATOR_D10);
1197 const bool operatorSupported = (operatorType == OPERATOR_VALUE) || (operatorType == OPERATOR_GRAD) || (operatorType == OPERATOR_CURL) || (operatorType == OPERATOR_DIV) || operatorIsDk;
1198 INTREPID2_TEST_FOR_EXCEPTION(!operatorSupported, std::invalid_argument,
"operator is not supported by allocateBasisValues");
1201 const int spaceDim = this->getDomainDimension();
1202 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");
1205 ordinal_type numBasisComponents = tensorComponents_.size();
1211 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.");
1212 const ordinal_type numPoints = points.
extent_int(0);
1213 auto outputView = this->allocateOutputView(numPoints, operatorType);
1220 INTREPID2_TEST_FOR_EXCEPTION(numBasisComponents > points.
numTensorComponents(), std::invalid_argument,
"points must have at least as many tensorial components as basis.");
1224 ordinal_type numVectorComponents = opDecomposition.numVectorComponents();
1225 const bool useVectorData = numVectorComponents > 1;
1227 std::vector<ordinal_type> componentPointCounts(numBasisComponents);
1228 ordinal_type pointComponentNumber = 0;
1229 for (ordinal_type r=0; r<numBasisComponents; r++)
1231 const ordinal_type compSpaceDim = tensorComponents_[r]->getDomainDimension();
1232 ordinal_type dimsSoFar = 0;
1233 ordinal_type numPointsForBasisComponent = 1;
1234 while (dimsSoFar < compSpaceDim)
1236 INTREPID2_TEST_FOR_EXCEPTION(pointComponentNumber >= points.
numTensorComponents(), std::invalid_argument,
"Error in processing points container; perhaps it is mis-sized?");
1238 const int numComponentDims = points.
getTensorComponent(pointComponentNumber).extent_int(1);
1239 numPointsForBasisComponent *= numComponentPoints;
1240 dimsSoFar += numComponentDims;
1241 INTREPID2_TEST_FOR_EXCEPTION(dimsSoFar > points.
numTensorComponents(), std::invalid_argument,
"Error in processing points container; perhaps it is mis-sized?");
1242 pointComponentNumber++;
1244 componentPointCounts[r] = numPointsForBasisComponent;
1249 const int numFamilies = 1;
1252 const int familyOrdinal = 0;
1253 for (ordinal_type vectorComponentOrdinal=0; vectorComponentOrdinal<numVectorComponents; vectorComponentOrdinal++)
1255 if (!opDecomposition.identicallyZeroComponent(vectorComponentOrdinal))
1257 std::vector< Data<OutputValueType,DeviceType> > componentData;
1258 for (ordinal_type r=0; r<numBasisComponents; r++)
1260 const int numComponentPoints = componentPointCounts[r];
1261 const EOperator op = opDecomposition.op(vectorComponentOrdinal, r);
1262 auto componentView = tensorComponents_[r]->allocateOutputView(numComponentPoints, op);
1274 std::vector< Data<OutputValueType,DeviceType> > componentData;
1276 const ordinal_type vectorComponentOrdinal = 0;
1277 for (ordinal_type r=0; r<numBasisComponents; r++)
1279 const int numComponentPoints = componentPointCounts[r];
1280 const EOperator op = opDecomposition.op(vectorComponentOrdinal, r);
1281 auto componentView = tensorComponents_[r]->allocateOutputView(numComponentPoints, op);
1286 const Kokkos::Array<int,7> extents {componentView.extent_int(0), componentView.extent_int(1), 1,1,1,1,1};
1287 Kokkos::Array<DataVariationType,7> variationType {GENERAL, GENERAL, CONSTANT, CONSTANT, CONSTANT, CONSTANT, CONSTANT };
1293 std::vector< TensorData<OutputValueType,DeviceType> > tensorDataEntries {tensorData};
1301 using BasisBase::getValues;
1315 PointViewType & inputPoints1, PointViewType & inputPoints2,
bool &tensorDecompositionSucceeded)
const
1317 INTREPID2_TEST_FOR_EXCEPTION(attemptTensorDecomposition, std::invalid_argument,
"tensor decomposition not yet supported");
1331 int spaceDim1 = basis1_->getDomainDimension();
1332 int spaceDim2 = basis2_->getDomainDimension();
1334 int totalSpaceDim = inputPoints.extent_int(1);
1336 TEUCHOS_ASSERT(spaceDim1 + spaceDim2 == totalSpaceDim);
1340 inputPoints1 = Kokkos::subview(inputPoints,Kokkos::ALL(),std::make_pair(0,spaceDim1));
1341 inputPoints2 = Kokkos::subview(inputPoints,Kokkos::ALL(),std::make_pair(spaceDim1,totalSpaceDim));
1347 tensorDecompositionSucceeded =
false;
1358 virtual void getDofCoords(
typename BasisBase::ScalarViewType dofCoords )
const override
1360 int spaceDim1 = basis1_->getBaseCellTopology().getDimension();
1361 int spaceDim2 = basis2_->getBaseCellTopology().getDimension();
1363 using ValueType =
typename BasisBase::ScalarViewType::value_type;
1364 using ResultLayout =
typename DeduceLayout< typename BasisBase::ScalarViewType >::result_layout;
1365 using ViewType = Kokkos::DynRankView<ValueType, ResultLayout, DeviceType >;
1367 const ordinal_type basisCardinality1 = basis1_->getCardinality();
1368 const ordinal_type basisCardinality2 = basis2_->getCardinality();
1370 ViewType dofCoords1(
"dofCoords1",basisCardinality1,spaceDim1);
1371 ViewType dofCoords2(
"dofCoords2",basisCardinality2,spaceDim2);
1373 basis1_->getDofCoords(dofCoords1);
1374 basis2_->getDofCoords(dofCoords2);
1376 Kokkos::RangePolicy<ExecutionSpace> policy(0, basisCardinality2);
1377 Kokkos::parallel_for(policy, KOKKOS_LAMBDA (
const int fieldOrdinal2)
1379 for (
int fieldOrdinal1=0; fieldOrdinal1<basisCardinality1; fieldOrdinal1++)
1381 const ordinal_type fieldOrdinal = fieldOrdinal1 + fieldOrdinal2 * basisCardinality1;
1382 for (
int d1=0; d1<spaceDim1; d1++)
1384 dofCoords(fieldOrdinal,d1) = dofCoords1(fieldOrdinal1,d1);
1386 for (
int d2=0; d2<spaceDim2; d2++)
1388 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 const ordinal_type basisCardinality1 = basis1_->getCardinality();
1411 const ordinal_type basisCardinality2 = basis2_->getCardinality();
1413 bool isVectorBasis1 = getFieldRank(basis1_->getFunctionSpace()) == 1;
1414 bool isVectorBasis2 = getFieldRank(basis2_->getFunctionSpace()) == 1;
1416 INTREPID2_TEST_FOR_EXCEPTION(isVectorBasis1 && isVectorBasis2, std::invalid_argument,
"the case in which basis1 and basis2 are vector bases is not supported");
1418 int basisDim1 = isVectorBasis1 ? basis1_->getBaseCellTopology().getDimension() : 1;
1419 int basisDim2 = isVectorBasis2 ? basis2_->getBaseCellTopology().getDimension() : 1;
1421 auto dofCoeffs1 = isVectorBasis1 ? ViewType(
"dofCoeffs1",basis1_->getCardinality(), basisDim1) : ViewType(
"dofCoeffs1",basis1_->getCardinality());
1422 auto dofCoeffs2 = isVectorBasis2 ? ViewType(
"dofCoeffs2",basis2_->getCardinality(), basisDim2) : ViewType(
"dofCoeffs2",basis2_->getCardinality());
1424 basis1_->getDofCoeffs(dofCoeffs1);
1425 basis2_->getDofCoeffs(dofCoeffs2);
1427 Kokkos::RangePolicy<ExecutionSpace> policy(0, basisCardinality2);
1428 Kokkos::parallel_for(policy, KOKKOS_LAMBDA (
const int fieldOrdinal2)
1430 for (
int fieldOrdinal1=0; fieldOrdinal1<basisCardinality1; fieldOrdinal1++)
1432 const ordinal_type fieldOrdinal = fieldOrdinal1 + fieldOrdinal2 * basisCardinality1;
1433 for (
int d1 = 0; d1 <basisDim1; d1++) {
1434 for (
int d2 = 0; d2 <basisDim2; d2++) {
1435 dofCoeffs.access(fieldOrdinal,d1+d2) = dofCoeffs1.access(fieldOrdinal1,d1);
1436 dofCoeffs.access(fieldOrdinal,d1+d2) *= dofCoeffs2.access(fieldOrdinal2,d2);
1450 return name_.c_str();
1453 std::vector<BasisPtr> getTensorBasisComponents()
const
1455 return tensorComponents_;
1473 const EOperator operatorType = OPERATOR_VALUE )
const override
1475 const ordinal_type numTensorComponents = tensorComponents_.size();
1479 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" );
1480 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" );
1481 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" );
1483 OutputViewType outputView = outputValues.
tensorData().getTensorComponent(0).getUnderlyingView();
1485 this->
getValues(outputView, pointView, operatorType);
1491 const ordinal_type numVectorComponents = operatorDecomposition.numVectorComponents();
1492 const bool useVectorData = numVectorComponents > 1;
1493 const ordinal_type numBasisComponents = operatorDecomposition.numBasisComponents();
1495 for (ordinal_type vectorComponentOrdinal=0; vectorComponentOrdinal<numVectorComponents; vectorComponentOrdinal++)
1497 const double weight = operatorDecomposition.weight(vectorComponentOrdinal);
1498 ordinal_type pointComponentOrdinal = 0;
1499 for (ordinal_type basisOrdinal=0; basisOrdinal<numBasisComponents; basisOrdinal++, pointComponentOrdinal++)
1501 const EOperator op = operatorDecomposition.op(vectorComponentOrdinal, basisOrdinal);
1503 if (op != OPERATOR_MAX)
1505 const int vectorFamily = 0;
1506 auto tensorData = useVectorData ? outputValues.
vectorData().getComponent(vectorFamily,vectorComponentOrdinal) : outputValues.
tensorData();
1507 INTREPID2_TEST_FOR_EXCEPTION( ! tensorData.getTensorComponent(basisOrdinal).isValid(), std::invalid_argument,
"Invalid output component encountered");
1513 const ordinal_type basisDomainDimension = tensorComponents_[basisOrdinal]->getDomainDimension();
1514 if (pointView.extent_int(1) == basisDomainDimension)
1516 tensorComponents_[basisOrdinal]->getValues(basisValueView, pointView, op);
1523 ordinal_type dimsSoFar = 0;
1524 std::vector< ScalarView<PointValueType,DeviceType> > basisPointComponents;
1525 while (dimsSoFar < basisDomainDimension)
1527 INTREPID2_TEST_FOR_EXCEPTION(pointComponentOrdinal >= inputPoints.
numTensorComponents(), std::invalid_argument,
"Error in processing points container; perhaps it is mis-sized?");
1529 const ordinal_type numComponentDims = pointComponent.extent_int(1);
1530 dimsSoFar += numComponentDims;
1531 INTREPID2_TEST_FOR_EXCEPTION(dimsSoFar > inputPoints.
numTensorComponents(), std::invalid_argument,
"Error in processing points container; perhaps it is mis-sized?");
1532 basisPointComponents.push_back(pointComponent);
1533 if (dimsSoFar < basisDomainDimension)
1536 pointComponentOrdinal++;
1542 bool useVectorData2 = (basisValueView.rank() == 3);
1556 tensorComponents_[basisOrdinal]->getValues(basisValues, basisPoints, op);
1564 const bool spansXY = (vectorComponentOrdinal == 0) && (basisValueView.extent_int(2) == 2);
1568 auto policy = Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<2>>({0,0},{basisValueView.extent_int(0),basisValueView.extent_int(1)});
1569 Kokkos::parallel_for(
"rotateXYNinetyDegrees", policy,
1570 KOKKOS_LAMBDA (
const int &fieldOrdinal,
const int &pointOrdinal) {
1571 const auto f_x = basisValueView(fieldOrdinal,pointOrdinal,0);
1572 const auto &f_y = basisValueView(fieldOrdinal,pointOrdinal,1);
1573 basisValueView(fieldOrdinal,pointOrdinal,0) = -f_y;
1574 basisValueView(fieldOrdinal,pointOrdinal,1) = f_x;
1580 if ((weight != 1.0) && (basisOrdinal == 0))
1582 if (basisValueView.rank() == 2)
1584 auto policy = Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<2>>({0,0},{basisValueView.extent_int(0),basisValueView.extent_int(1)});
1585 Kokkos::parallel_for(
"multiply basisValueView by weight", policy,
1586 KOKKOS_LAMBDA (
const int &fieldOrdinal,
const int &pointOrdinal) {
1587 basisValueView(fieldOrdinal,pointOrdinal) *= weight;
1590 else if (basisValueView.rank() == 3)
1592 auto policy = Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<3>>({0,0,0},{basisValueView.extent_int(0),basisValueView.extent_int(1),basisValueView.extent_int(2)});
1593 Kokkos::parallel_for(
"multiply basisValueView by weight", policy,
1594 KOKKOS_LAMBDA (
const int &fieldOrdinal,
const int &pointOrdinal,
const int &d) {
1595 basisValueView(fieldOrdinal,pointOrdinal,d) *= weight;
1600 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported rank for basisValueView");
1626 void getValues( OutputViewType outputValues,
const PointViewType inputPoints,
1627 const EOperator operatorType = OPERATOR_VALUE )
const override
1630 bool attemptTensorDecomposition =
false;
1631 PointViewType inputPoints1, inputPoints2;
1632 getComponentPoints(inputPoints, attemptTensorDecomposition, inputPoints1, inputPoints2, tensorPoints);
1634 const auto functionSpace = this->getFunctionSpace();
1636 if ((functionSpace == FUNCTION_SPACE_HVOL) || (functionSpace == FUNCTION_SPACE_HGRAD))
1639 switch (operatorType)
1641 case OPERATOR_VALUE:
1654 auto opOrder = getOperatorOrder(operatorType);
1657 for (
int derivativeCountComp2=0; derivativeCountComp2<=opOrder; derivativeCountComp2++)
1659 int derivativeCountComp1=opOrder-derivativeCountComp2;
1660 EOperator op1 = (derivativeCountComp1 == 0) ? OPERATOR_VALUE : EOperator(OPERATOR_D1 + (derivativeCountComp1 - 1));
1661 EOperator op2 = (derivativeCountComp2 == 0) ? OPERATOR_VALUE : EOperator(OPERATOR_D1 + (derivativeCountComp2 - 1));
1663 int spaceDim1 = inputPoints1.extent_int(1);
1664 int spaceDim2 = inputPoints2.extent_int(1);
1666 int dkCardinality1 = (op1 != OPERATOR_VALUE) ? getDkCardinality(op1, spaceDim1) : 1;
1667 int dkCardinality2 = (op2 != OPERATOR_VALUE) ? getDkCardinality(op2, spaceDim2) : 1;
1669 int basisCardinality1 = basis1_->getCardinality();
1670 int basisCardinality2 = basis2_->getCardinality();
1672 int totalPointCount = tensorPoints ? inputPoints1.extent_int(0) * inputPoints2.extent_int(0) : inputPoints1.extent_int(0);
1674 int pointCount1, pointCount2;
1677 pointCount1 = inputPoints1.extent_int(0);
1678 pointCount2 = inputPoints2.extent_int(0);
1682 pointCount1 = totalPointCount;
1683 pointCount2 = totalPointCount;
1686 OutputViewType outputValues1, outputValues2;
1687 if (op1 == OPERATOR_VALUE)
1688 outputValues1 = getMatchingViewWithLabel(outputValues,
"output values - basis 1",basisCardinality1,pointCount1);
1690 outputValues1 = getMatchingViewWithLabel(outputValues,
"output values - basis 1",basisCardinality1,pointCount1,dkCardinality1);
1692 if (op2 == OPERATOR_VALUE)
1693 outputValues2 = getMatchingViewWithLabel(outputValues,
"output values - basis 2",basisCardinality2,pointCount2);
1695 outputValues2 = getMatchingViewWithLabel(outputValues,
"output values - basis 2",basisCardinality2,pointCount2,dkCardinality2);
1697 basis1_->getValues(outputValues1,inputPoints1,op1);
1698 basis2_->getValues(outputValues2,inputPoints2,op2);
1700 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputValueType>();
1701 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointValueType>();
1702 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
1704 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(basisCardinality1,Kokkos::AUTO(),vectorSize);
1706 double weight = 1.0;
1709 for (
int dkEnum1=0; dkEnum1<dkCardinality1; dkEnum1++)
1711 auto outputValues1_dkEnum1 = (op1 != OPERATOR_VALUE) ? Kokkos::subview(outputValues1,Kokkos::ALL(),Kokkos::ALL(),dkEnum1)
1712 : Kokkos::subview(outputValues1,Kokkos::ALL(),Kokkos::ALL());
1713 for (
int dkEnum2=0; dkEnum2<dkCardinality2; dkEnum2++)
1715 auto outputValues2_dkEnum2 = (op2 != OPERATOR_VALUE) ? Kokkos::subview(outputValues2,Kokkos::ALL(),Kokkos::ALL(),dkEnum2)
1716 : Kokkos::subview(outputValues2,Kokkos::ALL(),Kokkos::ALL());
1718 ordinal_type dkTensorIndex =
getTensorDkEnumeration(dkEnum1, derivativeCountComp1, dkEnum2, derivativeCountComp2);
1719 auto outputValues_dkTensor = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),dkTensorIndex);
1725 FunctorType functor(outputValues_dkTensor, outputValues1_dkEnum1, outputValues2_dkEnum2, tensorPoints, weight);
1726 Kokkos::parallel_for(
"TensorViewFunctor", policy , functor);
1733 this->
getValues(outputValues, operatorType, inputPoints1, inputPoints2, tensorPoints);
1739 this->
getValues(outputValues, operatorType, inputPoints1, inputPoints2, tensorPoints);
1768 virtual void getValues(OutputViewType outputValues,
const EOperator operatorType,
1769 const PointViewType inputPoints1,
const PointViewType inputPoints2,
1770 bool tensorPoints)
const
1772 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"one-operator, two-inputPoints getValues should be overridden by TensorBasis subclasses");
1799 const PointViewType inputPoints1,
const EOperator operatorType1,
1800 const PointViewType inputPoints2,
const EOperator operatorType2,
1801 bool tensorPoints,
double weight=1.0)
const
1803 int basisCardinality1 = basis1_->getCardinality();
1804 int basisCardinality2 = basis2_->getCardinality();
1806 int totalPointCount = tensorPoints ? inputPoints1.extent_int(0) * inputPoints2.extent_int(0) : inputPoints1.extent_int(0);
1808 int pointCount1, pointCount2;
1811 pointCount1 = inputPoints1.extent_int(0);
1812 pointCount2 = inputPoints2.extent_int(0);
1816 pointCount1 = totalPointCount;
1817 pointCount2 = totalPointCount;
1820 const ordinal_type spaceDim1 = inputPoints1.extent_int(1);
1821 const ordinal_type spaceDim2 = inputPoints2.extent_int(1);
1823 INTREPID2_TEST_FOR_EXCEPTION(!tensorPoints && (totalPointCount != inputPoints2.extent_int(0)),
1824 std::invalid_argument,
"If tensorPoints is false, the point counts must match!");
1826 const ordinal_type opRank1 = getOperatorRank(basis1_->getFunctionSpace(), operatorType1, spaceDim1);
1827 const ordinal_type opRank2 = getOperatorRank(basis2_->getFunctionSpace(), operatorType2, spaceDim2);
1829 const ordinal_type outputRank1 = opRank1 + getFieldRank(basis1_->getFunctionSpace());
1830 const ordinal_type outputRank2 = opRank2 + getFieldRank(basis2_->getFunctionSpace());
1832 OutputViewType outputValues1, outputValues2;
1833 if (outputRank1 == 0)
1835 outputValues1 = getMatchingViewWithLabel(outputValues,
"output values - basis 1",basisCardinality1,pointCount1);
1837 else if (outputRank1 == 1)
1839 outputValues1 = getMatchingViewWithLabel(outputValues,
"output values - basis 1",basisCardinality1,pointCount1,spaceDim1);
1843 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported opRank1");
1846 if (outputRank2 == 0)
1848 outputValues2 = getMatchingViewWithLabel(outputValues,
"output values - basis 2",basisCardinality2,pointCount2);
1850 else if (outputRank2 == 1)
1852 outputValues2 = getMatchingViewWithLabel(outputValues,
"output values - basis 2",basisCardinality2,pointCount2,spaceDim2);
1856 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unsupported opRank2");
1859 basis1_->getValues(outputValues1,inputPoints1,operatorType1);
1860 basis2_->getValues(outputValues2,inputPoints2,operatorType2);
1862 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputValueType>();
1863 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointValueType>();
1864 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
1866 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(basisCardinality1,Kokkos::AUTO(),vectorSize);
1870 FunctorType functor(outputValues, outputValues1, outputValues2, tensorPoints, weight);
1871 Kokkos::parallel_for(
"TensorViewFunctor", policy , functor);
1878 virtual HostBasisPtr<OutputValueType, PointValueType>
1880 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"TensorBasis subclasses must override getHostBasis");
1891 template<
class ExecutionSpace,
class OutputScalar,
class OutputFieldType>
1894 using ScratchSpace =
typename ExecutionSpace::scratch_memory_space;
1895 using OutputScratchView = Kokkos::View<OutputScalar*,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
1897 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
1898 using TeamMember =
typename TeamPolicy::member_type;
1900 OutputFieldType output_;
1901 OutputFieldType input1_;
1902 OutputFieldType input2_;
1903 OutputFieldType input3_;
1905 int numFields_, numPoints_;
1906 int numFields1_, numPoints1_;
1907 int numFields2_, numPoints2_;
1908 int numFields3_, numPoints3_;
1914 TensorBasis3_Functor(OutputFieldType output, OutputFieldType inputValues1, OutputFieldType inputValues2, OutputFieldType inputValues3,
1915 bool tensorPoints,
double weight)
1916 : output_(output), input1_(inputValues1), input2_(inputValues2), input3_(inputValues3), tensorPoints_(tensorPoints), weight_(weight)
1918 numFields_ = output.extent_int(0);
1919 numPoints_ = output.extent_int(1);
1921 numFields1_ = inputValues1.extent_int(0);
1922 numPoints1_ = inputValues1.extent_int(1);
1924 numFields2_ = inputValues2.extent_int(0);
1925 numPoints2_ = inputValues2.extent_int(1);
1927 numFields3_ = inputValues3.extent_int(0);
1928 numPoints3_ = inputValues3.extent_int(1);
1935 INTREPID2_TEST_FOR_EXCEPTION(inputValues1.rank() > 3, std::invalid_argument,
"ranks greater than 3 not yet supported");
1936 INTREPID2_TEST_FOR_EXCEPTION(inputValues2.rank() > 3, std::invalid_argument,
"ranks greater than 3 not yet supported");
1937 INTREPID2_TEST_FOR_EXCEPTION(inputValues3.rank() > 3, std::invalid_argument,
"ranks greater than 3 not yet supported");
1938 INTREPID2_TEST_FOR_EXCEPTION((inputValues1.rank() == 3) && (inputValues2.rank() == 3), std::invalid_argument,
"two vector-valued input ranks not yet supported");
1939 INTREPID2_TEST_FOR_EXCEPTION((inputValues1.rank() == 3) && (inputValues3.rank() == 3), std::invalid_argument,
"two vector-valued input ranks not yet supported");
1940 INTREPID2_TEST_FOR_EXCEPTION((inputValues2.rank() == 3) && (inputValues3.rank() == 3), std::invalid_argument,
"two vector-valued input ranks not yet supported");
1945 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != numPoints1_, std::invalid_argument,
"incompatible point counts");
1946 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != numPoints2_, std::invalid_argument,
"incompatible point counts");
1947 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != numPoints3_, std::invalid_argument,
"incompatible point counts");
1951 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != numPoints1_ * numPoints2_ * numPoints3_, std::invalid_argument,
"incompatible point counts");
1954 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != numFields1_ * numFields2_ * numFields3_, std::invalid_argument,
"incompatible field sizes");
1957 KOKKOS_INLINE_FUNCTION
1958 void operator()(
const TeamMember & teamMember )
const
1960 auto fieldOrdinal1 = teamMember.league_rank();
1964 if ((input1_.rank() == 2) && (input2_.rank() == 2) && (input3_.rank() == 2))
1966 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (
const int& fieldOrdinal2) {
1967 for (
int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
1969 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
1970 for (
int pointOrdinal=0; pointOrdinal<numPoints_; pointOrdinal++)
1972 output_(fieldOrdinal,pointOrdinal) = weight_ * input1_(fieldOrdinal1,pointOrdinal) * input2_(fieldOrdinal2,pointOrdinal) * input3_(fieldOrdinal3,pointOrdinal);
1977 else if (input1_.rank() == 3)
1979 int spaceDim = input1_.extent_int(2);
1980 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (
const int& fieldOrdinal2) {
1981 for (
int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
1983 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
1984 for (
int pointOrdinal=0; pointOrdinal<numPoints_; pointOrdinal++)
1986 for (
int d=0; d<spaceDim; d++)
1988 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal,d) * input2_(fieldOrdinal2,pointOrdinal) * input3_(fieldOrdinal3,pointOrdinal);
1994 else if (input2_.rank() == 3)
1996 int spaceDim = input2_.extent_int(2);
1997 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (
const int& fieldOrdinal2) {
1998 for (
int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
2000 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
2001 for (
int pointOrdinal=0; pointOrdinal<numPoints_; pointOrdinal++)
2003 for (
int d=0; d<spaceDim; d++)
2005 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal) * input2_(fieldOrdinal2,pointOrdinal,d) * input3_(fieldOrdinal3,pointOrdinal);
2011 else if (input3_.rank() == 3)
2013 int spaceDim = input3_.extent_int(2);
2014 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (
const int& fieldOrdinal2) {
2015 for (
int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
2017 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
2018 for (
int pointOrdinal=0; pointOrdinal<numPoints_; pointOrdinal++)
2020 for (
int d=0; d<spaceDim; d++)
2022 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal) * input2_(fieldOrdinal2,pointOrdinal) * input3_(fieldOrdinal3,pointOrdinal,d);
2035 if ((input1_.rank() == 2) && (input2_.rank() == 2) && (input3_.rank() == 2) )
2037 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (
const int& fieldOrdinal2) {
2038 for (
int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
2040 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
2041 for (
int pointOrdinal3=0; pointOrdinal3<numPoints3_; pointOrdinal3++)
2043 for (
int pointOrdinal2=0; pointOrdinal2<numPoints2_; pointOrdinal2++)
2045 for (
int pointOrdinal1=0; pointOrdinal1<numPoints1_; pointOrdinal1++)
2047 int pointOrdinal = (pointOrdinal3 * numPoints2_ + pointOrdinal2) * numPoints1_ + pointOrdinal1;
2048 output_(fieldOrdinal,pointOrdinal) = weight_ * input1_(fieldOrdinal1,pointOrdinal1) * input2_(fieldOrdinal2,pointOrdinal2) * input3_(fieldOrdinal3,pointOrdinal3);
2055 else if (input1_.rank() == 3)
2057 int spaceDim = input1_.extent_int(2);
2058 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (
const int& fieldOrdinal2) {
2059 for (
int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
2061 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
2062 for (
int pointOrdinal3=0; pointOrdinal3<numPoints3_; pointOrdinal3++)
2064 for (
int pointOrdinal2=0; pointOrdinal2<numPoints2_; pointOrdinal2++)
2066 for (
int pointOrdinal1=0; pointOrdinal1<numPoints1_; pointOrdinal1++)
2068 int pointOrdinal = (pointOrdinal3 * numPoints2_ + pointOrdinal2) * numPoints1_ + pointOrdinal1;
2069 for (
int d=0; d<spaceDim; d++)
2071 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal1,d) * input2_(fieldOrdinal2,pointOrdinal2) * input3_(fieldOrdinal3,pointOrdinal3);
2079 else if (input2_.rank() == 3)
2081 int spaceDim = input2_.extent_int(2);
2082 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (
const int& fieldOrdinal2) {
2083 for (
int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
2085 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
2086 for (
int pointOrdinal3=0; pointOrdinal3<numPoints3_; pointOrdinal3++)
2088 for (
int pointOrdinal2=0; pointOrdinal2<numPoints2_; pointOrdinal2++)
2090 for (
int pointOrdinal1=0; pointOrdinal1<numPoints1_; pointOrdinal1++)
2092 int pointOrdinal = (pointOrdinal3 * numPoints2_ + pointOrdinal2) * numPoints1_ + pointOrdinal1;
2093 for (
int d=0; d<spaceDim; d++)
2095 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal1) * input2_(fieldOrdinal2,pointOrdinal2,d) * input3_(fieldOrdinal3,pointOrdinal3);
2103 else if (input3_.rank() == 3)
2105 int spaceDim = input3_.extent_int(2);
2106 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (
const int& fieldOrdinal2) {
2107 for (
int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
2109 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
2110 for (
int pointOrdinal3=0; pointOrdinal3<numPoints3_; pointOrdinal3++)
2112 for (
int pointOrdinal2=0; pointOrdinal2<numPoints2_; pointOrdinal2++)
2114 for (
int pointOrdinal1=0; pointOrdinal1<numPoints1_; pointOrdinal1++)
2116 int pointOrdinal = (pointOrdinal3 * numPoints2_ + pointOrdinal2) * numPoints1_ + pointOrdinal1;
2117 for (
int d=0; d<spaceDim; d++)
2119 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal1) * input2_(fieldOrdinal2,pointOrdinal2) * input3_(fieldOrdinal3,pointOrdinal3,d);
2136 template<
typename BasisBaseClass =
void>
2140 using BasisBase = BasisBaseClass;
2143 using typename BasisBase::OutputViewType;
2144 using typename BasisBase::PointViewType;
2145 using typename BasisBase::ScalarViewType;
2147 using typename BasisBase::OutputValueType;
2148 using typename BasisBase::PointValueType;
2150 using typename BasisBase::ExecutionSpace;
2152 using BasisPtr = Teuchos::RCP<BasisBase>;
2158 Basis_TensorBasis3(BasisPtr basis1, BasisPtr basis2, BasisPtr basis3,
const bool useShardsCellTopologyAndTags =
false)
2160 TensorBasis(Teuchos::rcp(
new TensorBasis(basis1,basis2,FUNCTION_SPACE_MAX,useShardsCellTopologyAndTags)),
2162 FUNCTION_SPACE_MAX,useShardsCellTopologyAndTags),
2177 std::vector<BasisPtr> componentBases {basis1_, basis2_, basis3_};
2207 virtual void getValues(OutputViewType outputValues,
const EOperator operatorType,
2208 const PointViewType inputPoints12,
const PointViewType inputPoints3,
2209 bool tensorPoints)
const override
2213 int spaceDim1 = basis1_->getDomainDimension();
2214 int spaceDim2 = basis2_->getDomainDimension();
2216 int totalSpaceDim12 = inputPoints12.extent_int(1);
2218 TEUCHOS_ASSERT(spaceDim1 + spaceDim2 == totalSpaceDim12);
2222 auto inputPoints1 = Kokkos::subview(inputPoints12,Kokkos::ALL(),std::make_pair(0,spaceDim1));
2223 auto inputPoints2 = Kokkos::subview(inputPoints12,Kokkos::ALL(),std::make_pair(spaceDim1,totalSpaceDim12));
2225 this->
getValues(outputValues, operatorType, inputPoints1, inputPoints2, inputPoints3, tensorPoints);
2232 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"This method does not yet handle tensorPoints=true");
2263 virtual void getValues(OutputViewType outputValues,
const EOperator operatorType,
2264 const PointViewType inputPoints1,
const PointViewType inputPoints2,
const PointViewType inputPoints3,
2265 bool tensorPoints)
const = 0;
2295 const PointViewType inputPoints1,
const EOperator operatorType1,
2296 const PointViewType inputPoints2,
const EOperator operatorType2,
2297 const PointViewType inputPoints3,
const EOperator operatorType3,
2298 bool tensorPoints,
double weight=1.0)
const
2300 int basisCardinality1 = basis1_->getCardinality();
2301 int basisCardinality2 = basis2_->getCardinality();
2302 int basisCardinality3 = basis3_->getCardinality();
2304 int spaceDim1 = inputPoints1.extent_int(1);
2305 int spaceDim2 = inputPoints2.extent_int(1);
2306 int spaceDim3 = inputPoints3.extent_int(1);
2308 int totalPointCount;
2309 int pointCount1, pointCount2, pointCount3;
2312 pointCount1 = inputPoints1.extent_int(0);
2313 pointCount2 = inputPoints2.extent_int(0);
2314 pointCount3 = inputPoints3.extent_int(0);
2315 totalPointCount = pointCount1 * pointCount2 * pointCount3;
2319 totalPointCount = inputPoints1.extent_int(0);
2320 pointCount1 = totalPointCount;
2321 pointCount2 = totalPointCount;
2322 pointCount3 = totalPointCount;
2324 INTREPID2_TEST_FOR_EXCEPTION((totalPointCount != inputPoints2.extent_int(0)) || (totalPointCount != inputPoints3.extent_int(0)),
2325 std::invalid_argument,
"If tensorPoints is false, the point counts must match!");
2344 OutputViewType outputValues1, outputValues2, outputValues3;
2348 if ((spaceDim1 == 1) && (operatorType1 == OPERATOR_VALUE))
2351 outputValues1 = getMatchingViewWithLabel(outputValues,
"output values - basis 1",basisCardinality1,pointCount1);
2355 outputValues1 = getMatchingViewWithLabel(outputValues,
"output values - basis 1",basisCardinality1,pointCount1,spaceDim1);
2357 if ((spaceDim2 == 1) && (operatorType2 == OPERATOR_VALUE))
2360 outputValues2 = getMatchingViewWithLabel(outputValues,
"output values - basis 2",basisCardinality2,pointCount2);
2364 outputValues2 = getMatchingViewWithLabel(outputValues,
"output values - basis 2",basisCardinality2,pointCount2,spaceDim2);
2366 if ((spaceDim3 == 1) && (operatorType3 == OPERATOR_VALUE))
2369 outputValues3 = getMatchingViewWithLabel(outputValues,
"output values - basis 3",basisCardinality3,pointCount3);
2373 outputValues3 = getMatchingViewWithLabel(outputValues,
"output values - basis 3",basisCardinality3,pointCount3,spaceDim3);
2376 basis1_->getValues(outputValues1,inputPoints1,operatorType1);
2377 basis2_->getValues(outputValues2,inputPoints2,operatorType2);
2378 basis3_->getValues(outputValues3,inputPoints3,operatorType3);
2380 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputValueType>();
2381 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointValueType>();
2382 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
2384 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(basisCardinality1,Kokkos::AUTO(),vectorSize);
2387 FunctorType functor(outputValues, outputValues1, outputValues2, outputValues3, tensorPoints, weight);
2388 Kokkos::parallel_for(
"TensorBasis3_Functor", policy , functor);
2398 virtual void getDofCoeffs(
typename BasisBase::ScalarViewType dofCoeffs )
const override
2400 using ValueType =
typename BasisBase::ScalarViewType::value_type;
2401 using ResultLayout =
typename DeduceLayout< typename BasisBase::ScalarViewType >::result_layout;
2402 using ViewType = Kokkos::DynRankView<ValueType, ResultLayout, typename TensorBasis::DeviceType >;
2404 const ordinal_type basisCardinality1 = basis1_->getCardinality();
2405 const ordinal_type basisCardinality2 = basis2_->getCardinality();
2406 const ordinal_type basisCardinality3 = basis3_->getCardinality();
2408 auto dofCoeffs1 = ViewType(
"dofCoeffs1",basisCardinality1);
2409 auto dofCoeffs2 = ViewType(
"dofCoeffs2",basisCardinality2);
2410 auto dofCoeffs3 = ViewType(
"dofCoeffs3",basisCardinality3);
2412 basis1_->getDofCoeffs(dofCoeffs1);
2413 basis2_->getDofCoeffs(dofCoeffs2);
2414 basis3_->getDofCoeffs(dofCoeffs3);
2416 Kokkos::RangePolicy<ExecutionSpace> policy(0, basisCardinality3);
2417 Kokkos::parallel_for(policy, KOKKOS_LAMBDA (
const int fieldOrdinal3)
2419 for (
int fieldOrdinal2=0; fieldOrdinal2<basisCardinality2; fieldOrdinal2++)
2420 for (
int fieldOrdinal1=0; fieldOrdinal1<basisCardinality1; fieldOrdinal1++)
2422 const ordinal_type fieldOrdinal = fieldOrdinal1 + fieldOrdinal2 * basisCardinality1 + fieldOrdinal3 * (basisCardinality1*basisCardinality2);
2423 dofCoeffs(fieldOrdinal) = dofCoeffs1(fieldOrdinal1);
2424 dofCoeffs(fieldOrdinal) *= dofCoeffs2(fieldOrdinal2) * dofCoeffs3(fieldOrdinal3);
2433 virtual HostBasisPtr<OutputValueType, PointValueType>
2435 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"TensorBasis3 subclasses must override getHostBasis");
KOKKOS_INLINE_FUNCTION ScalarType getView1Entry()
Functor for computing values for the TensorBasis class.
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...
KOKKOS_INLINE_FUNCTION int numFamilies() const
For valid vectorData, returns the number of families in vectorData; otherwise, returns number of Tens...
Implements arbitrary-dimensional extrusion of a base shards::CellTopology.
View-like interface to tensor points; point components are stored separately; the appropriate coordin...
#define INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(test, x, msg)
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
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...
For a multi-component tensor basis, specifies the operators to be applied to the components to produc...
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.
Wrapper around a Kokkos::View that allows data that is constant or repeating in various logical dimen...
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...
void getComponentPoints(const PointViewType inputPoints, const bool attemptTensorDecomposition, PointViewType &inputPoints1, PointViewType &inputPoints2, bool &tensorDecompositionSucceeded) const
Method to extract component points from composite points.
Header function for Intrepid2::Util class and other utility functions.
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.
Implementation of an assert that can safely be called from device code.
const VectorDataType & vectorData() const
VectorData accessor.
Basis_TensorBasis(BasisPtr basis1, BasisPtr basis2, EFunctionSpace functionSpace=FUNCTION_SPACE_MAX, const bool useShardsCellTopologyAndTags=false)
Constructor.
The data containers in Intrepid2 that support sum factorization and other reduced-data optimizations ...
bool rotateXYNinetyDegrees() const
If true, this flag indicates that a vector component that spans the first two dimensions should be ro...
virtual OperatorTensorDecomposition getOperatorDecomposition(const EOperator operatorType) const
Returns a full decomposition of the specified operator. (Full meaning that all TensorBasis components...
KOKKOS_INLINE_FUNCTION void set(ScalarType value)
KOKKOS_INLINE_FUNCTION int increment()
KOKKOS_INLINE_FUNCTION void setLocation(const Kokkos::Array< int, 7 > location)
KOKKOS_INLINE_FUNCTION ordinal_type numTensorComponents() const
Returns the number of tensorial components.
KOKKOS_INLINE_FUNCTION ScalarType getView2Entry()
virtual HostBasisPtr< OutputValueType, PointValueType > getHostBasis() const override
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
virtual OperatorTensorDecomposition getSimpleOperatorDecomposition(const EOperator &operatorType) const
Returns a simple decomposition of the specified operator: what operator(s) should be applied to basis...
static CellTopoPtr cellTopology(const shards::CellTopology &shardsCellTopo, ordinal_type tensorialDegree=0)
static accessor that returns a CellTopoPtr; these are lazily constructed and cached.
TensorDataType & tensorData()
TensorData accessor for single-family scalar data.
virtual void getDofCoords(typename BasisBase::ScalarViewType dofCoords) const override
Fills in spatial locations (coordinates) of degrees of freedom (nodes) on the reference cell...
virtual const char * getName() const override
Returns basis name.
Implementation of support for traversing component views alongside a view that represents a combinati...
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...
KOKKOS_INLINE_FUNCTION int nextIncrementRank()
virtual void getDofCoeffs(typename BasisBase::ScalarViewType dofCoeffs) const override
Fills in coefficients of degrees of freedom on the reference cell.
ordinal_type componentPointCount(const ordinal_type &tensorComponentOrdinal) const
Returns the number of points in the indicated component.
virtual void getDofCoeffs(typename BasisBase::ScalarViewType dofCoeffs) const override
Fills in coefficients of degrees of freedom on the reference cell.
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...
Class that defines mappings from component cell topologies to their tensor product topologies...
virtual HostBasisPtr< OutputValueType, PointValueType > getHostBasis() const override
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
View-like interface to tensor data; tensor components are stored separately and multiplied together a...
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 OperatorTensorDecomposition getOperatorDecomposition(const EOperator operatorType) const override
Returns a full decomposition of the specified operator. (Full meaning that all TensorBasis components...
A helper class that allows iteration over three Kokkos Views simultaneously, according to tensor comb...
Basis defined as the tensor product of two component bases.
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...
KOKKOS_INLINE_FUNCTION ScalarView< PointScalar, DeviceType > getTensorComponent(const ordinal_type &r) const
Returns the requested tensor component.
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.
Reference-space field values for a basis, designed to support typical vector-valued bases...
Functor for computing values for the TensorBasis3 class.
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.
Header file for the abstract base class Intrepid2::Basis.