Compadre  1.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
GMLS_Tutorial.hpp
Go to the documentation of this file.
1 #ifndef _GMLS_TUTORIAL_HPP_
2 #define _GMLS_TUTORIAL_HPP_
3 
4 #include <Kokkos_Core.hpp>
6 #include <Compadre_GMLS.hpp>
7 
8 KOKKOS_INLINE_FUNCTION
9 double device_max(double d1, double d2) {
10  return (d1 > d2) ? d1 : d2;
11 }
12 
13 KOKKOS_INLINE_FUNCTION
14 double trueSolution(double x, double y, double z, int order, int dimension) {
15  double ans = 0;
16  for (int i=0; i<order+1; i++) {
17  for (int j=0; j<order+1; j++) {
18  for (int k=0; k<order+1; k++) {
19  if (i+j+k <= order) {
20  ans += std::pow(x,i)*std::pow(y,j)*std::pow(z,k);
21  }
22  }
23  }
24  }
25  return ans;
26 }
27 
28 KOKKOS_INLINE_FUNCTION
29 double trueLaplacian(double x, double y, double z, int order, int dimension) {
30  double ans = 0;
31  for (int i=0; i<order+1; i++) {
32  for (int j=0; j<order+1; j++) {
33  for (int k=0; k<order+1; k++) {
34  if (i+j+k <= order) {
35  ans += device_max(0,i-1)*device_max(0,i)*
36  std::pow(x,device_max(0,i-2))*std::pow(y,j)*std::pow(z,k);
37  }
38  }
39  }
40  }
41  if (dimension>1) {
42  for (int i=0; i<order+1; i++) {
43  for (int j=0; j<order+1; j++) {
44  for (int k=0; k<order+1; k++) {
45  if (i+j+k <= order) {
46  ans += device_max(0,j-1)*device_max(0,j)*
47  std::pow(x,i)*std::pow(y,device_max(0,j-2))*std::pow(z,k);
48  }
49  }
50  }
51  }
52  }
53  if (dimension>2) {
54  for (int i=0; i<order+1; i++) {
55  for (int j=0; j<order+1; j++) {
56  for (int k=0; k<order+1; k++) {
57  if (i+j+k <= order) {
58  ans += device_max(0,k-1)*device_max(0,k)*
59  std::pow(x,i)*std::pow(y,j)*std::pow(z,device_max(0,k-2));
60  }
61  }
62  }
63  }
64  }
65  return ans;
66 }
67 
68 KOKKOS_INLINE_FUNCTION
69 void trueGradient(double* ans, double x, double y, double z, int order, int dimension) {
70 
71  for (int i=0; i<order+1; i++) {
72  for (int j=0; j<order+1; j++) {
73  for (int k=0; k<order+1; k++) {
74  if (i+j+k <= order) {
75  ans[0] += device_max(0,i)*
76  std::pow(x,device_max(0,i-1))*std::pow(y,j)*std::pow(z,k);
77  }
78  }
79  }
80  }
81  if (dimension>1) {
82  for (int i=0; i<order+1; i++) {
83  for (int j=0; j<order+1; j++) {
84  for (int k=0; k<order+1; k++) {
85  if (i+j+k <= order) {
86  ans[1] += device_max(0,j)*
87  std::pow(x,i)*std::pow(y,device_max(0,j-1))*std::pow(z,k);
88  }
89  }
90  }
91  }
92  }
93  if (dimension>2) {
94  for (int i=0; i<order+1; i++) {
95  for (int j=0; j<order+1; j++) {
96  for (int k=0; k<order+1; k++) {
97  if (i+j+k <= order) {
98  ans[2] += device_max(0,k)*
99  std::pow(x,i)*std::pow(y,j)*std::pow(z,device_max(0,k-1));
100  }
101  }
102  }
103  }
104  }
105 }
106 
107 KOKKOS_INLINE_FUNCTION
108 double trueDivergence(double x, double y, double z, int order, int dimension) {
109  double ans = 0;
110  for (int i=0; i<order+1; i++) {
111  for (int j=0; j<order+1; j++) {
112  for (int k=0; k<order+1; k++) {
113  if (i+j+k <= order) {
114  ans += device_max(0,i)*
115  std::pow(x,device_max(0,i-1))*std::pow(y,j)*std::pow(z,k);
116  }
117  }
118  }
119  }
120  if (dimension>1) {
121  for (int i=0; i<order+1; i++) {
122  for (int j=0; j<order+1; j++) {
123  for (int k=0; k<order+1; k++) {
124  if (i+j+k <= order) {
125  ans += device_max(0,j)*
126  std::pow(x,i)*std::pow(y,device_max(0,j-1))*std::pow(z,k);
127  }
128  }
129  }
130  }
131  }
132  if (dimension>2) {
133  for (int i=0; i<order+1; i++) {
134  for (int j=0; j<order+1; j++) {
135  for (int k=0; k<order+1; k++) {
136  if (i+j+k <= order) {
137  ans += device_max(0,k)*
138  std::pow(x,i)*std::pow(y,j)*std::pow(z,device_max(0,k-1));
139  }
140  }
141  }
142  }
143  }
144  return ans;
145 }
146 
147 KOKKOS_INLINE_FUNCTION
148 void trueHessian(double* ans, double x, double y, double z, int order, int dimension) {
149  for (int i=0; i<order+1; i++) {
150  for (int j=0; j<order+1; j++) {
151  for (int k=0; k<order+1; k++) {
152  if (i+j+k <= order) {
153  // XX
154  ans[0] += device_max(0,i)*device_max(0,i-1)*
155  std::pow(x,device_max(0,i-2))*std::pow(y,j)*std::pow(z,k);
156  if (dimension>1) {
157  // XY
158  ans[1] += device_max(0,i)*device_max(0,j)*
159  std::pow(x,device_max(0,i-1))*std::pow(y,device_max(0,j-1))*std::pow(z,k);
160  // YX = XY
161  ans[1*dimension+0] = ans[1];
162  // YY
163  ans[1*dimension+1] += device_max(0,j)*device_max(0,j-1)*
164  std::pow(x,i)*std::pow(y,device_max(0,j-2))*std::pow(z,k);
165  }
166  if (dimension>2) {
167  // XZ
168  ans[2] += device_max(0,i)*device_max(0,k)*
169  std::pow(x,device_max(0,i-1))*std::pow(y,j)*std::pow(z,device_max(0,k-1));
170  // YZ
171  ans[1*dimension+2] += device_max(0,j)*device_max(0,k)*
172  std::pow(x,i)*std::pow(y,device_max(0,j-1))*std::pow(z,device_max(0,k-1));
173  // ZX = XZ
174  ans[2*dimension+0] = ans[2];
175  // ZY = YZ
176  ans[2*dimension+1] = ans[1*dimension+2];
177  // ZZ
178  ans[2*dimension+2] += device_max(0,k)*device_max(0,k-1)*
179  std::pow(x,i)*std::pow(y,j)*std::pow(z,device_max(0,k-2));
180  }
181  }
182  }
183  }
184  }
185 }
186 
187 KOKKOS_INLINE_FUNCTION
188 double divergenceTestSamples(double x, double y, double z, int component, int dimension) {
189  // solution can be captured exactly by at least 2rd order
190  switch (component) {
191  case 0:
192  return x*x + y*y - z*z;
193  case 1:
194  return 2*x +3*y + 4*z;
195  default:
196  return -21*x*y + 3*z*x*y + 4*z;
197  }
198 }
199 
200 KOKKOS_INLINE_FUNCTION
201 double divergenceTestSolution(double x, double y, double z, int dimension) {
202  switch (dimension) {
203  case 1:
204  // returns divergence of divergenceTestSamples
205  return 2*x;
206  case 2:
207  return 2*x + 3;
208  default:
209  return 2*x + 3 + 3*x*y + 4;
210  }
211 }
212 
213 KOKKOS_INLINE_FUNCTION
214 double curlTestSolution(double x, double y, double z, int component, int dimension) {
215  if (dimension==3) {
216  // returns curl of divergenceTestSamples
217  switch (component) {
218  case 0:
219  // u3,y- u2,z
220  return (-21*x + 3*z*x) - 4;
221  case 1:
222  // -u3,x + u1,z
223  return (21*y - 3*z*y) - 2*z;
224  default:
225  // u2,x - u1,y
226  return 2 - 2*y;
227  }
228  } else if (dimension==2) {
229  switch (component) {
230  case 0:
231  // u2,y
232  return 3;
233  default:
234  // -u1,x
235  return -2*x;
236  }
237  } else {
238  return 0;
239  }
240 }
241 
242 KOKKOS_INLINE_FUNCTION
243 double divfreeTestSolution_single_polynomial(double x, double y, double z, int component, int dimension) {
244  if (dimension==3) {
245  // returns divfreeTestSamples
246  switch (component) {
247  case 0:
248  return 6.0*x*x*y - 9.0*x*y + 7.0*x*z*z + 6.0*y*y*z;
249  case 1:
250  return 10.0*x*x*z - 7.0*y*z*z - 6.0*x*y*y;
251  default:
252  return -2.0*x*x*x + 9.0*x*y*y + 9.0*y*z;
253  }
254  } else if (dimension==2) {
255  switch (component) {
256  case 0:
257  return 6.0*x*x*y;
258  default:
259  return -6.0*x*y*y;
260  }
261  } else {
262  return 0.0;
263  }
264 }
265 
266 KOKKOS_INLINE_FUNCTION
267 double divfreeTestSolution_span_basis(double x, double y, double z, int component, int dimension, int exact_order) {
268  double val = 0.0;
270  if (dimension==3) {
271  for (int i=0; i<NP; ++i) {
273  val += basis_i[component];
274  }
275  } else {
276  for (int i=0; i<NP; ++i) {
278  val += basis_i[component];
279  }
280  }
281  return val;
282 }
283 
284 KOKKOS_INLINE_FUNCTION
285 double curldivfreeTestSolution_single_polynomial(double x, double y, double z, int component, int dimension) {
286  if (dimension==3) {
287  // returns curl of divergenceTestSamples
288  switch (component) {
289  case 0:
290  return -10.0*x*x + 18.0*x*y + 14.0*y*z + 9.0*z;
291  case 1:
292  return 6.0*x*x + 14.0*x*z - 3.0*y*y;
293  default:
294  return -6.0*x*x + 20.0*x*z + 9.0*x - 6.0*y*y - 12.0*y*z;
295  }
296  } else {
297  return 0;
298  }
299 }
300 
301 KOKKOS_INLINE_FUNCTION
302 double curldivfreeTestSolution_span_basis(double x, double y, double z, int component, int dimension, int exact_order) {
303  double val = 0.0;
305  if (dimension==3) {
306  if (component==0) {
307  for (int i=0; i<NP; ++i) {
310  val += grad_y[2] - grad_z[1];
311  }
312  } else if (component==1) {
313  for (int i=0; i<NP; ++i) {
316  val += -grad_x[2] + grad_z[0];
317  }
318  } else if (component==2) {
319  for (int i=0; i<NP; ++i) {
322  val += grad_x[1] - grad_y[0];
323  }
324  }
325  return val;
326  } else if (dimension==2) {
327  if (component==0) {
328  for (int i=0; i<NP; ++i) {
331  val += grad_x[1] - grad_y[0];
332  }
333  return val;
334  } else return 0;
335  } else {
336  return 0;
337  }
338 }
339 
340 KOKKOS_INLINE_FUNCTION
341 double curlcurldivfreeTestSolution_single_polynomial(double x, double y, double z, int component, int dimension) {
342  if (dimension==3) {
343  // returns curl of divergenceTestSamples
344  switch (component) {
345  case 0:
346  return -14.0*x - 12.0*y - 12.0*z;
347  case 1:
348  return 12.0*x + 14.0*y - 20.0*z;
349  default:
350  return -6.0*x;
351  }
352  } else if (dimension==2) {
353  switch (component) {
354  case 0:
355  return -12.0*y;
356  default:
357  return 12.0*x;
358  }
359  } else {
360  return 0;
361  }
362 }
363 
364 KOKKOS_INLINE_FUNCTION
365 double gradientdivfreeTestSolution_single_polynomial(double x, double y, double z, int component, int dimension) {
366  if (dimension==3) {
367  switch (component) {
368  case 0:
369  return 12.0*x*y - 9.0*y + 7.0*z*z;
370  case 1:
371  return 6.0*x*x - 9.0*x + 12.0*y*z;
372  case 2:
373  return 14.0*x*z + 6.0*y*y;
374  case 3:
375  return 20.0*x*z - 6.0*y*y;
376  case 4:
377  return -7.0*z*z - 12.0*x*y;
378  case 5:
379  return 10.0*x*x - 14.0*y*z;
380  case 6:
381  return -6.0*x*x + 9.0*y*y;
382  case 7:
383  return 18.0*x*y + 9.0*z;
384  case 8:
385  return 9.0*y;
386  }
387  }
388  if (dimension==2) {
389  switch (component) {
390  case 0:
391  return 12.0*x*y ;
392  case 1:
393  return 6.0*x*x;
394  case 2:
395  return -6.0*y*y;
396  case 3:
397  return -12.0*x*y;
398  }
399  }
400  return 0.0;
401 }
402 
403 KOKKOS_INLINE_FUNCTION
404 double gradientdivfreeTestSolution_span_basis(double x, double y, double z, int input_component, int output_component, int dimension, int exact_order) {
405  double val = 0.0;
407  if (dimension==3) {
408  for (int i=0; i<NP; ++i) {
409  Compadre::XYZ basis_i = Compadre::DivergenceFreePolynomialBasis::evaluatePartialDerivative(i, output_component, 1.0, x, y, z);
410  val += basis_i[input_component];
411  }
412  } else {
413  for (int i=0; i<NP; ++i) {
415  val += basis_i[input_component];
416  }
417  }
418  return val;
419 }
420 
421 /** Standard GMLS Example
422  *
423  * Exercises GMLS operator evaluation with data over various orders and numbers of targets for targets including point evaluation, Laplacian, divergence, curl, and gradient.
424  */
425 int main (int argc, char* args[]);
426 
427 /**
428  * \example "GMLS Tutorial" based on GMLS_Device.cpp
429  * \section ex GMLS Example with Device Views
430  *
431  * This tutorial sets up a batch of GMLS problems, solves the minimization problems, and applies the coefficients produced to data.
432  *
433  * \section ex1a Parse Command Line Arguments
434  * \snippet GMLS_Device.cpp Parse Command Line Arguments
435  *
436  * \section ex1b Setting Up The Point Cloud
437  * \snippet GMLS_Device.cpp Setting Up The Point Cloud
438  *
439  * \section ex1c Performing Neighbor Search
440  * \snippet GMLS_Device.cpp Performing Neighbor Search
441  *
442  * \section ex2 Creating The Data
443  * \snippet GMLS_Device.cpp Creating The Data
444  *
445  * \section ex3 Setting Up The GMLS Object
446  * \snippet GMLS_Device.cpp Setting Up The GMLS Object
447  *
448  * \section ex4 Apply GMLS Alphas To Data
449  * \snippet GMLS_Device.cpp Apply GMLS Alphas To Data
450  *
451  * \section ex5 Check That Solutions Are Correct
452  * \snippet GMLS_Device.cpp Check That Solutions Are Correct
453  *
454  * \section ex6 Finalize Program
455  * \snippet GMLS_Device.cpp Finalize Program
456  */
457 
458 #endif
Divergence-free vector polynomial basis.
KOKKOS_INLINE_FUNCTION double divfreeTestSolution_span_basis(double x, double y, double z, int component, int dimension, int exact_order)
KOKKOS_INLINE_FUNCTION double gradientdivfreeTestSolution_single_polynomial(double x, double y, double z, int component, int dimension)
KOKKOS_INLINE_FUNCTION void evaluate(const member_type &teamMember, double *delta, double *workspace, const int dimension, const int max_degree, const int component, const double h, const double x, const double y, const double z, const int starting_order=0, const double weight_of_original_value=0.0, const double weight_of_new_value=1.0)
Evaluates the divergence-free polynomial basis delta[j] = weight_of_original_value * delta[j] + weigh...
KOKKOS_INLINE_FUNCTION void trueHessian(double *ans, double x, double y, double z, int order, int dimension)
KOKKOS_INLINE_FUNCTION double curldivfreeTestSolution_span_basis(double x, double y, double z, int component, int dimension, int exact_order)
KOKKOS_INLINE_FUNCTION double curlcurldivfreeTestSolution_single_polynomial(double x, double y, double z, int component, int dimension)
KOKKOS_INLINE_FUNCTION double device_max(double d1, double d2)
int main(int argc, char *args[])
[Parse Command Line Arguments]
Definition: GMLS_Device.cpp:29
static KOKKOS_INLINE_FUNCTION int getNP(const int m, const int dimension=3, const ReconstructionSpace r_space=ReconstructionSpace::ScalarTaylorPolynomial)
Returns size of the basis for a given polynomial order and dimension General to dimension 1...
KOKKOS_INLINE_FUNCTION void trueGradient(double *ans, double x, double y, double z, int order, int dimension)
KOKKOS_INLINE_FUNCTION double divfreeTestSolution_single_polynomial(double x, double y, double z, int component, int dimension)
KOKKOS_INLINE_FUNCTION void evaluatePartialDerivative(const member_type &teamMember, double *delta, double *workspace, const int dimension, const int max_degree, const int component, const int partial_direction, const double h, const double x, const double y, const double z, const int starting_order=0, const double weight_of_original_value=0.0, const double weight_of_new_value=1.0)
Evaluates the first partial derivatives of the divergence-free polynomial basis delta[j] = weight_of_...
KOKKOS_INLINE_FUNCTION double divergenceTestSolution(double x, double y, double z, int dimension)
KOKKOS_INLINE_FUNCTION double trueLaplacian(double x, double y, double z, int order, int dimension)
KOKKOS_INLINE_FUNCTION double curlTestSolution(double x, double y, double z, int component, int dimension)
KOKKOS_INLINE_FUNCTION double divergenceTestSamples(double x, double y, double z, int component, int dimension)
KOKKOS_INLINE_FUNCTION double gradientdivfreeTestSolution_span_basis(double x, double y, double z, int input_component, int output_component, int dimension, int exact_order)
KOKKOS_INLINE_FUNCTION double curldivfreeTestSolution_single_polynomial(double x, double y, double z, int component, int dimension)
KOKKOS_INLINE_FUNCTION double trueSolution(double x, double y, double z, int order, int dimension)
KOKKOS_INLINE_FUNCTION double trueDivergence(double x, double y, double z, int order, int dimension)