CeedElemRestriction for H(curl)

AUTHORS

@@ -6,6 +6,7 @@ Jean-Sylvain Camier
 Veselin Dobrev
 Yohann Dudouit
 Leila Ghaffari
+Sebastian Grimberg
 Tzanio Kolev
 David Medina
 Will Pazner

backends/blocked/ceed-blocked-operator.c

@@ -41,29 +41,54 @@ static int CeedOperatorSetupFields_Blocked(CeedQFunction qf, CeedOperator op, bo
     CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode));
 
     if (eval_mode != CEED_EVAL_WEIGHT) {
+      Ceed ceed_rstr;
       CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &r));
-      CeedCallBackend(CeedElemRestrictionGetCeed(r, &ceed));
+      CeedCallBackend(CeedElemRestrictionGetCeed(r, &ceed_rstr));
       CeedSize l_size;
       CeedInt  num_elem, elem_size, comp_stride;
       CeedCallBackend(CeedElemRestrictionGetNumElements(r, &num_elem));
       CeedCallBackend(CeedElemRestrictionGetElementSize(r, &elem_size));
       CeedCallBackend(CeedElemRestrictionGetLVectorSize(r, &l_size));
       CeedCallBackend(CeedElemRestrictionGetNumComponents(r, &num_comp));
-
-      bool strided;
-      CeedCallBackend(CeedElemRestrictionIsStrided(r, &strided));
-      if (strided) {
-        CeedInt strides[3];
-        CeedCallBackend(CeedElemRestrictionGetStrides(r, &strides));
-        CeedCallBackend(
-            CeedElemRestrictionCreateBlockedStrided(ceed, num_elem, elem_size, blk_size, num_comp, l_size, strides, &blk_restr[i + start_e]));
-      } else {
-        const CeedInt *offsets = NULL;
-        CeedCallBackend(CeedElemRestrictionGetOffsets(r, CEED_MEM_HOST, &offsets));
-        CeedCallBackend(CeedElemRestrictionGetCompStride(r, &comp_stride));
-        CeedCallBackend(CeedElemRestrictionCreateBlocked(ceed, num_elem, elem_size, blk_size, num_comp, comp_stride, l_size, CEED_MEM_HOST,
-                                                         CEED_COPY_VALUES, offsets, &blk_restr[i + start_e]));
-        CeedCallBackend(CeedElemRestrictionRestoreOffsets(r, &offsets));
+      CeedCallBackend(CeedElemRestrictionGetCompStride(r, &comp_stride));
+
+      CeedRestrictionType rstr_type;
+      CeedCallBackend(CeedElemRestrictionGetType(r, &rstr_type));
+      switch (rstr_type) {
+        case CEED_RESTRICTION_STANDARD: {
+          const CeedInt *offsets = NULL;
+          CeedCallBackend(CeedElemRestrictionGetOffsets(r, CEED_MEM_HOST, &offsets));
+          CeedCallBackend(CeedElemRestrictionCreateBlocked(ceed_rstr, num_elem, elem_size, blk_size, num_comp, comp_stride, l_size, CEED_MEM_HOST,
+                                                           CEED_COPY_VALUES, offsets, &blk_restr[i + start_e]));
+          CeedCallBackend(CeedElemRestrictionRestoreOffsets(r, &offsets));
+        } break;
+        case CEED_RESTRICTION_ORIENTED: {
+          const CeedInt *offsets = NULL;
+          const bool    *orients = NULL;
+          CeedCallBackend(CeedElemRestrictionGetOffsets(r, CEED_MEM_HOST, &offsets));
+          CeedCallBackend(CeedElemRestrictionGetOrientations(r, CEED_MEM_HOST, &orients));
+          CeedCallBackend(CeedElemRestrictionCreateBlockedOriented(ceed_rstr, num_elem, elem_size, blk_size, num_comp, comp_stride, l_size,
+                                                                   CEED_MEM_HOST, CEED_COPY_VALUES, offsets, orients, &blk_restr[i + start_e]));
+          CeedCallBackend(CeedElemRestrictionRestoreOffsets(r, &offsets));
+          CeedCallBackend(CeedElemRestrictionRestoreOrientations(r, &orients));
+        } break;
+        case CEED_RESTRICTION_CURL_ORIENTED: {
+          const CeedInt  *offsets      = NULL;
+          const CeedInt8 *curl_orients = NULL;
+          CeedCallBackend(CeedElemRestrictionGetOffsets(r, CEED_MEM_HOST, &offsets));
+          CeedCallBackend(CeedElemRestrictionGetCurlOrientations(r, CEED_MEM_HOST, &curl_orients));
+          CeedCallBackend(CeedElemRestrictionCreateBlockedCurlOriented(ceed_rstr, num_elem, elem_size, blk_size, num_comp, comp_stride, l_size,
+                                                                       CEED_MEM_HOST, CEED_COPY_VALUES, offsets, curl_orients,
+                                                                       &blk_restr[i + start_e]));
+          CeedCallBackend(CeedElemRestrictionRestoreOffsets(r, &offsets));
+          CeedCallBackend(CeedElemRestrictionRestoreCurlOrientations(r, &curl_orients));
+        } break;
+        case CEED_RESTRICTION_STRIDED: {
+          CeedInt strides[3];
+          CeedCallBackend(CeedElemRestrictionGetStrides(r, &strides));
+          CeedCallBackend(
+              CeedElemRestrictionCreateBlockedStrided(ceed_rstr, num_elem, elem_size, blk_size, num_comp, l_size, strides, &blk_restr[i + start_e]));
+        } break;
       }
       CeedCallBackend(CeedElemRestrictionCreateVector(blk_restr[i + start_e], NULL, &e_vecs_full[i + start_e]));
     }

