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 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 373 if (t_mode == CEED_TRANSPOSE) { 374 switch (rstr_type) { 375 case CEED_RESTRICTION_STRIDED: 376 case CEED_RESTRICTION_STANDARD: 377 scale = 1; 378 break; 379 case CEED_RESTRICTION_ORIENTED: 380 scale = 2; 381 break; 382 case CEED_RESTRICTION_CURL_ORIENTED: 383 scale = 6; 384 break; 385 } 386 } else { 387 switch (rstr_type) { 388 case CEED_RESTRICTION_STRIDED: 389 case CEED_RESTRICTION_STANDARD: 390 scale = 0; 391 break; 392 case CEED_RESTRICTION_ORIENTED: 393 scale = 1; 394 break; 395 case CEED_RESTRICTION_CURL_ORIENTED: 396 scale = 5; 397 break; 398 } 399 } 400 *flops = e_size * scale; 401 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 /// Indicate that no CeedElemRestriction is provided by the user 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; 669 CeedInt num_block = (num_elem / block_size) + !!(num_elem % block_size); 670 671 if (!ceed->ElemRestrictionCreateBlocked) { 672 Ceed delegate; 673 674 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 675 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreateBlocked"); 676 CeedCall(CeedElemRestrictionCreateBlocked(delegate, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, 677 rstr)); 678 return CEED_ERROR_SUCCESS; 679 } 680 681 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 682 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 683 CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); 684 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); 685 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction component stride must be at least 1"); 686 687 CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_offsets)); 688 CeedCall(CeedPermutePadOffsets(offsets, block_offsets, num_block, num_elem, block_size, elem_size)); 689 690 CeedCall(CeedCalloc(1, rstr)); 691 CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); 692 (*rstr)->ref_count = 1; 693 (*rstr)->num_elem = num_elem; 694 (*rstr)->elem_size = elem_size; 695 (*rstr)->num_comp = num_comp; 696 (*rstr)->comp_stride = comp_stride; 697 (*rstr)->l_size = l_size; 698 (*rstr)->num_block = num_block; 699 (*rstr)->block_size = block_size; 700 (*rstr)->rstr_type = CEED_RESTRICTION_STANDARD; 701 CeedCall(ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, (const CeedInt *)block_offsets, NULL, NULL, *rstr)); 702 if (copy_mode == CEED_OWN_POINTER) { 703 CeedCall(CeedFree(&offsets)); 704 } 705 return CEED_ERROR_SUCCESS; 706 } 707 708 /** 709 @brief Create a blocked oriented CeedElemRestriction, typically only called by backends 710 711 @param[in] ceed Ceed object where the CeedElemRestriction will be created. 712 @param[in] num_elem Number of elements described in the @a offsets array. 713 @param[in] elem_size Size (number of unknowns) per element 714 @param[in] block_size Number of elements in a block 715 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 716 @param[in] comp_stride Stride between components for the same L-vector "node". 717 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. 718 @param[in] l_size The size of the L-vector. 719 This vector may be larger than the elements and fields given by this restriction. 720 @param[in] mem_type Memory type of the @a offsets array, see CeedMemType 721 @param[in] copy_mode Copy mode for the @a offsets array, see CeedCopyMode 722 @param[in] offsets Array of shape [@a num_elem, @a elem_size]. 723 Row i holds the ordered list of the offsets (into the input CeedVector) for the unknowns corresponding to element i, where 724 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 725 the blocksize, which is typically given by the backend. The default reordering is to interlace elements. 726 @param[in] orients Array of shape [@a num_elem, @a elem_size] with bool false for positively oriented and true to flip the orientation. 727 Will also be permuted and padded similarly to @a offsets. 728 @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored 729 730 @return An error code: 0 - success, otherwise - failure 731 732 @ref Backend 733 **/ 734 int CeedElemRestrictionCreateBlockedOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, 735 CeedInt comp_stride, CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, 736 const CeedInt *offsets, const bool *orients, CeedElemRestriction *rstr) { 737 CeedInt *block_offsets; 738 bool *block_orients; 739 CeedInt num_block = (num_elem / block_size) + !!(num_elem % block_size); 740 741 if (!ceed->ElemRestrictionCreateBlocked) { 742 Ceed delegate; 743 744 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 745 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreateBlocked"); 746 CeedCall(CeedElemRestrictionCreateBlockedOriented(delegate, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, mem_type, copy_mode, 747 offsets, orients, rstr)); 748 return CEED_ERROR_SUCCESS; 749 } 750 751 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 752 CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); 753 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); 754 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction component stride must be at least 1"); 755 756 CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_offsets)); 757 CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_orients)); 758 CeedCall(CeedPermutePadOffsets(offsets, block_offsets, num_block, num_elem, block_size, elem_size)); 759 CeedCall(CeedPermutePadOrients(orients, block_orients, num_block, num_elem, block_size, elem_size)); 760 761 CeedCall(CeedCalloc(1, rstr)); 762 CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); 763 (*rstr)->ref_count = 1; 764 (*rstr)->num_elem = num_elem; 765 (*rstr)->elem_size = elem_size; 766 (*rstr)->num_comp = num_comp; 767 (*rstr)->comp_stride = comp_stride; 768 (*rstr)->l_size = l_size; 769 (*rstr)->num_block = num_block; 770 (*rstr)->block_size = block_size; 771 (*rstr)->rstr_type = CEED_RESTRICTION_ORIENTED; 772 CeedCall( 773 ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, (const CeedInt *)block_offsets, (const bool *)block_orients, NULL, *rstr)); 774 if (copy_mode == CEED_OWN_POINTER) { 775 CeedCall(CeedFree(&offsets)); 776 } 777 return CEED_ERROR_SUCCESS; 778 } 779 780 /** 781 @brief Create a blocked curl-oriented CeedElemRestriction, typically only called by backends 782 783 @param[in] ceed Ceed object where the CeedElemRestriction will be created. 784 @param[in] num_elem Number of elements described in the @a offsets array. 785 @param[in] elem_size Size (number of unknowns) per element 786 @param[in] block_size Number of elements in a block 787 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 788 @param[in] comp_stride Stride between components for the same L-vector "node". 789 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. 790 @param[in] l_size The size of the L-vector. 791 This vector may be larger than the elements and fields given by this restriction. 792 @param[in] mem_type Memory type of the @a offsets array, see CeedMemType 793 @param[in] copy_mode Copy mode for the @a offsets array, see CeedCopyMode 794 @param[in] offsets Array of shape [@a num_elem, @a elem_size]. 795 Row i holds the ordered list of the offsets (into the input CeedVector) for the unknowns corresponding to element i, 796 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 797 for the blocksize, which is typically given by the backend. The default reordering is to interlace elements. 798 @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] 799 = 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 800 orientation matrix allows for pairs of face degrees of freedom on elements for H(curl) spaces to be coupled in the element restriction operation, 801 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 802 similarly to @a offsets. 803 @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored 804 805 @return An error code: 0 - success, otherwise - failure 806 807 @ref Backend 808 **/ 809 int CeedElemRestrictionCreateBlockedCurlOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, 810 CeedInt comp_stride, CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, 811 const CeedInt *offsets, const CeedInt8 *curl_orients, CeedElemRestriction *rstr) { 812 CeedInt *block_offsets; 813 CeedInt8 *block_curl_orients; 814 CeedInt num_block = (num_elem / block_size) + !!(num_elem % block_size); 815 816 if (!ceed->ElemRestrictionCreateBlocked) { 817 Ceed delegate; 818 819 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 820 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreateBlocked"); 821 CeedCall(CeedElemRestrictionCreateBlockedCurlOriented(delegate, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, mem_type, 822 copy_mode, offsets, curl_orients, rstr)); 823 return CEED_ERROR_SUCCESS; 824 } 825 826 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 827 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 828 CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); 829 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); 830 CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction component stride must be at least 1"); 831 832 CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_offsets)); 833 CeedCall(CeedCalloc(num_block * block_size * 3 * elem_size, &block_curl_orients)); 834 CeedCall(CeedPermutePadOffsets(offsets, block_offsets, num_block, num_elem, block_size, elem_size)); 835 CeedCall(CeedPermutePadCurlOrients(curl_orients, block_curl_orients, num_block, num_elem, block_size, 3 * elem_size)); 836 837 CeedCall(CeedCalloc(1, rstr)); 838 CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); 839 (*rstr)->ref_count = 1; 840 (*rstr)->num_elem = num_elem; 841 (*rstr)->elem_size = elem_size; 842 (*rstr)->num_comp = num_comp; 843 (*rstr)->comp_stride = comp_stride; 844 (*rstr)->l_size = l_size; 845 (*rstr)->num_block = num_block; 846 (*rstr)->block_size = block_size; 847 (*rstr)->rstr_type = CEED_RESTRICTION_CURL_ORIENTED; 848 CeedCall(ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, (const CeedInt *)block_offsets, NULL, 849 (const CeedInt8 *)block_curl_orients, *rstr)); 850 if (copy_mode == CEED_OWN_POINTER) { 851 CeedCall(CeedFree(&offsets)); 852 } 853 return CEED_ERROR_SUCCESS; 854 } 855 856 /** 857 @brief Create a blocked strided CeedElemRestriction, typically only called by backends 858 859 @param[in] ceed Ceed object where the CeedElemRestriction will be created 860 @param[in] num_elem Number of elements described by the restriction 861 @param[in] elem_size Size (number of "nodes") per element 862 @param[in] block_size Number of elements in a block 863 @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) 864 @param[in] l_size The size of the L-vector. 865 This vector may be larger than the elements and fields given by this restriction. 866 @param[in] strides Array for strides between [nodes, components, elements]. 867 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]. 868 @a CEED_STRIDES_BACKEND may be used with vectors created by a Ceed backend. 869 @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored 870 871 @return An error code: 0 - success, otherwise - failure 872 873 @ref User 874 **/ 875 int CeedElemRestrictionCreateBlockedStrided(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, CeedSize l_size, 876 const CeedInt strides[3], CeedElemRestriction *rstr) { 877 CeedInt num_block = (num_elem / block_size) + !!(num_elem % block_size); 878 879 if (!ceed->ElemRestrictionCreateBlocked) { 880 Ceed delegate; 881 882 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); 883 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreateBlocked"); 884 CeedCall(CeedElemRestrictionCreateBlockedStrided(delegate, num_elem, elem_size, block_size, num_comp, l_size, strides, rstr)); 885 return CEED_ERROR_SUCCESS; 886 } 887 888 CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); 889 CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); 890 CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); 891 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); 892 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"); 893 894 CeedCall(CeedCalloc(1, rstr)); 895 CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); 896 (*rstr)->ref_count = 1; 897 (*rstr)->num_elem = num_elem; 898 (*rstr)->elem_size = elem_size; 899 (*rstr)->num_comp = num_comp; 900 (*rstr)->l_size = l_size; 901 (*rstr)->num_block = num_block; 902 (*rstr)->block_size = block_size; 903 (*rstr)->rstr_type = CEED_RESTRICTION_STRIDED; 904 CeedCall(CeedMalloc(3, &(*rstr)->strides)); 905 for (CeedInt i = 0; i < 3; i++) (*rstr)->strides[i] = strides[i]; 906 CeedCall(ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, NULL, NULL, NULL, *rstr)); 907 return CEED_ERROR_SUCCESS; 908 } 909 910 /** 911 @brief Copy the pointer to a CeedElemRestriction and set `CeedElemRestrictionApply()` implementation to use the unsigned version. 912 913 Both pointers should be destroyed with `CeedElemRestrictionDestroy()`. 914 915 @param[in] rstr CeedElemRestriction to create unsigned reference to 916 @param[in,out] rstr_unsigned Variable to store unsigned CeedElemRestriction 917 918 @return An error code: 0 - success, otherwise - failure 919 920 @ref User 921 **/ 922 int CeedElemRestrictionCreateUnsignedCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_unsigned) { 923 CeedCall(CeedCalloc(1, rstr_unsigned)); 924 925 // Copy old rstr 926 memcpy(*rstr_unsigned, rstr, sizeof(struct CeedElemRestriction_private)); 927 (*rstr_unsigned)->ceed = NULL; 928 CeedCall(CeedReferenceCopy(rstr->ceed, &(*rstr_unsigned)->ceed)); 929 (*rstr_unsigned)->ref_count = 1; 930 (*rstr_unsigned)->strides = NULL; 931 if (rstr->strides) { 932 CeedCall(CeedMalloc(3, &(*rstr_unsigned)->strides)); 933 for (CeedInt i = 0; i < 3; i++) (*rstr_unsigned)->strides[i] = rstr->strides[i]; 934 } 935 CeedCall(CeedElemRestrictionReferenceCopy(rstr, &(*rstr_unsigned)->rstr_base)); 936 937 // Override Apply 938 (*rstr_unsigned)->Apply = rstr->ApplyUnsigned; 939 940 return CEED_ERROR_SUCCESS; 941 } 942 943 /** 944 @brief Copy the pointer to a CeedElemRestriction and set `CeedElemRestrictionApply()` implementation to use the unoriented version. 945 946 Both pointers should be destroyed with `CeedElemRestrictionDestroy()`. 947 948 @param[in] rstr CeedElemRestriction to create unoriented reference to 949 @param[in,out] rstr_unoriented Variable to store unoriented CeedElemRestriction 950 951 @return An error code: 0 - success, otherwise - failure 952 953 @ref User 954 **/ 955 int CeedElemRestrictionCreateUnorientedCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_unoriented) { 956 CeedCall(CeedCalloc(1, rstr_unoriented)); 957 958 // Copy old rstr 959 memcpy(*rstr_unoriented, rstr, sizeof(struct CeedElemRestriction_private)); 960 (*rstr_unoriented)->ceed = NULL; 961 CeedCall(CeedReferenceCopy(rstr->ceed, &(*rstr_unoriented)->ceed)); 962 (*rstr_unoriented)->ref_count = 1; 963 (*rstr_unoriented)->strides = NULL; 964 if (rstr->strides) { 965 CeedCall(CeedMalloc(3, &(*rstr_unoriented)->strides)); 966 for (CeedInt i = 0; i < 3; i++) (*rstr_unoriented)->strides[i] = rstr->strides[i]; 967 } 968 CeedCall(CeedElemRestrictionReferenceCopy(rstr, &(*rstr_unoriented)->rstr_base)); 969 970 // Override Apply 971 (*rstr_unoriented)->Apply = rstr->ApplyUnoriented; 972 973 return CEED_ERROR_SUCCESS; 974 } 975 976 /** 977 @brief Copy the pointer to a CeedElemRestriction. 978 979 Both pointers should be destroyed with `CeedElemRestrictionDestroy()`. 980 981 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. 982 This CeedElemRestriction will be destroyed if `rstr_copy` is the only reference to this CeedElemRestriction. 983 984 @param[in] rstr CeedElemRestriction to copy reference to 985 @param[in,out] rstr_copy Variable to store copied reference 986 987 @return An error code: 0 - success, otherwise - failure 988 989 @ref User 990 **/ 991 int CeedElemRestrictionReferenceCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_copy) { 992 if (rstr != CEED_ELEMRESTRICTION_NONE) CeedCall(CeedElemRestrictionReference(rstr)); 993 CeedCall(CeedElemRestrictionDestroy(rstr_copy)); 994 *rstr_copy = rstr; 995 return CEED_ERROR_SUCCESS; 996 } 997 998 /** 999 @brief Create CeedVectors associated with a CeedElemRestriction 1000 1001 @param[in] rstr CeedElemRestriction 1002 @param[out] l_vec The address of the L-vector to be created, or NULL 1003 @param[out] e_vec The address of the E-vector to be created, or NULL 1004 1005 @return An error code: 0 - success, otherwise - failure 1006 1007 @ref User 1008 **/ 1009 int CeedElemRestrictionCreateVector(CeedElemRestriction rstr, CeedVector *l_vec, CeedVector *e_vec) { 1010 CeedSize e_size, l_size; 1011 l_size = rstr->l_size; 1012 e_size = rstr->num_block * rstr->block_size * rstr->elem_size * rstr->num_comp; 1013 if (l_vec) CeedCall(CeedVectorCreate(rstr->ceed, l_size, l_vec)); 1014 if (e_vec) CeedCall(CeedVectorCreate(rstr->ceed, e_size, e_vec)); 1015 return CEED_ERROR_SUCCESS; 1016 } 1017 1018 /** 1019 @brief Restrict an L-vector to an E-vector or apply its transpose 1020 1021 @param[in] rstr CeedElemRestriction 1022 @param[in] t_mode Apply restriction or transpose 1023 @param[in] u Input vector (of size @a l_size when t_mode=@ref CEED_NOTRANSPOSE) 1024 @param[out] ru Output vector (of shape [@a num_elem * @a elem_size] when t_mode=@ref CEED_NOTRANSPOSE). 1025 Ordering of the e-vector is decided by the backend. 1026 @param[in] request Request or @ref CEED_REQUEST_IMMEDIATE 1027 1028 @return An error code: 0 - success, otherwise - failure 1029 1030 @ref User 1031 **/ 1032 int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, CeedRequest *request) { 1033 CeedInt m, n; 1034 1035 if (t_mode == CEED_NOTRANSPOSE) { 1036 m = rstr->num_block * rstr->block_size * rstr->elem_size * rstr->num_comp; 1037 n = rstr->l_size; 1038 } else { 1039 m = rstr->l_size; 1040 n = rstr->num_block * rstr->block_size * rstr->elem_size * rstr->num_comp; 1041 } 1042 CeedCheck(n == u->length, rstr->ceed, CEED_ERROR_DIMENSION, 1043 "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", u->length, m, n); 1044 CeedCheck(m == ru->length, rstr->ceed, CEED_ERROR_DIMENSION, 1045 "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", ru->length, m, n); 1046 if (rstr->num_elem > 0) CeedCall(rstr->Apply(rstr, t_mode, u, ru, request)); 1047 return CEED_ERROR_SUCCESS; 1048 } 1049 1050 /** 1051 @brief Restrict an L-vector to a block of an E-vector or apply its transpose 1052 1053 @param[in] rstr CeedElemRestriction 1054 @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 1055 [3*block_size : 4*block_size] 1056 @param[in] t_mode Apply restriction or transpose 1057 @param[in] u Input vector (of size @a l_size when t_mode=@ref CEED_NOTRANSPOSE) 1058 @param[out] ru Output vector (of shape [@a block_size * @a elem_size] when t_mode=@ref CEED_NOTRANSPOSE). 1059 Ordering of the e-vector is decided by the backend. 1060 @param[in] request Request or @ref CEED_REQUEST_IMMEDIATE 1061 1062 @return An error code: 0 - success, otherwise - failure 1063 1064 @ref Backend 1065 **/ 1066 int CeedElemRestrictionApplyBlock(CeedElemRestriction rstr, CeedInt block, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, 1067 CeedRequest *request) { 1068 CeedInt m, n; 1069 1070 CeedCheck(rstr->ApplyBlock, rstr->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionApplyBlock"); 1071 1072 if (t_mode == CEED_NOTRANSPOSE) { 1073 m = rstr->block_size * rstr->elem_size * rstr->num_comp; 1074 n = rstr->l_size; 1075 } else { 1076 m = rstr->l_size; 1077 n = rstr->block_size * rstr->elem_size * rstr->num_comp; 1078 } 1079 CeedCheck(n == u->length, rstr->ceed, CEED_ERROR_DIMENSION, 1080 "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", u->length, m, n); 1081 CeedCheck(m == ru->length, rstr->ceed, CEED_ERROR_DIMENSION, 1082 "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", ru->length, m, n); 1083 CeedCheck(rstr->block_size * block <= rstr->num_elem, rstr->ceed, CEED_ERROR_DIMENSION, 1084 "Cannot retrieve block %" CeedInt_FMT ", element %" CeedInt_FMT " > total elements %" CeedInt_FMT "", block, rstr->block_size * block, 1085 rstr->num_elem); 1086 CeedCall(rstr->ApplyBlock(rstr, block, t_mode, u, ru, request)); 1087 return CEED_ERROR_SUCCESS; 1088 } 1089 1090 /** 1091 @brief Get the Ceed associated with a CeedElemRestriction 1092 1093 @param[in] rstr CeedElemRestriction 1094 @param[out] ceed Variable to store Ceed 1095 1096 @return An error code: 0 - success, otherwise - failure 1097 1098 @ref Advanced 1099 **/ 1100 int CeedElemRestrictionGetCeed(CeedElemRestriction rstr, Ceed *ceed) { 1101 *ceed = rstr->ceed; 1102 return CEED_ERROR_SUCCESS; 1103 } 1104 1105 /** 1106 @brief Get the L-vector component stride 1107 1108 @param[in] rstr CeedElemRestriction 1109 @param[out] comp_stride Variable to store component stride 1110 1111 @return An error code: 0 - success, otherwise - failure 1112 1113 @ref Advanced 1114 **/ 1115 int CeedElemRestrictionGetCompStride(CeedElemRestriction rstr, CeedInt *comp_stride) { 1116 *comp_stride = rstr->comp_stride; 1117 return CEED_ERROR_SUCCESS; 1118 } 1119 1120 /** 1121 @brief Get the total number of elements in the range of a CeedElemRestriction 1122 1123 @param[in] rstr CeedElemRestriction 1124 @param[out] num_elem Variable to store number of elements 1125 1126 @return An error code: 0 - success, otherwise - failure 1127 1128 @ref Advanced 1129 **/ 1130 int CeedElemRestrictionGetNumElements(CeedElemRestriction rstr, CeedInt *num_elem) { 1131 *num_elem = rstr->num_elem; 1132 return CEED_ERROR_SUCCESS; 1133 } 1134 1135 /** 1136 @brief Get the size of elements in the CeedElemRestriction 1137 1138 @param[in] rstr CeedElemRestriction 1139 @param[out] elem_size Variable to store size of elements 1140 1141 @return An error code: 0 - success, otherwise - failure 1142 1143 @ref Advanced 1144 **/ 1145 int CeedElemRestrictionGetElementSize(CeedElemRestriction rstr, CeedInt *elem_size) { 1146 *elem_size = rstr->elem_size; 1147 return CEED_ERROR_SUCCESS; 1148 } 1149 1150 /** 1151 @brief Get the size of the l-vector for a CeedElemRestriction 1152 1153 @param[in] rstr CeedElemRestriction 1154 @param[out] l_size Variable to store number of nodes 1155 1156 @return An error code: 0 - success, otherwise - failure 1157 1158 @ref Advanced 1159 **/ 1160 int CeedElemRestrictionGetLVectorSize(CeedElemRestriction rstr, CeedSize *l_size) { 1161 *l_size = rstr->l_size; 1162 return CEED_ERROR_SUCCESS; 1163 } 1164 1165 /** 1166 @brief Get the number of components in the elements of a CeedElemRestriction 1167 1168 @param[in] rstr CeedElemRestriction 1169 @param[out] num_comp Variable to store number of components 1170 1171 @return An error code: 0 - success, otherwise - failure 1172 1173 @ref Advanced 1174 **/ 1175 int CeedElemRestrictionGetNumComponents(CeedElemRestriction rstr, CeedInt *num_comp) { 1176 *num_comp = rstr->num_comp; 1177 return CEED_ERROR_SUCCESS; 1178 } 1179 1180 /** 1181 @brief Get the number of blocks in a CeedElemRestriction 1182 1183 @param[in] rstr CeedElemRestriction 1184 @param[out] num_block Variable to store number of blocks 1185 1186 @return An error code: 0 - success, otherwise - failure 1187 1188 @ref Advanced 1189 **/ 1190 int CeedElemRestrictionGetNumBlocks(CeedElemRestriction rstr, CeedInt *num_block) { 1191 *num_block = rstr->num_block; 1192 return CEED_ERROR_SUCCESS; 1193 } 1194 1195 /** 1196 @brief Get the size of blocks in the CeedElemRestriction 1197 1198 @param[in] rstr CeedElemRestriction 1199 @param[out] block_size Variable to store size of blocks 1200 1201 @return An error code: 0 - success, otherwise - failure 1202 1203 @ref Advanced 1204 **/ 1205 int CeedElemRestrictionGetBlockSize(CeedElemRestriction rstr, CeedInt *block_size) { 1206 *block_size = rstr->block_size; 1207 return CEED_ERROR_SUCCESS; 1208 } 1209 1210 /** 1211 @brief Get the multiplicity of nodes in a CeedElemRestriction 1212 1213 @param[in] rstr CeedElemRestriction 1214 @param[out] mult Vector to store multiplicity (of size l_size) 1215 1216 @return An error code: 0 - success, otherwise - failure 1217 1218 @ref User 1219 **/ 1220 int CeedElemRestrictionGetMultiplicity(CeedElemRestriction rstr, CeedVector mult) { 1221 CeedVector e_vec; 1222 1223 // Create e_vec to hold intermediate computation in E^T (E 1) 1224 CeedCall(CeedElemRestrictionCreateVector(rstr, NULL, &e_vec)); 1225 1226 // Compute e_vec = E * 1 1227 CeedCall(CeedVectorSetValue(mult, 1.0)); 1228 CeedCall(CeedElemRestrictionApply(rstr, CEED_NOTRANSPOSE, mult, e_vec, CEED_REQUEST_IMMEDIATE)); 1229 // Compute multiplicity, mult = E^T * e_vec = E^T (E 1) 1230 CeedCall(CeedVectorSetValue(mult, 0.0)); 1231 CeedCall(CeedElemRestrictionApply(rstr, CEED_TRANSPOSE, e_vec, mult, CEED_REQUEST_IMMEDIATE)); 1232 // Cleanup 1233 CeedCall(CeedVectorDestroy(&e_vec)); 1234 return CEED_ERROR_SUCCESS; 1235 } 1236 1237 /** 1238 @brief View a CeedElemRestriction 1239 1240 @param[in] rstr CeedElemRestriction to view 1241 @param[in] stream Stream to write; typically stdout/stderr or a file 1242 1243 @return Error code: 0 - success, otherwise - failure 1244 1245 @ref User 1246 **/ 1247 int CeedElemRestrictionView(CeedElemRestriction rstr, FILE *stream) { 1248 char stridesstr[500]; 1249 if (rstr->strides) { 1250 sprintf(stridesstr, "[%" CeedInt_FMT ", %" CeedInt_FMT ", %" CeedInt_FMT "]", rstr->strides[0], rstr->strides[1], rstr->strides[2]); 1251 } else { 1252 sprintf(stridesstr, "%" CeedInt_FMT, rstr->comp_stride); 1253 } 1254 1255 fprintf(stream, "%sCeedElemRestriction from (%td, %" CeedInt_FMT ") to %" CeedInt_FMT " elements with %" CeedInt_FMT " nodes each and %s %s\n", 1256 rstr->block_size > 1 ? "Blocked " : "", rstr->l_size, rstr->num_comp, rstr->num_elem, rstr->elem_size, 1257 rstr->strides ? "strides" : "component stride", stridesstr); 1258 return CEED_ERROR_SUCCESS; 1259 } 1260 1261 /** 1262 @brief Destroy a CeedElemRestriction 1263 1264 @param[in,out] rstr CeedElemRestriction to destroy 1265 1266 @return An error code: 0 - success, otherwise - failure 1267 1268 @ref User 1269 **/ 1270 int CeedElemRestrictionDestroy(CeedElemRestriction *rstr) { 1271 if (!*rstr || *rstr == CEED_ELEMRESTRICTION_NONE || --(*rstr)->ref_count > 0) { 1272 *rstr = NULL; 1273 return CEED_ERROR_SUCCESS; 1274 } 1275 CeedCheck((*rstr)->num_readers == 0, (*rstr)->ceed, CEED_ERROR_ACCESS, 1276 "Cannot destroy CeedElemRestriction, a process has read access to the offset data"); 1277 1278 // Only destroy backend data once between rstr and unsigned copy 1279 if ((*rstr)->rstr_base) CeedCall(CeedElemRestrictionDestroy(&(*rstr)->rstr_base)); 1280 else if ((*rstr)->Destroy) CeedCall((*rstr)->Destroy(*rstr)); 1281 1282 CeedCall(CeedFree(&(*rstr)->strides)); 1283 CeedCall(CeedDestroy(&(*rstr)->ceed)); 1284 CeedCall(CeedFree(rstr)); 1285 return CEED_ERROR_SUCCESS; 1286 } 1287 1288 /// @} 1289