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