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
main(int argc,char ** argv)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