Intrepid2
Intrepid2_TensorData.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
50#ifndef Intrepid2_TensorData_h
51#define Intrepid2_TensorData_h
52
53#include "Intrepid2_Data.hpp"
54
55#include "Intrepid2_ScalarView.hpp"
56#include "Intrepid2_Types.hpp"
57
58namespace Intrepid2
59{
63 template<class Scalar, typename DeviceType>
64 class TensorData {
65 protected:
66 Kokkos::Array< Data<Scalar,DeviceType>, Parameters::MaxTensorComponents> tensorComponents_;
67 Kokkos::Array<ordinal_type, 7> extents_;
68 Kokkos::Array<Kokkos::Array<ordinal_type, Parameters::MaxTensorComponents>, 7> entryModulus_;
69 ordinal_type rank_;
70 bool separateFirstComponent_ = false; // true supported only for rank 1 components; uses a rank-2 operator() in that case, where the first argument corresponds to the index in the first component, while the second argument corresponds to the tensorially-multiplied remaining components
71 ordinal_type numTensorComponents_ = 0;
72
77 {
78 ordinal_type maxComponentRank = -1;
79 for (ordinal_type r=0; r<numTensorComponents_; r++)
80 {
81 const ordinal_type componentRank = tensorComponents_[r].rank();
82 maxComponentRank = (maxComponentRank > componentRank) ? maxComponentRank : componentRank;
83 }
84 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(separateFirstComponent_ && (maxComponentRank != 1), std::invalid_argument, "separateFirstComponent = true only supported if all components have rank 1");
85 ordinal_type rd_start; // where to begin the extents_ and entryModulus_ loops below
86 if ((maxComponentRank == 1) && separateFirstComponent_)
87 {
88 rank_ = 2;
89 rd_start = 1;
90 extents_[0] = tensorComponents_[0].extent_int(0);
91 entryModulus_[0][0] = extents_[0]; // should not be used
92 }
93 else
94 {
95 rank_ = maxComponentRank;
96 rd_start = 0;
97 }
98
99 for (ordinal_type d=rd_start; d<7; d++)
100 {
101 extents_[d] = 1;
102 for (ordinal_type r=rd_start; r<numTensorComponents_; r++)
103 {
104 extents_[d] *= tensorComponents_[r].extent_int(d-rd_start);
105 }
106 ordinal_type entryModulus = extents_[d];
107 for (ordinal_type r=rd_start; r<numTensorComponents_; r++)
108 {
109 entryModulus /= tensorComponents_[r].extent_int(d-rd_start);
110 entryModulus_[d][r] = entryModulus;
111 }
112 }
113 }
114 public:
124 template<size_t numTensorComponents>
125 TensorData(Kokkos::Array< Data<Scalar,DeviceType>, numTensorComponents> tensorComponents, bool separateFirstComponent = false)
126 :
127 separateFirstComponent_(separateFirstComponent),
128 numTensorComponents_(numTensorComponents)
129 {
130 for (size_t r=0; r< numTensorComponents; r++)
131 {
132 tensorComponents_[r] = tensorComponents[r];
133 }
134
135 initialize();
136 }
137
147 TensorData(std::vector< Data<Scalar,DeviceType> > tensorComponents, bool separateFirstComponent = false)
148 :
149 separateFirstComponent_(separateFirstComponent),
150 numTensorComponents_(tensorComponents.size())
151 {
152 for (ordinal_type r=0; r<numTensorComponents_; r++)
153 {
154 tensorComponents_[r] = tensorComponents[r];
155 }
156
157 initialize();
158 }
159
168 TensorData(const TensorData &first, const TensorData &second, bool separateFirstComponent = false)
169 :
170 separateFirstComponent_(separateFirstComponent),
171 numTensorComponents_(first.numTensorComponents() + second.numTensorComponents())
172 {
173 ordinal_type r = 0;
174 for (ordinal_type r1=0; r1<first.numTensorComponents(); r1++, r++)
175 {
176 tensorComponents_[r] = first.getTensorComponent(r1);
177 }
178 for (ordinal_type r2=0; r2<second.numTensorComponents(); r2++, r++)
179 {
180 tensorComponents_[r] = second.getTensorComponent(r2);
181 }
182
183 initialize();
184 }
185
193 :
194 TensorData(Kokkos::Array< Data<Scalar,DeviceType>, 1>({tensorComponent}), false)
195 {}
196
203 :
204 extents_({0,0,0,0,0,0,0}),
205 rank_(0)
206 {}
207
215 TensorData(TensorData otherTensorData, std::vector<int> whichComps)
216 :
217 numTensorComponents_(whichComps.size())
218 {
219 int r = 0;
220 for (const auto & componentOrdinal : whichComps)
221 {
222 tensorComponents_[r++] = otherTensorData.getTensorComponent(componentOrdinal);
223 }
224
225 initialize();
226 }
227
229 template<typename OtherDeviceType, class = typename std::enable_if< std::is_same<typename DeviceType::memory_space, typename OtherDeviceType::memory_space>::value>::type,
230 class = typename std::enable_if<!std::is_same<DeviceType,OtherDeviceType>::value>::type>
231 TensorData(const TensorData<Scalar,OtherDeviceType> &tensorData)
232 {
233 if (tensorData.isValid())
234 {
235 numTensorComponents_ = tensorData.numTensorComponents();
236 for (ordinal_type r=0; r<numTensorComponents_; r++)
237 {
238 Data<Scalar,OtherDeviceType> otherTensorComponent = tensorData.getTensorComponent(r);
239 tensorComponents_[r] = Data<Scalar,DeviceType>(otherTensorComponent);
240 }
241 initialize();
242 }
243 else
244 {
245 extents_ = Kokkos::Array<ordinal_type,7>{0,0,0,0,0,0,0};
246 rank_ = 0;
247 }
248 }
249
253 template<typename OtherDeviceType, class = typename std::enable_if<!std::is_same<typename DeviceType::memory_space, typename OtherDeviceType::memory_space>::value>::type>
255 {
256 if (tensorData.isValid())
257 {
258 numTensorComponents_ = tensorData.numTensorComponents();
259 for (ordinal_type r=0; r<numTensorComponents_; r++)
260 {
261 Data<Scalar,OtherDeviceType> otherTensorComponent = tensorData.getTensorComponent(r);
262 tensorComponents_[r] = Data<Scalar,DeviceType>(otherTensorComponent);
263 }
264 initialize();
265 }
266 else
267 {
268 extents_ = Kokkos::Array<ordinal_type,7>{0,0,0,0,0,0,0};
269 rank_ = 0;
270 }
271 }
272
278 KOKKOS_INLINE_FUNCTION
279 const Data<Scalar,DeviceType> & getTensorComponent(const ordinal_type &r) const
280 {
281 return tensorComponents_[r];
282 }
283
289 template <typename iType0>
290 KOKKOS_INLINE_FUNCTION typename std::enable_if<std::is_integral<iType0>::value, Scalar>::type
291 operator()(const iType0& tensorEntryIndex) const {
292#ifdef HAVE_INTREPID2_DEBUG
293 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(rank_ != 1, std::invalid_argument, "This method is only valid for rank 1 containers.");
294#endif
295 Scalar value = 1.0;
296 iType0 remainingEntryOrdinal = tensorEntryIndex;
297 for (ordinal_type r=0; r<numTensorComponents_; r++)
298 {
299 const ordinal_type componentEntryCount = tensorComponents_[r].extent_int(0);
300 const ordinal_type componentEntryOrdinal = remainingEntryOrdinal % componentEntryCount;
301 remainingEntryOrdinal /= componentEntryCount;
302
303 value *= tensorComponents_[r](componentEntryOrdinal);
304 }
305
306 return value;
307 }
308
314 template <typename iType0, ordinal_type numTensorComponents>
315 KOKKOS_INLINE_FUNCTION typename std::enable_if<std::is_integral<iType0>::value, Scalar>::type
316 operator()(const Kokkos::Array<iType0,numTensorComponents>& entryComponents) const {
317#ifdef HAVE_INTREPID2_DEBUG
318 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(rank_ != 1, std::invalid_argument, "This method is only valid for rank 1 containers.");
319 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(numTensorComponents_ != numTensorComponents, std::invalid_argument, "Tensorial component count mismatch");
320#endif
321 Scalar value = 1.0;
322 for (ordinal_type r=0; r<numTensorComponents; r++)
323 {
324 value *= tensorComponents_[r](entryComponents[r]);
325 }
326 return value;
327 }
328
339 template <typename iType0, typename iType1>
340 KOKKOS_INLINE_FUNCTION typename std::enable_if<
341 (std::is_integral<iType0>::value && std::is_integral<iType1>::value),
342 Scalar>::type
343 operator()(const iType0& tensorEntryIndex0, const iType1& tensorEntryIndex1) const {
344#ifdef HAVE_INTREPID2_DEBUG
345 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(rank_ != 2, std::invalid_argument, "This method is only valid for rank 2 containers.");
346#endif
347
348 if (numTensorComponents_ == 1)
349 {
350 return tensorComponents_[0](tensorEntryIndex0,tensorEntryIndex1);
351 }
352
353 if (!separateFirstComponent_)
354 {
355 Scalar value = 1.0;
356 iType0 remainingEntryOrdinal0 = tensorEntryIndex0;
357 iType1 remainingEntryOrdinal1 = tensorEntryIndex1;
358 for (ordinal_type r=0; r<numTensorComponents_; r++)
359 {
360 auto & component = tensorComponents_[r];
361 const ordinal_type componentEntryCount0 = component.extent_int(0);
362 const ordinal_type componentEntryCount1 = component.extent_int(1);
363 const iType0 componentEntryOrdinal0 = remainingEntryOrdinal0 % componentEntryCount0;
364 const iType1 componentEntryOrdinal1 = remainingEntryOrdinal1 % componentEntryCount1;
365 remainingEntryOrdinal0 /= componentEntryCount0;
366 remainingEntryOrdinal1 /= componentEntryCount1;
367
368 const ordinal_type componentRank = component.rank();
369
370 if (componentRank == 2)
371 {
372 value *= component(componentEntryOrdinal0,componentEntryOrdinal1);
373 }
374 else if (componentRank == 1)
375 {
376 value *= component(componentEntryOrdinal0);
377 }
378 else
379 {
380 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::logic_error, "unsupported component rank encountered");
381 }
382 }
383
384 return value;
385 }
386 else
387 {
388 Scalar value = tensorComponents_[0](tensorEntryIndex0);
389 iType0 remainingEntryOrdinal = tensorEntryIndex1;
390 for (ordinal_type r=1; r<numTensorComponents_; r++)
391 {
392 const ordinal_type componentEntryCount = tensorComponents_[r].extent_int(0);
393 const ordinal_type componentEntryOrdinal = remainingEntryOrdinal % componentEntryCount;
394 remainingEntryOrdinal /= componentEntryCount;
395
396 value *= tensorComponents_[r](componentEntryOrdinal);
397 }
398 return value;
399 }
400 }
401
403 KOKKOS_INLINE_FUNCTION
404 ordinal_type getTensorComponentIndex(const ordinal_type &tensorComponent, const ordinal_type &dim, const ordinal_type &enumerationIndex) const
405 {
406 ordinal_type remainingEntryOrdinal = enumerationIndex;
407 for (ordinal_type r=0; r<tensorComponent; r++)
408 {
409 const auto & component = tensorComponents_[r];
410 const ordinal_type & componentEntryCount = component.extent_int(dim);
411
412 remainingEntryOrdinal /= componentEntryCount;
413 }
414 return remainingEntryOrdinal % tensorComponents_[tensorComponent].extent_int(dim);
415 }
416
426 template <typename iType0, typename iType1, typename iType2>
427 KOKKOS_INLINE_FUNCTION typename std::enable_if<
428 (std::is_integral<iType0>::value && std::is_integral<iType1>::value && std::is_integral<iType2>::value),
429 Scalar>::type
430 operator()(const iType0& tensorEntryIndex0, const iType1& tensorEntryIndex1, const iType2& tensorEntryIndex2) const
431 {
432#ifdef HAVE_INTREPID2_DEBUG
433 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(rank_ != 3, std::invalid_argument, "This method is only valid for rank 3 containers.");
434 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(separateFirstComponent_, std::logic_error, "This method should never be called when separateFirstComponent is true");
435#endif
436
437 Scalar value = 1.0;
438 Kokkos::Array<ordinal_type,3> remainingEntryOrdinal {tensorEntryIndex0, tensorEntryIndex1, tensorEntryIndex2};
439 for (ordinal_type r=0; r<numTensorComponents_; r++)
440 {
441 auto & component = tensorComponents_[r];
442 const ordinal_type componentEntryCount0 = component.extent_int(0);
443 const ordinal_type componentEntryCount1 = component.extent_int(1);
444 const ordinal_type componentEntryCount2 = component.extent_int(2);
445 const ordinal_type componentEntryOrdinal0 = remainingEntryOrdinal[0] % componentEntryCount0;
446 const ordinal_type componentEntryOrdinal1 = remainingEntryOrdinal[1] % componentEntryCount1;
447 const ordinal_type componentEntryOrdinal2 = remainingEntryOrdinal[2] % componentEntryCount2;
448 remainingEntryOrdinal[0] /= componentEntryCount0;
449 remainingEntryOrdinal[1] /= componentEntryCount1;
450 remainingEntryOrdinal[2] /= componentEntryCount2;
451
452 const ordinal_type componentRank = component.rank();
453
454 if (componentRank == 3)
455 {
456 value *= component(componentEntryOrdinal0,componentEntryOrdinal1,componentEntryOrdinal2);
457 }
458 else if (componentRank == 2)
459 {
460 value *= component(componentEntryOrdinal0,componentEntryOrdinal1);
461 }
462 else if (componentRank == 1)
463 {
464 value *= component(componentEntryOrdinal0);
465 }
466 else
467 {
468 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::logic_error, "unsupported component rank encountered");
469 }
470 }
471 return value;
472 }
473
481 template <typename iType0, typename iType1, ordinal_type numTensorComponents>
482 KOKKOS_INLINE_FUNCTION typename std::enable_if<
483 (std::is_integral<iType0>::value && std::is_integral<iType1>::value),
484 Scalar>::type
485 operator()(const Kokkos::Array<iType0,numTensorComponents>& entryComponents0, const Kokkos::Array<iType1,numTensorComponents>& entryComponents1) const {
486#ifdef HAVE_INTREPID2_DEBUG
487 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(numTensorComponents_ != numTensorComponents, std::invalid_argument, "Tensorial component count mismatch");
488#endif
489 Scalar value = 1.0;
490 for (ordinal_type r=0; r<numTensorComponents; r++)
491 {
492 auto & component = tensorComponents_[r];
493 const ordinal_type componentRank = component.rank();
494 if (componentRank == 2)
495 {
496 value *= component(entryComponents0[r],entryComponents1[r]);
497 }
498 else if (componentRank == 1)
499 {
500 value *= component(entryComponents0[r]);
501 }
502 else
503 {
504 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::logic_error, "unsupported component rank encountered");
505 }
506 }
507 return value;
508 }
509
515 template <typename iType>
516 KOKKOS_INLINE_FUNCTION
517 typename std::enable_if<std::is_integral<iType>::value, ordinal_type>::type
518 extent_int(const iType& d) const {
519 return extents_[d];
520 }
521
527 template <typename iType>
528 KOKKOS_INLINE_FUNCTION constexpr
529 typename std::enable_if<std::is_integral<iType>::value, size_t>::type
530 extent(const iType& d) const {
531 return extents_[d];
532 }
533
535 KOKKOS_INLINE_FUNCTION constexpr bool isValid() const
536 {
537 return extents_[0] > 0;
538 }
539
541 KOKKOS_INLINE_FUNCTION
542 ordinal_type rank() const
543 {
544 return rank_;
545 }
546
548 KOKKOS_INLINE_FUNCTION
549 ordinal_type numTensorComponents() const
550 {
551 return numTensorComponents_;
552 }
553
555 KOKKOS_INLINE_FUNCTION
557 {
558 return separateFirstComponent_;
559 }
560
562 void setFirstComponentExtentInDimension0(const ordinal_type &newExtent)
563 {
564 INTREPID2_TEST_FOR_EXCEPTION(!separateFirstComponent_ && (numTensorComponents_ != 1), std::invalid_argument, "setFirstComponentExtent() is only allowed when separateFirstComponent_ is true, or there is only one component");
565 tensorComponents_[0].setExtent(0,newExtent);
566 extents_[0] = newExtent;
567 }
568 };
569}
570
571#endif /* Intrepid2_TensorData_h */
Defines the Data class, a wrapper around a Kokkos::View that allows data that is constant or repeatin...
Contains definitions of custom data types in Intrepid2.
#define INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(test, x, msg)
Wrapper around a Kokkos::View that allows data that is constant or repeating in various logical dimen...
static constexpr ordinal_type MaxTensorComponents
Maximum number of tensor/Cartesian products that can be taken: this allows hypercube basis in 7D to b...
View-like interface to tensor data; tensor components are stored separately and multiplied together a...
void initialize()
Initialize members based on constructor parameters.
TensorData()
Default constructor.
KOKKOS_INLINE_FUNCTION bool separateFirstComponent() const
Returns true if the first component is indexed separately; false if not.
KOKKOS_INLINE_FUNCTION ordinal_type getTensorComponentIndex(const ordinal_type &tensorComponent, const ordinal_type &dim, const ordinal_type &enumerationIndex) const
return the index into the specified tensorial component in the dimension specified corresponding to t...
KOKKOS_INLINE_FUNCTION const Data< Scalar, DeviceType > & getTensorComponent(const ordinal_type &r) const
Returns the requested tensor component.
KOKKOS_INLINE_FUNCTION std::enable_if< std::is_integral< iType0 >::value, Scalar >::type operator()(const Kokkos::Array< iType0, numTensorComponents > &entryComponents) const
Accessor that accepts a fixed-length array with entries corresponding to component indices.
TensorData(const TensorData< Scalar, OtherDeviceType > &tensorData)
Copy-like constructor for differing execution spaces. This performs a deep copy of the underlying dat...
KOKKOS_INLINE_FUNCTION std::enable_if< std::is_integral< iType0 >::value, Scalar >::type operator()(const iType0 &tensorEntryIndex) const
Accessor for rank-1 objects.
TensorData(Kokkos::Array< Data< Scalar, DeviceType >, numTensorComponents > tensorComponents, bool separateFirstComponent=false)
Constructor with fixed-length Kokkos::Array argument.
KOKKOS_INLINE_FUNCTION constexpr bool isValid() const
Returns true for containers that have data; false for those that don't (e.g., those that have been co...
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< std::is_integral< iType >::value, size_t >::type extent(const iType &d) const
Returns the logical extent in the requested dimension.
KOKKOS_INLINE_FUNCTION ordinal_type rank() const
Returns the rank of the container.
KOKKOS_INLINE_FUNCTION std::enable_if< std::is_integral< iType >::value, ordinal_type >::type extent_int(const iType &d) const
Returns the logical extent in the requested dimension.
void setFirstComponentExtentInDimension0(const ordinal_type &newExtent)
Sets the extent of the first component. Only valid when either there is only one component,...
TensorData(Data< Scalar, DeviceType > tensorComponent)
Simple constructor for the case of trivial tensor-product structure (single component)
TensorData(std::vector< Data< Scalar, DeviceType > > tensorComponents, bool separateFirstComponent=false)
Constructor with variable-length std::vector containing the components.
KOKKOS_INLINE_FUNCTION ordinal_type numTensorComponents() const
Return the number of tensorial components.
TensorData(TensorData otherTensorData, std::vector< int > whichComps)
Constructor that takes a subset of the tensorial components of another TensorData container.
TensorData(const TensorData &first, const TensorData &second, bool separateFirstComponent=false)
Constructor to combine two other TensorData objects.