Compadre  1.5.9
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Compadre_Manifold_Functions.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Compadre: COMpatible PArticle Discretization and REmap Toolkit
4 //
5 // Copyright 2018 NTESS and the Compadre contributors.
6 // SPDX-License-Identifier: BSD-2-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 namespace Compadre {
11 
12  //! Metric factor (det(G)) at any point in the local chart
13  KOKKOS_INLINE_FUNCTION
14  double MetricFactor(const scratch_vector_type a_, const double h, const double u1, const double u2) {
15  double det_g = 1.0;
16  double q1 = 0;
17  double q2 = 0;
18  switch (a_.extent(0)) {
19  case 45:
20  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);
21  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);
22  case 36:
23  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);
24  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);
25  case 28:
26  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);
27  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);
28  case 21:
29  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);
30  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);
31  case 15:
32  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);
33  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);
34  case 10:
35  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);
36  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);
37  case 6:
38  q1 += a_[3]*u1/std::pow(h, 2) + a_[4]*u2/std::pow(h, 2);
39  q2 += a_[4]*u1/std::pow(h, 2) + a_[5]*u2/std::pow(h, 2);
40  case 3:
41  q1 += a_[1]/h;
42  q2 += a_[2]/h;
43  case 1:
44  break;
45  default:
46  compadre_kernel_assert_release(false && "curvature polynomial order is greater than 8.");
47 
48  }
49  q1 = std::pow(q1, 2);
50  q2 = std::pow(q2, 2);
51  det_g += q1;
52  det_g += q2;
53  return det_g;
54  }
55 
56  //! Gaussian curvature K at any point in the local chart
57  KOKKOS_INLINE_FUNCTION
58  double GaussianCurvature(const scratch_vector_type a_, const double h, const double u1, const double u2) {
59 
60  double q1=0, q2=0, q3=0;
61  switch (a_.extent(0)) {
62  case 45:
63  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);
64  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);
65  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);
66  case 36:
67  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);
68  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);
69  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);
70  case 28:
71  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);
72  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);
73  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);
74  case 21:
75  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);
76  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);
77  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);
78  case 15:
79  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);
80  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);
81  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);
82  case 10:
83  q1 += a_[6]*u1/std::pow(h, 3) + a_[7]*u2/std::pow(h, 3);
84  q2 += a_[8]*u1/std::pow(h, 3) + a_[9]*u2/std::pow(h, 3);
85  q3 += a_[7]*u1/std::pow(h, 3) + a_[8]*u2/std::pow(h, 3);
86  case 6:
87  q1 += a_[3]/std::pow(h, 2);
88  q2 += a_[5]/std::pow(h, 2);
89  q3 += a_[4]/std::pow(h, 2);
90  // case 3: nothing to do here
91  case 1:
92  break;
93  default:
94  compadre_kernel_assert_release(false && "curvature polynomial order is greater than 8.");
95 
96  }
97  const double det_g = MetricFactor(a_, h, u1, u2);
98  return (q1*q2-std::pow(q3,2))/(det_g*det_g);
99 
100  }
101 
102  //! Surface curl at any point in the local chart
103  KOKKOS_INLINE_FUNCTION
104  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) {
105 
106  const double factorial[15] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
107 
108  const double det_g = MetricFactor(a_, h, u1, u2);
109 
110  int partial_direction = -1;
111  if (component == 0) { // comp 0 of surface curl, i.e. dy(p)
112  partial_direction = 1;
113  } else if (component == 1) { // comp 1 of surface curl, i.e. -dx(p)
114  partial_direction = 0;
115  } else {
116  compadre_kernel_assert_release(false);
117  }
118  const double sign_change = (component == 1) ? -1 : 1;
119  int n_x_pow = (partial_direction == 0) ? x_pow-1 : x_pow;
120  int n_y_pow = (partial_direction == 1) ? y_pow-1 : y_pow;
121 
122  double return_val = 0;
123  if (n_x_pow<0 || n_y_pow<0) {
124  return 0;
125  } else {
126  double alphaf = factorial[n_x_pow]*factorial[n_y_pow];
127  return_val = 1./h
128  *std::pow(u1/h,n_x_pow)
129  *std::pow(u2/h,n_y_pow)/alphaf;
130  }
131  return_val *= sign_change;
132  return_val /= sqrt(det_g);
133  return return_val;
134 
135  }
136 
137 } // Compadre namespace
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:162
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:58
scratch_vector_type
Kokkos::View< double *, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_vector_type
Definition: Compadre_Typedefs.hpp:71
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:104
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:14
Generated on Fri Nov 22 2024 09:25:51 for Compadre by   doxygen 1.8.5