From fc0567d9dd43c4a8f0744c44d73ac6beb4777ae6 Mon Sep 17 00:00:00 2001 From: rezgarshakeri Date: Tue, 1 Feb 2022 13:47:04 -0700 Subject: [PATCH 1/8] backends/ref: Added oriented element restriction. --- backends/ref/ceed-ref-restriction.c | 42 +++++++++++++++- backends/ref/ceed-ref.c | 3 ++ backends/ref/ceed-ref.h | 7 +++ include/ceed-impl.h | 3 ++ include/ceed/ceed.h | 4 ++ interface/ceed-elemrestriction.c | 74 +++++++++++++++++++++++++++++ interface/ceed.c | 1 + 7 files changed, 132 insertions(+), 2 deletions(-) diff --git a/backends/ref/ceed-ref-restriction.c b/backends/ref/ceed-ref-restriction.c index b125fc818a..23e82e2a15 100644 --- a/backends/ref/ceed-ref-restriction.c +++ b/backends/ref/ceed-ref-restriction.c @@ -91,7 +91,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] * + (impl->orient && impl->orient[i+elem_size*e] ? -1. : 1.); } } else { // Restriction from E-vector to L-vector @@ -137,7 +138,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] * + (impl->orient && impl->orient[j+e*elem_size] ? -1. : 1.); } } ierr = CeedVectorRestoreArrayRead(u, &uu); CeedChkBackend(ierr); @@ -456,4 +458,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 index 0c3dfcce41..671f4649c2 100644 --- 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 index 0b472ab22e..19e3e1870d 100644 --- 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 index c3be7fa134..2b8ccefba0 100644 --- 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 *, diff --git a/include/ceed/ceed.h b/include/ceed/ceed.h index 0d47b13618..2fff9be788 100644 --- 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 index a76dd19f2c..52bd36bf91 100644 --- a/interface/ceed-elemrestriction.c +++ b/interface/ceed-elemrestriction.c @@ -357,6 +357,80 @@ int CeedElemRestrictionCreate(Ceed ceed, CeedInt num_elem, CeedInt elem_size, 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 scale An scalar value that scales the dofs in assembly. + @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->ElemRestrictionCreate) { + Ceed delegate; + ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); + CeedChk(ierr); + + if (!delegate) + // LCOV_EXCL_START + return CeedError(ceed, CEED_ERROR_UNSUPPORTED, + "Backend does not support 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; + ierr = ceed->ElemRestrictionCreateOriented(mem_type, copy_mode, offsets, orient, + *rstr); + CeedChk(ierr); + return CEED_ERROR_SUCCESS; +} + /** @brief Create a strided CeedElemRestriction diff --git a/interface/ceed.c b/interface/ceed.c index a3e42daa31..690a8080fb 100644 --- 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), From cf6be9071e5bf8fe8883d9b8f62fff643818f7ab Mon Sep 17 00:00:00 2001 From: rezgarshakeri Date: Tue, 1 Feb 2022 15:06:46 -0700 Subject: [PATCH 2/8] tests: added t220-elemrestriction.c --- tests/t220-elemrestriction.c | 55 ++++++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) create mode 100644 tests/t220-elemrestriction.c diff --git a/tests/t220-elemrestriction.c b/tests/t220-elemrestriction.c new file mode 100644 index 0000000000..7f8d3a8347 --- /dev/null +++ b/tests/t220-elemrestriction.c @@ -0,0 +1,55 @@ +/// @file +/// Test creation, use, and destruction of an element restriction oriented +/// \test Test creation, use, and destruction of an element restriction oriented +#include + +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 Date: Wed, 2 Feb 2022 12:45:52 -0700 Subject: [PATCH 3/8] update ceed-elemrestriction.c: fixed formatting --- interface/ceed-elemrestriction.c | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/interface/ceed-elemrestriction.c b/interface/ceed-elemrestriction.c index 52bd36bf91..7492990b16 100644 --- a/interface/ceed-elemrestriction.c +++ b/interface/ceed-elemrestriction.c @@ -380,7 +380,6 @@ int CeedElemRestrictionCreate(Ceed ceed, CeedInt num_elem, CeedInt elem_size, [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 scale An scalar value that scales the dofs in assembly. @param[out] rstr Address of the variable where the newly created CeedElemRestriction will be stored @@ -408,9 +407,9 @@ int CeedElemRestrictionCreateOriented(Ceed ceed, CeedInt num_elem, // 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); + num_comp, comp_stride, l_size, + mem_type, copy_mode, offsets, + orient, rstr); CeedChk(ierr); return CEED_ERROR_SUCCESS; } @@ -425,9 +424,8 @@ int CeedElemRestrictionCreateOriented(Ceed ceed, CeedInt num_elem, (*rstr)->l_size = l_size; (*rstr)->num_blk = num_elem; (*rstr)->blk_size = 1; - ierr = ceed->ElemRestrictionCreateOriented(mem_type, copy_mode, offsets, orient, - *rstr); - CeedChk(ierr); + ierr = ceed->ElemRestrictionCreateOriented(mem_type, copy_mode, + offsets, orient, *rstr); CeedChk(ierr); return CEED_ERROR_SUCCESS; } From 61e7462c3a07c3c6a5221ea9e89d0c3cd40e49e5 Mon Sep 17 00:00:00 2001 From: Rezgar Shakeri <42816410+rezgarshakeri@users.noreply.github.com> Date: Wed, 2 Feb 2022 17:42:04 -0700 Subject: [PATCH 4/8] Update interface/ceed-elemrestriction.c Co-authored-by: Jeremy L Thompson --- interface/ceed-elemrestriction.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/interface/ceed-elemrestriction.c b/interface/ceed-elemrestriction.c index 7492990b16..a5cfa0ce59 100644 --- a/interface/ceed-elemrestriction.c +++ b/interface/ceed-elemrestriction.c @@ -403,7 +403,7 @@ int CeedElemRestrictionCreateOriented(Ceed ceed, CeedInt num_elem, if (!delegate) // LCOV_EXCL_START return CeedError(ceed, CEED_ERROR_UNSUPPORTED, - "Backend does not support ElemRestrictionCreateOriented"); + "Backend does not implement ElemRestrictionCreateOriented"); // LCOV_EXCL_STOP ierr = CeedElemRestrictionCreateOriented(delegate, num_elem, elem_size, From d4b88fd2eeefc92f28e90a50a0d82e97abcdcb7c Mon Sep 17 00:00:00 2001 From: rezgarshakeri Date: Thu, 3 Feb 2022 11:56:41 -0700 Subject: [PATCH 5/8] tests: deleted CeedVectorSetValue in restriction tests --- tests/t200-elemrestriction.c | 1 - tests/t201-elemrestriction.c | 1 - tests/t202-elemrestriction.c | 1 - tests/t203-elemrestriction.c | 1 - tests/t204-elemrestriction.c | 1 - tests/t205-elemrestriction.c | 1 - tests/t208-elemrestriction.c | 1 - tests/t213-elemrestriction.c | 1 - tests/t220-elemrestriction.c | 1 - 9 files changed, 9 deletions(-) diff --git a/tests/t200-elemrestriction.c b/tests/t200-elemrestriction.c index 15985962dd..f527b19a77 100644 --- 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 index 0981718f6a..3e50375499 100644 --- 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 index 0a11b126cb..e48a56e209 100644 --- 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 index 848ea8d2a1..de48c5890e 100644 --- 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 index e482cde984..ab37491408 100644 --- 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 index 099b7a0b24..48cd9cecd2 100644 --- 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 index 2fa483497e..d836eaf8ef 100644 --- 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 index 17bf192c9a..cae33ae689 100644 --- 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 index 7f8d3a8347..94367d1dc7 100644 --- a/tests/t220-elemrestriction.c +++ b/tests/t220-elemrestriction.c @@ -31,7 +31,6 @@ int main(int argc, char **argv) { num_elem+1, CEED_MEM_HOST, CEED_USE_POINTER, ind, orient, &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); From c7745053b0774c228ac7429b6b3e5eca9baa354c Mon Sep 17 00:00:00 2001 From: Rezgar Shakeri <42816410+rezgarshakeri@users.noreply.github.com> Date: Thu, 3 Feb 2022 17:14:40 -0700 Subject: [PATCH 6/8] Update interface/ceed-elemrestriction.c Co-authored-by: Jeremy L Thompson --- interface/ceed-elemrestriction.c | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/interface/ceed-elemrestriction.c b/interface/ceed-elemrestriction.c index a5cfa0ce59..65d5a8431d 100644 --- a/interface/ceed-elemrestriction.c +++ b/interface/ceed-elemrestriction.c @@ -395,7 +395,7 @@ int CeedElemRestrictionCreateOriented(Ceed ceed, CeedInt num_elem, CeedElemRestriction *rstr) { int ierr; - if (!ceed->ElemRestrictionCreate) { + if (!ceed->ElemRestrictionCreateOriented) { Ceed delegate; ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); CeedChk(ierr); @@ -407,8 +407,7 @@ int CeedElemRestrictionCreateOriented(Ceed ceed, CeedInt num_elem, // LCOV_EXCL_STOP ierr = CeedElemRestrictionCreateOriented(delegate, num_elem, elem_size, - num_comp, comp_stride, l_size, - mem_type, copy_mode, offsets, + num_comp, comp_stride, l_size, mem_type, copy_mode, offsets, orient, rstr); CeedChk(ierr); return CEED_ERROR_SUCCESS; } From b435c5a60b1833cef613d6281360673b9d9fef16 Mon Sep 17 00:00:00 2001 From: rezgarshakeri Date: Thu, 3 Feb 2022 19:53:08 -0700 Subject: [PATCH 7/8] Added CeedElemRestrictionIsOriented function --- backends/ref/ceed-ref-restriction.c | 52 ++++++++++++++++++++--------- include/ceed-impl.h | 1 + include/ceed/backend.h | 2 ++ interface/ceed-elemrestriction.c | 20 +++++++++++ 4 files changed, 60 insertions(+), 15 deletions(-) diff --git a/backends/ref/ceed-ref-restriction.c b/backends/ref/ceed-ref-restriction.c index 23e82e2a15..e9e70645ab 100644 --- 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 @@ -85,14 +87,24 @@ static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r, // Offsets provided, standard or blocked restriction // vv has shape [elem_size, num_comp, num_elem], row-major // uu has shape [nnodes, num_comp] - for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size) - CeedPragmaSIMD - for (CeedInt k = 0; k < num_comp; k++) + if (!is_oriented) { + for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size) + CeedPragmaSIMD + for (CeedInt k = 0; k < num_comp; k++) + 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]; + } else { + for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size) 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] * - (impl->orient && impl->orient[i+elem_size*e] ? -1. : 1.); + for (CeedInt k = 0; k < num_comp; k++) + 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] * + (impl->orient && impl->orient[i+elem_size*e] ? -1. : 1.); + } } } else { // Restriction from E-vector to L-vector @@ -132,14 +144,24 @@ static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r, // Offsets provided, standard or blocked restriction // uu has shape [elem_size, num_comp, num_elem] // vv has shape [nnodes, num_comp] - for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size) - for (CeedInt k = 0; k < num_comp; k++) - for (CeedInt i = 0; i < elem_size*blk_size; i+=blk_size) - // 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] * - (impl->orient && impl->orient[j+e*elem_size] ? -1. : 1.); + if (!is_oriented) { + for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size) + for (CeedInt k = 0; k < num_comp; k++) + for (CeedInt i = 0; i < elem_size*blk_size; i+=blk_size) + // 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]; + } else { + for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size) + for (CeedInt k = 0; k < num_comp; k++) + for (CeedInt i = 0; i < elem_size*blk_size; i+=blk_size) + // 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] * + (impl->orient && impl->orient[j+e*elem_size] ? -1. : 1.); + } } } ierr = CeedVectorRestoreArrayRead(u, &uu); CeedChkBackend(ierr); diff --git a/include/ceed-impl.h b/include/ceed-impl.h index 2b8ccefba0..510e523ec8 100644 --- a/include/ceed-impl.h +++ b/include/ceed-impl.h @@ -180,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 index 73f5db242d..78f15ae9c4 100644 --- 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/interface/ceed-elemrestriction.c b/interface/ceed-elemrestriction.c index 65d5a8431d..a18221faf4 100644 --- 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,6 +367,7 @@ 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; @@ -423,6 +439,7 @@ int CeedElemRestrictionCreateOriented(Ceed ceed, CeedInt num_elem, (*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; @@ -485,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]; @@ -571,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) { @@ -636,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]; From 000294e3464a97bfeded1398360f60f8eec27728 Mon Sep 17 00:00:00 2001 From: rezgarshakeri Date: Thu, 3 Feb 2022 20:32:20 -0700 Subject: [PATCH 8/8] updated ceed-ref-restriction.c --- backends/ref/ceed-ref-restriction.c | 50 +++++++++-------------------- 1 file changed, 15 insertions(+), 35 deletions(-) diff --git a/backends/ref/ceed-ref-restriction.c b/backends/ref/ceed-ref-restriction.c index e9e70645ab..fff1236672 100644 --- a/backends/ref/ceed-ref-restriction.c +++ b/backends/ref/ceed-ref-restriction.c @@ -87,24 +87,14 @@ static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r, // Offsets provided, standard or blocked restriction // vv has shape [elem_size, num_comp, num_elem], row-major // uu has shape [nnodes, num_comp] - if (!is_oriented) { - for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size) - CeedPragmaSIMD - for (CeedInt k = 0; k < num_comp; k++) - 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]; - } else { - for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size) + for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size) + CeedPragmaSIMD + for (CeedInt k = 0; k < num_comp; k++) CeedPragmaSIMD - for (CeedInt k = 0; k < num_comp; k++) - 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] * - (impl->orient && impl->orient[i+elem_size*e] ? -1. : 1.); - } + 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] * + (is_oriented && impl->orient[i+elem_size*e] ? -1. : 1.); } } else { // Restriction from E-vector to L-vector @@ -144,24 +134,14 @@ static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r, // Offsets provided, standard or blocked restriction // uu has shape [elem_size, num_comp, num_elem] // vv has shape [nnodes, num_comp] - if (!is_oriented) { - for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size) - for (CeedInt k = 0; k < num_comp; k++) - for (CeedInt i = 0; i < elem_size*blk_size; i+=blk_size) - // 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]; - } else { - for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size) - for (CeedInt k = 0; k < num_comp; k++) - for (CeedInt i = 0; i < elem_size*blk_size; i+=blk_size) - // 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] * - (impl->orient && impl->orient[j+e*elem_size] ? -1. : 1.); - } + for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size) + for (CeedInt k = 0; k < num_comp; k++) + for (CeedInt i = 0; i < elem_size*blk_size; i+=blk_size) + // 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] * + (is_oriented && impl->orient[j+e*elem_size] ? -1. : 1.); } } ierr = CeedVectorRestoreArrayRead(u, &uu); CeedChkBackend(ierr);