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