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-impl.h> 18 #include <ceed-backend.h> 19 20 /// @file 21 /// Implementation of public CeedElemRestriction interfaces 22 /// 23 /// @addtogroup CeedElemRestriction 24 /// @{ 25 26 /** 27 @brief Create a CeedElemRestriction 28 29 @param ceed A Ceed object where the CeedElemRestriction will be created 30 @param imode Ordering of the ncomp components, i.e. it specifies 31 the ordering of the components of the L-vector used 32 by this CeedElemRestriction. CEED_NONINTERLACED indicates 33 the component is the outermost index and CEED_INTERLACED 34 indicates the component is the innermost index in 35 ordering of the L-vector. 36 @param nelem Number of elements described in the @a indices array 37 @param elemsize Size (number of "nodes") per element 38 @param nnodes The number of nodes in the L-vector. The input CeedVector 39 to which the restriction will be applied is of size 40 @a nnodes * @a ncomp. This size may include data 41 used by other CeedElemRestriction objects describing 42 different types of elements. 43 @param ncomp Number of field components per interpolation node 44 (1 for scalar fields) 45 @param mtype Memory type of the @a indices array, see CeedMemType 46 @param cmode Copy mode for the @a indices array, see CeedCopyMode 47 @param indices Array of shape [@a nelem, @a elemsize]. Row i holds the 48 ordered list of the indices (into the input CeedVector) 49 for the unknowns corresponding to element i, where 50 0 <= i < @a nelem. All indices must be in the range 51 [0, @a nnodes - 1]. 52 @param[out] rstr Address of the variable where the newly created 53 CeedElemRestriction will be stored 54 55 @return An error code: 0 - success, otherwise - failure 56 57 @ref Basic 58 **/ 59 int CeedElemRestrictionCreate(Ceed ceed, CeedInterlaceMode imode, 60 CeedInt nelem, 61 CeedInt elemsize, CeedInt nnodes, CeedInt ncomp, 62 CeedMemType mtype, CeedCopyMode cmode, 63 const CeedInt *indices, 64 CeedElemRestriction *rstr) { 65 int ierr; 66 67 if (!ceed->ElemRestrictionCreate) { 68 Ceed delegate; 69 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 70 CeedChk(ierr); 71 72 if (!delegate) 73 // LCOV_EXCL_START 74 return CeedError(ceed, 1, "Backend does not support ElemRestrictionCreate"); 75 // LCOV_EXCL_STOP 76 77 ierr = CeedElemRestrictionCreate(delegate, imode, nelem, elemsize, 78 nnodes, ncomp, mtype, cmode, 79 indices, rstr); CeedChk(ierr); 80 return 0; 81 } 82 83 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 84 (*rstr)->ceed = ceed; 85 ceed->refcount++; 86 (*rstr)->refcount = 1; 87 (*rstr)->imode = imode; 88 (*rstr)->nelem = nelem; 89 (*rstr)->elemsize = elemsize; 90 (*rstr)->nnodes = nnodes; 91 (*rstr)->ncomp = ncomp; 92 (*rstr)->nblk = nelem; 93 (*rstr)->blksize = 1; 94 ierr = ceed->ElemRestrictionCreate(mtype, cmode, indices, *rstr); CeedChk(ierr); 95 return 0; 96 } 97 98 /** 99 @brief Create an identity CeedElemRestriction 100 101 @param ceed A Ceed object where the CeedElemRestriction will be created 102 @param imode Ordering of the ncomp components, i.e. it specifies 103 the ordering of the components of the L-vector used 104 by this CeedElemRestriction. CEED_NONINTERLACED indicates 105 the component is the outermost index and CEED_INTERLACED 106 indicates the component is the innermost index in 107 ordering of the L-vector. 108 @param nelem Number of elements described in the @a indices array 109 @param elemsize Size (number of "nodes") per element 110 @param nnodes The number of nodes in the L-vector. The input CeedVector 111 to which the restriction will be applied is of size 112 @a nnodes * @a ncomp. This size may include data 113 used by other CeedElemRestriction objects describing 114 different types of elements. 115 @param ncomp Number of field components per interpolation node 116 (1 for scalar fields) 117 @param rstr Address of the variable where the newly created 118 CeedElemRestriction will be stored 119 120 @return An error code: 0 - success, otherwise - failure 121 122 @ref Basic 123 **/ 124 int CeedElemRestrictionCreateIdentity(Ceed ceed, CeedInterlaceMode imode, 125 CeedInt nelem, CeedInt elemsize, 126 CeedInt nnodes, CeedInt ncomp, 127 CeedElemRestriction *rstr) { 128 int ierr; 129 130 if (!ceed->ElemRestrictionCreate) { 131 Ceed delegate; 132 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 133 CeedChk(ierr); 134 135 if (!delegate) 136 // LCOV_EXCL_START 137 return CeedError(ceed, 1, "Backend does not support ElemRestrictionCreate"); 138 // LCOV_EXCL_STOP 139 140 ierr = CeedElemRestrictionCreateIdentity(delegate, imode, nelem, elemsize, 141 nnodes, ncomp, rstr); CeedChk(ierr); 142 return 0; 143 } 144 145 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 146 (*rstr)->ceed = ceed; 147 ceed->refcount++; 148 (*rstr)->refcount = 1; 149 (*rstr)->imode = imode; 150 (*rstr)->nelem = nelem; 151 (*rstr)->elemsize = elemsize; 152 (*rstr)->nnodes = nnodes; 153 (*rstr)->ncomp = ncomp; 154 (*rstr)->nblk = nelem; 155 (*rstr)->blksize = 1; 156 ierr = ceed->ElemRestrictionCreate(CEED_MEM_HOST, CEED_OWN_POINTER, NULL, 157 *rstr); 158 CeedChk(ierr); 159 return 0; 160 } 161 162 /** 163 @brief Permute and pad indices for a blocked restriction 164 165 @param indices Array of shape [@a nelem, @a elemsize]. Row i holds the 166 ordered list of the indices (into the input CeedVector) 167 for the unknowns corresponding to element i, where 168 0 <= i < @a nelem. All indices must be in the range 169 [0, @a nnodes). 170 @param blkindices Array of permuted and padded indices of 171 shape [@a nblk, @a elemsize, @a blksize]. 172 @param nblk Number of blocks 173 @param nelem Number of elements 174 @param blksize Number of elements in a block 175 @param elemsize Size of each element 176 177 @return An error code: 0 - success, otherwise - failure 178 179 @ref Utility 180 **/ 181 int CeedPermutePadIndices(const CeedInt *indices, CeedInt *blkindices, 182 CeedInt nblk, CeedInt nelem, CeedInt blksize, 183 CeedInt elemsize) { 184 for (CeedInt e = 0; e < nblk*blksize; e+=blksize) 185 for (int j = 0; j < blksize; j++) 186 for (int k = 0; k < elemsize; k++) 187 blkindices[e*elemsize + k*blksize + j] 188 = indices[CeedIntMin(e+j,nelem-1)*elemsize + k]; 189 return 0; 190 } 191 192 /** 193 @brief Create a blocked CeedElemRestriction, typically only called by backends 194 195 @param ceed A Ceed object where the CeedElemRestriction will be created. 196 @param imode Ordering of the ncomp components, i.e. it specifies 197 the ordering of the components of the L-vector used 198 by this CeedElemRestriction. CEED_NONINTERLACED indicates 199 the component is the outermost index and CEED_INTERLACED 200 indicates the component is the innermost index in 201 ordering of the L-vector. 202 @param nelem Number of elements described in the @a indices array. 203 @param elemsize Size (number of unknowns) per element 204 @param blksize Number of elements in a block 205 @param nnodes The number of nodes in the L-vector. The input CeedVector 206 to which the restriction will be applied is of size 207 @a nnodes * @a ncomp. This size may include data 208 used by other CeedElemRestriction objects describing 209 different types of elements. 210 @param ncomp Number of field components per interpolation node 211 (1 for scalar fields) 212 @param mtype Memory type of the @a indices array, see CeedMemType 213 @param cmode Copy mode for the @a indices array, see CeedCopyMode 214 @param indices Array of shape [@a nelem, @a elemsize]. Row i holds the 215 ordered list of the indices (into the input CeedVector) 216 for the unknowns corresponding to element i, where 217 0 <= i < @a nelem. All indices must be in the range 218 [0, @a nnodes). The backend will permute and pad this 219 array to the desired ordering for the blocksize, which is 220 typically given by the backend. The default reordering is 221 to interlace elements. 222 @param rstr Address of the variable where the newly created 223 CeedElemRestriction will be stored 224 225 @return An error code: 0 - success, otherwise - failure 226 227 @ref Advanced 228 **/ 229 int CeedElemRestrictionCreateBlocked(Ceed ceed, CeedInterlaceMode imode, 230 CeedInt nelem, CeedInt elemsize, 231 CeedInt blksize, CeedInt nnodes, 232 CeedInt ncomp, CeedMemType mtype, 233 CeedCopyMode cmode, const CeedInt *indices, 234 CeedElemRestriction *rstr) { 235 int ierr; 236 CeedInt *blkindices; 237 CeedInt nblk = (nelem / blksize) + !!(nelem % blksize); 238 239 if (!ceed->ElemRestrictionCreateBlocked) { 240 Ceed delegate; 241 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 242 CeedChk(ierr); 243 244 if (!delegate) 245 // LCOV_EXCL_START 246 return CeedError(ceed, 1, "Backend does not support " 247 "ElemRestrictionCreateBlocked"); 248 // LCOV_EXCL_STOP 249 250 ierr = CeedElemRestrictionCreateBlocked(delegate, imode, nelem, elemsize, 251 blksize, nnodes, ncomp, mtype, cmode, 252 indices, rstr); CeedChk(ierr); 253 return 0; 254 } 255 256 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 257 258 if (indices) { 259 ierr = CeedCalloc(nblk*blksize*elemsize, &blkindices); CeedChk(ierr); 260 ierr = CeedPermutePadIndices(indices, blkindices, nblk, nelem, blksize, 261 elemsize); 262 CeedChk(ierr); 263 } else { 264 blkindices = NULL; 265 } 266 267 (*rstr)->ceed = ceed; 268 ceed->refcount++; 269 (*rstr)->refcount = 1; 270 (*rstr)->imode = imode; 271 (*rstr)->nelem = nelem; 272 (*rstr)->elemsize = elemsize; 273 (*rstr)->nnodes = nnodes; 274 (*rstr)->ncomp = ncomp; 275 (*rstr)->nblk = nblk; 276 (*rstr)->blksize = blksize; 277 ierr = ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, 278 (const CeedInt *) blkindices, *rstr); CeedChk(ierr); 279 280 if (cmode == CEED_OWN_POINTER) { 281 ierr = CeedFree(&indices); CeedChk(ierr); 282 } 283 284 return 0; 285 } 286 287 /** 288 @brief Create CeedVectors associated with a CeedElemRestriction 289 290 @param rstr CeedElemRestriction 291 @param lvec The address of the L-vector to be created, or NULL 292 @param evec The address of the E-vector to be created, or NULL 293 294 @return An error code: 0 - success, otherwise - failure 295 296 @ref Advanced 297 **/ 298 int CeedElemRestrictionCreateVector(CeedElemRestriction rstr, CeedVector *lvec, 299 CeedVector *evec) { 300 int ierr; 301 CeedInt n, m; 302 m = rstr->nnodes * rstr->ncomp; 303 n = rstr->nblk * rstr->blksize * rstr->elemsize * rstr->ncomp; 304 if (lvec) { 305 ierr = CeedVectorCreate(rstr->ceed, m, lvec); CeedChk(ierr); 306 } 307 if (evec) { 308 ierr = CeedVectorCreate(rstr->ceed, n, evec); CeedChk(ierr); 309 } 310 return 0; 311 } 312 313 /** 314 @brief Restrict an L-vector to an E-vector or apply its transpose 315 316 @param rstr CeedElemRestriction 317 @param tmode Apply restriction or transpose 318 @param u Input vector (of shape [@a nnodes, @a ncomp] when 319 tmode=CEED_NOTRANSPOSE, imode=CEED_INTERLACED) 320 @param ru Output vector (of shape [@a nelem * @a elemsize] when 321 tmode=CEED_NOTRANSPOSE). Ordering of the e-vector is decided 322 by the backend. 323 @param request Request or CEED_REQUEST_IMMEDIATE 324 325 @return An error code: 0 - success, otherwise - failure 326 327 @ref Advanced 328 **/ 329 int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode tmode, 330 CeedVector u, CeedVector ru, 331 CeedRequest *request) { 332 CeedInt m,n; 333 int ierr; 334 335 if (tmode == CEED_NOTRANSPOSE) { 336 m = rstr->nblk * rstr->blksize * rstr->elemsize * rstr->ncomp; 337 n = rstr->nnodes * rstr->ncomp; 338 } else { 339 m = rstr->nnodes * rstr->ncomp; 340 n = rstr->nblk * rstr->blksize * rstr->elemsize * rstr->ncomp; 341 } 342 if (n != u->length) 343 // LCOV_EXCL_START 344 return CeedError(rstr->ceed, 2, "Input vector size %d not compatible with " 345 "element restriction (%d, %d)", u->length, m, n); 346 // LCOV_EXCL_STOP 347 if (m != ru->length) 348 // LCOV_EXCL_START 349 return CeedError(rstr->ceed, 2, "Output vector size %d not compatible with " 350 "element restriction (%d, %d)", ru->length, m, n); 351 // LCOV_EXCL_STOP 352 ierr = rstr->Apply(rstr, tmode, u, ru, request); CeedChk(ierr); 353 354 return 0; 355 } 356 357 /** 358 @brief Restrict an L-vector to a block of an E-vector or apply its transpose 359 360 @param rstr CeedElemRestriction 361 @param block Block number to restrict to/from, i.e. block=0 will handle 362 elements [0 : blksize] and block=3 will handle elements 363 [3*blksize : 4*blksize] 364 @param tmode Apply restriction or transpose 365 @param u Input vector (of shape [@a nnodes, @a ncomp] when 366 tmode=CEED_NOTRANSPOSE, imode=CEED_INTERLACED) 367 @param ru Output vector (of shape [@a blksize * @a elemsize] when 368 tmode=CEED_NOTRANSPOSE). Ordering of the e-vector is decided 369 by the backend. 370 @param request Request or CEED_REQUEST_IMMEDIATE 371 372 @return An error code: 0 - success, otherwise - failure 373 374 @ref Advanced 375 **/ 376 int CeedElemRestrictionApplyBlock(CeedElemRestriction rstr, CeedInt block, 377 CeedTransposeMode tmode, CeedVector u, 378 CeedVector ru, CeedRequest *request) { 379 CeedInt m,n; 380 int ierr; 381 382 if (tmode == CEED_NOTRANSPOSE) { 383 m = rstr->blksize * rstr->elemsize * rstr->ncomp; 384 n = rstr->nnodes * rstr->ncomp; 385 } else { 386 m = rstr->nnodes * rstr->ncomp; 387 n = rstr->blksize * rstr->elemsize * rstr->ncomp; 388 } 389 if (n != u->length) 390 // LCOV_EXCL_START 391 return CeedError(rstr->ceed, 2, "Input vector size %d not compatible with " 392 "element restriction (%d, %d)", u->length, m, n); 393 // LCOV_EXCL_STOP 394 if (m != ru->length) 395 // LCOV_EXCL_START 396 return CeedError(rstr->ceed, 2, "Output vector size %d not compatible with " 397 "element restriction (%d, %d)", ru->length, m, n); 398 // LCOV_EXCL_STOP 399 if (rstr->blksize*block > rstr->nelem) 400 // LCOV_EXCL_START 401 return CeedError(rstr->ceed, 2, "Cannot retrieve block %d, element %d > " 402 "total elements %d", block, rstr->blksize*block, 403 rstr->nelem); 404 // LCOV_EXCL_STOP 405 ierr = rstr->ApplyBlock(rstr, block, tmode, u, ru, request); 406 CeedChk(ierr); 407 408 return 0; 409 } 410 411 /** 412 @brief Get the multiplicity of nodes in a CeedElemRestriction 413 414 @param rstr CeedElemRestriction 415 @param[out] mult Vector to store multiplicity (of size nnodes*ncomp) 416 417 @return An error code: 0 - success, otherwise - failure 418 419 @ref Advanced 420 **/ 421 int CeedElemRestrictionGetMultiplicity(CeedElemRestriction rstr, 422 CeedVector mult) { 423 int ierr; 424 CeedVector evec; 425 426 // Create and set evec 427 ierr = CeedElemRestrictionCreateVector(rstr, NULL, &evec); CeedChk(ierr); 428 ierr = CeedVectorSetValue(evec, 1.0); CeedChk(ierr); 429 ierr = CeedVectorSetValue(mult, 0.0); CeedChk(ierr); 430 431 // Apply to get multiplicity 432 ierr = CeedElemRestrictionApply(rstr, CEED_TRANSPOSE, evec, mult, 433 CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 434 435 // Cleanup 436 ierr = CeedVectorDestroy(&evec); CeedChk(ierr); 437 438 return 0; 439 } 440 441 /** 442 @brief Get the Ceed associated with a CeedElemRestriction 443 444 @param rstr CeedElemRestriction 445 @param[out] ceed Variable to store Ceed 446 447 @return An error code: 0 - success, otherwise - failure 448 449 @ref Advanced 450 **/ 451 int CeedElemRestrictionGetCeed(CeedElemRestriction rstr, Ceed *ceed) { 452 *ceed = rstr->ceed; 453 return 0; 454 } 455 456 /** 457 @brief Get the L-vector interlaced mode of a CeedElemRestriction 458 459 @param rstr CeedElemRestriction 460 @param[out] imode Variable to store imode 461 462 @return An error code: 0 - success, otherwise - failure 463 464 @ref Advanced 465 **/ 466 int CeedElemRestrictionGetIMode(CeedElemRestriction rstr, 467 CeedInterlaceMode *imode) { 468 *imode = rstr->imode; 469 return 0; 470 } 471 472 /** 473 @brief Get the total number of elements in the range of a CeedElemRestriction 474 475 @param rstr CeedElemRestriction 476 @param[out] numelem Variable to store number of elements 477 478 @return An error code: 0 - success, otherwise - failure 479 480 @ref Advanced 481 **/ 482 int CeedElemRestrictionGetNumElements(CeedElemRestriction rstr, 483 CeedInt *numelem) { 484 *numelem = rstr->nelem; 485 return 0; 486 } 487 488 /** 489 @brief Get the size of elements in the CeedElemRestriction 490 491 @param rstr CeedElemRestriction 492 @param[out] elemsize Variable to store size of elements 493 494 @return An error code: 0 - success, otherwise - failure 495 496 @ref Advanced 497 **/ 498 int CeedElemRestrictionGetElementSize(CeedElemRestriction rstr, 499 CeedInt *elemsize) { 500 *elemsize = rstr->elemsize; 501 return 0; 502 } 503 504 /** 505 @brief Get the number of degrees of freedom in the range of a 506 CeedElemRestriction 507 508 @param rstr CeedElemRestriction 509 @param[out] numnodes Variable to store number of nodes 510 511 @return An error code: 0 - success, otherwise - failure 512 513 @ref Advanced 514 **/ 515 int CeedElemRestrictionGetNumNodes(CeedElemRestriction rstr, 516 CeedInt *numnodes) { 517 *numnodes = rstr->nnodes; 518 return 0; 519 } 520 521 /** 522 @brief Get the number of components in the elements of a 523 CeedElemRestriction 524 525 @param rstr CeedElemRestriction 526 @param[out] numcomp Variable to store number of components 527 528 @return An error code: 0 - success, otherwise - failure 529 530 @ref Advanced 531 **/ 532 int CeedElemRestrictionGetNumComponents(CeedElemRestriction rstr, 533 CeedInt *numcomp) { 534 *numcomp = rstr->ncomp; 535 return 0; 536 } 537 538 /** 539 @brief Get the number of blocks in a CeedElemRestriction 540 541 @param rstr CeedElemRestriction 542 @param[out] numblock Variable to store number of blocks 543 544 @return An error code: 0 - success, otherwise - failure 545 546 @ref Advanced 547 **/ 548 int CeedElemRestrictionGetNumBlocks(CeedElemRestriction rstr, 549 CeedInt *numblock) { 550 *numblock = rstr->nblk; 551 return 0; 552 } 553 554 /** 555 @brief Get the size of blocks in the CeedElemRestriction 556 557 @param rstr CeedElemRestriction 558 @param[out] blksize Variable to store size of blocks 559 560 @return An error code: 0 - success, otherwise - failure 561 562 @ref Advanced 563 **/ 564 int CeedElemRestrictionGetBlockSize(CeedElemRestriction rstr, 565 CeedInt *blksize) { 566 *blksize = rstr->blksize; 567 return 0; 568 } 569 570 /** 571 @brief Get the backend data of a CeedElemRestriction 572 573 @param rstr CeedElemRestriction 574 @param[out] data Variable to store data 575 576 @return An error code: 0 - success, otherwise - failure 577 578 @ref Advanced 579 **/ 580 int CeedElemRestrictionGetData(CeedElemRestriction rstr, void **data) { 581 *data = rstr->data; 582 return 0; 583 } 584 585 /** 586 @brief Set the backend data of a CeedElemRestriction 587 588 @param[out] rstr CeedElemRestriction 589 @param data Data to set 590 591 @return An error code: 0 - success, otherwise - failure 592 593 @ref Advanced 594 **/ 595 int CeedElemRestrictionSetData(CeedElemRestriction rstr, void **data) { 596 rstr->data = *data; 597 return 0; 598 } 599 600 /** 601 @brief View a CeedElemRestriction 602 603 @param[in] rstr CeedElemRestriction to view 604 @param[in] stream Stream to write; typically stdout/stderr or a file 605 606 @return Error code: 0 - success, otherwise - failure 607 608 @ref Utility 609 **/ 610 int CeedElemRestrictionView(CeedElemRestriction rstr, FILE *stream) { 611 fprintf(stream, "CeedElemRestriction from (%d, %d) to %d elements with %d " 612 "nodes each and L-vector components %s\n", rstr->nnodes, rstr->ncomp, 613 rstr->nelem, rstr->elemsize, CeedInterlaceModes[rstr->imode]); 614 return 0; 615 } 616 617 /** 618 @brief Destroy a CeedElemRestriction 619 620 @param rstr CeedElemRestriction to destroy 621 622 @return An error code: 0 - success, otherwise - failure 623 624 @ref Basic 625 **/ 626 int CeedElemRestrictionDestroy(CeedElemRestriction *rstr) { 627 int ierr; 628 629 if (!*rstr || --(*rstr)->refcount > 0) 630 return 0; 631 if ((*rstr)->Destroy) { 632 ierr = (*rstr)->Destroy(*rstr); CeedChk(ierr); 633 } 634 ierr = CeedDestroy(&(*rstr)->ceed); CeedChk(ierr); 635 ierr = CeedFree(rstr); CeedChk(ierr); 636 return 0; 637 } 638 639 /// @} 640