Intrepid2
Intrepid2_TensorBasis.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid2 Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Kyungjoo Kim (kyukim@sandia.gov),
38// Mauro Perego (mperego@sandia.gov), or
39// Nate Roberts (nvrober@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
49#ifndef Intrepid2_TensorBasis_h
50#define Intrepid2_TensorBasis_h
51
52#include <Kokkos_DynRankView.hpp>
53
54#include <Teuchos_RCP.hpp>
55
56#include <Intrepid2_config.h>
57
58#include <map>
59#include <set>
60#include <vector>
61
62#include "Intrepid2_Basis.hpp"
66#include "Intrepid2_Utils.hpp" // defines FAD_VECTOR_SIZE, VECTOR_SIZE
67
69
70namespace Intrepid2
71{
72 template<ordinal_type spaceDim>
73 KOKKOS_INLINE_FUNCTION
74 ordinal_type getDkEnumeration(Kokkos::Array<int,spaceDim> &entries);
75
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);
79
80 template<>
81 KOKKOS_INLINE_FUNCTION
82 void getDkEnumerationInverse<1>(Kokkos::Array<int,1> &entries, const ordinal_type dkEnum, const ordinal_type operatorOrder)
83 {
84 entries[0] = operatorOrder;
85 }
86
87 template<>
88 KOKKOS_INLINE_FUNCTION
89 void getDkEnumerationInverse<2>(Kokkos::Array<int,2> &entries, const ordinal_type dkEnum, const ordinal_type operatorOrder)
90 {
91 entries[0] = operatorOrder - dkEnum;
92 entries[1] = dkEnum;
93 }
94
95 template<>
96 KOKKOS_INLINE_FUNCTION
97 void getDkEnumerationInverse<3>(Kokkos::Array<int,3> &entries, const ordinal_type dkEnum, const ordinal_type operatorOrder)
98 {
99 // formula is zMult + (yMult+zMult)*(yMult+zMult+1)/2; where xMult+yMult+zMult = operatorOrder
100 // it seems complicated to work out a formula that will invert this. For the present we just take a brute force approach,
101 // using getDkEnumeration() to check each possibility
102 for (ordinal_type yMult=0; yMult<=operatorOrder; yMult++)
103 {
104 for (ordinal_type zMult=0; zMult<=operatorOrder-yMult; zMult++)
105 {
106 const ordinal_type xMult = operatorOrder-(zMult+yMult);
107 if (dkEnum == getDkEnumeration<3>(xMult,yMult,zMult))
108 {
109 entries[0] = xMult;
110 entries[1] = yMult;
111 entries[2] = zMult;
112 return;
113 }
114 }
115 }
116 }
117
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)
121 {
122 // for operator order k, the recursive formula defining getDkEnumeration is:
123 // getDkEnumeration(k0,k1,…,k_{n-1}) = getDkCardinality(k - k0) + getDkEnumeration(k1,…,k_{n-1})
124 // The entries are in reverse lexicographic order. We search for k0, by incrementing k0 until getDkEnumeration(k0,0,…,0) <= dkEnum
125 // Then we recursively call getDkEnumerationInverse<spaceDim-1>({k1,…,k_{n-1}}, dkEnum - getDkEnumeration(k0,0,…,0) - 1)
126
127 for (int k0=0; k0<=operatorOrder; k0++)
128 {
129 entries[0] = k0;
130 for (int d=1; d<spaceDim-1; d++)
131 {
132 entries[d] = 0;
133 }
134 // sum of entries must be equal to operatorOrder
135 if (spaceDim > 1) entries[spaceDim-1] = operatorOrder - k0;
136 else if (k0 != operatorOrder) continue; // if spaceDim == 1, then the only way the sum of the entries is operatorOrder is if k0 == operatorOrder
137 const ordinal_type dkEnumFor_k0 = getDkEnumeration<spaceDim>(entries);
138
139 if (dkEnumFor_k0 > dkEnum) continue; // next k0
140 else if (dkEnumFor_k0 == dkEnum) return; // entries has (k0,0,…,0), and this has dkEnum as its enumeration value
141 else
142 {
143 // (k0,0,…,0) is prior to the dkEnum entry, which means that the dkEnum entry starts with k0-1.
144 entries[0] = k0 - 1;
145
146 // We determine the rest of the entries through a recursive call to getDkEnumerationInverse<spaceDim - 1>().
147
148 // ensure that we don't try to allocate an empty array…
149 constexpr ordinal_type sizeForSubArray = (spaceDim > 2) ? spaceDim - 1 : 1;
150 Kokkos::Array<int,sizeForSubArray> subEntries;
151
152 // the -1 in sub-entry enumeration value accounts for the fact that the entry is the one *after* (k0,0,…,0)
153 getDkEnumerationInverse<spaceDim-1>(subEntries, dkEnum - dkEnumFor_k0 - 1, operatorOrder - entries[0]);
154
155 for (int i=1; i<spaceDim; i++)
156 {
157 entries[i] = subEntries[i-1];
158 }
159 return;
160 }
161 }
162 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "entries corresponding to dkEnum not found");
163 }
164
165 template<>
166 KOKKOS_INLINE_FUNCTION
167 ordinal_type getDkEnumeration<1>(Kokkos::Array<int,1> &entries)
168 {
169 return getDkEnumeration<1>(entries[0]);
170 }
171
172 template<ordinal_type spaceDim>
173 KOKKOS_INLINE_FUNCTION
174 ordinal_type getDkEnumeration(Kokkos::Array<int,spaceDim> &entries)
175 {
176 ordinal_type k_minus_k0 = 0; // sum of all the entries but the first
177
178 // recursive formula in general is: getDkEnumeration(k0,k1,…,k_{n-1}) = getDkCardinality(k - k0) + getDkEnumeration(k1,…,k_{n-1})
179 // ensure that we don't try to allocate an empty array…
180 constexpr ordinal_type sizeForSubArray = (spaceDim > 2) ? spaceDim - 1 : 1;
181 Kokkos::Array<int,sizeForSubArray> remainingEntries;
182 for (int i=1; i<spaceDim; i++)
183 {
184 k_minus_k0 += entries[i];
185 remainingEntries[i-1] = entries[i];
186 }
187
188 if (k_minus_k0 == 0)
189 {
190 return 0;
191 }
192 else
193 {
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);
197 return dkEnum;
198 }
199 }
200
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)
205 {
206 Kokkos::Array<int,spaceDim1> entries1;
207 getDkEnumerationInverse<spaceDim1>(entries1, dkEnum1, operatorOrder1);
208
209 Kokkos::Array<int,spaceDim2> entries2;
210 getDkEnumerationInverse<spaceDim2>(entries2, dkEnum2, operatorOrder2);
211
212 const int spaceDim = spaceDim1 + spaceDim2;
213 Kokkos::Array<int,spaceDim> entries;
214
215 for (ordinal_type d=0; d<spaceDim1; d++)
216 {
217 entries[d] = entries1[d];
218 }
219
220 for (ordinal_type d=0; d<spaceDim2; d++)
221 {
222 entries[d+spaceDim1] = entries2[d];
223 }
224
225 return getDkEnumeration<spaceDim>(entries);
226 }
227
228template<typename BasisBase>
229class Basis_TensorBasis;
230
235{
236 // if we want to make this usable on device, we could switch to Kokkos::Array instead of std::vector. But this is not our immediate use case.
237 std::vector< std::vector<EOperator> > ops; // outer index: vector entry ordinal; inner index: basis component ordinal. (scalar-valued operators have a single entry in outer vector)
238 std::vector<double> weights; // weights for each vector entry
239 ordinal_type numBasisComponents_;
240 bool rotateXYNinetyDegrees_ = false; // if true, indicates that something that otherwise would have values (f_x, f_y, …) should be mapped to (-f_y, f_x, …). This is used for H(curl) wedges (specifically, for OPERATOR_CURL).
241
242 OperatorTensorDecomposition(const std::vector<EOperator> &opsBasis1, const std::vector<EOperator> &opsBasis2, const std::vector<double> vectorComponentWeights)
243 :
244 weights(vectorComponentWeights),
245 numBasisComponents_(2)
246 {
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()");
252
253 for (ordinal_type i=0; i<size; i++)
254 {
255 ops.push_back(std::vector<EOperator>{opsBasis1[i],opsBasis2[i]});
256 }
257 }
258
259 OperatorTensorDecomposition(const std::vector< std::vector<EOperator> > &vectorEntryOps, const std::vector<double> &vectorComponentWeights)
260 :
261 ops(vectorEntryOps),
262 weights(vectorComponentWeights)
263 {
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()");
267
268 INTREPID2_TEST_FOR_EXCEPTION(numVectorComponents == 0, std::invalid_argument, "must have at least one entry!");
269
270 ordinal_type numBases = 0;
271 for (ordinal_type i=0; i<numVectorComponents; i++)
272 {
273 if (numBases == 0)
274 {
275 numBases = ops[i].size();
276 }
277 else if (ops[i].size() != 0)
278 {
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");
281 }
282 }
283 INTREPID2_TEST_FOR_EXCEPTION(numBases == 0, std::invalid_argument, "at least one vectorEntryOps entry must be non-trivial");
284 numBasisComponents_ = numBases;
285 }
286
287 OperatorTensorDecomposition(const std::vector<EOperator> &basisOps, const double weight = 1.0)
288 :
289 ops({basisOps}),
290 weights({weight}),
291 numBasisComponents_(basisOps.size())
292 {}
293
294 OperatorTensorDecomposition(const EOperator &opBasis1, const EOperator &opBasis2, double weight = 1.0)
295 :
296 ops({ std::vector<EOperator>{opBasis1, opBasis2} }),
297 weights({weight}),
298 numBasisComponents_(2)
299 {}
300
301 OperatorTensorDecomposition(const EOperator &opBasis1, const EOperator &opBasis2, const EOperator &opBasis3, double weight = 1.0)
302 :
303 ops({ std::vector<EOperator>{opBasis1, opBasis2, opBasis3} }),
304 weights({weight}),
305 numBasisComponents_(3)
306 {}
307
308 ordinal_type numVectorComponents() const
309 {
310 return ops.size(); // will match weights.size()
311 }
312
313 ordinal_type numBasisComponents() const
314 {
315 return numBasisComponents_;
316 }
317
318 double weight(const ordinal_type &vectorComponentOrdinal) const
319 {
320 return weights[vectorComponentOrdinal];
321 }
322
323 bool identicallyZeroComponent(const ordinal_type &vectorComponentOrdinal) const
324 {
325 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(vectorComponentOrdinal < 0, std::invalid_argument, "vectorComponentOrdinal is out of bounds");
326 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(vectorComponentOrdinal >= numVectorComponents(), std::invalid_argument, "vectorComponentOrdinal is out of bounds");
327 return ops[vectorComponentOrdinal].size() == 0;
328 }
329
330 EOperator op(const ordinal_type &vectorComponentOrdinal, const ordinal_type &basisOrdinal) const
331 {
332 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(vectorComponentOrdinal < 0, std::invalid_argument, "vectorComponentOrdinal is out of bounds");
333 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(vectorComponentOrdinal >= numVectorComponents(), std::invalid_argument, "vectorComponentOrdinal is out of bounds");
334 if (identicallyZeroComponent(vectorComponentOrdinal))
335 {
336 return OPERATOR_MAX; // by convention: zero in this component
337 }
338 else
339 {
340 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(basisOrdinal < 0, std::invalid_argument, "basisOrdinal is out of bounds");
341 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(basisOrdinal >= numBasisComponents_, std::invalid_argument, "basisOrdinal is out of bounds");
342 return ops[vectorComponentOrdinal][basisOrdinal];
343 }
344 }
345
347 template<typename DeviceType, typename OutputValueType, class PointValueType>
349 {
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");
352
353 ordinal_type numExpandedBasisComponents = 0;
355 using TensorBasis = Basis_TensorBasis<BasisBase>;
356 std::vector<TensorBasis*> basesAsTensorBasis(numBasisComponents_);
357 for (ordinal_type basisComponentOrdinal=0; basisComponentOrdinal<numBasisComponents_; basisComponentOrdinal++)
358 {
359 TensorBasis* basisAsTensorBasis = dynamic_cast<TensorBasis*>(bases[basisComponentOrdinal].get());
360 basesAsTensorBasis[basisComponentOrdinal] = basisAsTensorBasis;
361 if (basisAsTensorBasis)
362 {
363 numExpandedBasisComponents += basisAsTensorBasis->getTensorBasisComponents().size();
364 }
365 else
366 {
367 numExpandedBasisComponents += 1;
368 }
369 }
370
371 std::vector< std::vector<EOperator> > expandedOps; // outer index: vector entry ordinal; inner index: basis component ordinal.
372 std::vector<double> expandedWeights;
373 const ordinal_type opsSize = ops.size();
374 for (ordinal_type simpleVectorEntryOrdinal=0; simpleVectorEntryOrdinal<opsSize; simpleVectorEntryOrdinal++)
375 {
376 if (identicallyZeroComponent(simpleVectorEntryOrdinal))
377 {
378 expandedOps.push_back(std::vector<EOperator>{});
379 expandedWeights.push_back(0.0);
380 continue;
381 }
382
383 std::vector< std::vector<EOperator> > expandedBasisOpsForSimpleVectorEntry(1); // start out with one outer entry; expands if a component is vector-valued
384
385 // this lambda appends an op to each of the vector components
386 auto addExpandedOp = [&expandedBasisOpsForSimpleVectorEntry](const EOperator &op)
387 {
388 const ordinal_type size = expandedBasisOpsForSimpleVectorEntry.size();
389 for (ordinal_type i=0; i<size; i++)
390 {
391 expandedBasisOpsForSimpleVectorEntry[i].push_back(op);
392 }
393 };
394
395 // this lambda takes a scalar-valued (single outer entry) expandedBasisOps and expands it
396 // according to the number of vector entries coming from the vector-valued component basis
397 auto vectorizeExpandedOps = [&expandedBasisOpsForSimpleVectorEntry](const int &numSubVectors)
398 {
399 // we require that this only gets called once per simpleVectorEntryOrdinal -- i.e., only one basis component gets to be vector-valued.
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++)
402 {
403 expandedBasisOpsForSimpleVectorEntry.push_back(expandedBasisOpsForSimpleVectorEntry[0]);
404 }
405 };
406
407 std::vector<EOperator> subVectorOps; // only used if one of the components is vector-valued
408 std::vector<double> subVectorWeights {weights[simpleVectorEntryOrdinal]};
409 for (ordinal_type basisComponentOrdinal=0; basisComponentOrdinal<numBasisComponents_; basisComponentOrdinal++)
410 {
411 const auto &op = ops[simpleVectorEntryOrdinal][basisComponentOrdinal];
412
413 if (! basesAsTensorBasis[basisComponentOrdinal])
414 {
415 addExpandedOp(op);
416 }
417 else
418 {
419 OperatorTensorDecomposition basisOpDecomposition = basesAsTensorBasis[basisComponentOrdinal]->getOperatorDecomposition(op);
420 if (basisOpDecomposition.numVectorComponents() > 1)
421 {
422 // We don't currently support a use case where we have multiple component bases that are vector-valued:
423 INTREPID2_TEST_FOR_EXCEPTION(subVectorWeights.size() > 1, std::invalid_argument, "Unhandled case: multiple component bases are vector-valued");
424 // We do support a single vector-valued case, though; this splits the current simpleVectorEntryOrdinal into an appropriate number of components that come in order in the expanded vector
425 ordinal_type numSubVectors = basisOpDecomposition.numVectorComponents();
426 vectorizeExpandedOps(numSubVectors);
427
428 double weightSoFar = subVectorWeights[0];
429 for (ordinal_type subVectorEntryOrdinal=1; subVectorEntryOrdinal<numSubVectors; subVectorEntryOrdinal++)
430 {
431 subVectorWeights.push_back(weightSoFar * basisOpDecomposition.weight(subVectorEntryOrdinal));
432 }
433 subVectorWeights[0] *= basisOpDecomposition.weight(0);
434 for (ordinal_type subVectorEntryOrdinal=0; subVectorEntryOrdinal<numSubVectors; subVectorEntryOrdinal++)
435 {
436 for (ordinal_type subComponentBasis=0; subComponentBasis<basisOpDecomposition.numBasisComponents(); subComponentBasis++)
437 {
438 const auto &basisOp = basisOpDecomposition.op(subVectorEntryOrdinal, subComponentBasis);
439 expandedBasisOpsForSimpleVectorEntry[subVectorEntryOrdinal].push_back(basisOp);
440 }
441 }
442 }
443 else
444 {
445 double componentWeight = basisOpDecomposition.weight(0);
446 const ordinal_type size = subVectorWeights.size();
447 for (ordinal_type i=0; i<size; i++)
448 {
449 subVectorWeights[i] *= componentWeight;
450 }
451 ordinal_type subVectorEntryOrdinal = 0;
452 const ordinal_type numBasisComponents = basisOpDecomposition.numBasisComponents();
453 for (ordinal_type subComponentBasis=0; subComponentBasis<numBasisComponents; subComponentBasis++)
454 {
455 const auto &basisOp = basisOpDecomposition.op(subVectorEntryOrdinal, basisComponentOrdinal);
456 addExpandedOp( basisOp );
457 }
458 }
459 }
460 }
461
462 // sanity check on the new expandedOps entries:
463 for (ordinal_type i=0; i<static_cast<ordinal_type>(expandedBasisOpsForSimpleVectorEntry.size()); i++)
464 {
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");
467 }
468
469 expandedOps.insert(expandedOps.end(), expandedBasisOpsForSimpleVectorEntry.begin(), expandedBasisOpsForSimpleVectorEntry.end());
470 expandedWeights.insert(expandedWeights.end(), subVectorWeights.begin(), subVectorWeights.end());
471 }
472 // check that vector lengths agree:
473 INTREPID2_TEST_FOR_EXCEPTION(expandedOps.size() != expandedWeights.size(), std::logic_error, "expandedWeights and expandedOps do not agree on the number of vector components");
474
475 OperatorTensorDecomposition result(expandedOps, expandedWeights);
476 result.setRotateXYNinetyDegrees(rotateXYNinetyDegrees_);
477 return result;
478 }
479
482 {
483 return rotateXYNinetyDegrees_;
484 }
485
486 void setRotateXYNinetyDegrees(const bool &value)
487 {
488 rotateXYNinetyDegrees_ = value;
489 }
490};
491
497 template<class ExecutionSpace, class OutputScalar, class OutputFieldType>
499 {
500 using ScratchSpace = typename ExecutionSpace::scratch_memory_space;
501 using OutputScratchView = Kokkos::View<OutputScalar*,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
502
503 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
504 using TeamMember = typename TeamPolicy::member_type;
505
507 using RankCombinationType = typename TensorViewIteratorType::RankCombinationType;
508
509 OutputFieldType output_; // F,P[,D…]
510 OutputFieldType input1_; // F1,P[,D…] or F1,P1[,D…]
511 OutputFieldType input2_; // F2,P[,D…] or F2,P2[,D…]
512
513 int numFields_, numPoints_;
514 int numFields1_, numPoints1_;
515 int numFields2_, numPoints2_;
516
517 bool tensorPoints_; // if true, input1 and input2 refer to values at decomposed points, and P = P1 * P2. If false, then the two inputs refer to points in the full-dimensional space, and their point lengths are the same as that of the final output.
518
519 using RankCombinationViewType = typename TensorViewIteratorType::RankCombinationViewType;
520 RankCombinationViewType rank_combinations_;// indicates the policy by which the input views will be combined in output view
521
522 double weight_;
523
524 public:
525
526 TensorViewFunctor(OutputFieldType output, OutputFieldType inputValues1, OutputFieldType inputValues2,
527 bool tensorPoints, double weight)
528 : output_(output), input1_(inputValues1), input2_(inputValues2), tensorPoints_(tensorPoints), weight_(weight)
529 {
530 numFields_ = output.extent_int(0);
531 numPoints_ = output.extent_int(1);
532
533 numFields1_ = inputValues1.extent_int(0);
534 numPoints1_ = inputValues1.extent_int(1);
535
536 numFields2_ = inputValues2.extent_int(0);
537 numPoints2_ = inputValues2.extent_int(1);
538
539 if (!tensorPoints_)
540 {
541 // then the point counts should all match
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");
544 }
545 else
546 {
547 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != numPoints1_ * numPoints2_, std::invalid_argument, "incompatible point counts");
548 }
549
550 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != numFields1_ * numFields2_, std::invalid_argument, "incompatible field sizes");
551
552 const ordinal_type max_rank = std::max(inputValues1.rank(),inputValues2.rank());
553 // at present, no supported case will result in an output rank greater than both input ranks
554
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_);
559
560 rank_combinations_host[0] = TensorViewIteratorType::TENSOR_PRODUCT; // field combination is always tensor product
561 rank_combinations_host[1] = tensorPoints ? TensorViewIteratorType::TENSOR_PRODUCT : TensorViewIteratorType::DIMENSION_MATCH; // tensorPoints controls interpretation of the point dimension
562 for (ordinal_type d=2; d<max_rank; d++)
563 {
564 // d >= 2 have the interpretation of spatial dimensions (gradients, etc.)
565 // we let the extents of the containers determine what we're doing here
566 if ((inputValues1.extent_int(d) == inputValues2.extent_int(d)) && (output.extent_int(d) == 1))
567 {
568 rank_combinations_host[d] = TensorViewIteratorType::TENSOR_CONTRACTION;
569 }
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))
572 )
573 {
574 // this looks like multiplication of a vector by a scalar, resulting in a vector
575 // this can be understood as a tensor product
576 rank_combinations_host[d] = TensorViewIteratorType::TENSOR_PRODUCT;
577 }
578 else if ((inputValues1.extent_int(d) == inputValues2.extent_int(d)) && (output.extent_int(d) == inputValues1.extent_int(d) * inputValues2.extent_int(d)))
579 {
580 // this is actually a generalization of the above case: a tensor product, something like a vector outer product
581 rank_combinations_host[d] = TensorViewIteratorType::TENSOR_PRODUCT;
582 }
583 else if ((inputValues1.extent_int(d) == inputValues2.extent_int(d)) && (output.extent_int(d) == inputValues1.extent_int(d)))
584 {
585 // it's a bit weird (I'm not aware of the use case, in the present context), but we can handle this case by adopting DIMENSION_MATCH here
586 // this is something like MATLAB's .* and .+ operators, which operate entry-wise
587 rank_combinations_host[d] = TensorViewIteratorType::DIMENSION_MATCH;
588 }
589 else
590 {
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");
595 }
596 }
597 Kokkos::deep_copy(rank_combinations_,rank_combinations_host);
598 }
599
600 KOKKOS_INLINE_FUNCTION
601 void operator()( const TeamMember & teamMember ) const
602 {
603 auto fieldOrdinal1 = teamMember.league_rank();
604
605 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (const int& fieldOrdinal2) {
606 TensorViewIteratorType it(output_,input1_,input2_,rank_combinations_);
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; // used to initialize prev_increment_rank at the start of the do/while loop. Notionally, we last incremented in the field ordinal rank to get to the {fieldOrdinal1,0,0,0,0,0,0},{fieldOrdinal2,0,0,0,0,0,0} location.
610 OutputScalar accumulator = 0;
611
612 do
613 {
614 accumulator += weight_ * it.getView1Entry() * it.getView2Entry();
615 next_increment_rank = it.nextIncrementRank();
616
617 if ((next_increment_rank < 0) || (rank_combinations_[next_increment_rank] != TensorViewIteratorType::TENSOR_CONTRACTION))
618 {
619 // then we've finished the accumulation and should set the value
620 it.set(accumulator);
621 // reset the accumulator:
622 accumulator = 0;
623 }
624 } while (it.increment() > FIELD_ORDINAL_DIMENSION);
625 });
626 }
627 };
628
642 template<typename BasisBaseClass = void>
644 :
645 public BasisBaseClass
646 {
647 public:
648 using BasisBase = BasisBaseClass;
649 using BasisPtr = Teuchos::RCP<BasisBase>;
650
651 protected:
652 BasisPtr basis1_;
653 BasisPtr basis2_;
654
655 std::vector<BasisPtr> tensorComponents_;
656
657 std::string name_; // name of the basis
658
659 int numTensorialExtrusions_; // relative to cell topo returned by getBaseCellTopology().
660 public:
661 using DeviceType = typename BasisBase::DeviceType;
662 using ExecutionSpace = typename BasisBase::ExecutionSpace;
663 using OutputValueType = typename BasisBase::OutputValueType;
664 using PointValueType = typename BasisBase::PointValueType;
665
666 using OrdinalTypeArray1DHost = typename BasisBase::OrdinalTypeArray1DHost;
667 using OrdinalTypeArray2DHost = typename BasisBase::OrdinalTypeArray2DHost;
668 using OutputViewType = typename BasisBase::OutputViewType;
669 using PointViewType = typename BasisBase::PointViewType;
671 public:
678 Basis_TensorBasis(BasisPtr basis1, BasisPtr basis2, EFunctionSpace functionSpace = FUNCTION_SPACE_MAX,
679 const bool useShardsCellTopologyAndTags = false)
680 :
681 basis1_(basis1),basis2_(basis2)
682 {
683 this->functionSpace_ = functionSpace;
684
685 Basis_TensorBasis* basis1AsTensor = dynamic_cast<Basis_TensorBasis*>(basis1_.get());
686 if (basis1AsTensor)
687 {
688 auto basis1Components = basis1AsTensor->getTensorBasisComponents();
689 tensorComponents_.insert(tensorComponents_.end(), basis1Components.begin(), basis1Components.end());
690 }
691 else
692 {
693 tensorComponents_.push_back(basis1_);
694 }
695
696 Basis_TensorBasis* basis2AsTensor = dynamic_cast<Basis_TensorBasis*>(basis2_.get());
697 if (basis2AsTensor)
698 {
699 auto basis2Components = basis2AsTensor->getTensorBasisComponents();
700 tensorComponents_.insert(tensorComponents_.end(), basis2Components.begin(), basis2Components.end());
701 }
702 else
703 {
704 tensorComponents_.push_back(basis2_);
705 }
706
707 this->basisCardinality_ = basis1->getCardinality() * basis2->getCardinality();
708 this->basisDegree_ = std::max(basis1->getDegree(), basis2->getDegree());
709
710 {
711 std::ostringstream basisName;
712 basisName << basis1->getName() << " x " << basis2->getName();
713 name_ = basisName.str();
714 }
715
716 // set cell topology
717 this->basisCellTopology_ = tensorComponents_[0]->getBaseCellTopology();
718 this->numTensorialExtrusions_ = tensorComponents_.size() - 1;
719
720 this->basisType_ = basis1_->getBasisType();
721 this->basisCoordinates_ = COORDINATES_CARTESIAN;
722
723 ordinal_type spaceDim1 = basis1_->getDomainDimension();
724 ordinal_type spaceDim2 = basis2_->getDomainDimension();
725
726 INTREPID2_TEST_FOR_EXCEPTION(spaceDim2 != 1, std::invalid_argument, "TensorBasis only supports 1D bases in basis2_ position");
727
728 if (this->getBasisType() == BASIS_FEM_HIERARCHICAL)
729 {
730 // fill in degree lookup:
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);
734
735 const ordinal_type basis1Cardinality = basis1_->getCardinality();
736 const ordinal_type basis2Cardinality = basis2_->getCardinality();
737
738 int degreeLengthField1 = basis1_->getPolynomialDegreeLength();
739 int degreeLengthField2 = basis2_->getPolynomialDegreeLength();
740
741 for (ordinal_type fieldOrdinal1 = 0; fieldOrdinal1 < basis1Cardinality; fieldOrdinal1++)
742 {
743 OrdinalTypeArray1DHost degreesField1 = basis1_->getPolynomialDegreeOfField(fieldOrdinal1);
744 OrdinalTypeArray1DHost h1DegreesField1 = basis1_->getH1PolynomialDegreeOfField(fieldOrdinal1);
745 for (ordinal_type fieldOrdinal2 = 0; fieldOrdinal2 < basis2Cardinality; fieldOrdinal2++)
746 {
747 OrdinalTypeArray1DHost degreesField2 = basis2_->getPolynomialDegreeOfField(fieldOrdinal2);
748 OrdinalTypeArray1DHost h1DegreesField2 = basis2_->getH1PolynomialDegreeOfField(fieldOrdinal2);
749 const ordinal_type tensorFieldOrdinal = fieldOrdinal2 * basis1Cardinality + fieldOrdinal1;
750
751 for (int d3=0; d3<degreeLengthField1; d3++)
752 {
753 this->fieldOrdinalPolynomialDegree_ (tensorFieldOrdinal,d3) = degreesField1(d3);
754 this->fieldOrdinalH1PolynomialDegree_(tensorFieldOrdinal,d3) = h1DegreesField1(d3);
755 }
756 for (int d3=0; d3<degreeLengthField2; d3++)
757 {
758 this->fieldOrdinalPolynomialDegree_ (tensorFieldOrdinal,d3+degreeLengthField1) = degreesField2(d3);
759 this->fieldOrdinalH1PolynomialDegree_(tensorFieldOrdinal,d3+degreeLengthField1) = h1DegreesField2(d3);
760 }
761 }
762 }
763 }
764
765 if (useShardsCellTopologyAndTags)
766 {
767 setShardsTopologyAndTags();
768 }
769 else
770 {
771 // we build tags recursively, making reference to basis1_ and basis2_'s tags to produce the tensor product tags.
772 // // initialize tags
773 const auto & cardinality = this->basisCardinality_;
774
775 // Basis-dependent initializations
776 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
777 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
778 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
779 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
780 const ordinal_type posDfCnt = 3; // position in the tag, counting from 0, of DoF count for the subcell
781
782 OrdinalTypeArray1DHost tagView("tag view", cardinality*tagSize);
783
784 // we assume that basis2_ is defined on a line, and that basis1_ is defined on a domain that is once-extruded in by that line.
785 auto cellTopo = CellTopology::cellTopology(this->basisCellTopology_, numTensorialExtrusions_);
786 auto basis1Topo = cellTopo->getTensorialComponent();
787
788 const ordinal_type spaceDim = spaceDim1 + spaceDim2;
789 const ordinal_type sideDim = spaceDim - 1;
790
791 const OrdinalTypeArray2DHost ordinalToTag1 = basis1_->getAllDofTags();
792 const OrdinalTypeArray2DHost ordinalToTag2 = basis2_->getAllDofTags();
793
794 for (int fieldOrdinal1=0; fieldOrdinal1<basis1_->getCardinality(); fieldOrdinal1++)
795 {
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++)
800 {
801 ordinal_type subcellDim2 = ordinalToTag2(fieldOrdinal2,posScDim);
802 ordinal_type subcellOrd2 = ordinalToTag2(fieldOrdinal2,posScOrd);
803 ordinal_type subcellDfCnt2 = ordinalToTag2(fieldOrdinal2,posDfCnt);
804
805 ordinal_type subcellDim = subcellDim1 + subcellDim2;
806 ordinal_type subcellOrd;
807 if (subcellDim2 == 0)
808 {
809 // vertex node in extrusion; the subcell is not extruded but belongs to one of the two "copies"
810 // of the basis1 topology
811 ordinal_type sideOrdinal = cellTopo->getTensorialComponentSideOrdinal(subcellOrd2); // subcellOrd2 is a "side" of the line topology
812 subcellOrd = CellTopology::getSubcellOrdinalMap(cellTopo, sideDim, sideOrdinal,
813 subcellDim1, subcellOrd1);
814 }
815 else
816 {
817 // line subcell in time; the subcell *is* extruded in final dimension
818 subcellOrd = cellTopo->getExtrudedSubcellOrdinal(subcellDim1, subcellOrd1);
819 if (subcellOrd == -1)
820 {
821 std::cout << "ERROR: -1 subcell ordinal.\n";
822 subcellOrd = cellTopo->getExtrudedSubcellOrdinal(subcellDim1, subcellOrd1);
823 }
824 }
825 ordinal_type tensorFieldOrdinal = fieldOrdinal2 * basis1_->getCardinality() + fieldOrdinal1;
826 // cout << "(" << fieldOrdinal1 << "," << fieldOrdinal2 << ") --> " << i << endl;
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; // subcellDim
832 tagView(tagSize*tensorFieldOrdinal + posScOrd) = subcellOrd; // subcell ordinal
833 tagView(tagSize*tensorFieldOrdinal + posDfOrd) = dofOffsetOrdinal; // ordinal of the specified DoF relative to the subcell
834 tagView(tagSize*tensorFieldOrdinal + posDfCnt) = subcellDfCnt1 * subcellDfCnt2; // total number of DoFs associated with the subcell
835 }
836 }
837
838 // // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
839 // // tags are constructed on host
840 this->setOrdinalTagData(this->tagToOrdinal_,
841 this->ordinalToTag_,
842 tagView,
843 this->basisCardinality_,
844 tagSize,
845 posScDim,
846 posScOrd,
847 posDfOrd);
848 }
849 }
850
851 void setShardsTopologyAndTags()
852 {
853 shards::CellTopology cellTopo1 = basis1_->getBaseCellTopology();
854 shards::CellTopology cellTopo2 = basis2_->getBaseCellTopology();
855
856 auto cellKey1 = basis1_->getBaseCellTopology().getKey();
857 auto cellKey2 = basis2_->getBaseCellTopology().getKey();
858
859 const int numTensorialExtrusions = basis1_->getNumTensorialExtrusions() + basis2_->getNumTensorialExtrusions();
860 if ((cellKey1 == shards::Line<2>::key) && (cellKey2 == shards::Line<2>::key) && (numTensorialExtrusions == 0))
861 {
862 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
863 }
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))
867 )
868 {
869 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
870 }
871 else if ((cellKey1 == shards::Triangle<3>::key) && (cellKey2 == shards::Line<2>::key))
872 {
873 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Wedge<6> >() );
874 }
875 else
876 {
877 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Cell topology combination not yet supported");
878 }
879
880 // numTensorialExtrusions_ is relative to the basisCellTopology_; what we've just done is found a cell topology of the same spatial dimension as the extruded topology, so now numTensorialExtrusions_ should be 0.
881 numTensorialExtrusions_ = 0;
882
883 // initialize tags
884 {
885 const auto & cardinality = this->basisCardinality_;
886
887 // Basis-dependent initializations
888 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
889 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
890 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
891 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
892
893 OrdinalTypeArray1DHost tagView("tag view", cardinality*tagSize);
894
895 shards::CellTopology cellTopo = this->basisCellTopology_;
896
897 ordinal_type tensorSpaceDim = cellTopo.getDimension();
898 ordinal_type spaceDim1 = cellTopo1.getDimension();
899 ordinal_type spaceDim2 = cellTopo2.getDimension();
900
901 TensorTopologyMap topoMap(cellTopo1, cellTopo2);
902
903 for (ordinal_type d=0; d<=tensorSpaceDim; d++) // d: tensorial dimension
904 {
905 ordinal_type d2_max = std::min(spaceDim2,d);
906 int subcellOffset = 0; // for this dimension of tensor subcells, how many subcells have we already counted with other d2/d1 combos?
907 for (ordinal_type d2=0; d2<=d2_max; d2++)
908 {
909 ordinal_type d1 = d-d2;
910 if (d1 > spaceDim1) continue;
911
912 ordinal_type subcellCount2 = cellTopo2.getSubcellCount(d2);
913 ordinal_type subcellCount1 = cellTopo1.getSubcellCount(d1);
914 for (ordinal_type subcellOrdinal2=0; subcellOrdinal2<subcellCount2; subcellOrdinal2++)
915 {
916 ordinal_type subcellDofCount2 = basis2_->getDofCount(d2, subcellOrdinal2);
917 for (ordinal_type subcellOrdinal1=0; subcellOrdinal1<subcellCount1; subcellOrdinal1++)
918 {
919 ordinal_type subcellDofCount1 = basis1_->getDofCount(d1, subcellOrdinal1);
920 ordinal_type tensorLocalDofCount = subcellDofCount1 * subcellDofCount2;
921 for (ordinal_type localDofID2 = 0; localDofID2<subcellDofCount2; localDofID2++)
922 {
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++)
927 {
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; // subcell dimension
932 tagView(tensorFieldOrdinal*tagSize+1) = topoMap.getCompositeSubcellOrdinal(d1, subcellOrdinal1, d2, subcellOrdinal2);
933 tagView(tensorFieldOrdinal*tagSize+2) = tensorLocalDofID;
934 tagView(tensorFieldOrdinal*tagSize+3) = tensorLocalDofCount;
935 } // localDofID1
936 } // localDofID2
937 } // subcellOrdinal1
938 } // subcellOrdinal2
939 subcellOffset += subcellCount1 * subcellCount2;
940 }
941 }
942
943 // // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
944 // // tags are constructed on host
945 this->setOrdinalTagData(this->tagToOrdinal_,
946 this->ordinalToTag_,
947 tagView,
948 this->basisCardinality_,
949 tagSize,
950 posScDim,
951 posScOrd,
952 posDfOrd);
953 }
954 }
955
956 virtual int getNumTensorialExtrusions() const override
957 {
958 return numTensorialExtrusions_;
959 }
960
969 ordinal_type getTensorDkEnumeration(ordinal_type dkEnum1, ordinal_type operatorOrder1,
970 ordinal_type dkEnum2, ordinal_type operatorOrder2) const
971 {
972 ordinal_type spaceDim1 = basis1_->getDomainDimension();
973 ordinal_type spaceDim2 = basis2_->getDomainDimension();
974
975 // We support total spaceDim <= 7.
976 switch (spaceDim1)
977 {
978 case 0:
979 {
980 INTREPID2_TEST_FOR_EXCEPTION(operatorOrder1 > 0, std::invalid_argument, "For spaceDim1 = 0, operatorOrder1 must be 0.");
981 return dkEnum2;
982 }
983 case 1:
984 switch (spaceDim2)
985 {
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);
992 default:
993 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported dimension combination");
994 }
995 case 2:
996 switch (spaceDim2)
997 {
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);
1003 default:
1004 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported dimension combination");
1005 }
1006 case 3:
1007 switch (spaceDim2)
1008 {
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);
1013 default:
1014 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported dimension combination");
1015 }
1016 case 4:
1017 switch (spaceDim2)
1018 {
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);
1022 default:
1023 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported dimension combination");
1024 }
1025 case 5:
1026 switch (spaceDim2)
1027 {
1028 case 1: return getDkTensorIndex<5, 1>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1029 case 2: return getDkTensorIndex<5, 2>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1030 default:
1031 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported dimension combination");
1032 }
1033 case 6:
1034 switch (spaceDim2)
1035 {
1036 case 1: return getDkTensorIndex<6, 1>(dkEnum1, operatorOrder1, dkEnum2, operatorOrder2);
1037 default:
1038 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported dimension combination");
1039 }
1040 default:
1041 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported dimension combination");
1042 }
1043 }
1044
1049 virtual OperatorTensorDecomposition getSimpleOperatorDecomposition(const EOperator &operatorType) const
1050 {
1051 const int spaceDim = this->getDomainDimension();
1052
1053 const EOperator VALUE = Intrepid2::OPERATOR_VALUE;
1054
1055 std::vector< std::vector<EOperator> > opsVALUE{{VALUE, VALUE}};
1056
1057 std::vector< std::vector<EOperator> > ops(spaceDim);
1058
1059 switch (operatorType)
1060 {
1061 case VALUE:
1062 ops = opsVALUE;
1063 break;
1064 case OPERATOR_DIV:
1065 case OPERATOR_CURL:
1066 // DIV and CURL are multi-family bases; subclasses are required to override
1067 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported operator type - TensorBasis subclass should override");
1068 break;
1069 case OPERATOR_GRAD:
1070 case OPERATOR_D1:
1071 case OPERATOR_D2:
1072 case OPERATOR_D3:
1073 case OPERATOR_D4:
1074 case OPERATOR_D5:
1075 case OPERATOR_D6:
1076 case OPERATOR_D7:
1077 case OPERATOR_D8:
1078 case OPERATOR_D9:
1079 case OPERATOR_D10:
1080 case OPERATOR_Dn:
1081 {
1082 auto opOrder = getOperatorOrder(operatorType); // number of derivatives that we take in total
1083 const int dkCardinality = getDkCardinality(operatorType, 2); // 2 because we have two tensor component bases, basis1_ and basis2_
1084
1085 ops = std::vector< std::vector<EOperator> >(dkCardinality);
1086
1087 // the Dk enumeration happens in lexicographic order (reading from left to right: x, y, z, etc.)
1088 // this governs the nesting order of the dkEnum1, dkEnum2 for loops below: dkEnum2 should increment fastest.
1089 for (int derivativeCountComp2=0; derivativeCountComp2<=opOrder; derivativeCountComp2++)
1090 {
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));
1094
1095 int dkCardinality1 = getDkCardinality(op1, 1); // use dim = 1 because this is a "simple" decomposition -- full decomposition will expand within the dimensions of basis1_
1096 int dkCardinality2 = getDkCardinality(op2, 1); // use dim = 1 because this is a "simple" decomposition -- full decomposition will expand within the dimensions of basis2_
1097
1098 for (int dkEnum1=0; dkEnum1<dkCardinality1; dkEnum1++)
1099 {
1100 for (int dkEnum2=0; dkEnum2<dkCardinality2; dkEnum2++)
1101 {
1102 ordinal_type dkTensorIndex = getDkTensorIndex<1, 1>(dkEnum1, derivativeCountComp1, dkEnum2, derivativeCountComp2);
1103 ops[dkTensorIndex] = std::vector<EOperator>{op1, op2};
1104 }
1105 }
1106 }
1107 }
1108 break;
1109 }
1110
1111 std::vector<double> weights(ops.size(), 1.0);
1112 return OperatorTensorDecomposition(ops, weights);
1113 }
1114
1117 virtual OperatorTensorDecomposition getOperatorDecomposition(const EOperator operatorType) const
1118 {
1119 if (((operatorType >= OPERATOR_D1) && (operatorType <= OPERATOR_D10)) || (operatorType == OPERATOR_GRAD))
1120 {
1121 // ordering of the operators is reverse-lexicographic, reading left to right (highest-dimension is fastest-moving).
1122 // first entry will be (operatorType, VALUE, …, VALUE)
1123 // next will be (operatorType - 1, OP_D1, VALUE, …, VALUE)
1124 // then (operatorType - 1, VALUE, OP_D1, …, VALUE)
1125
1126 ordinal_type numBasisComponents = tensorComponents_.size();
1127
1128 auto opOrder = getOperatorOrder(operatorType); // number of derivatives that we take in total
1129 const int dkCardinality = getDkCardinality(operatorType, numBasisComponents);
1130
1131 std::vector< std::vector<EOperator> > ops(dkCardinality);
1132
1133 std::vector<EOperator> prevEntry(numBasisComponents, OPERATOR_VALUE);
1134 prevEntry[0] = operatorType;
1135
1136 ops[0] = prevEntry;
1137
1138 for (ordinal_type dkOrdinal=1; dkOrdinal<dkCardinality; dkOrdinal++)
1139 {
1140 std::vector<EOperator> entry = prevEntry;
1141
1142 // decrement to follow reverse lexicographic ordering:
1143 /*
1144 How to tell when it is time to decrement the nth entry:
1145 1. Let a be the sum of the opOrders for entries 0 through n-1.
1146 2. Let b be the sum of the nth entry and the final entry.
1147 3. If opOrder == a + b, then the nth entry should be decremented.
1148 */
1149 ordinal_type cumulativeOpOrder = 0;
1150 ordinal_type finalOpOrder = getOperatorOrder(entry[numBasisComponents-1]);
1151 for (ordinal_type compOrdinal=0; compOrdinal<numBasisComponents; compOrdinal++)
1152 {
1153 const ordinal_type thisOpOrder = getOperatorOrder(entry[compOrdinal]);
1154 cumulativeOpOrder += thisOpOrder;
1155 if (cumulativeOpOrder + finalOpOrder == opOrder)
1156 {
1157 // decrement this
1158 EOperator decrementedOp;
1159 if (thisOpOrder == 1)
1160 {
1161 decrementedOp = OPERATOR_VALUE;
1162 }
1163 else
1164 {
1165 decrementedOp = static_cast<EOperator>(OPERATOR_D1 + ((thisOpOrder - 1) - 1));
1166 }
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++)
1171 {
1172 entry[i] = OPERATOR_VALUE;
1173 }
1174 break;
1175 }
1176 }
1177 ops[dkOrdinal] = entry;
1178 prevEntry = entry;
1179 }
1180 std::vector<double> weights(dkCardinality, 1.0);
1181
1182 return OperatorTensorDecomposition(ops, weights);
1183 }
1184 else
1185 {
1186 OperatorTensorDecomposition opSimpleDecomposition = this->getSimpleOperatorDecomposition(operatorType);
1187 std::vector<BasisPtr> componentBases {basis1_, basis2_};
1188 return opSimpleDecomposition.expandedDecomposition(componentBases);
1189 }
1190 }
1191
1196 virtual BasisValues<OutputValueType,DeviceType> allocateBasisValues( TensorPoints<PointValueType,DeviceType> points, const EOperator operatorType = OPERATOR_VALUE) const override
1197 {
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");
1201
1202 // check that points's spatial dimension matches the basis
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");
1205
1206 // check that points has enough tensor components
1207 ordinal_type numBasisComponents = tensorComponents_.size();
1208 if (numBasisComponents > points.numTensorComponents())
1209 {
1210 // Then we require points to have a trivial tensor structure. (Subclasses could be more sophisticated.)
1211 // (More sophisticated approaches are possible here, too, but likely the most common use case in which there is not a one-to-one correspondence
1212 // between basis components and point components will involve trivial tensor structure in the points...)
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);
1216
1217 Data<OutputValueType,DeviceType> outputData(outputView);
1218 TensorData<OutputValueType,DeviceType> outputTensorData(outputData);
1219
1220 return BasisValues<OutputValueType,DeviceType>(outputTensorData);
1221 }
1222 INTREPID2_TEST_FOR_EXCEPTION(numBasisComponents > points.numTensorComponents(), std::invalid_argument, "points must have at least as many tensorial components as basis.");
1223
1224 OperatorTensorDecomposition opDecomposition = getOperatorDecomposition(operatorType);
1225
1226 ordinal_type numVectorComponents = opDecomposition.numVectorComponents();
1227 const bool useVectorData = numVectorComponents > 1;
1228
1229 std::vector<ordinal_type> componentPointCounts(numBasisComponents);
1230 ordinal_type pointComponentNumber = 0;
1231 for (ordinal_type r=0; r<numBasisComponents; r++)
1232 {
1233 const ordinal_type compSpaceDim = tensorComponents_[r]->getDomainDimension();
1234 ordinal_type dimsSoFar = 0;
1235 ordinal_type numPointsForBasisComponent = 1;
1236 while (dimsSoFar < compSpaceDim)
1237 {
1238 INTREPID2_TEST_FOR_EXCEPTION(pointComponentNumber >= points.numTensorComponents(), std::invalid_argument, "Error in processing points container; perhaps it is mis-sized?");
1239 const int numComponentPoints = points.componentPointCount(pointComponentNumber);
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++;
1245 }
1246 componentPointCounts[r] = numPointsForBasisComponent;
1247 }
1248
1249 if (useVectorData)
1250 {
1251 const int numFamilies = 1;
1252 std::vector< std::vector<TensorData<OutputValueType,DeviceType> > > vectorComponents(numFamilies, std::vector<TensorData<OutputValueType,DeviceType> >(numVectorComponents));
1253
1254 const int familyOrdinal = 0;
1255 for (ordinal_type vectorComponentOrdinal=0; vectorComponentOrdinal<numVectorComponents; vectorComponentOrdinal++)
1256 {
1257 if (!opDecomposition.identicallyZeroComponent(vectorComponentOrdinal))
1258 {
1259 std::vector< Data<OutputValueType,DeviceType> > componentData;
1260 for (ordinal_type r=0; r<numBasisComponents; r++)
1261 {
1262 const int numComponentPoints = componentPointCounts[r];
1263 const EOperator op = opDecomposition.op(vectorComponentOrdinal, r);
1264 auto componentView = tensorComponents_[r]->allocateOutputView(numComponentPoints, op);
1265 componentData.push_back(Data<OutputValueType,DeviceType>(componentView));
1266 }
1267 vectorComponents[familyOrdinal][vectorComponentOrdinal] = TensorData<OutputValueType,DeviceType>(componentData);
1268 }
1269 }
1270 VectorData<OutputValueType,DeviceType> vectorData(vectorComponents);
1271 return BasisValues<OutputValueType,DeviceType>(vectorData);
1272 }
1273 else
1274 {
1275 // TensorData: single tensor product
1276 std::vector< Data<OutputValueType,DeviceType> > componentData;
1277
1278 const ordinal_type vectorComponentOrdinal = 0;
1279 for (ordinal_type r=0; r<numBasisComponents; r++)
1280 {
1281 const int numComponentPoints = componentPointCounts[r];
1282 const EOperator op = opDecomposition.op(vectorComponentOrdinal, r);
1283 auto componentView = tensorComponents_[r]->allocateOutputView(numComponentPoints, op);
1284
1285 const int rank = 2; // (F,P) -- TensorData-only BasisValues are always scalar-valued. Use VectorData for vector-valued.
1286 // (we need to be explicit about the rank argument because GRAD, even in 1D, elevates to rank 3), so e.g. DIV of HDIV uses a componentView that is rank 3;
1287 // we want Data to insulate us from that fact)
1288 const Kokkos::Array<int,7> extents {componentView.extent_int(0), componentView.extent_int(1), 1,1,1,1,1};
1289 Kokkos::Array<DataVariationType,7> variationType {GENERAL, GENERAL, CONSTANT, CONSTANT, CONSTANT, CONSTANT, CONSTANT };
1290 componentData.push_back(Data<OutputValueType,DeviceType>(componentView, rank, extents, variationType));
1291 }
1292
1293 TensorData<OutputValueType,DeviceType> tensorData(componentData);
1294
1295 std::vector< TensorData<OutputValueType,DeviceType> > tensorDataEntries {tensorData};
1296 return BasisValues<OutputValueType,DeviceType>(tensorDataEntries);
1297 }
1298 }
1299
1300 // since the getValues() below only overrides the FEM variant, we specify that
1301 // we use the base class's getValues(), which implements the FVD variant by throwing an exception.
1302 // (It's an error to use the FVD variant on this basis.)
1303 using BasisBase::getValues;
1304
1316 void getComponentPoints(const PointViewType inputPoints, const bool attemptTensorDecomposition,
1317 PointViewType & inputPoints1, PointViewType & inputPoints2, bool &tensorDecompositionSucceeded) const
1318 {
1319 INTREPID2_TEST_FOR_EXCEPTION(attemptTensorDecomposition, std::invalid_argument, "tensor decomposition not yet supported");
1320
1321 // for inputPoints that are actually tensor-product of component quadrature points (say),
1322 // having just the one input (which will have a lot of redundant point data) is suboptimal
1323 // The general case can have unique x/y/z coordinates at every point, though, so we have to support that
1324 // when this interface is used. But we may try detecting that the data is tensor-product and compressing
1325 // from there... Ultimately, we should also add a getValues() variant that takes multiple input point containers,
1326 // one for each tensorial dimension.
1327
1328 // this initial implementation is intended to simplify development of 2D and 3D bases, while also opening
1329 // the possibility of higher-dimensional bases. It is not necessarily optimized for speed/memory. There
1330 // are things we can do in this regard, which may become important for matrix-free computations wherein
1331 // basis values don't get stored but are computed dynamically.
1332
1333 int spaceDim1 = basis1_->getDomainDimension();
1334 int spaceDim2 = basis2_->getDomainDimension();
1335
1336 int totalSpaceDim = inputPoints.extent_int(1);
1337
1338 TEUCHOS_ASSERT(spaceDim1 + spaceDim2 == totalSpaceDim);
1339
1340 // first pass: just take subviews to get input points -- this will result in redundant computations when points are themselves tensor product (i.e., inputPoints itself contains redundant data)
1341
1342 inputPoints1 = Kokkos::subview(inputPoints,Kokkos::ALL(),std::make_pair(0,spaceDim1));
1343 inputPoints2 = Kokkos::subview(inputPoints,Kokkos::ALL(),std::make_pair(spaceDim1,totalSpaceDim));
1344
1345 // std::cout << "inputPoints : " << inputPoints.extent(0) << " x " << inputPoints.extent(1) << std::endl;
1346 // std::cout << "inputPoints1 : " << inputPoints1.extent(0) << " x " << inputPoints1.extent(1) << std::endl;
1347 // std::cout << "inputPoints2 : " << inputPoints2.extent(0) << " x " << inputPoints2.extent(1) << std::endl;
1348
1349 tensorDecompositionSucceeded = false;
1350 }
1351
1360 virtual void getDofCoords( typename BasisBase::ScalarViewType dofCoords ) const override
1361 {
1362 int spaceDim1 = basis1_->getBaseCellTopology().getDimension();
1363 int spaceDim2 = basis2_->getBaseCellTopology().getDimension();
1364
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 >;
1368
1369 const ordinal_type basisCardinality1 = basis1_->getCardinality();
1370 const ordinal_type basisCardinality2 = basis2_->getCardinality();
1371
1372 ViewType dofCoords1("dofCoords1",basisCardinality1,spaceDim1);
1373 ViewType dofCoords2("dofCoords2",basisCardinality2,spaceDim2);
1374
1375 basis1_->getDofCoords(dofCoords1);
1376 basis2_->getDofCoords(dofCoords2);
1377
1378 Kokkos::RangePolicy<ExecutionSpace> policy(0, basisCardinality2);
1379 Kokkos::parallel_for(policy, KOKKOS_LAMBDA (const int fieldOrdinal2)
1380 {
1381 for (int fieldOrdinal1=0; fieldOrdinal1<basisCardinality1; fieldOrdinal1++)
1382 {
1383 const ordinal_type fieldOrdinal = fieldOrdinal1 + fieldOrdinal2 * basisCardinality1;
1384 for (int d1=0; d1<spaceDim1; d1++)
1385 {
1386 dofCoords(fieldOrdinal,d1) = dofCoords1(fieldOrdinal1,d1);
1387 }
1388 for (int d2=0; d2<spaceDim2; d2++)
1389 {
1390 dofCoords(fieldOrdinal,spaceDim1+d2) = dofCoords2(fieldOrdinal2,d2);
1391 }
1392 }
1393 });
1394 }
1395
1396
1404 virtual void getDofCoeffs( typename BasisBase::ScalarViewType dofCoeffs ) const override
1405 {
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 >;
1409
1410 ViewType dofCoeffs1("dofCoeffs1",basis1_->getCardinality());
1411 ViewType dofCoeffs2("dofCoeffs2",basis2_->getCardinality());
1412
1413 basis1_->getDofCoeffs(dofCoeffs1);
1414 basis2_->getDofCoeffs(dofCoeffs2);
1415
1416 const ordinal_type basisCardinality1 = basis1_->getCardinality();
1417 const ordinal_type basisCardinality2 = basis2_->getCardinality();
1418
1419 Kokkos::RangePolicy<ExecutionSpace> policy(0, basisCardinality2);
1420 Kokkos::parallel_for(policy, KOKKOS_LAMBDA (const int fieldOrdinal2)
1421 {
1422 for (int fieldOrdinal1=0; fieldOrdinal1<basisCardinality1; fieldOrdinal1++)
1423 {
1424 const ordinal_type fieldOrdinal = fieldOrdinal1 + fieldOrdinal2 * basisCardinality1;
1425 dofCoeffs(fieldOrdinal) = dofCoeffs1(fieldOrdinal1);
1426 dofCoeffs(fieldOrdinal) *= dofCoeffs2(fieldOrdinal2);
1427 }
1428 });
1429 }
1430
1435 virtual
1436 const char*
1437 getName() const override {
1438 return name_.c_str();
1439 }
1440
1441 std::vector<BasisPtr> getTensorBasisComponents() const
1442 {
1443 return tensorComponents_;
1444 }
1445
1457 virtual
1458 void
1461 const EOperator operatorType = OPERATOR_VALUE ) const override
1462 {
1463 const ordinal_type numTensorComponents = tensorComponents_.size();
1464 if (inputPoints.numTensorComponents() < numTensorComponents)
1465 {
1466 // then we require that both inputPoints and outputValues trivial tensor structure
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" );
1470
1471 OutputViewType outputView = outputValues.tensorData().getTensorComponent(0).getUnderlyingView();
1472 PointViewType pointView = inputPoints.getTensorComponent(0);
1473 this->getValues(outputView, pointView, operatorType);
1474 return;
1475 }
1476
1477 OperatorTensorDecomposition operatorDecomposition = getOperatorDecomposition(operatorType);
1478
1479 const ordinal_type numVectorComponents = operatorDecomposition.numVectorComponents();
1480 const bool useVectorData = numVectorComponents > 1;
1481 const ordinal_type numBasisComponents = operatorDecomposition.numBasisComponents();
1482
1483 for (ordinal_type vectorComponentOrdinal=0; vectorComponentOrdinal<numVectorComponents; vectorComponentOrdinal++)
1484 {
1485 const double weight = operatorDecomposition.weight(vectorComponentOrdinal);
1486 ordinal_type pointComponentOrdinal = 0;
1487 for (ordinal_type basisOrdinal=0; basisOrdinal<numBasisComponents; basisOrdinal++, pointComponentOrdinal++)
1488 {
1489 const EOperator op = operatorDecomposition.op(vectorComponentOrdinal, basisOrdinal);
1490 // by convention, op == OPERATOR_MAX signals a zero component; skip
1491 if (op != OPERATOR_MAX)
1492 {
1493 const int vectorFamily = 0; // TensorBasis always has just a single family; multiple families arise in DirectSumBasis
1494 auto tensorData = useVectorData ? outputValues.vectorData().getComponent(vectorFamily,vectorComponentOrdinal) : outputValues.tensorData();
1495 INTREPID2_TEST_FOR_EXCEPTION( ! tensorData.getTensorComponent(basisOrdinal).isValid(), std::invalid_argument, "Invalid output component encountered");
1496
1497 const Data<OutputValueType,DeviceType> & outputData = tensorData.getTensorComponent(basisOrdinal);
1498
1499 auto basisValueView = outputData.getUnderlyingView();
1500 PointViewType pointView = inputPoints.getTensorComponent(pointComponentOrdinal);
1501 const ordinal_type basisDomainDimension = tensorComponents_[basisOrdinal]->getDomainDimension();
1502 if (pointView.extent_int(1) == basisDomainDimension)
1503 {
1504 tensorComponents_[basisOrdinal]->getValues(basisValueView, pointView, op);
1505 }
1506 else
1507 {
1508 // we need to wrap the basisValueView in a BasisValues container, and to wrap the point components in a TensorPoints container.
1509
1510 // combine point components to build up to basisDomainDimension
1511 ordinal_type dimsSoFar = 0;
1512 std::vector< ScalarView<PointValueType,DeviceType> > basisPointComponents;
1513 while (dimsSoFar < basisDomainDimension)
1514 {
1515 INTREPID2_TEST_FOR_EXCEPTION(pointComponentOrdinal >= inputPoints.numTensorComponents(), std::invalid_argument, "Error in processing points container; perhaps it is mis-sized?");
1516 const auto & pointComponent = inputPoints.getTensorComponent(pointComponentOrdinal);
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)
1522 {
1523 // we will pass through this loop again, so we should increment the point component ordinal
1524 pointComponentOrdinal++;
1525 }
1526 }
1527
1528 TensorPoints<PointValueType, DeviceType> basisPoints(basisPointComponents);
1529
1530 bool useVectorData2 = (basisValueView.rank() == 3);
1531
1533 if (useVectorData2)
1534 {
1535 VectorData<OutputValueType,DeviceType> vectorData(outputData);
1536 basisValues = BasisValues<OutputValueType,DeviceType>(vectorData);
1537 }
1538 else
1539 {
1540 TensorData<OutputValueType,DeviceType> tensorData2(outputData);
1541 basisValues = BasisValues<OutputValueType,DeviceType>(tensorData2);
1542 }
1543
1544 tensorComponents_[basisOrdinal]->getValues(basisValues, basisPoints, op);
1545 }
1546
1547 // op.rotateXYNinetyDegrees() is set to true for one of the H(curl) wedge families
1548 // (due to the fact that Intrepid2::EOperator does not allow us to extract individual vector components
1549 // via, e.g., OPERATOR_X, OPERATOR_Y, etc., we don't have a way of expressing the decomposition
1550 // just in terms of EOperator and component-wise scalar weights; we could also do this via component-wise
1551 // matrix weights, but this would involve a more intrusive change to the implementation).
1552 const bool spansXY = (vectorComponentOrdinal == 0) && (basisValueView.extent_int(2) == 2);
1553 if (spansXY && operatorDecomposition.rotateXYNinetyDegrees())
1554 {
1555 // map from (f_x,f_y) --> (-f_y,f_x)
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); // copy
1560 const auto &f_y = basisValueView(fieldOrdinal,pointOrdinal,1); // reference
1561 basisValueView(fieldOrdinal,pointOrdinal,0) = -f_y;
1562 basisValueView(fieldOrdinal,pointOrdinal,1) = f_x;
1563 });
1564 }
1565
1566 // if weight is non-trivial (not 1.0), then we need to multiply one of the component views by weight.
1567 // we do that for the first basisOrdinal's values
1568 if ((weight != 1.0) && (basisOrdinal == 0))
1569 {
1570 if (basisValueView.rank() == 2)
1571 {
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;
1576 });
1577 }
1578 else if (basisValueView.rank() == 3)
1579 {
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;
1584 });
1585 }
1586 else
1587 {
1588 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported rank for basisValueView");
1589 }
1590 }
1591 }
1592 }
1593 }
1594 }
1595
1614 void getValues( OutputViewType outputValues, const PointViewType inputPoints,
1615 const EOperator operatorType = OPERATOR_VALUE ) const override
1616 {
1617 bool tensorPoints; // true would mean that we take the tensor product of inputPoints1 and inputPoints2 (and that this would be equivalent to inputPoints as given -- i.e., inputPoints1 and inputPoints2 would be a tensor decomposition of inputPoints)
1618 bool attemptTensorDecomposition = false; // support for this not yet implemented
1619 PointViewType inputPoints1, inputPoints2;
1620 getComponentPoints(inputPoints, attemptTensorDecomposition, inputPoints1, inputPoints2, tensorPoints);
1621
1622 const auto functionSpace = this->getFunctionSpace();
1623
1624 if ((functionSpace == FUNCTION_SPACE_HVOL) || (functionSpace == FUNCTION_SPACE_HGRAD))
1625 {
1626 // then we can handle VALUE, GRAD, and Op_Dn without reference to subclass
1627 switch (operatorType)
1628 {
1629 case OPERATOR_VALUE:
1630 case OPERATOR_GRAD:
1631 case OPERATOR_D1:
1632 case OPERATOR_D2:
1633 case OPERATOR_D3:
1634 case OPERATOR_D4:
1635 case OPERATOR_D5:
1636 case OPERATOR_D6:
1637 case OPERATOR_D7:
1638 case OPERATOR_D8:
1639 case OPERATOR_D9:
1640 case OPERATOR_D10:
1641 {
1642 auto opOrder = getOperatorOrder(operatorType); // number of derivatives that we take in total
1643 // the Dk enumeration happens in lexicographic order (reading from left to right: x, y, z, etc.)
1644 // this governs the nesting order of the dkEnum1, dkEnum2 for loops below: dkEnum2 should increment fastest.
1645 for (int derivativeCountComp2=0; derivativeCountComp2<=opOrder; derivativeCountComp2++)
1646 {
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));
1650
1651 int spaceDim1 = inputPoints1.extent_int(1);
1652 int spaceDim2 = inputPoints2.extent_int(1);
1653
1654 int dkCardinality1 = (op1 != OPERATOR_VALUE) ? getDkCardinality(op1, spaceDim1) : 1;
1655 int dkCardinality2 = (op2 != OPERATOR_VALUE) ? getDkCardinality(op2, spaceDim2) : 1;
1656
1657 int basisCardinality1 = basis1_->getCardinality();
1658 int basisCardinality2 = basis2_->getCardinality();
1659
1660 int totalPointCount = tensorPoints ? inputPoints1.extent_int(0) * inputPoints2.extent_int(0) : inputPoints1.extent_int(0);
1661
1662 int pointCount1, pointCount2;
1663 if (tensorPoints)
1664 {
1665 pointCount1 = inputPoints1.extent_int(0);
1666 pointCount2 = inputPoints2.extent_int(0);
1667 }
1668 else
1669 {
1670 pointCount1 = totalPointCount;
1671 pointCount2 = totalPointCount;
1672 }
1673
1674 OutputViewType outputValues1, outputValues2;
1675 if (op1 == OPERATOR_VALUE)
1676 outputValues1 = getMatchingViewWithLabel(outputValues, "output values - basis 1",basisCardinality1,pointCount1);
1677 else
1678 outputValues1 = getMatchingViewWithLabel(outputValues, "output values - basis 1",basisCardinality1,pointCount1,dkCardinality1);
1679
1680 if (op2 == OPERATOR_VALUE)
1681 outputValues2 = getMatchingViewWithLabel(outputValues, "output values - basis 2",basisCardinality2,pointCount2);
1682 else
1683 outputValues2 = getMatchingViewWithLabel(outputValues, "output values - basis 2",basisCardinality2,pointCount2,dkCardinality2);
1684
1685 basis1_->getValues(outputValues1,inputPoints1,op1);
1686 basis2_->getValues(outputValues2,inputPoints2,op2);
1687
1688 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputValueType>();
1689 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointValueType>();
1690 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
1691
1692 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(basisCardinality1,Kokkos::AUTO(),vectorSize);
1693
1694 double weight = 1.0;
1696
1697 for (int dkEnum1=0; dkEnum1<dkCardinality1; dkEnum1++)
1698 {
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++)
1702 {
1703 auto outputValues2_dkEnum2 = (op2 != OPERATOR_VALUE) ? Kokkos::subview(outputValues2,Kokkos::ALL(),Kokkos::ALL(),dkEnum2)
1704 : Kokkos::subview(outputValues2,Kokkos::ALL(),Kokkos::ALL());
1705
1706 ordinal_type dkTensorIndex = getTensorDkEnumeration(dkEnum1, derivativeCountComp1, dkEnum2, derivativeCountComp2);
1707 auto outputValues_dkTensor = Kokkos::subview(outputValues,Kokkos::ALL(),Kokkos::ALL(),dkTensorIndex);
1708 // Note that there may be performance optimizations available here:
1709 // - could eliminate interior for loop in favor of having a vector-valued outputValues1_dk
1710 // - could add support to TensorViewFunctor (and probably TensorViewIterator) for this kind of tensor Dk type of traversal
1711 // (this would allow us to eliminate both for loops here)
1712 // At the moment, we defer such optimizations on the idea that this may not ever become a performance bottleneck.
1713 FunctorType functor(outputValues_dkTensor, outputValues1_dkEnum1, outputValues2_dkEnum2, tensorPoints, weight);
1714 Kokkos::parallel_for("TensorViewFunctor", policy , functor);
1715 }
1716 }
1717 }
1718 }
1719 break;
1720 default: // non-OPERATOR_Dn case must be handled by subclass.
1721 this->getValues(outputValues, operatorType, inputPoints1, inputPoints2, tensorPoints);
1722 }
1723 }
1724 else
1725 {
1726 // not HVOL or HGRAD; subclass must handle
1727 this->getValues(outputValues, operatorType, inputPoints1, inputPoints2, tensorPoints);
1728 }
1729 }
1730
1756 virtual void getValues(OutputViewType outputValues, const EOperator operatorType,
1757 const PointViewType inputPoints1, const PointViewType inputPoints2,
1758 bool tensorPoints) const
1759 {
1760 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "one-operator, two-inputPoints getValues should be overridden by TensorBasis subclasses");
1761 }
1762
1786 void getValues( OutputViewType outputValues,
1787 const PointViewType inputPoints1, const EOperator operatorType1,
1788 const PointViewType inputPoints2, const EOperator operatorType2,
1789 bool tensorPoints, double weight=1.0) const
1790 {
1791 int basisCardinality1 = basis1_->getCardinality();
1792 int basisCardinality2 = basis2_->getCardinality();
1793
1794 int totalPointCount = tensorPoints ? inputPoints1.extent_int(0) * inputPoints2.extent_int(0) : inputPoints1.extent_int(0);
1795
1796 int pointCount1, pointCount2;
1797 if (tensorPoints)
1798 {
1799 pointCount1 = inputPoints1.extent_int(0);
1800 pointCount2 = inputPoints2.extent_int(0);
1801 }
1802 else
1803 {
1804 pointCount1 = totalPointCount;
1805 pointCount2 = totalPointCount;
1806 }
1807
1808 const ordinal_type spaceDim1 = inputPoints1.extent_int(1);
1809 const ordinal_type spaceDim2 = inputPoints2.extent_int(1);
1810
1811 INTREPID2_TEST_FOR_EXCEPTION(!tensorPoints && (totalPointCount != inputPoints2.extent_int(0)),
1812 std::invalid_argument, "If tensorPoints is false, the point counts must match!");
1813
1814 const ordinal_type opRank1 = getOperatorRank(basis1_->getFunctionSpace(), operatorType1, spaceDim1);
1815 const ordinal_type opRank2 = getOperatorRank(basis2_->getFunctionSpace(), operatorType2, spaceDim2);
1816
1817 const ordinal_type outputRank1 = opRank1 + getFieldRank(basis1_->getFunctionSpace());
1818 const ordinal_type outputRank2 = opRank2 + getFieldRank(basis2_->getFunctionSpace());
1819
1820 OutputViewType outputValues1, outputValues2;
1821 if (outputRank1 == 0)
1822 {
1823 outputValues1 = getMatchingViewWithLabel(outputValues,"output values - basis 1",basisCardinality1,pointCount1);
1824 }
1825 else if (outputRank1 == 1)
1826 {
1827 outputValues1 = getMatchingViewWithLabel(outputValues,"output values - basis 1",basisCardinality1,pointCount1,spaceDim1);
1828 }
1829 else
1830 {
1831 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported opRank1");
1832 }
1833
1834 if (outputRank2 == 0)
1835 {
1836 outputValues2 = getMatchingViewWithLabel(outputValues,"output values - basis 2",basisCardinality2,pointCount2);
1837 }
1838 else if (outputRank2 == 1)
1839 {
1840 outputValues2 = getMatchingViewWithLabel(outputValues,"output values - basis 2",basisCardinality2,pointCount2,spaceDim2);
1841 }
1842 else
1843 {
1844 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported opRank2");
1845 }
1846
1847 basis1_->getValues(outputValues1,inputPoints1,operatorType1);
1848 basis2_->getValues(outputValues2,inputPoints2,operatorType2);
1849
1850 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputValueType>();
1851 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointValueType>();
1852 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
1853
1854 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(basisCardinality1,Kokkos::AUTO(),vectorSize);
1855
1857
1858 FunctorType functor(outputValues, outputValues1, outputValues2, tensorPoints, weight);
1859 Kokkos::parallel_for("TensorViewFunctor", policy , functor);
1860 }
1861
1867 getHostBasis() const override {
1868 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "TensorBasis subclasses must override getHostBasis");
1869 }
1870 }; // Basis_TensorBasis
1871
1879 template<class ExecutionSpace, class OutputScalar, class OutputFieldType>
1881 {
1882 using ScratchSpace = typename ExecutionSpace::scratch_memory_space;
1883 using OutputScratchView = Kokkos::View<OutputScalar*,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
1884
1885 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
1886 using TeamMember = typename TeamPolicy::member_type;
1887
1888 OutputFieldType output_; // F,P
1889 OutputFieldType input1_; // F1,P[,D] or F1,P1[,D]
1890 OutputFieldType input2_; // F2,P[,D] or F2,P2[,D]
1891 OutputFieldType input3_; // F2,P[,D] or F2,P2[,D]
1892
1893 int numFields_, numPoints_;
1894 int numFields1_, numPoints1_;
1895 int numFields2_, numPoints2_;
1896 int numFields3_, numPoints3_;
1897
1898 bool tensorPoints_; // if true, input1, input2, input3 refer to values at decomposed points, and P = P1 * P2 * P3. If false, then the three inputs refer to points in the full-dimensional space, and their point lengths are the same as that of the final output.
1899
1900 double weight_;
1901
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)
1905 {
1906 numFields_ = output.extent_int(0);
1907 numPoints_ = output.extent_int(1);
1908
1909 numFields1_ = inputValues1.extent_int(0);
1910 numPoints1_ = inputValues1.extent_int(1);
1911
1912 numFields2_ = inputValues2.extent_int(0);
1913 numPoints2_ = inputValues2.extent_int(1);
1914
1915 numFields3_ = inputValues3.extent_int(0);
1916 numPoints3_ = inputValues3.extent_int(1);
1917 /*
1918 We don't yet support tensor-valued bases here (only vector and scalar). The main design question is how the layouts
1919 of the input containers relates to the layout of the output container. The work we've done in TensorViewIterator basically
1920 shows the choices that can be made. It does appear that in most cases (at least (most of?) those supported by TensorViewIterator),
1921 we can infer from the dimensions of input/output containers what choice should be made in each dimension.
1922 */
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");
1929
1930 if (!tensorPoints_)
1931 {
1932 // then the point counts should all match
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");
1936 }
1937 else
1938 {
1939 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != numPoints1_ * numPoints2_ * numPoints3_, std::invalid_argument, "incompatible point counts");
1940 }
1941
1942 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != numFields1_ * numFields2_ * numFields3_, std::invalid_argument, "incompatible field sizes");
1943 }
1944
1945 KOKKOS_INLINE_FUNCTION
1946 void operator()( const TeamMember & teamMember ) const
1947 {
1948 auto fieldOrdinal1 = teamMember.league_rank();
1949
1950 if (!tensorPoints_)
1951 {
1952 if ((input1_.rank() == 2) && (input2_.rank() == 2) && (input3_.rank() == 2))
1953 {
1954 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (const int& fieldOrdinal2) {
1955 for (int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
1956 {
1957 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
1958 for (int pointOrdinal=0; pointOrdinal<numPoints_; pointOrdinal++)
1959 {
1960 output_(fieldOrdinal,pointOrdinal) = weight_ * input1_(fieldOrdinal1,pointOrdinal) * input2_(fieldOrdinal2,pointOrdinal) * input3_(fieldOrdinal3,pointOrdinal);
1961 }
1962 }
1963 });
1964 }
1965 else if (input1_.rank() == 3)
1966 {
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++)
1970 {
1971 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
1972 for (int pointOrdinal=0; pointOrdinal<numPoints_; pointOrdinal++)
1973 {
1974 for (int d=0; d<spaceDim; d++)
1975 {
1976 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal,d) * input2_(fieldOrdinal2,pointOrdinal) * input3_(fieldOrdinal3,pointOrdinal);
1977 }
1978 }
1979 }
1980 });
1981 }
1982 else if (input2_.rank() == 3)
1983 {
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++)
1987 {
1988 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
1989 for (int pointOrdinal=0; pointOrdinal<numPoints_; pointOrdinal++)
1990 {
1991 for (int d=0; d<spaceDim; d++)
1992 {
1993 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal) * input2_(fieldOrdinal2,pointOrdinal,d) * input3_(fieldOrdinal3,pointOrdinal);
1994 }
1995 }
1996 }
1997 });
1998 }
1999 else if (input3_.rank() == 3)
2000 {
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++)
2004 {
2005 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
2006 for (int pointOrdinal=0; pointOrdinal<numPoints_; pointOrdinal++)
2007 {
2008 for (int d=0; d<spaceDim; d++)
2009 {
2010 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal) * input2_(fieldOrdinal2,pointOrdinal) * input3_(fieldOrdinal3,pointOrdinal,d);
2011 }
2012 }
2013 }
2014 });
2015 }
2016 else
2017 {
2018 // unsupported rank combination -- enforced in constructor
2019 }
2020 }
2021 else
2022 {
2023 if ((input1_.rank() == 2) && (input2_.rank() == 2) && (input3_.rank() == 2) )
2024 {
2025 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,0,numFields2_), [&] (const int& fieldOrdinal2) {
2026 for (int fieldOrdinal3=0; fieldOrdinal3 < numFields3_; fieldOrdinal3++)
2027 {
2028 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
2029 for (int pointOrdinal3=0; pointOrdinal3<numPoints3_; pointOrdinal3++)
2030 {
2031 for (int pointOrdinal2=0; pointOrdinal2<numPoints2_; pointOrdinal2++)
2032 {
2033 for (int pointOrdinal1=0; pointOrdinal1<numPoints1_; pointOrdinal1++)
2034 {
2035 int pointOrdinal = (pointOrdinal3 * numPoints2_ + pointOrdinal2) * numPoints1_ + pointOrdinal1;
2036 output_(fieldOrdinal,pointOrdinal) = weight_ * input1_(fieldOrdinal1,pointOrdinal1) * input2_(fieldOrdinal2,pointOrdinal2) * input3_(fieldOrdinal3,pointOrdinal3);
2037 }
2038 }
2039 }
2040 }
2041 });
2042 }
2043 else if (input1_.rank() == 3) // based on constructor requirements, this means the others are rank 2
2044 {
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++)
2048 {
2049 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
2050 for (int pointOrdinal3=0; pointOrdinal3<numPoints3_; pointOrdinal3++)
2051 {
2052 for (int pointOrdinal2=0; pointOrdinal2<numPoints2_; pointOrdinal2++)
2053 {
2054 for (int pointOrdinal1=0; pointOrdinal1<numPoints1_; pointOrdinal1++)
2055 {
2056 int pointOrdinal = (pointOrdinal3 * numPoints2_ + pointOrdinal2) * numPoints1_ + pointOrdinal1;
2057 for (int d=0; d<spaceDim; d++)
2058 {
2059 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal1,d) * input2_(fieldOrdinal2,pointOrdinal2) * input3_(fieldOrdinal3,pointOrdinal3);
2060 }
2061 }
2062 }
2063 }
2064 }
2065 });
2066 }
2067 else if (input2_.rank() == 3) // based on constructor requirements, this means the others are rank 2
2068 {
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++)
2072 {
2073 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
2074 for (int pointOrdinal3=0; pointOrdinal3<numPoints3_; pointOrdinal3++)
2075 {
2076 for (int pointOrdinal2=0; pointOrdinal2<numPoints2_; pointOrdinal2++)
2077 {
2078 for (int pointOrdinal1=0; pointOrdinal1<numPoints1_; pointOrdinal1++)
2079 {
2080 int pointOrdinal = (pointOrdinal3 * numPoints2_ + pointOrdinal2) * numPoints1_ + pointOrdinal1;
2081 for (int d=0; d<spaceDim; d++)
2082 {
2083 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal1) * input2_(fieldOrdinal2,pointOrdinal2,d) * input3_(fieldOrdinal3,pointOrdinal3);
2084 }
2085 }
2086 }
2087 }
2088 }
2089 });
2090 }
2091 else if (input3_.rank() == 3) // based on constructor requirements, this means the others are rank 2
2092 {
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++)
2096 {
2097 int fieldOrdinal = (fieldOrdinal3 * numFields2_ + fieldOrdinal2) * numFields1_ + fieldOrdinal1;
2098 for (int pointOrdinal3=0; pointOrdinal3<numPoints3_; pointOrdinal3++)
2099 {
2100 for (int pointOrdinal2=0; pointOrdinal2<numPoints2_; pointOrdinal2++)
2101 {
2102 for (int pointOrdinal1=0; pointOrdinal1<numPoints1_; pointOrdinal1++)
2103 {
2104 int pointOrdinal = (pointOrdinal3 * numPoints2_ + pointOrdinal2) * numPoints1_ + pointOrdinal1;
2105 for (int d=0; d<spaceDim; d++)
2106 {
2107 output_(fieldOrdinal,pointOrdinal,d) = weight_ * input1_(fieldOrdinal1,pointOrdinal1) * input2_(fieldOrdinal2,pointOrdinal2) * input3_(fieldOrdinal3,pointOrdinal3,d);
2108 }
2109 }
2110 }
2111 }
2112 }
2113 });
2114 }
2115 else
2116 {
2117 // unsupported rank combination -- enforced in constructor
2118 }
2119 }
2120 }
2121 }; // TensorBasis3_Functor
2122
2123
2124 template<typename BasisBaseClass = void>
2126 : public Basis_TensorBasis<BasisBaseClass>
2127 {
2128 using BasisBase = BasisBaseClass;
2130 public:
2131 using typename BasisBase::OutputViewType;
2132 using typename BasisBase::PointViewType;
2133 using typename BasisBase::ScalarViewType;
2134
2135 using typename BasisBase::OutputValueType;
2136 using typename BasisBase::PointValueType;
2137
2138 using typename BasisBase::ExecutionSpace;
2139
2140 using BasisPtr = Teuchos::RCP<BasisBase>;
2141 protected:
2142 BasisPtr basis1_;
2143 BasisPtr basis2_;
2144 BasisPtr basis3_;
2145 public:
2146 Basis_TensorBasis3(BasisPtr basis1, BasisPtr basis2, BasisPtr basis3, const bool useShardsCellTopologyAndTags = false)
2147 :
2148 TensorBasis(Teuchos::rcp( new TensorBasis(basis1,basis2,FUNCTION_SPACE_MAX,useShardsCellTopologyAndTags)),
2149 basis3,
2150 FUNCTION_SPACE_MAX,useShardsCellTopologyAndTags),
2151 basis1_(basis1),
2152 basis2_(basis2),
2153 basis3_(basis3)
2154 {}
2155
2162 virtual OperatorTensorDecomposition getOperatorDecomposition(const EOperator operatorType) const override
2163 {
2164 OperatorTensorDecomposition opSimpleDecomposition = this->getSimpleOperatorDecomposition(operatorType);
2165 std::vector<BasisPtr> componentBases {basis1_, basis2_, basis3_};
2166 return opSimpleDecomposition.expandedDecomposition(componentBases);
2167 }
2168
2170
2195 virtual void getValues(OutputViewType outputValues, const EOperator operatorType,
2196 const PointViewType inputPoints12, const PointViewType inputPoints3,
2197 bool tensorPoints) const override
2198 {
2199 // TODO: rework this to use superclass's getComponentPoints.
2200
2201 int spaceDim1 = basis1_->getDomainDimension();
2202 int spaceDim2 = basis2_->getDomainDimension();
2203
2204 int totalSpaceDim12 = inputPoints12.extent_int(1);
2205
2206 TEUCHOS_ASSERT(spaceDim1 + spaceDim2 == totalSpaceDim12);
2207
2208 if (!tensorPoints)
2209 {
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));
2212
2213 this->getValues(outputValues, operatorType, inputPoints1, inputPoints2, inputPoints3, tensorPoints);
2214 }
2215 else
2216 {
2217 // superclass doesn't (yet) have a clever way to detect tensor points in a single container
2218 // we'd need something along those lines here to detect them in inputPoints12.
2219 // if we do add such a mechanism to superclass, it should be simple enough to call that from here
2220 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "This method does not yet handle tensorPoints=true");
2221 }
2222 }
2223
2251 virtual void getValues(OutputViewType outputValues, const EOperator operatorType,
2252 const PointViewType inputPoints1, const PointViewType inputPoints2, const PointViewType inputPoints3,
2253 bool tensorPoints) const = 0;
2254
2282 void getValues( OutputViewType outputValues,
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
2287 {
2288 int basisCardinality1 = basis1_->getCardinality();
2289 int basisCardinality2 = basis2_->getCardinality();
2290 int basisCardinality3 = basis3_->getCardinality();
2291
2292 int spaceDim1 = inputPoints1.extent_int(1);
2293 int spaceDim2 = inputPoints2.extent_int(1);
2294 int spaceDim3 = inputPoints3.extent_int(1);
2295
2296 int totalPointCount;
2297 int pointCount1, pointCount2, pointCount3;
2298 if (tensorPoints)
2299 {
2300 pointCount1 = inputPoints1.extent_int(0);
2301 pointCount2 = inputPoints2.extent_int(0);
2302 pointCount3 = inputPoints3.extent_int(0);
2303 totalPointCount = pointCount1 * pointCount2 * pointCount3;
2304 }
2305 else
2306 {
2307 totalPointCount = inputPoints1.extent_int(0);
2308 pointCount1 = totalPointCount;
2309 pointCount2 = totalPointCount;
2310 pointCount3 = totalPointCount;
2311
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!");
2314 }
2315
2316 // structure of this implementation:
2317 /*
2318 - allocate output1, output2, output3 containers
2319 - either:
2320 1. split off the tensor functor call into its own method in TensorBasis, and
2321 - call it once with output1, output2, placing these in another newly allocated output12, then
2322 - call it again with output12, output3
2323 OR
2324 2. create a 3-argument tensor functor and call it with output1,output2,output3
2325
2326 At the moment, the 3-argument functor seems like a better approach. It's likely more code, but somewhat
2327 more efficient and easier to understand/debug. And the code is fairly straightforward to produce.
2328 */
2329
2330 // copied from the 2-argument TensorBasis implementation:
2331
2332 OutputViewType outputValues1, outputValues2, outputValues3;
2333
2334 //Note: the gradient of HGRAD basis on a line has an output vector of rank 3, the last dimension being of size 1.
2335 // in particular this holds even when computing the divergence of an HDIV basis, which is scalar and has rank 2.
2336 if ((spaceDim1 == 1) && (operatorType1 == OPERATOR_VALUE))
2337 {
2338 // use a rank 2 container for basis1
2339 outputValues1 = getMatchingViewWithLabel(outputValues,"output values - basis 1",basisCardinality1,pointCount1);
2340 }
2341 else
2342 {
2343 outputValues1 = getMatchingViewWithLabel(outputValues,"output values - basis 1",basisCardinality1,pointCount1,spaceDim1);
2344 }
2345 if ((spaceDim2 == 1) && (operatorType2 == OPERATOR_VALUE))
2346 {
2347 // use a rank 2 container for basis2
2348 outputValues2 = getMatchingViewWithLabel(outputValues,"output values - basis 2",basisCardinality2,pointCount2);
2349 }
2350 else
2351 {
2352 outputValues2 = getMatchingViewWithLabel(outputValues,"output values - basis 2",basisCardinality2,pointCount2,spaceDim2);
2353 }
2354 if ((spaceDim3 == 1) && (operatorType3 == OPERATOR_VALUE))
2355 {
2356 // use a rank 2 container for basis2
2357 outputValues3 = getMatchingViewWithLabel(outputValues,"output values - basis 3",basisCardinality3,pointCount3);
2358 }
2359 else
2360 {
2361 outputValues3 = getMatchingViewWithLabel(outputValues,"output values - basis 3",basisCardinality3,pointCount3,spaceDim3);
2362 }
2363
2364 basis1_->getValues(outputValues1,inputPoints1,operatorType1);
2365 basis2_->getValues(outputValues2,inputPoints2,operatorType2);
2366 basis3_->getValues(outputValues3,inputPoints3,operatorType3);
2367
2368 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputValueType>();
2369 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointValueType>();
2370 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
2371
2372 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(basisCardinality1,Kokkos::AUTO(),vectorSize);
2373
2375 FunctorType functor(outputValues, outputValues1, outputValues2, outputValues3, tensorPoints, weight);
2376 Kokkos::parallel_for("TensorBasis3_Functor", policy , functor);
2377 }
2378
2384 getHostBasis() const override {
2385 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "TensorBasis3 subclasses must override getHostBasis");
2386 }
2387 };
2388} // end namespace Intrepid2
2389
2390#endif /* Intrepid2_TensorBasis_h */
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.
@ CONSTANT
does not vary
@ 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 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.