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 
57 namespace 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
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
Functor for mapping reference points to physical frame see Intrepid2::CellTools for more...
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...
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.
void CellTools_mapToPhysicalFrameArgs(const physPointViewType physPoints, const refPointViewType refPoints, const worksetCellViewType worksetCell, const shards::CellTopology cellTopo)
Validates arguments to Intrepid2::CellTools::mapToPhysicalFrame.