1a8de75f0Sjeremylt /// @file 252bfb9bbSJeremy L Thompson /// Test polynomial interpolation in 1D 352bfb9bbSJeremy L Thompson /// \test Test polynomial interpolation in 1D 4a8de75f0Sjeremylt #include <ceed.h> 5a8de75f0Sjeremylt #include <math.h> 6a8de75f0Sjeremylt 752bfb9bbSJeremy L Thompson #define ALEN(a) (sizeof(a) / sizeof((a)[0])) 852bfb9bbSJeremy L Thompson 952bfb9bbSJeremy L Thompson static CeedScalar PolyEval(CeedScalar x, CeedInt n, const CeedScalar *p) { 1052bfb9bbSJeremy L Thompson CeedScalar y = p[n - 1]; 1152bfb9bbSJeremy L Thompson for (CeedInt i = n - 2; i >= 0; i--) y = y * x + p[i]; 1252bfb9bbSJeremy L Thompson return y; 13a8de75f0Sjeremylt } 14a8de75f0Sjeremylt 15a8de75f0Sjeremylt int main(int argc, char **argv) { 16a8de75f0Sjeremylt Ceed ceed; 17d1d35e2fSjeremylt CeedVector X, X_q, U, U_q, W; 18d1d35e2fSjeremylt CeedBasis basis_x_lobatto, basis_x_gauss, basis_u_gauss; 1952bfb9bbSJeremy L Thompson CeedInt Q = 6; 2052bfb9bbSJeremy L Thompson const CeedScalar p[6] = {1, 2, 3, 4, 5, 6}; // 1 + 2x + 3x^2 + ... 2152bfb9bbSJeremy L Thompson const CeedScalar *xq, *uq, *w; 2252bfb9bbSJeremy L Thompson CeedScalar u[Q], x[2], sum, error, pint[ALEN(p) + 1]; 23a8de75f0Sjeremylt 24a8de75f0Sjeremylt CeedInit(argv[1], &ceed); 25aedaa0e5Sjeremylt 2652bfb9bbSJeremy L Thompson CeedVectorCreate(ceed, 2, &X); 27d1d35e2fSjeremylt CeedVectorCreate(ceed, Q, &X_q); 28d1d35e2fSjeremylt CeedVectorSetValue(X_q, 0); 2952bfb9bbSJeremy L Thompson CeedVectorCreate(ceed, Q, &U); 30d1d35e2fSjeremylt CeedVectorCreate(ceed, Q, &U_q); 31d1d35e2fSjeremylt CeedVectorSetValue(U_q, 0); 3252bfb9bbSJeremy L Thompson CeedVectorCreate(ceed, Q, &W); 3352bfb9bbSJeremy L Thompson CeedVectorSetValue(W, 0); 34a8de75f0Sjeremylt 35*2b730f8bSJeremy L Thompson CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS_LOBATTO, &basis_x_lobatto); 36aedaa0e5Sjeremylt 37*2b730f8bSJeremy L Thompson for (int i = 0; i < 2; i++) x[i] = CeedIntPow(-1, i + 1); 3852bfb9bbSJeremy L Thompson CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 39aedaa0e5Sjeremylt 40d1d35e2fSjeremylt CeedBasisApply(basis_x_lobatto, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, X_q); 41a8de75f0Sjeremylt 42d1d35e2fSjeremylt CeedVectorGetArrayRead(X_q, CEED_MEM_HOST, &xq); 43*2b730f8bSJeremy L Thompson for (CeedInt i = 0; i < Q; i++) u[i] = PolyEval(xq[i], ALEN(p), p); 44d1d35e2fSjeremylt CeedVectorRestoreArrayRead(X_q, &xq); 4552bfb9bbSJeremy L Thompson CeedVectorSetArray(U, CEED_MEM_HOST, CEED_USE_POINTER, u); 4652bfb9bbSJeremy L Thompson 47d1d35e2fSjeremylt CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS, &basis_x_gauss); 48d1d35e2fSjeremylt CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, Q, Q, CEED_GAUSS, &basis_u_gauss); 4952bfb9bbSJeremy L Thompson 50d1d35e2fSjeremylt CeedBasisApply(basis_x_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, X_q); 51d1d35e2fSjeremylt CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, U, U_q); 52*2b730f8bSJeremy L Thompson CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, W); 5352bfb9bbSJeremy L Thompson 5452bfb9bbSJeremy L Thompson CeedVectorGetArrayRead(W, CEED_MEM_HOST, &w); 55d1d35e2fSjeremylt CeedVectorGetArrayRead(U_q, CEED_MEM_HOST, &uq); 56a8de75f0Sjeremylt sum = 0; 57*2b730f8bSJeremy L Thompson for (CeedInt i = 0; i < Q; i++) sum += w[i] * uq[i]; 5852bfb9bbSJeremy L Thompson CeedVectorRestoreArrayRead(W, &w); 59d1d35e2fSjeremylt CeedVectorRestoreArrayRead(U_q, &uq); 60a8de75f0Sjeremylt 6152bfb9bbSJeremy L Thompson pint[0] = 0; 62*2b730f8bSJeremy L Thompson for (CeedInt i = 0; i < (int)ALEN(p); i++) pint[i + 1] = p[i] / (i + 1); 6352bfb9bbSJeremy L Thompson error = sum - PolyEval(1, ALEN(pint), pint) + PolyEval(-1, ALEN(pint), pint); 64*2b730f8bSJeremy L Thompson if (error > 100. * CEED_EPSILON) { 6552bfb9bbSJeremy L Thompson // LCOV_EXCL_START 66*2b730f8bSJeremy L Thompson printf("Error %e sum %g exact %g\n", error, sum, PolyEval(1, ALEN(pint), pint) - PolyEval(-1, ALEN(pint), pint)); 6752bfb9bbSJeremy L Thompson // LCOV_EXCL_STOP 68*2b730f8bSJeremy L Thompson } 6952bfb9bbSJeremy L Thompson 7052bfb9bbSJeremy L Thompson CeedVectorDestroy(&X); 71d1d35e2fSjeremylt CeedVectorDestroy(&X_q); 7252bfb9bbSJeremy L Thompson CeedVectorDestroy(&U); 73d1d35e2fSjeremylt CeedVectorDestroy(&U_q); 7452bfb9bbSJeremy L Thompson CeedVectorDestroy(&W); 75d1d35e2fSjeremylt CeedBasisDestroy(&basis_x_lobatto); 76d1d35e2fSjeremylt CeedBasisDestroy(&basis_x_gauss); 77d1d35e2fSjeremylt CeedBasisDestroy(&basis_u_gauss); 78a8de75f0Sjeremylt CeedDestroy(&ceed); 79a8de75f0Sjeremylt return 0; 80a8de75f0Sjeremylt } 81