1 // Copyright (c) 2017-2024, 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 `CeedQFunction` Field 26 27 @param[in] ceed `Ceed` object for error handling 28 @param[in] qf_field `CeedQFunction` Field matching `CeedOperator` Field 29 @param[in] rstr `CeedOperator` Field `CeedElemRestriction` 30 @param[in] basis `CeedOperator` Field `CeedBasis` 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 rstr, CeedBasis basis) { 37 const char *field_name; 38 CeedInt dim = 1, num_comp = 1, q_comp = 1, rstr_num_comp = 1, size; 39 CeedEvalMode eval_mode; 40 41 // Field data 42 CeedCall(CeedQFunctionFieldGetData(qf_field, &field_name, &size, &eval_mode)); 43 44 // Restriction 45 CeedCheck((rstr == CEED_ELEMRESTRICTION_NONE) == (eval_mode == CEED_EVAL_WEIGHT), ceed, CEED_ERROR_INCOMPATIBLE, 46 "CEED_ELEMRESTRICTION_NONE and CEED_EVAL_WEIGHT must be used together."); 47 if (rstr != CEED_ELEMRESTRICTION_NONE) { 48 CeedCall(CeedElemRestrictionGetNumComponents(rstr, &rstr_num_comp)); 49 } 50 // Basis 51 CeedCheck((basis == CEED_BASIS_NONE) == (eval_mode == CEED_EVAL_NONE), ceed, CEED_ERROR_INCOMPATIBLE, 52 "CEED_BASIS_NONE and CEED_EVAL_NONE must be used together."); 53 if (basis != CEED_BASIS_NONE) { 54 CeedCall(CeedBasisGetDimension(basis, &dim)); 55 CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); 56 CeedCall(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp)); 57 CeedCheck(rstr == CEED_ELEMRESTRICTION_NONE || rstr_num_comp == num_comp, ceed, CEED_ERROR_DIMENSION, 58 "Field '%s' of size %" CeedInt_FMT " and EvalMode %s: CeedElemRestriction has %" CeedInt_FMT 59 " components, but CeedBasis has %" CeedInt_FMT " components", 60 field_name, size, CeedEvalModes[eval_mode], rstr_num_comp, num_comp); 61 } 62 // Field size 63 switch (eval_mode) { 64 case CEED_EVAL_NONE: 65 CeedCheck(size == rstr_num_comp, ceed, CEED_ERROR_DIMENSION, 66 "Field '%s' of size %" CeedInt_FMT " and EvalMode %s: CeedElemRestriction has %" CeedInt_FMT " components", field_name, size, 67 CeedEvalModes[eval_mode], rstr_num_comp); 68 break; 69 case CEED_EVAL_INTERP: 70 case CEED_EVAL_GRAD: 71 case CEED_EVAL_DIV: 72 case CEED_EVAL_CURL: 73 CeedCheck(size == num_comp * q_comp, ceed, CEED_ERROR_DIMENSION, 74 "Field '%s' of size %" CeedInt_FMT " and EvalMode %s: CeedElemRestriction/Basis has %" CeedInt_FMT " components", field_name, size, 75 CeedEvalModes[eval_mode], num_comp * q_comp); 76 break; 77 case CEED_EVAL_WEIGHT: 78 // No additional checks required 79 break; 80 } 81 return CEED_ERROR_SUCCESS; 82 } 83 84 /** 85 @brief View a field of a `CeedOperator` 86 87 @param[in] op_field `CeedOperator` Field to view 88 @param[in] qf_field `CeedQFunction` Field (carries field name) 89 @param[in] field_number Number of field being viewed 90 @param[in] sub true indicates sub-operator, which increases indentation; false for top-level operator 91 @param[in] input true for an input field; false for output field 92 @param[in] stream Stream to view to, e.g., `stdout` 93 94 @return An error code: 0 - success, otherwise - failure 95 96 @ref Utility 97 **/ 98 static int CeedOperatorFieldView(CeedOperatorField op_field, CeedQFunctionField qf_field, CeedInt field_number, bool sub, bool input, FILE *stream) { 99 const char *pre = sub ? " " : ""; 100 const char *in_out = input ? "Input" : "Output"; 101 const char *field_name; 102 CeedInt size; 103 CeedEvalMode eval_mode; 104 CeedVector vec; 105 CeedBasis basis; 106 107 // Field data 108 CeedCall(CeedQFunctionFieldGetData(qf_field, &field_name, &size, &eval_mode)); 109 CeedCall(CeedOperatorFieldGetData(op_field, NULL, NULL, &basis, &vec)); 110 111 fprintf(stream, 112 "%s %s field %" CeedInt_FMT 113 ":\n" 114 "%s Name: \"%s\"\n", 115 pre, in_out, field_number, pre, field_name); 116 fprintf(stream, "%s Size: %" CeedInt_FMT "\n", pre, size); 117 fprintf(stream, "%s EvalMode: %s\n", pre, CeedEvalModes[eval_mode]); 118 if (basis == CEED_BASIS_NONE) fprintf(stream, "%s No basis\n", pre); 119 if (vec == CEED_VECTOR_ACTIVE) fprintf(stream, "%s Active vector\n", pre); 120 else if (vec == CEED_VECTOR_NONE) fprintf(stream, "%s No vector\n", pre); 121 return CEED_ERROR_SUCCESS; 122 } 123 124 /** 125 @brief View a single `CeedOperator` 126 127 @param[in] op `CeedOperator` to view 128 @param[in] sub Boolean flag for sub-operator 129 @param[in] stream Stream to write; typically `stdout` or a file 130 131 @return Error code: 0 - success, otherwise - failure 132 133 @ref Utility 134 **/ 135 int CeedOperatorSingleView(CeedOperator op, bool sub, FILE *stream) { 136 const char *pre = sub ? " " : ""; 137 CeedInt num_elem, num_qpts, total_fields = 0, num_input_fields, num_output_fields; 138 CeedQFunction qf; 139 CeedQFunctionField *qf_input_fields, *qf_output_fields; 140 CeedOperatorField *op_input_fields, *op_output_fields; 141 142 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 143 CeedCall(CeedOperatorGetNumQuadraturePoints(op, &num_qpts)); 144 CeedCall(CeedOperatorGetNumArgs(op, &total_fields)); 145 CeedCall(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 146 CeedCall(CeedOperatorGetQFunction(op, &qf)); 147 CeedCall(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 148 149 fprintf(stream, "%s %" CeedInt_FMT " elements with %" CeedInt_FMT " quadrature points each\n", pre, num_elem, num_qpts); 150 fprintf(stream, "%s %" CeedInt_FMT " field%s\n", pre, total_fields, total_fields > 1 ? "s" : ""); 151 fprintf(stream, "%s %" CeedInt_FMT " input field%s:\n", pre, num_input_fields, num_input_fields > 1 ? "s" : ""); 152 for (CeedInt i = 0; i < num_input_fields; i++) { 153 CeedCall(CeedOperatorFieldView(op_input_fields[i], qf_input_fields[i], i, sub, 1, stream)); 154 } 155 fprintf(stream, "%s %" CeedInt_FMT " output field%s:\n", pre, num_output_fields, num_output_fields > 1 ? "s" : ""); 156 for (CeedInt i = 0; i < num_output_fields; i++) { 157 CeedCall(CeedOperatorFieldView(op_output_fields[i], qf_output_fields[i], i, sub, 0, stream)); 158 } 159 return CEED_ERROR_SUCCESS; 160 } 161 162 /** 163 @brief Find the active input vector `CeedBasis` for a non-composite `CeedOperator` 164 165 @param[in] op `CeedOperator` to find active `CeedBasis` for 166 @param[out] active_basis `CeedBasis` for active input vector or `NULL` for composite operator 167 168 @return An error code: 0 - success, otherwise - failure 169 170 @ref Developer 171 **/ 172 int CeedOperatorGetActiveBasis(CeedOperator op, CeedBasis *active_basis) { 173 CeedCall(CeedOperatorGetActiveBases(op, active_basis, NULL)); 174 return CEED_ERROR_SUCCESS; 175 } 176 177 /** 178 @brief Find the active input and output vector `CeedBasis` for a non-composite `CeedOperator` 179 180 @param[in] op `CeedOperator` to find active `CeedBasis` for 181 @param[out] active_input_basis `CeedBasis` for active input vector or `NULL` for composite operator 182 @param[out] active_output_basis `CeedBasis` for active output vector or `NULL` for composite operator 183 184 @return An error code: 0 - success, otherwise - failure 185 186 @ref Developer 187 **/ 188 int CeedOperatorGetActiveBases(CeedOperator op, CeedBasis *active_input_basis, CeedBasis *active_output_basis) { 189 bool is_composite; 190 CeedInt num_input_fields, num_output_fields; 191 Ceed ceed; 192 CeedOperatorField *op_input_fields, *op_output_fields; 193 194 CeedCall(CeedOperatorGetCeed(op, &ceed)); 195 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 196 CeedCall(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 197 198 if (active_input_basis) { 199 *active_input_basis = NULL; 200 if (!is_composite) { 201 for (CeedInt i = 0; i < num_input_fields; i++) { 202 CeedVector vec; 203 204 CeedCall(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 205 if (vec == CEED_VECTOR_ACTIVE) { 206 CeedBasis basis; 207 208 CeedCall(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 209 CeedCheck(!*active_input_basis || *active_input_basis == basis, ceed, CEED_ERROR_MINOR, "Multiple active input CeedBases found"); 210 *active_input_basis = basis; 211 } 212 } 213 CeedCheck(*active_input_basis, ceed, CEED_ERROR_INCOMPLETE, "No active input CeedBasis found"); 214 } 215 } 216 if (active_output_basis) { 217 *active_output_basis = NULL; 218 if (!is_composite) { 219 for (CeedInt i = 0; i < num_output_fields; i++) { 220 CeedVector vec; 221 222 CeedCall(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 223 if (vec == CEED_VECTOR_ACTIVE) { 224 CeedBasis basis; 225 226 CeedCall(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 227 CeedCheck(!*active_output_basis || *active_output_basis == basis, ceed, CEED_ERROR_MINOR, "Multiple active output CeedBases found"); 228 *active_output_basis = basis; 229 } 230 } 231 CeedCheck(*active_output_basis, ceed, CEED_ERROR_INCOMPLETE, "No active output CeedBasis found"); 232 } 233 } 234 return CEED_ERROR_SUCCESS; 235 } 236 237 /** 238 @brief Find the active vector `CeedElemRestriction` for a non-composite `CeedOperator` 239 240 @param[in] op `CeedOperator` to find active `CeedElemRestriction` for 241 @param[out] active_rstr `CeedElemRestriction` for active input vector or NULL for composite operator 242 243 @return An error code: 0 - success, otherwise - failure 244 245 @ref Utility 246 **/ 247 int CeedOperatorGetActiveElemRestriction(CeedOperator op, CeedElemRestriction *active_rstr) { 248 CeedCall(CeedOperatorGetActiveElemRestrictions(op, active_rstr, NULL)); 249 return CEED_ERROR_SUCCESS; 250 } 251 252 /** 253 @brief Find the active input and output vector `CeedElemRestriction` for a non-composite `CeedOperator` 254 255 @param[in] op `CeedOperator` to find active `CeedElemRestriction` for 256 @param[out] active_input_rstr `CeedElemRestriction` for active input vector or NULL for composite operator 257 @param[out] active_output_rstr `CeedElemRestriction` for active output vector or NULL for composite operator 258 259 @return An error code: 0 - success, otherwise - failure 260 261 @ref Utility 262 **/ 263 int CeedOperatorGetActiveElemRestrictions(CeedOperator op, CeedElemRestriction *active_input_rstr, CeedElemRestriction *active_output_rstr) { 264 bool is_composite; 265 CeedInt num_input_fields, num_output_fields; 266 Ceed ceed; 267 CeedOperatorField *op_input_fields, *op_output_fields; 268 269 CeedCall(CeedOperatorGetCeed(op, &ceed)); 270 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 271 CeedCall(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 272 273 if (active_input_rstr) { 274 *active_input_rstr = NULL; 275 if (!is_composite) { 276 for (CeedInt i = 0; i < num_input_fields; i++) { 277 CeedVector vec; 278 279 CeedCall(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 280 if (vec == CEED_VECTOR_ACTIVE) { 281 CeedElemRestriction rstr; 282 283 CeedCall(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &rstr)); 284 CeedCheck(!*active_input_rstr || *active_input_rstr == rstr, ceed, CEED_ERROR_MINOR, "Multiple active input CeedElemRestrictions found"); 285 *active_input_rstr = rstr; 286 } 287 } 288 CeedCheck(*active_input_rstr, ceed, CEED_ERROR_INCOMPLETE, "No active input CeedElemRestriction found"); 289 } 290 } 291 if (active_output_rstr) { 292 *active_output_rstr = NULL; 293 if (!is_composite) { 294 for (CeedInt i = 0; i < num_output_fields; i++) { 295 CeedVector vec; 296 297 CeedCall(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 298 if (vec == CEED_VECTOR_ACTIVE) { 299 CeedElemRestriction rstr; 300 301 CeedCall(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &rstr)); 302 CeedCheck(!*active_output_rstr || *active_output_rstr == rstr, ceed, CEED_ERROR_MINOR, "Multiple active output CeedElemRestrictions found"); 303 *active_output_rstr = rstr; 304 } 305 } 306 CeedCheck(*active_output_rstr, ceed, CEED_ERROR_INCOMPLETE, "No active output CeedElemRestriction found"); 307 } 308 } 309 return CEED_ERROR_SUCCESS; 310 } 311 312 /** 313 @brief Set `CeedQFunctionContext` field values of the specified type. 314 315 For composite operators, the value is set in all sub-operator `CeedQFunctionContext` that have a matching `field_name`. 316 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 not have any field of a matching type. 317 318 @param[in,out] op `CeedOperator` 319 @param[in] field_label Label of field to set 320 @param[in] field_type Type of field to set 321 @param[in] values Values to set 322 323 @return An error code: 0 - success, otherwise - failure 324 325 @ref User 326 **/ 327 static int CeedOperatorContextSetGeneric(CeedOperator op, CeedContextFieldLabel field_label, CeedContextFieldType field_type, void *values) { 328 bool is_composite = false; 329 Ceed ceed; 330 331 CeedCall(CeedOperatorGetCeed(op, &ceed)); 332 CeedCheck(field_label, ceed, CEED_ERROR_UNSUPPORTED, "Invalid field label"); 333 334 // Check if field_label and op correspond 335 if (field_label->from_op) { 336 CeedInt index = -1; 337 338 for (CeedInt i = 0; i < op->num_context_labels; i++) { 339 if (op->context_labels[i] == field_label) index = i; 340 } 341 CeedCheck(index != -1, ceed, CEED_ERROR_UNSUPPORTED, "ContextFieldLabel does not correspond to the operator"); 342 } 343 344 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 345 if (is_composite) { 346 CeedInt num_sub; 347 CeedOperator *sub_operators; 348 349 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub)); 350 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 351 CeedCheck(num_sub == field_label->num_sub_labels, ceed, CEED_ERROR_UNSUPPORTED, "Composite operator modified after ContextFieldLabel created"); 352 353 for (CeedInt i = 0; i < num_sub; i++) { 354 CeedQFunction qf; 355 CeedQFunctionContext ctx; 356 357 CeedCall(CeedOperatorGetQFunction(sub_operators[i], &qf)); 358 CeedCall(CeedQFunctionGetContext(qf, &ctx)); 359 // Try every sub-operator, ok if some sub-operators do not have field 360 if (field_label->sub_labels[i] && ctx) { 361 CeedCall(CeedQFunctionContextSetGeneric(ctx, field_label->sub_labels[i], field_type, values)); 362 } 363 } 364 } else { 365 CeedQFunction qf; 366 CeedQFunctionContext ctx; 367 368 CeedCall(CeedOperatorGetQFunction(op, &qf)); 369 CeedCall(CeedQFunctionGetContext(qf, &ctx)); 370 CeedCheck(ctx, ceed, CEED_ERROR_UNSUPPORTED, "QFunction does not have context data"); 371 CeedCall(CeedQFunctionContextSetGeneric(ctx, field_label, field_type, values)); 372 } 373 CeedCall(CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(op, true)); 374 return CEED_ERROR_SUCCESS; 375 } 376 377 /** 378 @brief Get `CeedQFunctionContext` field values of the specified type, read-only. 379 380 For composite operators, the values retrieved are for the first sub-operator `CeedQFunctionContext` that have a matching `field_name`. 381 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 not have any field of a matching type. 382 383 @param[in,out] op `CeedOperator` 384 @param[in] field_label Label of field to set 385 @param[in] field_type Type of field to set 386 @param[out] num_values Number of values of type `field_type` in array `values` 387 @param[out] values Values in the label 388 389 @return An error code: 0 - success, otherwise - failure 390 391 @ref User 392 **/ 393 static int CeedOperatorContextGetGenericRead(CeedOperator op, CeedContextFieldLabel field_label, CeedContextFieldType field_type, size_t *num_values, 394 void *values) { 395 bool is_composite = false; 396 Ceed ceed; 397 398 CeedCall(CeedOperatorGetCeed(op, &ceed)); 399 CeedCheck(field_label, ceed, CEED_ERROR_UNSUPPORTED, "Invalid field label"); 400 401 *(void **)values = NULL; 402 *num_values = 0; 403 404 // Check if field_label and op correspond 405 if (field_label->from_op) { 406 CeedInt index = -1; 407 408 for (CeedInt i = 0; i < op->num_context_labels; i++) { 409 if (op->context_labels[i] == field_label) index = i; 410 } 411 CeedCheck(index != -1, ceed, CEED_ERROR_UNSUPPORTED, "ContextFieldLabel does not correspond to the operator"); 412 } 413 414 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 415 if (is_composite) { 416 CeedInt num_sub; 417 CeedOperator *sub_operators; 418 419 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub)); 420 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 421 CeedCheck(num_sub == field_label->num_sub_labels, ceed, CEED_ERROR_UNSUPPORTED, "Composite operator modified after ContextFieldLabel created"); 422 423 for (CeedInt i = 0; i < num_sub; i++) { 424 CeedQFunction qf; 425 CeedQFunctionContext ctx; 426 427 CeedCall(CeedOperatorGetQFunction(sub_operators[i], &qf)); 428 CeedCall(CeedQFunctionGetContext(qf, &ctx)); 429 // Try every sub-operator, ok if some sub-operators do not have field 430 if (field_label->sub_labels[i] && ctx) { 431 CeedCall(CeedQFunctionContextGetGenericRead(ctx, field_label->sub_labels[i], field_type, num_values, values)); 432 return CEED_ERROR_SUCCESS; 433 } 434 } 435 } else { 436 CeedQFunction qf; 437 CeedQFunctionContext ctx; 438 439 CeedCall(CeedOperatorGetQFunction(op, &qf)); 440 CeedCall(CeedQFunctionGetContext(qf, &ctx)); 441 CeedCheck(ctx, ceed, CEED_ERROR_UNSUPPORTED, "QFunction does not have context data"); 442 CeedCall(CeedQFunctionContextGetGenericRead(ctx, field_label, field_type, num_values, values)); 443 } 444 return CEED_ERROR_SUCCESS; 445 } 446 447 /** 448 @brief Restore `CeedQFunctionContext` field values of the specified type, read-only. 449 450 For composite operators, the values restored are for the first sub-operator `CeedQFunctionContext` that have a matching `field_name`. 451 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 not have any field of a matching type. 452 453 @param[in,out] op `CeedOperator` 454 @param[in] field_label Label of field to set 455 @param[in] field_type Type of field to set 456 @param[in] values Values array to restore 457 458 @return An error code: 0 - success, otherwise - failure 459 460 @ref User 461 **/ 462 static int CeedOperatorContextRestoreGenericRead(CeedOperator op, CeedContextFieldLabel field_label, CeedContextFieldType field_type, void *values) { 463 bool is_composite = false; 464 Ceed ceed; 465 466 CeedCall(CeedOperatorGetCeed(op, &ceed)); 467 CeedCheck(field_label, ceed, CEED_ERROR_UNSUPPORTED, "Invalid field label"); 468 469 // Check if field_label and op correspond 470 if (field_label->from_op) { 471 CeedInt index = -1; 472 473 for (CeedInt i = 0; i < op->num_context_labels; i++) { 474 if (op->context_labels[i] == field_label) index = i; 475 } 476 CeedCheck(index != -1, ceed, CEED_ERROR_UNSUPPORTED, "ContextFieldLabel does not correspond to the operator"); 477 } 478 479 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 480 if (is_composite) { 481 CeedInt num_sub; 482 CeedOperator *sub_operators; 483 484 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub)); 485 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 486 CeedCheck(num_sub == field_label->num_sub_labels, ceed, CEED_ERROR_UNSUPPORTED, "Composite operator modified after ContextFieldLabel created"); 487 488 for (CeedInt i = 0; i < num_sub; i++) { 489 CeedQFunction qf; 490 CeedQFunctionContext ctx; 491 492 CeedCall(CeedOperatorGetQFunction(sub_operators[i], &qf)); 493 CeedCall(CeedQFunctionGetContext(qf, &ctx)); 494 // Try every sub-operator, ok if some sub-operators do not have field 495 if (field_label->sub_labels[i] && ctx) { 496 CeedCall(CeedQFunctionContextRestoreGenericRead(ctx, field_label->sub_labels[i], field_type, values)); 497 return CEED_ERROR_SUCCESS; 498 } 499 } 500 } else { 501 CeedQFunction qf; 502 CeedQFunctionContext ctx; 503 504 CeedCall(CeedOperatorGetQFunction(op, &qf)); 505 CeedCall(CeedQFunctionGetContext(qf, &ctx)); 506 CeedCheck(ctx, ceed, CEED_ERROR_UNSUPPORTED, "QFunction does not have context data"); 507 CeedCall(CeedQFunctionContextRestoreGenericRead(ctx, field_label, field_type, values)); 508 } 509 return CEED_ERROR_SUCCESS; 510 } 511 512 /// @} 513 514 /// ---------------------------------------------------------------------------- 515 /// CeedOperator Backend API 516 /// ---------------------------------------------------------------------------- 517 /// @addtogroup CeedOperatorBackend 518 /// @{ 519 520 /** 521 @brief Get the number of arguments associated with a `CeedOperator` 522 523 @param[in] op `CeedOperator` 524 @param[out] num_args Variable to store vector number of arguments 525 526 @return An error code: 0 - success, otherwise - failure 527 528 @ref Backend 529 **/ 530 int CeedOperatorGetNumArgs(CeedOperator op, CeedInt *num_args) { 531 bool is_composite; 532 533 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 534 CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Not defined for composite operators"); 535 *num_args = op->num_fields; 536 return CEED_ERROR_SUCCESS; 537 } 538 539 /** 540 @brief Get the tensor product status of all bases for a `CeedOperator`. 541 542 `has_tensor_bases` is only set to `true` if every field uses a tensor-product basis. 543 544 @param[in] op `CeedOperator` 545 @param[out] has_tensor_bases Variable to store tensor bases status 546 547 @return An error code: 0 - success, otherwise - failure 548 549 @ref Backend 550 **/ 551 int CeedOperatorHasTensorBases(CeedOperator op, bool *has_tensor_bases) { 552 CeedInt num_inputs, num_outputs; 553 CeedOperatorField *input_fields, *output_fields; 554 555 CeedCall(CeedOperatorGetFields(op, &num_inputs, &input_fields, &num_outputs, &output_fields)); 556 *has_tensor_bases = true; 557 for (CeedInt i = 0; i < num_inputs; i++) { 558 bool is_tensor; 559 CeedBasis basis; 560 561 CeedCall(CeedOperatorFieldGetBasis(input_fields[i], &basis)); 562 if (basis != CEED_BASIS_NONE) { 563 CeedCall(CeedBasisIsTensor(basis, &is_tensor)); 564 *has_tensor_bases &= is_tensor; 565 } 566 } 567 for (CeedInt i = 0; i < num_outputs; i++) { 568 bool is_tensor; 569 CeedBasis basis; 570 571 CeedCall(CeedOperatorFieldGetBasis(output_fields[i], &basis)); 572 if (basis != CEED_BASIS_NONE) { 573 CeedCall(CeedBasisIsTensor(basis, &is_tensor)); 574 *has_tensor_bases &= is_tensor; 575 } 576 } 577 return CEED_ERROR_SUCCESS; 578 } 579 580 /** 581 @brief Get a boolean value indicating if the `CeedOperator` is immutable 582 583 @param[in] op `CeedOperator` 584 @param[out] is_immutable Variable to store immutability status 585 586 @return An error code: 0 - success, otherwise - failure 587 588 @ref Backend 589 **/ 590 int CeedOperatorIsImmutable(CeedOperator op, bool *is_immutable) { 591 *is_immutable = op->is_immutable; 592 return CEED_ERROR_SUCCESS; 593 } 594 595 /** 596 @brief Get the setup status of a `CeedOperator` 597 598 @param[in] op `CeedOperator` 599 @param[out] is_setup_done Variable to store setup status 600 601 @return An error code: 0 - success, otherwise - failure 602 603 @ref Backend 604 **/ 605 int CeedOperatorIsSetupDone(CeedOperator op, bool *is_setup_done) { 606 *is_setup_done = op->is_backend_setup; 607 return CEED_ERROR_SUCCESS; 608 } 609 610 /** 611 @brief Get the `CeedQFunction` associated with a `CeedOperator` 612 613 @param[in] op `CeedOperator` 614 @param[out] qf Variable to store `CeedQFunction` 615 616 @return An error code: 0 - success, otherwise - failure 617 618 @ref Backend 619 **/ 620 int CeedOperatorGetQFunction(CeedOperator op, CeedQFunction *qf) { 621 bool is_composite; 622 623 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 624 CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Not defined for composite operator"); 625 *qf = op->qf; 626 return CEED_ERROR_SUCCESS; 627 } 628 629 /** 630 @brief Get a boolean value indicating if the `CeedOperator` is composite 631 632 @param[in] op `CeedOperator` 633 @param[out] is_composite Variable to store composite status 634 635 @return An error code: 0 - success, otherwise - failure 636 637 @ref Backend 638 **/ 639 int CeedOperatorIsComposite(CeedOperator op, bool *is_composite) { 640 *is_composite = op->is_composite; 641 return CEED_ERROR_SUCCESS; 642 } 643 644 /** 645 @brief Get the backend data of a `CeedOperator` 646 647 @param[in] op `CeedOperator` 648 @param[out] data Variable to store data 649 650 @return An error code: 0 - success, otherwise - failure 651 652 @ref Backend 653 **/ 654 int CeedOperatorGetData(CeedOperator op, void *data) { 655 *(void **)data = op->data; 656 return CEED_ERROR_SUCCESS; 657 } 658 659 /** 660 @brief Set the backend data of a `CeedOperator` 661 662 @param[in,out] op `CeedOperator` 663 @param[in] data Data to set 664 665 @return An error code: 0 - success, otherwise - failure 666 667 @ref Backend 668 **/ 669 int CeedOperatorSetData(CeedOperator op, void *data) { 670 op->data = data; 671 return CEED_ERROR_SUCCESS; 672 } 673 674 /** 675 @brief Increment the reference counter for a `CeedOperator` 676 677 @param[in,out] op `CeedOperator` to increment the reference counter 678 679 @return An error code: 0 - success, otherwise - failure 680 681 @ref Backend 682 **/ 683 int CeedOperatorReference(CeedOperator op) { 684 op->ref_count++; 685 return CEED_ERROR_SUCCESS; 686 } 687 688 /** 689 @brief Set the setup flag of a `CeedOperator` to `true` 690 691 @param[in,out] op `CeedOperator` 692 693 @return An error code: 0 - success, otherwise - failure 694 695 @ref Backend 696 **/ 697 int CeedOperatorSetSetupDone(CeedOperator op) { 698 op->is_backend_setup = true; 699 return CEED_ERROR_SUCCESS; 700 } 701 702 /// @} 703 704 /// ---------------------------------------------------------------------------- 705 /// CeedOperator Public API 706 /// ---------------------------------------------------------------------------- 707 /// @addtogroup CeedOperatorUser 708 /// @{ 709 710 /** 711 @brief Create a `CeedOperator` and associate a `CeedQFunction`. 712 713 A `CeedBasis` and `CeedElemRestriction` can be associated with `CeedQFunction` fields with @ref CeedOperatorSetField(). 714 715 @param[in] ceed `Ceed` object used to create the `CeedOperator` 716 @param[in] qf `CeedQFunction` defining the action of the operator at quadrature points 717 @param[in] dqf `CeedQFunction` defining the action of the Jacobian of `qf` (or @ref CEED_QFUNCTION_NONE) 718 @param[in] dqfT `CeedQFunction` defining the action of the transpose of the Jacobian of `qf` (or @ref CEED_QFUNCTION_NONE) 719 @param[out] op Address of the variable where the newly created `CeedOperator` will be stored 720 721 @return An error code: 0 - success, otherwise - failure 722 723 @ref User 724 */ 725 int CeedOperatorCreate(Ceed ceed, CeedQFunction qf, CeedQFunction dqf, CeedQFunction dqfT, CeedOperator *op) { 726 if (!ceed->OperatorCreate) { 727 Ceed delegate; 728 729 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Operator")); 730 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedOperatorCreate"); 731 CeedCall(CeedOperatorCreate(delegate, qf, dqf, dqfT, op)); 732 return CEED_ERROR_SUCCESS; 733 } 734 735 CeedCheck(qf && qf != CEED_QFUNCTION_NONE, ceed, CEED_ERROR_MINOR, "Operator must have a valid CeedQFunction."); 736 737 CeedCall(CeedCalloc(1, op)); 738 CeedCall(CeedReferenceCopy(ceed, &(*op)->ceed)); 739 (*op)->ref_count = 1; 740 (*op)->input_size = -1; 741 (*op)->output_size = -1; 742 CeedCall(CeedQFunctionReferenceCopy(qf, &(*op)->qf)); 743 if (dqf && dqf != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionReferenceCopy(dqf, &(*op)->dqf)); 744 if (dqfT && dqfT != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionReferenceCopy(dqfT, &(*op)->dqfT)); 745 CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*op)->input_fields)); 746 CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*op)->output_fields)); 747 CeedCall(ceed->OperatorCreate(*op)); 748 return CEED_ERROR_SUCCESS; 749 } 750 751 /** 752 @brief Create a `CeedOperator` for evaluation at evaluation at arbitrary points in each element. 753 754 A `CeedBasis` and `CeedElemRestriction` can be associated with `CeedQFunction` fields with `CeedOperator` SetField. 755 The locations of each point are set with @ref CeedOperatorAtPointsSetPoints(). 756 757 @param[in] ceed `Ceed` object used to create the `CeedOperator` 758 @param[in] qf `CeedQFunction` defining the action of the operator at quadrature points 759 @param[in] dqf `CeedQFunction` defining the action of the Jacobian of @a qf (or @ref CEED_QFUNCTION_NONE) 760 @param[in] dqfT `CeedQFunction` defining the action of the transpose of the Jacobian of @a qf (or @ref CEED_QFUNCTION_NONE) 761 @param[out] op Address of the variable where the newly created CeedOperator will be stored 762 763 @return An error code: 0 - success, otherwise - failure 764 765 @ref User 766 */ 767 int CeedOperatorCreateAtPoints(Ceed ceed, CeedQFunction qf, CeedQFunction dqf, CeedQFunction dqfT, CeedOperator *op) { 768 if (!ceed->OperatorCreateAtPoints) { 769 Ceed delegate; 770 771 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Operator")); 772 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedOperatorCreateAtPoints"); 773 CeedCall(CeedOperatorCreateAtPoints(delegate, qf, dqf, dqfT, op)); 774 return CEED_ERROR_SUCCESS; 775 } 776 777 CeedCheck(qf && qf != CEED_QFUNCTION_NONE, ceed, CEED_ERROR_MINOR, "Operator must have a valid CeedQFunction."); 778 779 CeedCall(CeedCalloc(1, op)); 780 CeedCall(CeedReferenceCopy(ceed, &(*op)->ceed)); 781 (*op)->ref_count = 1; 782 (*op)->is_at_points = true; 783 (*op)->input_size = -1; 784 (*op)->output_size = -1; 785 CeedCall(CeedQFunctionReferenceCopy(qf, &(*op)->qf)); 786 if (dqf && dqf != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionReferenceCopy(dqf, &(*op)->dqf)); 787 if (dqfT && dqfT != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionReferenceCopy(dqfT, &(*op)->dqfT)); 788 CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*op)->input_fields)); 789 CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*op)->output_fields)); 790 CeedCall(ceed->OperatorCreateAtPoints(*op)); 791 return CEED_ERROR_SUCCESS; 792 } 793 794 /** 795 @brief Create a composite `CeedOperator` that composes the action of several `CeedOperator` 796 797 @param[in] ceed `Ceed` object used to create the `CeedOperator` 798 @param[out] op Address of the variable where the newly created composite `CeedOperator` will be stored 799 800 @return An error code: 0 - success, otherwise - failure 801 802 @ref User 803 */ 804 int CeedCompositeOperatorCreate(Ceed ceed, CeedOperator *op) { 805 if (!ceed->CompositeOperatorCreate) { 806 Ceed delegate; 807 808 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Operator")); 809 if (delegate) { 810 CeedCall(CeedCompositeOperatorCreate(delegate, op)); 811 return CEED_ERROR_SUCCESS; 812 } 813 } 814 815 CeedCall(CeedCalloc(1, op)); 816 CeedCall(CeedReferenceCopy(ceed, &(*op)->ceed)); 817 (*op)->ref_count = 1; 818 (*op)->is_composite = true; 819 CeedCall(CeedCalloc(CEED_COMPOSITE_MAX, &(*op)->sub_operators)); 820 (*op)->input_size = -1; 821 (*op)->output_size = -1; 822 823 if (ceed->CompositeOperatorCreate) CeedCall(ceed->CompositeOperatorCreate(*op)); 824 return CEED_ERROR_SUCCESS; 825 } 826 827 /** 828 @brief Copy the pointer to a `CeedOperator`. 829 830 Both pointers should be destroyed with @ref CeedOperatorDestroy(). 831 832 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`. 833 This `CeedOperator` will be destroyed if `*op_copy` is the only reference to this `CeedOperator`. 834 835 @param[in] op `CeedOperator` to copy reference to 836 @param[in,out] op_copy Variable to store copied reference 837 838 @return An error code: 0 - success, otherwise - failure 839 840 @ref User 841 **/ 842 int CeedOperatorReferenceCopy(CeedOperator op, CeedOperator *op_copy) { 843 CeedCall(CeedOperatorReference(op)); 844 CeedCall(CeedOperatorDestroy(op_copy)); 845 *op_copy = op; 846 return CEED_ERROR_SUCCESS; 847 } 848 849 /** 850 @brief Provide a field to a `CeedOperator` for use by its `CeedQFunction`. 851 852 This function is used to specify both active and passive fields to a `CeedOperator`. 853 For passive fields, a `CeedVector` `vec` must be provided. 854 Passive fields can inputs or outputs (updated in-place when operator is applied). 855 856 Active fields must be specified using this function, but their data (in a `CeedVector`) is passed in @ref CeedOperatorApply(). 857 There can be at most one active input `CeedVector` and at most one active output@ref CeedVector passed to @ref CeedOperatorApply(). 858 859 The number of quadrature points must agree across all points. 860 When using @ref CEED_BASIS_NONE, the number of quadrature points is determined by the element size of `rstr`. 861 862 @param[in,out] op `CeedOperator` on which to provide the field 863 @param[in] field_name Name of the field (to be matched with the name used by `CeedQFunction`) 864 @param[in] rstr `CeedElemRestriction` 865 @param[in] basis `CeedBasis` in which the field resides or @ref CEED_BASIS_NONE if collocated with quadrature points 866 @param[in] vec `CeedVector` to be used by CeedOperator or @ref CEED_VECTOR_ACTIVE if field is active or @ref CEED_VECTOR_NONE if using @ref CEED_EVAL_WEIGHT in the `CeedQFunction` 867 868 @return An error code: 0 - success, otherwise - failure 869 870 @ref User 871 **/ 872 int CeedOperatorSetField(CeedOperator op, const char *field_name, CeedElemRestriction rstr, CeedBasis basis, CeedVector vec) { 873 bool is_input = true, is_at_points, is_composite, is_immutable; 874 CeedInt num_elem = 0, num_qpts = 0, num_input_fields, num_output_fields; 875 Ceed ceed; 876 CeedQFunction qf; 877 CeedQFunctionField qf_field, *qf_input_fields, *qf_output_fields; 878 CeedOperatorField *op_field; 879 880 CeedCall(CeedOperatorGetCeed(op, &ceed)); 881 CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); 882 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 883 CeedCall(CeedOperatorIsImmutable(op, &is_immutable)); 884 CeedCheck(!is_composite, ceed, CEED_ERROR_INCOMPATIBLE, "Cannot add field to composite operator."); 885 CeedCheck(!is_immutable, ceed, CEED_ERROR_MAJOR, "Operator cannot be changed after set as immutable"); 886 CeedCheck(rstr, ceed, CEED_ERROR_INCOMPATIBLE, "CeedElemRestriction rstr for field \"%s\" must be non-NULL.", field_name); 887 CeedCheck(basis, ceed, CEED_ERROR_INCOMPATIBLE, "CeedBasis basis for field \"%s\" must be non-NULL.", field_name); 888 CeedCheck(vec, ceed, CEED_ERROR_INCOMPATIBLE, "CeedVector vec for field \"%s\" must be non-NULL.", field_name); 889 890 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 891 CeedCheck(rstr == CEED_ELEMRESTRICTION_NONE || !op->has_restriction || num_elem == op->num_elem, ceed, CEED_ERROR_DIMENSION, 892 "CeedElemRestriction with %" CeedInt_FMT " elements incompatible with prior %" CeedInt_FMT " elements", num_elem, op->num_elem); 893 { 894 CeedRestrictionType rstr_type; 895 896 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 897 if (rstr_type == CEED_RESTRICTION_POINTS) { 898 CeedCheck(is_at_points, ceed, CEED_ERROR_UNSUPPORTED, "CeedElemRestriction AtPoints not supported for standard operator fields"); 899 CeedCheck(basis == CEED_BASIS_NONE, ceed, CEED_ERROR_UNSUPPORTED, "CeedElemRestriction AtPoints must be used with CEED_BASIS_NONE"); 900 if (!op->first_points_rstr) { 901 CeedCall(CeedElemRestrictionReferenceCopy(rstr, &op->first_points_rstr)); 902 } else { 903 bool are_compatible; 904 905 CeedCall(CeedElemRestrictionAtPointsAreCompatible(op->first_points_rstr, rstr, &are_compatible)); 906 CeedCheck(are_compatible, ceed, CEED_ERROR_INCOMPATIBLE, 907 "CeedElemRestriction must have compatible offsets with previously set CeedElemRestriction"); 908 } 909 } 910 } 911 912 if (basis == CEED_BASIS_NONE) CeedCall(CeedElemRestrictionGetElementSize(rstr, &num_qpts)); 913 else CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); 914 CeedCheck(op->num_qpts == 0 || num_qpts == op->num_qpts, ceed, CEED_ERROR_DIMENSION, 915 "%s must correspond to the same number of quadrature points as previously added CeedBases. Found %" CeedInt_FMT 916 " quadrature points but expected %" CeedInt_FMT " quadrature points.", 917 basis == CEED_BASIS_NONE ? "CeedElemRestriction" : "CeedBasis", num_qpts, op->num_qpts); 918 919 CeedCall(CeedOperatorGetQFunction(op, &qf)); 920 CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &qf_input_fields, &num_output_fields, &qf_output_fields)); 921 for (CeedInt i = 0; i < num_input_fields; i++) { 922 const char *qf_field_name; 923 924 CeedCall(CeedQFunctionFieldGetName(qf_input_fields[i], &qf_field_name)); 925 if (!strcmp(field_name, qf_field_name)) { 926 qf_field = qf_input_fields[i]; 927 op_field = &op->input_fields[i]; 928 goto found; 929 } 930 } 931 is_input = false; 932 for (CeedInt i = 0; i < num_output_fields; i++) { 933 const char *qf_field_name; 934 935 CeedCall(CeedQFunctionFieldGetName(qf_output_fields[i], &qf_field_name)); 936 if (!strcmp(field_name, qf_field_name)) { 937 qf_field = qf_output_fields[i]; 938 op_field = &op->output_fields[i]; 939 goto found; 940 } 941 } 942 // LCOV_EXCL_START 943 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "CeedQFunction has no knowledge of field '%s'", field_name); 944 // LCOV_EXCL_STOP 945 found: 946 CeedCall(CeedOperatorCheckField(ceed, qf_field, rstr, basis)); 947 CeedCall(CeedCalloc(1, op_field)); 948 949 if (vec == CEED_VECTOR_ACTIVE) { 950 CeedSize l_size; 951 952 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &l_size)); 953 if (is_input) { 954 if (op->input_size == -1) op->input_size = l_size; 955 CeedCheck(l_size == op->input_size, ceed, CEED_ERROR_INCOMPATIBLE, 956 "LVector size %" CeedSize_FMT " does not match previous size %" CeedSize_FMT "", l_size, op->input_size); 957 } else { 958 if (op->output_size == -1) op->output_size = l_size; 959 CeedCheck(l_size == op->output_size, ceed, CEED_ERROR_INCOMPATIBLE, 960 "LVector size %" CeedSize_FMT " does not match previous size %" CeedSize_FMT "", l_size, op->output_size); 961 } 962 } 963 964 CeedCall(CeedVectorReferenceCopy(vec, &(*op_field)->vec)); 965 CeedCall(CeedElemRestrictionReferenceCopy(rstr, &(*op_field)->elem_rstr)); 966 if (rstr != CEED_ELEMRESTRICTION_NONE && !op->has_restriction) { 967 op->num_elem = num_elem; 968 op->has_restriction = true; // Restriction set, but num_elem may be 0 969 } 970 CeedCall(CeedBasisReferenceCopy(basis, &(*op_field)->basis)); 971 if (op->num_qpts == 0 && !is_at_points) op->num_qpts = num_qpts; // no consistent number of qpts for OperatorAtPoints 972 op->num_fields += 1; 973 CeedCall(CeedStringAllocCopy(field_name, (char **)&(*op_field)->field_name)); 974 return CEED_ERROR_SUCCESS; 975 } 976 977 /** 978 @brief Get the `CeedOperator` Field of a `CeedOperator`. 979 980 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 981 982 @param[in] op `CeedOperator` 983 @param[out] num_input_fields Variable to store number of input fields 984 @param[out] input_fields Variable to store input fields 985 @param[out] num_output_fields Variable to store number of output fields 986 @param[out] output_fields Variable to store output fields 987 988 @return An error code: 0 - success, otherwise - failure 989 990 @ref Advanced 991 **/ 992 int CeedOperatorGetFields(CeedOperator op, CeedInt *num_input_fields, CeedOperatorField **input_fields, CeedInt *num_output_fields, 993 CeedOperatorField **output_fields) { 994 bool is_composite; 995 CeedQFunction qf; 996 997 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 998 CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Not defined for composite operator"); 999 CeedCall(CeedOperatorCheckReady(op)); 1000 1001 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1002 CeedCall(CeedQFunctionGetFields(qf, num_input_fields, NULL, num_output_fields, NULL)); 1003 if (input_fields) *input_fields = op->input_fields; 1004 if (output_fields) *output_fields = op->output_fields; 1005 return CEED_ERROR_SUCCESS; 1006 } 1007 1008 /** 1009 @brief Set the arbitrary points in each element for a `CeedOperator` at points. 1010 1011 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 1012 1013 @param[in,out] op `CeedOperator` at points 1014 @param[in] rstr_points `CeedElemRestriction` for the coordinates of each point by element 1015 @param[in] point_coords `CeedVector` holding coordinates of each point 1016 1017 @return An error code: 0 - success, otherwise - failure 1018 1019 @ref Advanced 1020 **/ 1021 int CeedOperatorAtPointsSetPoints(CeedOperator op, CeedElemRestriction rstr_points, CeedVector point_coords) { 1022 bool is_at_points, is_immutable; 1023 Ceed ceed; 1024 1025 CeedCall(CeedOperatorGetCeed(op, &ceed)); 1026 CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); 1027 CeedCall(CeedOperatorIsImmutable(op, &is_immutable)); 1028 CeedCheck(is_at_points, ceed, CEED_ERROR_MINOR, "Only defined for operator at points"); 1029 CeedCheck(!is_immutable, ceed, CEED_ERROR_MAJOR, "Operator cannot be changed after set as immutable"); 1030 1031 if (!op->first_points_rstr) { 1032 CeedCall(CeedElemRestrictionReferenceCopy(rstr_points, &op->first_points_rstr)); 1033 } else { 1034 bool are_compatible; 1035 1036 CeedCall(CeedElemRestrictionAtPointsAreCompatible(op->first_points_rstr, rstr_points, &are_compatible)); 1037 CeedCheck(are_compatible, ceed, CEED_ERROR_INCOMPATIBLE, 1038 "CeedElemRestriction must have compatible offsets with previously set field CeedElemRestriction"); 1039 } 1040 1041 CeedCall(CeedElemRestrictionReferenceCopy(rstr_points, &op->rstr_points)); 1042 CeedCall(CeedVectorReferenceCopy(point_coords, &op->point_coords)); 1043 return CEED_ERROR_SUCCESS; 1044 } 1045 1046 /** 1047 @brief Get a boolean value indicating if the `CeedOperator` was created with `CeedOperatorCreateAtPoints` 1048 1049 @param[in] op `CeedOperator` 1050 @param[out] is_at_points Variable to store at points status 1051 1052 @return An error code: 0 - success, otherwise - failure 1053 1054 @ref User 1055 **/ 1056 int CeedOperatorIsAtPoints(CeedOperator op, bool *is_at_points) { 1057 *is_at_points = op->is_at_points; 1058 return CEED_ERROR_SUCCESS; 1059 } 1060 1061 /** 1062 @brief Get the arbitrary points in each element for a `CeedOperator` at points. 1063 1064 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 1065 1066 @param[in] op `CeedOperator` at points 1067 @param[out] rstr_points Variable to hold `CeedElemRestriction` for the coordinates of each point by element 1068 @param[out] point_coords Variable to hold `CeedVector` holding coordinates of each point 1069 1070 @return An error code: 0 - success, otherwise - failure 1071 1072 @ref Advanced 1073 **/ 1074 int CeedOperatorAtPointsGetPoints(CeedOperator op, CeedElemRestriction *rstr_points, CeedVector *point_coords) { 1075 bool is_at_points; 1076 1077 CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); 1078 CeedCheck(is_at_points, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Only defined for operator at points"); 1079 CeedCall(CeedOperatorCheckReady(op)); 1080 1081 if (rstr_points) CeedCall(CeedElemRestrictionReferenceCopy(op->rstr_points, rstr_points)); 1082 if (point_coords) CeedCall(CeedVectorReferenceCopy(op->point_coords, point_coords)); 1083 return CEED_ERROR_SUCCESS; 1084 } 1085 1086 /** 1087 @brief Get a `CeedOperator` Field of a `CeedOperator` from its name. 1088 1089 `op_field` is set to `NULL` if the field is not found. 1090 1091 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 1092 1093 @param[in] op `CeedOperator` 1094 @param[in] field_name Name of desired `CeedOperator` Field 1095 @param[out] op_field `CeedOperator` Field corresponding to the name 1096 1097 @return An error code: 0 - success, otherwise - failure 1098 1099 @ref Advanced 1100 **/ 1101 int CeedOperatorGetFieldByName(CeedOperator op, const char *field_name, CeedOperatorField *op_field) { 1102 const char *name; 1103 CeedInt num_input_fields, num_output_fields; 1104 CeedOperatorField *input_fields, *output_fields; 1105 1106 *op_field = NULL; 1107 CeedCall(CeedOperatorGetFields(op, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); 1108 for (CeedInt i = 0; i < num_input_fields; i++) { 1109 CeedCall(CeedOperatorFieldGetName(input_fields[i], &name)); 1110 if (!strcmp(name, field_name)) { 1111 *op_field = input_fields[i]; 1112 return CEED_ERROR_SUCCESS; 1113 } 1114 } 1115 for (CeedInt i = 0; i < num_output_fields; i++) { 1116 CeedCall(CeedOperatorFieldGetName(output_fields[i], &name)); 1117 if (!strcmp(name, field_name)) { 1118 *op_field = output_fields[i]; 1119 return CEED_ERROR_SUCCESS; 1120 } 1121 } 1122 return CEED_ERROR_SUCCESS; 1123 } 1124 1125 /** 1126 @brief Get the name of a `CeedOperator` Field 1127 1128 @param[in] op_field `CeedOperator` Field 1129 @param[out] field_name Variable to store the field name 1130 1131 @return An error code: 0 - success, otherwise - failure 1132 1133 @ref Advanced 1134 **/ 1135 int CeedOperatorFieldGetName(CeedOperatorField op_field, const char **field_name) { 1136 *field_name = op_field->field_name; 1137 return CEED_ERROR_SUCCESS; 1138 } 1139 1140 /** 1141 @brief Get the `CeedElemRestriction` of a `CeedOperator` Field 1142 1143 @param[in] op_field `CeedOperator` Field 1144 @param[out] rstr Variable to store `CeedElemRestriction` 1145 1146 @return An error code: 0 - success, otherwise - failure 1147 1148 @ref Advanced 1149 **/ 1150 int CeedOperatorFieldGetElemRestriction(CeedOperatorField op_field, CeedElemRestriction *rstr) { 1151 *rstr = op_field->elem_rstr; 1152 return CEED_ERROR_SUCCESS; 1153 } 1154 1155 /** 1156 @brief Get the `CeedBasis` of a `CeedOperator` Field 1157 1158 @param[in] op_field `CeedOperator` Field 1159 @param[out] basis Variable to store `CeedBasis` 1160 1161 @return An error code: 0 - success, otherwise - failure 1162 1163 @ref Advanced 1164 **/ 1165 int CeedOperatorFieldGetBasis(CeedOperatorField op_field, CeedBasis *basis) { 1166 *basis = op_field->basis; 1167 return CEED_ERROR_SUCCESS; 1168 } 1169 1170 /** 1171 @brief Get the `CeedVector` of a `CeedOperator` Field 1172 1173 @param[in] op_field `CeedOperator` Field 1174 @param[out] vec Variable to store `CeedVector` 1175 1176 @return An error code: 0 - success, otherwise - failure 1177 1178 @ref Advanced 1179 **/ 1180 int CeedOperatorFieldGetVector(CeedOperatorField op_field, CeedVector *vec) { 1181 *vec = op_field->vec; 1182 return CEED_ERROR_SUCCESS; 1183 } 1184 1185 /** 1186 @brief Get the data of a `CeedOperator` Field. 1187 1188 Any arguments set as `NULL` are ignored. 1189 1190 @param[in] op_field `CeedOperator` Field 1191 @param[out] field_name Variable to store the field name 1192 @param[out] rstr Variable to store `CeedElemRestriction` 1193 @param[out] basis Variable to store `CeedBasis` 1194 @param[out] vec Variable to store `CeedVector` 1195 1196 @return An error code: 0 - success, otherwise - failure 1197 1198 @ref Advanced 1199 **/ 1200 int CeedOperatorFieldGetData(CeedOperatorField op_field, const char **field_name, CeedElemRestriction *rstr, CeedBasis *basis, CeedVector *vec) { 1201 if (field_name) CeedCall(CeedOperatorFieldGetName(op_field, field_name)); 1202 if (rstr) CeedCall(CeedOperatorFieldGetElemRestriction(op_field, rstr)); 1203 if (basis) CeedCall(CeedOperatorFieldGetBasis(op_field, basis)); 1204 if (vec) CeedCall(CeedOperatorFieldGetVector(op_field, vec)); 1205 return CEED_ERROR_SUCCESS; 1206 } 1207 1208 /** 1209 @brief Add a sub-operator to a composite `CeedOperator` 1210 1211 @param[in,out] composite_op Composite `CeedOperator` 1212 @param[in] sub_op Sub-operator `CeedOperator` 1213 1214 @return An error code: 0 - success, otherwise - failure 1215 1216 @ref User 1217 */ 1218 int CeedCompositeOperatorAddSub(CeedOperator composite_op, CeedOperator sub_op) { 1219 bool is_immutable; 1220 Ceed ceed; 1221 1222 CeedCall(CeedOperatorGetCeed(composite_op, &ceed)); 1223 CeedCheck(composite_op->is_composite, ceed, CEED_ERROR_MINOR, "CeedOperator is not a composite operator"); 1224 CeedCheck(composite_op->num_suboperators < CEED_COMPOSITE_MAX, ceed, CEED_ERROR_UNSUPPORTED, "Cannot add additional sub-operators"); 1225 CeedCall(CeedOperatorIsImmutable(composite_op, &is_immutable)); 1226 CeedCheck(!is_immutable, ceed, CEED_ERROR_MAJOR, "Operator cannot be changed after set as immutable"); 1227 1228 { 1229 CeedSize input_size, output_size; 1230 1231 CeedCall(CeedOperatorGetActiveVectorLengths(sub_op, &input_size, &output_size)); 1232 if (composite_op->input_size == -1) composite_op->input_size = input_size; 1233 if (composite_op->output_size == -1) composite_op->output_size = output_size; 1234 // Note, a size of -1 means no active vector restriction set, so no incompatibility 1235 CeedCheck((input_size == -1 || input_size == composite_op->input_size) && (output_size == -1 || output_size == composite_op->output_size), ceed, 1236 CEED_ERROR_MAJOR, 1237 "Sub-operators must have compatible dimensions; composite operator of shape (%" CeedSize_FMT ", %" CeedSize_FMT 1238 ") not compatible with sub-operator of " 1239 "shape (%" CeedSize_FMT ", %" CeedSize_FMT ")", 1240 composite_op->input_size, composite_op->output_size, input_size, output_size); 1241 } 1242 1243 composite_op->sub_operators[composite_op->num_suboperators] = sub_op; 1244 CeedCall(CeedOperatorReference(sub_op)); 1245 composite_op->num_suboperators++; 1246 return CEED_ERROR_SUCCESS; 1247 } 1248 1249 /** 1250 @brief Get the number of sub-operators associated with a `CeedOperator` 1251 1252 @param[in] op `CeedOperator` 1253 @param[out] num_suboperators Variable to store number of sub-operators 1254 1255 @return An error code: 0 - success, otherwise - failure 1256 1257 @ref Backend 1258 **/ 1259 int CeedCompositeOperatorGetNumSub(CeedOperator op, CeedInt *num_suboperators) { 1260 bool is_composite; 1261 1262 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1263 CeedCheck(is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Only defined for a composite operator"); 1264 *num_suboperators = op->num_suboperators; 1265 return CEED_ERROR_SUCCESS; 1266 } 1267 1268 /** 1269 @brief Get the list of sub-operators associated with a `CeedOperator` 1270 1271 @param[in] op `CeedOperator` 1272 @param[out] sub_operators Variable to store list of sub-operators 1273 1274 @return An error code: 0 - success, otherwise - failure 1275 1276 @ref Backend 1277 **/ 1278 int CeedCompositeOperatorGetSubList(CeedOperator op, CeedOperator **sub_operators) { 1279 bool is_composite; 1280 1281 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1282 CeedCheck(is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Only defined for a composite operator"); 1283 *sub_operators = op->sub_operators; 1284 return CEED_ERROR_SUCCESS; 1285 } 1286 1287 /** 1288 @brief Get a sub `CeedOperator` of a composite `CeedOperator` from its name. 1289 1290 `sub_op` is set to `NULL` if the sub operator is not found. 1291 1292 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 1293 1294 @param[in] op Composite `CeedOperator` 1295 @param[in] op_name Name of desired sub `CeedOperator` 1296 @param[out] sub_op Sub `CeedOperator` corresponding to the name 1297 1298 @return An error code: 0 - success, otherwise - failure 1299 1300 @ref Advanced 1301 **/ 1302 int CeedCompositeOperatorGetSubByName(CeedOperator op, const char *op_name, CeedOperator *sub_op) { 1303 bool is_composite; 1304 CeedInt num_sub_ops; 1305 CeedOperator *sub_ops; 1306 1307 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1308 CeedCheck(is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Only defined for a composite operator"); 1309 *sub_op = NULL; 1310 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub_ops)); 1311 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_ops)); 1312 for (CeedInt i = 0; i < num_sub_ops; i++) { 1313 if (sub_ops[i]->name && !strcmp(op_name, sub_ops[i]->name)) { 1314 *sub_op = sub_ops[i]; 1315 return CEED_ERROR_SUCCESS; 1316 } 1317 } 1318 return CEED_ERROR_SUCCESS; 1319 } 1320 1321 /** 1322 @brief Check if a `CeedOperator` is ready to be used. 1323 1324 @param[in] op `CeedOperator` to check 1325 1326 @return An error code: 0 - success, otherwise - failure 1327 1328 @ref User 1329 **/ 1330 int CeedOperatorCheckReady(CeedOperator op) { 1331 bool is_at_points, is_composite; 1332 Ceed ceed; 1333 CeedQFunction qf = NULL; 1334 1335 if (op->is_interface_setup) return CEED_ERROR_SUCCESS; 1336 1337 CeedCall(CeedOperatorGetCeed(op, &ceed)); 1338 CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); 1339 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1340 if (!is_composite) CeedCall(CeedOperatorGetQFunction(op, &qf)); 1341 if (is_composite) { 1342 CeedInt num_suboperators; 1343 1344 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1345 if (!num_suboperators) { 1346 // Empty operator setup 1347 op->input_size = 0; 1348 op->output_size = 0; 1349 } else { 1350 CeedOperator *sub_operators; 1351 1352 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1353 for (CeedInt i = 0; i < num_suboperators; i++) { 1354 CeedCall(CeedOperatorCheckReady(sub_operators[i])); 1355 } 1356 // Sub-operators could be modified after adding to composite operator 1357 // Need to verify no lvec incompatibility from any changes 1358 CeedSize input_size, output_size; 1359 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 1360 } 1361 } else { 1362 CeedInt num_input_fields, num_output_fields; 1363 1364 CeedCheck(op->num_fields > 0, ceed, CEED_ERROR_INCOMPLETE, "No operator fields set"); 1365 CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, NULL, &num_output_fields, NULL)); 1366 CeedCheck(op->num_fields == num_input_fields + num_output_fields, ceed, CEED_ERROR_INCOMPLETE, "Not all operator fields set"); 1367 CeedCheck(op->has_restriction, ceed, CEED_ERROR_INCOMPLETE, "At least one restriction required"); 1368 CeedCheck(op->num_qpts > 0 || is_at_points, ceed, CEED_ERROR_INCOMPLETE, 1369 "At least one non-collocated CeedBasis is required or the number of quadrature points must be set"); 1370 } 1371 1372 // Flag as immutable and ready 1373 op->is_interface_setup = true; 1374 if (qf && qf != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionSetImmutable(qf)); 1375 if (op->dqf && op->dqf != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionSetImmutable(op->dqf)); 1376 if (op->dqfT && op->dqfT != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionSetImmutable(op->dqfT)); 1377 return CEED_ERROR_SUCCESS; 1378 } 1379 1380 /** 1381 @brief Get vector lengths for the active input and/or output `CeedVector` of a `CeedOperator`. 1382 1383 Note: Lengths of `-1` indicate that the CeedOperator does not have an active input and/or output. 1384 1385 @param[in] op `CeedOperator` 1386 @param[out] input_size Variable to store active input vector length, or `NULL` 1387 @param[out] output_size Variable to store active output vector length, or `NULL` 1388 1389 @return An error code: 0 - success, otherwise - failure 1390 1391 @ref User 1392 **/ 1393 int CeedOperatorGetActiveVectorLengths(CeedOperator op, CeedSize *input_size, CeedSize *output_size) { 1394 bool is_composite; 1395 1396 if (input_size) *input_size = op->input_size; 1397 if (output_size) *output_size = op->output_size; 1398 1399 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1400 if (is_composite && (op->input_size == -1 || op->output_size == -1)) { 1401 CeedInt num_suboperators; 1402 CeedOperator *sub_operators; 1403 1404 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1405 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1406 for (CeedInt i = 0; i < num_suboperators; i++) { 1407 CeedSize sub_input_size, sub_output_size; 1408 1409 CeedCall(CeedOperatorGetActiveVectorLengths(sub_operators[i], &sub_input_size, &sub_output_size)); 1410 if (op->input_size == -1) op->input_size = sub_input_size; 1411 if (op->output_size == -1) op->output_size = sub_output_size; 1412 // Note, a size of -1 means no active vector restriction set, so no incompatibility 1413 CeedCheck((sub_input_size == -1 || sub_input_size == op->input_size) && (sub_output_size == -1 || sub_output_size == op->output_size), 1414 CeedOperatorReturnCeed(op), CEED_ERROR_MAJOR, 1415 "Sub-operators must have compatible dimensions; composite operator of shape (%" CeedSize_FMT ", %" CeedSize_FMT 1416 ") not compatible with sub-operator of " 1417 "shape (%" CeedSize_FMT ", %" CeedSize_FMT ")", 1418 op->input_size, op->output_size, input_size, output_size); 1419 } 1420 } 1421 return CEED_ERROR_SUCCESS; 1422 } 1423 1424 /** 1425 @brief Set reuse of `CeedQFunction` data in `CeedOperatorLinearAssemble*()` functions. 1426 1427 When `reuse_assembly_data = false` (default), the `CeedQFunction` associated with this `CeedOperator` is re-assembled every time a `CeedOperatorLinearAssemble*()` function is called. 1428 When `reuse_assembly_data = true`, the `CeedQFunction` associated with this `CeedOperator` is reused between calls to @ref CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(). 1429 1430 @param[in] op `CeedOperator` 1431 @param[in] reuse_assembly_data Boolean flag setting assembly data reuse 1432 1433 @return An error code: 0 - success, otherwise - failure 1434 1435 @ref Advanced 1436 **/ 1437 int CeedOperatorSetQFunctionAssemblyReuse(CeedOperator op, bool reuse_assembly_data) { 1438 bool is_composite; 1439 1440 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1441 if (is_composite) { 1442 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1443 CeedCall(CeedOperatorSetQFunctionAssemblyReuse(op->sub_operators[i], reuse_assembly_data)); 1444 } 1445 } else { 1446 CeedQFunctionAssemblyData data; 1447 1448 CeedCall(CeedOperatorGetQFunctionAssemblyData(op, &data)); 1449 CeedCall(CeedQFunctionAssemblyDataSetReuse(data, reuse_assembly_data)); 1450 } 1451 return CEED_ERROR_SUCCESS; 1452 } 1453 1454 /** 1455 @brief Mark `CeedQFunction` data as updated and the `CeedQFunction` as requiring re-assembly. 1456 1457 @param[in] op `CeedOperator` 1458 @param[in] needs_data_update Boolean flag setting assembly data reuse 1459 1460 @return An error code: 0 - success, otherwise - failure 1461 1462 @ref Advanced 1463 **/ 1464 int CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(CeedOperator op, bool needs_data_update) { 1465 bool is_composite; 1466 1467 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1468 if (is_composite) { 1469 CeedInt num_suboperators; 1470 CeedOperator *sub_operators; 1471 1472 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1473 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1474 for (CeedInt i = 0; i < num_suboperators; i++) { 1475 CeedCall(CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(sub_operators[i], needs_data_update)); 1476 } 1477 } else { 1478 CeedQFunctionAssemblyData data; 1479 1480 CeedCall(CeedOperatorGetQFunctionAssemblyData(op, &data)); 1481 CeedCall(CeedQFunctionAssemblyDataSetUpdateNeeded(data, needs_data_update)); 1482 } 1483 return CEED_ERROR_SUCCESS; 1484 } 1485 1486 /** 1487 @brief Set name of `CeedOperator` for @ref CeedOperatorView() output 1488 1489 @param[in,out] op `CeedOperator` 1490 @param[in] name Name to set, or NULL to remove previously set name 1491 1492 @return An error code: 0 - success, otherwise - failure 1493 1494 @ref User 1495 **/ 1496 int CeedOperatorSetName(CeedOperator op, const char *name) { 1497 char *name_copy; 1498 size_t name_len = name ? strlen(name) : 0; 1499 1500 CeedCall(CeedFree(&op->name)); 1501 if (name_len > 0) { 1502 CeedCall(CeedCalloc(name_len + 1, &name_copy)); 1503 memcpy(name_copy, name, name_len); 1504 op->name = name_copy; 1505 } 1506 return CEED_ERROR_SUCCESS; 1507 } 1508 1509 /** 1510 @brief Core logic for viewing a `CeedOperator` 1511 1512 @param[in] op `CeedOperator` to view brief summary 1513 @param[in] stream Stream to write; typically `stdout` or a file 1514 @param[in] is_full Whether to write full operator view or terse 1515 1516 @return Error code: 0 - success, otherwise - failure 1517 **/ 1518 static int CeedOperatorView_Core(CeedOperator op, FILE *stream, bool is_full) { 1519 bool has_name = op->name, is_composite; 1520 1521 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1522 if (is_composite) { 1523 CeedInt num_suboperators; 1524 CeedOperator *sub_operators; 1525 1526 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1527 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1528 fprintf(stream, "Composite CeedOperator%s%s\n", has_name ? " - " : "", has_name ? op->name : ""); 1529 1530 for (CeedInt i = 0; i < num_suboperators; i++) { 1531 has_name = sub_operators[i]->name; 1532 fprintf(stream, " SubOperator %" CeedInt_FMT "%s%s%s\n", i, has_name ? " - " : "", has_name ? sub_operators[i]->name : "", is_full ? ":" : ""); 1533 if (is_full) CeedCall(CeedOperatorSingleView(sub_operators[i], 1, stream)); 1534 } 1535 } else { 1536 fprintf(stream, "CeedOperator%s%s\n", has_name ? " - " : "", has_name ? op->name : ""); 1537 if (is_full) CeedCall(CeedOperatorSingleView(op, 0, stream)); 1538 } 1539 return CEED_ERROR_SUCCESS; 1540 } 1541 1542 /** 1543 @brief View a `CeedOperator` 1544 1545 @param[in] op `CeedOperator` to view 1546 @param[in] stream Stream to write; typically `stdout` or a file 1547 1548 @return Error code: 0 - success, otherwise - failure 1549 1550 @ref User 1551 **/ 1552 int CeedOperatorView(CeedOperator op, FILE *stream) { 1553 CeedCall(CeedOperatorView_Core(op, stream, true)); 1554 return CEED_ERROR_SUCCESS; 1555 } 1556 1557 /** 1558 @brief View a brief summary `CeedOperator` 1559 1560 @param[in] op `CeedOperator` to view brief summary 1561 @param[in] stream Stream to write; typically `stdout` or a file 1562 1563 @return Error code: 0 - success, otherwise - failure 1564 1565 @ref User 1566 **/ 1567 int CeedOperatorViewTerse(CeedOperator op, FILE *stream) { 1568 CeedCall(CeedOperatorView_Core(op, stream, false)); 1569 return CEED_ERROR_SUCCESS; 1570 } 1571 1572 /** 1573 @brief Get the `Ceed` associated with a `CeedOperator` 1574 1575 @param[in] op `CeedOperator` 1576 @param[out] ceed Variable to store `Ceed` 1577 1578 @return An error code: 0 - success, otherwise - failure 1579 1580 @ref Advanced 1581 **/ 1582 int CeedOperatorGetCeed(CeedOperator op, Ceed *ceed) { 1583 *ceed = CeedOperatorReturnCeed(op); 1584 return CEED_ERROR_SUCCESS; 1585 } 1586 1587 /** 1588 @brief Return the `Ceed` associated with a `CeedOperator` 1589 1590 @param[in] op `CeedOperator` 1591 1592 @return `Ceed` associated with the `op` 1593 1594 @ref Advanced 1595 **/ 1596 Ceed CeedOperatorReturnCeed(CeedOperator op) { return op->ceed; } 1597 1598 /** 1599 @brief Get the number of elements associated with a `CeedOperator` 1600 1601 @param[in] op `CeedOperator` 1602 @param[out] num_elem Variable to store number of elements 1603 1604 @return An error code: 0 - success, otherwise - failure 1605 1606 @ref Advanced 1607 **/ 1608 int CeedOperatorGetNumElements(CeedOperator op, CeedInt *num_elem) { 1609 bool is_composite; 1610 1611 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1612 CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Not defined for composite operator"); 1613 *num_elem = op->num_elem; 1614 return CEED_ERROR_SUCCESS; 1615 } 1616 1617 /** 1618 @brief Get the number of quadrature points associated with a `CeedOperator` 1619 1620 @param[in] op `CeedOperator` 1621 @param[out] num_qpts Variable to store vector number of quadrature points 1622 1623 @return An error code: 0 - success, otherwise - failure 1624 1625 @ref Advanced 1626 **/ 1627 int CeedOperatorGetNumQuadraturePoints(CeedOperator op, CeedInt *num_qpts) { 1628 bool is_composite; 1629 1630 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1631 CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Not defined for composite operator"); 1632 *num_qpts = op->num_qpts; 1633 return CEED_ERROR_SUCCESS; 1634 } 1635 1636 /** 1637 @brief Estimate number of FLOPs required to apply `CeedOperator` on the active `CeedVector` 1638 1639 @param[in] op `CeedOperator` to estimate FLOPs for 1640 @param[out] flops Address of variable to hold FLOPs estimate 1641 1642 @ref Backend 1643 **/ 1644 int CeedOperatorGetFlopsEstimate(CeedOperator op, CeedSize *flops) { 1645 bool is_composite; 1646 1647 CeedCall(CeedOperatorCheckReady(op)); 1648 1649 *flops = 0; 1650 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1651 if (is_composite) { 1652 CeedInt num_suboperators; 1653 1654 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1655 CeedOperator *sub_operators; 1656 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1657 1658 // FLOPs for each suboperator 1659 for (CeedInt i = 0; i < num_suboperators; i++) { 1660 CeedSize suboperator_flops; 1661 1662 CeedCall(CeedOperatorGetFlopsEstimate(sub_operators[i], &suboperator_flops)); 1663 *flops += suboperator_flops; 1664 } 1665 } else { 1666 CeedInt num_input_fields, num_output_fields, num_elem = 0; 1667 CeedQFunction qf; 1668 CeedQFunctionField *qf_input_fields, *qf_output_fields; 1669 CeedOperatorField *op_input_fields, *op_output_fields; 1670 1671 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1672 CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &qf_input_fields, &num_output_fields, &qf_output_fields)); 1673 CeedCall(CeedOperatorGetFields(op, NULL, &op_input_fields, NULL, &op_output_fields)); 1674 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 1675 1676 // Input FLOPs 1677 for (CeedInt i = 0; i < num_input_fields; i++) { 1678 CeedVector vec; 1679 1680 CeedCall(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1681 if (vec == CEED_VECTOR_ACTIVE) { 1682 CeedEvalMode eval_mode; 1683 CeedSize rstr_flops, basis_flops; 1684 CeedElemRestriction rstr; 1685 CeedBasis basis; 1686 1687 CeedCall(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &rstr)); 1688 CeedCall(CeedElemRestrictionGetFlopsEstimate(rstr, CEED_NOTRANSPOSE, &rstr_flops)); 1689 *flops += rstr_flops; 1690 CeedCall(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 1691 CeedCall(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 1692 CeedCall(CeedBasisGetFlopsEstimate(basis, CEED_NOTRANSPOSE, eval_mode, &basis_flops)); 1693 *flops += basis_flops * num_elem; 1694 } 1695 } 1696 // QF FLOPs 1697 { 1698 CeedInt num_qpts; 1699 CeedSize qf_flops; 1700 CeedQFunction qf; 1701 1702 CeedCall(CeedOperatorGetNumQuadraturePoints(op, &num_qpts)); 1703 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1704 CeedCall(CeedQFunctionGetFlopsEstimate(qf, &qf_flops)); 1705 CeedCheck(qf_flops > -1, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE, 1706 "Must set CeedQFunction FLOPs estimate with CeedQFunctionSetUserFlopsEstimate"); 1707 *flops += num_elem * num_qpts * qf_flops; 1708 } 1709 1710 // Output FLOPs 1711 for (CeedInt i = 0; i < num_output_fields; i++) { 1712 CeedVector vec; 1713 1714 CeedCall(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 1715 if (vec == CEED_VECTOR_ACTIVE) { 1716 CeedEvalMode eval_mode; 1717 CeedSize rstr_flops, basis_flops; 1718 CeedElemRestriction rstr; 1719 CeedBasis basis; 1720 1721 CeedCall(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &rstr)); 1722 CeedCall(CeedElemRestrictionGetFlopsEstimate(rstr, CEED_TRANSPOSE, &rstr_flops)); 1723 *flops += rstr_flops; 1724 CeedCall(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 1725 CeedCall(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 1726 CeedCall(CeedBasisGetFlopsEstimate(basis, CEED_TRANSPOSE, eval_mode, &basis_flops)); 1727 *flops += basis_flops * num_elem; 1728 } 1729 } 1730 } 1731 return CEED_ERROR_SUCCESS; 1732 } 1733 1734 /** 1735 @brief Get `CeedQFunction` global context for a `CeedOperator`. 1736 1737 The caller is responsible for destroying `ctx` returned from this function via @ref CeedQFunctionContextDestroy(). 1738 1739 Note: If the value of `ctx` passed into this function is non-`NULL`, then it is assumed that `ctx` is a pointer to a `CeedQFunctionContext`. 1740 This `CeedQFunctionContext` will be destroyed if `ctx` is the only reference to this `CeedQFunctionContext`. 1741 1742 @param[in] op `CeedOperator` 1743 @param[out] ctx Variable to store `CeedQFunctionContext` 1744 1745 @return An error code: 0 - success, otherwise - failure 1746 1747 @ref Advanced 1748 **/ 1749 int CeedOperatorGetContext(CeedOperator op, CeedQFunctionContext *ctx) { 1750 bool is_composite; 1751 CeedQFunction qf; 1752 CeedQFunctionContext qf_ctx; 1753 1754 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1755 CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE, "Cannot retrieve CeedQFunctionContext for composite operator"); 1756 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1757 CeedCall(CeedQFunctionGetInnerContext(qf, &qf_ctx)); 1758 if (qf_ctx) CeedCall(CeedQFunctionContextReferenceCopy(qf_ctx, ctx)); 1759 return CEED_ERROR_SUCCESS; 1760 } 1761 1762 /** 1763 @brief Get label for a registered `CeedQFunctionContext` field, or `NULL` if no field has been registered with this `field_name`. 1764 1765 Fields are registered via `CeedQFunctionContextRegister*()` functions (eg. @ref CeedQFunctionContextRegisterDouble()). 1766 1767 @param[in] op `CeedOperator` 1768 @param[in] field_name Name of field to retrieve label 1769 @param[out] field_label Variable to field label 1770 1771 @return An error code: 0 - success, otherwise - failure 1772 1773 @ref User 1774 **/ 1775 int CeedOperatorGetContextFieldLabel(CeedOperator op, const char *field_name, CeedContextFieldLabel *field_label) { 1776 bool is_composite, field_found = false; 1777 1778 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1779 1780 if (is_composite) { 1781 // Composite operator 1782 // -- Check if composite label already created 1783 for (CeedInt i = 0; i < op->num_context_labels; i++) { 1784 if (!strcmp(op->context_labels[i]->name, field_name)) { 1785 *field_label = op->context_labels[i]; 1786 return CEED_ERROR_SUCCESS; 1787 } 1788 } 1789 1790 // -- Create composite label if needed 1791 CeedInt num_sub; 1792 CeedOperator *sub_operators; 1793 CeedContextFieldLabel new_field_label; 1794 1795 CeedCall(CeedCalloc(1, &new_field_label)); 1796 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub)); 1797 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1798 CeedCall(CeedCalloc(num_sub, &new_field_label->sub_labels)); 1799 new_field_label->num_sub_labels = num_sub; 1800 1801 for (CeedInt i = 0; i < num_sub; i++) { 1802 if (sub_operators[i]->qf->ctx) { 1803 CeedContextFieldLabel new_field_label_i; 1804 1805 CeedCall(CeedQFunctionContextGetFieldLabel(sub_operators[i]->qf->ctx, field_name, &new_field_label_i)); 1806 if (new_field_label_i) { 1807 field_found = true; 1808 new_field_label->sub_labels[i] = new_field_label_i; 1809 new_field_label->name = new_field_label_i->name; 1810 new_field_label->description = new_field_label_i->description; 1811 if (new_field_label->type && new_field_label->type != new_field_label_i->type) { 1812 // LCOV_EXCL_START 1813 CeedCall(CeedFree(&new_field_label)); 1814 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE, "Incompatible field types on sub-operator contexts. %s != %s", 1815 CeedContextFieldTypes[new_field_label->type], CeedContextFieldTypes[new_field_label_i->type]); 1816 // LCOV_EXCL_STOP 1817 } else { 1818 new_field_label->type = new_field_label_i->type; 1819 } 1820 if (new_field_label->num_values != 0 && new_field_label->num_values != new_field_label_i->num_values) { 1821 // LCOV_EXCL_START 1822 CeedCall(CeedFree(&new_field_label)); 1823 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE, 1824 "Incompatible field number of values on sub-operator contexts. %zu != %zu", new_field_label->num_values, 1825 new_field_label_i->num_values); 1826 // LCOV_EXCL_STOP 1827 } else { 1828 new_field_label->num_values = new_field_label_i->num_values; 1829 } 1830 } 1831 } 1832 } 1833 // -- Cleanup if field was found 1834 if (field_found) { 1835 *field_label = new_field_label; 1836 } else { 1837 // LCOV_EXCL_START 1838 CeedCall(CeedFree(&new_field_label->sub_labels)); 1839 CeedCall(CeedFree(&new_field_label)); 1840 *field_label = NULL; 1841 // LCOV_EXCL_STOP 1842 } 1843 } else { 1844 CeedQFunction qf; 1845 CeedQFunctionContext ctx; 1846 1847 // Single, non-composite operator 1848 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1849 CeedCall(CeedQFunctionGetInnerContext(qf, &ctx)); 1850 if (ctx) { 1851 CeedCall(CeedQFunctionContextGetFieldLabel(ctx, field_name, field_label)); 1852 } else { 1853 *field_label = NULL; 1854 } 1855 } 1856 1857 // Set label in operator 1858 if (*field_label) { 1859 (*field_label)->from_op = true; 1860 1861 // Move new composite label to operator 1862 if (op->num_context_labels == 0) { 1863 CeedCall(CeedCalloc(1, &op->context_labels)); 1864 op->max_context_labels = 1; 1865 } else if (op->num_context_labels == op->max_context_labels) { 1866 CeedCall(CeedRealloc(2 * op->num_context_labels, &op->context_labels)); 1867 op->max_context_labels *= 2; 1868 } 1869 op->context_labels[op->num_context_labels] = *field_label; 1870 op->num_context_labels++; 1871 } 1872 return CEED_ERROR_SUCCESS; 1873 } 1874 1875 /** 1876 @brief Set `CeedQFunctionContext` field holding double precision values. 1877 1878 For composite operators, the values are set in all sub-operator `CeedQFunctionContext` that have a matching `field_name`. 1879 1880 @param[in,out] op `CeedOperator` 1881 @param[in] field_label Label of field to set 1882 @param[in] values Values to set 1883 1884 @return An error code: 0 - success, otherwise - failure 1885 1886 @ref User 1887 **/ 1888 int CeedOperatorSetContextDouble(CeedOperator op, CeedContextFieldLabel field_label, double *values) { 1889 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, values); 1890 } 1891 1892 /** 1893 @brief Get `CeedQFunctionContext` field holding double precision values, read-only. 1894 1895 For composite operators, the values correspond to the first sub-operator `CeedQFunctionContext` that has a matching `field_name`. 1896 1897 @param[in] op `CeedOperator` 1898 @param[in] field_label Label of field to get 1899 @param[out] num_values Number of values in the field label 1900 @param[out] values Pointer to context values 1901 1902 @return An error code: 0 - success, otherwise - failure 1903 1904 @ref User 1905 **/ 1906 int CeedOperatorGetContextDoubleRead(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const double **values) { 1907 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, num_values, values); 1908 } 1909 1910 /** 1911 @brief Restore `CeedQFunctionContext` field holding double precision values, read-only. 1912 1913 @param[in] op `CeedOperator` 1914 @param[in] field_label Label of field to restore 1915 @param[out] values Pointer to context values 1916 1917 @return An error code: 0 - success, otherwise - failure 1918 1919 @ref User 1920 **/ 1921 int CeedOperatorRestoreContextDoubleRead(CeedOperator op, CeedContextFieldLabel field_label, const double **values) { 1922 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, values); 1923 } 1924 1925 /** 1926 @brief Set `CeedQFunctionContext` field holding `int32` values. 1927 1928 For composite operators, the values are set in all sub-operator `CeedQFunctionContext` that have a matching `field_name`. 1929 1930 @param[in,out] op `CeedOperator` 1931 @param[in] field_label Label of field to set 1932 @param[in] values Values to set 1933 1934 @return An error code: 0 - success, otherwise - failure 1935 1936 @ref User 1937 **/ 1938 int CeedOperatorSetContextInt32(CeedOperator op, CeedContextFieldLabel field_label, int32_t *values) { 1939 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_INT32, values); 1940 } 1941 1942 /** 1943 @brief Get `CeedQFunctionContext` field holding `int32` values, read-only. 1944 1945 For composite operators, the values correspond to the first sub-operator `CeedQFunctionContext` that has a matching `field_name`. 1946 1947 @param[in] op `CeedOperator` 1948 @param[in] field_label Label of field to get 1949 @param[out] num_values Number of `int32` values in `values` 1950 @param[out] values Pointer to context values 1951 1952 @return An error code: 0 - success, otherwise - failure 1953 1954 @ref User 1955 **/ 1956 int CeedOperatorGetContextInt32Read(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const int32_t **values) { 1957 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_INT32, num_values, values); 1958 } 1959 1960 /** 1961 @brief Restore `CeedQFunctionContext` field holding `int32` values, read-only. 1962 1963 @param[in] op `CeedOperator` 1964 @param[in] field_label Label of field to get 1965 @param[out] values Pointer to context values 1966 1967 @return An error code: 0 - success, otherwise - failure 1968 1969 @ref User 1970 **/ 1971 int CeedOperatorRestoreContextInt32Read(CeedOperator op, CeedContextFieldLabel field_label, const int32_t **values) { 1972 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_INT32, values); 1973 } 1974 1975 /** 1976 @brief Set `CeedQFunctionContext` field holding boolean values. 1977 1978 For composite operators, the values are set in all sub-operator `CeedQFunctionContext` that have a matching `field_name`. 1979 1980 @param[in,out] op `CeedOperator` 1981 @param[in] field_label Label of field to set 1982 @param[in] values Values to set 1983 1984 @return An error code: 0 - success, otherwise - failure 1985 1986 @ref User 1987 **/ 1988 int CeedOperatorSetContextBoolean(CeedOperator op, CeedContextFieldLabel field_label, bool *values) { 1989 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_BOOL, values); 1990 } 1991 1992 /** 1993 @brief Get `CeedQFunctionContext` field holding boolean values, read-only. 1994 1995 For composite operators, the values correspond to the first sub-operator `CeedQFunctionContext` that has a matching `field_name`. 1996 1997 @param[in] op `CeedOperator` 1998 @param[in] field_label Label of field to get 1999 @param[out] num_values Number of boolean values in `values` 2000 @param[out] values Pointer to context values 2001 2002 @return An error code: 0 - success, otherwise - failure 2003 2004 @ref User 2005 **/ 2006 int CeedOperatorGetContextBooleanRead(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const bool **values) { 2007 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_BOOL, num_values, values); 2008 } 2009 2010 /** 2011 @brief Restore `CeedQFunctionContext` field holding boolean values, read-only. 2012 2013 @param[in] op `CeedOperator` 2014 @param[in] field_label Label of field to get 2015 @param[out] values Pointer to context values 2016 2017 @return An error code: 0 - success, otherwise - failure 2018 2019 @ref User 2020 **/ 2021 int CeedOperatorRestoreContextBooleanRead(CeedOperator op, CeedContextFieldLabel field_label, const bool **values) { 2022 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_BOOL, values); 2023 } 2024 2025 /** 2026 @brief Apply `CeedOperator` to a `CeedVector`. 2027 2028 This computes the action of the operator on the specified (active) input, yielding its (active) output. 2029 All inputs and outputs must be specified using @ref CeedOperatorSetField(). 2030 2031 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2032 2033 @param[in] op `CeedOperator` to apply 2034 @param[in] in `CeedVector` containing input state or @ref CEED_VECTOR_NONE if there are no active inputs 2035 @param[out] out `CeedVector` to store result of applying operator (must be distinct from `in`) or @ref CEED_VECTOR_NONE if there are no active outputs 2036 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2037 2038 @return An error code: 0 - success, otherwise - failure 2039 2040 @ref User 2041 **/ 2042 int CeedOperatorApply(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) { 2043 bool is_composite; 2044 2045 CeedCall(CeedOperatorCheckReady(op)); 2046 2047 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2048 if (is_composite) { 2049 // Composite Operator 2050 if (op->ApplyComposite) { 2051 CeedCall(op->ApplyComposite(op, in, out, request)); 2052 } else { 2053 CeedInt num_suboperators; 2054 CeedOperator *sub_operators; 2055 2056 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 2057 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 2058 2059 // Zero all output vectors 2060 if (out != CEED_VECTOR_NONE) CeedCall(CeedVectorSetValue(out, 0.0)); 2061 for (CeedInt i = 0; i < num_suboperators; i++) { 2062 CeedInt num_output_fields; 2063 CeedOperatorField *output_fields; 2064 2065 CeedCall(CeedOperatorGetFields(sub_operators[i], NULL, NULL, &num_output_fields, &output_fields)); 2066 for (CeedInt j = 0; j < num_output_fields; j++) { 2067 CeedVector vec; 2068 2069 CeedCall(CeedOperatorFieldGetVector(output_fields[j], &vec)); 2070 if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) { 2071 CeedCall(CeedVectorSetValue(vec, 0.0)); 2072 } 2073 } 2074 } 2075 // Apply 2076 for (CeedInt i = 0; i < num_suboperators; i++) { 2077 CeedCall(CeedOperatorApplyAdd(sub_operators[i], in, out, request)); 2078 } 2079 } 2080 } else { 2081 // Standard Operator 2082 if (op->Apply) { 2083 CeedCall(op->Apply(op, in, out, request)); 2084 } else { 2085 CeedInt num_output_fields; 2086 CeedOperatorField *output_fields; 2087 2088 CeedCall(CeedOperatorGetFields(op, NULL, NULL, &num_output_fields, &output_fields)); 2089 // Zero all output vectors 2090 for (CeedInt i = 0; i < num_output_fields; i++) { 2091 CeedVector vec; 2092 2093 CeedCall(CeedOperatorFieldGetVector(output_fields[i], &vec)); 2094 if (vec == CEED_VECTOR_ACTIVE) vec = out; 2095 if (vec != CEED_VECTOR_NONE) CeedCall(CeedVectorSetValue(vec, 0.0)); 2096 } 2097 // Apply 2098 if (op->num_elem > 0) CeedCall(op->ApplyAdd(op, in, out, request)); 2099 } 2100 } 2101 return CEED_ERROR_SUCCESS; 2102 } 2103 2104 /** 2105 @brief Apply `CeedOperator` to a `CeedVector` and add result to output `CeedVector`. 2106 2107 This computes the action of the operator on the specified (active) input, yielding its (active) output. 2108 All inputs and outputs must be specified using @ref CeedOperatorSetField(). 2109 2110 @param[in] op `CeedOperator` to apply 2111 @param[in] in `CeedVector` containing input state or @ref CEED_VECTOR_NONE if there are no active inputs 2112 @param[out] out `CeedVector` to sum in result of applying operator (must be distinct from `in`) or @ref CEED_VECTOR_NONE if there are no active outputs 2113 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2114 2115 @return An error code: 0 - success, otherwise - failure 2116 2117 @ref User 2118 **/ 2119 int CeedOperatorApplyAdd(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) { 2120 bool is_composite; 2121 2122 CeedCall(CeedOperatorCheckReady(op)); 2123 2124 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2125 if (is_composite) { 2126 // Composite Operator 2127 if (op->ApplyAddComposite) { 2128 CeedCall(op->ApplyAddComposite(op, in, out, request)); 2129 } else { 2130 CeedInt num_suboperators; 2131 CeedOperator *sub_operators; 2132 2133 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 2134 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 2135 for (CeedInt i = 0; i < num_suboperators; i++) { 2136 CeedCall(CeedOperatorApplyAdd(sub_operators[i], in, out, request)); 2137 } 2138 } 2139 } else if (op->num_elem > 0) { 2140 // Standard Operator 2141 CeedCall(op->ApplyAdd(op, in, out, request)); 2142 } 2143 return CEED_ERROR_SUCCESS; 2144 } 2145 2146 /** 2147 @brief Destroy temporary assembly data associated with a `CeedOperator` 2148 2149 @param[in,out] op `CeedOperator` whose assembly data to destroy 2150 2151 @return An error code: 0 - success, otherwise - failure 2152 2153 @ref User 2154 **/ 2155 int CeedOperatorAssemblyDataStrip(CeedOperator op) { 2156 bool is_composite; 2157 2158 CeedCall(CeedQFunctionAssemblyDataDestroy(&op->qf_assembled)); 2159 CeedCall(CeedOperatorAssemblyDataDestroy(&op->op_assembled)); 2160 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2161 if (is_composite) { 2162 CeedInt num_suboperators; 2163 CeedOperator *sub_operators; 2164 2165 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 2166 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 2167 for (CeedInt i = 0; i < num_suboperators; i++) { 2168 CeedCall(CeedQFunctionAssemblyDataDestroy(&sub_operators[i]->qf_assembled)); 2169 CeedCall(CeedOperatorAssemblyDataDestroy(&sub_operators[i]->op_assembled)); 2170 } 2171 } 2172 return CEED_ERROR_SUCCESS; 2173 } 2174 2175 /** 2176 @brief Destroy a `CeedOperator` 2177 2178 @param[in,out] op `CeedOperator` to destroy 2179 2180 @return An error code: 0 - success, otherwise - failure 2181 2182 @ref User 2183 **/ 2184 int CeedOperatorDestroy(CeedOperator *op) { 2185 if (!*op || --(*op)->ref_count > 0) { 2186 *op = NULL; 2187 return CEED_ERROR_SUCCESS; 2188 } 2189 // Backend destroy 2190 if ((*op)->Destroy) { 2191 CeedCall((*op)->Destroy(*op)); 2192 } 2193 // Free fields 2194 for (CeedInt i = 0; i < (*op)->num_fields; i++) { 2195 if ((*op)->input_fields[i]) { 2196 if ((*op)->input_fields[i]->elem_rstr != CEED_ELEMRESTRICTION_NONE) { 2197 CeedCall(CeedElemRestrictionDestroy(&(*op)->input_fields[i]->elem_rstr)); 2198 } 2199 if ((*op)->input_fields[i]->basis != CEED_BASIS_NONE) { 2200 CeedCall(CeedBasisDestroy(&(*op)->input_fields[i]->basis)); 2201 } 2202 if ((*op)->input_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->input_fields[i]->vec != CEED_VECTOR_NONE) { 2203 CeedCall(CeedVectorDestroy(&(*op)->input_fields[i]->vec)); 2204 } 2205 CeedCall(CeedFree(&(*op)->input_fields[i]->field_name)); 2206 CeedCall(CeedFree(&(*op)->input_fields[i])); 2207 } 2208 } 2209 for (CeedInt i = 0; i < (*op)->num_fields; i++) { 2210 if ((*op)->output_fields[i]) { 2211 CeedCall(CeedElemRestrictionDestroy(&(*op)->output_fields[i]->elem_rstr)); 2212 if ((*op)->output_fields[i]->basis != CEED_BASIS_NONE) { 2213 CeedCall(CeedBasisDestroy(&(*op)->output_fields[i]->basis)); 2214 } 2215 if ((*op)->output_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->output_fields[i]->vec != CEED_VECTOR_NONE) { 2216 CeedCall(CeedVectorDestroy(&(*op)->output_fields[i]->vec)); 2217 } 2218 CeedCall(CeedFree(&(*op)->output_fields[i]->field_name)); 2219 CeedCall(CeedFree(&(*op)->output_fields[i])); 2220 } 2221 } 2222 CeedCall(CeedFree(&(*op)->input_fields)); 2223 CeedCall(CeedFree(&(*op)->output_fields)); 2224 // Destroy AtPoints data 2225 CeedCall(CeedVectorDestroy(&(*op)->point_coords)); 2226 CeedCall(CeedElemRestrictionDestroy(&(*op)->rstr_points)); 2227 CeedCall(CeedElemRestrictionDestroy(&(*op)->first_points_rstr)); 2228 // Destroy assembly data (must happen before destroying sub_operators) 2229 CeedCall(CeedOperatorAssemblyDataStrip(*op)); 2230 // Destroy sub_operators 2231 for (CeedInt i = 0; i < (*op)->num_suboperators; i++) { 2232 if ((*op)->sub_operators[i]) { 2233 CeedCall(CeedOperatorDestroy(&(*op)->sub_operators[i])); 2234 } 2235 } 2236 CeedCall(CeedFree(&(*op)->sub_operators)); 2237 CeedCall(CeedQFunctionDestroy(&(*op)->qf)); 2238 CeedCall(CeedQFunctionDestroy(&(*op)->dqf)); 2239 CeedCall(CeedQFunctionDestroy(&(*op)->dqfT)); 2240 // Destroy any composite labels 2241 if ((*op)->is_composite) { 2242 for (CeedInt i = 0; i < (*op)->num_context_labels; i++) { 2243 CeedCall(CeedFree(&(*op)->context_labels[i]->sub_labels)); 2244 CeedCall(CeedFree(&(*op)->context_labels[i])); 2245 } 2246 } 2247 CeedCall(CeedFree(&(*op)->context_labels)); 2248 2249 // Destroy fallback 2250 CeedCall(CeedOperatorDestroy(&(*op)->op_fallback)); 2251 2252 CeedCall(CeedFree(&(*op)->name)); 2253 CeedCall(CeedDestroy(&(*op)->ceed)); 2254 CeedCall(CeedFree(op)); 2255 return CEED_ERROR_SUCCESS; 2256 } 2257 2258 /// @} 2259