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 <string.h> 12 13 #include "ceed-ref.h" 14 15 //------------------------------------------------------------------------------ 16 // Core ElemRestriction Apply Code 17 //------------------------------------------------------------------------------ 18 static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, 19 CeedInt start, CeedInt stop, bool use_orients, CeedTransposeMode t_mode, CeedVector u, 20 CeedVector v, CeedRequest *request) { 21 CeedElemRestriction_Ref *impl; 22 CeedCallBackend(CeedElemRestrictionGetData(r, &impl)); 23 const CeedScalar *uu; 24 CeedScalar *vv; 25 CeedInt num_elem, elem_size, v_offset; 26 CeedCallBackend(CeedElemRestrictionGetNumElements(r, &num_elem)); 27 CeedCallBackend(CeedElemRestrictionGetElementSize(r, &elem_size)); 28 v_offset = start * blk_size * elem_size * num_comp; 29 bool is_oriented, is_curl_oriented; 30 is_oriented = (impl->orients != NULL); 31 is_curl_oriented = (impl->curl_orients != NULL); 32 33 CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_HOST, &uu)); 34 if (t_mode == CEED_TRANSPOSE) { 35 // Sum into for transpose mode, e-vec to l-vec 36 CeedCallBackend(CeedVectorGetArray(v, CEED_MEM_HOST, &vv)); 37 } else { 38 // Overwrite for notranspose mode, l-vec to e-vec 39 CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_HOST, &vv)); 40 } 41 // Restriction from L-vector to E-vector 42 // Perform: v = r * u 43 if (t_mode == CEED_NOTRANSPOSE) { 44 // No offsets provided, Identity Restriction 45 if (!impl->offsets) { 46 bool has_backend_strides; 47 CeedCallBackend(CeedElemRestrictionHasBackendStrides(r, &has_backend_strides)); 48 if (has_backend_strides) { 49 // CPU backend strides are {1, elem_size, elem_size*num_comp} 50 // This if branch is left separate to allow better inlining 51 for (CeedInt e = start * blk_size; e < stop * blk_size; e += blk_size) { 52 CeedPragmaSIMD for (CeedInt k = 0; k < num_comp; k++) { 53 CeedPragmaSIMD for (CeedInt n = 0; n < elem_size; n++) { 54 CeedPragmaSIMD for (CeedInt j = 0; j < blk_size; j++) { 55 vv[e * elem_size * num_comp + (k * elem_size + n) * blk_size + j - v_offset] = 56 uu[n + k * elem_size + CeedIntMin(e + j, num_elem - 1) * elem_size * num_comp]; 57 } 58 } 59 } 60 } 61 } else { 62 // User provided strides 63 CeedInt strides[3]; 64 CeedCallBackend(CeedElemRestrictionGetStrides(r, &strides)); 65 for (CeedInt e = start * blk_size; e < stop * blk_size; e += blk_size) { 66 CeedPragmaSIMD for (CeedInt k = 0; k < num_comp; k++) { 67 CeedPragmaSIMD for (CeedInt n = 0; n < elem_size; n++) { 68 CeedPragmaSIMD for (CeedInt j = 0; j < blk_size; j++) { 69 vv[e * elem_size * num_comp + (k * elem_size + n) * blk_size + j - v_offset] = 70 uu[n * strides[0] + k * strides[1] + CeedIntMin(e + j, num_elem - 1) * strides[2]]; 71 } 72 } 73 } 74 } 75 } 76 } else { 77 // Offsets provided, standard or blocked restriction 78 // vv has shape [elem_size, num_comp, num_elem], row-major 79 // uu has shape [nnodes, num_comp] 80 for (CeedInt e = start * blk_size; e < stop * blk_size; e += blk_size) { 81 CeedPragmaSIMD for (CeedInt k = 0; k < num_comp; k++) { 82 CeedPragmaSIMD for (CeedInt i = 0; i < elem_size * blk_size; i++) { 83 if (!use_orients || (!is_oriented && !is_curl_oriented)) { 84 // Unsigned restriction 85 vv[elem_size * (k * blk_size + num_comp * e) + i - v_offset] = uu[impl->offsets[i + elem_size * e] + k * comp_stride]; 86 } else if (!is_curl_oriented) { 87 // Signed restriction 88 vv[elem_size * (k * blk_size + num_comp * e) + i - v_offset] = 89 uu[impl->offsets[i + elem_size * e] + k * comp_stride] * (impl->orients[i + elem_size * e] ? -1.0 : 1.0); 90 } else { 91 // Restriction with tridiagonal transformation 92 CeedInt ii = i % elem_size; 93 if (ii == 0) { 94 vv[elem_size * (k * blk_size + num_comp * e) + i - v_offset] = 95 uu[impl->offsets[i + 0 + elem_size * e] + k * comp_stride] * impl->curl_orients[3 * i + 1 + 3 * elem_size * e] + 96 uu[impl->offsets[i + 1 + elem_size * e] + k * comp_stride] * impl->curl_orients[3 * i + 2 + 3 * elem_size * e]; 97 } else if (ii == elem_size - 1) { 98 vv[elem_size * (k * blk_size + num_comp * e) + i - v_offset] = 99 uu[impl->offsets[i - 1 + elem_size * e] + k * comp_stride] * impl->curl_orients[3 * i + 0 + 3 * elem_size * e] + 100 uu[impl->offsets[i + 0 + elem_size * e] + k * comp_stride] * impl->curl_orients[3 * i + 1 + 3 * elem_size * e]; 101 } else { 102 vv[elem_size * (k * blk_size + num_comp * e) + i - v_offset] = 103 uu[impl->offsets[i - 1 + elem_size * e] + k * comp_stride] * impl->curl_orients[3 * i + 0 + 3 * elem_size * e] + 104 uu[impl->offsets[i + 0 + elem_size * e] + k * comp_stride] * impl->curl_orients[3 * i + 1 + 3 * elem_size * e] + 105 uu[impl->offsets[i + 1 + elem_size * e] + k * comp_stride] * impl->curl_orients[3 * i + 2 + 3 * elem_size * e]; 106 } 107 } 108 } 109 } 110 } 111 } 112 } else { 113 // Restriction from E-vector to L-vector 114 // Performing v += r^T * u 115 // No offsets provided, Identity Restriction 116 if (!impl->offsets) { 117 bool has_backend_strides; 118 CeedCallBackend(CeedElemRestrictionHasBackendStrides(r, &has_backend_strides)); 119 if (has_backend_strides) { 120 // CPU backend strides are {1, elem_size, elem_size*num_comp} 121 // This if brach is left separate to allow better inlining 122 for (CeedInt e = start * blk_size; e < stop * blk_size; e += blk_size) { 123 CeedPragmaSIMD for (CeedInt k = 0; k < num_comp; k++) { 124 CeedPragmaSIMD for (CeedInt n = 0; n < elem_size; n++) { 125 CeedPragmaSIMD for (CeedInt j = 0; j < CeedIntMin(blk_size, num_elem - e); j++) { 126 vv[n + k * elem_size + (e + j) * elem_size * num_comp] += 127 uu[e * elem_size * num_comp + (k * elem_size + n) * blk_size + j - v_offset]; 128 } 129 } 130 } 131 } 132 } else { 133 // User provided strides 134 CeedInt strides[3]; 135 CeedCallBackend(CeedElemRestrictionGetStrides(r, &strides)); 136 for (CeedInt e = start * blk_size; e < stop * blk_size; e += blk_size) { 137 CeedPragmaSIMD for (CeedInt k = 0; k < num_comp; k++) { 138 CeedPragmaSIMD for (CeedInt n = 0; n < elem_size; n++) { 139 CeedPragmaSIMD for (CeedInt j = 0; j < CeedIntMin(blk_size, num_elem - e); j++) { 140 vv[n * strides[0] + k * strides[1] + (e + j) * strides[2]] += 141 uu[e * elem_size * num_comp + (k * elem_size + n) * blk_size + j - v_offset]; 142 } 143 } 144 } 145 } 146 } 147 } else { 148 // Offsets provided, standard or blocked restriction 149 // uu has shape [elem_size, num_comp, num_elem] 150 // vv has shape [nnodes, num_comp] 151 for (CeedInt e = start * blk_size; e < stop * blk_size; e += blk_size) { 152 for (CeedInt k = 0; k < num_comp; k++) { 153 for (CeedInt i = 0; i < elem_size * blk_size; i += blk_size) { 154 // Iteration bound set to discard padding elements 155 for (CeedInt j = i; j < i + CeedIntMin(blk_size, num_elem - e); j++) { 156 if (!use_orients || (!is_oriented && !is_curl_oriented)) { 157 // Unsigned restriction 158 vv[impl->offsets[j + e * elem_size] + k * comp_stride] += uu[elem_size * (k * blk_size + num_comp * e) + j - v_offset]; 159 } else if (!is_curl_oriented) { 160 // Signed restriction 161 vv[impl->offsets[j + e * elem_size] + k * comp_stride] += 162 uu[elem_size * (k * blk_size + num_comp * e) + j - v_offset] * (impl->orients[j + e * elem_size] ? -1.0 : 1.0); 163 } else { 164 // Restriction with tridiagonal transformation 165 CeedInt jj = j % elem_size; 166 if (jj == 0) { 167 vv[impl->offsets[j + e * elem_size] + k * comp_stride] += 168 uu[elem_size * (k * blk_size + num_comp * e) + j + 0 - v_offset] * impl->curl_orients[(j + 0) * 3 + 1 + e * 3 * elem_size] + 169 uu[elem_size * (k * blk_size + num_comp * e) + j + 1 - v_offset] * impl->curl_orients[(j + 1) * 3 + 0 + e * 3 * elem_size]; 170 } else if (jj == elem_size - 1) { 171 vv[impl->offsets[j + e * elem_size] + k * comp_stride] += 172 uu[elem_size * (k * blk_size + num_comp * e) + j - 1 - v_offset] * impl->curl_orients[(j - 1) * 3 + 2 + e * 3 * elem_size] + 173 uu[elem_size * (k * blk_size + num_comp * e) + j + 0 - v_offset] * impl->curl_orients[(j + 0) * 3 + 1 + e * 3 * elem_size]; 174 } else { 175 vv[impl->offsets[j + e * elem_size] + k * comp_stride] += 176 uu[elem_size * (k * blk_size + num_comp * e) + j - 1 - v_offset] * impl->curl_orients[(j - 1) * 3 + 2 + e * 3 * elem_size] + 177 uu[elem_size * (k * blk_size + num_comp * e) + j + 0 - v_offset] * impl->curl_orients[(j + 0) * 3 + 1 + e * 3 * elem_size] + 178 uu[elem_size * (k * blk_size + num_comp * e) + j + 1 - v_offset] * impl->curl_orients[(j + 1) * 3 + 0 + e * 3 * elem_size]; 179 } 180 } 181 } 182 } 183 } 184 } 185 } 186 } 187 CeedCallBackend(CeedVectorRestoreArrayRead(u, &uu)); 188 CeedCallBackend(CeedVectorRestoreArray(v, &vv)); 189 if (request != CEED_REQUEST_IMMEDIATE && request != CEED_REQUEST_ORDERED) *request = NULL; 190 return CEED_ERROR_SUCCESS; 191 } 192 193 //------------------------------------------------------------------------------ 194 // ElemRestriction Apply - Common Sizes 195 //------------------------------------------------------------------------------ 196 static int CeedElemRestrictionApply_Ref_110(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, 197 CeedInt start, CeedInt stop, bool use_orients, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 198 CeedRequest *request) { 199 return CeedElemRestrictionApply_Ref_Core(r, 1, 1, comp_stride, start, stop, use_orients, t_mode, u, v, request); 200 } 201 202 static int CeedElemRestrictionApply_Ref_111(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, 203 CeedInt start, CeedInt stop, bool use_orients, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 204 CeedRequest *request) { 205 return CeedElemRestrictionApply_Ref_Core(r, 1, 1, 1, start, stop, use_orients, t_mode, u, v, request); 206 } 207 208 static int CeedElemRestrictionApply_Ref_180(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, 209 CeedInt start, CeedInt stop, bool use_orients, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 210 CeedRequest *request) { 211 return CeedElemRestrictionApply_Ref_Core(r, 1, 8, comp_stride, start, stop, use_orients, t_mode, u, v, request); 212 } 213 214 static int CeedElemRestrictionApply_Ref_181(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, 215 CeedInt start, CeedInt stop, bool use_orients, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 216 CeedRequest *request) { 217 return CeedElemRestrictionApply_Ref_Core(r, 1, 8, 1, start, stop, use_orients, t_mode, u, v, request); 218 } 219 220 static int CeedElemRestrictionApply_Ref_310(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, 221 CeedInt start, CeedInt stop, bool use_orients, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 222 CeedRequest *request) { 223 return CeedElemRestrictionApply_Ref_Core(r, 3, 1, comp_stride, start, stop, use_orients, t_mode, u, v, request); 224 } 225 226 static int CeedElemRestrictionApply_Ref_311(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, 227 CeedInt start, CeedInt stop, bool use_orients, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 228 CeedRequest *request) { 229 return CeedElemRestrictionApply_Ref_Core(r, 3, 1, 1, start, stop, use_orients, t_mode, u, v, request); 230 } 231 232 static int CeedElemRestrictionApply_Ref_380(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, 233 CeedInt start, CeedInt stop, bool use_orients, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 234 CeedRequest *request) { 235 return CeedElemRestrictionApply_Ref_Core(r, 3, 8, comp_stride, start, stop, use_orients, t_mode, u, v, request); 236 } 237 238 static int CeedElemRestrictionApply_Ref_381(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, 239 CeedInt start, CeedInt stop, bool use_orients, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 240 CeedRequest *request) { 241 return CeedElemRestrictionApply_Ref_Core(r, 3, 8, 1, start, stop, use_orients, t_mode, u, v, request); 242 } 243 244 // LCOV_EXCL_START 245 static int CeedElemRestrictionApply_Ref_510(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, 246 CeedInt start, CeedInt stop, bool use_orients, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 247 CeedRequest *request) { 248 return CeedElemRestrictionApply_Ref_Core(r, 5, 1, comp_stride, start, stop, use_orients, t_mode, u, v, request); 249 } 250 // LCOV_EXCL_STOP 251 252 static int CeedElemRestrictionApply_Ref_511(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, 253 CeedInt start, CeedInt stop, bool use_orients, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 254 CeedRequest *request) { 255 return CeedElemRestrictionApply_Ref_Core(r, 5, 1, 1, start, stop, use_orients, t_mode, u, v, request); 256 } 257 258 // LCOV_EXCL_START 259 static int CeedElemRestrictionApply_Ref_580(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, 260 CeedInt start, CeedInt stop, bool use_orients, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 261 CeedRequest *request) { 262 return CeedElemRestrictionApply_Ref_Core(r, 5, 8, comp_stride, start, stop, use_orients, t_mode, u, v, request); 263 } 264 // LCOV_EXCL_STOP 265 266 static int CeedElemRestrictionApply_Ref_581(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, 267 CeedInt start, CeedInt stop, bool use_orients, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 268 CeedRequest *request) { 269 return CeedElemRestrictionApply_Ref_Core(r, 5, 8, 1, start, stop, use_orients, t_mode, u, v, request); 270 } 271 272 //------------------------------------------------------------------------------ 273 // ElemRestriction Apply 274 //------------------------------------------------------------------------------ 275 static int CeedElemRestrictionApply_Ref(CeedElemRestriction r, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { 276 CeedInt num_blk, blk_size, num_comp, comp_stride; 277 CeedCallBackend(CeedElemRestrictionGetNumBlocks(r, &num_blk)); 278 CeedCallBackend(CeedElemRestrictionGetBlockSize(r, &blk_size)); 279 CeedCallBackend(CeedElemRestrictionGetNumComponents(r, &num_comp)); 280 CeedCallBackend(CeedElemRestrictionGetCompStride(r, &comp_stride)); 281 CeedElemRestriction_Ref *impl; 282 CeedCallBackend(CeedElemRestrictionGetData(r, &impl)); 283 284 return impl->Apply(r, num_comp, blk_size, comp_stride, 0, num_blk, true, t_mode, u, v, request); 285 } 286 287 //------------------------------------------------------------------------------ 288 // ElemRestriction Apply Unsigned 289 //------------------------------------------------------------------------------ 290 static int CeedElemRestrictionApplyUnsigned_Ref(CeedElemRestriction r, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { 291 CeedInt num_blk, blk_size, num_comp, comp_stride; 292 CeedCallBackend(CeedElemRestrictionGetNumBlocks(r, &num_blk)); 293 CeedCallBackend(CeedElemRestrictionGetBlockSize(r, &blk_size)); 294 CeedCallBackend(CeedElemRestrictionGetNumComponents(r, &num_comp)); 295 CeedCallBackend(CeedElemRestrictionGetCompStride(r, &comp_stride)); 296 CeedElemRestriction_Ref *impl; 297 CeedCallBackend(CeedElemRestrictionGetData(r, &impl)); 298 299 return impl->Apply(r, num_comp, blk_size, comp_stride, 0, num_blk, false, t_mode, u, v, request); 300 } 301 302 //------------------------------------------------------------------------------ 303 // ElemRestriction Apply Block 304 //------------------------------------------------------------------------------ 305 static int CeedElemRestrictionApplyBlock_Ref(CeedElemRestriction r, CeedInt block, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 306 CeedRequest *request) { 307 CeedInt blk_size, num_comp, comp_stride; 308 CeedCallBackend(CeedElemRestrictionGetBlockSize(r, &blk_size)); 309 CeedCallBackend(CeedElemRestrictionGetNumComponents(r, &num_comp)); 310 CeedCallBackend(CeedElemRestrictionGetCompStride(r, &comp_stride)); 311 CeedElemRestriction_Ref *impl; 312 CeedCallBackend(CeedElemRestrictionGetData(r, &impl)); 313 314 return impl->Apply(r, num_comp, blk_size, comp_stride, block, block + 1, true, t_mode, u, v, request); 315 } 316 317 //------------------------------------------------------------------------------ 318 // ElemRestriction Get Offsets 319 //------------------------------------------------------------------------------ 320 static int CeedElemRestrictionGetOffsets_Ref(CeedElemRestriction rstr, CeedMemType mem_type, const CeedInt **offsets) { 321 CeedElemRestriction_Ref *impl; 322 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 323 Ceed ceed; 324 CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed)); 325 326 CeedCheck(mem_type == CEED_MEM_HOST, ceed, CEED_ERROR_BACKEND, "Can only provide to HOST memory"); 327 328 *offsets = impl->offsets; 329 return CEED_ERROR_SUCCESS; 330 } 331 332 //------------------------------------------------------------------------------ 333 // ElemRestriction Get Orientations 334 //------------------------------------------------------------------------------ 335 static int CeedElemRestrictionGetOrientations_Ref(CeedElemRestriction rstr, CeedMemType mem_type, const bool **orients) { 336 CeedElemRestriction_Ref *impl; 337 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 338 Ceed ceed; 339 CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed)); 340 341 if (mem_type != CEED_MEM_HOST) { 342 // LCOV_EXCL_START 343 return CeedError(ceed, CEED_ERROR_BACKEND, "Can only provide to HOST memory"); 344 // LCOV_EXCL_STOP 345 } 346 347 *orients = impl->orients; 348 return CEED_ERROR_SUCCESS; 349 } 350 351 //------------------------------------------------------------------------------ 352 // ElemRestriction Get Curl-Conforming Orientations 353 //------------------------------------------------------------------------------ 354 static int CeedElemRestrictionGetCurlOrientations_Ref(CeedElemRestriction rstr, CeedMemType mem_type, const CeedInt **curl_orients) { 355 CeedElemRestriction_Ref *impl; 356 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 357 Ceed ceed; 358 CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed)); 359 360 if (mem_type != CEED_MEM_HOST) { 361 // LCOV_EXCL_START 362 return CeedError(ceed, CEED_ERROR_BACKEND, "Can only provide to HOST memory"); 363 // LCOV_EXCL_STOP 364 } 365 366 *curl_orients = impl->curl_orients; 367 return CEED_ERROR_SUCCESS; 368 } 369 370 //------------------------------------------------------------------------------ 371 // ElemRestriction Destroy 372 //------------------------------------------------------------------------------ 373 static int CeedElemRestrictionDestroy_Ref(CeedElemRestriction r) { 374 CeedElemRestriction_Ref *impl; 375 CeedCallBackend(CeedElemRestrictionGetData(r, &impl)); 376 377 CeedCallBackend(CeedFree(&impl->offsets_allocated)); 378 CeedCallBackend(CeedFree(&impl->orients_allocated)); 379 CeedCallBackend(CeedFree(&impl->curl_orients_allocated)); 380 CeedCallBackend(CeedFree(&impl)); 381 return CEED_ERROR_SUCCESS; 382 } 383 384 //------------------------------------------------------------------------------ 385 // ElemRestriction Create 386 //------------------------------------------------------------------------------ 387 int CeedElemRestrictionCreate_Ref(CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, CeedElemRestriction r) { 388 CeedElemRestriction_Ref *impl; 389 CeedInt num_elem, elem_size, num_blk, blk_size, num_comp, comp_stride; 390 CeedCallBackend(CeedElemRestrictionGetNumElements(r, &num_elem)); 391 CeedCallBackend(CeedElemRestrictionGetElementSize(r, &elem_size)); 392 CeedCallBackend(CeedElemRestrictionGetNumBlocks(r, &num_blk)); 393 CeedCallBackend(CeedElemRestrictionGetBlockSize(r, &blk_size)); 394 CeedCallBackend(CeedElemRestrictionGetNumComponents(r, &num_comp)); 395 CeedCallBackend(CeedElemRestrictionGetCompStride(r, &comp_stride)); 396 Ceed ceed; 397 CeedCallBackend(CeedElemRestrictionGetCeed(r, &ceed)); 398 399 CeedCheck(mem_type == CEED_MEM_HOST, ceed, CEED_ERROR_BACKEND, "Only MemType = HOST supported"); 400 CeedCallBackend(CeedCalloc(1, &impl)); 401 402 // Offsets data 403 bool is_strided; 404 CeedCallBackend(CeedElemRestrictionIsStrided(r, &is_strided)); 405 if (!is_strided) { 406 // Check indices for ref or memcheck backends 407 Ceed parent_ceed = ceed, curr_ceed = NULL; 408 while (parent_ceed != curr_ceed) { 409 curr_ceed = parent_ceed; 410 CeedCallBackend(CeedGetParent(curr_ceed, &parent_ceed)); 411 } 412 const char *resource; 413 CeedCallBackend(CeedGetResource(parent_ceed, &resource)); 414 if (!strcmp(resource, "/cpu/self/ref/serial") || !strcmp(resource, "/cpu/self/ref/blocked") || !strcmp(resource, "/cpu/self/memcheck/serial") || 415 !strcmp(resource, "/cpu/self/memcheck/blocked")) { 416 CeedSize l_size; 417 CeedCallBackend(CeedElemRestrictionGetLVectorSize(r, &l_size)); 418 419 for (CeedInt i = 0; i < num_elem * elem_size; i++) { 420 CeedCheck(offsets[i] >= 0 && offsets[i] + (num_comp - 1) * comp_stride < l_size, ceed, CEED_ERROR_BACKEND, 421 "Restriction offset %" CeedInt_FMT " (%" CeedInt_FMT ") out of range [0, %" CeedInt_FMT "]", i, offsets[i], l_size); 422 } 423 } 424 425 // Copy data 426 switch (copy_mode) { 427 case CEED_COPY_VALUES: 428 CeedCallBackend(CeedMalloc(num_elem * elem_size, &impl->offsets_allocated)); 429 memcpy(impl->offsets_allocated, offsets, num_elem * elem_size * sizeof(offsets[0])); 430 impl->offsets = impl->offsets_allocated; 431 break; 432 case CEED_OWN_POINTER: 433 impl->offsets_allocated = (CeedInt *)offsets; 434 impl->offsets = impl->offsets_allocated; 435 break; 436 case CEED_USE_POINTER: 437 impl->offsets = offsets; 438 } 439 } 440 441 CeedCallBackend(CeedElemRestrictionSetData(r, impl)); 442 CeedInt layout[3] = {1, elem_size, elem_size * num_comp}; 443 CeedCallBackend(CeedElemRestrictionSetELayout(r, layout)); 444 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "Apply", CeedElemRestrictionApply_Ref)); 445 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "ApplyUnsigned", CeedElemRestrictionApplyUnsigned_Ref)); 446 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "ApplyBlock", CeedElemRestrictionApplyBlock_Ref)); 447 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "GetOffsets", CeedElemRestrictionGetOffsets_Ref)); 448 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "GetOrientations", CeedElemRestrictionGetOrientations_Ref)); 449 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "GetCurlOrientations", CeedElemRestrictionGetCurlOrientations_Ref)); 450 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "Destroy", CeedElemRestrictionDestroy_Ref)); 451 452 // Set apply function based upon num_comp, blk_size, and comp_stride 453 CeedInt idx = -1; 454 if (blk_size < 10) idx = 100 * num_comp + 10 * blk_size + (comp_stride == 1); 455 switch (idx) { 456 case 110: 457 impl->Apply = CeedElemRestrictionApply_Ref_110; 458 break; 459 case 111: 460 impl->Apply = CeedElemRestrictionApply_Ref_111; 461 break; 462 case 180: 463 impl->Apply = CeedElemRestrictionApply_Ref_180; 464 break; 465 case 181: 466 impl->Apply = CeedElemRestrictionApply_Ref_181; 467 break; 468 case 310: 469 impl->Apply = CeedElemRestrictionApply_Ref_310; 470 break; 471 case 311: 472 impl->Apply = CeedElemRestrictionApply_Ref_311; 473 break; 474 case 380: 475 impl->Apply = CeedElemRestrictionApply_Ref_380; 476 break; 477 case 381: 478 impl->Apply = CeedElemRestrictionApply_Ref_381; 479 break; 480 // LCOV_EXCL_START 481 case 510: 482 impl->Apply = CeedElemRestrictionApply_Ref_510; 483 break; 484 // LCOV_EXCL_STOP 485 case 511: 486 impl->Apply = CeedElemRestrictionApply_Ref_511; 487 break; 488 // LCOV_EXCL_START 489 case 580: 490 impl->Apply = CeedElemRestrictionApply_Ref_580; 491 break; 492 // LCOV_EXCL_STOP 493 case 581: 494 impl->Apply = CeedElemRestrictionApply_Ref_581; 495 break; 496 default: 497 impl->Apply = CeedElemRestrictionApply_Ref_Core; 498 break; 499 } 500 501 return CEED_ERROR_SUCCESS; 502 } 503 504 //------------------------------------------------------------------------------ 505 // ElemRestriction Create Oriented 506 //------------------------------------------------------------------------------ 507 int CeedElemRestrictionCreateOriented_Ref(CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, const bool *orients, 508 CeedElemRestriction r) { 509 // Set up for normal restriction with explicit offsets. This sets up dispatch to 510 // CeedElemRestrictionApply_Ref_* and manages the impl->offsets array copy/allocation. 511 CeedCallBackend(CeedElemRestrictionCreate_Ref(mem_type, copy_mode, offsets, r)); 512 513 CeedElemRestriction_Ref *impl; 514 CeedInt num_elem, elem_size; 515 CeedCallBackend(CeedElemRestrictionGetData(r, &impl)); 516 CeedCallBackend(CeedElemRestrictionGetNumElements(r, &num_elem)); 517 CeedCallBackend(CeedElemRestrictionGetElementSize(r, &elem_size)); 518 519 // Copy data 520 switch (copy_mode) { 521 case CEED_COPY_VALUES: 522 CeedCallBackend(CeedMalloc(num_elem * elem_size, &impl->orients_allocated)); 523 memcpy(impl->orients_allocated, orients, num_elem * elem_size * sizeof(orients[0])); 524 impl->orients = impl->orients_allocated; 525 break; 526 case CEED_OWN_POINTER: 527 impl->orients_allocated = (bool *)orients; 528 impl->orients = impl->orients_allocated; 529 break; 530 case CEED_USE_POINTER: 531 impl->orients = orients; 532 } 533 534 return CEED_ERROR_SUCCESS; 535 } 536 537 //------------------------------------------------------------------------------ 538 // ElemRestriction Create Curl-Conforming Oriented 539 //------------------------------------------------------------------------------ 540 int CeedElemRestrictionCreateCurlOriented_Ref(CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, const CeedInt *curl_orients, 541 CeedElemRestriction r) { 542 // Set up for normal restriction with explicit offsets. This sets up dispatch to 543 // CeedElemRestrictionApply_Ref_* and manages the impl->offsets array copy/allocation. 544 CeedCallBackend(CeedElemRestrictionCreate_Ref(mem_type, copy_mode, offsets, r)); 545 546 CeedElemRestriction_Ref *impl; 547 CeedInt num_elem, elem_size; 548 CeedCallBackend(CeedElemRestrictionGetData(r, &impl)); 549 CeedCallBackend(CeedElemRestrictionGetNumElements(r, &num_elem)); 550 CeedCallBackend(CeedElemRestrictionGetElementSize(r, &elem_size)); 551 552 // Orientation data 553 switch (copy_mode) { 554 case CEED_COPY_VALUES: 555 CeedCallBackend(CeedMalloc(num_elem * 3 * elem_size, &impl->curl_orients_allocated)); 556 memcpy(impl->curl_orients_allocated, curl_orients, num_elem * 3 * elem_size * sizeof(curl_orients[0])); 557 impl->curl_orients = impl->curl_orients_allocated; 558 break; 559 case CEED_OWN_POINTER: 560 impl->curl_orients_allocated = (CeedInt *)curl_orients; 561 impl->curl_orients = impl->curl_orients_allocated; 562 break; 563 case CEED_USE_POINTER: 564 impl->curl_orients = curl_orients; 565 } 566 567 return CEED_ERROR_SUCCESS; 568 } 569 570 //------------------------------------------------------------------------------ 571