xref: /libCEED/tests/t206-elemrestriction.c (revision 1ea55a34129e65adf46cc9aad604e9be02563d89)
1 /// @file
2 /// Test creation, transpose use, and destruction of a multi-component element restriction
3 /// \test Test creation, transpose use, and destruction of a multi-component 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 = 5;
12   CeedInt             ind[2 * num_elem];
13   CeedInt             layout[3];
14   CeedScalar          mult;
15   CeedElemRestriction elem_restriction;
16 
17   CeedInit(argv[1], &ceed);
18 
19   // Setup
20   CeedVectorCreate(ceed, 2 * (num_elem * 2), &x);
21   CeedVectorCreate(ceed, 2 * (num_elem + 1), &y);
22   CeedVectorSetValue(y, 0.0);  // Allocates array, transpose mode sums into vector
23 
24   for (CeedInt i = 0; i < num_elem; i++) {
25     ind[2 * i + 0] = i;
26     ind[2 * i + 1] = i + 1;
27   }
28   CeedElemRestrictionCreate(ceed, num_elem, 2, 2, num_elem + 1, 2 * (num_elem + 1), CEED_MEM_HOST, CEED_USE_POINTER, ind, &elem_restriction);
29 
30   // Set x data in backend E-layout
31   CeedElemRestrictionGetELayout(elem_restriction, &layout);
32   {
33     CeedScalar x_array[2 * (num_elem * 2)];
34 
35     for (CeedInt i = 0; i < 2; i++) {             // Node
36       for (CeedInt j = 0; j < 2; j++) {           // Component
37         for (CeedInt k = 0; k < num_elem; k++) {  // Element
38           x_array[i * layout[0] + j * layout[1] + k * layout[2]] = 10 * j + (2 * k + i + 1) / 2;
39         }
40       }
41     }
42     CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
43   }
44 
45   // Restrict
46   CeedElemRestrictionApply(elem_restriction, CEED_TRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE);
47 
48   // Check
49   {
50     const CeedScalar *y_array;
51 
52     CeedVectorGetArrayRead(y, CEED_MEM_HOST, &y_array);
53     for (CeedInt i = 0; i < num_elem + 1; i++) {
54       mult = i > 0 && i < num_elem ? 2 : 1;
55       if (y_array[i] != i * mult) printf("Error in restricted array y[%" CeedInt_FMT "] = %f != %f\n", i, (CeedScalar)y_array[i], i * mult);
56       if (y_array[i + num_elem + 1] != (10 + i) * mult) {
57         // LCOV_EXCL_START
58         printf("Error in restricted array y[%" CeedInt_FMT "] = %f != %f\n", i + num_elem + 1, (CeedScalar)y_array[i + num_elem + 1],
59                (10. + i) * mult);
60         // LCOV_EXCL_STOP
61       }
62     }
63     CeedVectorRestoreArrayRead(y, &y_array);
64   }
65 
66   CeedVectorDestroy(&x);
67   CeedVectorDestroy(&y);
68   CeedElemRestrictionDestroy(&elem_restriction);
69   CeedDestroy(&ceed);
70   return 0;
71 }
72