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 <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, CeedVector *e_vecs_full, CeedVector *e_vecs, 20 CeedVector *q_vecs, CeedInt start_e, CeedInt num_fields, CeedInt Q) { 21 Ceed ceed; 22 CeedSize e_size, q_size; 23 CeedInt num_comp, size, P; 24 CeedQFunctionField *qf_fields; 25 CeedOperatorField *op_fields; 26 27 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 28 if (is_input) { 29 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 30 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 31 } else { 32 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 33 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 34 } 35 36 // Loop over fields 37 for (CeedInt i = 0; i < num_fields; i++) { 38 CeedEvalMode eval_mode; 39 CeedElemRestriction elem_rstr; 40 CeedBasis basis; 41 42 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 43 if (eval_mode != CEED_EVAL_WEIGHT) { 44 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr)); 45 CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs_full[i + start_e])); 46 } 47 48 switch (eval_mode) { 49 case CEED_EVAL_NONE: 50 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 51 q_size = (CeedSize)Q * size; 52 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 53 break; 54 case CEED_EVAL_INTERP: 55 case CEED_EVAL_GRAD: 56 case CEED_EVAL_DIV: 57 case CEED_EVAL_CURL: 58 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 59 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 60 CeedCallBackend(CeedBasisGetNumNodes(basis, &P)); 61 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 62 e_size = (CeedSize)P * num_comp; 63 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); 64 q_size = (CeedSize)Q * size; 65 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 66 break; 67 case CEED_EVAL_WEIGHT: // Only on input fields 68 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 69 q_size = (CeedSize)Q; 70 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 71 CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i])); 72 break; 73 } 74 } 75 return CEED_ERROR_SUCCESS; 76 } 77 78 //------------------------------------------------------------------------------ 79 // Setup Operator 80 //------------------------------------------------------------------------------/* 81 static int CeedOperatorSetup_Ref(CeedOperator op) { 82 bool is_setup_done; 83 Ceed ceed; 84 CeedInt Q, num_input_fields, num_output_fields; 85 CeedQFunctionField *qf_input_fields, *qf_output_fields; 86 CeedQFunction qf; 87 CeedOperatorField *op_input_fields, *op_output_fields; 88 CeedOperator_Ref *impl; 89 90 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 91 if (is_setup_done) return CEED_ERROR_SUCCESS; 92 93 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 94 CeedCallBackend(CeedOperatorGetData(op, &impl)); 95 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 96 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 97 CeedCallBackend(CeedQFunctionIsIdentity(qf, &impl->is_identity_qf)); 98 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 99 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 100 101 // Allocate 102 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs_full)); 103 104 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states)); 105 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_in)); 106 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_out)); 107 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in)); 108 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out)); 109 110 impl->num_inputs = num_input_fields; 111 impl->num_outputs = num_output_fields; 112 113 // Set up infield and outfield e_vecs and q_vecs 114 // Infields 115 CeedCallBackend(CeedOperatorSetupFields_Ref(qf, op, true, impl->e_vecs_full, impl->e_vecs_in, impl->q_vecs_in, 0, num_input_fields, Q)); 116 // Outfields 117 CeedCallBackend( 118 CeedOperatorSetupFields_Ref(qf, op, false, impl->e_vecs_full, impl->e_vecs_out, impl->q_vecs_out, num_input_fields, num_output_fields, Q)); 119 120 // Identity QFunctions 121 if (impl->is_identity_qf) { 122 CeedEvalMode in_mode, out_mode; 123 CeedQFunctionField *in_fields, *out_fields; 124 125 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &in_fields, NULL, &out_fields)); 126 CeedCallBackend(CeedQFunctionFieldGetEvalMode(in_fields[0], &in_mode)); 127 CeedCallBackend(CeedQFunctionFieldGetEvalMode(out_fields[0], &out_mode)); 128 129 if (in_mode == CEED_EVAL_NONE && out_mode == CEED_EVAL_NONE) { 130 impl->is_identity_rstr_op = true; 131 } else { 132 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->q_vecs_out[0])); 133 } 134 } 135 136 CeedCallBackend(CeedOperatorSetSetupDone(op)); 137 return CEED_ERROR_SUCCESS; 138 } 139 140 //------------------------------------------------------------------------------ 141 // Setup Operator Inputs 142 //------------------------------------------------------------------------------ 143 static inline int CeedOperatorSetupInputs_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 144 CeedVector in_vec, const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], 145 CeedOperator_Ref *impl, CeedRequest *request) { 146 for (CeedInt i = 0; i < num_input_fields; i++) { 147 uint64_t state; 148 CeedEvalMode eval_mode; 149 CeedVector vec; 150 CeedElemRestriction elem_rstr; 151 152 // Get input vector 153 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 154 if (vec == CEED_VECTOR_ACTIVE) { 155 if (skip_active) continue; 156 else vec = in_vec; 157 } 158 159 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 160 // Restrict and Evec 161 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 162 } else { 163 // Restrict 164 CeedCallBackend(CeedVectorGetState(vec, &state)); 165 // Skip restriction if input is unchanged 166 if (state != impl->input_states[i] || vec == in_vec) { 167 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 168 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_NOTRANSPOSE, vec, impl->e_vecs_full[i], request)); 169 impl->input_states[i] = state; 170 } 171 // Get evec 172 CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs_full[i], CEED_MEM_HOST, (const CeedScalar **)&e_data_full[i])); 173 } 174 } 175 return CEED_ERROR_SUCCESS; 176 } 177 178 //------------------------------------------------------------------------------ 179 // Input Basis Action 180 //------------------------------------------------------------------------------ 181 static inline int CeedOperatorInputBasis_Ref(CeedInt e, CeedInt Q, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 182 CeedInt num_input_fields, const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], 183 CeedOperator_Ref *impl) { 184 for (CeedInt i = 0; i < num_input_fields; i++) { 185 CeedInt elem_size, size, num_comp; 186 CeedEvalMode eval_mode; 187 CeedElemRestriction elem_rstr; 188 CeedBasis basis; 189 190 // Skip active input 191 if (skip_active) { 192 CeedVector vec; 193 194 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 195 if (vec == CEED_VECTOR_ACTIVE) continue; 196 } 197 // Get elem_size, eval_mode, size 198 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 199 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 200 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 201 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 202 // Basis action 203 switch (eval_mode) { 204 case CEED_EVAL_NONE: 205 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * Q * size])); 206 break; 207 case CEED_EVAL_INTERP: 208 case CEED_EVAL_GRAD: 209 case CEED_EVAL_DIV: 210 case CEED_EVAL_CURL: 211 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 212 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 213 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * elem_size * num_comp])); 214 CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, eval_mode, impl->e_vecs_in[i], impl->q_vecs_in[i])); 215 break; 216 case CEED_EVAL_WEIGHT: 217 break; // No action 218 } 219 } 220 return CEED_ERROR_SUCCESS; 221 } 222 223 //------------------------------------------------------------------------------ 224 // Output Basis Action 225 //------------------------------------------------------------------------------ 226 static inline int CeedOperatorOutputBasis_Ref(CeedInt e, CeedInt Q, CeedQFunctionField *qf_output_fields, CeedOperatorField *op_output_fields, 227 CeedInt num_input_fields, CeedInt num_output_fields, CeedOperator op, 228 CeedScalar *e_data_full[2 * CEED_FIELD_MAX], CeedOperator_Ref *impl) { 229 for (CeedInt i = 0; i < num_output_fields; i++) { 230 CeedInt elem_size, num_comp; 231 CeedEvalMode eval_mode; 232 CeedElemRestriction elem_rstr; 233 CeedBasis basis; 234 235 // Get elem_size, eval_mode 236 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 237 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 238 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 239 // Basis action 240 switch (eval_mode) { 241 case CEED_EVAL_NONE: 242 break; // No action 243 case CEED_EVAL_INTERP: 244 case CEED_EVAL_GRAD: 245 case CEED_EVAL_DIV: 246 case CEED_EVAL_CURL: 247 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 248 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 249 CeedCallBackend( 250 CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * elem_size * num_comp])); 251 CeedCallBackend(CeedBasisApply(basis, 1, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs_out[i])); 252 break; 253 // LCOV_EXCL_START 254 case CEED_EVAL_WEIGHT: { 255 Ceed ceed; 256 257 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 258 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 259 // LCOV_EXCL_STOP 260 } 261 } 262 } 263 return CEED_ERROR_SUCCESS; 264 } 265 266 //------------------------------------------------------------------------------ 267 // Restore Input Vectors 268 //------------------------------------------------------------------------------ 269 static inline int CeedOperatorRestoreInputs_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 270 const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], CeedOperator_Ref *impl) { 271 for (CeedInt i = 0; i < num_input_fields; i++) { 272 CeedEvalMode eval_mode; 273 274 // Skip active inputs 275 if (skip_active) { 276 CeedVector vec; 277 278 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 279 if (vec == CEED_VECTOR_ACTIVE) continue; 280 } 281 // Restore input 282 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 283 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 284 } else { 285 CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs_full[i], (const CeedScalar **)&e_data_full[i])); 286 } 287 } 288 return CEED_ERROR_SUCCESS; 289 } 290 291 //------------------------------------------------------------------------------ 292 // Operator Apply 293 //------------------------------------------------------------------------------ 294 static int CeedOperatorApplyAdd_Ref(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) { 295 CeedInt Q, num_elem, num_input_fields, num_output_fields, size; 296 CeedEvalMode eval_mode; 297 CeedScalar *e_data_full[2 * CEED_FIELD_MAX] = {NULL}; 298 CeedQFunctionField *qf_input_fields, *qf_output_fields; 299 CeedQFunction qf; 300 CeedOperatorField *op_input_fields, *op_output_fields; 301 CeedOperator_Ref *impl; 302 303 CeedCallBackend(CeedOperatorGetData(op, &impl)); 304 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 305 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 306 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 307 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 308 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 309 310 // Setup 311 CeedCallBackend(CeedOperatorSetup_Ref(op)); 312 313 // Restriction only operator 314 if (impl->is_identity_rstr_op) { 315 CeedElemRestriction elem_rstr; 316 317 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[0], &elem_rstr)); 318 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_full[0], request)); 319 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[0], &elem_rstr)); 320 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs_full[0], out_vec, request)); 321 return CEED_ERROR_SUCCESS; 322 } 323 324 // Input Evecs and Restriction 325 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, in_vec, false, e_data_full, impl, request)); 326 327 // Output Evecs 328 for (CeedInt i = 0; i < num_output_fields; i++) { 329 CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs_full[i + impl->num_inputs], CEED_MEM_HOST, &e_data_full[i + num_input_fields])); 330 } 331 332 // Loop through elements 333 for (CeedInt e = 0; e < num_elem; e++) { 334 // Output pointers 335 for (CeedInt i = 0; i < num_output_fields; i++) { 336 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 337 if (eval_mode == CEED_EVAL_NONE) { 338 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size)); 339 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * Q * size])); 340 } 341 } 342 343 // Input basis apply 344 CeedCallBackend(CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, num_input_fields, false, e_data_full, impl)); 345 346 // Q function 347 if (!impl->is_identity_qf) { 348 CeedCallBackend(CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out)); 349 } 350 351 // Output basis apply 352 CeedCallBackend( 353 CeedOperatorOutputBasis_Ref(e, Q, qf_output_fields, op_output_fields, num_input_fields, num_output_fields, op, e_data_full, impl)); 354 } 355 356 // Output restriction 357 for (CeedInt i = 0; i < num_output_fields; i++) { 358 CeedVector vec; 359 CeedElemRestriction elem_rstr; 360 361 // Restore Evec 362 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_full[i + impl->num_inputs], &e_data_full[i + num_input_fields])); 363 // Get output vector 364 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 365 // Active 366 if (vec == CEED_VECTOR_ACTIVE) vec = out_vec; 367 // Restrict 368 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 369 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs_full[i + impl->num_inputs], vec, request)); 370 } 371 372 // Restore input arrays 373 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, false, e_data_full, impl)); 374 return CEED_ERROR_SUCCESS; 375 } 376 377 //------------------------------------------------------------------------------ 378 // Core code for assembling linear QFunction 379 //------------------------------------------------------------------------------ 380 static inline int CeedOperatorLinearAssembleQFunctionCore_Ref(CeedOperator op, bool build_objects, CeedVector *assembled, CeedElemRestriction *rstr, 381 CeedRequest *request) { 382 Ceed ceed, ceed_ceed_parent; 383 CeedSize q_size; 384 CeedInt num_active_in, num_active_out, Q, num_elem, num_input_fields, num_output_fields, size; 385 CeedScalar *assembled_array, *e_data_full[2 * CEED_FIELD_MAX] = {NULL}; 386 CeedVector *active_in; 387 CeedQFunctionField *qf_input_fields, *qf_output_fields; 388 CeedQFunction qf; 389 CeedOperatorField *op_input_fields, *op_output_fields; 390 CeedOperator_Ref *impl; 391 392 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 393 CeedCallBackend(CeedOperatorGetFallbackParentCeed(op, &ceed_ceed_parent)); 394 CeedCallBackend(CeedOperatorGetData(op, &impl)); 395 active_in = impl->qf_active_in; 396 num_active_in = impl->num_active_in, num_active_out = impl->num_active_out; 397 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 398 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 399 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 400 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 401 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 402 403 // Setup 404 CeedCallBackend(CeedOperatorSetup_Ref(op)); 405 406 // Check for restriction only operator 407 CeedCheck(!impl->is_identity_rstr_op, ceed, CEED_ERROR_BACKEND, "Assembling restriction only operators is not supported"); 408 409 // Input Evecs and Restriction 410 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data_full, impl, request)); 411 412 // Count number of active input fields 413 if (!num_active_in) { 414 for (CeedInt i = 0; i < num_input_fields; i++) { 415 CeedScalar *q_vec_array; 416 CeedVector vec; 417 418 // Get input vector 419 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 420 // Check if active input 421 if (vec == CEED_VECTOR_ACTIVE) { 422 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 423 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 424 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &q_vec_array)); 425 CeedCallBackend(CeedRealloc(num_active_in + size, &active_in)); 426 for (CeedInt field = 0; field < size; field++) { 427 q_size = (CeedSize)Q; 428 CeedCallBackend(CeedVectorCreate(ceed, q_size, &active_in[num_active_in + field])); 429 CeedCallBackend(CeedVectorSetArray(active_in[num_active_in + field], CEED_MEM_HOST, CEED_USE_POINTER, &q_vec_array[field * Q])); 430 } 431 num_active_in += size; 432 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &q_vec_array)); 433 } 434 } 435 impl->num_active_in = num_active_in; 436 impl->qf_active_in = active_in; 437 } 438 439 // Count number of active output fields 440 if (!num_active_out) { 441 for (CeedInt i = 0; i < num_output_fields; i++) { 442 CeedVector vec; 443 444 // Get output vector 445 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 446 // Check if active output 447 if (vec == CEED_VECTOR_ACTIVE) { 448 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size)); 449 num_active_out += size; 450 } 451 } 452 impl->num_active_out = num_active_out; 453 } 454 455 // Check sizes 456 CeedCheck(num_active_in > 0 && num_active_out > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 457 458 // Build objects if needed 459 if (build_objects) { 460 const CeedSize l_size = (CeedSize)num_elem * Q * num_active_in * num_active_out; 461 CeedInt strides[3] = {1, Q, num_active_in * num_active_out * Q}; /* *NOPAD* */ 462 463 // Create output restriction 464 CeedCallBackend(CeedElemRestrictionCreateStrided(ceed_ceed_parent, num_elem, Q, num_active_in * num_active_out, 465 num_active_in * num_active_out * num_elem * Q, strides, rstr)); 466 // Create assembled vector 467 CeedCallBackend(CeedVectorCreate(ceed_ceed_parent, l_size, assembled)); 468 } 469 // Clear output vector 470 CeedCallBackend(CeedVectorSetValue(*assembled, 0.0)); 471 CeedCallBackend(CeedVectorGetArray(*assembled, CEED_MEM_HOST, &assembled_array)); 472 473 // Loop through elements 474 for (CeedInt e = 0; e < num_elem; e++) { 475 // Input basis apply 476 CeedCallBackend(CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, num_input_fields, true, e_data_full, impl)); 477 478 // Assemble QFunction 479 for (CeedInt in = 0; in < num_active_in; in++) { 480 // Set Inputs 481 CeedCallBackend(CeedVectorSetValue(active_in[in], 1.0)); 482 if (num_active_in > 1) { 483 CeedCallBackend(CeedVectorSetValue(active_in[(in + num_active_in - 1) % num_active_in], 0.0)); 484 } 485 if (!impl->is_identity_qf) { 486 // Set Outputs 487 for (CeedInt out = 0; out < num_output_fields; out++) { 488 CeedVector vec; 489 490 // Get output vector 491 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 492 // Check if active output 493 if (vec == CEED_VECTOR_ACTIVE) { 494 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, CEED_USE_POINTER, assembled_array)); 495 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &size)); 496 assembled_array += size * Q; // Advance the pointer by the size of the output 497 } 498 } 499 // Apply QFunction 500 CeedCallBackend(CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out)); 501 } else { 502 const CeedScalar *q_vec_array; 503 504 // Copy Identity Outputs 505 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[0], &size)); 506 CeedCallBackend(CeedVectorGetArrayRead(impl->q_vecs_out[0], CEED_MEM_HOST, &q_vec_array)); 507 for (CeedInt i = 0; i < size * Q; i++) assembled_array[i] = q_vec_array[i]; 508 CeedCallBackend(CeedVectorRestoreArrayRead(impl->q_vecs_out[0], &q_vec_array)); 509 assembled_array += size * Q; 510 } 511 } 512 } 513 514 // Un-set output Qvecs to prevent accidental overwrite of Assembled 515 if (!impl->is_identity_qf) { 516 for (CeedInt out = 0; out < num_output_fields; out++) { 517 CeedVector vec; 518 519 // Get output vector 520 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 521 // Check if active output 522 if (vec == CEED_VECTOR_ACTIVE && num_elem > 0) { 523 CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_HOST, NULL)); 524 } 525 } 526 } 527 528 // Restore input arrays 529 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data_full, impl)); 530 531 // Restore output 532 CeedCallBackend(CeedVectorRestoreArray(*assembled, &assembled_array)); 533 return CEED_ERROR_SUCCESS; 534 } 535 536 //------------------------------------------------------------------------------ 537 // Assemble Linear QFunction 538 //------------------------------------------------------------------------------ 539 static int CeedOperatorLinearAssembleQFunction_Ref(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 540 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, true, assembled, rstr, request); 541 } 542 543 //------------------------------------------------------------------------------ 544 // Update Assembled Linear QFunction 545 //------------------------------------------------------------------------------ 546 static int CeedOperatorLinearAssembleQFunctionUpdate_Ref(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) { 547 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, false, &assembled, &rstr, request); 548 } 549 550 //------------------------------------------------------------------------------ 551 // Setup Input/Output Fields 552 //------------------------------------------------------------------------------ 553 static int CeedOperatorSetupFieldsAtPoints_Ref(CeedQFunction qf, CeedOperator op, bool is_input, CeedVector *e_vecs_full, CeedVector *e_vecs, 554 CeedVector *q_vecs, CeedInt start_e, CeedInt num_fields, CeedInt Q) { 555 Ceed ceed; 556 CeedSize e_size, q_size; 557 CeedInt max_num_points, num_comp, size, P; 558 CeedQFunctionField *qf_fields; 559 CeedOperatorField *op_fields; 560 561 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 562 if (is_input) { 563 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 564 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 565 } else { 566 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 567 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 568 } 569 570 // Get max number of points 571 { 572 CeedInt dim; 573 CeedElemRestriction rstr_points = NULL; 574 CeedOperator_Ref *impl; 575 576 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL)); 577 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points)); 578 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_points, &dim)); 579 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 580 CeedCallBackend(CeedOperatorGetData(op, &impl)); 581 if (is_input) CeedCallBackend(CeedVectorCreate(ceed, dim * max_num_points, &impl->point_coords_elem)); 582 } 583 584 // Loop over fields 585 for (CeedInt i = 0; i < num_fields; i++) { 586 CeedEvalMode eval_mode; 587 CeedBasis basis; 588 589 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 590 if (eval_mode != CEED_EVAL_WEIGHT) { 591 CeedRestrictionType rstr_type; 592 CeedElemRestriction elem_rstr; 593 594 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr)); 595 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 596 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 597 CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs_full[i + start_e])); 598 } 599 600 switch (eval_mode) { 601 case CEED_EVAL_NONE: { 602 CeedVector vec; 603 604 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 605 e_size = (CeedSize)max_num_points * size; 606 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); 607 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 608 if (vec == CEED_VECTOR_ACTIVE || !is_input) { 609 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &q_vecs[i])); 610 } else { 611 q_size = (CeedSize)max_num_points * size; 612 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 613 } 614 break; 615 } 616 case CEED_EVAL_INTERP: 617 case CEED_EVAL_GRAD: 618 case CEED_EVAL_DIV: 619 case CEED_EVAL_CURL: 620 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 621 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 622 CeedCallBackend(CeedBasisGetNumNodes(basis, &P)); 623 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 624 e_size = (CeedSize)P * num_comp; 625 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); 626 q_size = (CeedSize)max_num_points * size; 627 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 628 break; 629 case CEED_EVAL_WEIGHT: // Only on input fields 630 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 631 q_size = (CeedSize)max_num_points; 632 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 633 CeedCallBackend( 634 CeedBasisApplyAtPoints(basis, max_num_points, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, CEED_VECTOR_NONE, q_vecs[i])); 635 break; 636 } 637 if (is_input && e_vecs[i]) { 638 CeedCallBackend(CeedVectorSetArray(e_vecs[i], CEED_MEM_HOST, CEED_COPY_VALUES, NULL)); 639 } 640 } 641 return CEED_ERROR_SUCCESS; 642 } 643 644 //------------------------------------------------------------------------------ 645 // Setup Operator 646 //------------------------------------------------------------------------------ 647 static int CeedOperatorSetupAtPoints_Ref(CeedOperator op) { 648 bool is_setup_done; 649 Ceed ceed; 650 CeedInt Q, num_input_fields, num_output_fields; 651 CeedQFunctionField *qf_input_fields, *qf_output_fields; 652 CeedQFunction qf; 653 CeedOperatorField *op_input_fields, *op_output_fields; 654 CeedOperator_Ref *impl; 655 656 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 657 if (is_setup_done) return CEED_ERROR_SUCCESS; 658 659 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 660 CeedCallBackend(CeedOperatorGetData(op, &impl)); 661 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 662 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 663 CeedCallBackend(CeedQFunctionIsIdentity(qf, &impl->is_identity_qf)); 664 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 665 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 666 667 // Allocate 668 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs_full)); 669 670 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states)); 671 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_in)); 672 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_out)); 673 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in)); 674 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out)); 675 676 impl->num_inputs = num_input_fields; 677 impl->num_outputs = num_output_fields; 678 679 // Set up infield and outfield pointer arrays 680 // Infields 681 CeedCallBackend(CeedOperatorSetupFieldsAtPoints_Ref(qf, op, true, impl->e_vecs_full, impl->e_vecs_in, impl->q_vecs_in, 0, num_input_fields, Q)); 682 // Outfields 683 CeedCallBackend(CeedOperatorSetupFieldsAtPoints_Ref(qf, op, false, impl->e_vecs_full, impl->e_vecs_out, impl->q_vecs_out, num_input_fields, 684 num_output_fields, Q)); 685 686 // Identity QFunctions 687 if (impl->is_identity_qf) { 688 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->q_vecs_out[0])); 689 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->e_vecs_out[0])); 690 } 691 692 CeedCallBackend(CeedOperatorSetSetupDone(op)); 693 return CEED_ERROR_SUCCESS; 694 } 695 696 //------------------------------------------------------------------------------ 697 // Setup Input Fields 698 //------------------------------------------------------------------------------ 699 static inline int CeedOperatorSetupInputsAtPoints_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, 700 CeedOperatorField *op_input_fields, CeedVector in_vec, CeedScalar *e_data[2 * CEED_FIELD_MAX], 701 CeedOperator_Ref *impl, CeedRequest *request) { 702 for (CeedInt i = 0; i < num_input_fields; i++) { 703 uint64_t state; 704 CeedEvalMode eval_mode; 705 CeedVector vec; 706 CeedElemRestriction elem_rstr; 707 708 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 709 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 710 } else { 711 // Get input vector 712 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 713 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 714 if (vec != CEED_VECTOR_ACTIVE) { 715 // Restrict 716 CeedCallBackend(CeedVectorGetState(vec, &state)); 717 if (state != impl->input_states[i]) { 718 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_NOTRANSPOSE, vec, impl->e_vecs_full[i], request)); 719 impl->input_states[i] = state; 720 } 721 // Get evec 722 CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs_full[i], CEED_MEM_HOST, (const CeedScalar **)&e_data[i])); 723 } 724 } 725 } 726 return CEED_ERROR_SUCCESS; 727 } 728 729 //------------------------------------------------------------------------------ 730 // Input Basis Action 731 //------------------------------------------------------------------------------ 732 static inline int CeedOperatorInputBasisAtPoints_Ref(CeedInt e, CeedInt num_points_offset, CeedInt num_points, CeedQFunctionField *qf_input_fields, 733 CeedOperatorField *op_input_fields, CeedInt num_input_fields, CeedVector in_vec, 734 CeedVector point_coords_elem, bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX], 735 CeedOperator_Ref *impl, CeedRequest *request) { 736 for (CeedInt i = 0; i < num_input_fields; i++) { 737 bool is_active_input = false; 738 CeedInt elem_size, size, num_comp; 739 CeedRestrictionType rstr_type; 740 CeedEvalMode eval_mode; 741 CeedVector vec; 742 CeedElemRestriction elem_rstr; 743 CeedBasis basis; 744 745 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 746 // Skip active input 747 is_active_input = vec == CEED_VECTOR_ACTIVE; 748 if (skip_active && is_active_input) continue; 749 750 // Get elem_size, eval_mode, size 751 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 752 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 753 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 754 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 755 // Restrict block active input 756 if (is_active_input) { 757 if (rstr_type == CEED_RESTRICTION_POINTS) { 758 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[i], request)); 759 } else { 760 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[i], request)); 761 } 762 } 763 // Basis action 764 switch (eval_mode) { 765 case CEED_EVAL_NONE: 766 if (!is_active_input) { 767 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][num_points_offset * size])); 768 } 769 break; 770 // Note - these basis eval modes require FEM fields 771 case CEED_EVAL_INTERP: 772 case CEED_EVAL_GRAD: 773 case CEED_EVAL_DIV: 774 case CEED_EVAL_CURL: 775 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 776 if (!is_active_input) { 777 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 778 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 779 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][e * elem_size * num_comp])); 780 } 781 CeedCallBackend( 782 CeedBasisApplyAtPoints(basis, num_points, CEED_NOTRANSPOSE, eval_mode, point_coords_elem, impl->e_vecs_in[i], impl->q_vecs_in[i])); 783 break; 784 case CEED_EVAL_WEIGHT: 785 break; // No action 786 } 787 } 788 return CEED_ERROR_SUCCESS; 789 } 790 791 //------------------------------------------------------------------------------ 792 // Output Basis Action 793 //------------------------------------------------------------------------------ 794 static inline int CeedOperatorOutputBasisAtPoints_Ref(CeedInt e, CeedInt num_points_offset, CeedInt num_points, CeedQFunctionField *qf_output_fields, 795 CeedOperatorField *op_output_fields, CeedInt num_input_fields, CeedInt num_output_fields, 796 CeedOperator op, CeedVector out_vec, CeedVector point_coords_elem, CeedOperator_Ref *impl, 797 CeedRequest *request) { 798 for (CeedInt i = 0; i < num_output_fields; i++) { 799 CeedRestrictionType rstr_type; 800 CeedEvalMode eval_mode; 801 CeedVector vec; 802 CeedElemRestriction elem_rstr; 803 CeedBasis basis; 804 805 // Get elem_size, eval_mode, size 806 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 807 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 808 // Basis action 809 switch (eval_mode) { 810 case CEED_EVAL_NONE: 811 break; // No action 812 case CEED_EVAL_INTERP: 813 case CEED_EVAL_GRAD: 814 case CEED_EVAL_DIV: 815 case CEED_EVAL_CURL: 816 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 817 CeedCallBackend( 818 CeedBasisApplyAtPoints(basis, num_points, CEED_TRANSPOSE, eval_mode, point_coords_elem, impl->q_vecs_out[i], impl->e_vecs_out[i])); 819 break; 820 // LCOV_EXCL_START 821 case CEED_EVAL_WEIGHT: { 822 Ceed ceed; 823 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 824 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 825 // LCOV_EXCL_STOP 826 } 827 } 828 // Restrict output block 829 // Get output vector 830 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 831 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 832 if (vec == CEED_VECTOR_ACTIVE) vec = out_vec; 833 // Restrict 834 if (rstr_type == CEED_RESTRICTION_POINTS) { 835 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], vec, request)); 836 } else { 837 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], vec, request)); 838 } 839 } 840 return CEED_ERROR_SUCCESS; 841 } 842 843 //------------------------------------------------------------------------------ 844 // Restore Input Vectors 845 //------------------------------------------------------------------------------ 846 static inline int CeedOperatorRestoreInputsAtPoints_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, 847 CeedOperatorField *op_input_fields, CeedScalar *e_data[2 * CEED_FIELD_MAX], 848 CeedOperator_Ref *impl) { 849 for (CeedInt i = 0; i < num_input_fields; i++) { 850 CeedEvalMode eval_mode; 851 CeedVector vec; 852 853 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 854 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 855 if (eval_mode != CEED_EVAL_WEIGHT && vec != CEED_VECTOR_ACTIVE) { 856 CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs_full[i], (const CeedScalar **)&e_data[i])); 857 } 858 } 859 return CEED_ERROR_SUCCESS; 860 } 861 862 //------------------------------------------------------------------------------ 863 // Operator Apply 864 //------------------------------------------------------------------------------ 865 static int CeedOperatorApplyAddAtPoints_Ref(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) { 866 Ceed ceed; 867 CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem; 868 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0}; 869 CeedVector point_coords = NULL; 870 CeedElemRestriction rstr_points = NULL; 871 CeedQFunctionField *qf_input_fields, *qf_output_fields; 872 CeedQFunction qf; 873 CeedOperatorField *op_input_fields, *op_output_fields; 874 CeedOperator_Ref *impl; 875 876 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 877 CeedCallBackend(CeedOperatorGetData(op, &impl)); 878 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 879 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 880 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 881 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 882 883 // Setup 884 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op)); 885 886 // Point coordinates 887 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 888 889 // Input Evecs and Restriction 890 CeedCallBackend(CeedOperatorSetupInputsAtPoints_Ref(num_input_fields, qf_input_fields, op_input_fields, in_vec, e_data, impl, request)); 891 892 // Loop through elements 893 for (CeedInt e = 0; e < num_elem; e++) { 894 CeedInt num_points; 895 896 // Setup points for element 897 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 898 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points)); 899 900 // Input basis apply 901 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, in_vec, 902 impl->point_coords_elem, false, e_data, impl, request)); 903 904 // Q function 905 if (!impl->is_identity_qf) { 906 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out)); 907 } 908 909 // Output basis apply and restriction 910 CeedCallBackend(CeedOperatorOutputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_output_fields, op_output_fields, num_input_fields, 911 num_output_fields, op, out_vec, impl->point_coords_elem, impl, request)); 912 913 num_points_offset += num_points; 914 } 915 916 // Restore input arrays 917 CeedCallBackend(CeedOperatorRestoreInputsAtPoints_Ref(num_input_fields, qf_input_fields, op_input_fields, e_data, impl)); 918 919 // Cleanup point coordinates 920 CeedCallBackend(CeedVectorDestroy(&point_coords)); 921 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 922 return CEED_ERROR_SUCCESS; 923 } 924 925 //------------------------------------------------------------------------------ 926 // Operator Destroy 927 //------------------------------------------------------------------------------ 928 static int CeedOperatorDestroy_Ref(CeedOperator op) { 929 CeedOperator_Ref *impl; 930 931 CeedCallBackend(CeedOperatorGetData(op, &impl)); 932 for (CeedInt i = 0; i < impl->num_inputs + impl->num_outputs; i++) { 933 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_full[i])); 934 } 935 CeedCallBackend(CeedFree(&impl->e_vecs_full)); 936 CeedCallBackend(CeedFree(&impl->input_states)); 937 938 for (CeedInt i = 0; i < impl->num_inputs; i++) { 939 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_in[i])); 940 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_in[i])); 941 } 942 CeedCallBackend(CeedFree(&impl->e_vecs_in)); 943 CeedCallBackend(CeedFree(&impl->q_vecs_in)); 944 945 for (CeedInt i = 0; i < impl->num_outputs; i++) { 946 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_out[i])); 947 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_out[i])); 948 } 949 CeedCallBackend(CeedFree(&impl->e_vecs_out)); 950 CeedCallBackend(CeedFree(&impl->q_vecs_out)); 951 CeedCallBackend(CeedVectorDestroy(&impl->point_coords_elem)); 952 953 // QFunction assembly 954 for (CeedInt i = 0; i < impl->num_active_in; i++) { 955 CeedCallBackend(CeedVectorDestroy(&impl->qf_active_in[i])); 956 } 957 CeedCallBackend(CeedFree(&impl->qf_active_in)); 958 959 CeedCallBackend(CeedFree(&impl)); 960 return CEED_ERROR_SUCCESS; 961 } 962 963 //------------------------------------------------------------------------------ 964 // Operator Create 965 //------------------------------------------------------------------------------ 966 int CeedOperatorCreate_Ref(CeedOperator op) { 967 Ceed ceed; 968 CeedOperator_Ref *impl; 969 970 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 971 CeedCallBackend(CeedCalloc(1, &impl)); 972 CeedCallBackend(CeedOperatorSetData(op, impl)); 973 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunction_Ref)); 974 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionUpdate_Ref)); 975 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAdd_Ref)); 976 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref)); 977 return CEED_ERROR_SUCCESS; 978 } 979 980 //------------------------------------------------------------------------------ 981 // Operator Create At Points 982 //------------------------------------------------------------------------------ 983 int CeedOperatorCreateAtPoints_Ref(CeedOperator op) { 984 Ceed ceed; 985 CeedOperator_Ref *impl; 986 987 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 988 CeedCallBackend(CeedCalloc(1, &impl)); 989 CeedCallBackend(CeedOperatorSetData(op, impl)); 990 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAddAtPoints_Ref)); 991 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref)); 992 return CEED_ERROR_SUCCESS; 993 } 994 995 //------------------------------------------------------------------------------ 996