Intrepid2
Intrepid2_CellToolsDefRefToPhys.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Intrepid2 Package
4 //
5 // Copyright 2007 NTESS and the Intrepid2 contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 
16 #ifndef __INTREPID2_CELLTOOLS_DEF_REF_TO_PHYS_HPP__
17 #define __INTREPID2_CELLTOOLS_DEF_REF_TO_PHYS_HPP__
18 
19 // disable clang warnings
20 #if defined (__clang__) && !defined (__INTEL_COMPILER)
21 #pragma clang system_header
22 #endif
23 
24 namespace Intrepid2 {
25 
26  //============================================================================================//
27  // //
28  // Reference-to-physical frame mapping and its inverse //
29  // //
30  //============================================================================================//
31 
32  namespace FunctorCellTools {
36  template<typename physPointViewType,
37  typename worksetCellType,
38  typename basisValType>
39  struct F_mapToPhysicalFrame {
40  physPointViewType _physPoints;
41  const worksetCellType _worksetCells;
42  const basisValType _basisVals;
43 
44  KOKKOS_INLINE_FUNCTION
45  F_mapToPhysicalFrame( physPointViewType physPoints_,
46  worksetCellType worksetCells_,
47  basisValType basisVals_ )
48  : _physPoints(physPoints_), _worksetCells(worksetCells_), _basisVals(basisVals_) {}
49 
50  KOKKOS_INLINE_FUNCTION
51  void operator()(const size_type iter) const {
52  size_type cell, pt;
53  unrollIndex( cell, pt,
54  _physPoints.extent(0),
55  _physPoints.extent(1),
56  iter );
57  auto phys = Kokkos::subdynrankview( _physPoints, cell, pt, Kokkos::ALL());
58 
59  const auto valRank = rank(_basisVals);
60  const auto val = ( valRank == 2 ? Kokkos::subdynrankview( _basisVals, Kokkos::ALL(), pt) :
61  Kokkos::subdynrankview( _basisVals, cell, Kokkos::ALL(), pt));
62 
63  const ordinal_type dim = phys.extent(0);
64  const ordinal_type cardinality = val.extent(0);
65 
66  for (ordinal_type i=0;i<dim;++i) {
67  phys(i) = 0;
68  for (ordinal_type bf=0;bf<cardinality;++bf)
69  phys(i) += _worksetCells(cell, bf, i)*val(bf);
70  }
71  }
72  };
73 
74 
75  template<typename refSubcellViewType,
76  typename paramPointsViewType,
77  typename subcellMapViewType>
78  struct F_mapReferenceSubcell {
79  refSubcellViewType refSubcellPoints_;
80  const paramPointsViewType paramPoints_;
81  const subcellMapViewType subcellMap_;
82  const ordinal_type subcellOrd_;
83  const ordinal_type sideDim_;
84 
85  KOKKOS_INLINE_FUNCTION
86  F_mapReferenceSubcell( refSubcellViewType refSubcellPoints,
87  const paramPointsViewType paramPoints,
88  const subcellMapViewType subcellMap,
89  ordinal_type subcellOrd)
90  : refSubcellPoints_(refSubcellPoints), paramPoints_(paramPoints), subcellMap_(subcellMap),
91  subcellOrd_(subcellOrd), sideDim_(paramPoints_.extent(1)){};
92 
93  KOKKOS_INLINE_FUNCTION
94  void operator()(const size_type pt, const size_type d) const {
95 
96  refSubcellPoints_(pt, d) = subcellMap_(subcellOrd_, d, 0);
97  for(ordinal_type k=0; k<sideDim_; ++k)
98  refSubcellPoints_(pt, d) += subcellMap_(subcellOrd_, d, k+1)*paramPoints_(pt, k);
99  }
100  };
101 
102  template<typename refSubcellViewType,
103  typename paramPointsViewType,
104  typename subcellMapViewType,
105  typename ordViewType>
106  struct F_mapReferenceSubcellBatch {
107  refSubcellViewType refSubcellPoints_;
108  const paramPointsViewType paramPoints_;
109  const subcellMapViewType subcellMap_;
110  const ordViewType subcellOrd_;
111  const ordinal_type sideDim_;
112 
113  KOKKOS_INLINE_FUNCTION
114  F_mapReferenceSubcellBatch( refSubcellViewType refSubcellPoints,
115  const paramPointsViewType paramPoints,
116  const subcellMapViewType subcellMap,
117  ordViewType subcellOrd)
118  : refSubcellPoints_(refSubcellPoints), paramPoints_(paramPoints), subcellMap_(subcellMap),
119  subcellOrd_(subcellOrd), sideDim_(paramPoints_.extent(1)){};
120 
121  KOKKOS_INLINE_FUNCTION
122  void operator()(const size_type isc, const size_type pt, const size_type d) const {
123 
124  refSubcellPoints_(isc, pt, d) = subcellMap_(subcellOrd_(isc), d, 0);
125  for(ordinal_type k=0; k<sideDim_; ++k)
126  refSubcellPoints_(isc, pt, d) += subcellMap_(subcellOrd_(isc), d, k+1)*paramPoints_(pt, k);
127  }
128  };
129  }
130 
131  template<typename DeviceType>
132  template<typename PhysPointViewType,
133  typename RefPointViewType,
134  typename WorksetType,
135  typename HGradBasisPtrType>
136  void
137  CellTools<DeviceType>::
138  mapToPhysicalFrame( PhysPointViewType physPoints,
139  const RefPointViewType refPoints,
140  const WorksetType worksetCell,
141  const HGradBasisPtrType basis ) {
142 #ifdef HAVE_INTREPID2_DEBUG
143  CellTools_mapToPhysicalFrameArgs( physPoints, refPoints, worksetCell, basis->getBaseCellTopology() );
144 #endif
145  constexpr bool are_accessible =
146  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
147  typename decltype(physPoints)::memory_space>::accessible &&
148  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
149  typename decltype(refPoints)::memory_space>::accessible;
150 
151  static_assert(are_accessible, "CellTools<DeviceType>::mapToPhysicalFrame(..): input/output views' memory spaces are not compatible with DeviceType");
152 
153  const auto cellTopo = basis->getBaseCellTopology();
154  const auto numCells = worksetCell.extent(0);
155 
156  //points can be rank-2 (P,D), or rank-3 (C,P,D)
157  const auto refPointRank = refPoints.rank();
158  const auto numPoints = (refPointRank == 2 ? refPoints.extent(0) : refPoints.extent(1));
159  const auto basisCardinality = basis->getCardinality();
160  auto vcprop = Kokkos::common_view_alloc_prop(physPoints);
161 
162  using valViewType = Kokkos::DynRankView<decltype(basis->getDummyOutputValue()),DeviceType>;
163 
164  valViewType vals;
165 
166  switch (refPointRank) {
167  case 2: {
168  // refPoints is (P,D): single set of ref. points is mapped to one or multiple physical cells
169  vals = valViewType(Kokkos::view_alloc("CellTools::mapToPhysicalFrame::vals", vcprop), basisCardinality, numPoints);
170  basis->getValues(vals,
171  refPoints,
172  OPERATOR_VALUE);
173  break;
174  }
175  case 3: {
176  // refPoints is (C,P,D): multiple sets of ref. points are mapped to matching number of physical cells.
177  //vals = valViewType("CellTools::mapToPhysicalFrame::vals", numCells, basisCardinality, numPoints);
178  vals = valViewType(Kokkos::view_alloc("CellTools::mapToPhysicalFrame::vals", vcprop), numCells, basisCardinality, numPoints);
179  for (size_type cell=0;cell<numCells;++cell)
180  basis->getValues(Kokkos::subdynrankview( vals, cell, Kokkos::ALL(), Kokkos::ALL() ),
181  Kokkos::subdynrankview( refPoints, cell, Kokkos::ALL(), Kokkos::ALL() ),
182  OPERATOR_VALUE);
183  break;
184  }
185  }
186 
187  using FunctorType = FunctorCellTools::F_mapToPhysicalFrame<PhysPointViewType,WorksetType,valViewType>;
188 
189  const auto loopSize = physPoints.extent(0)*physPoints.extent(1);
190  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
191  Kokkos::parallel_for( policy, FunctorType(physPoints, worksetCell, vals) );
192  }
193 
194  template<typename DeviceType>
195  template<typename refSubcellViewType, typename paramViewType>
196  void
197  CellTools<DeviceType>::
198  mapToReferenceSubcell( refSubcellViewType refSubcellPoints,
199  const paramViewType paramPoints,
200  const ordinal_type subcellDim,
201  const ordinal_type subcellOrd,
202  const shards::CellTopology parentCell ) {
203 #ifdef HAVE_INTREPID2_DEBUG
204  ordinal_type parentCellDim = parentCell.getDimension();
205  INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(parentCell), std::invalid_argument,
206  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): the specified cell topology does not have a reference cell.");
207 
208  INTREPID2_TEST_FOR_EXCEPTION( subcellDim < 1 ||
209  subcellDim > parentCellDim-1, std::invalid_argument,
210  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): method defined only for subcells with dimension greater than 0 and less than the cell dimension");
211 
212  INTREPID2_TEST_FOR_EXCEPTION( subcellOrd < 0 ||
213  subcellOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(subcellDim)), std::invalid_argument,
214  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): subcell ordinal out of range.");
215 
216  // refSubcellPoints is rank-2 (P,D1), D1 = cell dimension
217  INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.rank() != 2, std::invalid_argument,
218  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints must have rank 2.");
219  INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.extent(1) != parentCell.getDimension(), std::invalid_argument,
220  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints dimension (1) does not match to parent cell dimension.");
221 
222  // paramPoints is rank-2 (P,D2) with D2 = subcell dimension
223  INTREPID2_TEST_FOR_EXCEPTION( paramPoints.rank() != 2, std::invalid_argument,
224  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): paramPoints must have rank 2.");
225  INTREPID2_TEST_FOR_EXCEPTION( static_cast<ordinal_type>(paramPoints.extent(1)) != subcellDim, std::invalid_argument,
226  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): paramPoints dimension (1) does not match to subcell dimension.");
227 
228  // cross check: refSubcellPoints and paramPoints: dimension 0 must match
229  INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.extent(0) != paramPoints.extent(0), std::invalid_argument,
230  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints dimension (0) does not match to paramPoints dimension(0).");
231 #endif
232 
233  // Get the subcell map, i.e., the coefficients of the parametrization function for the subcell
234  const auto subcellMap = RefSubcellParametrization<DeviceType>::get(subcellDim, parentCell.getKey());
235 
236  mapToReferenceSubcell(refSubcellPoints, paramPoints, subcellMap, subcellOrd);
237  }
238 
239  template<typename DeviceType>
240  template<typename refSubcellViewType, typename paramViewType>
241  void
242  CellTools<DeviceType>::
243  mapToReferenceSubcell( refSubcellViewType refSubcellPoints,
244  const paramViewType paramPoints,
245  const typename RefSubcellParametrization<DeviceType>::ConstViewType subcellParametrization,
246  const ordinal_type subcellOrd) {
247 
248  constexpr bool are_accessible =
249  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
250  typename decltype(refSubcellPoints)::memory_space>::accessible &&
251  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
252  typename decltype(paramPoints)::memory_space>::accessible;
253  static_assert(are_accessible, "CellTools<DeviceType>::mapToReferenceSubcell(..): input/output views' memory spaces are not compatible with DeviceType");
254 
255 
256  #ifdef HAVE_INTREPID2_DEBUG
257  const bool ranks_and_dims_compatible = (refSubcellPoints.rank() == 2) && (paramPoints.rank() == 2) && (subcellParametrization.rank() == 3) &&
258  (refSubcellPoints.extent(0) == paramPoints.extent(0)) &&
259  (refSubcellPoints.extent(1) == subcellParametrization.extent(1)) &&
260  (paramPoints.extent(1) == subcellParametrization.extent(2)-1);
261 
262  INTREPID2_TEST_FOR_EXCEPTION(!ranks_and_dims_compatible, std::invalid_argument, "CellTools<DeviceType>::mapToReferenceSubcell(..): input/output views' ranks and dimensions are not compatible");
263  #endif
264 
265  const ordinal_type parentCellDim = subcellParametrization.extent(1);
266  const ordinal_type numPts = paramPoints.extent(0);
267 
268  //Note: this function has several template parameters and the compiler gets confused if using a lambda function
269  using FunctorType = FunctorCellTools::F_mapReferenceSubcell<decltype(refSubcellPoints), decltype(paramPoints), decltype(subcellParametrization)>;
270 
271  Kokkos::MDRangePolicy<ExecSpaceType,Kokkos::Rank<2>> rangePolicy({0,0},{numPts,parentCellDim});
272 
273  // Apply the parametrization map to every point in parameter domain
274  Kokkos::parallel_for( rangePolicy, FunctorType(refSubcellPoints, paramPoints, subcellParametrization, subcellOrd) );
275  }
276 
277  template<typename DeviceType>
278  template<typename refSubcellViewType, typename paramViewType, typename ordViewType>
279  void
280  CellTools<DeviceType>::
281  mapToReferenceSubcell( refSubcellViewType refSubcellPoints,
282  const paramViewType paramPoints,
283  const typename RefSubcellParametrization<DeviceType>::ConstViewType subcellParametrization,
284  const ordViewType subcellOrd) {
285 
286  constexpr bool are_accessible =
287  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
288  typename decltype(refSubcellPoints)::memory_space>::accessible &&
289  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
290  typename decltype(paramPoints)::memory_space>::accessible &&
291  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
292  typename decltype(subcellOrd)::memory_space>::accessible;
293  static_assert(are_accessible, "CellTools<DeviceType>::mapToReferenceSubcell(..): input/output views' memory spaces are not compatible with DeviceType");
294 
295 
296 #ifdef HAVE_INTREPID2_DEBUG
297  const bool ranks_and_dims_compatible = (refSubcellPoints.rank() == 3) && (paramPoints.rank() == 2) && (subcellParametrization.rank() == 3) &&
298  (refSubcellPoints.extent(0) == subcellOrd.extent(0)) &&
299  (refSubcellPoints.extent(1) == paramPoints.extent(0)) &&
300  (refSubcellPoints.extent(2) == subcellParametrization.extent(1)) &&
301  (paramPoints.extent(1) == subcellParametrization.extent(2)-1);
302 
303  INTREPID2_TEST_FOR_EXCEPTION(!ranks_and_dims_compatible, std::invalid_argument, "CellTools<DeviceType>::mapToReferenceSubcell(..): input/output views' ranks and dimensions are not compatible");
304 #endif
305 
306  const ordinal_type numSubCells = refSubcellPoints.extent(0);
307  const ordinal_type parentCellDim = subcellParametrization.extent(1);
308  const ordinal_type numPts = paramPoints.extent(0);
309 
310 
311  //Note: this function has several template parameters and the compiler gets confused if using a lambda function
312  using FunctorType = FunctorCellTools::F_mapReferenceSubcellBatch<decltype(refSubcellPoints), decltype(paramPoints), decltype(subcellParametrization), decltype(subcellOrd)>;
313 
314  Kokkos::MDRangePolicy<ExecSpaceType,Kokkos::Rank<3>> rangePolicy({0,0,0},{numSubCells,numPts,parentCellDim});
315 
316  // Apply the parametrization map to every point in parameter domain
317  Kokkos::parallel_for( rangePolicy, FunctorType(refSubcellPoints, paramPoints, subcellParametrization, subcellOrd) );
318  }
319 
320 }
321 
322 #endif
Intrepid2::FunctorCellTools::F_mapToPhysicalFrame
Functor for mapping reference points to physical frame see Intrepid2::CellTools for more...
Definition: Intrepid2_CellToolsDefRefToPhys.hpp:39
Intrepid2::RefSubcellParametrization::get
static ConstViewType get(const ordinal_type subcellDim, const unsigned parentCellKey)
Returns a Kokkos view with the coefficients of the parametrization maps for the edges or faces of a r...
Definition: Intrepid2_CellDataDef.hpp:69
Intrepid2::FunctorCellTools::F_mapReferenceSubcell
Definition: Intrepid2_CellToolsDefRefToPhys.hpp:78
Intrepid2::FunctorCellTools::F_mapReferenceSubcellBatch
Definition: Intrepid2_CellToolsDefRefToPhys.hpp:106
Intrepid2::CellTools::mapToReferenceSubcell
static void mapToReferenceSubcell(refSubcellViewType refSubcellPoints, const paramPointViewType paramPoints, const ordinal_type subcellDim, const ordinal_type subcellOrd, const shards::CellTopology parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
Intrepid2::CellTools::mapToPhysicalFrame
static void mapToPhysicalFrame(PhysPointValueType physPoints, const RefPointValueType refPoints, const WorksetType worksetCell, const HGradBasisPtrType basis)
Computes F, the reference-to-physical frame map.
Generated by   doxygen 1.8.5