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