Intrepid2
Intrepid2_CellToolsDefRefToPhys.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), or
38// Mauro Perego (mperego@sandia.gov)
39//
40// ************************************************************************
41// @HEADER
42
43
49#ifndef __INTREPID2_CELLTOOLS_DEF_REF_TO_PHYS_HPP__
50#define __INTREPID2_CELLTOOLS_DEF_REF_TO_PHYS_HPP__
51
52// disable clang warnings
53#if defined (__clang__) && !defined (__INTEL_COMPILER)
54#pragma clang system_header
55#endif
56
57namespace Intrepid2 {
58
59 //============================================================================================//
60 // //
61 // Reference-to-physical frame mapping and its inverse //
62 // //
63 //============================================================================================//
64
65 namespace FunctorCellTools {
69 template<typename physPointViewType,
70 typename worksetCellType,
71 typename basisValType>
73 physPointViewType _physPoints;
74 const worksetCellType _worksetCells;
75 const basisValType _basisVals;
76
77 KOKKOS_INLINE_FUNCTION
78 F_mapToPhysicalFrame( physPointViewType physPoints_,
79 worksetCellType worksetCells_,
80 basisValType basisVals_ )
81 : _physPoints(physPoints_), _worksetCells(worksetCells_), _basisVals(basisVals_) {}
82
83 KOKKOS_INLINE_FUNCTION
84 void operator()(const size_type iter) const {
85 size_type cell, pt;
86 unrollIndex( cell, pt,
87 _physPoints.extent(0),
88 _physPoints.extent(1),
89 iter );
90 auto phys = Kokkos::subdynrankview( _physPoints, cell, pt, Kokkos::ALL());
91
92 const auto valRank = _basisVals.rank();
93 const auto val = ( valRank == 2 ? Kokkos::subdynrankview( _basisVals, Kokkos::ALL(), pt) :
94 Kokkos::subdynrankview( _basisVals, cell, Kokkos::ALL(), pt));
95
96 const ordinal_type dim = phys.extent(0);
97 const ordinal_type cardinality = val.extent(0);
98
99 for (ordinal_type i=0;i<dim;++i) {
100 phys(i) = 0;
101 for (ordinal_type bf=0;bf<cardinality;++bf)
102 phys(i) += _worksetCells(cell, bf, i)*val(bf);
103 }
104 }
105 };
106
107
108 template<typename refSubcellViewType,
109 typename paramPointsViewType,
110 typename subcellMapViewType>
112 refSubcellViewType refSubcellPoints_;
113 const paramPointsViewType paramPoints_;
114 const subcellMapViewType subcellMap_;
115 ordinal_type subcellOrd_, dim_;
116
117 KOKKOS_INLINE_FUNCTION
118 F_mapReferenceSubcell2( refSubcellViewType refSubcellPoints,
119 const paramPointsViewType paramPoints,
120 const subcellMapViewType subcellMap,
121 ordinal_type subcellOrd,
122 ordinal_type dim)
123 : refSubcellPoints_(refSubcellPoints), paramPoints_(paramPoints), subcellMap_(subcellMap),
124 subcellOrd_(subcellOrd), dim_(dim){};
125
126 KOKKOS_INLINE_FUNCTION
127 void operator()(const size_type pt) const {
128
129 const auto u = paramPoints_(pt, 0);
130 const auto v = paramPoints_(pt, 1);
131
132 // map_dim(u,v) = c_0(dim) + c_1(dim)*u + c_2(dim)*v because both Quad and Tri ref faces are affine!
133 for (ordinal_type i=0;i<dim_;++i)
134 refSubcellPoints_(pt, i) = subcellMap_(subcellOrd_, i, 0) + ( subcellMap_(subcellOrd_, i, 1)*u +
135 subcellMap_(subcellOrd_, i, 2)*v );
136 }
137 };
138
139 template<typename refSubcellViewType,
140 typename paramPointsViewType,
141 typename subcellMapViewType>
143 refSubcellViewType refSubcellPoints_;
144 const paramPointsViewType paramPoints_;
145 const subcellMapViewType subcellMap_;
146 ordinal_type subcellOrd_, dim_;
147
148 KOKKOS_INLINE_FUNCTION
149 F_mapReferenceSubcell1( refSubcellViewType refSubcellPoints,
150 const paramPointsViewType paramPoints,
151 const subcellMapViewType subcellMap,
152 ordinal_type subcellOrd,
153 ordinal_type dim)
154 : refSubcellPoints_(refSubcellPoints), paramPoints_(paramPoints), subcellMap_(subcellMap),
155 subcellOrd_(subcellOrd), dim_(dim){};
156
157 KOKKOS_INLINE_FUNCTION
158 void operator()(const size_type pt) const {
159 const auto u = paramPoints_(pt, 0);
160 for (ordinal_type i=0;i<dim_;++i)
161 refSubcellPoints_(pt, i) = subcellMap_(subcellOrd_, i, 0) + ( subcellMap_(subcellOrd_, i, 1)*u );
162 }
163 };
164 }
165
166
167/*
168 template<typename DeviceType>
169 template<typename physPointValueType, class ...physPointProperties,
170 typename refPointValueType, class ...refPointProperties,
171 typename worksetCellValueType, class ...worksetCellProperties>
172 void
173 CellTools<DeviceType>::
174 mapToPhysicalFrame( Kokkos::DynRankView<physPointValueType,physPointProperties...> physPoints,
175 const Kokkos::DynRankView<refPointValueType,refPointProperties...> refPoints,
176 const Kokkos::DynRankView<worksetCellValueType,worksetCellProperties...> worksetCell,
177 const shards::CellTopology cellTopo ) {
178
179 auto basis = createHGradBasis<refPointValueType,refPointValueType>(cellTopo);
180 mapToPhysicalFrame(physPoints,
181 refPoints,
182 worksetCell,
183 basis);
184 }
185*/
186
187 template<typename DeviceType>
188 template<typename physPointValueType, class ...physPointProperties,
189 typename refPointValueType, class ...refPointProperties,
190 typename WorksetType,
191 typename HGradBasisPtrType>
192 void
194 mapToPhysicalFrame( Kokkos::DynRankView<physPointValueType,physPointProperties...> physPoints,
195 const Kokkos::DynRankView<refPointValueType,refPointProperties...> refPoints,
196 const WorksetType worksetCell,
197 const HGradBasisPtrType basis ) {
198#ifdef HAVE_INTREPID2_DEBUG
199 CellTools_mapToPhysicalFrameArgs( physPoints, refPoints, worksetCell, basis->getBaseCellTopology() );
200#endif
201 constexpr bool are_accessible =
202 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
203 typename decltype(physPoints)::memory_space>::accessible &&
204 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
205 typename decltype(refPoints)::memory_space>::accessible;
206
207 static_assert(are_accessible, "CellTools<DeviceType>::mapToPhysicalFrame(..): input/output views' memory spaces are not compatible with DeviceType");
208
209 const auto cellTopo = basis->getBaseCellTopology();
210 const auto numCells = worksetCell.extent(0);
211
212 //points can be rank-2 (P,D), or rank-3 (C,P,D)
213 const auto refPointRank = refPoints.rank();
214 const auto numPoints = (refPointRank == 2 ? refPoints.extent(0) : refPoints.extent(1));
215 const auto basisCardinality = basis->getCardinality();
216 auto vcprop = Kokkos::common_view_alloc_prop(physPoints);
217
218 using physPointViewType =Kokkos::DynRankView<physPointValueType,physPointProperties...>;
219 using valViewType = Kokkos::DynRankView<decltype(basis->getDummyOutputValue()),DeviceType>;
220
221 valViewType vals;
222
223 switch (refPointRank) {
224 case 2: {
225 // refPoints is (P,D): single set of ref. points is mapped to one or multiple physical cells
226 vals = valViewType(Kokkos::view_alloc("CellTools::mapToPhysicalFrame::vals", vcprop), basisCardinality, numPoints);
227 basis->getValues(vals,
228 refPoints,
229 OPERATOR_VALUE);
230 break;
231 }
232 case 3: {
233 // refPoints is (C,P,D): multiple sets of ref. points are mapped to matching number of physical cells.
234 //vals = valViewType("CellTools::mapToPhysicalFrame::vals", numCells, basisCardinality, numPoints);
235 vals = valViewType(Kokkos::view_alloc("CellTools::mapToPhysicalFrame::vals", vcprop), numCells, basisCardinality, numPoints);
236 for (size_type cell=0;cell<numCells;++cell)
237 basis->getValues(Kokkos::subdynrankview( vals, cell, Kokkos::ALL(), Kokkos::ALL() ),
238 Kokkos::subdynrankview( refPoints, cell, Kokkos::ALL(), Kokkos::ALL() ),
239 OPERATOR_VALUE);
240 break;
241 }
242 }
243
245
246 const auto loopSize = physPoints.extent(0)*physPoints.extent(1);
247 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
248 Kokkos::parallel_for( policy, FunctorType(physPoints, worksetCell, vals) );
249 }
250
251 template<typename DeviceType>
252 template<typename refSubcellPointValueType, class ...refSubcellPointProperties,
253 typename paramPointValueType, class ...paramPointProperties>
254 void
256 mapToReferenceSubcell( Kokkos::DynRankView<refSubcellPointValueType,refSubcellPointProperties...> refSubcellPoints,
257 const Kokkos::DynRankView<paramPointValueType,paramPointProperties...> paramPoints,
258 const ordinal_type subcellDim,
259 const ordinal_type subcellOrd,
260 const shards::CellTopology parentCell ) {
261 ordinal_type parentCellDim = parentCell.getDimension();
262#ifdef HAVE_INTREPID2_DEBUG
263 INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(parentCell), std::invalid_argument,
264 ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): the specified cell topology does not have a reference cell.");
265
266 INTREPID2_TEST_FOR_EXCEPTION( subcellDim < 1 ||
267 subcellDim > parentCellDim-1, std::invalid_argument,
268 ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): method defined only for subcells with dimension greater than 0 and less than the cell dimension");
269
270 INTREPID2_TEST_FOR_EXCEPTION( subcellOrd < 0 ||
271 subcellOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(subcellDim)), std::invalid_argument,
272 ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): subcell ordinal out of range.");
273
274 // refSubcellPoints is rank-2 (P,D1), D1 = cell dimension
275 INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.rank() != 2, std::invalid_argument,
276 ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints must have rank 2.");
277 INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.extent(1) != parentCell.getDimension(), std::invalid_argument,
278 ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints dimension (1) does not match to parent cell dimension.");
279
280 // paramPoints is rank-2 (P,D2) with D2 = subcell dimension
281 INTREPID2_TEST_FOR_EXCEPTION( paramPoints.rank() != 2, std::invalid_argument,
282 ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): paramPoints must have rank 2.");
283 INTREPID2_TEST_FOR_EXCEPTION( static_cast<ordinal_type>(paramPoints.extent(1)) != subcellDim, std::invalid_argument,
284 ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): paramPoints dimension (1) does not match to subcell dimension.");
285
286 // cross check: refSubcellPoints and paramPoints: dimension 0 must match
287 INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.extent(0) < paramPoints.extent(0), std::invalid_argument,
288 ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints dimension (0) does not match to paramPoints dimension(0).");
289#endif
290
291 constexpr bool are_accessible =
292 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
293 typename decltype(refSubcellPoints)::memory_space>::accessible &&
294 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
295 typename decltype(paramPoints)::memory_space>::accessible;
296
297 static_assert(are_accessible, "CellTools<DeviceType>::mapToReferenceSubcell(..): input/output views' memory spaces are not compatible with DeviceType");
298
299 // Get the subcell map, i.e., the coefficients of the parametrization function for the subcell
300 const auto subcellMap = RefSubcellParametrization<DeviceType>::get(subcellDim, parentCell.getKey());
301
302 const ordinal_type numPts = paramPoints.extent(0);
303
304 //Note: this function has several template parameters and the compiler gets confused if using a lambda function
305 using FunctorType1 = FunctorCellTools::F_mapReferenceSubcell1<decltype(refSubcellPoints), decltype(paramPoints), decltype(subcellMap)>;
306 using FunctorType2 = FunctorCellTools::F_mapReferenceSubcell2<decltype(refSubcellPoints), decltype(paramPoints), decltype(subcellMap)>;
307
308 Kokkos::RangePolicy<ExecSpaceType> policy(0, numPts);
309 // Apply the parametrization map to every point in parameter domain
310 switch (subcellDim) {
311 case 2: {
312 Kokkos::parallel_for( policy, FunctorType2(refSubcellPoints, paramPoints, subcellMap, subcellOrd, parentCellDim) );
313 break;
314 }
315 case 1: {
316 Kokkos::parallel_for( policy, FunctorType1(refSubcellPoints, paramPoints, subcellMap, subcellOrd, parentCellDim) );
317 break;
318 }
319 default: {}
320 }
321 }
322}
323
324#endif
void CellTools_mapToPhysicalFrameArgs(const physPointViewType physPoints, const refPointViewType refPoints, const worksetCellViewType worksetCell, const shards::CellTopology cellTopo)
Validates arguments to Intrepid2::CellTools::mapToPhysicalFrame.
static void mapToPhysicalFrame(Kokkos::DynRankView< physPointValueType, physPointProperties... > physPoints, const Kokkos::DynRankView< refPointValueType, refPointProperties... > refPoints, const WorksetType worksetCell, const HGradBasisPtrType basis)
Computes F, the reference-to-physical frame map.
static void mapToReferenceSubcell(Kokkos::DynRankView< refSubcellPointValueType, refSubcellPointProperties... > refSubcellPoints, const Kokkos::DynRankView< paramPointValueType, paramPointProperties... > paramPoints, const ordinal_type subcellDim, const ordinal_type subcellOrd, const shards::CellTopology parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
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...
Functor for mapping reference points to physical frame see Intrepid2::CellTools for more.