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 
58 namespace 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 
76  void initialize()
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 
167  :
168  TensorData(Kokkos::Array< Data<Scalar,DeviceType>, 1>({tensorComponent}), false)
169  {}
170 
177  :
178  extents_({0,0,0,0,0,0,0}),
179  rank_(0)
180  {}
181 
183  template<typename OtherDeviceType, class = typename std::enable_if< std::is_same<typename DeviceType::memory_space, typename OtherDeviceType::memory_space>::value>::type,
184  class = typename std::enable_if<!std::is_same<DeviceType,OtherDeviceType>::value>::type>
185  TensorData(const TensorData<Scalar,OtherDeviceType> &tensorData)
186  {
187  if (tensorData.isValid())
188  {
189  numTensorComponents_ = tensorData.numTensorComponents();
190  for (ordinal_type r=0; r<numTensorComponents_; r++)
191  {
192  Data<Scalar,OtherDeviceType> otherTensorComponent = tensorData.getTensorComponent(r);
193  tensorComponents_[r] = Data<Scalar,DeviceType>(otherTensorComponent);
194  }
195  initialize();
196  }
197  else
198  {
199  extents_ = Kokkos::Array<ordinal_type,7>{0,0,0,0,0,0,0};
200  rank_ = 0;
201  }
202  }
203 
207  template<typename OtherDeviceType, class = typename std::enable_if<!std::is_same<typename DeviceType::memory_space, typename OtherDeviceType::memory_space>::value>::type>
209  {
210  if (tensorData.isValid())
211  {
212  numTensorComponents_ = tensorData.numTensorComponents();
213  for (ordinal_type r=0; r<numTensorComponents_; r++)
214  {
215  Data<Scalar,OtherDeviceType> otherTensorComponent = tensorData.getTensorComponent(r);
216  tensorComponents_[r] = Data<Scalar,DeviceType>(otherTensorComponent);
217  }
218  initialize();
219  }
220  else
221  {
222  extents_ = Kokkos::Array<ordinal_type,7>{0,0,0,0,0,0,0};
223  rank_ = 0;
224  }
225  }
226 
232  KOKKOS_INLINE_FUNCTION
233  const Data<Scalar,DeviceType> & getTensorComponent(const ordinal_type &r) const
234  {
235  return tensorComponents_[r];
236  }
237 
243  template <typename iType0>
244  KOKKOS_INLINE_FUNCTION typename std::enable_if<std::is_integral<iType0>::value, Scalar>::type
245  operator()(const iType0& tensorEntryIndex) const {
246 #ifdef HAVE_INTREPID2_DEBUG
247  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(rank_ != 1, std::invalid_argument, "This method is only valid for rank 1 containers.");
248 #endif
249  Scalar value = 1.0;
250  iType0 remainingEntryOrdinal = tensorEntryIndex;
251  for (ordinal_type r=0; r<numTensorComponents_; r++)
252  {
253  const ordinal_type componentEntryCount = tensorComponents_[r].extent_int(0);
254  const ordinal_type componentEntryOrdinal = remainingEntryOrdinal % componentEntryCount;
255  remainingEntryOrdinal /= componentEntryCount;
256 
257  value *= tensorComponents_[r](componentEntryOrdinal);
258  }
259 
260  return value;
261  }
262 
268  template <typename iType0, ordinal_type numTensorComponents>
269  KOKKOS_INLINE_FUNCTION typename std::enable_if<std::is_integral<iType0>::value, Scalar>::type
270  operator()(const Kokkos::Array<iType0,numTensorComponents>& entryComponents) const {
271 #ifdef HAVE_INTREPID2_DEBUG
272  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(rank_ != 1, std::invalid_argument, "This method is only valid for rank 1 containers.");
273  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(numTensorComponents_ != numTensorComponents, std::invalid_argument, "Tensorial component count mismatch");
274 #endif
275  Scalar value = 1.0;
276  for (ordinal_type r=0; r<numTensorComponents; r++)
277  {
278  value *= tensorComponents_[r](entryComponents[r]);
279  }
280  return value;
281  }
282 
293  template <typename iType0, typename iType1>
294  KOKKOS_INLINE_FUNCTION typename std::enable_if<
295  (std::is_integral<iType0>::value && std::is_integral<iType1>::value),
296  Scalar>::type
297  operator()(const iType0& tensorEntryIndex0, const iType1& tensorEntryIndex1) const {
298 #ifdef HAVE_INTREPID2_DEBUG
299  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(rank_ != 2, std::invalid_argument, "This method is only valid for rank 2 containers.");
300 #endif
301 
302  if (numTensorComponents_ == 1)
303  {
304  return tensorComponents_[0](tensorEntryIndex0,tensorEntryIndex1);
305  }
306 
307  if (!separateFirstComponent_)
308  {
309  Scalar value = 1.0;
310  iType0 remainingEntryOrdinal0 = tensorEntryIndex0;
311  iType1 remainingEntryOrdinal1 = tensorEntryIndex1;
312  for (ordinal_type r=0; r<numTensorComponents_; r++)
313  {
314  auto & component = tensorComponents_[r];
315  const ordinal_type componentEntryCount0 = component.extent_int(0);
316  const ordinal_type componentEntryCount1 = component.extent_int(1);
317  const iType0 componentEntryOrdinal0 = remainingEntryOrdinal0 % componentEntryCount0;
318  const iType1 componentEntryOrdinal1 = remainingEntryOrdinal1 % componentEntryCount1;
319  remainingEntryOrdinal0 /= componentEntryCount0;
320  remainingEntryOrdinal1 /= componentEntryCount1;
321 
322  const ordinal_type componentRank = component.rank();
323 
324  if (componentRank == 2)
325  {
326  value *= component(componentEntryOrdinal0,componentEntryOrdinal1);
327  }
328  else if (componentRank == 1)
329  {
330  value *= component(componentEntryOrdinal0);
331  }
332  else
333  {
334  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::logic_error, "unsupported component rank encountered");
335  }
336  }
337 
338  return value;
339  }
340  else
341  {
342  Scalar value = tensorComponents_[0](tensorEntryIndex0);
343  iType0 remainingEntryOrdinal = tensorEntryIndex1;
344  for (ordinal_type r=1; r<numTensorComponents_; r++)
345  {
346  const ordinal_type componentEntryCount = tensorComponents_[r].extent_int(0);
347  const ordinal_type componentEntryOrdinal = remainingEntryOrdinal % componentEntryCount;
348  remainingEntryOrdinal /= componentEntryCount;
349 
350  value *= tensorComponents_[r](componentEntryOrdinal);
351  }
352  return value;
353  }
354  }
355 
357  KOKKOS_INLINE_FUNCTION
358  ordinal_type getTensorComponentIndex(const ordinal_type &tensorComponent, const ordinal_type &dim, const ordinal_type &enumerationIndex) const
359  {
360  ordinal_type remainingEntryOrdinal = enumerationIndex;
361  for (ordinal_type r=0; r<tensorComponent; r++)
362  {
363  const auto & component = tensorComponents_[r];
364  const ordinal_type & componentEntryCount = component.extent_int(dim);
365 
366  remainingEntryOrdinal /= componentEntryCount;
367  }
368  return remainingEntryOrdinal % tensorComponents_[tensorComponent].extent_int(dim);
369  }
370 
380  template <typename iType0, typename iType1, typename iType2>
381  KOKKOS_INLINE_FUNCTION typename std::enable_if<
382  (std::is_integral<iType0>::value && std::is_integral<iType1>::value && std::is_integral<iType2>::value),
383  Scalar>::type
384  operator()(const iType0& tensorEntryIndex0, const iType1& tensorEntryIndex1, const iType2& tensorEntryIndex2) const
385  {
386 #ifdef HAVE_INTREPID2_DEBUG
387  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(rank_ != 3, std::invalid_argument, "This method is only valid for rank 3 containers.");
388  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(separateFirstComponent_, std::logic_error, "This method should never be called when separateFirstComponent is true");
389 #endif
390 
391  Scalar value = 1.0;
392  Kokkos::Array<ordinal_type,3> remainingEntryOrdinal {tensorEntryIndex0, tensorEntryIndex1, tensorEntryIndex2};
393  for (ordinal_type r=0; r<numTensorComponents_; r++)
394  {
395  auto & component = tensorComponents_[r];
396  const ordinal_type componentEntryCount0 = component.extent_int(0);
397  const ordinal_type componentEntryCount1 = component.extent_int(1);
398  const ordinal_type componentEntryCount2 = component.extent_int(2);
399  const ordinal_type componentEntryOrdinal0 = remainingEntryOrdinal[0] % componentEntryCount0;
400  const ordinal_type componentEntryOrdinal1 = remainingEntryOrdinal[1] % componentEntryCount1;
401  const ordinal_type componentEntryOrdinal2 = remainingEntryOrdinal[2] % componentEntryCount2;
402  remainingEntryOrdinal[0] /= componentEntryCount0;
403  remainingEntryOrdinal[1] /= componentEntryCount1;
404  remainingEntryOrdinal[2] /= componentEntryCount2;
405 
406  const ordinal_type componentRank = component.rank();
407 
408  if (componentRank == 3)
409  {
410  value *= component(componentEntryOrdinal0,componentEntryOrdinal1,componentEntryOrdinal2);
411  }
412  else if (componentRank == 2)
413  {
414  value *= component(componentEntryOrdinal0,componentEntryOrdinal1);
415  }
416  else if (componentRank == 1)
417  {
418  value *= component(componentEntryOrdinal0);
419  }
420  else
421  {
422  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::logic_error, "unsupported component rank encountered");
423  }
424  }
425  return value;
426  }
427 
435  template <typename iType0, typename iType1, ordinal_type numTensorComponents>
436  KOKKOS_INLINE_FUNCTION typename std::enable_if<
437  (std::is_integral<iType0>::value && std::is_integral<iType1>::value),
438  Scalar>::type
439  operator()(const Kokkos::Array<iType0,numTensorComponents>& entryComponents0, const Kokkos::Array<iType1,numTensorComponents>& entryComponents1) const {
440 #ifdef HAVE_INTREPID2_DEBUG
441  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(numTensorComponents_ != numTensorComponents, std::invalid_argument, "Tensorial component count mismatch");
442 #endif
443  Scalar value = 1.0;
444  for (ordinal_type r=0; r<numTensorComponents; r++)
445  {
446  auto & component = tensorComponents_[r];
447  const ordinal_type componentRank = component.rank();
448  if (componentRank == 2)
449  {
450  value *= component(entryComponents0[r],entryComponents1[r]);
451  }
452  else if (componentRank == 1)
453  {
454  value *= component(entryComponents0[r]);
455  }
456  else
457  {
458  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::logic_error, "unsupported component rank encountered");
459  }
460  }
461  return value;
462  }
463 
469  template <typename iType>
470  KOKKOS_INLINE_FUNCTION
471  typename std::enable_if<std::is_integral<iType>::value, ordinal_type>::type
472  extent_int(const iType& d) const {
473  return extents_[d];
474  }
475 
481  template <typename iType>
482  KOKKOS_INLINE_FUNCTION constexpr
483  typename std::enable_if<std::is_integral<iType>::value, size_t>::type
484  extent(const iType& d) const {
485  return extents_[d];
486  }
487 
489  KOKKOS_INLINE_FUNCTION constexpr bool isValid() const
490  {
491  return extents_[0] > 0;
492  }
493 
495  KOKKOS_INLINE_FUNCTION
496  ordinal_type rank() const
497  {
498  return rank_;
499  }
500 
502  KOKKOS_INLINE_FUNCTION
503  ordinal_type numTensorComponents() const
504  {
505  return numTensorComponents_;
506  }
507 
509  KOKKOS_INLINE_FUNCTION
511  {
512  return separateFirstComponent_;
513  }
514 
516  void setFirstComponentExtentInDimension0(const ordinal_type &newExtent)
517  {
518  INTREPID2_TEST_FOR_EXCEPTION(!separateFirstComponent_ && (numTensorComponents_ != 1), std::invalid_argument, "setFirstComponentExtent() is only allowed when separateFirstComponent_ is true, or there is only one component");
519  tensorComponents_[0].setExtent(0,newExtent);
520  }
521  };
522 }
523 
524 #endif /* Intrepid2_TensorData_h */
TensorData(Kokkos::Array< Data< Scalar, DeviceType >, numTensorComponents > tensorComponents, bool separateFirstComponent=false)
Constructor with fixed-length Kokkos::Array argument.
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...
Defines the Data class, a wrapper around a Kokkos::View that allows data that is constant or repeatin...
KOKKOS_INLINE_FUNCTION constexpr bool isValid() const
Returns true for containers that have data; false for those that don&#39;t (e.g., those that have been co...
void initialize()
Initialize members based on constructor parameters.
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 ordinal_type rank() const
Returns the rank of the container.
TensorData(std::vector< Data< Scalar, DeviceType > > tensorComponents, bool separateFirstComponent=false)
Constructor with variable-length std::vector containing the components.
KOKKOS_INLINE_FUNCTION bool separateFirstComponent() const
Returns true if the first component is indexed separately; false if not.
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.
KOKKOS_INLINE_FUNCTION const Data< Scalar, DeviceType > & getTensorComponent(const ordinal_type &r) const
Returns the requested tensor component.
Contains definitions of custom data types in Intrepid2.
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()
Default constructor.
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 numTensorComponents() const
Return the number of tensorial components.
void setFirstComponentExtentInDimension0(const ordinal_type &newExtent)
Sets the extent of the first component. Only valid when either there is only one component, or when separateFirstComponent() returns true. The intended use case is when the 0 dimension in first component represents a cell index, and the container is resized to match a workset size that does not evenly divide the number of cells.
KOKKOS_INLINE_FUNCTION std::enable_if< std::is_integral< iType0 >::value, Scalar >::type operator()(const iType0 &tensorEntryIndex) const
Accessor for rank-1 objects.
View-like interface to tensor data; tensor components are stored separately and multiplied together a...
TensorData(Data< Scalar, DeviceType > tensorComponent)
Simple constructor for the case of trivial tensor-product structure (single component) ...
static constexpr ordinal_type MaxTensorComponents
Maximum number of tensor/Cartesian products that can be taken: this allows hypercube basis in 7D to b...