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-impl.h> 9 #include <ceed.h> 10 #include <ceed/backend.h> 11 #include <stdbool.h> 12 #include <stdio.h> 13 #include <string.h> 14 15 /// @file 16 /// Implementation of CeedOperator interfaces 17 18 /// ---------------------------------------------------------------------------- 19 /// CeedOperator Library Internal Functions 20 /// ---------------------------------------------------------------------------- 21 /// @addtogroup CeedOperatorDeveloper 22 /// @{ 23 24 /** 25 @brief Check if a CeedOperator Field matches the QFunction Field 26 27 @param[in] ceed Ceed object for error handling 28 @param[in] qf_field QFunction Field matching Operator Field 29 @param[in] r Operator Field ElemRestriction 30 @param[in] b Operator Field Basis 31 32 @return An error code: 0 - success, otherwise - failure 33 34 @ref Developer 35 **/ 36 static int CeedOperatorCheckField(Ceed ceed, CeedQFunctionField qf_field, CeedElemRestriction r, CeedBasis b) { 37 CeedEvalMode eval_mode = qf_field->eval_mode; 38 CeedInt dim = 1, num_comp = 1, Q_comp = 1, restr_num_comp = 1, size = qf_field->size; 39 40 // Restriction 41 if (r != CEED_ELEMRESTRICTION_NONE) { 42 if (eval_mode == CEED_EVAL_WEIGHT) { 43 // LCOV_EXCL_START 44 return CeedError(ceed, CEED_ERROR_INCOMPATIBLE, "CEED_ELEMRESTRICTION_NONE should be used for a field with eval mode CEED_EVAL_WEIGHT"); 45 // LCOV_EXCL_STOP 46 } 47 CeedCall(CeedElemRestrictionGetNumComponents(r, &restr_num_comp)); 48 } 49 if ((r == CEED_ELEMRESTRICTION_NONE) != (eval_mode == CEED_EVAL_WEIGHT)) { 50 // LCOV_EXCL_START 51 return CeedError(ceed, CEED_ERROR_INCOMPATIBLE, "CEED_ELEMRESTRICTION_NONE and CEED_EVAL_WEIGHT must be used together."); 52 // LCOV_EXCL_STOP 53 } 54 // Basis 55 if (b != CEED_BASIS_COLLOCATED) { 56 if (eval_mode == CEED_EVAL_NONE) { 57 // LCOV_EXCL_START 58 return CeedError(ceed, CEED_ERROR_INCOMPATIBLE, "Field '%s' configured with CEED_EVAL_NONE must be used with CEED_BASIS_COLLOCATED", 59 qf_field->field_name); 60 // LCOV_EXCL_STOP 61 } 62 CeedCall(CeedBasisGetDimension(b, &dim)); 63 CeedCall(CeedBasisGetNumComponents(b, &num_comp)); 64 CeedCall(CeedBasisGetNumQuadratureComponents(b, &Q_comp)); 65 if (r != CEED_ELEMRESTRICTION_NONE && restr_num_comp != num_comp) { 66 // LCOV_EXCL_START 67 return CeedError(ceed, CEED_ERROR_DIMENSION, 68 "Field '%s' of size %" CeedInt_FMT " and EvalMode %s: ElemRestriction has %" CeedInt_FMT 69 " components, but Basis has %" CeedInt_FMT " components", 70 qf_field->field_name, qf_field->size, CeedEvalModes[qf_field->eval_mode], restr_num_comp, num_comp); 71 // LCOV_EXCL_STOP 72 } 73 } else if (eval_mode != CEED_EVAL_NONE) { 74 // LCOV_EXCL_START 75 return CeedError(ceed, CEED_ERROR_INCOMPATIBLE, "Field '%s' configured with %s cannot be used with CEED_BASIS_COLLOCATED", qf_field->field_name, 76 CeedEvalModes[eval_mode]); 77 // LCOV_EXCL_STOP 78 } 79 // Field size 80 switch (eval_mode) { 81 case CEED_EVAL_NONE: 82 if (size != restr_num_comp) { 83 // LCOV_EXCL_START 84 return CeedError(ceed, CEED_ERROR_DIMENSION, 85 "Field '%s' of size %" CeedInt_FMT " and EvalMode %s: ElemRestriction has %" CeedInt_FMT " components", qf_field->field_name, 86 qf_field->size, CeedEvalModes[qf_field->eval_mode], restr_num_comp); 87 // LCOV_EXCL_STOP 88 } 89 break; 90 case CEED_EVAL_INTERP: 91 if (size != num_comp * Q_comp) { 92 // LCOV_EXCL_START 93 return CeedError(ceed, CEED_ERROR_DIMENSION, 94 "Field '%s' of size %" CeedInt_FMT " and EvalMode %s: ElemRestriction/Basis has %" CeedInt_FMT " components", 95 qf_field->field_name, qf_field->size, CeedEvalModes[qf_field->eval_mode], num_comp * Q_comp); 96 // LCOV_EXCL_STOP 97 } 98 break; 99 case CEED_EVAL_GRAD: 100 if (size != num_comp * dim) { 101 // LCOV_EXCL_START 102 return CeedError(ceed, CEED_ERROR_DIMENSION, 103 "Field '%s' of size %" CeedInt_FMT " and EvalMode %s in %" CeedInt_FMT " dimensions: ElemRestriction/Basis has %" CeedInt_FMT 104 " components", 105 qf_field->field_name, qf_field->size, CeedEvalModes[qf_field->eval_mode], dim, num_comp); 106 // LCOV_EXCL_STOP 107 } 108 break; 109 case CEED_EVAL_WEIGHT: 110 // No additional checks required 111 break; 112 case CEED_EVAL_DIV: 113 if (size != num_comp) { 114 // LCOV_EXCL_START 115 return CeedError(ceed, CEED_ERROR_DIMENSION, 116 "Field '%s' of size %" CeedInt_FMT " and EvalMode %s: ElemRestriction/Basis has %" CeedInt_FMT " components", 117 qf_field->field_name, qf_field->size, CeedEvalModes[qf_field->eval_mode], num_comp); 118 // LCOV_EXCL_STOP 119 } 120 break; 121 case CEED_EVAL_CURL: 122 // Not implemented 123 break; 124 } 125 return CEED_ERROR_SUCCESS; 126 } 127 128 /** 129 @brief View a field of a CeedOperator 130 131 @param[in] field Operator field to view 132 @param[in] qf_field QFunction field (carries field name) 133 @param[in] field_number Number of field being viewed 134 @param[in] sub true indicates sub-operator, which increases indentation; false for top-level operator 135 @param[in] input true for an input field; false for output field 136 @param[in] stream Stream to view to, e.g., stdout 137 138 @return An error code: 0 - success, otherwise - failure 139 140 @ref Utility 141 **/ 142 static int CeedOperatorFieldView(CeedOperatorField field, CeedQFunctionField qf_field, CeedInt field_number, bool sub, bool input, FILE *stream) { 143 const char *pre = sub ? " " : ""; 144 const char *in_out = input ? "Input" : "Output"; 145 146 fprintf(stream, 147 "%s %s field %" CeedInt_FMT 148 ":\n" 149 "%s Name: \"%s\"\n", 150 pre, in_out, field_number, pre, qf_field->field_name); 151 152 fprintf(stream, "%s Size: %" CeedInt_FMT "\n", pre, qf_field->size); 153 154 fprintf(stream, "%s EvalMode: %s\n", pre, CeedEvalModes[qf_field->eval_mode]); 155 156 if (field->basis == CEED_BASIS_COLLOCATED) fprintf(stream, "%s Collocated basis\n", pre); 157 158 if (field->vec == CEED_VECTOR_ACTIVE) fprintf(stream, "%s Active vector\n", pre); 159 else if (field->vec == CEED_VECTOR_NONE) fprintf(stream, "%s No vector\n", pre); 160 161 return CEED_ERROR_SUCCESS; 162 } 163 164 /** 165 @brief View a single CeedOperator 166 167 @param[in] op CeedOperator to view 168 @param[in] sub Boolean flag for sub-operator 169 @param[in] stream Stream to write; typically stdout/stderr or a file 170 171 @return Error code: 0 - success, otherwise - failure 172 173 @ref Utility 174 **/ 175 int CeedOperatorSingleView(CeedOperator op, bool sub, FILE *stream) { 176 const char *pre = sub ? " " : ""; 177 178 CeedInt num_elem, num_qpts; 179 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 180 CeedCall(CeedOperatorGetNumQuadraturePoints(op, &num_qpts)); 181 182 CeedInt total_fields = 0; 183 CeedCall(CeedOperatorGetNumArgs(op, &total_fields)); 184 fprintf(stream, "%s %" CeedInt_FMT " elements with %" CeedInt_FMT " quadrature points each\n", pre, num_elem, num_qpts); 185 186 fprintf(stream, "%s %" CeedInt_FMT " field%s\n", pre, total_fields, total_fields > 1 ? "s" : ""); 187 188 fprintf(stream, "%s %" CeedInt_FMT " input field%s:\n", pre, op->qf->num_input_fields, op->qf->num_input_fields > 1 ? "s" : ""); 189 for (CeedInt i = 0; i < op->qf->num_input_fields; i++) { 190 CeedCall(CeedOperatorFieldView(op->input_fields[i], op->qf->input_fields[i], i, sub, 1, stream)); 191 } 192 193 fprintf(stream, "%s %" CeedInt_FMT " output field%s:\n", pre, op->qf->num_output_fields, op->qf->num_output_fields > 1 ? "s" : ""); 194 for (CeedInt i = 0; i < op->qf->num_output_fields; i++) { 195 CeedCall(CeedOperatorFieldView(op->output_fields[i], op->qf->output_fields[i], i, sub, 0, stream)); 196 } 197 return CEED_ERROR_SUCCESS; 198 } 199 200 /** 201 @brief Find the active vector basis for a non-composite CeedOperator 202 203 @param[in] op CeedOperator to find active basis for 204 @param[out] active_basis Basis for active input vector or NULL for composite operator 205 206 @return An error code: 0 - success, otherwise - failure 207 208 @ref Developer 209 **/ 210 int CeedOperatorGetActiveBasis(CeedOperator op, CeedBasis *active_basis) { 211 Ceed ceed; 212 CeedCall(CeedOperatorGetCeed(op, &ceed)); 213 214 *active_basis = NULL; 215 if (op->is_composite) return CEED_ERROR_SUCCESS; 216 for (CeedInt i = 0; i < op->qf->num_input_fields; i++) { 217 if (op->input_fields[i]->vec == CEED_VECTOR_ACTIVE) { 218 if (*active_basis && *active_basis != op->input_fields[i]->basis) { 219 // LCOV_EXCL_START 220 return CeedError(ceed, CEED_ERROR_MINOR, "Multiple active CeedBases found"); 221 // LCOV_EXCL_STOP 222 } 223 *active_basis = op->input_fields[i]->basis; 224 } 225 } 226 227 if (!*active_basis) { 228 // LCOV_EXCL_START 229 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "No active CeedBasis found"); 230 // LCOV_EXCL_STOP 231 } 232 return CEED_ERROR_SUCCESS; 233 } 234 235 /** 236 @brief Find the active vector ElemRestriction for a non-composite CeedOperator 237 238 @param[in] op CeedOperator to find active basis for 239 @param[out] active_rstr ElemRestriction for active input vector or NULL for composite operator 240 241 @return An error code: 0 - success, otherwise - failure 242 243 @ref Utility 244 **/ 245 int CeedOperatorGetActiveElemRestriction(CeedOperator op, CeedElemRestriction *active_rstr) { 246 Ceed ceed; 247 CeedCall(CeedOperatorGetCeed(op, &ceed)); 248 249 *active_rstr = NULL; 250 if (op->is_composite) return CEED_ERROR_SUCCESS; 251 for (CeedInt i = 0; i < op->qf->num_input_fields; i++) { 252 if (op->input_fields[i]->vec == CEED_VECTOR_ACTIVE) { 253 if (*active_rstr && *active_rstr != op->input_fields[i]->elem_rstr) { 254 // LCOV_EXCL_START 255 return CeedError(ceed, CEED_ERROR_MINOR, "Multiple active CeedElemRestrictions found"); 256 // LCOV_EXCL_STOP 257 } 258 *active_rstr = op->input_fields[i]->elem_rstr; 259 } 260 } 261 262 if (!*active_rstr) { 263 // LCOV_EXCL_START 264 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "No active CeedElemRestriction found"); 265 // LCOV_EXCL_STOP 266 } 267 return CEED_ERROR_SUCCESS; 268 } 269 270 /** 271 @brief Set QFunctionContext field values of the specified type. 272 For composite operators, the value is set in all sub-operator QFunctionContexts that have a matching `field_name`. 273 A non-zero error code is returned for single operators that do not have a matching field of the same type or composite operators that do 274 not have any field of a matching type. 275 276 @param[in,out] op CeedOperator 277 @param[in] field_label Label of field to set 278 @param[in] field_type Type of field to set 279 @param[in] values Values to set 280 281 @return An error code: 0 - success, otherwise - failure 282 283 @ref User 284 **/ 285 static int CeedOperatorContextSetGeneric(CeedOperator op, CeedContextFieldLabel field_label, CeedContextFieldType field_type, void *values) { 286 if (!field_label) { 287 // LCOV_EXCL_START 288 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "Invalid field label"); 289 // LCOV_EXCL_STOP 290 } 291 292 // Check if field_label and op correspond 293 if (field_label->from_op) { 294 CeedInt index = -1; 295 296 for (CeedInt i = 0; i < op->num_context_labels; i++) { 297 if (op->context_labels[i] == field_label) index = i; 298 } 299 if (index == -1) { 300 // LCOV_EXCL_START 301 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "ContextFieldLabel does not correspond to the operator"); 302 // LCOV_EXCL_STOP 303 } 304 } 305 306 bool is_composite = false; 307 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 308 if (is_composite) { 309 CeedInt num_sub; 310 CeedOperator *sub_operators; 311 312 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub)); 313 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 314 if (num_sub != field_label->num_sub_labels) { 315 // LCOV_EXCL_START 316 CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "Composite operator modified after ContextFieldLabel created"); 317 // LCOV_EXCL_STOP 318 } 319 320 for (CeedInt i = 0; i < num_sub; i++) { 321 // Try every sub-operator, ok if some sub-operators do not have field 322 if (field_label->sub_labels[i] && sub_operators[i]->qf->ctx) { 323 CeedCall(CeedQFunctionContextSetGeneric(sub_operators[i]->qf->ctx, field_label->sub_labels[i], field_type, values)); 324 } 325 } 326 } else { 327 if (!op->qf->ctx) { 328 // LCOV_EXCL_START 329 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "QFunction does not have context data"); 330 // LCOV_EXCL_STOP 331 } 332 333 CeedCall(CeedQFunctionContextSetGeneric(op->qf->ctx, field_label, field_type, values)); 334 } 335 336 return CEED_ERROR_SUCCESS; 337 } 338 339 /** 340 @brief Get QFunctionContext field values of the specified type, read-only. 341 For composite operators, the values retrieved are for the first sub-operator QFunctionContext that have a matching `field_name`. 342 A non-zero error code is returned for single operators that do not have a matching field of the same type or composite operators that 343 do not have any field of a matching type. 344 345 @param[in,out] op CeedOperator 346 @param[in] field_label Label of field to set 347 @param[in] field_type Type of field to set 348 @param[out] num_values Number of values of type `field_type` in array `values` 349 @param[out] values Values in the label 350 351 @return An error code: 0 - success, otherwise - failure 352 353 @ref User 354 **/ 355 static int CeedOperatorContextGetGenericRead(CeedOperator op, CeedContextFieldLabel field_label, CeedContextFieldType field_type, size_t *num_values, 356 void *values) { 357 if (!field_label) { 358 // LCOV_EXCL_START 359 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "Invalid field label"); 360 // LCOV_EXCL_STOP 361 } 362 363 *(void **)values = NULL; 364 *num_values = 0; 365 366 // Check if field_label and op correspond 367 if (field_label->from_op) { 368 CeedInt index = -1; 369 370 for (CeedInt i = 0; i < op->num_context_labels; i++) { 371 if (op->context_labels[i] == field_label) index = i; 372 } 373 if (index == -1) { 374 // LCOV_EXCL_START 375 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "ContextFieldLabel does not correspond to the operator"); 376 // LCOV_EXCL_STOP 377 } 378 } 379 380 bool is_composite = false; 381 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 382 if (is_composite) { 383 CeedInt num_sub; 384 CeedOperator *sub_operators; 385 386 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub)); 387 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 388 if (num_sub != field_label->num_sub_labels) { 389 // LCOV_EXCL_START 390 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "Composite operator modified after ContextFieldLabel created"); 391 // LCOV_EXCL_STOP 392 } 393 394 for (CeedInt i = 0; i < num_sub; i++) { 395 // Try every sub-operator, ok if some sub-operators do not have field 396 if (field_label->sub_labels[i] && sub_operators[i]->qf->ctx) { 397 CeedCall(CeedQFunctionContextGetGenericRead(sub_operators[i]->qf->ctx, field_label->sub_labels[i], field_type, num_values, values)); 398 return CEED_ERROR_SUCCESS; 399 } 400 } 401 } else { 402 if (!op->qf->ctx) { 403 // LCOV_EXCL_START 404 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "QFunction does not have context data"); 405 // LCOV_EXCL_STOP 406 } 407 408 CeedCall(CeedQFunctionContextGetGenericRead(op->qf->ctx, field_label, field_type, num_values, values)); 409 } 410 411 return CEED_ERROR_SUCCESS; 412 } 413 414 /** 415 @brief Restore QFunctionContext field values of the specified type, read-only. 416 For composite operators, the values restored are for the first sub-operator QFunctionContext that have a matching `field_name`. 417 A non-zero error code is returned for single operators that do not have a matching field of the same type or composite operators 418 that do not have any field of a matching type. 419 420 @param[in,out] op CeedOperator 421 @param[in] field_label Label of field to set 422 @param[in] field_type Type of field to set 423 @param[in] values Values array to restore 424 425 @return An error code: 0 - success, otherwise - failure 426 427 @ref User 428 **/ 429 static int CeedOperatorContextRestoreGenericRead(CeedOperator op, CeedContextFieldLabel field_label, CeedContextFieldType field_type, void *values) { 430 if (!field_label) { 431 // LCOV_EXCL_START 432 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "Invalid field label"); 433 // LCOV_EXCL_STOP 434 } 435 436 // Check if field_label and op correspond 437 if (field_label->from_op) { 438 CeedInt index = -1; 439 440 for (CeedInt i = 0; i < op->num_context_labels; i++) { 441 if (op->context_labels[i] == field_label) index = i; 442 } 443 if (index == -1) { 444 // LCOV_EXCL_START 445 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "ContextFieldLabel does not correspond to the operator"); 446 // LCOV_EXCL_STOP 447 } 448 } 449 450 bool is_composite = false; 451 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 452 if (is_composite) { 453 CeedInt num_sub; 454 CeedOperator *sub_operators; 455 456 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub)); 457 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 458 if (num_sub != field_label->num_sub_labels) { 459 // LCOV_EXCL_START 460 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "Composite operator modified after ContextFieldLabel created"); 461 // LCOV_EXCL_STOP 462 } 463 464 for (CeedInt i = 0; i < num_sub; i++) { 465 // Try every sub-operator, ok if some sub-operators do not have field 466 if (field_label->sub_labels[i] && sub_operators[i]->qf->ctx) { 467 CeedCall(CeedQFunctionContextRestoreGenericRead(sub_operators[i]->qf->ctx, field_label->sub_labels[i], field_type, values)); 468 return CEED_ERROR_SUCCESS; 469 } 470 } 471 } else { 472 if (!op->qf->ctx) { 473 // LCOV_EXCL_START 474 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "QFunction does not have context data"); 475 // LCOV_EXCL_STOP 476 } 477 478 CeedCall(CeedQFunctionContextRestoreGenericRead(op->qf->ctx, field_label, field_type, values)); 479 } 480 481 return CEED_ERROR_SUCCESS; 482 } 483 484 /// @} 485 486 /// ---------------------------------------------------------------------------- 487 /// CeedOperator Backend API 488 /// ---------------------------------------------------------------------------- 489 /// @addtogroup CeedOperatorBackend 490 /// @{ 491 492 /** 493 @brief Get the number of arguments associated with a CeedOperator 494 495 @param[in] op CeedOperator 496 @param[out] num_args Variable to store vector number of arguments 497 498 @return An error code: 0 - success, otherwise - failure 499 500 @ref Backend 501 **/ 502 503 int CeedOperatorGetNumArgs(CeedOperator op, CeedInt *num_args) { 504 if (op->is_composite) { 505 // LCOV_EXCL_START 506 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not defined for composite operators"); 507 // LCOV_EXCL_STOP 508 } 509 *num_args = op->num_fields; 510 return CEED_ERROR_SUCCESS; 511 } 512 513 /** 514 @brief Get the setup status of a CeedOperator 515 516 @param[in] op CeedOperator 517 @param[out] is_setup_done Variable to store setup status 518 519 @return An error code: 0 - success, otherwise - failure 520 521 @ref Backend 522 **/ 523 524 int CeedOperatorIsSetupDone(CeedOperator op, bool *is_setup_done) { 525 *is_setup_done = op->is_backend_setup; 526 return CEED_ERROR_SUCCESS; 527 } 528 529 /** 530 @brief Get the QFunction associated with a CeedOperator 531 532 @param[in] op CeedOperator 533 @param[out] qf Variable to store QFunction 534 535 @return An error code: 0 - success, otherwise - failure 536 537 @ref Backend 538 **/ 539 540 int CeedOperatorGetQFunction(CeedOperator op, CeedQFunction *qf) { 541 if (op->is_composite) { 542 // LCOV_EXCL_START 543 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not defined for composite operator"); 544 // LCOV_EXCL_STOP 545 } 546 *qf = op->qf; 547 return CEED_ERROR_SUCCESS; 548 } 549 550 /** 551 @brief Get a boolean value indicating if the CeedOperator is composite 552 553 @param[in] op CeedOperator 554 @param[out] is_composite Variable to store composite status 555 556 @return An error code: 0 - success, otherwise - failure 557 558 @ref Backend 559 **/ 560 561 int CeedOperatorIsComposite(CeedOperator op, bool *is_composite) { 562 *is_composite = op->is_composite; 563 return CEED_ERROR_SUCCESS; 564 } 565 566 /** 567 @brief Get the backend data of a CeedOperator 568 569 @param[in] op CeedOperator 570 @param[out] data Variable to store data 571 572 @return An error code: 0 - success, otherwise - failure 573 574 @ref Backend 575 **/ 576 577 int CeedOperatorGetData(CeedOperator op, void *data) { 578 *(void **)data = op->data; 579 return CEED_ERROR_SUCCESS; 580 } 581 582 /** 583 @brief Set the backend data of a CeedOperator 584 585 @param[in,out] op CeedOperator 586 @param[in] data Data to set 587 588 @return An error code: 0 - success, otherwise - failure 589 590 @ref Backend 591 **/ 592 593 int CeedOperatorSetData(CeedOperator op, void *data) { 594 op->data = data; 595 return CEED_ERROR_SUCCESS; 596 } 597 598 /** 599 @brief Increment the reference counter for a CeedOperator 600 601 @param[in,out] op CeedOperator to increment the reference counter 602 603 @return An error code: 0 - success, otherwise - failure 604 605 @ref Backend 606 **/ 607 int CeedOperatorReference(CeedOperator op) { 608 op->ref_count++; 609 return CEED_ERROR_SUCCESS; 610 } 611 612 /** 613 @brief Set the setup flag of a CeedOperator to True 614 615 @param[in,out] op CeedOperator 616 617 @return An error code: 0 - success, otherwise - failure 618 619 @ref Backend 620 **/ 621 622 int CeedOperatorSetSetupDone(CeedOperator op) { 623 op->is_backend_setup = true; 624 return CEED_ERROR_SUCCESS; 625 } 626 627 /// @} 628 629 /// ---------------------------------------------------------------------------- 630 /// CeedOperator Public API 631 /// ---------------------------------------------------------------------------- 632 /// @addtogroup CeedOperatorUser 633 /// @{ 634 635 /** 636 @brief Create a CeedOperator and associate a CeedQFunction. 637 A CeedBasis and CeedElemRestriction can be associated with CeedQFunction fields with \ref CeedOperatorSetField. 638 639 @param[in] ceed Ceed object where the CeedOperator will be created 640 @param[in] qf QFunction defining the action of the operator at quadrature points 641 @param[in] dqf QFunction defining the action of the Jacobian of @a qf (or @ref CEED_QFUNCTION_NONE) 642 @param[in] dqfT QFunction defining the action of the transpose of the Jacobian of @a qf (or @ref CEED_QFUNCTION_NONE) 643 @param[out] op Address of the variable where the newly created CeedOperator will be stored 644 645 @return An error code: 0 - success, otherwise - failure 646 647 @ref User 648 */ 649 int CeedOperatorCreate(Ceed ceed, CeedQFunction qf, CeedQFunction dqf, CeedQFunction dqfT, CeedOperator *op) { 650 if (!ceed->OperatorCreate) { 651 Ceed delegate; 652 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Operator")); 653 654 if (!delegate) { 655 // LCOV_EXCL_START 656 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support OperatorCreate"); 657 // LCOV_EXCL_STOP 658 } 659 660 CeedCall(CeedOperatorCreate(delegate, qf, dqf, dqfT, op)); 661 return CEED_ERROR_SUCCESS; 662 } 663 664 if (!qf || qf == CEED_QFUNCTION_NONE) { 665 // LCOV_EXCL_START 666 return CeedError(ceed, CEED_ERROR_MINOR, "Operator must have a valid QFunction."); 667 // LCOV_EXCL_STOP 668 } 669 CeedCall(CeedCalloc(1, op)); 670 (*op)->ceed = ceed; 671 CeedCall(CeedReference(ceed)); 672 (*op)->ref_count = 1; 673 (*op)->qf = qf; 674 (*op)->input_size = -1; 675 (*op)->output_size = -1; 676 CeedCall(CeedQFunctionReference(qf)); 677 if (dqf && dqf != CEED_QFUNCTION_NONE) { 678 (*op)->dqf = dqf; 679 CeedCall(CeedQFunctionReference(dqf)); 680 } 681 if (dqfT && dqfT != CEED_QFUNCTION_NONE) { 682 (*op)->dqfT = dqfT; 683 CeedCall(CeedQFunctionReference(dqfT)); 684 } 685 CeedCall(CeedQFunctionAssemblyDataCreate(ceed, &(*op)->qf_assembled)); 686 CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*op)->input_fields)); 687 CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*op)->output_fields)); 688 CeedCall(ceed->OperatorCreate(*op)); 689 return CEED_ERROR_SUCCESS; 690 } 691 692 /** 693 @brief Create an operator that composes the action of several operators 694 695 @param[in] ceed Ceed object where the CeedOperator will be created 696 @param[out] op Address of the variable where the newly created Composite CeedOperator will be stored 697 698 @return An error code: 0 - success, otherwise - failure 699 700 @ref User 701 */ 702 int CeedCompositeOperatorCreate(Ceed ceed, CeedOperator *op) { 703 if (!ceed->CompositeOperatorCreate) { 704 Ceed delegate; 705 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Operator")); 706 707 if (delegate) { 708 CeedCall(CeedCompositeOperatorCreate(delegate, op)); 709 return CEED_ERROR_SUCCESS; 710 } 711 } 712 713 CeedCall(CeedCalloc(1, op)); 714 (*op)->ceed = ceed; 715 CeedCall(CeedReference(ceed)); 716 (*op)->ref_count = 1; 717 (*op)->is_composite = true; 718 CeedCall(CeedCalloc(CEED_COMPOSITE_MAX, &(*op)->sub_operators)); 719 (*op)->input_size = -1; 720 (*op)->output_size = -1; 721 722 if (ceed->CompositeOperatorCreate) { 723 CeedCall(ceed->CompositeOperatorCreate(*op)); 724 } 725 return CEED_ERROR_SUCCESS; 726 } 727 728 /** 729 @brief Copy the pointer to a CeedOperator. 730 Both pointers should be destroyed with `CeedOperatorDestroy()`. 731 732 Note: If the value of `op_copy` passed to this function is non-NULL, then it is assumed that `op_copy` is a pointer to a CeedOperator. 733 This CeedOperator will be destroyed if `op_copy` is the only reference to this CeedOperator. 734 735 @param[in] op CeedOperator to copy reference to 736 @param[in,out] op_copy Variable to store copied reference 737 738 @return An error code: 0 - success, otherwise - failure 739 740 @ref User 741 **/ 742 int CeedOperatorReferenceCopy(CeedOperator op, CeedOperator *op_copy) { 743 CeedCall(CeedOperatorReference(op)); 744 CeedCall(CeedOperatorDestroy(op_copy)); 745 *op_copy = op; 746 return CEED_ERROR_SUCCESS; 747 } 748 749 /** 750 @brief Provide a field to a CeedOperator for use by its CeedQFunction. 751 752 This function is used to specify both active and passive fields to a CeedOperator. 753 For passive fields, a vector @arg v must be provided. 754 Passive fields can inputs or outputs (updated in-place when operator is applied). 755 756 Active fields must be specified using this function, but their data (in a CeedVector) is passed in CeedOperatorApply(). 757 There can be at most one active input CeedVector and at most one active output CeedVector passed to CeedOperatorApply(). 758 759 @param[in,out] op CeedOperator on which to provide the field 760 @param[in] field_name Name of the field (to be matched with the name used by CeedQFunction) 761 @param[in] r CeedElemRestriction 762 @param[in] b CeedBasis in which the field resides or @ref CEED_BASIS_COLLOCATED if collocated with quadrature points 763 @param[in] v CeedVector to be used by CeedOperator or @ref CEED_VECTOR_ACTIVE if field is active or @ref CEED_VECTOR_NONE 764 if using @ref CEED_EVAL_WEIGHT in the QFunction 765 766 @return An error code: 0 - success, otherwise - failure 767 768 @ref User 769 **/ 770 int CeedOperatorSetField(CeedOperator op, const char *field_name, CeedElemRestriction r, CeedBasis b, CeedVector v) { 771 if (op->is_composite) { 772 // LCOV_EXCL_START 773 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "Cannot add field to composite operator."); 774 // LCOV_EXCL_STOP 775 } 776 if (op->is_immutable) { 777 // LCOV_EXCL_START 778 return CeedError(op->ceed, CEED_ERROR_MAJOR, "Operator cannot be changed after set as immutable"); 779 // LCOV_EXCL_STOP 780 } 781 if (!r) { 782 // LCOV_EXCL_START 783 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "ElemRestriction r for field \"%s\" must be non-NULL.", field_name); 784 // LCOV_EXCL_STOP 785 } 786 if (!b) { 787 // LCOV_EXCL_START 788 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "Basis b for field \"%s\" must be non-NULL.", field_name); 789 // LCOV_EXCL_STOP 790 } 791 if (!v) { 792 // LCOV_EXCL_START 793 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "Vector v for field \"%s\" must be non-NULL.", field_name); 794 // LCOV_EXCL_STOP 795 } 796 797 CeedInt num_elem; 798 CeedCall(CeedElemRestrictionGetNumElements(r, &num_elem)); 799 if (r != CEED_ELEMRESTRICTION_NONE && op->has_restriction && op->num_elem != num_elem) { 800 // LCOV_EXCL_START 801 return CeedError(op->ceed, CEED_ERROR_DIMENSION, 802 "ElemRestriction with %" CeedInt_FMT " elements incompatible with prior %" CeedInt_FMT " elements", num_elem, op->num_elem); 803 // LCOV_EXCL_STOP 804 } 805 806 CeedInt num_qpts = 0; 807 if (b != CEED_BASIS_COLLOCATED) { 808 CeedCall(CeedBasisGetNumQuadraturePoints(b, &num_qpts)); 809 if (op->num_qpts && op->num_qpts != num_qpts) { 810 // LCOV_EXCL_START 811 return CeedError(op->ceed, CEED_ERROR_DIMENSION, 812 "Basis with %" CeedInt_FMT " quadrature points incompatible with prior %" CeedInt_FMT " points", num_qpts, op->num_qpts); 813 // LCOV_EXCL_STOP 814 } 815 } 816 CeedQFunctionField qf_field; 817 CeedOperatorField *op_field; 818 bool is_input = true; 819 for (CeedInt i = 0; i < op->qf->num_input_fields; i++) { 820 if (!strcmp(field_name, (*op->qf->input_fields[i]).field_name)) { 821 qf_field = op->qf->input_fields[i]; 822 op_field = &op->input_fields[i]; 823 goto found; 824 } 825 } 826 is_input = false; 827 for (CeedInt i = 0; i < op->qf->num_output_fields; i++) { 828 if (!strcmp(field_name, (*op->qf->output_fields[i]).field_name)) { 829 qf_field = op->qf->output_fields[i]; 830 op_field = &op->output_fields[i]; 831 goto found; 832 } 833 } 834 // LCOV_EXCL_START 835 return CeedError(op->ceed, CEED_ERROR_INCOMPLETE, "QFunction has no knowledge of field '%s'", field_name); 836 // LCOV_EXCL_STOP 837 found: 838 CeedCall(CeedOperatorCheckField(op->ceed, qf_field, r, b)); 839 CeedCall(CeedCalloc(1, op_field)); 840 841 if (v == CEED_VECTOR_ACTIVE) { 842 CeedSize l_size; 843 CeedCall(CeedElemRestrictionGetLVectorSize(r, &l_size)); 844 if (is_input) { 845 if (op->input_size == -1) op->input_size = l_size; 846 if (l_size != op->input_size) { 847 // LCOV_EXCL_START 848 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "LVector size %td does not match previous size %td", l_size, op->input_size); 849 // LCOV_EXCL_STOP 850 } 851 } else { 852 if (op->output_size == -1) op->output_size = l_size; 853 if (l_size != op->output_size) { 854 // LCOV_EXCL_START 855 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "LVector size %td does not match previous size %td", l_size, op->output_size); 856 // LCOV_EXCL_STOP 857 } 858 } 859 } 860 861 (*op_field)->vec = v; 862 if (v != CEED_VECTOR_ACTIVE && v != CEED_VECTOR_NONE) { 863 CeedCall(CeedVectorReference(v)); 864 } 865 866 (*op_field)->elem_rstr = r; 867 CeedCall(CeedElemRestrictionReference(r)); 868 if (r != CEED_ELEMRESTRICTION_NONE) { 869 op->num_elem = num_elem; 870 op->has_restriction = true; // Restriction set, but num_elem may be 0 871 } 872 873 (*op_field)->basis = b; 874 if (b != CEED_BASIS_COLLOCATED) { 875 if (!op->num_qpts) { 876 CeedCall(CeedOperatorSetNumQuadraturePoints(op, num_qpts)); 877 } 878 CeedCall(CeedBasisReference(b)); 879 } 880 881 op->num_fields += 1; 882 CeedCall(CeedStringAllocCopy(field_name, (char **)&(*op_field)->field_name)); 883 return CEED_ERROR_SUCCESS; 884 } 885 886 /** 887 @brief Get the CeedOperatorFields of a CeedOperator 888 889 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 890 891 @param[in] op CeedOperator 892 @param[out] num_input_fields Variable to store number of input fields 893 @param[out] input_fields Variable to store input_fields 894 @param[out] num_output_fields Variable to store number of output fields 895 @param[out] output_fields Variable to store output_fields 896 897 @return An error code: 0 - success, otherwise - failure 898 899 @ref Advanced 900 **/ 901 int CeedOperatorGetFields(CeedOperator op, CeedInt *num_input_fields, CeedOperatorField **input_fields, CeedInt *num_output_fields, 902 CeedOperatorField **output_fields) { 903 if (op->is_composite) { 904 // LCOV_EXCL_START 905 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not defined for composite operator"); 906 // LCOV_EXCL_STOP 907 } 908 CeedCall(CeedOperatorCheckReady(op)); 909 910 if (num_input_fields) *num_input_fields = op->qf->num_input_fields; 911 if (input_fields) *input_fields = op->input_fields; 912 if (num_output_fields) *num_output_fields = op->qf->num_output_fields; 913 if (output_fields) *output_fields = op->output_fields; 914 return CEED_ERROR_SUCCESS; 915 } 916 917 /** 918 @brief Get a CeedOperatorField of an CeedOperator from its name 919 920 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 921 922 @param[in] op CeedOperator 923 @param[in] field_name Name of desired CeedOperatorField 924 @param[out] op_field CeedOperatorField corresponding to the name 925 926 @return An error code: 0 - success, otherwise - failure 927 928 @ref Advanced 929 **/ 930 int CeedOperatorGetFieldByName(CeedOperator op, const char *field_name, CeedOperatorField *op_field) { 931 CeedInt num_input_fields, num_output_fields; 932 CeedOperatorField *input_fields, *output_fields; 933 char *name; 934 CeedCall(CeedOperatorGetFields(op, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); 935 936 for (CeedInt i = 0; i < num_input_fields; i++) { 937 CeedCall(CeedOperatorFieldGetName(input_fields[i], &name)); 938 if (!strcmp(name, field_name)) { 939 *op_field = input_fields[i]; 940 return CEED_ERROR_SUCCESS; 941 } 942 } 943 944 for (CeedInt i = 0; i < num_output_fields; i++) { 945 CeedCall(CeedOperatorFieldGetName(output_fields[i], &name)); 946 if (!strcmp(name, field_name)) { 947 *op_field = output_fields[i]; 948 return CEED_ERROR_SUCCESS; 949 } 950 } 951 952 bool has_name = op->name; 953 // LCOV_EXCL_START 954 return CeedError(op->ceed, CEED_ERROR_MINOR, "The field \"%s\" not found in CeedOperator%s%s%s.\n", field_name, has_name ? " \"" : "", 955 has_name ? op->name : "", has_name ? "\"" : ""); 956 // LCOV_EXCL_STOP 957 } 958 959 /** 960 @brief Get the name of a CeedOperatorField 961 962 @param[in] op_field CeedOperatorField 963 @param[out] field_name Variable to store the field name 964 965 @return An error code: 0 - success, otherwise - failure 966 967 @ref Advanced 968 **/ 969 int CeedOperatorFieldGetName(CeedOperatorField op_field, char **field_name) { 970 *field_name = (char *)op_field->field_name; 971 return CEED_ERROR_SUCCESS; 972 } 973 974 /** 975 @brief Get the CeedElemRestriction of a CeedOperatorField 976 977 @param[in] op_field CeedOperatorField 978 @param[out] rstr Variable to store CeedElemRestriction 979 980 @return An error code: 0 - success, otherwise - failure 981 982 @ref Advanced 983 **/ 984 int CeedOperatorFieldGetElemRestriction(CeedOperatorField op_field, CeedElemRestriction *rstr) { 985 *rstr = op_field->elem_rstr; 986 return CEED_ERROR_SUCCESS; 987 } 988 989 /** 990 @brief Get the CeedBasis of a CeedOperatorField 991 992 @param[in] op_field CeedOperatorField 993 @param[out] basis Variable to store CeedBasis 994 995 @return An error code: 0 - success, otherwise - failure 996 997 @ref Advanced 998 **/ 999 int CeedOperatorFieldGetBasis(CeedOperatorField op_field, CeedBasis *basis) { 1000 *basis = op_field->basis; 1001 return CEED_ERROR_SUCCESS; 1002 } 1003 1004 /** 1005 @brief Get the CeedVector of a CeedOperatorField 1006 1007 @param[in] op_field CeedOperatorField 1008 @param[out] vec Variable to store CeedVector 1009 1010 @return An error code: 0 - success, otherwise - failure 1011 1012 @ref Advanced 1013 **/ 1014 int CeedOperatorFieldGetVector(CeedOperatorField op_field, CeedVector *vec) { 1015 *vec = op_field->vec; 1016 return CEED_ERROR_SUCCESS; 1017 } 1018 1019 /** 1020 @brief Add a sub-operator to a composite CeedOperator 1021 1022 @param[in,out] composite_op Composite CeedOperator 1023 @param[in] sub_op Sub-operator CeedOperator 1024 1025 @return An error code: 0 - success, otherwise - failure 1026 1027 @ref User 1028 */ 1029 int CeedCompositeOperatorAddSub(CeedOperator composite_op, CeedOperator sub_op) { 1030 if (!composite_op->is_composite) { 1031 // LCOV_EXCL_START 1032 return CeedError(composite_op->ceed, CEED_ERROR_MINOR, "CeedOperator is not a composite operator"); 1033 // LCOV_EXCL_STOP 1034 } 1035 1036 if (composite_op->num_suboperators == CEED_COMPOSITE_MAX) { 1037 // LCOV_EXCL_START 1038 return CeedError(composite_op->ceed, CEED_ERROR_UNSUPPORTED, "Cannot add additional sub-operators"); 1039 // LCOV_EXCL_STOP 1040 } 1041 if (composite_op->is_immutable) { 1042 // LCOV_EXCL_START 1043 return CeedError(composite_op->ceed, CEED_ERROR_MAJOR, "Operator cannot be changed after set as immutable"); 1044 // LCOV_EXCL_STOP 1045 } 1046 1047 { 1048 CeedSize input_size, output_size; 1049 CeedCall(CeedOperatorGetActiveVectorLengths(sub_op, &input_size, &output_size)); 1050 if (composite_op->input_size == -1) composite_op->input_size = input_size; 1051 if (composite_op->output_size == -1) composite_op->output_size = output_size; 1052 // Note, a size of -1 means no active vector restriction set, so no incompatibility 1053 if ((input_size != -1 && input_size != composite_op->input_size) || (output_size != -1 && output_size != composite_op->output_size)) { 1054 // LCOV_EXCL_START 1055 return CeedError(composite_op->ceed, CEED_ERROR_MAJOR, 1056 "Sub-operators must have compatible dimensions; composite operator of shape (%td, %td) not compatible with sub-operator of " 1057 "shape (%td, %td)", 1058 composite_op->input_size, composite_op->output_size, input_size, output_size); 1059 // LCOV_EXCL_STOP 1060 } 1061 } 1062 1063 composite_op->sub_operators[composite_op->num_suboperators] = sub_op; 1064 CeedCall(CeedOperatorReference(sub_op)); 1065 composite_op->num_suboperators++; 1066 return CEED_ERROR_SUCCESS; 1067 } 1068 1069 /** 1070 @brief Get the number of sub_operators associated with a CeedOperator 1071 1072 @param[in] op CeedOperator 1073 @param[out] num_suboperators Variable to store number of sub_operators 1074 1075 @return An error code: 0 - success, otherwise - failure 1076 1077 @ref Backend 1078 **/ 1079 1080 int CeedCompositeOperatorGetNumSub(CeedOperator op, CeedInt *num_suboperators) { 1081 if (!op->is_composite) { 1082 // LCOV_EXCL_START 1083 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not a composite operator"); 1084 // LCOV_EXCL_STOP 1085 } 1086 *num_suboperators = op->num_suboperators; 1087 return CEED_ERROR_SUCCESS; 1088 } 1089 1090 /** 1091 @brief Get the list of sub_operators associated with a CeedOperator 1092 1093 @param op CeedOperator 1094 @param[out] sub_operators Variable to store list of sub_operators 1095 1096 @return An error code: 0 - success, otherwise - failure 1097 1098 @ref Backend 1099 **/ 1100 1101 int CeedCompositeOperatorGetSubList(CeedOperator op, CeedOperator **sub_operators) { 1102 if (!op->is_composite) { 1103 // LCOV_EXCL_START 1104 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not a composite operator"); 1105 // LCOV_EXCL_STOP 1106 } 1107 *sub_operators = op->sub_operators; 1108 return CEED_ERROR_SUCCESS; 1109 } 1110 1111 /** 1112 @brief Check if a CeedOperator is ready to be used. 1113 1114 @param[in] op CeedOperator to check 1115 1116 @return An error code: 0 - success, otherwise - failure 1117 1118 @ref User 1119 **/ 1120 int CeedOperatorCheckReady(CeedOperator op) { 1121 Ceed ceed; 1122 CeedCall(CeedOperatorGetCeed(op, &ceed)); 1123 1124 if (op->is_interface_setup) return CEED_ERROR_SUCCESS; 1125 1126 CeedQFunction qf = op->qf; 1127 if (op->is_composite) { 1128 if (!op->num_suboperators) { 1129 // Empty operator setup 1130 op->input_size = 0; 1131 op->output_size = 0; 1132 } else { 1133 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1134 CeedCall(CeedOperatorCheckReady(op->sub_operators[i])); 1135 } 1136 // Sub-operators could be modified after adding to composite operator 1137 // Need to verify no lvec incompatibility from any changes 1138 CeedSize input_size, output_size; 1139 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 1140 } 1141 } else { 1142 if (op->num_fields == 0) { 1143 // LCOV_EXCL_START 1144 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "No operator fields set"); 1145 // LCOV_EXCL_STOP 1146 } 1147 if (op->num_fields < qf->num_input_fields + qf->num_output_fields) { 1148 // LCOV_EXCL_START 1149 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "Not all operator fields set"); 1150 // LCOV_EXCL_STOP 1151 } 1152 if (!op->has_restriction) { 1153 // LCOV_EXCL_START 1154 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "At least one restriction required"); 1155 // LCOV_EXCL_STOP 1156 } 1157 if (op->num_qpts == 0) { 1158 // LCOV_EXCL_START 1159 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "At least one non-collocated basis is required or the number of quadrature points must be set"); 1160 // LCOV_EXCL_STOP 1161 } 1162 } 1163 1164 // Flag as immutable and ready 1165 op->is_interface_setup = true; 1166 if (op->qf && op->qf != CEED_QFUNCTION_NONE) op->qf->is_immutable = true; 1167 if (op->dqf && op->dqf != CEED_QFUNCTION_NONE) op->dqf->is_immutable = true; 1168 if (op->dqfT && op->dqfT != CEED_QFUNCTION_NONE) op->dqfT->is_immutable = true; 1169 return CEED_ERROR_SUCCESS; 1170 } 1171 1172 /** 1173 @brief Get vector lengths for the active input and/or output vectors of a CeedOperator. 1174 Note: Lengths of -1 indicate that the CeedOperator does not have an active input and/or output. 1175 1176 @param[in] op CeedOperator 1177 @param[out] input_size Variable to store active input vector length, or NULL 1178 @param[out] output_size Variable to store active output vector length, or NULL 1179 1180 @return An error code: 0 - success, otherwise - failure 1181 1182 @ref User 1183 **/ 1184 int CeedOperatorGetActiveVectorLengths(CeedOperator op, CeedSize *input_size, CeedSize *output_size) { 1185 bool is_composite; 1186 1187 if (input_size) *input_size = op->input_size; 1188 if (output_size) *output_size = op->output_size; 1189 1190 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1191 if (is_composite && (op->input_size == -1 || op->output_size == -1)) { 1192 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1193 CeedSize sub_input_size, sub_output_size; 1194 CeedCall(CeedOperatorGetActiveVectorLengths(op->sub_operators[i], &sub_input_size, &sub_output_size)); 1195 if (op->input_size == -1) op->input_size = sub_input_size; 1196 if (op->output_size == -1) op->output_size = sub_output_size; 1197 // Note, a size of -1 means no active vector restriction set, so no incompatibility 1198 if ((sub_input_size != -1 && sub_input_size != op->input_size) || (sub_output_size != -1 && sub_output_size != op->output_size)) { 1199 // LCOV_EXCL_START 1200 return CeedError(op->ceed, CEED_ERROR_MAJOR, 1201 "Sub-operators must have compatible dimensions; composite operator of shape (%td, %td) not compatible with sub-operator of " 1202 "shape (%td, %td)", 1203 op->input_size, op->output_size, input_size, output_size); 1204 // LCOV_EXCL_STOP 1205 } 1206 } 1207 } 1208 1209 return CEED_ERROR_SUCCESS; 1210 } 1211 1212 /** 1213 @brief Set reuse of CeedQFunction data in CeedOperatorLinearAssemble* functions. 1214 When `reuse_assembly_data = false` (default), the CeedQFunction associated with this CeedOperator is re-assembled every time a 1215 `CeedOperatorLinearAssemble*` function is called. When `reuse_assembly_data = true`, the CeedQFunction associated with this CeedOperator is reused 1216 between calls to `CeedOperatorSetQFunctionAssemblyDataUpdated`. 1217 1218 @param[in] op CeedOperator 1219 @param[in] reuse_assembly_data Boolean flag setting assembly data reuse 1220 1221 @return An error code: 0 - success, otherwise - failure 1222 1223 @ref Advanced 1224 **/ 1225 int CeedOperatorSetQFunctionAssemblyReuse(CeedOperator op, bool reuse_assembly_data) { 1226 bool is_composite; 1227 1228 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1229 if (is_composite) { 1230 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1231 CeedCall(CeedOperatorSetQFunctionAssemblyReuse(op->sub_operators[i], reuse_assembly_data)); 1232 } 1233 } else { 1234 CeedCall(CeedQFunctionAssemblyDataSetReuse(op->qf_assembled, reuse_assembly_data)); 1235 } 1236 1237 return CEED_ERROR_SUCCESS; 1238 } 1239 1240 /** 1241 @brief Mark CeedQFunction data as updated and the CeedQFunction as requiring re-assembly. 1242 1243 @param[in] op CeedOperator 1244 @param[in] needs_data_update Boolean flag setting assembly data reuse 1245 1246 @return An error code: 0 - success, otherwise - failure 1247 1248 @ref Advanced 1249 **/ 1250 int CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(CeedOperator op, bool needs_data_update) { 1251 bool is_composite; 1252 1253 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1254 if (is_composite) { 1255 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1256 CeedCall(CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(op->sub_operators[i], needs_data_update)); 1257 } 1258 } else { 1259 CeedCall(CeedQFunctionAssemblyDataSetUpdateNeeded(op->qf_assembled, needs_data_update)); 1260 } 1261 1262 return CEED_ERROR_SUCCESS; 1263 } 1264 1265 /** 1266 @brief Set the number of quadrature points associated with a CeedOperator. 1267 This should be used when creating a CeedOperator where every field has a collocated basis. 1268 This function cannot be used for composite CeedOperators. 1269 1270 @param[in,out] op CeedOperator 1271 @param[in] num_qpts Number of quadrature points to set 1272 1273 @return An error code: 0 - success, otherwise - failure 1274 1275 @ref Advanced 1276 **/ 1277 int CeedOperatorSetNumQuadraturePoints(CeedOperator op, CeedInt num_qpts) { 1278 if (op->is_composite) { 1279 // LCOV_EXCL_START 1280 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not defined for composite operator"); 1281 // LCOV_EXCL_STOP 1282 } 1283 if (op->num_qpts) { 1284 // LCOV_EXCL_START 1285 return CeedError(op->ceed, CEED_ERROR_MINOR, "Number of quadrature points already defined"); 1286 // LCOV_EXCL_STOP 1287 } 1288 if (op->is_immutable) { 1289 // LCOV_EXCL_START 1290 return CeedError(op->ceed, CEED_ERROR_MAJOR, "Operator cannot be changed after set as immutable"); 1291 // LCOV_EXCL_STOP 1292 } 1293 op->num_qpts = num_qpts; 1294 return CEED_ERROR_SUCCESS; 1295 } 1296 1297 /** 1298 @brief Set name of CeedOperator for CeedOperatorView output 1299 1300 @param[in,out] op CeedOperator 1301 @param[in] name Name to set, or NULL to remove previously set name 1302 1303 @return An error code: 0 - success, otherwise - failure 1304 1305 @ref User 1306 **/ 1307 int CeedOperatorSetName(CeedOperator op, const char *name) { 1308 char *name_copy; 1309 size_t name_len = name ? strlen(name) : 0; 1310 1311 CeedCall(CeedFree(&op->name)); 1312 if (name_len > 0) { 1313 CeedCall(CeedCalloc(name_len + 1, &name_copy)); 1314 memcpy(name_copy, name, name_len); 1315 op->name = name_copy; 1316 } 1317 1318 return CEED_ERROR_SUCCESS; 1319 } 1320 1321 /** 1322 @brief View a CeedOperator 1323 1324 @param[in] op CeedOperator to view 1325 @param[in] stream Stream to write; typically stdout/stderr or a file 1326 1327 @return Error code: 0 - success, otherwise - failure 1328 1329 @ref User 1330 **/ 1331 int CeedOperatorView(CeedOperator op, FILE *stream) { 1332 bool has_name = op->name; 1333 1334 if (op->is_composite) { 1335 fprintf(stream, "Composite CeedOperator%s%s\n", has_name ? " - " : "", has_name ? op->name : ""); 1336 1337 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1338 has_name = op->sub_operators[i]->name; 1339 fprintf(stream, " SubOperator %" CeedInt_FMT "%s%s:\n", i, has_name ? " - " : "", has_name ? op->sub_operators[i]->name : ""); 1340 CeedCall(CeedOperatorSingleView(op->sub_operators[i], 1, stream)); 1341 } 1342 } else { 1343 fprintf(stream, "CeedOperator%s%s\n", has_name ? " - " : "", has_name ? op->name : ""); 1344 CeedCall(CeedOperatorSingleView(op, 0, stream)); 1345 } 1346 return CEED_ERROR_SUCCESS; 1347 } 1348 1349 /** 1350 @brief Get the Ceed associated with a CeedOperator 1351 1352 @param[in] op CeedOperator 1353 @param[out] ceed Variable to store Ceed 1354 1355 @return An error code: 0 - success, otherwise - failure 1356 1357 @ref Advanced 1358 **/ 1359 int CeedOperatorGetCeed(CeedOperator op, Ceed *ceed) { 1360 *ceed = op->ceed; 1361 return CEED_ERROR_SUCCESS; 1362 } 1363 1364 /** 1365 @brief Get the number of elements associated with a CeedOperator 1366 1367 @param[in] op CeedOperator 1368 @param[out] num_elem Variable to store number of elements 1369 1370 @return An error code: 0 - success, otherwise - failure 1371 1372 @ref Advanced 1373 **/ 1374 int CeedOperatorGetNumElements(CeedOperator op, CeedInt *num_elem) { 1375 if (op->is_composite) { 1376 // LCOV_EXCL_START 1377 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not defined for composite operator"); 1378 // LCOV_EXCL_STOP 1379 } 1380 *num_elem = op->num_elem; 1381 return CEED_ERROR_SUCCESS; 1382 } 1383 1384 /** 1385 @brief Get the number of quadrature points associated with a CeedOperator 1386 1387 @param[in] op CeedOperator 1388 @param[out] num_qpts Variable to store vector number of quadrature points 1389 1390 @return An error code: 0 - success, otherwise - failure 1391 1392 @ref Advanced 1393 **/ 1394 int CeedOperatorGetNumQuadraturePoints(CeedOperator op, CeedInt *num_qpts) { 1395 if (op->is_composite) { 1396 // LCOV_EXCL_START 1397 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not defined for composite operator"); 1398 // LCOV_EXCL_STOP 1399 } 1400 *num_qpts = op->num_qpts; 1401 return CEED_ERROR_SUCCESS; 1402 } 1403 1404 /** 1405 @brief Estimate number of FLOPs required to apply CeedOperator on the active vector 1406 1407 @param[in] op CeedOperator to estimate FLOPs for 1408 @param[out] flops Address of variable to hold FLOPs estimate 1409 1410 @ref Backend 1411 **/ 1412 int CeedOperatorGetFlopsEstimate(CeedOperator op, CeedSize *flops) { 1413 bool is_composite; 1414 CeedCall(CeedOperatorCheckReady(op)); 1415 1416 *flops = 0; 1417 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1418 if (is_composite) { 1419 CeedInt num_suboperators; 1420 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1421 CeedOperator *sub_operators; 1422 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1423 1424 // FLOPs for each suboperator 1425 for (CeedInt i = 0; i < num_suboperators; i++) { 1426 CeedSize suboperator_flops; 1427 CeedCall(CeedOperatorGetFlopsEstimate(sub_operators[i], &suboperator_flops)); 1428 *flops += suboperator_flops; 1429 } 1430 } else { 1431 CeedInt num_input_fields, num_output_fields; 1432 CeedOperatorField *input_fields, *output_fields; 1433 1434 CeedCall(CeedOperatorGetFields(op, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); 1435 1436 CeedInt num_elem = 0; 1437 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 1438 // Input FLOPs 1439 for (CeedInt i = 0; i < num_input_fields; i++) { 1440 if (input_fields[i]->vec == CEED_VECTOR_ACTIVE) { 1441 CeedSize restr_flops, basis_flops; 1442 1443 CeedCall(CeedElemRestrictionGetFlopsEstimate(input_fields[i]->elem_rstr, CEED_NOTRANSPOSE, &restr_flops)); 1444 *flops += restr_flops; 1445 CeedCall(CeedBasisGetFlopsEstimate(input_fields[i]->basis, CEED_NOTRANSPOSE, op->qf->input_fields[i]->eval_mode, &basis_flops)); 1446 *flops += basis_flops * num_elem; 1447 } 1448 } 1449 // QF FLOPs 1450 CeedInt num_qpts; 1451 CeedSize qf_flops; 1452 CeedCall(CeedOperatorGetNumQuadraturePoints(op, &num_qpts)); 1453 CeedCall(CeedQFunctionGetFlopsEstimate(op->qf, &qf_flops)); 1454 *flops += num_elem * num_qpts * qf_flops; 1455 // Output FLOPs 1456 for (CeedInt i = 0; i < num_output_fields; i++) { 1457 if (output_fields[i]->vec == CEED_VECTOR_ACTIVE) { 1458 CeedSize restr_flops, basis_flops; 1459 1460 CeedCall(CeedElemRestrictionGetFlopsEstimate(output_fields[i]->elem_rstr, CEED_TRANSPOSE, &restr_flops)); 1461 *flops += restr_flops; 1462 CeedCall(CeedBasisGetFlopsEstimate(output_fields[i]->basis, CEED_TRANSPOSE, op->qf->output_fields[i]->eval_mode, &basis_flops)); 1463 *flops += basis_flops * num_elem; 1464 } 1465 } 1466 } 1467 1468 return CEED_ERROR_SUCCESS; 1469 } 1470 1471 /** 1472 @brief Get CeedQFunction global context for a CeedOperator. 1473 The caller is responsible for destroying `ctx` returned from this function via `CeedQFunctionContextDestroy()`. 1474 1475 Note: If the value of `ctx` passed into this function is non-NULL, then it is assumed that `ctx` is a pointer to a 1476 CeedQFunctionContext. This CeedQFunctionContext will be destroyed if `ctx` is the only reference to this CeedQFunctionContext. 1477 1478 @param[in] op CeedOperator 1479 @param[out] ctx Variable to store CeedQFunctionContext 1480 1481 @return An error code: 0 - success, otherwise - failure 1482 1483 @ref Advanced 1484 **/ 1485 int CeedOperatorGetContext(CeedOperator op, CeedQFunctionContext *ctx) { 1486 bool is_composite; 1487 1488 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1489 if (is_composite) return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "Cannot retrieve QFunctionContext for composite operator"); 1490 1491 if (op->qf->ctx) CeedCall(CeedQFunctionContextReferenceCopy(op->qf->ctx, ctx)); 1492 else *ctx = NULL; 1493 return CEED_ERROR_SUCCESS; 1494 } 1495 1496 /** 1497 @brief Get label for a registered QFunctionContext field, or `NULL` if no field has been registered with this `field_name`. 1498 1499 Fields are registered via `CeedQFunctionContextRegister*()` functions (eg. `CeedQFunctionContextRegisterDouble()`). 1500 1501 @param[in] op CeedOperator 1502 @param[in] field_name Name of field to retrieve label 1503 @param[out] field_label Variable to field label 1504 1505 @return An error code: 0 - success, otherwise - failure 1506 1507 @ref User 1508 **/ 1509 int CeedOperatorGetContextFieldLabel(CeedOperator op, const char *field_name, CeedContextFieldLabel *field_label) { 1510 bool is_composite; 1511 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1512 1513 if (is_composite) { 1514 // Check if composite label already created 1515 for (CeedInt i = 0; i < op->num_context_labels; i++) { 1516 if (!strcmp(op->context_labels[i]->name, field_name)) { 1517 *field_label = op->context_labels[i]; 1518 return CEED_ERROR_SUCCESS; 1519 } 1520 } 1521 1522 // Create composite label if needed 1523 CeedInt num_sub; 1524 CeedOperator *sub_operators; 1525 CeedContextFieldLabel new_field_label; 1526 1527 CeedCall(CeedCalloc(1, &new_field_label)); 1528 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub)); 1529 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1530 CeedCall(CeedCalloc(num_sub, &new_field_label->sub_labels)); 1531 new_field_label->num_sub_labels = num_sub; 1532 1533 bool label_found = false; 1534 for (CeedInt i = 0; i < num_sub; i++) { 1535 if (sub_operators[i]->qf->ctx) { 1536 CeedContextFieldLabel new_field_label_i; 1537 CeedCall(CeedQFunctionContextGetFieldLabel(sub_operators[i]->qf->ctx, field_name, &new_field_label_i)); 1538 if (new_field_label_i) { 1539 label_found = true; 1540 new_field_label->sub_labels[i] = new_field_label_i; 1541 new_field_label->name = new_field_label_i->name; 1542 new_field_label->description = new_field_label_i->description; 1543 if (new_field_label->type && new_field_label->type != new_field_label_i->type) { 1544 // LCOV_EXCL_START 1545 CeedCall(CeedFree(&new_field_label)); 1546 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "Incompatible field types on sub-operator contexts. %s != %s", 1547 CeedContextFieldTypes[new_field_label->type], CeedContextFieldTypes[new_field_label_i->type]); 1548 // LCOV_EXCL_STOP 1549 } else { 1550 new_field_label->type = new_field_label_i->type; 1551 } 1552 if (new_field_label->num_values != 0 && new_field_label->num_values != new_field_label_i->num_values) { 1553 // LCOV_EXCL_START 1554 CeedCall(CeedFree(&new_field_label)); 1555 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "Incompatible field number of values on sub-operator contexts. %ld != %ld", 1556 new_field_label->num_values, new_field_label_i->num_values); 1557 // LCOV_EXCL_STOP 1558 } else { 1559 new_field_label->num_values = new_field_label_i->num_values; 1560 } 1561 } 1562 } 1563 } 1564 if (!label_found) { 1565 // LCOV_EXCL_START 1566 CeedCall(CeedFree(&new_field_label->sub_labels)); 1567 CeedCall(CeedFree(&new_field_label)); 1568 *field_label = NULL; 1569 // LCOV_EXCL_STOP 1570 } else { 1571 *field_label = new_field_label; 1572 } 1573 } else { 1574 CeedCall(CeedQFunctionContextGetFieldLabel(op->qf->ctx, field_name, field_label)); 1575 } 1576 1577 // Set label in operator 1578 if (*field_label) { 1579 (*field_label)->from_op = true; 1580 1581 // Move new composite label to operator 1582 if (op->num_context_labels == 0) { 1583 CeedCall(CeedCalloc(1, &op->context_labels)); 1584 op->max_context_labels = 1; 1585 } else if (op->num_context_labels == op->max_context_labels) { 1586 CeedCall(CeedRealloc(2 * op->num_context_labels, &op->context_labels)); 1587 op->max_context_labels *= 2; 1588 } 1589 op->context_labels[op->num_context_labels] = *field_label; 1590 op->num_context_labels++; 1591 } 1592 1593 return CEED_ERROR_SUCCESS; 1594 } 1595 1596 /** 1597 @brief Set QFunctionContext field holding double precision values. 1598 For composite operators, the values are set in all sub-operator QFunctionContexts that have a matching `field_name`. 1599 1600 @param[in,out] op CeedOperator 1601 @param[in] field_label Label of field to set 1602 @param[in] values Values to set 1603 1604 @return An error code: 0 - success, otherwise - failure 1605 1606 @ref User 1607 **/ 1608 int CeedOperatorSetContextDouble(CeedOperator op, CeedContextFieldLabel field_label, double *values) { 1609 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, values); 1610 } 1611 1612 /** 1613 @brief Get QFunctionContext field holding double precision values, read-only. 1614 For composite operators, the values correspond to the first sub-operator QFunctionContexts that has a matching `field_name`. 1615 1616 @param[in] op CeedOperator 1617 @param[in] field_label Label of field to get 1618 @param[out] num_values Number of values in the field label 1619 @param[out] values Pointer to context values 1620 1621 @return An error code: 0 - success, otherwise - failure 1622 1623 @ref User 1624 **/ 1625 int CeedOperatorGetContextDoubleRead(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const double **values) { 1626 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, num_values, values); 1627 } 1628 1629 /** 1630 @brief Restore QFunctionContext field holding double precision values, read-only. 1631 1632 @param[in] op CeedOperator 1633 @param[in] field_label Label of field to restore 1634 @param[out] values Pointer to context values 1635 1636 @return An error code: 0 - success, otherwise - failure 1637 1638 @ref User 1639 **/ 1640 int CeedOperatorRestoreContextDoubleRead(CeedOperator op, CeedContextFieldLabel field_label, const double **values) { 1641 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, values); 1642 } 1643 1644 /** 1645 @brief Set QFunctionContext field holding int32 values. 1646 For composite operators, the values are set in all sub-operator QFunctionContexts that have a matching `field_name`. 1647 1648 @param[in,out] op CeedOperator 1649 @param[in] field_label Label of field to set 1650 @param[in] values Values to set 1651 1652 @return An error code: 0 - success, otherwise - failure 1653 1654 @ref User 1655 **/ 1656 int CeedOperatorSetContextInt32(CeedOperator op, CeedContextFieldLabel field_label, int *values) { 1657 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_INT32, values); 1658 } 1659 1660 /** 1661 @brief Get QFunctionContext field holding int32 values, read-only. 1662 For composite operators, the values correspond to the first sub-operator QFunctionContexts that has a matching `field_name`. 1663 1664 @param[in] op CeedOperator 1665 @param[in] field_label Label of field to get 1666 @param[out] num_values Number of int32 values in `values` 1667 @param[out] values Pointer to context values 1668 1669 @return An error code: 0 - success, otherwise - failure 1670 1671 @ref User 1672 **/ 1673 int CeedOperatorGetContextInt32Read(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const int **values) { 1674 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_INT32, num_values, values); 1675 } 1676 1677 /** 1678 @brief Restore QFunctionContext field holding int32 values, read-only. 1679 1680 @param[in] op CeedOperator 1681 @param[in] field_label Label of field to get 1682 @param[out] values Pointer to context values 1683 1684 @return An error code: 0 - success, otherwise - failure 1685 1686 @ref User 1687 **/ 1688 int CeedOperatorRestoreContextInt32Read(CeedOperator op, CeedContextFieldLabel field_label, const int **values) { 1689 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_INT32, values); 1690 } 1691 1692 /** 1693 @brief Apply CeedOperator to a vector 1694 1695 This computes the action of the operator on the specified (active) input, yielding its (active) output. 1696 All inputs and outputs must be specified using CeedOperatorSetField(). 1697 1698 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 1699 1700 @param[in] op CeedOperator to apply 1701 @param[in] in CeedVector containing input state or @ref CEED_VECTOR_NONE if there are no active inputs 1702 @param[out] out CeedVector to store result of applying operator (must be distinct from @a in) or @ref CEED_VECTOR_NONE if there are no 1703 active outputs 1704 @param[in] request Address of CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 1705 1706 @return An error code: 0 - success, otherwise - failure 1707 1708 @ref User 1709 **/ 1710 int CeedOperatorApply(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) { 1711 CeedCall(CeedOperatorCheckReady(op)); 1712 1713 if (op->num_elem) { 1714 // Standard Operator 1715 if (op->Apply) { 1716 CeedCall(op->Apply(op, in, out, request)); 1717 } else { 1718 // Zero all output vectors 1719 CeedQFunction qf = op->qf; 1720 for (CeedInt i = 0; i < qf->num_output_fields; i++) { 1721 CeedVector vec = op->output_fields[i]->vec; 1722 if (vec == CEED_VECTOR_ACTIVE) vec = out; 1723 if (vec != CEED_VECTOR_NONE) { 1724 CeedCall(CeedVectorSetValue(vec, 0.0)); 1725 } 1726 } 1727 // Apply 1728 CeedCall(op->ApplyAdd(op, in, out, request)); 1729 } 1730 } else if (op->is_composite) { 1731 // Composite Operator 1732 if (op->ApplyComposite) { 1733 CeedCall(op->ApplyComposite(op, in, out, request)); 1734 } else { 1735 CeedInt num_suboperators; 1736 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1737 CeedOperator *sub_operators; 1738 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1739 1740 // Zero all output vectors 1741 if (out != CEED_VECTOR_NONE) { 1742 CeedCall(CeedVectorSetValue(out, 0.0)); 1743 } 1744 for (CeedInt i = 0; i < num_suboperators; i++) { 1745 for (CeedInt j = 0; j < sub_operators[i]->qf->num_output_fields; j++) { 1746 CeedVector vec = sub_operators[i]->output_fields[j]->vec; 1747 if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) { 1748 CeedCall(CeedVectorSetValue(vec, 0.0)); 1749 } 1750 } 1751 } 1752 // Apply 1753 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1754 CeedCall(CeedOperatorApplyAdd(op->sub_operators[i], in, out, request)); 1755 } 1756 } 1757 } 1758 return CEED_ERROR_SUCCESS; 1759 } 1760 1761 /** 1762 @brief Apply CeedOperator to a vector and add result to output vector 1763 1764 This computes the action of the operator on the specified (active) input, yielding its (active) output. 1765 All inputs and outputs must be specified using CeedOperatorSetField(). 1766 1767 @param[in] op CeedOperator to apply 1768 @param[in] in CeedVector containing input state or NULL if there are no active inputs 1769 @param[out] out CeedVector to sum in result of applying operator (must be distinct from @a in) or NULL if there are no active outputs 1770 @param[in] request Address of CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 1771 1772 @return An error code: 0 - success, otherwise - failure 1773 1774 @ref User 1775 **/ 1776 int CeedOperatorApplyAdd(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) { 1777 CeedCall(CeedOperatorCheckReady(op)); 1778 1779 if (op->num_elem) { 1780 // Standard Operator 1781 CeedCall(op->ApplyAdd(op, in, out, request)); 1782 } else if (op->is_composite) { 1783 // Composite Operator 1784 if (op->ApplyAddComposite) { 1785 CeedCall(op->ApplyAddComposite(op, in, out, request)); 1786 } else { 1787 CeedInt num_suboperators; 1788 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1789 CeedOperator *sub_operators; 1790 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1791 1792 for (CeedInt i = 0; i < num_suboperators; i++) { 1793 CeedCall(CeedOperatorApplyAdd(sub_operators[i], in, out, request)); 1794 } 1795 } 1796 } 1797 return CEED_ERROR_SUCCESS; 1798 } 1799 1800 /** 1801 @brief Destroy a CeedOperator 1802 1803 @param[in,out] op CeedOperator to destroy 1804 1805 @return An error code: 0 - success, otherwise - failure 1806 1807 @ref User 1808 **/ 1809 int CeedOperatorDestroy(CeedOperator *op) { 1810 if (!*op || --(*op)->ref_count > 0) { 1811 *op = NULL; 1812 return CEED_ERROR_SUCCESS; 1813 } 1814 if ((*op)->Destroy) CeedCall((*op)->Destroy(*op)); 1815 CeedCall(CeedDestroy(&(*op)->ceed)); 1816 // Free fields 1817 for (CeedInt i = 0; i < (*op)->num_fields; i++) { 1818 if ((*op)->input_fields[i]) { 1819 if ((*op)->input_fields[i]->elem_rstr != CEED_ELEMRESTRICTION_NONE) { 1820 CeedCall(CeedElemRestrictionDestroy(&(*op)->input_fields[i]->elem_rstr)); 1821 } 1822 if ((*op)->input_fields[i]->basis != CEED_BASIS_COLLOCATED) { 1823 CeedCall(CeedBasisDestroy(&(*op)->input_fields[i]->basis)); 1824 } 1825 if ((*op)->input_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->input_fields[i]->vec != CEED_VECTOR_NONE) { 1826 CeedCall(CeedVectorDestroy(&(*op)->input_fields[i]->vec)); 1827 } 1828 CeedCall(CeedFree(&(*op)->input_fields[i]->field_name)); 1829 CeedCall(CeedFree(&(*op)->input_fields[i])); 1830 } 1831 } 1832 for (CeedInt i = 0; i < (*op)->num_fields; i++) { 1833 if ((*op)->output_fields[i]) { 1834 CeedCall(CeedElemRestrictionDestroy(&(*op)->output_fields[i]->elem_rstr)); 1835 if ((*op)->output_fields[i]->basis != CEED_BASIS_COLLOCATED) { 1836 CeedCall(CeedBasisDestroy(&(*op)->output_fields[i]->basis)); 1837 } 1838 if ((*op)->output_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->output_fields[i]->vec != CEED_VECTOR_NONE) { 1839 CeedCall(CeedVectorDestroy(&(*op)->output_fields[i]->vec)); 1840 } 1841 CeedCall(CeedFree(&(*op)->output_fields[i]->field_name)); 1842 CeedCall(CeedFree(&(*op)->output_fields[i])); 1843 } 1844 } 1845 // Destroy sub_operators 1846 for (CeedInt i = 0; i < (*op)->num_suboperators; i++) { 1847 if ((*op)->sub_operators[i]) { 1848 CeedCall(CeedOperatorDestroy(&(*op)->sub_operators[i])); 1849 } 1850 } 1851 CeedCall(CeedQFunctionDestroy(&(*op)->qf)); 1852 CeedCall(CeedQFunctionDestroy(&(*op)->dqf)); 1853 CeedCall(CeedQFunctionDestroy(&(*op)->dqfT)); 1854 // Destroy any composite labels 1855 if ((*op)->is_composite) { 1856 for (CeedInt i = 0; i < (*op)->num_context_labels; i++) { 1857 CeedCall(CeedFree(&(*op)->context_labels[i]->sub_labels)); 1858 CeedCall(CeedFree(&(*op)->context_labels[i])); 1859 } 1860 } 1861 CeedCall(CeedFree(&(*op)->context_labels)); 1862 1863 // Destroy fallback 1864 CeedCall(CeedOperatorDestroy(&(*op)->op_fallback)); 1865 1866 // Destroy assembly data 1867 CeedCall(CeedQFunctionAssemblyDataDestroy(&(*op)->qf_assembled)); 1868 CeedCall(CeedOperatorAssemblyDataDestroy(&(*op)->op_assembled)); 1869 1870 CeedCall(CeedFree(&(*op)->input_fields)); 1871 CeedCall(CeedFree(&(*op)->output_fields)); 1872 CeedCall(CeedFree(&(*op)->sub_operators)); 1873 CeedCall(CeedFree(&(*op)->name)); 1874 CeedCall(CeedFree(op)); 1875 return CEED_ERROR_SUCCESS; 1876 } 1877 1878 /// @} 1879