Intrepid2
Intrepid2_CellToolsDefInclusion.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_INCLUSION_HPP__
50#define __INTREPID2_CELLTOOLS_DEF_INCLUSION_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 // Inclusion tests //
62 // //
63 //============================================================================================//
64
65
66 template<typename DeviceType>
67 template<typename pointValueType, class ...pointProperties>
68 bool
70 checkPointInclusion( const Kokkos::DynRankView<pointValueType,pointProperties...> point,
71 const shards::CellTopology cellTopo,
72 const double threshold) {
73#ifdef HAVE_INTREPID2_DEBUG
74 INTREPID2_TEST_FOR_EXCEPTION( point.rank() != 1, std::invalid_argument,
75 ">>> ERROR (Intrepid2::CellTools::checkPointInclusion): Point must have rank 1. ");
76 INTREPID2_TEST_FOR_EXCEPTION( point.extent(0) != cellTopo.getDimension(), std::invalid_argument,
77 ">>> ERROR (Intrepid2::CellTools::checkPointInclusion): Point and cell dimensions do not match. ");
78#endif
79 bool testResult = true;
80
81 const double t = threshold; //(threshold < 0 ? threshold() : threshold);
82
83 // Using these values in the tests effectievly inflates the reference element to a larger one
84 const double minus_one = -1.0 - t;
85 const double plus_one = 1.0 + t;
86 const double minus_zero = -t;
87
88 // A cell with extended topology has the same reference cell as a cell with base topology.
89 // => testing for inclusion in a reference Triangle<> and a reference Triangle<6> relies on
90 // on the same set of inequalities. To eliminate unnecessary cases we switch on the base topology
91 const auto key = cellTopo.getBaseKey();
92 switch (key) {
93
94 case shards::Line<>::key :
95 if( !(minus_one <= point(0) && point(0) <= plus_one)) testResult = false;
96 break;
97
98 case shards::Triangle<>::key : {
99 const auto distance = Util<pointValueType>::max( std::max( -point(0), -point(1) ), point(0) + point(1) - 1.0 );
100 if( distance > threshold ) testResult = false;
101 break;
102 }
103
104 case shards::Quadrilateral<>::key :
105 if(!( (minus_one <= point(0) && point(0) <= plus_one) &&
106 (minus_one <= point(1) && point(1) <= plus_one) ) ) testResult = false;
107 break;
108
109 case shards::Tetrahedron<>::key : {
110 const auto distance = Util<pointValueType>::max( Util<pointValueType>::max(-point(0),-point(1)),
111 Util<pointValueType>::max(-point(2), point(0) + point(1) + point(2) - 1) );
112 if( distance > threshold ) testResult = false;
113 break;
114 }
115
116 case shards::Hexahedron<>::key :
117 if(!((minus_one <= point(0) && point(0) <= plus_one) &&
118 (minus_one <= point(1) && point(1) <= plus_one) &&
119 (minus_one <= point(2) && point(2) <= plus_one)))
120 testResult = false;
121 break;
122
123 // The base of the reference prism is the same as the reference triangle => apply triangle test
124 // to X and Y coordinates and test whether Z is in [-1,1]
125 case shards::Wedge<>::key : {
126 const auto distance = Util<pointValueType>::max( Util<pointValueType>::max( -point(0), -point(1) ), point(0) + point(1) - 1 );
127 if( distance > threshold ||
128 point(2) < minus_one || point(2) > plus_one)
129 testResult = false;
130 break;
131 }
132
133 // The base of the reference pyramid is the same as the reference quad cell => a horizontal plane
134 // through a point P(x,y,z) intersects the pyramid at a quadrilateral that equals the base quad
135 // scaled by (1-z) => P(x,y,z) is inside the pyramid <=> (x,y) is in [-1+z,1-z]^2 && 0 <= Z <= 1
136 case shards::Pyramid<>::key : {
137 const auto left = minus_one + point(2);
138 const auto right = plus_one - point(2);
139 if(!( (left <= point(0) && point(0) <= right) && \
140 (left <= point(1) && point(1) <= right) &&
141 (minus_zero <= point(2) && point(2) <= plus_one) ) ) \
142 testResult = false;
143 break;
144 }
145
146 default:
147 INTREPID2_TEST_FOR_EXCEPTION( !( (key == shards::Line<>::key ) ||
148 (key == shards::Triangle<>::key) ||
149 (key == shards::Quadrilateral<>::key) ||
150 (key == shards::Tetrahedron<>::key) ||
151 (key == shards::Hexahedron<>::key) ||
152 (key == shards::Wedge<>::key) ||
153 (key == shards::Pyramid<>::key) ),
154 std::invalid_argument,
155 ">>> ERROR (Intrepid2::CellTools::checkPointInclusion): Invalid cell topology. ");
156 }
157 return testResult;
158 }
159
160
161
162// template<class Scalar>
163// template<class ArrayPoint>
164// ordinal_type CellTools<Scalar>::checkPointsetInclusion(const ArrayPoint& points,
165// const shards::CellTopology & cellTopo,
166// const double & threshold) {
167
168// ordinal_type rank = points.rank();
169
170// #ifdef HAVE_INTREPID2_DEBUG
171// INTREPID2_TEST_FOR_EXCEPTION( !( (1 <=getrank(points) ) && (getrank(points) <= 3) ), std::invalid_argument,
172// ">>> ERROR (Intrepid2::CellTools::checkPointsetInclusion): rank-1, 2 or 3 required for input points array. ");
173
174// // The last dimension of points array at (rank - 1) is the spatial dimension. Must equal the cell dimension.
175// INTREPID2_TEST_FOR_EXCEPTION( !((index_type) points.extent(rank - 1) == (index_type)cellTopo.getDimension() ), std::invalid_argument,
176// ">>> ERROR (Intrepid2::CellTools::checkPointsetInclusion): Point and cell dimensions do not match. ");
177// #endif
178
179// // create temp output array depending on the rank of the input array
180// FieldContainer<ordinal_type> inRefCell;
181// index_type dim0(0), dim1(0);
182// switch(rank) {
183// case 1:
184// inRefCell.resize(1);
185// break;
186// case 2:
187// dim0 = static_cast<index_type>(points.extent(0));
188// inRefCell.resize(dim0);
189// break;
190// case 3:
191// dim0 = static_cast<index_type>(points.extent(0));
192// dim1 = static_cast<index_type>(points.extent(1));
193// inRefCell.resize(dim0, dim1);
194// break;
195// }
196
197// // Call the inclusion method which returns inclusion results for all points
198// checkPointwiseInclusion(inRefCell, points, cellTopo, threshold);
199
200// // Check if any points were outside, return 0 after finding one
201
202// switch(rank) {
203// case 1:
204// if (inRefCell(0) == 0)
205// return 0;
206// break;
207// case 2:
208// for(index_type i = 0; i < dim0; i++ )
209// if (inRefCell(i) == 0)
210// return 0;
211// break;
212
213// case 3:
214// for(index_type i = 0; i < dim0; i++ )
215// for(index_type j = 0; j < dim1; j++ )
216// if (inRefCell(i,j) == 0)
217// return 0;
218// break;
219// }
220
221// return 1; //all points are inside
222// }
223
224
225 template<typename DeviceType>
226 template<typename inCellValueType, class ...inCellProperties,
227 typename pointValueType, class ...pointProperties>
228 void
230 checkPointwiseInclusion( Kokkos::DynRankView<inCellValueType,inCellProperties...> inCell,
231 const Kokkos::DynRankView<pointValueType,pointProperties...> points,
232 const shards::CellTopology cellTopo,
233 const double threshold ) {
234#ifdef HAVE_INTREPID2_DEBUG
235 {
236 INTREPID2_TEST_FOR_EXCEPTION( inCell.rank() != (points.rank()-1), std::invalid_argument,
237 ">>> ERROR (Intrepid2::CellTools::checkPointwiseInclusion): rank difference between inCell and points is 1.");
238 const ordinal_type iend = inCell.rank();
239 for (ordinal_type i=0;i<iend;++i) {
240 INTREPID2_TEST_FOR_EXCEPTION( inCell.extent(i) != points.extent(i), std::invalid_argument,
241 ">>> ERROR (Intrepid2::CellTools::checkPointwiseInclusion): dimension mismatch between inCell and points.");
242 }
243 }
244#endif
245
246 // do we really need to support 3 ranks ?
247 switch (points.rank()) {
248 case 2: {
249 const ordinal_type iend = points.extent(0);
250 for (ordinal_type i=0;i<iend;++i) {
251 const auto point = Kokkos::subview(points, i, Kokkos::ALL());
252 inCell(i) = checkPointInclusion(point, cellTopo, threshold);
253 }
254 break;
255 }
256 case 3: {
257 const ordinal_type
258 iend = points.extent(0),
259 jend = points.extent(1);
260 for (ordinal_type i=0;i<iend;++i)
261 for (ordinal_type j=0;j<jend;++j) {
262 const auto point = Kokkos::subview(points, i, j, Kokkos::ALL());
263 inCell(i, j) = checkPointInclusion(point, cellTopo, threshold);
264 }
265 break;
266 }
267 }
268 }
269
270 template<typename DeviceType>
271 template<typename inCellValueType, class ...inCellProperties,
272 typename pointValueType, class ...pointProperties,
273 typename cellWorksetValueType, class ...cellWorksetProperties>
274 void
276 checkPointwiseInclusion( Kokkos::DynRankView<inCellValueType,inCellProperties...> inCell,
277 const Kokkos::DynRankView<pointValueType,pointProperties...> points,
278 const Kokkos::DynRankView<cellWorksetValueType,cellWorksetProperties...> cellWorkset,
279 const shards::CellTopology cellTopo,
280 const double threshold ) {
281#ifdef HAVE_INTREPID2_DEBUG
282 {
283 const auto key = cellTopo.getBaseKey();
284 INTREPID2_TEST_FOR_EXCEPTION( key != shards::Line<>::key &&
285 key != shards::Triangle<>::key &&
286 key != shards::Quadrilateral<>::key &&
287 key != shards::Tetrahedron<>::key &&
288 key != shards::Hexahedron<>::key &&
289 key != shards::Wedge<>::key &&
290 key != shards::Pyramid<>::key,
291 std::invalid_argument,
292 ">>> ERROR (Intrepid2::CellTools::checkPointwiseInclusion): cell topology not supported");
293
294 INTREPID2_TEST_FOR_EXCEPTION( points.rank() != 3, std::invalid_argument,
295 ">>> ERROR (Intrepid2::CellTools::checkPointwiseInclusion): Points must have rank 3. ");
296 INTREPID2_TEST_FOR_EXCEPTION( cellWorkset.rank() != 3, std::invalid_argument,
297 ">>> ERROR (Intrepid2::CellTools::checkPointwiseInclusion): cellWorkset must have rank 3. ");
298 INTREPID2_TEST_FOR_EXCEPTION( points.extent(2) != cellTopo.getDimension(), std::invalid_argument,
299 ">>> ERROR (Intrepid2::CellTools::checkPointInclusion): Points and cell dimensions do not match. ");
300 INTREPID2_TEST_FOR_EXCEPTION( cellWorkset.extent(2) != cellTopo.getDimension(), std::invalid_argument,
301 ">>> ERROR (Intrepid2::CellTools::checkPointInclusion): cellWorkset and cell dimensions do not match. ");
302 INTREPID2_TEST_FOR_EXCEPTION( points.extent(0) != cellWorkset.extent(0) , std::invalid_argument,
303 ">>> ERROR (Intrepid2::CellTools::checkPointInclusion): cellWorkset and points dimension(0) does not match. ");
304 }
305#endif
306 const ordinal_type
307 numCells = cellWorkset.extent(0),
308 numPoints = points.extent(1),
309 spaceDim = cellTopo.getDimension();
310
311 using result_layout = typename DeduceLayout< decltype(points) >::result_layout;
312 auto vcprop = Kokkos::common_view_alloc_prop(points);
313 using common_value_type = typename decltype(vcprop)::value_type;
314 Kokkos::DynRankView< common_value_type, result_layout, DeviceType > refPoints ( Kokkos::view_alloc("CellTools::checkPointwiseInclusion::refPoints", vcprop), numCells, numPoints, spaceDim);
315
316 // expect refPoints(CPD), points(CPD), cellWorkset(CND)
317 mapToReferenceFrame(refPoints, points, cellWorkset, cellTopo);
318 checkPointwiseInclusion(inCell, refPoints, cellTopo, threshold);
319 }
320
321}
322
323#endif
static bool checkPointInclusion(const Kokkos::DynRankView< pointValueType, pointProperties... > point, const shards::CellTopology cellTopo, const double thres=threshold())
Checks if a point belongs to a reference cell.
static void checkPointwiseInclusion(Kokkos::DynRankView< inCellValueType, inCellProperties... > inCell, const Kokkos::DynRankView< pointValueType, pointProperties... > points, const shards::CellTopology cellTopo, const double thres=threshold())
Checks every point in a set for inclusion in a reference cell.
small utility functions
layout deduction (temporary meta-function)