1 // Copyright (c) 2017-2024, 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 <stddef.h> 12 #include <stdint.h> 13 14 #include "ceed-ref.h" 15 16 //------------------------------------------------------------------------------ 17 // Setup Input/Output Fields 18 //------------------------------------------------------------------------------ 19 static int CeedOperatorSetupFields_Ref(CeedQFunction qf, CeedOperator op, bool is_input, bool *skip_rstr, CeedInt *e_data_out_indices, 20 bool *apply_add_basis, CeedVector *e_vecs_full, CeedVector *e_vecs, CeedVector *q_vecs, CeedInt start_e, 21 CeedInt num_fields, CeedInt Q) { 22 Ceed ceed; 23 CeedSize e_size, q_size; 24 CeedInt num_comp, size, P; 25 CeedQFunctionField *qf_fields; 26 CeedOperatorField *op_fields; 27 28 { 29 Ceed ceed_parent; 30 31 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 32 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 33 if (ceed_parent) ceed = ceed_parent; 34 } 35 if (is_input) { 36 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 37 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 38 } else { 39 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 40 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 41 } 42 43 // Loop over fields 44 for (CeedInt i = 0; i < num_fields; i++) { 45 CeedEvalMode eval_mode; 46 CeedElemRestriction elem_rstr; 47 CeedBasis basis; 48 49 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 50 if (eval_mode != CEED_EVAL_WEIGHT) { 51 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr)); 52 CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs_full[i + start_e])); 53 } 54 55 switch (eval_mode) { 56 case CEED_EVAL_NONE: 57 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 58 q_size = (CeedSize)Q * size; 59 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 60 break; 61 case CEED_EVAL_INTERP: 62 case CEED_EVAL_GRAD: 63 case CEED_EVAL_DIV: 64 case CEED_EVAL_CURL: 65 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 66 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 67 CeedCallBackend(CeedBasisGetNumNodes(basis, &P)); 68 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 69 e_size = (CeedSize)P * num_comp; 70 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); 71 q_size = (CeedSize)Q * size; 72 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 73 break; 74 case CEED_EVAL_WEIGHT: // Only on input fields 75 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 76 q_size = (CeedSize)Q; 77 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 78 CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i])); 79 break; 80 } 81 } 82 // Drop duplicate restrictions 83 if (is_input) { 84 for (CeedInt i = 0; i < num_fields; i++) { 85 CeedVector vec_i; 86 CeedElemRestriction rstr_i; 87 88 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec_i)); 89 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr_i)); 90 for (CeedInt j = i + 1; j < num_fields; j++) { 91 CeedVector vec_j; 92 CeedElemRestriction rstr_j; 93 94 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[j], &vec_j)); 95 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[j], &rstr_j)); 96 if (vec_i == vec_j && rstr_i == rstr_j) { 97 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &e_vecs[j])); 98 CeedCallBackend(CeedVectorReferenceCopy(e_vecs_full[i + start_e], &e_vecs_full[j + start_e])); 99 skip_rstr[j] = true; 100 } 101 } 102 } 103 } else { 104 for (CeedInt i = num_fields - 1; i >= 0; i--) { 105 CeedVector vec_i; 106 CeedElemRestriction rstr_i; 107 108 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec_i)); 109 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr_i)); 110 for (CeedInt j = i - 1; j >= 0; j--) { 111 CeedVector vec_j; 112 CeedElemRestriction rstr_j; 113 114 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[j], &vec_j)); 115 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[j], &rstr_j)); 116 if (vec_i == vec_j && rstr_i == rstr_j) { 117 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &e_vecs[j])); 118 CeedCallBackend(CeedVectorReferenceCopy(e_vecs_full[i + start_e], &e_vecs_full[j + start_e])); 119 skip_rstr[j] = true; 120 apply_add_basis[i] = true; 121 e_data_out_indices[j] = i; 122 } 123 } 124 } 125 } 126 return CEED_ERROR_SUCCESS; 127 } 128 129 //------------------------------------------------------------------------------ 130 // Setup Operator 131 //------------------------------------------------------------------------------/* 132 static int CeedOperatorSetup_Ref(CeedOperator op) { 133 bool is_setup_done; 134 CeedInt Q, num_input_fields, num_output_fields; 135 CeedQFunctionField *qf_input_fields, *qf_output_fields; 136 CeedQFunction qf; 137 CeedOperatorField *op_input_fields, *op_output_fields; 138 CeedOperator_Ref *impl; 139 140 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 141 if (is_setup_done) return CEED_ERROR_SUCCESS; 142 143 CeedCallBackend(CeedOperatorGetData(op, &impl)); 144 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 145 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 146 CeedCallBackend(CeedQFunctionIsIdentity(qf, &impl->is_identity_qf)); 147 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 148 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 149 150 // Allocate 151 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs_full)); 152 153 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->skip_rstr_in)); 154 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->skip_rstr_out)); 155 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_data_out_indices)); 156 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->apply_add_basis_out)); 157 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states)); 158 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_in)); 159 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_out)); 160 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in)); 161 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out)); 162 163 impl->num_inputs = num_input_fields; 164 impl->num_outputs = num_output_fields; 165 166 // Set up infield and outfield e_vecs and q_vecs 167 // Infields 168 CeedCallBackend(CeedOperatorSetupFields_Ref(qf, op, true, impl->skip_rstr_in, NULL, NULL, impl->e_vecs_full, impl->e_vecs_in, impl->q_vecs_in, 0, 169 num_input_fields, Q)); 170 // Outfields 171 CeedCallBackend(CeedOperatorSetupFields_Ref(qf, op, false, impl->skip_rstr_out, impl->e_data_out_indices, impl->apply_add_basis_out, 172 impl->e_vecs_full, impl->e_vecs_out, impl->q_vecs_out, num_input_fields, num_output_fields, Q)); 173 174 // Identity QFunctions 175 if (impl->is_identity_qf) { 176 CeedEvalMode in_mode, out_mode; 177 CeedQFunctionField *in_fields, *out_fields; 178 179 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &in_fields, NULL, &out_fields)); 180 CeedCallBackend(CeedQFunctionFieldGetEvalMode(in_fields[0], &in_mode)); 181 CeedCallBackend(CeedQFunctionFieldGetEvalMode(out_fields[0], &out_mode)); 182 183 if (in_mode == CEED_EVAL_NONE && out_mode == CEED_EVAL_NONE) { 184 impl->is_identity_rstr_op = true; 185 } else { 186 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->q_vecs_out[0])); 187 } 188 } 189 190 CeedCallBackend(CeedOperatorSetSetupDone(op)); 191 return CEED_ERROR_SUCCESS; 192 } 193 194 //------------------------------------------------------------------------------ 195 // Setup Operator Inputs 196 //------------------------------------------------------------------------------ 197 static inline int CeedOperatorSetupInputs_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 198 CeedVector in_vec, const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], 199 CeedOperator_Ref *impl, CeedRequest *request) { 200 for (CeedInt i = 0; i < num_input_fields; i++) { 201 uint64_t state; 202 CeedEvalMode eval_mode; 203 CeedVector vec; 204 CeedElemRestriction elem_rstr; 205 206 // Get input vector 207 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 208 if (vec == CEED_VECTOR_ACTIVE) { 209 if (skip_active) continue; 210 else vec = in_vec; 211 } 212 213 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 214 // Restrict and Evec 215 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 216 } else { 217 // Restrict 218 CeedCallBackend(CeedVectorGetState(vec, &state)); 219 // Skip restriction if input is unchanged 220 if ((state != impl->input_states[i] || vec == in_vec) && !impl->skip_rstr_in[i]) { 221 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 222 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_NOTRANSPOSE, vec, impl->e_vecs_full[i], request)); 223 } 224 impl->input_states[i] = state; 225 // Get evec 226 CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs_full[i], CEED_MEM_HOST, (const CeedScalar **)&e_data_full[i])); 227 } 228 } 229 return CEED_ERROR_SUCCESS; 230 } 231 232 //------------------------------------------------------------------------------ 233 // Input Basis Action 234 //------------------------------------------------------------------------------ 235 static inline int CeedOperatorInputBasis_Ref(CeedInt e, CeedInt Q, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 236 CeedInt num_input_fields, const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], 237 CeedOperator_Ref *impl) { 238 for (CeedInt i = 0; i < num_input_fields; i++) { 239 CeedInt elem_size, size, num_comp; 240 CeedEvalMode eval_mode; 241 CeedElemRestriction elem_rstr; 242 CeedBasis basis; 243 244 // Skip active input 245 if (skip_active) { 246 CeedVector vec; 247 248 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 249 if (vec == CEED_VECTOR_ACTIVE) continue; 250 } 251 // Get elem_size, eval_mode, size 252 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 253 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 254 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 255 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 256 // Basis action 257 switch (eval_mode) { 258 case CEED_EVAL_NONE: 259 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][(CeedSize)e * Q * size])); 260 break; 261 case CEED_EVAL_INTERP: 262 case CEED_EVAL_GRAD: 263 case CEED_EVAL_DIV: 264 case CEED_EVAL_CURL: 265 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 266 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 267 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][(CeedSize)e * elem_size * num_comp])); 268 CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, eval_mode, impl->e_vecs_in[i], impl->q_vecs_in[i])); 269 break; 270 case CEED_EVAL_WEIGHT: 271 break; // No action 272 } 273 } 274 return CEED_ERROR_SUCCESS; 275 } 276 277 //------------------------------------------------------------------------------ 278 // Output Basis Action 279 //------------------------------------------------------------------------------ 280 static inline int CeedOperatorOutputBasis_Ref(CeedInt e, CeedInt Q, CeedQFunctionField *qf_output_fields, CeedOperatorField *op_output_fields, 281 CeedInt num_input_fields, CeedInt num_output_fields, bool *apply_add_basis, CeedOperator op, 282 CeedScalar *e_data_full[2 * CEED_FIELD_MAX], CeedOperator_Ref *impl) { 283 for (CeedInt i = 0; i < num_output_fields; i++) { 284 CeedInt elem_size, num_comp; 285 CeedEvalMode eval_mode; 286 CeedElemRestriction elem_rstr; 287 CeedBasis basis; 288 289 // Get elem_size, eval_mode 290 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 291 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 292 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 293 // Basis action 294 switch (eval_mode) { 295 case CEED_EVAL_NONE: 296 break; // No action 297 case CEED_EVAL_INTERP: 298 case CEED_EVAL_GRAD: 299 case CEED_EVAL_DIV: 300 case CEED_EVAL_CURL: 301 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 302 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 303 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, 304 &e_data_full[i + num_input_fields][(CeedSize)e * elem_size * num_comp])); 305 if (apply_add_basis[i]) { 306 CeedCallBackend(CeedBasisApplyAdd(basis, 1, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs_out[i])); 307 } else { 308 CeedCallBackend(CeedBasisApply(basis, 1, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs_out[i])); 309 } 310 break; 311 // LCOV_EXCL_START 312 case CEED_EVAL_WEIGHT: { 313 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 314 // LCOV_EXCL_STOP 315 } 316 } 317 } 318 return CEED_ERROR_SUCCESS; 319 } 320 321 //------------------------------------------------------------------------------ 322 // Restore Input Vectors 323 //------------------------------------------------------------------------------ 324 static inline int CeedOperatorRestoreInputs_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 325 const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], CeedOperator_Ref *impl) { 326 for (CeedInt i = 0; i < num_input_fields; i++) { 327 CeedEvalMode eval_mode; 328 329 // Skip active inputs 330 if (skip_active) { 331 CeedVector vec; 332 333 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 334 if (vec == CEED_VECTOR_ACTIVE) continue; 335 } 336 // Restore input 337 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 338 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 339 } else { 340 CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs_full[i], (const CeedScalar **)&e_data_full[i])); 341 } 342 } 343 return CEED_ERROR_SUCCESS; 344 } 345 346 //------------------------------------------------------------------------------ 347 // Operator Apply 348 //------------------------------------------------------------------------------ 349 static int CeedOperatorApplyAdd_Ref(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) { 350 CeedInt Q, num_elem, num_input_fields, num_output_fields, size; 351 CeedEvalMode eval_mode; 352 CeedScalar *e_data_full[2 * CEED_FIELD_MAX] = {NULL}; 353 CeedQFunctionField *qf_input_fields, *qf_output_fields; 354 CeedQFunction qf; 355 CeedOperatorField *op_input_fields, *op_output_fields; 356 CeedOperator_Ref *impl; 357 358 CeedCallBackend(CeedOperatorGetData(op, &impl)); 359 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 360 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 361 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 362 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 363 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 364 365 // Setup 366 CeedCallBackend(CeedOperatorSetup_Ref(op)); 367 368 // Restriction only operator 369 if (impl->is_identity_rstr_op) { 370 CeedElemRestriction elem_rstr; 371 372 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[0], &elem_rstr)); 373 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_full[0], request)); 374 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[0], &elem_rstr)); 375 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs_full[0], out_vec, request)); 376 return CEED_ERROR_SUCCESS; 377 } 378 379 // Input Evecs and Restriction 380 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, in_vec, false, e_data_full, impl, request)); 381 382 // Output Evecs 383 for (CeedInt i = num_output_fields - 1; i >= 0; i--) { 384 if (impl->skip_rstr_out[i]) { 385 e_data_full[i + num_input_fields] = e_data_full[impl->e_data_out_indices[i] + num_input_fields]; 386 } else { 387 CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs_full[i + impl->num_inputs], CEED_MEM_HOST, &e_data_full[i + num_input_fields])); 388 } 389 } 390 391 // Loop through elements 392 for (CeedInt e = 0; e < num_elem; e++) { 393 // Output pointers 394 for (CeedInt i = 0; i < num_output_fields; i++) { 395 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 396 if (eval_mode == CEED_EVAL_NONE) { 397 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size)); 398 CeedCallBackend( 399 CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][(CeedSize)e * Q * size])); 400 } 401 } 402 403 // Input basis apply 404 CeedCallBackend(CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, num_input_fields, false, e_data_full, impl)); 405 406 // Q function 407 if (!impl->is_identity_qf) { 408 CeedCallBackend(CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out)); 409 } 410 411 // Output basis apply 412 CeedCallBackend(CeedOperatorOutputBasis_Ref(e, Q, qf_output_fields, op_output_fields, num_input_fields, num_output_fields, 413 impl->apply_add_basis_out, op, e_data_full, impl)); 414 } 415 416 // Output restriction 417 for (CeedInt i = 0; i < num_output_fields; i++) { 418 CeedVector vec; 419 CeedElemRestriction elem_rstr; 420 421 if (impl->skip_rstr_out[i]) continue; 422 // Restore Evec 423 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_full[i + impl->num_inputs], &e_data_full[i + num_input_fields])); 424 // Get output vector 425 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 426 // Active 427 if (vec == CEED_VECTOR_ACTIVE) vec = out_vec; 428 // Restrict 429 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 430 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs_full[i + impl->num_inputs], vec, request)); 431 } 432 433 // Restore input arrays 434 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, false, e_data_full, impl)); 435 return CEED_ERROR_SUCCESS; 436 } 437 438 //------------------------------------------------------------------------------ 439 // Core code for assembling linear QFunction 440 //------------------------------------------------------------------------------ 441 static inline int CeedOperatorLinearAssembleQFunctionCore_Ref(CeedOperator op, bool build_objects, CeedVector *assembled, CeedElemRestriction *rstr, 442 CeedRequest *request) { 443 Ceed ceed, ceed_parent; 444 CeedInt qf_size_in, qf_size_out, Q, num_elem, num_input_fields, num_output_fields; 445 CeedScalar *assembled_array, *e_data_full[2 * CEED_FIELD_MAX] = {NULL}; 446 CeedQFunctionField *qf_input_fields, *qf_output_fields; 447 CeedQFunction qf; 448 CeedOperatorField *op_input_fields, *op_output_fields; 449 CeedOperator_Ref *impl; 450 451 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 452 CeedCallBackend(CeedOperatorGetFallbackParentCeed(op, &ceed_parent)); 453 CeedCallBackend(CeedOperatorGetData(op, &impl)); 454 qf_size_in = impl->qf_size_in; 455 qf_size_out = impl->qf_size_out; 456 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 457 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 458 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 459 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 460 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 461 462 // Setup 463 CeedCallBackend(CeedOperatorSetup_Ref(op)); 464 465 // Check for restriction only operator 466 CeedCheck(!impl->is_identity_rstr_op, ceed, CEED_ERROR_BACKEND, "Assembling restriction only operators is not supported"); 467 468 // Input Evecs and Restriction 469 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data_full, impl, request)); 470 471 // Count number of active input fields 472 if (qf_size_in == 0) { 473 for (CeedInt i = 0; i < num_input_fields; i++) { 474 CeedInt field_size; 475 CeedVector vec; 476 477 // Get input vector 478 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 479 // Check if active input 480 if (vec == CEED_VECTOR_ACTIVE) { 481 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size)); 482 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 483 qf_size_in += field_size; 484 } 485 } 486 CeedCheck(qf_size_in > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 487 impl->qf_size_in = qf_size_in; 488 } 489 490 // Count number of active output fields 491 if (qf_size_out == 0) { 492 for (CeedInt i = 0; i < num_output_fields; i++) { 493 CeedInt field_size; 494 CeedVector vec; 495 496 // Get output vector 497 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 498 // Check if active output 499 if (vec == CEED_VECTOR_ACTIVE) { 500 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &field_size)); 501 qf_size_out += field_size; 502 } 503 } 504 CeedCheck(qf_size_out > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 505 impl->qf_size_out = qf_size_out; 506 } 507 508 // Build objects if needed 509 if (build_objects) { 510 const CeedSize l_size = (CeedSize)num_elem * Q * qf_size_in * qf_size_out; 511 CeedInt strides[3] = {1, Q, qf_size_in * qf_size_out * Q}; /* *NOPAD* */ 512 513 // Create output restriction 514 CeedCallBackend(CeedElemRestrictionCreateStrided(ceed_parent, num_elem, Q, qf_size_in * qf_size_out, 515 (CeedSize)qf_size_in * (CeedSize)qf_size_out * (CeedSize)num_elem * (CeedSize)Q, strides, rstr)); 516 // Create assembled vector 517 CeedCallBackend(CeedVectorCreate(ceed_parent, l_size, assembled)); 518 } 519 // Clear output vector 520 CeedCallBackend(CeedVectorSetValue(*assembled, 0.0)); 521 CeedCallBackend(CeedVectorGetArray(*assembled, CEED_MEM_HOST, &assembled_array)); 522 523 // Loop through elements 524 for (CeedInt e = 0; e < num_elem; e++) { 525 // Input basis apply 526 CeedCallBackend(CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, num_input_fields, true, e_data_full, impl)); 527 528 // Assemble QFunction 529 530 for (CeedInt i = 0; i < num_input_fields; i++) { 531 CeedInt field_size; 532 CeedVector vec; 533 534 // Set Inputs 535 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 536 if (vec != CEED_VECTOR_ACTIVE) continue; 537 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size)); 538 for (CeedInt field = 0; field < field_size; field++) { 539 // Set current portion of input to 1.0 540 { 541 CeedScalar *array; 542 543 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array)); 544 for (CeedInt j = 0; j < Q; j++) array[field * Q + j] = 1.0; 545 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array)); 546 } 547 548 if (!impl->is_identity_qf) { 549 // Set Outputs 550 for (CeedInt out = 0; out < num_output_fields; out++) { 551 CeedVector vec; 552 553 // Get output vector 554 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 555 // Check if active output 556 if (vec == CEED_VECTOR_ACTIVE) { 557 CeedInt field_size; 558 559 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, CEED_USE_POINTER, assembled_array)); 560 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &field_size)); 561 assembled_array += field_size * Q; // Advance the pointer by the size of the output 562 } 563 } 564 // Apply QFunction 565 CeedCallBackend(CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out)); 566 } else { 567 CeedInt field_size; 568 const CeedScalar *array; 569 570 // Copy Identity Outputs 571 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[0], &field_size)); 572 CeedCallBackend(CeedVectorGetArrayRead(impl->q_vecs_out[0], CEED_MEM_HOST, &array)); 573 for (CeedInt j = 0; j < field_size * Q; j++) assembled_array[j] = array[j]; 574 CeedCallBackend(CeedVectorRestoreArrayRead(impl->q_vecs_out[0], &array)); 575 assembled_array += field_size * Q; 576 } 577 // Reset input to 0.0 578 { 579 CeedScalar *array; 580 581 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array)); 582 for (CeedInt j = 0; j < Q; j++) array[field * Q + j] = 0.0; 583 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array)); 584 } 585 } 586 } 587 } 588 589 // Un-set output Qvecs to prevent accidental overwrite of Assembled 590 if (!impl->is_identity_qf) { 591 for (CeedInt out = 0; out < num_output_fields; out++) { 592 CeedVector vec; 593 594 // Get output vector 595 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 596 // Check if active output 597 if (vec == CEED_VECTOR_ACTIVE && num_elem > 0) { 598 CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_HOST, NULL)); 599 } 600 } 601 } 602 603 // Restore input arrays 604 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data_full, impl)); 605 606 // Restore output 607 CeedCallBackend(CeedVectorRestoreArray(*assembled, &assembled_array)); 608 return CEED_ERROR_SUCCESS; 609 } 610 611 //------------------------------------------------------------------------------ 612 // Assemble Linear QFunction 613 //------------------------------------------------------------------------------ 614 static int CeedOperatorLinearAssembleQFunction_Ref(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 615 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, true, assembled, rstr, request); 616 } 617 618 //------------------------------------------------------------------------------ 619 // Update Assembled Linear QFunction 620 //------------------------------------------------------------------------------ 621 static int CeedOperatorLinearAssembleQFunctionUpdate_Ref(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) { 622 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, false, &assembled, &rstr, request); 623 } 624 625 //------------------------------------------------------------------------------ 626 // Setup Input/Output Fields 627 //------------------------------------------------------------------------------ 628 static int CeedOperatorSetupFieldsAtPoints_Ref(CeedQFunction qf, CeedOperator op, bool is_input, bool *skip_rstr, bool *apply_add_basis, 629 CeedVector *e_vecs_full, CeedVector *e_vecs, CeedVector *q_vecs, CeedInt start_e, CeedInt num_fields, 630 CeedInt Q) { 631 Ceed ceed; 632 CeedSize e_size, q_size; 633 CeedInt max_num_points, num_comp, size, P; 634 CeedQFunctionField *qf_fields; 635 CeedOperatorField *op_fields; 636 637 { 638 Ceed ceed_parent; 639 640 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 641 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 642 if (ceed_parent) ceed = ceed_parent; 643 } 644 if (is_input) { 645 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 646 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 647 } else { 648 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 649 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 650 } 651 652 // Get max number of points 653 { 654 CeedInt dim; 655 CeedElemRestriction rstr_points = NULL; 656 CeedOperator_Ref *impl; 657 658 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL)); 659 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points)); 660 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_points, &dim)); 661 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 662 CeedCallBackend(CeedOperatorGetData(op, &impl)); 663 if (is_input) { 664 CeedCallBackend(CeedVectorCreate(ceed, dim * max_num_points, &impl->point_coords_elem)); 665 CeedCallBackend(CeedVectorSetValue(impl->point_coords_elem, 0.0)); 666 } 667 } 668 669 // Loop over fields 670 for (CeedInt i = 0; i < num_fields; i++) { 671 CeedEvalMode eval_mode; 672 CeedBasis basis; 673 674 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 675 if (eval_mode != CEED_EVAL_WEIGHT) { 676 CeedElemRestriction elem_rstr; 677 678 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr)); 679 CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs_full[i + start_e])); 680 CeedCallBackend(CeedVectorSetValue(e_vecs_full[i + start_e], 0.0)); 681 } 682 683 switch (eval_mode) { 684 case CEED_EVAL_NONE: { 685 CeedVector vec; 686 687 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 688 e_size = (CeedSize)max_num_points * size; 689 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); 690 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 691 if (vec == CEED_VECTOR_ACTIVE || !is_input) { 692 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &q_vecs[i])); 693 } else { 694 q_size = (CeedSize)max_num_points * size; 695 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 696 } 697 break; 698 } 699 case CEED_EVAL_INTERP: 700 case CEED_EVAL_GRAD: 701 case CEED_EVAL_DIV: 702 case CEED_EVAL_CURL: 703 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 704 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 705 CeedCallBackend(CeedBasisGetNumNodes(basis, &P)); 706 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 707 e_size = (CeedSize)P * num_comp; 708 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); 709 q_size = (CeedSize)max_num_points * size; 710 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 711 break; 712 case CEED_EVAL_WEIGHT: // Only on input fields 713 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 714 q_size = (CeedSize)max_num_points; 715 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 716 CeedCallBackend( 717 CeedBasisApplyAtPoints(basis, 1, &max_num_points, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, CEED_VECTOR_NONE, q_vecs[i])); 718 break; 719 } 720 // Initialize full arrays for E-vectors and Q-vectors 721 if (e_vecs[i]) CeedCallBackend(CeedVectorSetValue(e_vecs[i], 0.0)); 722 if (eval_mode != CEED_EVAL_WEIGHT) CeedCallBackend(CeedVectorSetValue(q_vecs[i], 0.0)); 723 } 724 // Drop duplicate restrictions 725 if (is_input) { 726 for (CeedInt i = 0; i < num_fields; i++) { 727 CeedVector vec_i; 728 CeedElemRestriction rstr_i; 729 730 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec_i)); 731 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr_i)); 732 for (CeedInt j = i + 1; j < num_fields; j++) { 733 CeedVector vec_j; 734 CeedElemRestriction rstr_j; 735 736 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[j], &vec_j)); 737 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[j], &rstr_j)); 738 if (vec_i == vec_j && rstr_i == rstr_j) { 739 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &e_vecs[j])); 740 CeedCallBackend(CeedVectorReferenceCopy(e_vecs_full[i + start_e], &e_vecs_full[j + start_e])); 741 skip_rstr[j] = true; 742 } 743 } 744 } 745 } else { 746 for (CeedInt i = num_fields - 1; i >= 0; i--) { 747 CeedVector vec_i; 748 CeedElemRestriction rstr_i; 749 750 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec_i)); 751 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr_i)); 752 for (CeedInt j = i - 1; j >= 0; j--) { 753 CeedVector vec_j; 754 CeedElemRestriction rstr_j; 755 756 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[j], &vec_j)); 757 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[j], &rstr_j)); 758 if (vec_i == vec_j && rstr_i == rstr_j) { 759 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &e_vecs[j])); 760 CeedCallBackend(CeedVectorReferenceCopy(e_vecs_full[i + start_e], &e_vecs_full[j + start_e])); 761 skip_rstr[j] = true; 762 apply_add_basis[i] = true; 763 } 764 } 765 } 766 } 767 return CEED_ERROR_SUCCESS; 768 } 769 770 //------------------------------------------------------------------------------ 771 // Setup Operator 772 //------------------------------------------------------------------------------ 773 static int CeedOperatorSetupAtPoints_Ref(CeedOperator op) { 774 bool is_setup_done; 775 CeedInt Q, num_input_fields, num_output_fields; 776 CeedQFunctionField *qf_input_fields, *qf_output_fields; 777 CeedQFunction qf; 778 CeedOperatorField *op_input_fields, *op_output_fields; 779 CeedOperator_Ref *impl; 780 781 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 782 if (is_setup_done) return CEED_ERROR_SUCCESS; 783 784 CeedCallBackend(CeedOperatorGetData(op, &impl)); 785 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 786 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 787 CeedCallBackend(CeedQFunctionIsIdentity(qf, &impl->is_identity_qf)); 788 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 789 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 790 791 // Allocate 792 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs_full)); 793 794 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->skip_rstr_in)); 795 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->skip_rstr_out)); 796 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->apply_add_basis_out)); 797 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states)); 798 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_in)); 799 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_out)); 800 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in)); 801 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out)); 802 803 impl->num_inputs = num_input_fields; 804 impl->num_outputs = num_output_fields; 805 806 // Set up infield and outfield pointer arrays 807 // Infields 808 CeedCallBackend(CeedOperatorSetupFieldsAtPoints_Ref(qf, op, true, impl->skip_rstr_in, NULL, impl->e_vecs_full, impl->e_vecs_in, impl->q_vecs_in, 0, 809 num_input_fields, Q)); 810 // Outfields 811 CeedCallBackend(CeedOperatorSetupFieldsAtPoints_Ref(qf, op, false, impl->skip_rstr_out, impl->apply_add_basis_out, impl->e_vecs_full, 812 impl->e_vecs_out, impl->q_vecs_out, num_input_fields, num_output_fields, Q)); 813 814 // Identity QFunctions 815 if (impl->is_identity_qf) { 816 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->q_vecs_out[0])); 817 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->e_vecs_out[0])); 818 } 819 820 CeedCallBackend(CeedOperatorSetSetupDone(op)); 821 return CEED_ERROR_SUCCESS; 822 } 823 824 //------------------------------------------------------------------------------ 825 // Input Basis Action 826 //------------------------------------------------------------------------------ 827 static inline int CeedOperatorInputBasisAtPoints_Ref(CeedInt e, CeedInt num_points_offset, CeedInt num_points, CeedQFunctionField *qf_input_fields, 828 CeedOperatorField *op_input_fields, CeedInt num_input_fields, CeedVector in_vec, 829 CeedVector point_coords_elem, bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX], 830 CeedOperator_Ref *impl, CeedRequest *request) { 831 for (CeedInt i = 0; i < num_input_fields; i++) { 832 bool is_active_input = false; 833 CeedInt elem_size, size, num_comp; 834 CeedRestrictionType rstr_type; 835 CeedEvalMode eval_mode; 836 CeedVector vec; 837 CeedElemRestriction elem_rstr; 838 CeedBasis basis; 839 840 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 841 // Skip active input 842 is_active_input = vec == CEED_VECTOR_ACTIVE; 843 if (skip_active && is_active_input) continue; 844 845 // Get elem_size, eval_mode, size 846 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 847 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 848 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 849 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 850 // Restrict block active input 851 if (is_active_input && !impl->skip_rstr_in[i]) { 852 if (rstr_type == CEED_RESTRICTION_POINTS) { 853 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[i], request)); 854 } else { 855 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[i], request)); 856 } 857 } 858 // Basis action 859 switch (eval_mode) { 860 case CEED_EVAL_NONE: 861 if (!is_active_input) { 862 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][num_points_offset * size])); 863 } 864 break; 865 // Note - these basis eval modes require FEM fields 866 case CEED_EVAL_INTERP: 867 case CEED_EVAL_GRAD: 868 case CEED_EVAL_DIV: 869 case CEED_EVAL_CURL: 870 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 871 if (!is_active_input) { 872 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 873 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 874 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][(CeedSize)e * elem_size * num_comp])); 875 } 876 CeedCallBackend( 877 CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_NOTRANSPOSE, eval_mode, point_coords_elem, impl->e_vecs_in[i], impl->q_vecs_in[i])); 878 break; 879 case CEED_EVAL_WEIGHT: 880 break; // No action 881 } 882 } 883 return CEED_ERROR_SUCCESS; 884 } 885 886 //------------------------------------------------------------------------------ 887 // Output Basis Action 888 //------------------------------------------------------------------------------ 889 static inline int CeedOperatorOutputBasisAtPoints_Ref(CeedInt e, CeedInt num_points_offset, CeedInt num_points, CeedQFunctionField *qf_output_fields, 890 CeedOperatorField *op_output_fields, CeedInt num_input_fields, CeedInt num_output_fields, 891 bool *apply_add_basis, bool *skip_rstr, CeedOperator op, CeedVector out_vec, 892 CeedVector point_coords_elem, CeedOperator_Ref *impl, CeedRequest *request) { 893 for (CeedInt i = 0; i < num_output_fields; i++) { 894 CeedRestrictionType rstr_type; 895 CeedEvalMode eval_mode; 896 CeedVector vec; 897 CeedElemRestriction elem_rstr; 898 CeedBasis basis; 899 900 // Get elem_size, eval_mode, size 901 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 902 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 903 // Basis action 904 switch (eval_mode) { 905 case CEED_EVAL_NONE: 906 break; // No action 907 case CEED_EVAL_INTERP: 908 case CEED_EVAL_GRAD: 909 case CEED_EVAL_DIV: 910 case CEED_EVAL_CURL: 911 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 912 if (apply_add_basis[i]) { 913 CeedCallBackend(CeedBasisApplyAddAtPoints(basis, 1, &num_points, CEED_TRANSPOSE, eval_mode, point_coords_elem, impl->q_vecs_out[i], 914 impl->e_vecs_out[i])); 915 } else { 916 CeedCallBackend( 917 CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_TRANSPOSE, eval_mode, point_coords_elem, impl->q_vecs_out[i], impl->e_vecs_out[i])); 918 } 919 break; 920 // LCOV_EXCL_START 921 case CEED_EVAL_WEIGHT: { 922 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 923 // LCOV_EXCL_STOP 924 } 925 } 926 // Restrict output block 927 if (skip_rstr[i]) continue; 928 // Get output vector 929 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 930 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 931 if (vec == CEED_VECTOR_ACTIVE) vec = out_vec; 932 // Restrict 933 if (rstr_type == CEED_RESTRICTION_POINTS) { 934 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], vec, request)); 935 } else { 936 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], vec, request)); 937 } 938 } 939 return CEED_ERROR_SUCCESS; 940 } 941 942 //------------------------------------------------------------------------------ 943 // Operator Apply 944 //------------------------------------------------------------------------------ 945 static int CeedOperatorApplyAddAtPoints_Ref(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) { 946 CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem; 947 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0}; 948 CeedVector point_coords = NULL; 949 CeedElemRestriction rstr_points = NULL; 950 CeedQFunctionField *qf_input_fields, *qf_output_fields; 951 CeedQFunction qf; 952 CeedOperatorField *op_input_fields, *op_output_fields; 953 CeedOperator_Ref *impl; 954 955 CeedCallBackend(CeedOperatorGetData(op, &impl)); 956 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 957 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 958 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 959 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 960 961 // Setup 962 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op)); 963 964 // Point coordinates 965 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 966 967 // Input Evecs and Restriction 968 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request)); 969 970 // Loop through elements 971 for (CeedInt e = 0; e < num_elem; e++) { 972 CeedInt num_points; 973 974 // Setup points for element 975 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 976 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points)); 977 978 // Input basis apply 979 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, in_vec, 980 impl->point_coords_elem, false, e_data, impl, request)); 981 982 // Q function 983 if (!impl->is_identity_qf) { 984 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out)); 985 } 986 987 // Output basis apply and restriction 988 CeedCallBackend(CeedOperatorOutputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_output_fields, op_output_fields, num_input_fields, 989 num_output_fields, impl->apply_add_basis_out, impl->skip_rstr_out, op, out_vec, 990 impl->point_coords_elem, impl, request)); 991 992 num_points_offset += num_points; 993 } 994 995 // Restore input arrays 996 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl)); 997 998 // Cleanup point coordinates 999 CeedCallBackend(CeedVectorDestroy(&point_coords)); 1000 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 1001 return CEED_ERROR_SUCCESS; 1002 } 1003 1004 //------------------------------------------------------------------------------ 1005 // Core code for assembling linear QFunction 1006 //------------------------------------------------------------------------------ 1007 static inline int CeedOperatorLinearAssembleQFunctionAtPointsCore_Ref(CeedOperator op, bool build_objects, CeedVector *assembled, 1008 CeedElemRestriction *rstr, CeedRequest *request) { 1009 Ceed ceed; 1010 CeedInt qf_size_in, qf_size_out, max_num_points, num_elem, num_input_fields, num_output_fields, num_points_offset = 0; 1011 CeedScalar *assembled_array, *e_data_full[2 * CEED_FIELD_MAX] = {NULL}; 1012 CeedVector point_coords = NULL; 1013 CeedQFunctionField *qf_input_fields, *qf_output_fields; 1014 CeedQFunction qf; 1015 CeedOperatorField *op_input_fields, *op_output_fields; 1016 CeedOperator_Ref *impl; 1017 CeedElemRestriction rstr_points = NULL; 1018 1019 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1020 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1021 qf_size_in = impl->qf_size_in; 1022 qf_size_out = impl->qf_size_out; 1023 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 1024 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 1025 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 1026 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 1027 1028 // Setup 1029 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op)); 1030 1031 // Check for restriction only operator 1032 CeedCheck(!impl->is_identity_rstr_op, ceed, CEED_ERROR_BACKEND, "Assembling restriction only operators is not supported"); 1033 1034 // Point coordinates 1035 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 1036 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points)); 1037 1038 // Input Evecs and Restriction 1039 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data_full, impl, request)); 1040 1041 // Count number of active input fields 1042 if (qf_size_in == 0) { 1043 for (CeedInt i = 0; i < num_input_fields; i++) { 1044 CeedInt field_size; 1045 CeedVector vec; 1046 1047 // Get input vector 1048 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1049 // Check if active input 1050 if (vec == CEED_VECTOR_ACTIVE) { 1051 // Check that all active inputs are nodal fields 1052 { 1053 CeedElemRestriction elem_rstr; 1054 bool is_at_points = false; 1055 1056 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 1057 CeedCallBackend(CeedElemRestrictionIsAtPoints(elem_rstr, &is_at_points)); 1058 CeedCheck(!is_at_points, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction with active input at points"); 1059 } 1060 // Get size of active input 1061 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size)); 1062 qf_size_in += field_size; 1063 } 1064 } 1065 CeedCheck(qf_size_in, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 1066 impl->qf_size_in = qf_size_in; 1067 } 1068 1069 // Count number of active output fields 1070 if (qf_size_out == 0) { 1071 for (CeedInt i = 0; i < num_output_fields; i++) { 1072 CeedInt field_size; 1073 CeedVector vec; 1074 1075 // Get output vector 1076 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 1077 // Check if active output 1078 if (vec == CEED_VECTOR_ACTIVE) { 1079 // Check that all active inputs are nodal fields 1080 { 1081 CeedElemRestriction elem_rstr; 1082 bool is_at_points = false; 1083 1084 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 1085 CeedCallBackend(CeedElemRestrictionIsAtPoints(elem_rstr, &is_at_points)); 1086 CeedCheck(!is_at_points, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction with active input at points"); 1087 } 1088 // Get size of active output 1089 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &field_size)); 1090 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 1091 qf_size_out += field_size; 1092 } 1093 } 1094 CeedCheck(qf_size_out > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 1095 impl->qf_size_out = qf_size_out; 1096 } 1097 1098 // Build objects if needed 1099 if (build_objects) { 1100 CeedInt num_points_total; 1101 const CeedInt *offsets; 1102 1103 CeedCallBackend(CeedElemRestrictionGetNumPoints(rstr_points, &num_points_total)); 1104 1105 // Create output restriction (at points) 1106 CeedCallBackend(CeedElemRestrictionGetOffsets(rstr_points, CEED_MEM_HOST, &offsets)); 1107 CeedCallBackend(CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points_total, qf_size_in * qf_size_out, 1108 qf_size_in * qf_size_out * num_points_total, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, rstr)); 1109 CeedCallBackend(CeedElemRestrictionRestoreOffsets(rstr_points, &offsets)); 1110 1111 // Create assembled vector 1112 CeedCallBackend(CeedElemRestrictionCreateVector(*rstr, assembled, NULL)); 1113 } 1114 // Clear output vector 1115 CeedCallBackend(CeedVectorSetValue(*assembled, 0.0)); 1116 CeedCallBackend(CeedVectorGetArray(*assembled, CEED_MEM_HOST, &assembled_array)); 1117 1118 // Loop through elements 1119 for (CeedInt e = 0; e < num_elem; e++) { 1120 CeedInt num_points; 1121 1122 // Setup points for element 1123 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 1124 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points)); 1125 1126 // Input basis apply 1127 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, NULL, 1128 impl->point_coords_elem, true, e_data_full, impl, request)); 1129 1130 // Assemble QFunction 1131 for (CeedInt i = 0; i < num_input_fields; i++) { 1132 CeedInt field_size; 1133 CeedVector vec; 1134 1135 // Get input vector 1136 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1137 // Check if active input 1138 if (vec != CEED_VECTOR_ACTIVE) continue; 1139 // Get size of active input 1140 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size)); 1141 for (CeedInt field = 0; field < field_size; field++) { 1142 // Set current portion of input to 1.0 1143 { 1144 CeedScalar *array; 1145 1146 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array)); 1147 for (CeedInt j = 0; j < num_points; j++) array[field * num_points + j] = 1.0; 1148 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array)); 1149 } 1150 1151 if (!impl->is_identity_qf) { 1152 // Set Outputs 1153 for (CeedInt out = 0; out < num_output_fields; out++) { 1154 CeedVector vec; 1155 CeedInt field_size; 1156 1157 // Get output vector 1158 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 1159 // Check if active output 1160 if (vec == CEED_VECTOR_ACTIVE) { 1161 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, CEED_USE_POINTER, assembled_array)); 1162 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &field_size)); 1163 assembled_array += field_size * num_points; // Advance the pointer by the size of the output 1164 } 1165 } 1166 // Apply QFunction 1167 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out)); 1168 } else { 1169 const CeedScalar *array; 1170 CeedInt field_size; 1171 1172 // Copy Identity Outputs 1173 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[0], &field_size)); 1174 CeedCallBackend(CeedVectorGetArrayRead(impl->q_vecs_out[0], CEED_MEM_HOST, &array)); 1175 for (CeedInt j = 0; j < field_size * num_points; j++) assembled_array[j] = array[j]; 1176 CeedCallBackend(CeedVectorRestoreArrayRead(impl->q_vecs_out[0], &array)); 1177 assembled_array += field_size * num_points; 1178 } 1179 // Reset input to 0.0 1180 { 1181 CeedScalar *array; 1182 1183 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array)); 1184 for (CeedInt j = 0; j < num_points; j++) array[field * num_points + j] = 0.0; 1185 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array)); 1186 } 1187 } 1188 } 1189 num_points_offset += num_points; 1190 } 1191 1192 // Un-set output Qvecs to prevent accidental overwrite of Assembled 1193 if (!impl->is_identity_qf) { 1194 for (CeedInt out = 0; out < num_output_fields; out++) { 1195 CeedVector vec; 1196 1197 // Get output vector 1198 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 1199 // Check if active output 1200 if (vec == CEED_VECTOR_ACTIVE && num_elem > 0) { 1201 CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_HOST, NULL)); 1202 } 1203 } 1204 } 1205 1206 // Restore input arrays 1207 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data_full, impl)); 1208 1209 // Restore output 1210 CeedCallBackend(CeedVectorRestoreArray(*assembled, &assembled_array)); 1211 1212 // Cleanup 1213 CeedCallBackend(CeedVectorDestroy(&point_coords)); 1214 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 1215 return CEED_ERROR_SUCCESS; 1216 } 1217 1218 //------------------------------------------------------------------------------ 1219 // Assemble Linear QFunction 1220 //------------------------------------------------------------------------------ 1221 static int CeedOperatorLinearAssembleQFunctionAtPoints_Ref(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 1222 return CeedOperatorLinearAssembleQFunctionAtPointsCore_Ref(op, true, assembled, rstr, request); 1223 } 1224 1225 //------------------------------------------------------------------------------ 1226 // Update Assembled Linear QFunction 1227 //------------------------------------------------------------------------------ 1228 static int CeedOperatorLinearAssembleQFunctionAtPointsUpdate_Ref(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr, 1229 CeedRequest *request) { 1230 return CeedOperatorLinearAssembleQFunctionAtPointsCore_Ref(op, false, &assembled, &rstr, request); 1231 } 1232 1233 //------------------------------------------------------------------------------ 1234 // Assemble Operator Diagonal AtPoints 1235 //------------------------------------------------------------------------------ 1236 static int CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref(CeedOperator op, CeedVector assembled, CeedRequest *request) { 1237 CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem, num_comp_active = 1; 1238 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0}; 1239 Ceed ceed; 1240 CeedVector point_coords = NULL, in_vec, out_vec; 1241 CeedElemRestriction rstr_points = NULL; 1242 CeedQFunctionField *qf_input_fields, *qf_output_fields; 1243 CeedQFunction qf; 1244 CeedOperatorField *op_input_fields, *op_output_fields; 1245 CeedOperator_Ref *impl; 1246 1247 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1248 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 1249 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 1250 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 1251 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 1252 1253 // Setup 1254 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op)); 1255 1256 // Ceed 1257 { 1258 Ceed ceed_parent; 1259 1260 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1261 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 1262 if (ceed_parent) ceed = ceed_parent; 1263 } 1264 1265 // Point coordinates 1266 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 1267 1268 // Input and output vectors 1269 { 1270 CeedSize input_size, output_size; 1271 1272 CeedCallBackend(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 1273 CeedCallBackend(CeedVectorCreate(ceed, input_size, &in_vec)); 1274 CeedCallBackend(CeedVectorCreate(ceed, output_size, &out_vec)); 1275 CeedCallBackend(CeedVectorSetValue(out_vec, 0.0)); 1276 } 1277 1278 // Clear input Qvecs 1279 for (CeedInt i = 0; i < num_input_fields; i++) { 1280 CeedVector vec; 1281 1282 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1283 if (vec != CEED_VECTOR_ACTIVE) continue; 1284 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 1285 } 1286 1287 // Input Evecs and Restriction 1288 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request)); 1289 1290 // Loop through elements 1291 for (CeedInt e = 0; e < num_elem; e++) { 1292 CeedInt num_points, e_vec_size = 0; 1293 1294 // Setup points for element 1295 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 1296 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points)); 1297 1298 // Input basis apply for non-active bases 1299 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, in_vec, 1300 impl->point_coords_elem, true, e_data, impl, request)); 1301 1302 // Loop over points on element 1303 for (CeedInt i = 0; i < num_input_fields; i++) { 1304 bool is_active_at_points = true; 1305 CeedInt elem_size_active = 1; 1306 CeedRestrictionType rstr_type; 1307 CeedVector vec; 1308 CeedElemRestriction elem_rstr; 1309 1310 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1311 // -- Skip non-active input 1312 if (vec != CEED_VECTOR_ACTIVE) continue; 1313 1314 // -- Get active restriction type 1315 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 1316 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 1317 is_active_at_points = rstr_type == CEED_RESTRICTION_POINTS; 1318 if (!is_active_at_points) CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size_active)); 1319 else elem_size_active = num_points; 1320 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp_active)); 1321 1322 e_vec_size = elem_size_active * num_comp_active; 1323 for (CeedInt s = 0; s < e_vec_size; s++) { 1324 CeedEvalMode eval_mode; 1325 CeedBasis basis; 1326 1327 // -- Update unit vector 1328 { 1329 CeedScalar *array; 1330 1331 if (s == 0) CeedCallBackend(CeedVectorSetValue(impl->e_vecs_in[i], 0.0)); 1332 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_in[i], CEED_MEM_HOST, &array)); 1333 array[s] = 1.0; 1334 if (s > 0) array[s - 1] = 0.0; 1335 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_in[i], &array)); 1336 } 1337 // -- Basis action 1338 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 1339 switch (eval_mode) { 1340 case CEED_EVAL_NONE: 1341 break; 1342 // Note - these basis eval modes require FEM fields 1343 case CEED_EVAL_INTERP: 1344 case CEED_EVAL_GRAD: 1345 case CEED_EVAL_DIV: 1346 case CEED_EVAL_CURL: 1347 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 1348 CeedCallBackend(CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_NOTRANSPOSE, eval_mode, impl->point_coords_elem, impl->e_vecs_in[i], 1349 impl->q_vecs_in[i])); 1350 break; 1351 case CEED_EVAL_WEIGHT: 1352 break; // No action 1353 } 1354 1355 // -- Q function 1356 if (!impl->is_identity_qf) { 1357 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out)); 1358 } 1359 1360 // -- Output basis apply and restriction 1361 CeedCallBackend(CeedOperatorOutputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_output_fields, op_output_fields, num_input_fields, 1362 num_output_fields, impl->apply_add_basis_out, impl->skip_rstr_out, op, out_vec, 1363 impl->point_coords_elem, impl, request)); 1364 1365 // -- Grab diagonal value 1366 for (CeedInt j = 0; j < num_output_fields; j++) { 1367 bool is_active_output = false; 1368 CeedInt elem_size = 0; 1369 CeedRestrictionType rstr_type; 1370 CeedEvalMode eval_mode; 1371 CeedVector vec; 1372 CeedElemRestriction elem_rstr; 1373 CeedBasis basis; 1374 1375 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[j], &vec)); 1376 // ---- Skip non-active output 1377 is_active_output = vec == CEED_VECTOR_ACTIVE; 1378 if (!is_active_output) continue; 1379 1380 // ---- Check if elem size matches 1381 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[j], &elem_rstr)); 1382 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 1383 if (is_active_at_points && rstr_type != CEED_RESTRICTION_POINTS) continue; 1384 if (rstr_type == CEED_RESTRICTION_POINTS) { 1385 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(elem_rstr, e, &elem_size)); 1386 } else { 1387 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 1388 } 1389 { 1390 CeedInt num_comp = 0; 1391 1392 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 1393 if (e_vec_size != num_comp * elem_size) continue; 1394 } 1395 1396 // ---- Basis action 1397 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[j], &eval_mode)); 1398 switch (eval_mode) { 1399 case CEED_EVAL_NONE: 1400 break; // No action 1401 case CEED_EVAL_INTERP: 1402 case CEED_EVAL_GRAD: 1403 case CEED_EVAL_DIV: 1404 case CEED_EVAL_CURL: 1405 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[j], &basis)); 1406 CeedCallBackend(CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_TRANSPOSE, eval_mode, impl->point_coords_elem, impl->q_vecs_out[j], 1407 impl->e_vecs_out[j])); 1408 break; 1409 // LCOV_EXCL_START 1410 case CEED_EVAL_WEIGHT: { 1411 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 1412 // LCOV_EXCL_STOP 1413 } 1414 } 1415 // ---- Update output vector 1416 { 1417 CeedScalar *array, current_value = 0.0; 1418 1419 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_out[j], CEED_MEM_HOST, &array)); 1420 current_value = array[s]; 1421 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_out[j], &array)); 1422 CeedCallBackend(CeedVectorSetValue(impl->e_vecs_out[j], 0.0)); 1423 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_out[j], CEED_MEM_HOST, &array)); 1424 array[s] = current_value; 1425 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_out[j], &array)); 1426 } 1427 // ---- Restrict output block 1428 if (rstr_type == CEED_RESTRICTION_POINTS) { 1429 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[j], assembled, request)); 1430 } else { 1431 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[j], assembled, request)); 1432 } 1433 } 1434 // -- Reset unit vector 1435 if (s == e_vec_size - 1) CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 1436 } 1437 } 1438 num_points_offset += num_points; 1439 } 1440 1441 // Restore input arrays 1442 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl)); 1443 1444 // Cleanup 1445 CeedCallBackend(CeedVectorDestroy(&in_vec)); 1446 CeedCallBackend(CeedVectorDestroy(&out_vec)); 1447 CeedCallBackend(CeedVectorDestroy(&point_coords)); 1448 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 1449 return CEED_ERROR_SUCCESS; 1450 } 1451 1452 //------------------------------------------------------------------------------ 1453 // Operator Destroy 1454 //------------------------------------------------------------------------------ 1455 static int CeedOperatorDestroy_Ref(CeedOperator op) { 1456 CeedOperator_Ref *impl; 1457 1458 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1459 CeedCallBackend(CeedFree(&impl->skip_rstr_in)); 1460 CeedCallBackend(CeedFree(&impl->skip_rstr_out)); 1461 CeedCallBackend(CeedFree(&impl->e_data_out_indices)); 1462 CeedCallBackend(CeedFree(&impl->apply_add_basis_out)); 1463 for (CeedInt i = 0; i < impl->num_inputs + impl->num_outputs; i++) { 1464 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_full[i])); 1465 } 1466 CeedCallBackend(CeedFree(&impl->e_vecs_full)); 1467 CeedCallBackend(CeedFree(&impl->input_states)); 1468 1469 for (CeedInt i = 0; i < impl->num_inputs; i++) { 1470 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_in[i])); 1471 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_in[i])); 1472 } 1473 CeedCallBackend(CeedFree(&impl->e_vecs_in)); 1474 CeedCallBackend(CeedFree(&impl->q_vecs_in)); 1475 1476 for (CeedInt i = 0; i < impl->num_outputs; i++) { 1477 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_out[i])); 1478 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_out[i])); 1479 } 1480 CeedCallBackend(CeedFree(&impl->e_vecs_out)); 1481 CeedCallBackend(CeedFree(&impl->q_vecs_out)); 1482 CeedCallBackend(CeedVectorDestroy(&impl->point_coords_elem)); 1483 1484 CeedCallBackend(CeedFree(&impl)); 1485 return CEED_ERROR_SUCCESS; 1486 } 1487 1488 //------------------------------------------------------------------------------ 1489 // Operator Create 1490 //------------------------------------------------------------------------------ 1491 int CeedOperatorCreate_Ref(CeedOperator op) { 1492 Ceed ceed; 1493 CeedOperator_Ref *impl; 1494 1495 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1496 CeedCallBackend(CeedCalloc(1, &impl)); 1497 CeedCallBackend(CeedOperatorSetData(op, impl)); 1498 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunction_Ref)); 1499 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionUpdate_Ref)); 1500 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAdd_Ref)); 1501 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref)); 1502 return CEED_ERROR_SUCCESS; 1503 } 1504 1505 //------------------------------------------------------------------------------ 1506 // Operator Create At Points 1507 //------------------------------------------------------------------------------ 1508 int CeedOperatorCreateAtPoints_Ref(CeedOperator op) { 1509 Ceed ceed; 1510 CeedOperator_Ref *impl; 1511 1512 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1513 CeedCallBackend(CeedCalloc(1, &impl)); 1514 CeedCallBackend(CeedOperatorSetData(op, impl)); 1515 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunctionAtPoints_Ref)); 1516 CeedCallBackend( 1517 CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionAtPointsUpdate_Ref)); 1518 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref)); 1519 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAddAtPoints_Ref)); 1520 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref)); 1521 return CEED_ERROR_SUCCESS; 1522 } 1523 1524 //------------------------------------------------------------------------------ 1525