1 // Copyright (c) 2017-2024, 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, bool *skip_rstr, 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 // Drop duplicate input restrictions 82 if (is_input) { 83 for (CeedInt i = 0; i < num_fields; i++) { 84 CeedVector vec_i; 85 CeedElemRestriction rstr_i; 86 87 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec_i)); 88 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr_i)); 89 for (CeedInt j = i + 1; j < num_fields; j++) { 90 CeedVector vec_j; 91 CeedElemRestriction rstr_j; 92 93 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[j], &vec_j)); 94 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[j], &rstr_j)); 95 if (vec_i == vec_j && rstr_i == rstr_j) { 96 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &e_vecs[j])); 97 CeedCallBackend(CeedVectorReferenceCopy(e_vecs_full[i], &e_vecs_full[j])); 98 skip_rstr[j] = true; 99 } 100 } 101 } 102 } 103 return CEED_ERROR_SUCCESS; 104 } 105 106 //------------------------------------------------------------------------------ 107 // Setup Operator 108 //------------------------------------------------------------------------------/* 109 static int CeedOperatorSetup_Ref(CeedOperator op) { 110 bool is_setup_done; 111 CeedInt Q, num_input_fields, num_output_fields; 112 CeedQFunctionField *qf_input_fields, *qf_output_fields; 113 CeedQFunction qf; 114 CeedOperatorField *op_input_fields, *op_output_fields; 115 CeedOperator_Ref *impl; 116 117 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 118 if (is_setup_done) return CEED_ERROR_SUCCESS; 119 120 CeedCallBackend(CeedOperatorGetData(op, &impl)); 121 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 122 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 123 CeedCallBackend(CeedQFunctionIsIdentity(qf, &impl->is_identity_qf)); 124 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 125 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 126 127 // Allocate 128 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs_full)); 129 130 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->skip_rstr_in)); 131 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states)); 132 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_in)); 133 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_out)); 134 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in)); 135 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out)); 136 137 impl->num_inputs = num_input_fields; 138 impl->num_outputs = num_output_fields; 139 140 // Set up infield and outfield e_vecs and q_vecs 141 // Infields 142 CeedCallBackend( 143 CeedOperatorSetupFields_Ref(qf, op, true, impl->skip_rstr_in, impl->e_vecs_full, impl->e_vecs_in, impl->q_vecs_in, 0, num_input_fields, Q)); 144 // Outfields 145 CeedCallBackend(CeedOperatorSetupFields_Ref(qf, op, false, NULL, impl->e_vecs_full, impl->e_vecs_out, impl->q_vecs_out, num_input_fields, 146 num_output_fields, Q)); 147 148 // Identity QFunctions 149 if (impl->is_identity_qf) { 150 CeedEvalMode in_mode, out_mode; 151 CeedQFunctionField *in_fields, *out_fields; 152 153 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &in_fields, NULL, &out_fields)); 154 CeedCallBackend(CeedQFunctionFieldGetEvalMode(in_fields[0], &in_mode)); 155 CeedCallBackend(CeedQFunctionFieldGetEvalMode(out_fields[0], &out_mode)); 156 157 if (in_mode == CEED_EVAL_NONE && out_mode == CEED_EVAL_NONE) { 158 impl->is_identity_rstr_op = true; 159 } else { 160 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->q_vecs_out[0])); 161 } 162 } 163 164 CeedCallBackend(CeedOperatorSetSetupDone(op)); 165 return CEED_ERROR_SUCCESS; 166 } 167 168 //------------------------------------------------------------------------------ 169 // Setup Operator Inputs 170 //------------------------------------------------------------------------------ 171 static inline int CeedOperatorSetupInputs_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 172 CeedVector in_vec, const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], 173 CeedOperator_Ref *impl, CeedRequest *request) { 174 for (CeedInt i = 0; i < num_input_fields; i++) { 175 uint64_t state; 176 CeedEvalMode eval_mode; 177 CeedVector vec; 178 CeedElemRestriction elem_rstr; 179 180 // Get input vector 181 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 182 if (vec == CEED_VECTOR_ACTIVE) { 183 if (skip_active) continue; 184 else vec = in_vec; 185 } 186 187 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 188 // Restrict and Evec 189 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 190 } else { 191 // Restrict 192 CeedCallBackend(CeedVectorGetState(vec, &state)); 193 // Skip restriction if input is unchanged 194 if ((state != impl->input_states[i] || vec == in_vec) && !impl->skip_rstr_in[i]) { 195 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 196 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_NOTRANSPOSE, vec, impl->e_vecs_full[i], request)); 197 } 198 impl->input_states[i] = state; 199 // Get evec 200 CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs_full[i], CEED_MEM_HOST, (const CeedScalar **)&e_data_full[i])); 201 } 202 } 203 return CEED_ERROR_SUCCESS; 204 } 205 206 //------------------------------------------------------------------------------ 207 // Input Basis Action 208 //------------------------------------------------------------------------------ 209 static inline int CeedOperatorInputBasis_Ref(CeedInt e, CeedInt Q, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 210 CeedInt num_input_fields, const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], 211 CeedOperator_Ref *impl) { 212 for (CeedInt i = 0; i < num_input_fields; i++) { 213 CeedInt elem_size, size, num_comp; 214 CeedEvalMode eval_mode; 215 CeedElemRestriction elem_rstr; 216 CeedBasis basis; 217 218 // Skip active input 219 if (skip_active) { 220 CeedVector vec; 221 222 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 223 if (vec == CEED_VECTOR_ACTIVE) continue; 224 } 225 // Get elem_size, eval_mode, size 226 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 227 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 228 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 229 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 230 // Basis action 231 switch (eval_mode) { 232 case CEED_EVAL_NONE: 233 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][(CeedSize)e * Q * size])); 234 break; 235 case CEED_EVAL_INTERP: 236 case CEED_EVAL_GRAD: 237 case CEED_EVAL_DIV: 238 case CEED_EVAL_CURL: 239 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 240 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 241 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][(CeedSize)e * elem_size * num_comp])); 242 CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, eval_mode, impl->e_vecs_in[i], impl->q_vecs_in[i])); 243 break; 244 case CEED_EVAL_WEIGHT: 245 break; // No action 246 } 247 } 248 return CEED_ERROR_SUCCESS; 249 } 250 251 //------------------------------------------------------------------------------ 252 // Output Basis Action 253 //------------------------------------------------------------------------------ 254 static inline int CeedOperatorOutputBasis_Ref(CeedInt e, CeedInt Q, CeedQFunctionField *qf_output_fields, CeedOperatorField *op_output_fields, 255 CeedInt num_input_fields, CeedInt num_output_fields, CeedOperator op, 256 CeedScalar *e_data_full[2 * CEED_FIELD_MAX], CeedOperator_Ref *impl) { 257 for (CeedInt i = 0; i < num_output_fields; i++) { 258 CeedInt elem_size, num_comp; 259 CeedEvalMode eval_mode; 260 CeedElemRestriction elem_rstr; 261 CeedBasis basis; 262 263 // Get elem_size, eval_mode 264 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 265 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 266 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 267 // Basis action 268 switch (eval_mode) { 269 case CEED_EVAL_NONE: 270 break; // No action 271 case CEED_EVAL_INTERP: 272 case CEED_EVAL_GRAD: 273 case CEED_EVAL_DIV: 274 case CEED_EVAL_CURL: 275 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 276 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 277 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, 278 &e_data_full[i + num_input_fields][(CeedSize)e * elem_size * num_comp])); 279 CeedCallBackend(CeedBasisApply(basis, 1, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs_out[i])); 280 break; 281 // LCOV_EXCL_START 282 case CEED_EVAL_WEIGHT: { 283 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 284 // LCOV_EXCL_STOP 285 } 286 } 287 } 288 return CEED_ERROR_SUCCESS; 289 } 290 291 //------------------------------------------------------------------------------ 292 // Restore Input Vectors 293 //------------------------------------------------------------------------------ 294 static inline int CeedOperatorRestoreInputs_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 295 const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], CeedOperator_Ref *impl) { 296 for (CeedInt i = 0; i < num_input_fields; i++) { 297 CeedEvalMode eval_mode; 298 299 // Skip active inputs 300 if (skip_active) { 301 CeedVector vec; 302 303 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 304 if (vec == CEED_VECTOR_ACTIVE) continue; 305 } 306 // Restore input 307 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 308 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 309 } else { 310 CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs_full[i], (const CeedScalar **)&e_data_full[i])); 311 } 312 } 313 return CEED_ERROR_SUCCESS; 314 } 315 316 //------------------------------------------------------------------------------ 317 // Operator Apply 318 //------------------------------------------------------------------------------ 319 static int CeedOperatorApplyAdd_Ref(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) { 320 CeedInt Q, num_elem, num_input_fields, num_output_fields, size; 321 CeedEvalMode eval_mode; 322 CeedScalar *e_data_full[2 * CEED_FIELD_MAX] = {NULL}; 323 CeedQFunctionField *qf_input_fields, *qf_output_fields; 324 CeedQFunction qf; 325 CeedOperatorField *op_input_fields, *op_output_fields; 326 CeedOperator_Ref *impl; 327 328 CeedCallBackend(CeedOperatorGetData(op, &impl)); 329 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 330 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 331 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 332 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 333 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 334 335 // Setup 336 CeedCallBackend(CeedOperatorSetup_Ref(op)); 337 338 // Restriction only operator 339 if (impl->is_identity_rstr_op) { 340 CeedElemRestriction elem_rstr; 341 342 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[0], &elem_rstr)); 343 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_full[0], request)); 344 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[0], &elem_rstr)); 345 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs_full[0], out_vec, request)); 346 return CEED_ERROR_SUCCESS; 347 } 348 349 // Input Evecs and Restriction 350 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, in_vec, false, e_data_full, impl, request)); 351 352 // Output Evecs 353 for (CeedInt i = 0; i < num_output_fields; i++) { 354 CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs_full[i + impl->num_inputs], CEED_MEM_HOST, &e_data_full[i + num_input_fields])); 355 } 356 357 // Loop through elements 358 for (CeedInt e = 0; e < num_elem; e++) { 359 // Output pointers 360 for (CeedInt i = 0; i < num_output_fields; i++) { 361 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 362 if (eval_mode == CEED_EVAL_NONE) { 363 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size)); 364 CeedCallBackend( 365 CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][(CeedSize)e * Q * size])); 366 } 367 } 368 369 // Input basis apply 370 CeedCallBackend(CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, num_input_fields, false, e_data_full, impl)); 371 372 // Q function 373 if (!impl->is_identity_qf) { 374 CeedCallBackend(CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out)); 375 } 376 377 // Output basis apply 378 CeedCallBackend( 379 CeedOperatorOutputBasis_Ref(e, Q, qf_output_fields, op_output_fields, num_input_fields, num_output_fields, op, e_data_full, impl)); 380 } 381 382 // Output restriction 383 for (CeedInt i = 0; i < num_output_fields; i++) { 384 CeedVector vec; 385 CeedElemRestriction elem_rstr; 386 387 // Restore Evec 388 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_full[i + impl->num_inputs], &e_data_full[i + num_input_fields])); 389 // Get output vector 390 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 391 // Active 392 if (vec == CEED_VECTOR_ACTIVE) vec = out_vec; 393 // Restrict 394 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 395 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs_full[i + impl->num_inputs], vec, request)); 396 } 397 398 // Restore input arrays 399 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, false, e_data_full, impl)); 400 return CEED_ERROR_SUCCESS; 401 } 402 403 //------------------------------------------------------------------------------ 404 // Core code for assembling linear QFunction 405 //------------------------------------------------------------------------------ 406 static inline int CeedOperatorLinearAssembleQFunctionCore_Ref(CeedOperator op, bool build_objects, CeedVector *assembled, CeedElemRestriction *rstr, 407 CeedRequest *request) { 408 Ceed ceed, ceed_parent; 409 CeedInt qf_size_in, qf_size_out, Q, num_elem, num_input_fields, num_output_fields; 410 CeedScalar *assembled_array, *e_data_full[2 * CEED_FIELD_MAX] = {NULL}; 411 CeedQFunctionField *qf_input_fields, *qf_output_fields; 412 CeedQFunction qf; 413 CeedOperatorField *op_input_fields, *op_output_fields; 414 CeedOperator_Ref *impl; 415 416 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 417 CeedCallBackend(CeedOperatorGetFallbackParentCeed(op, &ceed_parent)); 418 CeedCallBackend(CeedOperatorGetData(op, &impl)); 419 qf_size_in = impl->qf_size_in; 420 qf_size_out = impl->qf_size_out; 421 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 422 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 423 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 424 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 425 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 426 427 // Setup 428 CeedCallBackend(CeedOperatorSetup_Ref(op)); 429 430 // Check for restriction only operator 431 CeedCheck(!impl->is_identity_rstr_op, ceed, CEED_ERROR_BACKEND, "Assembling restriction only operators is not supported"); 432 433 // Input Evecs and Restriction 434 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data_full, impl, request)); 435 436 // Count number of active input fields 437 if (qf_size_in == 0) { 438 for (CeedInt i = 0; i < num_input_fields; i++) { 439 CeedInt field_size; 440 CeedVector vec; 441 442 // Get input vector 443 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 444 // Check if active input 445 if (vec == CEED_VECTOR_ACTIVE) { 446 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size)); 447 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 448 qf_size_in += field_size; 449 } 450 } 451 CeedCheck(qf_size_in > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 452 impl->qf_size_in = qf_size_in; 453 } 454 455 // Count number of active output fields 456 if (qf_size_out == 0) { 457 for (CeedInt i = 0; i < num_output_fields; i++) { 458 CeedInt field_size; 459 CeedVector vec; 460 461 // Get output vector 462 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 463 // Check if active output 464 if (vec == CEED_VECTOR_ACTIVE) { 465 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &field_size)); 466 qf_size_out += field_size; 467 } 468 } 469 CeedCheck(qf_size_out > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 470 impl->qf_size_out = qf_size_out; 471 } 472 473 // Build objects if needed 474 if (build_objects) { 475 const CeedSize l_size = (CeedSize)num_elem * Q * qf_size_in * qf_size_out; 476 CeedInt strides[3] = {1, Q, qf_size_in * qf_size_out * Q}; /* *NOPAD* */ 477 478 // Create output restriction 479 CeedCallBackend(CeedElemRestrictionCreateStrided(ceed_parent, num_elem, Q, qf_size_in * qf_size_out, 480 (CeedSize)qf_size_in * (CeedSize)qf_size_out * (CeedSize)num_elem * (CeedSize)Q, strides, rstr)); 481 // Create assembled vector 482 CeedCallBackend(CeedVectorCreate(ceed_parent, l_size, assembled)); 483 } 484 // Clear output vector 485 CeedCallBackend(CeedVectorSetValue(*assembled, 0.0)); 486 CeedCallBackend(CeedVectorGetArray(*assembled, CEED_MEM_HOST, &assembled_array)); 487 488 // Loop through elements 489 for (CeedInt e = 0; e < num_elem; e++) { 490 // Input basis apply 491 CeedCallBackend(CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, num_input_fields, true, e_data_full, impl)); 492 493 // Assemble QFunction 494 495 for (CeedInt i = 0; i < num_input_fields; i++) { 496 CeedInt field_size; 497 CeedVector vec; 498 499 // Set Inputs 500 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 501 if (vec != CEED_VECTOR_ACTIVE) continue; 502 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size)); 503 for (CeedInt field = 0; field < field_size; field++) { 504 // Set current portion of input to 1.0 505 { 506 CeedScalar *array; 507 508 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array)); 509 for (CeedInt j = 0; j < Q; j++) array[field * Q + j] = 1.0; 510 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array)); 511 } 512 513 if (!impl->is_identity_qf) { 514 // Set Outputs 515 for (CeedInt out = 0; out < num_output_fields; out++) { 516 CeedVector vec; 517 518 // Get output vector 519 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 520 // Check if active output 521 if (vec == CEED_VECTOR_ACTIVE) { 522 CeedInt field_size; 523 524 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, CEED_USE_POINTER, assembled_array)); 525 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &field_size)); 526 assembled_array += field_size * Q; // Advance the pointer by the size of the output 527 } 528 } 529 // Apply QFunction 530 CeedCallBackend(CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out)); 531 } else { 532 CeedInt field_size; 533 const CeedScalar *array; 534 535 // Copy Identity Outputs 536 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[0], &field_size)); 537 CeedCallBackend(CeedVectorGetArrayRead(impl->q_vecs_out[0], CEED_MEM_HOST, &array)); 538 for (CeedInt j = 0; j < field_size * Q; j++) assembled_array[j] = array[j]; 539 CeedCallBackend(CeedVectorRestoreArrayRead(impl->q_vecs_out[0], &array)); 540 assembled_array += field_size * Q; 541 } 542 // Reset input to 0.0 543 { 544 CeedScalar *array; 545 546 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array)); 547 for (CeedInt j = 0; j < Q; j++) array[field * Q + j] = 0.0; 548 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array)); 549 } 550 } 551 } 552 } 553 554 // Un-set output Qvecs to prevent accidental overwrite of Assembled 555 if (!impl->is_identity_qf) { 556 for (CeedInt out = 0; out < num_output_fields; out++) { 557 CeedVector vec; 558 559 // Get output vector 560 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 561 // Check if active output 562 if (vec == CEED_VECTOR_ACTIVE && num_elem > 0) { 563 CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_HOST, NULL)); 564 } 565 } 566 } 567 568 // Restore input arrays 569 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data_full, impl)); 570 571 // Restore output 572 CeedCallBackend(CeedVectorRestoreArray(*assembled, &assembled_array)); 573 return CEED_ERROR_SUCCESS; 574 } 575 576 //------------------------------------------------------------------------------ 577 // Assemble Linear QFunction 578 //------------------------------------------------------------------------------ 579 static int CeedOperatorLinearAssembleQFunction_Ref(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 580 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, true, assembled, rstr, request); 581 } 582 583 //------------------------------------------------------------------------------ 584 // Update Assembled Linear QFunction 585 //------------------------------------------------------------------------------ 586 static int CeedOperatorLinearAssembleQFunctionUpdate_Ref(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) { 587 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, false, &assembled, &rstr, request); 588 } 589 590 //------------------------------------------------------------------------------ 591 // Setup Input/Output Fields 592 //------------------------------------------------------------------------------ 593 static int CeedOperatorSetupFieldsAtPoints_Ref(CeedQFunction qf, CeedOperator op, bool is_input, bool *skip_rstr, CeedVector *e_vecs_full, 594 CeedVector *e_vecs, CeedVector *q_vecs, CeedInt start_e, CeedInt num_fields, CeedInt Q) { 595 Ceed ceed; 596 CeedSize e_size, q_size; 597 CeedInt max_num_points, num_comp, size, P; 598 CeedQFunctionField *qf_fields; 599 CeedOperatorField *op_fields; 600 601 { 602 Ceed ceed_parent; 603 604 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 605 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 606 if (ceed_parent) ceed = ceed_parent; 607 } 608 if (is_input) { 609 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 610 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 611 } else { 612 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 613 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 614 } 615 616 // Get max number of points 617 { 618 CeedInt dim; 619 CeedElemRestriction rstr_points = NULL; 620 CeedOperator_Ref *impl; 621 622 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL)); 623 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points)); 624 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_points, &dim)); 625 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 626 CeedCallBackend(CeedOperatorGetData(op, &impl)); 627 if (is_input) { 628 CeedCallBackend(CeedVectorCreate(ceed, dim * max_num_points, &impl->point_coords_elem)); 629 CeedCallBackend(CeedVectorSetValue(impl->point_coords_elem, 0.0)); 630 } 631 } 632 633 // Loop over fields 634 for (CeedInt i = 0; i < num_fields; i++) { 635 CeedEvalMode eval_mode; 636 CeedBasis basis; 637 638 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 639 if (eval_mode != CEED_EVAL_WEIGHT) { 640 CeedElemRestriction elem_rstr; 641 642 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr)); 643 CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs_full[i + start_e])); 644 CeedCallBackend(CeedVectorSetValue(e_vecs_full[i + start_e], 0.0)); 645 } 646 647 switch (eval_mode) { 648 case CEED_EVAL_NONE: { 649 CeedVector vec; 650 651 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 652 e_size = (CeedSize)max_num_points * size; 653 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); 654 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 655 if (vec == CEED_VECTOR_ACTIVE || !is_input) { 656 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &q_vecs[i])); 657 } else { 658 q_size = (CeedSize)max_num_points * size; 659 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 660 } 661 break; 662 } 663 case CEED_EVAL_INTERP: 664 case CEED_EVAL_GRAD: 665 case CEED_EVAL_DIV: 666 case CEED_EVAL_CURL: 667 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 668 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 669 CeedCallBackend(CeedBasisGetNumNodes(basis, &P)); 670 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 671 e_size = (CeedSize)P * num_comp; 672 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); 673 q_size = (CeedSize)max_num_points * size; 674 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 675 break; 676 case CEED_EVAL_WEIGHT: // Only on input fields 677 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 678 q_size = (CeedSize)max_num_points; 679 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 680 CeedCallBackend( 681 CeedBasisApplyAtPoints(basis, 1, &max_num_points, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, CEED_VECTOR_NONE, q_vecs[i])); 682 break; 683 } 684 // Initialize full arrays for E-vectors and Q-vectors 685 if (e_vecs[i]) CeedCallBackend(CeedVectorSetValue(e_vecs[i], 0.0)); 686 if (eval_mode != CEED_EVAL_WEIGHT) CeedCallBackend(CeedVectorSetValue(q_vecs[i], 0.0)); 687 } 688 // Drop duplicate input restrictions 689 if (is_input) { 690 for (CeedInt i = 0; i < num_fields; i++) { 691 CeedVector vec_i; 692 CeedElemRestriction rstr_i; 693 694 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec_i)); 695 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr_i)); 696 for (CeedInt j = i + 1; j < num_fields; j++) { 697 CeedVector vec_j; 698 CeedElemRestriction rstr_j; 699 700 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[j], &vec_j)); 701 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[j], &rstr_j)); 702 if (vec_i == vec_j && rstr_i == rstr_j) { 703 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &e_vecs[j])); 704 skip_rstr[j] = true; 705 } 706 } 707 } 708 } 709 return CEED_ERROR_SUCCESS; 710 } 711 712 //------------------------------------------------------------------------------ 713 // Setup Operator 714 //------------------------------------------------------------------------------ 715 static int CeedOperatorSetupAtPoints_Ref(CeedOperator op) { 716 bool is_setup_done; 717 CeedInt Q, num_input_fields, num_output_fields; 718 CeedQFunctionField *qf_input_fields, *qf_output_fields; 719 CeedQFunction qf; 720 CeedOperatorField *op_input_fields, *op_output_fields; 721 CeedOperator_Ref *impl; 722 723 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 724 if (is_setup_done) return CEED_ERROR_SUCCESS; 725 726 CeedCallBackend(CeedOperatorGetData(op, &impl)); 727 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 728 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 729 CeedCallBackend(CeedQFunctionIsIdentity(qf, &impl->is_identity_qf)); 730 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 731 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 732 733 // Allocate 734 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs_full)); 735 736 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->skip_rstr_in)); 737 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states)); 738 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_in)); 739 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_out)); 740 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in)); 741 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out)); 742 743 impl->num_inputs = num_input_fields; 744 impl->num_outputs = num_output_fields; 745 746 // Set up infield and outfield pointer arrays 747 // Infields 748 CeedCallBackend(CeedOperatorSetupFieldsAtPoints_Ref(qf, op, true, impl->skip_rstr_in, impl->e_vecs_full, impl->e_vecs_in, impl->q_vecs_in, 0, 749 num_input_fields, Q)); 750 // Outfields 751 CeedCallBackend(CeedOperatorSetupFieldsAtPoints_Ref(qf, op, false, NULL, impl->e_vecs_full, impl->e_vecs_out, impl->q_vecs_out, num_input_fields, 752 num_output_fields, Q)); 753 754 // Identity QFunctions 755 if (impl->is_identity_qf) { 756 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->q_vecs_out[0])); 757 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->e_vecs_out[0])); 758 } 759 760 CeedCallBackend(CeedOperatorSetSetupDone(op)); 761 return CEED_ERROR_SUCCESS; 762 } 763 764 //------------------------------------------------------------------------------ 765 // Input Basis Action 766 //------------------------------------------------------------------------------ 767 static inline int CeedOperatorInputBasisAtPoints_Ref(CeedInt e, CeedInt num_points_offset, CeedInt num_points, CeedQFunctionField *qf_input_fields, 768 CeedOperatorField *op_input_fields, CeedInt num_input_fields, CeedVector in_vec, 769 CeedVector point_coords_elem, bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX], 770 CeedOperator_Ref *impl, CeedRequest *request) { 771 for (CeedInt i = 0; i < num_input_fields; i++) { 772 bool is_active_input = false; 773 CeedInt elem_size, size, num_comp; 774 CeedRestrictionType rstr_type; 775 CeedEvalMode eval_mode; 776 CeedVector vec; 777 CeedElemRestriction elem_rstr; 778 CeedBasis basis; 779 780 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 781 // Skip active input 782 is_active_input = vec == CEED_VECTOR_ACTIVE; 783 if (skip_active && is_active_input) continue; 784 785 // Get elem_size, eval_mode, size 786 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 787 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 788 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 789 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 790 // Restrict block active input 791 if (is_active_input && !impl->skip_rstr_in[i]) { 792 if (rstr_type == CEED_RESTRICTION_POINTS) { 793 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[i], request)); 794 } else { 795 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[i], request)); 796 } 797 } 798 // Basis action 799 switch (eval_mode) { 800 case CEED_EVAL_NONE: 801 if (!is_active_input) { 802 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][num_points_offset * size])); 803 } 804 break; 805 // Note - these basis eval modes require FEM fields 806 case CEED_EVAL_INTERP: 807 case CEED_EVAL_GRAD: 808 case CEED_EVAL_DIV: 809 case CEED_EVAL_CURL: 810 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 811 if (!is_active_input) { 812 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 813 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 814 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][(CeedSize)e * elem_size * num_comp])); 815 } 816 CeedCallBackend( 817 CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_NOTRANSPOSE, eval_mode, point_coords_elem, impl->e_vecs_in[i], impl->q_vecs_in[i])); 818 break; 819 case CEED_EVAL_WEIGHT: 820 break; // No action 821 } 822 } 823 return CEED_ERROR_SUCCESS; 824 } 825 826 //------------------------------------------------------------------------------ 827 // Output Basis Action 828 //------------------------------------------------------------------------------ 829 static inline int CeedOperatorOutputBasisAtPoints_Ref(CeedInt e, CeedInt num_points_offset, CeedInt num_points, CeedQFunctionField *qf_output_fields, 830 CeedOperatorField *op_output_fields, CeedInt num_input_fields, CeedInt num_output_fields, 831 CeedOperator op, CeedVector out_vec, CeedVector point_coords_elem, CeedOperator_Ref *impl, 832 CeedRequest *request) { 833 for (CeedInt i = 0; i < num_output_fields; i++) { 834 CeedRestrictionType rstr_type; 835 CeedEvalMode eval_mode; 836 CeedVector vec; 837 CeedElemRestriction elem_rstr; 838 CeedBasis basis; 839 840 // Get elem_size, eval_mode, size 841 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 842 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 843 // Basis action 844 switch (eval_mode) { 845 case CEED_EVAL_NONE: 846 break; // No action 847 case CEED_EVAL_INTERP: 848 case CEED_EVAL_GRAD: 849 case CEED_EVAL_DIV: 850 case CEED_EVAL_CURL: 851 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 852 CeedCallBackend( 853 CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_TRANSPOSE, eval_mode, point_coords_elem, impl->q_vecs_out[i], impl->e_vecs_out[i])); 854 break; 855 // LCOV_EXCL_START 856 case CEED_EVAL_WEIGHT: { 857 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 858 // LCOV_EXCL_STOP 859 } 860 } 861 // Restrict output block 862 // Get output vector 863 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 864 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 865 if (vec == CEED_VECTOR_ACTIVE) vec = out_vec; 866 // Restrict 867 if (rstr_type == CEED_RESTRICTION_POINTS) { 868 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], vec, request)); 869 } else { 870 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], vec, request)); 871 } 872 } 873 return CEED_ERROR_SUCCESS; 874 } 875 876 //------------------------------------------------------------------------------ 877 // Operator Apply 878 //------------------------------------------------------------------------------ 879 static int CeedOperatorApplyAddAtPoints_Ref(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) { 880 CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem; 881 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0}; 882 CeedVector point_coords = NULL; 883 CeedElemRestriction rstr_points = NULL; 884 CeedQFunctionField *qf_input_fields, *qf_output_fields; 885 CeedQFunction qf; 886 CeedOperatorField *op_input_fields, *op_output_fields; 887 CeedOperator_Ref *impl; 888 889 CeedCallBackend(CeedOperatorGetData(op, &impl)); 890 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 891 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 892 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 893 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 894 895 // Setup 896 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op)); 897 898 // Point coordinates 899 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 900 901 // Input Evecs and Restriction 902 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request)); 903 904 // Loop through elements 905 for (CeedInt e = 0; e < num_elem; e++) { 906 CeedInt num_points; 907 908 // Setup points for element 909 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 910 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points)); 911 912 // Input basis apply 913 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, in_vec, 914 impl->point_coords_elem, false, e_data, impl, request)); 915 916 // Q function 917 if (!impl->is_identity_qf) { 918 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out)); 919 } 920 921 // Output basis apply and restriction 922 CeedCallBackend(CeedOperatorOutputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_output_fields, op_output_fields, num_input_fields, 923 num_output_fields, op, out_vec, impl->point_coords_elem, impl, request)); 924 925 num_points_offset += num_points; 926 } 927 928 // Restore input arrays 929 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl)); 930 931 // Cleanup point coordinates 932 CeedCallBackend(CeedVectorDestroy(&point_coords)); 933 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 934 return CEED_ERROR_SUCCESS; 935 } 936 937 //------------------------------------------------------------------------------ 938 // Core code for assembling linear QFunction 939 //------------------------------------------------------------------------------ 940 static inline int CeedOperatorLinearAssembleQFunctionAtPointsCore_Ref(CeedOperator op, bool build_objects, CeedVector *assembled, 941 CeedElemRestriction *rstr, CeedRequest *request) { 942 Ceed ceed; 943 CeedInt qf_size_in, qf_size_out, max_num_points, num_elem, num_input_fields, num_output_fields, num_points_offset = 0; 944 CeedScalar *assembled_array, *e_data_full[2 * CEED_FIELD_MAX] = {NULL}; 945 CeedVector point_coords = NULL; 946 CeedQFunctionField *qf_input_fields, *qf_output_fields; 947 CeedQFunction qf; 948 CeedOperatorField *op_input_fields, *op_output_fields; 949 CeedOperator_Ref *impl; 950 CeedElemRestriction rstr_points = NULL; 951 952 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 953 CeedCallBackend(CeedOperatorGetData(op, &impl)); 954 qf_size_in = impl->qf_size_in; 955 qf_size_out = impl->qf_size_out; 956 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 957 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 958 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 959 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 960 961 // Setup 962 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op)); 963 964 // Check for restriction only operator 965 CeedCheck(!impl->is_identity_rstr_op, ceed, CEED_ERROR_BACKEND, "Assembling restriction only operators is not supported"); 966 967 // Point coordinates 968 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 969 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points)); 970 971 // Input Evecs and Restriction 972 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data_full, impl, request)); 973 974 // Count number of active input fields 975 if (qf_size_in == 0) { 976 for (CeedInt i = 0; i < num_input_fields; i++) { 977 CeedInt field_size; 978 CeedVector vec; 979 980 // Get input vector 981 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 982 // Check if active input 983 if (vec == CEED_VECTOR_ACTIVE) { 984 // Check that all active inputs are nodal fields 985 { 986 CeedElemRestriction elem_rstr; 987 bool is_at_points = false; 988 989 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 990 CeedCallBackend(CeedElemRestrictionIsAtPoints(elem_rstr, &is_at_points)); 991 CeedCheck(!is_at_points, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction with active input at points"); 992 } 993 // Get size of active input 994 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size)); 995 qf_size_in += field_size; 996 } 997 } 998 CeedCheck(qf_size_in, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 999 impl->qf_size_in = qf_size_in; 1000 } 1001 1002 // Count number of active output fields 1003 if (qf_size_out == 0) { 1004 for (CeedInt i = 0; i < num_output_fields; i++) { 1005 CeedInt field_size; 1006 CeedVector vec; 1007 1008 // Get output vector 1009 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 1010 // Check if active output 1011 if (vec == CEED_VECTOR_ACTIVE) { 1012 // Check that all active inputs are nodal fields 1013 { 1014 CeedElemRestriction elem_rstr; 1015 bool is_at_points = false; 1016 1017 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 1018 CeedCallBackend(CeedElemRestrictionIsAtPoints(elem_rstr, &is_at_points)); 1019 CeedCheck(!is_at_points, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction with active input at points"); 1020 } 1021 // Get size of active output 1022 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &field_size)); 1023 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 1024 qf_size_out += field_size; 1025 } 1026 } 1027 CeedCheck(qf_size_out > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 1028 impl->qf_size_out = qf_size_out; 1029 } 1030 1031 // Build objects if needed 1032 if (build_objects) { 1033 CeedInt num_points_total; 1034 const CeedInt *offsets; 1035 1036 CeedCallBackend(CeedElemRestrictionGetNumPoints(rstr_points, &num_points_total)); 1037 1038 // Create output restriction (at points) 1039 CeedCallBackend(CeedElemRestrictionGetOffsets(rstr_points, CEED_MEM_HOST, &offsets)); 1040 CeedCallBackend(CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points_total, qf_size_in * qf_size_out, 1041 qf_size_in * qf_size_out * num_points_total, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, rstr)); 1042 CeedCallBackend(CeedElemRestrictionRestoreOffsets(rstr_points, &offsets)); 1043 1044 // Create assembled vector 1045 CeedCallBackend(CeedElemRestrictionCreateVector(*rstr, assembled, NULL)); 1046 } 1047 // Clear output vector 1048 CeedCallBackend(CeedVectorSetValue(*assembled, 0.0)); 1049 CeedCallBackend(CeedVectorGetArray(*assembled, CEED_MEM_HOST, &assembled_array)); 1050 1051 // Loop through elements 1052 for (CeedInt e = 0; e < num_elem; e++) { 1053 CeedInt num_points; 1054 1055 // Setup points for element 1056 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 1057 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points)); 1058 1059 // Input basis apply 1060 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, NULL, 1061 impl->point_coords_elem, true, e_data_full, impl, request)); 1062 1063 // Assemble QFunction 1064 for (CeedInt i = 0; i < num_input_fields; i++) { 1065 CeedInt field_size; 1066 CeedVector vec; 1067 1068 // Get input vector 1069 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1070 // Check if active input 1071 if (vec != CEED_VECTOR_ACTIVE) continue; 1072 // Get size of active input 1073 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size)); 1074 for (CeedInt field = 0; field < field_size; field++) { 1075 // Set current portion of input to 1.0 1076 { 1077 CeedScalar *array; 1078 1079 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array)); 1080 for (CeedInt j = 0; j < num_points; j++) array[field * num_points + j] = 1.0; 1081 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array)); 1082 } 1083 1084 if (!impl->is_identity_qf) { 1085 // Set Outputs 1086 for (CeedInt out = 0; out < num_output_fields; out++) { 1087 CeedVector vec; 1088 CeedInt field_size; 1089 1090 // Get output vector 1091 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 1092 // Check if active output 1093 if (vec == CEED_VECTOR_ACTIVE) { 1094 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, CEED_USE_POINTER, assembled_array)); 1095 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &field_size)); 1096 assembled_array += field_size * num_points; // Advance the pointer by the size of the output 1097 } 1098 } 1099 // Apply QFunction 1100 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out)); 1101 } else { 1102 const CeedScalar *array; 1103 CeedInt field_size; 1104 1105 // Copy Identity Outputs 1106 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[0], &field_size)); 1107 CeedCallBackend(CeedVectorGetArrayRead(impl->q_vecs_out[0], CEED_MEM_HOST, &array)); 1108 for (CeedInt j = 0; j < field_size * num_points; j++) assembled_array[j] = array[j]; 1109 CeedCallBackend(CeedVectorRestoreArrayRead(impl->q_vecs_out[0], &array)); 1110 assembled_array += field_size * num_points; 1111 } 1112 // Reset input to 0.0 1113 { 1114 CeedScalar *array; 1115 1116 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array)); 1117 for (CeedInt j = 0; j < num_points; j++) array[field * num_points + j] = 0.0; 1118 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array)); 1119 } 1120 } 1121 } 1122 num_points_offset += num_points; 1123 } 1124 1125 // Un-set output Qvecs to prevent accidental overwrite of Assembled 1126 if (!impl->is_identity_qf) { 1127 for (CeedInt out = 0; out < num_output_fields; out++) { 1128 CeedVector vec; 1129 1130 // Get output vector 1131 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 1132 // Check if active output 1133 if (vec == CEED_VECTOR_ACTIVE && num_elem > 0) { 1134 CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_HOST, NULL)); 1135 } 1136 } 1137 } 1138 1139 // Restore input arrays 1140 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data_full, impl)); 1141 1142 // Restore output 1143 CeedCallBackend(CeedVectorRestoreArray(*assembled, &assembled_array)); 1144 1145 // Cleanup 1146 CeedCallBackend(CeedVectorDestroy(&point_coords)); 1147 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 1148 return CEED_ERROR_SUCCESS; 1149 } 1150 1151 //------------------------------------------------------------------------------ 1152 // Assemble Linear QFunction 1153 //------------------------------------------------------------------------------ 1154 static int CeedOperatorLinearAssembleQFunctionAtPoints_Ref(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 1155 return CeedOperatorLinearAssembleQFunctionAtPointsCore_Ref(op, true, assembled, rstr, request); 1156 } 1157 1158 //------------------------------------------------------------------------------ 1159 // Update Assembled Linear QFunction 1160 //------------------------------------------------------------------------------ 1161 static int CeedOperatorLinearAssembleQFunctionAtPointsUpdate_Ref(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr, 1162 CeedRequest *request) { 1163 return CeedOperatorLinearAssembleQFunctionAtPointsCore_Ref(op, false, &assembled, &rstr, request); 1164 } 1165 1166 //------------------------------------------------------------------------------ 1167 // Assemble Operator Diagonal AtPoints 1168 //------------------------------------------------------------------------------ 1169 static int CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref(CeedOperator op, CeedVector assembled, CeedRequest *request) { 1170 CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem, num_comp_active = 1; 1171 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0}; 1172 Ceed ceed; 1173 CeedVector point_coords = NULL, in_vec, out_vec; 1174 CeedElemRestriction rstr_points = NULL; 1175 CeedQFunctionField *qf_input_fields, *qf_output_fields; 1176 CeedQFunction qf; 1177 CeedOperatorField *op_input_fields, *op_output_fields; 1178 CeedOperator_Ref *impl; 1179 1180 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1181 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 1182 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 1183 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 1184 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 1185 1186 // Setup 1187 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op)); 1188 1189 // Ceed 1190 { 1191 Ceed ceed_parent; 1192 1193 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1194 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 1195 if (ceed_parent) ceed = ceed_parent; 1196 } 1197 1198 // Point coordinates 1199 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 1200 1201 // Input and output vectors 1202 { 1203 CeedSize input_size, output_size; 1204 1205 CeedCallBackend(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 1206 CeedCallBackend(CeedVectorCreate(ceed, input_size, &in_vec)); 1207 CeedCallBackend(CeedVectorCreate(ceed, output_size, &out_vec)); 1208 CeedCallBackend(CeedVectorSetValue(out_vec, 0.0)); 1209 } 1210 1211 // Clear input Qvecs 1212 for (CeedInt i = 0; i < num_input_fields; i++) { 1213 CeedVector vec; 1214 1215 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1216 if (vec != CEED_VECTOR_ACTIVE) continue; 1217 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 1218 } 1219 1220 // Input Evecs and Restriction 1221 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request)); 1222 1223 // Loop through elements 1224 for (CeedInt e = 0; e < num_elem; e++) { 1225 CeedInt num_points, e_vec_size = 0; 1226 1227 // Setup points for element 1228 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 1229 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points)); 1230 1231 // Input basis apply for non-active bases 1232 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, in_vec, 1233 impl->point_coords_elem, true, e_data, impl, request)); 1234 1235 // Loop over points on element 1236 for (CeedInt i = 0; i < num_input_fields; i++) { 1237 bool is_active_at_points = true; 1238 CeedInt elem_size_active = 1; 1239 CeedRestrictionType rstr_type; 1240 CeedVector vec; 1241 CeedElemRestriction elem_rstr; 1242 1243 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1244 // -- Skip non-active input 1245 if (vec != CEED_VECTOR_ACTIVE) continue; 1246 1247 // -- Get active restriction type 1248 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 1249 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 1250 is_active_at_points = rstr_type == CEED_RESTRICTION_POINTS; 1251 if (!is_active_at_points) CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size_active)); 1252 else elem_size_active = num_points; 1253 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp_active)); 1254 1255 e_vec_size = elem_size_active * num_comp_active; 1256 for (CeedInt s = 0; s < e_vec_size; s++) { 1257 CeedEvalMode eval_mode; 1258 CeedBasis basis; 1259 1260 // -- Update unit vector 1261 { 1262 CeedScalar *array; 1263 1264 if (s == 0) CeedCallBackend(CeedVectorSetValue(impl->e_vecs_in[i], 0.0)); 1265 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_in[i], CEED_MEM_HOST, &array)); 1266 array[s] = 1.0; 1267 if (s > 0) array[s - 1] = 0.0; 1268 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_in[i], &array)); 1269 } 1270 // -- Basis action 1271 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 1272 switch (eval_mode) { 1273 case CEED_EVAL_NONE: 1274 break; 1275 // Note - these basis eval modes require FEM fields 1276 case CEED_EVAL_INTERP: 1277 case CEED_EVAL_GRAD: 1278 case CEED_EVAL_DIV: 1279 case CEED_EVAL_CURL: 1280 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 1281 CeedCallBackend(CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_NOTRANSPOSE, eval_mode, impl->point_coords_elem, impl->e_vecs_in[i], 1282 impl->q_vecs_in[i])); 1283 break; 1284 case CEED_EVAL_WEIGHT: 1285 break; // No action 1286 } 1287 1288 // -- Q function 1289 if (!impl->is_identity_qf) { 1290 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out)); 1291 } 1292 1293 // -- Output basis apply and restriction 1294 CeedCallBackend(CeedOperatorOutputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_output_fields, op_output_fields, num_input_fields, 1295 num_output_fields, op, out_vec, impl->point_coords_elem, impl, request)); 1296 1297 // -- Grab diagonal value 1298 for (CeedInt j = 0; j < num_output_fields; j++) { 1299 bool is_active_output = false; 1300 CeedInt elem_size = 0; 1301 CeedRestrictionType rstr_type; 1302 CeedEvalMode eval_mode; 1303 CeedVector vec; 1304 CeedElemRestriction elem_rstr; 1305 CeedBasis basis; 1306 1307 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[j], &vec)); 1308 // ---- Skip non-active output 1309 is_active_output = vec == CEED_VECTOR_ACTIVE; 1310 if (!is_active_output) continue; 1311 1312 // ---- Check if elem size matches 1313 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[j], &elem_rstr)); 1314 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 1315 if (is_active_at_points && rstr_type != CEED_RESTRICTION_POINTS) continue; 1316 if (rstr_type == CEED_RESTRICTION_POINTS) { 1317 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(elem_rstr, e, &elem_size)); 1318 } else { 1319 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 1320 } 1321 { 1322 CeedInt num_comp = 0; 1323 1324 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 1325 if (e_vec_size != num_comp * elem_size) continue; 1326 } 1327 1328 // ---- Basis action 1329 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[j], &eval_mode)); 1330 switch (eval_mode) { 1331 case CEED_EVAL_NONE: 1332 break; // No action 1333 case CEED_EVAL_INTERP: 1334 case CEED_EVAL_GRAD: 1335 case CEED_EVAL_DIV: 1336 case CEED_EVAL_CURL: 1337 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[j], &basis)); 1338 CeedCallBackend(CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_TRANSPOSE, eval_mode, impl->point_coords_elem, impl->q_vecs_out[j], 1339 impl->e_vecs_out[j])); 1340 break; 1341 // LCOV_EXCL_START 1342 case CEED_EVAL_WEIGHT: { 1343 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 1344 // LCOV_EXCL_STOP 1345 } 1346 } 1347 // ---- Update output vector 1348 { 1349 CeedScalar *array, current_value = 0.0; 1350 1351 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_out[j], CEED_MEM_HOST, &array)); 1352 current_value = array[s]; 1353 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_out[j], &array)); 1354 CeedCallBackend(CeedVectorSetValue(impl->e_vecs_out[j], 0.0)); 1355 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_out[j], CEED_MEM_HOST, &array)); 1356 array[s] = current_value; 1357 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_out[j], &array)); 1358 } 1359 // ---- Restrict output block 1360 if (rstr_type == CEED_RESTRICTION_POINTS) { 1361 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[j], assembled, request)); 1362 } else { 1363 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[j], assembled, request)); 1364 } 1365 } 1366 // -- Reset unit vector 1367 if (s == e_vec_size - 1) CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 1368 } 1369 } 1370 num_points_offset += num_points; 1371 } 1372 1373 // Restore input arrays 1374 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl)); 1375 1376 // Cleanup 1377 CeedCallBackend(CeedVectorDestroy(&in_vec)); 1378 CeedCallBackend(CeedVectorDestroy(&out_vec)); 1379 CeedCallBackend(CeedVectorDestroy(&point_coords)); 1380 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 1381 return CEED_ERROR_SUCCESS; 1382 } 1383 1384 //------------------------------------------------------------------------------ 1385 // Operator Destroy 1386 //------------------------------------------------------------------------------ 1387 static int CeedOperatorDestroy_Ref(CeedOperator op) { 1388 CeedOperator_Ref *impl; 1389 1390 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1391 CeedCallBackend(CeedFree(&impl->skip_rstr_in)); 1392 for (CeedInt i = 0; i < impl->num_inputs + impl->num_outputs; i++) { 1393 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_full[i])); 1394 } 1395 CeedCallBackend(CeedFree(&impl->e_vecs_full)); 1396 CeedCallBackend(CeedFree(&impl->input_states)); 1397 1398 for (CeedInt i = 0; i < impl->num_inputs; i++) { 1399 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_in[i])); 1400 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_in[i])); 1401 } 1402 CeedCallBackend(CeedFree(&impl->e_vecs_in)); 1403 CeedCallBackend(CeedFree(&impl->q_vecs_in)); 1404 1405 for (CeedInt i = 0; i < impl->num_outputs; i++) { 1406 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_out[i])); 1407 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_out[i])); 1408 } 1409 CeedCallBackend(CeedFree(&impl->e_vecs_out)); 1410 CeedCallBackend(CeedFree(&impl->q_vecs_out)); 1411 CeedCallBackend(CeedVectorDestroy(&impl->point_coords_elem)); 1412 1413 CeedCallBackend(CeedFree(&impl)); 1414 return CEED_ERROR_SUCCESS; 1415 } 1416 1417 //------------------------------------------------------------------------------ 1418 // Operator Create 1419 //------------------------------------------------------------------------------ 1420 int CeedOperatorCreate_Ref(CeedOperator op) { 1421 Ceed ceed; 1422 CeedOperator_Ref *impl; 1423 1424 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1425 CeedCallBackend(CeedCalloc(1, &impl)); 1426 CeedCallBackend(CeedOperatorSetData(op, impl)); 1427 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunction_Ref)); 1428 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionUpdate_Ref)); 1429 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAdd_Ref)); 1430 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref)); 1431 return CEED_ERROR_SUCCESS; 1432 } 1433 1434 //------------------------------------------------------------------------------ 1435 // Operator Create At Points 1436 //------------------------------------------------------------------------------ 1437 int CeedOperatorCreateAtPoints_Ref(CeedOperator op) { 1438 Ceed ceed; 1439 CeedOperator_Ref *impl; 1440 1441 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1442 CeedCallBackend(CeedCalloc(1, &impl)); 1443 CeedCallBackend(CeedOperatorSetData(op, impl)); 1444 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunctionAtPoints_Ref)); 1445 CeedCallBackend( 1446 CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionAtPointsUpdate_Ref)); 1447 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref)); 1448 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAddAtPoints_Ref)); 1449 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref)); 1450 return CEED_ERROR_SUCCESS; 1451 } 1452 1453 //------------------------------------------------------------------------------ 1454