Add CeedElemRestrictionCreateUnsignedCopy

include/ceed-impl.h

@@ -149,7 +149,8 @@ struct CeedVector_private {
 };
 
 struct CeedElemRestriction_private {
-  Ceed ceed;
+  Ceed                ceed;
+  CeedElemRestriction rstr_signed;
   int (*Apply)(CeedElemRestriction, CeedTransposeMode, CeedVector, CeedVector, CeedRequest *);
   int (*ApplyUnsigned)(CeedElemRestriction, CeedTransposeMode, CeedVector, CeedVector, CeedRequest *);
   int (*ApplyBlock)(CeedElemRestriction, CeedInt, CeedTransposeMode, CeedVector, CeedVector, CeedRequest *);

include/ceed/backend.h

@@ -164,6 +164,8 @@ CEED_EXTERN int CeedVectorGetData(CeedVector vec, void *data);
 CEED_EXTERN int CeedVectorSetData(CeedVector vec, void *data);
 CEED_EXTERN int CeedVectorReference(CeedVector vec);
 
+CEED_EXTERN int CeedElemRestrictionApplyUnsigned(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector ru,
+                                                 CeedRequest *request);
 CEED_EXTERN int CeedElemRestrictionGetStrides(CeedElemRestriction rstr, CeedInt (*strides)[3]);
 CEED_EXTERN int CeedElemRestrictionGetOffsets(CeedElemRestriction rstr, CeedMemType mem_type, const CeedInt **offsets);
 CEED_EXTERN int CeedElemRestrictionRestoreOffsets(CeedElemRestriction rstr, const CeedInt **offsets);

include/ceed/ceed.h

@@ -240,11 +240,10 @@ CEED_EXTERN int CeedElemRestrictionCreateBlocked(Ceed ceed, CeedInt num_elem, Ce
                                                  const CeedInt *offsets, CeedElemRestriction *rstr);
 CEED_EXTERN int CeedElemRestrictionCreateBlockedStrided(Ceed ceed, CeedInt num_elem, CeedInt elem_size, CeedInt blk_size, CeedInt num_comp,
                                                         CeedSize l_size, const CeedInt strides[3], CeedElemRestriction *rstr);
+CEED_EXTERN int CeedElemRestrictionCreateUnsignedCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_unsigned);
 CEED_EXTERN int CeedElemRestrictionReferenceCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_copy);
 CEED_EXTERN int CeedElemRestrictionCreateVector(CeedElemRestriction rstr, CeedVector *lvec, CeedVector *evec);
 CEED_EXTERN int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, CeedRequest *request);
-CEED_EXTERN int CeedElemRestrictionApplyUnsigned(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector ru,
-                                                 CeedRequest *request);
 CEED_EXTERN int CeedElemRestrictionApplyBlock(CeedElemRestriction rstr, CeedInt block, CeedTransposeMode t_mode, CeedVector u, CeedVector ru,
                                               CeedRequest *request);
 CEED_EXTERN int CeedElemRestrictionGetCeed(CeedElemRestriction rstr, Ceed *ceed);

interface/ceed-elemrestriction.c

@@ -10,6 +10,7 @@
 #include <ceed/backend.h>
 #include <stdbool.h>
 #include <stdio.h>
+#include <string.h>
 
 /// @file
 /// Implementation of CeedElemRestriction interfaces
