Merge pull request #889 from CEED/rezgar/oriented-restr

jedbrown · web-flow · commit c6e1a279e1a1 · 2022-02-07T06:47:11.000-07:00
Element Restriction Oriented
diff --git a/backends/ref/ceed-ref-restriction.c b/backends/ref/ceed-ref-restriction.c
@@ -37,6 +37,8 @@ static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r,
   ierr = CeedElemRestrictionGetElementSize(r, &elem_size); CeedChkBackend(ierr);
   v_offset = start*blk_size*elem_size*num_comp;
 
+  bool is_oriented;
+  ierr = CeedElemRestrictionIsOriented(r, &is_oriented); CeedChkBackend(ierr);
   ierr = CeedVectorGetArrayRead(u, CEED_MEM_HOST, &uu); CeedChkBackend(ierr);
   if (t_mode == CEED_TRANSPOSE) {
     // Sum into for transpose mode, e-vec to l-vec
@@ -91,7 +93,8 @@ static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r,
           CeedPragmaSIMD
           for (CeedInt i = 0; i < elem_size*blk_size; i++)
             vv[elem_size*(k*blk_size+num_comp*e) + i - v_offset]
-              = uu[impl->offsets[i+elem_size*e] + k*comp_stride];
+              = uu[impl->offsets[i+elem_size*e] + k*comp_stride] *
+                (is_oriented && impl->orient[i+elem_size*e] ? -1. : 1.);
     }
   } else {
     // Restriction from E-vector to L-vector
@@ -137,7 +140,8 @@ static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r,
             // Iteration bound set to discard padding elements
             for (CeedInt j = i; j < i+CeedIntMin(blk_size, num_elem-e); j++)
               vv[impl->offsets[j+e*elem_size] + k*comp_stride]
-              += uu[elem_size*(k*blk_size+num_comp*e) + j - v_offset];
+              += uu[elem_size*(k*blk_size+num_comp*e) + j - v_offset] *
+                 (is_oriented && impl->orient[j+e*elem_size] ? -1. : 1.);
     }
   }
   ierr = CeedVectorRestoreArrayRead(u, &uu); CeedChkBackend(ierr);
@@ -456,4 +460,40 @@ int CeedElemRestrictionCreate_Ref(CeedMemType mem_type, CeedCopyMode copy_mode,
 
   return CEED_ERROR_SUCCESS;
 }
+
+//------------------------------------------------------------------------------
+// ElemRestriction Create Oriented
+//------------------------------------------------------------------------------
+int CeedElemRestrictionCreateOriented_Ref(CeedMemType mem_type,
+    CeedCopyMode copy_mode,
+    const CeedInt *offsets, const bool *orient,
+    CeedElemRestriction r) {
+  int ierr;
+  CeedElemRestriction_Ref *impl;
+  CeedInt num_elem, elem_size;
+  // Set up for normal restriction with explicit offsets. This sets up dispatch to
+  // CeedElemRestrictionApply_Ref_* and manages the impl->offsets array copy/allocation.
+  ierr = CeedElemRestrictionCreate_Ref(mem_type, copy_mode, offsets, r);
+  CeedChkBackend(ierr);
+
+  ierr = CeedElemRestrictionGetData(r, &impl); CeedChkBackend(ierr);
+  ierr = CeedElemRestrictionGetNumElements(r, &num_elem); CeedChkBackend(ierr);
+  ierr = CeedElemRestrictionGetElementSize(r, &elem_size); CeedChkBackend(ierr);
+  switch (copy_mode) {
+  case CEED_COPY_VALUES:
+    ierr = CeedMalloc(num_elem * elem_size, &impl->orient_allocated);
+    CeedChkBackend(ierr);
+    memcpy(impl->orient_allocated, orient,
+           num_elem * elem_size * sizeof(orient[0]));
+    impl->orient = impl->orient_allocated;
+    break;
+  case CEED_OWN_POINTER:
+    impl->orient_allocated = (bool *)orient;
+    impl->orient = impl->orient_allocated;
+    break;
+  case CEED_USE_POINTER:
+    impl->orient = orient;
+  }
+  return CEED_ERROR_SUCCESS;
+}
 //------------------------------------------------------------------------------
