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]. 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, 165 CeedInt blksize, 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); 252 CeedChk(ierr); 253 254 if (cmode == CEED_OWN_POINTER) { 255 ierr = CeedFree(&indices); CeedChk(ierr); 256 } 257 258 return 0; 259 } 260 261 /** 262 @brief Create CeedVectors associated with a CeedElemRestriction 263 264 @param rstr CeedElemRestriction 265 @param lvec The address of the L-vector to be created, or NULL 266 @param evec The address of the E-vector to be created, or NULL 267 268 @return An error code: 0 - success, otherwise - failure 269 270 @ref Advanced 271 **/ 272 int CeedElemRestrictionCreateVector(CeedElemRestriction rstr, CeedVector *lvec, 273 CeedVector *evec) { 274 int ierr; 275 CeedInt n, m; 276 m = rstr->nnodes * rstr->ncomp; 277 n = rstr->nblk * rstr->blksize * rstr->elemsize * rstr->ncomp; 278 if (lvec) { 279 ierr = CeedVectorCreate(rstr->ceed, m, lvec); CeedChk(ierr); 280 } 281 if (evec) { 282 ierr = CeedVectorCreate(rstr->ceed, n, evec); CeedChk(ierr); 283 } 284 return 0; 285 } 286 287 /** 288 @brief Restrict an L-vector to an E-vector or apply its transpose 289 290 @param rstr CeedElemRestriction 291 @param tmode Apply restriction or transpose 292 @param lmode Ordering of the ncomp components, i.e. it specifies 293 the ordering of the components of the l-vector used 294 by this CeedElemRestriction. CEED_NOTRANSPOSE indicates 295 the component is the outermost index and CEED_TRANSPOSE 296 indicates the component is the innermost index in 297 ordering of the l-vector 298 @param u Input vector (of size @a nnodes * @a ncomp when 299 tmode=CEED_NOTRANSPOSE) 300 @param v Output vector (of size @a nelem * @a elemsize when 301 tmode=CEED_NOTRANSPOSE) 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, 310 CeedVector u, 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 v Output vector (of size @a nelem * @a elemsize when 352 tmode=CEED_NOTRANSPOSE) 353 @param request Request or CEED_REQUEST_IMMEDIATE 354 355 @return An error code: 0 - success, otherwise - failure 356 357 @ref Advanced 358 **/ 359 int CeedElemRestrictionApplyBlock(CeedElemRestriction rstr, CeedInt block, 360 CeedTransposeMode tmode, 361 CeedTransposeMode lmode, 362 CeedVector u, CeedVector v, 363 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