@@ -55,6 +56,41 @@ int CeedPermutePadOffsets(const CeedInt *offsets, CeedInt *blk_offsets, CeedInt
 /// @addtogroup CeedElemRestrictionBackend
 /// @{
 
+/**
+  @brief Restrict an L-vector to an E-vector or apply its transpose ignoring any
+         provided orientations
+
+  @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 CeedElemRestrictionApplyUnsigned(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, CeedRequest *request) {
+  CeedInt m, n;
+
+  CeedCheck(rstr->ApplyUnsigned, rstr->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support ElemRestrictionApplyUnsigned");
+
+  if (t_mode == CEED_NOTRANSPOSE) {
+    m = rstr->num_blk * rstr->blk_size * rstr->elem_size * rstr->num_comp;
+    n = rstr->l_size;
+  } else {
+    m = rstr->l_size;
+    n = rstr->num_blk * rstr->blk_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);
+  if (rstr->num_elem > 0) CeedCall(rstr->ApplyUnsigned(rstr, t_mode, u, ru, request));
+  return CEED_ERROR_SUCCESS;
+}
+
 /**
 
   @brief Get the strides of a strided CeedElemRestriction
@@ -85,9 +121,13 @@ int CeedElemRestrictionGetStrides(CeedElemRestriction rstr, CeedInt (*strides)[3
   @ref User
 **/
 int CeedElemRestrictionGetOffsets(CeedElemRestriction rstr, CeedMemType mem_type, const CeedInt **offsets) {
-  CeedCheck(rstr->GetOffsets, rstr->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support GetOffsets");
-  CeedCall(rstr->GetOffsets(rstr, mem_type, offsets));
-  rstr->num_readers++;
+  if (rstr->rstr_signed) {
+    CeedCall(CeedElemRestrictionGetOffsets(rstr->rstr_signed, 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;
 }
 
@@ -102,8 +142,12 @@ int CeedElemRestrictionGetOffsets(CeedElemRestriction rstr, CeedMemType mem_type
   @ref User
 **/
 int CeedElemRestrictionRestoreOffsets(CeedElemRestriction rstr, const CeedInt **offsets) {
-  *offsets = NULL;
-  rstr->num_readers--;
+  if (rstr->rstr_signed) {
+    CeedCall(CeedElemRestrictionRestoreOffsets(rstr->rstr_signed, offsets));
+  } else {
+    *offsets = NULL;
+    rstr->num_readers--;
+  }
   return CEED_ERROR_SUCCESS;
 }
 
@@ -558,6 +602,39 @@ int CeedElemRestrictionCreateBlockedStrided(Ceed ceed, CeedInt num_elem, CeedInt
   return CEED_ERROR_SUCCESS;
 }
 
+/**
+  @brief Copy the pointer to a CeedElemRestriction and set `CeedElemRestrictionApply()` implementation to `CeedElemRestrictionApplyUnsigned()`.
+
+  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_signed));
+
+  // Override Apply
+  (*rstr_unsigned)->Apply = rstr->ApplyUnsigned;
+
+  return CEED_ERROR_SUCCESS;
+}
+
 /**
   @brief Copy the pointer to a CeedElemRestriction.
 
@@ -632,41 +709,6 @@ int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode t_mode,
   return CEED_ERROR_SUCCESS;
 }
 
-/**
-  @brief Restrict an L-vector to an E-vector or apply its transpose ignoring any
-         provided orientations
-
-  @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 CeedElemRestrictionApplyUnsigned(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, CeedRequest *request) {
-  CeedInt m, n;
-
-  CeedCheck(rstr->ApplyUnsigned, rstr->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support ElemRestrictionApplyUnsigned");
-
-  if (t_mode == CEED_NOTRANSPOSE) {
-    m = rstr->num_blk * rstr->blk_size * rstr->elem_size * rstr->num_comp;
-    n = rstr->l_size;
-  } else {
-    m = rstr->l_size;
-    n = rstr->num_blk * rstr->blk_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);
-  if (rstr->num_elem > 0) CeedCall(rstr->ApplyUnsigned(rstr, 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
 
@@ -894,7 +936,11 @@ int CeedElemRestrictionDestroy(CeedElemRestriction *rstr) {
   }
   CeedCheck((*rstr)->num_readers == 0, (*rstr)->ceed, CEED_ERROR_ACCESS,
             "Cannot destroy CeedElemRestriction, a process has read access to the offset data");
-  if ((*rstr)->Destroy) CeedCall((*rstr)->Destroy(*rstr));
+
+  // Only destroy backend data once between rstr and unsigned copy
+  if ((*rstr)->rstr_signed) CeedCall(CeedElemRestrictionDestroy(&(*rstr)->rstr_signed));
+  else if ((*rstr)->Destroy) CeedCall((*rstr)->Destroy(*rstr));
+
   CeedCall(CeedFree(&(*rstr)->strides));
   CeedCall(CeedDestroy(&(*rstr)->ceed));
   CeedCall(CeedFree(rstr));

interface/ceed-preconditioning.c

@@ -702,8 +702,9 @@ static int CeedSingleOperatorAssemblyCountEntries(CeedOperator op, CeedInt *num_
 **/
 static int CeedSingleOperatorMultigridLevel(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse,
                                             CeedBasis basis_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong, CeedOperator *op_restrict) {
-  Ceed       ceed;
-  CeedVector mult_vec = NULL;
+  Ceed                ceed;
+  CeedVector          mult_vec         = NULL;
+  CeedElemRestriction rstr_p_mult_fine = NULL;
   CeedCall(CeedOperatorGetCeed(op_fine, &ceed));
 
   // Check for composite operator
@@ -740,12 +741,13 @@ static int CeedSingleOperatorMultigridLevel(CeedOperator op_fine, CeedVector p_m
   if (op_restrict || op_prolong) {
     CeedVector mult_e_vec;
 
+    CeedCall(CeedElemRestrictionCreateUnsignedCopy(rstr_fine, &rstr_p_mult_fine));
     CeedCheck(p_mult_fine, ceed, CEED_ERROR_INCOMPATIBLE, "Prolongation or restriction operator creation requires fine grid multiplicity vector");
     CeedCall(CeedElemRestrictionCreateVector(rstr_fine, &mult_vec, &mult_e_vec));
     CeedCall(CeedVectorSetValue(mult_e_vec, 0.0));
-    CeedCall(CeedElemRestrictionApply(rstr_fine, CEED_NOTRANSPOSE, p_mult_fine, mult_e_vec, CEED_REQUEST_IMMEDIATE));
+    CeedCall(CeedElemRestrictionApply(rstr_p_mult_fine, CEED_NOTRANSPOSE, p_mult_fine, mult_e_vec, CEED_REQUEST_IMMEDIATE));
     CeedCall(CeedVectorSetValue(mult_vec, 0.0));
-    CeedCall(CeedElemRestrictionApply(rstr_fine, CEED_TRANSPOSE, mult_e_vec, mult_vec, CEED_REQUEST_IMMEDIATE));
+    CeedCall(CeedElemRestrictionApply(rstr_p_mult_fine, CEED_TRANSPOSE, mult_e_vec, mult_vec, CEED_REQUEST_IMMEDIATE));
     CeedCall(CeedVectorDestroy(&mult_e_vec));
     CeedCall(CeedVectorReciprocal(mult_vec));
   }
@@ -783,7 +785,7 @@ static int CeedSingleOperatorMultigridLevel(CeedOperator op_fine, CeedVector p_m
 
     CeedCall(CeedOperatorCreate(ceed, qf_restrict, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, op_restrict));
     CeedCall(CeedOperatorSetField(*op_restrict, "input", rstr_fine, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE));
-    CeedCall(CeedOperatorSetField(*op_restrict, "scale", rstr_fine, CEED_BASIS_COLLOCATED, mult_vec));
+    CeedCall(CeedOperatorSetField(*op_restrict, "scale", rstr_p_mult_fine, CEED_BASIS_COLLOCATED, mult_vec));
     CeedCall(CeedOperatorSetField(*op_restrict, "output", rstr_coarse, basis_c_to_f, CEED_VECTOR_ACTIVE));
 
     // Set name
@@ -820,7 +822,7 @@ static int CeedSingleOperatorMultigridLevel(CeedOperator op_fine, CeedVector p_m
 
     CeedCall(CeedOperatorCreate(ceed, qf_prolong, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, op_prolong));
     CeedCall(CeedOperatorSetField(*op_prolong, "input", rstr_coarse, basis_c_to_f, CEED_VECTOR_ACTIVE));
-    CeedCall(CeedOperatorSetField(*op_prolong, "scale", rstr_fine, CEED_BASIS_COLLOCATED, mult_vec));
+    CeedCall(CeedOperatorSetField(*op_prolong, "scale", rstr_p_mult_fine, CEED_BASIS_COLLOCATED, mult_vec));
     CeedCall(CeedOperatorSetField(*op_prolong, "output", rstr_fine, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE));
 
     // Set name
@@ -842,6 +844,7 @@ static int CeedSingleOperatorMultigridLevel(CeedOperator op_fine, CeedVector p_m
 
   // Cleanup
   CeedCall(CeedVectorDestroy(&mult_vec));
+  CeedCall(CeedElemRestrictionDestroy(&rstr_p_mult_fine));
   CeedCall(CeedBasisDestroy(&basis_c_to_f));
 
   return CEED_ERROR_SUCCESS;

tests/t221-elemrestriction.c

@@ -0,0 +1,81 @@
+/// @file
+/// Test creation, use, and destruction of an element restriction oriented with unsigned application
+/// \test Test creation, use, and destruction of an element restriction oriented with unsigned application
+#include <ceed.h>
+#include <ceed/backend.h>
+#include <math.h>
+#include <stdio.h>
+
+int main(int argc, char **argv) {
+  Ceed                ceed;
+  CeedVector          x, y_oriented, y_unsigned, y_unsigned_copy;
+  CeedInt             num_elem = 6, p = 2, dim = 1;
+  CeedInt             ind[p * num_elem];
+  bool                orient[p * num_elem];
+  CeedScalar          x_array[num_elem + 1];
+  CeedElemRestriction elem_restriction, elem_restriction_copy;
+
+  CeedInit(argv[1], &ceed);
+
+  CeedVectorCreate(ceed, num_elem + 1, &x);
+  for (CeedInt i = 0; i < num_elem + 1; i++) x_array[i] = 10 + i;
+  CeedVectorSetArray(x, CEED_MEM_HOST, CEED_USE_POINTER, x_array);
+  CeedVectorCreate(ceed, num_elem * 2, &y_oriented);
+  CeedVectorCreate(ceed, num_elem * 2, &y_unsigned);
+  CeedVectorCreate(ceed, num_elem * 2, &y_unsigned_copy);
+
+  for (CeedInt i = 0; i < num_elem; i++) {
+    ind[2 * i + 0] = i;
+    ind[2 * i + 1] = i + 1;
+    // flip the dofs on element 1,3,...
+    orient[2 * i + 0] = (i % (2)) * -1 < 0;
+    orient[2 * i + 1] = (i % (2)) * -1 < 0;
+  }
+  CeedElemRestrictionCreateOriented(ceed, num_elem, p, dim, 1, num_elem + 1, CEED_MEM_HOST, CEED_USE_POINTER, ind, orient, &elem_restriction);
+  CeedElemRestrictionCreateUnsignedCopy(elem_restriction, &elem_restriction_copy);
+
+  CeedElemRestrictionApply(elem_restriction, CEED_NOTRANSPOSE, x, y_oriented, CEED_REQUEST_IMMEDIATE);
+  CeedElemRestrictionApplyUnsigned(elem_restriction, CEED_NOTRANSPOSE, x, y_unsigned, CEED_REQUEST_IMMEDIATE);
+  CeedElemRestrictionApply(elem_restriction_copy, CEED_NOTRANSPOSE, x, y_unsigned_copy, CEED_REQUEST_IMMEDIATE);
+
+  {
+    const CeedScalar *y_oriented_array, *y_unsigned_array, *y_unsigned_copy_array;
+
+    CeedVectorGetArrayRead(y_oriented, CEED_MEM_HOST, &y_oriented_array);
+    CeedVectorGetArrayRead(y_unsigned, CEED_MEM_HOST, &y_unsigned_array);
+    CeedVectorGetArrayRead(y_unsigned_copy, CEED_MEM_HOST, &y_unsigned_copy_array);
+    for (CeedInt i = 0; i < num_elem; i++) {
+      for (CeedInt j = 0; j < p; j++) {
+        CeedInt k = j + p * i;
+        // unsigned application should match oriented application, but with no sign change
+        if (y_oriented_array[k] * CeedIntPow(-1, i % 2) != 10 + (k + 1) / 2) {
+          // LCOV_EXCL_START
+          printf("Error in oriented restricted array y[%" CeedInt_FMT "] = %f\n", k, (CeedScalar)y_oriented_array[k]);
+          // LCOV_EXCL_STOP
+        }
+        if (y_unsigned_array[k] != 10 + (k + 1) / 2) {
+          // LCOV_EXCL_START
+          printf("Error in unsigned restricted array y[%" CeedInt_FMT "] = %f\n", k, (CeedScalar)y_unsigned_array[k]);
+          // LCOV_EXCL_STOP
+        }
+        if (y_unsigned_array[k] != y_unsigned_copy_array[k]) {
+          // LCOV_EXCL_START
+          printf("Error in copy restricted array y[%" CeedInt_FMT "] = %f\n", k, (CeedScalar)y_unsigned_copy_array[k]);
+          // LCOV_EXCL_STOP
+        }
+      }
+    }
+    CeedVectorRestoreArrayRead(y_oriented, &y_oriented_array);
+    CeedVectorRestoreArrayRead(y_unsigned, &y_unsigned_array);
+    CeedVectorRestoreArrayRead(y_unsigned_copy, &y_unsigned_copy_array);
+  }
+
+  CeedVectorDestroy(&x);
+  CeedVectorDestroy(&y_oriented);
+  CeedVectorDestroy(&y_unsigned);
+  CeedVectorDestroy(&y_unsigned_copy);
+  CeedElemRestrictionDestroy(&elem_restriction);
+  CeedElemRestrictionDestroy(&elem_restriction_copy);
+  CeedDestroy(&ceed);
+  return 0;
+}

tests/t552-operator.c

@@ -164,6 +164,8 @@ int main(int argc, char **argv) {
   CeedElemRestrictionDestroy(&elem_restriction_q_data);
   CeedBasisDestroy(&basis_u);
   CeedBasisDestroy(&basis_x);
+  CeedBasisDestroy(&basis_u_coarse);
+  CeedBasisDestroy(&basis_coarse_to_fine);
   CeedQFunctionDestroy(&qf_setup);
   CeedQFunctionDestroy(&qf_mass);
   CeedOperatorDestroy(&op_setup);