1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 #include <ceed-impl.h> 9 #include <ceed.h> 10 #include <ceed/backend.h> 11 #include <stdbool.h> 12 #include <stdio.h> 13 #include <string.h> 14 15 /// @file 16 /// Implementation of CeedOperator interfaces 17 18 /// ---------------------------------------------------------------------------- 19 /// CeedOperator Library Internal Functions 20 /// ---------------------------------------------------------------------------- 21 /// @addtogroup CeedOperatorDeveloper 22 /// @{ 23 24 /** 25 @brief Check if a `CeedOperator` Field matches the `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, op->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, op->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, op->ceed, CEED_ERROR_UNSUPPORTED, 487 "Composite operator modified after ContextFieldLabel created"); 488 489 for (CeedInt i = 0; i < num_sub; i++) { 490 CeedQFunction qf; 491 CeedQFunctionContext ctx; 492 493 CeedCall(CeedOperatorGetQFunction(sub_operators[i], &qf)); 494 CeedCall(CeedQFunctionGetContext(qf, &ctx)); 495 // Try every sub-operator, ok if some sub-operators do not have field 496 if (field_label->sub_labels[i] && ctx) { 497 CeedCall(CeedQFunctionContextRestoreGenericRead(ctx, field_label->sub_labels[i], field_type, values)); 498 return CEED_ERROR_SUCCESS; 499 } 500 } 501 } else { 502 CeedQFunction qf; 503 CeedQFunctionContext ctx; 504 505 CeedCall(CeedOperatorGetQFunction(op, &qf)); 506 CeedCall(CeedQFunctionGetContext(qf, &ctx)); 507 CeedCheck(ctx, ceed, CEED_ERROR_UNSUPPORTED, "QFunction does not have context data"); 508 CeedCall(CeedQFunctionContextRestoreGenericRead(ctx, field_label, field_type, values)); 509 } 510 return CEED_ERROR_SUCCESS; 511 } 512 513 /// @} 514 515 /// ---------------------------------------------------------------------------- 516 /// CeedOperator Backend API 517 /// ---------------------------------------------------------------------------- 518 /// @addtogroup CeedOperatorBackend 519 /// @{ 520 521 /** 522 @brief Get the number of arguments associated with a `CeedOperator` 523 524 @param[in] op `CeedOperator` 525 @param[out] num_args Variable to store vector number of arguments 526 527 @return An error code: 0 - success, otherwise - failure 528 529 @ref Backend 530 **/ 531 int CeedOperatorGetNumArgs(CeedOperator op, CeedInt *num_args) { 532 bool is_composite; 533 Ceed ceed; 534 535 CeedCall(CeedOperatorGetCeed(op, &ceed)); 536 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 537 CeedCheck(!is_composite, ceed, CEED_ERROR_MINOR, "Not defined for composite operators"); 538 *num_args = op->num_fields; 539 return CEED_ERROR_SUCCESS; 540 } 541 542 /** 543 @brief Get the tensor product status of all bases for a `CeedOperator`. 544 545 `has_tensor_bases` is only set to `true` if every field uses a tensor-product basis. 546 547 @param[in] op `CeedOperator` 548 @param[out] has_tensor_bases Variable to store tensor bases status 549 550 @return An error code: 0 - success, otherwise - failure 551 552 @ref Backend 553 **/ 554 int CeedOperatorHasTensorBases(CeedOperator op, bool *has_tensor_bases) { 555 CeedInt num_inputs, num_outputs; 556 CeedOperatorField *input_fields, *output_fields; 557 558 CeedCall(CeedOperatorGetFields(op, &num_inputs, &input_fields, &num_outputs, &output_fields)); 559 *has_tensor_bases = true; 560 for (CeedInt i = 0; i < num_inputs; i++) { 561 bool is_tensor; 562 CeedBasis basis; 563 564 CeedCall(CeedOperatorFieldGetBasis(input_fields[i], &basis)); 565 if (basis != CEED_BASIS_NONE) { 566 CeedCall(CeedBasisIsTensor(basis, &is_tensor)); 567 *has_tensor_bases &= is_tensor; 568 } 569 } 570 for (CeedInt i = 0; i < num_outputs; i++) { 571 bool is_tensor; 572 CeedBasis basis; 573 574 CeedCall(CeedOperatorFieldGetBasis(output_fields[i], &basis)); 575 if (basis != CEED_BASIS_NONE) { 576 CeedCall(CeedBasisIsTensor(basis, &is_tensor)); 577 *has_tensor_bases &= is_tensor; 578 } 579 } 580 return CEED_ERROR_SUCCESS; 581 } 582 583 /** 584 @brief Get a boolean value indicating if the `CeedOperator` is immutable 585 586 @param[in] op `CeedOperator` 587 @param[out] is_immutable Variable to store immutability status 588 589 @return An error code: 0 - success, otherwise - failure 590 591 @ref Backend 592 **/ 593 int CeedOperatorIsImmutable(CeedOperator op, bool *is_immutable) { 594 *is_immutable = op->is_immutable; 595 return CEED_ERROR_SUCCESS; 596 } 597 598 /** 599 @brief Get the setup status of a `CeedOperator` 600 601 @param[in] op `CeedOperator` 602 @param[out] is_setup_done Variable to store setup status 603 604 @return An error code: 0 - success, otherwise - failure 605 606 @ref Backend 607 **/ 608 int CeedOperatorIsSetupDone(CeedOperator op, bool *is_setup_done) { 609 *is_setup_done = op->is_backend_setup; 610 return CEED_ERROR_SUCCESS; 611 } 612 613 /** 614 @brief Get the `CeedQFunction` associated with a `CeedOperator` 615 616 @param[in] op `CeedOperator` 617 @param[out] qf Variable to store `CeedQFunction` 618 619 @return An error code: 0 - success, otherwise - failure 620 621 @ref Backend 622 **/ 623 int CeedOperatorGetQFunction(CeedOperator op, CeedQFunction *qf) { 624 bool is_composite; 625 Ceed ceed; 626 627 CeedCall(CeedOperatorGetCeed(op, &ceed)); 628 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 629 CeedCheck(!is_composite, ceed, CEED_ERROR_MINOR, "Not defined for composite operator"); 630 *qf = op->qf; 631 return CEED_ERROR_SUCCESS; 632 } 633 634 /** 635 @brief Get a boolean value indicating if the `CeedOperator` is composite 636 637 @param[in] op `CeedOperator` 638 @param[out] is_composite Variable to store composite status 639 640 @return An error code: 0 - success, otherwise - failure 641 642 @ref Backend 643 **/ 644 int CeedOperatorIsComposite(CeedOperator op, bool *is_composite) { 645 *is_composite = op->is_composite; 646 return CEED_ERROR_SUCCESS; 647 } 648 649 /** 650 @brief Get the backend data of a `CeedOperator` 651 652 @param[in] op `CeedOperator` 653 @param[out] data Variable to store data 654 655 @return An error code: 0 - success, otherwise - failure 656 657 @ref Backend 658 **/ 659 int CeedOperatorGetData(CeedOperator op, void *data) { 660 *(void **)data = op->data; 661 return CEED_ERROR_SUCCESS; 662 } 663 664 /** 665 @brief Set the backend data of a `CeedOperator` 666 667 @param[in,out] op `CeedOperator` 668 @param[in] data Data to set 669 670 @return An error code: 0 - success, otherwise - failure 671 672 @ref Backend 673 **/ 674 int CeedOperatorSetData(CeedOperator op, void *data) { 675 op->data = data; 676 return CEED_ERROR_SUCCESS; 677 } 678 679 /** 680 @brief Increment the reference counter for a `CeedOperator` 681 682 @param[in,out] op `CeedOperator` to increment the reference counter 683 684 @return An error code: 0 - success, otherwise - failure 685 686 @ref Backend 687 **/ 688 int CeedOperatorReference(CeedOperator op) { 689 op->ref_count++; 690 return CEED_ERROR_SUCCESS; 691 } 692 693 /** 694 @brief Set the setup flag of a `CeedOperator` to `true` 695 696 @param[in,out] op `CeedOperator` 697 698 @return An error code: 0 - success, otherwise - failure 699 700 @ref Backend 701 **/ 702 int CeedOperatorSetSetupDone(CeedOperator op) { 703 op->is_backend_setup = true; 704 return CEED_ERROR_SUCCESS; 705 } 706 707 /// @} 708 709 /// ---------------------------------------------------------------------------- 710 /// CeedOperator Public API 711 /// ---------------------------------------------------------------------------- 712 /// @addtogroup CeedOperatorUser 713 /// @{ 714 715 /** 716 @brief Create a `CeedOperator` and associate a `CeedQFunction`. 717 718 A `CeedBasis` and `CeedElemRestriction` can be associated with `CeedQFunction` fields with @ref CeedOperatorSetField(). 719 720 @param[in] ceed `Ceed` object used to create the `CeedOperator` 721 @param[in] qf `CeedQFunction` defining the action of the operator at quadrature points 722 @param[in] dqf `CeedQFunction` defining the action of the Jacobian of `qf` (or @ref CEED_QFUNCTION_NONE) 723 @param[in] dqfT `CeedQFunction` defining the action of the transpose of the Jacobian of `qf` (or @ref CEED_QFUNCTION_NONE) 724 @param[out] op Address of the variable where the newly created `CeedOperator` will be stored 725 726 @return An error code: 0 - success, otherwise - failure 727 728 @ref User 729 */ 730 int CeedOperatorCreate(Ceed ceed, CeedQFunction qf, CeedQFunction dqf, CeedQFunction dqfT, CeedOperator *op) { 731 if (!ceed->OperatorCreate) { 732 Ceed delegate; 733 734 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Operator")); 735 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support CeedOperatorCreate"); 736 CeedCall(CeedOperatorCreate(delegate, qf, dqf, dqfT, op)); 737 return CEED_ERROR_SUCCESS; 738 } 739 740 CeedCheck(qf && qf != CEED_QFUNCTION_NONE, ceed, CEED_ERROR_MINOR, "Operator must have a valid CeedQFunction."); 741 742 CeedCall(CeedCalloc(1, op)); 743 CeedCall(CeedReferenceCopy(ceed, &(*op)->ceed)); 744 (*op)->ref_count = 1; 745 (*op)->input_size = -1; 746 (*op)->output_size = -1; 747 CeedCall(CeedQFunctionReferenceCopy(qf, &(*op)->qf)); 748 if (dqf && dqf != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionReferenceCopy(dqf, &(*op)->dqf)); 749 if (dqfT && dqfT != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionReferenceCopy(dqfT, &(*op)->dqfT)); 750 CeedCall(CeedQFunctionAssemblyDataCreate(ceed, &(*op)->qf_assembled)); 751 CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*op)->input_fields)); 752 CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*op)->output_fields)); 753 CeedCall(ceed->OperatorCreate(*op)); 754 return CEED_ERROR_SUCCESS; 755 } 756 757 /** 758 @brief Create a `CeedOperator` for evaluation at evaluation at arbitrary points in each element. 759 760 A `CeedBasis` and `CeedElemRestriction` can be associated with `CeedQFunction` fields with `CeedOperator` SetField. 761 The locations of each point are set with @ref CeedOperatorAtPointsSetPoints(). 762 763 @param[in] ceed `Ceed` object used to create the `CeedOperator` 764 @param[in] qf `CeedQFunction` defining the action of the operator at quadrature points 765 @param[in] dqf `CeedQFunction` defining the action of the Jacobian of @a qf (or @ref CEED_QFUNCTION_NONE) 766 @param[in] dqfT `CeedQFunction` defining the action of the transpose of the Jacobian of @a qf (or @ref CEED_QFUNCTION_NONE) 767 @param[out] op Address of the variable where the newly created CeedOperator will be stored 768 769 @return An error code: 0 - success, otherwise - failure 770 771 @ref User 772 */ 773 int CeedOperatorCreateAtPoints(Ceed ceed, CeedQFunction qf, CeedQFunction dqf, CeedQFunction dqfT, CeedOperator *op) { 774 if (!ceed->OperatorCreateAtPoints) { 775 Ceed delegate; 776 777 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Operator")); 778 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support CeedOperatorCreateAtPoints"); 779 CeedCall(CeedOperatorCreateAtPoints(delegate, qf, dqf, dqfT, op)); 780 return CEED_ERROR_SUCCESS; 781 } 782 783 CeedCheck(qf && qf != CEED_QFUNCTION_NONE, ceed, CEED_ERROR_MINOR, "Operator must have a valid CeedQFunction."); 784 785 CeedCall(CeedCalloc(1, op)); 786 CeedCall(CeedReferenceCopy(ceed, &(*op)->ceed)); 787 (*op)->ref_count = 1; 788 (*op)->is_at_points = true; 789 (*op)->input_size = -1; 790 (*op)->output_size = -1; 791 CeedCall(CeedQFunctionReferenceCopy(qf, &(*op)->qf)); 792 if (dqf && dqf != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionReferenceCopy(dqf, &(*op)->dqf)); 793 if (dqfT && dqfT != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionReferenceCopy(dqfT, &(*op)->dqfT)); 794 CeedCall(CeedQFunctionAssemblyDataCreate(ceed, &(*op)->qf_assembled)); 795 CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*op)->input_fields)); 796 CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*op)->output_fields)); 797 CeedCall(ceed->OperatorCreateAtPoints(*op)); 798 return CEED_ERROR_SUCCESS; 799 } 800 801 /** 802 @brief Create a composite `CeedOperator` that composes the action of several `CeedOperator` 803 804 @param[in] ceed `Ceed` object used to create the `CeedOperator` 805 @param[out] op Address of the variable where the newly created composite `CeedOperator` will be stored 806 807 @return An error code: 0 - success, otherwise - failure 808 809 @ref User 810 */ 811 int CeedCompositeOperatorCreate(Ceed ceed, CeedOperator *op) { 812 if (!ceed->CompositeOperatorCreate) { 813 Ceed delegate; 814 815 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Operator")); 816 if (delegate) { 817 CeedCall(CeedCompositeOperatorCreate(delegate, op)); 818 return CEED_ERROR_SUCCESS; 819 } 820 } 821 822 CeedCall(CeedCalloc(1, op)); 823 CeedCall(CeedReferenceCopy(ceed, &(*op)->ceed)); 824 (*op)->ref_count = 1; 825 (*op)->is_composite = true; 826 CeedCall(CeedCalloc(CEED_COMPOSITE_MAX, &(*op)->sub_operators)); 827 (*op)->input_size = -1; 828 (*op)->output_size = -1; 829 830 if (ceed->CompositeOperatorCreate) CeedCall(ceed->CompositeOperatorCreate(*op)); 831 return CEED_ERROR_SUCCESS; 832 } 833 834 /** 835 @brief Copy the pointer to a `CeedOperator`. 836 837 Both pointers should be destroyed with @ref CeedOperatorDestroy(). 838 839 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`. 840 This `CeedOperator` will be destroyed if `*op_copy` is the only reference to this `CeedOperator`. 841 842 @param[in] op `CeedOperator` to copy reference to 843 @param[in,out] op_copy Variable to store copied reference 844 845 @return An error code: 0 - success, otherwise - failure 846 847 @ref User 848 **/ 849 int CeedOperatorReferenceCopy(CeedOperator op, CeedOperator *op_copy) { 850 CeedCall(CeedOperatorReference(op)); 851 CeedCall(CeedOperatorDestroy(op_copy)); 852 *op_copy = op; 853 return CEED_ERROR_SUCCESS; 854 } 855 856 /** 857 @brief Provide a field to a `CeedOperator` for use by its `CeedQFunction`. 858 859 This function is used to specify both active and passive fields to a `CeedOperator`. 860 For passive fields, a `CeedVector` `vec` must be provided. 861 Passive fields can inputs or outputs (updated in-place when operator is applied). 862 863 Active fields must be specified using this function, but their data (in a `CeedVector`) is passed in @ref CeedOperatorApply(). 864 There can be at most one active input `CeedVector` and at most one active output@ref CeedVector passed to @ref CeedOperatorApply(). 865 866 The number of quadrature points must agree across all points. 867 When using @ref CEED_BASIS_NONE, the number of quadrature points is determined by the element size of `rstr`. 868 869 @param[in,out] op `CeedOperator` on which to provide the field 870 @param[in] field_name Name of the field (to be matched with the name used by `CeedQFunction`) 871 @param[in] rstr `CeedElemRestriction` 872 @param[in] basis `CeedBasis` in which the field resides or @ref CEED_BASIS_NONE if collocated with quadrature points 873 @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` 874 875 @return An error code: 0 - success, otherwise - failure 876 877 @ref User 878 **/ 879 int CeedOperatorSetField(CeedOperator op, const char *field_name, CeedElemRestriction rstr, CeedBasis basis, CeedVector vec) { 880 bool is_input = true, is_at_points, is_composite, is_immutable; 881 CeedInt num_elem = 0, num_qpts = 0, num_input_fields, num_output_fields; 882 Ceed ceed; 883 CeedQFunction qf; 884 CeedQFunctionField qf_field, *qf_input_fields, *qf_output_fields; 885 CeedOperatorField *op_field; 886 887 CeedCall(CeedOperatorGetCeed(op, &ceed)); 888 CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); 889 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 890 CeedCall(CeedOperatorIsImmutable(op, &is_immutable)); 891 CeedCheck(!is_composite, ceed, CEED_ERROR_INCOMPATIBLE, "Cannot add field to composite operator."); 892 CeedCheck(!is_immutable, ceed, CEED_ERROR_MAJOR, "Operator cannot be changed after set as immutable"); 893 CeedCheck(rstr, ceed, CEED_ERROR_INCOMPATIBLE, "CeedElemRestriction rstr for field \"%s\" must be non-NULL.", field_name); 894 CeedCheck(basis, ceed, CEED_ERROR_INCOMPATIBLE, "CeedBasis basis for field \"%s\" must be non-NULL.", field_name); 895 CeedCheck(vec, ceed, CEED_ERROR_INCOMPATIBLE, "CeedVector vec for field \"%s\" must be non-NULL.", field_name); 896 897 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 898 CeedCheck(rstr == CEED_ELEMRESTRICTION_NONE || !op->has_restriction || num_elem == op->num_elem, ceed, CEED_ERROR_DIMENSION, 899 "CeedElemRestriction with %" CeedInt_FMT " elements incompatible with prior %" CeedInt_FMT " elements", num_elem, op->num_elem); 900 { 901 CeedRestrictionType rstr_type; 902 903 CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type)); 904 if (rstr_type == CEED_RESTRICTION_POINTS) { 905 CeedCheck(is_at_points, ceed, CEED_ERROR_UNSUPPORTED, "CeedElemRestriction AtPoints not supported for standard operator fields"); 906 CeedCheck(basis == CEED_BASIS_NONE, ceed, CEED_ERROR_UNSUPPORTED, "CeedElemRestriction AtPoints must be used with CEED_BASIS_NONE"); 907 if (!op->first_points_rstr) { 908 CeedCall(CeedElemRestrictionReferenceCopy(rstr, &op->first_points_rstr)); 909 } else { 910 bool are_compatible; 911 912 CeedCall(CeedElemRestrictionAtPointsAreCompatible(op->first_points_rstr, rstr, &are_compatible)); 913 CeedCheck(are_compatible, ceed, CEED_ERROR_INCOMPATIBLE, 914 "CeedElemRestriction must have compatible offsets with previously set CeedElemRestriction"); 915 } 916 } 917 } 918 919 if (basis == CEED_BASIS_NONE) CeedCall(CeedElemRestrictionGetElementSize(rstr, &num_qpts)); 920 else CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); 921 CeedCheck(op->num_qpts == 0 || num_qpts == op->num_qpts, ceed, CEED_ERROR_DIMENSION, 922 "%s must correspond to the same number of quadrature points as previously added CeedBases. Found %" CeedInt_FMT 923 " quadrature points but expected %" CeedInt_FMT " quadrature points.", 924 basis == CEED_BASIS_NONE ? "CeedElemRestriction" : "CeedBasis", num_qpts, op->num_qpts); 925 926 CeedCall(CeedOperatorGetQFunction(op, &qf)); 927 CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &qf_input_fields, &num_output_fields, &qf_output_fields)); 928 for (CeedInt i = 0; i < num_input_fields; i++) { 929 const char *qf_field_name; 930 931 CeedCall(CeedQFunctionFieldGetName(qf_input_fields[i], &qf_field_name)); 932 if (!strcmp(field_name, qf_field_name)) { 933 qf_field = qf_input_fields[i]; 934 op_field = &op->input_fields[i]; 935 goto found; 936 } 937 } 938 is_input = false; 939 for (CeedInt i = 0; i < num_output_fields; i++) { 940 const char *qf_field_name; 941 942 CeedCall(CeedQFunctionFieldGetName(qf_output_fields[i], &qf_field_name)); 943 if (!strcmp(field_name, qf_field_name)) { 944 qf_field = qf_output_fields[i]; 945 op_field = &op->output_fields[i]; 946 goto found; 947 } 948 } 949 // LCOV_EXCL_START 950 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "CeedQFunction has no knowledge of field '%s'", field_name); 951 // LCOV_EXCL_STOP 952 found: 953 CeedCall(CeedOperatorCheckField(ceed, qf_field, rstr, basis)); 954 CeedCall(CeedCalloc(1, op_field)); 955 956 if (vec == CEED_VECTOR_ACTIVE) { 957 CeedSize l_size; 958 959 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &l_size)); 960 if (is_input) { 961 if (op->input_size == -1) op->input_size = l_size; 962 CeedCheck(l_size == op->input_size, ceed, CEED_ERROR_INCOMPATIBLE, 963 "LVector size %" CeedSize_FMT " does not match previous size %" CeedSize_FMT "", l_size, op->input_size); 964 } else { 965 if (op->output_size == -1) op->output_size = l_size; 966 CeedCheck(l_size == op->output_size, ceed, CEED_ERROR_INCOMPATIBLE, 967 "LVector size %" CeedSize_FMT " does not match previous size %" CeedSize_FMT "", l_size, op->output_size); 968 } 969 } 970 971 CeedCall(CeedVectorReferenceCopy(vec, &(*op_field)->vec)); 972 CeedCall(CeedElemRestrictionReferenceCopy(rstr, &(*op_field)->elem_rstr)); 973 if (rstr != CEED_ELEMRESTRICTION_NONE && !op->has_restriction) { 974 op->num_elem = num_elem; 975 op->has_restriction = true; // Restriction set, but num_elem may be 0 976 } 977 CeedCall(CeedBasisReferenceCopy(basis, &(*op_field)->basis)); 978 if (op->num_qpts == 0 && !is_at_points) op->num_qpts = num_qpts; // no consistent number of qpts for OperatorAtPoints 979 op->num_fields += 1; 980 CeedCall(CeedStringAllocCopy(field_name, (char **)&(*op_field)->field_name)); 981 return CEED_ERROR_SUCCESS; 982 } 983 984 /** 985 @brief Get the `CeedOperator` Field of a `CeedOperator`. 986 987 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 988 989 @param[in] op `CeedOperator` 990 @param[out] num_input_fields Variable to store number of input fields 991 @param[out] input_fields Variable to store input fields 992 @param[out] num_output_fields Variable to store number of output fields 993 @param[out] output_fields Variable to store output fields 994 995 @return An error code: 0 - success, otherwise - failure 996 997 @ref Advanced 998 **/ 999 int CeedOperatorGetFields(CeedOperator op, CeedInt *num_input_fields, CeedOperatorField **input_fields, CeedInt *num_output_fields, 1000 CeedOperatorField **output_fields) { 1001 bool is_composite; 1002 Ceed ceed; 1003 CeedQFunction qf; 1004 1005 CeedCall(CeedOperatorGetCeed(op, &ceed)); 1006 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1007 CeedCheck(!is_composite, ceed, CEED_ERROR_MINOR, "Not defined for composite operator"); 1008 CeedCall(CeedOperatorCheckReady(op)); 1009 1010 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1011 CeedCall(CeedQFunctionGetFields(qf, num_input_fields, NULL, num_output_fields, NULL)); 1012 if (input_fields) *input_fields = op->input_fields; 1013 if (output_fields) *output_fields = op->output_fields; 1014 return CEED_ERROR_SUCCESS; 1015 } 1016 1017 /** 1018 @brief Set the arbitrary points in each element for a `CeedOperator` at points. 1019 1020 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 1021 1022 @param[in,out] op `CeedOperator` at points 1023 @param[in] rstr_points `CeedElemRestriction` for the coordinates of each point by element 1024 @param[in] point_coords `CeedVector` holding coordinates of each point 1025 1026 @return An error code: 0 - success, otherwise - failure 1027 1028 @ref Advanced 1029 **/ 1030 int CeedOperatorAtPointsSetPoints(CeedOperator op, CeedElemRestriction rstr_points, CeedVector point_coords) { 1031 bool is_at_points, is_immutable; 1032 Ceed ceed; 1033 1034 CeedCall(CeedOperatorGetCeed(op, &ceed)); 1035 CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); 1036 CeedCall(CeedOperatorIsImmutable(op, &is_immutable)); 1037 CeedCheck(is_at_points, ceed, CEED_ERROR_MINOR, "Only defined for operator at points"); 1038 CeedCheck(!is_immutable, ceed, CEED_ERROR_MAJOR, "Operator cannot be changed after set as immutable"); 1039 1040 if (!op->first_points_rstr) { 1041 CeedCall(CeedElemRestrictionReferenceCopy(rstr_points, &op->first_points_rstr)); 1042 } else { 1043 bool are_compatible; 1044 1045 CeedCall(CeedElemRestrictionAtPointsAreCompatible(op->first_points_rstr, rstr_points, &are_compatible)); 1046 CeedCheck(are_compatible, ceed, CEED_ERROR_INCOMPATIBLE, 1047 "CeedElemRestriction must have compatible offsets with previously set field CeedElemRestriction"); 1048 } 1049 1050 CeedCall(CeedElemRestrictionReferenceCopy(rstr_points, &op->rstr_points)); 1051 CeedCall(CeedVectorReferenceCopy(point_coords, &op->point_coords)); 1052 return CEED_ERROR_SUCCESS; 1053 } 1054 1055 /** 1056 @brief Get a boolean value indicating if the `CeedOperator` was created with `CeedOperatorCreateAtPoints` 1057 1058 @param[in] op `CeedOperator` 1059 @param[out] is_at_points Variable to store at points status 1060 1061 @return An error code: 0 - success, otherwise - failure 1062 1063 @ref User 1064 **/ 1065 int CeedOperatorIsAtPoints(CeedOperator op, bool *is_at_points) { 1066 *is_at_points = op->is_at_points; 1067 return CEED_ERROR_SUCCESS; 1068 } 1069 1070 /** 1071 @brief Get the arbitrary points in each element for a `CeedOperator` at points. 1072 1073 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 1074 1075 @param[in] op `CeedOperator` at points 1076 @param[out] rstr_points Variable to hold `CeedElemRestriction` for the coordinates of each point by element 1077 @param[out] point_coords Variable to hold `CeedVector` holding coordinates of each point 1078 1079 @return An error code: 0 - success, otherwise - failure 1080 1081 @ref Advanced 1082 **/ 1083 int CeedOperatorAtPointsGetPoints(CeedOperator op, CeedElemRestriction *rstr_points, CeedVector *point_coords) { 1084 bool is_at_points; 1085 Ceed ceed; 1086 1087 CeedCall(CeedOperatorGetCeed(op, &ceed)); 1088 CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); 1089 CeedCheck(is_at_points, ceed, CEED_ERROR_MINOR, "Only defined for operator at points"); 1090 CeedCall(CeedOperatorCheckReady(op)); 1091 1092 if (rstr_points) CeedCall(CeedElemRestrictionReferenceCopy(op->rstr_points, rstr_points)); 1093 if (point_coords) CeedCall(CeedVectorReferenceCopy(op->point_coords, point_coords)); 1094 return CEED_ERROR_SUCCESS; 1095 } 1096 1097 /** 1098 @brief Get a `CeedOperator` Field of a `CeedOperator` from its name. 1099 1100 `op_field` is set to `NULL` if the field is not found. 1101 1102 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 1103 1104 @param[in] op `CeedOperator` 1105 @param[in] field_name Name of desired `CeedOperator` Field 1106 @param[out] op_field `CeedOperator` Field corresponding to the name 1107 1108 @return An error code: 0 - success, otherwise - failure 1109 1110 @ref Advanced 1111 **/ 1112 int CeedOperatorGetFieldByName(CeedOperator op, const char *field_name, CeedOperatorField *op_field) { 1113 const char *name; 1114 CeedInt num_input_fields, num_output_fields; 1115 CeedOperatorField *input_fields, *output_fields; 1116 1117 *op_field = NULL; 1118 CeedCall(CeedOperatorGetFields(op, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); 1119 for (CeedInt i = 0; i < num_input_fields; i++) { 1120 CeedCall(CeedOperatorFieldGetName(input_fields[i], &name)); 1121 if (!strcmp(name, field_name)) { 1122 *op_field = input_fields[i]; 1123 return CEED_ERROR_SUCCESS; 1124 } 1125 } 1126 for (CeedInt i = 0; i < num_output_fields; i++) { 1127 CeedCall(CeedOperatorFieldGetName(output_fields[i], &name)); 1128 if (!strcmp(name, field_name)) { 1129 *op_field = output_fields[i]; 1130 return CEED_ERROR_SUCCESS; 1131 } 1132 } 1133 return CEED_ERROR_SUCCESS; 1134 } 1135 1136 /** 1137 @brief Get the name of a `CeedOperator` Field 1138 1139 @param[in] op_field `CeedOperator` Field 1140 @param[out] field_name Variable to store the field name 1141 1142 @return An error code: 0 - success, otherwise - failure 1143 1144 @ref Advanced 1145 **/ 1146 int CeedOperatorFieldGetName(CeedOperatorField op_field, const char **field_name) { 1147 *field_name = op_field->field_name; 1148 return CEED_ERROR_SUCCESS; 1149 } 1150 1151 /** 1152 @brief Get the `CeedElemRestriction` of a `CeedOperator` Field 1153 1154 @param[in] op_field `CeedOperator` Field 1155 @param[out] rstr Variable to store `CeedElemRestriction` 1156 1157 @return An error code: 0 - success, otherwise - failure 1158 1159 @ref Advanced 1160 **/ 1161 int CeedOperatorFieldGetElemRestriction(CeedOperatorField op_field, CeedElemRestriction *rstr) { 1162 *rstr = op_field->elem_rstr; 1163 return CEED_ERROR_SUCCESS; 1164 } 1165 1166 /** 1167 @brief Get the `CeedBasis` of a `CeedOperator` Field 1168 1169 @param[in] op_field `CeedOperator` Field 1170 @param[out] basis Variable to store `CeedBasis` 1171 1172 @return An error code: 0 - success, otherwise - failure 1173 1174 @ref Advanced 1175 **/ 1176 int CeedOperatorFieldGetBasis(CeedOperatorField op_field, CeedBasis *basis) { 1177 *basis = op_field->basis; 1178 return CEED_ERROR_SUCCESS; 1179 } 1180 1181 /** 1182 @brief Get the `CeedVector` of a `CeedOperator` Field 1183 1184 @param[in] op_field `CeedOperator` Field 1185 @param[out] vec Variable to store `CeedVector` 1186 1187 @return An error code: 0 - success, otherwise - failure 1188 1189 @ref Advanced 1190 **/ 1191 int CeedOperatorFieldGetVector(CeedOperatorField op_field, CeedVector *vec) { 1192 *vec = op_field->vec; 1193 return CEED_ERROR_SUCCESS; 1194 } 1195 1196 /** 1197 @brief Get the data of a `CeedOperator` Field. 1198 1199 Any arguments set as `NULL` are ignored. 1200 1201 @param[in] op_field `CeedOperator` Field 1202 @param[out] field_name Variable to store the field name 1203 @param[out] rstr Variable to store `CeedElemRestriction` 1204 @param[out] basis Variable to store `CeedBasis` 1205 @param[out] vec Variable to store `CeedVector` 1206 1207 @return An error code: 0 - success, otherwise - failure 1208 1209 @ref Advanced 1210 **/ 1211 int CeedOperatorFieldGetData(CeedOperatorField op_field, const char **field_name, CeedElemRestriction *rstr, CeedBasis *basis, CeedVector *vec) { 1212 if (field_name) CeedCall(CeedOperatorFieldGetName(op_field, field_name)); 1213 if (rstr) CeedCall(CeedOperatorFieldGetElemRestriction(op_field, rstr)); 1214 if (basis) CeedCall(CeedOperatorFieldGetBasis(op_field, basis)); 1215 if (vec) CeedCall(CeedOperatorFieldGetVector(op_field, vec)); 1216 return CEED_ERROR_SUCCESS; 1217 } 1218 1219 /** 1220 @brief Add a sub-operator to a composite `CeedOperator` 1221 1222 @param[in,out] composite_op Composite `CeedOperator` 1223 @param[in] sub_op Sub-operator `CeedOperator` 1224 1225 @return An error code: 0 - success, otherwise - failure 1226 1227 @ref User 1228 */ 1229 int CeedCompositeOperatorAddSub(CeedOperator composite_op, CeedOperator sub_op) { 1230 bool is_immutable; 1231 Ceed ceed; 1232 1233 CeedCall(CeedOperatorGetCeed(composite_op, &ceed)); 1234 CeedCheck(composite_op->is_composite, ceed, CEED_ERROR_MINOR, "CeedOperator is not a composite operator"); 1235 CeedCheck(composite_op->num_suboperators < CEED_COMPOSITE_MAX, ceed, CEED_ERROR_UNSUPPORTED, "Cannot add additional sub-operators"); 1236 CeedCall(CeedOperatorIsImmutable(composite_op, &is_immutable)); 1237 CeedCheck(!is_immutable, ceed, CEED_ERROR_MAJOR, "Operator cannot be changed after set as immutable"); 1238 1239 { 1240 CeedSize input_size, output_size; 1241 1242 CeedCall(CeedOperatorGetActiveVectorLengths(sub_op, &input_size, &output_size)); 1243 if (composite_op->input_size == -1) composite_op->input_size = input_size; 1244 if (composite_op->output_size == -1) composite_op->output_size = output_size; 1245 // Note, a size of -1 means no active vector restriction set, so no incompatibility 1246 CeedCheck((input_size == -1 || input_size == composite_op->input_size) && (output_size == -1 || output_size == composite_op->output_size), ceed, 1247 CEED_ERROR_MAJOR, 1248 "Sub-operators must have compatible dimensions; composite operator of shape (%" CeedSize_FMT ", %" CeedSize_FMT 1249 ") not compatible with sub-operator of " 1250 "shape (%" CeedSize_FMT ", %" CeedSize_FMT ")", 1251 composite_op->input_size, composite_op->output_size, input_size, output_size); 1252 } 1253 1254 composite_op->sub_operators[composite_op->num_suboperators] = sub_op; 1255 CeedCall(CeedOperatorReference(sub_op)); 1256 composite_op->num_suboperators++; 1257 return CEED_ERROR_SUCCESS; 1258 } 1259 1260 /** 1261 @brief Get the number of sub-operators associated with a `CeedOperator` 1262 1263 @param[in] op `CeedOperator` 1264 @param[out] num_suboperators Variable to store number of sub-operators 1265 1266 @return An error code: 0 - success, otherwise - failure 1267 1268 @ref Backend 1269 **/ 1270 int CeedCompositeOperatorGetNumSub(CeedOperator op, CeedInt *num_suboperators) { 1271 bool is_composite; 1272 Ceed ceed; 1273 1274 CeedCall(CeedOperatorGetCeed(op, &ceed)); 1275 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1276 CeedCheck(is_composite, ceed, CEED_ERROR_MINOR, "Only defined for a composite operator"); 1277 *num_suboperators = op->num_suboperators; 1278 return CEED_ERROR_SUCCESS; 1279 } 1280 1281 /** 1282 @brief Get the list of sub-operators associated with a `CeedOperator` 1283 1284 @param[in] op `CeedOperator` 1285 @param[out] sub_operators Variable to store list of sub-operators 1286 1287 @return An error code: 0 - success, otherwise - failure 1288 1289 @ref Backend 1290 **/ 1291 int CeedCompositeOperatorGetSubList(CeedOperator op, CeedOperator **sub_operators) { 1292 bool is_composite; 1293 Ceed ceed; 1294 1295 CeedCall(CeedOperatorGetCeed(op, &ceed)); 1296 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1297 CeedCheck(is_composite, ceed, CEED_ERROR_MINOR, "Only defined for a composite operator"); 1298 *sub_operators = op->sub_operators; 1299 return CEED_ERROR_SUCCESS; 1300 } 1301 1302 /** 1303 @brief Check if a `CeedOperator` is ready to be used. 1304 1305 @param[in] op `CeedOperator` to check 1306 1307 @return An error code: 0 - success, otherwise - failure 1308 1309 @ref User 1310 **/ 1311 int CeedOperatorCheckReady(CeedOperator op) { 1312 bool is_at_points, is_composite; 1313 Ceed ceed; 1314 CeedQFunction qf = NULL; 1315 1316 if (op->is_interface_setup) return CEED_ERROR_SUCCESS; 1317 1318 CeedCall(CeedOperatorGetCeed(op, &ceed)); 1319 CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); 1320 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1321 if (!is_composite) CeedCall(CeedOperatorGetQFunction(op, &qf)); 1322 if (is_composite) { 1323 CeedInt num_suboperators; 1324 1325 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1326 if (!num_suboperators) { 1327 // Empty operator setup 1328 op->input_size = 0; 1329 op->output_size = 0; 1330 } else { 1331 CeedOperator *sub_operators; 1332 1333 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1334 for (CeedInt i = 0; i < num_suboperators; i++) { 1335 CeedCall(CeedOperatorCheckReady(sub_operators[i])); 1336 } 1337 // Sub-operators could be modified after adding to composite operator 1338 // Need to verify no lvec incompatibility from any changes 1339 CeedSize input_size, output_size; 1340 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 1341 } 1342 } else { 1343 CeedInt num_input_fields, num_output_fields; 1344 1345 CeedCheck(op->num_fields > 0, ceed, CEED_ERROR_INCOMPLETE, "No operator fields set"); 1346 CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, NULL, &num_output_fields, NULL)); 1347 CeedCheck(op->num_fields == num_input_fields + num_output_fields, ceed, CEED_ERROR_INCOMPLETE, "Not all operator fields set"); 1348 CeedCheck(op->has_restriction, ceed, CEED_ERROR_INCOMPLETE, "At least one restriction required"); 1349 CeedCheck(op->num_qpts > 0 || is_at_points, ceed, CEED_ERROR_INCOMPLETE, 1350 "At least one non-collocated CeedBasis is required or the number of quadrature points must be set"); 1351 } 1352 1353 // Flag as immutable and ready 1354 op->is_interface_setup = true; 1355 if (qf && qf != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionSetImmutable(qf)); 1356 if (op->dqf && op->dqf != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionSetImmutable(op->dqf)); 1357 if (op->dqfT && op->dqfT != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionSetImmutable(op->dqfT)); 1358 return CEED_ERROR_SUCCESS; 1359 } 1360 1361 /** 1362 @brief Get vector lengths for the active input and/or output `CeedVector` of a `CeedOperator`. 1363 1364 Note: Lengths of `-1` indicate that the CeedOperator does not have an active input and/or output. 1365 1366 @param[in] op `CeedOperator` 1367 @param[out] input_size Variable to store active input vector length, or `NULL` 1368 @param[out] output_size Variable to store active output vector length, or `NULL` 1369 1370 @return An error code: 0 - success, otherwise - failure 1371 1372 @ref User 1373 **/ 1374 int CeedOperatorGetActiveVectorLengths(CeedOperator op, CeedSize *input_size, CeedSize *output_size) { 1375 bool is_composite; 1376 1377 if (input_size) *input_size = op->input_size; 1378 if (output_size) *output_size = op->output_size; 1379 1380 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1381 if (is_composite && (op->input_size == -1 || op->output_size == -1)) { 1382 CeedInt num_suboperators; 1383 CeedOperator *sub_operators; 1384 1385 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1386 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1387 for (CeedInt i = 0; i < num_suboperators; i++) { 1388 CeedSize sub_input_size, sub_output_size; 1389 1390 CeedCall(CeedOperatorGetActiveVectorLengths(sub_operators[i], &sub_input_size, &sub_output_size)); 1391 if (op->input_size == -1) op->input_size = sub_input_size; 1392 if (op->output_size == -1) op->output_size = sub_output_size; 1393 // Note, a size of -1 means no active vector restriction set, so no incompatibility 1394 CeedCheck((sub_input_size == -1 || sub_input_size == op->input_size) && (sub_output_size == -1 || sub_output_size == op->output_size), op->ceed, 1395 CEED_ERROR_MAJOR, 1396 "Sub-operators must have compatible dimensions; composite operator of shape (%" CeedSize_FMT ", %" CeedSize_FMT 1397 ") not compatible with sub-operator of " 1398 "shape (%" CeedSize_FMT ", %" CeedSize_FMT ")", 1399 op->input_size, op->output_size, input_size, output_size); 1400 } 1401 } 1402 return CEED_ERROR_SUCCESS; 1403 } 1404 1405 /** 1406 @brief Set reuse of `CeedQFunction` data in `CeedOperatorLinearAssemble*()` functions. 1407 1408 When `reuse_assembly_data = false` (default), the `CeedQFunction` associated with this `CeedOperator` is re-assembled every time a `CeedOperatorLinearAssemble*()` function is called. 1409 When `reuse_assembly_data = true`, the `CeedQFunction` associated with this `CeedOperator` is reused between calls to @ref CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(). 1410 1411 @param[in] op `CeedOperator` 1412 @param[in] reuse_assembly_data Boolean flag setting assembly data reuse 1413 1414 @return An error code: 0 - success, otherwise - failure 1415 1416 @ref Advanced 1417 **/ 1418 int CeedOperatorSetQFunctionAssemblyReuse(CeedOperator op, bool reuse_assembly_data) { 1419 bool is_composite; 1420 1421 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1422 if (is_composite) { 1423 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1424 CeedCall(CeedOperatorSetQFunctionAssemblyReuse(op->sub_operators[i], reuse_assembly_data)); 1425 } 1426 } else { 1427 CeedCall(CeedQFunctionAssemblyDataSetReuse(op->qf_assembled, reuse_assembly_data)); 1428 } 1429 return CEED_ERROR_SUCCESS; 1430 } 1431 1432 /** 1433 @brief Mark `CeedQFunction` data as updated and the `CeedQFunction` as requiring re-assembly. 1434 1435 @param[in] op `CeedOperator` 1436 @param[in] needs_data_update Boolean flag setting assembly data reuse 1437 1438 @return An error code: 0 - success, otherwise - failure 1439 1440 @ref Advanced 1441 **/ 1442 int CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(CeedOperator op, bool needs_data_update) { 1443 bool is_composite; 1444 1445 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1446 if (is_composite) { 1447 CeedInt num_suboperators; 1448 CeedOperator *sub_operators; 1449 1450 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1451 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1452 for (CeedInt i = 0; i < num_suboperators; i++) { 1453 CeedCall(CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(sub_operators[i], needs_data_update)); 1454 } 1455 } else { 1456 CeedCall(CeedQFunctionAssemblyDataSetUpdateNeeded(op->qf_assembled, needs_data_update)); 1457 } 1458 return CEED_ERROR_SUCCESS; 1459 } 1460 1461 /** 1462 @brief Set name of `CeedOperator` for @ref CeedOperatorView() output 1463 1464 @param[in,out] op `CeedOperator` 1465 @param[in] name Name to set, or NULL to remove previously set name 1466 1467 @return An error code: 0 - success, otherwise - failure 1468 1469 @ref User 1470 **/ 1471 int CeedOperatorSetName(CeedOperator op, const char *name) { 1472 char *name_copy; 1473 size_t name_len = name ? strlen(name) : 0; 1474 1475 CeedCall(CeedFree(&op->name)); 1476 if (name_len > 0) { 1477 CeedCall(CeedCalloc(name_len + 1, &name_copy)); 1478 memcpy(name_copy, name, name_len); 1479 op->name = name_copy; 1480 } 1481 return CEED_ERROR_SUCCESS; 1482 } 1483 1484 /** 1485 @brief View a `CeedOperator` 1486 1487 @param[in] op `CeedOperator` to view 1488 @param[in] stream Stream to write; typically `stdout` or a file 1489 1490 @return Error code: 0 - success, otherwise - failure 1491 1492 @ref User 1493 **/ 1494 int CeedOperatorView(CeedOperator op, FILE *stream) { 1495 bool has_name = op->name, is_composite; 1496 1497 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1498 if (is_composite) { 1499 CeedInt num_suboperators; 1500 CeedOperator *sub_operators; 1501 1502 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1503 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1504 fprintf(stream, "Composite CeedOperator%s%s\n", has_name ? " - " : "", has_name ? op->name : ""); 1505 1506 for (CeedInt i = 0; i < num_suboperators; i++) { 1507 has_name = sub_operators[i]->name; 1508 fprintf(stream, " SubOperator %" CeedInt_FMT "%s%s:\n", i, has_name ? " - " : "", has_name ? sub_operators[i]->name : ""); 1509 CeedCall(CeedOperatorSingleView(sub_operators[i], 1, stream)); 1510 } 1511 } else { 1512 fprintf(stream, "CeedOperator%s%s\n", has_name ? " - " : "", has_name ? op->name : ""); 1513 CeedCall(CeedOperatorSingleView(op, 0, stream)); 1514 } 1515 return CEED_ERROR_SUCCESS; 1516 } 1517 1518 /** 1519 @brief Get the `Ceed` associated with a `CeedOperator` 1520 1521 @param[in] op `CeedOperator` 1522 @param[out] ceed Variable to store `Ceed` 1523 1524 @return An error code: 0 - success, otherwise - failure 1525 1526 @ref Advanced 1527 **/ 1528 int CeedOperatorGetCeed(CeedOperator op, Ceed *ceed) { 1529 *ceed = op->ceed; 1530 return CEED_ERROR_SUCCESS; 1531 } 1532 1533 /** 1534 @brief Get the number of elements associated with a `CeedOperator` 1535 1536 @param[in] op `CeedOperator` 1537 @param[out] num_elem Variable to store number of elements 1538 1539 @return An error code: 0 - success, otherwise - failure 1540 1541 @ref Advanced 1542 **/ 1543 int CeedOperatorGetNumElements(CeedOperator op, CeedInt *num_elem) { 1544 bool is_composite; 1545 Ceed ceed; 1546 1547 CeedCall(CeedOperatorGetCeed(op, &ceed)); 1548 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1549 CeedCheck(!is_composite, ceed, CEED_ERROR_MINOR, "Not defined for composite operator"); 1550 *num_elem = op->num_elem; 1551 return CEED_ERROR_SUCCESS; 1552 } 1553 1554 /** 1555 @brief Get the number of quadrature points associated with a `CeedOperator` 1556 1557 @param[in] op `CeedOperator` 1558 @param[out] num_qpts Variable to store vector number of quadrature points 1559 1560 @return An error code: 0 - success, otherwise - failure 1561 1562 @ref Advanced 1563 **/ 1564 int CeedOperatorGetNumQuadraturePoints(CeedOperator op, CeedInt *num_qpts) { 1565 bool is_composite; 1566 Ceed ceed; 1567 1568 CeedCall(CeedOperatorGetCeed(op, &ceed)); 1569 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1570 CeedCheck(!is_composite, ceed, CEED_ERROR_MINOR, "Not defined for composite operator"); 1571 *num_qpts = op->num_qpts; 1572 return CEED_ERROR_SUCCESS; 1573 } 1574 1575 /** 1576 @brief Estimate number of FLOPs required to apply `CeedOperator` on the active `CeedVector` 1577 1578 @param[in] op `CeedOperator` to estimate FLOPs for 1579 @param[out] flops Address of variable to hold FLOPs estimate 1580 1581 @ref Backend 1582 **/ 1583 int CeedOperatorGetFlopsEstimate(CeedOperator op, CeedSize *flops) { 1584 bool is_composite; 1585 1586 CeedCall(CeedOperatorCheckReady(op)); 1587 1588 *flops = 0; 1589 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1590 if (is_composite) { 1591 CeedInt num_suboperators; 1592 1593 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1594 CeedOperator *sub_operators; 1595 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1596 1597 // FLOPs for each suboperator 1598 for (CeedInt i = 0; i < num_suboperators; i++) { 1599 CeedSize suboperator_flops; 1600 1601 CeedCall(CeedOperatorGetFlopsEstimate(sub_operators[i], &suboperator_flops)); 1602 *flops += suboperator_flops; 1603 } 1604 } else { 1605 CeedInt num_input_fields, num_output_fields, num_elem = 0; 1606 CeedQFunction qf; 1607 CeedQFunctionField *qf_input_fields, *qf_output_fields; 1608 CeedOperatorField *op_input_fields, *op_output_fields; 1609 1610 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1611 CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &qf_input_fields, &num_output_fields, &qf_output_fields)); 1612 CeedCall(CeedOperatorGetFields(op, NULL, &op_input_fields, NULL, &op_output_fields)); 1613 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 1614 1615 // Input FLOPs 1616 for (CeedInt i = 0; i < num_input_fields; i++) { 1617 CeedVector vec; 1618 1619 CeedCall(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1620 if (vec == CEED_VECTOR_ACTIVE) { 1621 CeedEvalMode eval_mode; 1622 CeedSize rstr_flops, basis_flops; 1623 CeedElemRestriction rstr; 1624 CeedBasis basis; 1625 1626 CeedCall(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &rstr)); 1627 CeedCall(CeedElemRestrictionGetFlopsEstimate(rstr, CEED_NOTRANSPOSE, &rstr_flops)); 1628 *flops += rstr_flops; 1629 CeedCall(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 1630 CeedCall(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 1631 CeedCall(CeedBasisGetFlopsEstimate(basis, CEED_NOTRANSPOSE, eval_mode, &basis_flops)); 1632 *flops += basis_flops * num_elem; 1633 } 1634 } 1635 // QF FLOPs 1636 { 1637 CeedInt num_qpts; 1638 CeedSize qf_flops; 1639 Ceed ceed; 1640 CeedQFunction qf; 1641 1642 CeedCall(CeedOperatorGetNumQuadraturePoints(op, &num_qpts)); 1643 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1644 CeedCall(CeedQFunctionGetFlopsEstimate(qf, &qf_flops)); 1645 CeedCall(CeedOperatorGetCeed(op, &ceed)); 1646 CeedCheck(qf_flops > -1, ceed, CEED_ERROR_INCOMPLETE, "Must set CeedQFunction FLOPs estimate with CeedQFunctionSetUserFlopsEstimate"); 1647 *flops += num_elem * num_qpts * qf_flops; 1648 } 1649 1650 // Output FLOPs 1651 for (CeedInt i = 0; i < num_output_fields; i++) { 1652 CeedVector vec; 1653 1654 CeedCall(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 1655 if (vec == CEED_VECTOR_ACTIVE) { 1656 CeedEvalMode eval_mode; 1657 CeedSize rstr_flops, basis_flops; 1658 CeedElemRestriction rstr; 1659 CeedBasis basis; 1660 1661 CeedCall(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &rstr)); 1662 CeedCall(CeedElemRestrictionGetFlopsEstimate(rstr, CEED_TRANSPOSE, &rstr_flops)); 1663 *flops += rstr_flops; 1664 CeedCall(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 1665 CeedCall(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 1666 CeedCall(CeedBasisGetFlopsEstimate(basis, CEED_TRANSPOSE, eval_mode, &basis_flops)); 1667 *flops += basis_flops * num_elem; 1668 } 1669 } 1670 } 1671 return CEED_ERROR_SUCCESS; 1672 } 1673 1674 /** 1675 @brief Get `CeedQFunction` global context for a `CeedOperator`. 1676 1677 The caller is responsible for destroying `ctx` returned from this function via @ref CeedQFunctionContextDestroy(). 1678 1679 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`. 1680 This `CeedQFunctionContext` will be destroyed if `ctx` is the only reference to this `CeedQFunctionContext`. 1681 1682 @param[in] op `CeedOperator` 1683 @param[out] ctx Variable to store `CeedQFunctionContext` 1684 1685 @return An error code: 0 - success, otherwise - failure 1686 1687 @ref Advanced 1688 **/ 1689 int CeedOperatorGetContext(CeedOperator op, CeedQFunctionContext *ctx) { 1690 bool is_composite; 1691 Ceed ceed; 1692 CeedQFunction qf; 1693 CeedQFunctionContext qf_ctx; 1694 1695 CeedCall(CeedOperatorGetCeed(op, &ceed)); 1696 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1697 CeedCheck(!is_composite, ceed, CEED_ERROR_INCOMPATIBLE, "Cannot retrieve CeedQFunctionContext for composite operator"); 1698 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1699 CeedCall(CeedQFunctionGetInnerContext(qf, &qf_ctx)); 1700 if (qf_ctx) CeedCall(CeedQFunctionContextReferenceCopy(qf_ctx, ctx)); 1701 return CEED_ERROR_SUCCESS; 1702 } 1703 1704 /** 1705 @brief Get label for a registered `CeedQFunctionContext` field, or `NULL` if no field has been registered with this `field_name`. 1706 1707 Fields are registered via `CeedQFunctionContextRegister*()` functions (eg. @ref CeedQFunctionContextRegisterDouble()). 1708 1709 @param[in] op `CeedOperator` 1710 @param[in] field_name Name of field to retrieve label 1711 @param[out] field_label Variable to field label 1712 1713 @return An error code: 0 - success, otherwise - failure 1714 1715 @ref User 1716 **/ 1717 int CeedOperatorGetContextFieldLabel(CeedOperator op, const char *field_name, CeedContextFieldLabel *field_label) { 1718 bool is_composite, field_found = false; 1719 1720 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1721 1722 if (is_composite) { 1723 // Composite operator 1724 // -- Check if composite label already created 1725 for (CeedInt i = 0; i < op->num_context_labels; i++) { 1726 if (!strcmp(op->context_labels[i]->name, field_name)) { 1727 *field_label = op->context_labels[i]; 1728 return CEED_ERROR_SUCCESS; 1729 } 1730 } 1731 1732 // -- Create composite label if needed 1733 CeedInt num_sub; 1734 CeedOperator *sub_operators; 1735 CeedContextFieldLabel new_field_label; 1736 1737 CeedCall(CeedCalloc(1, &new_field_label)); 1738 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub)); 1739 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1740 CeedCall(CeedCalloc(num_sub, &new_field_label->sub_labels)); 1741 new_field_label->num_sub_labels = num_sub; 1742 1743 for (CeedInt i = 0; i < num_sub; i++) { 1744 if (sub_operators[i]->qf->ctx) { 1745 CeedContextFieldLabel new_field_label_i; 1746 1747 CeedCall(CeedQFunctionContextGetFieldLabel(sub_operators[i]->qf->ctx, field_name, &new_field_label_i)); 1748 if (new_field_label_i) { 1749 field_found = true; 1750 new_field_label->sub_labels[i] = new_field_label_i; 1751 new_field_label->name = new_field_label_i->name; 1752 new_field_label->description = new_field_label_i->description; 1753 if (new_field_label->type && new_field_label->type != new_field_label_i->type) { 1754 // LCOV_EXCL_START 1755 CeedCall(CeedFree(&new_field_label)); 1756 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, "Incompatible field types on sub-operator contexts. %s != %s", 1757 CeedContextFieldTypes[new_field_label->type], CeedContextFieldTypes[new_field_label_i->type]); 1758 // LCOV_EXCL_STOP 1759 } else { 1760 new_field_label->type = new_field_label_i->type; 1761 } 1762 if (new_field_label->num_values != 0 && new_field_label->num_values != new_field_label_i->num_values) { 1763 // LCOV_EXCL_START 1764 Ceed ceed; 1765 1766 CeedCall(CeedFree(&new_field_label)); 1767 CeedCall(CeedOperatorGetCeed(op, &ceed)); 1768 return CeedError(ceed, CEED_ERROR_INCOMPATIBLE, "Incompatible field number of values on sub-operator contexts. %zu != %zu", 1769 new_field_label->num_values, new_field_label_i->num_values); 1770 // LCOV_EXCL_STOP 1771 } else { 1772 new_field_label->num_values = new_field_label_i->num_values; 1773 } 1774 } 1775 } 1776 } 1777 // -- Cleanup if field was found 1778 if (field_found) { 1779 *field_label = new_field_label; 1780 } else { 1781 // LCOV_EXCL_START 1782 CeedCall(CeedFree(&new_field_label->sub_labels)); 1783 CeedCall(CeedFree(&new_field_label)); 1784 *field_label = NULL; 1785 // LCOV_EXCL_STOP 1786 } 1787 } else { 1788 CeedQFunction qf; 1789 CeedQFunctionContext ctx; 1790 1791 // Single, non-composite operator 1792 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1793 CeedCall(CeedQFunctionGetInnerContext(qf, &ctx)); 1794 if (ctx) { 1795 CeedCall(CeedQFunctionContextGetFieldLabel(ctx, field_name, field_label)); 1796 } else { 1797 *field_label = NULL; 1798 } 1799 } 1800 1801 // Set label in operator 1802 if (*field_label) { 1803 (*field_label)->from_op = true; 1804 1805 // Move new composite label to operator 1806 if (op->num_context_labels == 0) { 1807 CeedCall(CeedCalloc(1, &op->context_labels)); 1808 op->max_context_labels = 1; 1809 } else if (op->num_context_labels == op->max_context_labels) { 1810 CeedCall(CeedRealloc(2 * op->num_context_labels, &op->context_labels)); 1811 op->max_context_labels *= 2; 1812 } 1813 op->context_labels[op->num_context_labels] = *field_label; 1814 op->num_context_labels++; 1815 } 1816 return CEED_ERROR_SUCCESS; 1817 } 1818 1819 /** 1820 @brief Set `CeedQFunctionContext` field holding double precision values. 1821 1822 For composite operators, the values are set in all sub-operator `CeedQFunctionContext` that have a matching `field_name`. 1823 1824 @param[in,out] op `CeedOperator` 1825 @param[in] field_label Label of field to set 1826 @param[in] values Values to set 1827 1828 @return An error code: 0 - success, otherwise - failure 1829 1830 @ref User 1831 **/ 1832 int CeedOperatorSetContextDouble(CeedOperator op, CeedContextFieldLabel field_label, double *values) { 1833 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, values); 1834 } 1835 1836 /** 1837 @brief Get `CeedQFunctionContext` field holding double precision values, read-only. 1838 1839 For composite operators, the values correspond to the first sub-operator `CeedQFunctionContext` that has a matching `field_name`. 1840 1841 @param[in] op `CeedOperator` 1842 @param[in] field_label Label of field to get 1843 @param[out] num_values Number of values in the field label 1844 @param[out] values Pointer to context values 1845 1846 @return An error code: 0 - success, otherwise - failure 1847 1848 @ref User 1849 **/ 1850 int CeedOperatorGetContextDoubleRead(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const double **values) { 1851 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, num_values, values); 1852 } 1853 1854 /** 1855 @brief Restore `CeedQFunctionContext` field holding double precision values, read-only. 1856 1857 @param[in] op `CeedOperator` 1858 @param[in] field_label Label of field to restore 1859 @param[out] values Pointer to context values 1860 1861 @return An error code: 0 - success, otherwise - failure 1862 1863 @ref User 1864 **/ 1865 int CeedOperatorRestoreContextDoubleRead(CeedOperator op, CeedContextFieldLabel field_label, const double **values) { 1866 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, values); 1867 } 1868 1869 /** 1870 @brief Set `CeedQFunctionContext` field holding `int32` values. 1871 1872 For composite operators, the values are set in all sub-operator `CeedQFunctionContext` that have a matching `field_name`. 1873 1874 @param[in,out] op `CeedOperator` 1875 @param[in] field_label Label of field to set 1876 @param[in] values Values to set 1877 1878 @return An error code: 0 - success, otherwise - failure 1879 1880 @ref User 1881 **/ 1882 int CeedOperatorSetContextInt32(CeedOperator op, CeedContextFieldLabel field_label, int32_t *values) { 1883 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_INT32, values); 1884 } 1885 1886 /** 1887 @brief Get `CeedQFunctionContext` field holding `int32` values, read-only. 1888 1889 For composite operators, the values correspond to the first sub-operator `CeedQFunctionContext` that has a matching `field_name`. 1890 1891 @param[in] op `CeedOperator` 1892 @param[in] field_label Label of field to get 1893 @param[out] num_values Number of `int32` values in `values` 1894 @param[out] values Pointer to context values 1895 1896 @return An error code: 0 - success, otherwise - failure 1897 1898 @ref User 1899 **/ 1900 int CeedOperatorGetContextInt32Read(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const int32_t **values) { 1901 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_INT32, num_values, values); 1902 } 1903 1904 /** 1905 @brief Restore `CeedQFunctionContext` field holding `int32` values, read-only. 1906 1907 @param[in] op `CeedOperator` 1908 @param[in] field_label Label of field to get 1909 @param[out] values Pointer to context values 1910 1911 @return An error code: 0 - success, otherwise - failure 1912 1913 @ref User 1914 **/ 1915 int CeedOperatorRestoreContextInt32Read(CeedOperator op, CeedContextFieldLabel field_label, const int32_t **values) { 1916 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_INT32, values); 1917 } 1918 1919 /** 1920 @brief Set `CeedQFunctionContext` field holding boolean values. 1921 1922 For composite operators, the values are set in all sub-operator `CeedQFunctionContext` that have a matching `field_name`. 1923 1924 @param[in,out] op `CeedOperator` 1925 @param[in] field_label Label of field to set 1926 @param[in] values Values to set 1927 1928 @return An error code: 0 - success, otherwise - failure 1929 1930 @ref User 1931 **/ 1932 int CeedOperatorSetContextBoolean(CeedOperator op, CeedContextFieldLabel field_label, bool *values) { 1933 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_BOOL, values); 1934 } 1935 1936 /** 1937 @brief Get `CeedQFunctionContext` field holding boolean values, read-only. 1938 1939 For composite operators, the values correspond to the first sub-operator `CeedQFunctionContext` that has a matching `field_name`. 1940 1941 @param[in] op `CeedOperator` 1942 @param[in] field_label Label of field to get 1943 @param[out] num_values Number of boolean values in `values` 1944 @param[out] values Pointer to context values 1945 1946 @return An error code: 0 - success, otherwise - failure 1947 1948 @ref User 1949 **/ 1950 int CeedOperatorGetContextBooleanRead(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const bool **values) { 1951 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_BOOL, num_values, values); 1952 } 1953 1954 /** 1955 @brief Restore `CeedQFunctionContext` field holding boolean values, read-only. 1956 1957 @param[in] op `CeedOperator` 1958 @param[in] field_label Label of field to get 1959 @param[out] values Pointer to context values 1960 1961 @return An error code: 0 - success, otherwise - failure 1962 1963 @ref User 1964 **/ 1965 int CeedOperatorRestoreContextBooleanRead(CeedOperator op, CeedContextFieldLabel field_label, const bool **values) { 1966 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_BOOL, values); 1967 } 1968 1969 /** 1970 @brief Apply `CeedOperator` to a `CeedVector`. 1971 1972 This computes the action of the operator on the specified (active) input, yielding its (active) output. 1973 All inputs and outputs must be specified using @ref CeedOperatorSetField(). 1974 1975 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 1976 1977 @param[in] op `CeedOperator` to apply 1978 @param[in] in `CeedVector` containing input state or @ref CEED_VECTOR_NONE if there are no active inputs 1979 @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 1980 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 1981 1982 @return An error code: 0 - success, otherwise - failure 1983 1984 @ref User 1985 **/ 1986 int CeedOperatorApply(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) { 1987 bool is_composite; 1988 1989 CeedCall(CeedOperatorCheckReady(op)); 1990 1991 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1992 if (is_composite) { 1993 // Composite Operator 1994 if (op->ApplyComposite) { 1995 CeedCall(op->ApplyComposite(op, in, out, request)); 1996 } else { 1997 CeedInt num_suboperators; 1998 CeedOperator *sub_operators; 1999 2000 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 2001 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 2002 2003 // Zero all output vectors 2004 if (out != CEED_VECTOR_NONE) CeedCall(CeedVectorSetValue(out, 0.0)); 2005 for (CeedInt i = 0; i < num_suboperators; i++) { 2006 CeedInt num_output_fields; 2007 CeedOperatorField *output_fields; 2008 2009 CeedCall(CeedOperatorGetFields(sub_operators[i], NULL, NULL, &num_output_fields, &output_fields)); 2010 for (CeedInt j = 0; j < num_output_fields; j++) { 2011 CeedVector vec; 2012 2013 CeedCall(CeedOperatorFieldGetVector(output_fields[j], &vec)); 2014 if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) { 2015 CeedCall(CeedVectorSetValue(vec, 0.0)); 2016 } 2017 } 2018 } 2019 // Apply 2020 for (CeedInt i = 0; i < num_suboperators; i++) { 2021 CeedCall(CeedOperatorApplyAdd(sub_operators[i], in, out, request)); 2022 } 2023 } 2024 } else { 2025 // Standard Operator 2026 if (op->Apply) { 2027 CeedCall(op->Apply(op, in, out, request)); 2028 } else { 2029 CeedInt num_output_fields; 2030 CeedOperatorField *output_fields; 2031 2032 CeedCall(CeedOperatorGetFields(op, NULL, NULL, &num_output_fields, &output_fields)); 2033 // Zero all output vectors 2034 for (CeedInt i = 0; i < num_output_fields; i++) { 2035 CeedVector vec; 2036 2037 CeedCall(CeedOperatorFieldGetVector(output_fields[i], &vec)); 2038 if (vec == CEED_VECTOR_ACTIVE) vec = out; 2039 if (vec != CEED_VECTOR_NONE) CeedCall(CeedVectorSetValue(vec, 0.0)); 2040 } 2041 // Apply 2042 if (op->num_elem > 0) CeedCall(op->ApplyAdd(op, in, out, request)); 2043 } 2044 } 2045 return CEED_ERROR_SUCCESS; 2046 } 2047 2048 /** 2049 @brief Apply `CeedOperator` to a `CeedVector` and add result to output `CeedVector`. 2050 2051 This computes the action of the operator on the specified (active) input, yielding its (active) output. 2052 All inputs and outputs must be specified using @ref CeedOperatorSetField(). 2053 2054 @param[in] op `CeedOperator` to apply 2055 @param[in] in `CeedVector` containing input state or @ref CEED_VECTOR_NONE if there are no active inputs 2056 @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 2057 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2058 2059 @return An error code: 0 - success, otherwise - failure 2060 2061 @ref User 2062 **/ 2063 int CeedOperatorApplyAdd(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) { 2064 bool is_composite; 2065 2066 CeedCall(CeedOperatorCheckReady(op)); 2067 2068 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2069 if (is_composite) { 2070 // Composite Operator 2071 if (op->ApplyAddComposite) { 2072 CeedCall(op->ApplyAddComposite(op, in, out, request)); 2073 } else { 2074 CeedInt num_suboperators; 2075 CeedOperator *sub_operators; 2076 2077 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 2078 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 2079 for (CeedInt i = 0; i < num_suboperators; i++) { 2080 CeedCall(CeedOperatorApplyAdd(sub_operators[i], in, out, request)); 2081 } 2082 } 2083 } else if (op->num_elem > 0) { 2084 // Standard Operator 2085 CeedCall(op->ApplyAdd(op, in, out, request)); 2086 } 2087 return CEED_ERROR_SUCCESS; 2088 } 2089 2090 /** 2091 @brief Destroy a `CeedOperator` 2092 2093 @param[in,out] op `CeedOperator` to destroy 2094 2095 @return An error code: 0 - success, otherwise - failure 2096 2097 @ref User 2098 **/ 2099 int CeedOperatorDestroy(CeedOperator *op) { 2100 if (!*op || --(*op)->ref_count > 0) { 2101 *op = NULL; 2102 return CEED_ERROR_SUCCESS; 2103 } 2104 if ((*op)->Destroy) CeedCall((*op)->Destroy(*op)); 2105 CeedCall(CeedDestroy(&(*op)->ceed)); 2106 // Free fields 2107 for (CeedInt i = 0; i < (*op)->num_fields; i++) { 2108 if ((*op)->input_fields[i]) { 2109 if ((*op)->input_fields[i]->elem_rstr != CEED_ELEMRESTRICTION_NONE) { 2110 CeedCall(CeedElemRestrictionDestroy(&(*op)->input_fields[i]->elem_rstr)); 2111 } 2112 if ((*op)->input_fields[i]->basis != CEED_BASIS_NONE) { 2113 CeedCall(CeedBasisDestroy(&(*op)->input_fields[i]->basis)); 2114 } 2115 if ((*op)->input_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->input_fields[i]->vec != CEED_VECTOR_NONE) { 2116 CeedCall(CeedVectorDestroy(&(*op)->input_fields[i]->vec)); 2117 } 2118 CeedCall(CeedFree(&(*op)->input_fields[i]->field_name)); 2119 CeedCall(CeedFree(&(*op)->input_fields[i])); 2120 } 2121 } 2122 for (CeedInt i = 0; i < (*op)->num_fields; i++) { 2123 if ((*op)->output_fields[i]) { 2124 CeedCall(CeedElemRestrictionDestroy(&(*op)->output_fields[i]->elem_rstr)); 2125 if ((*op)->output_fields[i]->basis != CEED_BASIS_NONE) { 2126 CeedCall(CeedBasisDestroy(&(*op)->output_fields[i]->basis)); 2127 } 2128 if ((*op)->output_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->output_fields[i]->vec != CEED_VECTOR_NONE) { 2129 CeedCall(CeedVectorDestroy(&(*op)->output_fields[i]->vec)); 2130 } 2131 CeedCall(CeedFree(&(*op)->output_fields[i]->field_name)); 2132 CeedCall(CeedFree(&(*op)->output_fields[i])); 2133 } 2134 } 2135 // AtPoints data 2136 CeedCall(CeedVectorDestroy(&(*op)->point_coords)); 2137 CeedCall(CeedElemRestrictionDestroy(&(*op)->rstr_points)); 2138 CeedCall(CeedElemRestrictionDestroy(&(*op)->first_points_rstr)); 2139 // Destroy sub_operators 2140 for (CeedInt i = 0; i < (*op)->num_suboperators; i++) { 2141 if ((*op)->sub_operators[i]) { 2142 CeedCall(CeedOperatorDestroy(&(*op)->sub_operators[i])); 2143 } 2144 } 2145 CeedCall(CeedQFunctionDestroy(&(*op)->qf)); 2146 CeedCall(CeedQFunctionDestroy(&(*op)->dqf)); 2147 CeedCall(CeedQFunctionDestroy(&(*op)->dqfT)); 2148 // Destroy any composite labels 2149 if ((*op)->is_composite) { 2150 for (CeedInt i = 0; i < (*op)->num_context_labels; i++) { 2151 CeedCall(CeedFree(&(*op)->context_labels[i]->sub_labels)); 2152 CeedCall(CeedFree(&(*op)->context_labels[i])); 2153 } 2154 } 2155 CeedCall(CeedFree(&(*op)->context_labels)); 2156 2157 // Destroy fallback 2158 CeedCall(CeedOperatorDestroy(&(*op)->op_fallback)); 2159 2160 // Destroy assembly data 2161 CeedCall(CeedQFunctionAssemblyDataDestroy(&(*op)->qf_assembled)); 2162 CeedCall(CeedOperatorAssemblyDataDestroy(&(*op)->op_assembled)); 2163 2164 CeedCall(CeedFree(&(*op)->input_fields)); 2165 CeedCall(CeedFree(&(*op)->output_fields)); 2166 CeedCall(CeedFree(&(*op)->sub_operators)); 2167 CeedCall(CeedFree(&(*op)->name)); 2168 CeedCall(CeedFree(op)); 2169 return CEED_ERROR_SUCCESS; 2170 } 2171 2172 /// @} 2173