Intrepid2
Intrepid2_PyramidCoords.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 
15 #ifndef Intrepid2_Intrepid2_PyramidCoords_h
16 #define Intrepid2_Intrepid2_PyramidCoords_h
17 
18 #include <Kokkos_DynRankView.hpp>
19 
20 #include <Intrepid2_config.h>
21 
22 namespace Intrepid2
23 {
25  template<class PointScalar>
26  KOKKOS_INLINE_FUNCTION
27  void affinePyramid(Kokkos::Array<PointScalar,5> &lambda,
28  Kokkos::Array<Kokkos::Array<PointScalar,3>,5> &lambdaGrad,
29  Kokkos::Array<Kokkos::Array<PointScalar,3>,2> &mu,
30  Kokkos::Array<Kokkos::Array<Kokkos::Array<PointScalar,3>,3>,2> &muGrad,
31  Kokkos::Array<Kokkos::Array<PointScalar,2>,3> &nu,
32  Kokkos::Array<Kokkos::Array<Kokkos::Array<PointScalar,3>,2>,3> &nuGrad,
33  Kokkos::Array<PointScalar,3> &coords)
34  {
35  const auto & x = coords[0];
36  const auto & y = coords[1];
37  const auto & z = coords[2];
38  nu[0][0] = 1. - x - z; // nu_0^{\zeta,\xi_1}
39  nu[0][1] = 1. - y - z; // nu_0^{\zeta,\xi_2}
40  nu[1][0] = x ; // nu_1^{\zeta,\xi_1}
41  nu[1][1] = y ; // nu_1^{\zeta,\xi_2}
42  nu[2][0] = z ; // nu_2^{\zeta,\xi_1}
43  nu[2][1] = z ; // nu_2^{\zeta,\xi_2}
44 
45  nuGrad[0][0][0] = -1. ; // nu_0^{\zeta,\xi_1}_dxi_1
46  nuGrad[0][0][1] = 0. ; // nu_0^{\zeta,\xi_1}_dxi_2
47  nuGrad[0][0][2] = -1. ; // nu_0^{\zeta,\xi_1}_dzeta
48 
49  nuGrad[0][1][0] = 0. ; // nu_0^{\zeta,\xi_2}_dxi_1
50  nuGrad[0][1][1] = -1. ; // nu_0^{\zeta,\xi_2}_dxi_2
51  nuGrad[0][1][2] = -1. ; // nu_0^{\zeta,\xi_2}_dzeta
52 
53  nuGrad[1][0][0] = 1. ; // nu_1^{\zeta,\xi_1}_dxi_1
54  nuGrad[1][0][1] = 0. ; // nu_1^{\zeta,\xi_1}_dxi_2
55  nuGrad[1][0][2] = 0. ; // nu_1^{\zeta,\xi_1}_dzeta
56 
57  nuGrad[1][1][0] = 0. ; // nu_1^{\zeta,\xi_2}_dxi_1
58  nuGrad[1][1][1] = 1. ; // nu_1^{\zeta,\xi_2}_dxi_2
59  nuGrad[1][1][2] = 0. ; // nu_1^{\zeta,\xi_2}_dzeta
60 
61  nuGrad[2][0][0] = 0. ; // nu_2^{\zeta,\xi_1}_dxi_1
62  nuGrad[2][0][1] = 0. ; // nu_2^{\zeta,\xi_1}_dxi_2
63  nuGrad[2][0][2] = 1. ; // nu_2^{\zeta,\xi_1}_dzeta
64 
65  nuGrad[2][1][0] = 0. ; // nu_2^{\zeta,\xi_2}_dxi_1
66  nuGrad[2][1][1] = 0. ; // nu_2^{\zeta,\xi_2}_dxi_2
67  nuGrad[2][1][2] = 1. ; // nu_2^{\zeta,\xi_2}_dzeta
68 
69  // (1 - z) goes in denominator -- so we check for 1-z=0
70  auto & muZ_0 = mu[0][2];
71  auto & muZ_1 = mu[1][2];
72  const double epsilon = 1e-12;
73  muZ_0 = (fabs(1.-z) > epsilon) ? 1. - z : epsilon;
74  muZ_1 = (fabs(1.-z) > epsilon) ? z : 1. - epsilon;
75  PointScalar scaling = 1. / muZ_0;
76  mu[0][0] = 1. - x * scaling;
77  mu[0][1] = 1. - y * scaling;
78  mu[1][0] = x * scaling;
79  mu[1][1] = y * scaling;
80 
81  PointScalar scaling2 = scaling * scaling;
82  muGrad[0][0][0] = -scaling ;
83  muGrad[0][0][1] = 0. ;
84  muGrad[0][0][2] = - x * scaling2;
85 
86  muGrad[0][1][0] = 0. ;
87  muGrad[0][1][1] = -scaling ;
88  muGrad[0][1][2] = -y * scaling2;
89 
90  muGrad[0][2][0] = 0. ;
91  muGrad[0][2][1] = 0. ;
92  muGrad[0][2][2] = -1. ;
93 
94  muGrad[1][0][0] = scaling ;
95  muGrad[1][0][1] = 0. ;
96  muGrad[1][0][2] = x * scaling2;
97 
98  muGrad[1][1][0] = 0. ;
99  muGrad[1][1][1] = scaling ;
100  muGrad[1][1][2] = y * scaling2;
101 
102  muGrad[1][2][0] = 0. ;
103  muGrad[1][2][1] = 0. ;
104  muGrad[1][2][2] = 1. ;
105 
106  lambda[0] = nu[0][0] * mu[0][1];
107  lambda[1] = nu[0][1] * mu[1][0];
108  lambda[2] = nu[1][0] * mu[1][1];
109  lambda[3] = nu[1][1] * mu[0][0];
110  lambda[4] = z;
111 
112  for (int d=0; d<3; d++)
113  {
114  lambdaGrad[0][d] = nu[0][0] * muGrad[0][1][d] + nuGrad[0][0][d] * mu[0][1];
115  lambdaGrad[1][d] = nu[0][1] * muGrad[1][0][d] + nuGrad[0][1][d] * mu[1][0];
116  lambdaGrad[2][d] = nu[1][0] * muGrad[1][1][d] + nuGrad[1][0][d] * mu[1][1];
117  lambdaGrad[3][d] = nu[1][1] * muGrad[0][0][d] + nuGrad[1][1][d] * mu[0][0];
118  }
119  lambdaGrad[4][0] = 0;
120  lambdaGrad[4][1] = 0;
121  lambdaGrad[4][2] = 1;
122  }
123 
125  template<class PointScalar>
126  KOKKOS_INLINE_FUNCTION
127  void transformToESEASPyramid( PointScalar &x_eseas, PointScalar &y_eseas, PointScalar &z_eseas,
128  const PointScalar &x_int2, const PointScalar &y_int2, const PointScalar &z_int2)
129  {
130  x_eseas = (x_int2 + 1. - z_int2) / 2.;
131  y_eseas = (y_int2 + 1. - z_int2) / 2.;
132  z_eseas = z_int2;
133  }
134 
136  template<class OutputScalar>
137  KOKKOS_INLINE_FUNCTION
138  void transformFromESEASPyramidGradient( OutputScalar &dx_int2, OutputScalar &dy_int2, OutputScalar &dz_int2,
139  const OutputScalar &dx_eseas, const OutputScalar &dy_eseas, const OutputScalar &dz_eseas)
140  {
141  dx_int2 = dx_eseas / 2.;
142  dy_int2 = dy_eseas / 2.;
143  dz_int2 = dz_eseas - dx_int2 - dy_int2;
144  }
145 } // end namespace Intrepid2
146 
147 #endif /* Intrepid2_Intrepid2_PyramidCoords_h */
Generated by   doxygen 1.8.5