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 72 @brief Get the strides of a strided CeedElemRestriction 73 74 @param rstr CeedElemRestriction 75 @param[out] strides Variable to store strides array 76 77 @return An error code: 0 - success, otherwise - failure 78 79 @ref Backend 80 **/ 81 int CeedElemRestrictionGetStrides(CeedElemRestriction rstr, 82 CeedInt (*strides)[3]) { 83 if (!rstr->strides) 84 // LCOV_EXCL_START 85 return CeedError(rstr->ceed, CEED_ERROR_MINOR, 86 "ElemRestriction has no stride data"); 87 // LCOV_EXCL_STOP 88 89 for (int i=0; i<3; i++) 90 (*strides)[i] = rstr->strides[i]; 91 return CEED_ERROR_SUCCESS; 92 } 93 94 /** 95 @brief Get read-only access to a CeedElemRestriction offsets array by memtype 96 97 @param rstr CeedElemRestriction to retrieve offsets 98 @param mem_type Memory type on which to access the array. If the backend 99 uses a different memory type, this will perform a copy 100 (possibly cached). 101 @param[out] offsets Array on memory type mem_type 102 103 @return An error code: 0 - success, otherwise - failure 104 105 @ref User 106 **/ 107 int CeedElemRestrictionGetOffsets(CeedElemRestriction rstr, 108 CeedMemType mem_type, 109 const CeedInt **offsets) { 110 int ierr; 111 112 if (!rstr->GetOffsets) 113 // LCOV_EXCL_START 114 return CeedError(rstr->ceed, CEED_ERROR_UNSUPPORTED, 115 "Backend does not support GetOffsets"); 116 // LCOV_EXCL_STOP 117 118 ierr = rstr->GetOffsets(rstr, mem_type, offsets); CeedChk(ierr); 119 rstr->num_readers++; 120 return CEED_ERROR_SUCCESS; 121 } 122 123 /** 124 @brief Restore an offsets array obtained using CeedElemRestrictionGetOffsets() 125 126 @param rstr CeedElemRestriction to restore 127 @param offsets Array of offset data 128 129 @return An error code: 0 - success, otherwise - failure 130 131 @ref User 132 **/ 133 int CeedElemRestrictionRestoreOffsets(CeedElemRestriction rstr, 134 const CeedInt **offsets) { 135 *offsets = NULL; 136 rstr->num_readers--; 137 return CEED_ERROR_SUCCESS; 138 } 139 140 /** 141 @brief Get the strided status of a CeedElemRestriction 142 143 @param rstr CeedElemRestriction 144 @param[out] is_strided Variable to store strided status, 1 if strided else 0 145 146 @return An error code: 0 - success, otherwise - failure 147 148 @ref Backend 149 **/ 150 int CeedElemRestrictionIsStrided(CeedElemRestriction rstr, bool *is_strided) { 151 *is_strided = rstr->strides ? true : false; 152 return CEED_ERROR_SUCCESS; 153 } 154 155 /** 156 @brief Get the backend stride status of a CeedElemRestriction 157 158 @param rstr CeedElemRestriction 159 @param[out] has_backend_strides Variable to store stride status 160 161 @return An error code: 0 - success, otherwise - failure 162 163 @ref Backend 164 **/ 165 int CeedElemRestrictionHasBackendStrides(CeedElemRestriction rstr, 166 bool *has_backend_strides) { 167 if (!rstr->strides) 168 // LCOV_EXCL_START 169 return CeedError(rstr->ceed, CEED_ERROR_MINOR, 170 "ElemRestriction has no stride data"); 171 // LCOV_EXCL_STOP 172 173 *has_backend_strides = ((rstr->strides[0] == CEED_STRIDES_BACKEND[0]) && 174 (rstr->strides[1] == CEED_STRIDES_BACKEND[1]) && 175 (rstr->strides[2] == CEED_STRIDES_BACKEND[2])); 176 return CEED_ERROR_SUCCESS; 177 } 178 179 /** 180 181 @brief Get the E-vector layout of a CeedElemRestriction 182 183 @param rstr CeedElemRestriction 184 @param[out] layout Variable to store layout array, 185 stored as [nodes, components, elements]. 186 The data for node i, component j, element k in the 187 E-vector is given by 188 i*layout[0] + j*layout[1] + k*layout[2] 189 190 @return An error code: 0 - success, otherwise - failure 191 192 @ref Backend 193 **/ 194 int CeedElemRestrictionGetELayout(CeedElemRestriction rstr, 195 CeedInt (*layout)[3]) { 196 if (!rstr->layout[0]) 197 // LCOV_EXCL_START 198 return CeedError(rstr->ceed, CEED_ERROR_MINOR, 199 "ElemRestriction has no layout data"); 200 // LCOV_EXCL_STOP 201 202 for (int i=0; i<3; i++) 203 (*layout)[i] = rstr->layout[i]; 204 return CEED_ERROR_SUCCESS; 205 } 206 207 /** 208 209 @brief Set the E-vector layout of a CeedElemRestriction 210 211 @param rstr CeedElemRestriction 212 @param layout Variable to containing layout array, 213 stored as [nodes, components, elements]. 214 The data for node i, component j, element k in the 215 E-vector is given by 216 i*layout[0] + j*layout[1] + k*layout[2] 217 218 @return An error code: 0 - success, otherwise - failure 219 220 @ref Backend 221 **/ 222 int CeedElemRestrictionSetELayout(CeedElemRestriction rstr, 223 CeedInt layout[3]) { 224 for (int i = 0; i<3; i++) 225 rstr->layout[i] = layout[i]; 226 return CEED_ERROR_SUCCESS; 227 } 228 229 /** 230 @brief Get the backend data of a CeedElemRestriction 231 232 @param rstr CeedElemRestriction 233 @param[out] data Variable to store data 234 235 @return An error code: 0 - success, otherwise - failure 236 237 @ref Backend 238 **/ 239 int CeedElemRestrictionGetData(CeedElemRestriction rstr, void *data) { 240 *(void **)data = rstr->data; 241 return CEED_ERROR_SUCCESS; 242 } 243 244 /** 245 @brief Set the backend data of a CeedElemRestriction 246 247 @param[out] rstr CeedElemRestriction 248 @param data Data to set 249 250 @return An error code: 0 - success, otherwise - failure 251 252 @ref Backend 253 **/ 254 int CeedElemRestrictionSetData(CeedElemRestriction rstr, void *data) { 255 rstr->data = data; 256 return CEED_ERROR_SUCCESS; 257 } 258 259 /** 260 @brief Increment the reference counter for a CeedElemRestriction 261 262 @param rstr ElemRestriction to increment the reference counter 263 264 @return An error code: 0 - success, otherwise - failure 265 266 @ref Backend 267 **/ 268 int CeedElemRestrictionReference(CeedElemRestriction rstr) { 269 rstr->ref_count++; 270 return CEED_ERROR_SUCCESS; 271 } 272 273 /// @} 274 275 /// @cond DOXYGEN_SKIP 276 static struct CeedElemRestriction_private ceed_elemrestriction_none; 277 /// @endcond 278 279 /// ---------------------------------------------------------------------------- 280 /// CeedElemRestriction Public API 281 /// ---------------------------------------------------------------------------- 282 /// @addtogroup CeedElemRestrictionUser 283 /// @{ 284 285 /// Indicate that the stride is determined by the backend 286 const CeedInt CEED_STRIDES_BACKEND[3] = {0}; 287 288 /// Indicate that no CeedElemRestriction is provided by the user 289 const CeedElemRestriction CEED_ELEMRESTRICTION_NONE = 290 &ceed_elemrestriction_none; 291 292 /** 293 @brief Create a CeedElemRestriction 294 295 @param ceed A Ceed object where the CeedElemRestriction will be created 296 @param num_elem Number of elements described in the @a offsets array 297 @param elem_size Size (number of "nodes") per element 298 @param num_comp Number of field components per interpolation node 299 (1 for scalar fields) 300 @param comp_stride Stride between components for the same L-vector "node". 301 Data for node i, component j, element k can be found in 302 the L-vector at index 303 offsets[i + k*elem_size] + j*comp_stride. 304 @param l_size The size of the L-vector. This vector may be larger than 305 the elements and fields given by this restriction. 306 @param mem_type Memory type of the @a offsets array, see CeedMemType 307 @param copy_mode Copy mode for the @a offsets array, see CeedCopyMode 308 @param offsets Array of shape [@a num_elem, @a elem_size]. Row i holds the 309 ordered list of the offsets (into the input CeedVector) 310 for the unknowns corresponding to element i, where 311 0 <= i < @a num_elem. All offsets must be in the range 312 [0, @a l_size - 1]. 313 @param[out] rstr Address of the variable where the newly created 314 CeedElemRestriction will be stored 315 316 @return An error code: 0 - success, otherwise - failure 317 318 @ref User 319 **/ 320 int CeedElemRestrictionCreate(Ceed ceed, CeedInt num_elem, CeedInt elem_size, 321 CeedInt num_comp, CeedInt comp_stride, 322 CeedInt l_size, CeedMemType mem_type, 323 CeedCopyMode copy_mode, const CeedInt *offsets, 324 CeedElemRestriction *rstr) { 325 int ierr; 326 327 if (!ceed->ElemRestrictionCreate) { 328 Ceed delegate; 329 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 330 CeedChk(ierr); 331 332 if (!delegate) 333 // LCOV_EXCL_START 334 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 335 "Backend does not support ElemRestrictionCreate"); 336 // LCOV_EXCL_STOP 337 338 ierr = CeedElemRestrictionCreate(delegate, num_elem, elem_size, num_comp, 339 comp_stride, l_size, mem_type, copy_mode, 340 offsets, rstr); CeedChk(ierr); 341 return CEED_ERROR_SUCCESS; 342 } 343 344 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 345 (*rstr)->ceed = ceed; 346 ierr = CeedReference(ceed); CeedChk(ierr); 347 (*rstr)->ref_count = 1; 348 (*rstr)->num_elem = num_elem; 349 (*rstr)->elem_size = elem_size; 350 (*rstr)->num_comp = num_comp; 351 (*rstr)->comp_stride = comp_stride; 352 (*rstr)->l_size = l_size; 353 (*rstr)->num_blk = num_elem; 354 (*rstr)->blk_size = 1; 355 ierr = ceed->ElemRestrictionCreate(mem_type, copy_mode, offsets, *rstr); 356 CeedChk(ierr); 357 return CEED_ERROR_SUCCESS; 358 } 359 360 /** 361 @brief Create a strided CeedElemRestriction 362 363 @param ceed A Ceed object where the CeedElemRestriction will be created 364 @param num_elem Number of elements described by the restriction 365 @param elem_size Size (number of "nodes") per element 366 @param num_comp Number of field components per interpolation "node" 367 (1 for scalar fields) 368 @param l_size The size of the L-vector. This vector may be larger than 369 the elements and fields given by this restriction. 370 @param strides Array for strides between [nodes, components, elements]. 371 Data for node i, component j, element k can be found in 372 the L-vector at index 373 i*strides[0] + j*strides[1] + k*strides[2]. 374 @a CEED_STRIDES_BACKEND may be used with vectors created 375 by a Ceed backend. 376 @param rstr Address of the variable where the newly created 377 CeedElemRestriction will be stored 378 379 @return An error code: 0 - success, otherwise - failure 380 381 @ref User 382 **/ 383 int CeedElemRestrictionCreateStrided(Ceed ceed, CeedInt num_elem, 384 CeedInt elem_size, 385 CeedInt num_comp, CeedInt l_size, 386 const CeedInt strides[3], 387 CeedElemRestriction *rstr) { 388 int ierr; 389 390 if (!ceed->ElemRestrictionCreate) { 391 Ceed delegate; 392 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 393 CeedChk(ierr); 394 395 if (!delegate) 396 // LCOV_EXCL_START 397 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 398 "Backend does not support ElemRestrictionCreate"); 399 // LCOV_EXCL_STOP 400 401 ierr = CeedElemRestrictionCreateStrided(delegate, num_elem, elem_size, num_comp, 402 l_size, strides, rstr); 403 CeedChk(ierr); 404 return CEED_ERROR_SUCCESS; 405 } 406 407 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 408 (*rstr)->ceed = ceed; 409 ierr = CeedReference(ceed); CeedChk(ierr); 410 (*rstr)->ref_count = 1; 411 (*rstr)->num_elem = num_elem; 412 (*rstr)->elem_size = elem_size; 413 (*rstr)->num_comp = num_comp; 414 (*rstr)->l_size = l_size; 415 (*rstr)->num_blk = num_elem; 416 (*rstr)->blk_size = 1; 417 ierr = CeedMalloc(3, &(*rstr)->strides); CeedChk(ierr); 418 for (int i=0; i<3; i++) 419 (*rstr)->strides[i] = strides[i]; 420 ierr = ceed->ElemRestrictionCreate(CEED_MEM_HOST, CEED_OWN_POINTER, NULL, 421 *rstr); 422 CeedChk(ierr); 423 return CEED_ERROR_SUCCESS; 424 } 425 426 /** 427 @brief Create a blocked CeedElemRestriction, typically only called by backends 428 429 @param ceed A Ceed object where the CeedElemRestriction will be created. 430 @param num_elem Number of elements described in the @a offsets array. 431 @param elem_size Size (number of unknowns) per element 432 @param blk_size Number of elements in a block 433 @param num_comp Number of field components per interpolation node 434 (1 for scalar fields) 435 @param comp_stride Stride between components for the same L-vector "node". 436 Data for node i, component j, element k can be found in 437 the L-vector at index 438 offsets[i + k*elem_size] + j*comp_stride. 439 @param l_size The size of the L-vector. This vector may be larger than 440 the elements and fields given by this restriction. 441 @param mem_type Memory type of the @a offsets array, see CeedMemType 442 @param copy_mode Copy mode for the @a offsets array, see CeedCopyMode 443 @param offsets Array of shape [@a num_elem, @a elem_size]. Row i holds the 444 ordered list of the offsets (into the input CeedVector) 445 for the unknowns corresponding to element i, where 446 0 <= i < @a num_elem. All offsets must be in the range 447 [0, @a l_size - 1]. The backend will permute and pad this 448 array to the desired ordering for the blocksize, which is 449 typically given by the backend. The default reordering is 450 to interlace elements. 451 @param rstr Address of the variable where the newly created 452 CeedElemRestriction will be stored 453 454 @return An error code: 0 - success, otherwise - failure 455 456 @ref Backend 457 **/ 458 int CeedElemRestrictionCreateBlocked(Ceed ceed, CeedInt num_elem, 459 CeedInt elem_size, 460 CeedInt blk_size, CeedInt num_comp, 461 CeedInt comp_stride, CeedInt l_size, 462 CeedMemType mem_type, CeedCopyMode copy_mode, 463 const CeedInt *offsets, 464 CeedElemRestriction *rstr) { 465 int ierr; 466 CeedInt *blk_offsets; 467 CeedInt num_blk = (num_elem / blk_size) + !!(num_elem % blk_size); 468 469 if (!ceed->ElemRestrictionCreateBlocked) { 470 Ceed delegate; 471 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 472 CeedChk(ierr); 473 474 if (!delegate) 475 // LCOV_EXCL_START 476 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support " 477 "ElemRestrictionCreateBlocked"); 478 // LCOV_EXCL_STOP 479 480 ierr = CeedElemRestrictionCreateBlocked(delegate, num_elem, elem_size, blk_size, 481 num_comp, comp_stride, l_size, mem_type, 482 copy_mode, offsets, rstr); 483 CeedChk(ierr); 484 return CEED_ERROR_SUCCESS; 485 } 486 487 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 488 489 ierr = CeedCalloc(num_blk*blk_size*elem_size, &blk_offsets); CeedChk(ierr); 490 ierr = CeedPermutePadOffsets(offsets, blk_offsets, num_blk, num_elem, blk_size, 491 elem_size); CeedChk(ierr); 492 493 (*rstr)->ceed = ceed; 494 ierr = CeedReference(ceed); CeedChk(ierr); 495 (*rstr)->ref_count = 1; 496 (*rstr)->num_elem = num_elem; 497 (*rstr)->elem_size = elem_size; 498 (*rstr)->num_comp = num_comp; 499 (*rstr)->comp_stride = comp_stride; 500 (*rstr)->l_size = l_size; 501 (*rstr)->num_blk = num_blk; 502 (*rstr)->blk_size = blk_size; 503 ierr = ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, 504 (const CeedInt *) blk_offsets, *rstr); CeedChk(ierr); 505 if (copy_mode == CEED_OWN_POINTER) { 506 ierr = CeedFree(&offsets); CeedChk(ierr); 507 } 508 return CEED_ERROR_SUCCESS; 509 } 510 511 /** 512 @brief Create a blocked strided CeedElemRestriction 513 514 @param ceed A Ceed object where the CeedElemRestriction will be created 515 @param num_elem Number of elements described by the restriction 516 @param elem_size Size (number of "nodes") per element 517 @param blk_size Number of elements in a block 518 @param num_comp Number of field components per interpolation node 519 (1 for scalar fields) 520 @param l_size The size of the L-vector. This vector may be larger than 521 the elements and fields given by this restriction. 522 @param strides Array for strides between [nodes, components, elements]. 523 Data for node i, component j, element k can be found in 524 the L-vector at index 525 i*strides[0] + j*strides[1] + k*strides[2]. 526 @a CEED_STRIDES_BACKEND may be used with vectors created 527 by a Ceed backend. 528 @param rstr Address of the variable where the newly created 529 CeedElemRestriction will be stored 530 531 @return An error code: 0 - success, otherwise - failure 532 533 @ref User 534 **/ 535 int CeedElemRestrictionCreateBlockedStrided(Ceed ceed, CeedInt num_elem, 536 CeedInt elem_size, CeedInt blk_size, CeedInt num_comp, CeedInt l_size, 537 const CeedInt strides[3], CeedElemRestriction *rstr) { 538 int ierr; 539 CeedInt num_blk = (num_elem / blk_size) + !!(num_elem % blk_size); 540 541 if (!ceed->ElemRestrictionCreateBlocked) { 542 Ceed delegate; 543 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 544 CeedChk(ierr); 545 546 if (!delegate) 547 // LCOV_EXCL_START 548 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support " 549 "ElemRestrictionCreateBlocked"); 550 // LCOV_EXCL_STOP 551 552 ierr = CeedElemRestrictionCreateBlockedStrided(delegate, num_elem, elem_size, 553 blk_size, num_comp, l_size, strides, rstr); CeedChk(ierr); 554 return CEED_ERROR_SUCCESS; 555 } 556 557 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 558 559 (*rstr)->ceed = ceed; 560 ierr = CeedReference(ceed); CeedChk(ierr); 561 (*rstr)->ref_count = 1; 562 (*rstr)->num_elem = num_elem; 563 (*rstr)->elem_size = elem_size; 564 (*rstr)->num_comp = num_comp; 565 (*rstr)->l_size = l_size; 566 (*rstr)->num_blk = num_blk; 567 (*rstr)->blk_size = blk_size; 568 ierr = CeedMalloc(3, &(*rstr)->strides); CeedChk(ierr); 569 for (int i=0; i<3; i++) 570 (*rstr)->strides[i] = strides[i]; 571 ierr = ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, 572 NULL, *rstr); CeedChk(ierr); 573 return CEED_ERROR_SUCCESS; 574 } 575 576 /** 577 @brief Copy the pointer to a CeedElemRestriction. Both pointers should 578 be destroyed with `CeedElemRestrictionDestroy()`; 579 Note: If `*rstr_copy` is non-NULL, then it is assumed that 580 `*rstr_copy` is a pointer to a CeedElemRestriction. This 581 CeedElemRestriction will be destroyed if `*rstr_copy` is the 582 only reference to this CeedElemRestriction. 583 584 @param rstr CeedElemRestriction to copy reference to 585 @param[out] rstr_copy Variable to store copied reference 586 587 @return An error code: 0 - success, otherwise - failure 588 589 @ref User 590 **/ 591 int CeedElemRestrictionReferenceCopy(CeedElemRestriction rstr, 592 CeedElemRestriction *rstr_copy) { 593 int ierr; 594 595 ierr = CeedElemRestrictionReference(rstr); CeedChk(ierr); 596 ierr = CeedElemRestrictionDestroy(rstr_copy); CeedChk(ierr); 597 *rstr_copy = rstr; 598 return CEED_ERROR_SUCCESS; 599 } 600 601 /** 602 @brief Create CeedVectors associated with a CeedElemRestriction 603 604 @param rstr CeedElemRestriction 605 @param l_vec The address of the L-vector to be created, or NULL 606 @param e_vec The address of the E-vector to be created, or NULL 607 608 @return An error code: 0 - success, otherwise - failure 609 610 @ref User 611 **/ 612 int CeedElemRestrictionCreateVector(CeedElemRestriction rstr, CeedVector *l_vec, 613 CeedVector *e_vec) { 614 int ierr; 615 CeedInt e_size, l_size; 616 l_size = rstr->l_size; 617 e_size = rstr->num_blk * rstr->blk_size * rstr->elem_size * rstr->num_comp; 618 if (l_vec) { 619 ierr = CeedVectorCreate(rstr->ceed, l_size, l_vec); CeedChk(ierr); 620 } 621 if (e_vec) { 622 ierr = CeedVectorCreate(rstr->ceed, e_size, e_vec); CeedChk(ierr); 623 } 624 return CEED_ERROR_SUCCESS; 625 } 626 627 /** 628 @brief Restrict an L-vector to an E-vector or apply its transpose 629 630 @param rstr CeedElemRestriction 631 @param t_mode Apply restriction or transpose 632 @param u Input vector (of size @a l_size when t_mode=@ref CEED_NOTRANSPOSE) 633 @param ru Output vector (of shape [@a num_elem * @a elem_size] when 634 t_mode=@ref CEED_NOTRANSPOSE). Ordering of the e-vector is decided 635 by the backend. 636 @param request Request or @ref CEED_REQUEST_IMMEDIATE 637 638 @return An error code: 0 - success, otherwise - failure 639 640 @ref User 641 **/ 642 int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode t_mode, 643 CeedVector u, CeedVector ru, 644 CeedRequest *request) { 645 CeedInt m, n; 646 int ierr; 647 648 if (t_mode == CEED_NOTRANSPOSE) { 649 m = rstr->num_blk * rstr->blk_size * rstr->elem_size * rstr->num_comp; 650 n = rstr->l_size; 651 } else { 652 m = rstr->l_size; 653 n = rstr->num_blk * rstr->blk_size * rstr->elem_size * rstr->num_comp; 654 } 655 if (n != u->length) 656 // LCOV_EXCL_START 657 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 658 "Input vector size %d not compatible with " 659 "element restriction (%d, %d)", u->length, m, n); 660 // LCOV_EXCL_STOP 661 if (m != ru->length) 662 // LCOV_EXCL_START 663 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 664 "Output vector size %d not compatible with " 665 "element restriction (%d, %d)", ru->length, m, n); 666 // LCOV_EXCL_STOP 667 ierr = rstr->Apply(rstr, t_mode, u, ru, request); CeedChk(ierr); 668 return CEED_ERROR_SUCCESS; 669 } 670 671 /** 672 @brief Restrict an L-vector to a block of an E-vector or apply its transpose 673 674 @param rstr CeedElemRestriction 675 @param block Block number to restrict to/from, i.e. block=0 will handle 676 elements [0 : blk_size] and block=3 will handle elements 677 [3*blk_size : 4*blk_size] 678 @param t_mode Apply restriction or transpose 679 @param u Input vector (of size @a l_size when t_mode=@ref CEED_NOTRANSPOSE) 680 @param ru Output vector (of shape [@a blk_size * @a elem_size] when 681 t_mode=@ref CEED_NOTRANSPOSE). Ordering of the e-vector is decided 682 by the backend. 683 @param request Request or @ref CEED_REQUEST_IMMEDIATE 684 685 @return An error code: 0 - success, otherwise - failure 686 687 @ref Backend 688 **/ 689 int CeedElemRestrictionApplyBlock(CeedElemRestriction rstr, CeedInt block, 690 CeedTransposeMode t_mode, CeedVector u, 691 CeedVector ru, CeedRequest *request) { 692 CeedInt m, n; 693 int ierr; 694 695 if (t_mode == CEED_NOTRANSPOSE) { 696 m = rstr->blk_size * rstr->elem_size * rstr->num_comp; 697 n = rstr->l_size; 698 } else { 699 m = rstr->l_size; 700 n = rstr->blk_size * rstr->elem_size * rstr->num_comp; 701 } 702 if (n != u->length) 703 // LCOV_EXCL_START 704 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 705 "Input vector size %d not compatible with " 706 "element restriction (%d, %d)", u->length, m, n); 707 // LCOV_EXCL_STOP 708 if (m != ru->length) 709 // LCOV_EXCL_START 710 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 711 "Output vector size %d not compatible with " 712 "element restriction (%d, %d)", ru->length, m, n); 713 // LCOV_EXCL_STOP 714 if (rstr->blk_size*block > rstr->num_elem) 715 // LCOV_EXCL_START 716 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 717 "Cannot retrieve block %d, element %d > " 718 "total elements %d", block, rstr->blk_size*block, 719 rstr->num_elem); 720 // LCOV_EXCL_STOP 721 ierr = rstr->ApplyBlock(rstr, block, t_mode, u, ru, request); 722 CeedChk(ierr); 723 return CEED_ERROR_SUCCESS; 724 } 725 726 /** 727 @brief Get the Ceed associated with a CeedElemRestriction 728 729 @param rstr CeedElemRestriction 730 @param[out] ceed Variable to store Ceed 731 732 @return An error code: 0 - success, otherwise - failure 733 734 @ref Advanced 735 **/ 736 int CeedElemRestrictionGetCeed(CeedElemRestriction rstr, Ceed *ceed) { 737 *ceed = rstr->ceed; 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 Advanced 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 Advanced 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 Advanced 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 Advanced 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 Advanced 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 Advanced 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 Advanced 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