// Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. // // SPDX-License-Identifier: BSD-2-Clause // // This file is part of CEED: http://github.com/ceed #include #include #include #include #include #include /// @file /// Implementation of CeedElemRestriction interfaces /// ---------------------------------------------------------------------------- /// CeedElemRestriction Library Internal Functions /// ---------------------------------------------------------------------------- /// @addtogroup CeedElemRestrictionDeveloper /// @{ /** @brief Permute and pad offsets for a blocked restriction @param[in] offsets Array of shape [@a num_elem, @a elem_size]. @param[out] block_offsets Array of permuted and padded array values of shape [@a num_block, @a elem_size, @a block_size]. @param[in] num_block Number of blocks @param[in] num_elem Number of elements @param[in] block_size Number of elements in a block @param[in] elem_size Size of each element @return An error code: 0 - success, otherwise - failure @ref Utility **/ int CeedPermutePadOffsets(const CeedInt *offsets, CeedInt *block_offsets, CeedInt num_block, CeedInt num_elem, CeedInt block_size, CeedInt elem_size) { for (CeedInt e = 0; e < num_block * block_size; e += block_size) { for (CeedInt j = 0; j < block_size; j++) { for (CeedInt k = 0; k < elem_size; k++) { block_offsets[e * elem_size + k * block_size + j] = offsets[CeedIntMin(e + j, num_elem - 1) * elem_size + k]; } } } return CEED_ERROR_SUCCESS; } /** @brief Permute and pad orientations for a blocked restriction @param[in] orients Array of shape [@a num_elem, @a elem_size]. @param[out] block_orients Array of permuted and padded array values of shape [@a num_block, @a elem_size, @a block_size]. @param[in] num_block Number of blocks @param[in] num_elem Number of elements @param[in] block_size Number of elements in a block @param[in] elem_size Size of each element @return An error code: 0 - success, otherwise - failure @ref Utility **/ int CeedPermutePadOrients(const bool *orients, bool *block_orients, CeedInt num_block, CeedInt num_elem, CeedInt block_size, CeedInt elem_size) { for (CeedInt e = 0; e < num_block * block_size; e += block_size) { for (CeedInt j = 0; j < block_size; j++) { for (CeedInt k = 0; k < elem_size; k++) { block_orients[e * elem_size + k * block_size + j] = orients[CeedIntMin(e + j, num_elem - 1) * elem_size + k]; } } } return CEED_ERROR_SUCCESS; } /** @brief Permute and pad curl-conforming orientations for a blocked restriction @param[in] curl_orients Array of shape [@a num_elem, @a 3 * elem_size]. @param[out] block_curl_orients Array of permuted and padded array values of shape [@a num_block, @a elem_size, @a block_size]. @param[in] num_block Number of blocks @param[in] num_elem Number of elements @param[in] block_size Number of elements in a block @param[in] elem_size Size of each element @return An error code: 0 - success, otherwise - failure @ref Utility **/ int CeedPermutePadCurlOrients(const CeedInt8 *curl_orients, CeedInt8 *block_curl_orients, CeedInt num_block, CeedInt num_elem, CeedInt block_size, CeedInt elem_size) { for (CeedInt e = 0; e < num_block * block_size; e += block_size) { for (CeedInt j = 0; j < block_size; j++) { for (CeedInt k = 0; k < elem_size; k++) { block_curl_orients[e * elem_size + k * block_size + j] = curl_orients[CeedIntMin(e + j, num_elem - 1) * elem_size + k]; } } } return CEED_ERROR_SUCCESS; } /// @} /// ---------------------------------------------------------------------------- /// CeedElemRestriction Backend API /// ---------------------------------------------------------------------------- /// @addtogroup CeedElemRestrictionBackend /// @{ /** @brief Get the type of a CeedElemRestriction @param[in] rstr CeedElemRestriction @param[out] rstr_type Variable to store restriction type @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedElemRestrictionGetType(CeedElemRestriction rstr, CeedRestrictionType *rstr_type) { *rstr_type = rstr->rstr_type; return CEED_ERROR_SUCCESS; } /** @brief Get the strided status of a CeedElemRestriction @param[in] rstr CeedElemRestriction @param[out] is_strided Variable to store strided status, 1 if strided else 0 **/ int CeedElemRestrictionIsStrided(CeedElemRestriction rstr, bool *is_strided) { *is_strided = (rstr->rstr_type == CEED_RESTRICTION_STRIDED); return CEED_ERROR_SUCCESS; } /** @brief Get the points status of a CeedElemRestriction @param[in] rstr CeedElemRestriction @param[out] is_points Variable to store points status, 1 if points else 0 **/ int CeedElemRestrictionIsPoints(CeedElemRestriction rstr, bool *is_points) { *is_points = (rstr->rstr_type == CEED_RESTRICTION_POINTS); return CEED_ERROR_SUCCESS; } /** @brief Get the strides of a strided CeedElemRestriction @param[in] rstr CeedElemRestriction @param[out] strides Variable to store strides array @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedElemRestrictionGetStrides(CeedElemRestriction rstr, CeedInt (*strides)[3]) { CeedCheck(rstr->strides, rstr->ceed, CEED_ERROR_MINOR, "ElemRestriction has no stride data"); for (CeedInt i = 0; i < 3; i++) (*strides)[i] = rstr->strides[i]; return CEED_ERROR_SUCCESS; } /** @brief Get the backend stride status of a CeedElemRestriction @param[in] rstr CeedElemRestriction @param[out] has_backend_strides Variable to store stride status @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedElemRestrictionHasBackendStrides(CeedElemRestriction rstr, bool *has_backend_strides) { CeedCheck(rstr->strides, rstr->ceed, CEED_ERROR_MINOR, "ElemRestriction has no stride data"); *has_backend_strides = ((rstr->strides[0] == CEED_STRIDES_BACKEND[0]) && (rstr->strides[1] == CEED_STRIDES_BACKEND[1]) && (rstr->strides[2] == CEED_STRIDES_BACKEND[2])); return CEED_ERROR_SUCCESS; } /** @brief Get read-only access to a CeedElemRestriction offsets array by memtype @param[in] rstr CeedElemRestriction to retrieve offsets @param[in] mem_type Memory type on which to access the array. If the backend uses a different memory type, this will perform a copy (possibly cached). @param[out] offsets Array on memory type mem_type @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionGetOffsets(CeedElemRestriction rstr, CeedMemType mem_type, const CeedInt **offsets) { if (rstr->rstr_base) { CeedCall(CeedElemRestrictionGetOffsets(rstr->rstr_base, mem_type, offsets)); } else { CeedCheck(rstr->GetOffsets, rstr->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support GetOffsets"); CeedCall(rstr->GetOffsets(rstr, mem_type, offsets)); rstr->num_readers++; } return CEED_ERROR_SUCCESS; } /** @brief Restore an offsets array obtained using CeedElemRestrictionGetOffsets() @param[in] rstr CeedElemRestriction to restore @param[in] offsets Array of offset data @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionRestoreOffsets(CeedElemRestriction rstr, const CeedInt **offsets) { if (rstr->rstr_base) { CeedCall(CeedElemRestrictionRestoreOffsets(rstr->rstr_base, offsets)); } else { *offsets = NULL; rstr->num_readers--; } return CEED_ERROR_SUCCESS; } /** @brief Get read-only access to a CeedElemRestriction orientations array by memtype @param[in] rstr CeedElemRestriction to retrieve orientations @param[in] mem_type Memory type on which to access the array. If the backend uses a different memory type, this will perform a copy (possibly cached). @param[out] orients Array on memory type mem_type @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionGetOrientations(CeedElemRestriction rstr, CeedMemType mem_type, const bool **orients) { CeedCheck(rstr->GetOrientations, rstr->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support GetOrientations"); CeedCall(rstr->GetOrientations(rstr, mem_type, orients)); rstr->num_readers++; return CEED_ERROR_SUCCESS; } /** @brief Restore an orientations array obtained using CeedElemRestrictionGetOrientations() @param[in] rstr CeedElemRestriction to restore @param[in] orients Array of orientation data @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionRestoreOrientations(CeedElemRestriction rstr, const bool **orients) { *orients = NULL; rstr->num_readers--; return CEED_ERROR_SUCCESS; } /** @brief Get read-only access to a CeedElemRestriction curl-conforming orientations array by memtype @param[in] rstr CeedElemRestriction to retrieve curl-conforming orientations @param[in] mem_type Memory type on which to access the array. If the backend uses a different memory type, this will perform a copy (possibly cached). @param[out] curl_orients Array on memory type mem_type @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionGetCurlOrientations(CeedElemRestriction rstr, CeedMemType mem_type, const CeedInt8 **curl_orients) { CeedCheck(rstr->GetCurlOrientations, rstr->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support GetCurlOrientations"); CeedCall(rstr->GetCurlOrientations(rstr, mem_type, curl_orients)); rstr->num_readers++; return CEED_ERROR_SUCCESS; } /** @brief Restore an orientations array obtained using CeedElemRestrictionGetCurlOrientations() @param[in] rstr CeedElemRestriction to restore @param[in] curl_orients Array of orientation data @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionRestoreCurlOrientations(CeedElemRestriction rstr, const CeedInt8 **curl_orients) { *curl_orients = NULL; rstr->num_readers--; return CEED_ERROR_SUCCESS; } /** @brief Get the E-vector layout of a CeedElemRestriction @param[in] rstr CeedElemRestriction @param[out] layout Variable to store layout array, stored as [nodes, components, elements]. 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] @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedElemRestrictionGetELayout(CeedElemRestriction rstr, CeedInt (*layout)[3]) { CeedCheck(rstr->layout[0], rstr->ceed, CEED_ERROR_MINOR, "ElemRestriction has no layout data"); for (CeedInt i = 0; i < 3; i++) (*layout)[i] = rstr->layout[i]; return CEED_ERROR_SUCCESS; } /** @brief Set the E-vector layout of a CeedElemRestriction @param[in] rstr CeedElemRestriction @param[in] layout Variable to containing layout array, stored as [nodes, components, elements]. 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] @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedElemRestrictionSetELayout(CeedElemRestriction rstr, CeedInt layout[3]) { for (CeedInt i = 0; i < 3; i++) rstr->layout[i] = layout[i]; return CEED_ERROR_SUCCESS; } /** @brief Get the backend data of a CeedElemRestriction @param[in] rstr CeedElemRestriction @param[out] data Variable to store data @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedElemRestrictionGetData(CeedElemRestriction rstr, void *data) { *(void **)data = rstr->data; return CEED_ERROR_SUCCESS; } /** @brief Set the backend data of a CeedElemRestriction @param[in,out] rstr CeedElemRestriction @param[in] data Data to set @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedElemRestrictionSetData(CeedElemRestriction rstr, void *data) { rstr->data = data; return CEED_ERROR_SUCCESS; } /** @brief Increment the reference counter for a CeedElemRestriction @param[in,out] rstr ElemRestriction to increment the reference counter @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedElemRestrictionReference(CeedElemRestriction rstr) { rstr->ref_count++; return CEED_ERROR_SUCCESS; } /** @brief Estimate number of FLOPs required to apply CeedElemRestriction in t_mode @param[in] rstr ElemRestriction to estimate FLOPs for @param[in] t_mode Apply restriction or transpose @param[out] flops Address of variable to hold FLOPs estimate @ref Backend **/ int CeedElemRestrictionGetFlopsEstimate(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedSize *flops) { CeedInt e_size = rstr->num_block * rstr->block_size * rstr->elem_size * rstr->num_comp, scale = 0; CeedRestrictionType rstr_type; CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); if (rstr_type == CEED_RESTRICTION_POINTS) e_size = rstr->num_points * rstr->num_comp; if (t_mode == CEED_TRANSPOSE) { switch (rstr_type) { case CEED_RESTRICTION_POINTS: scale = 0; break; case CEED_RESTRICTION_STRIDED: case CEED_RESTRICTION_STANDARD: scale = 1; break; case CEED_RESTRICTION_ORIENTED: scale = 2; break; case CEED_RESTRICTION_CURL_ORIENTED: scale = 6; break; } } else { switch (rstr_type) { case CEED_RESTRICTION_STRIDED: case CEED_RESTRICTION_STANDARD: case CEED_RESTRICTION_POINTS: scale = 0; break; case CEED_RESTRICTION_ORIENTED: scale = 1; break; case CEED_RESTRICTION_CURL_ORIENTED: scale = 5; break; } } *flops = e_size * scale; return CEED_ERROR_SUCCESS; } /// @} /// @cond DOXYGEN_SKIP static struct CeedElemRestriction_private ceed_elemrestriction_none; /// @endcond /// ---------------------------------------------------------------------------- /// CeedElemRestriction Public API /// ---------------------------------------------------------------------------- /// @addtogroup CeedElemRestrictionUser /// @{ /// Indicate that the stride is determined by the backend const CeedInt CEED_STRIDES_BACKEND[3] = {0}; /// Argument for CeedOperatorSetField indicating that the field does not requre a CeedElemRestriction const CeedElemRestriction CEED_ELEMRESTRICTION_NONE = &ceed_elemrestriction_none; /** @brief Create a CeedElemRestriction @param[in] ceed Ceed object where the CeedElemRestriction will be created @param[in] num_elem Number of elements described in the @a offsets array @param[in] elem_size Size (number of "nodes") per element @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) @param[in] comp_stride Stride between components for the same L-vector "node". 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. @param[in] l_size The size of the L-vector. This vector may be larger than the elements and fields given by this restriction. @param[in] mem_type Memory type of the @a offsets array, see CeedMemType @param[in] copy_mode Copy mode for the @a offsets array, see CeedCopyMode @param[in] offsets Array of shape [@a num_elem, @a elem_size]. Row i holds the ordered list of the offsets (into the input CeedVector) for the unknowns corresponding to element i, where 0 <= i < @a num_elem. All offsets must be in the range [0, @a l_size - 1]. @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionCreate(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt num_comp, CeedInt comp_stride, CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, CeedElemRestriction *rstr) { if (!ceed->ElemRestrictionCreate) { Ceed delegate; CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreate"); CeedCall(CeedElemRestrictionCreate(delegate, num_elem, elem_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, rstr)); return CEED_ERROR_SUCCESS; } CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction component stride must be at least 1"); CeedCall(CeedCalloc(1, rstr)); CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); (*rstr)->ref_count = 1; (*rstr)->num_elem = num_elem; (*rstr)->elem_size = elem_size; (*rstr)->num_comp = num_comp; (*rstr)->comp_stride = comp_stride; (*rstr)->l_size = l_size; (*rstr)->e_size = num_elem * elem_size * num_comp; (*rstr)->num_block = num_elem; (*rstr)->block_size = 1; (*rstr)->rstr_type = CEED_RESTRICTION_STANDARD; CeedCall(ceed->ElemRestrictionCreate(mem_type, copy_mode, offsets, NULL, NULL, *rstr)); return CEED_ERROR_SUCCESS; } /** @brief Create a CeedElemRestriction with orientation signs @param[in] ceed Ceed object where the CeedElemRestriction will be created @param[in] num_elem Number of elements described in the @a offsets array @param[in] elem_size Size (number of "nodes") per element @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) @param[in] comp_stride Stride between components for the same L-vector "node". 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. @param[in] l_size The size of the L-vector. This vector may be larger than the elements and fields given by this restriction. @param[in] mem_type Memory type of the @a offsets array, see CeedMemType @param[in] copy_mode Copy mode for the @a offsets array, see CeedCopyMode @param[in] offsets Array of shape [@a num_elem, @a elem_size]. Row i holds the ordered list of the offsets (into the input CeedVector) for the unknowns corresponding to element i, where 0 <= i < @a num_elem. All offsets must be in the range [0, @a l_size - 1]. @param[in] orients Array of shape [@a num_elem, @a elem_size] with bool false for positively oriented and true to flip the orientation. @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionCreateOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt num_comp, CeedInt comp_stride, CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, const bool *orients, CeedElemRestriction *rstr) { if (!ceed->ElemRestrictionCreate) { Ceed delegate; CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreate"); CeedCall( CeedElemRestrictionCreateOriented(delegate, num_elem, elem_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, orients, rstr)); return CEED_ERROR_SUCCESS; } CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction component stride must be at least 1"); CeedCall(CeedCalloc(1, rstr)); CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); (*rstr)->ref_count = 1; (*rstr)->num_elem = num_elem; (*rstr)->elem_size = elem_size; (*rstr)->num_comp = num_comp; (*rstr)->comp_stride = comp_stride; (*rstr)->l_size = l_size; (*rstr)->e_size = num_elem * elem_size * num_comp; (*rstr)->num_block = num_elem; (*rstr)->block_size = 1; (*rstr)->rstr_type = CEED_RESTRICTION_ORIENTED; CeedCall(ceed->ElemRestrictionCreate(mem_type, copy_mode, offsets, orients, NULL, *rstr)); return CEED_ERROR_SUCCESS; } /** @brief Create a CeedElemRestriction with a general tridiagonal transformation matrix for curl-conforming elements @param[in] ceed Ceed object where the CeedElemRestriction will be created @param[in] num_elem Number of elements described in the @a offsets array @param[in] elem_size Size (number of "nodes") per element @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) @param[in] comp_stride Stride between components for the same L-vector "node". 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. @param[in] l_size The size of the L-vector. This vector may be larger than the elements and fields given by this restriction. @param[in] mem_type Memory type of the @a offsets array, see CeedMemType @param[in] copy_mode Copy mode for the @a offsets array, see CeedCopyMode @param[in] offsets Array of shape [@a num_elem, @a elem_size]. Row i holds the ordered list of the offsets (into the input CeedVector) for the unknowns corresponding to element i, where 0 <= i < @a num_elem. All offsets must be in the range [0, @a l_size - 1]. @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] = 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 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 to resolve face orientation issues for 3D meshes (https://dl.acm.org/doi/pdf/10.1145/3524456). @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionCreateCurlOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt num_comp, CeedInt comp_stride, CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, const CeedInt8 *curl_orients, CeedElemRestriction *rstr) { if (!ceed->ElemRestrictionCreate) { Ceed delegate; CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreate"); CeedCall(CeedElemRestrictionCreateCurlOriented(delegate, num_elem, elem_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, curl_orients, rstr)); return CEED_ERROR_SUCCESS; } CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction component stride must be at least 1"); CeedCall(CeedCalloc(1, rstr)); CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); (*rstr)->ref_count = 1; (*rstr)->num_elem = num_elem; (*rstr)->elem_size = elem_size; (*rstr)->num_comp = num_comp; (*rstr)->comp_stride = comp_stride; (*rstr)->l_size = l_size; (*rstr)->e_size = num_elem * elem_size * num_comp; (*rstr)->num_block = num_elem; (*rstr)->block_size = 1; (*rstr)->rstr_type = CEED_RESTRICTION_CURL_ORIENTED; CeedCall(ceed->ElemRestrictionCreate(mem_type, copy_mode, offsets, NULL, curl_orients, *rstr)); return CEED_ERROR_SUCCESS; } /** @brief Create a strided CeedElemRestriction @param[in] ceed Ceed object where the CeedElemRestriction will be created @param[in] num_elem Number of elements described by the restriction @param[in] elem_size Size (number of "nodes") per element @param[in] num_comp Number of field components per interpolation "node" (1 for scalar fields) @param[in] l_size The size of the L-vector. This vector may be larger than the elements and fields given by this restriction. @param[in] strides Array for strides between [nodes, components, elements]. 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]. @a CEED_STRIDES_BACKEND may be used with vectors created by a Ceed backend. @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionCreateStrided(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt num_comp, CeedSize l_size, const CeedInt strides[3], CeedElemRestriction *rstr) { if (!ceed->ElemRestrictionCreate) { Ceed delegate; CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreate"); CeedCall(CeedElemRestrictionCreateStrided(delegate, num_elem, elem_size, num_comp, l_size, strides, rstr)); return CEED_ERROR_SUCCESS; } CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); 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"); CeedCall(CeedCalloc(1, rstr)); CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); (*rstr)->ref_count = 1; (*rstr)->num_elem = num_elem; (*rstr)->elem_size = elem_size; (*rstr)->num_comp = num_comp; (*rstr)->l_size = l_size; (*rstr)->e_size = num_elem * elem_size * num_comp; (*rstr)->num_block = num_elem; (*rstr)->block_size = 1; (*rstr)->rstr_type = CEED_RESTRICTION_STRIDED; CeedCall(CeedMalloc(3, &(*rstr)->strides)); for (CeedInt i = 0; i < 3; i++) (*rstr)->strides[i] = strides[i]; CeedCall(ceed->ElemRestrictionCreate(CEED_MEM_HOST, CEED_OWN_POINTER, NULL, NULL, NULL, *rstr)); return CEED_ERROR_SUCCESS; } /** @brief Create a points CeedElemRestriction, for restricting for restricting from a all local points to the current element in which they are located. The offsets array is arranged as element_0_start_index element_1_start_index ... element_n_start_index element_n_stop_index element_0_point_0 element_0_point_1 ... @param[in] ceed Ceed object where the CeedElemRestriction will be created @param[in] num_elem Number of elements described in the @a offsets array @param[in] num_points Number of points described in the @a offsets array @param[in] num_comp Number of field components per interpolation node (1 for scalar fields). Components are assumed to be contiguous by point. @param[in] l_size The size of the L-vector. This vector may be larger than the elements and fields given by this restriction. @param[in] mem_type Memory type of the @a offsets array, see CeedMemType @param[in] copy_mode Copy mode for the @a offsets array, see CeedCopyMode @param[in] offsets Array of size num_elem + 1 + num_points. The first portion of the offsets array holds the ranges of indices corresponding to each element. The second portion holds the indices for each element. @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedElemRestrictionCreateAtPoints(Ceed ceed, CeedInt num_elem, CeedInt num_points, CeedInt num_comp, CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, CeedElemRestriction *rstr) { if (!ceed->ElemRestrictionCreateAtPoints) { Ceed delegate; CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreateAtPoints"); CeedCall(CeedElemRestrictionCreateAtPoints(delegate, num_elem, num_points, num_comp, l_size, mem_type, copy_mode, offsets, rstr)); return CEED_ERROR_SUCCESS; } CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); CeedCheck(num_points >= 0, ceed, CEED_ERROR_DIMENSION, "Number of points must be non-negative"); CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); CeedCheck(l_size >= num_points * num_comp, ceed, CEED_ERROR_DIMENSION, "L-vector must be at least num_points * num_comp"); CeedCall(CeedCalloc(1, rstr)); CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); (*rstr)->ref_count = 1; (*rstr)->num_elem = num_elem; (*rstr)->num_points = num_points; (*rstr)->num_comp = num_comp; (*rstr)->l_size = l_size; (*rstr)->e_size = num_points * num_comp; (*rstr)->num_block = num_elem; (*rstr)->block_size = 1; (*rstr)->rstr_type = CEED_RESTRICTION_POINTS; CeedCall(ceed->ElemRestrictionCreateAtPoints(mem_type, copy_mode, offsets, NULL, NULL, *rstr)); return CEED_ERROR_SUCCESS; } /** @brief Create a blocked CeedElemRestriction, typically only called by backends @param[in] ceed Ceed object where the CeedElemRestriction will be created @param[in] num_elem Number of elements described in the @a offsets array @param[in] elem_size Size (number of unknowns) per element @param[in] block_size Number of elements in a block @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) @param[in] comp_stride Stride between components for the same L-vector "node". 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. @param[in] l_size The size of the L-vector. This vector may be larger than the elements and fields given by this restriction. @param[in] mem_type Memory type of the @a offsets array, see CeedMemType @param[in] copy_mode Copy mode for the @a offsets array, see CeedCopyMode @param[in] offsets Array of shape [@a num_elem, @a elem_size]. Row i holds the ordered list of the offsets (into the input CeedVector) for the unknowns corresponding to element i, 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 for the blocksize, which is typically given by the backend. The default reordering is to interlace elements. @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedElemRestrictionCreateBlocked(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, CeedInt comp_stride, CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, CeedElemRestriction *rstr) { CeedInt *block_offsets, num_block = (num_elem / block_size) + !!(num_elem % block_size); if (!ceed->ElemRestrictionCreateBlocked) { Ceed delegate; CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreateBlocked"); CeedCall(CeedElemRestrictionCreateBlocked(delegate, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, rstr)); return CEED_ERROR_SUCCESS; } CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction component stride must be at least 1"); CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_offsets)); CeedCall(CeedPermutePadOffsets(offsets, block_offsets, num_block, num_elem, block_size, elem_size)); CeedCall(CeedCalloc(1, rstr)); CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); (*rstr)->ref_count = 1; (*rstr)->num_elem = num_elem; (*rstr)->elem_size = elem_size; (*rstr)->num_comp = num_comp; (*rstr)->comp_stride = comp_stride; (*rstr)->l_size = l_size; (*rstr)->e_size = num_block * block_size * elem_size * num_comp; (*rstr)->num_block = num_block; (*rstr)->block_size = block_size; (*rstr)->rstr_type = CEED_RESTRICTION_STANDARD; CeedCall(ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, (const CeedInt *)block_offsets, NULL, NULL, *rstr)); if (copy_mode == CEED_OWN_POINTER) CeedCall(CeedFree(&offsets)); return CEED_ERROR_SUCCESS; } /** @brief Create a blocked oriented CeedElemRestriction, typically only called by backends @param[in] ceed Ceed object where the CeedElemRestriction will be created. @param[in] num_elem Number of elements described in the @a offsets array. @param[in] elem_size Size (number of unknowns) per element @param[in] block_size Number of elements in a block @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) @param[in] comp_stride Stride between components for the same L-vector "node". 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. @param[in] l_size The size of the L-vector. This vector may be larger than the elements and fields given by this restriction. @param[in] mem_type Memory type of the @a offsets array, see CeedMemType @param[in] copy_mode Copy mode for the @a offsets array, see CeedCopyMode @param[in] offsets Array of shape [@a num_elem, @a elem_size]. Row i holds the ordered list of the offsets (into the input CeedVector) for the unknowns corresponding to element i, 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 for the blocksize, which is typically given by the backend. The default reordering is to interlace elements. @param[in] orients Array of shape [@a num_elem, @a elem_size] with bool false for positively oriented and true to flip the orientation. Will also be permuted and padded similarly to @a offsets. @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedElemRestrictionCreateBlockedOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, CeedInt comp_stride, CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, const bool *orients, CeedElemRestriction *rstr) { bool *block_orients; CeedInt *block_offsets, num_block = (num_elem / block_size) + !!(num_elem % block_size); if (!ceed->ElemRestrictionCreateBlocked) { Ceed delegate; CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreateBlocked"); CeedCall(CeedElemRestrictionCreateBlockedOriented(delegate, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, orients, rstr)); return CEED_ERROR_SUCCESS; } CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction component stride must be at least 1"); CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_offsets)); CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_orients)); CeedCall(CeedPermutePadOffsets(offsets, block_offsets, num_block, num_elem, block_size, elem_size)); CeedCall(CeedPermutePadOrients(orients, block_orients, num_block, num_elem, block_size, elem_size)); CeedCall(CeedCalloc(1, rstr)); CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); (*rstr)->ref_count = 1; (*rstr)->num_elem = num_elem; (*rstr)->elem_size = elem_size; (*rstr)->num_comp = num_comp; (*rstr)->comp_stride = comp_stride; (*rstr)->l_size = l_size; (*rstr)->e_size = num_block * block_size * elem_size * num_comp; (*rstr)->num_block = num_block; (*rstr)->block_size = block_size; (*rstr)->rstr_type = CEED_RESTRICTION_ORIENTED; CeedCall( ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, (const CeedInt *)block_offsets, (const bool *)block_orients, NULL, *rstr)); if (copy_mode == CEED_OWN_POINTER) CeedCall(CeedFree(&offsets)); return CEED_ERROR_SUCCESS; } /** @brief Create a blocked curl-oriented CeedElemRestriction, typically only called by backends @param[in] ceed Ceed object where the CeedElemRestriction will be created. @param[in] num_elem Number of elements described in the @a offsets array. @param[in] elem_size Size (number of unknowns) per element @param[in] block_size Number of elements in a block @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) @param[in] comp_stride Stride between components for the same L-vector "node". 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. @param[in] l_size The size of the L-vector. This vector may be larger than the elements and fields given by this restriction. @param[in] mem_type Memory type of the @a offsets array, see CeedMemType @param[in] copy_mode Copy mode for the @a offsets array, see CeedCopyMode @param[in] offsets Array of shape [@a num_elem, @a elem_size]. Row i holds the ordered list of the offsets (into the input CeedVector) for the unknowns corresponding to element i, 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 for the blocksize, which is typically given by the backend. The default reordering is to interlace elements. @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] = 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 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 to resolve face orientation issues for 3D meshes (https://dl.acm.org/doi/pdf/10.1145/3524456). Will also be permuted and padded similarly to @a offsets. @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedElemRestrictionCreateBlockedCurlOriented(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, CeedInt comp_stride, CeedSize l_size, CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, const CeedInt8 *curl_orients, CeedElemRestriction *rstr) { CeedInt8 *block_curl_orients; CeedInt *block_offsets, num_block = (num_elem / block_size) + !!(num_elem % block_size); if (!ceed->ElemRestrictionCreateBlocked) { Ceed delegate; CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreateBlocked"); CeedCall(CeedElemRestrictionCreateBlockedCurlOriented(delegate, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, curl_orients, rstr)); return CEED_ERROR_SUCCESS; } CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); CeedCheck(num_comp == 1 || comp_stride > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction component stride must be at least 1"); CeedCall(CeedCalloc(num_block * block_size * elem_size, &block_offsets)); CeedCall(CeedCalloc(num_block * block_size * 3 * elem_size, &block_curl_orients)); CeedCall(CeedPermutePadOffsets(offsets, block_offsets, num_block, num_elem, block_size, elem_size)); CeedCall(CeedPermutePadCurlOrients(curl_orients, block_curl_orients, num_block, num_elem, block_size, 3 * elem_size)); CeedCall(CeedCalloc(1, rstr)); CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); (*rstr)->ref_count = 1; (*rstr)->num_elem = num_elem; (*rstr)->elem_size = elem_size; (*rstr)->num_comp = num_comp; (*rstr)->comp_stride = comp_stride; (*rstr)->l_size = l_size; (*rstr)->e_size = num_block * block_size * elem_size * num_comp; (*rstr)->num_block = num_block; (*rstr)->block_size = block_size; (*rstr)->rstr_type = CEED_RESTRICTION_CURL_ORIENTED; CeedCall(ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, (const CeedInt *)block_offsets, NULL, (const CeedInt8 *)block_curl_orients, *rstr)); if (copy_mode == CEED_OWN_POINTER) CeedCall(CeedFree(&offsets)); return CEED_ERROR_SUCCESS; } /** @brief Create a blocked strided CeedElemRestriction, typically only called by backends @param[in] ceed Ceed object where the CeedElemRestriction will be created @param[in] num_elem Number of elements described by the restriction @param[in] elem_size Size (number of "nodes") per element @param[in] block_size Number of elements in a block @param[in] num_comp Number of field components per interpolation node (1 for scalar fields) @param[in] l_size The size of the L-vector. This vector may be larger than the elements and fields given by this restriction. @param[in] strides Array for strides between [nodes, components, elements]. 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]. @a CEED_STRIDES_BACKEND may be used with vectors created by a Ceed backend. @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionCreateBlockedStrided(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt block_size, CeedInt num_comp, CeedSize l_size, const CeedInt strides[3], CeedElemRestriction *rstr) { CeedInt num_block = (num_elem / block_size) + !!(num_elem % block_size); if (!ceed->ElemRestrictionCreateBlocked) { Ceed delegate; CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionCreateBlocked"); CeedCall(CeedElemRestrictionCreateBlockedStrided(delegate, num_elem, elem_size, block_size, num_comp, l_size, strides, rstr)); return CEED_ERROR_SUCCESS; } CeedCheck(num_elem >= 0, ceed, CEED_ERROR_DIMENSION, "Number of elements must be non-negative"); CeedCheck(elem_size > 0, ceed, CEED_ERROR_DIMENSION, "Element size must be at least 1"); CeedCheck(block_size > 0, ceed, CEED_ERROR_DIMENSION, "Block size must be at least 1"); CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "ElemRestriction must have at least 1 component"); 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"); CeedCall(CeedCalloc(1, rstr)); CeedCall(CeedReferenceCopy(ceed, &(*rstr)->ceed)); (*rstr)->ref_count = 1; (*rstr)->num_elem = num_elem; (*rstr)->elem_size = elem_size; (*rstr)->num_comp = num_comp; (*rstr)->l_size = l_size; (*rstr)->e_size = num_block * block_size * elem_size * num_comp; (*rstr)->num_block = num_block; (*rstr)->block_size = block_size; (*rstr)->rstr_type = CEED_RESTRICTION_STRIDED; CeedCall(CeedMalloc(3, &(*rstr)->strides)); for (CeedInt i = 0; i < 3; i++) (*rstr)->strides[i] = strides[i]; CeedCall(ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, NULL, NULL, NULL, *rstr)); return CEED_ERROR_SUCCESS; } /** @brief Copy the pointer to a CeedElemRestriction and set `CeedElemRestrictionApply()` implementation to use the unsigned version. Both pointers should be destroyed with `CeedElemRestrictionDestroy()`. @param[in] rstr CeedElemRestriction to create unsigned reference to @param[in,out] rstr_unsigned Variable to store unsigned CeedElemRestriction @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionCreateUnsignedCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_unsigned) { CeedCall(CeedCalloc(1, rstr_unsigned)); // Copy old rstr memcpy(*rstr_unsigned, rstr, sizeof(struct CeedElemRestriction_private)); (*rstr_unsigned)->ceed = NULL; CeedCall(CeedReferenceCopy(rstr->ceed, &(*rstr_unsigned)->ceed)); (*rstr_unsigned)->ref_count = 1; (*rstr_unsigned)->strides = NULL; if (rstr->strides) { CeedCall(CeedMalloc(3, &(*rstr_unsigned)->strides)); for (CeedInt i = 0; i < 3; i++) (*rstr_unsigned)->strides[i] = rstr->strides[i]; } CeedCall(CeedElemRestrictionReferenceCopy(rstr, &(*rstr_unsigned)->rstr_base)); // Override Apply (*rstr_unsigned)->Apply = rstr->ApplyUnsigned; return CEED_ERROR_SUCCESS; } /** @brief Copy the pointer to a CeedElemRestriction and set `CeedElemRestrictionApply()` implementation to use the unoriented version. Both pointers should be destroyed with `CeedElemRestrictionDestroy()`. @param[in] rstr CeedElemRestriction to create unoriented reference to @param[in,out] rstr_unoriented Variable to store unoriented CeedElemRestriction @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionCreateUnorientedCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_unoriented) { CeedCall(CeedCalloc(1, rstr_unoriented)); // Copy old rstr memcpy(*rstr_unoriented, rstr, sizeof(struct CeedElemRestriction_private)); (*rstr_unoriented)->ceed = NULL; CeedCall(CeedReferenceCopy(rstr->ceed, &(*rstr_unoriented)->ceed)); (*rstr_unoriented)->ref_count = 1; (*rstr_unoriented)->strides = NULL; if (rstr->strides) { CeedCall(CeedMalloc(3, &(*rstr_unoriented)->strides)); for (CeedInt i = 0; i < 3; i++) (*rstr_unoriented)->strides[i] = rstr->strides[i]; } CeedCall(CeedElemRestrictionReferenceCopy(rstr, &(*rstr_unoriented)->rstr_base)); // Override Apply (*rstr_unoriented)->Apply = rstr->ApplyUnoriented; return CEED_ERROR_SUCCESS; } /** @brief Copy the pointer to a CeedElemRestriction. Both pointers should be destroyed with `CeedElemRestrictionDestroy()`. 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. This CeedElemRestriction will be destroyed if `rstr_copy` is the only reference to this CeedElemRestriction. @param[in] rstr CeedElemRestriction to copy reference to @param[in,out] rstr_copy Variable to store copied reference @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionReferenceCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_copy) { if (rstr != CEED_ELEMRESTRICTION_NONE) CeedCall(CeedElemRestrictionReference(rstr)); CeedCall(CeedElemRestrictionDestroy(rstr_copy)); *rstr_copy = rstr; return CEED_ERROR_SUCCESS; } /** @brief Create CeedVectors associated with a CeedElemRestriction @param[in] rstr CeedElemRestriction @param[out] l_vec The address of the L-vector to be created, or NULL @param[out] e_vec The address of the E-vector to be created, or NULL @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionCreateVector(CeedElemRestriction rstr, CeedVector *l_vec, CeedVector *e_vec) { CeedSize e_size, l_size; l_size = rstr->l_size; e_size = rstr->num_block * rstr->block_size * rstr->elem_size * rstr->num_comp; if (l_vec) CeedCall(CeedVectorCreate(rstr->ceed, l_size, l_vec)); if (e_vec) CeedCall(CeedVectorCreate(rstr->ceed, e_size, e_vec)); return CEED_ERROR_SUCCESS; } /** @brief Restrict an L-vector to an E-vector or apply its transpose @param[in] rstr CeedElemRestriction @param[in] t_mode Apply restriction or transpose @param[in] u Input vector (of size @a l_size when t_mode=@ref CEED_NOTRANSPOSE) @param[out] ru Output vector (of shape [@a num_elem * @a elem_size] when t_mode=@ref CEED_NOTRANSPOSE). Ordering of the e-vector is decided by the backend. @param[in] request Request or @ref CEED_REQUEST_IMMEDIATE @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, CeedRequest *request) { CeedInt m, n; if (t_mode == CEED_NOTRANSPOSE) { m = rstr->e_size; n = rstr->l_size; } else { m = rstr->l_size; n = rstr->e_size; } CeedCheck(n <= u->length, rstr->ceed, CEED_ERROR_DIMENSION, "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", u->length, m, n); CeedCheck(m <= ru->length, rstr->ceed, CEED_ERROR_DIMENSION, "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", ru->length, m, n); if (rstr->num_elem > 0) CeedCall(rstr->Apply(rstr, t_mode, u, ru, request)); return CEED_ERROR_SUCCESS; } /** @brief Restrict an L-vector of points to a single element or apply its transpose @param[in] rstr CeedElemRestriction @param[in] elem Element number in range 0..@a num_elem @param[in] t_mode Apply restriction or transpose @param[in] u Input vector (of size @a l_size when t_mode=@ref CEED_NOTRANSPOSE) @param[out] ru Output vector (of shape [@a num_elem * @a elem_size] when t_mode=@ref CEED_NOTRANSPOSE). Ordering of the e-vector is decided by the backend. @param[in] request Request or @ref CEED_REQUEST_IMMEDIATE @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionApplyAtPointsInElement(CeedElemRestriction rstr, CeedInt elem, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, CeedRequest *request) { CeedInt m, n; if (t_mode == CEED_NOTRANSPOSE) { CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr, elem, &m)); n = rstr->l_size; } else { m = rstr->l_size; CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr, elem, &n)); } CeedCheck(n <= u->length, rstr->ceed, CEED_ERROR_DIMENSION, "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ") for element %" CeedInt_FMT, u->length, m, n, elem); CeedCheck(m <= ru->length, rstr->ceed, CEED_ERROR_DIMENSION, "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ") for element %" CeedInt_FMT, ru->length, m, n, elem); CeedCheck(elem < rstr->num_elem, rstr->ceed, CEED_ERROR_DIMENSION, "Cannot retrieve element %" CeedInt_FMT ", element %" CeedInt_FMT " > total elements %" CeedInt_FMT "", elem, elem, rstr->num_elem); if (rstr->num_elem > 0) CeedCall(rstr->ApplyAtPointsInElement(rstr, elem, t_mode, u, ru, request)); return CEED_ERROR_SUCCESS; } /** @brief Restrict an L-vector to a block of an E-vector or apply its transpose @param[in] rstr CeedElemRestriction @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 [3*block_size : 4*block_size] @param[in] t_mode Apply restriction or transpose @param[in] u Input vector (of size @a l_size when t_mode=@ref CEED_NOTRANSPOSE) @param[out] ru Output vector (of shape [@a block_size * @a elem_size] when t_mode=@ref CEED_NOTRANSPOSE). Ordering of the e-vector is decided by the backend. @param[in] request Request or @ref CEED_REQUEST_IMMEDIATE @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedElemRestrictionApplyBlock(CeedElemRestriction rstr, CeedInt block, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, CeedRequest *request) { CeedInt m, n; CeedCheck(rstr->ApplyBlock, rstr->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionApplyBlock"); if (t_mode == CEED_NOTRANSPOSE) { m = rstr->block_size * rstr->elem_size * rstr->num_comp; n = rstr->l_size; } else { m = rstr->l_size; n = rstr->block_size * rstr->elem_size * rstr->num_comp; } CeedCheck(n == u->length, rstr->ceed, CEED_ERROR_DIMENSION, "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", u->length, m, n); CeedCheck(m == ru->length, rstr->ceed, CEED_ERROR_DIMENSION, "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", ru->length, m, n); CeedCheck(rstr->block_size * block <= rstr->num_elem, rstr->ceed, CEED_ERROR_DIMENSION, "Cannot retrieve block %" CeedInt_FMT ", element %" CeedInt_FMT " > total elements %" CeedInt_FMT "", block, rstr->block_size * block, rstr->num_elem); CeedCall(rstr->ApplyBlock(rstr, block, t_mode, u, ru, request)); return CEED_ERROR_SUCCESS; } /** @brief Get the Ceed associated with a CeedElemRestriction @param[in] rstr CeedElemRestriction @param[out] ceed Variable to store Ceed @return An error code: 0 - success, otherwise - failure @ref Advanced **/ int CeedElemRestrictionGetCeed(CeedElemRestriction rstr, Ceed *ceed) { *ceed = rstr->ceed; return CEED_ERROR_SUCCESS; } /** @brief Get the L-vector component stride @param[in] rstr CeedElemRestriction @param[out] comp_stride Variable to store component stride @return An error code: 0 - success, otherwise - failure @ref Advanced **/ int CeedElemRestrictionGetCompStride(CeedElemRestriction rstr, CeedInt *comp_stride) { *comp_stride = rstr->comp_stride; return CEED_ERROR_SUCCESS; } /** @brief Get the total number of elements in the range of a CeedElemRestriction @param[in] rstr CeedElemRestriction @param[out] num_elem Variable to store number of elements @return An error code: 0 - success, otherwise - failure @ref Advanced **/ int CeedElemRestrictionGetNumElements(CeedElemRestriction rstr, CeedInt *num_elem) { *num_elem = rstr->num_elem; return CEED_ERROR_SUCCESS; } /** @brief Get the size of elements in the CeedElemRestriction @param[in] rstr CeedElemRestriction @param[out] elem_size Variable to store size of elements @return An error code: 0 - success, otherwise - failure @ref Advanced **/ int CeedElemRestrictionGetElementSize(CeedElemRestriction rstr, CeedInt *elem_size) { *elem_size = rstr->elem_size; return CEED_ERROR_SUCCESS; } /** @brief Get the number of points in the l-vector for a points CeedElemRestriction @param[in] rstr CeedElemRestriction @param[out] num_points The number of points in the l-vector @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionGetNumPoints(CeedElemRestriction rstr, CeedInt *num_points) { Ceed ceed; CeedCall(CeedElemRestrictionGetCeed(rstr, &ceed)); CeedCheck(rstr->rstr_type == CEED_RESTRICTION_POINTS, ceed, CEED_ERROR_INCOMPATIBLE, "Can only retrieve the number of points for a points CeedElemRestriction"); *num_points = rstr->num_points; return CEED_ERROR_SUCCESS; } /** @brief Get the number of points in an element of a points CeedElemRestriction @param[in] rstr CeedElemRestriction @param[in] elem Index number of element to retrieve the number of points for @param[out] num_points The number of points in the element at index elem @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionGetNumPointsInElement(CeedElemRestriction rstr, CeedInt elem, CeedInt *num_points) { Ceed ceed; const CeedInt *offsets; CeedCall(CeedElemRestrictionGetCeed(rstr, &ceed)); CeedCheck(rstr->rstr_type == CEED_RESTRICTION_POINTS, ceed, CEED_ERROR_INCOMPATIBLE, "Can only retrieve the number of points for a points CeedElemRestriction"); CeedCall(CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets)); *num_points = offsets[elem + 1] - offsets[elem]; CeedCall(CeedElemRestrictionRestoreOffsets(rstr, &offsets)); return CEED_ERROR_SUCCESS; } /** @brief Get the maximum number of points in an element for a CeedElemRestriction at points @param[in] rstr CeedElemRestriction @param[out] max_points Variable to store size of elements @return An error code: 0 - success, otherwise - failure @ref Advanced **/ int CeedElemRestrictionGetMaxPointsInElement(CeedElemRestriction rstr, CeedInt *max_points) { Ceed ceed; CeedInt num_elem; CeedRestrictionType rstr_type; CeedCall(CeedElemRestrictionGetCeed(rstr, &ceed)); CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); CeedCheck(rstr_type == CEED_RESTRICTION_POINTS, ceed, CEED_ERROR_INCOMPATIBLE, "Cannot compute max points for a CeedElemRestriction that does not use points"); CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); *max_points = 0; for (CeedInt e = 0; e < num_elem; e++) { CeedInt num_points; CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr, e, &num_points)); *max_points = CeedIntMax(num_points, *max_points); } return CEED_ERROR_SUCCESS; } /** @brief Get the size of the l-vector for a CeedElemRestriction @param[in] rstr CeedElemRestriction @param[out] l_size Variable to store number of nodes @return An error code: 0 - success, otherwise - failure @ref Advanced **/ int CeedElemRestrictionGetLVectorSize(CeedElemRestriction rstr, CeedSize *l_size) { *l_size = rstr->l_size; return CEED_ERROR_SUCCESS; } /** @brief Get the number of components in the elements of a CeedElemRestriction @param[in] rstr CeedElemRestriction @param[out] num_comp Variable to store number of components @return An error code: 0 - success, otherwise - failure @ref Advanced **/ int CeedElemRestrictionGetNumComponents(CeedElemRestriction rstr, CeedInt *num_comp) { *num_comp = rstr->num_comp; return CEED_ERROR_SUCCESS; } /** @brief Get the number of blocks in a CeedElemRestriction @param[in] rstr CeedElemRestriction @param[out] num_block Variable to store number of blocks @return An error code: 0 - success, otherwise - failure @ref Advanced **/ int CeedElemRestrictionGetNumBlocks(CeedElemRestriction rstr, CeedInt *num_block) { *num_block = rstr->num_block; return CEED_ERROR_SUCCESS; } /** @brief Get the size of blocks in the CeedElemRestriction @param[in] rstr CeedElemRestriction @param[out] block_size Variable to store size of blocks @return An error code: 0 - success, otherwise - failure @ref Advanced **/ int CeedElemRestrictionGetBlockSize(CeedElemRestriction rstr, CeedInt *block_size) { *block_size = rstr->block_size; return CEED_ERROR_SUCCESS; } /** @brief Get the multiplicity of nodes in a CeedElemRestriction @param[in] rstr CeedElemRestriction @param[out] mult Vector to store multiplicity (of size l_size) @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionGetMultiplicity(CeedElemRestriction rstr, CeedVector mult) { CeedVector e_vec; // Create e_vec to hold intermediate computation in E^T (E 1) CeedCall(CeedElemRestrictionCreateVector(rstr, NULL, &e_vec)); // Compute e_vec = E * 1 CeedCall(CeedVectorSetValue(mult, 1.0)); CeedCall(CeedElemRestrictionApply(rstr, CEED_NOTRANSPOSE, mult, e_vec, CEED_REQUEST_IMMEDIATE)); // Compute multiplicity, mult = E^T * e_vec = E^T (E 1) CeedCall(CeedVectorSetValue(mult, 0.0)); CeedCall(CeedElemRestrictionApply(rstr, CEED_TRANSPOSE, e_vec, mult, CEED_REQUEST_IMMEDIATE)); // Cleanup CeedCall(CeedVectorDestroy(&e_vec)); return CEED_ERROR_SUCCESS; } /** @brief View a CeedElemRestriction @param[in] rstr CeedElemRestriction to view @param[in] stream Stream to write; typically stdout/stderr or a file @return Error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionView(CeedElemRestriction rstr, FILE *stream) { CeedRestrictionType rstr_type; CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); if (rstr_type == CEED_RESTRICTION_POINTS) { CeedInt max_points; CeedCall(CeedElemRestrictionGetMaxPointsInElement(rstr, &max_points)); fprintf(stream, "CeedElemRestriction at points from (%td, %" CeedInt_FMT ") to %" CeedInt_FMT " elements with a maximum of %" CeedInt_FMT " points on an element\n", rstr->l_size, rstr->num_comp, rstr->num_elem, max_points); } else { char stridesstr[500]; if (rstr->strides) { sprintf(stridesstr, "[%" CeedInt_FMT ", %" CeedInt_FMT ", %" CeedInt_FMT "]", rstr->strides[0], rstr->strides[1], rstr->strides[2]); } else { sprintf(stridesstr, "%" CeedInt_FMT, rstr->comp_stride); } fprintf(stream, "%sCeedElemRestriction from (%td, %" CeedInt_FMT ") to %" CeedInt_FMT " elements with %" CeedInt_FMT " nodes each and %s %s\n", rstr->block_size > 1 ? "Blocked " : "", rstr->l_size, rstr->num_comp, rstr->num_elem, rstr->elem_size, rstr->strides ? "strides" : "component stride", stridesstr); } return CEED_ERROR_SUCCESS; } /** @brief Destroy a CeedElemRestriction @param[in,out] rstr CeedElemRestriction to destroy @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedElemRestrictionDestroy(CeedElemRestriction *rstr) { if (!*rstr || *rstr == CEED_ELEMRESTRICTION_NONE || --(*rstr)->ref_count > 0) { *rstr = NULL; return CEED_ERROR_SUCCESS; } CeedCheck((*rstr)->num_readers == 0, (*rstr)->ceed, CEED_ERROR_ACCESS, "Cannot destroy CeedElemRestriction, a process has read access to the offset data"); // Only destroy backend data once between rstr and unsigned copy if ((*rstr)->rstr_base) CeedCall(CeedElemRestrictionDestroy(&(*rstr)->rstr_base)); else if ((*rstr)->Destroy) CeedCall((*rstr)->Destroy(*rstr)); CeedCall(CeedFree(&(*rstr)->strides)); CeedCall(CeedDestroy(&(*rstr)->ceed)); CeedCall(CeedFree(rstr)); return CEED_ERROR_SUCCESS; } /// @}