1 // Test interpolation in multiple dimensions 2 #include <ceed.h> 3 #include <math.h> 4 5 static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) { 6 CeedScalar result = 1, center = 0.1; 7 for (CeedInt d=0; d<dim; d++) { 8 result *= tanh(x[d] - center); 9 center += 0.1; 10 } 11 return result; 12 } 13 14 int main(int argc, char **argv) { 15 Ceed ceed; 16 17 CeedInit(argv[1], &ceed); 18 for (CeedInt dim=1; dim<=3; dim++) { 19 CeedBasis bxl, bul, bxg, bug; 20 CeedInt Q = 10, Qdim = CeedPowInt(Q, dim), Xdim = CeedPowInt(2, dim); 21 CeedScalar x[Xdim*dim]; 22 CeedScalar xq[Qdim*dim], uq[Qdim], u[Qdim]; 23 24 for (CeedInt d=0; d<dim; d++) { 25 for (CeedInt i=0; i<Xdim; i++) { 26 x[d*Xdim + i] = (i % CeedPowInt(2, dim-d)) / CeedPowInt(2, dim-d-1) ? 1 : -1; 27 } 28 } 29 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, Q, CEED_GAUSS_LOBATTO, &bxl); 30 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, Q, Q, CEED_GAUSS_LOBATTO, &bul); 31 CeedBasisApply(bxl, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, xq); 32 for (CeedInt i=0; i<Qdim; i++) { 33 CeedScalar xx[dim]; 34 for (CeedInt d=0; d<dim; d++) xx[d] = xq[d*Qdim + i]; 35 uq[i] = Eval(dim, xx); 36 } 37 38 // This operation is the identity because the quadrature is collocated 39 CeedBasisApply(bul, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, uq, u); 40 41 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, Q, CEED_GAUSS, &bxg); 42 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, Q, Q, CEED_GAUSS, &bug); 43 CeedBasisApply(bxg, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, xq); 44 CeedBasisApply(bug, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, uq); 45 for (CeedInt i=0; i<Qdim; i++) { 46 CeedScalar xx[dim]; 47 for (CeedInt d=0; d<dim; d++) xx[d] = xq[d*Qdim + i]; 48 CeedScalar fx = Eval(dim, xx); 49 if (!(fabs(uq[i] - fx) < 1e-4)) { 50 printf("[%d] %f != %f=f(%f", dim, uq[i], fx, xx[0]); 51 for (CeedInt d=1; d<dim; d++) printf(",%f", xx[d]); 52 puts(")"); 53 } 54 } 55 56 CeedBasisDestroy(&bxl); 57 CeedBasisDestroy(&bul); 58 CeedBasisDestroy(&bxg); 59 CeedBasisDestroy(&bug); 60 } 61 CeedDestroy(&ceed); 62 return 0; 63 } 64