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[out] mult Vector to store multiplicity (of size nnodes) 405 406 @return An error code: 0 - success, otherwise - failure 407 408 @ref Advanced 409 **/ 410 int CeedElemRestrictionGetMultiplicity(CeedElemRestriction rstr, 411 CeedVector mult) { 412 int ierr; 413 CeedVector evec; 414 415 // Create and set evec 416 ierr = CeedElemRestrictionCreateVector(rstr, NULL, &evec); CeedChk(ierr); 417 ierr = CeedVectorSetValue(evec, 1.0); CeedChk(ierr); 418 419 // Apply to get multiplicity 420 ierr = CeedElemRestrictionApply(rstr, CEED_TRANSPOSE, CEED_NOTRANSPOSE, evec, 421 mult, CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 422 423 // Cleanup 424 ierr = CeedVectorDestroy(&evec); CeedChk(ierr); 425 426 return 0; 427 } 428 429 /** 430 @brief Get the Ceed associated with a CeedElemRestriction 431 432 @param rstr CeedElemRestriction 433 @param[out] ceed Variable to store Ceed 434 435 @return An error code: 0 - success, otherwise - failure 436 437 @ref Advanced 438 **/ 439 int CeedElemRestrictionGetCeed(CeedElemRestriction rstr, Ceed *ceed) { 440 *ceed = rstr->ceed; 441 return 0; 442 } 443 444 /** 445 @brief Get the total number of elements in the range of a CeedElemRestriction 446 447 @param rstr CeedElemRestriction 448 @param[out] numelem Variable to store number of elements 449 450 @return An error code: 0 - success, otherwise - failure 451 452 @ref Advanced 453 **/ 454 int CeedElemRestrictionGetNumElements(CeedElemRestriction rstr, 455 CeedInt *numelem) { 456 *numelem = rstr->nelem; 457 return 0; 458 } 459 460 /** 461 @brief Get the size of elements in the CeedElemRestriction 462 463 @param rstr CeedElemRestriction 464 @param[out] elemsize Variable to store size of elements 465 466 @return An error code: 0 - success, otherwise - failure 467 468 @ref Advanced 469 **/ 470 int CeedElemRestrictionGetElementSize(CeedElemRestriction rstr, 471 CeedInt *elemsize) { 472 *elemsize = rstr->elemsize; 473 return 0; 474 } 475 476 /** 477 @brief Get the number of degrees of freedom in the range of a 478 CeedElemRestriction 479 480 @param rstr CeedElemRestriction 481 @param[out] numnodes Variable to store number of nodes 482 483 @return An error code: 0 - success, otherwise - failure 484 485 @ref Advanced 486 **/ 487 int CeedElemRestrictionGetNumNodes(CeedElemRestriction rstr, 488 CeedInt *numnodes) { 489 *numnodes = rstr->nnodes; 490 return 0; 491 } 492 493 /** 494 @brief Get the number of components in the elements of a 495 CeedElemRestriction 496 497 @param rstr CeedElemRestriction 498 @param[out] numcomp Variable to store number of components 499 500 @return An error code: 0 - success, otherwise - failure 501 502 @ref Advanced 503 **/ 504 int CeedElemRestrictionGetNumComponents(CeedElemRestriction rstr, 505 CeedInt *numcomp) { 506 *numcomp = rstr->ncomp; 507 return 0; 508 } 509 510 /** 511 @brief Get the number of blocks in a CeedElemRestriction 512 513 @param rstr CeedElemRestriction 514 @param[out] numblock Variable to store number of blocks 515 516 @return An error code: 0 - success, otherwise - failure 517 518 @ref Advanced 519 **/ 520 int CeedElemRestrictionGetNumBlocks(CeedElemRestriction rstr, 521 CeedInt *numblock) { 522 *numblock = rstr->nblk; 523 return 0; 524 } 525 526 /** 527 @brief Get the size of blocks in the CeedElemRestriction 528 529 @param rstr CeedElemRestriction 530 @param[out] blksize Variable to store size of blocks 531 532 @return An error code: 0 - success, otherwise - failure 533 534 @ref Advanced 535 **/ 536 int CeedElemRestrictionGetBlockSize(CeedElemRestriction rstr, 537 CeedInt *blksize) { 538 *blksize = rstr->blksize; 539 return 0; 540 } 541 542 /** 543 @brief Get the backend data of a CeedElemRestriction 544 545 @param rstr CeedElemRestriction 546 @param[out] data Variable to store data 547 548 @return An error code: 0 - success, otherwise - failure 549 550 @ref Advanced 551 **/ 552 int CeedElemRestrictionGetData(CeedElemRestriction rstr, void **data) { 553 *data = rstr->data; 554 return 0; 555 } 556 557 /** 558 @brief Set the backend data of a CeedElemRestriction 559 560 @param[out] rstr CeedElemRestriction 561 @param data Data to set 562 563 @return An error code: 0 - success, otherwise - failure 564 565 @ref Advanced 566 **/ 567 int CeedElemRestrictionSetData(CeedElemRestriction rstr, void **data) { 568 rstr->data = *data; 569 return 0; 570 } 571 572 /** 573 @brief View a CeedElemRestriction 574 575 @param[in] rstr CeedElemRestriction to view 576 @param[in] stream Stream to write; typically stdout/stderr or a file 577 578 @return Error code: 0 - success, otherwise - failure 579 580 @ref Utility 581 **/ 582 int CeedElemRestrictionView(CeedElemRestriction rstr, FILE *stream) { 583 fprintf(stream, "CeedElemRestriction from (%d, %d) to %d elements with %d " 584 "nodes each\n", rstr->nnodes, rstr->ncomp, rstr->nelem, 585 rstr->elemsize); 586 return 0; 587 } 588 589 /** 590 @brief Destroy a CeedElemRestriction 591 592 @param rstr CeedElemRestriction to destroy 593 594 @return An error code: 0 - success, otherwise - failure 595 596 @ref Basic 597 **/ 598 int CeedElemRestrictionDestroy(CeedElemRestriction *rstr) { 599 int ierr; 600 601 if (!*rstr || --(*rstr)->refcount > 0) 602 return 0; 603 if ((*rstr)->Destroy) { 604 ierr = (*rstr)->Destroy(*rstr); CeedChk(ierr); 605 } 606 ierr = CeedDestroy(&(*rstr)->ceed); CeedChk(ierr); 607 ierr = CeedFree(rstr); CeedChk(ierr); 608 return 0; 609 } 610 611 /// @} 612