14411cf47Sjeremylt /// @file 24411cf47Sjeremylt /// Test polynomial interpolation in 1D 34411cf47Sjeremylt /// \test Test polynomial interpolation in 1D 457c64913Sjeremylt #include <ceed.h> 557c64913Sjeremylt #include <math.h> 657c64913Sjeremylt 757c64913Sjeremylt #define ALEN(a) (sizeof(a) / sizeof((a)[0])) 857c64913Sjeremylt 957c64913Sjeremylt static CeedScalar PolyEval(CeedScalar x, CeedInt n, const CeedScalar *p) { 1057c64913Sjeremylt CeedScalar y = p[n-1]; 1157c64913Sjeremylt for (CeedInt i=n-2; i>=0; i--) y = y*x + p[i]; 1257c64913Sjeremylt return y; 1357c64913Sjeremylt } 1457c64913Sjeremylt 1557c64913Sjeremylt int main(int argc, char **argv) { 1657c64913Sjeremylt Ceed ceed; 17*aedaa0e5Sjeremylt CeedVector X, Xq, U, Uq, W; 1857c64913Sjeremylt CeedBasis bxl, bxg, bug; 1957c64913Sjeremylt CeedInt Q = 6; 20*aedaa0e5Sjeremylt const CeedScalar p[6] = {1, 2, 3, 4, 5, 6}; // 1 + 2x + 3x^2 + ... 21*aedaa0e5Sjeremylt const CeedScalar *xq, *uq, *w; 22*aedaa0e5Sjeremylt CeedScalar u[Q], x[2], sum, error, pint[ALEN(p)+1]; 2357c64913Sjeremylt 2457c64913Sjeremylt CeedInit(argv[1], &ceed); 25*aedaa0e5Sjeremylt 26*aedaa0e5Sjeremylt CeedVectorCreate(ceed, 2, &X); 27*aedaa0e5Sjeremylt CeedVectorCreate(ceed, Q, &Xq); 28*aedaa0e5Sjeremylt CeedVectorSetValue(Xq, 0); 29*aedaa0e5Sjeremylt CeedVectorCreate(ceed, Q, &U); 30*aedaa0e5Sjeremylt CeedVectorCreate(ceed, Q, &Uq); 31*aedaa0e5Sjeremylt CeedVectorSetValue(Uq, 0); 32*aedaa0e5Sjeremylt CeedVectorCreate(ceed, Q, &W); 33*aedaa0e5Sjeremylt CeedVectorSetValue(W, 0); 34*aedaa0e5Sjeremylt 3557c64913Sjeremylt CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS_LOBATTO, &bxl); 36*aedaa0e5Sjeremylt 37*aedaa0e5Sjeremylt for (int i = 0; i < 2; i++) x[i] = CeedIntPow(-1, i+1); 38*aedaa0e5Sjeremylt CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, (CeedScalar *)&x); 39*aedaa0e5Sjeremylt 40*aedaa0e5Sjeremylt CeedBasisApply(bxl, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, Xq); 41*aedaa0e5Sjeremylt 42*aedaa0e5Sjeremylt CeedVectorGetArrayRead(Xq, CEED_MEM_HOST, &xq); 4357c64913Sjeremylt for (CeedInt i=0; i<Q; i++) u[i] = PolyEval(xq[i], ALEN(p), p); 44*aedaa0e5Sjeremylt CeedVectorRestoreArrayRead(Xq, &xq); 45*aedaa0e5Sjeremylt CeedVectorSetArray(U, CEED_MEM_HOST, CEED_USE_POINTER, (CeedScalar *)&u); 4657c64913Sjeremylt 4757c64913Sjeremylt CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS, &bxg); 4857c64913Sjeremylt CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, Q, Q, CEED_GAUSS, &bug); 49*aedaa0e5Sjeremylt 50*aedaa0e5Sjeremylt CeedBasisApply(bxg, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, Xq); 51*aedaa0e5Sjeremylt CeedBasisApply(bug, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, U, Uq); 52*aedaa0e5Sjeremylt CeedBasisApply(bug, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, NULL, W); 53*aedaa0e5Sjeremylt 54*aedaa0e5Sjeremylt CeedVectorGetArrayRead(W, CEED_MEM_HOST, &w); 55*aedaa0e5Sjeremylt CeedVectorGetArrayRead(Uq, CEED_MEM_HOST, &uq); 5657c64913Sjeremylt sum = 0; 5757c64913Sjeremylt for (CeedInt i=0; i<Q; i++) { 5857c64913Sjeremylt sum += w[i] * uq[i]; 5957c64913Sjeremylt } 60*aedaa0e5Sjeremylt CeedVectorRestoreArrayRead(W, &w); 61*aedaa0e5Sjeremylt CeedVectorRestoreArrayRead(Uq, &uq); 62*aedaa0e5Sjeremylt 6357c64913Sjeremylt pint[0] = 0; 64*aedaa0e5Sjeremylt for (CeedInt i=0; i<(int)ALEN(p); i++) pint[i+1] = p[i] / (i+1); 6557c64913Sjeremylt error = sum - PolyEval(1, ALEN(pint), pint) + PolyEval(-1, ALEN(pint), pint); 66*aedaa0e5Sjeremylt if (error > 1.e-10) 6757c64913Sjeremylt printf("Error %e sum %g exact %g\n", error, sum, 6857c64913Sjeremylt PolyEval(1, ALEN(pint), pint) - PolyEval(-1, ALEN(pint), pint)); 6957c64913Sjeremylt 70*aedaa0e5Sjeremylt CeedVectorDestroy(&X); 71*aedaa0e5Sjeremylt CeedVectorDestroy(&Xq); 72*aedaa0e5Sjeremylt CeedVectorDestroy(&U); 73*aedaa0e5Sjeremylt CeedVectorDestroy(&Uq); 74*aedaa0e5Sjeremylt CeedVectorDestroy(&W); 7557c64913Sjeremylt CeedBasisDestroy(&bxl); 7657c64913Sjeremylt CeedBasisDestroy(&bxg); 7757c64913Sjeremylt CeedBasisDestroy(&bug); 7857c64913Sjeremylt CeedDestroy(&ceed); 7957c64913Sjeremylt return 0; 8057c64913Sjeremylt } 81