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 254 For composite operators, the value is set in all sub-operator QFunctionContexts that have a matching `field_name`. 255 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 not have 256 any field of a matching type. 257 258 @param[in,out] op CeedOperator 259 @param[in] field_label Label of field to set 260 @param[in] field_type Type of field to set 261 @param[in] values Values to set 262 263 @return An error code: 0 - success, otherwise - failure 264 265 @ref User 266 **/ 267 static int CeedOperatorContextSetGeneric(CeedOperator op, CeedContextFieldLabel field_label, CeedContextFieldType field_type, void *values) { 268 if (!field_label) { 269 // LCOV_EXCL_START 270 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "Invalid field label"); 271 // LCOV_EXCL_STOP 272 } 273 274 // Check if field_label and op correspond 275 if (field_label->from_op) { 276 CeedInt index = -1; 277 278 for (CeedInt i = 0; i < op->num_context_labels; i++) { 279 if (op->context_labels[i] == field_label) index = i; 280 } 281 if (index == -1) { 282 // LCOV_EXCL_START 283 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "ContextFieldLabel does not correspond to the operator"); 284 // LCOV_EXCL_STOP 285 } 286 } 287 288 bool is_composite = false; 289 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 290 if (is_composite) { 291 CeedInt num_sub; 292 CeedOperator *sub_operators; 293 294 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub)); 295 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 296 if (num_sub != field_label->num_sub_labels) { 297 // LCOV_EXCL_START 298 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "Composite operator modified after ContextFieldLabel created"); 299 // LCOV_EXCL_STOP 300 } 301 302 for (CeedInt i = 0; i < num_sub; i++) { 303 // Try every sub-operator, ok if some sub-operators do not have field 304 if (field_label->sub_labels[i] && sub_operators[i]->qf->ctx) { 305 CeedCall(CeedQFunctionContextSetGeneric(sub_operators[i]->qf->ctx, field_label->sub_labels[i], field_type, values)); 306 } 307 } 308 } else { 309 if (!op->qf->ctx) { 310 // LCOV_EXCL_START 311 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "QFunction does not have context data"); 312 // LCOV_EXCL_STOP 313 } 314 315 CeedCall(CeedQFunctionContextSetGeneric(op->qf->ctx, field_label, field_type, values)); 316 } 317 318 return CEED_ERROR_SUCCESS; 319 } 320 321 /** 322 @brief Get QFunctionContext field values of the specified type, read-only. 323 324 For composite operators, the values retrieved are for the first sub-operator QFunctionContext that have a matching `field_name`. 325 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 not have 326 any field of a matching type. 327 328 @param[in,out] op CeedOperator 329 @param[in] field_label Label of field to set 330 @param[in] field_type Type of field to set 331 @param[out] num_values Number of values of type `field_type` in array `values` 332 @param[out] values Values in the label 333 334 @return An error code: 0 - success, otherwise - failure 335 336 @ref User 337 **/ 338 static int CeedOperatorContextGetGenericRead(CeedOperator op, CeedContextFieldLabel field_label, CeedContextFieldType field_type, size_t *num_values, 339 void *values) { 340 if (!field_label) { 341 // LCOV_EXCL_START 342 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "Invalid field label"); 343 // LCOV_EXCL_STOP 344 } 345 346 *(void **)values = NULL; 347 *num_values = 0; 348 349 // Check if field_label and op correspond 350 if (field_label->from_op) { 351 CeedInt index = -1; 352 353 for (CeedInt i = 0; i < op->num_context_labels; i++) { 354 if (op->context_labels[i] == field_label) index = i; 355 } 356 if (index == -1) { 357 // LCOV_EXCL_START 358 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "ContextFieldLabel does not correspond to the operator"); 359 // LCOV_EXCL_STOP 360 } 361 } 362 363 bool is_composite = false; 364 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 365 if (is_composite) { 366 CeedInt num_sub; 367 CeedOperator *sub_operators; 368 369 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub)); 370 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 371 if (num_sub != field_label->num_sub_labels) { 372 // LCOV_EXCL_START 373 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "Composite operator modified after ContextFieldLabel created"); 374 // LCOV_EXCL_STOP 375 } 376 377 for (CeedInt i = 0; i < num_sub; i++) { 378 // Try every sub-operator, ok if some sub-operators do not have field 379 if (field_label->sub_labels[i] && sub_operators[i]->qf->ctx) { 380 CeedCall(CeedQFunctionContextGetGenericRead(sub_operators[i]->qf->ctx, field_label->sub_labels[i], field_type, num_values, values)); 381 return CEED_ERROR_SUCCESS; 382 } 383 } 384 } else { 385 if (!op->qf->ctx) { 386 // LCOV_EXCL_START 387 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "QFunction does not have context data"); 388 // LCOV_EXCL_STOP 389 } 390 391 CeedCall(CeedQFunctionContextGetGenericRead(op->qf->ctx, field_label, field_type, num_values, values)); 392 } 393 394 return CEED_ERROR_SUCCESS; 395 } 396 397 /** 398 @brief Restore QFunctionContext field values of the specified type, read-only. 399 400 For composite operators, the values restored are for the first sub-operator QFunctionContext that have a matching `field_name`. 401 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 not have 402 any field of a matching type. 403 404 @param[in,out] op CeedOperator 405 @param[in] field_label Label of field to set 406 @param[in] field_type Type of field to set 407 @param[in] values Values array to restore 408 409 @return An error code: 0 - success, otherwise - failure 410 411 @ref User 412 **/ 413 static int CeedOperatorContextRestoreGenericRead(CeedOperator op, CeedContextFieldLabel field_label, CeedContextFieldType field_type, void *values) { 414 if (!field_label) { 415 // LCOV_EXCL_START 416 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "Invalid field label"); 417 // LCOV_EXCL_STOP 418 } 419 420 // Check if field_label and op correspond 421 if (field_label->from_op) { 422 CeedInt index = -1; 423 424 for (CeedInt i = 0; i < op->num_context_labels; i++) { 425 if (op->context_labels[i] == field_label) index = i; 426 } 427 if (index == -1) { 428 // LCOV_EXCL_START 429 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "ContextFieldLabel does not correspond to the operator"); 430 // LCOV_EXCL_STOP 431 } 432 } 433 434 bool is_composite = false; 435 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 436 if (is_composite) { 437 CeedInt num_sub; 438 CeedOperator *sub_operators; 439 440 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub)); 441 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 442 if (num_sub != field_label->num_sub_labels) { 443 // LCOV_EXCL_START 444 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "Composite operator modified after ContextFieldLabel created"); 445 // LCOV_EXCL_STOP 446 } 447 448 for (CeedInt i = 0; i < num_sub; i++) { 449 // Try every sub-operator, ok if some sub-operators do not have field 450 if (field_label->sub_labels[i] && sub_operators[i]->qf->ctx) { 451 CeedCall(CeedQFunctionContextRestoreGenericRead(sub_operators[i]->qf->ctx, field_label->sub_labels[i], field_type, values)); 452 return CEED_ERROR_SUCCESS; 453 } 454 } 455 } else { 456 if (!op->qf->ctx) { 457 // LCOV_EXCL_START 458 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "QFunction does not have context data"); 459 // LCOV_EXCL_STOP 460 } 461 462 CeedCall(CeedQFunctionContextRestoreGenericRead(op->qf->ctx, field_label, field_type, values)); 463 } 464 465 return CEED_ERROR_SUCCESS; 466 } 467 468 /// @} 469 470 /// ---------------------------------------------------------------------------- 471 /// CeedOperator Backend API 472 /// ---------------------------------------------------------------------------- 473 /// @addtogroup CeedOperatorBackend 474 /// @{ 475 476 /** 477 @brief Get the number of arguments associated with a CeedOperator 478 479 @param[in] op CeedOperator 480 @param[out] num_args Variable to store vector number of arguments 481 482 @return An error code: 0 - success, otherwise - failure 483 484 @ref Backend 485 **/ 486 int CeedOperatorGetNumArgs(CeedOperator op, CeedInt *num_args) { 487 if (op->is_composite) { 488 // LCOV_EXCL_START 489 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not defined for composite operators"); 490 // LCOV_EXCL_STOP 491 } 492 *num_args = op->num_fields; 493 return CEED_ERROR_SUCCESS; 494 } 495 496 /** 497 @brief Get the setup status of a CeedOperator 498 499 @param[in] op CeedOperator 500 @param[out] is_setup_done Variable to store setup status 501 502 @return An error code: 0 - success, otherwise - failure 503 504 @ref Backend 505 **/ 506 int CeedOperatorIsSetupDone(CeedOperator op, bool *is_setup_done) { 507 *is_setup_done = op->is_backend_setup; 508 return CEED_ERROR_SUCCESS; 509 } 510 511 /** 512 @brief Get the QFunction associated with a CeedOperator 513 514 @param[in] op CeedOperator 515 @param[out] qf Variable to store QFunction 516 517 @return An error code: 0 - success, otherwise - failure 518 519 @ref Backend 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 int CeedOperatorIsComposite(CeedOperator op, bool *is_composite) { 542 *is_composite = op->is_composite; 543 return CEED_ERROR_SUCCESS; 544 } 545 546 /** 547 @brief Get the backend data of a CeedOperator 548 549 @param[in] op CeedOperator 550 @param[out] data Variable to store data 551 552 @return An error code: 0 - success, otherwise - failure 553 554 @ref Backend 555 **/ 556 int CeedOperatorGetData(CeedOperator op, void *data) { 557 *(void **)data = op->data; 558 return CEED_ERROR_SUCCESS; 559 } 560 561 /** 562 @brief Set the backend data of a CeedOperator 563 564 @param[in,out] op CeedOperator 565 @param[in] data Data to set 566 567 @return An error code: 0 - success, otherwise - failure 568 569 @ref Backend 570 **/ 571 int CeedOperatorSetData(CeedOperator op, void *data) { 572 op->data = data; 573 return CEED_ERROR_SUCCESS; 574 } 575 576 /** 577 @brief Increment the reference counter for a CeedOperator 578 579 @param[in,out] op CeedOperator to increment the reference counter 580 581 @return An error code: 0 - success, otherwise - failure 582 583 @ref Backend 584 **/ 585 int CeedOperatorReference(CeedOperator op) { 586 op->ref_count++; 587 return CEED_ERROR_SUCCESS; 588 } 589 590 /** 591 @brief Set the setup flag of a CeedOperator to True 592 593 @param[in,out] op CeedOperator 594 595 @return An error code: 0 - success, otherwise - failure 596 597 @ref Backend 598 **/ 599 int CeedOperatorSetSetupDone(CeedOperator op) { 600 op->is_backend_setup = true; 601 return CEED_ERROR_SUCCESS; 602 } 603 604 /// @} 605 606 /// ---------------------------------------------------------------------------- 607 /// CeedOperator Public API 608 /// ---------------------------------------------------------------------------- 609 /// @addtogroup CeedOperatorUser 610 /// @{ 611 612 /** 613 @brief Create a CeedOperator and associate a CeedQFunction. 614 615 A CeedBasis and CeedElemRestriction can be associated with CeedQFunction fields with \ref CeedOperatorSetField. 616 617 @param[in] ceed Ceed object where the CeedOperator will be created 618 @param[in] qf QFunction defining the action of the operator at quadrature points 619 @param[in] dqf QFunction defining the action of the Jacobian of @a qf (or @ref CEED_QFUNCTION_NONE) 620 @param[in] dqfT QFunction defining the action of the transpose of the Jacobian of @a qf (or @ref CEED_QFUNCTION_NONE) 621 @param[out] op Address of the variable where the newly created CeedOperator will be stored 622 623 @return An error code: 0 - success, otherwise - failure 624 625 @ref User 626 */ 627 int CeedOperatorCreate(Ceed ceed, CeedQFunction qf, CeedQFunction dqf, CeedQFunction dqfT, CeedOperator *op) { 628 if (!ceed->OperatorCreate) { 629 Ceed delegate; 630 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Operator")); 631 632 if (!delegate) { 633 // LCOV_EXCL_START 634 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support OperatorCreate"); 635 // LCOV_EXCL_STOP 636 } 637 638 CeedCall(CeedOperatorCreate(delegate, qf, dqf, dqfT, op)); 639 return CEED_ERROR_SUCCESS; 640 } 641 642 if (!qf || qf == CEED_QFUNCTION_NONE) { 643 // LCOV_EXCL_START 644 return CeedError(ceed, CEED_ERROR_MINOR, "Operator must have a valid QFunction."); 645 // LCOV_EXCL_STOP 646 } 647 CeedCall(CeedCalloc(1, op)); 648 (*op)->ceed = ceed; 649 CeedCall(CeedReference(ceed)); 650 (*op)->ref_count = 1; 651 (*op)->qf = qf; 652 (*op)->input_size = -1; 653 (*op)->output_size = -1; 654 CeedCall(CeedQFunctionReference(qf)); 655 if (dqf && dqf != CEED_QFUNCTION_NONE) { 656 (*op)->dqf = dqf; 657 CeedCall(CeedQFunctionReference(dqf)); 658 } 659 if (dqfT && dqfT != CEED_QFUNCTION_NONE) { 660 (*op)->dqfT = dqfT; 661 CeedCall(CeedQFunctionReference(dqfT)); 662 } 663 CeedCall(CeedQFunctionAssemblyDataCreate(ceed, &(*op)->qf_assembled)); 664 CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*op)->input_fields)); 665 CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*op)->output_fields)); 666 CeedCall(ceed->OperatorCreate(*op)); 667 return CEED_ERROR_SUCCESS; 668 } 669 670 /** 671 @brief Create an operator that composes the action of several operators 672 673 @param[in] ceed Ceed object where the CeedOperator will be created 674 @param[out] op Address of the variable where the newly created Composite CeedOperator will be stored 675 676 @return An error code: 0 - success, otherwise - failure 677 678 @ref User 679 */ 680 int CeedCompositeOperatorCreate(Ceed ceed, CeedOperator *op) { 681 if (!ceed->CompositeOperatorCreate) { 682 Ceed delegate; 683 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Operator")); 684 685 if (delegate) { 686 CeedCall(CeedCompositeOperatorCreate(delegate, op)); 687 return CEED_ERROR_SUCCESS; 688 } 689 } 690 691 CeedCall(CeedCalloc(1, op)); 692 (*op)->ceed = ceed; 693 CeedCall(CeedReference(ceed)); 694 (*op)->ref_count = 1; 695 (*op)->is_composite = true; 696 CeedCall(CeedCalloc(CEED_COMPOSITE_MAX, &(*op)->sub_operators)); 697 (*op)->input_size = -1; 698 (*op)->output_size = -1; 699 700 if (ceed->CompositeOperatorCreate) { 701 CeedCall(ceed->CompositeOperatorCreate(*op)); 702 } 703 return CEED_ERROR_SUCCESS; 704 } 705 706 /** 707 @brief Copy the pointer to a CeedOperator. 708 709 Both pointers should be destroyed with `CeedOperatorDestroy()`. 710 711 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. 712 This CeedOperator will be destroyed if `op_copy` is the only reference to this CeedOperator. 713 714 @param[in] op CeedOperator to copy reference to 715 @param[in,out] op_copy Variable to store copied reference 716 717 @return An error code: 0 - success, otherwise - failure 718 719 @ref User 720 **/ 721 int CeedOperatorReferenceCopy(CeedOperator op, CeedOperator *op_copy) { 722 CeedCall(CeedOperatorReference(op)); 723 CeedCall(CeedOperatorDestroy(op_copy)); 724 *op_copy = op; 725 return CEED_ERROR_SUCCESS; 726 } 727 728 /** 729 @brief Provide a field to a CeedOperator for use by its CeedQFunction. 730 731 This function is used to specify both active and passive fields to a CeedOperator. 732 For passive fields, a vector @arg v must be provided. 733 Passive fields can inputs or outputs (updated in-place when operator is applied). 734 735 Active fields must be specified using this function, but their data (in a CeedVector) is passed in CeedOperatorApply(). 736 There can be at most one active input CeedVector and at most one active output CeedVector passed to CeedOperatorApply(). 737 738 @param[in,out] op CeedOperator on which to provide the field 739 @param[in] field_name Name of the field (to be matched with the name used by CeedQFunction) 740 @param[in] r CeedElemRestriction 741 @param[in] b CeedBasis in which the field resides or @ref CEED_BASIS_COLLOCATED if collocated with quadrature points 742 @param[in] v CeedVector to be used by CeedOperator or @ref CEED_VECTOR_ACTIVE if field is active or @ref CEED_VECTOR_NONE 743 if using @ref CEED_EVAL_WEIGHT in the QFunction 744 745 @return An error code: 0 - success, otherwise - failure 746 747 @ref User 748 **/ 749 int CeedOperatorSetField(CeedOperator op, const char *field_name, CeedElemRestriction r, CeedBasis b, CeedVector v) { 750 if (op->is_composite) { 751 // LCOV_EXCL_START 752 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "Cannot add field to composite operator."); 753 // LCOV_EXCL_STOP 754 } 755 if (op->is_immutable) { 756 // LCOV_EXCL_START 757 return CeedError(op->ceed, CEED_ERROR_MAJOR, "Operator cannot be changed after set as immutable"); 758 // LCOV_EXCL_STOP 759 } 760 if (!r) { 761 // LCOV_EXCL_START 762 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "ElemRestriction r for field \"%s\" must be non-NULL.", field_name); 763 // LCOV_EXCL_STOP 764 } 765 if (!b) { 766 // LCOV_EXCL_START 767 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "Basis b for field \"%s\" must be non-NULL.", field_name); 768 // LCOV_EXCL_STOP 769 } 770 if (!v) { 771 // LCOV_EXCL_START 772 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "Vector v for field \"%s\" must be non-NULL.", field_name); 773 // LCOV_EXCL_STOP 774 } 775 776 CeedInt num_elem; 777 CeedCall(CeedElemRestrictionGetNumElements(r, &num_elem)); 778 if (r != CEED_ELEMRESTRICTION_NONE && op->has_restriction && op->num_elem != num_elem) { 779 // LCOV_EXCL_START 780 return CeedError(op->ceed, CEED_ERROR_DIMENSION, 781 "ElemRestriction with %" CeedInt_FMT " elements incompatible with prior %" CeedInt_FMT " elements", num_elem, op->num_elem); 782 // LCOV_EXCL_STOP 783 } 784 785 CeedInt num_qpts = 0; 786 if (b != CEED_BASIS_COLLOCATED) { 787 CeedCall(CeedBasisGetNumQuadraturePoints(b, &num_qpts)); 788 if (op->num_qpts && op->num_qpts != num_qpts) { 789 // LCOV_EXCL_START 790 return CeedError(op->ceed, CEED_ERROR_DIMENSION, 791 "Basis with %" CeedInt_FMT " quadrature points incompatible with prior %" CeedInt_FMT " points", num_qpts, op->num_qpts); 792 // LCOV_EXCL_STOP 793 } 794 } 795 CeedQFunctionField qf_field; 796 CeedOperatorField *op_field; 797 bool is_input = true; 798 for (CeedInt i = 0; i < op->qf->num_input_fields; i++) { 799 if (!strcmp(field_name, (*op->qf->input_fields[i]).field_name)) { 800 qf_field = op->qf->input_fields[i]; 801 op_field = &op->input_fields[i]; 802 goto found; 803 } 804 } 805 is_input = false; 806 for (CeedInt i = 0; i < op->qf->num_output_fields; i++) { 807 if (!strcmp(field_name, (*op->qf->output_fields[i]).field_name)) { 808 qf_field = op->qf->output_fields[i]; 809 op_field = &op->output_fields[i]; 810 goto found; 811 } 812 } 813 // LCOV_EXCL_START 814 return CeedError(op->ceed, CEED_ERROR_INCOMPLETE, "QFunction has no knowledge of field '%s'", field_name); 815 // LCOV_EXCL_STOP 816 found: 817 CeedCall(CeedOperatorCheckField(op->ceed, qf_field, r, b)); 818 CeedCall(CeedCalloc(1, op_field)); 819 820 if (v == CEED_VECTOR_ACTIVE) { 821 CeedSize l_size; 822 CeedCall(CeedElemRestrictionGetLVectorSize(r, &l_size)); 823 if (is_input) { 824 if (op->input_size == -1) op->input_size = l_size; 825 if (l_size != op->input_size) { 826 // LCOV_EXCL_START 827 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "LVector size %td does not match previous size %td", l_size, op->input_size); 828 // LCOV_EXCL_STOP 829 } 830 } else { 831 if (op->output_size == -1) op->output_size = l_size; 832 if (l_size != op->output_size) { 833 // LCOV_EXCL_START 834 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "LVector size %td does not match previous size %td", l_size, op->output_size); 835 // LCOV_EXCL_STOP 836 } 837 } 838 } 839 840 CeedCall(CeedVectorReferenceCopy(v, &(*op_field)->vec)); 841 CeedCall(CeedElemRestrictionReferenceCopy(r, &(*op_field)->elem_rstr)); 842 if (r != CEED_ELEMRESTRICTION_NONE) { 843 op->num_elem = num_elem; 844 op->has_restriction = true; // Restriction set, but num_elem may be 0 845 } 846 CeedCall(CeedBasisReferenceCopy(b, &(*op_field)->basis)); 847 if (!op->num_qpts && b != CEED_BASIS_COLLOCATED) CeedCall(CeedOperatorSetNumQuadraturePoints(op, num_qpts)); 848 849 op->num_fields += 1; 850 CeedCall(CeedStringAllocCopy(field_name, (char **)&(*op_field)->field_name)); 851 return CEED_ERROR_SUCCESS; 852 } 853 854 /** 855 @brief Get the CeedOperatorFields of a CeedOperator 856 857 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 858 859 @param[in] op CeedOperator 860 @param[out] num_input_fields Variable to store number of input fields 861 @param[out] input_fields Variable to store input_fields 862 @param[out] num_output_fields Variable to store number of output fields 863 @param[out] output_fields Variable to store output_fields 864 865 @return An error code: 0 - success, otherwise - failure 866 867 @ref Advanced 868 **/ 869 int CeedOperatorGetFields(CeedOperator op, CeedInt *num_input_fields, CeedOperatorField **input_fields, CeedInt *num_output_fields, 870 CeedOperatorField **output_fields) { 871 if (op->is_composite) { 872 // LCOV_EXCL_START 873 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not defined for composite operator"); 874 // LCOV_EXCL_STOP 875 } 876 CeedCall(CeedOperatorCheckReady(op)); 877 878 if (num_input_fields) *num_input_fields = op->qf->num_input_fields; 879 if (input_fields) *input_fields = op->input_fields; 880 if (num_output_fields) *num_output_fields = op->qf->num_output_fields; 881 if (output_fields) *output_fields = op->output_fields; 882 return CEED_ERROR_SUCCESS; 883 } 884 885 /** 886 @brief Get a CeedOperatorField of an CeedOperator from its name 887 888 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 889 890 @param[in] op CeedOperator 891 @param[in] field_name Name of desired CeedOperatorField 892 @param[out] op_field CeedOperatorField corresponding to the name 893 894 @return An error code: 0 - success, otherwise - failure 895 896 @ref Advanced 897 **/ 898 int CeedOperatorGetFieldByName(CeedOperator op, const char *field_name, CeedOperatorField *op_field) { 899 CeedInt num_input_fields, num_output_fields; 900 CeedOperatorField *input_fields, *output_fields; 901 char *name; 902 CeedCall(CeedOperatorGetFields(op, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); 903 904 for (CeedInt i = 0; i < num_input_fields; i++) { 905 CeedCall(CeedOperatorFieldGetName(input_fields[i], &name)); 906 if (!strcmp(name, field_name)) { 907 *op_field = input_fields[i]; 908 return CEED_ERROR_SUCCESS; 909 } 910 } 911 912 for (CeedInt i = 0; i < num_output_fields; i++) { 913 CeedCall(CeedOperatorFieldGetName(output_fields[i], &name)); 914 if (!strcmp(name, field_name)) { 915 *op_field = output_fields[i]; 916 return CEED_ERROR_SUCCESS; 917 } 918 } 919 920 bool has_name = op->name; 921 // LCOV_EXCL_START 922 return CeedError(op->ceed, CEED_ERROR_MINOR, "The field \"%s\" not found in CeedOperator%s%s%s.\n", field_name, has_name ? " \"" : "", 923 has_name ? op->name : "", has_name ? "\"" : ""); 924 // LCOV_EXCL_STOP 925 } 926 927 /** 928 @brief Get the name of a CeedOperatorField 929 930 @param[in] op_field CeedOperatorField 931 @param[out] field_name Variable to store the field name 932 933 @return An error code: 0 - success, otherwise - failure 934 935 @ref Advanced 936 **/ 937 int CeedOperatorFieldGetName(CeedOperatorField op_field, char **field_name) { 938 *field_name = (char *)op_field->field_name; 939 return CEED_ERROR_SUCCESS; 940 } 941 942 /** 943 @brief Get the CeedElemRestriction of a CeedOperatorField 944 945 @param[in] op_field CeedOperatorField 946 @param[out] rstr Variable to store CeedElemRestriction 947 948 @return An error code: 0 - success, otherwise - failure 949 950 @ref Advanced 951 **/ 952 int CeedOperatorFieldGetElemRestriction(CeedOperatorField op_field, CeedElemRestriction *rstr) { 953 *rstr = op_field->elem_rstr; 954 return CEED_ERROR_SUCCESS; 955 } 956 957 /** 958 @brief Get the CeedBasis of a CeedOperatorField 959 960 @param[in] op_field CeedOperatorField 961 @param[out] basis Variable to store CeedBasis 962 963 @return An error code: 0 - success, otherwise - failure 964 965 @ref Advanced 966 **/ 967 int CeedOperatorFieldGetBasis(CeedOperatorField op_field, CeedBasis *basis) { 968 *basis = op_field->basis; 969 return CEED_ERROR_SUCCESS; 970 } 971 972 /** 973 @brief Get the CeedVector of a CeedOperatorField 974 975 @param[in] op_field CeedOperatorField 976 @param[out] vec Variable to store CeedVector 977 978 @return An error code: 0 - success, otherwise - failure 979 980 @ref Advanced 981 **/ 982 int CeedOperatorFieldGetVector(CeedOperatorField op_field, CeedVector *vec) { 983 *vec = op_field->vec; 984 return CEED_ERROR_SUCCESS; 985 } 986 987 /** 988 @brief Add a sub-operator to a composite CeedOperator 989 990 @param[in,out] composite_op Composite CeedOperator 991 @param[in] sub_op Sub-operator CeedOperator 992 993 @return An error code: 0 - success, otherwise - failure 994 995 @ref User 996 */ 997 int CeedCompositeOperatorAddSub(CeedOperator composite_op, CeedOperator sub_op) { 998 if (!composite_op->is_composite) { 999 // LCOV_EXCL_START 1000 return CeedError(composite_op->ceed, CEED_ERROR_MINOR, "CeedOperator is not a composite operator"); 1001 // LCOV_EXCL_STOP 1002 } 1003 1004 if (composite_op->num_suboperators == CEED_COMPOSITE_MAX) { 1005 // LCOV_EXCL_START 1006 return CeedError(composite_op->ceed, CEED_ERROR_UNSUPPORTED, "Cannot add additional sub-operators"); 1007 // LCOV_EXCL_STOP 1008 } 1009 if (composite_op->is_immutable) { 1010 // LCOV_EXCL_START 1011 return CeedError(composite_op->ceed, CEED_ERROR_MAJOR, "Operator cannot be changed after set as immutable"); 1012 // LCOV_EXCL_STOP 1013 } 1014 1015 { 1016 CeedSize input_size, output_size; 1017 CeedCall(CeedOperatorGetActiveVectorLengths(sub_op, &input_size, &output_size)); 1018 if (composite_op->input_size == -1) composite_op->input_size = input_size; 1019 if (composite_op->output_size == -1) composite_op->output_size = output_size; 1020 // Note, a size of -1 means no active vector restriction set, so no incompatibility 1021 if ((input_size != -1 && input_size != composite_op->input_size) || (output_size != -1 && output_size != composite_op->output_size)) { 1022 // LCOV_EXCL_START 1023 return CeedError(composite_op->ceed, CEED_ERROR_MAJOR, 1024 "Sub-operators must have compatible dimensions; composite operator of shape (%td, %td) not compatible with sub-operator of " 1025 "shape (%td, %td)", 1026 composite_op->input_size, composite_op->output_size, input_size, output_size); 1027 // LCOV_EXCL_STOP 1028 } 1029 } 1030 1031 composite_op->sub_operators[composite_op->num_suboperators] = sub_op; 1032 CeedCall(CeedOperatorReference(sub_op)); 1033 composite_op->num_suboperators++; 1034 return CEED_ERROR_SUCCESS; 1035 } 1036 1037 /** 1038 @brief Get the number of sub_operators associated with a CeedOperator 1039 1040 @param[in] op CeedOperator 1041 @param[out] num_suboperators Variable to store number of sub_operators 1042 1043 @return An error code: 0 - success, otherwise - failure 1044 1045 @ref Backend 1046 **/ 1047 int CeedCompositeOperatorGetNumSub(CeedOperator op, CeedInt *num_suboperators) { 1048 if (!op->is_composite) { 1049 // LCOV_EXCL_START 1050 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not a composite operator"); 1051 // LCOV_EXCL_STOP 1052 } 1053 *num_suboperators = op->num_suboperators; 1054 return CEED_ERROR_SUCCESS; 1055 } 1056 1057 /** 1058 @brief Get the list of sub_operators associated with a CeedOperator 1059 1060 @param op CeedOperator 1061 @param[out] sub_operators Variable to store list of sub_operators 1062 1063 @return An error code: 0 - success, otherwise - failure 1064 1065 @ref Backend 1066 **/ 1067 int CeedCompositeOperatorGetSubList(CeedOperator op, CeedOperator **sub_operators) { 1068 if (!op->is_composite) { 1069 // LCOV_EXCL_START 1070 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not a composite operator"); 1071 // LCOV_EXCL_STOP 1072 } 1073 *sub_operators = op->sub_operators; 1074 return CEED_ERROR_SUCCESS; 1075 } 1076 1077 /** 1078 @brief Check if a CeedOperator is ready to be used. 1079 1080 @param[in] op CeedOperator to check 1081 1082 @return An error code: 0 - success, otherwise - failure 1083 1084 @ref User 1085 **/ 1086 int CeedOperatorCheckReady(CeedOperator op) { 1087 Ceed ceed; 1088 CeedCall(CeedOperatorGetCeed(op, &ceed)); 1089 1090 if (op->is_interface_setup) return CEED_ERROR_SUCCESS; 1091 1092 CeedQFunction qf = op->qf; 1093 if (op->is_composite) { 1094 if (!op->num_suboperators) { 1095 // Empty operator setup 1096 op->input_size = 0; 1097 op->output_size = 0; 1098 } else { 1099 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1100 CeedCall(CeedOperatorCheckReady(op->sub_operators[i])); 1101 } 1102 // Sub-operators could be modified after adding to composite operator 1103 // Need to verify no lvec incompatibility from any changes 1104 CeedSize input_size, output_size; 1105 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 1106 } 1107 } else { 1108 if (op->num_fields == 0) { 1109 // LCOV_EXCL_START 1110 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "No operator fields set"); 1111 // LCOV_EXCL_STOP 1112 } 1113 if (op->num_fields < qf->num_input_fields + qf->num_output_fields) { 1114 // LCOV_EXCL_START 1115 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "Not all operator fields set"); 1116 // LCOV_EXCL_STOP 1117 } 1118 if (!op->has_restriction) { 1119 // LCOV_EXCL_START 1120 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "At least one restriction required"); 1121 // LCOV_EXCL_STOP 1122 } 1123 if (op->num_qpts == 0) { 1124 // LCOV_EXCL_START 1125 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "At least one non-collocated basis is required or the number of quadrature points must be set"); 1126 // LCOV_EXCL_STOP 1127 } 1128 } 1129 1130 // Flag as immutable and ready 1131 op->is_interface_setup = true; 1132 if (op->qf && op->qf != CEED_QFUNCTION_NONE) op->qf->is_immutable = true; 1133 if (op->dqf && op->dqf != CEED_QFUNCTION_NONE) op->dqf->is_immutable = true; 1134 if (op->dqfT && op->dqfT != CEED_QFUNCTION_NONE) op->dqfT->is_immutable = true; 1135 return CEED_ERROR_SUCCESS; 1136 } 1137 1138 /** 1139 @brief Get vector lengths for the active input and/or output vectors of a CeedOperator. 1140 1141 Note: Lengths of -1 indicate that the CeedOperator does not have an active input and/or output. 1142 1143 @param[in] op CeedOperator 1144 @param[out] input_size Variable to store active input vector length, or NULL 1145 @param[out] output_size Variable to store active output vector length, or NULL 1146 1147 @return An error code: 0 - success, otherwise - failure 1148 1149 @ref User 1150 **/ 1151 int CeedOperatorGetActiveVectorLengths(CeedOperator op, CeedSize *input_size, CeedSize *output_size) { 1152 bool is_composite; 1153 1154 if (input_size) *input_size = op->input_size; 1155 if (output_size) *output_size = op->output_size; 1156 1157 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1158 if (is_composite && (op->input_size == -1 || op->output_size == -1)) { 1159 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1160 CeedSize sub_input_size, sub_output_size; 1161 CeedCall(CeedOperatorGetActiveVectorLengths(op->sub_operators[i], &sub_input_size, &sub_output_size)); 1162 if (op->input_size == -1) op->input_size = sub_input_size; 1163 if (op->output_size == -1) op->output_size = sub_output_size; 1164 // Note, a size of -1 means no active vector restriction set, so no incompatibility 1165 if ((sub_input_size != -1 && sub_input_size != op->input_size) || (sub_output_size != -1 && sub_output_size != op->output_size)) { 1166 // LCOV_EXCL_START 1167 return CeedError(op->ceed, CEED_ERROR_MAJOR, 1168 "Sub-operators must have compatible dimensions; composite operator of shape (%td, %td) not compatible with sub-operator of " 1169 "shape (%td, %td)", 1170 op->input_size, op->output_size, input_size, output_size); 1171 // LCOV_EXCL_STOP 1172 } 1173 } 1174 } 1175 1176 return CEED_ERROR_SUCCESS; 1177 } 1178 1179 /** 1180 @brief Set reuse of CeedQFunction data in CeedOperatorLinearAssemble* functions. 1181 1182 When `reuse_assembly_data = false` (default), the CeedQFunction associated with this CeedOperator is re-assembled every time a 1183 `CeedOperatorLinearAssemble*` function is called. 1184 When `reuse_assembly_data = true`, the CeedQFunction associated with this CeedOperator is reused between calls to 1185 `CeedOperatorSetQFunctionAssemblyDataUpdated`. 1186 1187 @param[in] op CeedOperator 1188 @param[in] reuse_assembly_data Boolean flag setting assembly data reuse 1189 1190 @return An error code: 0 - success, otherwise - failure 1191 1192 @ref Advanced 1193 **/ 1194 int CeedOperatorSetQFunctionAssemblyReuse(CeedOperator op, bool reuse_assembly_data) { 1195 bool is_composite; 1196 1197 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1198 if (is_composite) { 1199 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1200 CeedCall(CeedOperatorSetQFunctionAssemblyReuse(op->sub_operators[i], reuse_assembly_data)); 1201 } 1202 } else { 1203 CeedCall(CeedQFunctionAssemblyDataSetReuse(op->qf_assembled, reuse_assembly_data)); 1204 } 1205 1206 return CEED_ERROR_SUCCESS; 1207 } 1208 1209 /** 1210 @brief Mark CeedQFunction data as updated and the CeedQFunction as requiring re-assembly. 1211 1212 @param[in] op CeedOperator 1213 @param[in] needs_data_update Boolean flag setting assembly data reuse 1214 1215 @return An error code: 0 - success, otherwise - failure 1216 1217 @ref Advanced 1218 **/ 1219 int CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(CeedOperator op, bool needs_data_update) { 1220 bool is_composite; 1221 1222 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1223 if (is_composite) { 1224 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1225 CeedCall(CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(op->sub_operators[i], needs_data_update)); 1226 } 1227 } else { 1228 CeedCall(CeedQFunctionAssemblyDataSetUpdateNeeded(op->qf_assembled, needs_data_update)); 1229 } 1230 1231 return CEED_ERROR_SUCCESS; 1232 } 1233 1234 /** 1235 @brief Set the number of quadrature points associated with a CeedOperator. 1236 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 CeedCall(CeedOperatorGetFields(op, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); 1404 CeedInt num_elem = 0; 1405 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 1406 1407 // Input FLOPs 1408 for (CeedInt i = 0; i < num_input_fields; i++) { 1409 if (input_fields[i]->vec == CEED_VECTOR_ACTIVE) { 1410 CeedSize restr_flops, basis_flops; 1411 1412 CeedCall(CeedElemRestrictionGetFlopsEstimate(input_fields[i]->elem_rstr, CEED_NOTRANSPOSE, &restr_flops)); 1413 *flops += restr_flops; 1414 CeedCall(CeedBasisGetFlopsEstimate(input_fields[i]->basis, CEED_NOTRANSPOSE, op->qf->input_fields[i]->eval_mode, &basis_flops)); 1415 *flops += basis_flops * num_elem; 1416 } 1417 } 1418 // QF FLOPs 1419 CeedInt num_qpts; 1420 CeedSize qf_flops; 1421 CeedCall(CeedOperatorGetNumQuadraturePoints(op, &num_qpts)); 1422 CeedCall(CeedQFunctionGetFlopsEstimate(op->qf, &qf_flops)); 1423 *flops += num_elem * num_qpts * qf_flops; 1424 // Output FLOPs 1425 for (CeedInt i = 0; i < num_output_fields; i++) { 1426 if (output_fields[i]->vec == CEED_VECTOR_ACTIVE) { 1427 CeedSize restr_flops, basis_flops; 1428 1429 CeedCall(CeedElemRestrictionGetFlopsEstimate(output_fields[i]->elem_rstr, CEED_TRANSPOSE, &restr_flops)); 1430 *flops += restr_flops; 1431 CeedCall(CeedBasisGetFlopsEstimate(output_fields[i]->basis, CEED_TRANSPOSE, op->qf->output_fields[i]->eval_mode, &basis_flops)); 1432 *flops += basis_flops * num_elem; 1433 } 1434 } 1435 } 1436 1437 return CEED_ERROR_SUCCESS; 1438 } 1439 1440 /** 1441 @brief Get CeedQFunction global context for a CeedOperator. 1442 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 CeedQFunctionContext. 1446 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 1569 For composite operators, the values are set in all sub-operator QFunctionContexts that have a matching `field_name`. 1570 1571 @param[in,out] op CeedOperator 1572 @param[in] field_label Label of field to set 1573 @param[in] values Values to set 1574 1575 @return An error code: 0 - success, otherwise - failure 1576 1577 @ref User 1578 **/ 1579 int CeedOperatorSetContextDouble(CeedOperator op, CeedContextFieldLabel field_label, double *values) { 1580 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, values); 1581 } 1582 1583 /** 1584 @brief Get QFunctionContext field holding double precision values, read-only. 1585 1586 For composite operators, the values correspond to the first sub-operator QFunctionContexts that has a matching `field_name`. 1587 1588 @param[in] op CeedOperator 1589 @param[in] field_label Label of field to get 1590 @param[out] num_values Number of values in the field label 1591 @param[out] values Pointer to context values 1592 1593 @return An error code: 0 - success, otherwise - failure 1594 1595 @ref User 1596 **/ 1597 int CeedOperatorGetContextDoubleRead(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const double **values) { 1598 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, num_values, values); 1599 } 1600 1601 /** 1602 @brief Restore QFunctionContext field holding double precision values, read-only. 1603 1604 @param[in] op CeedOperator 1605 @param[in] field_label Label of field to restore 1606 @param[out] values Pointer to context values 1607 1608 @return An error code: 0 - success, otherwise - failure 1609 1610 @ref User 1611 **/ 1612 int CeedOperatorRestoreContextDoubleRead(CeedOperator op, CeedContextFieldLabel field_label, const double **values) { 1613 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, values); 1614 } 1615 1616 /** 1617 @brief Set QFunctionContext field holding int32 values. 1618 1619 For composite operators, the values are set in all sub-operator QFunctionContexts that have a matching `field_name`. 1620 1621 @param[in,out] op CeedOperator 1622 @param[in] field_label Label of field to set 1623 @param[in] values Values to set 1624 1625 @return An error code: 0 - success, otherwise - failure 1626 1627 @ref User 1628 **/ 1629 int CeedOperatorSetContextInt32(CeedOperator op, CeedContextFieldLabel field_label, int *values) { 1630 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_INT32, values); 1631 } 1632 1633 /** 1634 @brief Get QFunctionContext field holding int32 values, read-only. 1635 1636 For composite operators, the values correspond to the first sub-operator QFunctionContexts that has a matching `field_name`. 1637 1638 @param[in] op CeedOperator 1639 @param[in] field_label Label of field to get 1640 @param[out] num_values Number of int32 values in `values` 1641 @param[out] values Pointer to context values 1642 1643 @return An error code: 0 - success, otherwise - failure 1644 1645 @ref User 1646 **/ 1647 int CeedOperatorGetContextInt32Read(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const int **values) { 1648 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_INT32, num_values, values); 1649 } 1650 1651 /** 1652 @brief Restore QFunctionContext field holding int32 values, read-only. 1653 1654 @param[in] op CeedOperator 1655 @param[in] field_label Label of field to get 1656 @param[out] values Pointer to context values 1657 1658 @return An error code: 0 - success, otherwise - failure 1659 1660 @ref User 1661 **/ 1662 int CeedOperatorRestoreContextInt32Read(CeedOperator op, CeedContextFieldLabel field_label, const int **values) { 1663 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_INT32, values); 1664 } 1665 1666 /** 1667 @brief Apply CeedOperator to a vector 1668 1669 This computes the action of the operator on the specified (active) input, yielding its (active) output. 1670 All inputs and outputs must be specified using CeedOperatorSetField(). 1671 1672 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 1673 1674 @param[in] op CeedOperator to apply 1675 @param[in] in CeedVector containing input state or @ref CEED_VECTOR_NONE if there are no active inputs 1676 @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 1677 active outputs 1678 @param[in] request Address of CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 1679 1680 @return An error code: 0 - success, otherwise - failure 1681 1682 @ref User 1683 **/ 1684 int CeedOperatorApply(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) { 1685 CeedCall(CeedOperatorCheckReady(op)); 1686 1687 if (op->num_elem) { 1688 // Standard Operator 1689 if (op->Apply) { 1690 CeedCall(op->Apply(op, in, out, request)); 1691 } else { 1692 // Zero all output vectors 1693 CeedQFunction qf = op->qf; 1694 for (CeedInt i = 0; i < qf->num_output_fields; i++) { 1695 CeedVector vec = op->output_fields[i]->vec; 1696 if (vec == CEED_VECTOR_ACTIVE) vec = out; 1697 if (vec != CEED_VECTOR_NONE) { 1698 CeedCall(CeedVectorSetValue(vec, 0.0)); 1699 } 1700 } 1701 // Apply 1702 CeedCall(op->ApplyAdd(op, in, out, request)); 1703 } 1704 } else if (op->is_composite) { 1705 // Composite Operator 1706 if (op->ApplyComposite) { 1707 CeedCall(op->ApplyComposite(op, in, out, request)); 1708 } else { 1709 CeedInt num_suboperators; 1710 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1711 CeedOperator *sub_operators; 1712 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1713 1714 // Zero all output vectors 1715 if (out != CEED_VECTOR_NONE) { 1716 CeedCall(CeedVectorSetValue(out, 0.0)); 1717 } 1718 for (CeedInt i = 0; i < num_suboperators; i++) { 1719 for (CeedInt j = 0; j < sub_operators[i]->qf->num_output_fields; j++) { 1720 CeedVector vec = sub_operators[i]->output_fields[j]->vec; 1721 if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) { 1722 CeedCall(CeedVectorSetValue(vec, 0.0)); 1723 } 1724 } 1725 } 1726 // Apply 1727 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1728 CeedCall(CeedOperatorApplyAdd(op->sub_operators[i], in, out, request)); 1729 } 1730 } 1731 } 1732 return CEED_ERROR_SUCCESS; 1733 } 1734 1735 /** 1736 @brief Apply CeedOperator to a vector and add result to output vector 1737 1738 This computes the action of the operator on the specified (active) input, yielding its (active) output. 1739 All inputs and outputs must be specified using CeedOperatorSetField(). 1740 1741 @param[in] op CeedOperator to apply 1742 @param[in] in CeedVector containing input state or NULL if there are no active inputs 1743 @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 1744 @param[in] request Address of CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 1745 1746 @return An error code: 0 - success, otherwise - failure 1747 1748 @ref User 1749 **/ 1750 int CeedOperatorApplyAdd(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) { 1751 CeedCall(CeedOperatorCheckReady(op)); 1752 1753 if (op->num_elem) { 1754 // Standard Operator 1755 CeedCall(op->ApplyAdd(op, in, out, request)); 1756 } else if (op->is_composite) { 1757 // Composite Operator 1758 if (op->ApplyAddComposite) { 1759 CeedCall(op->ApplyAddComposite(op, in, out, request)); 1760 } else { 1761 CeedInt num_suboperators; 1762 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1763 CeedOperator *sub_operators; 1764 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1765 1766 for (CeedInt i = 0; i < num_suboperators; i++) { 1767 CeedCall(CeedOperatorApplyAdd(sub_operators[i], in, out, request)); 1768 } 1769 } 1770 } 1771 return CEED_ERROR_SUCCESS; 1772 } 1773 1774 /** 1775 @brief Destroy a CeedOperator 1776 1777 @param[in,out] op CeedOperator to destroy 1778 1779 @return An error code: 0 - success, otherwise - failure 1780 1781 @ref User 1782 **/ 1783 int CeedOperatorDestroy(CeedOperator *op) { 1784 if (!*op || --(*op)->ref_count > 0) { 1785 *op = NULL; 1786 return CEED_ERROR_SUCCESS; 1787 } 1788 if ((*op)->Destroy) CeedCall((*op)->Destroy(*op)); 1789 CeedCall(CeedDestroy(&(*op)->ceed)); 1790 // Free fields 1791 for (CeedInt i = 0; i < (*op)->num_fields; i++) { 1792 if ((*op)->input_fields[i]) { 1793 if ((*op)->input_fields[i]->elem_rstr != CEED_ELEMRESTRICTION_NONE) { 1794 CeedCall(CeedElemRestrictionDestroy(&(*op)->input_fields[i]->elem_rstr)); 1795 } 1796 if ((*op)->input_fields[i]->basis != CEED_BASIS_COLLOCATED) { 1797 CeedCall(CeedBasisDestroy(&(*op)->input_fields[i]->basis)); 1798 } 1799 if ((*op)->input_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->input_fields[i]->vec != CEED_VECTOR_NONE) { 1800 CeedCall(CeedVectorDestroy(&(*op)->input_fields[i]->vec)); 1801 } 1802 CeedCall(CeedFree(&(*op)->input_fields[i]->field_name)); 1803 CeedCall(CeedFree(&(*op)->input_fields[i])); 1804 } 1805 } 1806 for (CeedInt i = 0; i < (*op)->num_fields; i++) { 1807 if ((*op)->output_fields[i]) { 1808 CeedCall(CeedElemRestrictionDestroy(&(*op)->output_fields[i]->elem_rstr)); 1809 if ((*op)->output_fields[i]->basis != CEED_BASIS_COLLOCATED) { 1810 CeedCall(CeedBasisDestroy(&(*op)->output_fields[i]->basis)); 1811 } 1812 if ((*op)->output_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->output_fields[i]->vec != CEED_VECTOR_NONE) { 1813 CeedCall(CeedVectorDestroy(&(*op)->output_fields[i]->vec)); 1814 } 1815 CeedCall(CeedFree(&(*op)->output_fields[i]->field_name)); 1816 CeedCall(CeedFree(&(*op)->output_fields[i])); 1817 } 1818 } 1819 // Destroy sub_operators 1820 for (CeedInt i = 0; i < (*op)->num_suboperators; i++) { 1821 if ((*op)->sub_operators[i]) { 1822 CeedCall(CeedOperatorDestroy(&(*op)->sub_operators[i])); 1823 } 1824 } 1825 CeedCall(CeedQFunctionDestroy(&(*op)->qf)); 1826 CeedCall(CeedQFunctionDestroy(&(*op)->dqf)); 1827 CeedCall(CeedQFunctionDestroy(&(*op)->dqfT)); 1828 // Destroy any composite labels 1829 if ((*op)->is_composite) { 1830 for (CeedInt i = 0; i < (*op)->num_context_labels; i++) { 1831 CeedCall(CeedFree(&(*op)->context_labels[i]->sub_labels)); 1832 CeedCall(CeedFree(&(*op)->context_labels[i])); 1833 } 1834 } 1835 CeedCall(CeedFree(&(*op)->context_labels)); 1836 1837 // Destroy fallback 1838 CeedCall(CeedOperatorDestroy(&(*op)->op_fallback)); 1839 1840 // Destroy assembly data 1841 CeedCall(CeedQFunctionAssemblyDataDestroy(&(*op)->qf_assembled)); 1842 CeedCall(CeedOperatorAssemblyDataDestroy(&(*op)->op_assembled)); 1843 1844 CeedCall(CeedFree(&(*op)->input_fields)); 1845 CeedCall(CeedFree(&(*op)->output_fields)); 1846 CeedCall(CeedFree(&(*op)->sub_operators)); 1847 CeedCall(CeedFree(&(*op)->name)); 1848 CeedCall(CeedFree(op)); 1849 return CEED_ERROR_SUCCESS; 1850 } 1851 1852 /// @} 1853