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