Intrepid2
Intrepid2_CellToolsDefInclusion.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_INCLUSION_HPP__
17 #define __INTREPID2_CELLTOOLS_DEF_INCLUSION_HPP__
18 
19 // disable clang warnings
20 #if defined (__clang__) && !defined (__INTEL_COMPILER)
21 #pragma clang system_header
22 #endif
23 
24 namespace Intrepid2 {
25 
26  //============================================================================================//
27  // //
28  // Inclusion tests //
29  // //
30  //============================================================================================//
31 
32 
33  template<typename DeviceType>
34  template<typename PointViewType>
35  bool
36  CellTools<DeviceType>::
37  checkPointInclusion( const PointViewType point,
38  const shards::CellTopology cellTopo,
39  const typename ScalarTraits<typename PointViewType::value_type>::scalar_type threshold) {
40 #ifdef HAVE_INTREPID2_DEBUG
41  INTREPID2_TEST_FOR_EXCEPTION( point.rank() != 1, std::invalid_argument,
42  ">>> ERROR (Intrepid2::CellTools::checkPointInclusion): Point must have rank 1. ");
43  INTREPID2_TEST_FOR_EXCEPTION( point.extent(0) != cellTopo.getDimension(), std::invalid_argument,
44  ">>> ERROR (Intrepid2::CellTools::checkPointInclusion): Point and cell dimensions do not match. ");
45 #endif
46  bool testResult = true;
47 
48  // A cell with extended topology has the same reference cell as a cell with base topology.
49  // => testing for inclusion in a reference Triangle<> and a reference Triangle<6> relies on
50  // on the same set of inequalities. To eliminate unnecessary cases we switch on the base topology
51  const auto key = cellTopo.getBaseKey();
52  switch (key) {
53 
54  case shards::Line<>::key :
55  testResult = PointInclusion<shards::Line<>::key>::check(point, threshold);
56  break;
57  case shards::Triangle<>::key :
58  testResult = PointInclusion<shards::Triangle<>::key>::check(point, threshold);
59  break;
60  case shards::Quadrilateral<>::key :
61  testResult = PointInclusion<shards::Quadrilateral<>::key>::check(point, threshold);
62  break;
63  case shards::Tetrahedron<>::key :
64  testResult = PointInclusion<shards::Tetrahedron<>::key>::check(point, threshold);
65  break;
66  case shards::Hexahedron<>::key :
67  testResult = PointInclusion<shards::Hexahedron<>::key>::check(point, threshold);
68  break;
69  case shards::Wedge<>::key :
70  testResult = PointInclusion<shards::Wedge<>::key>::check(point, threshold);
71  break;
72  case shards::Pyramid<>::key :
73  testResult = PointInclusion<shards::Pyramid<>::key>::check(point, threshold);
74  break;
75  default:
76  INTREPID2_TEST_FOR_EXCEPTION( !( (key == shards::Line<>::key ) ||
77  (key == shards::Triangle<>::key) ||
78  (key == shards::Quadrilateral<>::key) ||
79  (key == shards::Tetrahedron<>::key) ||
80  (key == shards::Hexahedron<>::key) ||
81  (key == shards::Wedge<>::key) ||
82  (key == shards::Pyramid<>::key) ),
83  std::invalid_argument,
84  ">>> ERROR (Intrepid2::CellTools::checkPointInclusion): Invalid cell topology. ");
85  }
86  return testResult;
87  }
88 
89 
90 
91  template<unsigned cellTopologyKey,
92  typename OutputViewType,
93  typename InputViewType>
94  struct checkPointInclusionFunctor {
95  OutputViewType output_;
96  InputViewType input_;
97  using ScalarType = typename ScalarTraits<typename InputViewType::value_type>::scalar_type;
98  ScalarType threshold_;
99 
100  KOKKOS_INLINE_FUNCTION
101  checkPointInclusionFunctor( OutputViewType output,
102  const InputViewType input,
103  const ScalarType threshold)
104  : output_(output),
105  input_(input),
106  threshold_(threshold) {}
107 
108  KOKKOS_INLINE_FUNCTION
109  void
110  operator()(const ordinal_type i) const {
111  const auto in = Kokkos::subview(input_,i,Kokkos::ALL());
112  const auto check = PointInclusion<cellTopologyKey>::check(in, threshold_);
113  output_(i) = check;
114  }
115 
116  KOKKOS_INLINE_FUNCTION
117  void
118  operator()(const ordinal_type i, const ordinal_type j) const {
119  const auto in = Kokkos::subview(input_,i,j,Kokkos::ALL());
120  const auto check = PointInclusion<cellTopologyKey>::check(in, threshold_);
121  output_(i,j) = check;
122  }
123  };
124 
125 
126  template<typename DeviceType>
127  template<unsigned cellTopologyKey,
128  typename OutputViewType,
129  typename InputViewType>
130  void CellTools<DeviceType>::
131  checkPointwiseInclusion( OutputViewType inCell,
132  const InputViewType points,
133  const typename ScalarTraits<typename InputViewType::value_type>::scalar_type threshold) {
134 
135  using FunctorType = checkPointInclusionFunctor<cellTopologyKey,decltype(inCell),decltype(points)>;
136  if (points.rank() == 2) { // inCell.rank() == 1
137  Kokkos::RangePolicy<typename DeviceType::execution_space> policy(0, points.extent(0));
138  Kokkos::parallel_for(policy, FunctorType(inCell, points, threshold));
139  } else { //points.rank() == 3, inCell.rank() == 2
140  Kokkos::MDRangePolicy<typename DeviceType::execution_space,Kokkos::Rank<2>> policy({0,0},{points.extent(0),points.extent(1)});
141  Kokkos::parallel_for(policy, FunctorType(inCell, points, threshold));
142  }
143  }
144 
145 
146  template<typename DeviceType>
147  template<typename InCellViewType,
148  typename InputViewType>
149  void
150  CellTools<DeviceType>::
151  checkPointwiseInclusion( InCellViewType inCell,
152  const InputViewType points,
153  const shards::CellTopology cellTopo,
154  const typename ScalarTraits<typename InputViewType::value_type>::scalar_type threshold ) {
155 #ifdef HAVE_INTREPID2_DEBUG
156  {
157  INTREPID2_TEST_FOR_EXCEPTION( (inCell.rank() != 1) && (inCell.rank() != 2), std::invalid_argument,
158  ">>> ERROR (Intrepid2::CellTools::checkPointwiseInclusion): InCell must have rank 1 or 2. ");
159  INTREPID2_TEST_FOR_EXCEPTION( points.extent(points.rank()-1) != cellTopo.getDimension(), std::invalid_argument,
160  ">>> ERROR (Intrepid2::CellTools::checkPointInclusion): Points and cell dimensions do not match. ");
161  INTREPID2_TEST_FOR_EXCEPTION( inCell.rank() != (points.rank()-1), std::invalid_argument,
162  ">>> ERROR (Intrepid2::CellTools::checkPointwiseInclusion): rank difference between inCell and points is 1.");
163  const ordinal_type iend = inCell.rank();
164 
165  for (ordinal_type i=0;i<iend;++i) {
166  INTREPID2_TEST_FOR_EXCEPTION( inCell.extent(i) != points.extent(i), std::invalid_argument,
167  ">>> ERROR (Intrepid2::CellTools::checkPointwiseInclusion): dimension mismatch between inCell and points. " );
168  }
169  }
170 #endif
171 
172  const auto key = cellTopo.getBaseKey();
173  switch (key) {
174 
175  case shards::Line<>::key :
176  checkPointwiseInclusion<shards::Line<>::key,decltype(inCell),decltype(points)>(inCell, points, threshold);
177  break;
178 
179  case shards::Triangle<>::key :
180  checkPointwiseInclusion<shards::Triangle<>::key,decltype(inCell),decltype(points)>(inCell, points, threshold);
181  break;
182 
183  case shards::Quadrilateral<>::key :
184  checkPointwiseInclusion<shards::Quadrilateral<>::key,decltype(inCell),decltype(points)>(inCell, points, threshold);
185  break;
186 
187  case shards::Tetrahedron<>::key :
188  checkPointwiseInclusion<shards::Tetrahedron<>::key,decltype(inCell),decltype(points)>(inCell, points, threshold);
189  break;
190 
191  case shards::Hexahedron<>::key :
192  checkPointwiseInclusion<shards::Hexahedron<>::key,decltype(inCell),decltype(points)>(inCell, points, threshold);
193  break;
194 
195  case shards::Wedge<>::key :
196  checkPointwiseInclusion<shards::Wedge<>::key,decltype(inCell),decltype(points)>(inCell, points, threshold);
197  break;
198 
199  case shards::Pyramid<>::key :
200  checkPointwiseInclusion<shards::Pyramid<>::key,decltype(inCell),decltype(points)>(inCell, points, threshold);
201  break;
202 
203  default:
204  INTREPID2_TEST_FOR_EXCEPTION( false,
205  std::invalid_argument,
206  ">>> ERROR (Intrepid2::CellTools::checkPointInclusion): Invalid cell topology. ");
207  }
208  }
209 
210  template<typename DeviceType>
211  template<typename inCellValueType, class ...inCellProperties,
212  typename pointValueType, class ...pointProperties,
213  typename cellWorksetValueType, class ...cellWorksetProperties>
214  void
215  CellTools<DeviceType>::
216  checkPointwiseInclusion( Kokkos::DynRankView<inCellValueType,inCellProperties...> inCell,
217  const Kokkos::DynRankView<pointValueType,pointProperties...> points,
218  const Kokkos::DynRankView<cellWorksetValueType,cellWorksetProperties...> cellWorkset,
219  const shards::CellTopology cellTopo,
220  const typename ScalarTraits<pointValueType>::scalar_type threshold ) {
221 #ifdef HAVE_INTREPID2_DEBUG
222  {
223  const auto key = cellTopo.getBaseKey();
224  INTREPID2_TEST_FOR_EXCEPTION( key != shards::Line<>::key &&
225  key != shards::Triangle<>::key &&
226  key != shards::Quadrilateral<>::key &&
227  key != shards::Tetrahedron<>::key &&
228  key != shards::Hexahedron<>::key &&
229  key != shards::Wedge<>::key &&
230  key != shards::Pyramid<>::key,
231  std::invalid_argument,
232  ">>> ERROR (Intrepid2::CellTools::checkPointwiseInclusion): cell topology not supported");
233  INTREPID2_TEST_FOR_EXCEPTION( inCell.rank() != 2, std::invalid_argument,
234  ">>> ERROR (Intrepid2::CellTools::checkPointwiseInclusion): InCell must have rank 2. ");
235  INTREPID2_TEST_FOR_EXCEPTION( cellWorkset.rank() != 3, std::invalid_argument,
236  ">>> ERROR (Intrepid2::CellTools::checkPointwiseInclusion): cellWorkset must have rank 3. ");
237  INTREPID2_TEST_FOR_EXCEPTION( cellWorkset.extent(2) != cellTopo.getDimension(), std::invalid_argument,
238  ">>> ERROR (Intrepid2::CellTools::checkPointInclusion): cellWorkset and points have incompatible dimensions. ");
239  INTREPID2_TEST_FOR_EXCEPTION( (points.rank() == 3) && (points.extent(0) != cellWorkset.extent(0)) , std::invalid_argument,
240  ">>> ERROR (Intrepid2::CellTools::checkPointInclusion): cellWorkset and points have incompatible dimensions. ");
241  }
242 #endif
243  const ordinal_type
244  numCells = cellWorkset.extent(0),
245  numPoints = points.extent(points.rank()-2),
246  spaceDim = cellTopo.getDimension();
247 
248  using result_layout = typename DeduceLayout< decltype(points) >::result_layout;
249  auto vcprop = Kokkos::common_view_alloc_prop(points);
250  using common_value_type = typename decltype(vcprop)::value_type;
251  Kokkos::DynRankView< common_value_type, result_layout, DeviceType > refPoints ( Kokkos::view_alloc("CellTools::checkPointwiseInclusion::refPoints", vcprop), numCells, numPoints, spaceDim);
252 
253  // expect refPoints(CPD), points (CPD or PD), cellWorkset(CND)
254  if(points.rank() == 3)
255  mapToReferenceFrame(refPoints, points, cellWorkset, cellTopo);
256  else { //points.rank() == 2
257  Kokkos::DynRankView< common_value_type, result_layout, DeviceType > cellPoints ( Kokkos::view_alloc("CellTools::checkPointwiseInclusion::physCellPoints", vcprop), numCells, numPoints, spaceDim);
258  RealSpaceTools<DeviceType>::clone(cellPoints,points);
259  mapToReferenceFrame(refPoints, cellPoints, cellWorkset, cellTopo);
260  }
261  checkPointwiseInclusion(inCell, refPoints, cellTopo, threshold);
262  }
263 
264 }
265 
266 #endif
Intrepid2::CellTools::checkPointwiseInclusion
static void checkPointwiseInclusion(OutputViewType inCell, const InputViewType points, const typename ScalarTraits< typename InputViewType::value_type >::scalar_type thresh=threshold< typename ScalarTraits< typename InputViewType::value_type >::scalar_type >())
Checks every point for inclusion in the reference cell of a given topology. The points can belong to ...
Definition: Intrepid2_CellToolsDefInclusion.hpp:131
Intrepid2::PointInclusion
This class implements a check function that determines whether a given point is inside or outside the...
Definition: Intrepid2_CellData.hpp:333
Intrepid2::RealSpaceTools::clone
static void clone(Kokkos::DynRankView< outputValueType, outputProperties...> output, const Kokkos::DynRankView< inputValueType, inputProperties...> input)
Clone input array.
Definition: Intrepid2_RealSpaceToolsDef.hpp:301
Intrepid2::ScalarTraits
scalar type traits
Definition: Intrepid2_Utils.hpp:120
Intrepid2::CellTools::checkPointInclusion
static bool checkPointInclusion(const PointViewType point, const shards::CellTopology cellTopo, const typename ScalarTraits< typename PointViewType::value_type >::scalar_type thres=threshold< typename ScalarTraits< typename PointViewType::value_type >::scalar_type >())
Checks if a point belongs to a reference cell.
Definition: Intrepid2_CellToolsDefInclusion.hpp:37
Generated by   doxygen 1.8.5