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>
72  struct F_mapToPhysicalFrame {
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 = rank(_basisVals);
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>
111  struct F_mapReferenceSubcell {
112  refSubcellViewType refSubcellPoints_;
113  const paramPointsViewType paramPoints_;
114  const subcellMapViewType subcellMap_;
115  const ordinal_type subcellOrd_;
116  const ordinal_type sideDim_;
117 
118  KOKKOS_INLINE_FUNCTION
119  F_mapReferenceSubcell( refSubcellViewType refSubcellPoints,
120  const paramPointsViewType paramPoints,
121  const subcellMapViewType subcellMap,
122  ordinal_type subcellOrd)
123  : refSubcellPoints_(refSubcellPoints), paramPoints_(paramPoints), subcellMap_(subcellMap),
124  subcellOrd_(subcellOrd), sideDim_(paramPoints_.extent(1)){};
125 
126  KOKKOS_INLINE_FUNCTION
127  void operator()(const size_type pt, const size_type d) const {
128 
129  refSubcellPoints_(pt, d) = subcellMap_(subcellOrd_, d, 0);
130  for(ordinal_type k=0; k<sideDim_; ++k)
131  refSubcellPoints_(pt, d) += subcellMap_(subcellOrd_, d, k+1)*paramPoints_(pt, k);
132  }
133  };
134 
135  template<typename refSubcellViewType,
136  typename paramPointsViewType,
137  typename subcellMapViewType,
138  typename ordViewType>
139  struct F_mapReferenceSubcellBatch {
140  refSubcellViewType refSubcellPoints_;
141  const paramPointsViewType paramPoints_;
142  const subcellMapViewType subcellMap_;
143  const ordViewType subcellOrd_;
144  const ordinal_type sideDim_;
145 
146  KOKKOS_INLINE_FUNCTION
147  F_mapReferenceSubcellBatch( refSubcellViewType refSubcellPoints,
148  const paramPointsViewType paramPoints,
149  const subcellMapViewType subcellMap,
150  ordViewType subcellOrd)
151  : refSubcellPoints_(refSubcellPoints), paramPoints_(paramPoints), subcellMap_(subcellMap),
152  subcellOrd_(subcellOrd), sideDim_(paramPoints_.extent(1)){};
153 
154  KOKKOS_INLINE_FUNCTION
155  void operator()(const size_type isc, const size_type pt, const size_type d) const {
156 
157  refSubcellPoints_(isc, pt, d) = subcellMap_(subcellOrd_(isc), d, 0);
158  for(ordinal_type k=0; k<sideDim_; ++k)
159  refSubcellPoints_(isc, pt, d) += subcellMap_(subcellOrd_(isc), d, k+1)*paramPoints_(pt, k);
160  }
161  };
162  }
163 
164  template<typename DeviceType>
165  template<typename PhysPointViewType,
166  typename RefPointViewType,
167  typename WorksetType,
168  typename HGradBasisPtrType>
169  void
170  CellTools<DeviceType>::
171  mapToPhysicalFrame( PhysPointViewType physPoints,
172  const RefPointViewType refPoints,
173  const WorksetType worksetCell,
174  const HGradBasisPtrType basis ) {
175 #ifdef HAVE_INTREPID2_DEBUG
176  CellTools_mapToPhysicalFrameArgs( physPoints, refPoints, worksetCell, basis->getBaseCellTopology() );
177 #endif
178  constexpr bool are_accessible =
179  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
180  typename decltype(physPoints)::memory_space>::accessible &&
181  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
182  typename decltype(refPoints)::memory_space>::accessible;
183 
184  static_assert(are_accessible, "CellTools<DeviceType>::mapToPhysicalFrame(..): input/output views' memory spaces are not compatible with DeviceType");
185 
186  const auto cellTopo = basis->getBaseCellTopology();
187  const auto numCells = worksetCell.extent(0);
188 
189  //points can be rank-2 (P,D), or rank-3 (C,P,D)
190  const auto refPointRank = refPoints.rank();
191  const auto numPoints = (refPointRank == 2 ? refPoints.extent(0) : refPoints.extent(1));
192  const auto basisCardinality = basis->getCardinality();
193  auto vcprop = Kokkos::common_view_alloc_prop(physPoints);
194 
195  using valViewType = Kokkos::DynRankView<decltype(basis->getDummyOutputValue()),DeviceType>;
196 
197  valViewType vals;
198 
199  switch (refPointRank) {
200  case 2: {
201  // refPoints is (P,D): single set of ref. points is mapped to one or multiple physical cells
202  vals = valViewType(Kokkos::view_alloc("CellTools::mapToPhysicalFrame::vals", vcprop), basisCardinality, numPoints);
203  basis->getValues(vals,
204  refPoints,
205  OPERATOR_VALUE);
206  break;
207  }
208  case 3: {
209  // refPoints is (C,P,D): multiple sets of ref. points are mapped to matching number of physical cells.
210  //vals = valViewType("CellTools::mapToPhysicalFrame::vals", numCells, basisCardinality, numPoints);
211  vals = valViewType(Kokkos::view_alloc("CellTools::mapToPhysicalFrame::vals", vcprop), numCells, basisCardinality, numPoints);
212  for (size_type cell=0;cell<numCells;++cell)
213  basis->getValues(Kokkos::subdynrankview( vals, cell, Kokkos::ALL(), Kokkos::ALL() ),
214  Kokkos::subdynrankview( refPoints, cell, Kokkos::ALL(), Kokkos::ALL() ),
215  OPERATOR_VALUE);
216  break;
217  }
218  }
219 
220  using FunctorType = FunctorCellTools::F_mapToPhysicalFrame<PhysPointViewType,WorksetType,valViewType>;
221 
222  const auto loopSize = physPoints.extent(0)*physPoints.extent(1);
223  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
224  Kokkos::parallel_for( policy, FunctorType(physPoints, worksetCell, vals) );
225  }
226 
227  template<typename DeviceType>
228  template<typename refSubcellViewType, typename paramViewType>
229  void
230  CellTools<DeviceType>::
231  mapToReferenceSubcell( refSubcellViewType refSubcellPoints,
232  const paramViewType paramPoints,
233  const ordinal_type subcellDim,
234  const ordinal_type subcellOrd,
235  const shards::CellTopology parentCell ) {
236 #ifdef HAVE_INTREPID2_DEBUG
237  ordinal_type parentCellDim = parentCell.getDimension();
238  INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(parentCell), std::invalid_argument,
239  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): the specified cell topology does not have a reference cell.");
240 
241  INTREPID2_TEST_FOR_EXCEPTION( subcellDim < 1 ||
242  subcellDim > parentCellDim-1, std::invalid_argument,
243  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): method defined only for subcells with dimension greater than 0 and less than the cell dimension");
244 
245  INTREPID2_TEST_FOR_EXCEPTION( subcellOrd < 0 ||
246  subcellOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(subcellDim)), std::invalid_argument,
247  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): subcell ordinal out of range.");
248 
249  // refSubcellPoints is rank-2 (P,D1), D1 = cell dimension
250  INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.rank() != 2, std::invalid_argument,
251  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints must have rank 2.");
252  INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.extent(1) != parentCell.getDimension(), std::invalid_argument,
253  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints dimension (1) does not match to parent cell dimension.");
254 
255  // paramPoints is rank-2 (P,D2) with D2 = subcell dimension
256  INTREPID2_TEST_FOR_EXCEPTION( paramPoints.rank() != 2, std::invalid_argument,
257  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): paramPoints must have rank 2.");
258  INTREPID2_TEST_FOR_EXCEPTION( static_cast<ordinal_type>(paramPoints.extent(1)) != subcellDim, std::invalid_argument,
259  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): paramPoints dimension (1) does not match to subcell dimension.");
260 
261  // cross check: refSubcellPoints and paramPoints: dimension 0 must match
262  INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.extent(0) != paramPoints.extent(0), std::invalid_argument,
263  ">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints dimension (0) does not match to paramPoints dimension(0).");
264 #endif
265 
266  // Get the subcell map, i.e., the coefficients of the parametrization function for the subcell
267  const auto subcellMap = RefSubcellParametrization<DeviceType>::get(subcellDim, parentCell.getKey());
268 
269  mapToReferenceSubcell(refSubcellPoints, paramPoints, subcellMap, subcellOrd);
270  }
271 
272  template<typename DeviceType>
273  template<typename refSubcellViewType, typename paramViewType>
274  void
275  CellTools<DeviceType>::
276  mapToReferenceSubcell( refSubcellViewType refSubcellPoints,
277  const paramViewType paramPoints,
278  const typename RefSubcellParametrization<DeviceType>::ConstViewType subcellParametrization,
279  const ordinal_type subcellOrd) {
280 
281  constexpr bool are_accessible =
282  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
283  typename decltype(refSubcellPoints)::memory_space>::accessible &&
284  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
285  typename decltype(paramPoints)::memory_space>::accessible;
286  static_assert(are_accessible, "CellTools<DeviceType>::mapToReferenceSubcell(..): input/output views' memory spaces are not compatible with DeviceType");
287 
288 
289  #ifdef HAVE_INTREPID2_DEBUG
290  const bool ranks_and_dims_compatible = (refSubcellPoints.rank() == 2) && (paramPoints.rank() == 2) && (subcellParametrization.rank() == 3) &&
291  (refSubcellPoints.extent(0) == paramPoints.extent(0)) &&
292  (refSubcellPoints.extent(1) == subcellParametrization.extent(1)) &&
293  (paramPoints.extent(1) == subcellParametrization.extent(2)-1);
294 
295  INTREPID2_TEST_FOR_EXCEPTION(!ranks_and_dims_compatible, std::invalid_argument, "CellTools<DeviceType>::mapToReferenceSubcell(..): input/output views' ranks and dimensions are not compatible");
296  #endif
297 
298  const ordinal_type parentCellDim = subcellParametrization.extent(1);
299  const ordinal_type numPts = paramPoints.extent(0);
300 
301  //Note: this function has several template parameters and the compiler gets confused if using a lambda function
302  using FunctorType = FunctorCellTools::F_mapReferenceSubcell<decltype(refSubcellPoints), decltype(paramPoints), decltype(subcellParametrization)>;
303 
304  Kokkos::MDRangePolicy<ExecSpaceType,Kokkos::Rank<2>> rangePolicy({0,0},{numPts,parentCellDim});
305 
306  // Apply the parametrization map to every point in parameter domain
307  Kokkos::parallel_for( rangePolicy, FunctorType(refSubcellPoints, paramPoints, subcellParametrization, subcellOrd) );
308  }
309 
310  template<typename DeviceType>
311  template<typename refSubcellViewType, typename paramViewType, typename ordViewType>
312  void
313  CellTools<DeviceType>::
314  mapToReferenceSubcell( refSubcellViewType refSubcellPoints,
315  const paramViewType paramPoints,
316  const typename RefSubcellParametrization<DeviceType>::ConstViewType subcellParametrization,
317  const ordViewType subcellOrd) {
318 
319  constexpr bool are_accessible =
320  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
321  typename decltype(refSubcellPoints)::memory_space>::accessible &&
322  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
323  typename decltype(paramPoints)::memory_space>::accessible &&
324  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
325  typename decltype(subcellOrd)::memory_space>::accessible;
326  static_assert(are_accessible, "CellTools<DeviceType>::mapToReferenceSubcell(..): input/output views' memory spaces are not compatible with DeviceType");
327 
328 
329 #ifdef HAVE_INTREPID2_DEBUG
330  const bool ranks_and_dims_compatible = (refSubcellPoints.rank() == 3) && (paramPoints.rank() == 2) && (subcellParametrization.rank() == 3) &&
331  (refSubcellPoints.extent(0) == subcellOrd.extent(0)) &&
332  (refSubcellPoints.extent(1) == paramPoints.extent(0)) &&
333  (refSubcellPoints.extent(2) == subcellParametrization.extent(1)) &&
334  (paramPoints.extent(1) == subcellParametrization.extent(2)-1);
335 
336  INTREPID2_TEST_FOR_EXCEPTION(!ranks_and_dims_compatible, std::invalid_argument, "CellTools<DeviceType>::mapToReferenceSubcell(..): input/output views' ranks and dimensions are not compatible");
337 #endif
338 
339  const ordinal_type numSubCells = refSubcellPoints.extent(0);
340  const ordinal_type parentCellDim = subcellParametrization.extent(1);
341  const ordinal_type numPts = paramPoints.extent(0);
342 
343 
344  //Note: this function has several template parameters and the compiler gets confused if using a lambda function
345  using FunctorType = FunctorCellTools::F_mapReferenceSubcellBatch<decltype(refSubcellPoints), decltype(paramPoints), decltype(subcellParametrization), decltype(subcellOrd)>;
346 
347  Kokkos::MDRangePolicy<ExecSpaceType,Kokkos::Rank<3>> rangePolicy({0,0,0},{numSubCells,numPts,parentCellDim});
348 
349  // Apply the parametrization map to every point in parameter domain
350  Kokkos::parallel_for( rangePolicy, FunctorType(refSubcellPoints, paramPoints, subcellParametrization, subcellOrd) );
351  }
352 
353 }
354 
355 #endif
Intrepid2::FunctorCellTools::F_mapToPhysicalFrame
Functor for mapping reference points to physical frame see Intrepid2::CellTools for more...
Definition: Intrepid2_CellToolsDefRefToPhys.hpp:72
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:103
Intrepid2::FunctorCellTools::F_mapReferenceSubcell
Definition: Intrepid2_CellToolsDefRefToPhys.hpp:111
Intrepid2::FunctorCellTools::F_mapReferenceSubcellBatch
Definition: Intrepid2_CellToolsDefRefToPhys.hpp:139
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::CellTools::mapToPhysicalFrame
static void mapToPhysicalFrame(PhysPointValueType physPoints, const RefPointValueType refPoints, const WorksetType worksetCell, const HGradBasisPtrType basis)
Computes F, the reference-to-physical frame map.
Generated by   doxygen 1.8.5