xref: /libCEED/tests/t203-elemrestriction.c (revision 0acb07cd915a46efccf10dc99d60f8eab5d4016d)
14411cf47Sjeremylt /// @file
24411cf47Sjeremylt /// Test creation, use, and destruction of a blocked element restriction with multiple components in the lvector
34411cf47Sjeremylt /// \test Test creation, use, and destruction of a blocked element restriction with multiple components in the lvector
457c64913Sjeremylt #include <ceed.h>
5*0acb07cdSJeremy L Thompson #include <ceed/backend.h>
657c64913Sjeremylt 
757c64913Sjeremylt int main(int argc, char **argv) {
857c64913Sjeremylt   Ceed ceed;
957c64913Sjeremylt   CeedVector x, y;
10d1d35e2fSjeremylt   CeedInt num_elem = 8;
11*0acb07cdSJeremy L Thompson   CeedInt elem_size = 2;
12*0acb07cdSJeremy L Thompson   CeedInt num_blk = 2;
13d1d35e2fSjeremylt   CeedInt blk_size = 5;
14d1d35e2fSjeremylt   CeedInt num_comp = 3;
15*0acb07cdSJeremy L Thompson   CeedInt ind[elem_size*num_elem];
16d1d35e2fSjeremylt   CeedScalar a[num_comp*(num_elem + 1)];
17*0acb07cdSJeremy L Thompson   const CeedScalar *xx, *yy;
18*0acb07cdSJeremy L Thompson   CeedInt layout[3];
1957c64913Sjeremylt   CeedElemRestriction r;
2057c64913Sjeremylt 
2157c64913Sjeremylt   CeedInit(argv[1], &ceed);
22288c0443SJeremy L Thompson 
23*0acb07cdSJeremy L Thompson   CeedVectorCreate(ceed, num_comp*(num_elem+1), &x);
24*0acb07cdSJeremy L Thompson   for (CeedInt i=0; i<num_elem+1; i++) {
25d1d35e2fSjeremylt     a[i+0*(num_elem+1)] = 10 + i;
26d1d35e2fSjeremylt     a[i+1*(num_elem+1)] = 20 + i;
27d1d35e2fSjeremylt     a[i+2*(num_elem+1)] = 30 + i;
2857c64913Sjeremylt   }
2957c64913Sjeremylt   CeedVectorSetArray(x, CEED_MEM_HOST, CEED_USE_POINTER, a);
30*0acb07cdSJeremy L Thompson 
31d1d35e2fSjeremylt   for (CeedInt i=0; i<num_elem; i++) {
3257c64913Sjeremylt     ind[2*i+0] = i;
3357c64913Sjeremylt     ind[2*i+1] = i+1;
3457c64913Sjeremylt   }
35*0acb07cdSJeremy L Thompson   CeedElemRestrictionCreateBlocked(ceed, num_elem, elem_size, blk_size, num_comp,
36*0acb07cdSJeremy L Thompson                                    num_elem+1, num_comp*(num_elem+1), CEED_MEM_HOST,
37d979a051Sjeremylt                                    CEED_USE_POINTER, ind, &r);
38*0acb07cdSJeremy L Thompson   CeedVectorCreate(ceed, num_comp*num_blk*blk_size*elem_size, &y);
3973d26085Sjeremylt   CeedVectorSetValue(y, 0); // Allocates array
4057c64913Sjeremylt 
4157c64913Sjeremylt   // NoTranspose
42a8d32208Sjeremylt   CeedElemRestrictionApply(r, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE);
43*0acb07cdSJeremy L Thompson   CeedVectorGetArrayRead(y, CEED_MEM_HOST, &yy);
44*0acb07cdSJeremy L Thompson   CeedElemRestrictionGetELayout(r, &layout);
45*0acb07cdSJeremy L Thompson   for (CeedInt i=0; i<elem_size; i++)      // Node
46*0acb07cdSJeremy L Thompson     for (CeedInt j=0; j<num_comp; j++)     // Component
47*0acb07cdSJeremy L Thompson       for (CeedInt k=0; k<num_elem; k++) { // Element
48*0acb07cdSJeremy L Thompson         CeedInt block = k / blk_size;
49*0acb07cdSJeremy L Thompson         CeedInt elem = k % blk_size;
50*0acb07cdSJeremy L Thompson         CeedInt index = (i*blk_size+elem)*layout[0] + j*layout[1]*blk_size +
51*0acb07cdSJeremy L Thompson                         block*layout[2]*blk_size;
52*0acb07cdSJeremy L Thompson         if (yy[index] != a[ind[k*elem_size + i]+j*(num_elem+1)])
53*0acb07cdSJeremy L Thompson           // LCOV_EXCL_START
54*0acb07cdSJeremy L Thompson           printf("Error in restricted array y[%d][%d][%d] = %f\n",
55*0acb07cdSJeremy L Thompson                  i, j, k, (double)yy[index]);
56*0acb07cdSJeremy L Thompson         // LCOV_EXCL_STOP
57*0acb07cdSJeremy L Thompson       }
58*0acb07cdSJeremy L Thompson   CeedVectorRestoreArrayRead(y, &yy);
5957c64913Sjeremylt 
6057c64913Sjeremylt   // Transpose
61*0acb07cdSJeremy L Thompson   CeedVectorSetValue(x, 0);
62a8d32208Sjeremylt   CeedElemRestrictionApply(r, CEED_TRANSPOSE, y, x, CEED_REQUEST_IMMEDIATE);
63*0acb07cdSJeremy L Thompson   CeedVectorGetArrayRead(x, CEED_MEM_HOST, &xx);
64*0acb07cdSJeremy L Thompson   for (CeedInt i=0; i<num_elem+1; i++) {
65*0acb07cdSJeremy L Thompson     for (CeedInt j=0; j<num_comp; j++) {
66*0acb07cdSJeremy L Thompson       if (xx[i+j*(num_elem+1)] != ((j+1)*10+i)*(i > 0 && i < num_elem ? 2.0 : 1.0))
67*0acb07cdSJeremy L Thompson         // LCOV_EXCL_START
68*0acb07cdSJeremy L Thompson         printf("Error in restricted array x[%d][%d] = %f\n",
69*0acb07cdSJeremy L Thompson                j, i, (double)xx[i+j*(num_elem+1)]);
70*0acb07cdSJeremy L Thompson       // LCOV_EXCL_STOP
71*0acb07cdSJeremy L Thompson     }
72*0acb07cdSJeremy L Thompson   }
73*0acb07cdSJeremy L Thompson   CeedVectorRestoreArrayRead(x, &xx);
7457c64913Sjeremylt 
7557c64913Sjeremylt   CeedVectorDestroy(&x);
7657c64913Sjeremylt   CeedVectorDestroy(&y);
7757c64913Sjeremylt   CeedElemRestrictionDestroy(&r);
7857c64913Sjeremylt   CeedDestroy(&ceed);
7957c64913Sjeremylt   return 0;
8057c64913Sjeremylt }
81