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 the name of a CeedOperatorField 877 878 @param[in] op_field CeedOperatorField 879 @param[out] field_name Variable to store the field name 880 881 @return An error code: 0 - success, otherwise - failure 882 883 @ref Advanced 884 **/ 885 int CeedOperatorFieldGetName(CeedOperatorField op_field, char **field_name) { 886 *field_name = (char *)op_field->field_name; 887 return CEED_ERROR_SUCCESS; 888 } 889 890 /** 891 @brief Get the CeedElemRestriction of a CeedOperatorField 892 893 @param[in] op_field CeedOperatorField 894 @param[out] rstr Variable to store CeedElemRestriction 895 896 @return An error code: 0 - success, otherwise - failure 897 898 @ref Advanced 899 **/ 900 int CeedOperatorFieldGetElemRestriction(CeedOperatorField op_field, CeedElemRestriction *rstr) { 901 *rstr = op_field->elem_restr; 902 return CEED_ERROR_SUCCESS; 903 } 904 905 /** 906 @brief Get the CeedBasis of a CeedOperatorField 907 908 @param[in] op_field CeedOperatorField 909 @param[out] basis Variable to store CeedBasis 910 911 @return An error code: 0 - success, otherwise - failure 912 913 @ref Advanced 914 **/ 915 int CeedOperatorFieldGetBasis(CeedOperatorField op_field, CeedBasis *basis) { 916 *basis = op_field->basis; 917 return CEED_ERROR_SUCCESS; 918 } 919 920 /** 921 @brief Get the CeedVector of a CeedOperatorField 922 923 @param[in] op_field CeedOperatorField 924 @param[out] vec Variable to store CeedVector 925 926 @return An error code: 0 - success, otherwise - failure 927 928 @ref Advanced 929 **/ 930 int CeedOperatorFieldGetVector(CeedOperatorField op_field, CeedVector *vec) { 931 *vec = op_field->vec; 932 return CEED_ERROR_SUCCESS; 933 } 934 935 /** 936 @brief Add a sub-operator to a composite CeedOperator 937 938 @param[in,out] composite_op Composite CeedOperator 939 @param[in] sub_op Sub-operator CeedOperator 940 941 @return An error code: 0 - success, otherwise - failure 942 943 @ref User 944 */ 945 int CeedCompositeOperatorAddSub(CeedOperator composite_op, CeedOperator sub_op) { 946 if (!composite_op->is_composite) { 947 // LCOV_EXCL_START 948 return CeedError(composite_op->ceed, CEED_ERROR_MINOR, "CeedOperator is not a composite operator"); 949 // LCOV_EXCL_STOP 950 } 951 952 if (composite_op->num_suboperators == CEED_COMPOSITE_MAX) { 953 // LCOV_EXCL_START 954 return CeedError(composite_op->ceed, CEED_ERROR_UNSUPPORTED, "Cannot add additional sub-operators"); 955 // LCOV_EXCL_STOP 956 } 957 if (composite_op->is_immutable) { 958 // LCOV_EXCL_START 959 return CeedError(composite_op->ceed, CEED_ERROR_MAJOR, "Operator cannot be changed after set as immutable"); 960 // LCOV_EXCL_STOP 961 } 962 963 { 964 CeedSize input_size, output_size; 965 CeedCall(CeedOperatorGetActiveVectorLengths(sub_op, &input_size, &output_size)); 966 if (composite_op->input_size == -1) composite_op->input_size = input_size; 967 if (composite_op->output_size == -1) composite_op->output_size = output_size; 968 // Note, a size of -1 means no active vector restriction set, so no incompatibility 969 if ((input_size != -1 && input_size != composite_op->input_size) || (output_size != -1 && output_size != composite_op->output_size)) { 970 // LCOV_EXCL_START 971 return CeedError(composite_op->ceed, CEED_ERROR_MAJOR, 972 "Sub-operators must have compatible dimensions; composite operator of shape (%td, %td) not compatible with sub-operator of " 973 "shape (%td, %td)", 974 composite_op->input_size, composite_op->output_size, input_size, output_size); 975 // LCOV_EXCL_STOP 976 } 977 } 978 979 composite_op->sub_operators[composite_op->num_suboperators] = sub_op; 980 CeedCall(CeedOperatorReference(sub_op)); 981 composite_op->num_suboperators++; 982 return CEED_ERROR_SUCCESS; 983 } 984 985 /** 986 @brief Get the number of sub_operators associated with a CeedOperator 987 988 @param[in] op CeedOperator 989 @param[out] num_suboperators Variable to store number of sub_operators 990 991 @return An error code: 0 - success, otherwise - failure 992 993 @ref Backend 994 **/ 995 996 int CeedCompositeOperatorGetNumSub(CeedOperator op, CeedInt *num_suboperators) { 997 if (!op->is_composite) { 998 // LCOV_EXCL_START 999 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not a composite operator"); 1000 // LCOV_EXCL_STOP 1001 } 1002 *num_suboperators = op->num_suboperators; 1003 return CEED_ERROR_SUCCESS; 1004 } 1005 1006 /** 1007 @brief Get the list of sub_operators associated with a CeedOperator 1008 1009 @param op CeedOperator 1010 @param[out] sub_operators Variable to store list of sub_operators 1011 1012 @return An error code: 0 - success, otherwise - failure 1013 1014 @ref Backend 1015 **/ 1016 1017 int CeedCompositeOperatorGetSubList(CeedOperator op, CeedOperator **sub_operators) { 1018 if (!op->is_composite) { 1019 // LCOV_EXCL_START 1020 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not a composite operator"); 1021 // LCOV_EXCL_STOP 1022 } 1023 *sub_operators = op->sub_operators; 1024 return CEED_ERROR_SUCCESS; 1025 } 1026 1027 /** 1028 @brief Check if a CeedOperator is ready to be used. 1029 1030 @param[in] op CeedOperator to check 1031 1032 @return An error code: 0 - success, otherwise - failure 1033 1034 @ref User 1035 **/ 1036 int CeedOperatorCheckReady(CeedOperator op) { 1037 Ceed ceed; 1038 CeedCall(CeedOperatorGetCeed(op, &ceed)); 1039 1040 if (op->is_interface_setup) return CEED_ERROR_SUCCESS; 1041 1042 CeedQFunction qf = op->qf; 1043 if (op->is_composite) { 1044 if (!op->num_suboperators) { 1045 // Empty operator setup 1046 op->input_size = 0; 1047 op->output_size = 0; 1048 } else { 1049 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1050 CeedCall(CeedOperatorCheckReady(op->sub_operators[i])); 1051 } 1052 // Sub-operators could be modified after adding to composite operator 1053 // Need to verify no lvec incompatibility from any changes 1054 CeedSize input_size, output_size; 1055 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 1056 } 1057 } else { 1058 if (op->num_fields == 0) { 1059 // LCOV_EXCL_START 1060 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "No operator fields set"); 1061 // LCOV_EXCL_STOP 1062 } 1063 if (op->num_fields < qf->num_input_fields + qf->num_output_fields) { 1064 // LCOV_EXCL_START 1065 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "Not all operator fields set"); 1066 // LCOV_EXCL_STOP 1067 } 1068 if (!op->has_restriction) { 1069 // LCOV_EXCL_START 1070 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "At least one restriction required"); 1071 // LCOV_EXCL_STOP 1072 } 1073 if (op->num_qpts == 0) { 1074 // LCOV_EXCL_START 1075 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "At least one non-collocated basis is required or the number of quadrature points must be set"); 1076 // LCOV_EXCL_STOP 1077 } 1078 } 1079 1080 // Flag as immutable and ready 1081 op->is_interface_setup = true; 1082 if (op->qf && op->qf != CEED_QFUNCTION_NONE) op->qf->is_immutable = true; 1083 if (op->dqf && op->dqf != CEED_QFUNCTION_NONE) op->dqf->is_immutable = true; 1084 if (op->dqfT && op->dqfT != CEED_QFUNCTION_NONE) op->dqfT->is_immutable = true; 1085 return CEED_ERROR_SUCCESS; 1086 } 1087 1088 /** 1089 @brief Get vector lengths for the active input and/or output vectors of a CeedOperator. 1090 Note: Lengths of -1 indicate that the CeedOperator does not have an active input and/or output. 1091 1092 @param[in] op CeedOperator 1093 @param[out] input_size Variable to store active input vector length, or NULL 1094 @param[out] output_size Variable to store active output vector length, or NULL 1095 1096 @return An error code: 0 - success, otherwise - failure 1097 1098 @ref User 1099 **/ 1100 int CeedOperatorGetActiveVectorLengths(CeedOperator op, CeedSize *input_size, CeedSize *output_size) { 1101 bool is_composite; 1102 1103 if (input_size) *input_size = op->input_size; 1104 if (output_size) *output_size = op->output_size; 1105 1106 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1107 if (is_composite && (op->input_size == -1 || op->output_size == -1)) { 1108 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1109 CeedSize sub_input_size, sub_output_size; 1110 CeedCall(CeedOperatorGetActiveVectorLengths(op->sub_operators[i], &sub_input_size, &sub_output_size)); 1111 if (op->input_size == -1) op->input_size = sub_input_size; 1112 if (op->output_size == -1) op->output_size = sub_output_size; 1113 // Note, a size of -1 means no active vector restriction set, so no incompatibility 1114 if ((sub_input_size != -1 && sub_input_size != op->input_size) || (sub_output_size != -1 && sub_output_size != op->output_size)) { 1115 // LCOV_EXCL_START 1116 return CeedError(op->ceed, CEED_ERROR_MAJOR, 1117 "Sub-operators must have compatible dimensions; composite operator of shape (%td, %td) not compatible with sub-operator of " 1118 "shape (%td, %td)", 1119 op->input_size, op->output_size, input_size, output_size); 1120 // LCOV_EXCL_STOP 1121 } 1122 } 1123 } 1124 1125 return CEED_ERROR_SUCCESS; 1126 } 1127 1128 /** 1129 @brief Set reuse of CeedQFunction data in CeedOperatorLinearAssemble* functions. 1130 When `reuse_assembly_data = false` (default), the CeedQFunction associated with this CeedOperator is re-assembled every time a 1131 `CeedOperatorLinearAssemble*` function is called. When `reuse_assembly_data = true`, the CeedQFunction associated with this CeedOperator is reused 1132 between calls to `CeedOperatorSetQFunctionAssemblyDataUpdated`. 1133 1134 @param[in] op CeedOperator 1135 @param[in] reuse_assembly_data Boolean flag setting assembly data reuse 1136 1137 @return An error code: 0 - success, otherwise - failure 1138 1139 @ref Advanced 1140 **/ 1141 int CeedOperatorSetQFunctionAssemblyReuse(CeedOperator op, bool reuse_assembly_data) { 1142 bool is_composite; 1143 1144 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1145 if (is_composite) { 1146 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1147 CeedCall(CeedOperatorSetQFunctionAssemblyReuse(op->sub_operators[i], reuse_assembly_data)); 1148 } 1149 } else { 1150 CeedCall(CeedQFunctionAssemblyDataSetReuse(op->qf_assembled, reuse_assembly_data)); 1151 } 1152 1153 return CEED_ERROR_SUCCESS; 1154 } 1155 1156 /** 1157 @brief Mark CeedQFunction data as updated and the CeedQFunction as requiring re-assembly. 1158 1159 @param[in] op CeedOperator 1160 @param[in] needs_data_update Boolean flag setting assembly data reuse 1161 1162 @return An error code: 0 - success, otherwise - failure 1163 1164 @ref Advanced 1165 **/ 1166 int CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(CeedOperator op, bool needs_data_update) { 1167 bool is_composite; 1168 1169 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1170 if (is_composite) { 1171 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1172 CeedCall(CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(op->sub_operators[i], needs_data_update)); 1173 } 1174 } else { 1175 CeedCall(CeedQFunctionAssemblyDataSetUpdateNeeded(op->qf_assembled, needs_data_update)); 1176 } 1177 1178 return CEED_ERROR_SUCCESS; 1179 } 1180 1181 /** 1182 @brief Set the number of quadrature points associated with a CeedOperator. 1183 This should be used when creating a CeedOperator where every field has a collocated basis. 1184 This function cannot be used for composite CeedOperators. 1185 1186 @param[in,out] op CeedOperator 1187 @param[in] num_qpts Number of quadrature points to set 1188 1189 @return An error code: 0 - success, otherwise - failure 1190 1191 @ref Advanced 1192 **/ 1193 int CeedOperatorSetNumQuadraturePoints(CeedOperator op, CeedInt num_qpts) { 1194 if (op->is_composite) { 1195 // LCOV_EXCL_START 1196 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not defined for composite operator"); 1197 // LCOV_EXCL_STOP 1198 } 1199 if (op->num_qpts) { 1200 // LCOV_EXCL_START 1201 return CeedError(op->ceed, CEED_ERROR_MINOR, "Number of quadrature points already defined"); 1202 // LCOV_EXCL_STOP 1203 } 1204 if (op->is_immutable) { 1205 // LCOV_EXCL_START 1206 return CeedError(op->ceed, CEED_ERROR_MAJOR, "Operator cannot be changed after set as immutable"); 1207 // LCOV_EXCL_STOP 1208 } 1209 op->num_qpts = num_qpts; 1210 return CEED_ERROR_SUCCESS; 1211 } 1212 1213 /** 1214 @brief Set name of CeedOperator for CeedOperatorView output 1215 1216 @param[in,out] op CeedOperator 1217 @param[in] name Name to set, or NULL to remove previously set name 1218 1219 @return An error code: 0 - success, otherwise - failure 1220 1221 @ref User 1222 **/ 1223 int CeedOperatorSetName(CeedOperator op, const char *name) { 1224 char *name_copy; 1225 size_t name_len = name ? strlen(name) : 0; 1226 1227 CeedCall(CeedFree(&op->name)); 1228 if (name_len > 0) { 1229 CeedCall(CeedCalloc(name_len + 1, &name_copy)); 1230 memcpy(name_copy, name, name_len); 1231 op->name = name_copy; 1232 } 1233 1234 return CEED_ERROR_SUCCESS; 1235 } 1236 1237 /** 1238 @brief View a CeedOperator 1239 1240 @param[in] op CeedOperator to view 1241 @param[in] stream Stream to write; typically stdout/stderr or a file 1242 1243 @return Error code: 0 - success, otherwise - failure 1244 1245 @ref User 1246 **/ 1247 int CeedOperatorView(CeedOperator op, FILE *stream) { 1248 bool has_name = op->name; 1249 1250 if (op->is_composite) { 1251 fprintf(stream, "Composite CeedOperator%s%s\n", has_name ? " - " : "", has_name ? op->name : ""); 1252 1253 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1254 has_name = op->sub_operators[i]->name; 1255 fprintf(stream, " SubOperator %" CeedInt_FMT "%s%s:\n", i, has_name ? " - " : "", has_name ? op->sub_operators[i]->name : ""); 1256 CeedCall(CeedOperatorSingleView(op->sub_operators[i], 1, stream)); 1257 } 1258 } else { 1259 fprintf(stream, "CeedOperator%s%s\n", has_name ? " - " : "", has_name ? op->name : ""); 1260 CeedCall(CeedOperatorSingleView(op, 0, stream)); 1261 } 1262 return CEED_ERROR_SUCCESS; 1263 } 1264 1265 /** 1266 @brief Get the Ceed associated with a CeedOperator 1267 1268 @param[in] op CeedOperator 1269 @param[out] ceed Variable to store Ceed 1270 1271 @return An error code: 0 - success, otherwise - failure 1272 1273 @ref Advanced 1274 **/ 1275 int CeedOperatorGetCeed(CeedOperator op, Ceed *ceed) { 1276 *ceed = op->ceed; 1277 return CEED_ERROR_SUCCESS; 1278 } 1279 1280 /** 1281 @brief Get the number of elements associated with a CeedOperator 1282 1283 @param[in] op CeedOperator 1284 @param[out] num_elem Variable to store number of elements 1285 1286 @return An error code: 0 - success, otherwise - failure 1287 1288 @ref Advanced 1289 **/ 1290 int CeedOperatorGetNumElements(CeedOperator op, CeedInt *num_elem) { 1291 if (op->is_composite) { 1292 // LCOV_EXCL_START 1293 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not defined for composite operator"); 1294 // LCOV_EXCL_STOP 1295 } 1296 *num_elem = op->num_elem; 1297 return CEED_ERROR_SUCCESS; 1298 } 1299 1300 /** 1301 @brief Get the number of quadrature points associated with a CeedOperator 1302 1303 @param[in] op CeedOperator 1304 @param[out] num_qpts Variable to store vector number of quadrature points 1305 1306 @return An error code: 0 - success, otherwise - failure 1307 1308 @ref Advanced 1309 **/ 1310 int CeedOperatorGetNumQuadraturePoints(CeedOperator op, CeedInt *num_qpts) { 1311 if (op->is_composite) { 1312 // LCOV_EXCL_START 1313 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not defined for composite operator"); 1314 // LCOV_EXCL_STOP 1315 } 1316 *num_qpts = op->num_qpts; 1317 return CEED_ERROR_SUCCESS; 1318 } 1319 1320 /** 1321 @brief Estimate number of FLOPs required to apply CeedOperator on the active vector 1322 1323 @param[in] op CeedOperator to estimate FLOPs for 1324 @param[out] flops Address of variable to hold FLOPs estimate 1325 1326 @ref Backend 1327 **/ 1328 int CeedOperatorGetFlopsEstimate(CeedOperator op, CeedSize *flops) { 1329 bool is_composite; 1330 CeedCall(CeedOperatorCheckReady(op)); 1331 1332 *flops = 0; 1333 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1334 if (is_composite) { 1335 CeedInt num_suboperators; 1336 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1337 CeedOperator *sub_operators; 1338 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1339 1340 // FLOPs for each suboperator 1341 for (CeedInt i = 0; i < num_suboperators; i++) { 1342 CeedSize suboperator_flops; 1343 CeedCall(CeedOperatorGetFlopsEstimate(sub_operators[i], &suboperator_flops)); 1344 *flops += suboperator_flops; 1345 } 1346 } else { 1347 CeedInt num_input_fields, num_output_fields; 1348 CeedOperatorField *input_fields, *output_fields; 1349 1350 CeedCall(CeedOperatorGetFields(op, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); 1351 1352 CeedInt num_elem = 0; 1353 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 1354 // Input FLOPs 1355 for (CeedInt i = 0; i < num_input_fields; i++) { 1356 if (input_fields[i]->vec == CEED_VECTOR_ACTIVE) { 1357 CeedSize restr_flops, basis_flops; 1358 1359 CeedCall(CeedElemRestrictionGetFlopsEstimate(input_fields[i]->elem_restr, CEED_NOTRANSPOSE, &restr_flops)); 1360 *flops += restr_flops; 1361 CeedCall(CeedBasisGetFlopsEstimate(input_fields[i]->basis, CEED_NOTRANSPOSE, op->qf->input_fields[i]->eval_mode, &basis_flops)); 1362 *flops += basis_flops * num_elem; 1363 } 1364 } 1365 // QF FLOPs 1366 CeedInt num_qpts; 1367 CeedSize qf_flops; 1368 CeedCall(CeedOperatorGetNumQuadraturePoints(op, &num_qpts)); 1369 CeedCall(CeedQFunctionGetFlopsEstimate(op->qf, &qf_flops)); 1370 *flops += num_elem * num_qpts * qf_flops; 1371 // Output FLOPs 1372 for (CeedInt i = 0; i < num_output_fields; i++) { 1373 if (output_fields[i]->vec == CEED_VECTOR_ACTIVE) { 1374 CeedSize restr_flops, basis_flops; 1375 1376 CeedCall(CeedElemRestrictionGetFlopsEstimate(output_fields[i]->elem_restr, CEED_TRANSPOSE, &restr_flops)); 1377 *flops += restr_flops; 1378 CeedCall(CeedBasisGetFlopsEstimate(output_fields[i]->basis, CEED_TRANSPOSE, op->qf->output_fields[i]->eval_mode, &basis_flops)); 1379 *flops += basis_flops * num_elem; 1380 } 1381 } 1382 } 1383 1384 return CEED_ERROR_SUCCESS; 1385 } 1386 1387 /** 1388 @brief Get label for a registered QFunctionContext field, or `NULL` if no field has been registered with this `field_name`. 1389 1390 @param[in] op CeedOperator 1391 @param[in] field_name Name of field to retrieve label 1392 @param[out] field_label Variable to field label 1393 1394 @return An error code: 0 - success, otherwise - failure 1395 1396 @ref User 1397 **/ 1398 int CeedOperatorContextGetFieldLabel(CeedOperator op, const char *field_name, CeedContextFieldLabel *field_label) { 1399 bool is_composite; 1400 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1401 1402 if (is_composite) { 1403 // Check if composite label already created 1404 for (CeedInt i = 0; i < op->num_context_labels; i++) { 1405 if (!strcmp(op->context_labels[i]->name, field_name)) { 1406 *field_label = op->context_labels[i]; 1407 return CEED_ERROR_SUCCESS; 1408 } 1409 } 1410 1411 // Create composite label if needed 1412 CeedInt num_sub; 1413 CeedOperator *sub_operators; 1414 CeedContextFieldLabel new_field_label; 1415 1416 CeedCall(CeedCalloc(1, &new_field_label)); 1417 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub)); 1418 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1419 CeedCall(CeedCalloc(num_sub, &new_field_label->sub_labels)); 1420 new_field_label->num_sub_labels = num_sub; 1421 1422 bool label_found = false; 1423 for (CeedInt i = 0; i < num_sub; i++) { 1424 if (sub_operators[i]->qf->ctx) { 1425 CeedContextFieldLabel new_field_label_i; 1426 CeedCall(CeedQFunctionContextGetFieldLabel(sub_operators[i]->qf->ctx, field_name, &new_field_label_i)); 1427 if (new_field_label_i) { 1428 label_found = true; 1429 new_field_label->sub_labels[i] = new_field_label_i; 1430 new_field_label->name = new_field_label_i->name; 1431 new_field_label->description = new_field_label_i->description; 1432 if (new_field_label->type && new_field_label->type != new_field_label_i->type) { 1433 // LCOV_EXCL_START 1434 CeedCall(CeedFree(&new_field_label)); 1435 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "Incompatible field types on sub-operator contexts. %s != %s", 1436 CeedContextFieldTypes[new_field_label->type], CeedContextFieldTypes[new_field_label_i->type]); 1437 // LCOV_EXCL_STOP 1438 } else { 1439 new_field_label->type = new_field_label_i->type; 1440 } 1441 if (new_field_label->num_values != 0 && new_field_label->num_values != new_field_label_i->num_values) { 1442 // LCOV_EXCL_START 1443 CeedCall(CeedFree(&new_field_label)); 1444 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "Incompatible field number of values on sub-operator contexts. %ld != %ld", 1445 new_field_label->num_values, new_field_label_i->num_values); 1446 // LCOV_EXCL_STOP 1447 } else { 1448 new_field_label->num_values = new_field_label_i->num_values; 1449 } 1450 } 1451 } 1452 } 1453 if (!label_found) { 1454 // LCOV_EXCL_START 1455 CeedCall(CeedFree(&new_field_label->sub_labels)); 1456 CeedCall(CeedFree(&new_field_label)); 1457 *field_label = NULL; 1458 // LCOV_EXCL_STOP 1459 } else { 1460 // Move new composite label to operator 1461 if (op->num_context_labels == 0) { 1462 CeedCall(CeedCalloc(1, &op->context_labels)); 1463 op->max_context_labels = 1; 1464 } else if (op->num_context_labels == op->max_context_labels) { 1465 CeedCall(CeedRealloc(2 * op->num_context_labels, &op->context_labels)); 1466 op->max_context_labels *= 2; 1467 } 1468 op->context_labels[op->num_context_labels] = new_field_label; 1469 *field_label = new_field_label; 1470 op->num_context_labels++; 1471 } 1472 1473 return CEED_ERROR_SUCCESS; 1474 } else { 1475 return CeedQFunctionContextGetFieldLabel(op->qf->ctx, field_name, field_label); 1476 } 1477 } 1478 1479 /** 1480 @brief Set QFunctionContext field holding double precision values. 1481 For composite operators, the values are set in all sub-operator QFunctionContexts that have a matching `field_name`. 1482 1483 @param[in,out] op CeedOperator 1484 @param[in] field_label Label of field to set 1485 @param[in] values Values to set 1486 1487 @return An error code: 0 - success, otherwise - failure 1488 1489 @ref User 1490 **/ 1491 int CeedOperatorContextSetDouble(CeedOperator op, CeedContextFieldLabel field_label, double *values) { 1492 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, values); 1493 } 1494 1495 /** 1496 @brief Get QFunctionContext field holding double precision values, read-only. 1497 For composite operators, the values correspond to the first sub-operator QFunctionContexts that has a matching `field_name`. 1498 1499 @param[in] op CeedOperator 1500 @param[in] field_label Label of field to get 1501 @param[out] num_values Number of values in the field label 1502 @param[out] values Pointer to context values 1503 1504 @return An error code: 0 - success, otherwise - failure 1505 1506 @ref User 1507 **/ 1508 int CeedOperatorContextGetDoubleRead(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const double **values) { 1509 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, num_values, values); 1510 } 1511 1512 /** 1513 @brief Restore QFunctionContext field holding double precision values, read-only. 1514 1515 @param[in] op CeedOperator 1516 @param[in] field_label Label of field to restore 1517 @param[out] values Pointer to context values 1518 1519 @return An error code: 0 - success, otherwise - failure 1520 1521 @ref User 1522 **/ 1523 int CeedOperatorContextRestoreDoubleRead(CeedOperator op, CeedContextFieldLabel field_label, const double **values) { 1524 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, values); 1525 } 1526 1527 /** 1528 @brief Set QFunctionContext field holding int32 values. 1529 For composite operators, the values are set in all sub-operator QFunctionContexts that have a matching `field_name`. 1530 1531 @param[in,out] op CeedOperator 1532 @param[in] field_label Label of field to set 1533 @param[in] values Values to set 1534 1535 @return An error code: 0 - success, otherwise - failure 1536 1537 @ref User 1538 **/ 1539 int CeedOperatorContextSetInt32(CeedOperator op, CeedContextFieldLabel field_label, int *values) { 1540 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_INT32, values); 1541 } 1542 1543 /** 1544 @brief Get QFunctionContext field holding int32 values, read-only. 1545 For composite operators, the values correspond to the first sub-operator QFunctionContexts that has a matching `field_name`. 1546 1547 @param[in] op CeedOperator 1548 @param[in] field_label Label of field to get 1549 @param[out] num_values Number of values in the field label 1550 @param[out] values Pointer to context values 1551 1552 @return An error code: 0 - success, otherwise - failure 1553 1554 @ref User 1555 **/ 1556 int CeedOperatorContextGetInt32Read(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const int **values) { 1557 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_INT32, num_values, values); 1558 } 1559 1560 /** 1561 @brief Restore QFunctionContext field holding int32 values, read-only. 1562 1563 @param[in] op CeedOperator 1564 @param[in] field_label Label of field to get 1565 @param[out] num_values Number of values in the field label 1566 @param[out] values Pointer to context values 1567 1568 @return An error code: 0 - success, otherwise - failure 1569 1570 @ref User 1571 **/ 1572 int CeedOperatorContextRestoreInt32Read(CeedOperator op, CeedContextFieldLabel field_label, const int **values) { 1573 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_INT32, values); 1574 } 1575 1576 /** 1577 @brief Apply CeedOperator to a vector 1578 1579 This computes the action of the operator on the specified (active) input, yielding its (active) output. 1580 All inputs and outputs must be specified using CeedOperatorSetField(). 1581 1582 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 1583 1584 @param[in] op CeedOperator to apply 1585 @param[in] in CeedVector containing input state or @ref CEED_VECTOR_NONE if there are no active inputs 1586 @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 1587 outputs 1588 @param[in] request Address of CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 1589 1590 @return An error code: 0 - success, otherwise - failure 1591 1592 @ref User 1593 **/ 1594 int CeedOperatorApply(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) { 1595 CeedCall(CeedOperatorCheckReady(op)); 1596 1597 if (op->num_elem) { 1598 // Standard Operator 1599 if (op->Apply) { 1600 CeedCall(op->Apply(op, in, out, request)); 1601 } else { 1602 // Zero all output vectors 1603 CeedQFunction qf = op->qf; 1604 for (CeedInt i = 0; i < qf->num_output_fields; i++) { 1605 CeedVector vec = op->output_fields[i]->vec; 1606 if (vec == CEED_VECTOR_ACTIVE) vec = out; 1607 if (vec != CEED_VECTOR_NONE) { 1608 CeedCall(CeedVectorSetValue(vec, 0.0)); 1609 } 1610 } 1611 // Apply 1612 CeedCall(op->ApplyAdd(op, in, out, request)); 1613 } 1614 } else if (op->is_composite) { 1615 // Composite Operator 1616 if (op->ApplyComposite) { 1617 CeedCall(op->ApplyComposite(op, in, out, request)); 1618 } else { 1619 CeedInt num_suboperators; 1620 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1621 CeedOperator *sub_operators; 1622 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1623 1624 // Zero all output vectors 1625 if (out != CEED_VECTOR_NONE) { 1626 CeedCall(CeedVectorSetValue(out, 0.0)); 1627 } 1628 for (CeedInt i = 0; i < num_suboperators; i++) { 1629 for (CeedInt j = 0; j < sub_operators[i]->qf->num_output_fields; j++) { 1630 CeedVector vec = sub_operators[i]->output_fields[j]->vec; 1631 if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) { 1632 CeedCall(CeedVectorSetValue(vec, 0.0)); 1633 } 1634 } 1635 } 1636 // Apply 1637 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1638 CeedCall(CeedOperatorApplyAdd(op->sub_operators[i], in, out, request)); 1639 } 1640 } 1641 } 1642 return CEED_ERROR_SUCCESS; 1643 } 1644 1645 /** 1646 @brief Apply CeedOperator to a vector and add result to output vector 1647 1648 This computes the action of the operator on the specified (active) input, yielding its (active) output. 1649 All inputs and outputs must be specified using CeedOperatorSetField(). 1650 1651 @param[in] op CeedOperator to apply 1652 @param[in] in CeedVector containing input state or NULL if there are no active inputs 1653 @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 1654 @param[in] request Address of CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 1655 1656 @return An error code: 0 - success, otherwise - failure 1657 1658 @ref User 1659 **/ 1660 int CeedOperatorApplyAdd(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) { 1661 CeedCall(CeedOperatorCheckReady(op)); 1662 1663 if (op->num_elem) { 1664 // Standard Operator 1665 CeedCall(op->ApplyAdd(op, in, out, request)); 1666 } else if (op->is_composite) { 1667 // Composite Operator 1668 if (op->ApplyAddComposite) { 1669 CeedCall(op->ApplyAddComposite(op, in, out, request)); 1670 } else { 1671 CeedInt num_suboperators; 1672 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1673 CeedOperator *sub_operators; 1674 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1675 1676 for (CeedInt i = 0; i < num_suboperators; i++) { 1677 CeedCall(CeedOperatorApplyAdd(sub_operators[i], in, out, request)); 1678 } 1679 } 1680 } 1681 return CEED_ERROR_SUCCESS; 1682 } 1683 1684 /** 1685 @brief Destroy a CeedOperator 1686 1687 @param[in,out] op CeedOperator to destroy 1688 1689 @return An error code: 0 - success, otherwise - failure 1690 1691 @ref User 1692 **/ 1693 int CeedOperatorDestroy(CeedOperator *op) { 1694 if (!*op || --(*op)->ref_count > 0) return CEED_ERROR_SUCCESS; 1695 if ((*op)->Destroy) CeedCall((*op)->Destroy(*op)); 1696 CeedCall(CeedDestroy(&(*op)->ceed)); 1697 // Free fields 1698 for (CeedInt i = 0; i < (*op)->num_fields; i++) { 1699 if ((*op)->input_fields[i]) { 1700 if ((*op)->input_fields[i]->elem_restr != CEED_ELEMRESTRICTION_NONE) { 1701 CeedCall(CeedElemRestrictionDestroy(&(*op)->input_fields[i]->elem_restr)); 1702 } 1703 if ((*op)->input_fields[i]->basis != CEED_BASIS_COLLOCATED) { 1704 CeedCall(CeedBasisDestroy(&(*op)->input_fields[i]->basis)); 1705 } 1706 if ((*op)->input_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->input_fields[i]->vec != CEED_VECTOR_NONE) { 1707 CeedCall(CeedVectorDestroy(&(*op)->input_fields[i]->vec)); 1708 } 1709 CeedCall(CeedFree(&(*op)->input_fields[i]->field_name)); 1710 CeedCall(CeedFree(&(*op)->input_fields[i])); 1711 } 1712 } 1713 for (CeedInt i = 0; i < (*op)->num_fields; i++) { 1714 if ((*op)->output_fields[i]) { 1715 CeedCall(CeedElemRestrictionDestroy(&(*op)->output_fields[i]->elem_restr)); 1716 if ((*op)->output_fields[i]->basis != CEED_BASIS_COLLOCATED) { 1717 CeedCall(CeedBasisDestroy(&(*op)->output_fields[i]->basis)); 1718 } 1719 if ((*op)->output_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->output_fields[i]->vec != CEED_VECTOR_NONE) { 1720 CeedCall(CeedVectorDestroy(&(*op)->output_fields[i]->vec)); 1721 } 1722 CeedCall(CeedFree(&(*op)->output_fields[i]->field_name)); 1723 CeedCall(CeedFree(&(*op)->output_fields[i])); 1724 } 1725 } 1726 // Destroy sub_operators 1727 for (CeedInt i = 0; i < (*op)->num_suboperators; i++) { 1728 if ((*op)->sub_operators[i]) { 1729 CeedCall(CeedOperatorDestroy(&(*op)->sub_operators[i])); 1730 } 1731 } 1732 CeedCall(CeedQFunctionDestroy(&(*op)->qf)); 1733 CeedCall(CeedQFunctionDestroy(&(*op)->dqf)); 1734 CeedCall(CeedQFunctionDestroy(&(*op)->dqfT)); 1735 // Destroy any composite labels 1736 for (CeedInt i = 0; i < (*op)->num_context_labels; i++) { 1737 CeedCall(CeedFree(&(*op)->context_labels[i]->sub_labels)); 1738 CeedCall(CeedFree(&(*op)->context_labels[i])); 1739 } 1740 CeedCall(CeedFree(&(*op)->context_labels)); 1741 1742 // Destroy fallback 1743 CeedCall(CeedOperatorDestroy(&(*op)->op_fallback)); 1744 1745 // Destroy assembly data 1746 CeedCall(CeedQFunctionAssemblyDataDestroy(&(*op)->qf_assembled)); 1747 CeedCall(CeedOperatorAssemblyDataDestroy(&(*op)->op_assembled)); 1748 1749 CeedCall(CeedFree(&(*op)->input_fields)); 1750 CeedCall(CeedFree(&(*op)->output_fields)); 1751 CeedCall(CeedFree(&(*op)->sub_operators)); 1752 CeedCall(CeedFree(&(*op)->name)); 1753 CeedCall(CeedFree(op)); 1754 return CEED_ERROR_SUCCESS; 1755 } 1756 1757 /// @} 1758