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