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 Set whether the sub-operators of the composite `CeedOperator` must be run sequentially. 1395 1396 Note: This value currently only affects the GPU `/gpu/cuda/gen` and `/gpu/hip/gen` backends. 1397 1398 @param[in] op Composite `CeedOperator` 1399 @param[in] is_sequential Flag value to set, if `true`, forces the composite `CeedOperator` to execute sequentially 1400 1401 @return An error code: 0 - success, otherwise - failure 1402 1403 @ref Advanced 1404 **/ 1405 int CeedOperatorCompositeSetSequential(CeedOperator op, bool is_sequential) { 1406 bool is_composite; 1407 1408 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1409 CeedCheck(is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Only defined for a composite operator"); 1410 op->is_sequential = is_sequential; 1411 return CEED_ERROR_SUCCESS; 1412 } 1413 1414 /** 1415 @brief Get whether the sub-operators of the composite `CeedOperator` must be run sequentially. 1416 1417 Note: This value currently only affects the GPU `/gpu/cuda/gen` and `/gpu/hip/gen` backends. 1418 1419 @param[in] op Composite `CeedOperator` 1420 @param[out] is_sequential Variable to store sequential status 1421 1422 @return An error code: 0 - success, otherwise - failure 1423 1424 @ref Advanced 1425 **/ 1426 int CeedOperatorCompositeIsSequential(CeedOperator op, bool *is_sequential) { 1427 bool is_composite; 1428 1429 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1430 CeedCheck(is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Only defined for a composite operator"); 1431 *is_sequential = op->is_sequential; 1432 return CEED_ERROR_SUCCESS; 1433 } 1434 1435 /** 1436 @brief Check if a `CeedOperator` is ready to be used. 1437 1438 @param[in] op `CeedOperator` to check 1439 1440 @return An error code: 0 - success, otherwise - failure 1441 1442 @ref User 1443 **/ 1444 int CeedOperatorCheckReady(CeedOperator op) { 1445 bool is_at_points, is_composite; 1446 CeedQFunction qf = NULL; 1447 1448 if (op->is_interface_setup) return CEED_ERROR_SUCCESS; 1449 1450 CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); 1451 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1452 if (!is_composite) CeedCall(CeedOperatorGetQFunction(op, &qf)); 1453 if (is_composite) { 1454 CeedInt num_suboperators; 1455 1456 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 1457 if (!num_suboperators) { 1458 // Empty operator setup 1459 op->input_size = 0; 1460 op->output_size = 0; 1461 } else { 1462 CeedOperator *sub_operators; 1463 1464 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 1465 for (CeedInt i = 0; i < num_suboperators; i++) { 1466 CeedCall(CeedOperatorCheckReady(sub_operators[i])); 1467 } 1468 // Sub-operators could be modified after adding to composite operator 1469 // Need to verify no lvec incompatibility from any changes 1470 CeedSize input_size, output_size; 1471 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 1472 } 1473 } else { 1474 CeedInt num_input_fields, num_output_fields; 1475 1476 CeedCheck(op->num_fields > 0, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE, "No operator fields set"); 1477 CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, NULL, &num_output_fields, NULL)); 1478 CeedCheck(op->num_fields == num_input_fields + num_output_fields, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE, 1479 "Not all operator fields set"); 1480 CeedCheck(op->has_restriction, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE, "At least one restriction required"); 1481 CeedCheck(op->num_qpts > 0 || is_at_points, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE, 1482 "At least one non-collocated CeedBasis is required or the number of quadrature points must be set"); 1483 } 1484 1485 // Flag as immutable and ready 1486 op->is_interface_setup = true; 1487 if (qf && qf != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionSetImmutable(qf)); 1488 CeedCall(CeedQFunctionDestroy(&qf)); 1489 if (op->dqf && op->dqf != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionSetImmutable(op->dqf)); 1490 if (op->dqfT && op->dqfT != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionSetImmutable(op->dqfT)); 1491 return CEED_ERROR_SUCCESS; 1492 } 1493 1494 /** 1495 @brief Get vector lengths for the active input and/or output `CeedVector` of a `CeedOperator`. 1496 1497 Note: Lengths of `-1` indicate that the CeedOperator does not have an active input and/or output. 1498 1499 @param[in] op `CeedOperator` 1500 @param[out] input_size Variable to store active input vector length, or `NULL` 1501 @param[out] output_size Variable to store active output vector length, or `NULL` 1502 1503 @return An error code: 0 - success, otherwise - failure 1504 1505 @ref User 1506 **/ 1507 int CeedOperatorGetActiveVectorLengths(CeedOperator op, CeedSize *input_size, CeedSize *output_size) { 1508 bool is_composite; 1509 1510 if (input_size) *input_size = op->input_size; 1511 if (output_size) *output_size = op->output_size; 1512 1513 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1514 if (is_composite && (op->input_size == -1 || op->output_size == -1)) { 1515 CeedInt num_suboperators; 1516 CeedOperator *sub_operators; 1517 1518 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 1519 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 1520 for (CeedInt i = 0; i < num_suboperators; i++) { 1521 CeedSize sub_input_size, sub_output_size; 1522 1523 CeedCall(CeedOperatorGetActiveVectorLengths(sub_operators[i], &sub_input_size, &sub_output_size)); 1524 if (op->input_size == -1) op->input_size = sub_input_size; 1525 if (op->output_size == -1) op->output_size = sub_output_size; 1526 // Note, a size of -1 means no active vector restriction set, so no incompatibility 1527 CeedCheck((sub_input_size == -1 || sub_input_size == op->input_size) && (sub_output_size == -1 || sub_output_size == op->output_size), 1528 CeedOperatorReturnCeed(op), CEED_ERROR_MAJOR, 1529 "Sub-operators must have compatible dimensions; composite operator of shape (%" CeedSize_FMT ", %" CeedSize_FMT 1530 ") not compatible with sub-operator of " 1531 "shape (%" CeedSize_FMT ", %" CeedSize_FMT ")", 1532 op->input_size, op->output_size, input_size, output_size); 1533 } 1534 } 1535 return CEED_ERROR_SUCCESS; 1536 } 1537 1538 /** 1539 @brief Set reuse of `CeedQFunction` data in `CeedOperatorLinearAssemble*()` functions. 1540 1541 When `reuse_assembly_data = false` (default), the `CeedQFunction` associated with this `CeedOperator` is re-assembled every time a `CeedOperatorLinearAssemble*()` function is called. 1542 When `reuse_assembly_data = true`, the `CeedQFunction` associated with this `CeedOperator` is reused between calls to @ref CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(). 1543 1544 @param[in] op `CeedOperator` 1545 @param[in] reuse_assembly_data Boolean flag setting assembly data reuse 1546 1547 @return An error code: 0 - success, otherwise - failure 1548 1549 @ref Advanced 1550 **/ 1551 int CeedOperatorSetQFunctionAssemblyReuse(CeedOperator op, bool reuse_assembly_data) { 1552 bool is_composite; 1553 1554 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1555 if (is_composite) { 1556 for (CeedInt i = 0; i < op->num_suboperators; i++) { 1557 CeedCall(CeedOperatorSetQFunctionAssemblyReuse(op->sub_operators[i], reuse_assembly_data)); 1558 } 1559 } else { 1560 CeedQFunctionAssemblyData data; 1561 1562 CeedCall(CeedOperatorGetQFunctionAssemblyData(op, &data)); 1563 CeedCall(CeedQFunctionAssemblyDataSetReuse(data, reuse_assembly_data)); 1564 } 1565 return CEED_ERROR_SUCCESS; 1566 } 1567 1568 /** 1569 @brief Mark `CeedQFunction` data as updated and the `CeedQFunction` as requiring re-assembly. 1570 1571 @param[in] op `CeedOperator` 1572 @param[in] needs_data_update Boolean flag setting assembly data reuse 1573 1574 @return An error code: 0 - success, otherwise - failure 1575 1576 @ref Advanced 1577 **/ 1578 int CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(CeedOperator op, bool needs_data_update) { 1579 bool is_composite; 1580 1581 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1582 if (is_composite) { 1583 CeedInt num_suboperators; 1584 CeedOperator *sub_operators; 1585 1586 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 1587 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 1588 for (CeedInt i = 0; i < num_suboperators; i++) { 1589 CeedCall(CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(sub_operators[i], needs_data_update)); 1590 } 1591 } else { 1592 CeedQFunctionAssemblyData data; 1593 1594 CeedCall(CeedOperatorGetQFunctionAssemblyData(op, &data)); 1595 CeedCall(CeedQFunctionAssemblyDataSetUpdateNeeded(data, needs_data_update)); 1596 } 1597 return CEED_ERROR_SUCCESS; 1598 } 1599 1600 /** 1601 @brief Set name of `CeedOperator` for @ref CeedOperatorView() output 1602 1603 @param[in,out] op `CeedOperator` 1604 @param[in] name Name to set, or NULL to remove previously set name 1605 1606 @return An error code: 0 - success, otherwise - failure 1607 1608 @ref User 1609 **/ 1610 int CeedOperatorSetName(CeedOperator op, const char *name) { 1611 char *name_copy; 1612 size_t name_len = name ? strlen(name) : 0; 1613 1614 CeedCall(CeedFree(&op->name)); 1615 if (name_len > 0) { 1616 CeedCall(CeedCalloc(name_len + 1, &name_copy)); 1617 memcpy(name_copy, name, name_len); 1618 op->name = name_copy; 1619 } 1620 return CEED_ERROR_SUCCESS; 1621 } 1622 1623 /** 1624 @brief Get name of `CeedOperator` 1625 1626 @param[in] op `CeedOperator` 1627 @param[in,out] name Address of variable to hold currently set name 1628 1629 @return An error code: 0 - success, otherwise - failure 1630 1631 @ref User 1632 **/ 1633 int CeedOperatorGetName(CeedOperator op, const char **name) { 1634 if (op->name) { 1635 *name = op->name; 1636 } else if (!op->is_composite) { 1637 CeedQFunction qf; 1638 1639 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1640 if (qf) CeedCall(CeedQFunctionGetName(qf, name)); 1641 CeedCall(CeedQFunctionDestroy(&qf)); 1642 } 1643 return CEED_ERROR_SUCCESS; 1644 } 1645 1646 /** 1647 @brief Core logic for viewing a `CeedOperator` 1648 1649 @param[in] op `CeedOperator` to view brief summary 1650 @param[in] stream Stream to write; typically `stdout` or a file 1651 @param[in] is_full Whether to write full operator view or terse 1652 1653 @return Error code: 0 - success, otherwise - failure 1654 1655 @ref Developer 1656 **/ 1657 static int CeedOperatorView_Core(CeedOperator op, FILE *stream, bool is_full) { 1658 bool has_name, is_composite, is_at_points; 1659 char *tabs = NULL; 1660 const char *name = NULL; 1661 CeedInt num_tabs = 0; 1662 1663 CeedCall(CeedOperatorGetName(op, &name)); 1664 has_name = name ? strlen(name) : false; 1665 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1666 CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); 1667 // Set tabs 1668 CeedCall(CeedOperatorGetNumViewTabs(op, &num_tabs)); 1669 CeedCall(CeedCalloc(CEED_TAB_WIDTH * (num_tabs + is_composite) + 1, &tabs)); 1670 for (CeedInt i = 0; i < CEED_TAB_WIDTH * num_tabs; i++) tabs[i] = ' '; 1671 if (is_composite) { 1672 CeedInt num_suboperators; 1673 CeedOperator *sub_operators; 1674 1675 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 1676 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 1677 fprintf(stream, "%s", tabs); 1678 fprintf(stream, "Composite CeedOperator%s%s\n", has_name ? " - " : "", has_name ? name : ""); 1679 for (CeedInt i = 0; i < CEED_TAB_WIDTH; i++) tabs[CEED_TAB_WIDTH * num_tabs + i] = ' '; 1680 for (CeedInt i = 0; i < num_suboperators; i++) { 1681 has_name = sub_operators[i]->name; 1682 fprintf(stream, "%s", tabs); 1683 fprintf(stream, "SubOperator%s %" CeedInt_FMT "%s%s%s\n", is_at_points ? " AtPoints" : "", i, has_name ? " - " : "", 1684 has_name ? sub_operators[i]->name : "", is_full ? ":" : ""); 1685 if (is_full) CeedCall(CeedOperatorSingleView(sub_operators[i], tabs, stream)); 1686 } 1687 } else { 1688 fprintf(stream, "%s", tabs); 1689 fprintf(stream, "CeedOperator%s%s%s\n", is_at_points ? " AtPoints" : "", has_name ? " - " : "", has_name ? name : ""); 1690 if (is_full) CeedCall(CeedOperatorSingleView(op, tabs, stream)); 1691 } 1692 CeedCall(CeedFree(&tabs)); 1693 return CEED_ERROR_SUCCESS; 1694 } 1695 1696 /** 1697 @brief Set the number of tabs to indent for @ref CeedOperatorView() output 1698 1699 @param[in] op `CeedOperator` to set the number of view tabs 1700 @param[in] num_tabs Number of view tabs to set 1701 1702 @return Error code: 0 - success, otherwise - failure 1703 1704 @ref User 1705 **/ 1706 int CeedOperatorSetNumViewTabs(CeedOperator op, CeedInt num_tabs) { 1707 CeedCall(CeedObjectSetNumViewTabs((CeedObject)op, num_tabs)); 1708 return CEED_ERROR_SUCCESS; 1709 } 1710 1711 /** 1712 @brief Get the number of tabs to indent for @ref CeedOperatorView() output 1713 1714 @param[in] op `CeedOperator` to get the number of view tabs 1715 @param[out] num_tabs Number of view tabs 1716 1717 @return Error code: 0 - success, otherwise - failure 1718 1719 @ref User 1720 **/ 1721 int CeedOperatorGetNumViewTabs(CeedOperator op, CeedInt *num_tabs) { 1722 CeedCall(CeedObjectGetNumViewTabs((CeedObject)op, num_tabs)); 1723 return CEED_ERROR_SUCCESS; 1724 } 1725 1726 /** 1727 @brief View a `CeedOperator` 1728 1729 @param[in] op `CeedOperator` to view 1730 @param[in] stream Stream to write; typically `stdout` or a file 1731 1732 @return Error code: 0 - success, otherwise - failure 1733 1734 @ref User 1735 **/ 1736 int CeedOperatorView(CeedOperator op, FILE *stream) { 1737 CeedCall(CeedOperatorView_Core(op, stream, true)); 1738 return CEED_ERROR_SUCCESS; 1739 } 1740 1741 /** 1742 @brief View a brief summary `CeedOperator` 1743 1744 @param[in] op `CeedOperator` to view brief summary 1745 @param[in] stream Stream to write; typically `stdout` or a file 1746 1747 @return Error code: 0 - success, otherwise - failure 1748 1749 @ref User 1750 **/ 1751 int CeedOperatorViewTerse(CeedOperator op, FILE *stream) { 1752 CeedCall(CeedOperatorView_Core(op, stream, false)); 1753 return CEED_ERROR_SUCCESS; 1754 } 1755 1756 /** 1757 @brief Get the `Ceed` associated with a `CeedOperator` 1758 1759 @param[in] op `CeedOperator` 1760 @param[out] ceed Variable to store `Ceed` 1761 1762 @return An error code: 0 - success, otherwise - failure 1763 1764 @ref Advanced 1765 **/ 1766 int CeedOperatorGetCeed(CeedOperator op, Ceed *ceed) { 1767 CeedCall(CeedObjectGetCeed((CeedObject)op, ceed)); 1768 return CEED_ERROR_SUCCESS; 1769 } 1770 1771 /** 1772 @brief Return the `Ceed` associated with a `CeedOperator` 1773 1774 @param[in] op `CeedOperator` 1775 1776 @return `Ceed` associated with the `op` 1777 1778 @ref Advanced 1779 **/ 1780 Ceed CeedOperatorReturnCeed(CeedOperator op) { return CeedObjectReturnCeed((CeedObject)op); } 1781 1782 /** 1783 @brief Get the number of elements associated with a `CeedOperator` 1784 1785 @param[in] op `CeedOperator` 1786 @param[out] num_elem Variable to store number of elements 1787 1788 @return An error code: 0 - success, otherwise - failure 1789 1790 @ref Advanced 1791 **/ 1792 int CeedOperatorGetNumElements(CeedOperator op, CeedInt *num_elem) { 1793 bool is_composite; 1794 1795 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1796 CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Not defined for composite operator"); 1797 *num_elem = op->num_elem; 1798 return CEED_ERROR_SUCCESS; 1799 } 1800 1801 /** 1802 @brief Get the number of quadrature points associated with a `CeedOperator` 1803 1804 @param[in] op `CeedOperator` 1805 @param[out] num_qpts Variable to store vector number of quadrature points 1806 1807 @return An error code: 0 - success, otherwise - failure 1808 1809 @ref Advanced 1810 **/ 1811 int CeedOperatorGetNumQuadraturePoints(CeedOperator op, CeedInt *num_qpts) { 1812 bool is_composite; 1813 1814 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1815 CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Not defined for composite operator"); 1816 *num_qpts = op->num_qpts; 1817 return CEED_ERROR_SUCCESS; 1818 } 1819 1820 /** 1821 @brief Estimate number of FLOPs required to apply `CeedOperator` on the active `CeedVector` 1822 1823 @param[in] op `CeedOperator` to estimate FLOPs for 1824 @param[out] flops Address of variable to hold FLOPs estimate 1825 1826 @ref Backend 1827 **/ 1828 int CeedOperatorGetFlopsEstimate(CeedOperator op, CeedSize *flops) { 1829 bool is_composite; 1830 1831 CeedCall(CeedOperatorCheckReady(op)); 1832 1833 *flops = 0; 1834 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1835 if (is_composite) { 1836 CeedInt num_suboperators; 1837 1838 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 1839 CeedOperator *sub_operators; 1840 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 1841 1842 // FLOPs for each suboperator 1843 for (CeedInt i = 0; i < num_suboperators; i++) { 1844 CeedSize suboperator_flops; 1845 1846 CeedCall(CeedOperatorGetFlopsEstimate(sub_operators[i], &suboperator_flops)); 1847 *flops += suboperator_flops; 1848 } 1849 } else { 1850 bool is_at_points; 1851 CeedInt num_input_fields, num_output_fields, num_elem = 0, num_points = 0; 1852 CeedQFunction qf; 1853 CeedQFunctionField *qf_input_fields, *qf_output_fields; 1854 CeedOperatorField *op_input_fields, *op_output_fields; 1855 1856 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 1857 if (num_elem == 0) return CEED_ERROR_SUCCESS; 1858 CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); 1859 if (is_at_points) { 1860 CeedMemType mem_type; 1861 CeedElemRestriction rstr_points = NULL; 1862 1863 CeedCall(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL)); 1864 CeedCall(CeedGetPreferredMemType(CeedOperatorReturnCeed(op), &mem_type)); 1865 if (mem_type == CEED_MEM_DEVICE) { 1866 // Device backends pad out to the same number of points per element 1867 CeedCall(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &num_points)); 1868 } else { 1869 num_points = 0; 1870 for (CeedInt i = 0; i < num_elem; i++) { 1871 CeedInt points_in_elem = 0; 1872 1873 CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr_points, i, &points_in_elem)); 1874 num_points += points_in_elem; 1875 } 1876 num_points = num_points / num_elem + (num_points % num_elem > 0); 1877 } 1878 CeedCall(CeedElemRestrictionDestroy(&rstr_points)); 1879 } 1880 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1881 CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &qf_input_fields, &num_output_fields, &qf_output_fields)); 1882 CeedCall(CeedQFunctionDestroy(&qf)); 1883 CeedCall(CeedOperatorGetFields(op, NULL, &op_input_fields, NULL, &op_output_fields)); 1884 1885 // Input FLOPs 1886 for (CeedInt i = 0; i < num_input_fields; i++) { 1887 CeedVector vec; 1888 1889 CeedCall(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1890 if (vec == CEED_VECTOR_ACTIVE) { 1891 CeedEvalMode eval_mode; 1892 CeedSize rstr_flops, basis_flops; 1893 CeedElemRestriction rstr; 1894 CeedBasis basis; 1895 1896 CeedCall(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &rstr)); 1897 CeedCall(CeedElemRestrictionGetFlopsEstimate(rstr, CEED_NOTRANSPOSE, &rstr_flops)); 1898 CeedCall(CeedElemRestrictionDestroy(&rstr)); 1899 *flops += rstr_flops; 1900 CeedCall(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 1901 CeedCall(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 1902 CeedCall(CeedBasisGetFlopsEstimate(basis, CEED_NOTRANSPOSE, eval_mode, is_at_points, num_points, &basis_flops)); 1903 CeedCall(CeedBasisDestroy(&basis)); 1904 *flops += basis_flops * num_elem; 1905 } 1906 CeedCall(CeedVectorDestroy(&vec)); 1907 } 1908 // QF FLOPs 1909 { 1910 CeedInt num_qpts; 1911 CeedSize qf_flops; 1912 CeedQFunction qf; 1913 1914 if (is_at_points) num_qpts = num_points; 1915 else CeedCall(CeedOperatorGetNumQuadraturePoints(op, &num_qpts)); 1916 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1917 CeedCall(CeedQFunctionGetFlopsEstimate(qf, &qf_flops)); 1918 CeedCall(CeedQFunctionDestroy(&qf)); 1919 CeedCheck(qf_flops > -1, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE, 1920 "Must set CeedQFunction FLOPs estimate with CeedQFunctionSetUserFlopsEstimate"); 1921 *flops += num_elem * num_qpts * qf_flops; 1922 } 1923 1924 // Output FLOPs 1925 for (CeedInt i = 0; i < num_output_fields; i++) { 1926 CeedVector vec; 1927 1928 CeedCall(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 1929 if (vec == CEED_VECTOR_ACTIVE) { 1930 CeedEvalMode eval_mode; 1931 CeedSize rstr_flops, basis_flops; 1932 CeedElemRestriction rstr; 1933 CeedBasis basis; 1934 1935 CeedCall(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &rstr)); 1936 CeedCall(CeedElemRestrictionGetFlopsEstimate(rstr, CEED_TRANSPOSE, &rstr_flops)); 1937 CeedCall(CeedElemRestrictionDestroy(&rstr)); 1938 *flops += rstr_flops; 1939 CeedCall(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 1940 CeedCall(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 1941 CeedCall(CeedBasisGetFlopsEstimate(basis, CEED_TRANSPOSE, eval_mode, is_at_points, num_points, &basis_flops)); 1942 CeedCall(CeedBasisDestroy(&basis)); 1943 *flops += basis_flops * num_elem; 1944 } 1945 CeedCall(CeedVectorDestroy(&vec)); 1946 } 1947 } 1948 return CEED_ERROR_SUCCESS; 1949 } 1950 1951 /** 1952 @brief Get `CeedQFunction` global context for a `CeedOperator`. 1953 1954 The caller is responsible for destroying `ctx` returned from this function via @ref CeedQFunctionContextDestroy(). 1955 1956 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`. 1957 This `CeedQFunctionContext` will be destroyed if `ctx` is the only reference to this `CeedQFunctionContext`. 1958 1959 @param[in] op `CeedOperator` 1960 @param[out] ctx Variable to store `CeedQFunctionContext` 1961 1962 @return An error code: 0 - success, otherwise - failure 1963 1964 @ref Advanced 1965 **/ 1966 int CeedOperatorGetContext(CeedOperator op, CeedQFunctionContext *ctx) { 1967 bool is_composite; 1968 CeedQFunction qf; 1969 CeedQFunctionContext qf_ctx; 1970 1971 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1972 CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE, "Cannot retrieve CeedQFunctionContext for composite operator"); 1973 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1974 CeedCall(CeedQFunctionGetInnerContext(qf, &qf_ctx)); 1975 CeedCall(CeedQFunctionDestroy(&qf)); 1976 *ctx = NULL; 1977 if (qf_ctx) CeedCall(CeedQFunctionContextReferenceCopy(qf_ctx, ctx)); 1978 return CEED_ERROR_SUCCESS; 1979 } 1980 1981 /** 1982 @brief Get label for a registered `CeedQFunctionContext` field, or `NULL` if no field has been registered with this `field_name`. 1983 1984 Fields are registered via `CeedQFunctionContextRegister*()` functions (eg. @ref CeedQFunctionContextRegisterDouble()). 1985 1986 @param[in] op `CeedOperator` 1987 @param[in] field_name Name of field to retrieve label 1988 @param[out] field_label Variable to field label 1989 1990 @return An error code: 0 - success, otherwise - failure 1991 1992 @ref User 1993 **/ 1994 int CeedOperatorGetContextFieldLabel(CeedOperator op, const char *field_name, CeedContextFieldLabel *field_label) { 1995 bool is_composite, field_found = false; 1996 1997 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1998 1999 if (is_composite) { 2000 // Composite operator 2001 // -- Check if composite label already created 2002 for (CeedInt i = 0; i < op->num_context_labels; i++) { 2003 if (!strcmp(op->context_labels[i]->name, field_name)) { 2004 *field_label = op->context_labels[i]; 2005 return CEED_ERROR_SUCCESS; 2006 } 2007 } 2008 2009 // -- Create composite label if needed 2010 CeedInt num_sub; 2011 CeedOperator *sub_operators; 2012 CeedContextFieldLabel new_field_label; 2013 2014 CeedCall(CeedCalloc(1, &new_field_label)); 2015 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_sub)); 2016 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 2017 CeedCall(CeedCalloc(num_sub, &new_field_label->sub_labels)); 2018 new_field_label->num_sub_labels = num_sub; 2019 2020 for (CeedInt i = 0; i < num_sub; i++) { 2021 if (sub_operators[i]->qf->ctx) { 2022 CeedContextFieldLabel new_field_label_i; 2023 2024 CeedCall(CeedQFunctionContextGetFieldLabel(sub_operators[i]->qf->ctx, field_name, &new_field_label_i)); 2025 if (new_field_label_i) { 2026 field_found = true; 2027 new_field_label->sub_labels[i] = new_field_label_i; 2028 new_field_label->name = new_field_label_i->name; 2029 new_field_label->description = new_field_label_i->description; 2030 if (new_field_label->type && new_field_label->type != new_field_label_i->type) { 2031 // LCOV_EXCL_START 2032 CeedCall(CeedFree(&new_field_label)); 2033 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE, "Incompatible field types on sub-operator contexts. %s != %s", 2034 CeedContextFieldTypes[new_field_label->type], CeedContextFieldTypes[new_field_label_i->type]); 2035 // LCOV_EXCL_STOP 2036 } else { 2037 new_field_label->type = new_field_label_i->type; 2038 } 2039 if (new_field_label->num_values != 0 && new_field_label->num_values != new_field_label_i->num_values) { 2040 // LCOV_EXCL_START 2041 CeedCall(CeedFree(&new_field_label)); 2042 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE, 2043 "Incompatible field number of values on sub-operator contexts. %zu != %zu", new_field_label->num_values, 2044 new_field_label_i->num_values); 2045 // LCOV_EXCL_STOP 2046 } else { 2047 new_field_label->num_values = new_field_label_i->num_values; 2048 } 2049 } 2050 } 2051 } 2052 // -- Cleanup if field was found 2053 if (field_found) { 2054 *field_label = new_field_label; 2055 } else { 2056 // LCOV_EXCL_START 2057 CeedCall(CeedFree(&new_field_label->sub_labels)); 2058 CeedCall(CeedFree(&new_field_label)); 2059 *field_label = NULL; 2060 // LCOV_EXCL_STOP 2061 } 2062 } else { 2063 CeedQFunction qf; 2064 CeedQFunctionContext ctx; 2065 2066 // Single, non-composite operator 2067 CeedCall(CeedOperatorGetQFunction(op, &qf)); 2068 CeedCall(CeedQFunctionGetInnerContext(qf, &ctx)); 2069 CeedCall(CeedQFunctionDestroy(&qf)); 2070 if (ctx) { 2071 CeedCall(CeedQFunctionContextGetFieldLabel(ctx, field_name, field_label)); 2072 } else { 2073 *field_label = NULL; 2074 } 2075 } 2076 2077 // Set label in operator 2078 if (*field_label) { 2079 (*field_label)->from_op = true; 2080 2081 // Move new composite label to operator 2082 if (op->num_context_labels == 0) { 2083 CeedCall(CeedCalloc(1, &op->context_labels)); 2084 op->max_context_labels = 1; 2085 } else if (op->num_context_labels == op->max_context_labels) { 2086 CeedCall(CeedRealloc(2 * op->num_context_labels, &op->context_labels)); 2087 op->max_context_labels *= 2; 2088 } 2089 op->context_labels[op->num_context_labels] = *field_label; 2090 op->num_context_labels++; 2091 } 2092 return CEED_ERROR_SUCCESS; 2093 } 2094 2095 /** 2096 @brief Set `CeedQFunctionContext` field holding double precision values. 2097 2098 For composite operators, the values are set in all sub-operator `CeedQFunctionContext` that have a matching `field_name`. 2099 2100 @param[in,out] op `CeedOperator` 2101 @param[in] field_label Label of field to set 2102 @param[in] values Values to set 2103 2104 @return An error code: 0 - success, otherwise - failure 2105 2106 @ref User 2107 **/ 2108 int CeedOperatorSetContextDouble(CeedOperator op, CeedContextFieldLabel field_label, double *values) { 2109 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, values); 2110 } 2111 2112 /** 2113 @brief Get `CeedQFunctionContext` field holding double precision values, read-only. 2114 2115 For composite operators, the values correspond to the first sub-operator `CeedQFunctionContext` that has a matching `field_name`. 2116 2117 @param[in] op `CeedOperator` 2118 @param[in] field_label Label of field to get 2119 @param[out] num_values Number of values in the field label 2120 @param[out] values Pointer to context values 2121 2122 @return An error code: 0 - success, otherwise - failure 2123 2124 @ref User 2125 **/ 2126 int CeedOperatorGetContextDoubleRead(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const double **values) { 2127 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, num_values, values); 2128 } 2129 2130 /** 2131 @brief Restore `CeedQFunctionContext` field holding double precision values, read-only. 2132 2133 @param[in] op `CeedOperator` 2134 @param[in] field_label Label of field to restore 2135 @param[out] values Pointer to context values 2136 2137 @return An error code: 0 - success, otherwise - failure 2138 2139 @ref User 2140 **/ 2141 int CeedOperatorRestoreContextDoubleRead(CeedOperator op, CeedContextFieldLabel field_label, const double **values) { 2142 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, values); 2143 } 2144 2145 /** 2146 @brief Set `CeedQFunctionContext` field holding `int32` values. 2147 2148 For composite operators, the values are set in all sub-operator `CeedQFunctionContext` that have a matching `field_name`. 2149 2150 @param[in,out] op `CeedOperator` 2151 @param[in] field_label Label of field to set 2152 @param[in] values Values to set 2153 2154 @return An error code: 0 - success, otherwise - failure 2155 2156 @ref User 2157 **/ 2158 int CeedOperatorSetContextInt32(CeedOperator op, CeedContextFieldLabel field_label, int32_t *values) { 2159 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_INT32, values); 2160 } 2161 2162 /** 2163 @brief Get `CeedQFunctionContext` field holding `int32` values, read-only. 2164 2165 For composite operators, the values correspond to the first sub-operator `CeedQFunctionContext` that has a matching `field_name`. 2166 2167 @param[in] op `CeedOperator` 2168 @param[in] field_label Label of field to get 2169 @param[out] num_values Number of `int32` values in `values` 2170 @param[out] values Pointer to context values 2171 2172 @return An error code: 0 - success, otherwise - failure 2173 2174 @ref User 2175 **/ 2176 int CeedOperatorGetContextInt32Read(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const int32_t **values) { 2177 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_INT32, num_values, values); 2178 } 2179 2180 /** 2181 @brief Restore `CeedQFunctionContext` field holding `int32` values, read-only. 2182 2183 @param[in] op `CeedOperator` 2184 @param[in] field_label Label of field to get 2185 @param[out] values Pointer to context values 2186 2187 @return An error code: 0 - success, otherwise - failure 2188 2189 @ref User 2190 **/ 2191 int CeedOperatorRestoreContextInt32Read(CeedOperator op, CeedContextFieldLabel field_label, const int32_t **values) { 2192 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_INT32, values); 2193 } 2194 2195 /** 2196 @brief Set `CeedQFunctionContext` field holding boolean values. 2197 2198 For composite operators, the values are set in all sub-operator `CeedQFunctionContext` that have a matching `field_name`. 2199 2200 @param[in,out] op `CeedOperator` 2201 @param[in] field_label Label of field to set 2202 @param[in] values Values to set 2203 2204 @return An error code: 0 - success, otherwise - failure 2205 2206 @ref User 2207 **/ 2208 int CeedOperatorSetContextBoolean(CeedOperator op, CeedContextFieldLabel field_label, bool *values) { 2209 return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_BOOL, values); 2210 } 2211 2212 /** 2213 @brief Get `CeedQFunctionContext` field holding boolean values, read-only. 2214 2215 For composite operators, the values correspond to the first sub-operator `CeedQFunctionContext` that has a matching `field_name`. 2216 2217 @param[in] op `CeedOperator` 2218 @param[in] field_label Label of field to get 2219 @param[out] num_values Number of boolean values in `values` 2220 @param[out] values Pointer to context values 2221 2222 @return An error code: 0 - success, otherwise - failure 2223 2224 @ref User 2225 **/ 2226 int CeedOperatorGetContextBooleanRead(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const bool **values) { 2227 return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_BOOL, num_values, values); 2228 } 2229 2230 /** 2231 @brief Restore `CeedQFunctionContext` field holding boolean values, read-only. 2232 2233 @param[in] op `CeedOperator` 2234 @param[in] field_label Label of field to get 2235 @param[out] values Pointer to context values 2236 2237 @return An error code: 0 - success, otherwise - failure 2238 2239 @ref User 2240 **/ 2241 int CeedOperatorRestoreContextBooleanRead(CeedOperator op, CeedContextFieldLabel field_label, const bool **values) { 2242 return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_BOOL, values); 2243 } 2244 2245 /** 2246 @brief Apply `CeedOperator` to a `CeedVector`. 2247 2248 This computes the action of the operator on the specified (active) input, yielding its (active) output. 2249 All inputs and outputs must be specified using @ref CeedOperatorSetField(). 2250 2251 @note Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 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 store 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 CeedOperatorApply(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 && op->ApplyComposite) { 2269 // Composite Operator 2270 CeedCall(op->ApplyComposite(op, in, out, request)); 2271 } else if (!is_composite && op->Apply) { 2272 // Standard Operator 2273 CeedCall(op->Apply(op, in, out, request)); 2274 } else { 2275 // Standard or composite, default to zeroing out and calling ApplyAddActive 2276 // Zero active output 2277 if (out != CEED_VECTOR_NONE) CeedCall(CeedVectorSetValue(out, 0.0)); 2278 2279 // ApplyAddActive 2280 CeedCall(CeedOperatorApplyAddActive(op, in, out, request)); 2281 } 2282 return CEED_ERROR_SUCCESS; 2283 } 2284 2285 /** 2286 @brief Apply `CeedOperator` to a `CeedVector` and add result to output `CeedVector`. 2287 2288 This computes the action of the operator on the specified (active) input, yielding its (active) output. 2289 All inputs and outputs must be specified using @ref CeedOperatorSetField(). 2290 2291 @note Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2292 @warning This function adds into ALL outputs, including passive outputs. To only add into the active output, use `CeedOperatorApplyAddActive()`. 2293 @see `CeedOperatorApplyAddActive()` 2294 2295 @param[in] op `CeedOperator` to apply 2296 @param[in] in `CeedVector` containing input state or @ref CEED_VECTOR_NONE if there are no active inputs 2297 @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 2298 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2299 2300 @return An error code: 0 - success, otherwise - failure 2301 2302 @ref User 2303 **/ 2304 int CeedOperatorApplyAdd(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) { 2305 bool is_composite; 2306 2307 CeedCall(CeedOperatorCheckReady(op)); 2308 2309 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2310 if (is_composite) { 2311 // Composite Operator 2312 if (op->ApplyAddComposite) { 2313 CeedCall(op->ApplyAddComposite(op, in, out, request)); 2314 } else { 2315 CeedInt num_suboperators; 2316 CeedOperator *sub_operators; 2317 2318 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 2319 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 2320 for (CeedInt i = 0; i < num_suboperators; i++) { 2321 CeedCall(CeedOperatorApplyAdd(sub_operators[i], in, out, request)); 2322 } 2323 } 2324 } else if (op->num_elem > 0) { 2325 // Standard Operator 2326 CeedCall(op->ApplyAdd(op, in, out, request)); 2327 } 2328 return CEED_ERROR_SUCCESS; 2329 } 2330 2331 /** 2332 @brief Apply `CeedOperator` to a `CeedVector` and add result to output `CeedVector`. Only sums into active outputs, overwrites passive outputs. 2333 2334 This computes the action of the operator on the specified (active) input, yielding its (active) output. 2335 All inputs and outputs must be specified using @ref CeedOperatorSetField(). 2336 2337 @note Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2338 2339 @param[in] op `CeedOperator` to apply 2340 @param[in] in `CeedVector` containing input state or @ref CEED_VECTOR_NONE if there are no active inputs 2341 @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 2342 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2343 2344 @return An error code: 0 - success, otherwise - failure 2345 2346 @ref User 2347 **/ 2348 int CeedOperatorApplyAddActive(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) { 2349 bool is_composite; 2350 2351 CeedCall(CeedOperatorCheckReady(op)); 2352 2353 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2354 if (is_composite) { 2355 // Composite Operator 2356 CeedInt num_suboperators; 2357 CeedOperator *sub_operators; 2358 2359 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 2360 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 2361 2362 // Zero all output vectors 2363 for (CeedInt i = 0; i < num_suboperators; i++) { 2364 CeedInt num_output_fields; 2365 CeedOperatorField *output_fields; 2366 2367 CeedCall(CeedOperatorGetFields(sub_operators[i], NULL, NULL, &num_output_fields, &output_fields)); 2368 for (CeedInt j = 0; j < num_output_fields; j++) { 2369 CeedVector vec; 2370 2371 CeedCall(CeedOperatorFieldGetVector(output_fields[j], &vec)); 2372 if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) CeedCall(CeedVectorSetValue(vec, 0.0)); 2373 CeedCall(CeedVectorDestroy(&vec)); 2374 } 2375 } 2376 // ApplyAdd 2377 CeedCall(CeedOperatorApplyAdd(op, in, out, request)); 2378 } else { 2379 // Standard Operator 2380 CeedInt num_output_fields; 2381 CeedOperatorField *output_fields; 2382 2383 CeedCall(CeedOperatorGetFields(op, NULL, NULL, &num_output_fields, &output_fields)); 2384 // Zero all output vectors 2385 for (CeedInt i = 0; i < num_output_fields; i++) { 2386 CeedVector vec; 2387 2388 CeedCall(CeedOperatorFieldGetVector(output_fields[i], &vec)); 2389 if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) CeedCall(CeedVectorSetValue(vec, 0.0)); 2390 CeedCall(CeedVectorDestroy(&vec)); 2391 } 2392 // ApplyAdd 2393 CeedCall(CeedOperatorApplyAdd(op, in, out, request)); 2394 } 2395 return CEED_ERROR_SUCCESS; 2396 } 2397 2398 /** 2399 @brief Destroy temporary assembly data associated with a `CeedOperator` 2400 2401 @param[in,out] op `CeedOperator` whose assembly data to destroy 2402 2403 @return An error code: 0 - success, otherwise - failure 2404 2405 @ref User 2406 **/ 2407 int CeedOperatorAssemblyDataStrip(CeedOperator op) { 2408 bool is_composite; 2409 2410 CeedCall(CeedQFunctionAssemblyDataDestroy(&op->qf_assembled)); 2411 CeedCall(CeedOperatorAssemblyDataDestroy(&op->op_assembled)); 2412 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2413 if (is_composite) { 2414 CeedInt num_suboperators; 2415 CeedOperator *sub_operators; 2416 2417 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 2418 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 2419 for (CeedInt i = 0; i < num_suboperators; i++) { 2420 CeedCall(CeedQFunctionAssemblyDataDestroy(&sub_operators[i]->qf_assembled)); 2421 CeedCall(CeedOperatorAssemblyDataDestroy(&sub_operators[i]->op_assembled)); 2422 } 2423 } 2424 return CEED_ERROR_SUCCESS; 2425 } 2426 2427 /** 2428 @brief Destroy a `CeedOperator` 2429 2430 @param[in,out] op `CeedOperator` to destroy 2431 2432 @return An error code: 0 - success, otherwise - failure 2433 2434 @ref User 2435 **/ 2436 int CeedOperatorDestroy(CeedOperator *op) { 2437 if (!*op || CeedObjectDereference((CeedObject)*op) > 0) { 2438 *op = NULL; 2439 return CEED_ERROR_SUCCESS; 2440 } 2441 // Backend destroy 2442 if ((*op)->Destroy) { 2443 CeedCall((*op)->Destroy(*op)); 2444 } 2445 // Free fields 2446 for (CeedInt i = 0; i < (*op)->num_fields; i++) { 2447 if ((*op)->input_fields[i]) { 2448 if ((*op)->input_fields[i]->elem_rstr != CEED_ELEMRESTRICTION_NONE) { 2449 CeedCall(CeedElemRestrictionDestroy(&(*op)->input_fields[i]->elem_rstr)); 2450 } 2451 if ((*op)->input_fields[i]->basis != CEED_BASIS_NONE) { 2452 CeedCall(CeedBasisDestroy(&(*op)->input_fields[i]->basis)); 2453 } 2454 if ((*op)->input_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->input_fields[i]->vec != CEED_VECTOR_NONE) { 2455 CeedCall(CeedVectorDestroy(&(*op)->input_fields[i]->vec)); 2456 } 2457 CeedCall(CeedFree(&(*op)->input_fields[i]->field_name)); 2458 CeedCall(CeedFree(&(*op)->input_fields[i])); 2459 } 2460 } 2461 for (CeedInt i = 0; i < (*op)->num_fields; i++) { 2462 if ((*op)->output_fields[i]) { 2463 CeedCall(CeedElemRestrictionDestroy(&(*op)->output_fields[i]->elem_rstr)); 2464 if ((*op)->output_fields[i]->basis != CEED_BASIS_NONE) { 2465 CeedCall(CeedBasisDestroy(&(*op)->output_fields[i]->basis)); 2466 } 2467 if ((*op)->output_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->output_fields[i]->vec != CEED_VECTOR_NONE) { 2468 CeedCall(CeedVectorDestroy(&(*op)->output_fields[i]->vec)); 2469 } 2470 CeedCall(CeedFree(&(*op)->output_fields[i]->field_name)); 2471 CeedCall(CeedFree(&(*op)->output_fields[i])); 2472 } 2473 } 2474 CeedCall(CeedFree(&(*op)->input_fields)); 2475 CeedCall(CeedFree(&(*op)->output_fields)); 2476 // Destroy AtPoints data 2477 CeedCall(CeedVectorDestroy(&(*op)->point_coords)); 2478 CeedCall(CeedElemRestrictionDestroy(&(*op)->rstr_points)); 2479 CeedCall(CeedElemRestrictionDestroy(&(*op)->first_points_rstr)); 2480 // Destroy assembly data (must happen before destroying sub_operators) 2481 CeedCall(CeedOperatorAssemblyDataStrip(*op)); 2482 // Destroy sub_operators 2483 for (CeedInt i = 0; i < (*op)->num_suboperators; i++) { 2484 if ((*op)->sub_operators[i]) { 2485 CeedCall(CeedOperatorDestroy(&(*op)->sub_operators[i])); 2486 } 2487 } 2488 CeedCall(CeedFree(&(*op)->sub_operators)); 2489 CeedCall(CeedQFunctionDestroy(&(*op)->qf)); 2490 CeedCall(CeedQFunctionDestroy(&(*op)->dqf)); 2491 CeedCall(CeedQFunctionDestroy(&(*op)->dqfT)); 2492 // Destroy any composite labels 2493 if ((*op)->is_composite) { 2494 for (CeedInt i = 0; i < (*op)->num_context_labels; i++) { 2495 CeedCall(CeedFree(&(*op)->context_labels[i]->sub_labels)); 2496 CeedCall(CeedFree(&(*op)->context_labels[i])); 2497 } 2498 } 2499 CeedCall(CeedFree(&(*op)->context_labels)); 2500 2501 // Destroy fallback 2502 CeedCall(CeedOperatorDestroy(&(*op)->op_fallback)); 2503 2504 CeedCall(CeedFree(&(*op)->name)); 2505 CeedCall(CeedObjectDestroy_Private(&(*op)->obj)); 2506 CeedCall(CeedFree(op)); 2507 return CEED_ERROR_SUCCESS; 2508 } 2509 2510 /// @} 2511