1 /// @file 2 /// Test high order interpolation in multiple dimensions 3 /// \test Test high order interpolation in multiple dimensions 4 #include <ceed.h> 5 #include <math.h> 6 7 static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) { 8 CeedScalar result = 1, center = 0.1; 9 for (CeedInt d = 0; d < dim; d++) { 10 result *= tanh(x[d] - center); 11 center += 0.1; 12 } 13 return result; 14 } 15 16 int main(int argc, char **argv) { 17 Ceed ceed; 18 19 CeedInit(argv[1], &ceed); 20 for (CeedInt dim = 1; dim <= 3; dim++) { 21 CeedVector X, X_q, U, U_q; 22 CeedBasis basis_x_lobatto, basis_u_lobatto, basis_x_gauss, basis_u_gauss; 23 CeedInt Q = 11, Q_dim = CeedIntPow(Q, dim), X_dim = CeedIntPow(2, dim); 24 CeedScalar x[X_dim * dim]; 25 const CeedScalar *xq, *u; 26 CeedScalar uq[Q_dim]; 27 28 for (CeedInt d = 0; d < dim; d++) 29 for (CeedInt i = 0; i < X_dim; i++) x[d * X_dim + i] = (i % CeedIntPow(2, dim - d)) / CeedIntPow(2, dim - d - 1) ? 1 : -1; 30 31 CeedVectorCreate(ceed, X_dim * dim, &X); 32 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 33 CeedVectorCreate(ceed, Q_dim * dim, &X_q); 34 CeedVectorSetValue(X_q, 0); 35 CeedVectorCreate(ceed, Q_dim, &U); 36 CeedVectorSetValue(U, 0); 37 CeedVectorCreate(ceed, Q_dim, &U_q); 38 39 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, Q, CEED_GAUSS_LOBATTO, &basis_x_lobatto); 40 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, Q, Q, CEED_GAUSS_LOBATTO, &basis_u_lobatto); 41 42 CeedBasisApply(basis_x_lobatto, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, X_q); 43 44 CeedVectorGetArrayRead(X_q, CEED_MEM_HOST, &xq); 45 for (CeedInt i = 0; i < Q_dim; i++) { 46 CeedScalar xx[dim]; 47 for (CeedInt d = 0; d < dim; d++) xx[d] = xq[d * Q_dim + i]; 48 uq[i] = Eval(dim, xx); 49 } 50 CeedVectorRestoreArrayRead(X_q, &xq); 51 CeedVectorSetArray(U_q, CEED_MEM_HOST, CEED_USE_POINTER, uq); 52 53 // This operation is the identity because the quadrature is collocated 54 CeedBasisApply(basis_u_lobatto, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, U_q, U); 55 56 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, Q, CEED_GAUSS, &basis_x_gauss); 57 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, Q, Q, CEED_GAUSS, &basis_u_gauss); 58 59 CeedBasisApply(basis_x_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, X_q); 60 CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, U, U_q); 61 62 CeedVectorGetArrayRead(X_q, CEED_MEM_HOST, &xq); 63 CeedVectorGetArrayRead(U_q, CEED_MEM_HOST, &u); 64 for (CeedInt i = 0; i < Q_dim; i++) { 65 CeedScalar xx[dim]; 66 for (CeedInt d = 0; d < dim; d++) xx[d] = xq[d * Q_dim + i]; 67 CeedScalar fx = Eval(dim, xx); 68 if (fabs(u[i] - fx) > 1E-4) { 69 // LCOV_EXCL_START 70 printf("[%" CeedInt_FMT "] %f != %f = f(%f", dim, u[i], fx, xx[0]); 71 for (CeedInt d = 1; d < dim; d++) printf(",%f", xx[d]); 72 puts(")"); 73 // LCOV_EXCL_STOP 74 } 75 } 76 CeedVectorRestoreArrayRead(X_q, &xq); 77 CeedVectorRestoreArrayRead(U_q, &u); 78 79 CeedVectorDestroy(&X); 80 CeedVectorDestroy(&X_q); 81 CeedVectorDestroy(&U); 82 CeedVectorDestroy(&U_q); 83 CeedBasisDestroy(&basis_x_lobatto); 84 CeedBasisDestroy(&basis_u_lobatto); 85 CeedBasisDestroy(&basis_x_gauss); 86 CeedBasisDestroy(&basis_u_gauss); 87 } 88 CeedDestroy(&ceed); 89 return 0; 90 } 91