1 // Copyright (c) 2017-2022, 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 restriction 26 27 @param[in] offsets Array of shape [@a num_elem, @a elem_size]. 28 @param[out] block_offsets Array of permuted and padded array values of shape [@a num_block, @a elem_size, @a 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 restriction 52 53 @param[in] orients Array of shape [@a num_elem, @a elem_size]. 54 @param[out] block_orients Array of permuted and padded array values of shape [@a num_block, @a elem_size, @a 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 restriction 77 78 @param[in] curl_orients Array of shape [@a num_elem, @a 3 * elem_size]. 79 @param[out] block_curl_orients Array of permuted and padded array values of shape [@a num_block, @a elem_size, @a 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 strided status of a CeedElemRestriction 126 127 @param[in] rstr CeedElemRestriction 128 @param[out] is_strided Variable to store strided status, 1 if strided else 0 129 **/ 130 int CeedElemRestrictionIsStrided(CeedElemRestriction rstr, bool *is_strided) { 131 *is_strided = (rstr->rstr_type == CEED_RESTRICTION_STRIDED); 132 return CEED_ERROR_SUCCESS; 133 } 134 135 /** 136 @brief Get the points status of a CeedElemRestriction 137 138 @param[in] rstr CeedElemRestriction 139 @param[out] is_points Variable to store points status, 1 if points else 0 140 **/ 141 int CeedElemRestrictionIsPoints(CeedElemRestriction rstr, bool *is_points) { 142 *is_points = (rstr->rstr_type == CEED_RESTRICTION_POINTS); 143 return CEED_ERROR_SUCCESS; 144 } 145 146 /** 147 @brief Get the strides of a strided CeedElemRestriction 148 149 @param[in] rstr CeedElemRestriction 150 @param[out] strides Variable to store strides array 151 152 @return An error code: 0 - success, otherwise - failure 153 154 @ref Backend 155 **/ 156 int CeedElemRestrictionGetStrides(CeedElemRestriction rstr, CeedInt (*strides)[3]) { 157 CeedCheck(rstr->strides, rstr->ceed, CEED_ERROR_MINOR, "ElemRestriction has no stride data"); 158 for (CeedInt i = 0; i < 3; i++) (*strides)[i] = rstr->strides[i]; 159 return CEED_ERROR_SUCCESS; 160 } 161 162 /** 163 @brief Get the backend stride status of a CeedElemRestriction 164 165 @param[in] rstr CeedElemRestriction 166 @param[out] has_backend_strides Variable to store stride status 167 168 @return An error code: 0 - success, otherwise - failure 169 170 @ref Backend 171 **/ 172 int CeedElemRestrictionHasBackendStrides(CeedElemRestriction rstr, bool *has_backend_strides) { 173 CeedCheck(rstr->strides, rstr->ceed, CEED_ERROR_MINOR, "ElemRestriction has no stride data"); 174 *has_backend_strides = ((rstr->strides[0] == CEED_STRIDES_BACKEND[0]) && (rstr->strides[1] == CEED_STRIDES_BACKEND[1]) && 175 (rstr->strides[2] == CEED_STRIDES_BACKEND[2])); 176 return CEED_ERROR_SUCCESS; 177 } 178 179 /** 180 @brief Get read-only access to a CeedElemRestriction offsets array by memtype 181 182 @param[in] rstr CeedElemRestriction to retrieve offsets 183 @param[in] mem_type Memory type on which to access the array. 184 If the backend uses a different memory type, this will perform a copy (possibly cached). 185 @param[out] offsets Array on memory type mem_type 186 187 @return An error code: 0 - success, otherwise - failure 188 189 @ref User 190 **/ 191 int CeedElemRestrictionGetOffsets(CeedElemRestriction rstr, CeedMemType mem_type, const CeedInt **offsets) { 192 if (rstr->rstr_base) { 193 CeedCall(CeedElemRestrictionGetOffsets(rstr->rstr_base, mem_type, offsets)); 194 } else { 195 CeedCheck(rstr->GetOffsets, rstr->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support GetOffsets"); 196 CeedCall(rstr->GetOffsets(rstr, mem_type, offsets)); 197 rstr->num_readers++; 198 } 199 return CEED_ERROR_SUCCESS; 200 } 201 202 /** 203 @brief Restore an offsets array obtained using CeedElemRestrictionGetOffsets() 204 205 @param[in] rstr CeedElemRestriction to restore 206 @param[in] offsets Array of offset data 207 208 @return An error code: 0 - success, otherwise - failure 209 210 @ref User 211 **/ 212 int CeedElemRestrictionRestoreOffsets(CeedElemRestriction rstr, const CeedInt **offsets) { 213 if (rstr->rstr_base) { 214 CeedCall(CeedElemRestrictionRestoreOffsets(rstr->rstr_base, offsets)); 215 } else { 216 *offsets = NULL; 217 rstr->num_readers--; 218 } 219 return CEED_ERROR_SUCCESS; 220 } 221 222 /** 223 @brief Get read-only access to a CeedElemRestriction orientations array by memtype 224 225 @param[in] rstr CeedElemRestriction to retrieve orientations 226 @param[in] mem_type Memory type on which to access the array. 227 If the backend uses a different memory type, this will perform a copy (possibly cached). 228 @param[out] orients Array on memory type mem_type 229 230 @return An error code: 0 - success, otherwise - failure 231 232 @ref User 233 **/ 234 int CeedElemRestrictionGetOrientations(CeedElemRestriction rstr, CeedMemType mem_type, const bool **orients) { 235 CeedCheck(rstr->GetOrientations, rstr->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support GetOrientations"); 236 CeedCall(rstr->GetOrientations(rstr, mem_type, orients)); 237 rstr->num_readers++; 238 return CEED_ERROR_SUCCESS; 239 } 240 241 /** 242 @brief Restore an orientations array obtained using CeedElemRestrictionGetOrientations() 243 244 @param[in] rstr CeedElemRestriction to restore 245 @param[in] orients Array of orientation data 246 247 @return An error code: 0 - success, otherwise - failure 248 249 @ref User 250 **/ 251 int CeedElemRestrictionRestoreOrientations(CeedElemRestriction rstr, const bool **orients) { 252 *orients = NULL; 253 rstr->num_readers--; 254 return CEED_ERROR_SUCCESS; 255 } 256 257 /** 258 @brief Get read-only access to a CeedElemRestriction curl-conforming orientations array by memtype 259 260 @param[in] rstr CeedElemRestriction to retrieve curl-conforming orientations 261 @param[in] mem_type Memory type on which to access the array. 262 If the backend uses a different memory type, this will perform a copy (possibly cached). 263 @param[out] curl_orients Array on memory type mem_type 264 265 @return An error code: 0 - success, otherwise - failure 266 267 @ref User 268 **/ 269 int CeedElemRestrictionGetCurlOrientations(CeedElemRestriction rstr, CeedMemType mem_type, const CeedInt8 **curl_orients) { 270 CeedCheck(rstr->GetCurlOrientations, rstr->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support GetCurlOrientations"); 271 CeedCall(rstr->GetCurlOrientations(rstr, mem_type, curl_orients)); 272 rstr->num_readers++; 273 return CEED_ERROR_SUCCESS; 274 } 275 276 /** 277 @brief Restore an orientations array obtained using CeedElemRestrictionGetCurlOrientations() 278 279 @param[in] rstr CeedElemRestriction to restore 280 @param[in] curl_orients Array of orientation data 281 282 @return An error code: 0 - success, otherwise - failure 283 284 @ref User 285 **/ 286 int CeedElemRestrictionRestoreCurlOrientations(CeedElemRestriction rstr, const CeedInt8 **curl_orients) { 287 *curl_orients = NULL; 288 rstr->num_readers--; 289 return CEED_ERROR_SUCCESS; 290 } 291 292 /** 293 294 @brief Get the number of points in an element of a points CeedElemRestriction 295 296 @param[in] rstr CeedElemRestriction 297 @param[in] elem Index number of element to retrieve the number of points for 298 @param[out] num_points The number of points in the element at index elem 299 300 @return An error code: 0 - success, otherwise - failure 301 302 @ref Backend 303 **/ 304 int CeedElemRestrictionGetNumPointsInElement(CeedElemRestriction rstr, CeedInt elem, CeedInt *num_points) { 305 Ceed ceed; 306 const CeedInt *offsets; 307 308 CeedCall(CeedElemRestrictionGetCeed(rstr, &ceed)); 309 CeedCheck(rstr->rstr_type == CEED_RESTRICTION_POINTS, ceed, CEED_ERROR_INCOMPATIBLE, 310 "Can only retrieve the number of points for a points CeedElemRestriction"); 311 312 CeedCall(CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets)); 313 *num_points = offsets[elem - 1] - offsets[elem]; 314 CeedCall(CeedElemRestrictionRestoreOffsets(rstr, &offsets)); 315 return CEED_ERROR_SUCCESS; 316 } 317 318 /** 319 320 @brief Get the E-vector layout of a CeedElemRestriction 321 322 @param[in] rstr CeedElemRestriction 323 @param[out] layout Variable to store layout array, stored as [nodes, components, elements]. 324 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] 325 326 @return An error code: 0 - success, otherwise - failure 327 328 @ref Backend 329 **/ 330 int CeedElemRestrictionGetELayout(CeedElemRestriction rstr, CeedInt (*layout)[3]) { 331 CeedCheck(rstr->layout[0], rstr->ceed, CEED_ERROR_MINOR, "ElemRestriction has no layout data"); 332 for (CeedInt i = 0; i < 3; i++) (*layout)[i] = rstr->layout[i]; 333 return CEED_ERROR_SUCCESS; 334 } 335 336 /** 337 338 @brief Set the E-vector layout of a CeedElemRestriction 339 340 @param[in] rstr CeedElemRestriction 341 @param[in] layout Variable to containing layout array, stored as [nodes, components, elements]. 342 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] 343 344 @return An error code: 0 - success, otherwise - failure 345 346 @ref Backend 347 **/ 348 int CeedElemRestrictionSetELayout(CeedElemRestriction rstr, CeedInt layout[3]) { 349 for (CeedInt i = 0; i < 3; i++) rstr->layout[i] = layout[i]; 350 return CEED_ERROR_SUCCESS; 351 } 352 353 /** 354 @brief Get the backend data of a CeedElemRestriction 355 356 @param[in] rstr CeedElemRestriction 357 @param[out] data Variable to store data 358 359 @return An error code: 0 - success, otherwise - failure 360 361 @ref Backend 362 **/ 363 int CeedElemRestrictionGetData(CeedElemRestriction rstr, void *data) { 364 *(void **)data = rstr->data; 365 return CEED_ERROR_SUCCESS; 366 } 367 368 /** 369 @brief Set the backend data of a CeedElemRestriction 370 371 @param[in,out] rstr CeedElemRestriction 372 @param[in] data Data to set 373 374 @return An error code: 0 - success, otherwise - failure 375 376 @ref Backend 377 **/ 378 int CeedElemRestrictionSetData(CeedElemRestriction rstr, void *data) { 379 rstr->data = data; 380 return CEED_ERROR_SUCCESS; 381 } 382 383 /** 384 @brief Increment the reference counter for a CeedElemRestriction 385 386 @param[in,out] rstr ElemRestriction to increment the reference counter 387 388 @return An error code: 0 - success, otherwise - failure 389 390 @ref Backend 391 **/ 392 int CeedElemRestrictionReference(CeedElemRestriction rstr) { 393 rstr->ref_count++; 394 return CEED_ERROR_SUCCESS; 395 } 396 397 /** 398 @brief Estimate number of FLOPs required to apply CeedElemRestriction in t_mode 399 400 @param[in] rstr ElemRestriction to estimate FLOPs for 401 @param[in] t_mode Apply restriction or transpose 402 @param[out] flops Address of variable to hold FLOPs estimate 403 404 @ref Backend 405 **/ 406 int CeedElemRestrictionGetFlopsEstimate(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedSize *flops) { 407 CeedInt e_size = rstr->num_block * rstr->block_size * rstr->elem_size * rstr->num_comp, scale = 0; 408 CeedRestrictionType rstr_type; 409 410 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 411 if (rstr_type == CEED_RESTRICTION_POINTS) e_size = rstr->num_points * rstr->num_comp; 412 if (t_mode == CEED_TRANSPOSE) { 413 switch (rstr_type) { 414 case CEED_RESTRICTION_POINTS: 415 scale = 0; 416 break; 417 case CEED_RESTRICTION_STRIDED: 418 case CEED_RESTRICTION_STANDARD: 419 scale = 1; 420 break; 421 case CEED_RESTRICTION_ORIENTED: 422 scale = 2; 423 break; 424 case CEED_RESTRICTION_CURL_ORIENTED: 425 scale = 6; 426 break; 427 } 428 } else { 429 switch (rstr_type) { 430 case CEED_RESTRICTION_STRIDED: 431 case CEED_RESTRICTION_STANDARD: 432 case CEED_RESTRICTION_POINTS: 433 scale = 0; 434 break; 435 case CEED_RESTRICTION_ORIENTED: 436 scale = 1; 437 break; 438 case CEED_RESTRICTION_CURL_ORIENTED: 439 scale = 5; 440 break; 441 } 442 } 443 *flops = e_size * scale; 444 return CEED_ERROR_SUCCESS; 445 } 446 447 /// @} 448 449 /// @cond DOXYGEN_SKIP 450 static struct CeedElemRestriction_private ceed_elemrestriction_none; 451 /// @endcond 452 453 /// ---------------------------------------------------------------------------- 454 /// CeedElemRestriction Public API 455 /// ---------------------------------------------------------------------------- 456 /// @addtogroup CeedElemRestrictionUser 457 /// @{ 458 459 /// Indicate that the stride is determined by the backend 460 const CeedInt CEED_STRIDES_BACKEND[3] = {0}; 461 462 /// Argument for CeedOperatorSetField indicating that the field does not requre a CeedElemRestriction 463 const CeedElemRestriction CEED_ELEMRESTRICTION_NONE = &ceed_elemrestriction_none; 464 465 /** 466 @brief Create a CeedElemRestriction 467 468 @param[in] ceed Ceed object where the CeedElemRestriction will be created 469 @param[in] num_elem Number of elements described in the @a offsets array 470 @param[in] elem_size Size (number of "nodes") per element 471 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 472 @param[in] comp_stride Stride between components for the same L-vector "node". 473 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. 474 @param[in] l_size The size of the L-vector. 475 This vector may be larger than the elements and fields given by this restriction. 476 @param[in] mem_type Memory type of the @a offsets array, see CeedMemType 477 @param[in] copy_mode Copy mode for the @a offsets array, see CeedCopyMode 478 @param[in] offsets Array of shape [@a num_elem, @a elem_size]. 479 Row i holds the ordered list of the offsets (into the input CeedVector) for the unknowns corresponding to element i, where 480 0 <= i < @a num_elem. All offsets must be in the range [0, @a l_size - 1]. 481 @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored 482 483 @return An error code: 0 - success, otherwise - failure 484 485 @ref User 486 **/ 487 int CeedElemRestrictionCreate(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt num_comp, CeedInt comp_stride, CeedSize l_size, 488 CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, CeedElemRestriction *rstr) { 489 if (!ceed->ElemRestrictionCreate) { 490 Ceed delegate; 491 492 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 493 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreate"); 494 CeedCall(CeedElemRestrictionCreate(delegate, num_elem, elem_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, rstr)); 495 return CEED_ERROR_SUCCESS; 496 } 497 498 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 499 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 500 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); 501 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction component stride must be at least 1"); 502 503 CeedCall(CeedCalloc(1, rstr)); 504 CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); 505 (*rstr)->ref_count = 1; 506 (*rstr)->num_elem = num_elem; 507 (*rstr)->elem_size = elem_size; 508 (*rstr)->num_comp = num_comp; 509 (*rstr)->comp_stride = comp_stride; 510 (*rstr)->l_size = l_size; 511 (*rstr)->num_block = num_elem; 512 (*rstr)->block_size = 1; 513 (*rstr)->rstr_type = CEED_RESTRICTION_STANDARD; 514 CeedCall(ceed->ElemRestrictionCreate(mem_type, copy_mode, offsets, NULL, NULL, *rstr)); 515 return CEED_ERROR_SUCCESS; 516 } 517 518 /** 519 @brief Create a CeedElemRestriction with orientation signs 520 521 @param[in] ceed Ceed object where the CeedElemRestriction will be created 522 @param[in] num_elem Number of elements described in the @a offsets array 523 @param[in] elem_size Size (number of "nodes") per element 524 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 525 @param[in] comp_stride Stride between components for the same L-vector "node". 526 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. 527 @param[in] l_size The size of the L-vector. 528 This vector may be larger than the elements and fields given by this restriction. 529 @param[in] mem_type Memory type of the @a offsets array, see CeedMemType 530 @param[in] copy_mode Copy mode for the @a offsets array, see CeedCopyMode 531 @param[in] offsets Array of shape [@a num_elem, @a elem_size]. 532 Row i holds the ordered list of the offsets (into the input CeedVector) for the unknowns corresponding to element i, where 533 0 <= i < @a num_elem. All offsets must be in the range [0, @a l_size - 1]. 534 @param[in] orients Array of shape [@a num_elem, @a elem_size] with bool false for positively oriented and true to flip the orientation. 535 @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored 536 537 @return An error code: 0 - success, otherwise - failure 538 539 @ref User 540 **/ 541 int CeedElemRestrictionCreateOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt num_comp, CeedInt comp_stride, CeedSize l_size, 542 CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, const bool *orients, 543 CeedElemRestriction *rstr) { 544 if (!ceed->ElemRestrictionCreate) { 545 Ceed delegate; 546 547 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 548 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreate"); 549 CeedCall( 550 CeedElemRestrictionCreateOriented(delegate, num_elem, elem_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, orients, rstr)); 551 return CEED_ERROR_SUCCESS; 552 } 553 554 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 555 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 556 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); 557 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction component stride must be at least 1"); 558 559 CeedCall(CeedCalloc(1, rstr)); 560 CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); 561 (*rstr)->ref_count = 1; 562 (*rstr)->num_elem = num_elem; 563 (*rstr)->elem_size = elem_size; 564 (*rstr)->num_comp = num_comp; 565 (*rstr)->comp_stride = comp_stride; 566 (*rstr)->l_size = l_size; 567 (*rstr)->num_block = num_elem; 568 (*rstr)->block_size = 1; 569 (*rstr)->rstr_type = CEED_RESTRICTION_ORIENTED; 570 CeedCall(ceed->ElemRestrictionCreate(mem_type, copy_mode, offsets, orients, NULL, *rstr)); 571 return CEED_ERROR_SUCCESS; 572 } 573 574 /** 575 @brief Create a CeedElemRestriction with a general tridiagonal transformation matrix for curl-conforming elements 576 577 @param[in] ceed Ceed object where the CeedElemRestriction will be created 578 @param[in] num_elem Number of elements described in the @a offsets array 579 @param[in] elem_size Size (number of "nodes") per element 580 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 581 @param[in] comp_stride Stride between components for the same L-vector "node". 582 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. 583 @param[in] l_size The size of the L-vector. 584 This vector may be larger than the elements and fields given by this restriction. 585 @param[in] mem_type Memory type of the @a offsets array, see CeedMemType 586 @param[in] copy_mode Copy mode for the @a offsets array, see CeedCopyMode 587 @param[in] offsets Array of shape [@a num_elem, @a elem_size]. 588 Row i holds the ordered list of the offsets (into the input CeedVector) for the unknowns corresponding to element i, 589 where 0 <= i < @a num_elem. All offsets must be in the range [0, @a l_size - 1]. 590 @param[in] curl_orients Array of shape [@a num_elem, @a 3 * elem_size] representing a row-major tridiagonal matrix (curl_orients[i * 3 * elem_size] 591 = curl_orients[(i + 1) * 3 * elem_size - 1] = 0, where 0 <= i < @a num_elem) which is applied to the element unknowns upon restriction. This 592 orientation matrix allows for pairs of face degrees of freedom on elements for H(curl) spaces to be coupled in the element restriction operation, 593 which is a way to resolve face orientation issues for 3D meshes (https://dl.acm.org/doi/pdf/10.1145/3524456). 594 @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored 595 596 @return An error code: 0 - success, otherwise - failure 597 598 @ref User 599 **/ 600 int CeedElemRestrictionCreateCurlOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt num_comp, CeedInt comp_stride, CeedSize l_size, 601 CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, const CeedInt8 *curl_orients, 602 CeedElemRestriction *rstr) { 603 if (!ceed->ElemRestrictionCreate) { 604 Ceed delegate; 605 606 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 607 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreate"); 608 CeedCall(CeedElemRestrictionCreateCurlOriented(delegate, num_elem, elem_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, 609 curl_orients, rstr)); 610 return CEED_ERROR_SUCCESS; 611 } 612 613 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 614 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 615 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); 616 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction component stride must be at least 1"); 617 618 CeedCall(CeedCalloc(1, rstr)); 619 CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); 620 (*rstr)->ref_count = 1; 621 (*rstr)->num_elem = num_elem; 622 (*rstr)->elem_size = elem_size; 623 (*rstr)->num_comp = num_comp; 624 (*rstr)->comp_stride = comp_stride; 625 (*rstr)->l_size = l_size; 626 (*rstr)->num_block = num_elem; 627 (*rstr)->block_size = 1; 628 (*rstr)->rstr_type = CEED_RESTRICTION_CURL_ORIENTED; 629 CeedCall(ceed->ElemRestrictionCreate(mem_type, copy_mode, offsets, NULL, curl_orients, *rstr)); 630 return CEED_ERROR_SUCCESS; 631 } 632 633 /** 634 @brief Create a strided CeedElemRestriction 635 636 @param[in] ceed Ceed object where the CeedElemRestriction will be created 637 @param[in] num_elem Number of elements described by the restriction 638 @param[in] elem_size Size (number of "nodes") per element 639 @param[in] num_comp Number of field components per interpolation "node" (1 for scalar fields) 640 @param[in] l_size The size of the L-vector. 641 This vector may be larger than the elements and fields given by this restriction. 642 @param[in] strides Array for strides between [nodes, components, elements]. 643 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]. 644 @a CEED_STRIDES_BACKEND may be used with vectors created by a Ceed backend. 645 @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored 646 647 @return An error code: 0 - success, otherwise - failure 648 649 @ref User 650 **/ 651 int CeedElemRestrictionCreateStrided(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt num_comp, CeedSize l_size, const CeedInt strides[3], 652 CeedElemRestriction *rstr) { 653 if (!ceed->ElemRestrictionCreate) { 654 Ceed delegate; 655 656 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 657 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreate"); 658 CeedCall(CeedElemRestrictionCreateStrided(delegate, num_elem, elem_size, num_comp, l_size, strides, rstr)); 659 return CEED_ERROR_SUCCESS; 660 } 661 662 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 663 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 664 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); 665 CeedCheck(l_size >= num_elem * elem_size * num_comp, ceed, CEED_ERROR_DIMENSION, "L-vector size must be at least num_elem * elem_size * num_comp"); 666 667 CeedCall(CeedCalloc(1, rstr)); 668 CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); 669 (*rstr)->ref_count = 1; 670 (*rstr)->num_elem = num_elem; 671 (*rstr)->elem_size = elem_size; 672 (*rstr)->num_comp = num_comp; 673 (*rstr)->l_size = l_size; 674 (*rstr)->num_block = num_elem; 675 (*rstr)->block_size = 1; 676 (*rstr)->rstr_type = CEED_RESTRICTION_STRIDED; 677 CeedCall(CeedMalloc(3, &(*rstr)->strides)); 678 for (CeedInt i = 0; i < 3; i++) (*rstr)->strides[i] = strides[i]; 679 CeedCall(ceed->ElemRestrictionCreate(CEED_MEM_HOST, CEED_OWN_POINTER, NULL, NULL, NULL, *rstr)); 680 return CEED_ERROR_SUCCESS; 681 } 682 683 /** 684 @brief Create a points CeedElemRestriction, for restricting for restricting from a all local points to the current element in which they are 685 located. 686 687 The offsets array is arranged as 688 689 element_0_start_index 690 element_1_start_index 691 ... 692 element_n_start_index 693 element_n_stop_index 694 element_0_point_0 695 element_0_point_1 696 ... 697 698 @param[in] ceed Ceed object where the CeedElemRestriction will be created 699 @param[in] num_elem Number of elements described in the @a offsets array 700 @param[in] num_points Number of points described in the @a offsets array 701 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields). 702 Components are assumed to be contiguous by point. 703 @param[in] l_size The size of the L-vector. 704 This vector may be larger than the elements and fields given by this restriction. 705 @param[in] mem_type Memory type of the @a offsets array, see CeedMemType 706 @param[in] copy_mode Copy mode for the @a offsets array, see CeedCopyMode 707 @param[in] offsets Array of size num_elem + 1 + num_points. 708 The first portion of the offsets array holds the ranges of indices corresponding to each element. 709 The second portion holds the indices for each element. 710 @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored 711 712 @return An error code: 0 - success, otherwise - failure 713 714 @ref Backend 715 **/ 716 int CeedElemRestrictionCreateAtPoints(Ceed ceed, CeedInt num_elem, CeedInt num_points, CeedInt num_comp, CeedSize l_size, CeedMemType mem_type, 717 CeedCopyMode copy_mode, const CeedInt *offsets, CeedElemRestriction *rstr) { 718 if (!ceed->ElemRestrictionCreateAtPoints) { 719 Ceed delegate; 720 721 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 722 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreateAtPoints"); 723 CeedCall(CeedElemRestrictionCreateAtPoints(delegate, num_elem, num_points, num_comp, l_size, mem_type, copy_mode, offsets, rstr)); 724 return CEED_ERROR_SUCCESS; 725 } 726 727 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 728 CeedCheck(num_points >= 0, ceed, CEED_ERROR_DIMENSION, "Number of points must be non-negative"); 729 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); 730 CeedCheck(l_size >= num_points * num_comp, ceed, CEED_ERROR_DIMENSION, "L-vector must be at least num_points * num_comp"); 731 732 CeedCall(CeedCalloc(1, rstr)); 733 CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); 734 (*rstr)->ref_count = 1; 735 (*rstr)->num_elem = num_elem; 736 (*rstr)->num_points = num_points; 737 (*rstr)->num_comp = num_comp; 738 (*rstr)->l_size = l_size; 739 (*rstr)->block_size = 1; 740 (*rstr)->rstr_type = CEED_RESTRICTION_POINTS; 741 CeedCall(ceed->ElemRestrictionCreateAtPoints(mem_type, copy_mode, offsets, NULL, NULL, *rstr)); 742 return CEED_ERROR_SUCCESS; 743 } 744 745 /** 746 @brief Create a blocked CeedElemRestriction, typically only called by backends 747 748 @param[in] ceed Ceed object where the CeedElemRestriction will be created 749 @param[in] num_elem Number of elements described in the @a offsets array 750 @param[in] elem_size Size (number of unknowns) per element 751 @param[in] block_size Number of elements in a block 752 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 753 @param[in] comp_stride Stride between components for the same L-vector "node". 754 Data for node i, component j, element k can be found in the L-vector at index offsets[i + k*elem_size] + j*comp_stride. 755 @param[in] l_size The size of the L-vector. 756 This vector may be larger than the elements and fields given by this restriction. 757 @param[in] mem_type Memory type of the @a offsets array, see CeedMemType 758 @param[in] copy_mode Copy mode for the @a offsets array, see CeedCopyMode 759 @param[in] offsets Array of shape [@a num_elem, @a elem_size]. 760 Row i holds the ordered list of the offsets (into the input CeedVector) for the unknowns corresponding to element i, 761 where 0 <= i < @a num_elem. All offsets must be in the range [0, @a l_size - 1]. The backend will permute and pad this array to the desired ordering 762 for the blocksize, which is typically given by the backend. The default reordering is to interlace elements. 763 @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored 764 765 @return An error code: 0 - success, otherwise - failure 766 767 @ref Backend 768 **/ 769 int CeedElemRestrictionCreateBlocked(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, CeedInt comp_stride, 770 CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, 771 CeedElemRestriction *rstr) { 772 CeedInt *block_offsets, num_block = (num_elem / block_size) + !!(num_elem % block_size); 773 774 if (!ceed->ElemRestrictionCreateBlocked) { 775 Ceed delegate; 776 777 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 778 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreateBlocked"); 779 CeedCall(CeedElemRestrictionCreateBlocked(delegate, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, 780 rstr)); 781 return CEED_ERROR_SUCCESS; 782 } 783 784 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 785 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 786 CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); 787 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); 788 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction component stride must be at least 1"); 789 790 CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_offsets)); 791 CeedCall(CeedPermutePadOffsets(offsets, block_offsets, num_block, num_elem, block_size, elem_size)); 792 793 CeedCall(CeedCalloc(1, rstr)); 794 CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); 795 (*rstr)->ref_count = 1; 796 (*rstr)->num_elem = num_elem; 797 (*rstr)->elem_size = elem_size; 798 (*rstr)->num_comp = num_comp; 799 (*rstr)->comp_stride = comp_stride; 800 (*rstr)->l_size = l_size; 801 (*rstr)->num_block = num_block; 802 (*rstr)->block_size = block_size; 803 (*rstr)->rstr_type = CEED_RESTRICTION_STANDARD; 804 CeedCall(ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, (const CeedInt *)block_offsets, NULL, NULL, *rstr)); 805 if (copy_mode == CEED_OWN_POINTER) CeedCall(CeedFree(&offsets)); 806 return CEED_ERROR_SUCCESS; 807 } 808 809 /** 810 @brief Create a blocked oriented CeedElemRestriction, typically only called by backends 811 812 @param[in] ceed Ceed object where the CeedElemRestriction will be created. 813 @param[in] num_elem Number of elements described in the @a offsets array. 814 @param[in] elem_size Size (number of unknowns) per element 815 @param[in] block_size Number of elements in a block 816 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 817 @param[in] comp_stride Stride between components for the same L-vector "node". 818 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. 819 @param[in] l_size The size of the L-vector. 820 This vector may be larger than the elements and fields given by this restriction. 821 @param[in] mem_type Memory type of the @a offsets array, see CeedMemType 822 @param[in] copy_mode Copy mode for the @a offsets array, see CeedCopyMode 823 @param[in] offsets Array of shape [@a num_elem, @a elem_size]. 824 Row i holds the ordered list of the offsets (into the input CeedVector) for the unknowns corresponding to element i, where 825 0 <= i < @a num_elem. All offsets must be in the range [0, @a l_size - 1]. The backend will permute and pad this array to the desired ordering for 826 the blocksize, which is typically given by the backend. The default reordering is to interlace elements. 827 @param[in] orients Array of shape [@a num_elem, @a elem_size] with bool false for positively oriented and true to flip the orientation. 828 Will also be permuted and padded similarly to @a offsets. 829 @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored 830 831 @return An error code: 0 - success, otherwise - failure 832 833 @ref Backend 834 **/ 835 int CeedElemRestrictionCreateBlockedOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, 836 CeedInt comp_stride, CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, 837 const CeedInt *offsets, const bool *orients, CeedElemRestriction *rstr) { 838 bool *block_orients; 839 CeedInt *block_offsets, num_block = (num_elem / block_size) + !!(num_elem % block_size); 840 841 if (!ceed->ElemRestrictionCreateBlocked) { 842 Ceed delegate; 843 844 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 845 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreateBlocked"); 846 CeedCall(CeedElemRestrictionCreateBlockedOriented(delegate, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, mem_type, copy_mode, 847 offsets, orients, rstr)); 848 return CEED_ERROR_SUCCESS; 849 } 850 851 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 852 CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); 853 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); 854 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction component stride must be at least 1"); 855 856 CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_offsets)); 857 CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_orients)); 858 CeedCall(CeedPermutePadOffsets(offsets, block_offsets, num_block, num_elem, block_size, elem_size)); 859 CeedCall(CeedPermutePadOrients(orients, block_orients, num_block, num_elem, block_size, elem_size)); 860 861 CeedCall(CeedCalloc(1, rstr)); 862 CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); 863 (*rstr)->ref_count = 1; 864 (*rstr)->num_elem = num_elem; 865 (*rstr)->elem_size = elem_size; 866 (*rstr)->num_comp = num_comp; 867 (*rstr)->comp_stride = comp_stride; 868 (*rstr)->l_size = l_size; 869 (*rstr)->num_block = num_block; 870 (*rstr)->block_size = block_size; 871 (*rstr)->rstr_type = CEED_RESTRICTION_ORIENTED; 872 CeedCall( 873 ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, (const CeedInt *)block_offsets, (const bool *)block_orients, NULL, *rstr)); 874 if (copy_mode == CEED_OWN_POINTER) CeedCall(CeedFree(&offsets)); 875 return CEED_ERROR_SUCCESS; 876 } 877 878 /** 879 @brief Create a blocked curl-oriented CeedElemRestriction, typically only called by backends 880 881 @param[in] ceed Ceed object where the CeedElemRestriction will be created. 882 @param[in] num_elem Number of elements described in the @a offsets array. 883 @param[in] elem_size Size (number of unknowns) per element 884 @param[in] block_size Number of elements in a block 885 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 886 @param[in] comp_stride Stride between components for the same L-vector "node". 887 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. 888 @param[in] l_size The size of the L-vector. 889 This vector may be larger than the elements and fields given by this restriction. 890 @param[in] mem_type Memory type of the @a offsets array, see CeedMemType 891 @param[in] copy_mode Copy mode for the @a offsets array, see CeedCopyMode 892 @param[in] offsets Array of shape [@a num_elem, @a elem_size]. 893 Row i holds the ordered list of the offsets (into the input CeedVector) for the unknowns corresponding to element i, 894 where 0 <= i < @a num_elem. All offsets must be in the range [0, @a l_size - 1]. The backend will permute and pad this array to the desired ordering 895 for the blocksize, which is typically given by the backend. The default reordering is to interlace elements. 896 @param[in] curl_orients Array of shape [@a num_elem, @a 3 * elem_size] representing a row-major tridiagonal matrix (curl_orients[i * 3 * elem_size] 897 = curl_orients[(i + 1) * 3 * elem_size - 1] = 0, where 0 <= i < @a num_elem) which is applied to the element unknowns upon restriction. This 898 orientation matrix allows for pairs of face degrees of freedom on elements for H(curl) spaces to be coupled in the element restriction operation, 899 which is a way to resolve face orientation issues for 3D meshes (https://dl.acm.org/doi/pdf/10.1145/3524456). Will also be permuted and padded 900 similarly to @a offsets. 901 @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored 902 903 @return An error code: 0 - success, otherwise - failure 904 905 @ref Backend 906 **/ 907 int CeedElemRestrictionCreateBlockedCurlOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, 908 CeedInt comp_stride, CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, 909 const CeedInt *offsets, const CeedInt8 *curl_orients, CeedElemRestriction *rstr) { 910 CeedInt8 *block_curl_orients; 911 CeedInt *block_offsets, num_block = (num_elem / block_size) + !!(num_elem % block_size); 912 913 if (!ceed->ElemRestrictionCreateBlocked) { 914 Ceed delegate; 915 916 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 917 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreateBlocked"); 918 CeedCall(CeedElemRestrictionCreateBlockedCurlOriented(delegate, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, mem_type, 919 copy_mode, offsets, curl_orients, rstr)); 920 return CEED_ERROR_SUCCESS; 921 } 922 923 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 924 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 925 CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); 926 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); 927 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction component stride must be at least 1"); 928 929 CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_offsets)); 930 CeedCall(CeedCalloc(num_block * block_size * 3 * elem_size, &block_curl_orients)); 931 CeedCall(CeedPermutePadOffsets(offsets, block_offsets, num_block, num_elem, block_size, elem_size)); 932 CeedCall(CeedPermutePadCurlOrients(curl_orients, block_curl_orients, num_block, num_elem, block_size, 3 * elem_size)); 933 934 CeedCall(CeedCalloc(1, rstr)); 935 CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); 936 (*rstr)->ref_count = 1; 937 (*rstr)->num_elem = num_elem; 938 (*rstr)->elem_size = elem_size; 939 (*rstr)->num_comp = num_comp; 940 (*rstr)->comp_stride = comp_stride; 941 (*rstr)->l_size = l_size; 942 (*rstr)->num_block = num_block; 943 (*rstr)->block_size = block_size; 944 (*rstr)->rstr_type = CEED_RESTRICTION_CURL_ORIENTED; 945 CeedCall(ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, (const CeedInt *)block_offsets, NULL, 946 (const CeedInt8 *)block_curl_orients, *rstr)); 947 if (copy_mode == CEED_OWN_POINTER) CeedCall(CeedFree(&offsets)); 948 return CEED_ERROR_SUCCESS; 949 } 950 951 /** 952 @brief Create a blocked strided CeedElemRestriction, typically only called by backends 953 954 @param[in] ceed Ceed object where the CeedElemRestriction will be created 955 @param[in] num_elem Number of elements described by the restriction 956 @param[in] elem_size Size (number of "nodes") per element 957 @param[in] block_size Number of elements in a block 958 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 959 @param[in] l_size The size of the L-vector. 960 This vector may be larger than the elements and fields given by this restriction. 961 @param[in] strides Array for strides between [nodes, components, elements]. 962 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]. 963 @a CEED_STRIDES_BACKEND may be used with vectors created by a Ceed backend. 964 @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored 965 966 @return An error code: 0 - success, otherwise - failure 967 968 @ref User 969 **/ 970 int CeedElemRestrictionCreateBlockedStrided(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, CeedSize l_size, 971 const CeedInt strides[3], CeedElemRestriction *rstr) { 972 CeedInt num_block = (num_elem / block_size) + !!(num_elem % block_size); 973 974 if (!ceed->ElemRestrictionCreateBlocked) { 975 Ceed delegate; 976 977 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 978 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreateBlocked"); 979 CeedCall(CeedElemRestrictionCreateBlockedStrided(delegate, num_elem, elem_size, block_size, num_comp, l_size, strides, rstr)); 980 return CEED_ERROR_SUCCESS; 981 } 982 983 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 984 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 985 CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); 986 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); 987 CeedCheck(l_size >= num_elem * elem_size * num_comp, ceed, CEED_ERROR_DIMENSION, "L-vector size must be at least num_elem * elem_size * num_comp"); 988 989 CeedCall(CeedCalloc(1, rstr)); 990 CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); 991 (*rstr)->ref_count = 1; 992 (*rstr)->num_elem = num_elem; 993 (*rstr)->elem_size = elem_size; 994 (*rstr)->num_comp = num_comp; 995 (*rstr)->l_size = l_size; 996 (*rstr)->num_block = num_block; 997 (*rstr)->block_size = block_size; 998 (*rstr)->rstr_type = CEED_RESTRICTION_STRIDED; 999 CeedCall(CeedMalloc(3, &(*rstr)->strides)); 1000 for (CeedInt i = 0; i < 3; i++) (*rstr)->strides[i] = strides[i]; 1001 CeedCall(ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, NULL, NULL, NULL, *rstr)); 1002 return CEED_ERROR_SUCCESS; 1003 } 1004 1005 /** 1006 @brief Copy the pointer to a CeedElemRestriction and set `CeedElemRestrictionApply()` implementation to use the unsigned version. 1007 1008 Both pointers should be destroyed with `CeedElemRestrictionDestroy()`. 1009 1010 @param[in] rstr CeedElemRestriction to create unsigned reference to 1011 @param[in,out] rstr_unsigned Variable to store unsigned CeedElemRestriction 1012 1013 @return An error code: 0 - success, otherwise - failure 1014 1015 @ref User 1016 **/ 1017 int CeedElemRestrictionCreateUnsignedCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_unsigned) { 1018 CeedCall(CeedCalloc(1, rstr_unsigned)); 1019 1020 // Copy old rstr 1021 memcpy(*rstr_unsigned, rstr, sizeof(struct CeedElemRestriction_private)); 1022 (*rstr_unsigned)->ceed = NULL; 1023 CeedCall(CeedReferenceCopy(rstr->ceed, &(*rstr_unsigned)->ceed)); 1024 (*rstr_unsigned)->ref_count = 1; 1025 (*rstr_unsigned)->strides = NULL; 1026 if (rstr->strides) { 1027 CeedCall(CeedMalloc(3, &(*rstr_unsigned)->strides)); 1028 for (CeedInt i = 0; i < 3; i++) (*rstr_unsigned)->strides[i] = rstr->strides[i]; 1029 } 1030 CeedCall(CeedElemRestrictionReferenceCopy(rstr, &(*rstr_unsigned)->rstr_base)); 1031 1032 // Override Apply 1033 (*rstr_unsigned)->Apply = rstr->ApplyUnsigned; 1034 return CEED_ERROR_SUCCESS; 1035 } 1036 1037 /** 1038 @brief Copy the pointer to a CeedElemRestriction and set `CeedElemRestrictionApply()` implementation to use the unoriented version. 1039 1040 Both pointers should be destroyed with `CeedElemRestrictionDestroy()`. 1041 1042 @param[in] rstr CeedElemRestriction to create unoriented reference to 1043 @param[in,out] rstr_unoriented Variable to store unoriented CeedElemRestriction 1044 1045 @return An error code: 0 - success, otherwise - failure 1046 1047 @ref User 1048 **/ 1049 int CeedElemRestrictionCreateUnorientedCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_unoriented) { 1050 CeedCall(CeedCalloc(1, rstr_unoriented)); 1051 1052 // Copy old rstr 1053 memcpy(*rstr_unoriented, rstr, sizeof(struct CeedElemRestriction_private)); 1054 (*rstr_unoriented)->ceed = NULL; 1055 CeedCall(CeedReferenceCopy(rstr->ceed, &(*rstr_unoriented)->ceed)); 1056 (*rstr_unoriented)->ref_count = 1; 1057 (*rstr_unoriented)->strides = NULL; 1058 if (rstr->strides) { 1059 CeedCall(CeedMalloc(3, &(*rstr_unoriented)->strides)); 1060 for (CeedInt i = 0; i < 3; i++) (*rstr_unoriented)->strides[i] = rstr->strides[i]; 1061 } 1062 CeedCall(CeedElemRestrictionReferenceCopy(rstr, &(*rstr_unoriented)->rstr_base)); 1063 1064 // Override Apply 1065 (*rstr_unoriented)->Apply = rstr->ApplyUnoriented; 1066 return CEED_ERROR_SUCCESS; 1067 } 1068 1069 /** 1070 @brief Copy the pointer to a CeedElemRestriction. 1071 1072 Both pointers should be destroyed with `CeedElemRestrictionDestroy()`. 1073 1074 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. 1075 This CeedElemRestriction will be destroyed if `rstr_copy` is the only reference to this CeedElemRestriction. 1076 1077 @param[in] rstr CeedElemRestriction to copy reference to 1078 @param[in,out] rstr_copy Variable to store copied reference 1079 1080 @return An error code: 0 - success, otherwise - failure 1081 1082 @ref User 1083 **/ 1084 int CeedElemRestrictionReferenceCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_copy) { 1085 if (rstr != CEED_ELEMRESTRICTION_NONE) CeedCall(CeedElemRestrictionReference(rstr)); 1086 CeedCall(CeedElemRestrictionDestroy(rstr_copy)); 1087 *rstr_copy = rstr; 1088 return CEED_ERROR_SUCCESS; 1089 } 1090 1091 /** 1092 @brief Create CeedVectors associated with a CeedElemRestriction 1093 1094 @param[in] rstr CeedElemRestriction 1095 @param[out] l_vec The address of the L-vector to be created, or NULL 1096 @param[out] e_vec The address of the E-vector to be created, or NULL 1097 1098 @return An error code: 0 - success, otherwise - failure 1099 1100 @ref User 1101 **/ 1102 int CeedElemRestrictionCreateVector(CeedElemRestriction rstr, CeedVector *l_vec, CeedVector *e_vec) { 1103 CeedSize e_size, l_size; 1104 l_size = rstr->l_size; 1105 e_size = rstr->num_block * rstr->block_size * rstr->elem_size * rstr->num_comp; 1106 if (l_vec) CeedCall(CeedVectorCreate(rstr->ceed, l_size, l_vec)); 1107 if (e_vec) CeedCall(CeedVectorCreate(rstr->ceed, e_size, e_vec)); 1108 return CEED_ERROR_SUCCESS; 1109 } 1110 1111 /** 1112 @brief Restrict an L-vector to an E-vector or apply its transpose 1113 1114 @param[in] rstr CeedElemRestriction 1115 @param[in] t_mode Apply restriction or transpose 1116 @param[in] u Input vector (of size @a l_size when t_mode=@ref CEED_NOTRANSPOSE) 1117 @param[out] ru Output vector (of shape [@a num_elem * @a elem_size] when t_mode=@ref CEED_NOTRANSPOSE). 1118 Ordering of the e-vector is decided by the backend. 1119 @param[in] request Request or @ref CEED_REQUEST_IMMEDIATE 1120 1121 @return An error code: 0 - success, otherwise - failure 1122 1123 @ref User 1124 **/ 1125 int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, CeedRequest *request) { 1126 CeedInt m, n; 1127 1128 if (t_mode == CEED_NOTRANSPOSE) { 1129 m = rstr->num_block * rstr->block_size * rstr->elem_size * rstr->num_comp; 1130 n = rstr->l_size; 1131 } else { 1132 m = rstr->l_size; 1133 n = rstr->num_block * rstr->block_size * rstr->elem_size * rstr->num_comp; 1134 } 1135 CeedCheck(n == u->length, rstr->ceed, CEED_ERROR_DIMENSION, 1136 "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", u->length, m, n); 1137 CeedCheck(m == ru->length, rstr->ceed, CEED_ERROR_DIMENSION, 1138 "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", ru->length, m, n); 1139 if (rstr->num_elem > 0) CeedCall(rstr->Apply(rstr, t_mode, u, ru, request)); 1140 return CEED_ERROR_SUCCESS; 1141 } 1142 1143 /** 1144 @brief Restrict an L-vector of points to a single element or apply its transpose 1145 1146 @param[in] rstr CeedElemRestriction 1147 @param[in] t_mode Apply restriction or transpose 1148 @param[in] u Input vector (of size @a l_size when t_mode=@ref CEED_NOTRANSPOSE) 1149 @param[out] ru Output vector (of shape [@a num_elem * @a elem_size] when t_mode=@ref CEED_NOTRANSPOSE). 1150 Ordering of the e-vector is decided by the backend. 1151 @param[in] request Request or @ref CEED_REQUEST_IMMEDIATE 1152 1153 @return An error code: 0 - success, otherwise - failure 1154 1155 @ref User 1156 **/ 1157 int CeedElemRestrictionApplyAtPoints(CeedElemRestriction rstr, CeedInt elem, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, 1158 CeedRequest *request) { 1159 CeedInt m, n; 1160 1161 if (t_mode == CEED_NOTRANSPOSE) { 1162 CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr, elem, &m)); 1163 n = rstr->l_size; 1164 } else { 1165 m = rstr->l_size; 1166 CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr, elem, &n)); 1167 } 1168 CeedCheck(n <= u->length, rstr->ceed, CEED_ERROR_DIMENSION, 1169 "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT 1170 ") for element %" CeedInt_FMT, 1171 u->length, m, n, elem); 1172 CeedCheck(m <= ru->length, rstr->ceed, CEED_ERROR_DIMENSION, 1173 "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT 1174 ") for element %" CeedInt_FMT, 1175 ru->length, m, n, elem); 1176 CeedCheck(elem <= rstr->num_elem, rstr->ceed, CEED_ERROR_DIMENSION, 1177 "Cannot retrieve element %" CeedInt_FMT ", element %" CeedInt_FMT " > total elements %" CeedInt_FMT "", elem, elem, rstr->num_elem); 1178 if (rstr->num_elem > 0) CeedCall(rstr->ApplyAtPoints(rstr, elem, t_mode, u, ru, request)); 1179 return CEED_ERROR_SUCCESS; 1180 } 1181 1182 /** 1183 @brief Restrict an L-vector to a block of an E-vector or apply its transpose 1184 1185 @param[in] rstr CeedElemRestriction 1186 @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 1187 [3*block_size : 4*block_size] 1188 @param[in] t_mode Apply restriction or transpose 1189 @param[in] u Input vector (of size @a l_size when t_mode=@ref CEED_NOTRANSPOSE) 1190 @param[out] ru Output vector (of shape [@a block_size * @a elem_size] when t_mode=@ref CEED_NOTRANSPOSE). 1191 Ordering of the e-vector is decided by the backend. 1192 @param[in] request Request or @ref CEED_REQUEST_IMMEDIATE 1193 1194 @return An error code: 0 - success, otherwise - failure 1195 1196 @ref Backend 1197 **/ 1198 int CeedElemRestrictionApplyBlock(CeedElemRestriction rstr, CeedInt block, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, 1199 CeedRequest *request) { 1200 CeedInt m, n; 1201 1202 CeedCheck(rstr->ApplyBlock, rstr->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionApplyBlock"); 1203 1204 if (t_mode == CEED_NOTRANSPOSE) { 1205 m = rstr->block_size * rstr->elem_size * rstr->num_comp; 1206 n = rstr->l_size; 1207 } else { 1208 m = rstr->l_size; 1209 n = rstr->block_size * rstr->elem_size * rstr->num_comp; 1210 } 1211 CeedCheck(n == u->length, rstr->ceed, CEED_ERROR_DIMENSION, 1212 "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", u->length, m, n); 1213 CeedCheck(m == ru->length, rstr->ceed, CEED_ERROR_DIMENSION, 1214 "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", ru->length, m, n); 1215 CeedCheck(rstr->block_size * block <= rstr->num_elem, rstr->ceed, CEED_ERROR_DIMENSION, 1216 "Cannot retrieve block %" CeedInt_FMT ", element %" CeedInt_FMT " > total elements %" CeedInt_FMT "", block, rstr->block_size * block, 1217 rstr->num_elem); 1218 CeedCall(rstr->ApplyBlock(rstr, block, t_mode, u, ru, request)); 1219 return CEED_ERROR_SUCCESS; 1220 } 1221 1222 /** 1223 @brief Get the Ceed associated with a CeedElemRestriction 1224 1225 @param[in] rstr CeedElemRestriction 1226 @param[out] ceed Variable to store Ceed 1227 1228 @return An error code: 0 - success, otherwise - failure 1229 1230 @ref Advanced 1231 **/ 1232 int CeedElemRestrictionGetCeed(CeedElemRestriction rstr, Ceed *ceed) { 1233 *ceed = rstr->ceed; 1234 return CEED_ERROR_SUCCESS; 1235 } 1236 1237 /** 1238 @brief Get the L-vector component stride 1239 1240 @param[in] rstr CeedElemRestriction 1241 @param[out] comp_stride Variable to store component stride 1242 1243 @return An error code: 0 - success, otherwise - failure 1244 1245 @ref Advanced 1246 **/ 1247 int CeedElemRestrictionGetCompStride(CeedElemRestriction rstr, CeedInt *comp_stride) { 1248 *comp_stride = rstr->comp_stride; 1249 return CEED_ERROR_SUCCESS; 1250 } 1251 1252 /** 1253 @brief Get the total number of elements in the range of a CeedElemRestriction 1254 1255 @param[in] rstr CeedElemRestriction 1256 @param[out] num_elem Variable to store number of elements 1257 1258 @return An error code: 0 - success, otherwise - failure 1259 1260 @ref Advanced 1261 **/ 1262 int CeedElemRestrictionGetNumElements(CeedElemRestriction rstr, CeedInt *num_elem) { 1263 *num_elem = rstr->num_elem; 1264 return CEED_ERROR_SUCCESS; 1265 } 1266 1267 /** 1268 @brief Get the size of elements in the CeedElemRestriction 1269 1270 @param[in] rstr CeedElemRestriction 1271 @param[out] elem_size Variable to store size of elements 1272 1273 @return An error code: 0 - success, otherwise - failure 1274 1275 @ref Advanced 1276 **/ 1277 int CeedElemRestrictionGetElementSize(CeedElemRestriction rstr, CeedInt *elem_size) { 1278 *elem_size = rstr->elem_size; 1279 return CEED_ERROR_SUCCESS; 1280 } 1281 1282 /** 1283 @brief Get the size of the l-vector for a CeedElemRestriction 1284 1285 @param[in] rstr CeedElemRestriction 1286 @param[out] l_size Variable to store number of nodes 1287 1288 @return An error code: 0 - success, otherwise - failure 1289 1290 @ref Advanced 1291 **/ 1292 int CeedElemRestrictionGetLVectorSize(CeedElemRestriction rstr, CeedSize *l_size) { 1293 *l_size = rstr->l_size; 1294 return CEED_ERROR_SUCCESS; 1295 } 1296 1297 /** 1298 @brief Get the number of components in the elements of a CeedElemRestriction 1299 1300 @param[in] rstr CeedElemRestriction 1301 @param[out] num_comp Variable to store number of components 1302 1303 @return An error code: 0 - success, otherwise - failure 1304 1305 @ref Advanced 1306 **/ 1307 int CeedElemRestrictionGetNumComponents(CeedElemRestriction rstr, CeedInt *num_comp) { 1308 *num_comp = rstr->num_comp; 1309 return CEED_ERROR_SUCCESS; 1310 } 1311 1312 /** 1313 @brief Get the number of blocks in a CeedElemRestriction 1314 1315 @param[in] rstr CeedElemRestriction 1316 @param[out] num_block Variable to store number of blocks 1317 1318 @return An error code: 0 - success, otherwise - failure 1319 1320 @ref Advanced 1321 **/ 1322 int CeedElemRestrictionGetNumBlocks(CeedElemRestriction rstr, CeedInt *num_block) { 1323 *num_block = rstr->num_block; 1324 return CEED_ERROR_SUCCESS; 1325 } 1326 1327 /** 1328 @brief Get the size of blocks in the CeedElemRestriction 1329 1330 @param[in] rstr CeedElemRestriction 1331 @param[out] block_size Variable to store size of blocks 1332 1333 @return An error code: 0 - success, otherwise - failure 1334 1335 @ref Advanced 1336 **/ 1337 int CeedElemRestrictionGetBlockSize(CeedElemRestriction rstr, CeedInt *block_size) { 1338 *block_size = rstr->block_size; 1339 return CEED_ERROR_SUCCESS; 1340 } 1341 1342 /** 1343 @brief Get the multiplicity of nodes in a CeedElemRestriction 1344 1345 @param[in] rstr CeedElemRestriction 1346 @param[out] mult Vector to store multiplicity (of size l_size) 1347 1348 @return An error code: 0 - success, otherwise - failure 1349 1350 @ref User 1351 **/ 1352 int CeedElemRestrictionGetMultiplicity(CeedElemRestriction rstr, CeedVector mult) { 1353 CeedVector e_vec; 1354 1355 // Create e_vec to hold intermediate computation in E^T (E 1) 1356 CeedCall(CeedElemRestrictionCreateVector(rstr, NULL, &e_vec)); 1357 1358 // Compute e_vec = E * 1 1359 CeedCall(CeedVectorSetValue(mult, 1.0)); 1360 CeedCall(CeedElemRestrictionApply(rstr, CEED_NOTRANSPOSE, mult, e_vec, CEED_REQUEST_IMMEDIATE)); 1361 // Compute multiplicity, mult = E^T * e_vec = E^T (E 1) 1362 CeedCall(CeedVectorSetValue(mult, 0.0)); 1363 CeedCall(CeedElemRestrictionApply(rstr, CEED_TRANSPOSE, e_vec, mult, CEED_REQUEST_IMMEDIATE)); 1364 // Cleanup 1365 CeedCall(CeedVectorDestroy(&e_vec)); 1366 return CEED_ERROR_SUCCESS; 1367 } 1368 1369 /** 1370 @brief View a CeedElemRestriction 1371 1372 @param[in] rstr CeedElemRestriction to view 1373 @param[in] stream Stream to write; typically stdout/stderr or a file 1374 1375 @return Error code: 0 - success, otherwise - failure 1376 1377 @ref User 1378 **/ 1379 int CeedElemRestrictionView(CeedElemRestriction rstr, FILE *stream) { 1380 char stridesstr[500]; 1381 1382 if (rstr->strides) { 1383 sprintf(stridesstr, "[%" CeedInt_FMT ", %" CeedInt_FMT ", %" CeedInt_FMT "]", rstr->strides[0], rstr->strides[1], rstr->strides[2]); 1384 } else { 1385 sprintf(stridesstr, "%" CeedInt_FMT, rstr->comp_stride); 1386 } 1387 1388 fprintf(stream, "%sCeedElemRestriction from (%td, %" CeedInt_FMT ") to %" CeedInt_FMT " elements with %" CeedInt_FMT " nodes each and %s %s\n", 1389 rstr->block_size > 1 ? "Blocked " : "", rstr->l_size, rstr->num_comp, rstr->num_elem, rstr->elem_size, 1390 rstr->strides ? "strides" : "component stride", stridesstr); 1391 return CEED_ERROR_SUCCESS; 1392 } 1393 1394 /** 1395 @brief Destroy a CeedElemRestriction 1396 1397 @param[in,out] rstr CeedElemRestriction to destroy 1398 1399 @return An error code: 0 - success, otherwise - failure 1400 1401 @ref User 1402 **/ 1403 int CeedElemRestrictionDestroy(CeedElemRestriction *rstr) { 1404 if (!*rstr || *rstr == CEED_ELEMRESTRICTION_NONE || --(*rstr)->ref_count > 0) { 1405 *rstr = NULL; 1406 return CEED_ERROR_SUCCESS; 1407 } 1408 CeedCheck((*rstr)->num_readers == 0, (*rstr)->ceed, CEED_ERROR_ACCESS, 1409 "Cannot destroy CeedElemRestriction, a process has read access to the offset data"); 1410 1411 // Only destroy backend data once between rstr and unsigned copy 1412 if ((*rstr)->rstr_base) CeedCall(CeedElemRestrictionDestroy(&(*rstr)->rstr_base)); 1413 else if ((*rstr)->Destroy) CeedCall((*rstr)->Destroy(*rstr)); 1414 1415 CeedCall(CeedFree(&(*rstr)->strides)); 1416 CeedCall(CeedDestroy(&(*rstr)->ceed)); 1417 CeedCall(CeedFree(rstr)); 1418 return CEED_ERROR_SUCCESS; 1419 } 1420 1421 /// @} 1422