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