xref: /libCEED/rust/libceed-sys/c-src/backends/ref/ceed-ref-restriction.c (revision 60224bc541fd55b59daf56d957e93315ca03a274)
1 // Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC.
2 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707.
3 // All Rights reserved. See files LICENSE and NOTICE for details.
4 //
5 // This file is part of CEED, a collection of benchmarks, miniapps, software
6 // libraries and APIs for efficient high-order finite element and spectral
7 // element discretizations for exascale applications. For more information and
8 // source code availability see http://github.com/ceed.
9 //
10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
11 // a collaborative effort of two U.S. Department of Energy organizations (Office
12 // of Science and the National Nuclear Security Administration) responsible for
13 // the planning and preparation of a capable exascale ecosystem, including
14 // software, applications, hardware, advanced system engineering and early
15 // testbed platforms, in support of the nation's exascale computing imperative.
16 
17 #include <ceed/ceed.h>
18 #include <ceed/backend.h>
19 #include <stdbool.h>
20 #include <string.h>
21 #include "ceed-ref.h"
22 
23 //------------------------------------------------------------------------------
24 // Core ElemRestriction Apply Code
25 //------------------------------------------------------------------------------
26 static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r,
27     const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride,
28     CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u,
29     CeedVector v, CeedRequest *request) {
30   int ierr;
31   CeedElemRestriction_Ref *impl;
32   ierr = CeedElemRestrictionGetData(r, &impl); CeedChkBackend(ierr);
33   const CeedScalar *uu;
34   CeedScalar *vv;
35   CeedInt num_elem, elem_size, v_offset;
36   ierr = CeedElemRestrictionGetNumElements(r, &num_elem); CeedChkBackend(ierr);
37   ierr = CeedElemRestrictionGetElementSize(r, &elem_size); CeedChkBackend(ierr);
38   v_offset = start*blk_size*elem_size*num_comp;
39 
40   bool is_oriented;
41   ierr = CeedElemRestrictionIsOriented(r, &is_oriented); CeedChkBackend(ierr);
42   ierr = CeedVectorGetArrayRead(u, CEED_MEM_HOST, &uu); CeedChkBackend(ierr);
43   if (t_mode == CEED_TRANSPOSE) {
44     // Sum into for transpose mode, e-vec to l-vec
45     ierr = CeedVectorGetArray(v, CEED_MEM_HOST, &vv); CeedChkBackend(ierr);
46   } else {
47     // Overwrite for notranspose mode, l-vec to e-vec
48     ierr = CeedVectorGetArrayWrite(v, CEED_MEM_HOST, &vv); CeedChkBackend(ierr);
49   }
50   // Restriction from L-vector to E-vector
51   // Perform: v = r * u
52   if (t_mode == CEED_NOTRANSPOSE) {
53     // No offsets provided, Identity Restriction
54     if (!impl->offsets) {
55       bool has_backend_strides;
56       ierr = CeedElemRestrictionHasBackendStrides(r, &has_backend_strides);
57       CeedChkBackend(ierr);
58       if (has_backend_strides) {
59         // CPU backend strides are {1, elem_size, elem_size*num_comp}
60         // This if branch is left separate to allow better inlining
61         for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size)
62           CeedPragmaSIMD
63           for (CeedInt k = 0; k < num_comp; k++)
64             CeedPragmaSIMD
65             for (CeedInt n = 0; n < elem_size; n++)
66               CeedPragmaSIMD
67               for (CeedInt j = 0; j < blk_size; j++)
68                 vv[e*elem_size*num_comp + (k*elem_size+n)*blk_size + j - v_offset]
69                   = uu[n + k*elem_size +
70                          CeedIntMin(e+j, num_elem-1)*elem_size*num_comp];
71       } else {
72         // User provided strides
73         CeedInt strides[3];
74         ierr = CeedElemRestrictionGetStrides(r, &strides); CeedChkBackend(ierr);
75         for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size)
76           CeedPragmaSIMD
77           for (CeedInt k = 0; k < num_comp; k++)
78             CeedPragmaSIMD
79             for (CeedInt n = 0; n < elem_size; n++)
80               CeedPragmaSIMD
81               for (CeedInt j = 0; j < blk_size; j++)
82                 vv[e*elem_size*num_comp + (k*elem_size+n)*blk_size + j - v_offset]
83                   = uu[n*strides[0] + k*strides[1] +
84                                     CeedIntMin(e+j, num_elem-1)*strides[2]];
85       }
86     } else {
87       // Offsets provided, standard or blocked restriction
88       // vv has shape [elem_size, num_comp, num_elem], row-major
89       // uu has shape [nnodes, num_comp]
90       for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size)
91         CeedPragmaSIMD
92         for (CeedInt k = 0; k < num_comp; k++)
93           CeedPragmaSIMD
94           for (CeedInt i = 0; i < elem_size*blk_size; i++)
95             vv[elem_size*(k*blk_size+num_comp*e) + i - v_offset]
96               = uu[impl->offsets[i+elem_size*e] + k*comp_stride] *
97                 (is_oriented && impl->orient[i+elem_size*e] ? -1. : 1.);
98     }
99   } else {
100     // Restriction from E-vector to L-vector
101     // Performing v += r^T * u
102     // No offsets provided, Identity Restriction
103     if (!impl->offsets) {
104       bool has_backend_strides;
105       ierr = CeedElemRestrictionHasBackendStrides(r, &has_backend_strides);
106       CeedChkBackend(ierr);
107       if (has_backend_strides) {
108         // CPU backend strides are {1, elem_size, elem_size*num_comp}
109         // This if brach is left separate to allow better inlining
110         for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size)
111           CeedPragmaSIMD
112           for (CeedInt k = 0; k < num_comp; k++)
113             CeedPragmaSIMD
114             for (CeedInt n = 0; n < elem_size; n++)
115               CeedPragmaSIMD
116               for (CeedInt j = 0; j < CeedIntMin(blk_size, num_elem-e); j++)
117                 vv[n + k*elem_size + (e+j)*elem_size*num_comp]
118                 += uu[e*elem_size*num_comp + (k*elem_size+n)*blk_size + j - v_offset];
119       } else {
120         // User provided strides
121         CeedInt strides[3];
122         ierr = CeedElemRestrictionGetStrides(r, &strides); CeedChkBackend(ierr);
123         for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size)
124           CeedPragmaSIMD
125           for (CeedInt k = 0; k < num_comp; k++)
126             CeedPragmaSIMD
127             for (CeedInt n = 0; n < elem_size; n++)
128               CeedPragmaSIMD
129               for (CeedInt j = 0; j < CeedIntMin(blk_size, num_elem-e); j++)
130                 vv[n*strides[0] + k*strides[1] + (e+j)*strides[2]]
131                 += uu[e*elem_size*num_comp + (k*elem_size+n)*blk_size + j - v_offset];
132       }
133     } else {
134       // Offsets provided, standard or blocked restriction
135       // uu has shape [elem_size, num_comp, num_elem]
136       // vv has shape [nnodes, num_comp]
137       for (CeedInt e = start*blk_size; e < stop*blk_size; e+=blk_size)
138         for (CeedInt k = 0; k < num_comp; k++)
139           for (CeedInt i = 0; i < elem_size*blk_size; i+=blk_size)
140             // Iteration bound set to discard padding elements
141             for (CeedInt j = i; j < i+CeedIntMin(blk_size, num_elem-e); j++)
142               vv[impl->offsets[j+e*elem_size] + k*comp_stride]
143               += uu[elem_size*(k*blk_size+num_comp*e) + j - v_offset] *
144                  (is_oriented && impl->orient[j+e*elem_size] ? -1. : 1.);
145     }
146   }
147   ierr = CeedVectorRestoreArrayRead(u, &uu); CeedChkBackend(ierr);
148   ierr = CeedVectorRestoreArray(v, &vv); CeedChkBackend(ierr);
149   if (request != CEED_REQUEST_IMMEDIATE && request != CEED_REQUEST_ORDERED)
150     *request = NULL;
151   return CEED_ERROR_SUCCESS;
152 }
153 
154 //------------------------------------------------------------------------------
155 // ElemRestriction Apply - Common Sizes
156 //------------------------------------------------------------------------------
157 static int CeedElemRestrictionApply_Ref_110(CeedElemRestriction r,
158     const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride,
159     CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u,
160     CeedVector v, CeedRequest *request) {
161   return CeedElemRestrictionApply_Ref_Core(r, 1, 1, comp_stride, start, stop,
162          t_mode, u, v, request);
163 }
164 
165 static int CeedElemRestrictionApply_Ref_111(CeedElemRestriction r,
166     const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride,
167     CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u,
168     CeedVector v, CeedRequest *request) {
169   return CeedElemRestrictionApply_Ref_Core(r, 1, 1, 1, start, stop, t_mode,
170          u, v, request);
171 }
172 
173 static int CeedElemRestrictionApply_Ref_180(CeedElemRestriction r,
174     const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride,
175     CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u,
176     CeedVector v, CeedRequest *request) {
177   return CeedElemRestrictionApply_Ref_Core(r, 1, 8, comp_stride, start, stop,
178          t_mode, u, v, request);
179 }
180 
181 static int CeedElemRestrictionApply_Ref_181(CeedElemRestriction r,
182     const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride,
183     CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u,
184     CeedVector v, CeedRequest *request) {
185   return CeedElemRestrictionApply_Ref_Core(r, 1, 8, 1, start, stop, t_mode,
186          u, v, request);
187 }
188 
189 static int CeedElemRestrictionApply_Ref_310(CeedElemRestriction r,
190     const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride,
191     CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u,
192     CeedVector v, CeedRequest *request) {
193   return CeedElemRestrictionApply_Ref_Core(r, 3, 1, comp_stride, start, stop,
194          t_mode, u, v, request);
195 }
196 
197 static int CeedElemRestrictionApply_Ref_311(CeedElemRestriction r,
198     const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride,
199     CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u,
200     CeedVector v, CeedRequest *request) {
201   return CeedElemRestrictionApply_Ref_Core(r, 3, 1, 1, start, stop, t_mode,
202          u, v, request);
203 }
204 
205 static int CeedElemRestrictionApply_Ref_380(CeedElemRestriction r,
206     const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride,
207     CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u,
208     CeedVector v, CeedRequest *request) {
209   return CeedElemRestrictionApply_Ref_Core(r, 3, 8, comp_stride, start, stop,
210          t_mode, u, v, request);
211 }
212 
213 static int CeedElemRestrictionApply_Ref_381(CeedElemRestriction r,
214     const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride,
215     CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u,
216     CeedVector v, CeedRequest *request) {
217   return CeedElemRestrictionApply_Ref_Core(r, 3, 8, 1, start, stop, t_mode,
218          u, v, request);
219 }
220 
221 // LCOV_EXCL_START
222 static int CeedElemRestrictionApply_Ref_510(CeedElemRestriction r,
223     const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride,
224     CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u,
225     CeedVector v, CeedRequest *request) {
226   return CeedElemRestrictionApply_Ref_Core(r, 5, 1, comp_stride, start, stop,
227          t_mode, u, v, request);
228 }
229 // LCOV_EXCL_STOP
230 
231 static int CeedElemRestrictionApply_Ref_511(CeedElemRestriction r,
232     const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride,
233     CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u,
234     CeedVector v, CeedRequest *request) {
235   return CeedElemRestrictionApply_Ref_Core(r, 5, 1, 1, start, stop, t_mode,
236          u, v, request);
237 }
238 
239 // LCOV_EXCL_START
240 static int CeedElemRestrictionApply_Ref_580(CeedElemRestriction r,
241     const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride,
242     CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u,
243     CeedVector v, CeedRequest *request) {
244   return CeedElemRestrictionApply_Ref_Core(r, 5, 8, comp_stride, start, stop,
245          t_mode, u, v, request);
246 }
247 // LCOV_EXCL_STOP
248 
249 static int CeedElemRestrictionApply_Ref_581(CeedElemRestriction r,
250     const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride,
251     CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u,
252     CeedVector v, CeedRequest *request) {
253   return CeedElemRestrictionApply_Ref_Core(r, 5, 8, 1, start, stop, t_mode,
254          u, v, request);
255 }
256 
257 //------------------------------------------------------------------------------
258 // ElemRestriction Apply
259 //------------------------------------------------------------------------------
260 static int CeedElemRestrictionApply_Ref(CeedElemRestriction r,
261                                         CeedTransposeMode t_mode, CeedVector u,
262                                         CeedVector v, CeedRequest *request) {
263   int ierr;
264   CeedInt num_blk, blk_size, num_comp, comp_stride;
265   ierr = CeedElemRestrictionGetNumBlocks(r, &num_blk); CeedChkBackend(ierr);
266   ierr = CeedElemRestrictionGetBlockSize(r, &blk_size); CeedChkBackend(ierr);
267   ierr = CeedElemRestrictionGetNumComponents(r, &num_comp); CeedChkBackend(ierr);
268   ierr = CeedElemRestrictionGetCompStride(r, &comp_stride); CeedChkBackend(ierr);
269   CeedElemRestriction_Ref *impl;
270   ierr = CeedElemRestrictionGetData(r, &impl); CeedChkBackend(ierr);
271 
272   return impl->Apply(r, num_comp, blk_size, comp_stride, 0, num_blk, t_mode, u, v,
273                      request);
274 }
275 
276 //------------------------------------------------------------------------------
277 // ElemRestriction Apply Block
278 //------------------------------------------------------------------------------
279 static int CeedElemRestrictionApplyBlock_Ref(CeedElemRestriction r,
280     CeedInt block, CeedTransposeMode t_mode, CeedVector u, CeedVector v,
281     CeedRequest *request) {
282   int ierr;
283   CeedInt blk_size, num_comp, comp_stride;
284   ierr = CeedElemRestrictionGetBlockSize(r, &blk_size); CeedChkBackend(ierr);
285   ierr = CeedElemRestrictionGetNumComponents(r, &num_comp); CeedChkBackend(ierr);
286   ierr = CeedElemRestrictionGetCompStride(r, &comp_stride); CeedChkBackend(ierr);
287   CeedElemRestriction_Ref *impl;
288   ierr = CeedElemRestrictionGetData(r, &impl); CeedChkBackend(ierr);
289 
290   return impl->Apply(r, num_comp, blk_size, comp_stride, block, block+1, t_mode,
291                      u, v, request);
292 }
293 
294 //------------------------------------------------------------------------------
295 // ElemRestriction Get Offsets
296 //------------------------------------------------------------------------------
297 static int CeedElemRestrictionGetOffsets_Ref(CeedElemRestriction rstr,
298     CeedMemType mem_type, const CeedInt **offsets) {
299   int ierr;
300   CeedElemRestriction_Ref *impl;
301   ierr = CeedElemRestrictionGetData(rstr, &impl); CeedChkBackend(ierr);
302   Ceed ceed;
303   ierr = CeedElemRestrictionGetCeed(rstr, &ceed); CeedChkBackend(ierr);
304 
305   if (mem_type != CEED_MEM_HOST)
306     // LCOV_EXCL_START
307     return CeedError(ceed, CEED_ERROR_BACKEND, "Can only provide to HOST memory");
308   // LCOV_EXCL_STOP
309 
310   *offsets = impl->offsets;
311   return CEED_ERROR_SUCCESS;
312 }
313 
314 //------------------------------------------------------------------------------
315 // ElemRestriction Destroy
316 //------------------------------------------------------------------------------
317 static int CeedElemRestrictionDestroy_Ref(CeedElemRestriction r) {
318   int ierr;
319   CeedElemRestriction_Ref *impl;
320   ierr = CeedElemRestrictionGetData(r, &impl); CeedChkBackend(ierr);
321 
322   ierr = CeedFree(&impl->offsets_allocated); CeedChkBackend(ierr);
323   ierr = CeedFree(&impl); CeedChkBackend(ierr);
324   return CEED_ERROR_SUCCESS;
325 }
326 
327 //------------------------------------------------------------------------------
328 // ElemRestriction Create
329 //------------------------------------------------------------------------------
330 int CeedElemRestrictionCreate_Ref(CeedMemType mem_type, CeedCopyMode copy_mode,
331                                   const CeedInt *offsets,
332                                   CeedElemRestriction r) {
333   int ierr;
334   CeedElemRestriction_Ref *impl;
335   CeedInt num_elem, elem_size, num_blk, blk_size, num_comp, comp_stride;
336   ierr = CeedElemRestrictionGetNumElements(r, &num_elem); CeedChkBackend(ierr);
337   ierr = CeedElemRestrictionGetElementSize(r, &elem_size); CeedChkBackend(ierr);
338   ierr = CeedElemRestrictionGetNumBlocks(r, &num_blk); CeedChkBackend(ierr);
339   ierr = CeedElemRestrictionGetBlockSize(r, &blk_size); CeedChkBackend(ierr);
340   ierr = CeedElemRestrictionGetNumComponents(r, &num_comp); CeedChkBackend(ierr);
341   ierr = CeedElemRestrictionGetCompStride(r, &comp_stride); CeedChkBackend(ierr);
342   Ceed ceed;
343   ierr = CeedElemRestrictionGetCeed(r, &ceed); CeedChkBackend(ierr);
344 
345   if (mem_type != CEED_MEM_HOST)
346     // LCOV_EXCL_START
347     return CeedError(ceed, CEED_ERROR_BACKEND, "Only MemType = HOST supported");
348   // LCOV_EXCL_STOP
349   ierr = CeedCalloc(1, &impl); CeedChkBackend(ierr);
350 
351   // Offsets data
352   bool is_strided;
353   ierr = CeedElemRestrictionIsStrided(r, &is_strided); CeedChkBackend(ierr);
354   if (!is_strided) {
355     // Check indices for ref or memcheck backends
356     Ceed parent_ceed = ceed, curr_ceed = NULL;
357     while (parent_ceed != curr_ceed) {
358       curr_ceed = parent_ceed;
359       ierr = CeedGetParent(curr_ceed, &parent_ceed); CeedChkBackend(ierr);
360     }
361     const char *resource;
362     ierr = CeedGetResource(parent_ceed, &resource); CeedChkBackend(ierr);
363     if (!strcmp(resource, "/cpu/self/ref/serial") ||
364         !strcmp(resource, "/cpu/self/ref/blocked") ||
365         !strcmp(resource, "/cpu/self/memcheck/serial") ||
366         !strcmp(resource, "/cpu/self/memcheck/blocked")) {
367       CeedSize l_size;
368       ierr = CeedElemRestrictionGetLVectorSize(r, &l_size); CeedChkBackend(ierr);
369 
370       for (CeedInt i = 0; i < num_elem*elem_size; i++)
371         if (offsets[i] < 0 || l_size <= offsets[i] + (num_comp - 1) * comp_stride)
372           // LCOV_EXCL_START
373           return CeedError(ceed, CEED_ERROR_BACKEND,
374                            "Restriction offset %d (%d) out of range "
375                            "[0, %d]", i, offsets[i], l_size);
376       // LCOV_EXCL_STOP
377     }
378 
379     // Copy data
380     switch (copy_mode) {
381     case CEED_COPY_VALUES:
382       ierr = CeedMalloc(num_elem*elem_size, &impl->offsets_allocated);
383       CeedChkBackend(ierr);
384       memcpy(impl->offsets_allocated, offsets,
385              num_elem * elem_size * sizeof(offsets[0]));
386       impl->offsets = impl->offsets_allocated;
387       break;
388     case CEED_OWN_POINTER:
389       impl->offsets_allocated = (CeedInt *)offsets;
390       impl->offsets = impl->offsets_allocated;
391       break;
392     case CEED_USE_POINTER:
393       impl->offsets = offsets;
394     }
395   }
396 
397   ierr = CeedElemRestrictionSetData(r, impl); CeedChkBackend(ierr);
398   CeedInt layout[3] = {1, elem_size, elem_size*num_comp};
399   ierr = CeedElemRestrictionSetELayout(r, layout); CeedChkBackend(ierr);
400   ierr = CeedSetBackendFunction(ceed, "ElemRestriction", r, "Apply",
401                                 CeedElemRestrictionApply_Ref); CeedChkBackend(ierr);
402   ierr = CeedSetBackendFunction(ceed, "ElemRestriction", r, "ApplyBlock",
403                                 CeedElemRestrictionApplyBlock_Ref);
404   CeedChkBackend(ierr);
405   ierr = CeedSetBackendFunction(ceed, "ElemRestriction", r, "GetOffsets",
406                                 CeedElemRestrictionGetOffsets_Ref);
407   CeedChkBackend(ierr);
408   ierr = CeedSetBackendFunction(ceed, "ElemRestriction", r, "Destroy",
409                                 CeedElemRestrictionDestroy_Ref); CeedChkBackend(ierr);
410 
411   // Set apply function based upon num_comp, blk_size, and comp_stride
412   CeedInt idx = -1;
413   if (blk_size < 10)
414     idx = 100*num_comp + 10*blk_size + (comp_stride == 1);
415   switch (idx) {
416   case 110:
417     impl->Apply = CeedElemRestrictionApply_Ref_110;
418     break;
419   case 111:
420     impl->Apply = CeedElemRestrictionApply_Ref_111;
421     break;
422   case 180:
423     impl->Apply = CeedElemRestrictionApply_Ref_180;
424     break;
425   case 181:
426     impl->Apply = CeedElemRestrictionApply_Ref_181;
427     break;
428   case 310:
429     impl->Apply = CeedElemRestrictionApply_Ref_310;
430     break;
431   case 311:
432     impl->Apply = CeedElemRestrictionApply_Ref_311;
433     break;
434   case 380:
435     impl->Apply = CeedElemRestrictionApply_Ref_380;
436     break;
437   case 381:
438     impl->Apply = CeedElemRestrictionApply_Ref_381;
439     break;
440   // LCOV_EXCL_START
441   case 510:
442     impl->Apply = CeedElemRestrictionApply_Ref_510;
443     break;
444   // LCOV_EXCL_STOP
445   case 511:
446     impl->Apply = CeedElemRestrictionApply_Ref_511;
447     break;
448   // LCOV_EXCL_START
449   case 580:
450     impl->Apply = CeedElemRestrictionApply_Ref_580;
451     break;
452   // LCOV_EXCL_STOP
453   case 581:
454     impl->Apply = CeedElemRestrictionApply_Ref_581;
455     break;
456   default:
457     impl->Apply = CeedElemRestrictionApply_Ref_Core;
458     break;
459   }
460 
461   return CEED_ERROR_SUCCESS;
462 }
463 
464 //------------------------------------------------------------------------------
465 // ElemRestriction Create Oriented
466 //------------------------------------------------------------------------------
467 int CeedElemRestrictionCreateOriented_Ref(CeedMemType mem_type,
468     CeedCopyMode copy_mode,
469     const CeedInt *offsets, const bool *orient,
470     CeedElemRestriction r) {
471   int ierr;
472   CeedElemRestriction_Ref *impl;
473   CeedInt num_elem, elem_size;
474   // Set up for normal restriction with explicit offsets. This sets up dispatch to
475   // CeedElemRestrictionApply_Ref_* and manages the impl->offsets array copy/allocation.
476   ierr = CeedElemRestrictionCreate_Ref(mem_type, copy_mode, offsets, r);
477   CeedChkBackend(ierr);
478 
479   ierr = CeedElemRestrictionGetData(r, &impl); CeedChkBackend(ierr);
480   ierr = CeedElemRestrictionGetNumElements(r, &num_elem); CeedChkBackend(ierr);
481   ierr = CeedElemRestrictionGetElementSize(r, &elem_size); CeedChkBackend(ierr);
482   switch (copy_mode) {
483   case CEED_COPY_VALUES:
484     ierr = CeedMalloc(num_elem * elem_size, &impl->orient_allocated);
485     CeedChkBackend(ierr);
486     memcpy(impl->orient_allocated, orient,
487            num_elem * elem_size * sizeof(orient[0]));
488     impl->orient = impl->orient_allocated;
489     break;
490   case CEED_OWN_POINTER:
491     impl->orient_allocated = (bool *)orient;
492     impl->orient = impl->orient_allocated;
493     break;
494   case CEED_USE_POINTER:
495     impl->orient = orient;
496   }
497   return CEED_ERROR_SUCCESS;
498 }
499 //------------------------------------------------------------------------------
500