Intrepid2
Intrepid2_OrientationToolsDefModifyPoints.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
48#ifndef __INTREPID2_ORIENTATIONTOOLS_DEF_MODIFY_POINTS_HPP__
49#define __INTREPID2_ORIENTATIONTOOLS_DEF_MODIFY_POINTS_HPP__
50
51// disable clang warnings
52#if defined (__clang__) && !defined (__INTEL_COMPILER)
53#pragma clang system_header
54#endif
55
56namespace Intrepid2 {
57
58 namespace Impl {
59
60 // ------------------------------------------------------------------------------------
61 // Modified points according to orientations
62 //
63 //
64 template<typename VT>
65 KOKKOS_INLINE_FUNCTION
66 void
69 const VT pt,
70 const ordinal_type ort) {
71#ifdef HAVE_INTREPID2_DEBUG
72 INTREPID2_TEST_FOR_ABORT( !( -1.0 <= pt && pt <= 1.0 ),
73 ">>> ERROR (Intrepid::OrientationTools::getModifiedLinePoint): "
74 "Input point is out of range [-1, 1].");
75#endif
76
77 switch (ort) {
78 case 0: ot = pt; break;
79 case 1: ot = -pt; break;
80 default:
81 INTREPID2_TEST_FOR_ABORT( true,
82 ">>> ERROR (Intrepid2::OrientationTools::getModifiedLinePoint): "
83 "Orientation is invalid (0--1)." );
84 }
85 }
86
87 template<typename JacobianViewType>
88 KOKKOS_INLINE_FUNCTION
89 void
91 getLineJacobian(JacobianViewType jacobian, const ordinal_type ort) {
92#ifdef HAVE_INTREPID2_DEBUG
93 INTREPID2_TEST_FOR_ABORT( ort <0 || ort >1,
94 ">>> ERROR (Intrepid2::OrientationTools::getLineJacobian): " \
95 "Orientation is invalid (0--1)." );
96#endif
97
98 ordinal_type jac[2] = { 1, -1 };
99
100 jacobian(0,0) = jac[ort];
101 }
102
103 template<typename VT>
105 void
108 VT &ot1,
109 const VT pt0,
110 const VT pt1,
111 const ordinal_type ort) {
112 const VT lambda[3] = { 1.0 - pt0 - pt1,
113 pt0,
114 pt1 };
115
116#ifdef HAVE_INTREPID2_DEBUG
117 // hard-coded epsilon to avoid CUDA issue with calling host code. Would be nice to have a portable way to define this, but probably not worth the effort at the moment.
118 VT eps = 1e-14; // = 10.0*std::numeric_limits<VT>::epsilon();
119 INTREPID2_TEST_FOR_ABORT( !( -eps <= lambda[0] && lambda[0] <= 1.0+eps ),
120 ">>> ERROR (Intrepid::OrientationTools::getModifiedTrianglePoint): " \
121 "Computed bicentric coordinate (lamba[0]) is out of range [0, 1].");
122
123 INTREPID2_TEST_FOR_ABORT( !( -eps <= lambda[1] && lambda[1] <= 1.0+eps ),
124 ">>> ERROR (Intrepid::OrientationTools::getModifiedTrianglePoint): " \
125 "Computed bicentric coordinate (lamba[1]) is out of range [0, 1].");
126
127 INTREPID2_TEST_FOR_ABORT( !( -eps <= lambda[2] && lambda[2] <= 1.0+eps ),
128 ">>> ERROR (Intrepid::OrientationTools::getModifiedTrianglePoint): "
129 "Computed bicentric coordinate (lamba[2]) is out of range [0, 1].");
130#endif
131
132 switch (ort) {
133 case 0: ot0 = lambda[1]; ot1 = lambda[2]; break;
134 case 1: ot0 = lambda[0]; ot1 = lambda[1]; break;
135 case 2: ot0 = lambda[2]; ot1 = lambda[0]; break;
136
137 case 3: ot0 = lambda[2]; ot1 = lambda[1]; break;
138 case 4: ot0 = lambda[0]; ot1 = lambda[2]; break;
139 case 5: ot0 = lambda[1]; ot1 = lambda[0]; break;
140 default:
141 INTREPID2_TEST_FOR_ABORT( true,
142 ">>> ERROR (Intrepid2::OrientationTools::getModifiedTrianglePoint): " \
143 "Orientation is invalid (0--5)." );
144 }
145 }
146
147 template<typename JacobianViewType>
148 KOKKOS_INLINE_FUNCTION
149 void
152#ifdef HAVE_INTREPID2_DEBUG
153 INTREPID2_TEST_FOR_ABORT( ort <0 || ort >5,
154 ">>> ERROR (Intrepid2::OrientationTools::getTriangleJacobian): " \
155 "Orientation is invalid (0--5)." );
156#endif
157
158 ordinal_type jac[6][2][2] = { { { 1, 0 },
159 { 0, 1 } }, // 0
160 { { -1, -1 },
161 { 1, 0 } }, // 1
162 { { 0, 1 },
163 { -1, -1 } }, // 2
164 { { 0, 1 },
165 { 1, 0 } }, // 3
166 { { -1, -1 },
167 { 0, 1 } }, // 4
168 { { 1, 0 },
169 { -1, -1 } } }; // 5
170
171 jacobian(0,0) = jac[ort][0][0];
172 jacobian(0,1) = jac[ort][0][1];
173 jacobian(1,0) = jac[ort][1][0];
174 jacobian(1,1) = jac[ort][1][1];
175 }
176
177 template<typename VT>
179 void
182 VT &ot1,
183 const VT pt0,
184 const VT pt1,
185 const ordinal_type ort) {
186#ifdef HAVE_INTREPID2_DEBUG
187 INTREPID2_TEST_FOR_ABORT( !( -1.0 <= pt0 && pt0 <= 1.0 ),
188 ">>> ERROR (Intrepid::OrientationTools::getModifiedQuadrilateralPoint): " \
189 "Input point(0) is out of range [-1, 1].");
190
191 INTREPID2_TEST_FOR_ABORT( !( -1.0 <= pt1 && pt1 <= 1.0 ),
192 ">>> ERROR (Intrepid::OrientationTools::getModifiedQuadrilateralPoint): " \
193 "Input point(1) is out of range [-1, 1].");
194#endif
195
196 const VT lambda[2][2] = { { pt0, -pt0 },
197 { pt1, -pt1 } };
198
199 switch (ort) {
200 case 0: ot0 = lambda[0][0]; ot1 = lambda[1][0]; break;
201 case 1: ot0 = lambda[1][1]; ot1 = lambda[0][0]; break;
202 case 2: ot0 = lambda[0][1]; ot1 = lambda[1][1]; break;
203 case 3: ot0 = lambda[1][0]; ot1 = lambda[0][1]; break;
204 // flip
205 case 4: ot0 = lambda[1][0]; ot1 = lambda[0][0]; break;
206 case 5: ot0 = lambda[0][1]; ot1 = lambda[1][0]; break;
207 case 6: ot0 = lambda[1][1]; ot1 = lambda[0][1]; break;
208 case 7: ot0 = lambda[0][0]; ot1 = lambda[1][1]; break;
209 default:
210 INTREPID2_TEST_FOR_ABORT( true,
211 ">>> ERROR (Intrepid2::OrientationTools::getModifiedQuadrilateralPoint): " \
212 "Orientation is invalid (0--7)." );
213 }
214 }
215
216 template<typename JacobianViewType>
217 KOKKOS_INLINE_FUNCTION
218 void
221#ifdef HAVE_INTREPID2_DEBUG
222 INTREPID2_TEST_FOR_ABORT( ort <0 || ort >7,
223 ">>> ERROR (Intrepid2::OrientationTools::getQuadrilateralJacobian): " \
224 "Orientation is invalid (0--7)." );
225#endif
226
227 ordinal_type jac[8][2][2] = { { { 1, 0 },
228 { 0, 1 } }, // 0
229 { { 0, -1 },
230 { 1, 0 } }, // 1
231 { { -1, 0 },
232 { 0, -1 } }, // 2
233 { { 0, 1 },
234 { -1, 0 } }, // 3
235 { { 0, 1 },
236 { 1, 0 } }, // 4
237 { { -1, 0 },
238 { 0, 1 } }, // 5
239 { { 0, -1 },
240 { -1, 0 } }, // 6
241 { { 1, 0 },
242 { 0, -1 } } }; // 7
243
244 jacobian(0,0) = jac[ort][0][0];
245 jacobian(0,1) = jac[ort][0][1];
246 jacobian(1,0) = jac[ort][1][0];
247 jacobian(1,1) = jac[ort][1][1];
248 }
249
250 template<typename outPointViewType,
251 typename refPointViewType>
252 inline
253 void
257 const shards::CellTopology cellTopo,
258 const ordinal_type cellOrt) {
259#ifdef HAVE_INTREPID2_DEBUG
260 {
261 const auto cellDim = cellTopo.getDimension();
262 INTREPID2_TEST_FOR_EXCEPTION( !( (1 <= cellDim) && (cellDim <= 2 ) ), std::invalid_argument,
263 ">>> ERROR (Intrepid::OrientationTools::mapToModifiedReference): " \
264 "Method defined only for 1 and 2-dimensional subcells.");
265
266 INTREPID2_TEST_FOR_EXCEPTION( !( outPoints.extent(0) == refPoints.extent(0) ), std::invalid_argument,
267 ">>> ERROR (Intrepid::OrientationTools::mapToModifiedReference): " \
268 "Size of input and output point arrays does not match each other.");
269 }
270#endif
271
273 }
274
275
276 template<typename outPointViewType,
277 typename refPointViewType>
279 void
283 const unsigned cellTopoKey,
284 const ordinal_type cellOrt) {
285 // Apply the parametrization map to every point in parameter domain
286 const ordinal_type numPts = outPoints.extent(0);
287 switch (cellTopoKey) {
288 case shards::Line<>::key : {
289 for (ordinal_type pt=0;pt<numPts;++pt)
291 refPoints(pt, 0),
292 cellOrt);
293 break;
294 }
295 case shards::Triangle<>::key : {
296 for (ordinal_type pt=0;pt<numPts;++pt)
298 refPoints(pt, 0), refPoints(pt, 1),
299 cellOrt);
300 break;
301 }
302 case shards::Quadrilateral<>::key : {
303 for (ordinal_type pt=0;pt<numPts;++pt)
305 refPoints(pt, 0), refPoints(pt, 1),
306 cellOrt);
307 break;
308 }
309 default: {
310 INTREPID2_TEST_FOR_ABORT( true,
311 ">>> ERROR (Intrepid2::OrientationTools::mapToModifiedReference): " \
312 "Invalid cell topology key." );
313 break;
314 }
315 }
316 }
317
318
319 template<typename outPointViewType>
320 inline
321 void
324 const shards::CellTopology cellTopo,
325 const ordinal_type cellOrt) {
326#ifdef HAVE_INTREPID2_DEBUG
327 {
328 const auto cellDim = cellTopo.getDimension();
329 INTREPID2_TEST_FOR_EXCEPTION( !( (1 <= cellDim) && (cellDim <= 2 ) ), std::invalid_argument,
330 ">>> ERROR (Intrepid::OrientationTools::getJacobianOfOrientationMap): " \
331 "Method defined only for 1 and 2-dimensional subcells.");
332
333 INTREPID2_TEST_FOR_ABORT( jacobian.rank() != 2,
334 ">>> ERROR (Intrepid2::OrientationTools::getJacobianOfOrientationMap): " \
335 "Jacobian should have rank 2" );
336
337 INTREPID2_TEST_FOR_EXCEPTION( ((jacobian.extent(0) != cellDim) || (jacobian.extent(1) != cellDim)), std::invalid_argument,
338 ">>> ERROR (Intrepid::OrientationTools::getJacobianOfOrientationMap): " \
339 "Size of jacobian is not compatible with cell topology.");
340 }
341#endif
342
343 return getJacobianOfOrientationMap(jacobian, cellTopo.getKey(), cellOrt);
344 }
345
346 template<typename outPointViewType>
347 KOKKOS_INLINE_FUNCTION
348 void
350 getJacobianOfOrientationMap(outPointViewType jacobian,
351 const unsigned cellTopoKey,
352 const ordinal_type cellOrt) {
353 switch (cellTopoKey) {
354 case shards::Line<>::key :
355 getLineJacobian(jacobian, cellOrt);
356 break;
357 case shards::Triangle<>::key :
358 getTriangleJacobian(jacobian, cellOrt);
359 break;
360 case shards::Quadrilateral<>::key :
361 getQuadrilateralJacobian(jacobian, cellOrt);
362 break;
363 default: {
364 INTREPID2_TEST_FOR_ABORT( true,
365 ">>> ERROR (Intrepid2::OrientationTools::mapToModifiedReference): " \
366 "Invalid cell topology key." );
367 break;
368 }
369 }
370 }
371
372
373 template<typename TanViewType, typename ParamViewType>
374 KOKKOS_INLINE_FUNCTION
375 void
378 const unsigned subcellTopoKey,
379 const ordinal_type subcellOrd,
380 const ordinal_type ort){
381 typename ParamViewType::non_const_value_type data[4];
382 typename ParamViewType::non_const_type jac(data, 2, 2);
383
384 ordinal_type cellDim = subcellParametrization.extent(1);
385 ordinal_type numTans = subcellParametrization.extent(2)-1;
386
388 for(ordinal_type d=0; d<cellDim; ++d)
389 for(ordinal_type j=0; j<numTans; ++j) {
390 tangents(j,d) = 0;
391 for(ordinal_type k=0; k<numTans; ++k)
393 }
394 }
395
396
397 template<typename TanNormViewType, typename ParamViewType>
399 void
402 const unsigned subcellTopoKey,
403 const ordinal_type subcellOrd,
404 const ordinal_type ort){
405 ordinal_type cellDim = subcellParametrization.extent(1);
406 auto range = Kokkos::pair<ordinal_type, ordinal_type>(0,cellDim-1);
407 auto tangents = Kokkos::subview(tangentsAndNormal,range,Kokkos::ALL());
409
410 //compute normal
411 if(cellDim==2){
412 tangentsAndNormal(1,0) = tangents(0,1);
413 tangentsAndNormal(1,1) = -tangents(0,0);
414 } else { // cellDim==3
415 //cross-product
416 tangentsAndNormal(2,0) = tangents(0,1)*tangents(1,2) - tangents(0,2)*tangents(1,1);
417 tangentsAndNormal(2,1) = tangents(0,2)*tangents(1,0) - tangents(0,0)*tangents(1,2);
418 tangentsAndNormal(2,2) = tangents(0,0)*tangents(1,1) - tangents(0,1)*tangents(1,0);
419 }
420 }
421
422 template<typename coordsViewType, typename subcellCoordsViewType, typename ParamViewType>
424 void
428 const unsigned subcellTopoKey,
429 const ordinal_type subcellOrd,
430 const ordinal_type ort){
431
432 ordinal_type cellDim = subcellParametrization.extent(1);
433 ordinal_type numCoords = subcellCoords.extent(0);
434 ordinal_type subcellDim = subcellCoords.extent(1);
435 auto range = Kokkos::pair<ordinal_type, ordinal_type>(0,subcellDim);
436 auto modSubcellCoords = Kokkos::subview(cellCoords, Kokkos::ALL(),range);
438 typename coordsViewType::value_type subCoord[2];
439
440 for(ordinal_type i=0; i<numCoords; ++i) {
441 for(ordinal_type d=0; d<subcellDim; ++d)
443
444 for(ordinal_type d=0; d<cellDim; ++d) {
446 for(ordinal_type k=0; k<subcellDim; ++k)
448 }
449 }
450 }
451 }
452}
453
454#endif
static KOKKOS_INLINE_FUNCTION void getModifiedTrianglePoint(ValueType &ot0, ValueType &ot1, const ValueType pt0, const ValueType pt1, const ordinal_type ort)
Computes modified point for triangle.
static KOKKOS_INLINE_FUNCTION void getModifiedQuadrilateralPoint(ValueType &ot0, ValueType &ot1, const ValueType pt0, const ValueType pt1, const ordinal_type ort)
Computes modified point for quadrilateral.
static void getJacobianOfOrientationMap(JacobianViewType jacobian, const shards::CellTopology cellTopo, const ordinal_type cellOrt)
Computes jacobian of the parameterization maps of 1- and 2-subcells with orientation.
static KOKKOS_INLINE_FUNCTION void getRefSubcellTangents(TanViewType tangents, const ParamViewType subCellParametrization, const unsigned subcellTopoKey, const ordinal_type subCellOrd, const ordinal_type ort)
Computes the (oriented) subCell tangents.
static KOKKOS_INLINE_FUNCTION void getQuadrilateralJacobian(JacobianViewType jacobian, const ordinal_type ort)
Computes Jacobian of orientation map for quadrilateral.
static KOKKOS_INLINE_FUNCTION void getLineJacobian(JacobianViewType jacobian, const ordinal_type ort)
Computes Jacobian of orientation map for line segment.
static void mapToModifiedReference(outPointViewType outPoints, const refPointViewType refPoints, const shards::CellTopology cellTopo, const ordinal_type cellOrt=0)
Computes modified parameterization maps of 1- and 2-subcells with orientation.
static KOKKOS_INLINE_FUNCTION void getRefSideTangentsAndNormal(TanNormViewType tangentsAndNormal, const ParamViewType subCellParametrization, const unsigned subcellTopoKey, const ordinal_type subCellOrd, const ordinal_type ort)
Computes the (oriented) tangents and normal of the side subCell The normals are only defined for side...
static KOKKOS_INLINE_FUNCTION void mapSubcellCoordsToRefCell(coordsViewType cellCoords, const subcellCoordsViewType subCellCoords, const ParamViewType subcellParametrization, const unsigned subcellTopoKey, const ordinal_type subCellOrd, const ordinal_type ort)
Maps points defined on the subCell manifold into the parent Cell.
static KOKKOS_INLINE_FUNCTION void getModifiedLinePoint(ValueType &ot, const ValueType pt, const ordinal_type ort)
Computes modified point for line segment.
static KOKKOS_INLINE_FUNCTION void getTriangleJacobian(JacobianViewType jacobian, const ordinal_type ort)
Computes Jacobian of orientation map for triangle.