Intrepid2
Intrepid2_VectorData.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_VectorData_h
51 #define Intrepid2_VectorData_h
52 
53 namespace Intrepid2 {
64  template<class Scalar, typename DeviceType>
65  class VectorData
66  {
67  public:
68  using VectorArray = Kokkos::Array< TensorData<Scalar,DeviceType>, Parameters::MaxTensorComponents>; // for axis-aligned case, these correspond entry-wise to the axis with which the vector values align
69  using FamilyVectorArray = Kokkos::Array< VectorArray, Parameters::MaxTensorComponents>;
70 
71  FamilyVectorArray vectorComponents_; // outer: family ordinal; inner: component/spatial dimension ordinal
72  bool axialComponents_; // if true, each entry in vectorComponents_ is an axial component vector; for 3D: (f1,0,0); (0,f2,0); (0,0,f3). The 0s are represented by trivial/invalid TensorData objects. In this case, numComponents_ == numFamilies_.
73 
74  int totalDimension_;
75  Kokkos::Array<int,Parameters::MaxTensorComponents> dimToComponent_;
76  Kokkos::Array<int,Parameters::MaxTensorComponents> dimToComponentDim_;
77 
78  Kokkos::Array<int,Parameters::MaxTensorComponents> numDimsForComponent_;
79 
80  Kokkos::Array<int,Parameters::MaxTensorComponents> familyFieldUpperBound_; // one higher than the end of family indicated
81 
82  unsigned numFamilies_; // number of valid entries in vectorComponents_
83  unsigned numComponents_; // number of valid entries in each entry of vectorComponents_
84  unsigned numPoints_; // the second dimension of each (valid) TensorData entry
85 
89  void initialize()
90  {
91  int lastFieldUpperBound = 0;
92  int numPoints = 0;
93  axialComponents_ = true; // will set to false if we find any valid entries that are not on the "diagonal" (like position for family/component)
94  for (unsigned i=0; i<numFamilies_; i++)
95  {
96  bool validEntryFoundForFamily = false;
97  int numFieldsInFamily = 0;
98  for (unsigned j=0; j<numComponents_; j++)
99  {
100  if (vectorComponents_[i][j].isValid())
101  {
102  if (!validEntryFoundForFamily)
103  {
104  numFieldsInFamily = vectorComponents_[i][j].extent_int(0); // (F,P[,D])
105  validEntryFoundForFamily = true;
106  }
107  else
108  {
109  INTREPID2_TEST_FOR_EXCEPTION(numFieldsInFamily != vectorComponents_[i][j].extent_int(0), std::invalid_argument, "Each valid TensorData entry within a family must agree with the others on the number of fields in the family");
110  }
111  if (numPoints == 0)
112  {
113  numPoints = vectorComponents_[i][j].extent_int(1); // (F,P[,D])
114  }
115  else
116  {
117  INTREPID2_TEST_FOR_EXCEPTION(numPoints != vectorComponents_[i][j].extent_int(1), std::invalid_argument, "Each valid TensorData entry must agree with the others on the number of points");
118  }
119  if (i != j)
120  {
121  // valid entry found that is not on the "diagonal": axialComponents is false
122  axialComponents_ = false;
123  }
124  }
125  }
126  lastFieldUpperBound += numFieldsInFamily;
127  familyFieldUpperBound_[i] = lastFieldUpperBound;
128  INTREPID2_TEST_FOR_EXCEPTION(!validEntryFoundForFamily, std::invalid_argument, "Each family must have at least one valid TensorData entry");
129  }
130 
131  int currentDim = 0;
132  for (unsigned j=0; j<numComponents_; j++)
133  {
134  bool validEntryFoundForComponent = false;
135  int numDimsForComponent = 0;
136  for (unsigned i=0; i<numFamilies_; i++)
137  {
138  if (vectorComponents_[i][j].isValid())
139  {
140  if (!validEntryFoundForComponent)
141  {
142  validEntryFoundForComponent = true;
143  numDimsForComponent = vectorComponents_[i][j].extent_int(2); // (F,P,D) container or (F,P) container
144  for (int dim=0; dim<numDimsForComponent; dim++)
145  {
146  dimToComponent_[currentDim+dim] = j;
147  dimToComponentDim_[currentDim+dim] = dim;
148  }
149  }
150  else
151  {
152  INTREPID2_TEST_FOR_EXCEPTION(numDimsForComponent != vectorComponents_[i][j].extent_int(2), std::invalid_argument, "Components in like positions must agree across families on the number of dimensions spanned by that component position");
153  }
154  }
155  }
156  if (!validEntryFoundForComponent)
157  {
158  // assume that the component takes up exactly one space dim
160  dimToComponent_[currentDim] = j;
161  dimToComponentDim_[currentDim] = 0;
162  }
163 
164  numDimsForComponent_[j] = numDimsForComponent;
165 
166  currentDim += numDimsForComponent;
167  }
168  numPoints_ = numPoints;
169  totalDimension_ = currentDim;
170  }
171  public:
178  template<size_t numFamilies, size_t numComponents>
179  VectorData(Kokkos::Array< Kokkos::Array<TensorData<Scalar,DeviceType>, numComponents>, numFamilies> vectorComponents)
180  :
181  numFamilies_(numFamilies),
182  numComponents_(numComponents)
183  {
184  static_assert(numFamilies <= Parameters::MaxTensorComponents, "numFamilies must be less than Parameters::MaxTensorComponents");
185  static_assert(numComponents <= Parameters::MaxTensorComponents, "numComponents must be less than Parameters::MaxTensorComponents");
186  for (unsigned i=0; i<numFamilies; i++)
187  {
188  for (unsigned j=0; j<numComponents; j++)
189  {
190  vectorComponents_[i][j] = vectorComponents[i][j];
191  }
192  }
193  initialize();
194  }
195 
202  VectorData(const std::vector< std::vector<TensorData<Scalar,DeviceType> > > &vectorComponents)
203  {
204  numFamilies_ = vectorComponents.size();
205  INTREPID2_TEST_FOR_EXCEPTION(numFamilies_ <= 0, std::invalid_argument, "numFamilies must be at least 1");
206  numComponents_ = vectorComponents[0].size();
207  for (unsigned i=1; i<numFamilies_; i++)
208  {
209  INTREPID2_TEST_FOR_EXCEPTION(vectorComponents[i].size() != numComponents_, std::invalid_argument, "each family must have the same number of components");
210  }
211 
212  INTREPID2_TEST_FOR_EXCEPTION(numFamilies_ > Parameters::MaxTensorComponents, std::invalid_argument, "numFamilies must be less than Parameters::MaxTensorComponents");
213  INTREPID2_TEST_FOR_EXCEPTION(numComponents_ > Parameters::MaxTensorComponents, std::invalid_argument, "numComponents must be less than Parameters::MaxTensorComponents");
214  for (unsigned i=0; i<numFamilies_; i++)
215  {
216  for (unsigned j=0; j<numComponents_; j++)
217  {
218  vectorComponents_[i][j] = vectorComponents[i][j];
219  }
220  }
221  initialize();
222  }
223 
231  template<size_t numComponents>
233  {
234  if (axialComponents)
235  {
236  numFamilies_ = numComponents;
237  numComponents_ = numComponents;
238  for (unsigned d=0; d<numComponents_; d++)
239  {
240  vectorComponents_[d][d] = vectorComponents[d];
241  }
242  }
243  else
244  {
245  numFamilies_ = 1;
246  numComponents_ = numComponents;
247  for (unsigned d=0; d<numComponents_; d++)
248  {
249  vectorComponents_[0][d] = vectorComponents[d];
250  }
251  }
252  initialize();
253  }
254 
262  VectorData(std::vector< TensorData<Scalar,DeviceType> > vectorComponents, bool axialComponents)
263  : numComponents_(vectorComponents.size())
264  {
265  if (axialComponents)
266  {
267  numFamilies_ = numComponents_;
268  for (unsigned d=0; d<numComponents_; d++)
269  {
270  vectorComponents_[d][d] = vectorComponents[d];
271  }
272  }
273  else
274  {
275  numFamilies_ = 1;
276  for (unsigned d=0; d<numComponents_; d++)
277  {
278  vectorComponents_[0][d] = vectorComponents[d];
279  }
280  }
281  initialize();
282  }
283 
285  template<typename OtherDeviceType, class = typename std::enable_if< std::is_same<typename DeviceType::memory_space, typename OtherDeviceType::memory_space>::value>::type,
286  class = typename std::enable_if<!std::is_same<DeviceType,OtherDeviceType>::value>::type>
287  VectorData(const VectorData<Scalar,OtherDeviceType> &vectorData)
288  :
289  numFamilies_(vectorData.numFamilies()),
290  numComponents_(vectorData.numComponents())
291  {
292  if (vectorData.isValid())
293  {
294  for (unsigned i=0; i<numFamilies_; i++)
295  {
296  for (unsigned j=0; j<numComponents_; j++)
297  {
298  vectorComponents_[i][j] = vectorData.getComponent(i, j);
299  }
300  }
301  initialize();
302  }
303  }
304 
306  template<typename OtherDeviceType, class = typename std::enable_if<!std::is_same<typename DeviceType::memory_space, typename OtherDeviceType::memory_space>::value>::type>
308  :
309  numFamilies_(vectorData.numFamilies()),
310  numComponents_(vectorData.numComponents())
311  {
312  if (vectorData.isValid())
313  {
314  for (unsigned i=0; i<numFamilies_; i++)
315  {
316  for (unsigned j=0; j<numComponents_; j++)
317  {
318  vectorComponents_[i][j] = vectorData.getComponent(i, j);
319  }
320  }
321  initialize();
322  }
323  }
324 
327  :
328  VectorData(Kokkos::Array< TensorData<Scalar,DeviceType>, 1>(data), true)
329  {}
330 
333  :
334  VectorData(Kokkos::Array< TensorData<Scalar,DeviceType>, 1>({TensorData<Scalar,DeviceType>(data)}), true)
335  {}
336 
339  :
340  numFamilies_(0), numComponents_(0)
341  {}
342 
344  KOKKOS_INLINE_FUNCTION
345  bool axialComponents() const
346  {
347  return axialComponents_;
348  }
349 
351  KOKKOS_INLINE_FUNCTION
352  int numDimsForComponent(int &componentOrdinal) const
353  {
354  return numDimsForComponent_[componentOrdinal];
355  }
356 
358  KOKKOS_INLINE_FUNCTION
359  int numFields() const
360  {
361  return familyFieldUpperBound_[numFamilies_-1];
362  }
363 
365  KOKKOS_INLINE_FUNCTION
366  int familyFieldOrdinalOffset(const int &familyOrdinal) const
367  {
368  return (familyOrdinal == 0) ? 0 : familyFieldUpperBound_[familyOrdinal-1];
369  }
370 
372  KOKKOS_INLINE_FUNCTION
373  int numPoints() const
374  {
375  return numPoints_;
376  }
377 
379  KOKKOS_INLINE_FUNCTION
380  int spaceDim() const
381  {
382  return totalDimension_;
383  }
384 
386  KOKKOS_INLINE_FUNCTION
387  Scalar operator()(const int &fieldOrdinal, const int &pointOrdinal, const int &dim) const
388  {
389  int fieldOrdinalInFamily = fieldOrdinal;
390  int familyForField = 0;
391  if (numFamilies_ > 1)
392  {
393  familyForField = -1;
394  int previousFamilyEnd = -1;
395  int fieldAdjustment = 0;
396  // this loop is written in such a way as to avoid branching for CUDA performance
397  for (unsigned family=0; family<numFamilies_; family++)
398  {
399  const bool fieldInRange = (fieldOrdinal > previousFamilyEnd) && (fieldOrdinal < familyFieldUpperBound_[family]);
400  familyForField = fieldInRange ? family : familyForField;
401  fieldAdjustment = fieldInRange ? previousFamilyEnd + 1 : fieldAdjustment;
402  previousFamilyEnd = familyFieldUpperBound_[family] - 1;
403  }
404 #ifdef HAVE_INTREPID2_DEBUG
405  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(familyForField == -1, std::invalid_argument, "family for field not found");
406 #endif
407 
408  fieldOrdinalInFamily = fieldOrdinal - fieldAdjustment;
409  }
410 
411  const int componentOrdinal = dimToComponent_[dim];
412 
413  const auto &component = vectorComponents_[familyForField][componentOrdinal];
414  if (component.isValid())
415  {
416  const int componentRank = component.rank();
417  if (componentRank == 2) // (F,P) container
418  {
419  return component(fieldOrdinalInFamily,pointOrdinal);
420  }
421  else if (componentRank == 3) // (F,P,D)
422  {
423  return component(fieldOrdinalInFamily,pointOrdinal,dimToComponentDim_[dim]);
424  }
425  else
426  {
427  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::logic_error, "Unsupported component rank");
428  return -1; // unreachable, but compilers complain otherwise...
429  }
430  }
431  else // invalid component: placeholder means 0
432  {
433  return 0;
434  }
435  }
436 
441  KOKKOS_INLINE_FUNCTION
442  const TensorData<Scalar,DeviceType> & getComponent(const int &componentOrdinal) const
443  {
444  if (axialComponents_)
445  {
446  return vectorComponents_[componentOrdinal][componentOrdinal];
447  }
448  else if (numFamilies_ == 1)
449  {
450  return vectorComponents_[0][componentOrdinal];
451  }
452  else
453  {
454  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "Ambiguous component request; use the two-argument getComponent()");
455  }
456  // nvcc warns here about a missing return.
457  return vectorComponents_[6][6]; // likely this is an empty container, but anyway it's an unreachable line...
458  }
459 
465  KOKKOS_INLINE_FUNCTION
466  const TensorData<Scalar,DeviceType> & getComponent(const int &familyOrdinal, const int &componentOrdinal) const
467  {
468  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(familyOrdinal < 0, std::invalid_argument, "familyOrdinal must be non-negative");
469  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(static_cast<unsigned>(familyOrdinal) >= numFamilies_, std::invalid_argument, "familyOrdinal out of bounds");
470  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(componentOrdinal < 0, std::invalid_argument, "componentOrdinal must be non-negative");
471  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(static_cast<unsigned>(componentOrdinal) >= numComponents_, std::invalid_argument, "componentOrdinal out of bounds");
472 
473  return vectorComponents_[familyOrdinal][componentOrdinal];
474  }
475 
477  KOKKOS_INLINE_FUNCTION
478  int extent_int(const int &r) const
479  {
480  // logically (F,P,D) container
481  if (r == 0) return numFields();
482  else if (r == 1) return numPoints();
483  else if (r == 2) return totalDimension_;
484  else if (r > 2) return 1;
485 
486  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "Unsupported rank");
487  return -1; // unreachable; return here included to avoid compiler warnings.
488  }
489 
491  KOKKOS_INLINE_FUNCTION
492  unsigned rank() const
493  {
494  // logically (F,P,D) container
495  return 3;
496  }
497 
499  KOKKOS_INLINE_FUNCTION int numComponents() const
500  {
501  return numComponents_;
502  }
503 
505  KOKKOS_INLINE_FUNCTION int numFamilies() const
506  {
507  return numFamilies_;
508  }
509 
511  KOKKOS_INLINE_FUNCTION int familyForFieldOrdinal(const int &fieldOrdinal) const
512  {
513  int matchingFamily = -1;
514  int fieldsSoFar = 0;
515  // logic here is a little bit more complex to avoid branch divergence
516  for (int i=0; i<numFamilies_; i++)
517  {
518  const bool fieldIsBeyondPreviousFamily = (fieldOrdinal >= fieldsSoFar);
519  fieldsSoFar += numFieldsInFamily(i);
520  const bool fieldIsBeforeCurrentFamily = (fieldOrdinal < fieldsSoFar);
521  const bool fieldMatchesFamily = fieldIsBeyondPreviousFamily && fieldIsBeforeCurrentFamily;
522  matchingFamily = fieldMatchesFamily ? i : matchingFamily;
523  }
524  return matchingFamily;
525  }
526 
528  KOKKOS_INLINE_FUNCTION int numFieldsInFamily(const unsigned &familyOrdinal) const
529  {
530  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(familyOrdinal >= numFamilies_, std::invalid_argument, "familyOrdinal out of bounds");
531  int numFields = -1;
532  for (unsigned componentOrdinal=0; componentOrdinal<numComponents_; componentOrdinal++)
533  {
534  numFields = vectorComponents_[familyOrdinal][componentOrdinal].isValid() ? vectorComponents_[familyOrdinal][componentOrdinal].extent_int(0) : numFields;
535  }
536  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(numFields < 0, std::logic_error, "numFields was not properly initialized");
537  return numFields;
538  }
539 
541  KOKKOS_INLINE_FUNCTION constexpr bool isValid() const
542  {
543  return numComponents_ > 0;
544  }
545  };
546 }
547 
548 #endif /* Intrepid2_VectorData_h */
KOKKOS_INLINE_FUNCTION int spaceDim() const
Returns the spatial dimension; corresponds to the third dimension of this container.
VectorData(Kokkos::Array< Kokkos::Array< TensorData< Scalar, DeviceType >, numComponents >, numFamilies > vectorComponents)
Standard constructor for the arbitrary case, accepting a fixed-length array argument.
KOKKOS_INLINE_FUNCTION int numPoints() const
Returns the number of points; corresponds to the second dimension of this container.
KOKKOS_INLINE_FUNCTION int extent_int(const int &r) const
Returns the extent in the specified dimension as an int.
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...
VectorData(const std::vector< std::vector< TensorData< Scalar, DeviceType > > > &vectorComponents)
Standard constructor for the arbitrary case, accepting a variable-length std::vector argument...
VectorData(Kokkos::Array< TensorData< Scalar, DeviceType >, numComponents > vectorComponents, bool axialComponents)
Simplified constructor for gradients of HGRAD, and values of HDIV and HCURL vector bases...
KOKKOS_INLINE_FUNCTION const TensorData< Scalar, DeviceType > & getComponent(const int &familyOrdinal, const int &componentOrdinal) const
General component accessor.
void initialize()
Initialize members based on constructor parameters; all constructors should call this after populatin...
KOKKOS_INLINE_FUNCTION int numFieldsInFamily(const unsigned &familyOrdinal) const
returns the number of fields in the specified family
KOKKOS_INLINE_FUNCTION Scalar operator()(const int &fieldOrdinal, const int &pointOrdinal, const int &dim) const
Accessor for the container, which has shape (F,P,D).
VectorData(Data< Scalar, DeviceType > data)
Simple 1-argument constructor for the case of trivial tensor product structure. The argument should h...
VectorData()
default constructor; results in an invalid container.
KOKKOS_INLINE_FUNCTION int numComponents() const
returns the number of components
KOKKOS_INLINE_FUNCTION unsigned rank() const
Returns the rank of this container, which is 3.
VectorData(const VectorData< Scalar, OtherDeviceType > &vectorData)
copy-like constructor for differing execution spaces. This does a deep copy of underlying views...
KOKKOS_INLINE_FUNCTION int familyFieldOrdinalOffset(const int &familyOrdinal) const
Returns the field ordinal offset for the specified family.
KOKKOS_INLINE_FUNCTION int numFields() const
Returns the total number of fields; corresponds to the first dimension of this container.
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...
KOKKOS_INLINE_FUNCTION int numDimsForComponent(int &componentOrdinal) const
Returns the number of dimensions corresponding to the specified component.
KOKKOS_INLINE_FUNCTION int numFamilies() const
returns the number of families
VectorData(TensorData< Scalar, DeviceType > data)
Simple 1-argument constructor for the case of trivial tensor product structure. The argument should h...
KOKKOS_INLINE_FUNCTION bool axialComponents() const
Returns true only if the families are so structured that the first family has nonzeros only in the x ...
VectorData(std::vector< TensorData< Scalar, DeviceType > > vectorComponents, bool axialComponents)
Simplified constructor for gradients of HGRAD, and values of HDIV and HCURL vector bases...
View-like interface to tensor data; tensor components are stored separately and multiplied together a...
KOKKOS_INLINE_FUNCTION int familyForFieldOrdinal(const int &fieldOrdinal) const
Returns the family ordinal corresponding to the indicated field ordinal.
Reference-space field values for a basis, designed to support typical vector-valued bases...
static constexpr ordinal_type MaxTensorComponents
Maximum number of tensor/Cartesian products that can be taken: this allows hypercube basis in 7D to b...