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, w;
19 CeedBasis basis_x_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 CeedScalar sum, error, pint[ALEN(p) + 1];
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 CeedVectorCreate(ceed, q, &u_q);
31 CeedVectorSetValue(u_q, 0);
32 CeedVectorCreate(ceed, q, &w);
33 CeedVectorSetValue(w, 0);
34
35 {
36 CeedScalar x_array[2];
37
38 for (int i = 0; i < 2; i++) x_array[i] = CeedIntPow(-1, i + 1);
39 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
40 }
41
42 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, q, CEED_GAUSS_LOBATTO, &basis_x_lobatto);
43
44 CeedBasisApply(basis_x_lobatto, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_q);
45 {
46 const CeedScalar *x_q_array;
47 CeedScalar u_array[q];
48
49 CeedVectorGetArrayRead(x_q, CEED_MEM_HOST, &x_q_array);
50 for (CeedInt i = 0; i < q; i++) u_array[i] = Eval(x_q_array[i], ALEN(p), p);
51 CeedVectorRestoreArrayRead(x_q, &x_q_array);
52 CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, u_array);
53 }
54
55 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, q, CEED_GAUSS, &basis_x_gauss);
56 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, q, q, CEED_GAUSS, &basis_u_gauss);
57
58 CeedBasisApply(basis_x_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_q);
59 CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, u_q);
60 CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, w);
61
62 {
63 const CeedScalar *w_array, *u_q_array;
64
65 CeedVectorGetArrayRead(w, CEED_MEM_HOST, &w_array);
66 CeedVectorGetArrayRead(u_q, CEED_MEM_HOST, &u_q_array);
67 sum = 0;
68 for (CeedInt i = 0; i < q; i++) sum += w_array[i] * u_q_array[i];
69 CeedVectorRestoreArrayRead(w, &w_array);
70 CeedVectorRestoreArrayRead(u_q, &u_q_array);
71 }
72
73 pint[0] = 0;
74 for (CeedInt i = 0; i < (int)ALEN(p); i++) pint[i + 1] = p[i] / (i + 1);
75 error = sum - Eval(1, ALEN(pint), pint) + Eval(-1, ALEN(pint), pint);
76 if (error > 100. * CEED_EPSILON) {
77 // LCOV_EXCL_START
78 printf("Error %e sum %g exact %g\n", error, sum, Eval(1, ALEN(pint), pint) - Eval(-1, ALEN(pint), pint));
79 // LCOV_EXCL_STOP
80 }
81
82 CeedVectorDestroy(&x);
83 CeedVectorDestroy(&x_q);
84 CeedVectorDestroy(&u);
85 CeedVectorDestroy(&u_q);
86 CeedVectorDestroy(&w);
87 CeedBasisDestroy(&basis_x_lobatto);
88 CeedBasisDestroy(&basis_x_gauss);
89 CeedBasisDestroy(&basis_u_gauss);
90 CeedDestroy(&ceed);
91 return 0;
92 }
93