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