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