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 nelements. 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 nelements. 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 nelements. 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 size @a nnodes * @a ncomp when 298 tmode=CEED_NOTRANSPOSE) 299 @param v Output vector (of size @a nelem * @a elemsize when 300 tmode=CEED_NOTRANSPOSE) 301 @param request Request or CEED_REQUEST_IMMEDIATE 302 303 @return An error code: 0 - success, otherwise - failure 304 305 @ref Advanced 306 **/ 307 int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode tmode, 308 CeedTransposeMode lmode, CeedVector u, 309 CeedVector v, CeedRequest *request) { 310 CeedInt m,n; 311 int ierr; 312 313 if (tmode == CEED_NOTRANSPOSE) { 314 m = rstr->nblk * rstr->blksize * rstr->elemsize * rstr->ncomp; 315 n = rstr->nnodes * rstr->ncomp; 316 } else { 317 m = rstr->nnodes * rstr->ncomp; 318 n = rstr->nblk * rstr->blksize * rstr->elemsize * rstr->ncomp; 319 } 320 if (n != u->length) 321 // LCOV_EXCL_START 322 return CeedError(rstr->ceed, 2, "Input vector size %d not compatible with " 323 "element restriction (%d, %d)", u->length, m, n); 324 // LCOV_EXCL_STOP 325 if (m != v->length) 326 // LCOV_EXCL_START 327 return CeedError(rstr->ceed, 2, "Output vector size %d not compatible with " 328 "element restriction (%d, %d)", v->length, m, n); 329 // LCOV_EXCL_STOP 330 ierr = rstr->Apply(rstr, tmode, lmode, u, v, request); CeedChk(ierr); 331 332 return 0; 333 } 334 335 /** 336 @brief Restrict an L-vector to a block of an E-vector or apply its transpose 337 338 @param rstr CeedElemRestriction 339 @param block Block number to restrict to/from, i.e. block=0 will handle 340 elements [0 : blksize] and block=3 will handle elements 341 [3*blksize : 4*blksize] 342 @param tmode Apply restriction or transpose 343 @param lmode Ordering of the ncomp components, i.e. it specifies 344 the ordering of the components of the l-vector used 345 by this CeedElemRestriction. CEED_NOTRANSPOSE indicates 346 the component is the outermost index and CEED_TRANSPOSE 347 indicates the component is the innermost index in 348 ordering of the l-vector 349 tmode=CEED_NOTRANSPOSE) 350 @param u Input vector (of size @a nnodes * @a ncomp when 351 tmode=CEED_NOTRANSPOSE) 352 @param v Output vector (of size @a blksize * @a elemsize when 353 tmode=CEED_NOTRANSPOSE) 354 @param request Request or CEED_REQUEST_IMMEDIATE 355 356 @return An error code: 0 - success, otherwise - failure 357 358 @ref Advanced 359 **/ 360 int CeedElemRestrictionApplyBlock(CeedElemRestriction rstr, CeedInt block, 361 CeedTransposeMode tmode, 362 CeedTransposeMode lmode, CeedVector u, 363 CeedVector v, CeedRequest *request) { 364 CeedInt m,n; 365 int ierr; 366 367 if (tmode == CEED_NOTRANSPOSE) { 368 m = rstr->blksize * rstr->elemsize * rstr->ncomp; 369 n = rstr->nnodes * rstr->ncomp; 370 } else { 371 m = rstr->nnodes * rstr->ncomp; 372 n = rstr->blksize * rstr->elemsize * rstr->ncomp; 373 } 374 if (n != u->length) 375 // LCOV_EXCL_START 376 return CeedError(rstr->ceed, 2, "Input vector size %d not compatible with " 377 "element restriction (%d, %d)", u->length, m, n); 378 // LCOV_EXCL_STOP 379 if (m != v->length) 380 // LCOV_EXCL_START 381 return CeedError(rstr->ceed, 2, "Output vector size %d not compatible with " 382 "element restriction (%d, %d)", v->length, m, n); 383 // LCOV_EXCL_STOP 384 if (rstr->blksize*block > rstr->nelem) 385 // LCOV_EXCL_START 386 return CeedError(rstr->ceed, 2, "Cannot retrieve block %d, element %d > " 387 "total elements %d", block, rstr->blksize*block, 388 rstr->nelem); 389 // LCOV_EXCL_STOP 390 ierr = rstr->ApplyBlock(rstr, block, tmode, lmode, u, v, request); 391 CeedChk(ierr); 392 393 return 0; 394 } 395 396 /** 397 @brief Get the multiplicity of nodes in a CeedElemRestriction 398 399 @param rstr CeedElemRestriction 400 @param[out] mult Vector to store multiplicity (of size nnodes) 401 402 @return An error code: 0 - success, otherwise - failure 403 404 @ref Advanced 405 **/ 406 int CeedElemRestrictionGetMultiplicity(CeedElemRestriction rstr, 407 CeedVector mult) { 408 int ierr; 409 CeedVector evec; 410 411 // Create and set evec 412 ierr = CeedElemRestrictionCreateVector(rstr, NULL, &evec); CeedChk(ierr); 413 ierr = CeedVectorSetValue(evec, 1.0); CeedChk(ierr); 414 415 // Apply to get multiplicity 416 ierr = CeedElemRestrictionApply(rstr, CEED_TRANSPOSE, CEED_NOTRANSPOSE, evec, 417 mult, CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 418 419 // Cleanup 420 ierr = CeedVectorDestroy(&evec); CeedChk(ierr); 421 422 return 0; 423 } 424 425 /** 426 @brief Get the Ceed associated with a CeedElemRestriction 427 428 @param rstr CeedElemRestriction 429 @param[out] ceed Variable to store Ceed 430 431 @return An error code: 0 - success, otherwise - failure 432 433 @ref Advanced 434 **/ 435 int CeedElemRestrictionGetCeed(CeedElemRestriction rstr, Ceed *ceed) { 436 *ceed = rstr->ceed; 437 return 0; 438 } 439 440 /** 441 @brief Get the total number of elements in the range of a CeedElemRestriction 442 443 @param rstr CeedElemRestriction 444 @param[out] numelem Variable to store number of elements 445 446 @return An error code: 0 - success, otherwise - failure 447 448 @ref Advanced 449 **/ 450 int CeedElemRestrictionGetNumElements(CeedElemRestriction rstr, 451 CeedInt *numelem) { 452 *numelem = rstr->nelem; 453 return 0; 454 } 455 456 /** 457 @brief Get the size of elements in the CeedElemRestriction 458 459 @param rstr CeedElemRestriction 460 @param[out] elemsize Variable to store size of elements 461 462 @return An error code: 0 - success, otherwise - failure 463 464 @ref Advanced 465 **/ 466 int CeedElemRestrictionGetElementSize(CeedElemRestriction rstr, 467 CeedInt *elemsize) { 468 *elemsize = rstr->elemsize; 469 return 0; 470 } 471 472 /** 473 @brief Get the number of degrees of freedom in the range of a 474 CeedElemRestriction 475 476 @param rstr CeedElemRestriction 477 @param[out] numnodes Variable to store number of nodes 478 479 @return An error code: 0 - success, otherwise - failure 480 481 @ref Advanced 482 **/ 483 int CeedElemRestrictionGetNumNodes(CeedElemRestriction rstr, 484 CeedInt *numnodes) { 485 *numnodes = rstr->nnodes; 486 return 0; 487 } 488 489 /** 490 @brief Get the number of components in the elements of a 491 CeedElemRestriction 492 493 @param rstr CeedElemRestriction 494 @param[out] numcomp Variable to store number of components 495 496 @return An error code: 0 - success, otherwise - failure 497 498 @ref Advanced 499 **/ 500 int CeedElemRestrictionGetNumComponents(CeedElemRestriction rstr, 501 CeedInt *numcomp) { 502 *numcomp = rstr->ncomp; 503 return 0; 504 } 505 506 /** 507 @brief Get the number of blocks in a CeedElemRestriction 508 509 @param rstr CeedElemRestriction 510 @param[out] numblock Variable to store number of blocks 511 512 @return An error code: 0 - success, otherwise - failure 513 514 @ref Advanced 515 **/ 516 int CeedElemRestrictionGetNumBlocks(CeedElemRestriction rstr, 517 CeedInt *numblock) { 518 *numblock = rstr->nblk; 519 return 0; 520 } 521 522 /** 523 @brief Get the size of blocks in the CeedElemRestriction 524 525 @param rstr CeedElemRestriction 526 @param[out] blksize Variable to store size of blocks 527 528 @return An error code: 0 - success, otherwise - failure 529 530 @ref Advanced 531 **/ 532 int CeedElemRestrictionGetBlockSize(CeedElemRestriction rstr, 533 CeedInt *blksize) { 534 *blksize = rstr->blksize; 535 return 0; 536 } 537 538 /** 539 @brief Get the backend data of a CeedElemRestriction 540 541 @param rstr CeedElemRestriction 542 @param[out] data Variable to store data 543 544 @return An error code: 0 - success, otherwise - failure 545 546 @ref Advanced 547 **/ 548 int CeedElemRestrictionGetData(CeedElemRestriction rstr, void **data) { 549 *data = rstr->data; 550 return 0; 551 } 552 553 /** 554 @brief Set the backend data of a CeedElemRestriction 555 556 @param[out] rstr CeedElemRestriction 557 @param data Data to set 558 559 @return An error code: 0 - success, otherwise - failure 560 561 @ref Advanced 562 **/ 563 int CeedElemRestrictionSetData(CeedElemRestriction rstr, void **data) { 564 rstr->data = *data; 565 return 0; 566 } 567 568 /** 569 @brief View a CeedElemRestriction 570 571 @param[in] rstr CeedElemRestriction to view 572 @param[in] stream Stream to write; typically stdout/stderr or a file 573 574 @return Error code: 0 - success, otherwise - failure 575 576 @ref Utility 577 **/ 578 int CeedElemRestrictionView(CeedElemRestriction rstr, FILE *stream) { 579 fprintf(stream, "CeedElemRestriction from (%d, %d) to %d elements with %d " 580 "nodes each\n", rstr->nnodes, rstr->ncomp, rstr->nelem, 581 rstr->elemsize); 582 return 0; 583 } 584 585 /** 586 @brief Destroy a CeedElemRestriction 587 588 @param rstr CeedElemRestriction to destroy 589 590 @return An error code: 0 - success, otherwise - failure 591 592 @ref Basic 593 **/ 594 int CeedElemRestrictionDestroy(CeedElemRestriction *rstr) { 595 int ierr; 596 597 if (!*rstr || --(*rstr)->refcount > 0) 598 return 0; 599 if ((*rstr)->Destroy) { 600 ierr = (*rstr)->Destroy(*rstr); CeedChk(ierr); 601 } 602 ierr = CeedDestroy(&(*rstr)->ceed); CeedChk(ierr); 603 ierr = CeedFree(rstr); CeedChk(ierr); 604 return 0; 605 } 606 607 /// @} 608