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