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