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