1 /// @file 2 /// Test interpolation with a 2D Quad non-tensor H(div) basis 3 /// \test Test interpolation with 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 CeedVector u, v; 13 const CeedInt p = 8, q = 3, dim = 2, num_qpts = q * q; 14 CeedBasis basis; 15 CeedScalar q_ref[dim * num_qpts], q_weights[num_qpts]; 16 CeedScalar interp[dim * p * num_qpts], div[p * num_qpts]; 17 18 CeedInit(argv[1], &ceed); 19 20 BuildHdivQuadrilateral(q, q_ref, q_weights, interp, div, CEED_GAUSS); 21 CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_QUAD, 1, p, num_qpts, interp, div, q_ref, q_weights, &basis); 22 23 // Test interpolation for H(div) 24 { 25 const CeedScalar *interp_in_basis; 26 27 CeedBasisGetInterp(basis, &interp_in_basis); 28 for (CeedInt i = 0; i < dim * p * num_qpts; i++) { 29 if (fabs(interp[i] - interp_in_basis[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", interp[i], interp_in_basis[i]); 30 } 31 } 32 33 CeedVectorCreate(ceed, p, &u); 34 CeedVectorSetValue(u, 1.0); 35 CeedVectorCreate(ceed, dim * num_qpts, &v); 36 CeedVectorSetValue(v, 0.0); 37 38 CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, v); 39 40 { 41 const CeedScalar *v_array; 42 43 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); 44 for (CeedInt i = 0; i < dim * num_qpts; i++) { 45 if (fabs(q_ref[i] - v_array[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", q_ref[i], v_array[i]); 46 } 47 CeedVectorRestoreArrayRead(v, &v_array); 48 } 49 50 CeedVectorSetValue(v, 1.0); 51 CeedVectorSetValue(u, 0.0); 52 53 CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, v, u); 54 55 { 56 const CeedScalar *u_array; 57 58 CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array); 59 CeedScalar sum = 0.; 60 for (CeedInt i = 0; i < p; i++) { 61 sum += u_array[i]; 62 } 63 if (fabs(sum) > 100. * CEED_EPSILON) printf("sum of array %f != %f\n", sum, 0.0); 64 CeedVectorRestoreArrayRead(u, &u_array); 65 } 66 67 CeedBasisDestroy(&basis); 68 CeedVectorDestroy(&u); 69 CeedVectorDestroy(&v); 70 CeedDestroy(&ceed); 71 return 0; 72 } 73