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/ceed.h> 18 #include <ceed/backend.h> 19 #include <ceed-impl.h> 20 #include <stdbool.h> 21 #include <stdio.h> 22 23 /// @file 24 /// Implementation of CeedElemRestriction interfaces 25 26 /// ---------------------------------------------------------------------------- 27 /// CeedElemRestriction Library Internal Functions 28 /// ---------------------------------------------------------------------------- 29 /// @addtogroup CeedElemRestrictionDeveloper 30 /// @{ 31 32 /** 33 @brief Permute and pad offsets for a blocked restriction 34 35 @param offsets Array of shape [@a num_elem, @a elem_size]. Row i holds the 36 ordered list of the offsets (into the input CeedVector) 37 for the unknowns corresponding to element i, where 38 0 <= i < @a num_elem. All offsets must be in the range 39 [0, @a l_size - 1]. 40 @param blk_offsets Array of permuted and padded offsets of 41 shape [@a num_blk, @a elem_size, @a blk_size]. 42 @param num_blk Number of blocks 43 @param num_elem Number of elements 44 @param blk_size Number of elements in a block 45 @param elem_size Size of each element 46 47 @return An error code: 0 - success, otherwise - failure 48 49 @ref Utility 50 **/ 51 int CeedPermutePadOffsets(const CeedInt *offsets, CeedInt *blk_offsets, 52 CeedInt num_blk, CeedInt num_elem, CeedInt blk_size, 53 CeedInt elem_size) { 54 for (CeedInt e=0; e<num_blk*blk_size; e+=blk_size) 55 for (int j=0; j<blk_size; j++) 56 for (int k=0; k<elem_size; k++) 57 blk_offsets[e*elem_size + k*blk_size + j] 58 = offsets[CeedIntMin(e+j,num_elem-1)*elem_size + k]; 59 return CEED_ERROR_SUCCESS; 60 } 61 62 /// @} 63 64 /// ---------------------------------------------------------------------------- 65 /// CeedElemRestriction Backend API 66 /// ---------------------------------------------------------------------------- 67 /// @addtogroup CeedElemRestrictionBackend 68 /// @{ 69 70 /** 71 @brief Get the Ceed associated with a CeedElemRestriction 72 73 @param rstr CeedElemRestriction 74 @param[out] ceed Variable to store Ceed 75 76 @return An error code: 0 - success, otherwise - failure 77 78 @ref Backend 79 **/ 80 int CeedElemRestrictionGetCeed(CeedElemRestriction rstr, Ceed *ceed) { 81 *ceed = rstr->ceed; 82 return CEED_ERROR_SUCCESS; 83 } 84 85 /** 86 87 @brief Get the strides of a strided CeedElemRestriction 88 89 @param rstr CeedElemRestriction 90 @param[out] strides Variable to store strides array 91 92 @return An error code: 0 - success, otherwise - failure 93 94 @ref Backend 95 **/ 96 int CeedElemRestrictionGetStrides(CeedElemRestriction rstr, 97 CeedInt (*strides)[3]) { 98 if (!rstr->strides) 99 // LCOV_EXCL_START 100 return CeedError(rstr->ceed, CEED_ERROR_MINOR, 101 "ElemRestriction has no stride data"); 102 // LCOV_EXCL_STOP 103 104 for (int i=0; i<3; i++) 105 (*strides)[i] = rstr->strides[i]; 106 return CEED_ERROR_SUCCESS; 107 } 108 109 /** 110 @brief Get read-only access to a CeedElemRestriction offsets array by memtype 111 112 @param rstr CeedElemRestriction to retrieve offsets 113 @param mem_type Memory type on which to access the array. If the backend 114 uses a different memory type, this will perform a copy 115 (possibly cached). 116 @param[out] offsets Array on memory type mem_type 117 118 @return An error code: 0 - success, otherwise - failure 119 120 @ref User 121 **/ 122 int CeedElemRestrictionGetOffsets(CeedElemRestriction rstr, 123 CeedMemType mem_type, 124 const CeedInt **offsets) { 125 int ierr; 126 127 if (!rstr->GetOffsets) 128 // LCOV_EXCL_START 129 return CeedError(rstr->ceed, CEED_ERROR_UNSUPPORTED, 130 "Backend does not support GetOffsets"); 131 // LCOV_EXCL_STOP 132 133 ierr = rstr->GetOffsets(rstr, mem_type, offsets); CeedChk(ierr); 134 rstr->num_readers++; 135 return CEED_ERROR_SUCCESS; 136 } 137 138 /** 139 @brief Restore an offsets array obtained using CeedElemRestrictionGetOffsets() 140 141 @param rstr CeedElemRestriction to restore 142 @param offsets Array of offset data 143 144 @return An error code: 0 - success, otherwise - failure 145 146 @ref User 147 **/ 148 int CeedElemRestrictionRestoreOffsets(CeedElemRestriction rstr, 149 const CeedInt **offsets) { 150 *offsets = NULL; 151 rstr->num_readers--; 152 return CEED_ERROR_SUCCESS; 153 } 154 155 /** 156 @brief Get the strided status of a CeedElemRestriction 157 158 @param rstr CeedElemRestriction 159 @param[out] is_strided Variable to store strided status, 1 if strided else 0 160 161 @return An error code: 0 - success, otherwise - failure 162 163 @ref Backend 164 **/ 165 int CeedElemRestrictionIsStrided(CeedElemRestriction rstr, bool *is_strided) { 166 *is_strided = rstr->strides ? true : false; 167 return CEED_ERROR_SUCCESS; 168 } 169 170 /** 171 @brief Get the backend stride status of a CeedElemRestriction 172 173 @param rstr CeedElemRestriction 174 @param[out] has_backend_strides Variable to store stride status 175 176 @return An error code: 0 - success, otherwise - failure 177 178 @ref Backend 179 **/ 180 int CeedElemRestrictionHasBackendStrides(CeedElemRestriction rstr, 181 bool *has_backend_strides) { 182 if (!rstr->strides) 183 // LCOV_EXCL_START 184 return CeedError(rstr->ceed, CEED_ERROR_MINOR, 185 "ElemRestriction has no stride data"); 186 // LCOV_EXCL_STOP 187 188 *has_backend_strides = ((rstr->strides[0] == CEED_STRIDES_BACKEND[0]) && 189 (rstr->strides[1] == CEED_STRIDES_BACKEND[1]) && 190 (rstr->strides[2] == CEED_STRIDES_BACKEND[2])); 191 return CEED_ERROR_SUCCESS; 192 } 193 194 /** 195 196 @brief Get the E-vector layout of a CeedElemRestriction 197 198 @param rstr CeedElemRestriction 199 @param[out] layout Variable to store layout array, 200 stored as [nodes, components, elements]. 201 The data for node i, component j, element k in the 202 E-vector is given by 203 i*layout[0] + j*layout[1] + k*layout[2] 204 205 @return An error code: 0 - success, otherwise - failure 206 207 @ref Backend 208 **/ 209 int CeedElemRestrictionGetELayout(CeedElemRestriction rstr, 210 CeedInt (*layout)[3]) { 211 if (!rstr->layout[0]) 212 // LCOV_EXCL_START 213 return CeedError(rstr->ceed, CEED_ERROR_MINOR, 214 "ElemRestriction has no layout data"); 215 // LCOV_EXCL_STOP 216 217 for (int i=0; i<3; i++) 218 (*layout)[i] = rstr->layout[i]; 219 return CEED_ERROR_SUCCESS; 220 } 221 222 /** 223 224 @brief Set the E-vector layout of a CeedElemRestriction 225 226 @param rstr CeedElemRestriction 227 @param layout Variable to containing layout array, 228 stored as [nodes, components, elements]. 229 The data for node i, component j, element k in the 230 E-vector is given by 231 i*layout[0] + j*layout[1] + k*layout[2] 232 233 @return An error code: 0 - success, otherwise - failure 234 235 @ref Backend 236 **/ 237 int CeedElemRestrictionSetELayout(CeedElemRestriction rstr, 238 CeedInt layout[3]) { 239 for (int i = 0; i<3; i++) 240 rstr->layout[i] = layout[i]; 241 return CEED_ERROR_SUCCESS; 242 } 243 244 /** 245 @brief Get the backend data of a CeedElemRestriction 246 247 @param rstr CeedElemRestriction 248 @param[out] data Variable to store data 249 250 @return An error code: 0 - success, otherwise - failure 251 252 @ref Backend 253 **/ 254 int CeedElemRestrictionGetData(CeedElemRestriction rstr, void *data) { 255 *(void **)data = rstr->data; 256 return CEED_ERROR_SUCCESS; 257 } 258 259 /** 260 @brief Set the backend data of a CeedElemRestriction 261 262 @param[out] rstr CeedElemRestriction 263 @param data Data to set 264 265 @return An error code: 0 - success, otherwise - failure 266 267 @ref Backend 268 **/ 269 int CeedElemRestrictionSetData(CeedElemRestriction rstr, void *data) { 270 rstr->data = data; 271 return CEED_ERROR_SUCCESS; 272 } 273 274 /** 275 @brief Increment the reference counter for a CeedElemRestriction 276 277 @param rstr ElemRestriction to increment the reference counter 278 279 @return An error code: 0 - success, otherwise - failure 280 281 @ref Backend 282 **/ 283 int CeedElemRestrictionReference(CeedElemRestriction rstr) { 284 rstr->ref_count++; 285 return CEED_ERROR_SUCCESS; 286 } 287 288 /// @} 289 290 /// @cond DOXYGEN_SKIP 291 static struct CeedElemRestriction_private ceed_elemrestriction_none; 292 /// @endcond 293 294 /// ---------------------------------------------------------------------------- 295 /// CeedElemRestriction Public API 296 /// ---------------------------------------------------------------------------- 297 /// @addtogroup CeedElemRestrictionUser 298 /// @{ 299 300 /// Indicate that the stride is determined by the backend 301 const CeedInt CEED_STRIDES_BACKEND[3] = {}; 302 303 /// Indicate that no CeedElemRestriction is provided by the user 304 const CeedElemRestriction CEED_ELEMRESTRICTION_NONE = 305 &ceed_elemrestriction_none; 306 307 /** 308 @brief Create a CeedElemRestriction 309 310 @param ceed A Ceed object where the CeedElemRestriction will be created 311 @param num_elem Number of elements described in the @a offsets array 312 @param elem_size Size (number of "nodes") per element 313 @param num_comp Number of field components per interpolation node 314 (1 for scalar fields) 315 @param comp_stride Stride between components for the same L-vector "node". 316 Data for node i, component j, element k can be found in 317 the L-vector at index 318 offsets[i + k*elem_size] + j*comp_stride. 319 @param l_size The size of the L-vector. This vector may be larger than 320 the elements and fields given by this restriction. 321 @param mem_type Memory type of the @a offsets array, see CeedMemType 322 @param copy_mode Copy mode for the @a offsets array, see CeedCopyMode 323 @param offsets Array of shape [@a num_elem, @a elem_size]. Row i holds the 324 ordered list of the offsets (into the input CeedVector) 325 for the unknowns corresponding to element i, where 326 0 <= i < @a num_elem. All offsets must be in the range 327 [0, @a l_size - 1]. 328 @param[out] rstr Address of the variable where the newly created 329 CeedElemRestriction will be stored 330 331 @return An error code: 0 - success, otherwise - failure 332 333 @ref User 334 **/ 335 int CeedElemRestrictionCreate(Ceed ceed, CeedInt num_elem, CeedInt elem_size, 336 CeedInt num_comp, CeedInt comp_stride, 337 CeedInt l_size, CeedMemType mem_type, 338 CeedCopyMode copy_mode, const CeedInt *offsets, 339 CeedElemRestriction *rstr) { 340 int ierr; 341 342 if (!ceed->ElemRestrictionCreate) { 343 Ceed delegate; 344 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 345 CeedChk(ierr); 346 347 if (!delegate) 348 // LCOV_EXCL_START 349 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 350 "Backend does not support ElemRestrictionCreate"); 351 // LCOV_EXCL_STOP 352 353 ierr = CeedElemRestrictionCreate(delegate, num_elem, elem_size, num_comp, 354 comp_stride, l_size, mem_type, copy_mode, 355 offsets, rstr); CeedChk(ierr); 356 return CEED_ERROR_SUCCESS; 357 } 358 359 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 360 (*rstr)->ceed = ceed; 361 ierr = CeedReference(ceed); CeedChk(ierr); 362 (*rstr)->ref_count = 1; 363 (*rstr)->num_elem = num_elem; 364 (*rstr)->elem_size = elem_size; 365 (*rstr)->num_comp = num_comp; 366 (*rstr)->comp_stride = comp_stride; 367 (*rstr)->l_size = l_size; 368 (*rstr)->num_blk = num_elem; 369 (*rstr)->blk_size = 1; 370 ierr = ceed->ElemRestrictionCreate(mem_type, copy_mode, offsets, *rstr); 371 CeedChk(ierr); 372 return CEED_ERROR_SUCCESS; 373 } 374 375 /** 376 @brief Create a strided CeedElemRestriction 377 378 @param ceed A Ceed object where the CeedElemRestriction will be created 379 @param num_elem Number of elements described by the restriction 380 @param elem_size Size (number of "nodes") per element 381 @param num_comp Number of field components per interpolation "node" 382 (1 for scalar fields) 383 @param l_size The size of the L-vector. This vector may be larger than 384 the elements and fields given by this restriction. 385 @param strides Array for strides between [nodes, components, elements]. 386 Data for node i, component j, element k can be found in 387 the L-vector at index 388 i*strides[0] + j*strides[1] + k*strides[2]. 389 @a CEED_STRIDES_BACKEND may be used with vectors created 390 by a Ceed backend. 391 @param rstr Address of the variable where the newly created 392 CeedElemRestriction will be stored 393 394 @return An error code: 0 - success, otherwise - failure 395 396 @ref User 397 **/ 398 int CeedElemRestrictionCreateStrided(Ceed ceed, CeedInt num_elem, 399 CeedInt elem_size, 400 CeedInt num_comp, CeedInt l_size, 401 const CeedInt strides[3], 402 CeedElemRestriction *rstr) { 403 int ierr; 404 405 if (!ceed->ElemRestrictionCreate) { 406 Ceed delegate; 407 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 408 CeedChk(ierr); 409 410 if (!delegate) 411 // LCOV_EXCL_START 412 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 413 "Backend does not support ElemRestrictionCreate"); 414 // LCOV_EXCL_STOP 415 416 ierr = CeedElemRestrictionCreateStrided(delegate, num_elem, elem_size, num_comp, 417 l_size, strides, rstr); 418 CeedChk(ierr); 419 return CEED_ERROR_SUCCESS; 420 } 421 422 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 423 (*rstr)->ceed = ceed; 424 ierr = CeedReference(ceed); CeedChk(ierr); 425 (*rstr)->ref_count = 1; 426 (*rstr)->num_elem = num_elem; 427 (*rstr)->elem_size = elem_size; 428 (*rstr)->num_comp = num_comp; 429 (*rstr)->l_size = l_size; 430 (*rstr)->num_blk = num_elem; 431 (*rstr)->blk_size = 1; 432 ierr = CeedMalloc(3, &(*rstr)->strides); CeedChk(ierr); 433 for (int i=0; i<3; i++) 434 (*rstr)->strides[i] = strides[i]; 435 ierr = ceed->ElemRestrictionCreate(CEED_MEM_HOST, CEED_OWN_POINTER, NULL, 436 *rstr); 437 CeedChk(ierr); 438 return CEED_ERROR_SUCCESS; 439 } 440 441 /** 442 @brief Create a blocked CeedElemRestriction, typically only called by backends 443 444 @param ceed A Ceed object where the CeedElemRestriction will be created. 445 @param num_elem Number of elements described in the @a offsets array. 446 @param elem_size Size (number of unknowns) per element 447 @param blk_size Number of elements in a block 448 @param num_comp Number of field components per interpolation node 449 (1 for scalar fields) 450 @param comp_stride Stride between components for the same L-vector "node". 451 Data for node i, component j, element k can be found in 452 the L-vector at index 453 offsets[i + k*elem_size] + j*comp_stride. 454 @param l_size The size of the L-vector. This vector may be larger than 455 the elements and fields given by this restriction. 456 @param mem_type Memory type of the @a offsets array, see CeedMemType 457 @param copy_mode Copy mode for the @a offsets array, see CeedCopyMode 458 @param offsets Array of shape [@a num_elem, @a elem_size]. Row i holds the 459 ordered list of the offsets (into the input CeedVector) 460 for the unknowns corresponding to element i, where 461 0 <= i < @a num_elem. All offsets must be in the range 462 [0, @a l_size - 1]. The backend will permute and pad this 463 array to the desired ordering for the blocksize, which is 464 typically given by the backend. The default reordering is 465 to interlace elements. 466 @param rstr Address of the variable where the newly created 467 CeedElemRestriction will be stored 468 469 @return An error code: 0 - success, otherwise - failure 470 471 @ref Backend 472 **/ 473 int CeedElemRestrictionCreateBlocked(Ceed ceed, CeedInt num_elem, 474 CeedInt elem_size, 475 CeedInt blk_size, CeedInt num_comp, 476 CeedInt comp_stride, CeedInt l_size, 477 CeedMemType mem_type, CeedCopyMode copy_mode, 478 const CeedInt *offsets, 479 CeedElemRestriction *rstr) { 480 int ierr; 481 CeedInt *blk_offsets; 482 CeedInt num_blk = (num_elem / blk_size) + !!(num_elem % blk_size); 483 484 if (!ceed->ElemRestrictionCreateBlocked) { 485 Ceed delegate; 486 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 487 CeedChk(ierr); 488 489 if (!delegate) 490 // LCOV_EXCL_START 491 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support " 492 "ElemRestrictionCreateBlocked"); 493 // LCOV_EXCL_STOP 494 495 ierr = CeedElemRestrictionCreateBlocked(delegate, num_elem, elem_size, blk_size, 496 num_comp, comp_stride, l_size, mem_type, 497 copy_mode, offsets, rstr); 498 CeedChk(ierr); 499 return CEED_ERROR_SUCCESS; 500 } 501 502 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 503 504 ierr = CeedCalloc(num_blk*blk_size*elem_size, &blk_offsets); CeedChk(ierr); 505 ierr = CeedPermutePadOffsets(offsets, blk_offsets, num_blk, num_elem, blk_size, 506 elem_size); CeedChk(ierr); 507 508 (*rstr)->ceed = ceed; 509 ierr = CeedReference(ceed); CeedChk(ierr); 510 (*rstr)->ref_count = 1; 511 (*rstr)->num_elem = num_elem; 512 (*rstr)->elem_size = elem_size; 513 (*rstr)->num_comp = num_comp; 514 (*rstr)->comp_stride = comp_stride; 515 (*rstr)->l_size = l_size; 516 (*rstr)->num_blk = num_blk; 517 (*rstr)->blk_size = blk_size; 518 ierr = ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, 519 (const CeedInt *) blk_offsets, *rstr); CeedChk(ierr); 520 if (copy_mode == CEED_OWN_POINTER) { 521 ierr = CeedFree(&offsets); CeedChk(ierr); 522 } 523 return CEED_ERROR_SUCCESS; 524 } 525 526 /** 527 @brief Create a blocked strided CeedElemRestriction 528 529 @param ceed A Ceed object where the CeedElemRestriction will be created 530 @param num_elem Number of elements described by the restriction 531 @param elem_size Size (number of "nodes") per element 532 @param blk_size Number of elements in a block 533 @param num_comp Number of field components per interpolation node 534 (1 for scalar fields) 535 @param l_size The size of the L-vector. This vector may be larger than 536 the elements and fields given by this restriction. 537 @param strides Array for strides between [nodes, components, elements]. 538 Data for node i, component j, element k can be found in 539 the L-vector at index 540 i*strides[0] + j*strides[1] + k*strides[2]. 541 @a CEED_STRIDES_BACKEND may be used with vectors created 542 by a Ceed backend. 543 @param rstr Address of the variable where the newly created 544 CeedElemRestriction will be stored 545 546 @return An error code: 0 - success, otherwise - failure 547 548 @ref User 549 **/ 550 int CeedElemRestrictionCreateBlockedStrided(Ceed ceed, CeedInt num_elem, 551 CeedInt elem_size, CeedInt blk_size, CeedInt num_comp, CeedInt l_size, 552 const CeedInt strides[3], CeedElemRestriction *rstr) { 553 int ierr; 554 CeedInt num_blk = (num_elem / blk_size) + !!(num_elem % blk_size); 555 556 if (!ceed->ElemRestrictionCreateBlocked) { 557 Ceed delegate; 558 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 559 CeedChk(ierr); 560 561 if (!delegate) 562 // LCOV_EXCL_START 563 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support " 564 "ElemRestrictionCreateBlocked"); 565 // LCOV_EXCL_STOP 566 567 ierr = CeedElemRestrictionCreateBlockedStrided(delegate, num_elem, elem_size, 568 blk_size, num_comp, l_size, strides, rstr); CeedChk(ierr); 569 return CEED_ERROR_SUCCESS; 570 } 571 572 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 573 574 (*rstr)->ceed = ceed; 575 ierr = CeedReference(ceed); CeedChk(ierr); 576 (*rstr)->ref_count = 1; 577 (*rstr)->num_elem = num_elem; 578 (*rstr)->elem_size = elem_size; 579 (*rstr)->num_comp = num_comp; 580 (*rstr)->l_size = l_size; 581 (*rstr)->num_blk = num_blk; 582 (*rstr)->blk_size = blk_size; 583 ierr = CeedMalloc(3, &(*rstr)->strides); CeedChk(ierr); 584 for (int i=0; i<3; i++) 585 (*rstr)->strides[i] = strides[i]; 586 ierr = ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, 587 NULL, *rstr); CeedChk(ierr); 588 return CEED_ERROR_SUCCESS; 589 } 590 591 /** 592 @brief Copy the pointer to a CeedElemRestriction. Both pointers should 593 be destroyed with `CeedElemRestrictionDestroy()`; 594 Note: If `*rstr_copy` is non-NULL, then it is assumed that 595 `*rstr_copy` is a pointer to a CeedElemRestriction. This 596 CeedElemRestriction will be destroyed if `*rstr_copy` is the 597 only reference to this CeedElemRestriction. 598 599 @param rstr CeedElemRestriction to copy reference to 600 @param[out] rstr_copy Variable to store copied reference 601 602 @return An error code: 0 - success, otherwise - failure 603 604 @ref User 605 **/ 606 int CeedElemRestrictionReferenceCopy(CeedElemRestriction rstr, 607 CeedElemRestriction *rstr_copy) { 608 int ierr; 609 610 ierr = CeedElemRestrictionReference(rstr); CeedChk(ierr); 611 ierr = CeedElemRestrictionDestroy(rstr_copy); CeedChk(ierr); 612 *rstr_copy = rstr; 613 return CEED_ERROR_SUCCESS; 614 } 615 616 /** 617 @brief Create CeedVectors associated with a CeedElemRestriction 618 619 @param rstr CeedElemRestriction 620 @param l_vec The address of the L-vector to be created, or NULL 621 @param e_vec The address of the E-vector to be created, or NULL 622 623 @return An error code: 0 - success, otherwise - failure 624 625 @ref User 626 **/ 627 int CeedElemRestrictionCreateVector(CeedElemRestriction rstr, CeedVector *l_vec, 628 CeedVector *e_vec) { 629 int ierr; 630 CeedInt e_size, l_size; 631 l_size = rstr->l_size; 632 e_size = rstr->num_blk * rstr->blk_size * rstr->elem_size * rstr->num_comp; 633 if (l_vec) { 634 ierr = CeedVectorCreate(rstr->ceed, l_size, l_vec); CeedChk(ierr); 635 } 636 if (e_vec) { 637 ierr = CeedVectorCreate(rstr->ceed, e_size, e_vec); CeedChk(ierr); 638 } 639 return CEED_ERROR_SUCCESS; 640 } 641 642 /** 643 @brief Restrict an L-vector to an E-vector or apply its transpose 644 645 @param rstr CeedElemRestriction 646 @param t_mode Apply restriction or transpose 647 @param u Input vector (of size @a l_size when t_mode=@ref CEED_NOTRANSPOSE) 648 @param ru Output vector (of shape [@a num_elem * @a elem_size] when 649 t_mode=@ref CEED_NOTRANSPOSE). Ordering of the e-vector is decided 650 by the backend. 651 @param request Request or @ref CEED_REQUEST_IMMEDIATE 652 653 @return An error code: 0 - success, otherwise - failure 654 655 @ref User 656 **/ 657 int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode t_mode, 658 CeedVector u, CeedVector ru, 659 CeedRequest *request) { 660 CeedInt m, n; 661 int ierr; 662 663 if (t_mode == CEED_NOTRANSPOSE) { 664 m = rstr->num_blk * rstr->blk_size * rstr->elem_size * rstr->num_comp; 665 n = rstr->l_size; 666 } else { 667 m = rstr->l_size; 668 n = rstr->num_blk * rstr->blk_size * rstr->elem_size * rstr->num_comp; 669 } 670 if (n != u->length) 671 // LCOV_EXCL_START 672 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 673 "Input vector size %d not compatible with " 674 "element restriction (%d, %d)", u->length, m, n); 675 // LCOV_EXCL_STOP 676 if (m != ru->length) 677 // LCOV_EXCL_START 678 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 679 "Output vector size %d not compatible with " 680 "element restriction (%d, %d)", ru->length, m, n); 681 // LCOV_EXCL_STOP 682 ierr = rstr->Apply(rstr, t_mode, u, ru, request); CeedChk(ierr); 683 return CEED_ERROR_SUCCESS; 684 } 685 686 /** 687 @brief Restrict an L-vector to a block of an E-vector or apply its transpose 688 689 @param rstr CeedElemRestriction 690 @param block Block number to restrict to/from, i.e. block=0 will handle 691 elements [0 : blk_size] and block=3 will handle elements 692 [3*blk_size : 4*blk_size] 693 @param t_mode Apply restriction or transpose 694 @param u Input vector (of size @a l_size when t_mode=@ref CEED_NOTRANSPOSE) 695 @param ru Output vector (of shape [@a blk_size * @a elem_size] when 696 t_mode=@ref CEED_NOTRANSPOSE). Ordering of the e-vector is decided 697 by the backend. 698 @param request Request or @ref CEED_REQUEST_IMMEDIATE 699 700 @return An error code: 0 - success, otherwise - failure 701 702 @ref Backend 703 **/ 704 int CeedElemRestrictionApplyBlock(CeedElemRestriction rstr, CeedInt block, 705 CeedTransposeMode t_mode, CeedVector u, 706 CeedVector ru, CeedRequest *request) { 707 CeedInt m, n; 708 int ierr; 709 710 if (t_mode == CEED_NOTRANSPOSE) { 711 m = rstr->blk_size * rstr->elem_size * rstr->num_comp; 712 n = rstr->l_size; 713 } else { 714 m = rstr->l_size; 715 n = rstr->blk_size * rstr->elem_size * rstr->num_comp; 716 } 717 if (n != u->length) 718 // LCOV_EXCL_START 719 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 720 "Input vector size %d not compatible with " 721 "element restriction (%d, %d)", u->length, m, n); 722 // LCOV_EXCL_STOP 723 if (m != ru->length) 724 // LCOV_EXCL_START 725 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 726 "Output vector size %d not compatible with " 727 "element restriction (%d, %d)", ru->length, m, n); 728 // LCOV_EXCL_STOP 729 if (rstr->blk_size*block > rstr->num_elem) 730 // LCOV_EXCL_START 731 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 732 "Cannot retrieve block %d, element %d > " 733 "total elements %d", block, rstr->blk_size*block, 734 rstr->num_elem); 735 // LCOV_EXCL_STOP 736 ierr = rstr->ApplyBlock(rstr, block, t_mode, u, ru, request); 737 CeedChk(ierr); 738 return CEED_ERROR_SUCCESS; 739 } 740 741 /** 742 @brief Get the L-vector component stride 743 744 @param rstr CeedElemRestriction 745 @param[out] comp_stride Variable to store component stride 746 747 @return An error code: 0 - success, otherwise - failure 748 749 @ref Backend 750 **/ 751 int CeedElemRestrictionGetCompStride(CeedElemRestriction rstr, 752 CeedInt *comp_stride) { 753 *comp_stride = rstr->comp_stride; 754 return CEED_ERROR_SUCCESS; 755 } 756 757 /** 758 @brief Get the total number of elements in the range of a CeedElemRestriction 759 760 @param rstr CeedElemRestriction 761 @param[out] num_elem Variable to store number of elements 762 763 @return An error code: 0 - success, otherwise - failure 764 765 @ref Backend 766 **/ 767 int CeedElemRestrictionGetNumElements(CeedElemRestriction rstr, 768 CeedInt *num_elem) { 769 *num_elem = rstr->num_elem; 770 return CEED_ERROR_SUCCESS; 771 } 772 773 /** 774 @brief Get the size of elements in the CeedElemRestriction 775 776 @param rstr CeedElemRestriction 777 @param[out] elem_size Variable to store size of elements 778 779 @return An error code: 0 - success, otherwise - failure 780 781 @ref Backend 782 **/ 783 int CeedElemRestrictionGetElementSize(CeedElemRestriction rstr, 784 CeedInt *elem_size) { 785 *elem_size = rstr->elem_size; 786 return CEED_ERROR_SUCCESS; 787 } 788 789 /** 790 @brief Get the size of the l-vector for a CeedElemRestriction 791 792 @param rstr CeedElemRestriction 793 @param[out] l_size Variable to store number of nodes 794 795 @return An error code: 0 - success, otherwise - failure 796 797 @ref Backend 798 **/ 799 int CeedElemRestrictionGetLVectorSize(CeedElemRestriction rstr, 800 CeedInt *l_size) { 801 *l_size = rstr->l_size; 802 return CEED_ERROR_SUCCESS; 803 } 804 805 /** 806 @brief Get the number of components in the elements of a 807 CeedElemRestriction 808 809 @param rstr CeedElemRestriction 810 @param[out] num_comp Variable to store number of components 811 812 @return An error code: 0 - success, otherwise - failure 813 814 @ref Backend 815 **/ 816 int CeedElemRestrictionGetNumComponents(CeedElemRestriction rstr, 817 CeedInt *num_comp) { 818 *num_comp = rstr->num_comp; 819 return CEED_ERROR_SUCCESS; 820 } 821 822 /** 823 @brief Get the number of blocks in a CeedElemRestriction 824 825 @param rstr CeedElemRestriction 826 @param[out] num_block Variable to store number of blocks 827 828 @return An error code: 0 - success, otherwise - failure 829 830 @ref Backend 831 **/ 832 int CeedElemRestrictionGetNumBlocks(CeedElemRestriction rstr, 833 CeedInt *num_block) { 834 *num_block = rstr->num_blk; 835 return CEED_ERROR_SUCCESS; 836 } 837 838 /** 839 @brief Get the size of blocks in the CeedElemRestriction 840 841 @param rstr CeedElemRestriction 842 @param[out] blk_size Variable to store size of blocks 843 844 @return An error code: 0 - success, otherwise - failure 845 846 @ref Backend 847 **/ 848 int CeedElemRestrictionGetBlockSize(CeedElemRestriction rstr, 849 CeedInt *blk_size) { 850 *blk_size = rstr->blk_size; 851 return CEED_ERROR_SUCCESS; 852 } 853 854 /** 855 @brief Get the multiplicity of nodes in a CeedElemRestriction 856 857 @param rstr CeedElemRestriction 858 @param[out] mult Vector to store multiplicity (of size l_size) 859 860 @return An error code: 0 - success, otherwise - failure 861 862 @ref User 863 **/ 864 int CeedElemRestrictionGetMultiplicity(CeedElemRestriction rstr, 865 CeedVector mult) { 866 int ierr; 867 CeedVector e_vec; 868 869 // Create and set e_vec 870 ierr = CeedElemRestrictionCreateVector(rstr, NULL, &e_vec); CeedChk(ierr); 871 ierr = CeedVectorSetValue(e_vec, 1.0); CeedChk(ierr); 872 ierr = CeedVectorSetValue(mult, 0.0); CeedChk(ierr); 873 874 // Apply to get multiplicity 875 ierr = CeedElemRestrictionApply(rstr, CEED_TRANSPOSE, e_vec, mult, 876 CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 877 878 // Cleanup 879 ierr = CeedVectorDestroy(&e_vec); CeedChk(ierr); 880 return CEED_ERROR_SUCCESS; 881 } 882 883 /** 884 @brief View a CeedElemRestriction 885 886 @param[in] rstr CeedElemRestriction to view 887 @param[in] stream Stream to write; typically stdout/stderr or a file 888 889 @return Error code: 0 - success, otherwise - failure 890 891 @ref User 892 **/ 893 int CeedElemRestrictionView(CeedElemRestriction rstr, FILE *stream) { 894 char stridesstr[500]; 895 if (rstr->strides) 896 sprintf(stridesstr, "[%d, %d, %d]", rstr->strides[0], rstr->strides[1], 897 rstr->strides[2]); 898 else 899 sprintf(stridesstr, "%d", rstr->comp_stride); 900 901 fprintf(stream, "%sCeedElemRestriction from (%d, %d) to %d elements with %d " 902 "nodes each and %s %s\n", rstr->blk_size > 1 ? "Blocked " : "", 903 rstr->l_size, rstr->num_comp, rstr->num_elem, rstr->elem_size, 904 rstr->strides ? "strides" : "component stride", stridesstr); 905 return CEED_ERROR_SUCCESS; 906 } 907 908 /** 909 @brief Destroy a CeedElemRestriction 910 911 @param rstr CeedElemRestriction to destroy 912 913 @return An error code: 0 - success, otherwise - failure 914 915 @ref User 916 **/ 917 int CeedElemRestrictionDestroy(CeedElemRestriction *rstr) { 918 int ierr; 919 920 if (!*rstr || --(*rstr)->ref_count > 0) return CEED_ERROR_SUCCESS; 921 if ((*rstr)->num_readers) 922 // LCOV_EXCL_START 923 return CeedError((*rstr)->ceed, CEED_ERROR_ACCESS, 924 "Cannot destroy CeedElemRestriction, " 925 "a process has read access to the offset data"); 926 // LCOV_EXCL_STOP 927 if ((*rstr)->Destroy) { 928 ierr = (*rstr)->Destroy(*rstr); CeedChk(ierr); 929 } 930 ierr = CeedFree(&(*rstr)->strides); CeedChk(ierr); 931 ierr = CeedDestroy(&(*rstr)->ceed); CeedChk(ierr); 932 ierr = CeedFree(rstr); CeedChk(ierr); 933 return CEED_ERROR_SUCCESS; 934 } 935 936 /// @} 937