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; 17*d1d35e2fSjeremylt CeedVector X, X_q, U, U_q; 18*d1d35e2fSjeremylt CeedBasis basis_x_lobatto, basis_u_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, *uuq; 2252bfb9bbSJeremy L Thompson CeedScalar x[2], uq[Q]; 23a8de75f0Sjeremylt 24a8de75f0Sjeremylt CeedInit(argv[1], &ceed); 25aedaa0e5Sjeremylt 2652bfb9bbSJeremy L Thompson CeedVectorCreate(ceed, 2, &X); 27*d1d35e2fSjeremylt CeedVectorCreate(ceed, Q, &X_q); 28*d1d35e2fSjeremylt CeedVectorSetValue(X_q, 0); 2952bfb9bbSJeremy L Thompson CeedVectorCreate(ceed, Q, &U); 3052bfb9bbSJeremy L Thompson CeedVectorSetValue(U, 0); 31*d1d35e2fSjeremylt CeedVectorCreate(ceed, Q, &U_q); 32a8de75f0Sjeremylt 33*d1d35e2fSjeremylt CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS_LOBATTO, 34*d1d35e2fSjeremylt &basis_x_lobatto); 35*d1d35e2fSjeremylt CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, Q, Q, CEED_GAUSS_LOBATTO, 36*d1d35e2fSjeremylt &basis_u_lobatto); 37aedaa0e5Sjeremylt 3852bfb9bbSJeremy L Thompson for (int i = 0; i < 2; i++) 3952bfb9bbSJeremy L Thompson x[i] = CeedIntPow(-1, i+1); 4052bfb9bbSJeremy L Thompson CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, (CeedScalar *)&x); 41aedaa0e5Sjeremylt 42*d1d35e2fSjeremylt CeedBasisApply(basis_x_lobatto, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, X_q); 43a8de75f0Sjeremylt 44*d1d35e2fSjeremylt CeedVectorGetArrayRead(X_q, CEED_MEM_HOST, &xq); 4552bfb9bbSJeremy L Thompson for (CeedInt i=0; i<Q; i++) 4652bfb9bbSJeremy L Thompson uq[i] = PolyEval(xq[i], ALEN(p), p); 47*d1d35e2fSjeremylt CeedVectorRestoreArrayRead(X_q, &xq); 48*d1d35e2fSjeremylt CeedVectorSetArray(U_q, CEED_MEM_HOST, CEED_USE_POINTER, (CeedScalar *)&uq); 4952bfb9bbSJeremy L Thompson 5052bfb9bbSJeremy L Thompson // This operation is the identity because the quadrature is collocated 51*d1d35e2fSjeremylt CeedBasisApply(basis_u_lobatto, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, U_q, U); 5252bfb9bbSJeremy L Thompson 53*d1d35e2fSjeremylt CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS, &basis_x_gauss); 54*d1d35e2fSjeremylt CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, Q, Q, CEED_GAUSS, &basis_u_gauss); 5552bfb9bbSJeremy L Thompson 56*d1d35e2fSjeremylt CeedBasisApply(basis_x_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, X_q); 57*d1d35e2fSjeremylt CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, U, U_q); 5852bfb9bbSJeremy L Thompson 59*d1d35e2fSjeremylt CeedVectorGetArrayRead(X_q, CEED_MEM_HOST, &xq); 60*d1d35e2fSjeremylt CeedVectorGetArrayRead(U_q, CEED_MEM_HOST, &uuq); 6152bfb9bbSJeremy L Thompson for (CeedInt i=0; i<Q; i++) { 6252bfb9bbSJeremy L Thompson CeedScalar px = PolyEval(xq[i], ALEN(p), p); 63166c82afSvaleriabarra if (fabs(uuq[i] - px) > 1E-14) 64a2546046Sjeremylt // LCOV_EXCL_START 6552bfb9bbSJeremy L Thompson printf("%f != %f=p(%f)\n", uuq[i], px, xq[i]); 66de996c55Sjeremylt // LCOV_EXCL_STOP 67a8de75f0Sjeremylt } 68*d1d35e2fSjeremylt CeedVectorRestoreArrayRead(X_q, &xq); 69*d1d35e2fSjeremylt CeedVectorRestoreArrayRead(U_q, &uuq); 70a8de75f0Sjeremylt 7152bfb9bbSJeremy L Thompson CeedVectorDestroy(&X); 72*d1d35e2fSjeremylt CeedVectorDestroy(&X_q); 7352bfb9bbSJeremy L Thompson CeedVectorDestroy(&U); 74*d1d35e2fSjeremylt CeedVectorDestroy(&U_q); 75*d1d35e2fSjeremylt CeedBasisDestroy(&basis_x_lobatto); 76*d1d35e2fSjeremylt CeedBasisDestroy(&basis_u_lobatto); 77*d1d35e2fSjeremylt CeedBasisDestroy(&basis_x_gauss); 78*d1d35e2fSjeremylt CeedBasisDestroy(&basis_u_gauss); 79a8de75f0Sjeremylt CeedDestroy(&ceed); 80a8de75f0Sjeremylt return 0; 81a8de75f0Sjeremylt } 82