14411cf47Sjeremylt /// @file 24411cf47Sjeremylt /// Test interpolation in multiple dimensions 34411cf47Sjeremylt /// \test Test interpolation in multiple dimensions 457c64913Sjeremylt #include <ceed.h> 557c64913Sjeremylt #include <math.h> 657c64913Sjeremylt 757c64913Sjeremylt static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) { 857c64913Sjeremylt CeedScalar result = 1, center = 0.1; 957c64913Sjeremylt for (CeedInt d=0; d<dim; d++) { 1057c64913Sjeremylt result *= tanh(x[d] - center); 1157c64913Sjeremylt center += 0.1; 1257c64913Sjeremylt } 1357c64913Sjeremylt return result; 1457c64913Sjeremylt } 1557c64913Sjeremylt 1657c64913Sjeremylt int main(int argc, char **argv) { 1757c64913Sjeremylt Ceed ceed; 1857c64913Sjeremylt 1957c64913Sjeremylt CeedInit(argv[1], &ceed); 2057c64913Sjeremylt for (CeedInt dim=1; dim<=3; dim++) { 2157c64913Sjeremylt CeedBasis bxl, bul, bxg, bug; 22*b5cf12eeSjeremylt CeedInt Q = 10, Qdim = CeedIntPow(Q, dim), Xdim = CeedIntPow(2, dim); 2357c64913Sjeremylt CeedScalar x[Xdim*dim]; 2457c64913Sjeremylt CeedScalar xq[Qdim*dim], uq[Qdim], u[Qdim]; 2557c64913Sjeremylt 2657c64913Sjeremylt for (CeedInt d=0; d<dim; d++) { 2757c64913Sjeremylt for (CeedInt i=0; i<Xdim; i++) { 28*b5cf12eeSjeremylt x[d*Xdim + i] = (i % CeedIntPow(2, dim-d)) / CeedIntPow(2, dim-d-1) ? 1 : -1; 2957c64913Sjeremylt } 3057c64913Sjeremylt } 3157c64913Sjeremylt CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, Q, CEED_GAUSS_LOBATTO, &bxl); 3257c64913Sjeremylt CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, Q, Q, CEED_GAUSS_LOBATTO, &bul); 3357c64913Sjeremylt CeedBasisApply(bxl, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, xq); 3457c64913Sjeremylt for (CeedInt i=0; i<Qdim; i++) { 3557c64913Sjeremylt CeedScalar xx[dim]; 3657c64913Sjeremylt for (CeedInt d=0; d<dim; d++) xx[d] = xq[d*Qdim + i]; 3757c64913Sjeremylt uq[i] = Eval(dim, xx); 3857c64913Sjeremylt } 3957c64913Sjeremylt 4057c64913Sjeremylt // This operation is the identity because the quadrature is collocated 4157c64913Sjeremylt CeedBasisApply(bul, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, uq, u); 4257c64913Sjeremylt 4357c64913Sjeremylt CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, Q, CEED_GAUSS, &bxg); 4457c64913Sjeremylt CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, Q, Q, CEED_GAUSS, &bug); 4557c64913Sjeremylt CeedBasisApply(bxg, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, xq); 4657c64913Sjeremylt CeedBasisApply(bug, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, uq); 4757c64913Sjeremylt for (CeedInt i=0; i<Qdim; i++) { 4857c64913Sjeremylt CeedScalar xx[dim]; 4957c64913Sjeremylt for (CeedInt d=0; d<dim; d++) xx[d] = xq[d*Qdim + i]; 5057c64913Sjeremylt CeedScalar fx = Eval(dim, xx); 5157c64913Sjeremylt if (!(fabs(uq[i] - fx) < 1e-4)) { 5257c64913Sjeremylt printf("[%d] %f != %f=f(%f", dim, uq[i], fx, xx[0]); 5357c64913Sjeremylt for (CeedInt d=1; d<dim; d++) printf(",%f", xx[d]); 5457c64913Sjeremylt puts(")"); 5557c64913Sjeremylt } 5657c64913Sjeremylt } 5757c64913Sjeremylt 5857c64913Sjeremylt CeedBasisDestroy(&bxl); 5957c64913Sjeremylt CeedBasisDestroy(&bul); 6057c64913Sjeremylt CeedBasisDestroy(&bxg); 6157c64913Sjeremylt CeedBasisDestroy(&bug); 6257c64913Sjeremylt } 6357c64913Sjeremylt CeedDestroy(&ceed); 6457c64913Sjeremylt return 0; 6557c64913Sjeremylt } 66