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