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