1 /// @file 2 /// Test interpolation in multiple dimensions 3 /// \test Test 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 = 10, 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++) 30 x[d*X_dim + i] = (i % CeedIntPow(2, dim-d)) / CeedIntPow(2, dim-d-1) ? 1 : -1; 31 32 CeedVectorCreate(ceed, X_dim*dim, &X); 33 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 34 CeedVectorCreate(ceed, Q_dim*dim, &X_q); 35 CeedVectorSetValue(X_q, 0); 36 CeedVectorCreate(ceed, Q_dim, &U); 37 CeedVectorSetValue(U, 0); 38 CeedVectorCreate(ceed, Q_dim, &U_q); 39 40 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, Q, CEED_GAUSS_LOBATTO, 41 &basis_x_lobatto); 42 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, Q, Q, CEED_GAUSS_LOBATTO, 43 &basis_u_lobatto); 44 45 CeedBasisApply(basis_x_lobatto, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, X_q); 46 47 CeedVectorGetArrayRead(X_q, CEED_MEM_HOST, &xq); 48 for (CeedInt i=0; i<Q_dim; i++) { 49 CeedScalar xx[dim]; 50 for (CeedInt d=0; d<dim; d++) 51 xx[d] = xq[d*Q_dim + i]; 52 uq[i] = Eval(dim, xx); 53 } 54 CeedVectorRestoreArrayRead(X_q, &xq); 55 CeedVectorSetArray(U_q, CEED_MEM_HOST, CEED_USE_POINTER, uq); 56 57 // This operation is the identity because the quadrature is collocated 58 CeedBasisApply(basis_u_lobatto, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, U_q, U); 59 60 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, Q, CEED_GAUSS, 61 &basis_x_gauss); 62 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, Q, Q, CEED_GAUSS, &basis_u_gauss); 63 64 CeedBasisApply(basis_x_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, X_q); 65 CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, U, U_q); 66 67 CeedVectorGetArrayRead(X_q, CEED_MEM_HOST, &xq); 68 CeedVectorGetArrayRead(U_q, CEED_MEM_HOST, &u); 69 for (CeedInt i=0; i<Q_dim; i++) { 70 CeedScalar xx[dim]; 71 for (CeedInt d=0; d<dim; d++) 72 xx[d] = xq[d*Q_dim + i]; 73 CeedScalar fx = Eval(dim, xx); 74 if (fabs(u[i] - fx) > 1E-4) { 75 // LCOV_EXCL_START 76 printf("[%d] %f != %f=f(%f", dim, u[i], fx, xx[0]); 77 for (CeedInt d=1; d<dim; d++) printf(",%f", xx[d]); 78 puts(")"); 79 // LCOV_EXCL_STOP 80 } 81 } 82 CeedVectorRestoreArrayRead(X_q, &xq); 83 CeedVectorRestoreArrayRead(U_q, &u); 84 85 CeedVectorDestroy(&X); 86 CeedVectorDestroy(&X_q); 87 CeedVectorDestroy(&U); 88 CeedVectorDestroy(&U_q); 89 CeedBasisDestroy(&basis_x_lobatto); 90 CeedBasisDestroy(&basis_u_lobatto); 91 CeedBasisDestroy(&basis_x_gauss); 92 CeedBasisDestroy(&basis_u_gauss); 93 } 94 CeedDestroy(&ceed); 95 return 0; 96 } 97