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