150c301a5SRezgar Shakeri /// @file 2*b5404d5dSSebastian Grimberg /// Test div with a 2D Quad non-tensor H(div) basis 3*b5404d5dSSebastian Grimberg /// \test Test div with a 2D Quad non-tensor H(div) basis 450c301a5SRezgar Shakeri #include <ceed.h> 550c301a5SRezgar Shakeri #include <math.h> 649aac155SJeremy L Thompson #include <stdio.h> 72b730f8bSJeremy L Thompson 850c301a5SRezgar Shakeri #include "t330-basis.h" 950c301a5SRezgar Shakeri 1050c301a5SRezgar Shakeri int main(int argc, char **argv) { 1150c301a5SRezgar Shakeri Ceed ceed; 12*b5404d5dSSebastian Grimberg CeedVector u, v; 13*b5404d5dSSebastian Grimberg const CeedInt p = 8, q = 3, dim = 2, num_qpts = q * q; 144fee36f0SJeremy L Thompson CeedBasis basis; 1550c301a5SRezgar Shakeri CeedScalar q_ref[dim * num_qpts], q_weights[num_qpts]; 16*b5404d5dSSebastian Grimberg CeedScalar interp[dim * p * num_qpts], div[p * num_qpts]; 1750c301a5SRezgar Shakeri 1850c301a5SRezgar Shakeri CeedInit(argv[1], &ceed); 1950c301a5SRezgar Shakeri 204fee36f0SJeremy L Thompson BuildHdivQuadrilateral(q, q_ref, q_weights, interp, div, CEED_GAUSS); 21*b5404d5dSSebastian Grimberg CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_QUAD, 1, p, num_qpts, interp, div, q_ref, q_weights, &basis); 224fee36f0SJeremy L Thompson 23*b5404d5dSSebastian Grimberg // Test div for H(div) 244fee36f0SJeremy L Thompson { 25*b5404d5dSSebastian Grimberg const CeedScalar *div_in_basis; 264fee36f0SJeremy L Thompson 27*b5404d5dSSebastian Grimberg CeedBasisGetDiv(basis, &div_in_basis); 284fee36f0SJeremy L Thompson for (CeedInt i = 0; i < p * num_qpts; i++) { 29*b5404d5dSSebastian Grimberg if (fabs(div[i] - div_in_basis[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", div[i], div_in_basis[i]); 3050c301a5SRezgar Shakeri } 314fee36f0SJeremy L Thompson } 324fee36f0SJeremy L Thompson 334fee36f0SJeremy L Thompson CeedVectorCreate(ceed, p, &u); 34*b5404d5dSSebastian Grimberg CeedVectorSetValue(u, 1.0); 354fee36f0SJeremy L Thompson CeedVectorCreate(ceed, num_qpts, &v); 36*b5404d5dSSebastian Grimberg CeedVectorSetValue(v, 0.0); 374fee36f0SJeremy L Thompson 384fee36f0SJeremy L Thompson CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_DIV, u, v); 3950c301a5SRezgar Shakeri 404fee36f0SJeremy L Thompson { 414fee36f0SJeremy L Thompson const CeedScalar *v_array; 424fee36f0SJeremy L Thompson 434fee36f0SJeremy L Thompson CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); 4450c301a5SRezgar Shakeri for (CeedInt i = 0; i < num_qpts; i++) { 454fee36f0SJeremy L Thompson if (fabs(p * 0.25 - v_array[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", 2.0, v_array[i]); 4650c301a5SRezgar Shakeri } 474fee36f0SJeremy L Thompson CeedVectorRestoreArrayRead(v, &v_array); 484fee36f0SJeremy L Thompson } 4950c301a5SRezgar Shakeri 504fee36f0SJeremy L Thompson CeedVectorSetValue(v, 1.0); 514fee36f0SJeremy L Thompson CeedVectorSetValue(u, 0.0); 524fee36f0SJeremy L Thompson 534fee36f0SJeremy L Thompson CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_DIV, v, u); 5450c301a5SRezgar Shakeri 554fee36f0SJeremy L Thompson { 564fee36f0SJeremy L Thompson const CeedScalar *u_array; 574fee36f0SJeremy L Thompson 584fee36f0SJeremy L Thompson CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array); 594fee36f0SJeremy L Thompson for (CeedInt i = 0; i < p; i++) { 604fee36f0SJeremy L Thompson if (fabs(num_qpts * 0.25 - u_array[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", 2.0, u_array[i]); 6150c301a5SRezgar Shakeri } 624fee36f0SJeremy L Thompson CeedVectorRestoreArrayRead(u, &u_array); 634fee36f0SJeremy L Thompson } 6450c301a5SRezgar Shakeri 654fee36f0SJeremy L Thompson CeedBasisDestroy(&basis); 664fee36f0SJeremy L Thompson CeedVectorDestroy(&u); 674fee36f0SJeremy L Thompson CeedVectorDestroy(&v); 6850c301a5SRezgar Shakeri CeedDestroy(&ceed); 6950c301a5SRezgar Shakeri return 0; 7050c301a5SRezgar Shakeri } 71