1 /// @file 2 /// Test polynomial interpolation in 1D 3 /// \test Test polynomial 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, w; 19 CeedBasis basis_x_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 CeedScalar sum, error, pint[ALEN(p) + 1]; 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 CeedVectorCreate(ceed, q, &u_q); 31 CeedVectorSetValue(u_q, 0); 32 CeedVectorCreate(ceed, q, &w); 33 CeedVectorSetValue(w, 0); 34 35 { 36 CeedScalar x_array[2]; 37 38 for (int i = 0; i < 2; i++) x_array[i] = CeedIntPow(-1, i + 1); 39 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); 40 } 41 42 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, q, CEED_GAUSS_LOBATTO, &basis_x_lobatto); 43 44 CeedBasisApply(basis_x_lobatto, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_q); 45 { 46 const CeedScalar *x_q_array; 47 CeedScalar u_array[q]; 48 49 CeedVectorGetArrayRead(x_q, CEED_MEM_HOST, &x_q_array); 50 for (CeedInt i = 0; i < q; i++) u_array[i] = Eval(x_q_array[i], ALEN(p), p); 51 CeedVectorRestoreArrayRead(x_q, &x_q_array); 52 CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, u_array); 53 } 54 55 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, q, CEED_GAUSS, &basis_x_gauss); 56 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, q, q, CEED_GAUSS, &basis_u_gauss); 57 58 CeedBasisApply(basis_x_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_q); 59 CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, u_q); 60 CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, w); 61 62 { 63 const CeedScalar *w_array, *u_q_array; 64 65 CeedVectorGetArrayRead(w, CEED_MEM_HOST, &w_array); 66 CeedVectorGetArrayRead(u_q, CEED_MEM_HOST, &u_q_array); 67 sum = 0; 68 for (CeedInt i = 0; i < q; i++) sum += w_array[i] * u_q_array[i]; 69 CeedVectorRestoreArrayRead(w, &w_array); 70 CeedVectorRestoreArrayRead(u_q, &u_q_array); 71 } 72 73 pint[0] = 0; 74 for (CeedInt i = 0; i < (int)ALEN(p); i++) pint[i + 1] = p[i] / (i + 1); 75 error = sum - Eval(1, ALEN(pint), pint) + Eval(-1, ALEN(pint), pint); 76 if (error > 100. * CEED_EPSILON) { 77 // LCOV_EXCL_START 78 printf("Error %e sum %g exact %g\n", error, sum, Eval(1, ALEN(pint), pint) - Eval(-1, ALEN(pint), pint)); 79 // LCOV_EXCL_STOP 80 } 81 82 CeedVectorDestroy(&x); 83 CeedVectorDestroy(&x_q); 84 CeedVectorDestroy(&u); 85 CeedVectorDestroy(&u_q); 86 CeedVectorDestroy(&w); 87 CeedBasisDestroy(&basis_x_lobatto); 88 CeedBasisDestroy(&basis_x_gauss); 89 CeedBasisDestroy(&basis_u_gauss); 90 CeedDestroy(&ceed); 91 return 0; 92 } 93