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