1 /// @file 2 /// Test polynomial derivative interpolation in 1D 3 /// \test Test polynomial derivative 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, X_q, U, U_q; 18 CeedBasis basis_x_lobatto, basis_u_lobatto, basis_x_gauss, basis_u_gauss; 19 CeedInt Q = 6; 20 const CeedScalar p[6] = {1, 2, 3, 4, 5, 6}; // 1 + 2x + 3x^2 + ... 21 const CeedScalar dp[5] = {2, 6, 12, 20, 30}; // 2 + 6x + 12x^2 + ... 22 const CeedScalar *xq, *uuq; 23 CeedScalar x[2], uq[Q]; 24 25 CeedInit(argv[1], &ceed); 26 27 CeedVectorCreate(ceed, 2, &X); 28 CeedVectorCreate(ceed, Q, &X_q); 29 CeedVectorSetValue(X_q, 0); 30 CeedVectorCreate(ceed, Q, &U); 31 CeedVectorSetValue(U, 0); 32 CeedVectorCreate(ceed, Q, &U_q); 33 34 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS_LOBATTO, 35 &basis_x_lobatto); 36 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, Q, Q, CEED_GAUSS_LOBATTO, 37 &basis_u_lobatto); 38 39 for (int i = 0; i < 2; i++) 40 x[i] = CeedIntPow(-1, i+1); 41 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, (CeedScalar *)&x); 42 43 CeedBasisApply(basis_x_lobatto, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, X_q); 44 45 CeedVectorGetArrayRead(X_q, CEED_MEM_HOST, &xq); 46 for (CeedInt i=0; i<Q; i++) 47 uq[i] = PolyEval(xq[i], ALEN(p), p); 48 CeedVectorRestoreArrayRead(X_q, &xq); 49 CeedVectorSetArray(U_q, CEED_MEM_HOST, CEED_USE_POINTER, (CeedScalar *)&uq); 50 51 // This operation is the identity because the quadrature is collocated 52 CeedBasisApply(basis_u_lobatto, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, U_q, U); 53 54 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS, &basis_x_gauss); 55 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, Q, Q, CEED_GAUSS, &basis_u_gauss); 56 57 CeedBasisApply(basis_x_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, X_q); 58 CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, U, U_q); 59 60 CeedVectorGetArrayRead(X_q, CEED_MEM_HOST, &xq); 61 CeedVectorGetArrayRead(U_q, CEED_MEM_HOST, &uuq); 62 for (CeedInt i=0; i<Q; i++) { 63 CeedScalar px = PolyEval(xq[i], ALEN(dp), dp); 64 if (fabs(uuq[i] - px) > 1000.*CEED_EPSILON) 65 // LCOV_EXCL_START 66 printf("%f != %f = p(%f)\n", uuq[i], px, xq[i]); 67 // LCOV_EXCL_STOP 68 } 69 CeedVectorRestoreArrayRead(X_q, &xq); 70 CeedVectorRestoreArrayRead(U_q, &uuq); 71 72 CeedVectorDestroy(&X); 73 CeedVectorDestroy(&X_q); 74 CeedVectorDestroy(&U); 75 CeedVectorDestroy(&U_q); 76 CeedBasisDestroy(&basis_x_lobatto); 77 CeedBasisDestroy(&basis_u_lobatto); 78 CeedBasisDestroy(&basis_x_gauss); 79 CeedBasisDestroy(&basis_u_gauss); 80 CeedDestroy(&ceed); 81 return 0; 82 } 83