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