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
56namespace 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
Generated by doxygen 1.9.4