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