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#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::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:190
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