1 /// @file
2 /// Test polynomial interpolation in 1D
3 /// \test Test polynomial 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
23 CeedInit(argv[1], &ceed);
24
25 CeedVectorCreate(ceed, 2, &x);
26 CeedVectorCreate(ceed, q, &x_q);
27 CeedVectorSetValue(x_q, 0);
28 CeedVectorCreate(ceed, q, &u);
29 CeedVectorSetValue(u, 0);
30 CeedVectorCreate(ceed, q, &u_q);
31
32 {
33 CeedScalar x_array[2];
34
35 for (CeedInt i = 0; i < 2; i++) x_array[i] = CeedIntPow(-1, i + 1);
36 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
37 }
38
39 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, q, CEED_GAUSS_LOBATTO, &basis_x_lobatto);
40 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, q, q, CEED_GAUSS_LOBATTO, &basis_u_lobatto);
41
42 CeedBasisApply(basis_x_lobatto, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_q);
43 {
44 const CeedScalar *x_q_array;
45 CeedScalar u_array[q];
46
47 CeedVectorGetArrayRead(x_q, CEED_MEM_HOST, &x_q_array);
48 for (CeedInt i = 0; i < q; i++) u_array[i] = Eval(x_q_array[i], ALEN(p), p);
49 CeedVectorRestoreArrayRead(x_q, &x_q_array);
50 CeedVectorSetArray(u_q, CEED_MEM_HOST, CEED_COPY_VALUES, (CeedScalar *)&u_array);
51 }
52
53 // This operation is the identity because the quadrature is collocated
54 CeedBasisApply(basis_u_lobatto, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, u_q, u);
55
56 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, q, CEED_GAUSS, &basis_x_gauss);
57 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, q, q, CEED_GAUSS, &basis_u_gauss);
58
59 CeedBasisApply(basis_x_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_q);
60 CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, u_q);
61
62 {
63 const CeedScalar *x_q_array, *u_q_array;
64
65 CeedVectorGetArrayRead(x_q, CEED_MEM_HOST, &x_q_array);
66 CeedVectorGetArrayRead(u_q, CEED_MEM_HOST, &u_q_array);
67 for (CeedInt i = 0; i < q; i++) {
68 CeedScalar px = Eval(x_q_array[i], ALEN(p), p);
69 if (fabs(u_q_array[i] - px) > 100. * CEED_EPSILON) printf("%f != %f = p(%f)\n", u_q_array[i], px, x_q_array[i]);
70 }
71 CeedVectorRestoreArrayRead(x_q, &x_q_array);
72 CeedVectorRestoreArrayRead(u_q, &u_q_array);
73 }
74
75 CeedVectorDestroy(&x);
76 CeedVectorDestroy(&x_q);
77 CeedVectorDestroy(&u);
78 CeedVectorDestroy(&u_q);
79 CeedBasisDestroy(&basis_x_lobatto);
80 CeedBasisDestroy(&basis_u_lobatto);
81 CeedBasisDestroy(&basis_x_gauss);
82 CeedBasisDestroy(&basis_u_gauss);
83 CeedDestroy(&ceed);
84 return 0;
85 }
86