diff --git a/backends/ref/ceed-ref.c b/backends/ref/ceed-ref.c
@@ -44,6 +44,9 @@ static int CeedInit_Ref(const char *resource, Ceed ceed) {
                                 CeedTensorContractCreate_Ref); CeedChkBackend(ierr);
   ierr = CeedSetBackendFunction(ceed, "Ceed", ceed, "ElemRestrictionCreate",
                                 CeedElemRestrictionCreate_Ref); CeedChkBackend(ierr);
+  ierr = CeedSetBackendFunction(ceed, "Ceed", ceed,
+                                "ElemRestrictionCreateOriented",
+                                CeedElemRestrictionCreateOriented_Ref); CeedChkBackend(ierr);
   ierr = CeedSetBackendFunction(ceed, "Ceed", ceed,
                                 "ElemRestrictionCreateBlocked",
                                 CeedElemRestrictionCreate_Ref); CeedChkBackend(ierr);
diff --git a/backends/ref/ceed-ref.h b/backends/ref/ceed-ref.h
@@ -36,6 +36,9 @@ typedef struct {
 typedef struct {
   const CeedInt *offsets;
   CeedInt *offsets_allocated;
+  // Orientation, if it exists, is true when the face must be flipped (multiplies by -1.).
+  const bool *orient;
+  bool *orient_allocated;
   int (*Apply)(CeedElemRestriction, const CeedInt, const CeedInt,
                const CeedInt, CeedInt, CeedInt, CeedTransposeMode, CeedVector,
                CeedVector, CeedRequest *);
@@ -70,6 +73,10 @@ CEED_INTERN int CeedVectorCreate_Ref(CeedInt n, CeedVector vec);
 CEED_INTERN int CeedElemRestrictionCreate_Ref(CeedMemType mem_type,
     CeedCopyMode copy_mode, const CeedInt *indices, CeedElemRestriction r);
 
+CEED_INTERN int CeedElemRestrictionCreateOriented_Ref(CeedMemType mem_type,
+    CeedCopyMode copy_mode, const CeedInt *indices,
+    const bool *orient, CeedElemRestriction r);
+
 CEED_INTERN int CeedBasisCreateTensorH1_Ref(CeedInt dim, CeedInt P_1d,
     CeedInt Q_1d, const CeedScalar *interp_1d, const CeedScalar *grad_1d,
     const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis basis);
diff --git a/include/ceed-impl.h b/include/ceed-impl.h
@@ -105,6 +105,9 @@ struct Ceed_private {
   int (*VectorCreate)(CeedInt, CeedVector);
   int (*ElemRestrictionCreate)(CeedMemType, CeedCopyMode,
                                const CeedInt *, CeedElemRestriction);
+  int (*ElemRestrictionCreateOriented)(CeedMemType, CeedCopyMode,
+                                       const CeedInt *, const bool *,
+                                       CeedElemRestriction);
   int (*ElemRestrictionCreateBlocked)(CeedMemType, CeedCopyMode,
                                       const CeedInt *, CeedElemRestriction);
   int (*BasisCreateTensorH1)(CeedInt, CeedInt, CeedInt, const CeedScalar *,
@@ -177,6 +180,7 @@ struct CeedElemRestriction_private {
   CeedInt *strides;      /* strides between [nodes, components, elements] */
   CeedInt layout[3];     /* E-vector layout [nodes, components, elements] */
   uint64_t num_readers;  /* number of instances of offset read only access */
+  bool is_oriented;       /* flag for oriented restriction */
   void *data;            /* place for the backend to store any data */
 };
 
diff --git a/include/ceed/backend.h b/include/ceed/backend.h
@@ -148,6 +148,8 @@ CEED_EXTERN int CeedElemRestrictionRestoreOffsets(CeedElemRestriction rstr,
     const CeedInt **offsets);
 CEED_EXTERN int CeedElemRestrictionIsStrided(CeedElemRestriction rstr,
     bool *is_strided);
+CEED_EXTERN int CeedElemRestrictionIsOriented(CeedElemRestriction rstr,
+    bool *is_oriented);
 CEED_EXTERN int CeedElemRestrictionHasBackendStrides(CeedElemRestriction rstr,
     bool *has_backend_strides);
 CEED_EXTERN int CeedElemRestrictionGetELayout(CeedElemRestriction rstr,
diff --git a/include/ceed/ceed.h b/include/ceed/ceed.h
@@ -415,6 +415,10 @@ CEED_EXTERN int CeedElemRestrictionCreate(Ceed ceed, CeedInt num_elem,
     CeedInt elem_size, CeedInt num_comp, CeedInt comp_stride, CeedInt l_size,
     CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets,
     CeedElemRestriction *rstr);
+CEED_EXTERN int CeedElemRestrictionCreateOriented(Ceed ceed, CeedInt num_elem,
+    CeedInt elem_size, CeedInt num_comp, CeedInt comp_stride, CeedInt l_size,
+    CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets,
+    const bool *orient, CeedElemRestriction *rstr);
 CEED_EXTERN int CeedElemRestrictionCreateStrided(Ceed ceed,
     CeedInt num_elem, CeedInt elem_size, CeedInt num_comp, CeedInt l_size,
     const CeedInt strides[3], CeedElemRestriction *rstr);
diff --git a/interface/ceed-elemrestriction.c b/interface/ceed-elemrestriction.c
@@ -152,6 +152,21 @@ int CeedElemRestrictionIsStrided(CeedElemRestriction rstr, bool *is_strided) {
   return CEED_ERROR_SUCCESS;
 }
 
+/**
+  @brief Get oriented status of a CeedElemRestriction
+
+  @param rstr             CeedElemRestriction
+  @param[out] is_oriented  Variable to store oriented status, 1 if oriented else 0
+
+  @return An error code: 0 - success, otherwise - failure
+
+  @ref Backend
+**/
+int CeedElemRestrictionIsOriented(CeedElemRestriction rstr, bool *is_oriented) {
+  *is_oriented = rstr->is_oriented;
+  return CEED_ERROR_SUCCESS;
+}
+
 /**
   @brief Get the backend stride status of a CeedElemRestriction
 
@@ -352,11 +367,84 @@ int CeedElemRestrictionCreate(Ceed ceed, CeedInt num_elem, CeedInt elem_size,
   (*rstr)->l_size = l_size;
   (*rstr)->num_blk = num_elem;
   (*rstr)->blk_size = 1;
+  (*rstr)->is_oriented = 0;
   ierr = ceed->ElemRestrictionCreate(mem_type, copy_mode, offsets, *rstr);
   CeedChk(ierr);
   return CEED_ERROR_SUCCESS;
 }
 
+/**
+  @brief Create a CeedElemRestriction with orientation sign
+
+  @param ceed         A Ceed object where the CeedElemRestriction will be created
+  @param num_elem     Number of elements described in the @a offsets array
+  @param elem_size    Size (number of "nodes") per element
+  @param num_comp     Number of field components per interpolation node
+                        (1 for scalar fields)
+  @param 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 l_size       The size of the L-vector. This vector may be larger than
+                        the elements and fields given by this restriction.
+  @param mem_type     Memory type of the @a offsets array, see CeedMemType
+  @param copy_mode    Copy mode for the @a offsets array, see CeedCopyMode
+  @param 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 orient       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, CeedInt l_size,
+                                      CeedMemType mem_type, CeedCopyMode copy_mode,
+                                      const CeedInt *offsets, const bool *orient,
+                                      CeedElemRestriction *rstr) {
+  int ierr;
+
+  if (!ceed->ElemRestrictionCreateOriented) {
+    Ceed delegate;
+    ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction");
+    CeedChk(ierr);
+
+    if (!delegate)
+      // LCOV_EXCL_START
+      return CeedError(ceed, CEED_ERROR_UNSUPPORTED,
+                       "Backend does not implement ElemRestrictionCreateOriented");
+    // LCOV_EXCL_STOP
+
+    ierr = CeedElemRestrictionCreateOriented(delegate, num_elem, elem_size,
+           num_comp, comp_stride, l_size, mem_type, copy_mode, offsets,
+           orient, rstr); CeedChk(ierr);
+    return CEED_ERROR_SUCCESS;
+  }
+
+  ierr = CeedCalloc(1, rstr); CeedChk(ierr);
+  (*rstr)->ceed = ceed;
+  ierr = CeedReference(ceed); CeedChk(ierr);
+  (*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)->num_blk = num_elem;
+  (*rstr)->blk_size = 1;
+  (*rstr)->is_oriented = 1;
+  ierr = ceed->ElemRestrictionCreateOriented(mem_type, copy_mode,
+         offsets, orient, *rstr); CeedChk(ierr);
+  return CEED_ERROR_SUCCESS;
+}
+
 /**
   @brief Create a strided CeedElemRestriction
 
@@ -414,6 +502,7 @@ int CeedElemRestrictionCreateStrided(Ceed ceed, CeedInt num_elem,
   (*rstr)->l_size = l_size;
   (*rstr)->num_blk = num_elem;
   (*rstr)->blk_size = 1;
+  (*rstr)->is_oriented = 0;
   ierr = CeedMalloc(3, &(*rstr)->strides); CeedChk(ierr);
   for (int i=0; i<3; i++)
     (*rstr)->strides[i] = strides[i];
@@ -500,6 +589,7 @@ int CeedElemRestrictionCreateBlocked(Ceed ceed, CeedInt num_elem,
   (*rstr)->l_size = l_size;
   (*rstr)->num_blk = num_blk;
   (*rstr)->blk_size = blk_size;
+  (*rstr)->is_oriented = 0;
   ierr = ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER,
          (const CeedInt *) blk_offsets, *rstr); CeedChk(ierr);
   if (copy_mode == CEED_OWN_POINTER) {
@@ -565,6 +655,7 @@ int CeedElemRestrictionCreateBlockedStrided(Ceed ceed, CeedInt num_elem,
   (*rstr)->l_size = l_size;
   (*rstr)->num_blk = num_blk;
   (*rstr)->blk_size = blk_size;
+  (*rstr)->is_oriented = 0;
   ierr = CeedMalloc(3, &(*rstr)->strides); CeedChk(ierr);
   for (int i=0; i<3; i++)
     (*rstr)->strides[i] = strides[i];
diff --git a/interface/ceed.c b/interface/ceed.c
@@ -836,6 +836,7 @@ int CeedInit(const char *resource, Ceed *ceed) {
     CEED_FTABLE_ENTRY(Ceed, Destroy),
     CEED_FTABLE_ENTRY(Ceed, VectorCreate),
     CEED_FTABLE_ENTRY(Ceed, ElemRestrictionCreate),
+    CEED_FTABLE_ENTRY(Ceed, ElemRestrictionCreateOriented),
     CEED_FTABLE_ENTRY(Ceed, ElemRestrictionCreateBlocked),
     CEED_FTABLE_ENTRY(Ceed, BasisCreateTensorH1),
     CEED_FTABLE_ENTRY(Ceed, BasisCreateH1),
diff --git a/tests/t200-elemrestriction.c b/tests/t200-elemrestriction.c
@@ -26,7 +26,6 @@ int main(int argc, char **argv) {
   CeedElemRestrictionCreate(ceed, num_elem, 2, 1, 1, num_elem+1, CEED_MEM_HOST,
                             CEED_USE_POINTER, ind, &r);
   CeedVectorCreate(ceed, num_elem*2, &y);
-  CeedVectorSetValue(y, 0); // Allocates array
   CeedElemRestrictionApply(r, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE);
 
   CeedVectorGetArrayRead(y, CEED_MEM_HOST, &yy);
diff --git a/tests/t201-elemrestriction.c b/tests/t201-elemrestriction.c
@@ -23,7 +23,6 @@ int main(int argc, char **argv) {
 
   CeedElemRestrictionCreateStrided(ceed, num_elem, 2, 1, num_elem*2, strides, &r);
   CeedVectorCreate(ceed, num_elem*2, &y);
-  CeedVectorSetValue(y, 0); // Allocates array
   CeedElemRestrictionApply(r, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE);
 
   CeedVectorGetArrayRead(y, CEED_MEM_HOST, &yy);
diff --git a/tests/t202-elemrestriction.c b/tests/t202-elemrestriction.c
@@ -32,7 +32,6 @@ int main(int argc, char **argv) {
                                    num_elem+1, CEED_MEM_HOST, CEED_USE_POINTER,
                                    ind, &r);
   CeedVectorCreate(ceed, num_blk*blk_size*elem_size, &y);
-  CeedVectorSetValue(y, 0); // Allocates array
 
   // NoTranspose
   CeedElemRestrictionApply(r, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE);
diff --git a/tests/t203-elemrestriction.c b/tests/t203-elemrestriction.c
@@ -36,7 +36,6 @@ int main(int argc, char **argv) {
                                    num_elem+1, num_comp*(num_elem+1), CEED_MEM_HOST,
                                    CEED_USE_POINTER, ind, &r);
   CeedVectorCreate(ceed, num_comp*num_blk*blk_size*elem_size, &y);
-  CeedVectorSetValue(y, 0); // Allocates array
 
   // NoTranspose
   CeedElemRestrictionApply(r, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE);
diff --git a/tests/t204-elemrestriction.c b/tests/t204-elemrestriction.c
@@ -32,7 +32,6 @@ int main(int argc, char **argv) {
                             CEED_MEM_HOST,
                             CEED_USE_POINTER, ind, &r);
   CeedVectorCreate(ceed, 2*(num_elem*2), &y);
-  CeedVectorSetValue(y, 0); // Allocates array
 
   // Restrict
   CeedElemRestrictionApply(r, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE);
diff --git a/tests/t205-elemrestriction.c b/tests/t205-elemrestriction.c
@@ -32,7 +32,6 @@ int main(int argc, char **argv) {
                             CEED_MEM_HOST,
                             CEED_USE_POINTER, ind, &r);
   CeedVectorCreate(ceed, 2*(num_elem*2), &y);
-  CeedVectorSetValue(y, 0); // Allocates array
 
   // Restrict
   CeedElemRestrictionApply(r, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE);
diff --git a/tests/t208-elemrestriction.c b/tests/t208-elemrestriction.c
@@ -31,7 +31,6 @@ int main(int argc, char **argv) {
                                    num_elem+1, CEED_MEM_HOST, CEED_USE_POINTER,
                                    ind, &r);
   CeedVectorCreate(ceed, blk_size*elem_size, &y);
-  CeedVectorSetValue(y, 0); // Allocates array
 
   // NoTranspose
   CeedElemRestrictionApplyBlock(r, 1, CEED_NOTRANSPOSE, x, y,
diff --git a/tests/t213-elemrestriction.c b/tests/t213-elemrestriction.c
@@ -40,7 +40,6 @@ int main(int argc, char **argv) {
                                    num_elem+1, num_comp*(num_elem+1), CEED_MEM_HOST,
                                    CEED_OWN_POINTER, ceed_ind, &r);
   CeedVectorCreate(ceed, num_comp*num_blk*blk_size*elem_size, &y);
-  CeedVectorSetValue(y, 0); // Allocates array
 
   // NoTranspose
   CeedElemRestrictionApply(r, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE);
diff --git a/tests/t220-elemrestriction.c b/tests/t220-elemrestriction.c
@@ -0,0 +1,54 @@
+/// @file
+/// Test creation, use, and destruction of an element restriction oriented
+/// \test Test creation, use, and destruction of an element restriction oriented
+#include <ceed.h>
+
+int main(int argc, char **argv) {
+  Ceed ceed;
+  CeedVector x, y;
+  CeedInt num_elem = 6, P = 2, dim = 1;
+  CeedInt ind[P*num_elem];
+  bool orient[P*num_elem];
+  CeedScalar a[num_elem+1];
+  const CeedScalar *yy;
+  CeedElemRestriction r;
+
+  CeedInit(argv[1], &ceed);
+
+  CeedVectorCreate(ceed, num_elem+1, &x);
+  for (CeedInt i=0; i<num_elem+1; i++)
+    a[i] = 10 + i;
+  CeedVectorSetArray(x, CEED_MEM_HOST, CEED_USE_POINTER, a);
+
+  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, &r);
+  CeedVectorCreate(ceed, num_elem*2, &y);
+  CeedElemRestrictionApply(r, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE);
+
+  CeedVectorGetArrayRead(y, CEED_MEM_HOST, &yy);
+  for (CeedInt i=0; i<num_elem; i++) {
+    for (CeedInt j=0; j<P; j++) {
+      CeedInt k = j + P*i;
+      if (10+(k+1)/2 != yy[k] * CeedIntPow(-1, i%2))
+        // LCOV_EXCL_START
+        printf("Error in restricted array y[%d] = %f",
+               k, (CeedScalar)yy[k]);
+      // LCOV_EXCL_STOP
+    }
+  }
+  CeedVectorRestoreArrayRead(y, &yy);
+
+  CeedVectorDestroy(&x);
+  CeedVectorDestroy(&y);
+  CeedElemRestrictionDestroy(&r);
+  CeedDestroy(&ceed);
+  return 0;
+}
