1 /// @file 2 /// Test interpolation with a 2D Simplex non-tensor H1 basis 3 /// \test Test interpolaton with a 2D Simplex non-tensor H1 basis 4 #include <ceed.h> 5 #include <math.h> 6 #include "t320-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 CeedVector In, Out; 15 const CeedInt P = 6, Q = 4, dim = 2; 16 CeedBasis b; 17 CeedScalar q_ref[dim*Q], q_weight[Q]; 18 CeedScalar interp[P*Q], grad[dim*P*Q]; 19 CeedScalar xq[] = {0.2, 0.6, 1./3., 0.2, 0.2, 0.2, 1./3., 0.6}; 20 CeedScalar xr[] = {0., 0.5, 1., 0., 0.5, 0., 0., 0., 0., 0.5, 0.5, 1.}; 21 const CeedScalar *out; 22 CeedScalar in[P], value; 23 24 buildmats(q_ref, q_weight, interp, grad); 25 26 CeedInit(argv[1], &ceed); 27 28 CeedBasisCreateH1(ceed, CEED_TRIANGLE, 1, P, Q, interp, grad, q_ref, 29 q_weight, &b); 30 31 // Interpolate function to quadrature points 32 for (int i=0; i<P; i++) 33 in[i] = feval(xr[0*P+i], xr[1*P+i]); 34 35 CeedVectorCreate(ceed, P, &In); 36 CeedVectorSetArray(In, CEED_MEM_HOST, CEED_USE_POINTER, in); 37 CeedVectorCreate(ceed, Q, &Out); 38 CeedVectorSetValue(Out, 0); 39 40 CeedBasisApply(b, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, In, Out); 41 42 // Check values at quadrature points 43 CeedVectorGetArrayRead(Out, CEED_MEM_HOST, &out); 44 for (int i=0; i<Q; i++) { 45 value = feval(xq[0*Q+i], xq[1*Q+i]); 46 if (fabs(out[i] - value) > 1E-10) 47 // LCOV_EXCL_START 48 printf("[%d] %f != %f\n", i, out[i], value); 49 // LCOV_EXCL_STOP 50 } 51 CeedVectorRestoreArrayRead(Out, &out); 52 53 CeedVectorDestroy(&In); 54 CeedVectorDestroy(&Out); 55 CeedBasisDestroy(&b); 56 CeedDestroy(&ceed); 57 return 0; 58 } 59