Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_ConvertNormalToRotationMatrix.hpp
Go to the documentation of this file.
1#ifndef PANZER_CONVERT_NORMAL_TO_ROTATION_MATRIX_HPP
2#define PANZER_CONVERT_NORMAL_TO_ROTATION_MATRIX_HPP
3
4namespace panzer {
5
6template <typename Scalar>
7KOKKOS_INLINE_FUNCTION
8void
9convertNormalToRotationMatrix(const Scalar normal[3], Scalar transverse[3], Scalar binormal[3])
10{
11 using T = Scalar;
12
13 const T n = sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
14
15 // If this fails then the geometry for this cell is probably undefined
16 if(n > 0.){
17 // Make sure transverse is not parallel to normal within some margin of error
18 transverse[0]=0.;transverse[1]=1.;transverse[2]=0.;
19 if(Kokkos::fabs(normal[0]*transverse[0]+normal[1]*transverse[1])>0.9){
20 transverse[0]=1.;transverse[1]=0.;
21 }
22
23 const T nt = normal[0]*transverse[0]+normal[1]*transverse[1]+normal[2]*transverse[2];
24
25 // Note normal has unit length
26 const T mult = nt/(n*n); // = nt
27
28 // Remove normal projection from transverse
29 for(int dim=0;dim<3;++dim){
30 transverse[dim] = transverse[dim] - mult * normal[dim];
31 }
32
33 const T t = sqrt(transverse[0]*transverse[0]+transverse[1]*transverse[1]+transverse[2]*transverse[2]);
34 KOKKOS_ASSERT(t != 0.);
35 for(int dim=0;dim<3;++dim){
36 transverse[dim] /= t;
37 }
38
39 // We assume a right handed system such that b = n \times t
40 binormal[0] = (normal[1] * transverse[2] - normal[2] * transverse[1]);
41 binormal[1] = (normal[2] * transverse[0] - normal[0] * transverse[2]);
42 binormal[2] = (normal[0] * transverse[1] - normal[1] * transverse[0]);
43
44 // Normalize binormal
45 const T b = sqrt(binormal[0]*binormal[0]+binormal[1]*binormal[1]+binormal[2]*binormal[2]);
46 for(int dim=0;dim<3;++dim){
47 binormal[dim] /= b;
48 }
49 } else {
50 transverse[0] = 0.;
51 transverse[1] = 0.;
52 transverse[2] = 0.;
53 binormal[0] = 0.;
54 binormal[1] = 0.;
55 binormal[2] = 0.;
56 }
57}
58
59}
60
61#endif
KOKKOS_INLINE_FUNCTION void convertNormalToRotationMatrix(const Scalar normal[3], Scalar transverse[3], Scalar binormal[3])