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