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