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