1*d3f518cbSJeremy L Thompson /// @file 2*d3f518cbSJeremy L Thompson /// Test high order interpolation in multiple dimensions 3*d3f518cbSJeremy L Thompson /// \test Test high order interpolation in multiple dimensions 4*d3f518cbSJeremy L Thompson #include <ceed.h> 5*d3f518cbSJeremy L Thompson #include <math.h> 6*d3f518cbSJeremy L Thompson 7*d3f518cbSJeremy L Thompson static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) { 8*d3f518cbSJeremy L Thompson CeedScalar result = 1, center = 0.1; 9*d3f518cbSJeremy L Thompson for (CeedInt d=0; d<dim; d++) { 10*d3f518cbSJeremy L Thompson result *= tanh(x[d] - center); 11*d3f518cbSJeremy L Thompson center += 0.1; 12*d3f518cbSJeremy L Thompson } 13*d3f518cbSJeremy L Thompson return result; 14*d3f518cbSJeremy L Thompson } 15*d3f518cbSJeremy L Thompson 16*d3f518cbSJeremy L Thompson int main(int argc, char **argv) { 17*d3f518cbSJeremy L Thompson Ceed ceed; 18*d3f518cbSJeremy L Thompson 19*d3f518cbSJeremy L Thompson CeedInit(argv[1], &ceed); 20*d3f518cbSJeremy L Thompson for (CeedInt dim=1; dim<=3; dim++) { 21*d3f518cbSJeremy L Thompson CeedVector X, Xq, U, Uq; 22*d3f518cbSJeremy L Thompson CeedBasis bxl, bul, bxg, bug; 23*d3f518cbSJeremy L Thompson CeedInt Q = 11, Qdim = CeedIntPow(Q, dim), Xdim = CeedIntPow(2, dim); 24*d3f518cbSJeremy L Thompson CeedScalar x[Xdim*dim]; 25*d3f518cbSJeremy L Thompson const CeedScalar *xq, *u; 26*d3f518cbSJeremy L Thompson CeedScalar uq[Qdim]; 27*d3f518cbSJeremy L Thompson 28*d3f518cbSJeremy L Thompson for (CeedInt d=0; d<dim; d++) 29*d3f518cbSJeremy L Thompson for (CeedInt i=0; i<Xdim; i++) 30*d3f518cbSJeremy L Thompson x[d*Xdim + i] = (i % CeedIntPow(2, dim-d)) / CeedIntPow(2, dim-d-1) ? 1 : -1; 31*d3f518cbSJeremy L Thompson 32*d3f518cbSJeremy L Thompson CeedVectorCreate(ceed, Xdim*dim, &X); 33*d3f518cbSJeremy L Thompson CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 34*d3f518cbSJeremy L Thompson CeedVectorCreate(ceed, Qdim*dim, &Xq); 35*d3f518cbSJeremy L Thompson CeedVectorSetValue(Xq, 0); 36*d3f518cbSJeremy L Thompson CeedVectorCreate(ceed, Qdim, &U); 37*d3f518cbSJeremy L Thompson CeedVectorSetValue(U, 0); 38*d3f518cbSJeremy L Thompson CeedVectorCreate(ceed, Qdim, &Uq); 39*d3f518cbSJeremy L Thompson 40*d3f518cbSJeremy L Thompson CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, Q, CEED_GAUSS_LOBATTO, &bxl); 41*d3f518cbSJeremy L Thompson CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, Q, Q, CEED_GAUSS_LOBATTO, &bul); 42*d3f518cbSJeremy L Thompson 43*d3f518cbSJeremy L Thompson CeedBasisApply(bxl, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, Xq); 44*d3f518cbSJeremy L Thompson 45*d3f518cbSJeremy L Thompson CeedVectorGetArrayRead(Xq, CEED_MEM_HOST, &xq); 46*d3f518cbSJeremy L Thompson for (CeedInt i=0; i<Qdim; i++) { 47*d3f518cbSJeremy L Thompson CeedScalar xx[dim]; 48*d3f518cbSJeremy L Thompson for (CeedInt d=0; d<dim; d++) 49*d3f518cbSJeremy L Thompson xx[d] = xq[d*Qdim + i]; 50*d3f518cbSJeremy L Thompson uq[i] = Eval(dim, xx); 51*d3f518cbSJeremy L Thompson } 52*d3f518cbSJeremy L Thompson CeedVectorRestoreArrayRead(Xq, &xq); 53*d3f518cbSJeremy L Thompson CeedVectorSetArray(Uq, CEED_MEM_HOST, CEED_USE_POINTER, uq); 54*d3f518cbSJeremy L Thompson 55*d3f518cbSJeremy L Thompson // This operation is the identity because the quadrature is collocated 56*d3f518cbSJeremy L Thompson CeedBasisApply(bul, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, Uq, U); 57*d3f518cbSJeremy L Thompson 58*d3f518cbSJeremy L Thompson CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, Q, CEED_GAUSS, &bxg); 59*d3f518cbSJeremy L Thompson CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, Q, Q, CEED_GAUSS, &bug); 60*d3f518cbSJeremy L Thompson 61*d3f518cbSJeremy L Thompson CeedBasisApply(bxg, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, Xq); 62*d3f518cbSJeremy L Thompson CeedBasisApply(bug, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, U, Uq); 63*d3f518cbSJeremy L Thompson 64*d3f518cbSJeremy L Thompson CeedVectorGetArrayRead(Xq, CEED_MEM_HOST, &xq); 65*d3f518cbSJeremy L Thompson CeedVectorGetArrayRead(Uq, CEED_MEM_HOST, &u); 66*d3f518cbSJeremy L Thompson for (CeedInt i=0; i<Qdim; i++) { 67*d3f518cbSJeremy L Thompson CeedScalar xx[dim]; 68*d3f518cbSJeremy L Thompson for (CeedInt d=0; d<dim; d++) 69*d3f518cbSJeremy L Thompson xx[d] = xq[d*Qdim + i]; 70*d3f518cbSJeremy L Thompson CeedScalar fx = Eval(dim, xx); 71*d3f518cbSJeremy L Thompson if (fabs(u[i] - fx) > 1E-4) { 72*d3f518cbSJeremy L Thompson // LCOV_EXCL_START 73*d3f518cbSJeremy L Thompson printf("[%d] %f != %f=f(%f", dim, u[i], fx, xx[0]); 74*d3f518cbSJeremy L Thompson for (CeedInt d=1; d<dim; d++) printf(",%f", xx[d]); 75*d3f518cbSJeremy L Thompson puts(")"); 76*d3f518cbSJeremy L Thompson // LCOV_EXCL_STOP 77*d3f518cbSJeremy L Thompson } 78*d3f518cbSJeremy L Thompson } 79*d3f518cbSJeremy L Thompson CeedVectorRestoreArrayRead(Xq, &xq); 80*d3f518cbSJeremy L Thompson CeedVectorRestoreArrayRead(Uq, &u); 81*d3f518cbSJeremy L Thompson 82*d3f518cbSJeremy L Thompson CeedVectorDestroy(&X); 83*d3f518cbSJeremy L Thompson CeedVectorDestroy(&Xq); 84*d3f518cbSJeremy L Thompson CeedVectorDestroy(&U); 85*d3f518cbSJeremy L Thompson CeedVectorDestroy(&Uq); 86*d3f518cbSJeremy L Thompson CeedBasisDestroy(&bxl); 87*d3f518cbSJeremy L Thompson CeedBasisDestroy(&bul); 88*d3f518cbSJeremy L Thompson CeedBasisDestroy(&bxg); 89*d3f518cbSJeremy L Thompson CeedBasisDestroy(&bug); 90*d3f518cbSJeremy L Thompson } 91*d3f518cbSJeremy L Thompson CeedDestroy(&ceed); 92*d3f518cbSJeremy L Thompson return 0; 93*d3f518cbSJeremy L Thompson } 94