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