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