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] elem Element number in range 0..@a num_elem 1132 @param[in] t_mode Apply restriction or transpose 1133 @param[in] u Input vector (of size @a l_size when t_mode=@ref CEED_NOTRANSPOSE) 1134 @param[out] ru Output vector (of shape [@a num_elem * @a elem_size] when t_mode=@ref CEED_NOTRANSPOSE). 1135 Ordering of the e-vector is decided by the backend. 1136 @param[in] request Request or @ref CEED_REQUEST_IMMEDIATE 1137 1138 @return An error code: 0 - success, otherwise - failure 1139 1140 @ref User 1141 **/ 1142 int CeedElemRestrictionApplyAtPointsInElement(CeedElemRestriction rstr, CeedInt elem, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, 1143 CeedRequest *request) { 1144 CeedInt m, n; 1145 1146 if (t_mode == CEED_NOTRANSPOSE) { 1147 CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr, elem, &m)); 1148 n = rstr->l_size; 1149 } else { 1150 m = rstr->l_size; 1151 CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr, elem, &n)); 1152 } 1153 CeedCheck(n <= u->length, rstr->ceed, CEED_ERROR_DIMENSION, 1154 "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT 1155 ") for element %" CeedInt_FMT, 1156 u->length, m, n, elem); 1157 CeedCheck(m <= ru->length, rstr->ceed, CEED_ERROR_DIMENSION, 1158 "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT 1159 ") for element %" CeedInt_FMT, 1160 ru->length, m, n, elem); 1161 CeedCheck(elem < rstr->num_elem, rstr->ceed, CEED_ERROR_DIMENSION, 1162 "Cannot retrieve element %" CeedInt_FMT ", element %" CeedInt_FMT " > total elements %" CeedInt_FMT "", elem, elem, rstr->num_elem); 1163 if (rstr->num_elem > 0) CeedCall(rstr->ApplyAtPointsInElement(rstr, elem, t_mode, u, ru, request)); 1164 return CEED_ERROR_SUCCESS; 1165 } 1166 1167 /** 1168 @brief Restrict an L-vector to a block of an E-vector or apply its transpose 1169 1170 @param[in] rstr CeedElemRestriction 1171 @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 1172 [3*block_size : 4*block_size] 1173 @param[in] t_mode Apply restriction or transpose 1174 @param[in] u Input vector (of size @a l_size when t_mode=@ref CEED_NOTRANSPOSE) 1175 @param[out] ru Output vector (of shape [@a block_size * @a elem_size] when t_mode=@ref CEED_NOTRANSPOSE). 1176 Ordering of the e-vector is decided by the backend. 1177 @param[in] request Request or @ref CEED_REQUEST_IMMEDIATE 1178 1179 @return An error code: 0 - success, otherwise - failure 1180 1181 @ref Backend 1182 **/ 1183 int CeedElemRestrictionApplyBlock(CeedElemRestriction rstr, CeedInt block, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, 1184 CeedRequest *request) { 1185 CeedInt m, n; 1186 1187 CeedCheck(rstr->ApplyBlock, rstr->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionApplyBlock"); 1188 1189 if (t_mode == CEED_NOTRANSPOSE) { 1190 m = rstr->block_size * rstr->elem_size * rstr->num_comp; 1191 n = rstr->l_size; 1192 } else { 1193 m = rstr->l_size; 1194 n = rstr->block_size * rstr->elem_size * rstr->num_comp; 1195 } 1196 CeedCheck(n == u->length, rstr->ceed, CEED_ERROR_DIMENSION, 1197 "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", u->length, m, n); 1198 CeedCheck(m == ru->length, rstr->ceed, CEED_ERROR_DIMENSION, 1199 "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", ru->length, m, n); 1200 CeedCheck(rstr->block_size * block <= rstr->num_elem, rstr->ceed, CEED_ERROR_DIMENSION, 1201 "Cannot retrieve block %" CeedInt_FMT ", element %" CeedInt_FMT " > total elements %" CeedInt_FMT "", block, rstr->block_size * block, 1202 rstr->num_elem); 1203 CeedCall(rstr->ApplyBlock(rstr, block, t_mode, u, ru, request)); 1204 return CEED_ERROR_SUCCESS; 1205 } 1206 1207 /** 1208 @brief Get the Ceed associated with a CeedElemRestriction 1209 1210 @param[in] rstr CeedElemRestriction 1211 @param[out] ceed Variable to store Ceed 1212 1213 @return An error code: 0 - success, otherwise - failure 1214 1215 @ref Advanced 1216 **/ 1217 int CeedElemRestrictionGetCeed(CeedElemRestriction rstr, Ceed *ceed) { 1218 *ceed = rstr->ceed; 1219 return CEED_ERROR_SUCCESS; 1220 } 1221 1222 /** 1223 @brief Get the L-vector component stride 1224 1225 @param[in] rstr CeedElemRestriction 1226 @param[out] comp_stride Variable to store component stride 1227 1228 @return An error code: 0 - success, otherwise - failure 1229 1230 @ref Advanced 1231 **/ 1232 int CeedElemRestrictionGetCompStride(CeedElemRestriction rstr, CeedInt *comp_stride) { 1233 *comp_stride = rstr->comp_stride; 1234 return CEED_ERROR_SUCCESS; 1235 } 1236 1237 /** 1238 @brief Get the total number of elements in the range of a CeedElemRestriction 1239 1240 @param[in] rstr CeedElemRestriction 1241 @param[out] num_elem Variable to store number of elements 1242 1243 @return An error code: 0 - success, otherwise - failure 1244 1245 @ref Advanced 1246 **/ 1247 int CeedElemRestrictionGetNumElements(CeedElemRestriction rstr, CeedInt *num_elem) { 1248 *num_elem = rstr->num_elem; 1249 return CEED_ERROR_SUCCESS; 1250 } 1251 1252 /** 1253 @brief Get the size of elements in the CeedElemRestriction 1254 1255 @param[in] rstr CeedElemRestriction 1256 @param[out] elem_size Variable to store size of elements 1257 1258 @return An error code: 0 - success, otherwise - failure 1259 1260 @ref Advanced 1261 **/ 1262 int CeedElemRestrictionGetElementSize(CeedElemRestriction rstr, CeedInt *elem_size) { 1263 *elem_size = rstr->elem_size; 1264 return CEED_ERROR_SUCCESS; 1265 } 1266 1267 /** 1268 1269 @brief Get the number of points in the l-vector for a points CeedElemRestriction 1270 1271 @param[in] rstr CeedElemRestriction 1272 @param[out] num_points The number of points in the l-vector 1273 1274 @return An error code: 0 - success, otherwise - failure 1275 1276 @ref User 1277 **/ 1278 int CeedElemRestrictionGetNumPoints(CeedElemRestriction rstr, CeedInt *num_points) { 1279 Ceed ceed; 1280 1281 CeedCall(CeedElemRestrictionGetCeed(rstr, &ceed)); 1282 CeedCheck(rstr->rstr_type == CEED_RESTRICTION_POINTS, ceed, CEED_ERROR_INCOMPATIBLE, 1283 "Can only retrieve the number of points for a points CeedElemRestriction"); 1284 1285 *num_points = rstr->num_points; 1286 return CEED_ERROR_SUCCESS; 1287 } 1288 1289 /** 1290 1291 @brief Get the number of points in an element of a points CeedElemRestriction 1292 1293 @param[in] rstr CeedElemRestriction 1294 @param[in] elem Index number of element to retrieve the number of points for 1295 @param[out] num_points The number of points in the element at index elem 1296 1297 @return An error code: 0 - success, otherwise - failure 1298 1299 @ref User 1300 **/ 1301 int CeedElemRestrictionGetNumPointsInElement(CeedElemRestriction rstr, CeedInt elem, CeedInt *num_points) { 1302 Ceed ceed; 1303 const CeedInt *offsets; 1304 1305 CeedCall(CeedElemRestrictionGetCeed(rstr, &ceed)); 1306 CeedCheck(rstr->rstr_type == CEED_RESTRICTION_POINTS, ceed, CEED_ERROR_INCOMPATIBLE, 1307 "Can only retrieve the number of points for a points CeedElemRestriction"); 1308 1309 CeedCall(CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets)); 1310 *num_points = offsets[elem + 1] - offsets[elem]; 1311 CeedCall(CeedElemRestrictionRestoreOffsets(rstr, &offsets)); 1312 return CEED_ERROR_SUCCESS; 1313 } 1314 1315 /** 1316 @brief Get the maximum number of points in an element for a CeedElemRestriction at points 1317 1318 @param[in] rstr CeedElemRestriction 1319 @param[out] max_points Variable to store size of elements 1320 1321 @return An error code: 0 - success, otherwise - failure 1322 1323 @ref Advanced 1324 **/ 1325 int CeedElemRestrictionGetMaxPointsInElement(CeedElemRestriction rstr, CeedInt *max_points) { 1326 Ceed ceed; 1327 CeedInt num_elem; 1328 CeedRestrictionType rstr_type; 1329 1330 CeedCall(CeedElemRestrictionGetCeed(rstr, &ceed)); 1331 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 1332 CeedCheck(rstr_type == CEED_RESTRICTION_POINTS, ceed, CEED_ERROR_INCOMPATIBLE, 1333 "Cannot compute max points for a CeedElemRestriction that does not use points"); 1334 1335 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 1336 *max_points = 0; 1337 for (CeedInt e = 0; e < num_elem; e++) { 1338 CeedInt num_points; 1339 1340 CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr, e, &num_points)); 1341 *max_points = CeedIntMax(num_points, *max_points); 1342 } 1343 return CEED_ERROR_SUCCESS; 1344 } 1345 1346 /** 1347 @brief Get the size of the l-vector for a CeedElemRestriction 1348 1349 @param[in] rstr CeedElemRestriction 1350 @param[out] l_size Variable to store number of nodes 1351 1352 @return An error code: 0 - success, otherwise - failure 1353 1354 @ref Advanced 1355 **/ 1356 int CeedElemRestrictionGetLVectorSize(CeedElemRestriction rstr, CeedSize *l_size) { 1357 *l_size = rstr->l_size; 1358 return CEED_ERROR_SUCCESS; 1359 } 1360 1361 /** 1362 @brief Get the number of components in the elements of a CeedElemRestriction 1363 1364 @param[in] rstr CeedElemRestriction 1365 @param[out] num_comp Variable to store number of components 1366 1367 @return An error code: 0 - success, otherwise - failure 1368 1369 @ref Advanced 1370 **/ 1371 int CeedElemRestrictionGetNumComponents(CeedElemRestriction rstr, CeedInt *num_comp) { 1372 *num_comp = rstr->num_comp; 1373 return CEED_ERROR_SUCCESS; 1374 } 1375 1376 /** 1377 @brief Get the number of blocks in a CeedElemRestriction 1378 1379 @param[in] rstr CeedElemRestriction 1380 @param[out] num_block Variable to store number of blocks 1381 1382 @return An error code: 0 - success, otherwise - failure 1383 1384 @ref Advanced 1385 **/ 1386 int CeedElemRestrictionGetNumBlocks(CeedElemRestriction rstr, CeedInt *num_block) { 1387 *num_block = rstr->num_block; 1388 return CEED_ERROR_SUCCESS; 1389 } 1390 1391 /** 1392 @brief Get the size of blocks in the CeedElemRestriction 1393 1394 @param[in] rstr CeedElemRestriction 1395 @param[out] block_size Variable to store size of blocks 1396 1397 @return An error code: 0 - success, otherwise - failure 1398 1399 @ref Advanced 1400 **/ 1401 int CeedElemRestrictionGetBlockSize(CeedElemRestriction rstr, CeedInt *block_size) { 1402 *block_size = rstr->block_size; 1403 return CEED_ERROR_SUCCESS; 1404 } 1405 1406 /** 1407 @brief Get the multiplicity of nodes in a CeedElemRestriction 1408 1409 @param[in] rstr CeedElemRestriction 1410 @param[out] mult Vector to store multiplicity (of size l_size) 1411 1412 @return An error code: 0 - success, otherwise - failure 1413 1414 @ref User 1415 **/ 1416 int CeedElemRestrictionGetMultiplicity(CeedElemRestriction rstr, CeedVector mult) { 1417 CeedVector e_vec; 1418 1419 // Create e_vec to hold intermediate computation in E^T (E 1) 1420 CeedCall(CeedElemRestrictionCreateVector(rstr, NULL, &e_vec)); 1421 1422 // Compute e_vec = E * 1 1423 CeedCall(CeedVectorSetValue(mult, 1.0)); 1424 CeedCall(CeedElemRestrictionApply(rstr, CEED_NOTRANSPOSE, mult, e_vec, CEED_REQUEST_IMMEDIATE)); 1425 // Compute multiplicity, mult = E^T * e_vec = E^T (E 1) 1426 CeedCall(CeedVectorSetValue(mult, 0.0)); 1427 CeedCall(CeedElemRestrictionApply(rstr, CEED_TRANSPOSE, e_vec, mult, CEED_REQUEST_IMMEDIATE)); 1428 // Cleanup 1429 CeedCall(CeedVectorDestroy(&e_vec)); 1430 return CEED_ERROR_SUCCESS; 1431 } 1432 1433 /** 1434 @brief View a CeedElemRestriction 1435 1436 @param[in] rstr CeedElemRestriction to view 1437 @param[in] stream Stream to write; typically stdout/stderr or a file 1438 1439 @return Error code: 0 - success, otherwise - failure 1440 1441 @ref User 1442 **/ 1443 int CeedElemRestrictionView(CeedElemRestriction rstr, FILE *stream) { 1444 CeedRestrictionType rstr_type; 1445 1446 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 1447 1448 if (rstr_type == CEED_RESTRICTION_POINTS) { 1449 CeedInt max_points; 1450 1451 CeedCall(CeedElemRestrictionGetMaxPointsInElement(rstr, &max_points)); 1452 fprintf(stream, 1453 "CeedElemRestriction at points from (%td, %" CeedInt_FMT ") to %" CeedInt_FMT " elements with a maximum of %" CeedInt_FMT 1454 " points on an element\n", 1455 rstr->l_size, rstr->num_comp, rstr->num_elem, max_points); 1456 } else { 1457 char stridesstr[500]; 1458 1459 if (rstr->strides) { 1460 sprintf(stridesstr, "[%" CeedInt_FMT ", %" CeedInt_FMT ", %" CeedInt_FMT "]", rstr->strides[0], rstr->strides[1], rstr->strides[2]); 1461 } else { 1462 sprintf(stridesstr, "%" CeedInt_FMT, rstr->comp_stride); 1463 } 1464 fprintf(stream, "%sCeedElemRestriction from (%td, %" CeedInt_FMT ") to %" CeedInt_FMT " elements with %" CeedInt_FMT " nodes each and %s %s\n", 1465 rstr->block_size > 1 ? "Blocked " : "", rstr->l_size, rstr->num_comp, rstr->num_elem, rstr->elem_size, 1466 rstr->strides ? "strides" : "component stride", stridesstr); 1467 } 1468 return CEED_ERROR_SUCCESS; 1469 } 1470 1471 /** 1472 @brief Destroy a CeedElemRestriction 1473 1474 @param[in,out] rstr CeedElemRestriction to destroy 1475 1476 @return An error code: 0 - success, otherwise - failure 1477 1478 @ref User 1479 **/ 1480 int CeedElemRestrictionDestroy(CeedElemRestriction *rstr) { 1481 if (!*rstr || *rstr == CEED_ELEMRESTRICTION_NONE || --(*rstr)->ref_count > 0) { 1482 *rstr = NULL; 1483 return CEED_ERROR_SUCCESS; 1484 } 1485 CeedCheck((*rstr)->num_readers == 0, (*rstr)->ceed, CEED_ERROR_ACCESS, 1486 "Cannot destroy CeedElemRestriction, a process has read access to the offset data"); 1487 1488 // Only destroy backend data once between rstr and unsigned copy 1489 if ((*rstr)->rstr_base) CeedCall(CeedElemRestrictionDestroy(&(*rstr)->rstr_base)); 1490 else if ((*rstr)->Destroy) CeedCall((*rstr)->Destroy(*rstr)); 1491 1492 CeedCall(CeedFree(&(*rstr)->strides)); 1493 CeedCall(CeedDestroy(&(*rstr)->ceed)); 1494 CeedCall(CeedFree(rstr)); 1495 return CEED_ERROR_SUCCESS; 1496 } 1497 1498 /// @} 1499