xref: /libCEED/tests/t202-elemrestriction.c (revision 0acb07cd915a46efccf10dc99d60f8eab5d4016d)
14411cf47Sjeremylt /// @file
24411cf47Sjeremylt /// Test creation, use, and destruction of a blocked element restriction
34411cf47Sjeremylt /// \test Test creation, use, and destruction of a blocked element restriction
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;
14*0acb07cdSJeremy L Thompson   CeedInt ind[elem_size*num_elem];
15d1d35e2fSjeremylt   CeedScalar a[num_elem + 1];
16*0acb07cdSJeremy L Thompson   const CeedScalar *xx, *yy;
17*0acb07cdSJeremy L Thompson   CeedInt layout[3];
1857c64913Sjeremylt   CeedElemRestriction r;
1957c64913Sjeremylt 
2057c64913Sjeremylt   CeedInit(argv[1], &ceed);
21288c0443SJeremy L Thompson 
22d1d35e2fSjeremylt   CeedVectorCreate(ceed, num_elem+1, &x);
23d1d35e2fSjeremylt   for (CeedInt i=0; i<num_elem+1; i++)
24288c0443SJeremy L Thompson     a[i] = 10 + i;
2557c64913Sjeremylt   CeedVectorSetArray(x, CEED_MEM_HOST, CEED_USE_POINTER, a);
2657c64913Sjeremylt 
27d1d35e2fSjeremylt   for (CeedInt i=0; i<num_elem; i++) {
2857c64913Sjeremylt     ind[2*i+0] = i;
2957c64913Sjeremylt     ind[2*i+1] = i+1;
3057c64913Sjeremylt   }
31*0acb07cdSJeremy L Thompson   CeedElemRestrictionCreateBlocked(ceed, num_elem, elem_size, blk_size, 1, 1,
32*0acb07cdSJeremy L Thompson                                    num_elem+1, CEED_MEM_HOST, CEED_USE_POINTER,
33*0acb07cdSJeremy L Thompson                                    ind, &r);
34*0acb07cdSJeremy L Thompson   CeedVectorCreate(ceed, num_blk*blk_size*elem_size, &y);
3573d26085Sjeremylt   CeedVectorSetValue(y, 0); // Allocates array
3657c64913Sjeremylt 
3757c64913Sjeremylt   // NoTranspose
38a8d32208Sjeremylt   CeedElemRestrictionApply(r, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE);
39*0acb07cdSJeremy L Thompson   CeedVectorGetArrayRead(y, CEED_MEM_HOST, &yy);
40*0acb07cdSJeremy L Thompson   CeedElemRestrictionGetELayout(r, &layout);
41*0acb07cdSJeremy L Thompson   for (CeedInt i=0; i<elem_size; i++)      // Node
42*0acb07cdSJeremy L Thompson     for (CeedInt j=0; j<1; j++)            // Component
43*0acb07cdSJeremy L Thompson       for (CeedInt k=0; k<num_elem; k++) { // Element
44*0acb07cdSJeremy L Thompson         CeedInt block = k / blk_size;
45*0acb07cdSJeremy L Thompson         CeedInt elem = k % blk_size;
46*0acb07cdSJeremy L Thompson         CeedInt index = (i*blk_size+elem)*layout[0] + j*layout[1]*blk_size +
47*0acb07cdSJeremy L Thompson                         block*layout[2]*blk_size;
48*0acb07cdSJeremy L Thompson         if (yy[index] != a[ind[k*elem_size + i]])
49*0acb07cdSJeremy L Thompson           // LCOV_EXCL_START
50*0acb07cdSJeremy L Thompson           printf("Error in restricted array y[%d][%d][%d] = %f\n",
51*0acb07cdSJeremy L Thompson                  i, j, k, (double)yy[index]);
52*0acb07cdSJeremy L Thompson         // LCOV_EXCL_STOP
53*0acb07cdSJeremy L Thompson       }
54*0acb07cdSJeremy L Thompson   CeedVectorRestoreArrayRead(y, &yy);
5557c64913Sjeremylt 
5657c64913Sjeremylt   // Transpose
57*0acb07cdSJeremy L Thompson   CeedVectorSetValue(x, 0);
58a8d32208Sjeremylt   CeedElemRestrictionApply(r, CEED_TRANSPOSE, y, x, CEED_REQUEST_IMMEDIATE);
59*0acb07cdSJeremy L Thompson   CeedVectorGetArrayRead(x, CEED_MEM_HOST, &xx);
60*0acb07cdSJeremy L Thompson   for (CeedInt i=0; i<num_elem+1; i++)
61*0acb07cdSJeremy L Thompson     if (xx[i] != (10+i)*(i > 0 && i < num_elem ? 2.0 : 1.0))
62*0acb07cdSJeremy L Thompson       // LCOV_EXCL_START
63*0acb07cdSJeremy L Thompson       printf("Error in restricted array x[%d] = %f\n",
64*0acb07cdSJeremy L Thompson              i, (double)xx[i]);
65*0acb07cdSJeremy L Thompson   // LCOV_EXCL_STOP
66*0acb07cdSJeremy L Thompson   CeedVectorRestoreArrayRead(x, &xx);
6757c64913Sjeremylt 
6857c64913Sjeremylt   CeedVectorDestroy(&x);
6957c64913Sjeremylt   CeedVectorDestroy(&y);
7057c64913Sjeremylt   CeedElemRestrictionDestroy(&r);
7157c64913Sjeremylt   CeedDestroy(&ceed);
7257c64913Sjeremylt   return 0;
7357c64913Sjeremylt }
74