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 118 /// ---------------------------------------------------------------------------- 119 /// CeedElemRestriction Backend API 120 /// ---------------------------------------------------------------------------- 121 /// @addtogroup CeedElemRestrictionBackend 122 /// @{ 123 124 /** 125 @brief Get the type of a `CeedElemRestriction` 126 127 @param[in] rstr `CeedElemRestriction` 128 @param[out] rstr_type Variable to store restriction type 129 130 @return An error code: 0 - success, otherwise - failure 131 132 @ref Backend 133 **/ 134 int CeedElemRestrictionGetType(CeedElemRestriction rstr, CeedRestrictionType *rstr_type) { 135 *rstr_type = rstr->rstr_type; 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 CeedCall(CeedObjectReference((CeedObject)rstr)); 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(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, &(*rstr)->obj)); 662 (*rstr)->num_elem = num_elem; 663 (*rstr)->elem_size = elem_size; 664 (*rstr)->num_comp = num_comp; 665 (*rstr)->comp_stride = comp_stride; 666 (*rstr)->l_size = l_size; 667 (*rstr)->e_size = (CeedSize)num_elem * (CeedSize)elem_size * (CeedSize)num_comp; 668 (*rstr)->num_block = num_elem; 669 (*rstr)->block_size = 1; 670 (*rstr)->rstr_type = CEED_RESTRICTION_STANDARD; 671 CeedCall(ceed->ElemRestrictionCreate(mem_type, copy_mode, offsets, NULL, NULL, *rstr)); 672 return CEED_ERROR_SUCCESS; 673 } 674 675 /** 676 @brief Create a `CeedElemRestriction` with orientation signs 677 678 @param[in] ceed `Ceed` context used to create the `CeedElemRestriction` 679 @param[in] num_elem Number of elements described in the `offsets` array 680 @param[in] elem_size Size (number of "nodes") per element 681 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 682 @param[in] comp_stride Stride between components for the same L-vector "node". 683 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`. 684 @param[in] l_size The size of the L-vector. 685 This vector may be larger than the elements and fields given by this restriction. 686 @param[in] mem_type Memory type of the `offsets` array, see @ref CeedMemType 687 @param[in] copy_mode Copy mode for the `offsets` array, see @ref CeedCopyMode 688 @param[in] offsets Array of shape `[num_elem, elem_size]`. 689 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`. 690 All offsets must be in the range `[0, l_size - 1]`. 691 @param[in] orients Boolean array of shape `[num_elem, elem_size]` with `false` for positively oriented and `true` to flip the orientation 692 @param[out] rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 693 694 @return An error code: 0 - success, otherwise - failure 695 696 @ref User 697 **/ 698 int CeedElemRestrictionCreateOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt num_comp, CeedInt comp_stride, CeedSize l_size, 699 CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, const bool *orients, 700 CeedElemRestriction *rstr) { 701 if (!ceed->ElemRestrictionCreate) { 702 Ceed delegate; 703 704 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 705 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedElemRestrictionCreateOriented"); 706 CeedCall(CeedElemRestrictionCreateOriented(delegate, num_elem, elem_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, orients, 707 rstr)); 708 CeedCall(CeedDestroy(&delegate)); 709 return CEED_ERROR_SUCCESS; 710 } 711 712 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 713 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 714 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction must have at least 1 component"); 715 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction component stride must be at least 1"); 716 717 CeedCall(CeedCalloc(1, rstr)); 718 CeedCall(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, &(*rstr)->obj)); 719 (*rstr)->num_elem = num_elem; 720 (*rstr)->elem_size = elem_size; 721 (*rstr)->num_comp = num_comp; 722 (*rstr)->comp_stride = comp_stride; 723 (*rstr)->l_size = l_size; 724 (*rstr)->e_size = (CeedSize)num_elem * (CeedSize)elem_size * (CeedSize)num_comp; 725 (*rstr)->num_block = num_elem; 726 (*rstr)->block_size = 1; 727 (*rstr)->rstr_type = CEED_RESTRICTION_ORIENTED; 728 CeedCall(ceed->ElemRestrictionCreate(mem_type, copy_mode, offsets, orients, NULL, *rstr)); 729 return CEED_ERROR_SUCCESS; 730 } 731 732 /** 733 @brief Create a `CeedElemRestriction` with a general tridiagonal transformation matrix for curl-conforming elements 734 735 @param[in] ceed `Ceed` context used to create the `CeedElemRestriction` 736 @param[in] num_elem Number of elements described in the `offsets` array 737 @param[in] elem_size Size (number of "nodes") per element 738 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 739 @param[in] comp_stride Stride between components for the same L-vector "node". 740 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`. 741 @param[in] l_size The size of the L-vector. 742 This vector may be larger than the elements and fields given by this restriction. 743 @param[in] mem_type Memory type of the `offsets` array, see @ref CeedMemType 744 @param[in] copy_mode Copy mode for the `offsets` array, see @ref CeedCopyMode 745 @param[in] offsets Array of shape `[num_elem, elem_size]`. 746 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`. 747 All offsets must be in the range `[0, l_size - 1]`. 748 @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. 749 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). 750 @param[out] rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 751 752 @return An error code: 0 - success, otherwise - failure 753 754 @ref User 755 **/ 756 int CeedElemRestrictionCreateCurlOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt num_comp, CeedInt comp_stride, CeedSize l_size, 757 CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, const CeedInt8 *curl_orients, 758 CeedElemRestriction *rstr) { 759 if (!ceed->ElemRestrictionCreate) { 760 Ceed delegate; 761 762 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 763 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedElemRestrictionCreateCurlOriented"); 764 CeedCall(CeedElemRestrictionCreateCurlOriented(delegate, num_elem, elem_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, 765 curl_orients, rstr)); 766 CeedCall(CeedDestroy(&delegate)); 767 return CEED_ERROR_SUCCESS; 768 } 769 770 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 771 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 772 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction must have at least 1 component"); 773 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction component stride must be at least 1"); 774 775 CeedCall(CeedCalloc(1, rstr)); 776 CeedCall(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, &(*rstr)->obj)); 777 (*rstr)->num_elem = num_elem; 778 (*rstr)->elem_size = elem_size; 779 (*rstr)->num_comp = num_comp; 780 (*rstr)->comp_stride = comp_stride; 781 (*rstr)->l_size = l_size; 782 (*rstr)->e_size = (CeedSize)num_elem * (CeedSize)elem_size * (CeedSize)num_comp; 783 (*rstr)->num_block = num_elem; 784 (*rstr)->block_size = 1; 785 (*rstr)->rstr_type = CEED_RESTRICTION_CURL_ORIENTED; 786 CeedCall(ceed->ElemRestrictionCreate(mem_type, copy_mode, offsets, NULL, curl_orients, *rstr)); 787 return CEED_ERROR_SUCCESS; 788 } 789 790 /** 791 @brief Create a strided `CeedElemRestriction` 792 793 @param[in] ceed `Ceed` context used to create the `CeedElemRestriction` 794 @param[in] num_elem Number of elements described by the restriction 795 @param[in] elem_size Size (number of "nodes") per element 796 @param[in] num_comp Number of field components per interpolation "node" (1 for scalar fields) 797 @param[in] l_size The size of the L-vector. 798 This vector may be larger than the elements and fields given by this restriction. 799 @param[in] strides Array for strides between `[nodes, components, elements]`. 800 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]`. 801 @ref CEED_STRIDES_BACKEND may be used for `CeedVector` ordered by the same `Ceed` backend. 802 `CEED_STRIDES_BACKEND` should only be used pass data between `CeedOperator` created with the same `Ceed` backend. 803 The L-vector layout will, in general, be different between `Ceed` backends. 804 @param[out] rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 805 806 @return An error code: 0 - success, otherwise - failure 807 808 @ref User 809 **/ 810 int CeedElemRestrictionCreateStrided(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt num_comp, CeedSize l_size, const CeedInt strides[3], 811 CeedElemRestriction *rstr) { 812 if (!ceed->ElemRestrictionCreate) { 813 Ceed delegate; 814 815 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 816 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedElemRestrictionCreateStrided"); 817 CeedCall(CeedElemRestrictionCreateStrided(delegate, num_elem, elem_size, num_comp, l_size, strides, rstr)); 818 CeedCall(CeedDestroy(&delegate)); 819 return CEED_ERROR_SUCCESS; 820 } 821 822 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 823 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 824 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction must have at least 1 component"); 825 CeedCheck(l_size >= (CeedSize)num_elem * (CeedSize)elem_size * (CeedSize)num_comp, ceed, CEED_ERROR_DIMENSION, 826 "L-vector size must be at least num_elem * elem_size * num_comp. Expected: > %" CeedSize_FMT " Found: %" CeedSize_FMT, 827 (CeedSize)num_elem * (CeedSize)elem_size * (CeedSize)num_comp, l_size); 828 829 CeedCall(CeedCalloc(1, rstr)); 830 CeedCall(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, &(*rstr)->obj)); 831 (*rstr)->num_elem = num_elem; 832 (*rstr)->elem_size = elem_size; 833 (*rstr)->num_comp = num_comp; 834 (*rstr)->l_size = l_size; 835 (*rstr)->e_size = (CeedSize)num_elem * (CeedSize)elem_size * (CeedSize)num_comp; 836 (*rstr)->num_block = num_elem; 837 (*rstr)->block_size = 1; 838 (*rstr)->rstr_type = CEED_RESTRICTION_STRIDED; 839 CeedCall(CeedMalloc(3, &(*rstr)->strides)); 840 for (CeedInt i = 0; i < 3; i++) (*rstr)->strides[i] = strides[i]; 841 CeedCall(ceed->ElemRestrictionCreate(CEED_MEM_HOST, CEED_OWN_POINTER, NULL, NULL, NULL, *rstr)); 842 return CEED_ERROR_SUCCESS; 843 } 844 845 /** 846 @brief Create a points `CeedElemRestriction`, for restricting for restricting from a all local points to the current element in which they are located. 847 848 The offsets array is arranged as 849 850 element_0_start_index 851 element_1_start_index 852 ... 853 element_n_start_index 854 element_n_stop_index 855 element_0_point_0 856 element_0_point_1 857 ... 858 859 @param[in] ceed `Ceed` context used to create the `CeedElemRestriction` 860 @param[in] num_elem Number of elements described in the `offsets` array 861 @param[in] num_points Number of points described in the `offsets` array 862 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields). 863 Components are assumed to be contiguous by point. 864 @param[in] l_size The size of the L-vector. 865 This vector may be larger than the elements and fields given by this restriction. 866 @param[in] mem_type Memory type of the `offsets` array, see @ref CeedMemType 867 @param[in] copy_mode Copy mode for the `offsets` array, see @ref CeedCopyMode 868 @param[in] offsets Array of size `num_elem + 1 + num_points`. 869 The first portion of the offsets array holds the ranges of indices corresponding to each element. 870 The second portion holds the indices for each element. 871 @param[out] rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 872 873 @return An error code: 0 - success, otherwise - failure 874 875 @ref Backend 876 **/ 877 int CeedElemRestrictionCreateAtPoints(Ceed ceed, CeedInt num_elem, CeedInt num_points, CeedInt num_comp, CeedSize l_size, CeedMemType mem_type, 878 CeedCopyMode copy_mode, const CeedInt *offsets, CeedElemRestriction *rstr) { 879 if (!ceed->ElemRestrictionCreateAtPoints) { 880 Ceed delegate; 881 882 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 883 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedElemRestrictionCreateAtPoints"); 884 CeedCall(CeedElemRestrictionCreateAtPoints(delegate, num_elem, num_points, num_comp, l_size, mem_type, copy_mode, offsets, rstr)); 885 CeedCall(CeedDestroy(&delegate)); 886 return CEED_ERROR_SUCCESS; 887 } 888 889 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 890 CeedCheck(num_points >= 0, ceed, CEED_ERROR_DIMENSION, "Number of points must be non-negative"); 891 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction must have at least 1 component"); 892 CeedCheck(l_size >= (CeedSize)num_points * num_comp, ceed, CEED_ERROR_DIMENSION, 893 "L-vector must be at least num_points * num_comp. Expected: > %" CeedSize_FMT " Found: %" CeedSize_FMT, (CeedSize)num_points * num_comp, 894 l_size); 895 896 CeedCall(CeedCalloc(1, rstr)); 897 CeedCall(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, &(*rstr)->obj)); 898 (*rstr)->num_elem = num_elem; 899 (*rstr)->num_points = num_points; 900 (*rstr)->num_comp = num_comp; 901 (*rstr)->comp_stride = 1; 902 (*rstr)->l_size = l_size; 903 (*rstr)->e_size = (CeedSize)num_points * (CeedSize)num_comp; 904 (*rstr)->num_block = num_elem; 905 (*rstr)->block_size = 1; 906 (*rstr)->rstr_type = CEED_RESTRICTION_POINTS; 907 CeedCall(ceed->ElemRestrictionCreateAtPoints(mem_type, copy_mode, offsets, NULL, NULL, *rstr)); 908 return CEED_ERROR_SUCCESS; 909 } 910 911 /** 912 @brief Create a blocked `CeedElemRestriction`, typically only used by backends 913 914 @param[in] ceed `Ceed` context used to create the `CeedElemRestriction` 915 @param[in] num_elem Number of elements described in the `offsets` array 916 @param[in] elem_size Size (number of unknowns) per element 917 @param[in] block_size Number of elements in a block 918 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 919 @param[in] comp_stride Stride between components for the same L-vector "node". 920 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`. 921 @param[in] l_size The size of the L-vector. 922 This vector may be larger than the elements and fields given by this restriction. 923 @param[in] mem_type Memory type of the `offsets` array, see @ref CeedMemType 924 @param[in] copy_mode Copy mode for the `offsets` array, see @ref CeedCopyMode 925 @param[in] offsets Array of shape `[num_elem, elem_size]`. 926 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`. 927 All offsets must be in the range `[0, l_size - 1]`. 928 The backend will permute and pad this array to the desired ordering for the blocksize, which is typically given by the backend. 929 The default reordering is to interlace elements. 930 @param[out] rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 931 932 @return An error code: 0 - success, otherwise - failure 933 934 @ref Backend 935 **/ 936 int CeedElemRestrictionCreateBlocked(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, CeedInt comp_stride, 937 CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, 938 CeedElemRestriction *rstr) { 939 CeedInt *block_offsets, num_block = (num_elem / block_size) + !!(num_elem % block_size); 940 941 if (!ceed->ElemRestrictionCreateBlocked) { 942 Ceed delegate; 943 944 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 945 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedElemRestrictionCreateBlocked"); 946 CeedCall(CeedElemRestrictionCreateBlocked(delegate, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, 947 rstr)); 948 CeedCall(CeedDestroy(&delegate)); 949 return CEED_ERROR_SUCCESS; 950 } 951 952 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 953 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 954 CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); 955 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction must have at least 1 component"); 956 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction component stride must be at least 1"); 957 958 CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_offsets)); 959 CeedCall(CeedPermutePadOffsets(offsets, block_offsets, num_block, num_elem, block_size, elem_size)); 960 961 CeedCall(CeedCalloc(1, rstr)); 962 CeedCall(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, &(*rstr)->obj)); 963 (*rstr)->num_elem = num_elem; 964 (*rstr)->elem_size = elem_size; 965 (*rstr)->num_comp = num_comp; 966 (*rstr)->comp_stride = comp_stride; 967 (*rstr)->l_size = l_size; 968 (*rstr)->e_size = (CeedSize)num_block * (CeedSize)block_size * (CeedSize)elem_size * (CeedSize)num_comp; 969 (*rstr)->num_block = num_block; 970 (*rstr)->block_size = block_size; 971 (*rstr)->rstr_type = CEED_RESTRICTION_STANDARD; 972 CeedCall(ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, (const CeedInt *)block_offsets, NULL, NULL, *rstr)); 973 if (copy_mode == CEED_OWN_POINTER) CeedCall(CeedFree(&offsets)); 974 return CEED_ERROR_SUCCESS; 975 } 976 977 /** 978 @brief Create a blocked oriented `CeedElemRestriction`, typically only used by backends 979 980 @param[in] ceed `Ceed` context used to create the `CeedElemRestriction` 981 @param[in] num_elem Number of elements described in the `offsets` array. 982 @param[in] elem_size Size (number of unknowns) per element 983 @param[in] block_size Number of elements in a block 984 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 985 @param[in] comp_stride Stride between components for the same L-vector "node". 986 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`. 987 @param[in] l_size The size of the L-vector. 988 This vector may be larger than the elements and fields given by this restriction. 989 @param[in] mem_type Memory type of the `offsets` array, see @ref CeedMemType 990 @param[in] copy_mode Copy mode for the `offsets` array, see @ref CeedCopyMode 991 @param[in] offsets Array of shape `[num_elem, elem_size]`. 992 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`. 993 All offsets must be in the range `[0, l_size - 1]`. 994 The backend will permute and pad this array to the desired ordering for the blocksize, which is typically given by the backend. 995 The default reordering is to interlace elements. 996 @param[in] orients Boolean array of shape `[num_elem, elem_size]` with `false` for positively oriented and `true` to flip the orientation. 997 Will also be permuted and padded similarly to `offsets`. 998 @param[out] rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 999 1000 @return An error code: 0 - success, otherwise - failure 1001 1002 @ref Backend 1003 **/ 1004 int CeedElemRestrictionCreateBlockedOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, 1005 CeedInt comp_stride, CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, 1006 const CeedInt *offsets, const bool *orients, CeedElemRestriction *rstr) { 1007 bool *block_orients; 1008 CeedInt *block_offsets, num_block = (num_elem / block_size) + !!(num_elem % block_size); 1009 1010 if (!ceed->ElemRestrictionCreateBlocked) { 1011 Ceed delegate; 1012 1013 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 1014 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedElemRestrictionCreateBlockedOriented"); 1015 CeedCall(CeedElemRestrictionCreateBlockedOriented(delegate, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, mem_type, copy_mode, 1016 offsets, orients, rstr)); 1017 CeedCall(CeedDestroy(&delegate)); 1018 return CEED_ERROR_SUCCESS; 1019 } 1020 1021 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 1022 CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); 1023 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction must have at least 1 component"); 1024 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction component stride must be at least 1"); 1025 1026 CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_offsets)); 1027 CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_orients)); 1028 CeedCall(CeedPermutePadOffsets(offsets, block_offsets, num_block, num_elem, block_size, elem_size)); 1029 CeedCall(CeedPermutePadOrients(orients, block_orients, num_block, num_elem, block_size, elem_size)); 1030 1031 CeedCall(CeedCalloc(1, rstr)); 1032 CeedCall(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, &(*rstr)->obj)); 1033 (*rstr)->num_elem = num_elem; 1034 (*rstr)->elem_size = elem_size; 1035 (*rstr)->num_comp = num_comp; 1036 (*rstr)->comp_stride = comp_stride; 1037 (*rstr)->l_size = l_size; 1038 (*rstr)->e_size = (CeedSize)num_block * (CeedSize)block_size * (CeedSize)elem_size * (CeedSize)num_comp; 1039 (*rstr)->num_block = num_block; 1040 (*rstr)->block_size = block_size; 1041 (*rstr)->rstr_type = CEED_RESTRICTION_ORIENTED; 1042 CeedCall(ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, (const CeedInt *)block_offsets, (const bool *)block_orients, NULL, 1043 *rstr)); 1044 if (copy_mode == CEED_OWN_POINTER) CeedCall(CeedFree(&offsets)); 1045 return CEED_ERROR_SUCCESS; 1046 } 1047 1048 /** 1049 @brief Create a blocked curl-oriented `CeedElemRestriction`, typically only used by backends 1050 1051 @param[in] ceed `Ceed` context used to create the `CeedElemRestriction` 1052 @param[in] num_elem Number of elements described in the `offsets` array. 1053 @param[in] elem_size Size (number of unknowns) per element 1054 @param[in] block_size Number of elements in a block 1055 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 1056 @param[in] comp_stride Stride between components for the same L-vector "node". 1057 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`. 1058 @param[in] l_size The size of the L-vector. 1059 This vector may be larger than the elements and fields given by this restriction. 1060 @param[in] mem_type Memory type of the `offsets` array, see @ref CeedMemType 1061 @param[in] copy_mode Copy mode for the `offsets` array, see @ref CeedCopyMode 1062 @param[in] offsets Array of shape `[num_elem, elem_size]`. 1063 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`. 1064 All offsets must be in the range `[0, l_size - 1]`. 1065 The backend will permute and pad this array to the desired ordering for the blocksize, which is typically given by the backend. 1066 The default reordering is to interlace elements. 1067 @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. 1068 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). 1069 Will also be permuted and padded similarly to offsets. 1070 @param[out] rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 1071 1072 @return An error code: 0 - success, otherwise - failure 1073 1074 @ref Backend 1075 **/ 1076 int CeedElemRestrictionCreateBlockedCurlOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, 1077 CeedInt comp_stride, CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, 1078 const CeedInt *offsets, const CeedInt8 *curl_orients, CeedElemRestriction *rstr) { 1079 CeedInt8 *block_curl_orients; 1080 CeedInt *block_offsets, num_block = (num_elem / block_size) + !!(num_elem % block_size); 1081 1082 if (!ceed->ElemRestrictionCreateBlocked) { 1083 Ceed delegate; 1084 1085 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 1086 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedElemRestrictionCreateBlockedCurlOriented"); 1087 CeedCall(CeedElemRestrictionCreateBlockedCurlOriented(delegate, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, mem_type, 1088 copy_mode, offsets, curl_orients, rstr)); 1089 CeedCall(CeedDestroy(&delegate)); 1090 return CEED_ERROR_SUCCESS; 1091 } 1092 1093 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 1094 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 1095 CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); 1096 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction must have at least 1 component"); 1097 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction component stride must be at least 1"); 1098 1099 CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_offsets)); 1100 CeedCall(CeedCalloc(num_block * block_size * 3 * elem_size, &block_curl_orients)); 1101 CeedCall(CeedPermutePadOffsets(offsets, block_offsets, num_block, num_elem, block_size, elem_size)); 1102 CeedCall(CeedPermutePadCurlOrients(curl_orients, block_curl_orients, num_block, num_elem, block_size, 3 * elem_size)); 1103 1104 CeedCall(CeedCalloc(1, rstr)); 1105 CeedCall(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, &(*rstr)->obj)); 1106 (*rstr)->num_elem = num_elem; 1107 (*rstr)->elem_size = elem_size; 1108 (*rstr)->num_comp = num_comp; 1109 (*rstr)->comp_stride = comp_stride; 1110 (*rstr)->l_size = l_size; 1111 (*rstr)->e_size = (CeedSize)num_block * (CeedSize)block_size * (CeedSize)elem_size * (CeedSize)num_comp; 1112 (*rstr)->num_block = num_block; 1113 (*rstr)->block_size = block_size; 1114 (*rstr)->rstr_type = CEED_RESTRICTION_CURL_ORIENTED; 1115 CeedCall(ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, (const CeedInt *)block_offsets, NULL, 1116 (const CeedInt8 *)block_curl_orients, *rstr)); 1117 if (copy_mode == CEED_OWN_POINTER) CeedCall(CeedFree(&offsets)); 1118 return CEED_ERROR_SUCCESS; 1119 } 1120 1121 /** 1122 @brief Create a blocked strided `CeedElemRestriction`, typically only used by backends 1123 1124 @param[in] ceed `Ceed` context used to create the `CeedElemRestriction` 1125 @param[in] num_elem Number of elements described by the restriction 1126 @param[in] elem_size Size (number of "nodes") per element 1127 @param[in] block_size Number of elements in a block 1128 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 1129 @param[in] l_size The size of the L-vector. 1130 This vector may be larger than the elements and fields given by this restriction. 1131 @param[in] strides Array for strides between `[nodes, components, elements]`. 1132 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]`. 1133 @ref CEED_STRIDES_BACKEND may be used for `CeedVector` ordered by the same `Ceed` backend. 1134 @param[out] rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 1135 1136 @return An error code: 0 - success, otherwise - failure 1137 1138 @ref User 1139 **/ 1140 int CeedElemRestrictionCreateBlockedStrided(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, CeedSize l_size, 1141 const CeedInt strides[3], CeedElemRestriction *rstr) { 1142 CeedInt num_block = (num_elem / block_size) + !!(num_elem % block_size); 1143 1144 if (!ceed->ElemRestrictionCreateBlocked) { 1145 Ceed delegate; 1146 1147 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 1148 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedElemRestrictionCreateBlockedStrided"); 1149 CeedCall(CeedElemRestrictionCreateBlockedStrided(delegate, num_elem, elem_size, block_size, num_comp, l_size, strides, rstr)); 1150 CeedCall(CeedDestroy(&delegate)); 1151 return CEED_ERROR_SUCCESS; 1152 } 1153 1154 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 1155 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 1156 CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); 1157 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedElemRestriction must have at least 1 component"); 1158 CeedCheck(l_size >= (CeedSize)num_elem * (CeedSize)elem_size * (CeedSize)num_comp, ceed, CEED_ERROR_DIMENSION, 1159 "L-vector size must be at least num_elem * elem_size * num_comp. Expected: > %" CeedSize_FMT " Found: %" CeedSize_FMT, 1160 (CeedSize)num_elem * (CeedSize)elem_size * (CeedSize)num_comp, l_size); 1161 1162 CeedCall(CeedCalloc(1, rstr)); 1163 CeedCall(CeedObjectCreate(ceed, CeedElemRestrictionView_Object, &(*rstr)->obj)); 1164 (*rstr)->num_elem = num_elem; 1165 (*rstr)->elem_size = elem_size; 1166 (*rstr)->num_comp = num_comp; 1167 (*rstr)->l_size = l_size; 1168 (*rstr)->e_size = (CeedSize)num_block * (CeedSize)block_size * (CeedSize)elem_size * (CeedSize)num_comp; 1169 (*rstr)->num_block = num_block; 1170 (*rstr)->block_size = block_size; 1171 (*rstr)->rstr_type = CEED_RESTRICTION_STRIDED; 1172 CeedCall(CeedMalloc(3, &(*rstr)->strides)); 1173 for (CeedInt i = 0; i < 3; i++) (*rstr)->strides[i] = strides[i]; 1174 CeedCall(ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, NULL, NULL, NULL, *rstr)); 1175 return CEED_ERROR_SUCCESS; 1176 } 1177 1178 /** 1179 @brief Copy the pointer to a `CeedElemRestriction` and set @ref CeedElemRestrictionApply() implementation to use the unsigned version. 1180 1181 Both pointers should be destroyed with @ref CeedElemRestrictionDestroy(). 1182 1183 @param[in] rstr `CeedElemRestriction` to create unsigned reference to 1184 @param[in,out] rstr_unsigned Variable to store unsigned `CeedElemRestriction` 1185 1186 @return An error code: 0 - success, otherwise - failure 1187 1188 @ref User 1189 **/ 1190 int CeedElemRestrictionCreateUnsignedCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_unsigned) { 1191 CeedCall(CeedCalloc(1, rstr_unsigned)); 1192 1193 // Copy old rstr 1194 memcpy(*rstr_unsigned, rstr, sizeof(struct CeedElemRestriction_private)); 1195 CeedCall(CeedObjectCreate(CeedElemRestrictionReturnCeed(rstr), CeedElemRestrictionView_Object, &(*rstr_unsigned)->obj)); 1196 (*rstr_unsigned)->strides = NULL; 1197 if (rstr->strides) { 1198 CeedCall(CeedMalloc(3, &(*rstr_unsigned)->strides)); 1199 for (CeedInt i = 0; i < 3; i++) (*rstr_unsigned)->strides[i] = rstr->strides[i]; 1200 } 1201 CeedCall(CeedElemRestrictionReferenceCopy(rstr, &(*rstr_unsigned)->rstr_base)); 1202 1203 // Override Apply 1204 (*rstr_unsigned)->Apply = rstr->ApplyUnsigned; 1205 return CEED_ERROR_SUCCESS; 1206 } 1207 1208 /** 1209 @brief Copy the pointer to a `CeedElemRestriction` and set @ref CeedElemRestrictionApply() implementation to use the unoriented version. 1210 1211 Both pointers should be destroyed with @ref CeedElemRestrictionDestroy(). 1212 1213 @param[in] rstr `CeedElemRestriction` to create unoriented reference to 1214 @param[in,out] rstr_unoriented Variable to store unoriented `CeedElemRestriction` 1215 1216 @return An error code: 0 - success, otherwise - failure 1217 1218 @ref User 1219 **/ 1220 int CeedElemRestrictionCreateUnorientedCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_unoriented) { 1221 CeedCall(CeedCalloc(1, rstr_unoriented)); 1222 1223 // Copy old rstr 1224 memcpy(*rstr_unoriented, rstr, sizeof(struct CeedElemRestriction_private)); 1225 CeedCall(CeedObjectCreate(CeedElemRestrictionReturnCeed(rstr), CeedElemRestrictionView_Object, &(*rstr_unoriented)->obj)); 1226 (*rstr_unoriented)->strides = NULL; 1227 if (rstr->strides) { 1228 CeedCall(CeedMalloc(3, &(*rstr_unoriented)->strides)); 1229 for (CeedInt i = 0; i < 3; i++) (*rstr_unoriented)->strides[i] = rstr->strides[i]; 1230 } 1231 CeedCall(CeedElemRestrictionReferenceCopy(rstr, &(*rstr_unoriented)->rstr_base)); 1232 1233 // Override Apply 1234 (*rstr_unoriented)->Apply = rstr->ApplyUnoriented; 1235 return CEED_ERROR_SUCCESS; 1236 } 1237 1238 /** 1239 @brief Copy the pointer to a `CeedElemRestriction`. 1240 1241 Both pointers should be destroyed with @ref CeedElemRestrictionDestroy(). 1242 1243 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`. 1244 This `CeedElemRestriction` will be destroyed if `*rstr_copy` is the only reference to this `CeedElemRestriction`. 1245 1246 @param[in] rstr `CeedElemRestriction` to copy reference to 1247 @param[in,out] rstr_copy Variable to store copied reference 1248 1249 @return An error code: 0 - success, otherwise - failure 1250 1251 @ref User 1252 **/ 1253 int CeedElemRestrictionReferenceCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_copy) { 1254 if (rstr != CEED_ELEMRESTRICTION_NONE) CeedCall(CeedElemRestrictionReference(rstr)); 1255 CeedCall(CeedElemRestrictionDestroy(rstr_copy)); 1256 *rstr_copy = rstr; 1257 return CEED_ERROR_SUCCESS; 1258 } 1259 1260 /** 1261 @brief Create `CeedVector` associated with a `CeedElemRestriction` 1262 1263 @param[in] rstr `CeedElemRestriction` 1264 @param[out] l_vec The address of the L-vector to be created, or `NULL` 1265 @param[out] e_vec The address of the E-vector to be created, or `NULL` 1266 1267 @return An error code: 0 - success, otherwise - failure 1268 1269 @ref User 1270 **/ 1271 int CeedElemRestrictionCreateVector(CeedElemRestriction rstr, CeedVector *l_vec, CeedVector *e_vec) { 1272 CeedSize e_size, l_size; 1273 Ceed ceed; 1274 1275 CeedCall(CeedElemRestrictionGetCeed(rstr, &ceed)); 1276 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &l_size)); 1277 CeedCall(CeedElemRestrictionGetEVectorSize(rstr, &e_size)); 1278 if (l_vec) CeedCall(CeedVectorCreate(ceed, l_size, l_vec)); 1279 if (e_vec) CeedCall(CeedVectorCreate(ceed, e_size, e_vec)); 1280 CeedCall(CeedDestroy(&ceed)); 1281 return CEED_ERROR_SUCCESS; 1282 } 1283 1284 /** 1285 @brief Restrict an L-vector to an E-vector or apply its transpose 1286 1287 @param[in] rstr `CeedElemRestriction` 1288 @param[in] t_mode Apply restriction or transpose 1289 @param[in] u Input vector (of size `l_size` when `t_mode` = @ref CEED_NOTRANSPOSE) 1290 @param[out] ru Output vector (of shape `[num_elem * elem_size]` when `t_mode` = @ref CEED_NOTRANSPOSE). 1291 Ordering of the e-vector is decided by the backend. 1292 @param[in] request Request or @ref CEED_REQUEST_IMMEDIATE 1293 1294 @return An error code: 0 - success, otherwise - failure 1295 1296 @ref User 1297 **/ 1298 int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, CeedRequest *request) { 1299 CeedSize min_u_len, min_ru_len, len; 1300 CeedInt num_elem; 1301 1302 if (t_mode == CEED_NOTRANSPOSE) { 1303 CeedCall(CeedElemRestrictionGetEVectorSize(rstr, &min_ru_len)); 1304 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &min_u_len)); 1305 } else { 1306 CeedCall(CeedElemRestrictionGetEVectorSize(rstr, &min_u_len)); 1307 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &min_ru_len)); 1308 } 1309 CeedCall(CeedVectorGetLength(u, &len)); 1310 CeedCheck(min_u_len <= len, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_DIMENSION, 1311 "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", len, min_ru_len, 1312 min_u_len); 1313 CeedCall(CeedVectorGetLength(ru, &len)); 1314 CeedCheck(min_ru_len <= len, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_DIMENSION, 1315 "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", len, min_u_len, 1316 min_ru_len); 1317 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 1318 if (num_elem > 0) CeedCall(rstr->Apply(rstr, t_mode, u, ru, request)); 1319 return CEED_ERROR_SUCCESS; 1320 } 1321 1322 /** 1323 @brief Restrict an L-vector of points to a single element or apply its transpose 1324 1325 @param[in] rstr `CeedElemRestriction` 1326 @param[in] elem Element number in range `[0, num_elem)` 1327 @param[in] t_mode Apply restriction or transpose 1328 @param[in] u Input vector (of size `l_size` when `t_mode` = @ref CEED_NOTRANSPOSE) 1329 @param[out] ru Output vector (of shape [`num_points * num_comp]` when `t_mode` = @ref CEED_NOTRANSPOSE). 1330 Ordering of the e-vector is decided by the backend. 1331 @param[in] request Request or @ref CEED_REQUEST_IMMEDIATE 1332 1333 @return An error code: 0 - success, otherwise - failure 1334 1335 @ref User 1336 **/ 1337 int CeedElemRestrictionApplyAtPointsInElement(CeedElemRestriction rstr, CeedInt elem, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, 1338 CeedRequest *request) { 1339 CeedSize min_u_len, min_ru_len, len; 1340 CeedInt num_elem; 1341 1342 CeedCheck(rstr->ApplyAtPointsInElement, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_UNSUPPORTED, 1343 "Backend does not implement CeedElemRestrictionApplyAtPointsInElement"); 1344 1345 if (t_mode == CEED_NOTRANSPOSE) { 1346 CeedInt num_points, num_comp; 1347 1348 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &min_u_len)); 1349 CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr, elem, &num_points)); 1350 CeedCall(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 1351 min_ru_len = (CeedSize)num_points * (CeedSize)num_comp; 1352 } else { 1353 CeedInt num_points, num_comp; 1354 1355 CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr, elem, &num_points)); 1356 CeedCall(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 1357 min_u_len = (CeedSize)num_points * (CeedSize)num_comp; 1358 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &min_ru_len)); 1359 } 1360 CeedCall(CeedVectorGetLength(u, &len)); 1361 CeedCheck(min_u_len <= len, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_DIMENSION, 1362 "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT 1363 ") for element %" CeedInt_FMT, 1364 len, min_ru_len, min_u_len, elem); 1365 CeedCall(CeedVectorGetLength(ru, &len)); 1366 CeedCheck(min_ru_len <= len, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_DIMENSION, 1367 "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT 1368 ") for element %" CeedInt_FMT, 1369 len, min_ru_len, min_u_len, elem); 1370 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 1371 CeedCheck(elem < num_elem, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_DIMENSION, 1372 "Cannot retrieve element %" CeedInt_FMT ", element %" CeedInt_FMT " > total elements %" CeedInt_FMT "", elem, elem, num_elem); 1373 if (num_elem > 0) CeedCall(rstr->ApplyAtPointsInElement(rstr, elem, t_mode, u, ru, request)); 1374 return CEED_ERROR_SUCCESS; 1375 } 1376 1377 /** 1378 @brief Restrict an L-vector to a block of an E-vector or apply its transpose 1379 1380 @param[in] rstr `CeedElemRestriction` 1381 @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]` 1382 @param[in] t_mode Apply restriction or transpose 1383 @param[in] u Input vector (of size `l_size` when `t_mode` = @ref CEED_NOTRANSPOSE) 1384 @param[out] ru Output vector (of shape `[block_size * elem_size]` when `t_mode` = @ref CEED_NOTRANSPOSE). 1385 Ordering of the e-vector is decided by the backend. 1386 @param[in] request Request or @ref CEED_REQUEST_IMMEDIATE 1387 1388 @return An error code: 0 - success, otherwise - failure 1389 1390 @ref Backend 1391 **/ 1392 int CeedElemRestrictionApplyBlock(CeedElemRestriction rstr, CeedInt block, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, 1393 CeedRequest *request) { 1394 CeedSize min_u_len, min_ru_len, len; 1395 CeedInt block_size, num_elem; 1396 1397 CeedCheck(rstr->ApplyBlock, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_UNSUPPORTED, 1398 "Backend does not implement CeedElemRestrictionApplyBlock"); 1399 1400 CeedCall(CeedElemRestrictionGetBlockSize(rstr, &block_size)); 1401 if (t_mode == CEED_NOTRANSPOSE) { 1402 CeedInt elem_size, num_comp; 1403 1404 CeedCall(CeedElemRestrictionGetElementSize(rstr, &elem_size)); 1405 CeedCall(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 1406 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &min_u_len)); 1407 min_ru_len = (CeedSize)block_size * (CeedSize)elem_size * (CeedSize)num_comp; 1408 } else { 1409 CeedInt elem_size, num_comp; 1410 1411 CeedCall(CeedElemRestrictionGetElementSize(rstr, &elem_size)); 1412 CeedCall(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 1413 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &min_ru_len)); 1414 min_u_len = (CeedSize)block_size * (CeedSize)elem_size * (CeedSize)num_comp; 1415 } 1416 CeedCall(CeedVectorGetLength(u, &len)); 1417 CeedCheck(min_u_len == len, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_DIMENSION, 1418 "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", len, min_u_len, 1419 min_ru_len); 1420 CeedCall(CeedVectorGetLength(ru, &len)); 1421 CeedCheck(min_ru_len == len, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_DIMENSION, 1422 "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", len, min_ru_len, 1423 min_u_len); 1424 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 1425 CeedCheck(block_size * block <= num_elem, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_DIMENSION, 1426 "Cannot retrieve block %" CeedInt_FMT ", element %" CeedInt_FMT " > total elements %" CeedInt_FMT "", block, block_size * block, 1427 num_elem); 1428 CeedCall(rstr->ApplyBlock(rstr, block, t_mode, u, ru, request)); 1429 return CEED_ERROR_SUCCESS; 1430 } 1431 1432 /** 1433 @brief Get the `Ceed` associated with a `CeedElemRestriction` 1434 1435 @param[in] rstr `CeedElemRestriction` 1436 @param[out] ceed Variable to store `Ceed` 1437 1438 @return An error code: 0 - success, otherwise - failure 1439 1440 @ref Advanced 1441 **/ 1442 int CeedElemRestrictionGetCeed(CeedElemRestriction rstr, Ceed *ceed) { 1443 CeedCall(CeedObjectGetCeed((CeedObject)rstr, ceed)); 1444 return CEED_ERROR_SUCCESS; 1445 } 1446 1447 /** 1448 @brief Return the `Ceed` associated with a `CeedElemRestriction` 1449 1450 @param[in] rstr `CeedElemRestriction` 1451 1452 @return `Ceed` associated with the `rstr` 1453 1454 @ref Advanced 1455 **/ 1456 Ceed CeedElemRestrictionReturnCeed(CeedElemRestriction rstr) { return CeedObjectReturnCeed((CeedObject)rstr); } 1457 1458 /** 1459 @brief Get the L-vector component stride 1460 1461 @param[in] rstr `CeedElemRestriction` 1462 @param[out] comp_stride Variable to store component stride 1463 1464 @return An error code: 0 - success, otherwise - failure 1465 1466 @ref Advanced 1467 **/ 1468 int CeedElemRestrictionGetCompStride(CeedElemRestriction rstr, CeedInt *comp_stride) { 1469 *comp_stride = rstr->comp_stride; 1470 return CEED_ERROR_SUCCESS; 1471 } 1472 1473 /** 1474 @brief Get the total number of elements in the range of a `CeedElemRestriction` 1475 1476 @param[in] rstr `CeedElemRestriction` 1477 @param[out] num_elem Variable to store number of elements 1478 1479 @return An error code: 0 - success, otherwise - failure 1480 1481 @ref Advanced 1482 **/ 1483 int CeedElemRestrictionGetNumElements(CeedElemRestriction rstr, CeedInt *num_elem) { 1484 *num_elem = rstr->num_elem; 1485 return CEED_ERROR_SUCCESS; 1486 } 1487 1488 /** 1489 @brief Get the size of elements in the `CeedElemRestriction` 1490 1491 @param[in] rstr `CeedElemRestriction` 1492 @param[out] elem_size Variable to store size of elements 1493 1494 @return An error code: 0 - success, otherwise - failure 1495 1496 @ref Advanced 1497 **/ 1498 int CeedElemRestrictionGetElementSize(CeedElemRestriction rstr, CeedInt *elem_size) { 1499 *elem_size = rstr->elem_size; 1500 return CEED_ERROR_SUCCESS; 1501 } 1502 1503 /** 1504 1505 @brief Get the number of points in the offsets array for a points `CeedElemRestriction` 1506 1507 @param[in] rstr `CeedElemRestriction` 1508 @param[out] num_points The number of points in the offsets array 1509 1510 @return An error code: 0 - success, otherwise - failure 1511 1512 @ref User 1513 **/ 1514 int CeedElemRestrictionGetNumPoints(CeedElemRestriction rstr, CeedInt *num_points) { 1515 CeedRestrictionType rstr_type; 1516 1517 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 1518 CeedCheck(rstr_type == CEED_RESTRICTION_POINTS, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_INCOMPATIBLE, 1519 "Can only retrieve the number of points for a points CeedElemRestriction"); 1520 1521 *num_points = rstr->num_points; 1522 return CEED_ERROR_SUCCESS; 1523 } 1524 1525 /** 1526 1527 @brief Get the number of points in an element of a `CeedElemRestriction` at points 1528 1529 @param[in] rstr `CeedElemRestriction` 1530 @param[in] elem Index number of element to retrieve the number of points for 1531 @param[out] num_points The number of points in the element at index elem 1532 1533 @return An error code: 0 - success, otherwise - failure 1534 1535 @ref User 1536 **/ 1537 int CeedElemRestrictionGetNumPointsInElement(CeedElemRestriction rstr, CeedInt elem, CeedInt *num_points) { 1538 CeedRestrictionType rstr_type; 1539 const CeedInt *offsets; 1540 1541 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 1542 CeedCheck(rstr_type == CEED_RESTRICTION_POINTS, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_INCOMPATIBLE, 1543 "Can only retrieve the number of points for a points CeedElemRestriction"); 1544 1545 CeedCall(CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets)); 1546 *num_points = offsets[elem + 1] - offsets[elem]; 1547 CeedCall(CeedElemRestrictionRestoreOffsets(rstr, &offsets)); 1548 return CEED_ERROR_SUCCESS; 1549 } 1550 1551 /** 1552 @brief Get the minimum and/or maximum number of points in an element for a `CeedElemRestriction` at points 1553 1554 @param[in] rstr `CeedElemRestriction` 1555 @param[out] min_points Variable to minimum number of points in an element, or `NULL` 1556 @param[out] max_points Variable to maximum number of points in an element, or `NULL` 1557 1558 @return An error code: 0 - success, otherwise - failure 1559 1560 @ref Advanced 1561 **/ 1562 int CeedElemRestrictionGetMinMaxPointsInElement(CeedElemRestriction rstr, CeedInt *min_points, CeedInt *max_points) { 1563 CeedInt num_elem, num_points; 1564 CeedRestrictionType rstr_type; 1565 1566 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 1567 CeedCheck(rstr_type == CEED_RESTRICTION_POINTS, CeedElemRestrictionReturnCeed(rstr), CEED_ERROR_INCOMPATIBLE, 1568 "Cannot compute min/max points for a CeedElemRestriction that does not use points"); 1569 1570 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 1571 1572 // Exit early if there are no elements 1573 if (num_elem == 0) { 1574 if (min_points) *min_points = 0; 1575 if (max_points) *max_points = 0; 1576 return CEED_ERROR_SUCCESS; 1577 } 1578 1579 // Initialize to the number of points in the first element 1580 CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr, 0, &num_points)); 1581 if (min_points) *min_points = num_points; 1582 if (max_points) *max_points = num_points; 1583 for (CeedInt e = 1; e < num_elem; e++) { 1584 CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr, e, &num_points)); 1585 if (min_points) *min_points = CeedIntMin(num_points, *min_points); 1586 if (max_points) *max_points = CeedIntMax(num_points, *max_points); 1587 } 1588 return CEED_ERROR_SUCCESS; 1589 } 1590 1591 /** 1592 @brief Get the maximum number of points in an element for a `CeedElemRestriction` at points 1593 1594 @param[in] rstr `CeedElemRestriction` 1595 @param[out] max_points Variable to store maximum number of points in an element 1596 1597 @return An error code: 0 - success, otherwise - failure 1598 1599 @ref User 1600 1601 @see CeedElemRestrictionGetMinMaxPointsInElement() 1602 **/ 1603 int CeedElemRestrictionGetMaxPointsInElement(CeedElemRestriction rstr, CeedInt *max_points) { 1604 return CeedElemRestrictionGetMinMaxPointsInElement(rstr, NULL, max_points); 1605 } 1606 1607 /** 1608 @brief Get the minimum number of points in an element for a `CeedElemRestriction` at points 1609 1610 @param[in] rstr `CeedElemRestriction` 1611 @param[out] min_points Variable to store minimum 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 CeedElemRestrictionGetMinPointsInElement(CeedElemRestriction rstr, CeedInt *min_points) { 1620 return CeedElemRestrictionGetMinMaxPointsInElement(rstr, min_points, NULL); 1621 } 1622 1623 /** 1624 @brief Get the size of the l-vector for a `CeedElemRestriction` 1625 1626 @param[in] rstr `CeedElemRestriction` 1627 @param[out] l_size Variable to store l-vector size 1628 1629 @return An error code: 0 - success, otherwise - failure 1630 1631 @ref Advanced 1632 **/ 1633 int CeedElemRestrictionGetLVectorSize(CeedElemRestriction rstr, CeedSize *l_size) { 1634 *l_size = rstr->l_size; 1635 return CEED_ERROR_SUCCESS; 1636 } 1637 1638 /** 1639 @brief Get the size of the e-vector for a `CeedElemRestriction` 1640 1641 @param[in] rstr `CeedElemRestriction` 1642 @param[out] e_size Variable to store e-vector size 1643 1644 @return An error code: 0 - success, otherwise - failure 1645 1646 @ref Advanced 1647 **/ 1648 int CeedElemRestrictionGetEVectorSize(CeedElemRestriction rstr, CeedSize *e_size) { 1649 *e_size = rstr->e_size; 1650 return CEED_ERROR_SUCCESS; 1651 } 1652 1653 /** 1654 @brief Get the number of components in the elements of a `CeedElemRestriction` 1655 1656 @param[in] rstr `CeedElemRestriction` 1657 @param[out] num_comp Variable to store number of components 1658 1659 @return An error code: 0 - success, otherwise - failure 1660 1661 @ref Advanced 1662 **/ 1663 int CeedElemRestrictionGetNumComponents(CeedElemRestriction rstr, CeedInt *num_comp) { 1664 *num_comp = rstr->num_comp; 1665 return CEED_ERROR_SUCCESS; 1666 } 1667 1668 /** 1669 @brief Get the number of blocks in a `CeedElemRestriction` 1670 1671 @param[in] rstr `CeedElemRestriction` 1672 @param[out] num_block Variable to store number of blocks 1673 1674 @return An error code: 0 - success, otherwise - failure 1675 1676 @ref Advanced 1677 **/ 1678 int CeedElemRestrictionGetNumBlocks(CeedElemRestriction rstr, CeedInt *num_block) { 1679 *num_block = rstr->num_block; 1680 return CEED_ERROR_SUCCESS; 1681 } 1682 1683 /** 1684 @brief Get the size of blocks in the `CeedElemRestriction` 1685 1686 @param[in] rstr `CeedElemRestriction` 1687 @param[out] block_size Variable to store size of blocks 1688 1689 @return An error code: 0 - success, otherwise - failure 1690 1691 @ref Advanced 1692 **/ 1693 int CeedElemRestrictionGetBlockSize(CeedElemRestriction rstr, CeedInt *block_size) { 1694 *block_size = rstr->block_size; 1695 return CEED_ERROR_SUCCESS; 1696 } 1697 1698 /** 1699 @brief Get the multiplicity of nodes in a `CeedElemRestriction` 1700 1701 @param[in] rstr `CeedElemRestriction` 1702 @param[out] mult Vector to store multiplicity (of size `l_size`) 1703 1704 @return An error code: 0 - success, otherwise - failure 1705 1706 @ref User 1707 **/ 1708 int CeedElemRestrictionGetMultiplicity(CeedElemRestriction rstr, CeedVector mult) { 1709 CeedVector e_vec; 1710 1711 // Create e_vec to hold intermediate computation in E^T (E 1) 1712 CeedCall(CeedElemRestrictionCreateVector(rstr, NULL, &e_vec)); 1713 1714 // Compute e_vec = E * 1 1715 CeedCall(CeedVectorSetValue(mult, 1.0)); 1716 CeedCall(CeedElemRestrictionApply(rstr, CEED_NOTRANSPOSE, mult, e_vec, CEED_REQUEST_IMMEDIATE)); 1717 // Compute multiplicity, mult = E^T * e_vec = E^T (E 1) 1718 CeedCall(CeedVectorSetValue(mult, 0.0)); 1719 CeedCall(CeedElemRestrictionApply(rstr, CEED_TRANSPOSE, e_vec, mult, CEED_REQUEST_IMMEDIATE)); 1720 // Cleanup 1721 CeedCall(CeedVectorDestroy(&e_vec)); 1722 return CEED_ERROR_SUCCESS; 1723 } 1724 1725 /** 1726 @brief Set the number of tabs to indent for @ref CeedElemRestrictionView() output 1727 1728 @param[in] rstr `CeedElemRestriction` to set the number of view tabs 1729 @param[in] num_tabs Number of view tabs to set 1730 1731 @return Error code: 0 - success, otherwise - failure 1732 1733 @ref User 1734 **/ 1735 int CeedElemRestrictionSetNumViewTabs(CeedElemRestriction rstr, CeedInt num_tabs) { 1736 CeedCall(CeedObjectSetNumViewTabs((CeedObject)rstr, num_tabs)); 1737 return CEED_ERROR_SUCCESS; 1738 } 1739 1740 /** 1741 @brief Get the number of tabs to indent for @ref CeedElemRestrictionView() output 1742 1743 @param[in] rstr `CeedElemRestriction` to get the number of view tabs 1744 @param[out] num_tabs Number of view tabs 1745 1746 @return Error code: 0 - success, otherwise - failure 1747 1748 @ref User 1749 **/ 1750 int CeedElemRestrictionGetNumViewTabs(CeedElemRestriction rstr, CeedInt *num_tabs) { 1751 CeedCall(CeedObjectGetNumViewTabs((CeedObject)rstr, 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 || CeedObjectDereference((CeedObject)*rstr) > 0) { 1815 *rstr = NULL; 1816 return CEED_ERROR_SUCCESS; 1817 } 1818 CeedCheck((*rstr)->num_readers == 0, CeedElemRestrictionReturnCeed(*rstr), 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(CeedObjectDestroy(&(*rstr)->obj)); 1827 CeedCall(CeedFree(rstr)); 1828 return CEED_ERROR_SUCCESS; 1829 } 1830 1831 /// @} 1832