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