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