xref: /libCEED/backends/ref/ceed-ref-restriction.c (revision 2c7e74134068de5221afb7e7255d8033016bfd12)
1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED:  http://github.com/ceed
7 
8 #include <ceed.h>
9 #include <ceed/backend.h>
10 #include <stdbool.h>
11 #include <stdlib.h>
12 #include <string.h>
13 
14 #include "ceed-ref.h"
15 
16 //------------------------------------------------------------------------------
17 // Core ElemRestriction Apply Code
18 //------------------------------------------------------------------------------
19 static inline int CeedElemRestrictionApplyStridedNoTranspose_Ref_Core(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size,
20                                                                       CeedInt start, CeedInt stop, CeedInt num_elem, CeedInt elem_size,
21                                                                       CeedInt v_offset, const CeedScalar *uu, CeedScalar *vv) {
22   // No offsets provided, identity restriction
23   bool has_backend_strides;
24 
25   CeedCallBackend(CeedElemRestrictionHasBackendStrides(rstr, &has_backend_strides));
26   if (has_backend_strides) {
27     // CPU backend strides are {1, elem_size, elem_size*num_comp}
28     // This if branch is left separate to allow better inlining
29     for (CeedInt e = start * block_size; e < stop * block_size; e += block_size) {
30       CeedPragmaSIMD for (CeedInt k = 0; k < num_comp; k++) {
31         CeedPragmaSIMD for (CeedInt n = 0; n < elem_size; n++) {
32           CeedPragmaSIMD for (CeedInt j = 0; j < block_size; j++) {
33             vv[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] =
34                 uu[n + k * elem_size + CeedIntMin(e + j, num_elem - 1) * elem_size * num_comp];
35           }
36         }
37       }
38     }
39   } else {
40     // User provided strides
41     CeedInt strides[3];
42 
43     CeedCallBackend(CeedElemRestrictionGetStrides(rstr, &strides));
44     for (CeedInt e = start * block_size; e < stop * block_size; e += block_size) {
45       CeedPragmaSIMD for (CeedInt k = 0; k < num_comp; k++) {
46         CeedPragmaSIMD for (CeedInt n = 0; n < elem_size; n++) {
47           CeedPragmaSIMD for (CeedInt j = 0; j < block_size; j++) {
48             vv[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] =
49                 uu[n * strides[0] + k * strides[1] + CeedIntMin(e + j, num_elem - 1) * strides[2]];
50           }
51         }
52       }
53     }
54   }
55   return CEED_ERROR_SUCCESS;
56 }
57 
58 static inline int CeedElemRestrictionApplyStandardNoTranspose_Ref_Core(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size,
59                                                                        const CeedInt comp_stride, CeedInt start, CeedInt stop, CeedInt num_elem,
60                                                                        CeedInt elem_size, CeedInt v_offset, const CeedScalar *uu, CeedScalar *vv) {
61   // Default restriction with offsets
62   CeedElemRestriction_Ref *impl;
63 
64   CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl));
65   for (CeedInt e = start * block_size; e < stop * block_size; e += block_size) {
66     CeedPragmaSIMD for (CeedInt k = 0; k < num_comp; k++) {
67       CeedPragmaSIMD for (CeedInt i = 0; i < elem_size * block_size; i++) {
68         vv[elem_size * (k * block_size + e * num_comp) + i - v_offset] = uu[impl->offsets[i + e * elem_size] + k * comp_stride];
69       }
70     }
71   }
72   return CEED_ERROR_SUCCESS;
73 }
74 
75 static inline int CeedElemRestrictionApplyOrientedNoTranspose_Ref_Core(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size,
76                                                                        const CeedInt comp_stride, CeedInt start, CeedInt stop, CeedInt num_elem,
77                                                                        CeedInt elem_size, CeedInt v_offset, const CeedScalar *uu, CeedScalar *vv) {
78   // Restriction with orientations
79   CeedElemRestriction_Ref *impl;
80 
81   CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl));
82   for (CeedInt e = start * block_size; e < stop * block_size; e += block_size) {
83     CeedPragmaSIMD for (CeedInt k = 0; k < num_comp; k++) {
84       CeedPragmaSIMD for (CeedInt i = 0; i < elem_size * block_size; i++) {
85         vv[elem_size * (k * block_size + e * num_comp) + i - v_offset] =
86             uu[impl->offsets[i + e * elem_size] + k * comp_stride] * (impl->orients[i + e * elem_size] ? -1.0 : 1.0);
87       }
88     }
89   }
90   return CEED_ERROR_SUCCESS;
91 }
92 
93 static inline int CeedElemRestrictionApplyCurlOrientedNoTranspose_Ref_Core(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size,
94                                                                            const CeedInt comp_stride, CeedInt start, CeedInt stop, CeedInt num_elem,
95                                                                            CeedInt elem_size, CeedInt v_offset, const CeedScalar *uu,
96                                                                            CeedScalar *vv) {
97   // Restriction with tridiagonal transformation
98   CeedElemRestriction_Ref *impl;
99 
100   CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl));
101   for (CeedInt e = start * block_size; e < stop * block_size; e += block_size) {
102     CeedPragmaSIMD for (CeedInt k = 0; k < num_comp; k++) {
103       CeedInt n = 0;
104       CeedPragmaSIMD for (CeedInt j = 0; j < block_size; j++) {
105         vv[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] =
106             uu[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] *
107                 impl->curl_orients[j + (3 * n + 1) * block_size + e * 3 * elem_size] +
108             uu[impl->offsets[j + (n + 1) * block_size + e * elem_size] + k * comp_stride] *
109                 impl->curl_orients[j + (3 * n + 2) * block_size + e * 3 * elem_size];
110       }
111       for (n = 1; n < elem_size - 1; n++) {
112         CeedPragmaSIMD for (CeedInt j = 0; j < block_size; j++) {
113           vv[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] =
114               uu[impl->offsets[j + (n - 1) * block_size + e * elem_size] + k * comp_stride] *
115                   impl->curl_orients[j + (3 * n + 0) * block_size + e * 3 * elem_size] +
116               uu[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] *
117                   impl->curl_orients[j + (3 * n + 1) * block_size + e * 3 * elem_size] +
118               uu[impl->offsets[j + (n + 1) * block_size + e * elem_size] + k * comp_stride] *
119                   impl->curl_orients[j + (3 * n + 2) * block_size + e * 3 * elem_size];
120         }
121       }
122       CeedPragmaSIMD for (CeedInt j = 0; j < block_size; j++) {
123         vv[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] =
124             uu[impl->offsets[j + (n - 1) * block_size + e * elem_size] + k * comp_stride] *
125                 impl->curl_orients[j + (3 * n + 0) * block_size + e * 3 * elem_size] +
126             uu[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] *
127                 impl->curl_orients[j + (3 * n + 1) * block_size + e * 3 * elem_size];
128       }
129     }
130   }
131   return CEED_ERROR_SUCCESS;
132 }
133 
134 static inline int CeedElemRestrictionApplyCurlOrientedUnsignedNoTranspose_Ref_Core(CeedElemRestriction rstr, const CeedInt num_comp,
135                                                                                    const CeedInt block_size, const CeedInt comp_stride, CeedInt start,
136                                                                                    CeedInt stop, CeedInt num_elem, CeedInt elem_size,
137                                                                                    CeedInt v_offset, const CeedScalar *uu, CeedScalar *vv) {
138   // Restriction with (unsigned) tridiagonal transformation
139   CeedElemRestriction_Ref *impl;
140 
141   CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl));
142   for (CeedInt e = start * block_size; e < stop * block_size; e += block_size) {
143     CeedPragmaSIMD for (CeedInt k = 0; k < num_comp; k++) {
144       CeedInt n = 0;
145 
146       CeedPragmaSIMD for (CeedInt j = 0; j < block_size; j++) {
147         vv[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] =
148             uu[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] *
149                 abs(impl->curl_orients[j + (3 * n + 1) * block_size + e * 3 * elem_size]) +
150             uu[impl->offsets[j + (n + 1) * block_size + e * elem_size] + k * comp_stride] *
151                 abs(impl->curl_orients[j + (3 * n + 2) * block_size + e * 3 * elem_size]);
152       }
153       for (n = 1; n < elem_size - 1; n++) {
154         CeedPragmaSIMD for (CeedInt j = 0; j < block_size; j++) {
155           vv[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] =
156               uu[impl->offsets[j + (n - 1) * block_size + e * elem_size] + k * comp_stride] *
157                   abs(impl->curl_orients[j + (3 * n + 0) * block_size + e * 3 * elem_size]) +
158               uu[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] *
159                   abs(impl->curl_orients[j + (3 * n + 1) * block_size + e * 3 * elem_size]) +
160               uu[impl->offsets[j + (n + 1) * block_size + e * elem_size] + k * comp_stride] *
161                   abs(impl->curl_orients[j + (3 * n + 2) * block_size + e * 3 * elem_size]);
162         }
163       }
164       CeedPragmaSIMD for (CeedInt j = 0; j < block_size; j++) {
165         vv[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] =
166             uu[impl->offsets[j + (n - 1) * block_size + e * elem_size] + k * comp_stride] *
167                 abs(impl->curl_orients[j + (3 * n + 0) * block_size + e * 3 * elem_size]) +
168             uu[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] *
169                 abs(impl->curl_orients[j + (3 * n + 1) * block_size + e * 3 * elem_size]);
170       }
171     }
172   }
173   return CEED_ERROR_SUCCESS;
174 }
175 
176 static inline int CeedElemRestrictionApplyStridedTranspose_Ref_Core(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size,
177                                                                     CeedInt start, CeedInt stop, CeedInt num_elem, CeedInt elem_size,
178                                                                     CeedInt v_offset, const CeedScalar *uu, CeedScalar *vv) {
179   // No offsets provided, identity restriction
180   bool has_backend_strides;
181 
182   CeedCallBackend(CeedElemRestrictionHasBackendStrides(rstr, &has_backend_strides));
183   if (has_backend_strides) {
184     // CPU backend strides are {1, elem_size, elem_size*num_comp}
185     // This if brach is left separate to allow better inlining
186     for (CeedInt e = start * block_size; e < stop * block_size; e += block_size) {
187       CeedPragmaSIMD for (CeedInt k = 0; k < num_comp; k++) {
188         CeedPragmaSIMD for (CeedInt n = 0; n < elem_size; n++) {
189           CeedPragmaSIMD for (CeedInt j = 0; j < CeedIntMin(block_size, num_elem - e); j++) {
190             vv[n + k * elem_size + (e + j) * elem_size * num_comp] += uu[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset];
191           }
192         }
193       }
194     }
195   } else {
196     // User provided strides
197     CeedInt strides[3];
198 
199     CeedCallBackend(CeedElemRestrictionGetStrides(rstr, &strides));
200     for (CeedInt e = start * block_size; e < stop * block_size; e += block_size) {
201       CeedPragmaSIMD for (CeedInt k = 0; k < num_comp; k++) {
202         CeedPragmaSIMD for (CeedInt n = 0; n < elem_size; n++) {
203           CeedPragmaSIMD for (CeedInt j = 0; j < CeedIntMin(block_size, num_elem - e); j++) {
204             vv[n * strides[0] + k * strides[1] + (e + j) * strides[2]] +=
205                 uu[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset];
206           }
207         }
208       }
209     }
210   }
211   return CEED_ERROR_SUCCESS;
212 }
213 
214 static inline int CeedElemRestrictionApplyStandardTranspose_Ref_Core(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size,
215                                                                      const CeedInt comp_stride, CeedInt start, CeedInt stop, CeedInt num_elem,
216                                                                      CeedInt elem_size, CeedInt v_offset, const CeedScalar *uu, CeedScalar *vv) {
217   // Default restriction with offsets
218   CeedElemRestriction_Ref *impl;
219 
220   CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl));
221   for (CeedInt e = start * block_size; e < stop * block_size; e += block_size) {
222     for (CeedInt k = 0; k < num_comp; k++) {
223       for (CeedInt i = 0; i < elem_size * block_size; i += block_size) {
224         // Iteration bound set to discard padding elements
225         for (CeedInt j = i; j < i + CeedIntMin(block_size, num_elem - e); j++) {
226           vv[impl->offsets[j + e * elem_size] + k * comp_stride] += uu[elem_size * (k * block_size + e * num_comp) + j - v_offset];
227         }
228       }
229     }
230   }
231   return CEED_ERROR_SUCCESS;
232 }
233 
234 static inline int CeedElemRestrictionApplyOrientedTranspose_Ref_Core(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size,
235                                                                      const CeedInt comp_stride, CeedInt start, CeedInt stop, CeedInt num_elem,
236                                                                      CeedInt elem_size, CeedInt v_offset, const CeedScalar *uu, CeedScalar *vv) {
237   // Restriction with orientations
238   CeedElemRestriction_Ref *impl;
239 
240   CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl));
241   for (CeedInt e = start * block_size; e < stop * block_size; e += block_size) {
242     for (CeedInt k = 0; k < num_comp; k++) {
243       for (CeedInt i = 0; i < elem_size * block_size; i += block_size) {
244         // Iteration bound set to discard padding elements
245         for (CeedInt j = i; j < i + CeedIntMin(block_size, num_elem - e); j++) {
246           vv[impl->offsets[j + e * elem_size] + k * comp_stride] +=
247               uu[elem_size * (k * block_size + e * num_comp) + j - v_offset] * (impl->orients[j + e * elem_size] ? -1.0 : 1.0);
248         }
249       }
250     }
251   }
252   return CEED_ERROR_SUCCESS;
253 }
254 
255 static inline int CeedElemRestrictionApplyCurlOrientedTranspose_Ref_Core(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size,
256                                                                          const CeedInt comp_stride, CeedInt start, CeedInt stop, CeedInt num_elem,
257                                                                          CeedInt elem_size, CeedInt v_offset, const CeedScalar *uu, CeedScalar *vv) {
258   // Restriction with tridiagonal transformation
259   CeedElemRestriction_Ref *impl;
260 
261   CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl));
262   for (CeedInt e = start * block_size; e < stop * block_size; e += block_size) {
263     for (CeedInt k = 0; k < num_comp; k++) {
264       // Iteration bound set to discard padding elements
265       CeedInt block_end = CeedIntMin(block_size, num_elem - e), n = 0;
266       for (CeedInt j = 0; j < block_end; j++) {
267         vv[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] +=
268             uu[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] *
269                 impl->curl_orients[j + (3 * n + 1) * block_size + e * 3 * elem_size] +
270             uu[e * elem_size * num_comp + (k * elem_size + n + 1) * block_size + j - v_offset] *
271                 impl->curl_orients[j + (3 * n + 3) * block_size + e * 3 * elem_size];
272       }
273       for (n = 1; n < elem_size - 1; n++) {
274         for (CeedInt j = 0; j < block_end; j++) {
275           vv[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] +=
276               uu[e * elem_size * num_comp + (k * elem_size + n - 1) * block_size + j - v_offset] *
277                   impl->curl_orients[j + (3 * n - 1) * block_size + e * 3 * elem_size] +
278               uu[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] *
279                   impl->curl_orients[j + (3 * n + 1) * block_size + e * 3 * elem_size] +
280               uu[e * elem_size * num_comp + (k * elem_size + n + 1) * block_size + j - v_offset] *
281                   impl->curl_orients[j + (3 * n + 3) * block_size + e * 3 * elem_size];
282         }
283       }
284       for (CeedInt j = 0; j < block_end; j++) {
285         vv[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] +=
286             uu[e * elem_size * num_comp + (k * elem_size + n - 1) * block_size + j - v_offset] *
287                 impl->curl_orients[j + (3 * n - 1) * block_size + e * 3 * elem_size] +
288             uu[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] *
289                 impl->curl_orients[j + (3 * n + 1) * block_size + e * 3 * elem_size];
290       }
291     }
292   }
293   return CEED_ERROR_SUCCESS;
294 }
295 
296 static inline int CeedElemRestrictionApplyCurlOrientedUnsignedTranspose_Ref_Core(CeedElemRestriction rstr, const CeedInt num_comp,
297                                                                                  const CeedInt block_size, const CeedInt comp_stride, CeedInt start,
298                                                                                  CeedInt stop, CeedInt num_elem, CeedInt elem_size, CeedInt v_offset,
299                                                                                  const CeedScalar *uu, CeedScalar *vv) {
300   // Restriction with (unsigned) tridiagonal transformation
301   CeedElemRestriction_Ref *impl;
302 
303   CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl));
304   for (CeedInt e = start * block_size; e < stop * block_size; e += block_size) {
305     for (CeedInt k = 0; k < num_comp; k++) {
306       // Iteration bound set to discard padding elements
307       CeedInt       n         = 0;
308       const CeedInt block_end = CeedIntMin(block_size, num_elem - e);
309 
310       for (CeedInt j = 0; j < block_end; j++) {
311         vv[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] +=
312             uu[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] *
313                 abs(impl->curl_orients[j + (3 * n + 1) * block_size + e * 3 * elem_size]) +
314             uu[e * elem_size * num_comp + (k * elem_size + n + 1) * block_size + j - v_offset] *
315                 abs(impl->curl_orients[j + (3 * n + 3) * block_size + e * 3 * elem_size]);
316       }
317       for (n = 1; n < elem_size - 1; n++) {
318         for (CeedInt j = 0; j < block_end; j++) {
319           vv[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] +=
320               uu[e * elem_size * num_comp + (k * elem_size + n - 1) * block_size + j - v_offset] *
321                   abs(impl->curl_orients[j + (3 * n - 1) * block_size + e * 3 * elem_size]) +
322               uu[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] *
323                   abs(impl->curl_orients[j + (3 * n + 1) * block_size + e * 3 * elem_size]) +
324               uu[e * elem_size * num_comp + (k * elem_size + n + 1) * block_size + j - v_offset] *
325                   abs(impl->curl_orients[j + (3 * n + 3) * block_size + e * 3 * elem_size]);
326         }
327       }
328       for (CeedInt j = 0; j < block_end; j++) {
329         vv[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] +=
330             uu[e * elem_size * num_comp + (k * elem_size + n - 1) * block_size + j - v_offset] *
331                 abs(impl->curl_orients[j + (3 * n - 1) * block_size + e * 3 * elem_size]) +
332             uu[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] *
333                 abs(impl->curl_orients[j + (3 * n + 1) * block_size + e * 3 * elem_size]);
334       }
335     }
336   }
337   return CEED_ERROR_SUCCESS;
338 }
339 
340 static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size,
341                                                     const CeedInt comp_stride, CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs,
342                                                     bool use_orients, CeedVector u, CeedVector v, CeedRequest *request) {
343   CeedInt             num_elem, elem_size, v_offset;
344   CeedRestrictionType rstr_type;
345   const CeedScalar   *uu;
346   CeedScalar         *vv;
347 
348   CeedCallBackend(CeedElemRestrictionGetNumElements(rstr, &num_elem));
349   CeedCallBackend(CeedElemRestrictionGetElementSize(rstr, &elem_size));
350   v_offset = start * block_size * elem_size * num_comp;
351   CeedCallBackend(CeedElemRestrictionGetType(rstr, &rstr_type));
352   CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_HOST, &uu));
353 
354   if (t_mode == CEED_TRANSPOSE) {
355     // Sum into for transpose mode, E-vector to L-vector
356     CeedCallBackend(CeedVectorGetArray(v, CEED_MEM_HOST, &vv));
357   } else {
358     // Overwrite for notranspose mode, L-vector to E-vector
359     CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_HOST, &vv));
360   }
361   if (t_mode == CEED_TRANSPOSE) {
362     // Restriction from E-vector to L-vector
363     // Performing v += r^T * u
364     // uu has shape [elem_size, num_comp, num_elem], row-major
365     // vv has shape [nnodes, num_comp]
366     // Sum into for transpose mode
367     switch (rstr_type) {
368       case CEED_RESTRICTION_STRIDED:
369         CeedCallBackend(
370             CeedElemRestrictionApplyStridedTranspose_Ref_Core(rstr, num_comp, block_size, start, stop, num_elem, elem_size, v_offset, uu, vv));
371         break;
372       case CEED_RESTRICTION_STANDARD:
373         CeedCallBackend(CeedElemRestrictionApplyStandardTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem, elem_size,
374                                                                            v_offset, uu, vv));
375         break;
376       case CEED_RESTRICTION_ORIENTED:
377         if (use_signs) {
378           CeedCallBackend(CeedElemRestrictionApplyOrientedTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem,
379                                                                              elem_size, v_offset, uu, vv));
380         } else {
381           CeedCallBackend(CeedElemRestrictionApplyStandardTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem,
382                                                                              elem_size, v_offset, uu, vv));
383         }
384         break;
385       case CEED_RESTRICTION_CURL_ORIENTED:
386         if (use_signs && use_orients) {
387           CeedCallBackend(CeedElemRestrictionApplyCurlOrientedTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem,
388                                                                                  elem_size, v_offset, uu, vv));
389         } else if (use_orients) {
390           CeedCallBackend(CeedElemRestrictionApplyCurlOrientedUnsignedTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop,
391                                                                                          num_elem, elem_size, v_offset, uu, vv));
392         } else {
393           CeedCallBackend(CeedElemRestrictionApplyStandardTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem,
394                                                                              elem_size, v_offset, uu, vv));
395         }
396         break;
397       case CEED_RESTRICTION_POINTS:
398         // Empty case - won't occur
399         break;
400     }
401   } else {
402     // Restriction from L-vector to E-vector
403     // Perform: v = r * u
404     // vv has shape [elem_size, num_comp, num_elem], row-major
405     // uu has shape [nnodes, num_comp]
406     // Overwrite for notranspose mode
407     switch (rstr_type) {
408       case CEED_RESTRICTION_STRIDED:
409         CeedCallBackend(
410             CeedElemRestrictionApplyStridedNoTranspose_Ref_Core(rstr, num_comp, block_size, start, stop, num_elem, elem_size, v_offset, uu, vv));
411         break;
412       case CEED_RESTRICTION_STANDARD:
413         CeedCallBackend(CeedElemRestrictionApplyStandardNoTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem,
414                                                                              elem_size, v_offset, uu, vv));
415         break;
416       case CEED_RESTRICTION_ORIENTED:
417         if (use_signs) {
418           CeedCallBackend(CeedElemRestrictionApplyOrientedNoTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem,
419                                                                                elem_size, v_offset, uu, vv));
420         } else {
421           CeedCallBackend(CeedElemRestrictionApplyStandardNoTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem,
422                                                                                elem_size, v_offset, uu, vv));
423         }
424         break;
425       case CEED_RESTRICTION_CURL_ORIENTED:
426         if (use_signs && use_orients) {
427           CeedCallBackend(CeedElemRestrictionApplyCurlOrientedNoTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem,
428                                                                                    elem_size, v_offset, uu, vv));
429         } else if (use_orients) {
430           CeedCallBackend(CeedElemRestrictionApplyCurlOrientedUnsignedNoTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop,
431                                                                                            num_elem, elem_size, v_offset, uu, vv));
432         } else {
433           CeedCallBackend(CeedElemRestrictionApplyStandardNoTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem,
434                                                                                elem_size, v_offset, uu, vv));
435         }
436         break;
437       case CEED_RESTRICTION_POINTS:
438         // Empty case - won't occur
439         break;
440     }
441   }
442   CeedCallBackend(CeedVectorRestoreArrayRead(u, &uu));
443   CeedCallBackend(CeedVectorRestoreArray(v, &vv));
444   if (request != CEED_REQUEST_IMMEDIATE && request != CEED_REQUEST_ORDERED) *request = NULL;
445   return CEED_ERROR_SUCCESS;
446 }
447 
448 //------------------------------------------------------------------------------
449 // ElemRestriction Apply - Common Sizes
450 //------------------------------------------------------------------------------
451 static int CeedElemRestrictionApply_Ref_110(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride,
452                                             CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u,
453                                             CeedVector v, CeedRequest *request) {
454   return CeedElemRestrictionApply_Ref_Core(rstr, 1, 1, comp_stride, start, stop, t_mode, use_signs, use_orients, u, v, request);
455 }
456 
457 static int CeedElemRestrictionApply_Ref_111(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride,
458                                             CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u,
459                                             CeedVector v, CeedRequest *request) {
460   return CeedElemRestrictionApply_Ref_Core(rstr, 1, 1, 1, start, stop, t_mode, use_signs, use_orients, u, v, request);
461 }
462 
463 static int CeedElemRestrictionApply_Ref_180(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride,
464                                             CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u,
465                                             CeedVector v, CeedRequest *request) {
466   return CeedElemRestrictionApply_Ref_Core(rstr, 1, 8, comp_stride, start, stop, t_mode, use_signs, use_orients, u, v, request);
467 }
468 
469 static int CeedElemRestrictionApply_Ref_181(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride,
470                                             CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u,
471                                             CeedVector v, CeedRequest *request) {
472   return CeedElemRestrictionApply_Ref_Core(rstr, 1, 8, 1, start, stop, t_mode, use_signs, use_orients, u, v, request);
473 }
474 
475 static int CeedElemRestrictionApply_Ref_310(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride,
476                                             CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u,
477                                             CeedVector v, CeedRequest *request) {
478   return CeedElemRestrictionApply_Ref_Core(rstr, 3, 1, comp_stride, start, stop, t_mode, use_signs, use_orients, u, v, request);
479 }
480 
481 static int CeedElemRestrictionApply_Ref_311(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride,
482                                             CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u,
483                                             CeedVector v, CeedRequest *request) {
484   return CeedElemRestrictionApply_Ref_Core(rstr, 3, 1, 1, start, stop, t_mode, use_signs, use_orients, u, v, request);
485 }
486 
487 static int CeedElemRestrictionApply_Ref_380(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride,
488                                             CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u,
489                                             CeedVector v, CeedRequest *request) {
490   return CeedElemRestrictionApply_Ref_Core(rstr, 3, 8, comp_stride, start, stop, t_mode, use_signs, use_orients, u, v, request);
491 }
492 
493 static int CeedElemRestrictionApply_Ref_381(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride,
494                                             CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u,
495                                             CeedVector v, CeedRequest *request) {
496   return CeedElemRestrictionApply_Ref_Core(rstr, 3, 8, 1, start, stop, t_mode, use_signs, use_orients, u, v, request);
497 }
498 
499 // LCOV_EXCL_START
500 static int CeedElemRestrictionApply_Ref_510(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride,
501                                             CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u,
502                                             CeedVector v, CeedRequest *request) {
503   return CeedElemRestrictionApply_Ref_Core(rstr, 5, 1, comp_stride, start, stop, t_mode, use_signs, use_orients, u, v, request);
504 }
505 // LCOV_EXCL_STOP
506 
507 static int CeedElemRestrictionApply_Ref_511(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride,
508                                             CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u,
509                                             CeedVector v, CeedRequest *request) {
510   return CeedElemRestrictionApply_Ref_Core(rstr, 5, 1, 1, start, stop, t_mode, use_signs, use_orients, u, v, request);
511 }
512 
513 // LCOV_EXCL_START
514 static int CeedElemRestrictionApply_Ref_580(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride,
515                                             CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u,
516                                             CeedVector v, CeedRequest *request) {
517   return CeedElemRestrictionApply_Ref_Core(rstr, 5, 8, comp_stride, start, stop, t_mode, use_signs, use_orients, u, v, request);
518 }
519 // LCOV_EXCL_STOP
520 
521 static int CeedElemRestrictionApply_Ref_581(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride,
522                                             CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u,
523                                             CeedVector v, CeedRequest *request) {
524   return CeedElemRestrictionApply_Ref_Core(rstr, 5, 8, 1, start, stop, t_mode, use_signs, use_orients, u, v, request);
525 }
526 
527 //------------------------------------------------------------------------------
528 // ElemRestriction Apply
529 //------------------------------------------------------------------------------
530 static int CeedElemRestrictionApply_Ref(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) {
531   CeedInt                  num_block, block_size, num_comp, comp_stride;
532   CeedElemRestriction_Ref *impl;
533 
534   CeedCallBackend(CeedElemRestrictionGetNumBlocks(rstr, &num_block));
535   CeedCallBackend(CeedElemRestrictionGetBlockSize(rstr, &block_size));
536   CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr, &num_comp));
537   CeedCallBackend(CeedElemRestrictionGetCompStride(rstr, &comp_stride));
538   CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl));
539   CeedCallBackend(impl->Apply(rstr, num_comp, block_size, comp_stride, 0, num_block, t_mode, true, true, u, v, request));
540   return CEED_ERROR_SUCCESS;
541 }
542 
543 //------------------------------------------------------------------------------
544 // ElemRestriction Apply Unsigned
545 //------------------------------------------------------------------------------
546 static int CeedElemRestrictionApplyUnsigned_Ref(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector v,
547                                                 CeedRequest *request) {
548   CeedInt                  num_block, block_size, num_comp, comp_stride;
549   CeedElemRestriction_Ref *impl;
550 
551   CeedCallBackend(CeedElemRestrictionGetNumBlocks(rstr, &num_block));
552   CeedCallBackend(CeedElemRestrictionGetBlockSize(rstr, &block_size));
553   CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr, &num_comp));
554   CeedCallBackend(CeedElemRestrictionGetCompStride(rstr, &comp_stride));
555   CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl));
556   CeedCallBackend(impl->Apply(rstr, num_comp, block_size, comp_stride, 0, num_block, t_mode, false, true, u, v, request));
557   return CEED_ERROR_SUCCESS;
558 }
559 
560 //------------------------------------------------------------------------------
561 // ElemRestriction Apply Unoriented
562 //------------------------------------------------------------------------------
563 static int CeedElemRestrictionApplyUnoriented_Ref(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector v,
564                                                   CeedRequest *request) {
565   CeedInt                  num_block, block_size, num_comp, comp_stride;
566   CeedElemRestriction_Ref *impl;
567 
568   CeedCallBackend(CeedElemRestrictionGetNumBlocks(rstr, &num_block));
569   CeedCallBackend(CeedElemRestrictionGetBlockSize(rstr, &block_size));
570   CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr, &num_comp));
571   CeedCallBackend(CeedElemRestrictionGetCompStride(rstr, &comp_stride));
572   CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl));
573   CeedCallBackend(impl->Apply(rstr, num_comp, block_size, comp_stride, 0, num_block, t_mode, false, false, u, v, request));
574   return CEED_ERROR_SUCCESS;
575 }
576 
577 //------------------------------------------------------------------------------
578 // ElemRestriction Apply Points
579 //------------------------------------------------------------------------------
580 static int CeedElemRestrictionApplyAtPoints_Ref(CeedElemRestriction r, CeedInt elem, CeedTransposeMode t_mode, CeedVector u, CeedVector v,
581                                                 CeedRequest *request) {
582   CeedInt                  num_points, num_comp, l_vec_offset;
583   CeedScalar              *vv;
584   const CeedScalar        *uu;
585   CeedElemRestriction_Ref *impl;
586 
587   CeedCallBackend(CeedElemRestrictionGetData(r, &impl));
588   CeedCallBackend(CeedElemRestrictionGetNumComponents(r, &num_comp));
589   CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(r, elem, &num_points));
590 
591   CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_HOST, &uu));
592   if (t_mode == CEED_TRANSPOSE) {
593     // Sum into for transpose mode, E-vector to L-vector
594     CeedCallBackend(CeedVectorGetArray(v, CEED_MEM_HOST, &vv));
595   } else {
596     // Overwrite for notranspose mode, L-vector to E-vector
597     CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_HOST, &vv));
598   }
599 
600   l_vec_offset = impl->offsets[elem + 1] - impl->offsets[elem];
601   if (t_mode == CEED_NOTRANSPOSE) {
602     for (CeedInt i = 0; i < num_points; i++) {
603       for (CeedInt j = 0; j < num_comp; j++) vv[i * num_comp + j] = uu[impl->offsets[i + l_vec_offset] * num_comp + j];
604     }
605   } else {
606     for (CeedInt i = 0; i < num_points; i++) {
607       for (CeedInt j = 0; j < num_comp; j++) vv[impl->offsets[i + l_vec_offset] * num_comp + j] = uu[i * num_comp + j];
608     }
609   }
610 
611   return CEED_ERROR_SUCCESS;
612 }
613 
614 //------------------------------------------------------------------------------
615 // ElemRestriction Apply Block
616 //------------------------------------------------------------------------------
617 static int CeedElemRestrictionApplyBlock_Ref(CeedElemRestriction rstr, CeedInt block, CeedTransposeMode t_mode, CeedVector u, CeedVector v,
618                                              CeedRequest *request) {
619   CeedInt                  block_size, num_comp, comp_stride;
620   CeedElemRestriction_Ref *impl;
621 
622   CeedCallBackend(CeedElemRestrictionGetBlockSize(rstr, &block_size));
623   CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr, &num_comp));
624   CeedCallBackend(CeedElemRestrictionGetCompStride(rstr, &comp_stride));
625   CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl));
626   CeedCallBackend(impl->Apply(rstr, num_comp, block_size, comp_stride, block, block + 1, t_mode, true, true, u, v, request));
627   return CEED_ERROR_SUCCESS;
628 }
629 
630 //------------------------------------------------------------------------------
631 // ElemRestriction Get Offsets
632 //------------------------------------------------------------------------------
633 static int CeedElemRestrictionGetOffsets_Ref(CeedElemRestriction rstr, CeedMemType mem_type, const CeedInt **offsets) {
634   Ceed                     ceed;
635   CeedElemRestriction_Ref *impl;
636 
637   CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl));
638   CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed));
639 
640   CeedCheck(mem_type == CEED_MEM_HOST, ceed, CEED_ERROR_BACKEND, "Can only provide to HOST memory");
641 
642   *offsets = impl->offsets;
643   return CEED_ERROR_SUCCESS;
644 }
645 
646 //------------------------------------------------------------------------------
647 // ElemRestriction Get Orientations
648 //------------------------------------------------------------------------------
649 static int CeedElemRestrictionGetOrientations_Ref(CeedElemRestriction rstr, CeedMemType mem_type, const bool **orients) {
650   Ceed                     ceed;
651   CeedElemRestriction_Ref *impl;
652 
653   CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl));
654   CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed));
655 
656   CeedCheck(mem_type == CEED_MEM_HOST, ceed, CEED_ERROR_BACKEND, "Can only provide to HOST memory");
657 
658   *orients = impl->orients;
659   return CEED_ERROR_SUCCESS;
660 }
661 
662 //------------------------------------------------------------------------------
663 // ElemRestriction Get Curl-Conforming Orientations
664 //------------------------------------------------------------------------------
665 static int CeedElemRestrictionGetCurlOrientations_Ref(CeedElemRestriction rstr, CeedMemType mem_type, const CeedInt8 **curl_orients) {
666   Ceed                     ceed;
667   CeedElemRestriction_Ref *impl;
668 
669   CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl));
670   CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed));
671 
672   CeedCheck(mem_type == CEED_MEM_HOST, ceed, CEED_ERROR_BACKEND, "Can only provide to HOST memory");
673 
674   *curl_orients = impl->curl_orients;
675   return CEED_ERROR_SUCCESS;
676 }
677 
678 //------------------------------------------------------------------------------
679 // ElemRestriction Destroy
680 //------------------------------------------------------------------------------
681 static int CeedElemRestrictionDestroy_Ref(CeedElemRestriction rstr) {
682   CeedElemRestriction_Ref *impl;
683 
684   CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl));
685   CeedCallBackend(CeedFree(&impl->offsets_allocated));
686   CeedCallBackend(CeedFree(&impl->orients_allocated));
687   CeedCallBackend(CeedFree(&impl->curl_orients_allocated));
688   CeedCallBackend(CeedFree(&impl));
689   return CEED_ERROR_SUCCESS;
690 }
691 
692 //------------------------------------------------------------------------------
693 // ElemRestriction Create
694 //------------------------------------------------------------------------------
695 int CeedElemRestrictionCreate_Ref(CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, const bool *orients,
696                                   const CeedInt8 *curl_orients, CeedElemRestriction rstr) {
697   Ceed                     ceed;
698   CeedInt                  num_elem, elem_size, num_block, block_size, num_comp, comp_stride;
699   CeedRestrictionType      rstr_type;
700   CeedElemRestriction_Ref *impl;
701 
702   CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed));
703   CeedCallBackend(CeedElemRestrictionGetNumElements(rstr, &num_elem));
704   CeedCallBackend(CeedElemRestrictionGetElementSize(rstr, &elem_size));
705   CeedCallBackend(CeedElemRestrictionGetNumBlocks(rstr, &num_block));
706   CeedCallBackend(CeedElemRestrictionGetBlockSize(rstr, &block_size));
707   CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr, &num_comp));
708   CeedCallBackend(CeedElemRestrictionGetCompStride(rstr, &comp_stride));
709   CeedInt layout[3] = {1, elem_size, elem_size * num_comp};
710 
711   CeedCheck(mem_type == CEED_MEM_HOST, ceed, CEED_ERROR_BACKEND, "Only MemType = HOST supported");
712   CeedCallBackend(CeedCalloc(1, &impl));
713 
714   // Offsets data
715   CeedCallBackend(CeedElemRestrictionGetType(rstr, &rstr_type));
716   if (rstr_type != CEED_RESTRICTION_STRIDED) {
717     const char *resource;
718 
719     // Check indices for ref or memcheck backends
720     CeedCallBackend(CeedGetResource(ceed, &resource));
721     if (!strcmp(resource, "/cpu/self/ref/serial") || !strcmp(resource, "/cpu/self/ref/blocked") || !strcmp(resource, "/cpu/self/memcheck/serial") ||
722         !strcmp(resource, "/cpu/self/memcheck/blocked")) {
723       CeedSize l_size;
724 
725       CeedCallBackend(CeedElemRestrictionGetLVectorSize(rstr, &l_size));
726       for (CeedInt i = 0; i < num_elem * elem_size; i++) {
727         CeedCheck(offsets[i] >= 0 && offsets[i] + (num_comp - 1) * comp_stride < l_size, ceed, CEED_ERROR_BACKEND,
728                   "Restriction offset %" CeedInt_FMT " (%" CeedInt_FMT ") out of range [0, %" CeedInt_FMT "]", i, offsets[i], l_size);
729       }
730     }
731 
732     // Copy data
733     switch (copy_mode) {
734       case CEED_COPY_VALUES:
735         CeedCallBackend(CeedMalloc(num_elem * elem_size, &impl->offsets_allocated));
736         memcpy(impl->offsets_allocated, offsets, num_elem * elem_size * sizeof(offsets[0]));
737         impl->offsets = impl->offsets_allocated;
738         break;
739       case CEED_OWN_POINTER:
740         impl->offsets_allocated = (CeedInt *)offsets;
741         impl->offsets           = impl->offsets_allocated;
742         break;
743       case CEED_USE_POINTER:
744         impl->offsets = offsets;
745     }
746 
747     // Orientation data
748     if (rstr_type == CEED_RESTRICTION_ORIENTED) {
749       CeedCheck(orients != NULL, ceed, CEED_ERROR_BACKEND, "No orients array provided for oriented restriction");
750       switch (copy_mode) {
751         case CEED_COPY_VALUES:
752           CeedCallBackend(CeedMalloc(num_elem * elem_size, &impl->orients_allocated));
753           memcpy(impl->orients_allocated, orients, num_elem * elem_size * sizeof(orients[0]));
754           impl->orients = impl->orients_allocated;
755           break;
756         case CEED_OWN_POINTER:
757           impl->orients_allocated = (bool *)orients;
758           impl->orients           = impl->orients_allocated;
759           break;
760         case CEED_USE_POINTER:
761           impl->orients = orients;
762       }
763     } else if (rstr_type == CEED_RESTRICTION_CURL_ORIENTED) {
764       CeedCheck(curl_orients != NULL, ceed, CEED_ERROR_BACKEND, "No curl_orients array provided for oriented restriction");
765       switch (copy_mode) {
766         case CEED_COPY_VALUES:
767           CeedCallBackend(CeedMalloc(num_elem * 3 * elem_size, &impl->curl_orients_allocated));
768           memcpy(impl->curl_orients_allocated, curl_orients, num_elem * 3 * elem_size * sizeof(curl_orients[0]));
769           impl->curl_orients = impl->curl_orients_allocated;
770           break;
771         case CEED_OWN_POINTER:
772           impl->curl_orients_allocated = (CeedInt8 *)curl_orients;
773           impl->curl_orients           = impl->curl_orients_allocated;
774           break;
775         case CEED_USE_POINTER:
776           impl->curl_orients = curl_orients;
777       }
778     }
779   }
780 
781   CeedCallBackend(CeedElemRestrictionSetData(rstr, impl));
782   CeedCallBackend(CeedElemRestrictionSetELayout(rstr, layout));
783   CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "Apply", CeedElemRestrictionApply_Ref));
784   if (rstr_type == CEED_RESTRICTION_POINTS) {
785     CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "ApplyAtPoints", CeedElemRestrictionApplyAtPoints_Ref));
786   } else {
787     CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "ApplyUnsigned", CeedElemRestrictionApplyUnsigned_Ref));
788     CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "ApplyUnoriented", CeedElemRestrictionApplyUnoriented_Ref));
789     CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "ApplyBlock", CeedElemRestrictionApplyBlock_Ref));
790     CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "GetOrientations", CeedElemRestrictionGetOrientations_Ref));
791     CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "GetCurlOrientations", CeedElemRestrictionGetCurlOrientations_Ref));
792   }
793   CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "GetOffsets", CeedElemRestrictionGetOffsets_Ref));
794   CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "Destroy", CeedElemRestrictionDestroy_Ref));
795 
796   // Set apply function based upon num_comp, block_size, and comp_stride
797   CeedInt index = -1;
798 
799   if (block_size < 10) index = 100 * num_comp + 10 * block_size + (comp_stride == 1);
800   switch (index) {
801     case 110:
802       impl->Apply = CeedElemRestrictionApply_Ref_110;
803       break;
804     case 111:
805       impl->Apply = CeedElemRestrictionApply_Ref_111;
806       break;
807     case 180:
808       impl->Apply = CeedElemRestrictionApply_Ref_180;
809       break;
810     case 181:
811       impl->Apply = CeedElemRestrictionApply_Ref_181;
812       break;
813     case 310:
814       impl->Apply = CeedElemRestrictionApply_Ref_310;
815       break;
816     case 311:
817       impl->Apply = CeedElemRestrictionApply_Ref_311;
818       break;
819     case 380:
820       impl->Apply = CeedElemRestrictionApply_Ref_380;
821       break;
822     case 381:
823       impl->Apply = CeedElemRestrictionApply_Ref_381;
824       break;
825     // LCOV_EXCL_START
826     case 510:
827       impl->Apply = CeedElemRestrictionApply_Ref_510;
828       break;
829     // LCOV_EXCL_STOP
830     case 511:
831       impl->Apply = CeedElemRestrictionApply_Ref_511;
832       break;
833     // LCOV_EXCL_START
834     case 580:
835       impl->Apply = CeedElemRestrictionApply_Ref_580;
836       break;
837     // LCOV_EXCL_STOP
838     case 581:
839       impl->Apply = CeedElemRestrictionApply_Ref_581;
840       break;
841     default:
842       impl->Apply = CeedElemRestrictionApply_Ref_Core;
843       break;
844   }
845   return CEED_ERROR_SUCCESS;
846 }
847 
848 //------------------------------------------------------------------------------
849