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