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