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