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