backends/cuda-ref/ceed-cuda-ref-restriction.c

@@ -209,7 +209,8 @@ static int CeedElemRestrictionOffset_Cuda(const CeedElemRestriction r, const Cee
 //------------------------------------------------------------------------------
 // Create restriction
 //------------------------------------------------------------------------------
-int CeedElemRestrictionCreate_Cuda(CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *indices, CeedElemRestriction r) {
+int CeedElemRestrictionCreate_Cuda(CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *indices, const bool *orients,
+                                   const CeedInt8 *curl_orients, CeedElemRestriction r) {
   Ceed ceed;
   CeedCallBackend(CeedElemRestrictionGetCeed(r, &ceed));
   CeedElemRestriction_Cuda *impl;
@@ -222,6 +223,11 @@ int CeedElemRestrictionCreate_Cuda(CeedMemType mem_type, CeedCopyMode copy_mode,
   CeedInt strides[3]  = {1, size, elem_size};
   CeedInt comp_stride = 1;
 
+  CeedRestrictionType rstr_type;
+  CeedCallBackend(CeedElemRestrictionGetType(r, &rstr_type));
+  CeedCheck(rstr_type != CEED_RESTRICTION_ORIENTED && rstr_type != CEED_RESTRICTION_CURL_ORIENTED, ceed, CEED_ERROR_BACKEND,
+            "Backend does not implement CeedElemRestrictionCreateOriented or CeedElemRestrictionCreateCurlOriented");
+
   // Stride data
   bool is_strided;
   CeedCallBackend(CeedElemRestrictionIsStrided(r, &is_strided));
@@ -323,6 +329,7 @@ int CeedElemRestrictionCreate_Cuda(CeedMemType mem_type, CeedCopyMode copy_mode,
   // Register backend functions
   CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "Apply", CeedElemRestrictionApply_Cuda));
   CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "ApplyUnsigned", CeedElemRestrictionApply_Cuda));
+  CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "ApplyUnoriented", CeedElemRestrictionApply_Cuda));
   CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "GetOffsets", CeedElemRestrictionGetOffsets_Cuda));
   CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "Destroy", CeedElemRestrictionDestroy_Cuda));
   return CEED_ERROR_SUCCESS;

backends/cuda-ref/ceed-cuda-ref.h

@@ -113,7 +113,8 @@ CEED_INTERN int CeedGetCublasHandle_Cuda(Ceed ceed, cublasHandle_t *handle);
 
 CEED_INTERN int CeedVectorCreate_Cuda(CeedSize n, CeedVector vec);
 
