1 /// @file 2 /// Test polynomial derivative interpolation in 1D 3 /// \test Test polynomial derivative interpolation 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 *p) { 11 CeedScalar y = p[n - 1]; 12 for (CeedInt i = n - 2; i >= 0; i--) y = y * x + p[i]; 13 return y; 14 } 15 16 int main(int argc, char **argv) { 17 Ceed ceed; 18 CeedVector x, x_q, u, u_q; 19 CeedBasis basis_x_lobatto, basis_u_lobatto, basis_x_gauss, basis_u_gauss; 20 CeedInt q = 6; 21 const CeedScalar p[6] = {1, 2, 3, 4, 5, 6}; // 1 + 2x + 3x^2 + ... 22 const CeedScalar dp[5] = {2, 6, 12, 20, 30}; // 2 + 6x + 12x^2 + ... 23 24 CeedInit(argv[1], &ceed); 25 26 CeedVectorCreate(ceed, 2, &x); 27 CeedVectorCreate(ceed, q, &x_q); 28 CeedVectorSetValue(x_q, 0); 29 CeedVectorCreate(ceed, q, &u); 30 CeedVectorSetValue(u, 0); 31 CeedVectorCreate(ceed, q, &u_q); 32 33 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, q, CEED_GAUSS_LOBATTO, &basis_x_lobatto); 34 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, q, q, CEED_GAUSS_LOBATTO, &basis_u_lobatto); 35 36 { 37 CeedScalar x_array[2]; 38 39 for (int i = 0; i < 2; i++) x_array[i] = CeedIntPow(-1, i + 1); 40 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); 41 } 42 CeedBasisApply(basis_x_lobatto, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_q); 43 44 { 45 const CeedScalar *x_q_array; 46 CeedScalar u_q_array[q]; 47 48 CeedVectorGetArrayRead(x_q, CEED_MEM_HOST, &x_q_array); 49 for (CeedInt i = 0; i < q; i++) u_q_array[i] = Eval(x_q_array[i], ALEN(p), p); 50 CeedVectorRestoreArrayRead(x_q, &x_q_array); 51 CeedVectorSetArray(u_q, CEED_MEM_HOST, CEED_COPY_VALUES, u_q_array); 52 } 53 54 // This operation is the identity because the quadrature is collocated 55 CeedBasisApply(basis_u_lobatto, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, u_q, u); 56 57 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, q, CEED_GAUSS, &basis_x_gauss); 58 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, q, q, CEED_GAUSS, &basis_u_gauss); 59 60 CeedBasisApply(basis_x_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_q); 61 CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, u, u_q); 62 63 { 64 const CeedScalar *x_q_array, *u_q_array; 65 66 CeedVectorGetArrayRead(x_q, CEED_MEM_HOST, &x_q_array); 67 CeedVectorGetArrayRead(u_q, CEED_MEM_HOST, &u_q_array); 68 for (CeedInt i = 0; i < q; i++) { 69 CeedScalar px = Eval(x_q_array[i], ALEN(dp), dp); 70 if (fabs(u_q_array[i] - px) > 1000. * CEED_EPSILON) printf("%f != %f = p(%f)\n", u_q_array[i], px, x_q_array[i]); 71 } 72 CeedVectorRestoreArrayRead(x_q, &x_q_array); 73 CeedVectorRestoreArrayRead(u_q, &u_q_array); 74 } 75 76 CeedVectorDestroy(&x); 77 CeedVectorDestroy(&x_q); 78 CeedVectorDestroy(&u); 79 CeedVectorDestroy(&u_q); 80 CeedBasisDestroy(&basis_x_lobatto); 81 CeedBasisDestroy(&basis_u_lobatto); 82 CeedBasisDestroy(&basis_x_gauss); 83 CeedBasisDestroy(&basis_u_gauss); 84 CeedDestroy(&ceed); 85 return 0; 86 } 87