xref: /libCEED/tests/t218-elemrestriction.c (revision bd403d52388e064913153593f96f75d758a54e86)
1 /// @file
2 /// Test creation, use, and destruction of a curl-conforming oriented element restriction
3 /// \test Test creation, use, and destruction of a curl-conforming oriented element restriction
4 #include <ceed.h>
5 #include <stdio.h>
6 
7 int main(int argc, char **argv) {
8   Ceed                ceed;
9   CeedVector          x, y;
10   CeedInt             num_elem = 6, elem_size = 2;
11   CeedInt             ind[elem_size * num_elem], curl_orients[3 * elem_size * num_elem];
12   CeedScalar          x_array[num_elem + 1];
13   CeedElemRestriction elem_restriction;
14 
15   CeedInit(argv[1], &ceed);
16 
17   CeedVectorCreate(ceed, num_elem + 1, &x);
18   for (CeedInt i = 0; i < num_elem + 1; i++) x_array[i] = 10 + i;
19   CeedVectorSetArray(x, CEED_MEM_HOST, CEED_USE_POINTER, x_array);
20   CeedVectorCreate(ceed, num_elem * elem_size, &y);
21 
22   for (CeedInt i = 0; i < num_elem; i++) {
23     ind[2 * i + 0]          = i;
24     ind[2 * i + 1]          = i + 1;
25     curl_orients[3 * 2 * i] = curl_orients[3 * 2 * (i + 1) - 1] = 0;
26     if (i % 2 > 0) {
27       // T = [0  -1]
28       //     [-1  0]
29       curl_orients[3 * 2 * i + 1] = 0;
30       curl_orients[3 * 2 * i + 2] = -1;
31       curl_orients[3 * 2 * i + 3] = -1;
32       curl_orients[3 * 2 * i + 4] = 0;
33     } else {
34       // T = I
35       curl_orients[3 * 2 * i + 1] = 1;
36       curl_orients[3 * 2 * i + 2] = 0;
37       curl_orients[3 * 2 * i + 3] = 0;
38       curl_orients[3 * 2 * i + 4] = 1;
39     }
40   }
41   CeedElemRestrictionCreateCurlOriented(ceed, num_elem, elem_size, 1, 1, num_elem + 1, CEED_MEM_HOST, CEED_USE_POINTER, ind, curl_orients,
42                                         &elem_restriction);
43 
44   // NoTranspose
45   CeedElemRestrictionApply(elem_restriction, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE);
46   {
47     const CeedScalar *y_array;
48 
49     CeedVectorGetArrayRead(y, CEED_MEM_HOST, &y_array);
50     for (CeedInt i = 0; i < num_elem; i++) {
51       for (CeedInt j = 0; j < elem_size; j++) {
52         CeedInt k = j + elem_size * i;
53         if (i % 2 > 0) {
54           if (j == 0 && 10 + i + 1 != -y_array[k]) {
55             // LCOV_EXCL_START
56             printf("Error in restricted array y[%" CeedInt_FMT "] = %f\n", k, (CeedScalar)y_array[k]);
57             // LCOV_EXCL_STOP
58           } else if (j == 1 && 10 + i != -y_array[k]) {
59             // LCOV_EXCL_START
60             printf("Error in restricted array y[%" CeedInt_FMT "] = %f\n", k, (CeedScalar)y_array[k]);
61             // LCOV_EXCL_STOP
62           }
63         } else {
64           if (10 + (k + 1) / 2 != y_array[k]) {
65             // LCOV_EXCL_START
66             printf("Error in restricted array y[%" CeedInt_FMT "] = %f\n", k, (CeedScalar)y_array[k]);
67             // LCOV_EXCL_STOP
68           }
69         }
70       }
71     }
72     CeedVectorRestoreArrayRead(y, &y_array);
73   }
74 
75   // Transpose
76   CeedVectorSetValue(x, 0);
77   CeedElemRestrictionApply(elem_restriction, CEED_TRANSPOSE, y, x, CEED_REQUEST_IMMEDIATE);
78   {
79     const CeedScalar *x_array;
80 
81     CeedVectorGetArrayRead(x, CEED_MEM_HOST, &x_array);
82     for (CeedInt i = 0; i < num_elem + 1; i++) {
83       if (x_array[i] != (10 + i) * (i > 0 && i < num_elem ? 2.0 : 1.0))
84         printf("Error in restricted array x[%" CeedInt_FMT "] = %f\n", i, (CeedScalar)x_array[i]);
85     }
86     CeedVectorRestoreArrayRead(x, &x_array);
87   }
88 
89   CeedVectorDestroy(&x);
90   CeedVectorDestroy(&y);
91   CeedElemRestrictionDestroy(&elem_restriction);
92   CeedDestroy(&ceed);
93   return 0;
94 }
95