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