1 /// @file 2 /// Test GetDiv and BasisApply for a 2D Quad non-tensor H(div) basis 3 /// \test Test GetDiv 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 GetDiv 26 { 27 const CeedScalar *div_2; 28 29 CeedBasisGetDiv(basis, &div_2); 30 for (CeedInt i = 0; i < p * num_qpts; i++) { 31 if (fabs(div[i] - div_2[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", div[i], div_2[i]); 32 } 33 } 34 35 CeedVectorCreate(ceed, p, &u); 36 CeedVectorSetValue(u, 1); 37 CeedVectorCreate(ceed, num_qpts, &v); 38 CeedVectorSetValue(v, 0); 39 40 // BasisApply for H(div): CEED_EVAL_DIV, NOTRANSPOSE case 41 CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_DIV, u, v); 42 43 { 44 const CeedScalar *v_array; 45 46 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); 47 for (CeedInt i = 0; i < num_qpts; i++) { 48 if (fabs(p * 0.25 - v_array[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", 2.0, v_array[i]); 49 } 50 CeedVectorRestoreArrayRead(v, &v_array); 51 } 52 53 CeedVectorSetValue(v, 1.0); 54 CeedVectorSetValue(u, 0.0); 55 56 // BasisApply for H(div): CEED_EVAL_DIV, TRANSPOSE case 57 CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_DIV, v, u); 58 59 { 60 const CeedScalar *u_array; 61 62 CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array); 63 for (CeedInt i = 0; i < p; i++) { 64 if (fabs(num_qpts * 0.25 - u_array[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", 2.0, u_array[i]); 65 } 66 CeedVectorRestoreArrayRead(u, &u_array); 67 } 68 69 CeedBasisDestroy(&basis); 70 CeedVectorDestroy(&u); 71 CeedVectorDestroy(&v); 72 CeedDestroy(&ceed); 73 return 0; 74 } 75