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