Intrepid2
Intrepid2_OrientationDef.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_ORIENTATION_DEF_HPP__
49 #define __INTREPID2_ORIENTATION_DEF_HPP__
50 
51 // disable clang warnings
52 #if defined (__clang__) && !defined (__INTEL_COMPILER)
53 #pragma clang system_header
54 #endif
55 
56 namespace Intrepid2 {
57 
58  // ------------------------------------------------------------------------------------
59  // Orientation
60  //
61  //
62  template<typename elemNodeViewType>
63  inline
64  void
65  Orientation::getElementNodeMap(typename elemNodeViewType::non_const_value_type *subCellVerts,
66  ordinal_type &numVerts,
67  const shards::CellTopology cellTopo,
68  const elemNodeViewType elemNodes,
69  const ordinal_type subCellDim,
70  const ordinal_type subCellOrd) {
71  static_assert(Kokkos::Impl::MemorySpaceAccess
72  <Kokkos::HostSpace,typename elemNodeViewType::device_type::memory_space>::accessible,
73  "host space cannot access elemNodeViewType");
74  switch (subCellDim) {
75  case 0: {
76  numVerts = 1;
77  subCellVerts[0] = elemNodes(subCellOrd);
78  break;
79  }
80  default: {
81  numVerts = cellTopo.getVertexCount(subCellDim, subCellOrd);
82  for (ordinal_type i=0;i<numVerts;++i)
83  subCellVerts[i] = elemNodes(cellTopo.getNodeMap(subCellDim, subCellOrd, i));
84  break;
85  }
86  }
87  }
88 
89  template<typename subCellVertType>
90  inline
91  ordinal_type
92  Orientation::getOrientation(const subCellVertType subCellVerts[],
93  const ordinal_type numVerts) {
94  ordinal_type ort = 0;
95  switch (numVerts) {
96  case 2: {// edge
97 #ifdef HAVE_INTREPID2_DEBUG
98  INTREPID2_TEST_FOR_ABORT( ( subCellVerts[0] == subCellVerts[1] ),
99  ">>> ERROR (Intrepid::Orientation::getOrientation): " \
100  "Invalid subCellVerts, same vertex ids are repeated");
101 #endif
102  ort = (subCellVerts[0] > subCellVerts[1]);
103  break;
104  }
105  case 3: {
106 #ifdef HAVE_INTREPID2_DEBUG
107  INTREPID2_TEST_FOR_ABORT( ( subCellVerts[0] == subCellVerts[1] ||
108  subCellVerts[0] == subCellVerts[2] ||
109  subCellVerts[1] == subCellVerts[2] ),
110  ">>> ERROR (Intrepid::Orientation::getOrientation): " \
111  "Invalid subCellVerts, same vertex ids are repeated");
112 #endif
113  ordinal_type rotation = 0; // find smallest vertex id
114  for (ordinal_type i=1;i<3;++i)
115  rotation = ( subCellVerts[i] < subCellVerts[rotation] ? i : rotation );
116 
117  const ordinal_type axes[][2] = { {1,2}, {2,0}, {0,1} };
118  const ordinal_type flip = (subCellVerts[axes[rotation][0]] > subCellVerts[axes[rotation][1]]);
119 
120  ort = flip*3 + rotation;
121  break;
122  }
123  case 4: {
124 #ifdef HAVE_INTREPID2_DEBUG
125  INTREPID2_TEST_FOR_ABORT( ( subCellVerts[0] == subCellVerts[1] ||
126  subCellVerts[0] == subCellVerts[2] ||
127  subCellVerts[0] == subCellVerts[3] ||
128  subCellVerts[1] == subCellVerts[2] ||
129  subCellVerts[1] == subCellVerts[3] ||
130  subCellVerts[2] == subCellVerts[3] ),
131  ">>> ERROR (Intrepid::Orientation::getGlobalVertexNodes): " \
132  "Invalid subCellVerts, same vertex ids are repeated");
133 #endif
134  ordinal_type rotation = 0; // find smallest vertex id
135  for (ordinal_type i=1;i<4;++i)
136  rotation = ( subCellVerts[i] < subCellVerts[rotation] ? i : rotation );
137 
138  const ordinal_type axes[][2] = { {1,3}, {2,0}, {3,1}, {0,2} };
139  const ordinal_type flip = (subCellVerts[axes[rotation][0]] > subCellVerts[axes[rotation][1]]);
140 
141  ort = flip*4 + rotation;
142  break;
143  }
144  default: {
145  INTREPID2_TEST_FOR_ABORT( true,
146  ">>> ERROR (Intrepid::Orientation::getOrientation): " \
147  "Invalid numVerts (2 (edge),3 (triangle) and 4 (quadrilateral) are allowed)");
148  break;
149  }
150  }
151  return ort;
152  }
153 
154  template<typename elemNodeViewType>
155  inline
156  Orientation
157  Orientation::getOrientation(const shards::CellTopology cellTopo,
158  const elemNodeViewType elemNodes) {
159  static_assert(Kokkos::Impl::MemorySpaceAccess
160  <Kokkos::HostSpace,typename elemNodeViewType::device_type::memory_space>::accessible,
161  "host space cannot access elemNodeViewType");
162 
163  Orientation ort;
164  const ordinal_type nedge = cellTopo.getEdgeCount();
165 
166  if (nedge > 0) {
167  typename elemNodeViewType::non_const_value_type vertsSubCell[2];
168  ordinal_type orts[12], nvertSubCell;
169  for (ordinal_type i=0;i<nedge;++i) {
170  Orientation::getElementNodeMap(vertsSubCell,
171  nvertSubCell,
172  cellTopo,
173  elemNodes,
174  1, i);
175  orts[i] = Orientation::getOrientation(vertsSubCell, nvertSubCell);
176  }
177  ort.setEdgeOrientation(nedge, orts);
178  }
179  const ordinal_type nface = cellTopo.getFaceCount();
180  if (nface > 0) {
181  typename elemNodeViewType::non_const_value_type vertsSubCell[4];
182  ordinal_type orts[6], nvertSubCell;
183  for (ordinal_type i=0;i<nface;++i) {
184  Orientation::getElementNodeMap(vertsSubCell,
185  nvertSubCell,
186  cellTopo,
187  elemNodes,
188  2, i);
189  orts[i] = Orientation::getOrientation(vertsSubCell, nvertSubCell);
190  }
191  ort.setFaceOrientation(nface, orts);
192  }
193  return ort;
194  }
195 
196  inline
197  ordinal_type
198  Orientation::getEdgeOrdinalOfFace(const ordinal_type subsubcellOrd,
199  const ordinal_type subcellOrd,
200  const shards::CellTopology cellTopo) {
201  ordinal_type r_val = -1;
202 
203  const auto cellBaseKey = cellTopo.getBaseKey();
204  if (cellBaseKey == shards::Hexahedron<>::key) {
205  INTREPID2_TEST_FOR_EXCEPTION( !(subcellOrd < 6) &&
206  !(subsubcellOrd < 4),
207  std::logic_error,
208  "subcell and subsubcell information are not correct" );
209  const int quad_to_hex_edges[6][4] = { { 0, 9, 4, 8 },
210  { 1,10, 5, 9 },
211  { 2,11, 6,10 },
212  { 8, 7,11, 3 },
213  { 3, 2, 1, 0 },
214  { 4, 5, 6, 7 } };
215  r_val = quad_to_hex_edges[subcellOrd][subsubcellOrd];
216  } else if (cellBaseKey == shards::Tetrahedron<>::key) {
217  INTREPID2_TEST_FOR_EXCEPTION( !(subcellOrd < 4) &&
218  !(subsubcellOrd < 3),
219  std::logic_error,
220  "subcell and subsubcell information are not correct" );
221  const ordinal_type tri_to_tet_edges[4][3] = { { 0, 4, 3 },
222  { 1, 5, 4 },
223  { 3, 5, 2 },
224  { 2, 1, 0 } };
225  r_val = tri_to_tet_edges[subcellOrd][subsubcellOrd];
226  } else {
227  INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
228  "cellTopo is not supported: try TET and HEX" );
229  }
230  return r_val;
231  }
232 
233  /*
234  template<typename refTanType>
235  inline
236  void
237  Orientation::getReferenceEdgeTangent(const refTanType &tanE,
238  const ordinal_type subcellOrd,
239  const shards::CellTopology cellTopo,
240  const ordinal_type ort,
241  const bool is_normalize) {
242  const auto cellBaseKey = cellTopo.getBaseKey();
243  INTREPID2_TEST_FOR_EXCEPTION( !(cellBaseKey == shards::Hexahedron<>::key && subcellOrd < 12) &&
244  !(cellBaseKey == shards::Tetrahedron<>::key && subcellOrd < 6) &&
245  !(cellBaseKey == shards::Quadrilateral<>::key && subcellOrd < 4) &&
246  !(cellBaseKey == shards::Triangle<>::key && subcellOrd < 3),
247  std::logic_error,
248  "subcell information are not correct" );
249  const ordinal_type i[2][2] = { { 0, 1 },
250  { 1, 0 } };
251  const unsigned int v[2] = { cellTopo.getNodeMap(1, subcellOrd, 0),
252  cellTopo.getNodeMap(1, subcellOrd, 1) };
253 
254  auto normalize = [&](double *vv, ordinal_type iend) {
255  double norm = 0.0;
256  for (ordinal_type ii=0;ii<iend;++ii)
257  norm += vv[ii]*vv[ii];
258  norm = std::sqrt(norm);
259  for (ordinal_type ii=0;ii<iend;++ii)
260  vv[ii] /= norm;
261  };
262 
263  auto assign_tangent = [&](refTanType t, double *vv, ordinal_type iend) {
264  for (ordinal_type ii=0;ii<iend;++ii)
265  t(ii) = vv[ii];
266  };
267 
268  double t[3] = {};
269  const int cell_dim = cellTopo.getDimension();
270  if (cellBaseKey == shards::Hexahedron<>::key) {
271  const double hex_verts[8][3] = { { -1.0, -1.0, -1.0 },
272  { 1.0, -1.0, -1.0 },
273  { 1.0, 1.0, -1.0 },
274  { -1.0, 1.0, -1.0 },
275  //
276  { -1.0, -1.0, 1.0 },
277  { 1.0, -1.0, 1.0 },
278  { 1.0, 1.0, 1.0 },
279  { -1.0, 1.0, 1.0 } };
280  for (ordinal_type k=0;k<3;++k) {
281  const ordinal_type *ii = &i[ort][0];
282  t[k] = hex_verts[v[ii[1]]][k] - hex_verts[v[ii[0]]][k];
283  }
284  } else if (cellBaseKey == shards::Tetrahedron<>::key) {
285  const double tet_verts[4][3] = { { 0.0, 0.0, 0.0 },
286  { 1.0, 0.0, 0.0 },
287  { 0.0, 1.0, 0.0 },
288  { 0.0, 0.0, 1.0 } };
289  for (ordinal_type k=0;k<3;++k) {
290  const ordinal_type *ii = &i[ort][0];
291  t[k] = tet_verts[v[ii[1]]][k] - tet_verts[v[ii[0]]][k];
292  }
293  } else if (cellBaseKey == shards::Quadrilateral<>::key) {
294  const double quad_verts[8][3] = { { -1.0, -1.0 },
295  { 1.0, -1.0 },
296  { 1.0, 1.0 },
297  { -1.0, 1.0 } };
298  for (ordinal_type k=0;k<2;++k) {
299  const ordinal_type *ii = &i[ort][0];
300  t[k] = quad_verts[v[ii[1]]][k] - quad_verts[v[ii[0]]][k];
301  }
302  } else if (cellBaseKey == shards::Triangle<>::key) {
303  const double tri_verts[4][3] = { { 0.0, 0.0 },
304  { 1.0, 0.0 },
305  { 0.0, 1.0 } };
306  for (ordinal_type k=0;k<2;++k) {
307  const ordinal_type *ii = &i[ort][0];
308  t[k] = tri_verts[v[ii[1]]][k] - tri_verts[v[ii[0]]][k];
309  }
310  } else {
311  INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
312  "cellTopo is not supported: try TET and HEX" );
313  }
314 
315  if (is_normalize) normalize(t, cell_dim);
316  assign_tangent(tanE, t, cell_dim);
317  }
318 
319  template<typename refTanType>
320  inline
321  void
322  Orientation::getReferenceFaceTangents(const refTanType &tanU,
323  const refTanType &tanV,
324  const ordinal_type subcellOrd,
325  const shards::CellTopology cellTopo,
326  const ordinal_type ort,
327  const bool is_normalize) {
328  const auto cellBaseKey = cellTopo.getBaseKey();
329 
330  auto normalize = [&](double *v, ordinal_type iend) {
331  double norm = 0.0;
332  for (ordinal_type i=0;i<iend;++i)
333  norm += v[i]*v[i];
334  norm = std::sqrt(norm);
335  for (ordinal_type i=0;i<iend;++i)
336  v[i] /= norm;
337  };
338 
339  auto assign_tangent = [&](refTanType t, double *v, ordinal_type iend) {
340  for (ordinal_type i=0;i<iend;++i)
341  t(i) = v[i];
342  };
343 
344  double tu[3], tv[3];
345  if (cellBaseKey == shards::Hexahedron<>::key) {
346  INTREPID2_TEST_FOR_EXCEPTION( !(subcellOrd < 6),
347  std::logic_error,
348  "subcell information are not correct" );
349  const double hex_verts[8][3] = { { -1.0, -1.0, -1.0 },
350  { 1.0, -1.0, -1.0 },
351  { 1.0, 1.0, -1.0 },
352  { -1.0, 1.0, -1.0 },
353  //
354  { -1.0, -1.0, 1.0 },
355  { 1.0, -1.0, 1.0 },
356  { 1.0, 1.0, 1.0 },
357  { -1.0, 1.0, 1.0 } };
358  const unsigned int v[4] = { cellTopo.getNodeMap(2, subcellOrd, 0),
359  cellTopo.getNodeMap(2, subcellOrd, 1),
360  cellTopo.getNodeMap(2, subcellOrd, 2),
361  cellTopo.getNodeMap(2, subcellOrd, 3) };
362  const ordinal_type i[8][4] = { { 0, 1, 2, 3 },
363  { 1, 2, 3, 0 },
364  { 2, 3, 0, 1 },
365  { 3, 0, 1, 2 },
366  //
367  { 0, 3, 2, 1 },
368  { 1, 0, 3, 2 },
369  { 2, 1, 0, 3 },
370  { 3, 2, 1, 0 } };
371  for (ordinal_type k=0;k<3;++k) {
372  const ordinal_type *ii = &i[ort][0];
373 
374  tu[k] = hex_verts[v[ii[1]]][k] - hex_verts[v[ii[0]]][k];
375  tv[k] = hex_verts[v[ii[3]]][k] - hex_verts[v[ii[0]]][k];
376  }
377 
378  } else if (cellBaseKey == shards::Tetrahedron<>::key) {
379  INTREPID2_TEST_FOR_EXCEPTION( !(subcellOrd < 4),
380  std::logic_error,
381  "subcell information are not correct" );
382  const double tet_verts[4][3] = { { 0.0, 0.0, 0.0 },
383  { 1.0, 0.0, 0.0 },
384  { 0.0, 1.0, 0.0 },
385  { 0.0, 0.0, 1.0 } };
386  const unsigned int v[4] = { cellTopo.getNodeMap(2, subcellOrd, 0),
387  cellTopo.getNodeMap(2, subcellOrd, 1),
388  cellTopo.getNodeMap(2, subcellOrd, 2) };
389  const ordinal_type i[6][3] = { { 0, 1, 2 },
390  { 1, 2, 0 },
391  { 2, 0, 1 },
392  //
393  { 0, 2, 1 },
394  { 1, 0, 2 },
395  { 2, 1, 0 } };
396  for (ordinal_type k=0;k<3;++k) {
397  const ordinal_type *ii = &i[ort][0];
398 
399  tu[k] = tet_verts[v[ii[1]]][k] - tet_verts[v[ii[0]]][k];
400  tv[k] = tet_verts[v[ii[2]]][k] - tet_verts[v[ii[0]]][k];
401  }
402 
403  } else {
404  INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
405  "cellTopo is not supported: try TET and HEX" );
406  }
407 
408  if (is_normalize) {
409  normalize(tu, 3);
410  normalize(tv, 3);
411  }
412 
413  assign_tangent(tanU, tu, 3);
414  assign_tangent(tanV, tv, 3);
415  }
416 
417  template<typename refNormalType>
418  inline
419  void
420  Orientation::getReferenceFaceNormal(const refNormalType &normalV,
421  const ordinal_type subcellOrd,
422  const shards::CellTopology cellTopo,
423  const ordinal_type ort,
424  const bool is_normalize) {
425  const auto cellBaseKey = cellTopo.getBaseKey();
426 
427  auto normalize = [&](double *v, ordinal_type iend) {
428  double norm = 0.0;
429  for (ordinal_type i=0;i<iend;++i)
430  norm += v[i]*v[i];
431  norm = std::sqrt(norm);
432  for (ordinal_type i=0;i<iend;++i)
433  v[i] /= norm;
434  };
435 
436  auto assign_normal = [&](refNormalType n, double *v, ordinal_type iend) {
437  for (ordinal_type i=0;i<iend;++i)
438  n(i) = v[i];
439  };
440 
441  double buf[2][3];
442  Kokkos::View<double*,Kokkos::HostSpace> tanU(&buf[0][0], 3);
443  Kokkos::View<double*,Kokkos::HostSpace> tanV(&buf[1][0], 3);
444 
445  getReferenceFaceTangents(tanU, tanV,
446  subcellOrd,
447  cellTopo,
448  ort,
449  false);
450 
451  // cross product
452  double v[3];
453  v[0] = tanU(1)*tanV(2) - tanU(2)*tanV(1);
454  v[1] = tanU(2)*tanV(0) - tanU(0)*tanV(2);
455  v[2] = tanU(0)*tanV(1) - tanU(1)*tanV(0);
456 
457  if (is_normalize) normalize(v, 3);
458  assign_normal(normalV, v, 3);
459  }
460  */
461 
462  KOKKOS_INLINE_FUNCTION
463  Orientation::Orientation()
464  : _edgeOrt(0), _faceOrt(0) {}
465 
466  KOKKOS_INLINE_FUNCTION
467  bool
468  Orientation::isAlignedToReference() const {
469  return (_edgeOrt == 0 && _faceOrt == 0);
470  }
471 
472  KOKKOS_INLINE_FUNCTION
473  void
474  Orientation::setEdgeOrientation(const ordinal_type numEdge, const ordinal_type edgeOrt[]) {
475 #ifdef HAVE_INTREPID2_DEBUG
476  INTREPID2_TEST_FOR_ABORT( !( 3 <= numEdge && numEdge <= 12 ),
477  ">>> ERROR (Intrepid::Orientation::setEdgeOrientation): " \
478  "Invalid numEdge (3--12)");
479 #endif
480  _edgeOrt = 0;
481  for (ordinal_type i=0;i<numEdge;++i)
482  _edgeOrt |= (edgeOrt[i] & 1) << i;
483  }
484 
485  KOKKOS_INLINE_FUNCTION
486  void
487  Orientation::getEdgeOrientation(ordinal_type *edgeOrt, const ordinal_type numEdge) const {
488 #ifdef HAVE_INTREPID2_DEBUG
489  INTREPID2_TEST_FOR_ABORT( !( 3 <= numEdge && numEdge <= 12 ),
490  ">>> ERROR (Intrepid::Orientation::setEdgeOrientation): " \
491  "Invalid numEdge (3--12)");
492 #endif
493  for (ordinal_type i=0;i<numEdge;++i)
494  edgeOrt[i] = (_edgeOrt & (1 << i)) >> i;
495  }
496 
497  KOKKOS_INLINE_FUNCTION
498  void
499  Orientation::setFaceOrientation(const ordinal_type numFace, const ordinal_type faceOrt[]) {
500 #ifdef HAVE_INTREPID2_DEBUG
501  INTREPID2_TEST_FOR_ABORT( !( 4 <= numFace && numFace <= 6 ),
502  ">>> ERROR (Intrepid::Orientation::setFaceOrientation): "
503  "Invalid numFace (4--6)");
504 #endif
505  _faceOrt = 0;
506  for (ordinal_type i=0;i<numFace;++i) {
507  const ordinal_type s = i*3;
508  _faceOrt |= (faceOrt[i] & 7) << s;
509  }
510  }
511 
512  KOKKOS_INLINE_FUNCTION
513  void
514  Orientation::getFaceOrientation(ordinal_type *faceOrt, const ordinal_type numFace) const {
515 #ifdef HAVE_INTREPID2_DEBUG
516  INTREPID2_TEST_FOR_ABORT( !( 4 <= numFace && numFace <= 6 ),
517  ">>> ERROR (Intrepid::Orientation::setEdgeOrientation): "
518  "Invalid numFace (4--6)");
519 #endif
520  for (ordinal_type i=0;i<numFace;++i) {
521  const ordinal_type s = i*3;
522  faceOrt[i] = (_faceOrt & (7 << s)) >> s;
523  }
524  }
525 
526  inline std::string Orientation::to_string() const {
527  return "Orientation{ face: " + std::to_string(_faceOrt) + "; edge: " + std::to_string(_edgeOrt) + " }";
528  }
529 }
530 
531 #endif