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