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 
55 namespace Intrepid2 {
59  template<class PointScalar, typename DeviceType>
60  class TensorPoints {
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 
77  void initialize()
78  {
79  totalPointCount_ = 1;
80  totalDimension_ = 0;
81  for (ordinal_type r=0; r<numTensorComponents_; r++)
82  {
83  totalPointCount_ *= pointTensorComponents_[r].extent_int(0);
84  totalDimension_ += pointTensorComponents_[r].extent_int(1);
85  }
86  ordinal_type pointDivisor = 1;
87  for (ordinal_type r=0; r<numTensorComponents_; r++)
88  {
89  pointModulus_[r] = pointTensorComponents_[r].extent_int(0);
90  pointDivisor_[r] = pointDivisor;
91  pointDivisor *= pointTensorComponents_[r].extent_int(0);
92  }
93  dimToComponent_ = Kokkos::View<ordinal_type*, DeviceType>("dimToComponent_",totalDimension_);
94  dimToComponentDim_ = Kokkos::View<ordinal_type*, DeviceType>("dimToComponentDim_",totalDimension_);
95  ordinal_type d=0;
96  ordinal_type dimsSoFar = 0;
97 
98  auto dimToComponentHost = Kokkos::create_mirror_view(dimToComponent_);
99  auto dimToComponentDimHost = Kokkos::create_mirror_view(dimToComponentDim_);
100  for (ordinal_type r=0; r<numTensorComponents_; r++)
101  {
102  const int componentDim = pointTensorComponents_[r].extent_int(1);
103  for (int i=0; i<componentDim; i++)
104  {
105  dimToComponentHost[d] = r;
106  dimToComponentDimHost[d] = d - dimsSoFar;
107  d++;
108  }
109  dimsSoFar += componentDim;
110  }
111  Kokkos::deep_copy(dimToComponent_,dimToComponentHost);
112  Kokkos::deep_copy(dimToComponentDim_,dimToComponentDimHost);
113  }
114  public:
121  template<size_t numTensorComponents>
122  TensorPoints(Kokkos::Array< ScalarView<PointScalar,DeviceType>, numTensorComponents> pointTensorComponents)
123  :
124  numTensorComponents_(numTensorComponents),
125  isValid_(true)
126  {
127  for (unsigned r=0; r<numTensorComponents; r++)
128  {
129  pointTensorComponents_[r] = pointTensorComponents[r];
130  }
131 
132  initialize();
133  }
134 
141  TensorPoints(std::vector< ScalarView<PointScalar,DeviceType>> pointTensorComponents)
142  :
143  numTensorComponents_(pointTensorComponents.size()),
144  isValid_(true)
145  {
146  for (ordinal_type r=0; r<numTensorComponents_; r++)
147  {
148  pointTensorComponents_[r] = pointTensorComponents[r];
149  }
150 
151  initialize();
152  }
153 
160  TensorPoints(ScalarView<PointScalar,DeviceType> points)
161  :
162  numTensorComponents_(1),
163  isValid_(true)
164  {
165  pointTensorComponents_[0] = points;
166  initialize();
167  }
168 
174  template<class OtherPointsContainer>
175  void copyPointsContainer(ScalarView<PointScalar,DeviceType> toPoints, OtherPointsContainer fromPoints)
176  {
177  const int numPoints = fromPoints.extent_int(0);
178  const int numDims = fromPoints.extent_int(1);
179  using ExecutionSpace = typename DeviceType::execution_space;
180  auto policy = Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<2>>({0,0},{numPoints,numDims});
181  Kokkos::parallel_for("copy points", policy,
182  KOKKOS_LAMBDA (const int &i0, const int &i1) {
183  toPoints(i0,i1) = fromPoints(i0,i1);
184  });
185  }
186 
188  template<typename OtherDeviceType, class = typename std::enable_if< std::is_same<typename DeviceType::memory_space, typename OtherDeviceType::memory_space>::value>::type,
189  class = typename std::enable_if<!std::is_same<DeviceType,OtherDeviceType>::value>::type>
190  TensorPoints(const TensorPoints<PointScalar,OtherDeviceType> &tensorPoints)
191  :
192  numTensorComponents_(tensorPoints.numTensorComponents()),
193  isValid_(tensorPoints.isValid())
194  {
195  if (isValid_)
196  {
197  for (ordinal_type r=0; r<numTensorComponents_; r++)
198  {
199  pointTensorComponents_[r] = tensorPoints.getTensorComponent(r);
200  }
201  initialize();
202  }
203  }
204 
206  template<typename OtherDeviceType, class = typename std::enable_if<!std::is_same<typename DeviceType::memory_space, typename OtherDeviceType::memory_space>::value>::type>
208  :
209  numTensorComponents_(tensorPoints.numTensorComponents()),
210  isValid_(tensorPoints.isValid())
211  {
212  if (isValid_)
213  {
214  for (ordinal_type r=0; r<numTensorComponents_; r++)
215  {
216  ScalarView<PointScalar,OtherDeviceType> otherPointComponent = tensorPoints.getTensorComponent(r);
217  const int numPoints = otherPointComponent.extent_int(0);
218  const int numDims = otherPointComponent.extent_int(1);
219  pointTensorComponents_[r] = ScalarView<PointScalar,DeviceType>("Intrepid2 point component", numPoints, numDims);
220 
221  using MemorySpace = typename DeviceType::memory_space;
222  auto pointComponentMirror = Kokkos::create_mirror_view_and_copy(MemorySpace(), otherPointComponent);
223 
224  copyPointsContainer(pointTensorComponents_[r], pointComponentMirror);
225  }
226  initialize();
227  }
228  }
229 
232  isValid_(false)
233  // when constructed with default arguments, TensorPoints should not be used…
234  // default constructor is only provided for things like CellGeometry, which has TensorPoints as a member,
235  // but only uses it in certain modes.
236  {}
237 
239  ordinal_type componentPointCount(const ordinal_type &tensorComponentOrdinal) const
240  {
241  return pointTensorComponents_[tensorComponentOrdinal].extent_int(0);
242  }
243 
252  template <typename iType0, typename iType1>
253  KOKKOS_INLINE_FUNCTION typename std::enable_if<
254  (std::is_integral<iType0>::value && std::is_integral<iType1>::value),
255  reference_type>::type
256  operator()(const iType0& tensorPointIndex, const iType1& dim) const {
257  const ordinal_type component = dimToComponent_[dim];
258  const ordinal_type d = dimToComponentDim_[dim];
259  const ordinal_type componentPointOrdinal = (tensorPointIndex / pointDivisor_[component]) % pointModulus_[component];
260  return pointTensorComponents_[component](componentPointOrdinal,d);
261  }
262 
271  template <typename iType0, typename iType1, size_t numTensorComponents>
272  KOKKOS_INLINE_FUNCTION typename std::enable_if<
273  (std::is_integral<iType0>::value && std::is_integral<iType1>::value),
274  reference_type>::type
275  operator()(const Kokkos::Array<iType0,numTensorComponents>& pointOrdinalComponents, const iType1& dim) const {
276  const ordinal_type component = dimToComponent_[dim];
277  const ordinal_type d = dimToComponentDim_[dim];
278  const ordinal_type componentPointOrdinal = pointOrdinalComponents[component];
279  return pointTensorComponents_[component](componentPointOrdinal,d);
280  }
281 
287  template <typename iType>
288  KOKKOS_INLINE_FUNCTION
289  typename std::enable_if<std::is_integral<iType>::value, int>::type
290  extent_int(const iType& r) const {
291  if (r == static_cast<iType>(0))
292  {
293  return totalPointCount_;
294  }
295  else if (r == static_cast<iType>(1))
296  {
297  return totalDimension_;
298  }
299  else
300  {
301  return 1;
302  }
303  }
304 
310  template <typename iType>
311  KOKKOS_INLINE_FUNCTION constexpr
312  typename std::enable_if<std::is_integral<iType>::value, size_t>::type
313  extent(const iType& r) const {
314  // C++ prior to 14 doesn't allow constexpr if statements; the compound ternary is here to keep things constexpr
315  return (r == static_cast<iType>(0)) ? totalPointCount_
316  :
317  (r == static_cast<iType>(1)) ? totalDimension_ : 1;
318  }
319 
321  ScalarView<PointScalar,DeviceType> allocateAndFillExpandedRawPointView() const
322  {
323  const int numPoints = this->extent_int(0);
324  const int spaceDim = this->extent_int(1);
325  ScalarView<PointScalar,DeviceType> expandedRawPoints("expanded raw points from TensorPoints", numPoints, spaceDim);
326  TensorPoints<PointScalar,DeviceType> tensorPoints(*this); // (shallow) copy for lambda capture
327  using ExecutionSpace = typename DeviceType::execution_space;
328  Kokkos::parallel_for(
329  Kokkos::MDRangePolicy<ExecutionSpace,Kokkos::Rank<2>>({0,0},{numPoints,spaceDim}),
330  KOKKOS_LAMBDA (const int &pointOrdinal, const int &d) {
331  expandedRawPoints(pointOrdinal,d) = tensorPoints(pointOrdinal,d);
332  });
333  return expandedRawPoints;
334  }
335 
341  KOKKOS_INLINE_FUNCTION
342  ScalarView<PointScalar,DeviceType> getTensorComponent(const ordinal_type &r) const
343  {
344  return pointTensorComponents_[r];
345  }
346 
348  KOKKOS_INLINE_FUNCTION
349  bool isValid() const
350  {
351  return isValid_;
352  }
353 
355  KOKKOS_INLINE_FUNCTION
356  ordinal_type numTensorComponents() const
357  {
358  return numTensorComponents_;
359  }
360 
362  KOKKOS_INLINE_FUNCTION
363  constexpr ordinal_type rank() const
364  {
365  return 2; // shape is (P,D)
366  }
367 
368  // NOTE: the extractTensorPoints() code commented out below appears not to work. We don't have tests against it, though.
369  // TODO: either delete this, or re-enable, add tests, and fix.
370 // template<class ViewType>
371 // static TensorPoints<PointScalar,DeviceType> extractTensorPoints( ViewType expandedPoints, const std::vector<ordinal_type> &dimensionExtents )
372 // {
373 // const ordinal_type numComponents = dimensionExtents.size();
374 // const ordinal_type numPoints = expandedPoints.extent_int(0);
375 // Kokkos::Array<ordinal_type,Parameters::MaxTensorComponents> componentPointCounts;
376 // Kokkos::Array<ordinal_type,Parameters::MaxTensorComponents> componentPointTensorStride;
377 // Kokkos::Array<ordinal_type,Parameters::MaxTensorComponents> componentPointDimOffsets;
378 //
379 // // for simplicity of implementation, we copy to host:
380 // auto hostExpandedPoints = Kokkos::create_mirror_view_and_copy(typename Kokkos::HostSpace::memory_space(), expandedPoints);
381 //
382 // ordinal_type dimOffset = 0;
383 // ordinal_type tensorPointStride = 1; // increases with componentOrdinal
384 //
385 // TensorPoints<PointScalar,DeviceType> invalidTensorData; // will be returned if extraction does not succeed.
386 //
387 // for (ordinal_type componentOrdinal=0; componentOrdinal<numComponents; componentOrdinal++)
388 // {
389 // componentPointDimOffsets[componentOrdinal] = dimOffset;
390 // componentPointTensorStride[componentOrdinal] = tensorPointStride;
391 // const ordinal_type numDimsForComponent = dimensionExtents[componentOrdinal];
392 // std::vector<PointScalar> firstPoint(numDimsForComponent);
393 // for (ordinal_type d=0; d<numDimsForComponent; d++)
394 // {
395 // firstPoint[d] = hostExpandedPoints(0,d+dimOffset);
396 // }
397 //
398 // // we assume that once we see the same point twice, we've found the cycle length:
399 // componentPointCounts[componentOrdinal] = -1;
400 // for (ordinal_type pointOrdinal=1; pointOrdinal<numPoints; pointOrdinal += tensorPointStride)
401 // {
402 // bool matches = true;
403 // for (ordinal_type d=0; d<numDimsForComponent; d++)
404 // {
405 // matches = matches && (firstPoint[d] == hostExpandedPoints(pointOrdinal,d+dimOffset));
406 // }
407 // if (matches)
408 // {
409 // componentPointCounts[componentOrdinal] = pointOrdinal;
410 // break;
411 // }
412 // }
413 // if (componentPointCounts[componentOrdinal] == -1)
414 // {
415 // // no matches found -> no tensor decomposition available
416 // return invalidTensorData;
417 // }
418 //
419 // // check that we got the cycle length correct:
420 // for (ordinal_type pointOrdinal=0; pointOrdinal<componentPointCounts[componentOrdinal]; pointOrdinal += tensorPointStride)
421 // {
422 // std::vector<PointScalar> point(numDimsForComponent);
423 // for (ordinal_type d=0; d<numDimsForComponent; d++)
424 // {
425 // point[d] = hostExpandedPoints(pointOrdinal,d+dimOffset);
426 // }
427 // // each of the following points should match:
428 // for (ordinal_type secondPointOrdinal=0; secondPointOrdinal<numPoints; secondPointOrdinal += tensorPointStride*componentPointCounts[componentOrdinal])
429 // {
430 // bool matches = true;
431 // for (ordinal_type d=0; d<numDimsForComponent; d++)
432 // {
433 // matches = matches && (point[d] == hostExpandedPoints(secondPointOrdinal,d+dimOffset));
434 // }
435 // if (!matches)
436 // {
437 // // fail:
438 // return invalidTensorData;
439 // }
440 // }
441 // }
442 //
443 // dimOffset += numDimsForComponent;
444 // tensorPointStride *= componentPointCounts[componentOrdinal];
445 // }
446 //
447 // std::vector< ScalarView<PointScalar,DeviceType> > componentPoints(numComponents);
448 // std::vector< ScalarView<PointScalar,Kokkos::HostSpace> > componentPointsHost(numComponents);
449 // for (ordinal_type componentOrdinal=0; componentOrdinal<numComponents; componentOrdinal++)
450 // {
451 // const ordinal_type numPointsForComponent = componentPointCounts[componentOrdinal];
452 // const ordinal_type dimForComponent = dimensionExtents[componentOrdinal];
453 // componentPoints[componentOrdinal] = ScalarView<PointScalar,DeviceType>("extracted tensor components", numPointsForComponent, dimForComponent);
454 // componentPointsHost[componentOrdinal] = Kokkos::create_mirror(componentPoints[componentOrdinal]);
455 // }
456 //
457 // for (ordinal_type componentOrdinal=0; componentOrdinal<numComponents; componentOrdinal++)
458 // {
459 // const ordinal_type numComponentPoints = componentPointCounts[componentOrdinal];
460 //
461 // auto hostView = componentPointsHost[componentOrdinal];
462 // auto deviceView = componentPoints[componentOrdinal];
463 // const ordinal_type tensorPointStride = componentPointTensorStride[componentOrdinal];
464 // const ordinal_type dimOffset = componentPointDimOffsets[componentOrdinal];
465 // const ordinal_type numDimsForComponent = dimensionExtents[componentOrdinal];
466 // for (ordinal_type componentPointOrdinal=0; componentPointOrdinal<numComponentPoints; componentPointOrdinal++)
467 // {
468 // const ordinal_type expandedPointOrdinal = componentPointOrdinal*tensorPointStride;
469 // for (ordinal_type d=0; d<numDimsForComponent; d++)
470 // {
471 // hostView(componentPointOrdinal,d) = hostExpandedPoints(expandedPointOrdinal,d+dimOffset);
472 // }
473 // }
474 // Kokkos::deep_copy(deviceView, hostView);
475 // }
476 //
477 // // prior to return, check all points agree in all dimensions with the input points
478 // // for convenience, we do this check on host, too
479 // TensorPoints<PointScalar,Kokkos::HostSpace> hostTensorPoints(componentPointsHost);
480 // const ordinal_type totalDim = expandedPoints.extent_int(1);
481 // bool matches = true;
482 // for (ordinal_type pointOrdinal=0; pointOrdinal<numPoints; pointOrdinal++)
483 // {
484 // for (ordinal_type d=0; d<totalDim; d++)
485 // {
486 // const auto &originalCoord = hostExpandedPoints(pointOrdinal,d);
487 // const auto &tensorCoord = hostTensorPoints(pointOrdinal,d);
488 // if (originalCoord != tensorCoord)
489 // {
490 // matches = false;
491 // }
492 // }
493 // }
494 // if (!matches)
495 // {
496 // return invalidTensorData;
497 // }
498 //
499 // return TensorPoints<PointScalar,DeviceType>(componentPoints);
500 // }
501  };
502 
503 }
504 
505 #endif /* Intrepid2_TensorPoints_h */
TensorPoints(ScalarView< PointScalar, DeviceType > points)
Constructor for point set with trivial tensor structure.
View-like interface to tensor points; point components are stored separately; the appropriate coordin...
void initialize()
Initialize members based on constructor parameters.
KOKKOS_INLINE_FUNCTION constexpr ordinal_type rank() const
Return the rank of the container, which is 2.
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 std::enable_if< std::is_integral< iType >::value, int >::type extent_int(const iType &r) const
Returns the logical extent in the requested dimension.
void copyPointsContainer(ScalarView< PointScalar, DeviceType > toPoints, OtherPointsContainer fromPoints)
Copy from one points container, which may be an arbitrary functor, to a DynRankView.
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 bool isValid() const
Returns true for containers that have data; false for those that don&#39;t (e.g., those that have been co...
TensorPoints(Kokkos::Array< ScalarView< PointScalar, DeviceType >, numTensorComponents > pointTensorComponents)
Constructor with fixed-length Kokkos::Array argument.
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.
TensorPoints()
Default constructor. TensorPoints::isValid() will return false.
KOKKOS_INLINE_FUNCTION ScalarView< PointScalar, DeviceType > getTensorComponent(const ordinal_type &r) const
Returns the requested tensor component.
KOKKOS_INLINE_FUNCTION ordinal_type numTensorComponents() const
Returns the number of tensorial components.
TensorPoints(std::vector< ScalarView< PointScalar, DeviceType >> pointTensorComponents)
Constructor with variable-length std::vector argument.
static constexpr ordinal_type MaxTensorComponents
Maximum number of tensor/Cartesian products that can be taken: this allows hypercube basis in 7D to b...