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