1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 #include <ceed.h> 18 #include <ceed-backend.h> 19 #include <ceed-impl.h> 20 #include <stdbool.h> 21 #include <stdio.h> 22 23 /// @file 24 /// Implementation of CeedElemRestriction interfaces 25 26 /// ---------------------------------------------------------------------------- 27 /// CeedElemRestriction Library Internal Functions 28 /// ---------------------------------------------------------------------------- 29 /// @addtogroup CeedElemRestrictionDeveloper 30 /// @{ 31 32 /** 33 @brief Permute and pad offsets for a blocked restriction 34 35 @param offsets Array of shape [@a nelem, @a elemsize]. Row i holds the 36 ordered list of the offsets (into the input CeedVector) 37 for the unknowns corresponding to element i, where 38 0 <= i < @a nelem. All offsets must be in the range 39 [0, @a lsize - 1]. 40 @param blkoffsets Array of permuted and padded offsets of 41 shape [@a nblk, @a elemsize, @a blksize]. 42 @param nblk Number of blocks 43 @param nelem Number of elements 44 @param blksize Number of elements in a block 45 @param elemsize Size of each element 46 47 @return An error code: 0 - success, otherwise - failure 48 49 @ref Utility 50 **/ 51 int CeedPermutePadOffsets(const CeedInt *offsets, CeedInt *blkoffsets, 52 CeedInt nblk, CeedInt nelem, CeedInt blksize, 53 CeedInt elemsize) { 54 for (CeedInt e=0; e<nblk*blksize; e+=blksize) 55 for (int j=0; j<blksize; j++) 56 for (int k=0; k<elemsize; k++) 57 blkoffsets[e*elemsize + k*blksize + j] 58 = offsets[CeedIntMin(e+j,nelem-1)*elemsize + k]; 59 return CEED_ERROR_SUCCESS; 60 } 61 62 /// @} 63 64 /// ---------------------------------------------------------------------------- 65 /// CeedElemRestriction Backend API 66 /// ---------------------------------------------------------------------------- 67 /// @addtogroup CeedElemRestrictionBackend 68 /// @{ 69 70 /** 71 @brief Get the Ceed associated with a CeedElemRestriction 72 73 @param rstr CeedElemRestriction 74 @param[out] ceed Variable to store Ceed 75 76 @return An error code: 0 - success, otherwise - failure 77 78 @ref Backend 79 **/ 80 int CeedElemRestrictionGetCeed(CeedElemRestriction rstr, Ceed *ceed) { 81 *ceed = rstr->ceed; 82 return CEED_ERROR_SUCCESS; 83 } 84 85 /** 86 87 @brief Get the strides of a strided CeedElemRestriction 88 89 @param rstr CeedElemRestriction 90 @param[out] strides Variable to store strides array 91 92 @return An error code: 0 - success, otherwise - failure 93 94 @ref Backend 95 **/ 96 int CeedElemRestrictionGetStrides(CeedElemRestriction rstr, 97 CeedInt (*strides)[3]) { 98 if (!rstr->strides) 99 // LCOV_EXCL_START 100 return CeedError(rstr->ceed, CEED_ERROR_MINOR, 101 "ElemRestriction has no stride data"); 102 // LCOV_EXCL_STOP 103 104 for (int i=0; i<3; i++) 105 (*strides)[i] = rstr->strides[i]; 106 return CEED_ERROR_SUCCESS; 107 } 108 109 /** 110 @brief Get read-only access to a CeedElemRestriction offsets array by memtype 111 112 @param rstr CeedElemRestriction to retrieve offsets 113 @param mtype Memory type on which to access the array. If the backend 114 uses a different memory type, this will perform a copy 115 (possibly cached). 116 @param[out] offsets Array on memory type mtype 117 118 @return An error code: 0 - success, otherwise - failure 119 120 @ref User 121 **/ 122 int CeedElemRestrictionGetOffsets(CeedElemRestriction rstr, CeedMemType mtype, 123 const CeedInt **offsets) { 124 int ierr; 125 126 if (!rstr->GetOffsets) 127 // LCOV_EXCL_START 128 return CeedError(rstr->ceed, CEED_ERROR_UNSUPPORTED, 129 "Backend does not support GetOffsets"); 130 // LCOV_EXCL_STOP 131 132 ierr = rstr->GetOffsets(rstr, mtype, offsets); CeedChk(ierr); 133 rstr->numreaders++; 134 return CEED_ERROR_SUCCESS; 135 } 136 137 /** 138 @brief Restore an offsets array obtained using CeedElemRestrictionGetOffsets() 139 140 @param rstr CeedElemRestriction to restore 141 @param offsets Array of offset data 142 143 @return An error code: 0 - success, otherwise - failure 144 145 @ref User 146 **/ 147 int CeedElemRestrictionRestoreOffsets(CeedElemRestriction rstr, 148 const CeedInt **offsets) { 149 *offsets = NULL; 150 rstr->numreaders--; 151 return CEED_ERROR_SUCCESS; 152 } 153 154 /** 155 @brief Get the strided status of a CeedElemRestriction 156 157 @param rstr CeedElemRestriction 158 @param[out] isstrided Variable to store strided status, 1 if strided else 0 159 160 @return An error code: 0 - success, otherwise - failure 161 162 @ref Backend 163 **/ 164 int CeedElemRestrictionIsStrided(CeedElemRestriction rstr, bool *isstrided) { 165 *isstrided = rstr->strides ? true : false; 166 return CEED_ERROR_SUCCESS; 167 } 168 169 /** 170 @brief Get the backend stride status of a CeedElemRestriction 171 172 @param rstr CeedElemRestriction 173 @param[out] status Variable to store stride status 174 175 @return An error code: 0 - success, otherwise - failure 176 177 @ref Backend 178 **/ 179 int CeedElemRestrictionHasBackendStrides(CeedElemRestriction rstr, 180 bool *hasbackendstrides) { 181 if (!rstr->strides) 182 // LCOV_EXCL_START 183 return CeedError(rstr->ceed, CEED_ERROR_MINOR, 184 "ElemRestriction has no stride data"); 185 // LCOV_EXCL_STOP 186 187 *hasbackendstrides = ((rstr->strides[0] == CEED_STRIDES_BACKEND[0]) && 188 (rstr->strides[1] == CEED_STRIDES_BACKEND[1]) && 189 (rstr->strides[2] == CEED_STRIDES_BACKEND[2])); 190 return CEED_ERROR_SUCCESS; 191 } 192 193 /** 194 195 @brief Get the E-vector layout of a CeedElemRestriction 196 197 @param rstr CeedElemRestriction 198 @param[out] layout Variable to store layout array, 199 stored as [nodes, components, elements]. 200 The data for node i, component j, element k in the 201 E-vector is given by 202 i*layout[0] + j*layout[1] + k*layout[2] 203 204 @return An error code: 0 - success, otherwise - failure 205 206 @ref Backend 207 **/ 208 int CeedElemRestrictionGetELayout(CeedElemRestriction rstr, 209 CeedInt (*layout)[3]) { 210 if (!rstr->layout[0]) 211 // LCOV_EXCL_START 212 return CeedError(rstr->ceed, CEED_ERROR_MINOR, 213 "ElemRestriction has no layout data"); 214 // LCOV_EXCL_STOP 215 216 for (int i=0; i<3; i++) 217 (*layout)[i] = rstr->layout[i]; 218 return CEED_ERROR_SUCCESS; 219 } 220 221 /** 222 223 @brief Set the E-vector layout of a CeedElemRestriction 224 225 @param rstr CeedElemRestriction 226 @param layout Variable to containing layout array, 227 stored as [nodes, components, elements]. 228 The data for node i, component j, element k in the 229 E-vector is given by 230 i*layout[0] + j*layout[1] + k*layout[2] 231 232 @return An error code: 0 - success, otherwise - failure 233 234 @ref Backend 235 **/ 236 int CeedElemRestrictionSetELayout(CeedElemRestriction rstr, 237 CeedInt layout[3]) { 238 for (int i = 0; i<3; i++) 239 rstr->layout[i] = layout[i]; 240 return CEED_ERROR_SUCCESS; 241 } 242 243 /** 244 @brief Get the backend data of a CeedElemRestriction 245 246 @param rstr CeedElemRestriction 247 @param[out] data Variable to store data 248 249 @return An error code: 0 - success, otherwise - failure 250 251 @ref Backend 252 **/ 253 int CeedElemRestrictionGetData(CeedElemRestriction rstr, void *data) { 254 *(void **)data = rstr->data; 255 return CEED_ERROR_SUCCESS; 256 } 257 258 /** 259 @brief Set the backend data of a CeedElemRestriction 260 261 @param[out] rstr CeedElemRestriction 262 @param data Data to set 263 264 @return An error code: 0 - success, otherwise - failure 265 266 @ref Backend 267 **/ 268 int CeedElemRestrictionSetData(CeedElemRestriction rstr, void *data) { 269 rstr->data = data; 270 return CEED_ERROR_SUCCESS; 271 } 272 273 /// @} 274 275 /// @cond DOXYGEN_SKIP 276 static struct CeedElemRestriction_private ceed_elemrestriction_none; 277 /// @endcond 278 279 /// ---------------------------------------------------------------------------- 280 /// CeedElemRestriction Public API 281 /// ---------------------------------------------------------------------------- 282 /// @addtogroup CeedElemRestrictionUser 283 /// @{ 284 285 /// Indicate that the stride is determined by the backend 286 const CeedInt CEED_STRIDES_BACKEND[3] = {}; 287 288 /// Indicate that no CeedElemRestriction is provided by the user 289 const CeedElemRestriction CEED_ELEMRESTRICTION_NONE = 290 &ceed_elemrestriction_none; 291 292 /** 293 @brief Create a CeedElemRestriction 294 295 @param ceed A Ceed object where the CeedElemRestriction will be created 296 @param nelem Number of elements described in the @a offsets array 297 @param elemsize Size (number of "nodes") per element 298 @param ncomp Number of field components per interpolation node 299 (1 for scalar fields) 300 @param compstride Stride between components for the same L-vector "node". 301 Data for node i, component j, element k can be found in 302 the L-vector at index 303 offsets[i + k*elemsize] + j*compstride. 304 @param lsize The size of the L-vector. This vector may be larger than 305 the elements and fields given by this restriction. 306 @param mtype Memory type of the @a offsets array, see CeedMemType 307 @param cmode Copy mode for the @a offsets array, see CeedCopyMode 308 @param offsets Array of shape [@a nelem, @a elemsize]. Row i holds the 309 ordered list of the offsets (into the input CeedVector) 310 for the unknowns corresponding to element i, where 311 0 <= i < @a nelem. All offsets must be in the range 312 [0, @a lsize - 1]. 313 @param[out] rstr Address of the variable where the newly created 314 CeedElemRestriction will be stored 315 316 @return An error code: 0 - success, otherwise - failure 317 318 @ref User 319 **/ 320 int CeedElemRestrictionCreate(Ceed ceed, CeedInt nelem, CeedInt elemsize, 321 CeedInt ncomp, CeedInt compstride, 322 CeedInt lsize, CeedMemType mtype, 323 CeedCopyMode cmode, const CeedInt *offsets, 324 CeedElemRestriction *rstr) { 325 int ierr; 326 327 if (!ceed->ElemRestrictionCreate) { 328 Ceed delegate; 329 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 330 CeedChk(ierr); 331 332 if (!delegate) 333 // LCOV_EXCL_START 334 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 335 "Backend does not support ElemRestrictionCreate"); 336 // LCOV_EXCL_STOP 337 338 ierr = CeedElemRestrictionCreate(delegate, nelem, elemsize, ncomp, 339 compstride, lsize, mtype, cmode, 340 offsets, rstr); CeedChk(ierr); 341 return CEED_ERROR_SUCCESS; 342 } 343 344 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 345 (*rstr)->ceed = ceed; 346 ceed->refcount++; 347 (*rstr)->refcount = 1; 348 (*rstr)->nelem = nelem; 349 (*rstr)->elemsize = elemsize; 350 (*rstr)->ncomp = ncomp; 351 (*rstr)->compstride = compstride; 352 (*rstr)->lsize = lsize; 353 (*rstr)->nblk = nelem; 354 (*rstr)->blksize = 1; 355 ierr = ceed->ElemRestrictionCreate(mtype, cmode, offsets, *rstr); 356 CeedChk(ierr); 357 return CEED_ERROR_SUCCESS; 358 } 359 360 /** 361 @brief Create a strided CeedElemRestriction 362 363 @param ceed A Ceed object where the CeedElemRestriction will be created 364 @param nelem Number of elements described by the restriction 365 @param elemsize Size (number of "nodes") per element 366 @param ncomp Number of field components per interpolation "node" 367 (1 for scalar fields) 368 @param lsize The size of the L-vector. This vector may be larger than 369 the elements and fields given by this restriction. 370 @param strides Array for strides between [nodes, components, elements]. 371 Data for node i, component j, element k can be found in 372 the L-vector at index 373 i*strides[0] + j*strides[1] + k*strides[2]. 374 @a CEED_STRIDES_BACKEND may be used with vectors created 375 by a Ceed backend. 376 @param rstr Address of the variable where the newly created 377 CeedElemRestriction will be stored 378 379 @return An error code: 0 - success, otherwise - failure 380 381 @ref User 382 **/ 383 int CeedElemRestrictionCreateStrided(Ceed ceed, CeedInt nelem, CeedInt elemsize, 384 CeedInt ncomp, CeedInt lsize, 385 const CeedInt strides[3], 386 CeedElemRestriction *rstr) { 387 int ierr; 388 389 if (!ceed->ElemRestrictionCreate) { 390 Ceed delegate; 391 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 392 CeedChk(ierr); 393 394 if (!delegate) 395 // LCOV_EXCL_START 396 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 397 "Backend does not support ElemRestrictionCreate"); 398 // LCOV_EXCL_STOP 399 400 ierr = CeedElemRestrictionCreateStrided(delegate, nelem, elemsize, ncomp, 401 lsize, strides, rstr); 402 CeedChk(ierr); 403 return CEED_ERROR_SUCCESS; 404 } 405 406 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 407 (*rstr)->ceed = ceed; 408 ceed->refcount++; 409 (*rstr)->refcount = 1; 410 (*rstr)->nelem = nelem; 411 (*rstr)->elemsize = elemsize; 412 (*rstr)->ncomp = ncomp; 413 (*rstr)->lsize = lsize; 414 (*rstr)->nblk = nelem; 415 (*rstr)->blksize = 1; 416 ierr = CeedMalloc(3, &(*rstr)->strides); CeedChk(ierr); 417 for (int i=0; i<3; i++) 418 (*rstr)->strides[i] = strides[i]; 419 ierr = ceed->ElemRestrictionCreate(CEED_MEM_HOST, CEED_OWN_POINTER, NULL, 420 *rstr); 421 CeedChk(ierr); 422 return CEED_ERROR_SUCCESS; 423 } 424 425 /** 426 @brief Create a blocked CeedElemRestriction, typically only called by backends 427 428 @param ceed A Ceed object where the CeedElemRestriction will be created. 429 @param nelem Number of elements described in the @a offsets array. 430 @param elemsize Size (number of unknowns) per element 431 @param blksize Number of elements in a block 432 @param ncomp Number of field components per interpolation node 433 (1 for scalar fields) 434 @param compstride Stride between components for the same L-vector "node". 435 Data for node i, component j, element k can be found in 436 the L-vector at index 437 offsets[i + k*elemsize] + j*compstride. 438 @param lsize The size of the L-vector. This vector may be larger than 439 the elements and fields given by this restriction. 440 @param mtype Memory type of the @a offsets array, see CeedMemType 441 @param cmode Copy mode for the @a offsets array, see CeedCopyMode 442 @param offsets Array of shape [@a nelem, @a elemsize]. Row i holds the 443 ordered list of the offsets (into the input CeedVector) 444 for the unknowns corresponding to element i, where 445 0 <= i < @a nelem. All offsets must be in the range 446 [0, @a lsize - 1]. The backend will permute and pad this 447 array to the desired ordering for the blocksize, which is 448 typically given by the backend. The default reordering is 449 to interlace elements. 450 @param rstr Address of the variable where the newly created 451 CeedElemRestriction will be stored 452 453 @return An error code: 0 - success, otherwise - failure 454 455 @ref Backend 456 **/ 457 int CeedElemRestrictionCreateBlocked(Ceed ceed, CeedInt nelem, CeedInt elemsize, 458 CeedInt blksize, CeedInt ncomp, 459 CeedInt compstride, CeedInt lsize, 460 CeedMemType mtype, CeedCopyMode cmode, 461 const CeedInt *offsets, 462 CeedElemRestriction *rstr) { 463 int ierr; 464 CeedInt *blkoffsets; 465 CeedInt nblk = (nelem / blksize) + !!(nelem % blksize); 466 467 if (!ceed->ElemRestrictionCreateBlocked) { 468 Ceed delegate; 469 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 470 CeedChk(ierr); 471 472 if (!delegate) 473 // LCOV_EXCL_START 474 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support " 475 "ElemRestrictionCreateBlocked"); 476 // LCOV_EXCL_STOP 477 478 ierr = CeedElemRestrictionCreateBlocked(delegate, nelem, elemsize, blksize, 479 ncomp, compstride, lsize, mtype, 480 cmode, offsets, rstr); 481 CeedChk(ierr); 482 return CEED_ERROR_SUCCESS; 483 } 484 485 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 486 487 ierr = CeedCalloc(nblk*blksize*elemsize, &blkoffsets); CeedChk(ierr); 488 ierr = CeedPermutePadOffsets(offsets, blkoffsets, nblk, nelem, blksize, 489 elemsize); 490 CeedChk(ierr); 491 492 (*rstr)->ceed = ceed; 493 ceed->refcount++; 494 (*rstr)->refcount = 1; 495 (*rstr)->nelem = nelem; 496 (*rstr)->elemsize = elemsize; 497 (*rstr)->ncomp = ncomp; 498 (*rstr)->compstride = compstride; 499 (*rstr)->lsize = lsize; 500 (*rstr)->nblk = nblk; 501 (*rstr)->blksize = blksize; 502 ierr = ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, 503 (const CeedInt *) blkoffsets, *rstr); CeedChk(ierr); 504 if (cmode == CEED_OWN_POINTER) { 505 ierr = CeedFree(&offsets); CeedChk(ierr); 506 } 507 return CEED_ERROR_SUCCESS; 508 } 509 510 /** 511 @brief Create a blocked strided CeedElemRestriction 512 513 @param ceed A Ceed object where the CeedElemRestriction will be created 514 @param nelem Number of elements described by the restriction 515 @param elemsize Size (number of "nodes") per element 516 @param blksize Number of elements in a block 517 @param ncomp Number of field components per interpolation node 518 (1 for scalar fields) 519 @param lsize The size of the L-vector. This vector may be larger than 520 the elements and fields given by this restriction. 521 @param strides Array for strides between [nodes, components, elements]. 522 Data for node i, component j, element k can be found in 523 the L-vector at index 524 i*strides[0] + j*strides[1] + k*strides[2]. 525 @a CEED_STRIDES_BACKEND may be used with vectors created 526 by a Ceed backend. 527 @param rstr Address of the variable where the newly created 528 CeedElemRestriction will be stored 529 530 @return An error code: 0 - success, otherwise - failure 531 532 @ref User 533 **/ 534 int CeedElemRestrictionCreateBlockedStrided(Ceed ceed, CeedInt nelem, 535 CeedInt elemsize, CeedInt blksize, CeedInt ncomp, CeedInt lsize, 536 const CeedInt strides[3], CeedElemRestriction *rstr) { 537 int ierr; 538 CeedInt nblk = (nelem / blksize) + !!(nelem % blksize); 539 540 if (!ceed->ElemRestrictionCreateBlocked) { 541 Ceed delegate; 542 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 543 CeedChk(ierr); 544 545 if (!delegate) 546 // LCOV_EXCL_START 547 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support " 548 "ElemRestrictionCreateBlocked"); 549 // LCOV_EXCL_STOP 550 551 ierr = CeedElemRestrictionCreateBlockedStrided(delegate, nelem, elemsize, 552 blksize, ncomp, lsize, strides, rstr); 553 CeedChk(ierr); 554 return CEED_ERROR_SUCCESS; 555 } 556 557 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 558 559 (*rstr)->ceed = ceed; 560 ceed->refcount++; 561 (*rstr)->refcount = 1; 562 (*rstr)->nelem = nelem; 563 (*rstr)->elemsize = elemsize; 564 (*rstr)->ncomp = ncomp; 565 (*rstr)->lsize = lsize; 566 (*rstr)->nblk = nblk; 567 (*rstr)->blksize = blksize; 568 ierr = CeedMalloc(3, &(*rstr)->strides); CeedChk(ierr); 569 for (int i=0; i<3; i++) 570 (*rstr)->strides[i] = strides[i]; 571 ierr = ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, 572 NULL, *rstr); CeedChk(ierr); 573 return CEED_ERROR_SUCCESS; 574 } 575 576 /** 577 @brief Create CeedVectors associated with a CeedElemRestriction 578 579 @param rstr CeedElemRestriction 580 @param lvec The address of the L-vector to be created, or NULL 581 @param evec The address of the E-vector to be created, or NULL 582 583 @return An error code: 0 - success, otherwise - failure 584 585 @ref User 586 **/ 587 int CeedElemRestrictionCreateVector(CeedElemRestriction rstr, CeedVector *lvec, 588 CeedVector *evec) { 589 int ierr; 590 CeedInt n, m; 591 m = rstr->lsize; 592 n = rstr->nblk * rstr->blksize * rstr->elemsize * rstr->ncomp; 593 if (lvec) { 594 ierr = CeedVectorCreate(rstr->ceed, m, lvec); CeedChk(ierr); 595 } 596 if (evec) { 597 ierr = CeedVectorCreate(rstr->ceed, n, evec); CeedChk(ierr); 598 } 599 return CEED_ERROR_SUCCESS; 600 } 601 602 /** 603 @brief Restrict an L-vector to an E-vector or apply its transpose 604 605 @param rstr CeedElemRestriction 606 @param tmode Apply restriction or transpose 607 @param u Input vector (of size @a lsize when tmode=@ref CEED_NOTRANSPOSE) 608 @param ru Output vector (of shape [@a nelem * @a elemsize] when 609 tmode=@ref CEED_NOTRANSPOSE). Ordering of the e-vector is decided 610 by the backend. 611 @param request Request or @ref CEED_REQUEST_IMMEDIATE 612 613 @return An error code: 0 - success, otherwise - failure 614 615 @ref User 616 **/ 617 int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode tmode, 618 CeedVector u, CeedVector ru, 619 CeedRequest *request) { 620 CeedInt m,n; 621 int ierr; 622 623 if (tmode == CEED_NOTRANSPOSE) { 624 m = rstr->nblk * rstr->blksize * rstr->elemsize * rstr->ncomp; 625 n = rstr->lsize; 626 } else { 627 m = rstr->lsize; 628 n = rstr->nblk * rstr->blksize * rstr->elemsize * rstr->ncomp; 629 } 630 if (n != u->length) 631 // LCOV_EXCL_START 632 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 633 "Input vector size %d not compatible with " 634 "element restriction (%d, %d)", u->length, m, n); 635 // LCOV_EXCL_STOP 636 if (m != ru->length) 637 // LCOV_EXCL_START 638 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 639 "Output vector size %d not compatible with " 640 "element restriction (%d, %d)", ru->length, m, n); 641 // LCOV_EXCL_STOP 642 ierr = rstr->Apply(rstr, tmode, u, ru, request); CeedChk(ierr); 643 return CEED_ERROR_SUCCESS; 644 } 645 646 /** 647 @brief Restrict an L-vector to a block of an E-vector or apply its transpose 648 649 @param rstr CeedElemRestriction 650 @param block Block number to restrict to/from, i.e. block=0 will handle 651 elements [0 : blksize] and block=3 will handle elements 652 [3*blksize : 4*blksize] 653 @param tmode Apply restriction or transpose 654 @param u Input vector (of size @a lsize when tmode=@ref CEED_NOTRANSPOSE) 655 @param ru Output vector (of shape [@a blksize * @a elemsize] when 656 tmode=@ref CEED_NOTRANSPOSE). Ordering of the e-vector is decided 657 by the backend. 658 @param request Request or @ref CEED_REQUEST_IMMEDIATE 659 660 @return An error code: 0 - success, otherwise - failure 661 662 @ref Backend 663 **/ 664 int CeedElemRestrictionApplyBlock(CeedElemRestriction rstr, CeedInt block, 665 CeedTransposeMode tmode, CeedVector u, 666 CeedVector ru, CeedRequest *request) { 667 CeedInt m,n; 668 int ierr; 669 670 if (tmode == CEED_NOTRANSPOSE) { 671 m = rstr->blksize * rstr->elemsize * rstr->ncomp; 672 n = rstr->lsize; 673 } else { 674 m = rstr->lsize; 675 n = rstr->blksize * rstr->elemsize * rstr->ncomp; 676 } 677 if (n != u->length) 678 // LCOV_EXCL_START 679 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 680 "Input vector size %d not compatible with " 681 "element restriction (%d, %d)", u->length, m, n); 682 // LCOV_EXCL_STOP 683 if (m != ru->length) 684 // LCOV_EXCL_START 685 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 686 "Output vector size %d not compatible with " 687 "element restriction (%d, %d)", ru->length, m, n); 688 // LCOV_EXCL_STOP 689 if (rstr->blksize*block > rstr->nelem) 690 // LCOV_EXCL_START 691 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 692 "Cannot retrieve block %d, element %d > " 693 "total elements %d", block, rstr->blksize*block, 694 rstr->nelem); 695 // LCOV_EXCL_STOP 696 ierr = rstr->ApplyBlock(rstr, block, tmode, u, ru, request); 697 CeedChk(ierr); 698 return CEED_ERROR_SUCCESS; 699 } 700 701 /** 702 @brief Get the L-vector component stride 703 704 @param rstr CeedElemRestriction 705 @param[out] compstride Variable to store component stride 706 707 @return An error code: 0 - success, otherwise - failure 708 709 @ref Backend 710 **/ 711 int CeedElemRestrictionGetCompStride(CeedElemRestriction rstr, 712 CeedInt *compstride) { 713 *compstride = rstr->compstride; 714 return CEED_ERROR_SUCCESS; 715 } 716 717 /** 718 @brief Get the total number of elements in the range of a CeedElemRestriction 719 720 @param rstr CeedElemRestriction 721 @param[out] numelem Variable to store number of elements 722 723 @return An error code: 0 - success, otherwise - failure 724 725 @ref Backend 726 **/ 727 int CeedElemRestrictionGetNumElements(CeedElemRestriction rstr, 728 CeedInt *numelem) { 729 *numelem = rstr->nelem; 730 return CEED_ERROR_SUCCESS; 731 } 732 733 /** 734 @brief Get the size of elements in the CeedElemRestriction 735 736 @param rstr CeedElemRestriction 737 @param[out] elemsize Variable to store size of elements 738 739 @return An error code: 0 - success, otherwise - failure 740 741 @ref Backend 742 **/ 743 int CeedElemRestrictionGetElementSize(CeedElemRestriction rstr, 744 CeedInt *elemsize) { 745 *elemsize = rstr->elemsize; 746 return CEED_ERROR_SUCCESS; 747 } 748 749 /** 750 @brief Get the size of the l-vector for a CeedElemRestriction 751 752 @param rstr CeedElemRestriction 753 @param[out] numnodes Variable to store number of nodes 754 755 @return An error code: 0 - success, otherwise - failure 756 757 @ref Backend 758 **/ 759 int CeedElemRestrictionGetLVectorSize(CeedElemRestriction rstr, 760 CeedInt *lsize) { 761 *lsize = rstr->lsize; 762 return CEED_ERROR_SUCCESS; 763 } 764 765 /** 766 @brief Get the number of components in the elements of a 767 CeedElemRestriction 768 769 @param rstr CeedElemRestriction 770 @param[out] numcomp Variable to store number of components 771 772 @return An error code: 0 - success, otherwise - failure 773 774 @ref Backend 775 **/ 776 int CeedElemRestrictionGetNumComponents(CeedElemRestriction rstr, 777 CeedInt *numcomp) { 778 *numcomp = rstr->ncomp; 779 return CEED_ERROR_SUCCESS; 780 } 781 782 /** 783 @brief Get the number of blocks in a CeedElemRestriction 784 785 @param rstr CeedElemRestriction 786 @param[out] numblock Variable to store number of blocks 787 788 @return An error code: 0 - success, otherwise - failure 789 790 @ref Backend 791 **/ 792 int CeedElemRestrictionGetNumBlocks(CeedElemRestriction rstr, 793 CeedInt *numblock) { 794 *numblock = rstr->nblk; 795 return CEED_ERROR_SUCCESS; 796 } 797 798 /** 799 @brief Get the size of blocks in the CeedElemRestriction 800 801 @param rstr CeedElemRestriction 802 @param[out] blksize Variable to store size of blocks 803 804 @return An error code: 0 - success, otherwise - failure 805 806 @ref Backend 807 **/ 808 int CeedElemRestrictionGetBlockSize(CeedElemRestriction rstr, 809 CeedInt *blksize) { 810 *blksize = rstr->blksize; 811 return CEED_ERROR_SUCCESS; 812 } 813 814 /** 815 @brief Get the multiplicity of nodes in a CeedElemRestriction 816 817 @param rstr CeedElemRestriction 818 @param[out] mult Vector to store multiplicity (of size lsize) 819 820 @return An error code: 0 - success, otherwise - failure 821 822 @ref User 823 **/ 824 int CeedElemRestrictionGetMultiplicity(CeedElemRestriction rstr, 825 CeedVector mult) { 826 int ierr; 827 CeedVector evec; 828 829 // Create and set evec 830 ierr = CeedElemRestrictionCreateVector(rstr, NULL, &evec); CeedChk(ierr); 831 ierr = CeedVectorSetValue(evec, 1.0); CeedChk(ierr); 832 ierr = CeedVectorSetValue(mult, 0.0); CeedChk(ierr); 833 834 // Apply to get multiplicity 835 ierr = CeedElemRestrictionApply(rstr, CEED_TRANSPOSE, evec, mult, 836 CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 837 838 // Cleanup 839 ierr = CeedVectorDestroy(&evec); CeedChk(ierr); 840 return CEED_ERROR_SUCCESS; 841 } 842 843 /** 844 @brief View a CeedElemRestriction 845 846 @param[in] rstr CeedElemRestriction to view 847 @param[in] stream Stream to write; typically stdout/stderr or a file 848 849 @return Error code: 0 - success, otherwise - failure 850 851 @ref User 852 **/ 853 int CeedElemRestrictionView(CeedElemRestriction rstr, FILE *stream) { 854 char stridesstr[500]; 855 if (rstr->strides) 856 sprintf(stridesstr, "[%d, %d, %d]", rstr->strides[0], rstr->strides[1], 857 rstr->strides[2]); 858 else 859 sprintf(stridesstr, "%d", rstr->compstride); 860 861 fprintf(stream, "%sCeedElemRestriction from (%d, %d) to %d elements with %d " 862 "nodes each and %s %s\n", rstr->blksize > 1 ? "Blocked " : "", 863 rstr->lsize, rstr->ncomp, rstr->nelem, rstr->elemsize, 864 rstr->strides ? "strides" : "component stride", stridesstr); 865 return CEED_ERROR_SUCCESS; 866 } 867 868 /** 869 @brief Destroy a CeedElemRestriction 870 871 @param rstr CeedElemRestriction to destroy 872 873 @return An error code: 0 - success, otherwise - failure 874 875 @ref User 876 **/ 877 int CeedElemRestrictionDestroy(CeedElemRestriction *rstr) { 878 int ierr; 879 880 if (!*rstr || --(*rstr)->refcount > 0) return CEED_ERROR_SUCCESS; 881 if ((*rstr)->numreaders) 882 return CeedError((*rstr)->ceed, CEED_ERROR_ACCESS, 883 "Cannot destroy CeedElemRestriction, " 884 "a process has read access to the offset data"); 885 if ((*rstr)->Destroy) { 886 ierr = (*rstr)->Destroy(*rstr); CeedChk(ierr); 887 } 888 ierr = CeedFree(&(*rstr)->strides); CeedChk(ierr); 889 ierr = CeedDestroy(&(*rstr)->ceed); CeedChk(ierr); 890 ierr = CeedFree(rstr); CeedChk(ierr); 891 return CEED_ERROR_SUCCESS; 892 } 893 894 /// @} 895