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] status 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 276 /// @cond DOXYGEN_SKIP 277 static struct CeedElemRestriction_private ceed_elemrestriction_none; 278 /// @endcond 279 280 /// ---------------------------------------------------------------------------- 281 /// CeedElemRestriction Public API 282 /// ---------------------------------------------------------------------------- 283 /// @addtogroup CeedElemRestrictionUser 284 /// @{ 285 286 /// Indicate that the stride is determined by the backend 287 const CeedInt CEED_STRIDES_BACKEND[3] = {}; 288 289 /// Indicate that no CeedElemRestriction is provided by the user 290 const CeedElemRestriction CEED_ELEMRESTRICTION_NONE = 291 &ceed_elemrestriction_none; 292 293 /** 294 @brief Create a CeedElemRestriction 295 296 @param ceed A Ceed object where the CeedElemRestriction will be created 297 @param num_elem Number of elements described in the @a offsets array 298 @param elem_size Size (number of "nodes") per element 299 @param num_comp Number of field components per interpolation node 300 (1 for scalar fields) 301 @param comp_stride Stride between components for the same L-vector "node". 302 Data for node i, component j, element k can be found in 303 the L-vector at index 304 offsets[i + k*elem_size] + j*comp_stride. 305 @param l_size The size of the L-vector. This vector may be larger than 306 the elements and fields given by this restriction. 307 @param mem_type Memory type of the @a offsets array, see CeedMemType 308 @param copy_mode Copy mode for the @a offsets array, see CeedCopyMode 309 @param offsets Array of shape [@a num_elem, @a elem_size]. Row i holds the 310 ordered list of the offsets (into the input CeedVector) 311 for the unknowns corresponding to element i, where 312 0 <= i < @a num_elem. All offsets must be in the range 313 [0, @a l_size - 1]. 314 @param[out] rstr Address of the variable where the newly created 315 CeedElemRestriction will be stored 316 317 @return An error code: 0 - success, otherwise - failure 318 319 @ref User 320 **/ 321 int CeedElemRestrictionCreate(Ceed ceed, CeedInt num_elem, CeedInt elem_size, 322 CeedInt num_comp, CeedInt comp_stride, 323 CeedInt l_size, CeedMemType mem_type, 324 CeedCopyMode copy_mode, const CeedInt *offsets, 325 CeedElemRestriction *rstr) { 326 int ierr; 327 328 if (!ceed->ElemRestrictionCreate) { 329 Ceed delegate; 330 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 331 CeedChk(ierr); 332 333 if (!delegate) 334 // LCOV_EXCL_START 335 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 336 "Backend does not support ElemRestrictionCreate"); 337 // LCOV_EXCL_STOP 338 339 ierr = CeedElemRestrictionCreate(delegate, num_elem, elem_size, num_comp, 340 comp_stride, l_size, mem_type, copy_mode, 341 offsets, rstr); CeedChk(ierr); 342 return CEED_ERROR_SUCCESS; 343 } 344 345 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 346 (*rstr)->ceed = ceed; 347 ceed->ref_count++; 348 (*rstr)->ref_count = 1; 349 (*rstr)->num_elem = num_elem; 350 (*rstr)->elem_size = elem_size; 351 (*rstr)->num_comp = num_comp; 352 (*rstr)->comp_stride = comp_stride; 353 (*rstr)->l_size = l_size; 354 (*rstr)->num_blk = num_elem; 355 (*rstr)->blk_size = 1; 356 ierr = ceed->ElemRestrictionCreate(mem_type, copy_mode, offsets, *rstr); 357 CeedChk(ierr); 358 return CEED_ERROR_SUCCESS; 359 } 360 361 /** 362 @brief Create a strided CeedElemRestriction 363 364 @param ceed A Ceed object where the CeedElemRestriction will be created 365 @param num_elem Number of elements described by the restriction 366 @param elem_size Size (number of "nodes") per element 367 @param num_comp Number of field components per interpolation "node" 368 (1 for scalar fields) 369 @param l_size The size of the L-vector. This vector may be larger than 370 the elements and fields given by this restriction. 371 @param strides Array for strides between [nodes, components, elements]. 372 Data for node i, component j, element k can be found in 373 the L-vector at index 374 i*strides[0] + j*strides[1] + k*strides[2]. 375 @a CEED_STRIDES_BACKEND may be used with vectors created 376 by a Ceed backend. 377 @param rstr Address of the variable where the newly created 378 CeedElemRestriction will be stored 379 380 @return An error code: 0 - success, otherwise - failure 381 382 @ref User 383 **/ 384 int CeedElemRestrictionCreateStrided(Ceed ceed, CeedInt num_elem, 385 CeedInt elem_size, 386 CeedInt num_comp, CeedInt l_size, 387 const CeedInt strides[3], 388 CeedElemRestriction *rstr) { 389 int ierr; 390 391 if (!ceed->ElemRestrictionCreate) { 392 Ceed delegate; 393 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 394 CeedChk(ierr); 395 396 if (!delegate) 397 // LCOV_EXCL_START 398 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 399 "Backend does not support ElemRestrictionCreate"); 400 // LCOV_EXCL_STOP 401 402 ierr = CeedElemRestrictionCreateStrided(delegate, num_elem, elem_size, num_comp, 403 l_size, strides, rstr); 404 CeedChk(ierr); 405 return CEED_ERROR_SUCCESS; 406 } 407 408 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 409 (*rstr)->ceed = ceed; 410 ceed->ref_count++; 411 (*rstr)->ref_count = 1; 412 (*rstr)->num_elem = num_elem; 413 (*rstr)->elem_size = elem_size; 414 (*rstr)->num_comp = num_comp; 415 (*rstr)->l_size = l_size; 416 (*rstr)->num_blk = num_elem; 417 (*rstr)->blk_size = 1; 418 ierr = CeedMalloc(3, &(*rstr)->strides); CeedChk(ierr); 419 for (int i=0; i<3; i++) 420 (*rstr)->strides[i] = strides[i]; 421 ierr = ceed->ElemRestrictionCreate(CEED_MEM_HOST, CEED_OWN_POINTER, NULL, 422 *rstr); 423 CeedChk(ierr); 424 return CEED_ERROR_SUCCESS; 425 } 426 427 /** 428 @brief Create a blocked CeedElemRestriction, typically only called by backends 429 430 @param ceed A Ceed object where the CeedElemRestriction will be created. 431 @param num_elem Number of elements described in the @a offsets array. 432 @param elem_size Size (number of unknowns) per element 433 @param blk_size Number of elements in a block 434 @param num_comp Number of field components per interpolation node 435 (1 for scalar fields) 436 @param comp_stride Stride between components for the same L-vector "node". 437 Data for node i, component j, element k can be found in 438 the L-vector at index 439 offsets[i + k*elem_size] + j*comp_stride. 440 @param l_size The size of the L-vector. This vector may be larger than 441 the elements and fields given by this restriction. 442 @param mem_type Memory type of the @a offsets array, see CeedMemType 443 @param copy_mode Copy mode for the @a offsets array, see CeedCopyMode 444 @param offsets Array of shape [@a num_elem, @a elem_size]. Row i holds the 445 ordered list of the offsets (into the input CeedVector) 446 for the unknowns corresponding to element i, where 447 0 <= i < @a num_elem. All offsets must be in the range 448 [0, @a l_size - 1]. The backend will permute and pad this 449 array to the desired ordering for the blocksize, which is 450 typically given by the backend. The default reordering is 451 to interlace elements. 452 @param rstr Address of the variable where the newly created 453 CeedElemRestriction will be stored 454 455 @return An error code: 0 - success, otherwise - failure 456 457 @ref Backend 458 **/ 459 int CeedElemRestrictionCreateBlocked(Ceed ceed, CeedInt num_elem, 460 CeedInt elem_size, 461 CeedInt blk_size, CeedInt num_comp, 462 CeedInt comp_stride, CeedInt l_size, 463 CeedMemType mem_type, CeedCopyMode copy_mode, 464 const CeedInt *offsets, 465 CeedElemRestriction *rstr) { 466 int ierr; 467 CeedInt *blk_offsets; 468 CeedInt num_blk = (num_elem / blk_size) + !!(num_elem % blk_size); 469 470 if (!ceed->ElemRestrictionCreateBlocked) { 471 Ceed delegate; 472 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 473 CeedChk(ierr); 474 475 if (!delegate) 476 // LCOV_EXCL_START 477 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support " 478 "ElemRestrictionCreateBlocked"); 479 // LCOV_EXCL_STOP 480 481 ierr = CeedElemRestrictionCreateBlocked(delegate, num_elem, elem_size, blk_size, 482 num_comp, comp_stride, l_size, mem_type, 483 copy_mode, offsets, rstr); 484 CeedChk(ierr); 485 return CEED_ERROR_SUCCESS; 486 } 487 488 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 489 490 ierr = CeedCalloc(num_blk*blk_size*elem_size, &blk_offsets); CeedChk(ierr); 491 ierr = CeedPermutePadOffsets(offsets, blk_offsets, num_blk, num_elem, blk_size, 492 elem_size); CeedChk(ierr); 493 494 (*rstr)->ceed = ceed; 495 ceed->ref_count++; 496 (*rstr)->ref_count = 1; 497 (*rstr)->num_elem = num_elem; 498 (*rstr)->elem_size = elem_size; 499 (*rstr)->num_comp = num_comp; 500 (*rstr)->comp_stride = comp_stride; 501 (*rstr)->l_size = l_size; 502 (*rstr)->num_blk = num_blk; 503 (*rstr)->blk_size = blk_size; 504 ierr = ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, 505 (const CeedInt *) blk_offsets, *rstr); CeedChk(ierr); 506 if (copy_mode == CEED_OWN_POINTER) { 507 ierr = CeedFree(&offsets); CeedChk(ierr); 508 } 509 return CEED_ERROR_SUCCESS; 510 } 511 512 /** 513 @brief Create a blocked strided CeedElemRestriction 514 515 @param ceed A Ceed object where the CeedElemRestriction will be created 516 @param num_elem Number of elements described by the restriction 517 @param elem_size Size (number of "nodes") per element 518 @param blk_size Number of elements in a block 519 @param num_comp Number of field components per interpolation node 520 (1 for scalar fields) 521 @param l_size The size of the L-vector. This vector may be larger than 522 the elements and fields given by this restriction. 523 @param strides Array for strides between [nodes, components, elements]. 524 Data for node i, component j, element k can be found in 525 the L-vector at index 526 i*strides[0] + j*strides[1] + k*strides[2]. 527 @a CEED_STRIDES_BACKEND may be used with vectors created 528 by a Ceed backend. 529 @param rstr Address of the variable where the newly created 530 CeedElemRestriction will be stored 531 532 @return An error code: 0 - success, otherwise - failure 533 534 @ref User 535 **/ 536 int CeedElemRestrictionCreateBlockedStrided(Ceed ceed, CeedInt num_elem, 537 CeedInt elem_size, CeedInt blk_size, CeedInt num_comp, CeedInt l_size, 538 const CeedInt strides[3], CeedElemRestriction *rstr) { 539 int ierr; 540 CeedInt num_blk = (num_elem / blk_size) + !!(num_elem % blk_size); 541 542 if (!ceed->ElemRestrictionCreateBlocked) { 543 Ceed delegate; 544 ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 545 CeedChk(ierr); 546 547 if (!delegate) 548 // LCOV_EXCL_START 549 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support " 550 "ElemRestrictionCreateBlocked"); 551 // LCOV_EXCL_STOP 552 553 ierr = CeedElemRestrictionCreateBlockedStrided(delegate, num_elem, elem_size, 554 blk_size, num_comp, l_size, strides, rstr); CeedChk(ierr); 555 return CEED_ERROR_SUCCESS; 556 } 557 558 ierr = CeedCalloc(1, rstr); CeedChk(ierr); 559 560 (*rstr)->ceed = ceed; 561 ceed->ref_count++; 562 (*rstr)->ref_count = 1; 563 (*rstr)->num_elem = num_elem; 564 (*rstr)->elem_size = elem_size; 565 (*rstr)->num_comp = num_comp; 566 (*rstr)->l_size = l_size; 567 (*rstr)->num_blk = num_blk; 568 (*rstr)->blk_size = blk_size; 569 ierr = CeedMalloc(3, &(*rstr)->strides); CeedChk(ierr); 570 for (int i=0; i<3; i++) 571 (*rstr)->strides[i] = strides[i]; 572 ierr = ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, 573 NULL, *rstr); CeedChk(ierr); 574 return CEED_ERROR_SUCCESS; 575 } 576 577 /** 578 @brief Create CeedVectors associated with a CeedElemRestriction 579 580 @param rstr CeedElemRestriction 581 @param l_vec The address of the L-vector to be created, or NULL 582 @param e_vec The address of the E-vector to be created, or NULL 583 584 @return An error code: 0 - success, otherwise - failure 585 586 @ref User 587 **/ 588 int CeedElemRestrictionCreateVector(CeedElemRestriction rstr, CeedVector *l_vec, 589 CeedVector *e_vec) { 590 int ierr; 591 CeedInt e_size, l_size; 592 l_size = rstr->l_size; 593 e_size = rstr->num_blk * rstr->blk_size * rstr->elem_size * rstr->num_comp; 594 if (l_vec) { 595 ierr = CeedVectorCreate(rstr->ceed, l_size, l_vec); CeedChk(ierr); 596 } 597 if (e_vec) { 598 ierr = CeedVectorCreate(rstr->ceed, e_size, e_vec); CeedChk(ierr); 599 } 600 return CEED_ERROR_SUCCESS; 601 } 602 603 /** 604 @brief Restrict an L-vector to an E-vector or apply its transpose 605 606 @param rstr CeedElemRestriction 607 @param t_mode Apply restriction or transpose 608 @param u Input vector (of size @a l_size when t_mode=@ref CEED_NOTRANSPOSE) 609 @param ru Output vector (of shape [@a num_elem * @a elem_size] when 610 t_mode=@ref CEED_NOTRANSPOSE). Ordering of the e-vector is decided 611 by the backend. 612 @param request Request or @ref CEED_REQUEST_IMMEDIATE 613 614 @return An error code: 0 - success, otherwise - failure 615 616 @ref User 617 **/ 618 int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode t_mode, 619 CeedVector u, CeedVector ru, 620 CeedRequest *request) { 621 CeedInt m, n; 622 int ierr; 623 624 if (t_mode == CEED_NOTRANSPOSE) { 625 m = rstr->num_blk * rstr->blk_size * rstr->elem_size * rstr->num_comp; 626 n = rstr->l_size; 627 } else { 628 m = rstr->l_size; 629 n = rstr->num_blk * rstr->blk_size * rstr->elem_size * rstr->num_comp; 630 } 631 if (n != u->length) 632 // LCOV_EXCL_START 633 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 634 "Input vector size %d not compatible with " 635 "element restriction (%d, %d)", u->length, m, n); 636 // LCOV_EXCL_STOP 637 if (m != ru->length) 638 // LCOV_EXCL_START 639 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 640 "Output vector size %d not compatible with " 641 "element restriction (%d, %d)", ru->length, m, n); 642 // LCOV_EXCL_STOP 643 ierr = rstr->Apply(rstr, t_mode, u, ru, request); CeedChk(ierr); 644 return CEED_ERROR_SUCCESS; 645 } 646 647 /** 648 @brief Restrict an L-vector to a block of an E-vector or apply its transpose 649 650 @param rstr CeedElemRestriction 651 @param block Block number to restrict to/from, i.e. block=0 will handle 652 elements [0 : blk_size] and block=3 will handle elements 653 [3*blk_size : 4*blk_size] 654 @param t_mode Apply restriction or transpose 655 @param u Input vector (of size @a l_size when t_mode=@ref CEED_NOTRANSPOSE) 656 @param ru Output vector (of shape [@a blk_size * @a elem_size] when 657 t_mode=@ref CEED_NOTRANSPOSE). Ordering of the e-vector is decided 658 by the backend. 659 @param request Request or @ref CEED_REQUEST_IMMEDIATE 660 661 @return An error code: 0 - success, otherwise - failure 662 663 @ref Backend 664 **/ 665 int CeedElemRestrictionApplyBlock(CeedElemRestriction rstr, CeedInt block, 666 CeedTransposeMode t_mode, CeedVector u, 667 CeedVector ru, CeedRequest *request) { 668 CeedInt m, n; 669 int ierr; 670 671 if (t_mode == CEED_NOTRANSPOSE) { 672 m = rstr->blk_size * rstr->elem_size * rstr->num_comp; 673 n = rstr->l_size; 674 } else { 675 m = rstr->l_size; 676 n = rstr->blk_size * rstr->elem_size * rstr->num_comp; 677 } 678 if (n != u->length) 679 // LCOV_EXCL_START 680 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 681 "Input vector size %d not compatible with " 682 "element restriction (%d, %d)", u->length, m, n); 683 // LCOV_EXCL_STOP 684 if (m != ru->length) 685 // LCOV_EXCL_START 686 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 687 "Output vector size %d not compatible with " 688 "element restriction (%d, %d)", ru->length, m, n); 689 // LCOV_EXCL_STOP 690 if (rstr->blk_size*block > rstr->num_elem) 691 // LCOV_EXCL_START 692 return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, 693 "Cannot retrieve block %d, element %d > " 694 "total elements %d", block, rstr->blk_size*block, 695 rstr->num_elem); 696 // LCOV_EXCL_STOP 697 ierr = rstr->ApplyBlock(rstr, block, t_mode, u, ru, request); 698 CeedChk(ierr); 699 return CEED_ERROR_SUCCESS; 700 } 701 702 /** 703 @brief Get the L-vector component stride 704 705 @param rstr CeedElemRestriction 706 @param[out] comp_stride Variable to store component stride 707 708 @return An error code: 0 - success, otherwise - failure 709 710 @ref Backend 711 **/ 712 int CeedElemRestrictionGetCompStride(CeedElemRestriction rstr, 713 CeedInt *comp_stride) { 714 *comp_stride = rstr->comp_stride; 715 return CEED_ERROR_SUCCESS; 716 } 717 718 /** 719 @brief Get the total number of elements in the range of a CeedElemRestriction 720 721 @param rstr CeedElemRestriction 722 @param[out] num_elem Variable to store number of elements 723 724 @return An error code: 0 - success, otherwise - failure 725 726 @ref Backend 727 **/ 728 int CeedElemRestrictionGetNumElements(CeedElemRestriction rstr, 729 CeedInt *num_elem) { 730 *num_elem = rstr->num_elem; 731 return CEED_ERROR_SUCCESS; 732 } 733 734 /** 735 @brief Get the size of elements in the CeedElemRestriction 736 737 @param rstr CeedElemRestriction 738 @param[out] elem_size Variable to store size of elements 739 740 @return An error code: 0 - success, otherwise - failure 741 742 @ref Backend 743 **/ 744 int CeedElemRestrictionGetElementSize(CeedElemRestriction rstr, 745 CeedInt *elem_size) { 746 *elem_size = rstr->elem_size; 747 return CEED_ERROR_SUCCESS; 748 } 749 750 /** 751 @brief Get the size of the l-vector for a CeedElemRestriction 752 753 @param rstr CeedElemRestriction 754 @param[out] l_size Variable to store number of nodes 755 756 @return An error code: 0 - success, otherwise - failure 757 758 @ref Backend 759 **/ 760 int CeedElemRestrictionGetLVectorSize(CeedElemRestriction rstr, 761 CeedInt *l_size) { 762 *l_size = rstr->l_size; 763 return CEED_ERROR_SUCCESS; 764 } 765 766 /** 767 @brief Get the number of components in the elements of a 768 CeedElemRestriction 769 770 @param rstr CeedElemRestriction 771 @param[out] num_comp Variable to store number of components 772 773 @return An error code: 0 - success, otherwise - failure 774 775 @ref Backend 776 **/ 777 int CeedElemRestrictionGetNumComponents(CeedElemRestriction rstr, 778 CeedInt *num_comp) { 779 *num_comp = rstr->num_comp; 780 return CEED_ERROR_SUCCESS; 781 } 782 783 /** 784 @brief Get the number of blocks in a CeedElemRestriction 785 786 @param rstr CeedElemRestriction 787 @param[out] num_block Variable to store number of blocks 788 789 @return An error code: 0 - success, otherwise - failure 790 791 @ref Backend 792 **/ 793 int CeedElemRestrictionGetNumBlocks(CeedElemRestriction rstr, 794 CeedInt *num_block) { 795 *num_block = rstr->num_blk; 796 return CEED_ERROR_SUCCESS; 797 } 798 799 /** 800 @brief Get the size of blocks in the CeedElemRestriction 801 802 @param rstr CeedElemRestriction 803 @param[out] blk_size Variable to store size of blocks 804 805 @return An error code: 0 - success, otherwise - failure 806 807 @ref Backend 808 **/ 809 int CeedElemRestrictionGetBlockSize(CeedElemRestriction rstr, 810 CeedInt *blk_size) { 811 *blk_size = rstr->blk_size; 812 return CEED_ERROR_SUCCESS; 813 } 814 815 /** 816 @brief Get the multiplicity of nodes in a CeedElemRestriction 817 818 @param rstr CeedElemRestriction 819 @param[out] mult Vector to store multiplicity (of size l_size) 820 821 @return An error code: 0 - success, otherwise - failure 822 823 @ref User 824 **/ 825 int CeedElemRestrictionGetMultiplicity(CeedElemRestriction rstr, 826 CeedVector mult) { 827 int ierr; 828 CeedVector e_vec; 829 830 // Create and set e_vec 831 ierr = CeedElemRestrictionCreateVector(rstr, NULL, &e_vec); CeedChk(ierr); 832 ierr = CeedVectorSetValue(e_vec, 1.0); CeedChk(ierr); 833 ierr = CeedVectorSetValue(mult, 0.0); CeedChk(ierr); 834 835 // Apply to get multiplicity 836 ierr = CeedElemRestrictionApply(rstr, CEED_TRANSPOSE, e_vec, mult, 837 CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 838 839 // Cleanup 840 ierr = CeedVectorDestroy(&e_vec); CeedChk(ierr); 841 return CEED_ERROR_SUCCESS; 842 } 843 844 /** 845 @brief View a CeedElemRestriction 846 847 @param[in] rstr CeedElemRestriction to view 848 @param[in] stream Stream to write; typically stdout/stderr or a file 849 850 @return Error code: 0 - success, otherwise - failure 851 852 @ref User 853 **/ 854 int CeedElemRestrictionView(CeedElemRestriction rstr, FILE *stream) { 855 char stridesstr[500]; 856 if (rstr->strides) 857 sprintf(stridesstr, "[%d, %d, %d]", rstr->strides[0], rstr->strides[1], 858 rstr->strides[2]); 859 else 860 sprintf(stridesstr, "%d", rstr->comp_stride); 861 862 fprintf(stream, "%sCeedElemRestriction from (%d, %d) to %d elements with %d " 863 "nodes each and %s %s\n", rstr->blk_size > 1 ? "Blocked " : "", 864 rstr->l_size, rstr->num_comp, rstr->num_elem, rstr->elem_size, 865 rstr->strides ? "strides" : "component stride", stridesstr); 866 return CEED_ERROR_SUCCESS; 867 } 868 869 /** 870 @brief Destroy a CeedElemRestriction 871 872 @param rstr CeedElemRestriction to destroy 873 874 @return An error code: 0 - success, otherwise - failure 875 876 @ref User 877 **/ 878 int CeedElemRestrictionDestroy(CeedElemRestriction *rstr) { 879 int ierr; 880 881 if (!*rstr || --(*rstr)->ref_count > 0) return CEED_ERROR_SUCCESS; 882 if ((*rstr)->num_readers) 883 return CeedError((*rstr)->ceed, CEED_ERROR_ACCESS, 884 "Cannot destroy CeedElemRestriction, " 885 "a process has read access to the offset data"); 886 if ((*rstr)->Destroy) { 887 ierr = (*rstr)->Destroy(*rstr); CeedChk(ierr); 888 } 889 ierr = CeedFree(&(*rstr)->strides); CeedChk(ierr); 890 ierr = CeedDestroy(&(*rstr)->ceed); CeedChk(ierr); 891 ierr = CeedFree(rstr); CeedChk(ierr); 892 return CEED_ERROR_SUCCESS; 893 } 894 895 /// @} 896