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