-CEED_INTERN int CeedElemRestrictionCreate_Cuda(CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *indices, CeedElemRestriction r);
+CEED_INTERN int CeedElemRestrictionCreate_Cuda(CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *indices, const bool *orients,
+                                               const CeedInt8 *curl_orients, CeedElemRestriction r);
 
 CEED_INTERN int CeedBasisCreateTensorH1_Cuda(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);

backends/hip-ref/ceed-hip-ref-restriction.c

@@ -207,7 +207,8 @@ static int CeedElemRestrictionOffset_Hip(const CeedElemRestriction r, const Ceed
 //------------------------------------------------------------------------------
 // Create restriction
 //------------------------------------------------------------------------------
-int CeedElemRestrictionCreate_Hip(CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *indices, CeedElemRestriction r) {
+int CeedElemRestrictionCreate_Hip(CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *indices, const bool *orients,
+                                  const CeedInt8 *curl_orients, CeedElemRestriction r) {
   Ceed ceed;
   CeedCallBackend(CeedElemRestrictionGetCeed(r, &ceed));
   CeedElemRestriction_Hip *impl;
@@ -220,6 +221,11 @@ int CeedElemRestrictionCreate_Hip(CeedMemType mem_type, CeedCopyMode copy_mode,
   CeedInt strides[3]  = {1, size, elem_size};
   CeedInt comp_stride = 1;
 
+  CeedRestrictionType rstr_type;
+  CeedCallBackend(CeedElemRestrictionGetType(r, &rstr_type));
+  CeedCheck(rstr_type != CEED_RESTRICTION_ORIENTED && rstr_type != CEED_RESTRICTION_CURL_ORIENTED, ceed, CEED_ERROR_BACKEND,
+            "Backend does not implement CeedElemRestrictionCreateOriented or CeedElemRestrictionCreateCurlOriented");
+
   // Stride data
   bool is_strided;
   CeedCallBackend(CeedElemRestrictionIsStrided(r, &is_strided));
@@ -321,6 +327,7 @@ int CeedElemRestrictionCreate_Hip(CeedMemType mem_type, CeedCopyMode copy_mode,
   // Register backend functions
   CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "Apply", CeedElemRestrictionApply_Hip));
   CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "ApplyUnsigned", CeedElemRestrictionApply_Hip));
+  CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "ApplyUnoriented", CeedElemRestrictionApply_Hip));
   CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "GetOffsets", CeedElemRestrictionGetOffsets_Hip));
   CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "Destroy", CeedElemRestrictionDestroy_Hip));
   return CEED_ERROR_SUCCESS;

backends/hip-ref/ceed-hip-ref.h

@@ -118,7 +118,8 @@ CEED_INTERN int CeedGetHipblasHandle_Hip(Ceed ceed, hipblasHandle_t *handle);
 
 CEED_INTERN int CeedVectorCreate_Hip(CeedSize n, CeedVector vec);
 
-CEED_INTERN int CeedElemRestrictionCreate_Hip(CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *indices, CeedElemRestriction r);
+CEED_INTERN int CeedElemRestrictionCreate_Hip(CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *indices, const bool *orients,
+                                              const CeedInt8 *curl_orients, CeedElemRestriction r);
 
 CEED_INTERN int CeedBasisCreateTensorH1_Hip(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);

backends/magma/ceed-magma-restriction.c

@@ -151,7 +151,8 @@ static int CeedElemRestrictionDestroy_Magma(CeedElemRestriction r) {
   return CEED_ERROR_SUCCESS;
 }
 
-int CeedElemRestrictionCreate_Magma(CeedMemType mtype, CeedCopyMode cmode, const CeedInt *offsets, CeedElemRestriction r) {
+int CeedElemRestrictionCreate_Magma(CeedMemType mtype, CeedCopyMode cmode, const CeedInt *offsets, const bool *orients, const CeedInt8 *curl_orients,
+                                    CeedElemRestriction r) {
   Ceed ceed;
   CeedCallBackend(CeedElemRestrictionGetCeed(r, &ceed));
 
@@ -163,6 +164,11 @@ int CeedElemRestrictionCreate_Magma(CeedMemType mtype, CeedCopyMode cmode, const
   CeedCallBackend(CeedElemRestrictionGetElementSize(r, &elemsize));
   CeedInt size = elemsize * nelem;
 
+  CeedRestrictionType rstr_type;
+  CeedCallBackend(CeedElemRestrictionGetType(r, &rstr_type));
+  CeedCheck(rstr_type != CEED_RESTRICTION_ORIENTED && rstr_type != CEED_RESTRICTION_CURL_ORIENTED, ceed, CEED_ERROR_BACKEND,
+            "Backend does not implement CeedElemRestrictionCreateOriented or CeedElemRestrictionCreateCurlOriented");
+
   CeedElemRestriction_Magma *impl;
   CeedCallBackend(CeedCalloc(1, &impl));
 
@@ -261,6 +267,7 @@ int CeedElemRestrictionCreate_Magma(CeedMemType mtype, CeedCopyMode cmode, const
   CeedCallBackend(CeedElemRestrictionSetELayout(r, layout));
   CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "Apply", CeedElemRestrictionApply_Magma));
   CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "ApplyUnsigned", CeedElemRestrictionApply_Magma));
