10036de2cSjeremylt /// @file 20036de2cSjeremylt /// Test creation, use, and destruction of a blocked element restriction with multiple components in the lvector, and libCEED owns ind pointer 30036de2cSjeremylt /// \test Test creation, use, and destruction of a blocked element restriction with multiple components in the lvector, and libCEED owns ind pointer 40036de2cSjeremylt #include <ceed.h> 5*0acb07cdSJeremy L Thompson #include <ceed/backend.h> 6*0acb07cdSJeremy L Thompson #include <stdlib.h> 7*0acb07cdSJeremy L Thompson #include <string.h> 80036de2cSjeremylt 90036de2cSjeremylt int main(int argc, char **argv) { 100036de2cSjeremylt Ceed ceed; 110036de2cSjeremylt CeedVector x, y; 12d1d35e2fSjeremylt CeedInt num_elem = 8; 13*0acb07cdSJeremy L Thompson CeedInt elem_size = 2; 14*0acb07cdSJeremy L Thompson CeedInt num_blk = 2; 15d1d35e2fSjeremylt CeedInt blk_size = 5; 16d1d35e2fSjeremylt CeedInt num_comp = 3; 17*0acb07cdSJeremy L Thompson CeedInt ind[elem_size*num_elem]; 18*0acb07cdSJeremy L Thompson CeedInt *ceed_ind = malloc(sizeof(CeedInt)*elem_size*num_elem); 19d1d35e2fSjeremylt CeedScalar a[num_comp*(num_elem + 1)]; 20*0acb07cdSJeremy L Thompson const CeedScalar *xx, *yy; 21*0acb07cdSJeremy L Thompson CeedInt layout[3]; 220036de2cSjeremylt CeedElemRestriction r; 230036de2cSjeremylt 240036de2cSjeremylt CeedInit(argv[1], &ceed); 250036de2cSjeremylt 26*0acb07cdSJeremy L Thompson CeedVectorCreate(ceed, num_comp*(num_elem+1), &x); 27*0acb07cdSJeremy L Thompson for (CeedInt i=0; i<num_elem+1; i++) { 28d1d35e2fSjeremylt a[i+0*(num_elem+1)] = 10 + i; 29d1d35e2fSjeremylt a[i+1*(num_elem+1)] = 20 + i; 30d1d35e2fSjeremylt a[i+2*(num_elem+1)] = 30 + i; 310036de2cSjeremylt } 320036de2cSjeremylt CeedVectorSetArray(x, CEED_MEM_HOST, CEED_USE_POINTER, a); 33*0acb07cdSJeremy L Thompson 34d1d35e2fSjeremylt for (CeedInt i=0; i<num_elem; i++) { 350036de2cSjeremylt ind[2*i+0] = i; 360036de2cSjeremylt ind[2*i+1] = i+1; 370036de2cSjeremylt } 38*0acb07cdSJeremy L Thompson memcpy(ceed_ind, ind, sizeof(CeedInt)*elem_size*num_elem); 39*0acb07cdSJeremy L Thompson CeedElemRestrictionCreateBlocked(ceed, num_elem, elem_size, blk_size, num_comp, 40*0acb07cdSJeremy L Thompson num_elem+1, num_comp*(num_elem+1), CEED_MEM_HOST, 41*0acb07cdSJeremy L Thompson CEED_OWN_POINTER, ceed_ind, &r); 42*0acb07cdSJeremy L Thompson CeedVectorCreate(ceed, num_comp*num_blk*blk_size*elem_size, &y); 430036de2cSjeremylt CeedVectorSetValue(y, 0); // Allocates array 440036de2cSjeremylt 450036de2cSjeremylt // NoTranspose 460036de2cSjeremylt CeedElemRestrictionApply(r, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE); 47*0acb07cdSJeremy L Thompson CeedVectorGetArrayRead(y, CEED_MEM_HOST, &yy); 48*0acb07cdSJeremy L Thompson CeedElemRestrictionGetELayout(r, &layout); 49*0acb07cdSJeremy L Thompson for (CeedInt i=0; i<elem_size; i++) // Node 50*0acb07cdSJeremy L Thompson for (CeedInt j=0; j<num_comp; j++) // Component 51*0acb07cdSJeremy L Thompson for (CeedInt k=0; k<num_elem; k++) { // Element 52*0acb07cdSJeremy L Thompson CeedInt block = k / blk_size; 53*0acb07cdSJeremy L Thompson CeedInt elem = k % blk_size; 54*0acb07cdSJeremy L Thompson CeedInt index = (i*blk_size+elem)*layout[0] + j*layout[1]*blk_size + 55*0acb07cdSJeremy L Thompson block*layout[2]*blk_size; 56*0acb07cdSJeremy L Thompson if (yy[index] != a[ind[k*elem_size + i]+j*(num_elem+1)]) 57*0acb07cdSJeremy L Thompson // LCOV_EXCL_START 58*0acb07cdSJeremy L Thompson printf("Error in restricted array y[%d][%d][%d] = %f\n", 59*0acb07cdSJeremy L Thompson i, j, k, (double)yy[index]); 60*0acb07cdSJeremy L Thompson // LCOV_EXCL_STOP 61*0acb07cdSJeremy L Thompson } 62*0acb07cdSJeremy L Thompson CeedVectorRestoreArrayRead(y, &yy); 630036de2cSjeremylt 640036de2cSjeremylt // Transpose 65*0acb07cdSJeremy L Thompson CeedVectorSetValue(x, 0); 660036de2cSjeremylt CeedElemRestrictionApply(r, CEED_TRANSPOSE, y, x, CEED_REQUEST_IMMEDIATE); 67*0acb07cdSJeremy L Thompson CeedVectorGetArrayRead(x, CEED_MEM_HOST, &xx); 68*0acb07cdSJeremy L Thompson for (CeedInt i=0; i<num_elem+1; i++) { 69*0acb07cdSJeremy L Thompson for (CeedInt j=0; j<num_comp; j++) { 70*0acb07cdSJeremy L Thompson if (xx[i+j*(num_elem+1)] != ((j+1)*10+i)*(i > 0 && i < num_elem ? 2.0 : 1.0)) 71*0acb07cdSJeremy L Thompson // LCOV_EXCL_START 72*0acb07cdSJeremy L Thompson printf("Error in restricted array x[%d][%d] = %f\n", 73*0acb07cdSJeremy L Thompson j, i, (double)xx[i+j*(num_elem+1)]); 74*0acb07cdSJeremy L Thompson // LCOV_EXCL_STOP 75*0acb07cdSJeremy L Thompson } 76*0acb07cdSJeremy L Thompson } 77*0acb07cdSJeremy L Thompson CeedVectorRestoreArrayRead(x, &xx); 780036de2cSjeremylt 790036de2cSjeremylt CeedVectorDestroy(&x); 800036de2cSjeremylt CeedVectorDestroy(&y); 810036de2cSjeremylt CeedElemRestrictionDestroy(&r); 820036de2cSjeremylt CeedDestroy(&ceed); 830036de2cSjeremylt return 0; 840036de2cSjeremylt } 85