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
53namespace Intrepid2 {
64 template<class Scalar, typename DeviceType>
66 {
67 public:
68 using VectorArray = Kokkos::Array< TensorData<Scalar,DeviceType>, Parameters::MaxVectorComponents >; // 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::MaxVectorComponents> dimToComponent_;
76 Kokkos::Array<int, Parameters::MaxVectorComponents> dimToComponentDim_;
77 Kokkos::Array<int, Parameters::MaxVectorComponents> numDimsForComponent_;
78
79 Kokkos::Array<int,Parameters::MaxTensorComponents> familyFieldUpperBound_; // one higher than the end of family indicated
80
81 unsigned numFamilies_; // number of valid entries in vectorComponents_
82 unsigned numComponents_; // number of valid entries in each entry of vectorComponents_
83 unsigned numPoints_; // the second dimension of each (valid) TensorData entry
84
89 {
90 int lastFieldUpperBound = 0;
91 int numPoints = 0;
92 axialComponents_ = true; // will set to false if we find any valid entries that are not on the "diagonal" (like position for family/component)
93 for (unsigned i=0; i<numFamilies_; i++)
94 {
95 bool validEntryFoundForFamily = false;
96 int numFieldsInFamily = 0;
97 for (unsigned j=0; j<numComponents_; j++)
98 {
99 if (vectorComponents_[i][j].isValid())
100 {
101 if (!validEntryFoundForFamily)
102 {
103 numFieldsInFamily = vectorComponents_[i][j].extent_int(0); // (F,P[,D])
104 validEntryFoundForFamily = true;
105 }
106 else
107 {
108 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");
109 }
110 if (numPoints == 0)
111 {
112 numPoints = vectorComponents_[i][j].extent_int(1); // (F,P[,D])
113 }
114 else
115 {
116 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");
117 }
118 if (i != j)
119 {
120 // valid entry found that is not on the "diagonal": axialComponents is false
121 axialComponents_ = false;
122 }
123 }
124 }
125 lastFieldUpperBound += numFieldsInFamily;
126 familyFieldUpperBound_[i] = lastFieldUpperBound;
127 INTREPID2_TEST_FOR_EXCEPTION(!validEntryFoundForFamily, std::invalid_argument, "Each family must have at least one valid TensorData entry");
128 }
129
130 // do a pass through components to determine total component dim (totalDimension_) and size lookups appropriately
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 }
145 else
146 {
147 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");
148 }
149 }
150 }
151 if (!validEntryFoundForComponent)
152 {
153 // assume that the component takes up exactly one space dim
155 }
156
157 numDimsForComponent_[j] = numDimsForComponent;
158
159 currentDim += numDimsForComponent;
160 }
161 totalDimension_ = currentDim;
162
163 currentDim = 0;
164 for (unsigned j=0; j<numComponents_; j++)
165 {
166 int numDimsForComponent = numDimsForComponent_[j];
167 for (int dim=0; dim<numDimsForComponent; dim++)
168 {
169 dimToComponent_[currentDim+dim] = j;
170 dimToComponentDim_[currentDim+dim] = dim;
171 }
172 currentDim += numDimsForComponent;
173 }
174 numPoints_ = numPoints;
175 }
176 public:
183 template<size_t numFamilies, size_t numComponents>
184 VectorData(Kokkos::Array< Kokkos::Array<TensorData<Scalar,DeviceType>, numComponents>, numFamilies> vectorComponents)
185 :
186 numFamilies_(numFamilies),
187 numComponents_(numComponents)
188 {
189 static_assert(numFamilies <= Parameters::MaxTensorComponents, "numFamilies must be less than Parameters::MaxTensorComponents");
190 static_assert(numComponents <= Parameters::MaxVectorComponents, "numComponents must be less than Parameters::MaxVectorComponents");
191 for (unsigned i=0; i<numFamilies; i++)
192 {
193 for (unsigned j=0; j<numComponents; j++)
194 {
195 vectorComponents_[i][j] = vectorComponents[i][j];
196 }
197 }
198 initialize();
199 }
200
207 VectorData(const std::vector< std::vector<TensorData<Scalar,DeviceType> > > &vectorComponents)
208 {
209 numFamilies_ = vectorComponents.size();
210 INTREPID2_TEST_FOR_EXCEPTION(numFamilies_ <= 0, std::invalid_argument, "numFamilies must be at least 1");
211 numComponents_ = vectorComponents[0].size();
212 for (unsigned i=1; i<numFamilies_; i++)
213 {
214 INTREPID2_TEST_FOR_EXCEPTION(vectorComponents[i].size() != numComponents_, std::invalid_argument, "each family must have the same number of components");
215 }
216
217 INTREPID2_TEST_FOR_EXCEPTION(numFamilies_ > Parameters::MaxTensorComponents, std::invalid_argument, "numFamilies must be at most Parameters::MaxTensorComponents");
218 INTREPID2_TEST_FOR_EXCEPTION(numComponents_ > Parameters::MaxVectorComponents, std::invalid_argument, "numComponents must be at most Parameters::MaxVectorComponents");
219 for (unsigned i=0; i<numFamilies_; i++)
220 {
221 for (unsigned j=0; j<numComponents_; j++)
222 {
223 vectorComponents_[i][j] = vectorComponents[i][j];
224 }
225 }
226 initialize();
227 }
228
236 template<size_t numComponents>
238 {
239 if (axialComponents)
240 {
241 numFamilies_ = numComponents;
242 numComponents_ = numComponents;
243 for (unsigned d=0; d<numComponents_; d++)
244 {
245 vectorComponents_[d][d] = vectorComponents[d];
246 }
247 }
248 else
249 {
250 numFamilies_ = 1;
251 numComponents_ = numComponents;
252 for (unsigned d=0; d<numComponents_; d++)
253 {
254 vectorComponents_[0][d] = vectorComponents[d];
255 }
256 }
257 initialize();
258 }
259
267 VectorData(std::vector< TensorData<Scalar,DeviceType> > vectorComponents, bool axialComponents)
268 : numComponents_(vectorComponents.size())
269 {
270 if (axialComponents)
271 {
272 numFamilies_ = numComponents_;
273 for (unsigned d=0; d<numComponents_; d++)
274 {
275 vectorComponents_[d][d] = vectorComponents[d];
276 }
277 }
278 else
279 {
280 numFamilies_ = 1;
281 for (unsigned d=0; d<numComponents_; d++)
282 {
283 vectorComponents_[0][d] = vectorComponents[d];
284 }
285 }
286 initialize();
287 }
288
290 template<typename OtherDeviceType, class = typename std::enable_if< std::is_same<typename DeviceType::memory_space, typename OtherDeviceType::memory_space>::value>::type,
291 class = typename std::enable_if<!std::is_same<DeviceType,OtherDeviceType>::value>::type>
292 VectorData(const VectorData<Scalar,OtherDeviceType> &vectorData)
293 :
294 numFamilies_(vectorData.numFamilies()),
295 numComponents_(vectorData.numComponents())
296 {
297 if (vectorData.isValid())
298 {
299 for (unsigned i=0; i<numFamilies_; i++)
300 {
301 for (unsigned j=0; j<numComponents_; j++)
302 {
303 vectorComponents_[i][j] = vectorData.getComponent(i, j);
304 }
305 }
306 initialize();
307 }
308 }
309
311 template<typename OtherDeviceType, class = typename std::enable_if<!std::is_same<typename DeviceType::memory_space, typename OtherDeviceType::memory_space>::value>::type>
313 :
314 numFamilies_(vectorData.numFamilies()),
315 numComponents_(vectorData.numComponents())
316 {
317 if (vectorData.isValid())
318 {
319 for (unsigned i=0; i<numFamilies_; i++)
320 {
321 for (unsigned j=0; j<numComponents_; j++)
322 {
323 vectorComponents_[i][j] = vectorData.getComponent(i, j);
324 }
325 }
326 initialize();
327 }
328 }
329
332 :
333 VectorData(Kokkos::Array< TensorData<Scalar,DeviceType>, 1>(data), true)
334 {}
335
338 :
339 VectorData(Kokkos::Array< TensorData<Scalar,DeviceType>, 1>({TensorData<Scalar,DeviceType>(data)}), true)
340 {}
341
344 :
345 numFamilies_(0), numComponents_(0)
346 {}
347
349 KOKKOS_INLINE_FUNCTION
350 bool axialComponents() const
351 {
352 return axialComponents_;
353 }
354
356 KOKKOS_INLINE_FUNCTION
357 int numDimsForComponent(int &componentOrdinal) const
358 {
359 return numDimsForComponent_[componentOrdinal];
360 }
361
363 KOKKOS_INLINE_FUNCTION
364 int numFields() const
365 {
366 return familyFieldUpperBound_[numFamilies_-1];
367 }
368
370 KOKKOS_INLINE_FUNCTION
371 int familyFieldOrdinalOffset(const int &familyOrdinal) const
372 {
373 return (familyOrdinal == 0) ? 0 : familyFieldUpperBound_[familyOrdinal-1];
374 }
375
377 KOKKOS_INLINE_FUNCTION
378 int numPoints() const
379 {
380 return numPoints_;
381 }
382
384 KOKKOS_INLINE_FUNCTION
385 int spaceDim() const
386 {
387 return totalDimension_;
388 }
389
391 KOKKOS_INLINE_FUNCTION
392 Scalar operator()(const int &fieldOrdinal, const int &pointOrdinal, const int &dim) const
393 {
394 int fieldOrdinalInFamily = fieldOrdinal;
395 int familyForField = 0;
396 if (numFamilies_ > 1)
397 {
398 familyForField = -1;
399 int previousFamilyEnd = -1;
400 int fieldAdjustment = 0;
401 // this loop is written in such a way as to avoid branching for CUDA performance
402 for (unsigned family=0; family<numFamilies_; family++)
403 {
404 const bool fieldInRange = (fieldOrdinal > previousFamilyEnd) && (fieldOrdinal < familyFieldUpperBound_[family]);
405 familyForField = fieldInRange ? family : familyForField;
406 fieldAdjustment = fieldInRange ? previousFamilyEnd + 1 : fieldAdjustment;
407 previousFamilyEnd = familyFieldUpperBound_[family] - 1;
408 }
409#ifdef HAVE_INTREPID2_DEBUG
410 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(familyForField == -1, std::invalid_argument, "family for field not found");
411#endif
412
413 fieldOrdinalInFamily = fieldOrdinal - fieldAdjustment;
414 }
415
416 const int componentOrdinal = dimToComponent_[dim];
417
418 const auto &component = vectorComponents_[familyForField][componentOrdinal];
419 if (component.isValid())
420 {
421 const int componentRank = component.rank();
422 if (componentRank == 2) // (F,P) container
423 {
424 return component(fieldOrdinalInFamily,pointOrdinal);
425 }
426 else if (componentRank == 3) // (F,P,D)
427 {
428 return component(fieldOrdinalInFamily,pointOrdinal,dimToComponentDim_[dim]);
429 }
430 else
431 {
432 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::logic_error, "Unsupported component rank");
433 return -1; // unreachable, but compilers complain otherwise...
434 }
435 }
436 else // invalid component: placeholder means 0
437 {
438 return 0;
439 }
440 }
441
446 KOKKOS_INLINE_FUNCTION
447 const TensorData<Scalar,DeviceType> & getComponent(const int &componentOrdinal) const
448 {
449 if (axialComponents_)
450 {
451 return vectorComponents_[componentOrdinal][componentOrdinal];
452 }
453 else if (numFamilies_ == 1)
454 {
455 return vectorComponents_[0][componentOrdinal];
456 }
457 else
458 {
459 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "Ambiguous component request; use the two-argument getComponent()");
460 }
461 // nvcc warns here about a missing return.
462 return vectorComponents_[6][6]; // likely this is an empty container, but anyway it's an unreachable line...
463 }
464
470 KOKKOS_INLINE_FUNCTION
471 const TensorData<Scalar,DeviceType> & getComponent(const int &familyOrdinal, const int &componentOrdinal) const
472 {
473 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(familyOrdinal < 0, std::invalid_argument, "familyOrdinal must be non-negative");
474 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(static_cast<unsigned>(familyOrdinal) >= numFamilies_, std::invalid_argument, "familyOrdinal out of bounds");
475 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(componentOrdinal < 0, std::invalid_argument, "componentOrdinal must be non-negative");
476 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(static_cast<unsigned>(componentOrdinal) >= numComponents_, std::invalid_argument, "componentOrdinal out of bounds");
477
478 return vectorComponents_[familyOrdinal][componentOrdinal];
479 }
480
482 KOKKOS_INLINE_FUNCTION
483 int extent_int(const int &r) const
484 {
485 // logically (F,P,D) container
486 if (r == 0) return numFields();
487 else if (r == 1) return numPoints();
488 else if (r == 2) return totalDimension_;
489 else if (r > 2) return 1;
490
491 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "Unsupported rank");
492 return -1; // unreachable; return here included to avoid compiler warnings.
493 }
494
496 KOKKOS_INLINE_FUNCTION
497 unsigned rank() const
498 {
499 // logically (F,P,D) container
500 return 3;
501 }
502
504 KOKKOS_INLINE_FUNCTION int numComponents() const
505 {
506 return numComponents_;
507 }
508
510 KOKKOS_INLINE_FUNCTION int numFamilies() const
511 {
512 return numFamilies_;
513 }
514
516 KOKKOS_INLINE_FUNCTION int familyForFieldOrdinal(const int &fieldOrdinal) const
517 {
518 int matchingFamily = -1;
519 int fieldsSoFar = 0;
520 // logic here is a little bit more complex to avoid branch divergence
521 for (int i=0; i<numFamilies_; i++)
522 {
523 const bool fieldIsBeyondPreviousFamily = (fieldOrdinal >= fieldsSoFar);
524 fieldsSoFar += numFieldsInFamily(i);
525 const bool fieldIsBeforeCurrentFamily = (fieldOrdinal < fieldsSoFar);
526 const bool fieldMatchesFamily = fieldIsBeyondPreviousFamily && fieldIsBeforeCurrentFamily;
527 matchingFamily = fieldMatchesFamily ? i : matchingFamily;
528 }
529 return matchingFamily;
530 }
531
533 KOKKOS_INLINE_FUNCTION int numFieldsInFamily(const unsigned &familyOrdinal) const
534 {
535 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(familyOrdinal >= numFamilies_, std::invalid_argument, "familyOrdinal out of bounds");
536 int numFields = -1;
537 for (unsigned componentOrdinal=0; componentOrdinal<numComponents_; componentOrdinal++)
538 {
539 numFields = vectorComponents_[familyOrdinal][componentOrdinal].isValid() ? vectorComponents_[familyOrdinal][componentOrdinal].extent_int(0) : numFields;
540 }
541 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(numFields < 0, std::logic_error, "numFields was not properly initialized");
542 return numFields;
543 }
544
546 KOKKOS_INLINE_FUNCTION constexpr bool isValid() const
547 {
548 return numComponents_ > 0;
549 }
550 };
551}
552
553#endif /* Intrepid2_VectorData_h */
#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...
static constexpr ordinal_type MaxVectorComponents
Maximum number of components that a VectorData object will store – 66 corresponds to OPERATOR_D10 on ...
View-like interface to tensor data; tensor components are stored separately and multiplied together a...
Reference-space field values for a basis, designed to support typical vector-valued bases.
KOKKOS_INLINE_FUNCTION int numPoints() const
Returns the number of points; corresponds to the second dimension of this container.
VectorData(const VectorData< Scalar, OtherDeviceType > &vectorData)
copy-like constructor for differing execution spaces. This does a deep copy of underlying views.
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 bool axialComponents() const
Returns true only if the families are so structured that the first family has nonzeros only in the x ...
VectorData()
default constructor; results in an invalid container.
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
VectorData(std::vector< TensorData< Scalar, DeviceType > > vectorComponents, bool axialComponents)
Simplified constructor for gradients of HGRAD, and values of HDIV and HCURL vector bases.
VectorData(const std::vector< std::vector< TensorData< Scalar, DeviceType > > > &vectorComponents)
Standard constructor for the arbitrary case, accepting a variable-length std::vector argument.
KOKKOS_INLINE_FUNCTION int familyForFieldOrdinal(const int &fieldOrdinal) const
Returns the family ordinal corresponding to the indicated field ordinal.
KOKKOS_INLINE_FUNCTION unsigned rank() const
Returns the rank of this container, which is 3.
KOKKOS_INLINE_FUNCTION const TensorData< Scalar, DeviceType > & getComponent(const int &familyOrdinal, const int &componentOrdinal) const
General component accessor.
VectorData(Data< Scalar, DeviceType > data)
Simple 1-argument constructor for the case of trivial tensor product structure. The argument should h...
VectorData(TensorData< Scalar, DeviceType > data)
Simple 1-argument constructor for the case of trivial tensor product structure. The argument should h...
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).
KOKKOS_INLINE_FUNCTION int numComponents() const
returns the number of components
KOKKOS_INLINE_FUNCTION int familyFieldOrdinalOffset(const int &familyOrdinal) const
Returns the field ordinal offset for the specified family.
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 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 int numFamilies() const
returns the number of families
KOKKOS_INLINE_FUNCTION int spaceDim() const
Returns the spatial dimension; corresponds to the third dimension of this container.
KOKKOS_INLINE_FUNCTION int numDimsForComponent(int &componentOrdinal) const
Returns the number of dimensions corresponding to the specified component.
KOKKOS_INLINE_FUNCTION int numFields() const
Returns the total number of fields; corresponds to the first dimension of this container.
KOKKOS_INLINE_FUNCTION int extent_int(const int &r) const
Returns the extent in the specified dimension as an int.
VectorData(Kokkos::Array< TensorData< Scalar, DeviceType >, numComponents > vectorComponents, bool axialComponents)
Simplified constructor for gradients of HGRAD, and values of HDIV and HCURL vector bases.