1 /// @file
2 /// Test polynomial derivative interpolation in 1D
3 /// \test Test polynomial derivative interpolation in 1D
4 #include <ceed.h>
5 #include <math.h>
6 #include <stdio.h>
7
8 #define ALEN(a) (sizeof(a) / sizeof((a)[0]))
9
Eval(CeedScalar x,CeedInt n,const CeedScalar * p)10 static CeedScalar Eval(CeedScalar x, CeedInt n, const CeedScalar *p) {
11 CeedScalar y = p[n - 1];
12 for (CeedInt i = n - 2; i >= 0; i--) y = y * x + p[i];
13 return y;
14 }
15
main(int argc,char ** argv)16 int main(int argc, char **argv) {
17 Ceed ceed;
18 CeedVector x, x_q, u, u_q;
19 CeedBasis basis_x_lobatto, basis_u_lobatto, basis_x_gauss, basis_u_gauss;
20 CeedInt q = 6;
21 const CeedScalar p[6] = {1, 2, 3, 4, 5, 6}; // 1 + 2x + 3x^2 + ...
22 const CeedScalar dp[5] = {2, 6, 12, 20, 30}; // 2 + 6x + 12x^2 + ...
23
24 CeedInit(argv[1], &ceed);
25
26 CeedVectorCreate(ceed, 2, &x);
27 CeedVectorCreate(ceed, q, &x_q);
28 CeedVectorSetValue(x_q, 0);
29 CeedVectorCreate(ceed, q, &u);
30 CeedVectorSetValue(u, 0);
31 CeedVectorCreate(ceed, q, &u_q);
32
33 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, q, CEED_GAUSS_LOBATTO, &basis_x_lobatto);
34 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, q, q, CEED_GAUSS_LOBATTO, &basis_u_lobatto);
35
36 {
37 CeedScalar x_array[2];
38
39 for (int i = 0; i < 2; i++) x_array[i] = CeedIntPow(-1, i + 1);
40 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
41 }
42 CeedBasisApply(basis_x_lobatto, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_q);
43
44 {
45 const CeedScalar *x_q_array;
46 CeedScalar u_q_array[q];
47
48 CeedVectorGetArrayRead(x_q, CEED_MEM_HOST, &x_q_array);
49 for (CeedInt i = 0; i < q; i++) u_q_array[i] = Eval(x_q_array[i], ALEN(p), p);
50 CeedVectorRestoreArrayRead(x_q, &x_q_array);
51 CeedVectorSetArray(u_q, CEED_MEM_HOST, CEED_COPY_VALUES, u_q_array);
52 }
53
54 // This operation is the identity because the quadrature is collocated
55 CeedBasisApply(basis_u_lobatto, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, u_q, u);
56
57 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, q, CEED_GAUSS, &basis_x_gauss);
58 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, q, q, CEED_GAUSS, &basis_u_gauss);
59
60 CeedBasisApply(basis_x_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_q);
61 CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, u, u_q);
62
63 {
64 const CeedScalar *x_q_array, *u_q_array;
65
66 CeedVectorGetArrayRead(x_q, CEED_MEM_HOST, &x_q_array);
67 CeedVectorGetArrayRead(u_q, CEED_MEM_HOST, &u_q_array);
68 for (CeedInt i = 0; i < q; i++) {
69 CeedScalar px = Eval(x_q_array[i], ALEN(dp), dp);
70 if (fabs(u_q_array[i] - px) > 1000. * CEED_EPSILON) printf("%f != %f = p(%f)\n", u_q_array[i], px, x_q_array[i]);
71 }
72 CeedVectorRestoreArrayRead(x_q, &x_q_array);
73 CeedVectorRestoreArrayRead(u_q, &u_q_array);
74 }
75
76 CeedVectorDestroy(&x);
77 CeedVectorDestroy(&x_q);
78 CeedVectorDestroy(&u);
79 CeedVectorDestroy(&u_q);
80 CeedBasisDestroy(&basis_x_lobatto);
81 CeedBasisDestroy(&basis_u_lobatto);
82 CeedBasisDestroy(&basis_x_gauss);
83 CeedBasisDestroy(&basis_u_gauss);
84 CeedDestroy(&ceed);
85 return 0;
86 }
87