+  CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "ApplyUnoriented", CeedElemRestrictionApply_Magma));
   CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "GetOffsets", CeedElemRestrictionGetOffsets_Magma));
   CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "Destroy", CeedElemRestrictionDestroy_Magma));
   CeedCallBackend(CeedFree(&restriction_kernel_path));

backends/magma/ceed-magma.h

@@ -134,7 +134,8 @@ CEED_INTERN int CeedBasisCreateTensorH1_Magma(CeedInt dim, CeedInt P1d, CeedInt
 CEED_INTERN int CeedBasisCreateH1_Magma(CeedElemTopology topo, CeedInt dim, CeedInt ndof, CeedInt nqpts, const CeedScalar *interp,
                                         const CeedScalar *grad, const CeedScalar *qref, const CeedScalar *qweight, CeedBasis basis);
 
-CEED_INTERN int CeedElemRestrictionCreate_Magma(CeedMemType mtype, CeedCopyMode cmode, const CeedInt *offsets, CeedElemRestriction r);
+CEED_INTERN int CeedElemRestrictionCreate_Magma(CeedMemType mtype, CeedCopyMode cmode, const CeedInt *indices, const bool *orients,
+                                                const CeedInt8 *curl_orients, CeedElemRestriction r);
 
 // comment the line below to use the default magma_is_devptr function
 #define magma_is_devptr magma_isdevptr

backends/occa/ceed-occa-elem-restriction.cpp

@@ -297,14 +297,21 @@ int ElemRestriction::registerCeedFunction(Ceed ceed, CeedElemRestriction r, cons
   return CeedSetBackendFunction(ceed, "ElemRestriction", r, fname, f);
 }
 
-int ElemRestriction::ceedCreate(CeedMemType memType, CeedCopyMode copyMode, const CeedInt *indicesInput, CeedElemRestriction r) {
+int ElemRestriction::ceedCreate(CeedMemType memType, CeedCopyMode copyMode, const CeedInt *indicesInput, const bool *orientsInput,
+                                const CeedInt8 *curlOrientsInput, CeedElemRestriction r) {
   Ceed ceed;
   CeedCallBackend(CeedElemRestrictionGetCeed(r, &ceed));
 
   if ((memType != CEED_MEM_DEVICE) && (memType != CEED_MEM_HOST)) {
     return staticCeedError("Only HOST and DEVICE CeedMemType supported");
   }
 
+  CeedRestrictionType rstr_type;
+  CeedCallBackend(CeedElemRestrictionGetType(r, &rstr_type));
+  if ((rstr_type == CEED_RESTRICTION_ORIENTED) || (rstr_type == CEED_RESTRICTION_CURL_ORIENTED)) {
+    return staticCeedError("(OCCA) Backend does not implement CeedElemRestrictionCreateOriented or CeedElemRestrictionCreateCurlOriented");
+  }
+
   ElemRestriction *elemRestriction = new ElemRestriction();
   CeedCallBackend(CeedElemRestrictionSetData(r, elemRestriction));
 
@@ -317,17 +324,14 @@ int ElemRestriction::ceedCreate(CeedMemType memType, CeedCopyMode copyMode, cons
 
   CeedOccaRegisterFunction(r, "Apply", ElemRestriction::ceedApply);
   CeedOccaRegisterFunction(r, "ApplyUnsigned", ElemRestriction::ceedApply);
+  CeedOccaRegisterFunction(r, "ApplyUnoriented", ElemRestriction::ceedApply);
   CeedOccaRegisterFunction(r, "ApplyBlock", ElemRestriction::ceedApplyBlock);
   CeedOccaRegisterFunction(r, "GetOffsets", ElemRestriction::ceedGetOffsets);
   CeedOccaRegisterFunction(r, "Destroy", ElemRestriction::ceedDestroy);
 
   return CEED_ERROR_SUCCESS;
 }
 
-int ElemRestriction::ceedCreateBlocked(CeedMemType memType, CeedCopyMode copyMode, const CeedInt *indicesInput, CeedElemRestriction r) {
-  return staticCeedError("(OCCA) Backend does not implement CeedElemRestrictionCreateBlocked");
-}
-
 int ElemRestriction::ceedApply(CeedElemRestriction r, CeedTransposeMode tmode, CeedVector u, CeedVector v, CeedRequest *request) {
   ElemRestriction *elemRestriction = ElemRestriction::from(r);
   Vector          *uVector         = Vector::from(u);

backends/occa/ceed-occa-elem-restriction.hpp

@@ -77,9 +77,8 @@ class ElemRestriction : public CeedObject {
   //---[ Ceed Callbacks ]-----------
   static int registerCeedFunction(Ceed ceed, CeedElemRestriction r, const char *fname, ceed::occa::ceedFunction f);
 
-  static int ceedCreate(CeedMemType memType, CeedCopyMode copyMode, const CeedInt *indicesInput, CeedElemRestriction r);
-
-  static int ceedCreateBlocked(CeedMemType memType, CeedCopyMode copyMode, const CeedInt *indicesInput, CeedElemRestriction r);
+  static int ceedCreate(CeedMemType memType, CeedCopyMode copyMode, const CeedInt *indicesInput, const bool *orientsInput,
+                        const CeedInt8 *curlOrientsInput, CeedElemRestriction r);
 
   static int ceedApply(CeedElemRestriction r, CeedTransposeMode tmode, CeedVector u, CeedVector v, CeedRequest *request);
 

backends/occa/ceed-occa.cpp

@@ -290,7 +290,6 @@ static int registerMethods(Ceed ceed) {
   CeedOccaRegisterBaseFunction("BasisCreateTensorH1", ceed::occa::TensorBasis::ceedCreate);
   CeedOccaRegisterBaseFunction("BasisCreateH1", ceed::occa::SimplexBasis::ceedCreate);
   CeedOccaRegisterBaseFunction("ElemRestrictionCreate", ceed::occa::ElemRestriction::ceedCreate);
-  CeedOccaRegisterBaseFunction("ElemRestrictionCreateBlocked", ceed::occa::ElemRestriction::ceedCreateBlocked);
   CeedOccaRegisterBaseFunction("QFunctionCreate", ceed::occa::QFunction::ceedCreate);
   CeedOccaRegisterBaseFunction("QFunctionContextCreate", ceed::occa::QFunctionContext::ceedCreate);
   CeedOccaRegisterBaseFunction("OperatorCreate", ceed::occa::Operator::ceedCreate);

backends/occa/ceed-occa.h

@@ -142,12 +142,8 @@ CEED_INTERN int CeedQFunctionCreate_Occa(CeedQFunction qf);
 CEED_INTERN int CeedQFunctionContextCreate_Occa(CeedQFunctionContext ctx);
 
 // *****************************************************************************
-CEED_INTERN int CeedElemRestrictionCreate_Occa(const CeedMemType mtype, const CeedCopyMode cmode, const CeedInt *indices,
-                                               const CeedElemRestriction res);
-
-// *****************************************************************************
-CEED_INTERN int CeedElemRestrictionCreateBlocked_Occa(const CeedMemType mtype, const CeedCopyMode cmode, const CeedInt *indices,
-                                                      const CeedElemRestriction res);
+CEED_INTERN int CeedElemRestrictionCreate_Occa(const CeedMemType mtype, const CeedCopyMode cmode, const CeedInt *indices, const bool *orients,
+                                               const CeedInt8 *curl_orients, const CeedElemRestriction res);
 
 // *****************************************************************************
 CEED_INTERN int CeedVectorCreate_Occa(CeedInt n, CeedVector vec);