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 CeedScalar uu_val; 227 228 uu_val = uu[elem_size * (k * block_size + e * num_comp) + j - v_offset]; 229 CeedPragmaAtomic vv[impl->offsets[j + e * elem_size] + k * comp_stride] += uu_val; 230 } 231 } 232 } 233 } 234 return CEED_ERROR_SUCCESS; 235 } 236 237 static inline int CeedElemRestrictionApplyOrientedTranspose_Ref_Core(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, 238 const CeedInt comp_stride, CeedInt start, CeedInt stop, CeedInt num_elem, 239 CeedInt elem_size, CeedInt v_offset, const CeedScalar *uu, CeedScalar *vv) { 240 // Restriction with orientations 241 CeedElemRestriction_Ref *impl; 242 243 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 244 for (CeedInt e = start * block_size; e < stop * block_size; e += block_size) { 245 for (CeedInt k = 0; k < num_comp; k++) { 246 for (CeedInt i = 0; i < elem_size * block_size; i += block_size) { 247 // Iteration bound set to discard padding elements 248 for (CeedInt j = i; j < i + CeedIntMin(block_size, num_elem - e); j++) { 249 CeedScalar uu_val; 250 251 uu_val = uu[elem_size * (k * block_size + e * num_comp) + j - v_offset] * (impl->orients[j + e * elem_size] ? -1.0 : 1.0); 252 CeedPragmaAtomic vv[impl->offsets[j + e * elem_size] + k * comp_stride] += uu_val; 253 } 254 } 255 } 256 } 257 return CEED_ERROR_SUCCESS; 258 } 259 260 static inline int CeedElemRestrictionApplyCurlOrientedTranspose_Ref_Core(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, 261 const CeedInt comp_stride, CeedInt start, CeedInt stop, CeedInt num_elem, 262 CeedInt elem_size, CeedInt v_offset, const CeedScalar *uu, CeedScalar *vv) { 263 // Restriction with tridiagonal transformation 264 CeedElemRestriction_Ref *impl; 265 266 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 267 for (CeedInt e = start * block_size; e < stop * block_size; e += block_size) { 268 for (CeedInt k = 0; k < num_comp; k++) { 269 // Iteration bound set to discard padding elements 270 const CeedInt block_end = CeedIntMin(block_size, num_elem - e); 271 CeedInt n = 0; 272 273 for (CeedInt j = 0; j < block_end; j++) { 274 CeedScalar uu_val; 275 276 uu_val = uu[e * elem_size * num_comp + (k * elem_size + n) * 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 + 1) * block_size + j - v_offset] * 279 impl->curl_orients[j + (3 * n + 3) * block_size + e * 3 * elem_size]; 280 CeedPragmaAtomic vv[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] += uu_val; 281 } 282 for (n = 1; n < elem_size - 1; n++) { 283 for (CeedInt j = 0; j < block_end; j++) { 284 CeedScalar uu_val; 285 286 uu_val = 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 uu[e * elem_size * num_comp + (k * elem_size + n + 1) * block_size + j - v_offset] * 291 impl->curl_orients[j + (3 * n + 3) * block_size + e * 3 * elem_size]; 292 CeedPragmaAtomic vv[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] += uu_val; 293 } 294 } 295 for (CeedInt j = 0; j < block_end; j++) { 296 CeedScalar uu_val; 297 298 uu_val = uu[e * elem_size * num_comp + (k * elem_size + n - 1) * block_size + j - v_offset] * 299 impl->curl_orients[j + (3 * n - 1) * block_size + e * 3 * elem_size] + 300 uu[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] * 301 impl->curl_orients[j + (3 * n + 1) * block_size + e * 3 * elem_size]; 302 CeedPragmaAtomic vv[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] += uu_val; 303 } 304 } 305 } 306 return CEED_ERROR_SUCCESS; 307 } 308 309 static inline int CeedElemRestrictionApplyCurlOrientedUnsignedTranspose_Ref_Core(CeedElemRestriction rstr, const CeedInt num_comp, 310 const CeedInt block_size, const CeedInt comp_stride, CeedInt start, 311 CeedInt stop, CeedInt num_elem, CeedInt elem_size, CeedInt v_offset, 312 const CeedScalar *uu, CeedScalar *vv) { 313 // Restriction with (unsigned) tridiagonal transformation 314 CeedElemRestriction_Ref *impl; 315 316 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 317 for (CeedInt e = start * block_size; e < stop * block_size; e += block_size) { 318 for (CeedInt k = 0; k < num_comp; k++) { 319 // Iteration bound set to discard padding elements 320 const CeedInt block_end = CeedIntMin(block_size, num_elem - e); 321 CeedInt n = 0; 322 323 for (CeedInt j = 0; j < block_end; j++) { 324 CeedScalar uu_val; 325 326 uu_val = uu[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] * 327 abs(impl->curl_orients[j + (3 * n + 1) * block_size + e * 3 * elem_size]) + 328 uu[e * elem_size * num_comp + (k * elem_size + n + 1) * block_size + j - v_offset] * 329 abs(impl->curl_orients[j + (3 * n + 3) * block_size + e * 3 * elem_size]); 330 CeedPragmaAtomic vv[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] += uu_val; 331 } 332 for (n = 1; n < elem_size - 1; n++) { 333 for (CeedInt j = 0; j < block_end; j++) { 334 CeedScalar uu_val; 335 336 uu_val = uu[e * elem_size * num_comp + (k * elem_size + n - 1) * block_size + j - v_offset] * 337 abs(impl->curl_orients[j + (3 * n - 1) * block_size + e * 3 * elem_size]) + 338 uu[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] * 339 abs(impl->curl_orients[j + (3 * n + 1) * block_size + e * 3 * elem_size]) + 340 uu[e * elem_size * num_comp + (k * elem_size + n + 1) * block_size + j - v_offset] * 341 abs(impl->curl_orients[j + (3 * n + 3) * block_size + e * 3 * elem_size]); 342 CeedPragmaAtomic vv[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] += uu_val; 343 } 344 } 345 for (CeedInt j = 0; j < block_end; j++) { 346 CeedScalar uu_val; 347 348 uu_val = uu[e * elem_size * num_comp + (k * elem_size + n - 1) * block_size + j - v_offset] * 349 abs(impl->curl_orients[j + (3 * n - 1) * block_size + e * 3 * elem_size]) + 350 uu[e * elem_size * num_comp + (k * elem_size + n) * block_size + j - v_offset] * 351 abs(impl->curl_orients[j + (3 * n + 1) * block_size + e * 3 * elem_size]); 352 CeedPragmaAtomic vv[impl->offsets[j + n * block_size + e * elem_size] + k * comp_stride] += uu_val; 353 } 354 } 355 } 356 return CEED_ERROR_SUCCESS; 357 } 358 359 static inline int CeedElemRestrictionApplyAtPointsInElement_Ref_Core(CeedElemRestriction rstr, const CeedInt num_comp, CeedInt start, CeedInt stop, 360 CeedTransposeMode t_mode, const CeedScalar *uu, CeedScalar *vv) { 361 CeedInt num_points, l_vec_offset, e_vec_offset = 0; 362 CeedElemRestriction_Ref *impl; 363 364 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 365 366 for (CeedInt e = start; e < stop; e++) { 367 l_vec_offset = impl->offsets[e]; 368 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr, e, &num_points)); 369 if (t_mode == CEED_NOTRANSPOSE) { 370 for (CeedInt i = 0; i < num_points; i++) { 371 for (CeedInt j = 0; j < num_comp; j++) vv[j * num_points + i + e_vec_offset] = uu[impl->offsets[i + l_vec_offset] * num_comp + j]; 372 } 373 } else { 374 for (CeedInt i = 0; i < num_points; i++) { 375 for (CeedInt j = 0; j < num_comp; j++) vv[impl->offsets[i + l_vec_offset] * num_comp + j] = uu[j * num_points + i + e_vec_offset]; 376 } 377 } 378 e_vec_offset += num_points * num_comp; 379 } 380 return CEED_ERROR_SUCCESS; 381 } 382 383 static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, 384 const CeedInt comp_stride, CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, 385 bool use_orients, CeedVector u, CeedVector v, CeedRequest *request) { 386 CeedInt num_elem, elem_size, v_offset; 387 CeedRestrictionType rstr_type; 388 const CeedScalar *uu; 389 CeedScalar *vv; 390 391 CeedCallBackend(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 392 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr, &elem_size)); 393 v_offset = start * block_size * elem_size * num_comp; 394 CeedCallBackend(CeedElemRestrictionGetType(rstr, &rstr_type)); 395 CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_HOST, &uu)); 396 397 if (t_mode == CEED_TRANSPOSE) { 398 // Sum into for transpose mode, E-vector to L-vector 399 CeedCallBackend(CeedVectorGetArray(v, CEED_MEM_HOST, &vv)); 400 } else { 401 // Overwrite for notranspose mode, L-vector to E-vector 402 CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_HOST, &vv)); 403 } 404 if (t_mode == CEED_TRANSPOSE) { 405 // Restriction from E-vector to L-vector 406 // Performing v += r^T * u 407 // uu has shape [elem_size, num_comp, num_elem], row-major 408 // vv has shape [nnodes, num_comp] 409 // Sum into for transpose mode 410 switch (rstr_type) { 411 case CEED_RESTRICTION_STRIDED: 412 CeedCallBackend( 413 CeedElemRestrictionApplyStridedTranspose_Ref_Core(rstr, num_comp, block_size, start, stop, num_elem, elem_size, v_offset, uu, vv)); 414 break; 415 case CEED_RESTRICTION_STANDARD: 416 CeedCallBackend(CeedElemRestrictionApplyStandardTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem, elem_size, 417 v_offset, uu, vv)); 418 break; 419 case CEED_RESTRICTION_ORIENTED: 420 if (use_signs) { 421 CeedCallBackend(CeedElemRestrictionApplyOrientedTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem, 422 elem_size, v_offset, uu, vv)); 423 } else { 424 CeedCallBackend(CeedElemRestrictionApplyStandardTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem, 425 elem_size, v_offset, uu, vv)); 426 } 427 break; 428 case CEED_RESTRICTION_CURL_ORIENTED: 429 if (use_signs && use_orients) { 430 CeedCallBackend(CeedElemRestrictionApplyCurlOrientedTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem, 431 elem_size, v_offset, uu, vv)); 432 } else if (use_orients) { 433 CeedCallBackend(CeedElemRestrictionApplyCurlOrientedUnsignedTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, 434 num_elem, elem_size, v_offset, uu, vv)); 435 } else { 436 CeedCallBackend(CeedElemRestrictionApplyStandardTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem, 437 elem_size, v_offset, uu, vv)); 438 } 439 break; 440 case CEED_RESTRICTION_POINTS: 441 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement_Ref_Core(rstr, num_comp, start, stop, t_mode, uu, vv)); 442 break; 443 } 444 } else { 445 // Restriction from L-vector to E-vector 446 // Perform: v = r * u 447 // vv has shape [elem_size, num_comp, num_elem], row-major 448 // uu has shape [nnodes, num_comp] 449 // Overwrite for notranspose mode 450 switch (rstr_type) { 451 case CEED_RESTRICTION_STRIDED: 452 CeedCallBackend( 453 CeedElemRestrictionApplyStridedNoTranspose_Ref_Core(rstr, num_comp, block_size, start, stop, num_elem, elem_size, v_offset, uu, vv)); 454 break; 455 case CEED_RESTRICTION_STANDARD: 456 CeedCallBackend(CeedElemRestrictionApplyStandardNoTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem, 457 elem_size, v_offset, uu, vv)); 458 break; 459 case CEED_RESTRICTION_ORIENTED: 460 if (use_signs) { 461 CeedCallBackend(CeedElemRestrictionApplyOrientedNoTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem, 462 elem_size, v_offset, uu, vv)); 463 } else { 464 CeedCallBackend(CeedElemRestrictionApplyStandardNoTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem, 465 elem_size, v_offset, uu, vv)); 466 } 467 break; 468 case CEED_RESTRICTION_CURL_ORIENTED: 469 if (use_signs && use_orients) { 470 CeedCallBackend(CeedElemRestrictionApplyCurlOrientedNoTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem, 471 elem_size, v_offset, uu, vv)); 472 } else if (use_orients) { 473 CeedCallBackend(CeedElemRestrictionApplyCurlOrientedUnsignedNoTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, 474 num_elem, elem_size, v_offset, uu, vv)); 475 } else { 476 CeedCallBackend(CeedElemRestrictionApplyStandardNoTranspose_Ref_Core(rstr, num_comp, block_size, comp_stride, start, stop, num_elem, 477 elem_size, v_offset, uu, vv)); 478 } 479 break; 480 case CEED_RESTRICTION_POINTS: 481 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement_Ref_Core(rstr, num_comp, start, stop, t_mode, uu, vv)); 482 break; 483 } 484 } 485 CeedCallBackend(CeedVectorRestoreArrayRead(u, &uu)); 486 CeedCallBackend(CeedVectorRestoreArray(v, &vv)); 487 if (request != CEED_REQUEST_IMMEDIATE && request != CEED_REQUEST_ORDERED) *request = NULL; 488 return CEED_ERROR_SUCCESS; 489 } 490 491 //------------------------------------------------------------------------------ 492 // ElemRestriction Apply - Common Sizes 493 //------------------------------------------------------------------------------ 494 static int CeedElemRestrictionApply_Ref_110(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride, 495 CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u, 496 CeedVector v, CeedRequest *request) { 497 return CeedElemRestrictionApply_Ref_Core(rstr, 1, 1, comp_stride, start, stop, t_mode, use_signs, use_orients, u, v, request); 498 } 499 500 static int CeedElemRestrictionApply_Ref_111(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, 1, 1, 1, start, stop, t_mode, use_signs, use_orients, u, v, request); 504 } 505 506 static int CeedElemRestrictionApply_Ref_180(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride, 507 CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u, 508 CeedVector v, CeedRequest *request) { 509 return CeedElemRestrictionApply_Ref_Core(rstr, 1, 8, comp_stride, start, stop, t_mode, use_signs, use_orients, u, v, request); 510 } 511 512 static int CeedElemRestrictionApply_Ref_181(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride, 513 CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u, 514 CeedVector v, CeedRequest *request) { 515 return CeedElemRestrictionApply_Ref_Core(rstr, 1, 8, 1, start, stop, t_mode, use_signs, use_orients, u, v, request); 516 } 517 518 static int CeedElemRestrictionApply_Ref_310(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride, 519 CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u, 520 CeedVector v, CeedRequest *request) { 521 return CeedElemRestrictionApply_Ref_Core(rstr, 3, 1, comp_stride, start, stop, t_mode, use_signs, use_orients, u, v, request); 522 } 523 524 static int CeedElemRestrictionApply_Ref_311(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride, 525 CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u, 526 CeedVector v, CeedRequest *request) { 527 return CeedElemRestrictionApply_Ref_Core(rstr, 3, 1, 1, start, stop, t_mode, use_signs, use_orients, u, v, request); 528 } 529 530 static int CeedElemRestrictionApply_Ref_380(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride, 531 CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u, 532 CeedVector v, CeedRequest *request) { 533 return CeedElemRestrictionApply_Ref_Core(rstr, 3, 8, comp_stride, start, stop, t_mode, use_signs, use_orients, u, v, request); 534 } 535 536 static int CeedElemRestrictionApply_Ref_381(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride, 537 CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u, 538 CeedVector v, CeedRequest *request) { 539 return CeedElemRestrictionApply_Ref_Core(rstr, 3, 8, 1, start, stop, t_mode, use_signs, use_orients, u, v, request); 540 } 541 542 // LCOV_EXCL_START 543 static int CeedElemRestrictionApply_Ref_510(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride, 544 CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u, 545 CeedVector v, CeedRequest *request) { 546 return CeedElemRestrictionApply_Ref_Core(rstr, 5, 1, comp_stride, start, stop, t_mode, use_signs, use_orients, u, v, request); 547 } 548 // LCOV_EXCL_STOP 549 550 static int CeedElemRestrictionApply_Ref_511(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride, 551 CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u, 552 CeedVector v, CeedRequest *request) { 553 return CeedElemRestrictionApply_Ref_Core(rstr, 5, 1, 1, start, stop, t_mode, use_signs, use_orients, u, v, request); 554 } 555 556 // LCOV_EXCL_START 557 static int CeedElemRestrictionApply_Ref_580(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride, 558 CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u, 559 CeedVector v, CeedRequest *request) { 560 return CeedElemRestrictionApply_Ref_Core(rstr, 5, 8, comp_stride, start, stop, t_mode, use_signs, use_orients, u, v, request); 561 } 562 // LCOV_EXCL_STOP 563 564 static int CeedElemRestrictionApply_Ref_581(CeedElemRestriction rstr, const CeedInt num_comp, const CeedInt block_size, const CeedInt comp_stride, 565 CeedInt start, CeedInt stop, CeedTransposeMode t_mode, bool use_signs, bool use_orients, CeedVector u, 566 CeedVector v, CeedRequest *request) { 567 return CeedElemRestrictionApply_Ref_Core(rstr, 5, 8, 1, start, stop, t_mode, use_signs, use_orients, u, v, request); 568 } 569 570 //------------------------------------------------------------------------------ 571 // ElemRestriction Apply 572 //------------------------------------------------------------------------------ 573 static int CeedElemRestrictionApply_Ref(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { 574 CeedInt num_block, block_size, num_comp, comp_stride; 575 CeedElemRestriction_Ref *impl; 576 577 CeedCallBackend(CeedElemRestrictionGetNumBlocks(rstr, &num_block)); 578 CeedCallBackend(CeedElemRestrictionGetBlockSize(rstr, &block_size)); 579 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 580 CeedCallBackend(CeedElemRestrictionGetCompStride(rstr, &comp_stride)); 581 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 582 CeedCallBackend(impl->Apply(rstr, num_comp, block_size, comp_stride, 0, num_block, t_mode, true, true, u, v, request)); 583 return CEED_ERROR_SUCCESS; 584 } 585 586 //------------------------------------------------------------------------------ 587 // ElemRestriction Apply Unsigned 588 //------------------------------------------------------------------------------ 589 static int CeedElemRestrictionApplyUnsigned_Ref(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 590 CeedRequest *request) { 591 CeedInt num_block, block_size, num_comp, comp_stride; 592 CeedElemRestriction_Ref *impl; 593 594 CeedCallBackend(CeedElemRestrictionGetNumBlocks(rstr, &num_block)); 595 CeedCallBackend(CeedElemRestrictionGetBlockSize(rstr, &block_size)); 596 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 597 CeedCallBackend(CeedElemRestrictionGetCompStride(rstr, &comp_stride)); 598 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 599 CeedCallBackend(impl->Apply(rstr, num_comp, block_size, comp_stride, 0, num_block, t_mode, false, true, u, v, request)); 600 return CEED_ERROR_SUCCESS; 601 } 602 603 //------------------------------------------------------------------------------ 604 // ElemRestriction Apply Unoriented 605 //------------------------------------------------------------------------------ 606 static int CeedElemRestrictionApplyUnoriented_Ref(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 607 CeedRequest *request) { 608 CeedInt num_block, block_size, num_comp, comp_stride; 609 CeedElemRestriction_Ref *impl; 610 611 CeedCallBackend(CeedElemRestrictionGetNumBlocks(rstr, &num_block)); 612 CeedCallBackend(CeedElemRestrictionGetBlockSize(rstr, &block_size)); 613 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 614 CeedCallBackend(CeedElemRestrictionGetCompStride(rstr, &comp_stride)); 615 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 616 CeedCallBackend(impl->Apply(rstr, num_comp, block_size, comp_stride, 0, num_block, t_mode, false, false, u, v, request)); 617 return CEED_ERROR_SUCCESS; 618 } 619 620 //------------------------------------------------------------------------------ 621 // ElemRestriction Apply Points 622 //------------------------------------------------------------------------------ 623 static int CeedElemRestrictionApplyAtPointsInElement_Ref(CeedElemRestriction r, CeedInt elem, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 624 CeedRequest *request) { 625 CeedInt num_comp; 626 CeedElemRestriction_Ref *impl; 627 628 CeedCallBackend(CeedElemRestrictionGetNumComponents(r, &num_comp)); 629 CeedCallBackend(CeedElemRestrictionGetData(r, &impl)); 630 return impl->Apply(r, num_comp, 0, 1, elem, elem + 1, t_mode, false, false, u, v, request); 631 } 632 633 //------------------------------------------------------------------------------ 634 // ElemRestriction Apply Block 635 //------------------------------------------------------------------------------ 636 static int CeedElemRestrictionApplyBlock_Ref(CeedElemRestriction rstr, CeedInt block, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 637 CeedRequest *request) { 638 CeedInt block_size, num_comp, comp_stride; 639 CeedElemRestriction_Ref *impl; 640 641 CeedCallBackend(CeedElemRestrictionGetBlockSize(rstr, &block_size)); 642 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 643 CeedCallBackend(CeedElemRestrictionGetCompStride(rstr, &comp_stride)); 644 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 645 CeedCallBackend(impl->Apply(rstr, num_comp, block_size, comp_stride, block, block + 1, t_mode, true, true, u, v, request)); 646 return CEED_ERROR_SUCCESS; 647 } 648 649 //------------------------------------------------------------------------------ 650 // ElemRestriction Get Offsets 651 //------------------------------------------------------------------------------ 652 static int CeedElemRestrictionGetOffsets_Ref(CeedElemRestriction rstr, CeedMemType mem_type, const CeedInt **offsets) { 653 Ceed ceed; 654 CeedElemRestriction_Ref *impl; 655 656 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 657 CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed)); 658 659 CeedCheck(mem_type == CEED_MEM_HOST, ceed, CEED_ERROR_BACKEND, "Can only provide to HOST memory"); 660 661 *offsets = impl->offsets; 662 return CEED_ERROR_SUCCESS; 663 } 664 665 //------------------------------------------------------------------------------ 666 // ElemRestriction Get Orientations 667 //------------------------------------------------------------------------------ 668 static int CeedElemRestrictionGetOrientations_Ref(CeedElemRestriction rstr, CeedMemType mem_type, const bool **orients) { 669 Ceed ceed; 670 CeedElemRestriction_Ref *impl; 671 672 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 673 CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed)); 674 675 CeedCheck(mem_type == CEED_MEM_HOST, ceed, CEED_ERROR_BACKEND, "Can only provide to HOST memory"); 676 677 *orients = impl->orients; 678 return CEED_ERROR_SUCCESS; 679 } 680 681 //------------------------------------------------------------------------------ 682 // ElemRestriction Get Curl-Conforming Orientations 683 //------------------------------------------------------------------------------ 684 static int CeedElemRestrictionGetCurlOrientations_Ref(CeedElemRestriction rstr, CeedMemType mem_type, const CeedInt8 **curl_orients) { 685 Ceed ceed; 686 CeedElemRestriction_Ref *impl; 687 688 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 689 CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed)); 690 691 CeedCheck(mem_type == CEED_MEM_HOST, ceed, CEED_ERROR_BACKEND, "Can only provide to HOST memory"); 692 693 *curl_orients = impl->curl_orients; 694 return CEED_ERROR_SUCCESS; 695 } 696 697 //------------------------------------------------------------------------------ 698 // ElemRestriction Destroy 699 //------------------------------------------------------------------------------ 700 static int CeedElemRestrictionDestroy_Ref(CeedElemRestriction rstr) { 701 CeedElemRestriction_Ref *impl; 702 703 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 704 CeedCallBackend(CeedFree(&impl->offsets_allocated)); 705 CeedCallBackend(CeedFree(&impl->orients_allocated)); 706 CeedCallBackend(CeedFree(&impl->curl_orients_allocated)); 707 CeedCallBackend(CeedFree(&impl)); 708 return CEED_ERROR_SUCCESS; 709 } 710 711 //------------------------------------------------------------------------------ 712 // ElemRestriction Create 713 //------------------------------------------------------------------------------ 714 int CeedElemRestrictionCreate_Ref(CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, const bool *orients, 715 const CeedInt8 *curl_orients, CeedElemRestriction rstr) { 716 Ceed ceed; 717 CeedInt num_elem, elem_size, num_block, block_size, num_comp, comp_stride, num_points = 0, num_offsets; 718 CeedRestrictionType rstr_type; 719 CeedElemRestriction_Ref *impl; 720 721 CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed)); 722 CeedCallBackend(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 723 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr, &elem_size)); 724 CeedCallBackend(CeedElemRestrictionGetNumBlocks(rstr, &num_block)); 725 CeedCallBackend(CeedElemRestrictionGetBlockSize(rstr, &block_size)); 726 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 727 CeedCallBackend(CeedElemRestrictionGetCompStride(rstr, &comp_stride)); 728 CeedInt layout[3] = {1, elem_size, elem_size * num_comp}; 729 730 CeedCheck(mem_type == CEED_MEM_HOST, ceed, CEED_ERROR_BACKEND, "Only MemType = HOST supported"); 731 CeedCallBackend(CeedCalloc(1, &impl)); 732 733 // Offsets data 734 CeedCallBackend(CeedElemRestrictionGetType(rstr, &rstr_type)); 735 if (rstr_type != CEED_RESTRICTION_STRIDED) { 736 const char *resource; 737 738 // Check indices for ref or memcheck backends 739 { 740 Ceed current = ceed, parent = NULL; 741 742 CeedCallBackend(CeedGetParent(current, &parent)); 743 while (current != parent) { 744 current = parent; 745 CeedCallBackend(CeedGetParent(current, &parent)); 746 } 747 CeedCallBackend(CeedGetResource(parent, &resource)); 748 } 749 if (!strcmp(resource, "/cpu/self/ref/serial") || !strcmp(resource, "/cpu/self/ref/blocked") || !strcmp(resource, "/cpu/self/memcheck/serial") || 750 !strcmp(resource, "/cpu/self/memcheck/blocked")) { 751 CeedSize l_size; 752 753 CeedCallBackend(CeedElemRestrictionGetLVectorSize(rstr, &l_size)); 754 for (CeedInt i = 0; i < num_elem * elem_size; i++) { 755 CeedCheck(offsets[i] >= 0 && offsets[i] + (num_comp - 1) * comp_stride < l_size, ceed, CEED_ERROR_BACKEND, 756 "Restriction offset %" CeedInt_FMT " (%" CeedInt_FMT ") out of range [0, %" CeedInt_FMT "]", i, offsets[i], l_size); 757 } 758 } 759 760 // Copy data 761 if (rstr_type == CEED_RESTRICTION_POINTS) CeedCallBackend(CeedElemRestrictionGetNumPoints(rstr, &num_points)); 762 num_offsets = rstr_type == CEED_RESTRICTION_POINTS ? (num_elem + 1 + num_points) : (num_elem * elem_size); 763 switch (copy_mode) { 764 case CEED_COPY_VALUES: 765 CeedCallBackend(CeedMalloc(num_offsets, &impl->offsets_allocated)); 766 memcpy(impl->offsets_allocated, offsets, num_offsets * sizeof(offsets[0])); 767 impl->offsets = impl->offsets_allocated; 768 break; 769 case CEED_OWN_POINTER: 770 impl->offsets_allocated = (CeedInt *)offsets; 771 impl->offsets = impl->offsets_allocated; 772 break; 773 case CEED_USE_POINTER: 774 impl->offsets = offsets; 775 } 776 777 // Orientation data 778 if (rstr_type == CEED_RESTRICTION_ORIENTED) { 779 CeedCheck(orients != NULL, ceed, CEED_ERROR_BACKEND, "No orients array provided for oriented restriction"); 780 switch (copy_mode) { 781 case CEED_COPY_VALUES: 782 CeedCallBackend(CeedMalloc(num_offsets, &impl->orients_allocated)); 783 memcpy(impl->orients_allocated, orients, num_offsets * sizeof(orients[0])); 784 impl->orients = impl->orients_allocated; 785 break; 786 case CEED_OWN_POINTER: 787 impl->orients_allocated = (bool *)orients; 788 impl->orients = impl->orients_allocated; 789 break; 790 case CEED_USE_POINTER: 791 impl->orients = orients; 792 } 793 } else if (rstr_type == CEED_RESTRICTION_CURL_ORIENTED) { 794 CeedCheck(curl_orients != NULL, ceed, CEED_ERROR_BACKEND, "No curl_orients array provided for oriented restriction"); 795 switch (copy_mode) { 796 case CEED_COPY_VALUES: 797 CeedCallBackend(CeedMalloc(3 * num_offsets, &impl->curl_orients_allocated)); 798 memcpy(impl->curl_orients_allocated, curl_orients, 3 * num_offsets * sizeof(curl_orients[0])); 799 impl->curl_orients = impl->curl_orients_allocated; 800 break; 801 case CEED_OWN_POINTER: 802 impl->curl_orients_allocated = (CeedInt8 *)curl_orients; 803 impl->curl_orients = impl->curl_orients_allocated; 804 break; 805 case CEED_USE_POINTER: 806 impl->curl_orients = curl_orients; 807 } 808 } 809 } 810 811 CeedCallBackend(CeedElemRestrictionSetData(rstr, impl)); 812 CeedCallBackend(CeedElemRestrictionSetELayout(rstr, layout)); 813 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "Apply", CeedElemRestrictionApply_Ref)); 814 if (rstr_type == CEED_RESTRICTION_POINTS) { 815 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "ApplyAtPointsInElement", CeedElemRestrictionApplyAtPointsInElement_Ref)); 816 } 817 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "ApplyUnsigned", CeedElemRestrictionApplyUnsigned_Ref)); 818 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "ApplyUnoriented", CeedElemRestrictionApplyUnoriented_Ref)); 819 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "ApplyBlock", CeedElemRestrictionApplyBlock_Ref)); 820 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "GetOrientations", CeedElemRestrictionGetOrientations_Ref)); 821 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "GetCurlOrientations", CeedElemRestrictionGetCurlOrientations_Ref)); 822 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "GetOffsets", CeedElemRestrictionGetOffsets_Ref)); 823 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "Destroy", CeedElemRestrictionDestroy_Ref)); 824 825 // Set apply function based upon num_comp, block_size, and comp_stride 826 CeedInt index = -1; 827 828 if (block_size < 10) index = 100 * num_comp + 10 * block_size + (comp_stride == 1); 829 switch (index) { 830 case 110: 831 impl->Apply = CeedElemRestrictionApply_Ref_110; 832 break; 833 case 111: 834 impl->Apply = CeedElemRestrictionApply_Ref_111; 835 break; 836 case 180: 837 impl->Apply = CeedElemRestrictionApply_Ref_180; 838 break; 839 case 181: 840 impl->Apply = CeedElemRestrictionApply_Ref_181; 841 break; 842 case 310: 843 impl->Apply = CeedElemRestrictionApply_Ref_310; 844 break; 845 case 311: 846 impl->Apply = CeedElemRestrictionApply_Ref_311; 847 break; 848 case 380: 849 impl->Apply = CeedElemRestrictionApply_Ref_380; 850 break; 851 case 381: 852 impl->Apply = CeedElemRestrictionApply_Ref_381; 853 break; 854 // LCOV_EXCL_START 855 case 510: 856 impl->Apply = CeedElemRestrictionApply_Ref_510; 857 break; 858 // LCOV_EXCL_STOP 859 case 511: 860 impl->Apply = CeedElemRestrictionApply_Ref_511; 861 break; 862 // LCOV_EXCL_START 863 case 580: 864 impl->Apply = CeedElemRestrictionApply_Ref_580; 865 break; 866 // LCOV_EXCL_STOP 867 case 581: 868 impl->Apply = CeedElemRestrictionApply_Ref_581; 869 break; 870 default: 871 impl->Apply = CeedElemRestrictionApply_Ref_Core; 872 break; 873 } 874 return CEED_ERROR_SUCCESS; 875 } 876 877 //------------------------------------------------------------------------------ 878