1 /// @file 2 /// Test GetInterp and BasisApply for a 2D Quad non-tensor H(div) basis 3 /// \test Test GetInterp and BasisApply for a 2D Quad non-tensor H(div) basis 4 #include <ceed.h> 5 #include <math.h> 6 #include "t330-basis.h" 7 8 int main(int argc, char **argv) { 9 Ceed ceed; 10 const CeedInt num_nodes = 4, Q = 3, dim = 2, num_qpts = Q*Q; 11 CeedInt num_comp = 1; // one vector componenet 12 CeedInt P = dim*num_nodes; // dof per element! 13 CeedBasis b; 14 CeedScalar q_ref[dim*num_qpts], q_weights[num_qpts]; 15 CeedScalar div[P*num_qpts], interp[P*dim*num_qpts]; 16 CeedVector X, Y; 17 const CeedScalar *y, *x; 18 19 CeedInit(argv[1], &ceed); 20 21 HdivBasisQuad(Q, q_ref, q_weights, interp, div, CEED_GAUSS); 22 CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_QUAD, num_comp, P, num_qpts, interp, 23 div, q_ref, q_weights, &b); 24 25 // Test GetInterp for H(div) 26 const CeedScalar *interp2; 27 CeedBasisGetInterp(b, &interp2); 28 for (CeedInt i=0; i<P*dim*num_qpts; i++) { 29 if (fabs(interp[i] - interp2[i]) > 100.*CEED_EPSILON) 30 // LCOV_EXCL_START 31 printf("%f != %f\n", interp[i], interp2[i]); 32 // LCOV_EXCL_STOP 33 } 34 35 CeedVectorCreate(ceed, P, &X); 36 CeedVectorSetValue(X, 1.0); 37 CeedVectorCreate(ceed, num_qpts*dim, &Y); 38 CeedVectorSetValue(Y, 0.); 39 // BasisApply for H(div): CEED_EVAL_INTERP, NOTRANSPOSE case 40 CeedBasisApply(b, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, Y); 41 42 CeedVectorGetArrayRead(Y, CEED_MEM_HOST, &y); 43 for (CeedInt i=0; i<dim*num_qpts; i++) { 44 if (fabs(q_ref[i] - y[i]) > 100.*CEED_EPSILON) 45 // LCOV_EXCL_START 46 printf("%f != %f\n", q_ref[i], y[i]); 47 // LCOV_EXCL_STOP 48 } 49 CeedVectorRestoreArrayRead(Y, &y); 50 51 CeedVectorSetValue(Y, 1.0); 52 CeedVectorSetValue(X, 0.0); 53 // BasisApply for Hdiv: CEED_EVAL_INTERP, TRANSPOSE case 54 CeedBasisApply(b, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, Y, X); 55 56 CeedVectorGetArrayRead(X, CEED_MEM_HOST, &x); 57 CeedScalar sum = 0.; 58 for (CeedInt i = 0; i<P; i++) { 59 sum += x[i]; 60 } 61 if (fabs(sum) > 100.*CEED_EPSILON) 62 // LCOV_EXCL_START 63 printf("sum of array %f != %f\n", sum, 0.0); 64 // LCOV_EXCL_STOP 65 CeedVectorRestoreArrayRead(X, &x); 66 67 CeedBasisDestroy(&b); 68 CeedVectorDestroy(&X); 69 CeedVectorDestroy(&Y); 70 CeedDestroy(&ceed); 71 return 0; 72 } 73