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