Compadre  1.5.7
Compadre_Manifold_Functions.hpp
Go to the documentation of this file.
1 
2 namespace Compadre {
3 
4  //! Metric factor (det(G)) at any point in the local chart
5  KOKKOS_INLINE_FUNCTION
6  double MetricFactor(const scratch_vector_type a_, const double h, const double u1, const double u2) {
7  double det_g = 1.0;
8  double q1 = 0;
9  double q2 = 0;
10  switch (a_.extent(0)) {
11  case 45:
12  q1 += (1.0/5040.0)*a_[36]*std::pow(u1, 7)/std::pow(h, 8) + (1.0/720.0)*a_[37]*std::pow(u1, 6)*u2/std::pow(h, 8) + (1.0/240.0)*a_[38]*std::pow(u1, 5)*std::pow(u2, 2)/std::pow(h, 8) + (1.0/144.0)*a_[39]*std::pow(u1, 4)*std::pow(u2, 3)/std::pow(h, 8) + (1.0/144.0)*a_[40]*std::pow(u1, 3)*std::pow(u2, 4)/std::pow(h, 8) + (1.0/240.0)*a_[41]*std::pow(u1, 2)*std::pow(u2, 5)/std::pow(h, 8) + (1.0/720.0)*a_[42]*u1*std::pow(u2, 6)/std::pow(h, 8) + (1.0/5040.0)*a_[43]*std::pow(u2, 7)/std::pow(h, 8);
13  q2 += (1.0/5040.0)*a_[37]*std::pow(u1, 7)/std::pow(h, 8) + (1.0/720.0)*a_[38]*std::pow(u1, 6)*u2/std::pow(h, 8) + (1.0/240.0)*a_[39]*std::pow(u1, 5)*std::pow(u2, 2)/std::pow(h, 8) + (1.0/144.0)*a_[40]*std::pow(u1, 4)*std::pow(u2, 3)/std::pow(h, 8) + (1.0/144.0)*a_[41]*std::pow(u1, 3)*std::pow(u2, 4)/std::pow(h, 8) + (1.0/240.0)*a_[42]*std::pow(u1, 2)*std::pow(u2, 5)/std::pow(h, 8) + (1.0/720.0)*a_[43]*u1*std::pow(u2, 6)/std::pow(h, 8) + (1.0/5040.0)*a_[44]*std::pow(u2, 7)/std::pow(h, 8);
14  case 36:
15  q1 += (1.0/720.0)*a_[28]*std::pow(u1, 6)/std::pow(h, 7) + (1.0/120.0)*a_[29]*std::pow(u1, 5)*u2/std::pow(h, 7) + (1.0/48.0)*a_[30]*std::pow(u1, 4)*std::pow(u2, 2)/std::pow(h, 7) + (1.0/36.0)*a_[31]*std::pow(u1, 3)*std::pow(u2, 3)/std::pow(h, 7) + (1.0/48.0)*a_[32]*std::pow(u1, 2)*std::pow(u2, 4)/std::pow(h, 7) + (1.0/120.0)*a_[33]*u1*std::pow(u2, 5)/std::pow(h, 7) + (1.0/720.0)*a_[34]*std::pow(u2, 6)/std::pow(h, 7);
16  q2 += (1.0/720.0)*a_[29]*std::pow(u1, 6)/std::pow(h, 7) + (1.0/120.0)*a_[30]*std::pow(u1, 5)*u2/std::pow(h, 7) + (1.0/48.0)*a_[31]*std::pow(u1, 4)*std::pow(u2, 2)/std::pow(h, 7) + (1.0/36.0)*a_[32]*std::pow(u1, 3)*std::pow(u2, 3)/std::pow(h, 7) + (1.0/48.0)*a_[33]*std::pow(u1, 2)*std::pow(u2, 4)/std::pow(h, 7) + (1.0/120.0)*a_[34]*u1*std::pow(u2, 5)/std::pow(h, 7) + (1.0/720.0)*a_[35]*std::pow(u2, 6)/std::pow(h, 7);
17  case 28:
18  q1 += (1.0/120.0)*a_[21]*std::pow(u1, 5)/std::pow(h, 6) + (1.0/24.0)*a_[22]*std::pow(u1, 4)*u2/std::pow(h, 6) + (1.0/12.0)*a_[23]*std::pow(u1, 3)*std::pow(u2, 2)/std::pow(h, 6) + (1.0/12.0)*a_[24]*std::pow(u1, 2)*std::pow(u2, 3)/std::pow(h, 6) + (1.0/24.0)*a_[25]*u1*std::pow(u2, 4)/std::pow(h, 6) + (1.0/120.0)*a_[26]*std::pow(u2, 5)/std::pow(h, 6);
19  q2 += (1.0/120.0)*a_[22]*std::pow(u1, 5)/std::pow(h, 6) + (1.0/24.0)*a_[23]*std::pow(u1, 4)*u2/std::pow(h, 6) + (1.0/12.0)*a_[24]*std::pow(u1, 3)*std::pow(u2, 2)/std::pow(h, 6) + (1.0/12.0)*a_[25]*std::pow(u1, 2)*std::pow(u2, 3)/std::pow(h, 6) + (1.0/24.0)*a_[26]*u1*std::pow(u2, 4)/std::pow(h, 6) + (1.0/120.0)*a_[27]*std::pow(u2, 5)/std::pow(h, 6);
20  case 21:
21  q1 += (1.0/24.0)*a_[15]*std::pow(u1, 4)/std::pow(h, 5) + (1.0/6.0)*a_[16]*std::pow(u1, 3)*u2/std::pow(h, 5) + (1.0/4.0)*a_[17]*std::pow(u1, 2)*std::pow(u2, 2)/std::pow(h, 5) + (1.0/6.0)*a_[18]*u1*std::pow(u2, 3)/std::pow(h, 5) + (1.0/24.0)*a_[19]*std::pow(u2, 4)/std::pow(h, 5);
22  q2 += (1.0/24.0)*a_[16]*std::pow(u1, 4)/std::pow(h, 5) + (1.0/6.0)*a_[17]*std::pow(u1, 3)*u2/std::pow(h, 5) + (1.0/4.0)*a_[18]*std::pow(u1, 2)*std::pow(u2, 2)/std::pow(h, 5) + (1.0/6.0)*a_[19]*u1*std::pow(u2, 3)/std::pow(h, 5) + (1.0/24.0)*a_[20]*std::pow(u2, 4)/std::pow(h, 5);
23  case 15:
24  q1 += (1.0/6.0)*a_[10]*std::pow(u1, 3)/std::pow(h, 4) + (1.0/2.0)*a_[11]*std::pow(u1, 2)*u2/std::pow(h, 4) + (1.0/2.0)*a_[12]*u1*std::pow(u2, 2)/std::pow(h, 4) + (1.0/6.0)*a_[13]*std::pow(u2, 3)/std::pow(h, 4);
25  q2 += (1.0/6.0)*a_[11]*std::pow(u1, 3)/std::pow(h, 4) + (1.0/2.0)*a_[12]*std::pow(u1, 2)*u2/std::pow(h, 4) + (1.0/2.0)*a_[13]*u1*std::pow(u2, 2)/std::pow(h, 4) + (1.0/6.0)*a_[14]*std::pow(u2, 3)/std::pow(h, 4);
26  case 10:
27  q1 += (1.0/2.0)*a_[6]*std::pow(u1, 2)/std::pow(h, 3) + a_[7]*u1*u2/std::pow(h, 3) + (1.0/2.0)*a_[8]*std::pow(u2, 2)/std::pow(h, 3);
28  q2 += (1.0/2.0)*a_[7]*std::pow(u1, 2)/std::pow(h, 3) + a_[8]*u1*u2/std::pow(h, 3) + (1.0/2.0)*a_[9]*std::pow(u2, 2)/std::pow(h, 3);
29  case 6:
30  q1 += a_[3]*u1/std::pow(h, 2) + a_[4]*u2/std::pow(h, 2);
31  q2 += a_[4]*u1/std::pow(h, 2) + a_[5]*u2/std::pow(h, 2);
32  case 3:
33  q1 += a_[1]/h;
34  q2 += a_[2]/h;
35  case 1:
36  break;
37  default:
38  compadre_kernel_assert_release(false && "curvature polynomial order is greater than 8.");
39 
40  }
41  q1 = std::pow(q1, 2);
42  q2 = std::pow(q2, 2);
43  det_g += q1;
44  det_g += q2;
45  return det_g;
46  }
47 
48  //! Gaussian curvature K at any point in the local chart
49  KOKKOS_INLINE_FUNCTION
50  double GaussianCurvature(const scratch_vector_type a_, const double h, const double u1, const double u2) {
51 
52  double q1=0, q2=0, q3=0;
53  switch (a_.extent(0)) {
54  case 45:
55  q1 += (1.0/720.0)*a_[36]*std::pow(u1, 6)/std::pow(h, 8) + (1.0/120.0)*a_[37]*std::pow(u1, 5)*u2/std::pow(h, 8) + (1.0/48.0)*a_[38]*std::pow(u1, 4)*std::pow(u2, 2)/std::pow(h, 8) + (1.0/36.0)*a_[39]*std::pow(u1, 3)*std::pow(u2, 3)/std::pow(h, 8) + (1.0/48.0)*a_[40]*std::pow(u1, 2)*std::pow(u2, 4)/std::pow(h, 8) + (1.0/120.0)*a_[41]*u1*std::pow(u2, 5)/std::pow(h, 8) + (1.0/720.0)*a_[42]*std::pow(u2, 6)/std::pow(h, 8);
56  q2 += (1.0/720.0)*a_[38]*std::pow(u1, 6)/std::pow(h, 8) + (1.0/120.0)*a_[39]*std::pow(u1, 5)*u2/std::pow(h, 8) + (1.0/48.0)*a_[40]*std::pow(u1, 4)*std::pow(u2, 2)/std::pow(h, 8) + (1.0/36.0)*a_[41]*std::pow(u1, 3)*std::pow(u2, 3)/std::pow(h, 8) + (1.0/48.0)*a_[42]*std::pow(u1, 2)*std::pow(u2, 4)/std::pow(h, 8) + (1.0/120.0)*a_[43]*u1*std::pow(u2, 5)/std::pow(h, 8) + (1.0/720.0)*a_[44]*std::pow(u2, 6)/std::pow(h, 8);
57  q3 += (1.0/720.0)*a_[37]*std::pow(u1, 6)/std::pow(h, 8) + (1.0/120.0)*a_[38]*std::pow(u1, 5)*u2/std::pow(h, 8) + (1.0/48.0)*a_[39]*std::pow(u1, 4)*std::pow(u2, 2)/std::pow(h, 8) + (1.0/36.0)*a_[40]*std::pow(u1, 3)*std::pow(u2, 3)/std::pow(h, 8) + (1.0/48.0)*a_[41]*std::pow(u1, 2)*std::pow(u2, 4)/std::pow(h, 8) + (1.0/120.0)*a_[42]*u1*std::pow(u2, 5)/std::pow(h, 8) + (1.0/720.0)*a_[43]*std::pow(u2, 6)/std::pow(h, 8);
58  case 36:
59  q1 += (1.0/120.0)*a_[28]*std::pow(u1, 5)/std::pow(h, 7) + (1.0/24.0)*a_[29]*std::pow(u1, 4)*u2/std::pow(h, 7) + (1.0/12.0)*a_[30]*std::pow(u1, 3)*std::pow(u2, 2)/std::pow(h, 7) + (1.0/12.0)*a_[31]*std::pow(u1, 2)*std::pow(u2, 3)/std::pow(h, 7) + (1.0/24.0)*a_[32]*u1*std::pow(u2, 4)/std::pow(h, 7) + (1.0/120.0)*a_[33]*std::pow(u2, 5)/std::pow(h, 7);
60  q2 += (1.0/120.0)*a_[30]*std::pow(u1, 5)/std::pow(h, 7) + (1.0/24.0)*a_[31]*std::pow(u1, 4)*u2/std::pow(h, 7) + (1.0/12.0)*a_[32]*std::pow(u1, 3)*std::pow(u2, 2)/std::pow(h, 7) + (1.0/12.0)*a_[33]*std::pow(u1, 2)*std::pow(u2, 3)/std::pow(h, 7) + (1.0/24.0)*a_[34]*u1*std::pow(u2, 4)/std::pow(h, 7) + (1.0/120.0)*a_[35]*std::pow(u2, 5)/std::pow(h, 7);
61  q3 += (1.0/120.0)*a_[29]*std::pow(u1, 5)/std::pow(h, 7) + (1.0/24.0)*a_[30]*std::pow(u1, 4)*u2/std::pow(h, 7) + (1.0/12.0)*a_[31]*std::pow(u1, 3)*std::pow(u2, 2)/std::pow(h, 7) + (1.0/12.0)*a_[32]*std::pow(u1, 2)*std::pow(u2, 3)/std::pow(h, 7) + (1.0/24.0)*a_[33]*u1*std::pow(u2, 4)/std::pow(h, 7) + (1.0/120.0)*a_[34]*std::pow(u2, 5)/std::pow(h, 7);
62  case 28:
63  q1 += (1.0/24.0)*a_[21]*std::pow(u1, 4)/std::pow(h, 6) + (1.0/6.0)*a_[22]*std::pow(u1, 3)*u2/std::pow(h, 6) + (1.0/4.0)*a_[23]*std::pow(u1, 2)*std::pow(u2, 2)/std::pow(h, 6) + (1.0/6.0)*a_[24]*u1*std::pow(u2, 3)/std::pow(h, 6) + (1.0/24.0)*a_[25]*std::pow(u2, 4)/std::pow(h, 6);
64  q2 += (1.0/24.0)*a_[23]*std::pow(u1, 4)/std::pow(h, 6) + (1.0/6.0)*a_[24]*std::pow(u1, 3)*u2/std::pow(h, 6) + (1.0/4.0)*a_[25]*std::pow(u1, 2)*std::pow(u2, 2)/std::pow(h, 6) + (1.0/6.0)*a_[26]*u1*std::pow(u2, 3)/std::pow(h, 6) + (1.0/24.0)*a_[27]*std::pow(u2, 4)/std::pow(h, 6);
65  q3 += (1.0/24.0)*a_[22]*std::pow(u1, 4)/std::pow(h, 6) + (1.0/6.0)*a_[23]*std::pow(u1, 3)*u2/std::pow(h, 6) + (1.0/4.0)*a_[24]*std::pow(u1, 2)*std::pow(u2, 2)/std::pow(h, 6) + (1.0/6.0)*a_[25]*u1*std::pow(u2, 3)/std::pow(h, 6) + (1.0/24.0)*a_[26]*std::pow(u2, 4)/std::pow(h, 6);
66  case 21:
67  q1 += (1.0/6.0)*a_[15]*std::pow(u1, 3)/std::pow(h, 5) + (1.0/2.0)*a_[16]*std::pow(u1, 2)*u2/std::pow(h, 5) + (1.0/2.0)*a_[17]*u1*std::pow(u2, 2)/std::pow(h, 5) + (1.0/6.0)*a_[18]*std::pow(u2, 3)/std::pow(h, 5);
68  q2 += (1.0/6.0)*a_[17]*std::pow(u1, 3)/std::pow(h, 5) + (1.0/2.0)*a_[18]*std::pow(u1, 2)*u2/std::pow(h, 5) + (1.0/2.0)*a_[19]*u1*std::pow(u2, 2)/std::pow(h, 5) + (1.0/6.0)*a_[20]*std::pow(u2, 3)/std::pow(h, 5);
69  q3 += (1.0/6.0)*a_[16]*std::pow(u1, 3)/std::pow(h, 5) + (1.0/2.0)*a_[17]*std::pow(u1, 2)*u2/std::pow(h, 5) + (1.0/2.0)*a_[18]*u1*std::pow(u2, 2)/std::pow(h, 5) + (1.0/6.0)*a_[19]*std::pow(u2, 3)/std::pow(h, 5);
70  case 15:
71  q1 += (1.0/2.0)*a_[10]*std::pow(u1, 2)/std::pow(h, 4) + a_[11]*u1*u2/std::pow(h, 4) + (1.0/2.0)*a_[12]*std::pow(u2, 2)/std::pow(h, 4);
72  q2 += (1.0/2.0)*a_[12]*std::pow(u1, 2)/std::pow(h, 4) + a_[13]*u1*u2/std::pow(h, 4) + (1.0/2.0)*a_[14]*std::pow(u2, 2)/std::pow(h, 4);
73  q3 += (1.0/2.0)*a_[11]*std::pow(u1, 2)/std::pow(h, 4) + a_[12]*u1*u2/std::pow(h, 4) + (1.0/2.0)*a_[13]*std::pow(u2, 2)/std::pow(h, 4);
74  case 10:
75  q1 += a_[6]*u1/std::pow(h, 3) + a_[7]*u2/std::pow(h, 3);
76  q2 += a_[8]*u1/std::pow(h, 3) + a_[9]*u2/std::pow(h, 3);
77  q3 += a_[7]*u1/std::pow(h, 3) + a_[8]*u2/std::pow(h, 3);
78  case 6:
79  q1 += a_[3]/std::pow(h, 2);
80  q2 += a_[5]/std::pow(h, 2);
81  q3 += a_[4]/std::pow(h, 2);
82  // case 3: nothing to do here
83  case 1:
84  break;
85  default:
86  compadre_kernel_assert_release(false && "curvature polynomial order is greater than 8.");
87 
88  }
89  const double det_g = MetricFactor(a_, h, u1, u2);
90  return (q1*q2-std::pow(q3,2))/(det_g*det_g);
91 
92  }
93 
94  //! Surface curl at any point in the local chart
95  KOKKOS_INLINE_FUNCTION
96  double SurfaceCurlOfScalar(const scratch_vector_type a_, const double h, const double u1, const double u2, int x_pow, int y_pow, const int component) {
97 
98  const double factorial[15] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
99 
100  const double det_g = MetricFactor(a_, h, u1, u2);
101 
102  int partial_direction = -1;
103  if (component == 0) { // comp 0 of surface curl, i.e. dy(p)
104  partial_direction = 1;
105  } else if (component == 1) { // comp 1 of surface curl, i.e. -dx(p)
106  partial_direction = 0;
107  } else {
108  compadre_kernel_assert_release(false);
109  }
110  const double sign_change = (component == 1) ? -1 : 1;
111  int n_x_pow = (partial_direction == 0) ? x_pow-1 : x_pow;
112  int n_y_pow = (partial_direction == 1) ? y_pow-1 : y_pow;
113 
114  double return_val = 0;
115  if (n_x_pow<0 || n_y_pow<0) {
116  return 0;
117  } else {
118  double alphaf = factorial[n_x_pow]*factorial[n_y_pow];
119  return_val = 1./h
120  *std::pow(u1/h,n_x_pow)
121  *std::pow(u2/h,n_y_pow)/alphaf;
122  }
123  return_val *= sign_change;
124  return_val /= sqrt(det_g);
125  return return_val;
126 
127  }
128 
129 } // Compadre namespace
scratch_vector_type
Kokkos::View< double *, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_vector_type
Definition: Compadre_Typedefs.hpp:63
compadre_kernel_assert_release
#define compadre_kernel_assert_release(condition)
compadre_kernel_assert_release is similar to compadre_assert_release, but is a call on the device,...
Definition: Compadre_Typedefs.hpp:154
Compadre::GaussianCurvature
KOKKOS_INLINE_FUNCTION double GaussianCurvature(const scratch_vector_type a_, const double h, const double u1, const double u2)
Gaussian curvature K at any point in the local chart.
Definition: Compadre_Manifold_Functions.hpp:50
Compadre::MetricFactor
KOKKOS_INLINE_FUNCTION double MetricFactor(const scratch_vector_type a_, const double h, const double u1, const double u2)
Metric factor (det(G)) at any point in the local chart.
Definition: Compadre_Manifold_Functions.hpp:6
Compadre::SurfaceCurlOfScalar
KOKKOS_INLINE_FUNCTION double SurfaceCurlOfScalar(const scratch_vector_type a_, const double h, const double u1, const double u2, int x_pow, int y_pow, const int component)
Surface curl at any point in the local chart.
Definition: Compadre_Manifold_Functions.hpp:96
Generated on Thu Aug 24 2023 11:50:34 for Compadre by doxygen 1.9.1