150c301a5SRezgar Shakeri /// @file 250c301a5SRezgar Shakeri /// Test GetDiv and BasisApply for a 2D Quad non-tensor H(div) basis 350c301a5SRezgar Shakeri /// \test Test GetDiv and BasisApply for a 2D Quad non-tensor H(div) basis 450c301a5SRezgar Shakeri #include <ceed.h> 550c301a5SRezgar Shakeri #include <math.h> 6*2b730f8bSJeremy L Thompson 750c301a5SRezgar Shakeri #include "t330-basis.h" 850c301a5SRezgar Shakeri 950c301a5SRezgar Shakeri int main(int argc, char **argv) { 1050c301a5SRezgar Shakeri Ceed ceed; 1150c301a5SRezgar Shakeri const CeedInt num_nodes = 4, Q = 3, dim = 2, num_qpts = Q * Q; 1250c301a5SRezgar Shakeri CeedInt num_comp = 1; // one vector componenet 1350c301a5SRezgar Shakeri CeedInt P = dim * num_nodes; // dof per element! 1450c301a5SRezgar Shakeri CeedBasis b; 1550c301a5SRezgar Shakeri CeedScalar q_ref[dim * num_qpts], q_weights[num_qpts]; 1650c301a5SRezgar Shakeri CeedScalar div[P * num_qpts], interp[P * dim * num_qpts]; 1750c301a5SRezgar Shakeri CeedVector X, Y; 1850c301a5SRezgar Shakeri const CeedScalar *y, *x; 1950c301a5SRezgar Shakeri 2050c301a5SRezgar Shakeri CeedInit(argv[1], &ceed); 2150c301a5SRezgar Shakeri 2250c301a5SRezgar Shakeri HdivBasisQuad(Q, q_ref, q_weights, interp, div, CEED_GAUSS); 23*2b730f8bSJeremy L Thompson CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_QUAD, num_comp, P, num_qpts, interp, div, q_ref, q_weights, &b); 2450c301a5SRezgar Shakeri // Test GetDiv 2550c301a5SRezgar Shakeri const CeedScalar *div2; 2650c301a5SRezgar Shakeri CeedBasisGetDiv(b, &div2); 2750c301a5SRezgar Shakeri for (CeedInt i = 0; i < P * num_qpts; i++) { 28*2b730f8bSJeremy L Thompson if (fabs(div[i] - div2[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", div[i], div2[i]); 2950c301a5SRezgar Shakeri } 3050c301a5SRezgar Shakeri CeedVectorCreate(ceed, P, &X); 3150c301a5SRezgar Shakeri CeedVectorSetValue(X, 1); 3250c301a5SRezgar Shakeri CeedVectorCreate(ceed, num_qpts, &Y); 3350c301a5SRezgar Shakeri CeedVectorSetValue(Y, 0); 3450c301a5SRezgar Shakeri // BasisApply for H(div): CEED_EVAL_DIV, NOTRANSPOSE case 3550c301a5SRezgar Shakeri CeedBasisApply(b, 1, CEED_NOTRANSPOSE, CEED_EVAL_DIV, X, Y); 3650c301a5SRezgar Shakeri 3750c301a5SRezgar Shakeri CeedVectorGetArrayRead(Y, CEED_MEM_HOST, &y); 3850c301a5SRezgar Shakeri for (CeedInt i = 0; i < num_qpts; i++) { 39*2b730f8bSJeremy L Thompson if (fabs(P * 0.25 - y[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", 2.0, y[i]); 4050c301a5SRezgar Shakeri } 4150c301a5SRezgar Shakeri CeedVectorRestoreArrayRead(Y, &y); 4250c301a5SRezgar Shakeri 4350c301a5SRezgar Shakeri CeedVectorSetValue(Y, 1.0); 4450c301a5SRezgar Shakeri CeedVectorSetValue(X, 0.0); 4550c301a5SRezgar Shakeri // BasisApply for H(div): CEED_EVAL_DIV, TRANSPOSE case 4650c301a5SRezgar Shakeri CeedBasisApply(b, 1, CEED_TRANSPOSE, CEED_EVAL_DIV, Y, X); 4750c301a5SRezgar Shakeri 4850c301a5SRezgar Shakeri CeedVectorGetArrayRead(X, CEED_MEM_HOST, &x); 4950c301a5SRezgar Shakeri for (CeedInt i = 0; i < P; i++) { 50*2b730f8bSJeremy L Thompson if (fabs(num_qpts * 0.25 - x[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", 2.0, x[i]); 5150c301a5SRezgar Shakeri } 5250c301a5SRezgar Shakeri CeedVectorRestoreArrayRead(X, &x); 5350c301a5SRezgar Shakeri 5450c301a5SRezgar Shakeri CeedBasisDestroy(&b); 5550c301a5SRezgar Shakeri CeedVectorDestroy(&X); 5650c301a5SRezgar Shakeri CeedVectorDestroy(&Y); 5750c301a5SRezgar Shakeri CeedDestroy(&ceed); 5850c301a5SRezgar Shakeri return 0; 5950c301a5SRezgar Shakeri } 60