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