1 /// @file 2 /// Test polynomial interpolation in 1D 3 /// \test Test polynomial interpolation in 1D 4 #include <ceed.h> 5 #include <math.h> 6 7 #define ALEN(a) (sizeof(a) / sizeof((a)[0])) 8 9 static CeedScalar PolyEval(CeedScalar x, CeedInt n, const CeedScalar *p) { 10 CeedScalar y = p[n-1]; 11 for (CeedInt i=n-2; i>=0; i--) y = y*x + p[i]; 12 return y; 13 } 14 15 int main(int argc, char **argv) { 16 Ceed ceed; 17 CeedVector X, X_q, U, U_q, W; 18 CeedBasis basis_x_lobatto, basis_x_gauss, basis_u_gauss; 19 CeedInt Q = 6; 20 const CeedScalar p[6] = {1, 2, 3, 4, 5, 6}; // 1 + 2x + 3x^2 + ... 21 const CeedScalar *xq, *uq, *w; 22 CeedScalar u[Q], x[2], 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 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS_LOBATTO, 36 &basis_x_lobatto); 37 38 for (int i = 0; i < 2; i++) 39 x[i] = CeedIntPow(-1, i+1); 40 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 41 42 CeedBasisApply(basis_x_lobatto, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, X_q); 43 44 CeedVectorGetArrayRead(X_q, CEED_MEM_HOST, &xq); 45 for (CeedInt i=0; i<Q; i++) 46 u[i] = PolyEval(xq[i], ALEN(p), p); 47 CeedVectorRestoreArrayRead(X_q, &xq); 48 CeedVectorSetArray(U, CEED_MEM_HOST, CEED_USE_POINTER, u); 49 50 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS, &basis_x_gauss); 51 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, Q, Q, CEED_GAUSS, &basis_u_gauss); 52 53 CeedBasisApply(basis_x_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, X_q); 54 CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, U, U_q); 55 CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, 56 CEED_VECTOR_NONE, W); 57 58 CeedVectorGetArrayRead(W, CEED_MEM_HOST, &w); 59 CeedVectorGetArrayRead(U_q, CEED_MEM_HOST, &uq); 60 sum = 0; 61 for (CeedInt i=0; i<Q; i++) 62 sum += w[i] * uq[i]; 63 CeedVectorRestoreArrayRead(W, &w); 64 CeedVectorRestoreArrayRead(U_q, &uq); 65 66 pint[0] = 0; 67 for (CeedInt i=0; i<(int)ALEN(p); i++) 68 pint[i+1] = p[i] / (i+1); 69 error = sum - PolyEval(1, ALEN(pint), pint) + PolyEval(-1, ALEN(pint), pint); 70 if (error > 1.E-10) 71 // LCOV_EXCL_START 72 printf("Error %e sum %g exact %g\n", error, sum, 73 PolyEval(1, ALEN(pint), pint) - PolyEval(-1, ALEN(pint), pint)); 74 // LCOV_EXCL_STOP 75 76 CeedVectorDestroy(&X); 77 CeedVectorDestroy(&X_q); 78 CeedVectorDestroy(&U); 79 CeedVectorDestroy(&U_q); 80 CeedVectorDestroy(&W); 81 CeedBasisDestroy(&basis_x_lobatto); 82 CeedBasisDestroy(&basis_x_gauss); 83 CeedBasisDestroy(&basis_u_gauss); 84 CeedDestroy(&ceed); 85 return 0; 86 } 87