1 /// @file
2 /// Test polynomial interpolation transpose from arbitrary points in 1D
3 /// \test Test polynomial interpolation transpose from arbitrary points 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 * c)10 static CeedScalar Eval(CeedScalar x, CeedInt n, const CeedScalar *c) {
11 CeedScalar y = c[n - 1];
12 for (CeedInt i = n - 2; i >= 0; i--) y = y * x + c[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_nodes, x_points, x_point, u, v, u_point, v_point;
19 CeedBasis basis_x, basis_u;
20 const CeedInt p = 5, q = 5, num_points = 4;
21 const CeedScalar c[4] = {1, 2, 3, 4}; // 1 + 2x + 3x^2 + ...
22
23 CeedInit(argv[1], &ceed);
24
25 CeedVectorCreate(ceed, 2, &x);
26 CeedVectorCreate(ceed, p, &x_nodes);
27 CeedVectorCreate(ceed, num_points, &x_points);
28 CeedVectorCreate(ceed, 1, &x_point);
29 CeedVectorCreate(ceed, p, &u);
30 CeedVectorCreate(ceed, num_points, &v);
31 CeedVectorCreate(ceed, p, &u_point);
32 CeedVectorCreate(ceed, 1, &v_point);
33 CeedVectorSetValue(v_point, 1.0);
34
35 // Get nodal coordinates
36 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, p, CEED_GAUSS_LOBATTO, &basis_x);
37 {
38 CeedScalar x_array[2];
39
40 for (CeedInt i = 0; i < 2; i++) x_array[i] = CeedIntPow(-1, i + 1);
41 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
42 }
43 CeedBasisApply(basis_x, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_nodes);
44
45 // Set values of u at nodes
46 {
47 const CeedScalar *x_array;
48 CeedScalar u_array[p];
49
50 CeedVectorGetArrayRead(x_nodes, CEED_MEM_HOST, &x_array);
51 for (CeedInt i = 0; i < p; i++) u_array[i] = Eval(x_array[i], ALEN(c), c);
52 CeedVectorRestoreArrayRead(x_nodes, &x_array);
53 CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, (CeedScalar *)&u_array);
54 }
55
56 // Interpolate to arbitrary points
57 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, q, CEED_GAUSS, &basis_u);
58 {
59 CeedScalar x_array[4] = {-0.33, -0.65, 0.16, 0.99};
60
61 CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
62 }
63 CeedBasisApplyAtPoints(basis_u, 1, &num_points, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_points, u, v);
64
65 for (CeedInt i = 0; i < num_points; i++) {
66 const CeedInt num_point[1] = {1};
67 CeedScalar fx = 0.0;
68 const CeedScalar *x_array, *u_array, *v_array, *u_point_array;
69
70 CeedVectorGetArrayRead(x_points, CEED_MEM_HOST, &x_array);
71 CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array);
72 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
73 CeedVectorSetValue(x_point, x_array[i]);
74 CeedBasisApplyAtPoints(basis_u, 1, num_point, CEED_TRANSPOSE, CEED_EVAL_INTERP, x_point, v_point, u_point);
75 CeedVectorGetArrayRead(u_point, CEED_MEM_HOST, &u_point_array);
76 for (CeedInt j = 0; j < p; j++) fx += u_array[j] * u_point_array[j];
77 if (fabs(v_array[i] - fx) > 100. * CEED_EPSILON) printf("%f != %f = f(%f)\n", v_array[i], fx, x_array[i]);
78 CeedVectorRestoreArrayRead(u_point, &u_point_array);
79 CeedVectorRestoreArrayRead(x_points, &x_array);
80 CeedVectorRestoreArrayRead(u, &u_array);
81 CeedVectorRestoreArrayRead(v, &v_array);
82 }
83
84 CeedVectorDestroy(&x);
85 CeedVectorDestroy(&x_nodes);
86 CeedVectorDestroy(&x_points);
87 CeedVectorDestroy(&x_point);
88 CeedVectorDestroy(&u);
89 CeedVectorDestroy(&v);
90 CeedVectorDestroy(&u_point);
91 CeedVectorDestroy(&v_point);
92 CeedBasisDestroy(&basis_x);
93 CeedBasisDestroy(&basis_u);
94 CeedDestroy(&ceed);
95 return 0;
96 }
97