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; 1857c64913Sjeremylt CeedBasis bxl, bul, 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, *uuq; 22*aedaa0e5Sjeremylt CeedScalar x[2], uq[Q]; 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 CeedVectorSetValue(U, 0); 31*aedaa0e5Sjeremylt CeedVectorCreate(ceed, Q, &Uq); 32*aedaa0e5Sjeremylt 3357c64913Sjeremylt CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS_LOBATTO, &bxl); 3457c64913Sjeremylt CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, Q, Q, CEED_GAUSS_LOBATTO, &bul); 35*aedaa0e5Sjeremylt 36*aedaa0e5Sjeremylt for (int i = 0; i < 2; i++) x[i] = CeedIntPow(-1, i+1); 37*aedaa0e5Sjeremylt CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, (CeedScalar *)&x); 38*aedaa0e5Sjeremylt 39*aedaa0e5Sjeremylt CeedBasisApply(bxl, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, Xq); 40*aedaa0e5Sjeremylt 41*aedaa0e5Sjeremylt CeedVectorGetArrayRead(Xq, CEED_MEM_HOST, &xq); 4257c64913Sjeremylt for (CeedInt i=0; i<Q; i++) uq[i] = PolyEval(xq[i], ALEN(p), p); 43*aedaa0e5Sjeremylt CeedVectorRestoreArrayRead(Xq, &xq); 44*aedaa0e5Sjeremylt CeedVectorSetArray(Uq, CEED_MEM_HOST, CEED_USE_POINTER, (CeedScalar *)&uq); 4557c64913Sjeremylt 4657c64913Sjeremylt // This operation is the identity because the quadrature is collocated 47*aedaa0e5Sjeremylt CeedBasisApply(bul, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, Uq, U); 4857c64913Sjeremylt 4957c64913Sjeremylt CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS, &bxg); 5057c64913Sjeremylt CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, Q, Q, CEED_GAUSS, &bug); 51*aedaa0e5Sjeremylt 52*aedaa0e5Sjeremylt CeedBasisApply(bxg, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, Xq); 53*aedaa0e5Sjeremylt CeedBasisApply(bug, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, U, Uq); 54*aedaa0e5Sjeremylt 55*aedaa0e5Sjeremylt CeedVectorGetArrayRead(Xq, CEED_MEM_HOST, &xq); 56*aedaa0e5Sjeremylt CeedVectorGetArrayRead(Uq, CEED_MEM_HOST, &uuq); 5757c64913Sjeremylt for (CeedInt i=0; i<Q; i++) { 5857c64913Sjeremylt CeedScalar px = PolyEval(xq[i], ALEN(p), p); 59*aedaa0e5Sjeremylt if ((fabs(uuq[i] - px) > 1e-14)) { 60*aedaa0e5Sjeremylt printf("%f != %f=p(%f)\n", uuq[i], px, xq[i]); 6157c64913Sjeremylt } 6257c64913Sjeremylt } 63*aedaa0e5Sjeremylt CeedVectorRestoreArrayRead(Xq, &xq); 64*aedaa0e5Sjeremylt CeedVectorRestoreArrayRead(Uq, &uuq); 6557c64913Sjeremylt 66*aedaa0e5Sjeremylt CeedVectorDestroy(&X); 67*aedaa0e5Sjeremylt CeedVectorDestroy(&Xq); 68*aedaa0e5Sjeremylt CeedVectorDestroy(&U); 69*aedaa0e5Sjeremylt CeedVectorDestroy(&Uq); 7057c64913Sjeremylt CeedBasisDestroy(&bxl); 7157c64913Sjeremylt CeedBasisDestroy(&bul); 7257c64913Sjeremylt CeedBasisDestroy(&bxg); 7357c64913Sjeremylt CeedBasisDestroy(&bug); 7457c64913Sjeremylt CeedDestroy(&ceed); 7557c64913Sjeremylt return 0; 7657c64913Sjeremylt } 77