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