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