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