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
24namespace 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 constexpr bool are_accessible =
143 Kokkos::SpaceAccessibility<MemSpaceType,
144 typename decltype(physPoints)::memory_space>::accessible &&
145 Kokkos::SpaceAccessibility<MemSpaceType,
146 typename decltype(refPoints)::memory_space>::accessible;
147
148 static_assert(are_accessible, "CellTools<DeviceType>::mapToPhysicalFrame(..): input/output views' memory spaces are not compatible with DeviceType");
149
150 const auto cellTopo = basis->getBaseCellTopology();
151 const auto numCells = worksetCell.extent(0);
152
153 //points can be rank-2 (P,D), or rank-3 (C,P,D)
154 const auto refPointRank = refPoints.rank();
155 const auto numPoints = (refPointRank == 2 ? refPoints.extent(0) : refPoints.extent(1));
156 const auto basisCardinality = basis->getCardinality();
157 const auto basisDimRange = Kokkos::pair<ordinal_type,ordinal_type>(0, basis->getDomainDimension());
158
159 using valViewType = Kokkos::DynRankView<decltype(basis->getDummyOutputValue()),DeviceType>;
160
161 valViewType vals;
162
163 switch (refPointRank) {
164 case 2: {
165 // refPoints is (P,D): single set of ref. points is mapped to one or multiple physical cells
166 vals = Impl::createMatchingView<valViewType>(physPoints, "CellTools::mapToPhysicalFrame::vals", basisCardinality, numPoints);
167 const auto refPts = Kokkos::subdynrankview(refPoints, Kokkos::ALL(), basisDimRange); //needed for shell topologies
168 basis->getValues(vals, refPts, OPERATOR_VALUE);
169 break;
170 }
171 case 3: {
172 // refPoints is (C,P,D): multiple sets of ref. points are mapped to matching number of physical cells.
173 vals = Impl::createMatchingView<valViewType>(physPoints, "CellTools::mapToPhysicalFrame::vals", numCells, basisCardinality, numPoints);
174 const auto refPts = Kokkos::subdynrankview(refPoints, Kokkos::ALL(), Kokkos::ALL(), basisDimRange); //needed for shell topologies
175 getHGradValues<OPERATOR_VALUE>(vals, refPts, basis.get());
176 break;
177 }
178 }
179
180 using FunctorType = FunctorCellTools::F_mapToPhysicalFrame<PhysPointViewType,WorksetType,valViewType>;
181
182 const auto loopSize = physPoints.extent(0)*physPoints.extent(1);
183 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
184 Kokkos::parallel_for( policy, FunctorType(physPoints, worksetCell, vals) );
185 }
186
187 template<typename DeviceType>
188 template<typename PhysPointViewType,
189 typename RefPointViewType,
190 typename WorksetCellViewType>
191 void
192 CellTools<DeviceType>::
193 mapToPhysicalFrame( PhysPointViewType physPoints,
194 const RefPointViewType refPoints,
195 const WorksetCellViewType worksetCell,
196 const shards::CellTopology cellTopo ) {
197
198 #ifdef HAVE_INTREPID2_DEBUG
199 CellTools_mapToPhysicalFrameArgs( physPoints, refPoints, worksetCell, cellTopo );
200 #endif
201 using nonConstRefPointValueType = typename RefPointViewType::non_const_value_type;
202 auto basis = createHGradBasis<nonConstRefPointValueType,nonConstRefPointValueType>(cellTopo);
203 mapToPhysicalFrame(physPoints,
204 refPoints,
205 worksetCell,
206 basis);
207 }
208
209 template<typename DeviceType>
210 template<typename refSubcellViewType, typename paramViewType>
211 void
212 CellTools<DeviceType>::
213 mapToReferenceSubcell( refSubcellViewType refSubcellPoints,
214 const paramViewType paramPoints,
215 const ordinal_type subcellDim,
216 const ordinal_type subcellOrd,
217 const shards::CellTopology parentCell ) {
218#ifdef HAVE_INTREPID2_DEBUG
219 ordinal_type parentCellDim = parentCell.getDimension();
220 INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(parentCell), std::invalid_argument,
221 ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): the specified cell topology does not have a reference cell.");
222
223 INTREPID2_TEST_FOR_EXCEPTION( subcellDim < 1 ||
224 subcellDim > parentCellDim-1, std::invalid_argument,
225 ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): method defined only for subcells with dimension greater than 0 and less than the cell dimension");
226
227 INTREPID2_TEST_FOR_EXCEPTION( subcellOrd < 0 ||
228 subcellOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(subcellDim)), std::invalid_argument,
229 ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): subcell ordinal out of range.");
230
231 // refSubcellPoints is rank-2 (P,D1), D1 = cell dimension
232 INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.rank() != 2, std::invalid_argument,
233 ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints must have rank 2.");
234 INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.extent(1) != parentCell.getDimension(), std::invalid_argument,
235 ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints dimension (1) does not match to parent cell dimension.");
236
237 // paramPoints is rank-2 (P,D2) with D2 = subcell dimension
238 INTREPID2_TEST_FOR_EXCEPTION( paramPoints.rank() != 2, std::invalid_argument,
239 ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): paramPoints must have rank 2.");
240 INTREPID2_TEST_FOR_EXCEPTION( static_cast<ordinal_type>(paramPoints.extent(1)) != subcellDim, std::invalid_argument,
241 ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): paramPoints dimension (1) does not match to subcell dimension.");
242
243 // cross check: refSubcellPoints and paramPoints: dimension 0 must match
244 INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.extent(0) != paramPoints.extent(0), std::invalid_argument,
245 ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints dimension (0) does not match to paramPoints dimension(0).");
246#endif
247
248 // Get the subcell map, i.e., the coefficients of the parametrization function for the subcell
249 const auto subcellMap = RefSubcellParametrization<DeviceType>::get(subcellDim, parentCell.getKey());
250
251 mapToReferenceSubcell(refSubcellPoints, paramPoints, subcellMap, subcellOrd);
252 }
253
254 template<typename DeviceType>
255 template<typename refSubcellViewType, typename paramViewType>
256 void
257 CellTools<DeviceType>::
258 mapToReferenceSubcell( refSubcellViewType refSubcellPoints,
259 const paramViewType paramPoints,
260 const typename RefSubcellParametrization<DeviceType>::ConstViewType subcellParametrization,
261 const ordinal_type subcellOrd) {
262
263 constexpr bool are_accessible =
264 Kokkos::SpaceAccessibility<MemSpaceType,
265 typename decltype(refSubcellPoints)::memory_space>::accessible &&
266 Kokkos::SpaceAccessibility<MemSpaceType,
267 typename decltype(paramPoints)::memory_space>::accessible;
268 static_assert(are_accessible, "CellTools<DeviceType>::mapToReferenceSubcell(..): input/output views' memory spaces are not compatible with DeviceType");
269
270
271 #ifdef HAVE_INTREPID2_DEBUG
272 const bool ranks_and_dims_compatible = (refSubcellPoints.rank() == 2) && (paramPoints.rank() == 2) && (subcellParametrization.rank() == 3) &&
273 (refSubcellPoints.extent(0) == paramPoints.extent(0)) &&
274 (refSubcellPoints.extent(1) == subcellParametrization.extent(1)) &&
275 (paramPoints.extent(1) == subcellParametrization.extent(2)-1);
276
277 INTREPID2_TEST_FOR_EXCEPTION(!ranks_and_dims_compatible, std::invalid_argument, "CellTools<DeviceType>::mapToReferenceSubcell(..): input/output views' ranks and dimensions are not compatible");
278 #endif
279
280 const ordinal_type parentCellDim = subcellParametrization.extent(1);
281 const ordinal_type numPts = paramPoints.extent(0);
282
283 //Note: this function has several template parameters and the compiler gets confused if using a lambda function
284 using FunctorType = FunctorCellTools::F_mapReferenceSubcell<decltype(refSubcellPoints), decltype(paramPoints), decltype(subcellParametrization)>;
285
286 Kokkos::MDRangePolicy<ExecSpaceType,Kokkos::Rank<2>> rangePolicy({0,0},{numPts,parentCellDim});
287
288 // Apply the parametrization map to every point in parameter domain
289 Kokkos::parallel_for( rangePolicy, FunctorType(refSubcellPoints, paramPoints, subcellParametrization, subcellOrd) );
290 }
291
292 template<typename DeviceType>
293 template<typename refSubcellViewType, typename paramViewType, typename ordViewType>
294 void
295 CellTools<DeviceType>::
296 mapToReferenceSubcell( refSubcellViewType refSubcellPoints,
297 const paramViewType paramPoints,
298 const typename RefSubcellParametrization<DeviceType>::ConstViewType subcellParametrization,
299 const ordViewType subcellOrd) {
300
301 constexpr bool are_accessible =
302 Kokkos::SpaceAccessibility<MemSpaceType,
303 typename decltype(refSubcellPoints)::memory_space>::accessible &&
304 Kokkos::SpaceAccessibility<MemSpaceType,
305 typename decltype(paramPoints)::memory_space>::accessible &&
306 Kokkos::SpaceAccessibility<MemSpaceType,
307 typename decltype(subcellOrd)::memory_space>::accessible;
308 static_assert(are_accessible, "CellTools<DeviceType>::mapToReferenceSubcell(..): input/output views' memory spaces are not compatible with DeviceType");
309
310
311#ifdef HAVE_INTREPID2_DEBUG
312 const bool ranks_and_dims_compatible = (refSubcellPoints.rank() == 3) && (paramPoints.rank() == 2) && (subcellParametrization.rank() == 3) &&
313 (refSubcellPoints.extent(0) == subcellOrd.extent(0)) &&
314 (refSubcellPoints.extent(1) == paramPoints.extent(0)) &&
315 (refSubcellPoints.extent(2) == subcellParametrization.extent(1)) &&
316 (paramPoints.extent(1) == subcellParametrization.extent(2)-1);
317
318 INTREPID2_TEST_FOR_EXCEPTION(!ranks_and_dims_compatible, std::invalid_argument, "CellTools<DeviceType>::mapToReferenceSubcell(..): input/output views' ranks and dimensions are not compatible");
319#endif
320
321 const ordinal_type numSubCells = refSubcellPoints.extent(0);
322 const ordinal_type parentCellDim = subcellParametrization.extent(1);
323 const ordinal_type numPts = paramPoints.extent(0);
324
325
326 //Note: this function has several template parameters and the compiler gets confused if using a lambda function
327 using FunctorType = FunctorCellTools::F_mapReferenceSubcellBatch<decltype(refSubcellPoints), decltype(paramPoints), decltype(subcellParametrization), decltype(subcellOrd)>;
328
329 Kokkos::MDRangePolicy<ExecSpaceType,Kokkos::Rank<3>> rangePolicy({0,0,0},{numSubCells,numPts,parentCellDim});
330
331 // Apply the parametrization map to every point in parameter domain
332 Kokkos::parallel_for( rangePolicy, FunctorType(refSubcellPoints, paramPoints, subcellParametrization, subcellOrd) );
333 }
334
335}
336
337#endif
Intrepid2::CellTools_mapToPhysicalFrameArgs
void CellTools_mapToPhysicalFrameArgs(const physPointViewType physPoints, const refPointViewType refPoints, const worksetCellViewType worksetCell, const shards::CellTopology cellTopo)
Validates arguments to Intrepid2::CellTools::mapToPhysicalFrame.
Definition Intrepid2_CellToolsDefValidateArguments.hpp:184
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.
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::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_mapReferenceSubcellBatch
Definition Intrepid2_CellToolsDefRefToPhys.hpp:106
Intrepid2::FunctorCellTools::F_mapReferenceSubcell
Definition Intrepid2_CellToolsDefRefToPhys.hpp:78
Intrepid2::FunctorCellTools::F_mapToPhysicalFrame
Functor for mapping reference points to physical frame see Intrepid2::CellTools for more.
Definition Intrepid2_CellToolsDefRefToPhys.hpp:39
Generated by doxygen 1.9.8