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