152bfb9bbSJeremy L Thompson /// @file 252bfb9bbSJeremy L Thompson /// Test collocated grad in multiple dimensions 352bfb9bbSJeremy L Thompson /// \test Test collocated grad in multiple dimensions 452bfb9bbSJeremy L Thompson #include <ceed.h> 552bfb9bbSJeremy L Thompson #include <math.h> 652bfb9bbSJeremy L Thompson 752bfb9bbSJeremy L Thompson static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) { 852bfb9bbSJeremy L Thompson CeedScalar result = tanh(x[0] + 0.1); 952bfb9bbSJeremy L Thompson if (dim > 1) result += atan(x[1] + 0.2); 1052bfb9bbSJeremy L Thompson if (dim > 2) result += exp(-(x[2] + 0.3)*(x[2] + 0.3)); 1152bfb9bbSJeremy L Thompson return result; 1252bfb9bbSJeremy L Thompson } 1352bfb9bbSJeremy L Thompson 1452bfb9bbSJeremy L Thompson int main(int argc, char **argv) { 1552bfb9bbSJeremy L Thompson Ceed ceed; 1652bfb9bbSJeremy L Thompson 1752bfb9bbSJeremy L Thompson CeedInit(argv[1], &ceed); 1852bfb9bbSJeremy L Thompson 1952bfb9bbSJeremy L Thompson for (CeedInt dim=1; dim<=3; dim++) { 2052bfb9bbSJeremy L Thompson CeedVector X, Xq, U, Uq, Ones, Gtposeones; 2152bfb9bbSJeremy L Thompson CeedBasis bxl, bug; 2252bfb9bbSJeremy L Thompson CeedInt P = 8, Q = 8, Pdim = CeedIntPow(P, dim), Qdim = CeedIntPow(Q, dim), 2352bfb9bbSJeremy L Thompson Xdim = CeedIntPow(2, dim); 2452bfb9bbSJeremy L Thompson CeedScalar x[Xdim*dim], u[Pdim]; 2552bfb9bbSJeremy L Thompson const CeedScalar *xq, *uq, *gtposeones; 2652bfb9bbSJeremy L Thompson CeedScalar sum1 = 0, sum2 = 0; 2752bfb9bbSJeremy L Thompson 2852bfb9bbSJeremy L Thompson for (CeedInt d=0; d<dim; d++) 2952bfb9bbSJeremy L Thompson for (CeedInt i=0; i<Xdim; i++) 3052bfb9bbSJeremy L Thompson x[d*Xdim + i] = (i % CeedIntPow(2, dim-d)) / 3152bfb9bbSJeremy L Thompson CeedIntPow(2, dim-d-1) ? 1 : -1; 3252bfb9bbSJeremy L Thompson 3352bfb9bbSJeremy L Thompson CeedVectorCreate(ceed, Xdim*dim, &X); 3452bfb9bbSJeremy L Thompson CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 3552bfb9bbSJeremy L Thompson CeedVectorCreate(ceed, Pdim*dim, &Xq); 3652bfb9bbSJeremy L Thompson CeedVectorSetValue(Xq, 0); 3752bfb9bbSJeremy L Thompson CeedVectorCreate(ceed, Pdim, &U); 3852bfb9bbSJeremy L Thompson CeedVectorCreate(ceed, Qdim*dim, &Uq); 3952bfb9bbSJeremy L Thompson CeedVectorSetValue(Uq, 0); 4052bfb9bbSJeremy L Thompson CeedVectorCreate(ceed, Qdim*dim, &Ones); 4152bfb9bbSJeremy L Thompson CeedVectorSetValue(Ones, 1); 4252bfb9bbSJeremy L Thompson CeedVectorCreate(ceed, Pdim, &Gtposeones); 4352bfb9bbSJeremy L Thompson CeedVectorSetValue(Gtposeones, 0); 4452bfb9bbSJeremy L Thompson 4552bfb9bbSJeremy L Thompson // Get function values at quadrature points 4652bfb9bbSJeremy L Thompson CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, P, 4752bfb9bbSJeremy L Thompson CEED_GAUSS_LOBATTO, &bxl); 4852bfb9bbSJeremy L Thompson CeedBasisApply(bxl, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, Xq); 4952bfb9bbSJeremy L Thompson 5052bfb9bbSJeremy L Thompson CeedVectorGetArrayRead(Xq, CEED_MEM_HOST, &xq); 5152bfb9bbSJeremy L Thompson for (CeedInt i=0; i<Pdim; i++) { 5252bfb9bbSJeremy L Thompson CeedScalar xx[dim]; 5352bfb9bbSJeremy L Thompson for (CeedInt d=0; d<dim; d++) 5452bfb9bbSJeremy L Thompson xx[d] = xq[d*Pdim + i]; 5552bfb9bbSJeremy L Thompson u[i] = Eval(dim, xx); 5652bfb9bbSJeremy L Thompson } 5752bfb9bbSJeremy L Thompson CeedVectorRestoreArrayRead(Xq, &xq); 5852bfb9bbSJeremy L Thompson CeedVectorSetArray(U, CEED_MEM_HOST, CEED_USE_POINTER, u); 5952bfb9bbSJeremy L Thompson 6052bfb9bbSJeremy L Thompson // Calculate G u at quadrature points, G' * 1 at dofs 6152bfb9bbSJeremy L Thompson CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P, Q, 6252bfb9bbSJeremy L Thompson CEED_GAUSS_LOBATTO, &bug); 6352bfb9bbSJeremy L Thompson CeedBasisApply(bug, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, U, Uq); 6452bfb9bbSJeremy L Thompson CeedBasisApply(bug, 1, CEED_TRANSPOSE, CEED_EVAL_GRAD, Ones, Gtposeones); 6552bfb9bbSJeremy L Thompson 6652bfb9bbSJeremy L Thompson // Check if 1' * G * u = u' * (G' * 1) 6752bfb9bbSJeremy L Thompson CeedVectorGetArrayRead(Gtposeones, CEED_MEM_HOST, >poseones); 6852bfb9bbSJeremy L Thompson CeedVectorGetArrayRead(Uq, CEED_MEM_HOST, &uq); 6952bfb9bbSJeremy L Thompson for (CeedInt i=0; i<Pdim; i++) 7052bfb9bbSJeremy L Thompson sum1 += gtposeones[i]*u[i]; 7152bfb9bbSJeremy L Thompson for (CeedInt i=0; i<dim*Qdim; i++) 7252bfb9bbSJeremy L Thompson sum2 += uq[i]; 7352bfb9bbSJeremy L Thompson CeedVectorRestoreArrayRead(Gtposeones, >poseones); 7452bfb9bbSJeremy L Thompson CeedVectorRestoreArrayRead(Uq, &uq); 7552bfb9bbSJeremy L Thompson if (fabs(sum1 - sum2) > 1e-10) 76*c042f62fSJeremy L Thompson // LCOV_EXCL_START 7752bfb9bbSJeremy L Thompson printf("[%d] %f != %f\n", dim, sum1, sum2); 78*c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 7952bfb9bbSJeremy L Thompson 8052bfb9bbSJeremy L Thompson CeedVectorDestroy(&X); 8152bfb9bbSJeremy L Thompson CeedVectorDestroy(&Xq); 8252bfb9bbSJeremy L Thompson CeedVectorDestroy(&U); 8352bfb9bbSJeremy L Thompson CeedVectorDestroy(&Uq); 8452bfb9bbSJeremy L Thompson CeedVectorDestroy(&Ones); 8552bfb9bbSJeremy L Thompson CeedVectorDestroy(&Gtposeones); 8652bfb9bbSJeremy L Thompson CeedBasisDestroy(&bxl); 8752bfb9bbSJeremy L Thompson CeedBasisDestroy(&bug); 8852bfb9bbSJeremy L Thompson } 8952bfb9bbSJeremy L Thompson CeedDestroy(&ceed); 9052bfb9bbSJeremy L Thompson return 0; 9152bfb9bbSJeremy L Thompson } 92