1 /// @file 2 /// Test integration with a 2D Simplex non-tensor H1 basis 3 /// \test Test integration 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 int main(int argc, char **argv) { 13 Ceed ceed; 14 const CeedInt P = 6, Q = 4, dim = 2; 15 CeedBasis b; 16 CeedScalar qref[dim*Q], qweight[Q]; 17 CeedScalar interp[P*Q], grad[dim*P*Q], weights[Q]; 18 CeedScalar xr[] = {0., 0.5, 1., 0., 0.5, 0., 0., 0., 0., 0.5, 0.5, 1.}; 19 CeedScalar in[P], out[Q], sum; 20 21 buildmats(qref, qweight, interp, grad); 22 23 CeedInit(argv[1], &ceed); 24 CeedBasisCreateH1(ceed, CEED_TRIANGLE, 1, P, Q, interp, grad, qref, qweight, 25 &b); 26 27 // Interpolate function to quadrature points 28 for (int i=0; i<P; i++) 29 in[i] = feval(xr[0*P+i], xr[1*P+i]); 30 CeedBasisApply(b, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, in, out); 31 CeedBasisApply(b, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, NULL, weights); 32 33 // Check values at quadrature points 34 sum = 0; 35 for (int i=0; i<Q; i++) 36 sum += out[i]*weights[i]; 37 if (fabs(sum - 17./24.) > 1e-10) 38 printf("%f != %f\n", sum, 17./24.); 39 40 CeedBasisDestroy(&b); 41 CeedDestroy(&ceed); 42 return 0; 43 } 44