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