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 <stdio.h> 7 8 #include "t330-basis.h" 9 10 int main(int argc, char **argv) { 11 Ceed ceed; 12 const CeedInt num_nodes = 4, q = 3, dim = 2, num_qpts = q * q; 13 CeedInt num_comp = 1; // one vector component 14 CeedInt p = dim * num_nodes; // DoF per element 15 CeedBasis basis; 16 CeedScalar q_ref[dim * num_qpts], q_weights[num_qpts]; 17 CeedScalar div[p * num_qpts], interp[p * dim * num_qpts]; 18 CeedVector u, v; 19 20 CeedInit(argv[1], &ceed); 21 22 BuildHdivQuadrilateral(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, &basis); 24 25 // Test GetInterp for H(div) 26 { 27 const CeedScalar *interp_2; 28 CeedBasisGetInterp(basis, &interp_2); 29 for (CeedInt i = 0; i < p * dim * num_qpts; i++) { 30 if (fabs(interp[i] - interp_2[i]) > 100. * CEED_EPSILON) { 31 // LCOV_EXCL_START 32 printf("%f != %f\n", interp[i], interp_2[i]); 33 // LCOV_EXCL_STOP 34 } 35 } 36 } 37 38 CeedVectorCreate(ceed, p, &u); 39 CeedVectorSetValue(u, 1.0); 40 CeedVectorCreate(ceed, num_qpts * dim, &v); 41 CeedVectorSetValue(v, 0.); 42 43 // BasisApply for H(div): CEED_EVAL_INTERP, NOTRANSPOSE case 44 CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, v); 45 46 { 47 const CeedScalar *v_array; 48 49 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); 50 for (CeedInt i = 0; i < dim * num_qpts; i++) { 51 if (fabs(q_ref[i] - v_array[i]) > 100. * CEED_EPSILON) { 52 // LCOV_EXCL_START 53 printf("%f != %f\n", q_ref[i], v_array[i]); 54 // LCOV_EXCL_STOP 55 } 56 } 57 CeedVectorRestoreArrayRead(v, &v_array); 58 } 59 60 CeedVectorSetValue(v, 1.0); 61 CeedVectorSetValue(u, 0.0); 62 63 // BasisApply for Hdiv: CEED_EVAL_INTERP, TRANSPOSE case 64 CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, v, u); 65 66 { 67 const CeedScalar *u_array; 68 69 CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array); 70 CeedScalar sum = 0.; 71 for (CeedInt i = 0; i < p; i++) { 72 sum += u_array[i]; 73 } 74 if (fabs(sum) > 100. * CEED_EPSILON) printf("sum of array %f != %f\n", sum, 0.0); 75 CeedVectorRestoreArrayRead(u, &u_array); 76 } 77 78 CeedBasisDestroy(&basis); 79 CeedVectorDestroy(&u); 80 CeedVectorDestroy(&v); 81 CeedDestroy(&ceed); 82 return 0; 83 } 84