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