xref: /libCEED/tests/t213-elemrestriction.c (revision bc7bbd5dca53df9bb3e78c945912120808c904be)
1 /// @file
2 /// Test creation, use, and destruction of a blocked element restriction with multiple components in the lvector, and libCEED owns ind pointer
3 /// \test Test creation, use, and destruction of a blocked element restriction with multiple components in the lvector, and libCEED owns ind pointer
4 #include <ceed.h>
5 #include <ceed/backend.h>
6 #include <stdlib.h>
7 #include <string.h>
8 
9 int main(int argc, char **argv) {
10   Ceed ceed;
11   CeedVector x, y;
12   CeedInt num_elem = 8;
13   CeedInt elem_size = 2;
14   CeedInt num_blk = 2;
15   CeedInt blk_size = 5;
16   CeedInt num_comp = 3;
17   CeedInt ind[elem_size*num_elem];
18   CeedInt *ceed_ind = malloc(sizeof(CeedInt)*elem_size*num_elem);
19   CeedScalar a[num_comp*(num_elem + 1)];
20   const CeedScalar *xx, *yy;
21   CeedInt layout[3];
22   CeedElemRestriction r;
23 
24   CeedInit(argv[1], &ceed);
25 
26   CeedVectorCreate(ceed, num_comp*(num_elem+1), &x);
27   for (CeedInt i=0; i<num_elem+1; i++) {
28     a[i+0*(num_elem+1)] = 10 + i;
29     a[i+1*(num_elem+1)] = 20 + i;
30     a[i+2*(num_elem+1)] = 30 + i;
31   }
32   CeedVectorSetArray(x, CEED_MEM_HOST, CEED_USE_POINTER, a);
33 
34   for (CeedInt i=0; i<num_elem; i++) {
35     ind[2*i+0] = i;
36     ind[2*i+1] = i+1;
37   }
38   memcpy(ceed_ind, ind, sizeof(CeedInt)*elem_size*num_elem);
39   CeedElemRestrictionCreateBlocked(ceed, num_elem, elem_size, blk_size, num_comp,
40                                    num_elem+1, num_comp*(num_elem+1), CEED_MEM_HOST,
41                                    CEED_OWN_POINTER, ceed_ind, &r);
42   CeedVectorCreate(ceed, num_comp*num_blk*blk_size*elem_size, &y);
43   CeedVectorSetValue(y, 0); // Allocates array
44 
45   // NoTranspose
46   CeedElemRestrictionApply(r, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE);
47   CeedVectorGetArrayRead(y, CEED_MEM_HOST, &yy);
48   CeedElemRestrictionGetELayout(r, &layout);
49   for (CeedInt i=0; i<elem_size; i++)      // Node
50     for (CeedInt j=0; j<num_comp; j++)     // Component
51       for (CeedInt k=0; k<num_elem; k++) { // Element
52         CeedInt block = k / blk_size;
53         CeedInt elem = k % blk_size;
54         CeedInt index = (i*blk_size+elem)*layout[0] + j*layout[1]*blk_size +
55                         block*layout[2]*blk_size;
56         if (yy[index] != a[ind[k*elem_size + i]+j*(num_elem+1)])
57           // LCOV_EXCL_START
58           printf("Error in restricted array y[%d][%d][%d] = %f\n",
59                  i, j, k, (double)yy[index]);
60         // LCOV_EXCL_STOP
61       }
62   CeedVectorRestoreArrayRead(y, &yy);
63 
64   // Transpose
65   CeedVectorSetValue(x, 0);
66   CeedElemRestrictionApply(r, CEED_TRANSPOSE, y, x, CEED_REQUEST_IMMEDIATE);
67   CeedVectorGetArrayRead(x, CEED_MEM_HOST, &xx);
68   for (CeedInt i=0; i<num_elem+1; i++) {
69     for (CeedInt j=0; j<num_comp; j++) {
70       if (xx[i+j*(num_elem+1)] != ((j+1)*10+i)*(i > 0 && i < num_elem ? 2.0 : 1.0))
71         // LCOV_EXCL_START
72         printf("Error in restricted array x[%d][%d] = %f\n",
73                j, i, (double)xx[i+j*(num_elem+1)]);
74       // LCOV_EXCL_STOP
75     }
76   }
77   CeedVectorRestoreArrayRead(x, &xx);
78 
79   CeedVectorDestroy(&x);
80   CeedVectorDestroy(&y);
81   CeedElemRestrictionDestroy(&r);
82   CeedDestroy(&ceed);
83   return 0;
84 }
85