1 /// @file 2 /// Test grad with a 2D Simplex non-tensor H1 basis 3 /// \test Test grad with a 2D Simplex non-tensor H1 basis 4 #include <ceed.h> 5 #include <math.h> 6 #include "t310-basis.h" 7 8 double feval(double x1, double x2) { 9 return x1*x1 + x2*x2 + x1*x2 + 1; 10 } 11 12 double dfeval(double x1, double x2) { 13 return 2*x1 + x2; 14 } 15 16 int main(int argc, char **argv) { 17 Ceed ceed; 18 const CeedInt P = 6, Q = 4, dim = 2; 19 CeedBasis b; 20 CeedScalar qref[dim*Q], qweight[Q]; 21 CeedScalar interp[P*Q], grad[dim*P*Q]; 22 CeedScalar xq[] = {0.2, 0.6, 1./3., 0.2, 0.2, 0.2, 1./3., 0.6}; 23 CeedScalar xr[] = {0., 0.5, 1., 0., 0.5, 0., 0., 0., 0., 0.5, 0.5, 1.}; 24 CeedScalar in[P], out[dim*Q], value; 25 26 buildmats(qref, qweight, interp, grad); 27 28 CeedInit(argv[1], &ceed); 29 CeedBasisCreateH1(ceed, CEED_TRIANGLE, 1, P, Q, interp, grad, qref, qweight, 30 &b); 31 32 // Interpolate function to quadrature points 33 for (int i=0; i<P; i++) 34 in[i] = feval(xr[0*P+i], xr[1*P+i]); 35 CeedBasisApply(b, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, in, out); 36 37 // Check values at quadrature points 38 for (int i=0; i<Q; i++) { 39 value = dfeval(xq[0*Q+i], xq[1*Q+i]); 40 if (fabs(out[0*Q+i] - value) > 1e-10) 41 printf("[%d] %f != %f\n", i, out[0*Q+i], value); 42 value = dfeval(xq[1*Q+i], xq[0*Q+i]); 43 if (fabs(out[1*Q+i] - value) > 1e-10) 44 printf("[%d] %f != %f\n", i, out[1*Q+i], value); 45 } 46 47 CeedBasisDestroy(&b); 48 CeedDestroy(&ceed); 49 return 0; 50 } 51