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