Compadre 1.5.5
Loading...
Searching...
No Matches
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
8KOKKOS_INLINE_FUNCTION
9double device_max(double d1, double d2) {
10 return (d1 > d2) ? d1 : d2;
11}
12
13KOKKOS_INLINE_FUNCTION
14double 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
28KOKKOS_INLINE_FUNCTION
29double 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
68KOKKOS_INLINE_FUNCTION
69void 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
107KOKKOS_INLINE_FUNCTION
108double 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
147KOKKOS_INLINE_FUNCTION
148void 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
187KOKKOS_INLINE_FUNCTION
188double 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
200KOKKOS_INLINE_FUNCTION
201double 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
213KOKKOS_INLINE_FUNCTION
214double 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
242KOKKOS_INLINE_FUNCTION
243double 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
266KOKKOS_INLINE_FUNCTION
267double 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
284KOKKOS_INLINE_FUNCTION
285double 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
301KOKKOS_INLINE_FUNCTION
302double 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
340KOKKOS_INLINE_FUNCTION
341double 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
364KOKKOS_INLINE_FUNCTION
365double 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
403KOKKOS_INLINE_FUNCTION
404double 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) {
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 */
425int 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
KOKKOS_INLINE_FUNCTION double divergenceTestSolution(double x, double y, double z, int dimension)
KOKKOS_INLINE_FUNCTION double divfreeTestSolution_single_polynomial(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_span_basis(double x, double y, double z, int component, int dimension, int exact_order)
KOKKOS_INLINE_FUNCTION double trueSolution(double x, double y, double z, int order, int dimension)
KOKKOS_INLINE_FUNCTION void trueGradient(double *ans, double x, double y, double z, int order, int dimension)
int main(int argc, char *args[])
Standard GMLS Example.
Definition: GMLS_Device.cpp:29
KOKKOS_INLINE_FUNCTION double divergenceTestSamples(double x, double y, double z, int component, int dimension)
KOKKOS_INLINE_FUNCTION double curldivfreeTestSolution_single_polynomial(double x, double y, double z, int component, int dimension)
KOKKOS_INLINE_FUNCTION double curlcurldivfreeTestSolution_single_polynomial(double x, double y, double z, int component, int dimension)
KOKKOS_INLINE_FUNCTION double gradientdivfreeTestSolution_single_polynomial(double x, double y, double z, int component, int dimension)
KOKKOS_INLINE_FUNCTION double divfreeTestSolution_span_basis(double x, double y, double z, int component, int dimension, int exact_order)
KOKKOS_INLINE_FUNCTION double trueDivergence(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 void trueHessian(double *ans, double x, double y, double z, int order, int dimension)
KOKKOS_INLINE_FUNCTION double device_max(double d1, double d2)
KOKKOS_INLINE_FUNCTION double trueLaplacian(double x, double y, double z, int order, int dimension)
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 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 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_...
@ DivergenceFreeVectorTaylorPolynomial
Divergence-free vector polynomial basis.