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