1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 #include <ceed-impl.h> 9 #include <ceed.h> 10 #include <ceed/backend.h> 11 #include <stdbool.h> 12 #include <stdio.h> 13 #include <string.h> 14 15 /// @file 16 /// Implementation of CeedElemRestriction interfaces 17 18 /// ---------------------------------------------------------------------------- 19 /// CeedElemRestriction Library Internal Functions 20 /// ---------------------------------------------------------------------------- 21 /// @addtogroup CeedElemRestrictionDeveloper 22 /// @{ 23 24 /** 25 @brief Permute and pad offsets for a blocked `CeedElemRestriction` 26 27 @param[in] offsets Array of shape `[num_elem, elem_size]` 28 @param[out] block_offsets Array of permuted and padded array values of shape `[num_block, elem_size, block_size]` 29 @param[in] num_block Number of blocks 30 @param[in] num_elem Number of elements 31 @param[in] block_size Number of elements in a block 32 @param[in] elem_size Size of each element 33 34 @return An error code: 0 - success, otherwise - failure 35 36 @ref Utility 37 **/ 38 int CeedPermutePadOffsets(const CeedInt *offsets, CeedInt *block_offsets, CeedInt num_block, CeedInt num_elem, CeedInt block_size, 39 CeedInt elem_size) { 40 for (CeedInt e = 0; e < num_block * block_size; e += block_size) { 41 for (CeedInt j = 0; j < block_size; j++) { 42 for (CeedInt k = 0; k < elem_size; k++) { 43 block_offsets[e * elem_size + k * block_size + j] = offsets[CeedIntMin(e + j, num_elem - 1) * elem_size + k]; 44 } 45 } 46 } 47 return CEED_ERROR_SUCCESS; 48 } 49 50 /** 51 @brief Permute and pad orientations for a blocked `CeedElemRestriction` 52 53 @param[in] orients Array of shape `[num_elem, elem_size]` 54 @param[out] block_orients Array of permuted and padded array values of shape `[num_block, elem_size, block_size]` 55 @param[in] num_block Number of blocks 56 @param[in] num_elem Number of elements 57 @param[in] block_size Number of elements in a block 58 @param[in] elem_size Size of each element 59 60 @return An error code: 0 - success, otherwise - failure 61 62 @ref Utility 63 **/ 64 int CeedPermutePadOrients(const bool *orients, bool *block_orients, CeedInt num_block, CeedInt num_elem, CeedInt block_size, CeedInt elem_size) { 65 for (CeedInt e = 0; e < num_block * block_size; e += block_size) { 66 for (CeedInt j = 0; j < block_size; j++) { 67 for (CeedInt k = 0; k < elem_size; k++) { 68 block_orients[e * elem_size + k * block_size + j] = orients[CeedIntMin(e + j, num_elem - 1) * elem_size + k]; 69 } 70 } 71 } 72 return CEED_ERROR_SUCCESS; 73 } 74 75 /** 76 @brief Permute and pad curl-conforming orientations for a blocked `CeedElemRestriction` 77 78 @param[in] curl_orients Array of shape `[num_elem, elem_size]` 79 @param[out] block_curl_orients Array of permuted and padded array values of shape `[num_block, elem_size, block_size]` 80 @param[in] num_block Number of blocks 81 @param[in] num_elem Number of elements 82 @param[in] block_size Number of elements in a block 83 @param[in] elem_size Size of each element 84 85 @return An error code: 0 - success, otherwise - failure 86 87 @ref Utility 88 **/ 89 int CeedPermutePadCurlOrients(const CeedInt8 *curl_orients, CeedInt8 *block_curl_orients, CeedInt num_block, CeedInt num_elem, CeedInt block_size, 90 CeedInt elem_size) { 91 for (CeedInt e = 0; e < num_block * block_size; e += block_size) { 92 for (CeedInt j = 0; j < block_size; j++) { 93 for (CeedInt k = 0; k < elem_size; k++) { 94 block_curl_orients[e * elem_size + k * block_size + j] = curl_orients[CeedIntMin(e + j, num_elem - 1) * elem_size + k]; 95 } 96 } 97 } 98 return CEED_ERROR_SUCCESS; 99 } 100 101 /** 102 @brief View a `CeedElemRestriction` passed as a `CeedObject` 103 104 @param[in] rstr `CeedElemRestriction` to view 105 @param[in] stream Filestream to write to 106 107 @return An error code: 0 - success, otherwise - failure 108 109 @ref Developer 110 **/ 111 static int CeedElemRestrictionView_Object(CeedObject rstr, FILE *stream) { 112 CeedCall(CeedElemRestrictionView((CeedElemRestriction)rstr, stream)); 113 return CEED_ERROR_SUCCESS; 114 } 115 116 /** 117 @brief Destroy a `CeedElemRestricton` passed as a `CeedObject` 118 119 @param[in,out] rstr Address of `CeedElemRestriction` to destroy 120 121 @return An error code: 0 - success, otherwise - failure 122 123 @ref Developer 124 **/ 125 static int CeedElemRestrictionDestroy_Object(CeedObject *rstr) { 126 CeedCall(CeedElemRestrictionDestroy((CeedElemRestriction *)rstr)); 127 return CEED_ERROR_SUCCESS; 128 } 129 130 /// @} 131 132 /// ---------------------------------------------------------------------------- 133 /// CeedElemRestriction Backend API 134 /// ---------------------------------------------------------------------------- 135 /// @addtogroup CeedElemRestrictionBackend 136 /// @{ 137 138 /** 139 @brief Get the type of a `CeedElemRestriction` 140 141 @param[in] rstr `CeedElemRestriction` 142 @param[out] rstr_type Variable to store restriction type 143 144 @return An error code: 0 - success, otherwise - failure 145 146 @ref Backend 147 **/ 148 int CeedElemRestrictionGetType(CeedElemRestriction rstr, CeedRestrictionType *rstr_type) { 149 *rstr_type = rstr->rstr_type; 150 return CEED_ERROR_SUCCESS; 151 } 152 153 /** 154 @brief Get the strided status of a `CeedElemRestriction` 155 156 @param[in] rstr `CeedElemRestriction` 157 @param[out] is_strided Variable to store strided status 158 159 @return An error code: 0 - success, otherwise - failure 160 161 @ref Backend 162 **/ 163 int CeedElemRestrictionIsStrided(CeedElemRestriction rstr, bool *is_strided) { 164 *is_strided = (rstr->rstr_type == CEED_RESTRICTION_STRIDED); 165 return CEED_ERROR_SUCCESS; 166 } 167 168 /** 169 @brief Get the points status of a `CeedElemRestriction` 170 171 @param[in] rstr `CeedElemRestriction` 172 @param[out] is_points Variable to store points status 173 174 @return An error code: 0 - success, otherwise - failure 175 176 @ref Backend 177 **/ 178 int CeedElemRestrictionIsAtPoints(CeedElemRestriction rstr, bool *is_points) { 179 *is_points = (rstr->rstr_type == CEED_RESTRICTION_POINTS); 180 return CEED_ERROR_SUCCESS; 181 } 182 183 /** 184 @brief Check if two `CeedElemRestriction` created with @ref CeedElemRestrictionCreateAtPoints() and use the same points per element 185 186 @param[in] rstr_a First `CeedElemRestriction` 187 @param[in] rstr_b Second `CeedElemRestriction` 188 @param[out] are_compatible Variable to store compatibility status 189 190 @return An error code: 0 - success, otherwise - failure 191 192 @ref Backend 193 **/ 194 int CeedElemRestrictionAtPointsAreCompatible(CeedElemRestriction rstr_a, CeedElemRestriction rstr_b, bool *are_compatible) { 195 CeedInt num_elem_a, num_elem_b, num_points_a, num_points_b; 196 197 // Cannot compare non-points restrictions 198 CeedCheck(rstr_a->rstr_type == CEED_RESTRICTION_POINTS, CeedElemRestrictionReturnCeed(rstr_a), CEED_ERROR_UNSUPPORTED, 199 "First CeedElemRestriction must be AtPoints"); 200 CeedCheck(rstr_b->rstr_type == CEED_RESTRICTION_POINTS, CeedElemRestrictionReturnCeed(rstr_a), CEED_ERROR_UNSUPPORTED, 201 "Second CeedElemRestriction must be AtPoints"); 202 203 CeedCall(CeedElemRestrictionGetNumElements(rstr_a, &num_elem_a)); 204 CeedCall(CeedElemRestrictionGetNumElements(rstr_b, &num_elem_b)); 205 CeedCall(CeedElemRestrictionGetNumPoints(rstr_a, &num_points_a)); 206 CeedCall(CeedElemRestrictionGetNumPoints(rstr_b, &num_points_b)); 207 208 // Check size and contents of offsets arrays 209 *are_compatible = true; 210 if (num_elem_a != num_elem_b) *are_compatible = false; 211 if (num_points_a != num_points_b) *are_compatible = false; 212 if (*are_compatible) { 213 const CeedInt *offsets_a, *offsets_b; 214 215 CeedCall(CeedElemRestrictionGetOffsets(rstr_a, CEED_MEM_HOST, &offsets_a)); 216 CeedCall(CeedElemRestrictionGetOffsets(rstr_b, CEED_MEM_HOST, &offsets_b)); 217 for (CeedInt i = 0; i < num_elem_a + 1 + num_points_a; i++) *are_compatible &= offsets_a[i] == offsets_b[i]; 218 CeedCall(CeedElemRestrictionRestoreOffsets(rstr_a, &offsets_a)); 219 CeedCall(CeedElemRestrictionRestoreOffsets(rstr_b, &offsets_b)); 220 } 221 return CEED_ERROR_SUCCESS; 222 } 223 224 /** 225 @brief Get the strides of a strided `CeedElemRestriction` 226 227 @param[in] rstr `CeedElemRestriction` 228 @param[out] strides Variable to store strides array 229 230 @return An error code: 0 - success, otherwise - failure 231 232 @ref Backend 233 **/ 234 int CeedElemRestrictionGetStrides(CeedElemRestriction rstr, CeedInt strides[3]) { 235 CeedCheck(rstr->strides, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_MINOR, "CeedElemRestriction has no stride data"); 236 for (CeedInt i = 0; i < 3; i++) strides[i] = rstr->strides[i]; 237 return CEED_ERROR_SUCCESS; 238 } 239 240 /** 241 @brief Get the backend stride status of a `CeedElemRestriction` 242 243 @param[in] rstr `CeedElemRestriction` 244 @param[out] has_backend_strides Variable to store stride status 245 246 @return An error code: 0 - success, otherwise - failure 247 248 @ref Backend 249 **/ 250 int CeedElemRestrictionHasBackendStrides(CeedElemRestriction rstr, bool *has_backend_strides) { 251 CeedCheck(rstr->strides, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_MINOR, "CeedElemRestriction has no stride data"); 252 *has_backend_strides = ((rstr->strides[0] == CEED_STRIDES_BACKEND[0]) && (rstr->strides[1] == CEED_STRIDES_BACKEND[1]) && 253 (rstr->strides[2] == CEED_STRIDES_BACKEND[2])); 254 return CEED_ERROR_SUCCESS; 255 } 256 257 /** 258 @brief Get read-only access to a `CeedElemRestriction` offsets array by @ref CeedMemType 259 260 @param[in] rstr `CeedElemRestriction` to retrieve offsets 261 @param[in] mem_type Memory type on which to access the array. 262 If the backend uses a different memory type, this will perform a copy (possibly cached). 263 @param[out] offsets Array on memory type `mem_type` 264 265 @return An error code: 0 - success, otherwise - failure 266 267 @ref User 268 **/ 269 int CeedElemRestrictionGetOffsets(CeedElemRestriction rstr, CeedMemType mem_type, const CeedInt **offsets) { 270 if (rstr->rstr_base) { 271 CeedCall(CeedElemRestrictionGetOffsets(rstr->rstr_base, mem_type, offsets)); 272 } else { 273 CeedCheck(rstr->GetOffsets, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_UNSUPPORTED, 274 "Backend does not implement CeedElemRestrictionGetOffsets"); 275 CeedCall(rstr->GetOffsets(rstr, mem_type, offsets)); 276 rstr->num_readers++; 277 } 278 return CEED_ERROR_SUCCESS; 279 } 280 281 /** 282 @brief Restore an offsets array obtained using @ref CeedElemRestrictionGetOffsets() 283 284 @param[in] rstr `CeedElemRestriction` to restore 285 @param[in] offsets Array of offset data 286 287 @return An error code: 0 - success, otherwise - failure 288 289 @ref User 290 **/ 291 int CeedElemRestrictionRestoreOffsets(CeedElemRestriction rstr, const CeedInt **offsets) { 292 if (rstr->rstr_base) { 293 CeedCall(CeedElemRestrictionRestoreOffsets(rstr->rstr_base, offsets)); 294 } else { 295 *offsets = NULL; 296 rstr->num_readers--; 297 } 298 return CEED_ERROR_SUCCESS; 299 } 300 301 /** 302 @brief Get read-only access to a `CeedElemRestriction` orientations array by @ref CeedMemType 303 304 @param[in] rstr `CeedElemRestriction` to retrieve orientations 305 @param[in] mem_type Memory type on which to access the array. 306 If the backend uses a different memory type, this will perform a copy (possibly cached). 307 @param[out] orients Array on memory type `mem_type` 308 309 @return An error code: 0 - success, otherwise - failure 310 311 @ref User 312 **/ 313 int CeedElemRestrictionGetOrientations(CeedElemRestriction rstr, CeedMemType mem_type, const bool **orients) { 314 CeedCheck(rstr->GetOrientations, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_UNSUPPORTED, 315 "Backend does not implement CeedElemRestrictionGetOrientations"); 316 CeedCall(rstr->GetOrientations(rstr, mem_type, orients)); 317 rstr->num_readers++; 318 return CEED_ERROR_SUCCESS; 319 } 320 321 /** 322 @brief Restore an orientations array obtained using @ref CeedElemRestrictionGetOrientations() 323 324 @param[in] rstr `CeedElemRestriction` to restore 325 @param[in] orients Array of orientation data 326 327 @return An error code: 0 - success, otherwise - failure 328 329 @ref User 330 **/ 331 int CeedElemRestrictionRestoreOrientations(CeedElemRestriction rstr, const bool **orients) { 332 *orients = NULL; 333 rstr->num_readers--; 334 return CEED_ERROR_SUCCESS; 335 } 336 337 /** 338 @brief Get read-only access to a `CeedElemRestriction` curl-conforming orientations array by @ref CeedMemType 339 340 @param[in] rstr `CeedElemRestriction` to retrieve curl-conforming orientations 341 @param[in] mem_type Memory type on which to access the array. 342 If the backend uses a different memory type, this will perform a copy (possibly cached). 343 @param[out] curl_orients Array on memory type `mem_type` 344 345 @return An error code: 0 - success, otherwise - failure 346 347 @ref User 348 **/ 349 int CeedElemRestrictionGetCurlOrientations(CeedElemRestriction rstr, CeedMemType mem_type, const CeedInt8 **curl_orients) { 350 CeedCheck(rstr->GetCurlOrientations, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_UNSUPPORTED, 351 "Backend does not implement CeedElemRestrictionGetCurlOrientations"); 352 CeedCall(rstr->GetCurlOrientations(rstr, mem_type, curl_orients)); 353 rstr->num_readers++; 354 return CEED_ERROR_SUCCESS; 355 } 356 357 /** 358 @brief Restore an orientations array obtained using @ref CeedElemRestrictionGetCurlOrientations() 359 360 @param[in] rstr `CeedElemRestriction` to restore 361 @param[in] curl_orients Array of orientation data 362 363 @return An error code: 0 - success, otherwise - failure 364 365 @ref User 366 **/ 367 int CeedElemRestrictionRestoreCurlOrientations(CeedElemRestriction rstr, const CeedInt8 **curl_orients) { 368 *curl_orients = NULL; 369 rstr->num_readers--; 370 return CEED_ERROR_SUCCESS; 371 } 372 373 /** 374 375 @brief Get the L-vector layout of a strided `CeedElemRestriction` 376 377 @param[in] rstr `CeedElemRestriction` 378 @param[out] layout Variable to store layout array, stored as `[nodes, components, elements]`. 379 The data for node `i`, component `j`, element `k` in the E-vector is given by `i*layout[0] + j*layout[1] + k*layout[2]`. 380 381 @return An error code: 0 - success, otherwise - failure 382 383 @ref Backend 384 **/ 385 int CeedElemRestrictionGetLLayout(CeedElemRestriction rstr, CeedInt layout[3]) { 386 bool has_backend_strides; 387 CeedRestrictionType rstr_type; 388 389 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 390 CeedCheck(rstr_type == CEED_RESTRICTION_STRIDED, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_MINOR, 391 "Only strided CeedElemRestriction have strided L-vector layout"); 392 CeedCall(CeedElemRestrictionHasBackendStrides(rstr, &has_backend_strides)); 393 if (has_backend_strides) { 394 CeedCheck(rstr->l_layout[0], CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_MINOR, "CeedElemRestriction has no L-vector layout data"); 395 for (CeedInt i = 0; i < 3; i++) layout[i] = rstr->l_layout[i]; 396 } else { 397 CeedCall(CeedElemRestrictionGetStrides(rstr, layout)); 398 } 399 return CEED_ERROR_SUCCESS; 400 } 401 402 /** 403 404 @brief Set the L-vector layout of a strided `CeedElemRestriction` 405 406 @param[in] rstr `CeedElemRestriction` 407 @param[in] layout Variable to containing layout array, stored as `[nodes, components, elements]`. 408 The data for node `i`, component `j`, element `k` in the E-vector is given by `i*layout[0] + j*layout[1] + k*layout[2]`. 409 410 @return An error code: 0 - success, otherwise - failure 411 412 @ref Backend 413 **/ 414 int CeedElemRestrictionSetLLayout(CeedElemRestriction rstr, CeedInt layout[3]) { 415 CeedRestrictionType rstr_type; 416 417 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 418 CeedCheck(rstr_type == CEED_RESTRICTION_STRIDED, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_MINOR, 419 "Only strided CeedElemRestriction have strided L-vector layout"); 420 for (CeedInt i = 0; i < 3; i++) rstr->l_layout[i] = layout[i]; 421 return CEED_ERROR_SUCCESS; 422 } 423 424 /** 425 426 @brief Get the E-vector layout of a `CeedElemRestriction` 427 428 @param[in] rstr `CeedElemRestriction` 429 @param[out] layout Variable to store layout array, stored as `[nodes, components, elements]`. 430 The data for node `i`, component `j`, element `k` in the E-vector is given by `i*layout[0] + j*layout[1] + k*layout[2]`. 431 432 @return An error code: 0 - success, otherwise - failure 433 434 @ref Backend 435 **/ 436 int CeedElemRestrictionGetELayout(CeedElemRestriction rstr, CeedInt layout[3]) { 437 CeedCheck(rstr->e_layout[0], CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_MINOR, "CeedElemRestriction has no E-vector layout data"); 438 for (CeedInt i = 0; i < 3; i++) layout[i] = rstr->e_layout[i]; 439 return CEED_ERROR_SUCCESS; 440 } 441 442 /** 443 444 @brief Set the E-vector layout of a `CeedElemRestriction` 445 446 @param[in] rstr `CeedElemRestriction` 447 @param[in] layout Variable to containing layout array, stored as `[nodes, components, elements]`. 448 The data for node `i`, component `j`, element `k` in the E-vector is given by `i*layout[0] + j*layout[1] + k*layout[2]`. 449 450 @return An error code: 0 - success, otherwise - failure 451 452 @ref Backend 453 **/ 454 int CeedElemRestrictionSetELayout(CeedElemRestriction rstr, CeedInt layout[3]) { 455 for (CeedInt i = 0; i < 3; i++) rstr->e_layout[i] = layout[i]; 456 return CEED_ERROR_SUCCESS; 457 } 458 459 /** 460 461 @brief Get the E-vector element offset of a `CeedElemRestriction` at points 462 463 @param[in] rstr `CeedElemRestriction` 464 @param[in] elem Element number index into E-vector for 465 @param[out] elem_offset Offset for element `elem` in the E-vector. 466 The data for point `i`, component `j`, element `elem` in the E-vector is given by `i*e_layout[0] + j*e_layout[1] + elem_offset`. 467 468 @return An error code: 0 - success, otherwise - failure 469 470 @ref Backend 471 **/ 472 int CeedElemRestrictionGetAtPointsElementOffset(CeedElemRestriction rstr, CeedInt elem, CeedSize *elem_offset) { 473 CeedInt num_comp; 474 CeedRestrictionType rstr_type; 475 476 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 477 CeedCheck(rstr_type == CEED_RESTRICTION_POINTS, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_INCOMPATIBLE, 478 "Can only compute offset for a points CeedElemRestriction"); 479 480 // Backend method 481 if (rstr->GetAtPointsElementOffset) { 482 CeedCall(rstr->GetAtPointsElementOffset(rstr, elem, elem_offset)); 483 return CEED_ERROR_SUCCESS; 484 } 485 486 // Default layout (CPU) 487 *elem_offset = 0; 488 CeedCall(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 489 for (CeedInt i = 0; i < elem; i++) { 490 CeedInt num_points; 491 492 CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr, i, &num_points)); 493 *elem_offset += num_points * num_comp; 494 } 495 return CEED_ERROR_SUCCESS; 496 } 497 498 /** 499 500 @brief Set the E-vector size of a `CeedElemRestriction` at points 501 502 @param[in,out] rstr `CeedElemRestriction` 503 @param[in] e_size New E-vector size; must be longer than the current E-vector size 504 505 @return An error code: 0 - success, otherwise - failure 506 507 @ref Backend 508 **/ 509 int CeedElemRestrictionSetAtPointsEVectorSize(CeedElemRestriction rstr, CeedSize e_size) { 510 CeedRestrictionType rstr_type; 511 512 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 513 CeedCheck(rstr_type == CEED_RESTRICTION_POINTS, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_INCOMPATIBLE, 514 "Can only compute offset for a points CeedElemRestriction"); 515 CeedCheck(e_size >= rstr->e_size, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_INCOMPATIBLE, 516 "Can only increase the size of the E-vector for the CeedElemRestriction." 517 " Current size: %" CeedSize_FMT " New size: %" CeedSize_FMT, 518 rstr->e_size, e_size); 519 rstr->e_size = e_size; 520 return CEED_ERROR_SUCCESS; 521 } 522 523 /** 524 @brief Get the backend data of a `CeedElemRestriction` 525 526 @param[in] rstr `CeedElemRestriction` 527 @param[out] data Variable to store data 528 529 @return An error code: 0 - success, otherwise - failure 530 531 @ref Backend 532 **/ 533 int CeedElemRestrictionGetData(CeedElemRestriction rstr, void *data) { 534 *(void **)data = rstr->data; 535 return CEED_ERROR_SUCCESS; 536 } 537 538 /** 539 @brief Set the backend data of a `CeedElemRestriction` 540 541 @param[in,out] rstr `CeedElemRestriction` 542 @param[in] data Data to set 543 544 @return An error code: 0 - success, otherwise - failure 545 546 @ref Backend 547 **/ 548 int CeedElemRestrictionSetData(CeedElemRestriction rstr, void *data) { 549 rstr->data = data; 550 return CEED_ERROR_SUCCESS; 551 } 552 553 /** 554 @brief Increment the reference counter for a `CeedElemRestriction` 555 556 @param[in,out] rstr `CeedElemRestriction` to increment the reference counter 557 558 @return An error code: 0 - success, otherwise - failure 559 560 @ref Backend 561 **/ 562 int CeedElemRestrictionReference(CeedElemRestriction rstr) { 563 CeedCall(CeedObjectReference((CeedObject)rstr)); 564 return CEED_ERROR_SUCCESS; 565 } 566 567 /** 568 @brief Estimate number of FLOPs required to apply `CeedElemRestriction` in `t_mode` 569 570 @param[in] rstr `CeedElemRestriction` to estimate FLOPs for 571 @param[in] t_mode Apply restriction or transpose 572 @param[out] flops Address of variable to hold FLOPs estimate 573 574 @ref Backend 575 **/ 576 int CeedElemRestrictionGetFlopsEstimate(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedSize *flops) { 577 CeedSize e_size, scale = 0; 578 CeedRestrictionType rstr_type; 579 580 CeedCall(CeedElemRestrictionGetEVectorSize(rstr, &e_size)); 581 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 582 if (t_mode == CEED_TRANSPOSE) { 583 switch (rstr_type) { 584 case CEED_RESTRICTION_POINTS: 585 scale = 0; 586 break; 587 case CEED_RESTRICTION_STRIDED: 588 case CEED_RESTRICTION_STANDARD: 589 scale = 1; 590 break; 591 case CEED_RESTRICTION_ORIENTED: 592 scale = 2; 593 break; 594 case CEED_RESTRICTION_CURL_ORIENTED: 595 scale = 6; 596 break; 597 } 598 } else { 599 switch (rstr_type) { 600 case CEED_RESTRICTION_STRIDED: 601 case CEED_RESTRICTION_STANDARD: 602 case CEED_RESTRICTION_POINTS: 603 scale = 0; 604 break; 605 case CEED_RESTRICTION_ORIENTED: 606 scale = 1; 607 break; 608 case CEED_RESTRICTION_CURL_ORIENTED: 609 scale = 5; 610 break; 611 } 612 } 613 *flops = e_size * scale; 614 return CEED_ERROR_SUCCESS; 615 } 616 617 /// @} 618 619 /// @cond DOXYGEN_SKIP 620 static struct CeedElemRestriction_private ceed_elemrestriction_none; 621 /// @endcond 622 623 /// ---------------------------------------------------------------------------- 624 /// CeedElemRestriction Public API 625 /// ---------------------------------------------------------------------------- 626 /// @addtogroup CeedElemRestrictionUser 627 /// @{ 628 629 /// Indicate that the stride is determined by the backend 630 const CeedInt CEED_STRIDES_BACKEND[3] = {0}; 631 632 /// Argument for @ref CeedOperatorSetField() indicating that the field does not require a `CeedElemRestriction` 633 const CeedElemRestriction CEED_ELEMRESTRICTION_NONE = &ceed_elemrestriction_none; 634 635 /** 636 @brief Create a `CeedElemRestriction` 637 638 @param[in] ceed `Ceed` context used to create the `CeedElemRestriction` 639 @param[in] num_elem Number of elements described in the `offsets` array 640 @param[in] elem_size Size (number of "nodes") per element 641 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 642 @param[in] comp_stride Stride between components for the same L-vector "node". 643 Data for node `i`, component `j`, element `k` can be found in the L-vector at index `offsets[i + k*elem_size] + j*comp_stride`. 644 @param[in] l_size The size of the L-vector. 645 This vector may be larger than the elements and fields given by this restriction. 646 @param[in] mem_type Memory type of the `offsets` array, see @ref CeedMemType 647 @param[in] copy_mode Copy mode for the `offsets` array, see @ref CeedCopyMode 648 @param[in] offsets Array of shape `[num_elem, elem_size]`. 649 Row `i` holds the ordered list of the offsets (into the input `CeedVector`) for the unknowns corresponding to element `i`, where 0 <= i < @a num_elem. 650 All offsets must be in the range `[0, l_size - 1]`. 651 @param[out] rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 652 653 @return An error code: 0 - success, otherwise - failure 654 655 @ref User 656 **/ 657 int CeedElemRestrictionCreate(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt num_comp, CeedInt comp_stride, CeedSize l_size, 658 CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, CeedElemRestriction *rstr) { 659 if (!ceed->ElemRestrictionCreate) { 660 Ceed delegate; 661 662 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 663 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedElemRestrictionCreate"); 664 CeedCall(CeedElemRestrictionCreate(delegate, num_elem, elem_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, rstr)); 665 CeedCall(CeedDestroy(&delegate)); 666 return CEED_ERROR_SUCCESS; 667 } 668 669 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 670 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 671 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction must have at least 1 component"); 672 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction component stride must be at least 1"); 673 674 CeedCall(CeedCalloc(1, rstr)); 675 CeedCall(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, CeedElemRestrictionDestroy_Object, &(*rstr)->obj)); 676 (*rstr)->num_elem = num_elem; 677 (*rstr)->elem_size = elem_size; 678 (*rstr)->num_comp = num_comp; 679 (*rstr)->comp_stride = comp_stride; 680 (*rstr)->l_size = l_size; 681 (*rstr)->e_size = (CeedSize)num_elem * (CeedSize)elem_size * (CeedSize)num_comp; 682 (*rstr)->num_block = num_elem; 683 (*rstr)->block_size = 1; 684 (*rstr)->rstr_type = CEED_RESTRICTION_STANDARD; 685 CeedCall(ceed->ElemRestrictionCreate(mem_type, copy_mode, offsets, NULL, NULL, *rstr)); 686 return CEED_ERROR_SUCCESS; 687 } 688 689 /** 690 @brief Create a `CeedElemRestriction` with orientation signs 691 692 @param[in] ceed `Ceed` context used to create the `CeedElemRestriction` 693 @param[in] num_elem Number of elements described in the `offsets` array 694 @param[in] elem_size Size (number of "nodes") per element 695 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 696 @param[in] comp_stride Stride between components for the same L-vector "node". 697 Data for node `i`, component `j`, element `k` can be found in the L-vector at index `offsets[i + k*elem_size] + j*comp_stride`. 698 @param[in] l_size The size of the L-vector. 699 This vector may be larger than the elements and fields given by this restriction. 700 @param[in] mem_type Memory type of the `offsets` array, see @ref CeedMemType 701 @param[in] copy_mode Copy mode for the `offsets` array, see @ref CeedCopyMode 702 @param[in] offsets Array of shape `[num_elem, elem_size]`. 703 Row i holds the ordered list of the offsets (into the input `CeedVector`) for the unknowns corresponding to element `i`, where `0 <= i < num_elem`. 704 All offsets must be in the range `[0, l_size - 1]`. 705 @param[in] orients Boolean array of shape `[num_elem, elem_size]` with `false` for positively oriented and `true` to flip the orientation 706 @param[out] rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 707 708 @return An error code: 0 - success, otherwise - failure 709 710 @ref User 711 **/ 712 int CeedElemRestrictionCreateOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt num_comp, CeedInt comp_stride, CeedSize l_size, 713 CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, const bool *orients, 714 CeedElemRestriction *rstr) { 715 if (!ceed->ElemRestrictionCreate) { 716 Ceed delegate; 717 718 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 719 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedElemRestrictionCreateOriented"); 720 CeedCall(CeedElemRestrictionCreateOriented(delegate, num_elem, elem_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, orients, 721 rstr)); 722 CeedCall(CeedDestroy(&delegate)); 723 return CEED_ERROR_SUCCESS; 724 } 725 726 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 727 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 728 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction must have at least 1 component"); 729 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction component stride must be at least 1"); 730 731 CeedCall(CeedCalloc(1, rstr)); 732 CeedCall(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, CeedElemRestrictionDestroy_Object, &(*rstr)->obj)); 733 (*rstr)->num_elem = num_elem; 734 (*rstr)->elem_size = elem_size; 735 (*rstr)->num_comp = num_comp; 736 (*rstr)->comp_stride = comp_stride; 737 (*rstr)->l_size = l_size; 738 (*rstr)->e_size = (CeedSize)num_elem * (CeedSize)elem_size * (CeedSize)num_comp; 739 (*rstr)->num_block = num_elem; 740 (*rstr)->block_size = 1; 741 (*rstr)->rstr_type = CEED_RESTRICTION_ORIENTED; 742 CeedCall(ceed->ElemRestrictionCreate(mem_type, copy_mode, offsets, orients, NULL, *rstr)); 743 return CEED_ERROR_SUCCESS; 744 } 745 746 /** 747 @brief Create a `CeedElemRestriction` with a general tridiagonal transformation matrix for curl-conforming elements 748 749 @param[in] ceed `Ceed` context used to create the `CeedElemRestriction` 750 @param[in] num_elem Number of elements described in the `offsets` array 751 @param[in] elem_size Size (number of "nodes") per element 752 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 753 @param[in] comp_stride Stride between components for the same L-vector "node". 754 Data for node `i`, component `j`, element `k` can be found in the L-vector at index `offsets[i + k*elem_size] + j*comp_stride`. 755 @param[in] l_size The size of the L-vector. 756 This vector may be larger than the elements and fields given by this restriction. 757 @param[in] mem_type Memory type of the `offsets` array, see @ref CeedMemType 758 @param[in] copy_mode Copy mode for the `offsets` array, see @ref CeedCopyMode 759 @param[in] offsets Array of shape `[num_elem, elem_size]`. 760 Row `i` holds the ordered list of the offsets (into the input `CeedVector`) for the unknowns corresponding to element `i`, where `0 <= i < num_elem`. 761 All offsets must be in the range `[0, l_size - 1]`. 762 @param[in] curl_orients Array of shape `[num_elem, 3 * elem_size]` representing a row-major tridiagonal matrix (`curl_orients[i * 3 * elem_size] = curl_orients[(i + 1) * 3 * elem_size - 1] = 0`, where `0 <= i < num_elem`) which is applied to the element unknowns upon restriction. 763 This orientation matrix allows for pairs of face degrees of freedom on elements for \f$H(\mathrm{curl})\f$ spaces to be coupled in the element restriction operation, which is a way to resolve face orientation issues for 3D meshes (https://dl.acm.org/doi/pdf/10.1145/3524456). 764 @param[out] rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 765 766 @return An error code: 0 - success, otherwise - failure 767 768 @ref User 769 **/ 770 int CeedElemRestrictionCreateCurlOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt num_comp, CeedInt comp_stride, CeedSize l_size, 771 CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, const CeedInt8 *curl_orients, 772 CeedElemRestriction *rstr) { 773 if (!ceed->ElemRestrictionCreate) { 774 Ceed delegate; 775 776 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 777 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedElemRestrictionCreateCurlOriented"); 778 CeedCall(CeedElemRestrictionCreateCurlOriented(delegate, num_elem, elem_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, 779 curl_orients, rstr)); 780 CeedCall(CeedDestroy(&delegate)); 781 return CEED_ERROR_SUCCESS; 782 } 783 784 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 785 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 786 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction must have at least 1 component"); 787 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction component stride must be at least 1"); 788 789 CeedCall(CeedCalloc(1, rstr)); 790 CeedCall(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, CeedElemRestrictionDestroy_Object, &(*rstr)->obj)); 791 (*rstr)->num_elem = num_elem; 792 (*rstr)->elem_size = elem_size; 793 (*rstr)->num_comp = num_comp; 794 (*rstr)->comp_stride = comp_stride; 795 (*rstr)->l_size = l_size; 796 (*rstr)->e_size = (CeedSize)num_elem * (CeedSize)elem_size * (CeedSize)num_comp; 797 (*rstr)->num_block = num_elem; 798 (*rstr)->block_size = 1; 799 (*rstr)->rstr_type = CEED_RESTRICTION_CURL_ORIENTED; 800 CeedCall(ceed->ElemRestrictionCreate(mem_type, copy_mode, offsets, NULL, curl_orients, *rstr)); 801 return CEED_ERROR_SUCCESS; 802 } 803 804 /** 805 @brief Create a strided `CeedElemRestriction` 806 807 @param[in] ceed `Ceed` context used to create the `CeedElemRestriction` 808 @param[in] num_elem Number of elements described by the restriction 809 @param[in] elem_size Size (number of "nodes") per element 810 @param[in] num_comp Number of field components per interpolation "node" (1 for scalar fields) 811 @param[in] l_size The size of the L-vector. 812 This vector may be larger than the elements and fields given by this restriction. 813 @param[in] strides Array for strides between `[nodes, components, elements]`. 814 Data for node `i`, component `j`, element `k` can be found in the L-vector at index `i*strides[0] + j*strides[1] + k*strides[2]`. 815 @ref CEED_STRIDES_BACKEND may be used for `CeedVector` ordered by the same `Ceed` backend. 816 `CEED_STRIDES_BACKEND` should only be used pass data between `CeedOperator` created with the same `Ceed` backend. 817 The L-vector layout will, in general, be different between `Ceed` backends. 818 @param[out] rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 819 820 @return An error code: 0 - success, otherwise - failure 821 822 @ref User 823 **/ 824 int CeedElemRestrictionCreateStrided(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt num_comp, CeedSize l_size, const CeedInt strides[3], 825 CeedElemRestriction *rstr) { 826 if (!ceed->ElemRestrictionCreate) { 827 Ceed delegate; 828 829 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 830 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedElemRestrictionCreateStrided"); 831 CeedCall(CeedElemRestrictionCreateStrided(delegate, num_elem, elem_size, num_comp, l_size, strides, rstr)); 832 CeedCall(CeedDestroy(&delegate)); 833 return CEED_ERROR_SUCCESS; 834 } 835 836 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 837 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 838 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction must have at least 1 component"); 839 CeedCheck(l_size >= (CeedSize)num_elem * (CeedSize)elem_size * (CeedSize)num_comp, ceed, CEED_ERROR_DIMENSION, 840 "L-vector size must be at least num_elem * elem_size * num_comp. Expected: > %" CeedSize_FMT " Found: %" CeedSize_FMT, 841 (CeedSize)num_elem * (CeedSize)elem_size * (CeedSize)num_comp, l_size); 842 843 CeedCall(CeedCalloc(1, rstr)); 844 CeedCall(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, CeedElemRestrictionDestroy_Object, &(*rstr)->obj)); 845 (*rstr)->num_elem = num_elem; 846 (*rstr)->elem_size = elem_size; 847 (*rstr)->num_comp = num_comp; 848 (*rstr)->l_size = l_size; 849 (*rstr)->e_size = (CeedSize)num_elem * (CeedSize)elem_size * (CeedSize)num_comp; 850 (*rstr)->num_block = num_elem; 851 (*rstr)->block_size = 1; 852 (*rstr)->rstr_type = CEED_RESTRICTION_STRIDED; 853 CeedCall(CeedMalloc(3, &(*rstr)->strides)); 854 for (CeedInt i = 0; i < 3; i++) (*rstr)->strides[i] = strides[i]; 855 CeedCall(ceed->ElemRestrictionCreate(CEED_MEM_HOST, CEED_OWN_POINTER, NULL, NULL, NULL, *rstr)); 856 return CEED_ERROR_SUCCESS; 857 } 858 859 /** 860 @brief Create a points `CeedElemRestriction`, for restricting for restricting from a all local points to the current element in which they are located. 861 862 The offsets array is arranged as 863 864 element_0_start_index 865 element_1_start_index 866 ... 867 element_n_start_index 868 element_n_stop_index 869 element_0_point_0 870 element_0_point_1 871 ... 872 873 @param[in] ceed `Ceed` context used to create the `CeedElemRestriction` 874 @param[in] num_elem Number of elements described in the `offsets` array 875 @param[in] num_points Number of points described in the `offsets` array 876 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields). 877 Components are assumed to be contiguous by point. 878 @param[in] l_size The size of the L-vector. 879 This vector may be larger than the elements and fields given by this restriction. 880 @param[in] mem_type Memory type of the `offsets` array, see @ref CeedMemType 881 @param[in] copy_mode Copy mode for the `offsets` array, see @ref CeedCopyMode 882 @param[in] offsets Array of size `num_elem + 1 + num_points`. 883 The first portion of the offsets array holds the ranges of indices corresponding to each element. 884 The second portion holds the indices for each element. 885 @param[out] rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 886 887 @return An error code: 0 - success, otherwise - failure 888 889 @ref Backend 890 **/ 891 int CeedElemRestrictionCreateAtPoints(Ceed ceed, CeedInt num_elem, CeedInt num_points, CeedInt num_comp, CeedSize l_size, CeedMemType mem_type, 892 CeedCopyMode copy_mode, const CeedInt *offsets, CeedElemRestriction *rstr) { 893 if (!ceed->ElemRestrictionCreateAtPoints) { 894 Ceed delegate; 895 896 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 897 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedElemRestrictionCreateAtPoints"); 898 CeedCall(CeedElemRestrictionCreateAtPoints(delegate, num_elem, num_points, num_comp, l_size, mem_type, copy_mode, offsets, rstr)); 899 CeedCall(CeedDestroy(&delegate)); 900 return CEED_ERROR_SUCCESS; 901 } 902 903 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 904 CeedCheck(num_points >= 0, ceed, CEED_ERROR_DIMENSION, "Number of points must be non-negative"); 905 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction must have at least 1 component"); 906 CeedCheck(l_size >= (CeedSize)num_points * num_comp, ceed, CEED_ERROR_DIMENSION, 907 "L-vector must be at least num_points * num_comp. Expected: > %" CeedSize_FMT " Found: %" CeedSize_FMT, (CeedSize)num_points * num_comp, 908 l_size); 909 910 CeedCall(CeedCalloc(1, rstr)); 911 CeedCall(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, CeedElemRestrictionDestroy_Object, &(*rstr)->obj)); 912 (*rstr)->num_elem = num_elem; 913 (*rstr)->num_points = num_points; 914 (*rstr)->num_comp = num_comp; 915 (*rstr)->comp_stride = 1; 916 (*rstr)->l_size = l_size; 917 (*rstr)->e_size = (CeedSize)num_points * (CeedSize)num_comp; 918 (*rstr)->num_block = num_elem; 919 (*rstr)->block_size = 1; 920 (*rstr)->rstr_type = CEED_RESTRICTION_POINTS; 921 CeedCall(ceed->ElemRestrictionCreateAtPoints(mem_type, copy_mode, offsets, NULL, NULL, *rstr)); 922 return CEED_ERROR_SUCCESS; 923 } 924 925 /** 926 @brief Create a blocked `CeedElemRestriction`, typically only used by backends 927 928 @param[in] ceed `Ceed` context used to create the `CeedElemRestriction` 929 @param[in] num_elem Number of elements described in the `offsets` array 930 @param[in] elem_size Size (number of unknowns) per element 931 @param[in] block_size Number of elements in a block 932 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 933 @param[in] comp_stride Stride between components for the same L-vector "node". 934 Data for node `i`, component `j`, element `k` can be found in the L-vector at index `offsets[i + k*elem_size] + j*comp_stride`. 935 @param[in] l_size The size of the L-vector. 936 This vector may be larger than the elements and fields given by this restriction. 937 @param[in] mem_type Memory type of the `offsets` array, see @ref CeedMemType 938 @param[in] copy_mode Copy mode for the `offsets` array, see @ref CeedCopyMode 939 @param[in] offsets Array of shape `[num_elem, elem_size]`. 940 Row `i` holds the ordered list of the offsets (into the input `CeedVector`) for the unknowns corresponding to element `i`, where `0 <= i < num_elem`. 941 All offsets must be in the range `[0, l_size - 1]`. 942 The backend will permute and pad this array to the desired ordering for the blocksize, which is typically given by the backend. 943 The default reordering is to interlace elements. 944 @param[out] rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 945 946 @return An error code: 0 - success, otherwise - failure 947 948 @ref Backend 949 **/ 950 int CeedElemRestrictionCreateBlocked(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, CeedInt comp_stride, 951 CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, 952 CeedElemRestriction *rstr) { 953 CeedInt *block_offsets, num_block = (num_elem / block_size) + !!(num_elem % block_size); 954 955 if (!ceed->ElemRestrictionCreateBlocked) { 956 Ceed delegate; 957 958 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 959 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedElemRestrictionCreateBlocked"); 960 CeedCall(CeedElemRestrictionCreateBlocked(delegate, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, 961 rstr)); 962 CeedCall(CeedDestroy(&delegate)); 963 return CEED_ERROR_SUCCESS; 964 } 965 966 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 967 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 968 CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); 969 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction must have at least 1 component"); 970 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction component stride must be at least 1"); 971 972 CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_offsets)); 973 CeedCall(CeedPermutePadOffsets(offsets, block_offsets, num_block, num_elem, block_size, elem_size)); 974 975 CeedCall(CeedCalloc(1, rstr)); 976 CeedCall(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, CeedElemRestrictionDestroy_Object, &(*rstr)->obj)); 977 (*rstr)->num_elem = num_elem; 978 (*rstr)->elem_size = elem_size; 979 (*rstr)->num_comp = num_comp; 980 (*rstr)->comp_stride = comp_stride; 981 (*rstr)->l_size = l_size; 982 (*rstr)->e_size = (CeedSize)num_block * (CeedSize)block_size * (CeedSize)elem_size * (CeedSize)num_comp; 983 (*rstr)->num_block = num_block; 984 (*rstr)->block_size = block_size; 985 (*rstr)->rstr_type = CEED_RESTRICTION_STANDARD; 986 CeedCall(ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, (const CeedInt *)block_offsets, NULL, NULL, *rstr)); 987 if (copy_mode == CEED_OWN_POINTER) CeedCall(CeedFree(&offsets)); 988 return CEED_ERROR_SUCCESS; 989 } 990 991 /** 992 @brief Create a blocked oriented `CeedElemRestriction`, typically only used by backends 993 994 @param[in] ceed `Ceed` context used to create the `CeedElemRestriction` 995 @param[in] num_elem Number of elements described in the `offsets` array. 996 @param[in] elem_size Size (number of unknowns) per element 997 @param[in] block_size Number of elements in a block 998 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 999 @param[in] comp_stride Stride between components for the same L-vector "node". 1000 Data for node `i`, component `j`, element `k` can be found in the L-vector at index `offsets[i + k*elem_size] + j*comp_stride`. 1001 @param[in] l_size The size of the L-vector. 1002 This vector may be larger than the elements and fields given by this restriction. 1003 @param[in] mem_type Memory type of the `offsets` array, see @ref CeedMemType 1004 @param[in] copy_mode Copy mode for the `offsets` array, see @ref CeedCopyMode 1005 @param[in] offsets Array of shape `[num_elem, elem_size]`. 1006 Row `i` holds the ordered list of the offsets (into the input `CeedVector`) for the unknowns corresponding to element `i`, where `0 <= i < num_elem`. 1007 All offsets must be in the range `[0, l_size - 1]`. 1008 The backend will permute and pad this array to the desired ordering for the blocksize, which is typically given by the backend. 1009 The default reordering is to interlace elements. 1010 @param[in] orients Boolean array of shape `[num_elem, elem_size]` with `false` for positively oriented and `true` to flip the orientation. 1011 Will also be permuted and padded similarly to `offsets`. 1012 @param[out] rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 1013 1014 @return An error code: 0 - success, otherwise - failure 1015 1016 @ref Backend 1017 **/ 1018 int CeedElemRestrictionCreateBlockedOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, 1019 CeedInt comp_stride, CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, 1020 const CeedInt *offsets, const bool *orients, CeedElemRestriction *rstr) { 1021 bool *block_orients; 1022 CeedInt *block_offsets, num_block = (num_elem / block_size) + !!(num_elem % block_size); 1023 1024 if (!ceed->ElemRestrictionCreateBlocked) { 1025 Ceed delegate; 1026 1027 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 1028 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedElemRestrictionCreateBlockedOriented"); 1029 CeedCall(CeedElemRestrictionCreateBlockedOriented(delegate, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, mem_type, copy_mode, 1030 offsets, orients, rstr)); 1031 CeedCall(CeedDestroy(&delegate)); 1032 return CEED_ERROR_SUCCESS; 1033 } 1034 1035 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 1036 CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); 1037 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction must have at least 1 component"); 1038 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction component stride must be at least 1"); 1039 1040 CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_offsets)); 1041 CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_orients)); 1042 CeedCall(CeedPermutePadOffsets(offsets, block_offsets, num_block, num_elem, block_size, elem_size)); 1043 CeedCall(CeedPermutePadOrients(orients, block_orients, num_block, num_elem, block_size, elem_size)); 1044 1045 CeedCall(CeedCalloc(1, rstr)); 1046 CeedCall(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, CeedElemRestrictionDestroy_Object, &(*rstr)->obj)); 1047 (*rstr)->num_elem = num_elem; 1048 (*rstr)->elem_size = elem_size; 1049 (*rstr)->num_comp = num_comp; 1050 (*rstr)->comp_stride = comp_stride; 1051 (*rstr)->l_size = l_size; 1052 (*rstr)->e_size = (CeedSize)num_block * (CeedSize)block_size * (CeedSize)elem_size * (CeedSize)num_comp; 1053 (*rstr)->num_block = num_block; 1054 (*rstr)->block_size = block_size; 1055 (*rstr)->rstr_type = CEED_RESTRICTION_ORIENTED; 1056 CeedCall(ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, (const CeedInt *)block_offsets, (const bool *)block_orients, NULL, 1057 *rstr)); 1058 if (copy_mode == CEED_OWN_POINTER) CeedCall(CeedFree(&offsets)); 1059 return CEED_ERROR_SUCCESS; 1060 } 1061 1062 /** 1063 @brief Create a blocked curl-oriented `CeedElemRestriction`, typically only used by backends 1064 1065 @param[in] ceed `Ceed` context used to create the `CeedElemRestriction` 1066 @param[in] num_elem Number of elements described in the `offsets` array. 1067 @param[in] elem_size Size (number of unknowns) per element 1068 @param[in] block_size Number of elements in a block 1069 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 1070 @param[in] comp_stride Stride between components for the same L-vector "node". 1071 Data for node `i`, component `j`, element `k` can be found in the L-vector at index `offsets[i + k*elem_size] + j*comp_stride`. 1072 @param[in] l_size The size of the L-vector. 1073 This vector may be larger than the elements and fields given by this restriction. 1074 @param[in] mem_type Memory type of the `offsets` array, see @ref CeedMemType 1075 @param[in] copy_mode Copy mode for the `offsets` array, see @ref CeedCopyMode 1076 @param[in] offsets Array of shape `[num_elem, elem_size]`. 1077 Row `i` holds the ordered list of the offsets (into the input `CeedVector`) for the unknowns corresponding to element `i`, where `0 <= i < num_elem`. 1078 All offsets must be in the range `[0, l_size - 1]`. 1079 The backend will permute and pad this array to the desired ordering for the blocksize, which is typically given by the backend. 1080 The default reordering is to interlace elements. 1081 @param[in] curl_orients Array of shape `[num_elem, 3 * elem_size]` representing a row-major tridiagonal matrix (`curl_orients[i * 3 * elem_size] = curl_orients[(i + 1) * 3 * elem_size - 1] = 0`, where `0 <= i < num_elem`) which is applied to the element unknowns upon restriction. 1082 This orientation matrix allows for pairs of face degrees of freedom on elements for \f$H(\mathrm{curl})\f$ spaces to be coupled in the element restriction operation, which is a way to resolve face orientation issues for 3D meshes (https://dl.acm.org/doi/pdf/10.1145/3524456). 1083 Will also be permuted and padded similarly to offsets. 1084 @param[out] rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 1085 1086 @return An error code: 0 - success, otherwise - failure 1087 1088 @ref Backend 1089 **/ 1090 int CeedElemRestrictionCreateBlockedCurlOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, 1091 CeedInt comp_stride, CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, 1092 const CeedInt *offsets, const CeedInt8 *curl_orients, CeedElemRestriction *rstr) { 1093 CeedInt8 *block_curl_orients; 1094 CeedInt *block_offsets, num_block = (num_elem / block_size) + !!(num_elem % block_size); 1095 1096 if (!ceed->ElemRestrictionCreateBlocked) { 1097 Ceed delegate; 1098 1099 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 1100 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedElemRestrictionCreateBlockedCurlOriented"); 1101 CeedCall(CeedElemRestrictionCreateBlockedCurlOriented(delegate, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, mem_type, 1102 copy_mode, offsets, curl_orients, rstr)); 1103 CeedCall(CeedDestroy(&delegate)); 1104 return CEED_ERROR_SUCCESS; 1105 } 1106 1107 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 1108 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 1109 CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); 1110 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction must have at least 1 component"); 1111 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction component stride must be at least 1"); 1112 1113 CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_offsets)); 1114 CeedCall(CeedCalloc(num_block * block_size * 3 * elem_size, &block_curl_orients)); 1115 CeedCall(CeedPermutePadOffsets(offsets, block_offsets, num_block, num_elem, block_size, elem_size)); 1116 CeedCall(CeedPermutePadCurlOrients(curl_orients, block_curl_orients, num_block, num_elem, block_size, 3 * elem_size)); 1117 1118 CeedCall(CeedCalloc(1, rstr)); 1119 CeedCall(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, CeedElemRestrictionDestroy_Object, &(*rstr)->obj)); 1120 (*rstr)->num_elem = num_elem; 1121 (*rstr)->elem_size = elem_size; 1122 (*rstr)->num_comp = num_comp; 1123 (*rstr)->comp_stride = comp_stride; 1124 (*rstr)->l_size = l_size; 1125 (*rstr)->e_size = (CeedSize)num_block * (CeedSize)block_size * (CeedSize)elem_size * (CeedSize)num_comp; 1126 (*rstr)->num_block = num_block; 1127 (*rstr)->block_size = block_size; 1128 (*rstr)->rstr_type = CEED_RESTRICTION_CURL_ORIENTED; 1129 CeedCall(ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, (const CeedInt *)block_offsets, NULL, 1130 (const CeedInt8 *)block_curl_orients, *rstr)); 1131 if (copy_mode == CEED_OWN_POINTER) CeedCall(CeedFree(&offsets)); 1132 return CEED_ERROR_SUCCESS; 1133 } 1134 1135 /** 1136 @brief Create a blocked strided `CeedElemRestriction`, typically only used by backends 1137 1138 @param[in] ceed `Ceed` context used to create the `CeedElemRestriction` 1139 @param[in] num_elem Number of elements described by the restriction 1140 @param[in] elem_size Size (number of "nodes") per element 1141 @param[in] block_size Number of elements in a block 1142 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 1143 @param[in] l_size The size of the L-vector. 1144 This vector may be larger than the elements and fields given by this restriction. 1145 @param[in] strides Array for strides between `[nodes, components, elements]`. 1146 Data for node `i`, component `j`, element `k` can be found in the L-vector at index `i*strides[0] + j*strides[1] +k*strides[2]`. 1147 @ref CEED_STRIDES_BACKEND may be used for `CeedVector` ordered by the same `Ceed` backend. 1148 @param[out] rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 1149 1150 @return An error code: 0 - success, otherwise - failure 1151 1152 @ref User 1153 **/ 1154 int CeedElemRestrictionCreateBlockedStrided(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, CeedSize l_size, 1155 const CeedInt strides[3], CeedElemRestriction *rstr) { 1156 CeedInt num_block = (num_elem / block_size) + !!(num_elem % block_size); 1157 1158 if (!ceed->ElemRestrictionCreateBlocked) { 1159 Ceed delegate; 1160 1161 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 1162 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedElemRestrictionCreateBlockedStrided"); 1163 CeedCall(CeedElemRestrictionCreateBlockedStrided(delegate, num_elem, elem_size, block_size, num_comp, l_size, strides, rstr)); 1164 CeedCall(CeedDestroy(&delegate)); 1165 return CEED_ERROR_SUCCESS; 1166 } 1167 1168 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 1169 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 1170 CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); 1171 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction must have at least 1 component"); 1172 CeedCheck(l_size >= (CeedSize)num_elem * (CeedSize)elem_size * (CeedSize)num_comp, ceed, CEED_ERROR_DIMENSION, 1173 "L-vector size must be at least num_elem * elem_size * num_comp. Expected: > %" CeedSize_FMT " Found: %" CeedSize_FMT, 1174 (CeedSize)num_elem * (CeedSize)elem_size * (CeedSize)num_comp, l_size); 1175 1176 CeedCall(CeedCalloc(1, rstr)); 1177 CeedCall(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, CeedElemRestrictionDestroy_Object, &(*rstr)->obj)); 1178 (*rstr)->num_elem = num_elem; 1179 (*rstr)->elem_size = elem_size; 1180 (*rstr)->num_comp = num_comp; 1181 (*rstr)->l_size = l_size; 1182 (*rstr)->e_size = (CeedSize)num_block * (CeedSize)block_size * (CeedSize)elem_size * (CeedSize)num_comp; 1183 (*rstr)->num_block = num_block; 1184 (*rstr)->block_size = block_size; 1185 (*rstr)->rstr_type = CEED_RESTRICTION_STRIDED; 1186 CeedCall(CeedMalloc(3, &(*rstr)->strides)); 1187 for (CeedInt i = 0; i < 3; i++) (*rstr)->strides[i] = strides[i]; 1188 CeedCall(ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, NULL, NULL, NULL, *rstr)); 1189 return CEED_ERROR_SUCCESS; 1190 } 1191 1192 /** 1193 @brief Copy the pointer to a `CeedElemRestriction` and set @ref CeedElemRestrictionApply() implementation to use the unsigned version. 1194 1195 Both pointers should be destroyed with @ref CeedElemRestrictionDestroy(). 1196 1197 @param[in] rstr `CeedElemRestriction` to create unsigned reference to 1198 @param[in,out] rstr_unsigned Variable to store unsigned `CeedElemRestriction` 1199 1200 @return An error code: 0 - success, otherwise - failure 1201 1202 @ref User 1203 **/ 1204 int CeedElemRestrictionCreateUnsignedCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_unsigned) { 1205 CeedCall(CeedCalloc(1, rstr_unsigned)); 1206 1207 // Copy old rstr 1208 memcpy(*rstr_unsigned, rstr, sizeof(struct CeedElemRestriction_private)); 1209 CeedCall(CeedObjectCreate(CeedElemRestrictionReturnCeed(rstr), CeedElemRestrictionView_Object, CeedElemRestrictionDestroy_Object, 1210 &(*rstr_unsigned)->obj)); 1211 (*rstr_unsigned)->strides = NULL; 1212 if (rstr->strides) { 1213 CeedCall(CeedMalloc(3, &(*rstr_unsigned)->strides)); 1214 for (CeedInt i = 0; i < 3; i++) (*rstr_unsigned)->strides[i] = rstr->strides[i]; 1215 } 1216 CeedCall(CeedElemRestrictionReferenceCopy(rstr, &(*rstr_unsigned)->rstr_base)); 1217 1218 // Override Apply 1219 (*rstr_unsigned)->Apply = rstr->ApplyUnsigned; 1220 return CEED_ERROR_SUCCESS; 1221 } 1222 1223 /** 1224 @brief Copy the pointer to a `CeedElemRestriction` and set @ref CeedElemRestrictionApply() implementation to use the unoriented version. 1225 1226 Both pointers should be destroyed with @ref CeedElemRestrictionDestroy(). 1227 1228 @param[in] rstr `CeedElemRestriction` to create unoriented reference to 1229 @param[in,out] rstr_unoriented Variable to store unoriented `CeedElemRestriction` 1230 1231 @return An error code: 0 - success, otherwise - failure 1232 1233 @ref User 1234 **/ 1235 int CeedElemRestrictionCreateUnorientedCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_unoriented) { 1236 CeedCall(CeedCalloc(1, rstr_unoriented)); 1237 1238 // Copy old rstr 1239 memcpy(*rstr_unoriented, rstr, sizeof(struct CeedElemRestriction_private)); 1240 CeedCall(CeedObjectCreate(CeedElemRestrictionReturnCeed(rstr), CeedElemRestrictionView_Object, CeedElemRestrictionDestroy_Object, 1241 &(*rstr_unoriented)->obj)); 1242 (*rstr_unoriented)->strides = NULL; 1243 if (rstr->strides) { 1244 CeedCall(CeedMalloc(3, &(*rstr_unoriented)->strides)); 1245 for (CeedInt i = 0; i < 3; i++) (*rstr_unoriented)->strides[i] = rstr->strides[i]; 1246 } 1247 CeedCall(CeedElemRestrictionReferenceCopy(rstr, &(*rstr_unoriented)->rstr_base)); 1248 1249 // Override Apply 1250 (*rstr_unoriented)->Apply = rstr->ApplyUnoriented; 1251 return CEED_ERROR_SUCCESS; 1252 } 1253 1254 /** 1255 @brief Copy the pointer to a `CeedElemRestriction`. 1256 1257 Both pointers should be destroyed with @ref CeedElemRestrictionDestroy(). 1258 1259 Note: If the value of `*rstr_copy` passed into this function is non-`NULL`, then it is assumed that `*rstr_copy` is a pointer to a `CeedElemRestriction`. 1260 This `CeedElemRestriction` will be destroyed if `*rstr_copy` is the only reference to this `CeedElemRestriction`. 1261 1262 @param[in] rstr `CeedElemRestriction` to copy reference to 1263 @param[in,out] rstr_copy Variable to store copied reference 1264 1265 @return An error code: 0 - success, otherwise - failure 1266 1267 @ref User 1268 **/ 1269 int CeedElemRestrictionReferenceCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_copy) { 1270 if (rstr != CEED_ELEMRESTRICTION_NONE) CeedCall(CeedElemRestrictionReference(rstr)); 1271 CeedCall(CeedElemRestrictionDestroy(rstr_copy)); 1272 *rstr_copy = rstr; 1273 return CEED_ERROR_SUCCESS; 1274 } 1275 1276 /** 1277 @brief Create `CeedVector` associated with a `CeedElemRestriction` 1278 1279 @param[in] rstr `CeedElemRestriction` 1280 @param[out] l_vec The address of the L-vector to be created, or `NULL` 1281 @param[out] e_vec The address of the E-vector to be created, or `NULL` 1282 1283 @return An error code: 0 - success, otherwise - failure 1284 1285 @ref User 1286 **/ 1287 int CeedElemRestrictionCreateVector(CeedElemRestriction rstr, CeedVector *l_vec, CeedVector *e_vec) { 1288 CeedSize e_size, l_size; 1289 Ceed ceed; 1290 1291 CeedCall(CeedElemRestrictionGetCeed(rstr, &ceed)); 1292 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &l_size)); 1293 CeedCall(CeedElemRestrictionGetEVectorSize(rstr, &e_size)); 1294 if (l_vec) CeedCall(CeedVectorCreate(ceed, l_size, l_vec)); 1295 if (e_vec) CeedCall(CeedVectorCreate(ceed, e_size, e_vec)); 1296 CeedCall(CeedDestroy(&ceed)); 1297 return CEED_ERROR_SUCCESS; 1298 } 1299 1300 /** 1301 @brief Restrict an L-vector to an E-vector or apply its transpose 1302 1303 @param[in] rstr `CeedElemRestriction` 1304 @param[in] t_mode Apply restriction or transpose 1305 @param[in] u Input vector (of size `l_size` when `t_mode` = @ref CEED_NOTRANSPOSE) 1306 @param[out] ru Output vector (of shape `[num_elem * elem_size]` when `t_mode` = @ref CEED_NOTRANSPOSE). 1307 Ordering of the e-vector is decided by the backend. 1308 @param[in] request Request or @ref CEED_REQUEST_IMMEDIATE 1309 1310 @return An error code: 0 - success, otherwise - failure 1311 1312 @ref User 1313 **/ 1314 int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, CeedRequest *request) { 1315 CeedSize min_u_len, min_ru_len, len; 1316 CeedInt num_elem; 1317 1318 if (t_mode == CEED_NOTRANSPOSE) { 1319 CeedCall(CeedElemRestrictionGetEVectorSize(rstr, &min_ru_len)); 1320 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &min_u_len)); 1321 } else { 1322 CeedCall(CeedElemRestrictionGetEVectorSize(rstr, &min_u_len)); 1323 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &min_ru_len)); 1324 } 1325 CeedCall(CeedVectorGetLength(u, &len)); 1326 CeedCheck(min_u_len <= len, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_DIMENSION, 1327 "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", len, min_ru_len, 1328 min_u_len); 1329 CeedCall(CeedVectorGetLength(ru, &len)); 1330 CeedCheck(min_ru_len <= len, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_DIMENSION, 1331 "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", len, min_u_len, 1332 min_ru_len); 1333 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 1334 if (num_elem > 0) CeedCall(rstr->Apply(rstr, t_mode, u, ru, request)); 1335 return CEED_ERROR_SUCCESS; 1336 } 1337 1338 /** 1339 @brief Restrict an L-vector of points to a single element or apply its transpose 1340 1341 @param[in] rstr `CeedElemRestriction` 1342 @param[in] elem Element number in range `[0, num_elem)` 1343 @param[in] t_mode Apply restriction or transpose 1344 @param[in] u Input vector (of size `l_size` when `t_mode` = @ref CEED_NOTRANSPOSE) 1345 @param[out] ru Output vector (of shape [`num_points * num_comp]` when `t_mode` = @ref CEED_NOTRANSPOSE). 1346 Ordering of the e-vector is decided by the backend. 1347 @param[in] request Request or @ref CEED_REQUEST_IMMEDIATE 1348 1349 @return An error code: 0 - success, otherwise - failure 1350 1351 @ref User 1352 **/ 1353 int CeedElemRestrictionApplyAtPointsInElement(CeedElemRestriction rstr, CeedInt elem, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, 1354 CeedRequest *request) { 1355 CeedSize min_u_len, min_ru_len, len; 1356 CeedInt num_elem; 1357 1358 CeedCheck(rstr->ApplyAtPointsInElement, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_UNSUPPORTED, 1359 "Backend does not implement CeedElemRestrictionApplyAtPointsInElement"); 1360 1361 if (t_mode == CEED_NOTRANSPOSE) { 1362 CeedInt num_points, num_comp; 1363 1364 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &min_u_len)); 1365 CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr, elem, &num_points)); 1366 CeedCall(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 1367 min_ru_len = (CeedSize)num_points * (CeedSize)num_comp; 1368 } else { 1369 CeedInt num_points, num_comp; 1370 1371 CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr, elem, &num_points)); 1372 CeedCall(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 1373 min_u_len = (CeedSize)num_points * (CeedSize)num_comp; 1374 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &min_ru_len)); 1375 } 1376 CeedCall(CeedVectorGetLength(u, &len)); 1377 CeedCheck(min_u_len <= len, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_DIMENSION, 1378 "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT 1379 ") for element %" CeedInt_FMT, 1380 len, min_ru_len, min_u_len, elem); 1381 CeedCall(CeedVectorGetLength(ru, &len)); 1382 CeedCheck(min_ru_len <= len, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_DIMENSION, 1383 "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT 1384 ") for element %" CeedInt_FMT, 1385 len, min_ru_len, min_u_len, elem); 1386 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 1387 CeedCheck(elem < num_elem, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_DIMENSION, 1388 "Cannot retrieve element %" CeedInt_FMT ", element %" CeedInt_FMT " > total elements %" CeedInt_FMT "", elem, elem, num_elem); 1389 if (num_elem > 0) CeedCall(rstr->ApplyAtPointsInElement(rstr, elem, t_mode, u, ru, request)); 1390 return CEED_ERROR_SUCCESS; 1391 } 1392 1393 /** 1394 @brief Restrict an L-vector to a block of an E-vector or apply its transpose 1395 1396 @param[in] rstr `CeedElemRestriction` 1397 @param[in] block Block number to restrict to/from, i.e. `block = 0` will handle elements `[0 : block_size]` and `block = 3` will handle elements `[3*block_size : 4*block_size]` 1398 @param[in] t_mode Apply restriction or transpose 1399 @param[in] u Input vector (of size `l_size` when `t_mode` = @ref CEED_NOTRANSPOSE) 1400 @param[out] ru Output vector (of shape `[block_size * elem_size]` when `t_mode` = @ref CEED_NOTRANSPOSE). 1401 Ordering of the e-vector is decided by the backend. 1402 @param[in] request Request or @ref CEED_REQUEST_IMMEDIATE 1403 1404 @return An error code: 0 - success, otherwise - failure 1405 1406 @ref Backend 1407 **/ 1408 int CeedElemRestrictionApplyBlock(CeedElemRestriction rstr, CeedInt block, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, 1409 CeedRequest *request) { 1410 CeedSize min_u_len, min_ru_len, len; 1411 CeedInt block_size, num_elem; 1412 1413 CeedCheck(rstr->ApplyBlock, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_UNSUPPORTED, 1414 "Backend does not implement CeedElemRestrictionApplyBlock"); 1415 1416 CeedCall(CeedElemRestrictionGetBlockSize(rstr, &block_size)); 1417 if (t_mode == CEED_NOTRANSPOSE) { 1418 CeedInt elem_size, num_comp; 1419 1420 CeedCall(CeedElemRestrictionGetElementSize(rstr, &elem_size)); 1421 CeedCall(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 1422 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &min_u_len)); 1423 min_ru_len = (CeedSize)block_size * (CeedSize)elem_size * (CeedSize)num_comp; 1424 } else { 1425 CeedInt elem_size, num_comp; 1426 1427 CeedCall(CeedElemRestrictionGetElementSize(rstr, &elem_size)); 1428 CeedCall(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 1429 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &min_ru_len)); 1430 min_u_len = (CeedSize)block_size * (CeedSize)elem_size * (CeedSize)num_comp; 1431 } 1432 CeedCall(CeedVectorGetLength(u, &len)); 1433 CeedCheck(min_u_len == len, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_DIMENSION, 1434 "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", len, min_u_len, 1435 min_ru_len); 1436 CeedCall(CeedVectorGetLength(ru, &len)); 1437 CeedCheck(min_ru_len == len, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_DIMENSION, 1438 "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", len, min_ru_len, 1439 min_u_len); 1440 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 1441 CeedCheck(block_size * block <= num_elem, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_DIMENSION, 1442 "Cannot retrieve block %" CeedInt_FMT ", element %" CeedInt_FMT " > total elements %" CeedInt_FMT "", block, block_size * block, 1443 num_elem); 1444 CeedCall(rstr->ApplyBlock(rstr, block, t_mode, u, ru, request)); 1445 return CEED_ERROR_SUCCESS; 1446 } 1447 1448 /** 1449 @brief Get the `Ceed` associated with a `CeedElemRestriction` 1450 1451 @param[in] rstr `CeedElemRestriction` 1452 @param[out] ceed Variable to store `Ceed` 1453 1454 @return An error code: 0 - success, otherwise - failure 1455 1456 @ref Advanced 1457 **/ 1458 int CeedElemRestrictionGetCeed(CeedElemRestriction rstr, Ceed *ceed) { 1459 CeedCall(CeedObjectGetCeed((CeedObject)rstr, ceed)); 1460 return CEED_ERROR_SUCCESS; 1461 } 1462 1463 /** 1464 @brief Return the `Ceed` associated with a `CeedElemRestriction` 1465 1466 @param[in] rstr `CeedElemRestriction` 1467 1468 @return `Ceed` associated with the `rstr` 1469 1470 @ref Advanced 1471 **/ 1472 Ceed CeedElemRestrictionReturnCeed(CeedElemRestriction rstr) { return CeedObjectReturnCeed((CeedObject)rstr); } 1473 1474 /** 1475 @brief Get the L-vector component stride 1476 1477 @param[in] rstr `CeedElemRestriction` 1478 @param[out] comp_stride Variable to store component stride 1479 1480 @return An error code: 0 - success, otherwise - failure 1481 1482 @ref Advanced 1483 **/ 1484 int CeedElemRestrictionGetCompStride(CeedElemRestriction rstr, CeedInt *comp_stride) { 1485 *comp_stride = rstr->comp_stride; 1486 return CEED_ERROR_SUCCESS; 1487 } 1488 1489 /** 1490 @brief Get the total number of elements in the range of a `CeedElemRestriction` 1491 1492 @param[in] rstr `CeedElemRestriction` 1493 @param[out] num_elem Variable to store number of elements 1494 1495 @return An error code: 0 - success, otherwise - failure 1496 1497 @ref Advanced 1498 **/ 1499 int CeedElemRestrictionGetNumElements(CeedElemRestriction rstr, CeedInt *num_elem) { 1500 *num_elem = rstr->num_elem; 1501 return CEED_ERROR_SUCCESS; 1502 } 1503 1504 /** 1505 @brief Get the size of elements in the `CeedElemRestriction` 1506 1507 @param[in] rstr `CeedElemRestriction` 1508 @param[out] elem_size Variable to store size of elements 1509 1510 @return An error code: 0 - success, otherwise - failure 1511 1512 @ref Advanced 1513 **/ 1514 int CeedElemRestrictionGetElementSize(CeedElemRestriction rstr, CeedInt *elem_size) { 1515 *elem_size = rstr->elem_size; 1516 return CEED_ERROR_SUCCESS; 1517 } 1518 1519 /** 1520 1521 @brief Get the number of points in the offsets array for a points `CeedElemRestriction` 1522 1523 @param[in] rstr `CeedElemRestriction` 1524 @param[out] num_points The number of points in the offsets array 1525 1526 @return An error code: 0 - success, otherwise - failure 1527 1528 @ref User 1529 **/ 1530 int CeedElemRestrictionGetNumPoints(CeedElemRestriction rstr, CeedInt *num_points) { 1531 CeedRestrictionType rstr_type; 1532 1533 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 1534 CeedCheck(rstr_type == CEED_RESTRICTION_POINTS, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_INCOMPATIBLE, 1535 "Can only retrieve the number of points for a points CeedElemRestriction"); 1536 1537 *num_points = rstr->num_points; 1538 return CEED_ERROR_SUCCESS; 1539 } 1540 1541 /** 1542 1543 @brief Get the number of points in an element of a `CeedElemRestriction` at points 1544 1545 @param[in] rstr `CeedElemRestriction` 1546 @param[in] elem Index number of element to retrieve the number of points for 1547 @param[out] num_points The number of points in the element at index elem 1548 1549 @return An error code: 0 - success, otherwise - failure 1550 1551 @ref User 1552 **/ 1553 int CeedElemRestrictionGetNumPointsInElement(CeedElemRestriction rstr, CeedInt elem, CeedInt *num_points) { 1554 CeedRestrictionType rstr_type; 1555 const CeedInt *offsets; 1556 1557 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 1558 CeedCheck(rstr_type == CEED_RESTRICTION_POINTS, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_INCOMPATIBLE, 1559 "Can only retrieve the number of points for a points CeedElemRestriction"); 1560 1561 CeedCall(CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets)); 1562 *num_points = offsets[elem + 1] - offsets[elem]; 1563 CeedCall(CeedElemRestrictionRestoreOffsets(rstr, &offsets)); 1564 return CEED_ERROR_SUCCESS; 1565 } 1566 1567 /** 1568 @brief Get the minimum and/or maximum number of points in an element for a `CeedElemRestriction` at points 1569 1570 @param[in] rstr `CeedElemRestriction` 1571 @param[out] min_points Variable to minimum number of points in an element, or `NULL` 1572 @param[out] max_points Variable to maximum number of points in an element, or `NULL` 1573 1574 @return An error code: 0 - success, otherwise - failure 1575 1576 @ref Advanced 1577 **/ 1578 int CeedElemRestrictionGetMinMaxPointsInElement(CeedElemRestriction rstr, CeedInt *min_points, CeedInt *max_points) { 1579 CeedInt num_elem, num_points; 1580 CeedRestrictionType rstr_type; 1581 1582 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 1583 CeedCheck(rstr_type == CEED_RESTRICTION_POINTS, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_INCOMPATIBLE, 1584 "Cannot compute min/max points for a CeedElemRestriction that does not use points"); 1585 1586 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 1587 1588 // Exit early if there are no elements 1589 if (num_elem == 0) { 1590 if (min_points) *min_points = 0; 1591 if (max_points) *max_points = 0; 1592 return CEED_ERROR_SUCCESS; 1593 } 1594 1595 // Initialize to the number of points in the first element 1596 CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr, 0, &num_points)); 1597 if (min_points) *min_points = num_points; 1598 if (max_points) *max_points = num_points; 1599 for (CeedInt e = 1; e < num_elem; e++) { 1600 CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr, e, &num_points)); 1601 if (min_points) *min_points = CeedIntMin(num_points, *min_points); 1602 if (max_points) *max_points = CeedIntMax(num_points, *max_points); 1603 } 1604 return CEED_ERROR_SUCCESS; 1605 } 1606 1607 /** 1608 @brief Get the maximum number of points in an element for a `CeedElemRestriction` at points 1609 1610 @param[in] rstr `CeedElemRestriction` 1611 @param[out] max_points Variable to store maximum number of points in an element 1612 1613 @return An error code: 0 - success, otherwise - failure 1614 1615 @ref User 1616 1617 @see CeedElemRestrictionGetMinMaxPointsInElement() 1618 **/ 1619 int CeedElemRestrictionGetMaxPointsInElement(CeedElemRestriction rstr, CeedInt *max_points) { 1620 return CeedElemRestrictionGetMinMaxPointsInElement(rstr, NULL, max_points); 1621 } 1622 1623 /** 1624 @brief Get the minimum number of points in an element for a `CeedElemRestriction` at points 1625 1626 @param[in] rstr `CeedElemRestriction` 1627 @param[out] min_points Variable to store minimum number of points in an element 1628 1629 @return An error code: 0 - success, otherwise - failure 1630 1631 @ref User 1632 1633 @see CeedElemRestrictionGetMinMaxPointsInElement() 1634 **/ 1635 int CeedElemRestrictionGetMinPointsInElement(CeedElemRestriction rstr, CeedInt *min_points) { 1636 return CeedElemRestrictionGetMinMaxPointsInElement(rstr, min_points, NULL); 1637 } 1638 1639 /** 1640 @brief Get the size of the l-vector for a `CeedElemRestriction` 1641 1642 @param[in] rstr `CeedElemRestriction` 1643 @param[out] l_size Variable to store l-vector size 1644 1645 @return An error code: 0 - success, otherwise - failure 1646 1647 @ref Advanced 1648 **/ 1649 int CeedElemRestrictionGetLVectorSize(CeedElemRestriction rstr, CeedSize *l_size) { 1650 *l_size = rstr->l_size; 1651 return CEED_ERROR_SUCCESS; 1652 } 1653 1654 /** 1655 @brief Get the size of the e-vector for a `CeedElemRestriction` 1656 1657 @param[in] rstr `CeedElemRestriction` 1658 @param[out] e_size Variable to store e-vector size 1659 1660 @return An error code: 0 - success, otherwise - failure 1661 1662 @ref Advanced 1663 **/ 1664 int CeedElemRestrictionGetEVectorSize(CeedElemRestriction rstr, CeedSize *e_size) { 1665 *e_size = rstr->e_size; 1666 return CEED_ERROR_SUCCESS; 1667 } 1668 1669 /** 1670 @brief Get the number of components in the elements of a `CeedElemRestriction` 1671 1672 @param[in] rstr `CeedElemRestriction` 1673 @param[out] num_comp Variable to store number of components 1674 1675 @return An error code: 0 - success, otherwise - failure 1676 1677 @ref Advanced 1678 **/ 1679 int CeedElemRestrictionGetNumComponents(CeedElemRestriction rstr, CeedInt *num_comp) { 1680 *num_comp = rstr->num_comp; 1681 return CEED_ERROR_SUCCESS; 1682 } 1683 1684 /** 1685 @brief Get the number of blocks in a `CeedElemRestriction` 1686 1687 @param[in] rstr `CeedElemRestriction` 1688 @param[out] num_block Variable to store number of blocks 1689 1690 @return An error code: 0 - success, otherwise - failure 1691 1692 @ref Advanced 1693 **/ 1694 int CeedElemRestrictionGetNumBlocks(CeedElemRestriction rstr, CeedInt *num_block) { 1695 *num_block = rstr->num_block; 1696 return CEED_ERROR_SUCCESS; 1697 } 1698 1699 /** 1700 @brief Get the size of blocks in the `CeedElemRestriction` 1701 1702 @param[in] rstr `CeedElemRestriction` 1703 @param[out] block_size Variable to store size of blocks 1704 1705 @return An error code: 0 - success, otherwise - failure 1706 1707 @ref Advanced 1708 **/ 1709 int CeedElemRestrictionGetBlockSize(CeedElemRestriction rstr, CeedInt *block_size) { 1710 *block_size = rstr->block_size; 1711 return CEED_ERROR_SUCCESS; 1712 } 1713 1714 /** 1715 @brief Get the multiplicity of nodes in a `CeedElemRestriction` 1716 1717 @param[in] rstr `CeedElemRestriction` 1718 @param[out] mult Vector to store multiplicity (of size `l_size`) 1719 1720 @return An error code: 0 - success, otherwise - failure 1721 1722 @ref User 1723 **/ 1724 int CeedElemRestrictionGetMultiplicity(CeedElemRestriction rstr, CeedVector mult) { 1725 CeedVector e_vec; 1726 1727 // Create e_vec to hold intermediate computation in E^T (E 1) 1728 CeedCall(CeedElemRestrictionCreateVector(rstr, NULL, &e_vec)); 1729 1730 // Compute e_vec = E * 1 1731 CeedCall(CeedVectorSetValue(mult, 1.0)); 1732 CeedCall(CeedElemRestrictionApply(rstr, CEED_NOTRANSPOSE, mult, e_vec, CEED_REQUEST_IMMEDIATE)); 1733 // Compute multiplicity, mult = E^T * e_vec = E^T (E 1) 1734 CeedCall(CeedVectorSetValue(mult, 0.0)); 1735 CeedCall(CeedElemRestrictionApply(rstr, CEED_TRANSPOSE, e_vec, mult, CEED_REQUEST_IMMEDIATE)); 1736 // Cleanup 1737 CeedCall(CeedVectorDestroy(&e_vec)); 1738 return CEED_ERROR_SUCCESS; 1739 } 1740 1741 /** 1742 @brief Set the number of tabs to indent for @ref CeedElemRestrictionView() output 1743 1744 @param[in] rstr `CeedElemRestriction` to set the number of view tabs 1745 @param[in] num_tabs Number of view tabs to set 1746 1747 @return Error code: 0 - success, otherwise - failure 1748 1749 @ref User 1750 **/ 1751 int CeedElemRestrictionSetNumViewTabs(CeedElemRestriction rstr, CeedInt num_tabs) { 1752 CeedCall(CeedObjectSetNumViewTabs((CeedObject)rstr, num_tabs)); 1753 return CEED_ERROR_SUCCESS; 1754 } 1755 1756 /** 1757 @brief Get the number of tabs to indent for @ref CeedElemRestrictionView() output 1758 1759 @param[in] rstr `CeedElemRestriction` to get the number of view tabs 1760 @param[out] num_tabs Number of view tabs 1761 1762 @return Error code: 0 - success, otherwise - failure 1763 1764 @ref User 1765 **/ 1766 int CeedElemRestrictionGetNumViewTabs(CeedElemRestriction rstr, CeedInt *num_tabs) { 1767 CeedCall(CeedObjectGetNumViewTabs((CeedObject)rstr, num_tabs)); 1768 return CEED_ERROR_SUCCESS; 1769 } 1770 1771 /** 1772 @brief View a `CeedElemRestriction` 1773 1774 @param[in] rstr `CeedElemRestriction` to view 1775 @param[in] stream Stream to write; typically `stdout` or a file 1776 1777 @return Error code: 0 - success, otherwise - failure 1778 1779 @ref User 1780 **/ 1781 int CeedElemRestrictionView(CeedElemRestriction rstr, FILE *stream) { 1782 char *tabs = NULL; 1783 CeedRestrictionType rstr_type; 1784 1785 { 1786 CeedInt num_tabs = 0; 1787 1788 CeedCall(CeedElemRestrictionGetNumViewTabs(rstr, &num_tabs)); 1789 CeedCall(CeedCalloc(CEED_TAB_WIDTH * num_tabs + 1, &tabs)); 1790 for (CeedInt i = 0; i < CEED_TAB_WIDTH * num_tabs; i++) tabs[i] = ' '; 1791 } 1792 1793 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 1794 if (rstr_type == CEED_RESTRICTION_POINTS) { 1795 CeedInt max_points; 1796 1797 CeedCall(CeedElemRestrictionGetMaxPointsInElement(rstr, &max_points)); 1798 fprintf(stream, 1799 "%sCeedElemRestriction at points from (%" CeedSize_FMT ", %" CeedInt_FMT ") to %" CeedInt_FMT " elements with a maximum of %" CeedInt_FMT 1800 " points on an element\n", 1801 tabs, rstr->l_size, rstr->num_comp, rstr->num_elem, max_points); 1802 } else { 1803 char strides_str[500]; 1804 1805 if (rstr->strides) { 1806 sprintf(strides_str, "[%" CeedInt_FMT ", %" CeedInt_FMT ", %" CeedInt_FMT "]", rstr->strides[0], rstr->strides[1], rstr->strides[2]); 1807 } else { 1808 sprintf(strides_str, "%" CeedInt_FMT, rstr->comp_stride); 1809 } 1810 fprintf(stream, 1811 "%s%sCeedElemRestriction from (%" CeedSize_FMT ", %" CeedInt_FMT ") to %" CeedInt_FMT " elements with %" CeedInt_FMT 1812 " nodes each and %s %s\n", 1813 tabs, rstr->block_size > 1 ? "Blocked " : "", rstr->l_size, rstr->num_comp, rstr->num_elem, rstr->elem_size, 1814 rstr->strides ? "strides" : "component stride", strides_str); 1815 } 1816 CeedCall(CeedFree(&tabs)); 1817 return CEED_ERROR_SUCCESS; 1818 } 1819 1820 /** 1821 @brief Destroy a `CeedElemRestriction` 1822 1823 @param[in,out] rstr `CeedElemRestriction` to destroy 1824 1825 @return An error code: 0 - success, otherwise - failure 1826 1827 @ref User 1828 **/ 1829 int CeedElemRestrictionDestroy(CeedElemRestriction *rstr) { 1830 if (!*rstr || *rstr == CEED_ELEMRESTRICTION_NONE || CeedObjectDereference((CeedObject)*rstr) > 0) { 1831 *rstr = NULL; 1832 return CEED_ERROR_SUCCESS; 1833 } 1834 CeedCheck((*rstr)->num_readers == 0, CeedElemRestrictionReturnCeed(*rstr), CEED_ERROR_ACCESS, 1835 "Cannot destroy CeedElemRestriction, a process has read access to the offset data"); 1836 1837 // Only destroy backend data once between rstr and unsigned copy 1838 if ((*rstr)->rstr_base) CeedCall(CeedElemRestrictionDestroy(&(*rstr)->rstr_base)); 1839 else if ((*rstr)->Destroy) CeedCall((*rstr)->Destroy(*rstr)); 1840 1841 CeedCall(CeedFree(&(*rstr)->strides)); 1842 CeedCall(CeedObjectDestroy_Private(&(*rstr)->obj)); 1843 CeedCall(CeedFree(rstr)); 1844 return CEED_ERROR_SUCCESS; 1845 } 1846 1847 /// @} 1848