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