1 // Test grad in multiple dimensions 2 #include <ceed.h> 3 #include <math.h> 4 5 static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) { 6 CeedScalar result = tanh(x[0] + 0.1); 7 if (dim > 1) result += atan(x[1] + 0.2); 8 if (dim > 2) result += exp(-(x[2] + 0.3)*(x[2] + 0.3)); 9 return result; 10 } 11 12 int main(int argc, char **argv) { 13 Ceed ceed; 14 15 CeedInit(argv[1], &ceed); 16 for (CeedInt dim=1; dim<=3; dim++) { 17 CeedBasis bxl, bug; 18 CeedInt P = 8, Q = 10, Pdim = CeedPowInt(P, dim), Qdim = CeedPowInt(Q, dim), 19 Xdim = CeedPowInt(2, dim); 20 CeedScalar x[Xdim*dim], ones[dim*Qdim], gtposeones[Pdim]; 21 CeedScalar xq[Pdim*dim], uq[dim*Qdim], u[Pdim], sum1 = 0, sum2 = 0; 22 23 for (CeedInt i=0; i<dim*Qdim; i++) ones[i] = 1; 24 for (CeedInt d=0; d<dim; d++) { 25 for (CeedInt i=0; i<Xdim; i++) { 26 x[d*Xdim + i] = (i % CeedPowInt(2, dim-d)) / CeedPowInt(2, dim-d-1) ? 1 : -1; 27 } 28 } 29 30 // Get function values at quadrature points 31 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, P, CEED_GAUSS_LOBATTO, &bxl); 32 CeedBasisApply(bxl, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, xq); 33 for (CeedInt i=0; i<Pdim; i++) { 34 CeedScalar xx[dim]; 35 for (CeedInt d=0; d<dim; d++) xx[d] = xq[d*Pdim + i]; 36 u[i] = Eval(dim, xx); 37 } 38 39 // Calculate G u at quadrature points, G' * 1 at dofs 40 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P, Q, CEED_GAUSS, &bug); 41 CeedBasisApply(bug, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, u, uq); 42 CeedBasisApply(bug, 1, CEED_TRANSPOSE, CEED_EVAL_GRAD, ones, gtposeones); 43 44 // Check if 1' * G * u = u' * (G' * 1) 45 for (CeedInt i=0; i<Pdim; i++) { 46 sum1 += gtposeones[i]*u[i]; 47 } 48 for (CeedInt i=0; i<dim*Qdim; i++) { 49 sum2 += uq[i]; 50 } 51 if (fabs(sum1 - sum2) > 1e-10) { 52 printf("[%d] %f != %f\n", dim, sum1, sum2); 53 } 54 55 CeedBasisDestroy(&bxl); 56 CeedBasisDestroy(&bug); 57 } 58 CeedDestroy(&ceed); 59 return 0; 60 } 61