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