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