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 Backend 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 CeedInt num_tabs = 0; 1609 1610 CeedCall(CeedOperatorGetName(op, &name)); 1611 has_name = name ? strlen(name) : false; 1612 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1613 CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); 1614 // Set tabs 1615 CeedCall(CeedOperatorGetNumViewTabs(op, &num_tabs)); 1616 CeedCall(CeedCalloc(CEED_TAB_WIDTH * (num_tabs + is_composite) + 1, &tabs)); 1617 for (CeedInt i = 0; i < CEED_TAB_WIDTH * num_tabs; i++) tabs[i] = ' '; 1618 if (is_composite) { 1619 CeedInt num_suboperators; 1620 CeedOperator *sub_operators; 1621 1622 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 1623 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 1624 fprintf(stream, "%s", tabs); 1625 fprintf(stream, "Composite CeedOperator%s%s\n", has_name ? " - " : "", has_name ? name : ""); 1626 for (CeedInt i = 0; i < CEED_TAB_WIDTH; i++) tabs[CEED_TAB_WIDTH * num_tabs + i] = ' '; 1627 for (CeedInt i = 0; i < num_suboperators; i++) { 1628 has_name = sub_operators[i]->name; 1629 fprintf(stream, "%s", tabs); 1630 fprintf(stream, "SubOperator%s %" CeedInt_FMT "%s%s%s\n", is_at_points ? " AtPoints" : "", i, has_name ? " - " : "", 1631 has_name ? sub_operators[i]->name : "", is_full ? ":" : ""); 1632 if (is_full) CeedCall(CeedOperatorSingleView(sub_operators[i], tabs, stream)); 1633 } 1634 } else { 1635 fprintf(stream, "%s", tabs); 1636 fprintf(stream, "CeedOperator%s%s%s\n", is_at_points ? " AtPoints" : "", has_name ? " - " : "", has_name ? name : ""); 1637 if (is_full) CeedCall(CeedOperatorSingleView(op, tabs, stream)); 1638 } 1639 CeedCall(CeedFree(&tabs)); 1640 return CEED_ERROR_SUCCESS; 1641 } 1642 1643 /** 1644 @brief Set the number of tabs to indent for @ref CeedOperatorView() output 1645 1646 @param[in] op `CeedOperator` to set the number of view tabs 1647 @param[in] num_tabs Number of view tabs to set 1648 1649 @return Error code: 0 - success, otherwise - failure 1650 1651 @ref User 1652 **/ 1653 int CeedOperatorSetNumViewTabs(CeedOperator op, CeedInt num_tabs) { 1654 CeedCheck(num_tabs >= 0, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Number of view tabs must be non-negative"); 1655 op->num_tabs = num_tabs; 1656 return CEED_ERROR_SUCCESS; 1657 } 1658 1659 /** 1660 @brief View a `CeedOperator` 1661 1662 @param[in] op `CeedOperator` to view 1663 @param[in] stream Stream to write; typically `stdout` or a file 1664 1665 @return Error code: 0 - success, otherwise - failure 1666 1667 @ref User 1668 **/ 1669 int CeedOperatorView(CeedOperator op, FILE *stream) { 1670 CeedCall(CeedOperatorView_Core(op, stream, true)); 1671 return CEED_ERROR_SUCCESS; 1672 } 1673 1674 /** 1675 @brief View a brief summary `CeedOperator` 1676 1677 @param[in] op `CeedOperator` to view brief summary 1678 @param[in] stream Stream to write; typically `stdout` or a file 1679 1680 @return Error code: 0 - success, otherwise - failure 1681 1682 @ref User 1683 **/ 1684 int CeedOperatorViewTerse(CeedOperator op, FILE *stream) { 1685 CeedCall(CeedOperatorView_Core(op, stream, false)); 1686 return CEED_ERROR_SUCCESS; 1687 } 1688 1689 /** 1690 @brief Get the `Ceed` associated with a `CeedOperator` 1691 1692 @param[in] op `CeedOperator` 1693 @param[out] ceed Variable to store `Ceed` 1694 1695 @return An error code: 0 - success, otherwise - failure 1696 1697 @ref Advanced 1698 **/ 1699 int CeedOperatorGetCeed(CeedOperator op, Ceed *ceed) { 1700 *ceed = NULL; 1701 CeedCall(CeedReferenceCopy(CeedOperatorReturnCeed(op), ceed)); 1702 return CEED_ERROR_SUCCESS; 1703 } 1704 1705 /** 1706 @brief Return the `Ceed` associated with a `CeedOperator` 1707 1708 @param[in] op `CeedOperator` 1709 1710 @return `Ceed` associated with the `op` 1711 1712 @ref Advanced 1713 **/ 1714 Ceed CeedOperatorReturnCeed(CeedOperator op) { return op->ceed; } 1715 1716 /** 1717 @brief Get the number of elements associated with a `CeedOperator` 1718 1719 @param[in] op `CeedOperator` 1720 @param[out] num_elem Variable to store number of elements 1721 1722 @return An error code: 0 - success, otherwise - failure 1723 1724 @ref Advanced 1725 **/ 1726 int CeedOperatorGetNumElements(CeedOperator op, CeedInt *num_elem) { 1727 bool is_composite; 1728 1729 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1730 CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Not defined for composite operator"); 1731 *num_elem = op->num_elem; 1732 return CEED_ERROR_SUCCESS; 1733 } 1734 1735 /** 1736 @brief Get the number of quadrature points associated with a `CeedOperator` 1737 1738 @param[in] op `CeedOperator` 1739 @param[out] num_qpts Variable to store vector number of quadrature points 1740 1741 @return An error code: 0 - success, otherwise - failure 1742 1743 @ref Advanced 1744 **/ 1745 int CeedOperatorGetNumQuadraturePoints(CeedOperator op, CeedInt *num_qpts) { 1746 bool is_composite; 1747 1748 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1749 CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Not defined for composite operator"); 1750 *num_qpts = op->num_qpts; 1751 return CEED_ERROR_SUCCESS; 1752 } 1753 1754 /** 1755 @brief Estimate number of FLOPs required to apply `CeedOperator` on the active `CeedVector` 1756 1757 @param[in] op `CeedOperator` to estimate FLOPs for 1758 @param[out] flops Address of variable to hold FLOPs estimate 1759 1760 @ref Backend 1761 **/ 1762 int CeedOperatorGetFlopsEstimate(CeedOperator op, CeedSize *flops) { 1763 bool is_composite; 1764 1765 CeedCall(CeedOperatorCheckReady(op)); 1766 1767 *flops = 0; 1768 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1769 if (is_composite) { 1770 CeedInt num_suboperators; 1771 1772 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 1773 CeedOperator *sub_operators; 1774 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 1775 1776 // FLOPs for each suboperator 1777 for (CeedInt i = 0; i < num_suboperators; i++) { 1778 CeedSize suboperator_flops; 1779 1780 CeedCall(CeedOperatorGetFlopsEstimate(sub_operators[i], &suboperator_flops)); 1781 *flops += suboperator_flops; 1782 } 1783 } else { 1784 bool is_at_points; 1785 CeedInt num_input_fields, num_output_fields, num_elem = 0, num_points = 0; 1786 CeedQFunction qf; 1787 CeedQFunctionField *qf_input_fields, *qf_output_fields; 1788 CeedOperatorField *op_input_fields, *op_output_fields; 1789 1790 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 1791 if (num_elem == 0) return CEED_ERROR_SUCCESS; 1792 CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); 1793 if (is_at_points) { 1794 CeedMemType mem_type; 1795 CeedElemRestriction rstr_points = NULL; 1796 1797 CeedCall(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL)); 1798 CeedCall(CeedGetPreferredMemType(CeedOperatorReturnCeed(op), &mem_type)); 1799 if (mem_type == CEED_MEM_DEVICE) { 1800 // Device backends pad out to the same number of points per element 1801 CeedCall(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &num_points)); 1802 } else { 1803 num_points = 0; 1804 for (CeedInt i = 0; i < num_elem; i++) { 1805 CeedInt points_in_elem = 0; 1806 1807 CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr_points, i, &points_in_elem)); 1808 num_points += points_in_elem; 1809 } 1810 num_points = num_points / num_elem + (num_points % num_elem > 0); 1811 } 1812 CeedCall(CeedElemRestrictionDestroy(&rstr_points)); 1813 } 1814 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1815 CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &qf_input_fields, &num_output_fields, &qf_output_fields)); 1816 CeedCall(CeedQFunctionDestroy(&qf)); 1817 CeedCall(CeedOperatorGetFields(op, NULL, &op_input_fields, NULL, &op_output_fields)); 1818 1819 // Input FLOPs 1820 for (CeedInt i = 0; i < num_input_fields; i++) { 1821 CeedVector vec; 1822 1823 CeedCall(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1824 if (vec == CEED_VECTOR_ACTIVE) { 1825 CeedEvalMode eval_mode; 1826 CeedSize rstr_flops, basis_flops; 1827 CeedElemRestriction rstr; 1828 CeedBasis basis; 1829 1830 CeedCall(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &rstr)); 1831 CeedCall(CeedElemRestrictionGetFlopsEstimate(rstr, CEED_NOTRANSPOSE, &rstr_flops)); 1832 CeedCall(CeedElemRestrictionDestroy(&rstr)); 1833 *flops += rstr_flops; 1834 CeedCall(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 1835 CeedCall(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 1836 CeedCall(CeedBasisGetFlopsEstimate(basis, CEED_NOTRANSPOSE, eval_mode, is_at_points, num_points, &basis_flops)); 1837 CeedCall(CeedBasisDestroy(&basis)); 1838 *flops += basis_flops * num_elem; 1839 } 1840 CeedCall(CeedVectorDestroy(&vec)); 1841 } 1842 // QF FLOPs 1843 { 1844 CeedInt num_qpts; 1845 CeedSize qf_flops; 1846 CeedQFunction qf; 1847 1848 if (is_at_points) num_qpts = num_points; 1849 else CeedCall(CeedOperatorGetNumQuadraturePoints(op, &num_qpts)); 1850 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1851 CeedCall(CeedQFunctionGetFlopsEstimate(qf, &qf_flops)); 1852 CeedCall(CeedQFunctionDestroy(&qf)); 1853 CeedCheck(qf_flops > -1, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE, 1854 "Must set CeedQFunction FLOPs estimate with CeedQFunctionSetUserFlopsEstimate"); 1855 *flops += num_elem * num_qpts * qf_flops; 1856 } 1857 1858 // Output FLOPs 1859 for (CeedInt i = 0; i < num_output_fields; i++) { 1860 CeedVector vec; 1861 1862 CeedCall(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 1863 if (vec == CEED_VECTOR_ACTIVE) { 1864 CeedEvalMode eval_mode; 1865 CeedSize rstr_flops, basis_flops; 1866 CeedElemRestriction rstr; 1867 CeedBasis basis; 1868 1869 CeedCall(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &rstr)); 1870 CeedCall(CeedElemRestrictionGetFlopsEstimate(rstr, CEED_TRANSPOSE, &rstr_flops)); 1871 CeedCall(CeedElemRestrictionDestroy(&rstr)); 1872 *flops += rstr_flops; 1873 CeedCall(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 1874 CeedCall(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 1875 CeedCall(CeedBasisGetFlopsEstimate(basis, CEED_TRANSPOSE, eval_mode, is_at_points, num_points, &basis_flops)); 1876 CeedCall(CeedBasisDestroy(&basis)); 1877 *flops += basis_flops * num_elem; 1878 } 1879 CeedCall(CeedVectorDestroy(&vec)); 1880 } 1881 } 1882 return CEED_ERROR_SUCCESS; 1883 } 1884 1885 /** 1886 @brief Get `CeedQFunction` global context for a `CeedOperator`. 1887 1888 The caller is responsible for destroying `ctx` returned from this function via @ref CeedQFunctionContextDestroy(). 1889 1890 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`. 1891 This `CeedQFunctionContext` will be destroyed if `ctx` is the only reference to this `CeedQFunctionContext`. 1892 1893 @param[in] op `CeedOperator` 1894 @param[out] ctx Variable to store `CeedQFunctionContext` 1895 1896 @return An error code: 0 - success, otherwise - failure 1897 1898 @ref Advanced 1899 **/ 1900 int CeedOperatorGetContext(CeedOperator op, CeedQFunctionContext *ctx) { 1901 bool is_composite; 1902 CeedQFunction qf; 1903 CeedQFunctionContext qf_ctx; 1904 1905 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1906 CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE, "Cannot retrieve CeedQFunctionContext for composite operator"); 1907 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1908 CeedCall(CeedQFunctionGetInnerContext(qf, &qf_ctx)); 1909 CeedCall(CeedQFunctionDestroy(&qf)); 1910 *ctx = NULL; 1911 if (qf_ctx) CeedCall(CeedQFunctionContextReferenceCopy(qf_ctx, ctx)); 1912 return CEED_ERROR_SUCCESS; 1913 } 1914 1915 /** 1916 @brief Get label for a registered `CeedQFunctionContext` field, or `NULL` if no field has been registered with this `field_name`. 1917 1918 Fields are registered via `CeedQFunctionContextRegister*()` functions (eg. @ref CeedQFunctionContextRegisterDouble()). 1919 1920 @param[in] op `CeedOperator` 1921 @param[in] field_name Name of field to retrieve label 1922 @param[out] field_label Variable to field label 1923 1924 @return An error code: 0 - success, otherwise - failure 1925 1926 @ref User 1927 **/ 1928 int CeedOperatorGetContextFieldLabel(CeedOperator op, const char *field_name, CeedContextFieldLabel *field_label) { 1929 bool is_composite, field_found = false; 1930 1931 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1932 1933 if (is_composite) { 1934 // Composite operator 1935 // -- Check if composite label already created 1936 for (CeedInt i = 0; i < op->num_context_labels; i++) { 1937 if (!strcmp(op->context_labels[i]->name, field_name)) { 1938 *field_label = op->context_labels[i]; 1939 return CEED_ERROR_SUCCESS; 1940 } 1941 } 1942 1943 // -- Create composite label if needed 1944 CeedInt num_sub; 1945 CeedOperator *sub_operators; 1946 CeedContextFieldLabel new_field_label; 1947 1948 CeedCall(CeedCalloc(1, &new_field_label)); 1949 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_sub)); 1950 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 1951 CeedCall(CeedCalloc(num_sub, &new_field_label->sub_labels)); 1952 new_field_label->num_sub_labels = num_sub; 1953 1954 for (CeedInt i = 0; i < num_sub; i++) { 1955 if (sub_operators[i]->qf->ctx) { 1956 CeedContextFieldLabel new_field_label_i; 1957 1958 CeedCall(CeedQFunctionContextGetFieldLabel(sub_operators[i]->qf->ctx, field_name, &new_field_label_i)); 1959 if (new_field_label_i) { 1960 field_found = true; 1961 new_field_label->sub_labels[i] = new_field_label_i; 1962 new_field_label->name = new_field_label_i->name; 1963 new_field_label->description = new_field_label_i->description; 1964 if (new_field_label->type && new_field_label->type != new_field_label_i->type) { 1965 // LCOV_EXCL_START 1966 CeedCall(CeedFree(&new_field_label)); 1967 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE, "Incompatible field types on sub-operator contexts. %s != %s", 1968 CeedContextFieldTypes[new_field_label->type], CeedContextFieldTypes[new_field_label_i->type]); 1969 // LCOV_EXCL_STOP 1970 } else { 1971 new_field_label->type = new_field_label_i->type; 1972 } 1973 if (new_field_label->num_values != 0 && new_field_label->num_values != new_field_label_i->num_values) { 1974 // LCOV_EXCL_START 1975 CeedCall(CeedFree(&new_field_label)); 1976 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE, 1977 "Incompatible field number of values on sub-operator contexts. %zu != %zu", new_field_label->num_values, 1978 new_field_label_i->num_values); 1979 // LCOV_EXCL_STOP 1980 } else { 1981 new_field_label->num_values = new_field_label_i->num_values; 1982 } 1983 } 1984 } 1985 } 1986 // -- Cleanup if field was found 1987 if (field_found) { 1988 *field_label = new_field_label; 1989 } else { 1990 // LCOV_EXCL_START 1991 CeedCall(CeedFree(&new_field_label->sub_labels)); 1992 CeedCall(CeedFree(&new_field_label)); 1993 *field_label = NULL; 1994 // LCOV_EXCL_STOP 1995 } 1996 } else { 1997 CeedQFunction qf; 1998 CeedQFunctionContext ctx; 1999 2000 // Single, non-composite operator 2001 CeedCall(CeedOperatorGetQFunction(op, &qf)); 2002 CeedCall(CeedQFunctionGetInnerContext(qf, &ctx)); 2003 CeedCall(CeedQFunctionDestroy(&qf)); 2004 if (ctx) { 2005 CeedCall(CeedQFunctionContextGetFieldLabel(ctx, field_name, field_label)); 2006 } else { 2007 *field_label = NULL; 2008 } 2009 } 2010 2011 // Set label in operator 2012 if (*field_label) { 2013 (*field_label)->from_op = true; 2014 2015 // Move new composite label to operator 2016 if (op->num_context_labels == 0) { 2017 CeedCall(CeedCalloc(1, &op->context_labels)); 2018 op->max_context_labels = 1; 2019 } else if (op->num_context_labels == op->max_context_labels) { 2020 CeedCall(CeedRealloc(2 * op->num_context_labels, &op->context_labels)); 2021 op->max_context_labels *= 2; 2022 } 2023 op->context_labels[op->num_context_labels] = *field_label; 2024 op->num_context_labels++; 2025 } 2026 return CEED_ERROR_SUCCESS; 2027 } 2028 2029 /** 2030 @brief Set `CeedQFunctionContext` field holding double precision values. 2031 2032 For composite operators, the values are set in all sub-operator `CeedQFunctionContext` that have a matching `field_name`. 2033 2034 @param[in,out] op `CeedOperator` 2035 @param[in] field_label Label of field to set 2036 @param[in] values Values to set 2037 2038 @return An error code: 0 - success, otherwise - failure 2039 2040 @ref User 2041 **/ 2042 int CeedOperatorSetContextDouble(CeedOperator op, CeedContextFieldLabel field_label, double *values) { 2043 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, values); 2044 } 2045 2046 /** 2047 @brief Get `CeedQFunctionContext` field holding double precision values, read-only. 2048 2049 For composite operators, the values correspond to the first sub-operator `CeedQFunctionContext` that has a matching `field_name`. 2050 2051 @param[in] op `CeedOperator` 2052 @param[in] field_label Label of field to get 2053 @param[out] num_values Number of values in the field label 2054 @param[out] values Pointer to context values 2055 2056 @return An error code: 0 - success, otherwise - failure 2057 2058 @ref User 2059 **/ 2060 int CeedOperatorGetContextDoubleRead(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const double **values) { 2061 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, num_values, values); 2062 } 2063 2064 /** 2065 @brief Restore `CeedQFunctionContext` field holding double precision values, read-only. 2066 2067 @param[in] op `CeedOperator` 2068 @param[in] field_label Label of field to restore 2069 @param[out] values Pointer to context values 2070 2071 @return An error code: 0 - success, otherwise - failure 2072 2073 @ref User 2074 **/ 2075 int CeedOperatorRestoreContextDoubleRead(CeedOperator op, CeedContextFieldLabel field_label, const double **values) { 2076 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, values); 2077 } 2078 2079 /** 2080 @brief Set `CeedQFunctionContext` field holding `int32` values. 2081 2082 For composite operators, the values are set in all sub-operator `CeedQFunctionContext` that have a matching `field_name`. 2083 2084 @param[in,out] op `CeedOperator` 2085 @param[in] field_label Label of field to set 2086 @param[in] values Values to set 2087 2088 @return An error code: 0 - success, otherwise - failure 2089 2090 @ref User 2091 **/ 2092 int CeedOperatorSetContextInt32(CeedOperator op, CeedContextFieldLabel field_label, int32_t *values) { 2093 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_INT32, values); 2094 } 2095 2096 /** 2097 @brief Get `CeedQFunctionContext` field holding `int32` values, read-only. 2098 2099 For composite operators, the values correspond to the first sub-operator `CeedQFunctionContext` that has a matching `field_name`. 2100 2101 @param[in] op `CeedOperator` 2102 @param[in] field_label Label of field to get 2103 @param[out] num_values Number of `int32` values in `values` 2104 @param[out] values Pointer to context values 2105 2106 @return An error code: 0 - success, otherwise - failure 2107 2108 @ref User 2109 **/ 2110 int CeedOperatorGetContextInt32Read(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const int32_t **values) { 2111 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_INT32, num_values, values); 2112 } 2113 2114 /** 2115 @brief Restore `CeedQFunctionContext` field holding `int32` values, read-only. 2116 2117 @param[in] op `CeedOperator` 2118 @param[in] field_label Label of field to get 2119 @param[out] values Pointer to context values 2120 2121 @return An error code: 0 - success, otherwise - failure 2122 2123 @ref User 2124 **/ 2125 int CeedOperatorRestoreContextInt32Read(CeedOperator op, CeedContextFieldLabel field_label, const int32_t **values) { 2126 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_INT32, values); 2127 } 2128 2129 /** 2130 @brief Set `CeedQFunctionContext` field holding boolean values. 2131 2132 For composite operators, the values are set in all sub-operator `CeedQFunctionContext` that have a matching `field_name`. 2133 2134 @param[in,out] op `CeedOperator` 2135 @param[in] field_label Label of field to set 2136 @param[in] values Values to set 2137 2138 @return An error code: 0 - success, otherwise - failure 2139 2140 @ref User 2141 **/ 2142 int CeedOperatorSetContextBoolean(CeedOperator op, CeedContextFieldLabel field_label, bool *values) { 2143 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_BOOL, values); 2144 } 2145 2146 /** 2147 @brief Get `CeedQFunctionContext` field holding boolean values, read-only. 2148 2149 For composite operators, the values correspond to the first sub-operator `CeedQFunctionContext` that has a matching `field_name`. 2150 2151 @param[in] op `CeedOperator` 2152 @param[in] field_label Label of field to get 2153 @param[out] num_values Number of boolean values in `values` 2154 @param[out] values Pointer to context values 2155 2156 @return An error code: 0 - success, otherwise - failure 2157 2158 @ref User 2159 **/ 2160 int CeedOperatorGetContextBooleanRead(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const bool **values) { 2161 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_BOOL, num_values, values); 2162 } 2163 2164 /** 2165 @brief Restore `CeedQFunctionContext` field holding boolean values, read-only. 2166 2167 @param[in] op `CeedOperator` 2168 @param[in] field_label Label of field to get 2169 @param[out] values Pointer to context values 2170 2171 @return An error code: 0 - success, otherwise - failure 2172 2173 @ref User 2174 **/ 2175 int CeedOperatorRestoreContextBooleanRead(CeedOperator op, CeedContextFieldLabel field_label, const bool **values) { 2176 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_BOOL, values); 2177 } 2178 2179 /** 2180 @brief Apply `CeedOperator` to a `CeedVector`. 2181 2182 This computes the action of the operator on the specified (active) input, yielding its (active) output. 2183 All inputs and outputs must be specified using @ref CeedOperatorSetField(). 2184 2185 @note Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2186 2187 @param[in] op `CeedOperator` to apply 2188 @param[in] in `CeedVector` containing input state or @ref CEED_VECTOR_NONE if there are no active inputs 2189 @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 2190 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2191 2192 @return An error code: 0 - success, otherwise - failure 2193 2194 @ref User 2195 **/ 2196 int CeedOperatorApply(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) { 2197 bool is_composite; 2198 2199 CeedCall(CeedOperatorCheckReady(op)); 2200 2201 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2202 if (is_composite && op->ApplyComposite) { 2203 // Composite Operator 2204 CeedCall(op->ApplyComposite(op, in, out, request)); 2205 } else if (!is_composite && op->Apply) { 2206 // Standard Operator 2207 CeedCall(op->Apply(op, in, out, request)); 2208 } else { 2209 // Standard or composite, default to zeroing out and calling ApplyAddActive 2210 // Zero active output 2211 if (out != CEED_VECTOR_NONE) CeedCall(CeedVectorSetValue(out, 0.0)); 2212 2213 // ApplyAddActive 2214 CeedCall(CeedOperatorApplyAddActive(op, in, out, request)); 2215 } 2216 return CEED_ERROR_SUCCESS; 2217 } 2218 2219 /** 2220 @brief Apply `CeedOperator` to a `CeedVector` and add result to output `CeedVector`. 2221 2222 This computes the action of the operator on the specified (active) input, yielding its (active) output. 2223 All inputs and outputs must be specified using @ref CeedOperatorSetField(). 2224 2225 @note Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2226 @warning This function adds into ALL outputs, including passive outputs. To only add into the active output, use `CeedOperatorApplyAddActive()`. 2227 @see `CeedOperatorApplyAddActive()` 2228 2229 @param[in] op `CeedOperator` to apply 2230 @param[in] in `CeedVector` containing input state or @ref CEED_VECTOR_NONE if there are no active inputs 2231 @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 2232 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2233 2234 @return An error code: 0 - success, otherwise - failure 2235 2236 @ref User 2237 **/ 2238 int CeedOperatorApplyAdd(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) { 2239 bool is_composite; 2240 2241 CeedCall(CeedOperatorCheckReady(op)); 2242 2243 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2244 if (is_composite) { 2245 // Composite Operator 2246 if (op->ApplyAddComposite) { 2247 CeedCall(op->ApplyAddComposite(op, in, out, request)); 2248 } else { 2249 CeedInt num_suboperators; 2250 CeedOperator *sub_operators; 2251 2252 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 2253 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 2254 for (CeedInt i = 0; i < num_suboperators; i++) { 2255 CeedCall(CeedOperatorApplyAdd(sub_operators[i], in, out, request)); 2256 } 2257 } 2258 } else if (op->num_elem > 0) { 2259 // Standard Operator 2260 CeedCall(op->ApplyAdd(op, in, out, request)); 2261 } 2262 return CEED_ERROR_SUCCESS; 2263 } 2264 2265 /** 2266 @brief Apply `CeedOperator` to a `CeedVector` and add result to output `CeedVector`. Only sums into active outputs, overwrites passive outputs. 2267 2268 This computes the action of the operator on the specified (active) input, yielding its (active) output. 2269 All inputs and outputs must be specified using @ref CeedOperatorSetField(). 2270 2271 @note Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2272 2273 @param[in] op `CeedOperator` to apply 2274 @param[in] in `CeedVector` containing input state or @ref CEED_VECTOR_NONE if there are no active inputs 2275 @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 2276 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2277 2278 @return An error code: 0 - success, otherwise - failure 2279 2280 @ref User 2281 **/ 2282 int CeedOperatorApplyAddActive(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) { 2283 bool is_composite; 2284 2285 CeedCall(CeedOperatorCheckReady(op)); 2286 2287 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2288 if (is_composite) { 2289 // Composite Operator 2290 CeedInt num_suboperators; 2291 CeedOperator *sub_operators; 2292 2293 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 2294 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 2295 2296 // Zero all output vectors 2297 for (CeedInt i = 0; i < num_suboperators; i++) { 2298 CeedInt num_output_fields; 2299 CeedOperatorField *output_fields; 2300 2301 CeedCall(CeedOperatorGetFields(sub_operators[i], NULL, NULL, &num_output_fields, &output_fields)); 2302 for (CeedInt j = 0; j < num_output_fields; j++) { 2303 CeedVector vec; 2304 2305 CeedCall(CeedOperatorFieldGetVector(output_fields[j], &vec)); 2306 if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) CeedCall(CeedVectorSetValue(vec, 0.0)); 2307 CeedCall(CeedVectorDestroy(&vec)); 2308 } 2309 } 2310 // ApplyAdd 2311 CeedCall(CeedOperatorApplyAdd(op, in, out, request)); 2312 } else { 2313 // Standard Operator 2314 CeedInt num_output_fields; 2315 CeedOperatorField *output_fields; 2316 2317 CeedCall(CeedOperatorGetFields(op, NULL, NULL, &num_output_fields, &output_fields)); 2318 // Zero all output vectors 2319 for (CeedInt i = 0; i < num_output_fields; i++) { 2320 CeedVector vec; 2321 2322 CeedCall(CeedOperatorFieldGetVector(output_fields[i], &vec)); 2323 if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) CeedCall(CeedVectorSetValue(vec, 0.0)); 2324 CeedCall(CeedVectorDestroy(&vec)); 2325 } 2326 // ApplyAdd 2327 CeedCall(CeedOperatorApplyAdd(op, in, out, request)); 2328 } 2329 return CEED_ERROR_SUCCESS; 2330 } 2331 2332 /** 2333 @brief Destroy temporary assembly data associated with a `CeedOperator` 2334 2335 @param[in,out] op `CeedOperator` whose assembly data to destroy 2336 2337 @return An error code: 0 - success, otherwise - failure 2338 2339 @ref User 2340 **/ 2341 int CeedOperatorAssemblyDataStrip(CeedOperator op) { 2342 bool is_composite; 2343 2344 CeedCall(CeedQFunctionAssemblyDataDestroy(&op->qf_assembled)); 2345 CeedCall(CeedOperatorAssemblyDataDestroy(&op->op_assembled)); 2346 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2347 if (is_composite) { 2348 CeedInt num_suboperators; 2349 CeedOperator *sub_operators; 2350 2351 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 2352 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 2353 for (CeedInt i = 0; i < num_suboperators; i++) { 2354 CeedCall(CeedQFunctionAssemblyDataDestroy(&sub_operators[i]->qf_assembled)); 2355 CeedCall(CeedOperatorAssemblyDataDestroy(&sub_operators[i]->op_assembled)); 2356 } 2357 } 2358 return CEED_ERROR_SUCCESS; 2359 } 2360 2361 /** 2362 @brief Destroy a `CeedOperator` 2363 2364 @param[in,out] op `CeedOperator` to destroy 2365 2366 @return An error code: 0 - success, otherwise - failure 2367 2368 @ref User 2369 **/ 2370 int CeedOperatorDestroy(CeedOperator *op) { 2371 if (!*op || --(*op)->ref_count > 0) { 2372 *op = NULL; 2373 return CEED_ERROR_SUCCESS; 2374 } 2375 // Backend destroy 2376 if ((*op)->Destroy) { 2377 CeedCall((*op)->Destroy(*op)); 2378 } 2379 // Free fields 2380 for (CeedInt i = 0; i < (*op)->num_fields; i++) { 2381 if ((*op)->input_fields[i]) { 2382 if ((*op)->input_fields[i]->elem_rstr != CEED_ELEMRESTRICTION_NONE) { 2383 CeedCall(CeedElemRestrictionDestroy(&(*op)->input_fields[i]->elem_rstr)); 2384 } 2385 if ((*op)->input_fields[i]->basis != CEED_BASIS_NONE) { 2386 CeedCall(CeedBasisDestroy(&(*op)->input_fields[i]->basis)); 2387 } 2388 if ((*op)->input_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->input_fields[i]->vec != CEED_VECTOR_NONE) { 2389 CeedCall(CeedVectorDestroy(&(*op)->input_fields[i]->vec)); 2390 } 2391 CeedCall(CeedFree(&(*op)->input_fields[i]->field_name)); 2392 CeedCall(CeedFree(&(*op)->input_fields[i])); 2393 } 2394 } 2395 for (CeedInt i = 0; i < (*op)->num_fields; i++) { 2396 if ((*op)->output_fields[i]) { 2397 CeedCall(CeedElemRestrictionDestroy(&(*op)->output_fields[i]->elem_rstr)); 2398 if ((*op)->output_fields[i]->basis != CEED_BASIS_NONE) { 2399 CeedCall(CeedBasisDestroy(&(*op)->output_fields[i]->basis)); 2400 } 2401 if ((*op)->output_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->output_fields[i]->vec != CEED_VECTOR_NONE) { 2402 CeedCall(CeedVectorDestroy(&(*op)->output_fields[i]->vec)); 2403 } 2404 CeedCall(CeedFree(&(*op)->output_fields[i]->field_name)); 2405 CeedCall(CeedFree(&(*op)->output_fields[i])); 2406 } 2407 } 2408 CeedCall(CeedFree(&(*op)->input_fields)); 2409 CeedCall(CeedFree(&(*op)->output_fields)); 2410 // Destroy AtPoints data 2411 CeedCall(CeedVectorDestroy(&(*op)->point_coords)); 2412 CeedCall(CeedElemRestrictionDestroy(&(*op)->rstr_points)); 2413 CeedCall(CeedElemRestrictionDestroy(&(*op)->first_points_rstr)); 2414 // Destroy assembly data (must happen before destroying sub_operators) 2415 CeedCall(CeedOperatorAssemblyDataStrip(*op)); 2416 // Destroy sub_operators 2417 for (CeedInt i = 0; i < (*op)->num_suboperators; i++) { 2418 if ((*op)->sub_operators[i]) { 2419 CeedCall(CeedOperatorDestroy(&(*op)->sub_operators[i])); 2420 } 2421 } 2422 CeedCall(CeedFree(&(*op)->sub_operators)); 2423 CeedCall(CeedQFunctionDestroy(&(*op)->qf)); 2424 CeedCall(CeedQFunctionDestroy(&(*op)->dqf)); 2425 CeedCall(CeedQFunctionDestroy(&(*op)->dqfT)); 2426 // Destroy any composite labels 2427 if ((*op)->is_composite) { 2428 for (CeedInt i = 0; i < (*op)->num_context_labels; i++) { 2429 CeedCall(CeedFree(&(*op)->context_labels[i]->sub_labels)); 2430 CeedCall(CeedFree(&(*op)->context_labels[i])); 2431 } 2432 } 2433 CeedCall(CeedFree(&(*op)->context_labels)); 2434 2435 // Destroy fallback 2436 CeedCall(CeedOperatorDestroy(&(*op)->op_fallback)); 2437 2438 CeedCall(CeedFree(&(*op)->name)); 2439 CeedCall(CeedDestroy(&(*op)->ceed)); 2440 CeedCall(CeedFree(op)); 2441 return CEED_ERROR_SUCCESS; 2442 } 2443 2444 /// @} 2445