1 /// @file 2 /// Test grad transpose with a 2D Simplex non-tensor H1 basis 3 /// \test Test grad transposewith a 2D Simplex non-tensor H1 basis 4 #include <ceed.h> 5 #include <math.h> 6 #include "t320-basis.h" 7 8 int main(int argc, char **argv) { 9 Ceed ceed; 10 CeedVector In, Out; 11 const CeedInt P = 6, Q = 4, dim = 2, ncomp = 3; 12 CeedBasis b; 13 CeedScalar qref[dim*Q], qweight[Q]; 14 CeedScalar interp[P*Q], grad[dim*P*Q]; 15 const CeedScalar *out; 16 CeedScalar colsum[P], *in; 17 18 buildmats(qref, qweight, interp, grad); 19 20 CeedInit(argv[1], &ceed); 21 22 for (int i=0; i<P; i++) { 23 colsum[i] = 0; 24 for (int j=0; j<Q*dim; j++) { 25 colsum[i] += grad[i+j*P]; 26 } 27 } 28 29 CeedBasisCreateH1(ceed, CEED_TRIANGLE, ncomp, P, Q, interp, grad, qref, 30 qweight, &b); 31 32 CeedVectorCreate(ceed, Q*dim*ncomp, &In); 33 CeedVectorGetArray(In, CEED_MEM_HOST, &in); 34 for (int d=0; d<dim; d++) 35 for (int n=0; n<ncomp; n++) 36 for (int q=0; q<Q; q++) 37 in[q+(n+d*ncomp)*Q] = n*1.0; 38 CeedVectorRestoreArray(In, &in); 39 CeedVectorCreate(ceed, P*ncomp, &Out); 40 CeedVectorSetValue(Out, 0); 41 42 CeedBasisApply(b, 1, CEED_TRANSPOSE, CEED_EVAL_GRAD, In, Out); 43 44 // Check values at quadrature points 45 CeedVectorGetArrayRead(Out, CEED_MEM_HOST, &out); 46 for (int p=0; p<P; p++) 47 for (int n=0; n<ncomp; n++) 48 if (fabs(n*colsum[p] - out[p+n*P]) > 1E-14) 49 // LCOV_EXCL_START 50 printf("[%d] %f != %f\n", p, out[p+n*P], n*colsum[p]); 51 // LCOV_EXCL_STOP 52 CeedVectorRestoreArrayRead(Out, &out); 53 54 CeedVectorDestroy(&In); 55 CeedVectorDestroy(&Out); 56 CeedBasisDestroy(&b); 57 CeedDestroy(&ceed); 58 return 0; 59 } 60