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 
23 namespace 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::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::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::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::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::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::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
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::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::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
Generated by   doxygen 1.8.5