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