1 /// @file 2 /// Test grad in multiple dimensions 3 /// \test Test grad in multiple dimensions 4 #include <ceed.h> 5 #include <math.h> 6 7 static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) { 8 CeedScalar result = tanh(x[0] + 0.1); 9 if (dim > 1) result += atan(x[1] + 0.2); 10 if (dim > 2) result += exp(-(x[2] + 0.3)*(x[2] + 0.3)); 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 CeedVector X, Xq, U, Uq, Ones, Gtposeones; 20 CeedBasis bxl, bug; 21 CeedInt P = 8, Q = 10, Pdim = CeedIntPow(P, dim), Qdim = CeedIntPow(Q, dim), 22 Xdim = CeedIntPow(2, dim); 23 CeedScalar x[Xdim*dim], u[Pdim]; 24 const CeedScalar *xq, *uq, *gtposeones; 25 CeedScalar sum1 = 0, sum2 = 0; 26 27 for (CeedInt d=0; d<dim; d++) { 28 for (CeedInt i=0; i<Xdim; i++) { 29 x[d*Xdim + i] = (i % CeedIntPow(2, dim-d)) / CeedIntPow(2, dim-d-1) ? 1 : -1; 30 } 31 } 32 33 CeedVectorCreate(ceed, Xdim*dim, &X); 34 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, (CeedScalar *)&x); 35 CeedVectorCreate(ceed, Pdim*dim, &Xq); 36 CeedVectorSetValue(Xq, 0); 37 CeedVectorCreate(ceed, Pdim, &U); 38 CeedVectorCreate(ceed, Qdim*dim, &Uq); 39 CeedVectorSetValue(Uq, 0); 40 CeedVectorCreate(ceed, Qdim*dim, &Ones); 41 CeedVectorSetValue(Ones, 1); 42 CeedVectorCreate(ceed, Pdim, &Gtposeones); 43 CeedVectorSetValue(Gtposeones, 0); 44 45 // Get function values at quadrature points 46 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, P, CEED_GAUSS_LOBATTO, &bxl); 47 CeedBasisApply(bxl, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, Xq); 48 49 CeedVectorGetArrayRead(Xq, CEED_MEM_HOST, &xq); 50 for (CeedInt i=0; i<Pdim; i++) { 51 CeedScalar xx[dim]; 52 for (CeedInt d=0; d<dim; d++) xx[d] = xq[d*Pdim + i]; 53 u[i] = Eval(dim, xx); 54 } 55 CeedVectorRestoreArrayRead(Xq, &xq); 56 CeedVectorSetArray(U, CEED_MEM_HOST, CEED_USE_POINTER, (CeedScalar *)&u); 57 58 // Calculate G u at quadrature points, G' * 1 at dofs 59 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P, Q, CEED_GAUSS, &bug); 60 CeedBasisApply(bug, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, U, Uq); 61 CeedBasisApply(bug, 1, CEED_TRANSPOSE, CEED_EVAL_GRAD, Ones, Gtposeones); 62 63 // Check if 1' * G * u = u' * (G' * 1) 64 CeedVectorGetArrayRead(Gtposeones, CEED_MEM_HOST, >poseones); 65 CeedVectorGetArrayRead(Uq, CEED_MEM_HOST, &uq); 66 for (CeedInt i=0; i<Pdim; i++) { 67 sum1 += gtposeones[i]*u[i]; 68 } 69 for (CeedInt i=0; i<dim*Qdim; i++) { 70 sum2 += uq[i]; 71 } 72 CeedVectorRestoreArrayRead(Gtposeones, >poseones); 73 CeedVectorRestoreArrayRead(Uq, &uq); 74 if (fabs(sum1 - sum2) > 1e-10) { 75 printf("[%d] %f != %f\n", dim, sum1, sum2); 76 } 77 78 CeedVectorDestroy(&X); 79 CeedVectorDestroy(&Xq); 80 CeedVectorDestroy(&U); 81 CeedVectorDestroy(&Uq); 82 CeedVectorDestroy(&Ones); 83 CeedVectorDestroy(&Gtposeones); 84 CeedBasisDestroy(&bxl); 85 CeedBasisDestroy(&bug); 86 } 87 CeedDestroy(&ceed); 88 return 0; 89 } 90