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