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