50 #ifndef Intrepid2_VectorData_h 51 #define Intrepid2_VectorData_h 64 template<
class Scalar,
typename DeviceType>
69 using FamilyVectorArray = Kokkos::Array< VectorArray, Parameters::MaxTensorComponents>;
71 FamilyVectorArray vectorComponents_;
72 bool axialComponents_;
75 Kokkos::Array<int,Parameters::MaxTensorComponents> dimToComponent_;
76 Kokkos::Array<int,Parameters::MaxTensorComponents> dimToComponentDim_;
78 Kokkos::Array<int,Parameters::MaxTensorComponents> numDimsForComponent_;
80 Kokkos::Array<int,Parameters::MaxTensorComponents> familyFieldUpperBound_;
82 unsigned numFamilies_;
83 unsigned numComponents_;
91 int lastFieldUpperBound = 0;
93 axialComponents_ =
true;
94 for (
unsigned i=0; i<numFamilies_; i++)
96 bool validEntryFoundForFamily =
false;
98 for (
unsigned j=0; j<numComponents_; j++)
100 if (vectorComponents_[i][j].
isValid())
102 if (!validEntryFoundForFamily)
105 validEntryFoundForFamily =
true;
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");
113 numPoints = vectorComponents_[i][j].extent_int(1);
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");
122 axialComponents_ =
false;
127 familyFieldUpperBound_[i] = lastFieldUpperBound;
128 INTREPID2_TEST_FOR_EXCEPTION(!validEntryFoundForFamily, std::invalid_argument,
"Each family must have at least one valid TensorData entry");
132 for (
unsigned j=0; j<numComponents_; j++)
134 bool validEntryFoundForComponent =
false;
136 for (
unsigned i=0; i<numFamilies_; i++)
138 if (vectorComponents_[i][j].
isValid())
140 if (!validEntryFoundForComponent)
142 validEntryFoundForComponent =
true;
146 dimToComponent_[currentDim+dim] = j;
147 dimToComponentDim_[currentDim+dim] = dim;
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");
156 if (!validEntryFoundForComponent)
160 dimToComponent_[currentDim] = j;
161 dimToComponentDim_[currentDim] = 0;
169 totalDimension_ = currentDim;
178 template<
size_t numFamilies,
size_t numComponents>
190 vectorComponents_[i][j] = vectorComponents[i][j];
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++)
209 INTREPID2_TEST_FOR_EXCEPTION(vectorComponents[i].size() != numComponents_, std::invalid_argument,
"each family must have the same number of components");
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++)
216 for (
unsigned j=0; j<numComponents_; j++)
218 vectorComponents_[i][j] = vectorComponents[i][j];
231 template<
size_t numComponents>
238 for (
unsigned d=0; d<numComponents_; d++)
240 vectorComponents_[d][d] = vectorComponents[d];
247 for (
unsigned d=0; d<numComponents_; d++)
249 vectorComponents_[0][d] = vectorComponents[d];
263 : numComponents_(vectorComponents.size())
267 numFamilies_ = numComponents_;
268 for (
unsigned d=0; d<numComponents_; d++)
270 vectorComponents_[d][d] = vectorComponents[d];
276 for (
unsigned d=0; d<numComponents_; d++)
278 vectorComponents_[0][d] = vectorComponents[d];
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)
289 numFamilies_(vectorData.numFamilies()),
290 numComponents_(vectorData.numComponents())
292 if (vectorData.isValid())
294 for (
unsigned i=0; i<numFamilies_; i++)
296 for (
unsigned j=0; j<numComponents_; j++)
298 vectorComponents_[i][j] = vectorData.
getComponent(i, j);
306 template<typename OtherDeviceType, class = typename std::enable_if<!std::is_same<typename DeviceType::memory_space, typename OtherDeviceType::memory_space>::value>::type>
314 for (
unsigned i=0; i<numFamilies_; i++)
316 for (
unsigned j=0; j<numComponents_; j++)
318 vectorComponents_[i][j] = vectorData.
getComponent(i, j);
340 numFamilies_(0), numComponents_(0)
344 KOKKOS_INLINE_FUNCTION
347 return axialComponents_;
351 KOKKOS_INLINE_FUNCTION
354 return numDimsForComponent_[componentOrdinal];
358 KOKKOS_INLINE_FUNCTION
361 return familyFieldUpperBound_[numFamilies_-1];
365 KOKKOS_INLINE_FUNCTION
368 return (familyOrdinal == 0) ? 0 : familyFieldUpperBound_[familyOrdinal-1];
372 KOKKOS_INLINE_FUNCTION
379 KOKKOS_INLINE_FUNCTION
382 return totalDimension_;
386 KOKKOS_INLINE_FUNCTION
387 Scalar
operator()(
const int &fieldOrdinal,
const int &pointOrdinal,
const int &dim)
const 389 int fieldOrdinalInFamily = fieldOrdinal;
390 int familyForField = 0;
391 if (numFamilies_ > 1)
394 int previousFamilyEnd = -1;
395 int fieldAdjustment = 0;
397 for (
unsigned family=0; family<numFamilies_; family++)
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;
404 #ifdef HAVE_INTREPID2_DEBUG 405 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(familyForField == -1, std::invalid_argument,
"family for field not found");
408 fieldOrdinalInFamily = fieldOrdinal - fieldAdjustment;
411 const int componentOrdinal = dimToComponent_[dim];
413 const auto &component = vectorComponents_[familyForField][componentOrdinal];
414 if (component.isValid())
416 const int componentRank = component.rank();
417 if (componentRank == 2)
419 return component(fieldOrdinalInFamily,pointOrdinal);
421 else if (componentRank == 3)
423 return component(fieldOrdinalInFamily,pointOrdinal,dimToComponentDim_[dim]);
427 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(
true, std::logic_error,
"Unsupported component rank");
441 KOKKOS_INLINE_FUNCTION
444 if (axialComponents_)
446 return vectorComponents_[componentOrdinal][componentOrdinal];
448 else if (numFamilies_ == 1)
450 return vectorComponents_[0][componentOrdinal];
454 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(
true, std::invalid_argument,
"Ambiguous component request; use the two-argument getComponent()");
457 return vectorComponents_[6][6];
465 KOKKOS_INLINE_FUNCTION
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");
473 return vectorComponents_[familyOrdinal][componentOrdinal];
477 KOKKOS_INLINE_FUNCTION
483 else if (r == 2)
return totalDimension_;
484 else if (r > 2)
return 1;
486 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(
true, std::invalid_argument,
"Unsupported rank");
491 KOKKOS_INLINE_FUNCTION
501 return numComponents_;
513 int matchingFamily = -1;
516 for (
int i=0; i<numFamilies_; i++)
518 const bool fieldIsBeyondPreviousFamily = (fieldOrdinal >= fieldsSoFar);
520 const bool fieldIsBeforeCurrentFamily = (fieldOrdinal < fieldsSoFar);
521 const bool fieldMatchesFamily = fieldIsBeyondPreviousFamily && fieldIsBeforeCurrentFamily;
522 matchingFamily = fieldMatchesFamily ? i : matchingFamily;
524 return matchingFamily;
530 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(familyOrdinal >= numFamilies_, std::invalid_argument,
"familyOrdinal out of bounds");
532 for (
unsigned componentOrdinal=0; componentOrdinal<numComponents_; componentOrdinal++)
534 numFields = vectorComponents_[familyOrdinal][componentOrdinal].isValid() ? vectorComponents_[familyOrdinal][componentOrdinal].extent_int(0) :
numFields;
536 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(
numFields < 0, std::logic_error,
"numFields was not properly initialized");
541 KOKKOS_INLINE_FUNCTION constexpr
bool isValid()
const 543 return numComponents_ > 0;
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'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...