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