Intrepid2
Intrepid2_TensorPoints.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_TensorPoints_h
51#define Intrepid2_TensorPoints_h
52
53#include <Kokkos_Vector.hpp>
54
55namespace Intrepid2 {
59 template<class PointScalar, typename DeviceType>
61 protected:
62 Kokkos::Array< ScalarView<PointScalar,DeviceType>, Parameters::MaxTensorComponents> pointTensorComponents_; // each component has dimensions (P,D)
63 ordinal_type numTensorComponents_;
64 ordinal_type totalPointCount_;
65 ordinal_type totalDimension_;
66 Kokkos::View<ordinal_type*, DeviceType> dimToComponent_;
67 Kokkos::View<ordinal_type*, DeviceType> dimToComponentDim_;
68 Kokkos::Array<ordinal_type, Parameters::MaxTensorComponents> pointModulus_;
69 Kokkos::Array<ordinal_type, Parameters::MaxTensorComponents> pointDivisor_;
70
71 bool isValid_;
72 using reference_type = typename ScalarView<PointScalar,DeviceType>::reference_type;
73
74 void TEST_VALID_POINT_COMPONENTS()
75 {
76#ifdef HAVE_INTREPID2_DEBUG
77 if (isValid_)
78 {
79 for (ordinal_type r=0; r<numTensorComponents_; r++)
80 {
81 const auto & pointComponent = pointTensorComponents_[r];
82 INTREPID2_TEST_FOR_EXCEPTION(2 != pointComponent.rank(), std::invalid_argument, "Each component must have shape (P,D)");
83 INTREPID2_TEST_FOR_EXCEPTION(pointComponent.extent_int(0) <= 0, std::invalid_argument, "Each component must have at least one point");
84 }
85 }
86#endif
87 }
88
93 {
94 TEST_VALID_POINT_COMPONENTS();
95
96 totalPointCount_ = 1;
97 totalDimension_ = 0;
98 for (ordinal_type r=0; r<numTensorComponents_; r++)
99 {
100 totalPointCount_ *= pointTensorComponents_[r].extent_int(0);
101 totalDimension_ += pointTensorComponents_[r].extent_int(1);
102 }
103 ordinal_type pointDivisor = 1;
104 for (ordinal_type r=0; r<numTensorComponents_; r++)
105 {
106 pointModulus_[r] = pointTensorComponents_[r].extent_int(0);
107 pointDivisor_[r] = pointDivisor;
108 pointDivisor *= pointTensorComponents_[r].extent_int(0);
109 }
110 dimToComponent_ = Kokkos::View<ordinal_type*, DeviceType>("dimToComponent_",totalDimension_);
111 dimToComponentDim_ = Kokkos::View<ordinal_type*, DeviceType>("dimToComponentDim_",totalDimension_);
112 ordinal_type d=0;
113 ordinal_type dimsSoFar = 0;
114
115 auto dimToComponentHost = Kokkos::create_mirror_view(dimToComponent_);
116 auto dimToComponentDimHost = Kokkos::create_mirror_view(dimToComponentDim_);
117 for (ordinal_type r=0; r<numTensorComponents_; r++)
118 {
119 const int componentDim = pointTensorComponents_[r].extent_int(1);
120 for (int i=0; i<componentDim; i++)
121 {
122 dimToComponentHost[d] = r;
123 dimToComponentDimHost[d] = d - dimsSoFar;
124 d++;
125 }
126 dimsSoFar += componentDim;
127 }
128 Kokkos::deep_copy(dimToComponent_,dimToComponentHost);
129 Kokkos::deep_copy(dimToComponentDim_,dimToComponentDimHost);
130 }
131 public:
138 template<size_t numTensorComponents>
139 TensorPoints(Kokkos::Array< ScalarView<PointScalar,DeviceType>, numTensorComponents> pointTensorComponents)
140 :
141 numTensorComponents_(numTensorComponents),
142 isValid_(true)
143 {
144 for (unsigned r=0; r<numTensorComponents; r++)
145 {
146 pointTensorComponents_[r] = pointTensorComponents[r];
147 }
148
149 initialize();
150 }
151
159 TensorPoints(TensorPoints otherPointsContainer, std::vector<int> whichDims)
160 :
161 numTensorComponents_(whichDims.size()),
162 isValid_(true)
163 {
164 int r = 0;
165 for (const auto & componentOrdinal : whichDims)
166 {
167 pointTensorComponents_[r++] = otherPointsContainer.getTensorComponent(componentOrdinal);
168 }
169
170 initialize();
171 }
172
179 TensorPoints(std::vector< ScalarView<PointScalar,DeviceType>> pointTensorComponents)
180 :
181 numTensorComponents_(pointTensorComponents.size()),
182 isValid_(true)
183 {
184 for (ordinal_type r=0; r<numTensorComponents_; r++)
185 {
186 pointTensorComponents_[r] = pointTensorComponents[r];
187 }
188
189 initialize();
190 }
191
198 TensorPoints(ScalarView<PointScalar,DeviceType> points)
199 :
200 numTensorComponents_(1),
201 isValid_(true)
202 {
203 pointTensorComponents_[0] = points;
204 initialize();
205 }
206
212 template<class OtherPointsContainer>
213 void copyPointsContainer(ScalarView<PointScalar,DeviceType> toPoints, OtherPointsContainer fromPoints)
214 {
215 const int numPoints = fromPoints.extent_int(0);
216 const int numDims = fromPoints.extent_int(1);
217 using ExecutionSpace = typename DeviceType::execution_space;
218 auto policy = Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<2>>({0,0},{numPoints,numDims});
219 Kokkos::parallel_for("copy points", policy,
220 KOKKOS_LAMBDA (const int &i0, const int &i1) {
221 toPoints(i0,i1) = fromPoints(i0,i1);
222 });
223 }
224
226 template<typename OtherDeviceType, class = typename std::enable_if< std::is_same<typename DeviceType::memory_space, typename OtherDeviceType::memory_space>::value>::type,
227 class = typename std::enable_if<!std::is_same<DeviceType,OtherDeviceType>::value>::type>
228 TensorPoints(const TensorPoints<PointScalar,OtherDeviceType> &tensorPoints)
229 :
230 numTensorComponents_(tensorPoints.numTensorComponents()),
231 isValid_(tensorPoints.isValid())
232 {
233 if (isValid_)
234 {
235 for (ordinal_type r=0; r<numTensorComponents_; r++)
236 {
237 pointTensorComponents_[r] = tensorPoints.getTensorComponent(r);
238 }
239 initialize();
240 }
241 }
242
244 template<typename OtherDeviceType, class = typename std::enable_if<!std::is_same<typename DeviceType::memory_space, typename OtherDeviceType::memory_space>::value>::type>
246 :
247 numTensorComponents_(tensorPoints.numTensorComponents()),
248 isValid_(tensorPoints.isValid())
249 {
250 if (isValid_)
251 {
252 for (ordinal_type r=0; r<numTensorComponents_; r++)
253 {
254 ScalarView<PointScalar,OtherDeviceType> otherPointComponent = tensorPoints.getTensorComponent(r);
255 const int numPoints = otherPointComponent.extent_int(0);
256 const int numDims = otherPointComponent.extent_int(1);
257 pointTensorComponents_[r] = ScalarView<PointScalar,DeviceType>("Intrepid2 point component", numPoints, numDims);
258
259 using MemorySpace = typename DeviceType::memory_space;
260 auto pointComponentMirror = Kokkos::create_mirror_view_and_copy(MemorySpace(), otherPointComponent);
261
262 copyPointsContainer(pointTensorComponents_[r], pointComponentMirror);
263 }
264 initialize();
265 }
266 }
267
270 isValid_(false)
271 // when constructed with default arguments, TensorPoints should not be used…
272 // default constructor is only provided for things like CellGeometry, which has TensorPoints as a member,
273 // but only uses it in certain modes.
274 {}
275
277 ordinal_type componentPointCount(const ordinal_type &tensorComponentOrdinal) const
278 {
279 return pointTensorComponents_[tensorComponentOrdinal].extent_int(0);
280 }
281
290 template <typename iType0, typename iType1>
291 KOKKOS_INLINE_FUNCTION typename std::enable_if<
292 (std::is_integral<iType0>::value && std::is_integral<iType1>::value),
293 reference_type>::type
294 operator()(const iType0& tensorPointIndex, const iType1& dim) const {
295 const ordinal_type component = dimToComponent_[dim];
296 const ordinal_type d = dimToComponentDim_[dim];
297 const ordinal_type componentPointOrdinal = (tensorPointIndex / pointDivisor_[component]) % pointModulus_[component];
298 return pointTensorComponents_[component](componentPointOrdinal,d);
299 }
300
309 template <typename iType0, typename iType1, size_t numTensorComponents>
310 KOKKOS_INLINE_FUNCTION typename std::enable_if<
311 (std::is_integral<iType0>::value && std::is_integral<iType1>::value),
312 reference_type>::type
313 operator()(const Kokkos::Array<iType0,numTensorComponents>& pointOrdinalComponents, const iType1& dim) const {
314 const ordinal_type component = dimToComponent_[dim];
315 const ordinal_type d = dimToComponentDim_[dim];
316 const ordinal_type componentPointOrdinal = pointOrdinalComponents[component];
317 return pointTensorComponents_[component](componentPointOrdinal,d);
318 }
319
325 template <typename iType>
326 KOKKOS_INLINE_FUNCTION
327 typename std::enable_if<std::is_integral<iType>::value, int>::type
328 extent_int(const iType& r) const {
329 if (r == static_cast<iType>(0))
330 {
331 return totalPointCount_;
332 }
333 else if (r == static_cast<iType>(1))
334 {
335 return totalDimension_;
336 }
337 else
338 {
339 return 1;
340 }
341 }
342
348 template <typename iType>
349 KOKKOS_INLINE_FUNCTION constexpr
350 typename std::enable_if<std::is_integral<iType>::value, size_t>::type
351 extent(const iType& r) const {
352 // C++ prior to 14 doesn't allow constexpr if statements; the compound ternary is here to keep things constexpr
353 return (r == static_cast<iType>(0)) ? totalPointCount_
354 :
355 (r == static_cast<iType>(1)) ? totalDimension_ : 1;
356 }
357
359 ScalarView<PointScalar,DeviceType> allocateAndFillExpandedRawPointView() const
360 {
361 const int numPoints = this->extent_int(0);
362 const int spaceDim = this->extent_int(1);
363 ScalarView<PointScalar,DeviceType> expandedRawPoints("expanded raw points from TensorPoints", numPoints, spaceDim);
364 TensorPoints<PointScalar,DeviceType> tensorPoints(*this); // (shallow) copy for lambda capture
365 using ExecutionSpace = typename DeviceType::execution_space;
366 Kokkos::parallel_for(
367 Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<2>>({0,0},{numPoints,spaceDim}),
368 KOKKOS_LAMBDA (const int &pointOrdinal, const int &d) {
369 expandedRawPoints(pointOrdinal,d) = tensorPoints(pointOrdinal,d);
370 });
371 return expandedRawPoints;
372 }
373
379 KOKKOS_INLINE_FUNCTION
380 ScalarView<PointScalar,DeviceType> getTensorComponent(const ordinal_type &r) const
381 {
382 return pointTensorComponents_[r];
383 }
384
386 KOKKOS_INLINE_FUNCTION
387 bool isValid() const
388 {
389 return isValid_;
390 }
391
393 KOKKOS_INLINE_FUNCTION
394 ordinal_type numTensorComponents() const
395 {
396 return numTensorComponents_;
397 }
398
400 KOKKOS_INLINE_FUNCTION
401 constexpr ordinal_type rank() const
402 {
403 return 2; // shape is (P,D)
404 }
405
406 // NOTE: the extractTensorPoints() code commented out below appears not to work. We don't have tests against it, though.
407 // TODO: either delete this, or re-enable, add tests, and fix.
408// template<class ViewType>
409// static TensorPoints<PointScalar,DeviceType> extractTensorPoints( ViewType expandedPoints, const std::vector<ordinal_type> &dimensionExtents )
410// {
411// const ordinal_type numComponents = dimensionExtents.size();
412// const ordinal_type numPoints = expandedPoints.extent_int(0);
413// Kokkos::Array<ordinal_type,Parameters::MaxTensorComponents> componentPointCounts;
414// Kokkos::Array<ordinal_type,Parameters::MaxTensorComponents> componentPointTensorStride;
415// Kokkos::Array<ordinal_type,Parameters::MaxTensorComponents> componentPointDimOffsets;
416//
417// // for simplicity of implementation, we copy to host:
418// auto hostExpandedPoints = Kokkos::create_mirror_view_and_copy(typename Kokkos::HostSpace::memory_space(), expandedPoints);
419//
420// ordinal_type dimOffset = 0;
421// ordinal_type tensorPointStride = 1; // increases with componentOrdinal
422//
423// TensorPoints<PointScalar,DeviceType> invalidTensorData; // will be returned if extraction does not succeed.
424//
425// for (ordinal_type componentOrdinal=0; componentOrdinal<numComponents; componentOrdinal++)
426// {
427// componentPointDimOffsets[componentOrdinal] = dimOffset;
428// componentPointTensorStride[componentOrdinal] = tensorPointStride;
429// const ordinal_type numDimsForComponent = dimensionExtents[componentOrdinal];
430// std::vector<PointScalar> firstPoint(numDimsForComponent);
431// for (ordinal_type d=0; d<numDimsForComponent; d++)
432// {
433// firstPoint[d] = hostExpandedPoints(0,d+dimOffset);
434// }
435//
436// // we assume that once we see the same point twice, we've found the cycle length:
437// componentPointCounts[componentOrdinal] = -1;
438// for (ordinal_type pointOrdinal=1; pointOrdinal<numPoints; pointOrdinal += tensorPointStride)
439// {
440// bool matches = true;
441// for (ordinal_type d=0; d<numDimsForComponent; d++)
442// {
443// matches = matches && (firstPoint[d] == hostExpandedPoints(pointOrdinal,d+dimOffset));
444// }
445// if (matches)
446// {
447// componentPointCounts[componentOrdinal] = pointOrdinal;
448// break;
449// }
450// }
451// if (componentPointCounts[componentOrdinal] == -1)
452// {
453// // no matches found -> no tensor decomposition available
454// return invalidTensorData;
455// }
456//
457// // check that we got the cycle length correct:
458// for (ordinal_type pointOrdinal=0; pointOrdinal<componentPointCounts[componentOrdinal]; pointOrdinal += tensorPointStride)
459// {
460// std::vector<PointScalar> point(numDimsForComponent);
461// for (ordinal_type d=0; d<numDimsForComponent; d++)
462// {
463// point[d] = hostExpandedPoints(pointOrdinal,d+dimOffset);
464// }
465// // each of the following points should match:
466// for (ordinal_type secondPointOrdinal=0; secondPointOrdinal<numPoints; secondPointOrdinal += tensorPointStride*componentPointCounts[componentOrdinal])
467// {
468// bool matches = true;
469// for (ordinal_type d=0; d<numDimsForComponent; d++)
470// {
471// matches = matches && (point[d] == hostExpandedPoints(secondPointOrdinal,d+dimOffset));
472// }
473// if (!matches)
474// {
475// // fail:
476// return invalidTensorData;
477// }
478// }
479// }
480//
481// dimOffset += numDimsForComponent;
482// tensorPointStride *= componentPointCounts[componentOrdinal];
483// }
484//
485// std::vector< ScalarView<PointScalar,DeviceType> > componentPoints(numComponents);
486// std::vector< ScalarView<PointScalar,Kokkos::HostSpace> > componentPointsHost(numComponents);
487// for (ordinal_type componentOrdinal=0; componentOrdinal<numComponents; componentOrdinal++)
488// {
489// const ordinal_type numPointsForComponent = componentPointCounts[componentOrdinal];
490// const ordinal_type dimForComponent = dimensionExtents[componentOrdinal];
491// componentPoints[componentOrdinal] = ScalarView<PointScalar,DeviceType>("extracted tensor components", numPointsForComponent, dimForComponent);
492// componentPointsHost[componentOrdinal] = Kokkos::create_mirror(componentPoints[componentOrdinal]);
493// }
494//
495// for (ordinal_type componentOrdinal=0; componentOrdinal<numComponents; componentOrdinal++)
496// {
497// const ordinal_type numComponentPoints = componentPointCounts[componentOrdinal];
498//
499// auto hostView = componentPointsHost[componentOrdinal];
500// auto deviceView = componentPoints[componentOrdinal];
501// const ordinal_type tensorPointStride = componentPointTensorStride[componentOrdinal];
502// const ordinal_type dimOffset = componentPointDimOffsets[componentOrdinal];
503// const ordinal_type numDimsForComponent = dimensionExtents[componentOrdinal];
504// for (ordinal_type componentPointOrdinal=0; componentPointOrdinal<numComponentPoints; componentPointOrdinal++)
505// {
506// const ordinal_type expandedPointOrdinal = componentPointOrdinal*tensorPointStride;
507// for (ordinal_type d=0; d<numDimsForComponent; d++)
508// {
509// hostView(componentPointOrdinal,d) = hostExpandedPoints(expandedPointOrdinal,d+dimOffset);
510// }
511// }
512// Kokkos::deep_copy(deviceView, hostView);
513// }
514//
515// // prior to return, check all points agree in all dimensions with the input points
516// // for convenience, we do this check on host, too
517// TensorPoints<PointScalar,Kokkos::HostSpace> hostTensorPoints(componentPointsHost);
518// const ordinal_type totalDim = expandedPoints.extent_int(1);
519// bool matches = true;
520// for (ordinal_type pointOrdinal=0; pointOrdinal<numPoints; pointOrdinal++)
521// {
522// for (ordinal_type d=0; d<totalDim; d++)
523// {
524// const auto &originalCoord = hostExpandedPoints(pointOrdinal,d);
525// const auto &tensorCoord = hostTensorPoints(pointOrdinal,d);
526// if (originalCoord != tensorCoord)
527// {
528// matches = false;
529// }
530// }
531// }
532// if (!matches)
533// {
534// return invalidTensorData;
535// }
536//
537// return TensorPoints<PointScalar,DeviceType>(componentPoints);
538// }
539 };
540
541}
542
543#endif /* Intrepid2_TensorPoints_h */
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 points; point components are stored separately; the appropriate coordin...
KOKKOS_INLINE_FUNCTION constexpr ordinal_type rank() const
Return the rank of the container, which is 2.
KOKKOS_INLINE_FUNCTION ordinal_type numTensorComponents() const
Returns the number of tensorial components.
TensorPoints(const TensorPoints< PointScalar, OtherDeviceType > &tensorPoints)
copy-like constructor for differing memory spaces. This does a deep_copy of the underlying view.
KOKKOS_INLINE_FUNCTION ScalarView< PointScalar, DeviceType > getTensorComponent(const ordinal_type &r) const
Returns the requested tensor component.
KOKKOS_INLINE_FUNCTION bool isValid() const
Returns true for containers that have data; false for those that don't (e.g., those that have been co...
TensorPoints()
Default constructor. TensorPoints::isValid() will return false.
void initialize()
Initialize members based on constructor parameters.
TensorPoints(ScalarView< PointScalar, DeviceType > points)
Constructor for point set with trivial tensor structure.
KOKKOS_INLINE_FUNCTION std::enable_if< std::is_integral< iType >::value, int >::type extent_int(const iType &r) const
Returns the logical extent in the requested dimension.
TensorPoints(TensorPoints otherPointsContainer, std::vector< int > whichDims)
Constructor that takes a subset of the tensorial components of another points container.
void copyPointsContainer(ScalarView< PointScalar, DeviceType > toPoints, OtherPointsContainer fromPoints)
Copy from one points container, which may be an arbitrary functor, to a DynRankView.
TensorPoints(std::vector< ScalarView< PointScalar, DeviceType > > pointTensorComponents)
Constructor with variable-length std::vector argument.
TensorPoints(Kokkos::Array< ScalarView< PointScalar, DeviceType >, numTensorComponents > pointTensorComponents)
Constructor with fixed-length Kokkos::Array argument.
ScalarView< PointScalar, DeviceType > allocateAndFillExpandedRawPointView() const
This method is for compatibility with existing methods that take raw point views. Note that in genera...
ordinal_type componentPointCount(const ordinal_type &tensorComponentOrdinal) const
Returns the number of points in the indicated component.
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< std::is_integral< iType >::value, size_t >::type extent(const iType &r) const
Returns the logical extent in the requested dimension.