1 /// @file
2 /// Test integration ApplyAdd with a 2D Simplex non-tensor H^1 basis
3 /// \test Test integration ApplyAdd with a 2D Simplex non-tensor H^1 basis
4 #include <ceed.h>
5 #include <math.h>
6 #include <stdio.h>
7
8 #include "t320-basis.h"
9
10 // main test
main(int argc,char ** argv)11 int main(int argc, char **argv) {
12 Ceed ceed;
13 CeedVector u, v, u_q, v_q, w_q;
14 const CeedInt p = 6, q = 4, dim = 2;
15 CeedBasis basis;
16 CeedScalar q_ref[dim * q], q_weight[q];
17 CeedScalar interp[p * q], grad[dim * p * q];
18
19 CeedInit(argv[1], &ceed);
20
21 CeedVectorCreate(ceed, p, &u);
22 CeedVectorCreate(ceed, p, &v);
23 CeedVectorSetValue(u, 1.0);
24 CeedVectorSetValue(v, 0.0);
25 CeedVectorCreate(ceed, q, &u_q);
26 CeedVectorCreate(ceed, q, &v_q);
27 CeedVectorCreate(ceed, q, &w_q);
28
29 Build2DSimplex(q_ref, q_weight, interp, grad);
30 CeedBasisCreateH1(ceed, CEED_TOPOLOGY_TRIANGLE, 1, p, q, interp, grad, q_ref, q_weight, &basis);
31
32 // Compute area
33 CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, u_q);
34 CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, w_q);
35 CeedVectorPointwiseMult(v_q, u_q, w_q);
36 CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, v_q, v);
37 // Double area computed
38 CeedBasisApplyAdd(basis, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, v_q, v);
39
40 // Check area
41 {
42 const CeedScalar *v_array;
43 CeedScalar area = 0.0;
44
45 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
46 for (CeedInt i = 0; i < p; i++) area += v_array[i];
47 if (fabs(area - 1.0) > 1E-6) printf("Incorrect area computed %f != %f\n", area, 1.0);
48 CeedVectorRestoreArrayRead(v, &v_array);
49 }
50
51 CeedVectorDestroy(&u);
52 CeedVectorDestroy(&v);
53 CeedVectorDestroy(&u_q);
54 CeedVectorDestroy(&v_q);
55 CeedVectorDestroy(&w_q);
56 CeedBasisDestroy(&basis);
57 CeedDestroy(&ceed);
58 return 0;
59 }
60