xref: /libCEED/tests/t231-elemrestriction.c (revision 1ef3a2a99284aa21c437e4d5e30bc6aa0d0f4b0f)
1 /// @file
2 /// Test creation, use, and destruction of an element restriction at points
3 /// \test Test creation, use, and destruction of an element restriction at points
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   CeedInt             num_elem = 3, num_points = num_elem * 2;
11   CeedInt             ind[(num_elem + 1) + num_points];
12   CeedVector          x, y;
13   CeedElemRestriction elem_restriction;
14 
15   CeedInit(argv[1], &ceed);
16 
17   {
18     CeedInt offset      = num_elem + 1;
19     CeedInt point_index = num_elem;
20 
21     for (CeedInt i = 0; i < num_elem; i++) {
22       CeedInt num_points_in_elem = (i + 1) % num_elem + 1;
23 
24       ind[i] = offset;
25       for (CeedInt j = 0; j < num_points_in_elem; j++) {
26         ind[offset + j] = point_index;
27         point_index     = (point_index + 1) % num_points;
28       }
29       offset += num_points_in_elem;
30     }
31     ind[num_elem] = offset;
32   }
33   CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, 1, num_points, CEED_MEM_HOST, CEED_USE_POINTER, ind, &elem_restriction);
34 
35   CeedElemRestrictionCreateVector(elem_restriction, &x, &y);
36   {
37     CeedInt    point_index = num_elem;
38     CeedScalar array[num_points];
39 
40     for (CeedInt i = 0; i < num_elem; i++) {
41       CeedInt num_points_in_elem = (i + 1) % num_elem + 1;
42 
43       for (CeedInt j = 0; j < num_points_in_elem; j++) {
44         array[point_index] = i;
45         point_index        = (point_index + 1) % num_points;
46       }
47     }
48     CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, array);
49   }
50 
51   CeedElemRestrictionApply(elem_restriction, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE);
52   {
53     CeedInt           e_layout[3];
54     const CeedScalar *read_array;
55 
56     CeedVectorGetArrayRead(y, CEED_MEM_HOST, &read_array);
57     CeedElemRestrictionGetELayout(elem_restriction, e_layout);
58 
59     for (CeedInt i = 0; i < num_elem; i++) {
60       CeedSize      elem_offset        = 0;
61       const CeedInt num_points_in_elem = (i + 1) % num_elem + 1;
62 
63       CeedElemRestrictionGetAtPointsElementOffset(elem_restriction, i, &elem_offset);
64       for (CeedInt j = 0; j < num_points_in_elem; j++) {
65         if (i != read_array[elem_offset + j * e_layout[0]]) {
66           // LCOV_EXCL_START
67           printf("Error in restricted array y[%" CeedInt_FMT "] = %f\n != %f\n", (CeedInt)elem_offset + j * e_layout[0],
68                  (CeedScalar)read_array[elem_offset + j * e_layout[0]], (CeedScalar)i);
69           // LCOV_EXCL_STOP
70         }
71       }
72     }
73     CeedVectorRestoreArrayRead(y, &read_array);
74   }
75 
76   CeedVectorDestroy(&x);
77   CeedVectorDestroy(&y);
78   CeedElemRestrictionDestroy(&elem_restriction);
79   CeedDestroy(&ceed);
80   return 0;
81 }
82