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