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