1 /// @file
2 /// Test curl with a 2D Simplex non-tensor H(curl) basis
3 /// \test Test curl with a 2D Simplex non-tensor H(curl) basis
4 #include <ceed.h>
5 #include <math.h>
6
7 #include "t340-basis.h"
8
main(int argc,char ** argv)9 int main(int argc, char **argv) {
10 Ceed ceed;
11 CeedVector u, v;
12 const CeedInt p = 8, q = 4, dim = 2;
13 CeedBasis basis;
14 CeedScalar q_ref[dim * q], q_weight[q];
15 CeedScalar interp[dim * p * q], curl[p * q];
16 CeedScalar row_sum[dim * q], column_sum[p];
17
18 CeedInit(argv[1], &ceed);
19
20 BuildHcurl2DSimplex(q_ref, q_weight, interp, curl);
21 CeedBasisCreateHcurl(ceed, CEED_TOPOLOGY_TRIANGLE, 1, p, q, interp, curl, q_ref, q_weight, &basis);
22
23 // Test curl for H(curl)
24 {
25 const CeedScalar *curl_in_basis;
26
27 CeedBasisGetCurl(basis, &curl_in_basis);
28 for (CeedInt i = 0; i < p * q; i++) {
29 if (fabs(curl[i] - curl_in_basis[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", curl[i], curl_in_basis[i]);
30 }
31 }
32
33 for (int i = 0; i < q; i++) {
34 row_sum[i] = 0;
35 for (int j = 0; j < p; j++) {
36 row_sum[i] += curl[j + i * p];
37 }
38 }
39 for (int i = 0; i < p; i++) {
40 column_sum[i] = 0;
41 for (int j = 0; j < q; j++) {
42 column_sum[i] += curl[i + j * p];
43 }
44 }
45
46 CeedVectorCreate(ceed, p, &u);
47 CeedVectorSetValue(u, 1.0);
48 CeedVectorCreate(ceed, q, &v);
49 CeedVectorSetValue(v, 0.0);
50
51 CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_CURL, u, v);
52
53 {
54 const CeedScalar *v_array;
55
56 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
57 for (CeedInt i = 0; i < q; i++) {
58 if (fabs(row_sum[i] - v_array[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", row_sum[i], v_array[i]);
59 }
60 CeedVectorRestoreArrayRead(v, &v_array);
61 }
62
63 CeedVectorSetValue(v, 1.0);
64 CeedVectorSetValue(u, 0.0);
65
66 CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_CURL, v, u);
67
68 {
69 const CeedScalar *u_array;
70
71 CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array);
72 for (CeedInt i = 0; i < p; i++) {
73 if (fabs(column_sum[i] - u_array[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", column_sum[i], u_array[i]);
74 }
75 CeedVectorRestoreArrayRead(u, &u_array);
76 }
77
78 CeedBasisDestroy(&basis);
79 CeedVectorDestroy(&u);
80 CeedVectorDestroy(&v);
81 CeedDestroy(&ceed);
82 return 0;
83 }
84