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