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 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 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 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