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