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 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 603 e_size = (CeedSize)max_num_points * size; 604 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); 605 q_size = (CeedSize)max_num_points * size; 606 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 607 break; 608 case CEED_EVAL_INTERP: 609 case CEED_EVAL_GRAD: 610 case CEED_EVAL_DIV: 611 case CEED_EVAL_CURL: 612 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 613 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 614 CeedCallBackend(CeedBasisGetNumNodes(basis, &P)); 615 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 616 e_size = (CeedSize)P * num_comp; 617 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); 618 q_size = (CeedSize)max_num_points * size; 619 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 620 break; 621 case CEED_EVAL_WEIGHT: // Only on input fields 622 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 623 q_size = (CeedSize)max_num_points; 624 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 625 CeedCallBackend( 626 CeedBasisApplyAtPoints(basis, max_num_points, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, CEED_VECTOR_NONE, q_vecs[i])); 627 break; 628 } 629 if (is_input && e_vecs[i]) { 630 CeedCallBackend(CeedVectorSetArray(e_vecs[i], CEED_MEM_HOST, CEED_COPY_VALUES, NULL)); 631 } 632 } 633 return CEED_ERROR_SUCCESS; 634 } 635 636 //------------------------------------------------------------------------------ 637 // Setup Operator 638 //------------------------------------------------------------------------------ 639 static int CeedOperatorSetupAtPoints_Ref(CeedOperator op) { 640 bool is_setup_done; 641 Ceed ceed; 642 CeedInt Q, num_input_fields, num_output_fields; 643 CeedQFunctionField *qf_input_fields, *qf_output_fields; 644 CeedQFunction qf; 645 CeedOperatorField *op_input_fields, *op_output_fields; 646 CeedOperator_Ref *impl; 647 648 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 649 if (is_setup_done) return CEED_ERROR_SUCCESS; 650 651 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 652 CeedCallBackend(CeedOperatorGetData(op, &impl)); 653 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 654 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 655 CeedCallBackend(CeedQFunctionIsIdentity(qf, &impl->is_identity_qf)); 656 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 657 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 658 659 // Allocate 660 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs_full)); 661 662 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states)); 663 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_in)); 664 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_out)); 665 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in)); 666 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out)); 667 668 impl->num_inputs = num_input_fields; 669 impl->num_outputs = num_output_fields; 670 671 // Set up infield and outfield pointer arrays 672 // Infields 673 CeedCallBackend(CeedOperatorSetupFieldsAtPoints_Ref(qf, op, true, impl->e_vecs_full, impl->e_vecs_in, impl->q_vecs_in, 0, num_input_fields, Q)); 674 // Outfields 675 CeedCallBackend(CeedOperatorSetupFieldsAtPoints_Ref(qf, op, false, impl->e_vecs_full, impl->e_vecs_out, impl->q_vecs_out, num_input_fields, 676 num_output_fields, Q)); 677 678 // Identity QFunctions 679 if (impl->is_identity_qf) { 680 CeedEvalMode in_mode, out_mode; 681 CeedQFunctionField *in_fields, *out_fields; 682 683 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &in_fields, NULL, &out_fields)); 684 CeedCallBackend(CeedQFunctionFieldGetEvalMode(in_fields[0], &in_mode)); 685 CeedCallBackend(CeedQFunctionFieldGetEvalMode(out_fields[0], &out_mode)); 686 687 if (in_mode == CEED_EVAL_NONE && out_mode == CEED_EVAL_NONE) { 688 impl->is_identity_rstr_op = true; 689 } else { 690 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->q_vecs_out[0])); 691 } 692 } 693 694 CeedCallBackend(CeedOperatorSetSetupDone(op)); 695 return CEED_ERROR_SUCCESS; 696 } 697 698 //------------------------------------------------------------------------------ 699 // Setup Input Fields 700 //------------------------------------------------------------------------------ 701 static inline int CeedOperatorSetupInputsAtPoints_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, 702 CeedOperatorField *op_input_fields, CeedVector in_vec, CeedScalar *e_data[2 * CEED_FIELD_MAX], 703 CeedOperator_Ref *impl, CeedRequest *request) { 704 for (CeedInt i = 0; i < num_input_fields; i++) { 705 uint64_t state; 706 CeedEvalMode eval_mode; 707 CeedVector vec; 708 CeedElemRestriction elem_rstr; 709 710 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 711 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 712 } else { 713 // Get input vector 714 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 715 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 716 if (vec != CEED_VECTOR_ACTIVE) { 717 // Restrict 718 CeedCallBackend(CeedVectorGetState(vec, &state)); 719 if (state != impl->input_states[i]) { 720 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_NOTRANSPOSE, vec, impl->e_vecs_full[i], request)); 721 impl->input_states[i] = state; 722 } 723 // Get evec 724 CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs_full[i], CEED_MEM_HOST, (const CeedScalar **)&e_data[i])); 725 } else { 726 // Set Qvec for CEED_EVAL_NONE 727 if (eval_mode == CEED_EVAL_NONE) { 728 CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs_in[i], CEED_MEM_HOST, (const CeedScalar **)&e_data[i])); 729 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, e_data[i])); 730 CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs_in[i], (const CeedScalar **)&e_data[i])); 731 } 732 } 733 } 734 } 735 return CEED_ERROR_SUCCESS; 736 } 737 738 //------------------------------------------------------------------------------ 739 // Input Basis Action 740 //------------------------------------------------------------------------------ 741 static inline int CeedOperatorInputBasisAtPoints_Ref(CeedInt e, CeedInt num_points_offset, CeedInt num_points, CeedQFunctionField *qf_input_fields, 742 CeedOperatorField *op_input_fields, CeedInt num_input_fields, CeedVector in_vec, 743 CeedVector point_coords_elem, bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX], 744 CeedOperator_Ref *impl, CeedRequest *request) { 745 for (CeedInt i = 0; i < num_input_fields; i++) { 746 bool is_active_input = false; 747 CeedInt elem_size, size, num_comp; 748 CeedRestrictionType rstr_type; 749 CeedEvalMode eval_mode; 750 CeedVector vec; 751 CeedElemRestriction elem_rstr; 752 CeedBasis basis; 753 754 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 755 // Skip active input 756 is_active_input = vec == CEED_VECTOR_ACTIVE; 757 if (skip_active && is_active_input) continue; 758 759 // Get elem_size, eval_mode, size 760 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 761 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 762 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 763 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 764 // Restrict block active input 765 if (is_active_input) { 766 if (rstr_type == CEED_RESTRICTION_POINTS) { 767 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[i], request)); 768 } else { 769 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[i], request)); 770 } 771 } 772 // Basis action 773 switch (eval_mode) { 774 case CEED_EVAL_NONE: 775 if (!is_active_input) { 776 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][num_points_offset * size])); 777 } 778 break; 779 // Note - these basis eval modes require FEM fields 780 case CEED_EVAL_INTERP: 781 case CEED_EVAL_GRAD: 782 case CEED_EVAL_DIV: 783 case CEED_EVAL_CURL: 784 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 785 if (!is_active_input) { 786 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 787 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 788 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][e * elem_size * num_comp])); 789 } 790 CeedCallBackend( 791 CeedBasisApplyAtPoints(basis, num_points, CEED_NOTRANSPOSE, eval_mode, point_coords_elem, impl->e_vecs_in[i], impl->q_vecs_in[i])); 792 break; 793 case CEED_EVAL_WEIGHT: 794 break; // No action 795 } 796 } 797 return CEED_ERROR_SUCCESS; 798 } 799 800 //------------------------------------------------------------------------------ 801 // Output Basis Action 802 //------------------------------------------------------------------------------ 803 static inline int CeedOperatorOutputBasisAtPoints_Ref(CeedInt e, CeedInt num_points_offset, CeedInt num_points, CeedQFunctionField *qf_output_fields, 804 CeedOperatorField *op_output_fields, CeedInt num_input_fields, CeedInt num_output_fields, 805 CeedOperator op, CeedVector out_vec, CeedVector point_coords_elem, CeedOperator_Ref *impl, 806 CeedRequest *request) { 807 for (CeedInt i = 0; i < num_output_fields; i++) { 808 CeedRestrictionType rstr_type; 809 CeedEvalMode eval_mode; 810 CeedVector vec; 811 CeedElemRestriction elem_rstr; 812 CeedBasis basis; 813 814 // Get elem_size, eval_mode, size 815 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 816 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 817 // Basis action 818 switch (eval_mode) { 819 case CEED_EVAL_NONE: 820 break; // No action 821 case CEED_EVAL_INTERP: 822 case CEED_EVAL_GRAD: 823 case CEED_EVAL_DIV: 824 case CEED_EVAL_CURL: 825 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 826 CeedCallBackend( 827 CeedBasisApplyAtPoints(basis, num_points, CEED_TRANSPOSE, eval_mode, point_coords_elem, impl->q_vecs_out[i], impl->e_vecs_out[i])); 828 break; 829 // LCOV_EXCL_START 830 case CEED_EVAL_WEIGHT: { 831 Ceed ceed; 832 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 833 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 834 // LCOV_EXCL_STOP 835 } 836 } 837 // Restrict output block 838 // Get output vector 839 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 840 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 841 if (vec == CEED_VECTOR_ACTIVE) vec = out_vec; 842 // Restrict 843 if (rstr_type == CEED_RESTRICTION_POINTS) { 844 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], vec, request)); 845 } else { 846 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], vec, request)); 847 } 848 } 849 return CEED_ERROR_SUCCESS; 850 } 851 852 //------------------------------------------------------------------------------ 853 // Restore Input Vectors 854 //------------------------------------------------------------------------------ 855 static inline int CeedOperatorRestoreInputsAtPoints_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, 856 CeedOperatorField *op_input_fields, CeedScalar *e_data[2 * CEED_FIELD_MAX], 857 CeedOperator_Ref *impl) { 858 for (CeedInt i = 0; i < num_input_fields; i++) { 859 CeedEvalMode eval_mode; 860 CeedVector vec; 861 862 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 863 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 864 if (eval_mode != CEED_EVAL_WEIGHT && vec != CEED_VECTOR_ACTIVE) { 865 CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs_full[i], (const CeedScalar **)&e_data[i])); 866 } 867 } 868 return CEED_ERROR_SUCCESS; 869 } 870 871 //------------------------------------------------------------------------------ 872 // Operator Apply 873 //------------------------------------------------------------------------------ 874 static int CeedOperatorApplyAddAtPoints_Ref(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) { 875 Ceed ceed; 876 CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem; 877 CeedEvalMode eval_mode; 878 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0}; 879 CeedVector point_coords = NULL; 880 CeedElemRestriction rstr_points = NULL; 881 CeedQFunctionField *qf_input_fields, *qf_output_fields; 882 CeedQFunction qf; 883 CeedOperatorField *op_input_fields, *op_output_fields; 884 CeedOperator_Ref *impl; 885 886 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 887 CeedCallBackend(CeedOperatorGetData(op, &impl)); 888 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 889 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 890 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 891 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 892 893 // Setup 894 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op)); 895 896 // Restriction only operator 897 if (impl->is_identity_rstr_op) { 898 // TODO: Fix this up 899 for (CeedInt e = 0; e < num_elem; e++) { 900 // CeedCallBackend(CeedElemRestrictionApplyBlock(impl->block_rstr[0], e, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[0], request)); 901 // CeedCallBackend(CeedElemRestrictionApplyBlock(impl->block_rstr[1], e, CEED_TRANSPOSE, impl->e_vecs_in[0], out_vec, request)); 902 } 903 return CEED_ERROR_SUCCESS; 904 } 905 906 // Point coordinates 907 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 908 909 // Input Evecs and Restriction 910 CeedCallBackend(CeedOperatorSetupInputsAtPoints_Ref(num_input_fields, qf_input_fields, op_input_fields, in_vec, e_data, impl, request)); 911 912 // Output Lvecs, Evecs, and Qvecs 913 for (CeedInt i = 0; i < num_output_fields; i++) { 914 // Set Qvec if needed 915 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 916 if (eval_mode == CEED_EVAL_NONE) { 917 // Set qvec to single block evec 918 CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs_out[i], CEED_MEM_HOST, &e_data[i + num_input_fields])); 919 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, e_data[i + num_input_fields])); 920 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_out[i], &e_data[i + num_input_fields])); 921 } 922 } 923 924 // Loop through elements 925 for (CeedInt e = 0; e < num_elem; e++) { 926 CeedInt num_points; 927 928 // Setup points for element 929 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 930 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points)); 931 932 // Input basis apply 933 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, in_vec, 934 impl->point_coords_elem, false, e_data, impl, request)); 935 936 // Q function 937 if (!impl->is_identity_qf) { 938 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out)); 939 } 940 941 // Output basis apply and restriction 942 CeedCallBackend(CeedOperatorOutputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_output_fields, op_output_fields, num_input_fields, 943 num_output_fields, op, out_vec, impl->point_coords_elem, impl, request)); 944 945 num_points_offset += num_points; 946 } 947 948 // Restore input arrays 949 CeedCallBackend(CeedOperatorRestoreInputsAtPoints_Ref(num_input_fields, qf_input_fields, op_input_fields, e_data, impl)); 950 951 // Cleanup point coordinates 952 CeedCallBackend(CeedVectorDestroy(&point_coords)); 953 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 954 return CEED_ERROR_SUCCESS; 955 } 956 957 //------------------------------------------------------------------------------ 958 // Operator Destroy 959 //------------------------------------------------------------------------------ 960 static int CeedOperatorDestroy_Ref(CeedOperator op) { 961 CeedOperator_Ref *impl; 962 963 CeedCallBackend(CeedOperatorGetData(op, &impl)); 964 for (CeedInt i = 0; i < impl->num_inputs + impl->num_outputs; i++) { 965 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_full[i])); 966 } 967 CeedCallBackend(CeedFree(&impl->e_vecs_full)); 968 CeedCallBackend(CeedFree(&impl->input_states)); 969 970 for (CeedInt i = 0; i < impl->num_inputs; i++) { 971 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_in[i])); 972 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_in[i])); 973 } 974 CeedCallBackend(CeedFree(&impl->e_vecs_in)); 975 CeedCallBackend(CeedFree(&impl->q_vecs_in)); 976 977 for (CeedInt i = 0; i < impl->num_outputs; i++) { 978 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_out[i])); 979 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_out[i])); 980 } 981 CeedCallBackend(CeedFree(&impl->e_vecs_out)); 982 CeedCallBackend(CeedFree(&impl->q_vecs_out)); 983 CeedCallBackend(CeedVectorDestroy(&impl->point_coords_elem)); 984 985 // QFunction assembly 986 for (CeedInt i = 0; i < impl->num_active_in; i++) { 987 CeedCallBackend(CeedVectorDestroy(&impl->qf_active_in[i])); 988 } 989 CeedCallBackend(CeedFree(&impl->qf_active_in)); 990 991 CeedCallBackend(CeedFree(&impl)); 992 return CEED_ERROR_SUCCESS; 993 } 994 995 //------------------------------------------------------------------------------ 996 // Operator Create 997 //------------------------------------------------------------------------------ 998 int CeedOperatorCreate_Ref(CeedOperator op) { 999 Ceed ceed; 1000 CeedOperator_Ref *impl; 1001 1002 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1003 CeedCallBackend(CeedCalloc(1, &impl)); 1004 CeedCallBackend(CeedOperatorSetData(op, impl)); 1005 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunction_Ref)); 1006 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionUpdate_Ref)); 1007 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAdd_Ref)); 1008 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref)); 1009 return CEED_ERROR_SUCCESS; 1010 } 1011 1012 //------------------------------------------------------------------------------ 1013 // Operator Create At Points 1014 //------------------------------------------------------------------------------ 1015 int CeedOperatorCreateAtPoints_Ref(CeedOperator op) { 1016 Ceed ceed; 1017 CeedOperator_Ref *impl; 1018 1019 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1020 CeedCallBackend(CeedCalloc(1, &impl)); 1021 CeedCallBackend(CeedOperatorSetData(op, impl)); 1022 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAddAtPoints_Ref)); 1023 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref)); 1024 return CEED_ERROR_SUCCESS; 1025 } 1026 1027 //------------------------------------------------------------------------------ 1028