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, Xq, U, Uq; 22 CeedBasis bxl, bul, bxg, bug; 23 CeedInt Q = 10, Qdim = CeedIntPow(Q, dim), Xdim = CeedIntPow(2, dim); 24 CeedScalar x[Xdim*dim]; 25 const CeedScalar *xq, *u; 26 CeedScalar uq[Qdim]; 27 28 for (CeedInt d=0; d<dim; d++) 29 for (CeedInt i=0; i<Xdim; i++) 30 x[d*Xdim + i] = (i % CeedIntPow(2, dim-d)) / CeedIntPow(2, dim-d-1) ? 1 : -1; 31 32 CeedVectorCreate(ceed, Xdim*dim, &X); 33 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 34 CeedVectorCreate(ceed, Qdim*dim, &Xq); 35 CeedVectorSetValue(Xq, 0); 36 CeedVectorCreate(ceed, Qdim, &U); 37 CeedVectorSetValue(U, 0); 38 CeedVectorCreate(ceed, Qdim, &Uq); 39 40 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, Q, CEED_GAUSS_LOBATTO, &bxl); 41 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, Q, Q, CEED_GAUSS_LOBATTO, &bul); 42 43 CeedBasisApply(bxl, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, Xq); 44 45 CeedVectorGetArrayRead(Xq, CEED_MEM_HOST, &xq); 46 for (CeedInt i=0; i<Qdim; i++) { 47 CeedScalar xx[dim]; 48 for (CeedInt d=0; d<dim; d++) 49 xx[d] = xq[d*Qdim + i]; 50 uq[i] = Eval(dim, xx); 51 } 52 CeedVectorRestoreArrayRead(Xq, &xq); 53 CeedVectorSetArray(Uq, CEED_MEM_HOST, CEED_USE_POINTER, uq); 54 55 // This operation is the identity because the quadrature is collocated 56 CeedBasisApply(bul, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, Uq, U); 57 58 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, Q, CEED_GAUSS, &bxg); 59 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, Q, Q, CEED_GAUSS, &bug); 60 61 CeedBasisApply(bxg, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, Xq); 62 CeedBasisApply(bug, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, U, Uq); 63 64 CeedVectorGetArrayRead(Xq, CEED_MEM_HOST, &xq); 65 CeedVectorGetArrayRead(Uq, CEED_MEM_HOST, &u); 66 for (CeedInt i=0; i<Qdim; i++) { 67 CeedScalar xx[dim]; 68 for (CeedInt d=0; d<dim; d++) 69 xx[d] = xq[d*Qdim + i]; 70 CeedScalar fx = Eval(dim, xx); 71 if (fabs(u[i] - fx) > 1e-4) { 72 // LCOV_EXCL_START 73 printf("[%d] %f != %f=f(%f", dim, u[i], fx, xx[0]); 74 for (CeedInt d=1; d<dim; d++) printf(",%f", xx[d]); 75 puts(")"); 76 // LCOV_EXCL_STOP 77 } 78 } 79 CeedVectorRestoreArrayRead(Xq, &xq); 80 CeedVectorRestoreArrayRead(Uq, &u); 81 82 CeedVectorDestroy(&X); 83 CeedVectorDestroy(&Xq); 84 CeedVectorDestroy(&U); 85 CeedVectorDestroy(&Uq); 86 CeedBasisDestroy(&bxl); 87 CeedBasisDestroy(&bul); 88 CeedBasisDestroy(&bxg); 89 CeedBasisDestroy(&bug); 90 } 91 CeedDestroy(&ceed); 92 return 0; 93 } 94