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 CeedBasis bxl, bul, bxg, bug; 22 CeedInt Q = 10, Qdim = CeedIntPow(Q, dim), Xdim = CeedIntPow(2, dim); 23 CeedScalar x[Xdim*dim]; 24 CeedScalar xq[Qdim*dim], uq[Qdim], u[Qdim]; 25 26 for (CeedInt d=0; d<dim; d++) { 27 for (CeedInt i=0; i<Xdim; i++) { 28 x[d*Xdim + i] = (i % CeedIntPow(2, dim-d)) / CeedIntPow(2, dim-d-1) ? 1 : -1; 29 } 30 } 31 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, Q, CEED_GAUSS_LOBATTO, &bxl); 32 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, Q, Q, CEED_GAUSS_LOBATTO, &bul); 33 CeedBasisApply(bxl, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, xq); 34 for (CeedInt i=0; i<Qdim; i++) { 35 CeedScalar xx[dim]; 36 for (CeedInt d=0; d<dim; d++) xx[d] = xq[d*Qdim + i]; 37 uq[i] = Eval(dim, xx); 38 } 39 40 // This operation is the identity because the quadrature is collocated 41 CeedBasisApply(bul, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, uq, u); 42 43 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, Q, CEED_GAUSS, &bxg); 44 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, Q, Q, CEED_GAUSS, &bug); 45 CeedBasisApply(bxg, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, xq); 46 CeedBasisApply(bug, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, uq); 47 for (CeedInt i=0; i<Qdim; i++) { 48 CeedScalar xx[dim]; 49 for (CeedInt d=0; d<dim; d++) xx[d] = xq[d*Qdim + i]; 50 CeedScalar fx = Eval(dim, xx); 51 if (!(fabs(uq[i] - fx) < 1e-4)) { 52 printf("[%d] %f != %f=f(%f", dim, uq[i], fx, xx[0]); 53 for (CeedInt d=1; d<dim; d++) printf(",%f", xx[d]); 54 puts(")"); 55 } 56 } 57 58 CeedBasisDestroy(&bxl); 59 CeedBasisDestroy(&bul); 60 CeedBasisDestroy(&bxg); 61 CeedBasisDestroy(&bug); 62 } 63 CeedDestroy(&ceed); 64 return 0; 65 } 66