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