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) { 2163 // Composite Operator 2164 if (op->ApplyComposite) { 2165 CeedCall(op->ApplyComposite(op, in, out, request)); 2166 } else { 2167 CeedInt num_suboperators; 2168 CeedOperator *sub_operators; 2169 2170 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 2171 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 2172 2173 // Zero all output vectors 2174 if (out != CEED_VECTOR_NONE) CeedCall(CeedVectorSetValue(out, 0.0)); 2175 for (CeedInt i = 0; i < num_suboperators; i++) { 2176 CeedInt num_output_fields; 2177 CeedOperatorField *output_fields; 2178 2179 CeedCall(CeedOperatorGetFields(sub_operators[i], NULL, NULL, &num_output_fields, &output_fields)); 2180 for (CeedInt j = 0; j < num_output_fields; j++) { 2181 CeedVector vec; 2182 2183 CeedCall(CeedOperatorFieldGetVector(output_fields[j], &vec)); 2184 if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) { 2185 CeedCall(CeedVectorSetValue(vec, 0.0)); 2186 } 2187 CeedCall(CeedVectorDestroy(&vec)); 2188 } 2189 } 2190 // ApplyAdd 2191 CeedCall(CeedOperatorApplyAdd(op, in, out, request)); 2192 } 2193 } else { 2194 // Standard Operator 2195 if (op->Apply) { 2196 CeedCall(op->Apply(op, in, out, request)); 2197 } else { 2198 CeedInt num_output_fields; 2199 CeedOperatorField *output_fields; 2200 2201 CeedCall(CeedOperatorGetFields(op, NULL, NULL, &num_output_fields, &output_fields)); 2202 // Zero all output vectors 2203 for (CeedInt i = 0; i < num_output_fields; i++) { 2204 bool is_active; 2205 CeedVector vec; 2206 2207 CeedCall(CeedOperatorFieldGetVector(output_fields[i], &vec)); 2208 is_active = vec == CEED_VECTOR_ACTIVE; 2209 if (is_active) vec = out; 2210 if (vec != CEED_VECTOR_NONE) CeedCall(CeedVectorSetValue(vec, 0.0)); 2211 if (!is_active) CeedCall(CeedVectorDestroy(&vec)); 2212 } 2213 // Apply 2214 if (op->num_elem > 0) CeedCall(op->ApplyAdd(op, in, out, request)); 2215 } 2216 } 2217 return CEED_ERROR_SUCCESS; 2218 } 2219 2220 /** 2221 @brief Apply `CeedOperator` to a `CeedVector` and add result to output `CeedVector`. 2222 2223 This computes the action of the operator on the specified (active) input, yielding its (active) output. 2224 All inputs and outputs must be specified using @ref CeedOperatorSetField(). 2225 2226 @param[in] op `CeedOperator` to apply 2227 @param[in] in `CeedVector` containing input state or @ref CEED_VECTOR_NONE if there are no active inputs 2228 @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 2229 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2230 2231 @return An error code: 0 - success, otherwise - failure 2232 2233 @ref User 2234 **/ 2235 int CeedOperatorApplyAdd(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) { 2236 bool is_composite; 2237 2238 CeedCall(CeedOperatorCheckReady(op)); 2239 2240 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2241 if (is_composite) { 2242 // Composite Operator 2243 if (op->ApplyAddComposite) { 2244 CeedCall(op->ApplyAddComposite(op, in, out, request)); 2245 } else { 2246 CeedInt num_suboperators; 2247 CeedOperator *sub_operators; 2248 2249 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 2250 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 2251 for (CeedInt i = 0; i < num_suboperators; i++) { 2252 CeedCall(CeedOperatorApplyAdd(sub_operators[i], in, out, request)); 2253 } 2254 } 2255 } else if (op->num_elem > 0) { 2256 // Standard Operator 2257 CeedCall(op->ApplyAdd(op, in, out, request)); 2258 } 2259 return CEED_ERROR_SUCCESS; 2260 } 2261 2262 /** 2263 @brief Destroy temporary assembly data associated with a `CeedOperator` 2264 2265 @param[in,out] op `CeedOperator` whose assembly data to destroy 2266 2267 @return An error code: 0 - success, otherwise - failure 2268 2269 @ref User 2270 **/ 2271 int CeedOperatorAssemblyDataStrip(CeedOperator op) { 2272 bool is_composite; 2273 2274 CeedCall(CeedQFunctionAssemblyDataDestroy(&op->qf_assembled)); 2275 CeedCall(CeedOperatorAssemblyDataDestroy(&op->op_assembled)); 2276 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2277 if (is_composite) { 2278 CeedInt num_suboperators; 2279 CeedOperator *sub_operators; 2280 2281 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 2282 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 2283 for (CeedInt i = 0; i < num_suboperators; i++) { 2284 CeedCall(CeedQFunctionAssemblyDataDestroy(&sub_operators[i]->qf_assembled)); 2285 CeedCall(CeedOperatorAssemblyDataDestroy(&sub_operators[i]->op_assembled)); 2286 } 2287 } 2288 return CEED_ERROR_SUCCESS; 2289 } 2290 2291 /** 2292 @brief Destroy a `CeedOperator` 2293 2294 @param[in,out] op `CeedOperator` to destroy 2295 2296 @return An error code: 0 - success, otherwise - failure 2297 2298 @ref User 2299 **/ 2300 int CeedOperatorDestroy(CeedOperator *op) { 2301 if (!*op || --(*op)->ref_count > 0) { 2302 *op = NULL; 2303 return CEED_ERROR_SUCCESS; 2304 } 2305 // Backend destroy 2306 if ((*op)->Destroy) { 2307 CeedCall((*op)->Destroy(*op)); 2308 } 2309 // Free fields 2310 for (CeedInt i = 0; i < (*op)->num_fields; i++) { 2311 if ((*op)->input_fields[i]) { 2312 if ((*op)->input_fields[i]->elem_rstr != CEED_ELEMRESTRICTION_NONE) { 2313 CeedCall(CeedElemRestrictionDestroy(&(*op)->input_fields[i]->elem_rstr)); 2314 } 2315 if ((*op)->input_fields[i]->basis != CEED_BASIS_NONE) { 2316 CeedCall(CeedBasisDestroy(&(*op)->input_fields[i]->basis)); 2317 } 2318 if ((*op)->input_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->input_fields[i]->vec != CEED_VECTOR_NONE) { 2319 CeedCall(CeedVectorDestroy(&(*op)->input_fields[i]->vec)); 2320 } 2321 CeedCall(CeedFree(&(*op)->input_fields[i]->field_name)); 2322 CeedCall(CeedFree(&(*op)->input_fields[i])); 2323 } 2324 } 2325 for (CeedInt i = 0; i < (*op)->num_fields; i++) { 2326 if ((*op)->output_fields[i]) { 2327 CeedCall(CeedElemRestrictionDestroy(&(*op)->output_fields[i]->elem_rstr)); 2328 if ((*op)->output_fields[i]->basis != CEED_BASIS_NONE) { 2329 CeedCall(CeedBasisDestroy(&(*op)->output_fields[i]->basis)); 2330 } 2331 if ((*op)->output_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->output_fields[i]->vec != CEED_VECTOR_NONE) { 2332 CeedCall(CeedVectorDestroy(&(*op)->output_fields[i]->vec)); 2333 } 2334 CeedCall(CeedFree(&(*op)->output_fields[i]->field_name)); 2335 CeedCall(CeedFree(&(*op)->output_fields[i])); 2336 } 2337 } 2338 CeedCall(CeedFree(&(*op)->input_fields)); 2339 CeedCall(CeedFree(&(*op)->output_fields)); 2340 // Destroy AtPoints data 2341 CeedCall(CeedVectorDestroy(&(*op)->point_coords)); 2342 CeedCall(CeedElemRestrictionDestroy(&(*op)->rstr_points)); 2343 CeedCall(CeedElemRestrictionDestroy(&(*op)->first_points_rstr)); 2344 // Destroy assembly data (must happen before destroying sub_operators) 2345 CeedCall(CeedOperatorAssemblyDataStrip(*op)); 2346 // Destroy sub_operators 2347 for (CeedInt i = 0; i < (*op)->num_suboperators; i++) { 2348 if ((*op)->sub_operators[i]) { 2349 CeedCall(CeedOperatorDestroy(&(*op)->sub_operators[i])); 2350 } 2351 } 2352 CeedCall(CeedFree(&(*op)->sub_operators)); 2353 CeedCall(CeedQFunctionDestroy(&(*op)->qf)); 2354 CeedCall(CeedQFunctionDestroy(&(*op)->dqf)); 2355 CeedCall(CeedQFunctionDestroy(&(*op)->dqfT)); 2356 // Destroy any composite labels 2357 if ((*op)->is_composite) { 2358 for (CeedInt i = 0; i < (*op)->num_context_labels; i++) { 2359 CeedCall(CeedFree(&(*op)->context_labels[i]->sub_labels)); 2360 CeedCall(CeedFree(&(*op)->context_labels[i])); 2361 } 2362 } 2363 CeedCall(CeedFree(&(*op)->context_labels)); 2364 2365 // Destroy fallback 2366 CeedCall(CeedOperatorDestroy(&(*op)->op_fallback)); 2367 2368 CeedCall(CeedFree(&(*op)->name)); 2369 CeedCall(CeedDestroy(&(*op)->ceed)); 2370 CeedCall(CeedFree(op)); 2371 return CEED_ERROR_SUCCESS; 2372 } 2373 2374 /// @} 2375