xref: /libCEED/interface/ceed-operator.c (revision ed094490f53e580908aa80e9fe815a6fd76d7526)

// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and other CEED contributors.
// All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
//
// SPDX-License-Identifier: BSD-2-Clause
//
// This file is part of CEED:  http://github.com/ceed

#include <ceed-impl.h>
#include <ceed.h>
#include <ceed/backend.h>
#include <stdbool.h>
#include <stdio.h>
#include <string.h>

/// @file
/// Implementation of CeedOperator interfaces

/// ----------------------------------------------------------------------------
/// CeedOperator Library Internal Functions
/// ----------------------------------------------------------------------------
/// @addtogroup CeedOperatorDeveloper
/// @{

/**
  @brief Check if a `CeedOperator` Field matches the `CeedQFunction` Field

  @param[in] ceed     `Ceed` object for error handling
  @param[in] qf_field `CeedQFunction` Field matching `CeedOperator` Field
  @param[in] rstr     `CeedOperator` Field `CeedElemRestriction`
  @param[in] basis    `CeedOperator` Field `CeedBasis`

  @return An error code: 0 - success, otherwise - failure

  @ref Developer
**/
static int CeedOperatorCheckField(Ceed ceed, CeedQFunctionField qf_field, CeedElemRestriction rstr, CeedBasis basis) {
  const char  *field_name;
  CeedInt      dim = 1, num_comp = 1, q_comp = 1, rstr_num_comp = 1, size;
  CeedEvalMode eval_mode;

  // Field data
  CeedCall(CeedQFunctionFieldGetData(qf_field, &field_name, &size, &eval_mode));

  // Restriction
  CeedCheck((rstr == CEED_ELEMRESTRICTION_NONE) == (eval_mode == CEED_EVAL_WEIGHT), ceed, CEED_ERROR_INCOMPATIBLE,
            "CEED_ELEMRESTRICTION_NONE and CEED_EVAL_WEIGHT must be used together.");
  if (rstr != CEED_ELEMRESTRICTION_NONE) {
    CeedCall(CeedElemRestrictionGetNumComponents(rstr, &rstr_num_comp));
  }
  // Basis
  CeedCheck((basis == CEED_BASIS_NONE) == (eval_mode == CEED_EVAL_NONE), ceed, CEED_ERROR_INCOMPATIBLE,
            "CEED_BASIS_NONE and CEED_EVAL_NONE must be used together.");
  if (basis != CEED_BASIS_NONE) {
    CeedCall(CeedBasisGetDimension(basis, &dim));
    CeedCall(CeedBasisGetNumComponents(basis, &num_comp));
    CeedCall(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp));
    CeedCheck(rstr == CEED_ELEMRESTRICTION_NONE || rstr_num_comp == num_comp, ceed, CEED_ERROR_DIMENSION,
              "Field '%s' of size %" CeedInt_FMT " and EvalMode %s: CeedElemRestriction has %" CeedInt_FMT
              " components, but CeedBasis has %" CeedInt_FMT " components",
              field_name, size, CeedEvalModes[eval_mode], rstr_num_comp, num_comp);
  }
  // Field size
  switch (eval_mode) {
    case CEED_EVAL_NONE:
      CeedCheck(size == rstr_num_comp, ceed, CEED_ERROR_DIMENSION,
                "Field '%s' of size %" CeedInt_FMT " and EvalMode %s: CeedElemRestriction has %" CeedInt_FMT " components", field_name, size,
                CeedEvalModes[eval_mode], rstr_num_comp);
      break;
    case CEED_EVAL_INTERP:
    case CEED_EVAL_GRAD:
    case CEED_EVAL_DIV:
    case CEED_EVAL_CURL:
      CeedCheck(size == num_comp * q_comp, ceed, CEED_ERROR_DIMENSION,
                "Field '%s' of size %" CeedInt_FMT " and EvalMode %s: CeedElemRestriction/Basis has %" CeedInt_FMT " components", field_name, size,
                CeedEvalModes[eval_mode], num_comp * q_comp);
      break;
    case CEED_EVAL_WEIGHT:
      // No additional checks required
      break;
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief View a field of a `CeedOperator`

  @param[in] op_field     `CeedOperator` Field to view
  @param[in] qf_field     `CeedQFunction` Field (carries field name)
  @param[in] field_number Number of field being viewed
  @param[in] sub          true indicates sub-operator, which increases indentation; false for top-level operator
  @param[in] input        true for an input field; false for output field
  @param[in] stream       Stream to view to, e.g., `stdout`

  @return An error code: 0 - success, otherwise - failure

  @ref Utility
**/
static int CeedOperatorFieldView(CeedOperatorField op_field, CeedQFunctionField qf_field, CeedInt field_number, bool sub, bool input, FILE *stream) {
  const char  *pre    = sub ? "  " : "";
  const char  *in_out = input ? "Input" : "Output";
  const char  *field_name;
  CeedInt      size;
  CeedEvalMode eval_mode;
  CeedVector   vec;
  CeedBasis    basis;

  // Field data
  CeedCall(CeedQFunctionFieldGetData(qf_field, &field_name, &size, &eval_mode));
  CeedCall(CeedOperatorFieldGetData(op_field, NULL, NULL, &basis, &vec));

  fprintf(stream,
          "%s    %s field %" CeedInt_FMT
          ":\n"
          "%s      Name: \"%s\"\n",
          pre, in_out, field_number, pre, field_name);
  fprintf(stream, "%s      Size: %" CeedInt_FMT "\n", pre, size);
  fprintf(stream, "%s      EvalMode: %s\n", pre, CeedEvalModes[eval_mode]);
  if (basis == CEED_BASIS_NONE) fprintf(stream, "%s      No basis\n", pre);
  if (vec == CEED_VECTOR_ACTIVE) fprintf(stream, "%s      Active vector\n", pre);
  else if (vec == CEED_VECTOR_NONE) fprintf(stream, "%s      No vector\n", pre);

  CeedCall(CeedVectorDestroy(&vec));
  CeedCall(CeedBasisDestroy(&basis));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief View a single `CeedOperator`

  @param[in] op     `CeedOperator` to view
  @param[in] sub    Boolean flag for sub-operator
  @param[in] stream Stream to write; typically `stdout` or a file

  @return Error code: 0 - success, otherwise - failure

  @ref Utility
**/
int CeedOperatorSingleView(CeedOperator op, bool sub, FILE *stream) {
  bool                is_at_points;
  const char         *pre = sub ? "  " : "";
  CeedInt             num_elem, num_qpts, total_fields = 0, num_input_fields, num_output_fields;
  CeedQFunction       qf;
  CeedQFunctionField *qf_input_fields, *qf_output_fields;
  CeedOperatorField  *op_input_fields, *op_output_fields;

  CeedCall(CeedOperatorIsAtPoints(op, &is_at_points));
  CeedCall(CeedOperatorGetNumElements(op, &num_elem));
  CeedCall(CeedOperatorGetNumQuadraturePoints(op, &num_qpts));
  CeedCall(CeedOperatorGetNumArgs(op, &total_fields));
  CeedCall(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
  CeedCall(CeedOperatorGetQFunction(op, &qf));
  CeedCall(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
  CeedCall(CeedQFunctionDestroy(&qf));

  if (is_at_points) {
    CeedInt             max_points = 0;
    CeedElemRestriction rstr_points;

    CeedCall(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL));
    CeedCall(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_points));
    fprintf(stream, "%s  %" CeedInt_FMT " elements with %" CeedInt_FMT " max points each\n", pre, num_elem, max_points);
    CeedCall(CeedElemRestrictionDestroy(&rstr_points));
  } else {
    fprintf(stream, "%s  %" CeedInt_FMT " elements with %" CeedInt_FMT " quadrature points each\n", pre, num_elem, num_qpts);
  }
  fprintf(stream, "%s  %" CeedInt_FMT " field%s\n", pre, total_fields, total_fields > 1 ? "s" : "");
  fprintf(stream, "%s  %" CeedInt_FMT " input field%s:\n", pre, num_input_fields, num_input_fields > 1 ? "s" : "");
  for (CeedInt i = 0; i < num_input_fields; i++) {
    CeedCall(CeedOperatorFieldView(op_input_fields[i], qf_input_fields[i], i, sub, 1, stream));
  }
  fprintf(stream, "%s  %" CeedInt_FMT " output field%s:\n", pre, num_output_fields, num_output_fields > 1 ? "s" : "");
  for (CeedInt i = 0; i < num_output_fields; i++) {
    CeedCall(CeedOperatorFieldView(op_output_fields[i], qf_output_fields[i], i, sub, 0, stream));
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Find the active input vector `CeedBasis` for a non-composite `CeedOperator`.

  Note: Caller is responsible for destroying the `active_basis` with @ref CeedBasisDestroy().

  @param[in]  op           `CeedOperator` to find active `CeedBasis` for
  @param[out] active_basis `CeedBasis` for active input vector or `NULL` for composite operator

  @return An error code: 0 - success, otherwise - failure

  @ref Developer
**/
int CeedOperatorGetActiveBasis(CeedOperator op, CeedBasis *active_basis) {
  CeedCall(CeedOperatorGetActiveBases(op, active_basis, NULL));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Find the active input and output vector `CeedBasis` for a non-composite `CeedOperator`.

  Note: Caller is responsible for destroying the bases with @ref CeedBasisDestroy().

  @param[in]  op                  `CeedOperator` to find active `CeedBasis` for
  @param[out] active_input_basis  `CeedBasis` for active input vector or `NULL` for composite operator
  @param[out] active_output_basis `CeedBasis` for active output vector or `NULL` for composite operator

  @return An error code: 0 - success, otherwise - failure

  @ref Developer
**/
int CeedOperatorGetActiveBases(CeedOperator op, CeedBasis *active_input_basis, CeedBasis *active_output_basis) {
  bool               is_composite;
  CeedInt            num_input_fields, num_output_fields;
  CeedOperatorField *op_input_fields, *op_output_fields;

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCall(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));

  if (active_input_basis) {
    *active_input_basis = NULL;
    if (!is_composite) {
      for (CeedInt i = 0; i < num_input_fields; i++) {
        CeedVector vec;

        CeedCall(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
        if (vec == CEED_VECTOR_ACTIVE) {
          CeedBasis basis;

          CeedCall(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
          CeedCheck(!*active_input_basis || *active_input_basis == basis, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR,
                    "Multiple active input CeedBases found");
          if (!*active_input_basis) CeedCall(CeedBasisReferenceCopy(basis, active_input_basis));
          CeedCall(CeedBasisDestroy(&basis));
        }
        CeedCall(CeedVectorDestroy(&vec));
      }
      CeedCheck(*active_input_basis, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE, "No active input CeedBasis found");
    }
  }
  if (active_output_basis) {
    *active_output_basis = NULL;
    if (!is_composite) {
      for (CeedInt i = 0; i < num_output_fields; i++) {
        CeedVector vec;

        CeedCall(CeedOperatorFieldGetVector(op_output_fields[i], &vec));
        if (vec == CEED_VECTOR_ACTIVE) {
          CeedBasis basis;

          CeedCall(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
          CeedCheck(!*active_output_basis || *active_output_basis == basis, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR,
                    "Multiple active output CeedBases found");
          if (!*active_output_basis) CeedCall(CeedBasisReferenceCopy(basis, active_output_basis));
          CeedCall(CeedBasisDestroy(&basis));
        }
        CeedCall(CeedVectorDestroy(&vec));
      }
      CeedCheck(*active_output_basis, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE, "No active output CeedBasis found");
    }
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Find the active vector `CeedElemRestriction` for a non-composite `CeedOperator`.

  Note: Caller is responsible for destroying the `active_rstr` with @ref CeedElemRestrictionDestroy().

  @param[in]  op          `CeedOperator` to find active `CeedElemRestriction` for
  @param[out] active_rstr `CeedElemRestriction` for active input vector or NULL for composite operator

  @return An error code: 0 - success, otherwise - failure

  @ref Utility
**/
int CeedOperatorGetActiveElemRestriction(CeedOperator op, CeedElemRestriction *active_rstr) {
  CeedCall(CeedOperatorGetActiveElemRestrictions(op, active_rstr, NULL));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Find the active input and output vector `CeedElemRestriction` for a non-composite `CeedOperator`.

  Note: Caller is responsible for destroying the restrictions with @ref CeedElemRestrictionDestroy().

  @param[in]  op                 `CeedOperator` to find active `CeedElemRestriction` for
  @param[out] active_input_rstr  `CeedElemRestriction` for active input vector or NULL for composite operator
  @param[out] active_output_rstr `CeedElemRestriction` for active output vector or NULL for composite operator

  @return An error code: 0 - success, otherwise - failure

  @ref Utility
**/
int CeedOperatorGetActiveElemRestrictions(CeedOperator op, CeedElemRestriction *active_input_rstr, CeedElemRestriction *active_output_rstr) {
  bool               is_composite;
  CeedInt            num_input_fields, num_output_fields;
  CeedOperatorField *op_input_fields, *op_output_fields;

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCall(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));

  if (active_input_rstr) {
    *active_input_rstr = NULL;
    if (!is_composite) {
      for (CeedInt i = 0; i < num_input_fields; i++) {
        CeedVector vec;

        CeedCall(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
        if (vec == CEED_VECTOR_ACTIVE) {
          CeedElemRestriction rstr;

          CeedCall(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &rstr));
          CeedCheck(!*active_input_rstr || *active_input_rstr == rstr, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR,
                    "Multiple active input CeedElemRestrictions found");
          if (!*active_input_rstr) CeedCall(CeedElemRestrictionReferenceCopy(rstr, active_input_rstr));
          CeedCall(CeedElemRestrictionDestroy(&rstr));
        }
        CeedCall(CeedVectorDestroy(&vec));
      }
      CeedCheck(*active_input_rstr, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE, "No active input CeedElemRestriction found");
    }
  }
  if (active_output_rstr) {
    *active_output_rstr = NULL;
    if (!is_composite) {
      for (CeedInt i = 0; i < num_output_fields; i++) {
        CeedVector vec;

        CeedCall(CeedOperatorFieldGetVector(op_output_fields[i], &vec));
        if (vec == CEED_VECTOR_ACTIVE) {
          CeedElemRestriction rstr;

          CeedCall(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &rstr));
          CeedCheck(!*active_output_rstr || *active_output_rstr == rstr, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR,
                    "Multiple active output CeedElemRestrictions found");
          if (!*active_output_rstr) CeedCall(CeedElemRestrictionReferenceCopy(rstr, active_output_rstr));
          CeedCall(CeedElemRestrictionDestroy(&rstr));
        }
        CeedCall(CeedVectorDestroy(&vec));
      }
      CeedCheck(*active_output_rstr, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE, "No active output CeedElemRestriction found");
    }
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Set `CeedQFunctionContext` field values of the specified type.

  For composite operators, the value is set in all sub-operator `CeedQFunctionContext` that have a matching `field_name`.
  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.

  @param[in,out] op          `CeedOperator`
  @param[in]     field_label Label of field to set
  @param[in]     field_type  Type of field to set
  @param[in]     values      Values to set

  @return An error code: 0 - success, otherwise - failure

  @ref Developer
**/
static int CeedOperatorContextSetGeneric(CeedOperator op, CeedContextFieldLabel field_label, CeedContextFieldType field_type, void *values) {
  bool is_composite = false;

  CeedCheck(field_label, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, "Invalid field label");

  // Check if field_label and op correspond
  if (field_label->from_op) {
    CeedInt index = -1;

    for (CeedInt i = 0; i < op->num_context_labels; i++) {
      if (op->context_labels[i] == field_label) index = i;
    }
    CeedCheck(index != -1, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, "ContextFieldLabel does not correspond to the operator");
  }

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  if (is_composite) {
    CeedInt       num_sub;
    CeedOperator *sub_operators;

    CeedCall(CeedOperatorCompositeGetNumSub(op, &num_sub));
    CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators));
    CeedCheck(num_sub == field_label->num_sub_labels, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED,
              "Composite operator modified after ContextFieldLabel created");

    for (CeedInt i = 0; i < num_sub; i++) {
      CeedQFunctionContext ctx;

      CeedCall(CeedOperatorGetContext(sub_operators[i], &ctx));
      // Try every sub-operator, ok if some sub-operators do not have field
      if (ctx && field_label->sub_labels[i]) {
        CeedCall(CeedQFunctionContextSetGeneric(ctx, field_label->sub_labels[i], field_type, values));
      }
      CeedCall(CeedQFunctionContextDestroy(&ctx));
    }
  } else {
    CeedQFunctionContext ctx;

    CeedCall(CeedOperatorGetContext(op, &ctx));
    CeedCheck(ctx, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, "QFunction does not have context data");
    CeedCall(CeedQFunctionContextSetGeneric(ctx, field_label, field_type, values));
    CeedCall(CeedQFunctionContextDestroy(&ctx));
  }
  CeedCall(CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(op, true));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get `CeedQFunctionContext` field values of the specified type, read-only.

  For composite operators, the values retrieved are for the first sub-operator `CeedQFunctionContext` that have a matching `field_name`.
  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.

  @param[in,out] op          `CeedOperator`
  @param[in]     field_label Label of field to set
  @param[in]     field_type  Type of field to set
  @param[out]    num_values  Number of values of type `field_type` in array `values`
  @param[out]    values      Values in the label

  @return An error code: 0 - success, otherwise - failure

  @ref Developer
**/
static int CeedOperatorContextGetGenericRead(CeedOperator op, CeedContextFieldLabel field_label, CeedContextFieldType field_type, size_t *num_values,
                                             void *values) {
  bool is_composite = false;

  CeedCheck(field_label, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, "Invalid field label");

  *(void **)values = NULL;
  *num_values      = 0;

  // Check if field_label and op correspond
  if (field_label->from_op) {
    CeedInt index = -1;

    for (CeedInt i = 0; i < op->num_context_labels; i++) {
      if (op->context_labels[i] == field_label) index = i;
    }
    CeedCheck(index != -1, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, "ContextFieldLabel does not correspond to the operator");
  }

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  if (is_composite) {
    CeedInt       num_sub;
    CeedOperator *sub_operators;

    CeedCall(CeedOperatorCompositeGetNumSub(op, &num_sub));
    CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators));
    CeedCheck(num_sub == field_label->num_sub_labels, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED,
              "Composite operator modified after ContextFieldLabel created");

    for (CeedInt i = 0; i < num_sub; i++) {
      CeedQFunctionContext ctx;

      CeedCall(CeedOperatorGetContext(sub_operators[i], &ctx));
      // Try every sub-operator, ok if some sub-operators do not have field
      if (ctx && field_label->sub_labels[i]) {
        CeedCall(CeedQFunctionContextGetGenericRead(ctx, field_label->sub_labels[i], field_type, num_values, values));
        CeedCall(CeedQFunctionContextDestroy(&ctx));
        return CEED_ERROR_SUCCESS;
      }
      CeedCall(CeedQFunctionContextDestroy(&ctx));
    }
  } else {
    CeedQFunctionContext ctx;

    CeedCall(CeedOperatorGetContext(op, &ctx));
    CeedCheck(ctx, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, "QFunction does not have context data");
    CeedCall(CeedQFunctionContextGetGenericRead(ctx, field_label, field_type, num_values, values));
    CeedCall(CeedQFunctionContextDestroy(&ctx));
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Restore `CeedQFunctionContext` field values of the specified type, read-only.

  For composite operators, the values restored are for the first sub-operator `CeedQFunctionContext` that have a matching `field_name`.
  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.

  @param[in,out] op          `CeedOperator`
  @param[in]     field_label Label of field to set
  @param[in]     field_type  Type of field to set
  @param[in]     values      Values array to restore

  @return An error code: 0 - success, otherwise - failure

  @ref Developer
**/
static int CeedOperatorContextRestoreGenericRead(CeedOperator op, CeedContextFieldLabel field_label, CeedContextFieldType field_type, void *values) {
  bool is_composite = false;

  CeedCheck(field_label, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, "Invalid field label");

  // Check if field_label and op correspond
  if (field_label->from_op) {
    CeedInt index = -1;

    for (CeedInt i = 0; i < op->num_context_labels; i++) {
      if (op->context_labels[i] == field_label) index = i;
    }
    CeedCheck(index != -1, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, "ContextFieldLabel does not correspond to the operator");
  }

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  if (is_composite) {
    CeedInt       num_sub;
    CeedOperator *sub_operators;

    CeedCall(CeedOperatorCompositeGetNumSub(op, &num_sub));
    CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators));
    CeedCheck(num_sub == field_label->num_sub_labels, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED,
              "Composite operator modified after ContextFieldLabel created");

    for (CeedInt i = 0; i < num_sub; i++) {
      CeedQFunctionContext ctx;

      CeedCall(CeedOperatorGetContext(sub_operators[i], &ctx));
      // Try every sub-operator, ok if some sub-operators do not have field
      if (ctx && field_label->sub_labels[i]) {
        CeedCall(CeedQFunctionContextRestoreGenericRead(ctx, field_label->sub_labels[i], field_type, values));
        CeedCall(CeedQFunctionContextDestroy(&ctx));
        return CEED_ERROR_SUCCESS;
      }
      CeedCall(CeedQFunctionContextDestroy(&ctx));
    }
  } else {
    CeedQFunctionContext ctx;

    CeedCall(CeedOperatorGetContext(op, &ctx));
    CeedCheck(ctx, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, "QFunction does not have context data");
    CeedCall(CeedQFunctionContextRestoreGenericRead(ctx, field_label, field_type, values));
    CeedCall(CeedQFunctionContextDestroy(&ctx));
  }
  return CEED_ERROR_SUCCESS;
}

/// @}

/// ----------------------------------------------------------------------------
/// CeedOperator Backend API
/// ----------------------------------------------------------------------------
/// @addtogroup CeedOperatorBackend
/// @{

/**
  @brief Get the number of arguments associated with a `CeedOperator`

  @param[in]  op        `CeedOperator`
  @param[out] num_args  Variable to store vector number of arguments

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedOperatorGetNumArgs(CeedOperator op, CeedInt *num_args) {
  bool is_composite;

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Not defined for composite operators");
  *num_args = op->num_fields;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the tensor product status of all bases for a `CeedOperator`.

  `has_tensor_bases` is only set to `true` if every field uses a tensor-product basis.

  @param[in]  op               `CeedOperator`
  @param[out] has_tensor_bases Variable to store tensor bases status

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedOperatorHasTensorBases(CeedOperator op, bool *has_tensor_bases) {
  CeedInt            num_inputs, num_outputs;
  CeedOperatorField *input_fields, *output_fields;

  CeedCall(CeedOperatorGetFields(op, &num_inputs, &input_fields, &num_outputs, &output_fields));
  *has_tensor_bases = true;
  for (CeedInt i = 0; i < num_inputs; i++) {
    bool      is_tensor;
    CeedBasis basis;

    CeedCall(CeedOperatorFieldGetBasis(input_fields[i], &basis));
    if (basis != CEED_BASIS_NONE) {
      CeedCall(CeedBasisIsTensor(basis, &is_tensor));
      *has_tensor_bases = *has_tensor_bases & is_tensor;
    }
    CeedCall(CeedBasisDestroy(&basis));
  }
  for (CeedInt i = 0; i < num_outputs; i++) {
    bool      is_tensor;
    CeedBasis basis;

    CeedCall(CeedOperatorFieldGetBasis(output_fields[i], &basis));
    if (basis != CEED_BASIS_NONE) {
      CeedCall(CeedBasisIsTensor(basis, &is_tensor));
      *has_tensor_bases = *has_tensor_bases & is_tensor;
    }
    CeedCall(CeedBasisDestroy(&basis));
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get a boolean value indicating if the `CeedOperator` is immutable

  @param[in]  op           `CeedOperator`
  @param[out] is_immutable Variable to store immutability status

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedOperatorIsImmutable(CeedOperator op, bool *is_immutable) {
  *is_immutable = op->is_immutable;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the setup status of a `CeedOperator`

  @param[in]  op            `CeedOperator`
  @param[out] is_setup_done Variable to store setup status

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedOperatorIsSetupDone(CeedOperator op, bool *is_setup_done) {
  *is_setup_done = op->is_backend_setup;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the `CeedQFunction` associated with a `CeedOperator`

  @param[in]  op `CeedOperator`
  @param[out] qf Variable to store `CeedQFunction`

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedOperatorGetQFunction(CeedOperator op, CeedQFunction *qf) {
  bool is_composite;

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Not defined for composite operator");
  *qf = NULL;
  CeedCall(CeedQFunctionReferenceCopy(op->qf, qf));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get a boolean value indicating if the `CeedOperator` is composite

  @param[in]  op           `CeedOperator`
  @param[out] is_composite Variable to store composite status

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedOperatorIsComposite(CeedOperator op, bool *is_composite) {
  *is_composite = op->is_composite;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the backend data of a `CeedOperator`

  @param[in]  op   `CeedOperator`
  @param[out] data Variable to store data

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedOperatorGetData(CeedOperator op, void *data) {
  *(void **)data = op->data;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Set the backend data of a `CeedOperator`

  @param[in,out] op   `CeedOperator`
  @param[in]     data Data to set

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedOperatorSetData(CeedOperator op, void *data) {
  op->data = data;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Increment the reference counter for a `CeedOperator`

  @param[in,out] op `CeedOperator` to increment the reference counter

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedOperatorReference(CeedOperator op) {
  op->ref_count++;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Set the setup flag of a `CeedOperator` to `true`

  @param[in,out] op `CeedOperator`

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedOperatorSetSetupDone(CeedOperator op) {
  op->is_backend_setup = true;
  return CEED_ERROR_SUCCESS;
}

/// @}

/// ----------------------------------------------------------------------------
/// CeedOperator Public API
/// ----------------------------------------------------------------------------
/// @addtogroup CeedOperatorUser
/// @{

/**
  @brief Create a `CeedOperator` and associate a `CeedQFunction`.

  A `CeedBasis` and `CeedElemRestriction` can be associated with `CeedQFunction` fields with @ref CeedOperatorSetField().

  @param[in]  ceed `Ceed` object used to create the `CeedOperator`
  @param[in]  qf   `CeedQFunction` defining the action of the operator at quadrature points
  @param[in]  dqf  `CeedQFunction` defining the action of the Jacobian of `qf` (or @ref CEED_QFUNCTION_NONE)
  @param[in]  dqfT `CeedQFunction` defining the action of the transpose of the Jacobian of `qf` (or @ref CEED_QFUNCTION_NONE)
  @param[out] op   Address of the variable where the newly created `CeedOperator` will be stored

  @return An error code: 0 - success, otherwise - failure

  @ref User
 */
int CeedOperatorCreate(Ceed ceed, CeedQFunction qf, CeedQFunction dqf, CeedQFunction dqfT, CeedOperator *op) {
  if (!ceed->OperatorCreate) {
    Ceed delegate;

    CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Operator"));
    CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedOperatorCreate");
    CeedCall(CeedOperatorCreate(delegate, qf, dqf, dqfT, op));
    CeedCall(CeedDestroy(&delegate));
    return CEED_ERROR_SUCCESS;
  }

  CeedCheck(qf && qf != CEED_QFUNCTION_NONE, ceed, CEED_ERROR_MINOR, "Operator must have a valid CeedQFunction.");

  CeedCall(CeedCalloc(1, op));
  CeedCall(CeedReferenceCopy(ceed, &(*op)->ceed));
  (*op)->ref_count   = 1;
  (*op)->input_size  = -1;
  (*op)->output_size = -1;
  CeedCall(CeedQFunctionReferenceCopy(qf, &(*op)->qf));
  if (dqf && dqf != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionReferenceCopy(dqf, &(*op)->dqf));
  if (dqfT && dqfT != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionReferenceCopy(dqfT, &(*op)->dqfT));
  CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*op)->input_fields));
  CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*op)->output_fields));
  CeedCall(ceed->OperatorCreate(*op));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Create a `CeedOperator` for evaluation at evaluation at arbitrary points in each element.

  A `CeedBasis` and `CeedElemRestriction` can be associated with `CeedQFunction` fields with `CeedOperator` SetField.
  The locations of each point are set with @ref CeedOperatorAtPointsSetPoints().

  @param[in]  ceed `Ceed` object used to create the `CeedOperator`
  @param[in]  qf   `CeedQFunction` defining the action of the operator at quadrature points
  @param[in]  dqf  `CeedQFunction` defining the action of the Jacobian of @a qf (or @ref CEED_QFUNCTION_NONE)
  @param[in]  dqfT `CeedQFunction` defining the action of the transpose of the Jacobian of @a qf (or @ref CEED_QFUNCTION_NONE)
  @param[out] op   Address of the variable where the newly created CeedOperator will be stored

  @return An error code: 0 - success, otherwise - failure

  @ref User
 */
int CeedOperatorCreateAtPoints(Ceed ceed, CeedQFunction qf, CeedQFunction dqf, CeedQFunction dqfT, CeedOperator *op) {
  if (!ceed->OperatorCreateAtPoints) {
    Ceed delegate;

    CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Operator"));
    CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedOperatorCreateAtPoints");
    CeedCall(CeedOperatorCreateAtPoints(delegate, qf, dqf, dqfT, op));
    CeedCall(CeedDestroy(&delegate));
    return CEED_ERROR_SUCCESS;
  }

  CeedCheck(qf && qf != CEED_QFUNCTION_NONE, ceed, CEED_ERROR_MINOR, "Operator must have a valid CeedQFunction.");

  CeedCall(CeedCalloc(1, op));
  CeedCall(CeedReferenceCopy(ceed, &(*op)->ceed));
  (*op)->ref_count    = 1;
  (*op)->is_at_points = true;
  (*op)->input_size   = -1;
  (*op)->output_size  = -1;
  CeedCall(CeedQFunctionReferenceCopy(qf, &(*op)->qf));
  if (dqf && dqf != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionReferenceCopy(dqf, &(*op)->dqf));
  if (dqfT && dqfT != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionReferenceCopy(dqfT, &(*op)->dqfT));
  CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*op)->input_fields));
  CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*op)->output_fields));
  CeedCall(ceed->OperatorCreateAtPoints(*op));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Create a composite `CeedOperator` that composes the action of several `CeedOperator`

  @param[in]  ceed `Ceed` object used to create the `CeedOperator`
  @param[out] op   Address of the variable where the newly created composite `CeedOperator` will be stored

  @return An error code: 0 - success, otherwise - failure

  @ref User
 */
int CeedOperatorCreateComposite(Ceed ceed, CeedOperator *op) {
  if (!ceed->CompositeOperatorCreate) {
    Ceed delegate;

    CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Operator"));
    if (delegate) {
      CeedCall(CeedOperatorCreateComposite(delegate, op));
      CeedCall(CeedDestroy(&delegate));
      return CEED_ERROR_SUCCESS;
    }
  }

  CeedCall(CeedCalloc(1, op));
  CeedCall(CeedReferenceCopy(ceed, &(*op)->ceed));
  (*op)->ref_count    = 1;
  (*op)->is_composite = true;
  CeedCall(CeedCalloc(CEED_COMPOSITE_MAX, &(*op)->sub_operators));
  (*op)->input_size  = -1;
  (*op)->output_size = -1;

  if (ceed->CompositeOperatorCreate) CeedCall(ceed->CompositeOperatorCreate(*op));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Copy the pointer to a `CeedOperator`.

  Both pointers should be destroyed with @ref CeedOperatorDestroy().

  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`.
        This `CeedOperator` will be destroyed if `*op_copy` is the only reference to this `CeedOperator`.

  @param[in]     op      `CeedOperator` to copy reference to
  @param[in,out] op_copy Variable to store copied reference

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorReferenceCopy(CeedOperator op, CeedOperator *op_copy) {
  CeedCall(CeedOperatorReference(op));
  CeedCall(CeedOperatorDestroy(op_copy));
  *op_copy = op;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Provide a field to a `CeedOperator` for use by its `CeedQFunction`.

  This function is used to specify both active and passive fields to a `CeedOperator`.
  For passive fields, a `CeedVector` `vec` must be provided.
  Passive fields can inputs or outputs (updated in-place when operator is applied).

  Active fields must be specified using this function, but their data (in a `CeedVector`) is passed in @ref CeedOperatorApply().
  There can be at most one active input `CeedVector` and at most one active output@ref  CeedVector passed to @ref CeedOperatorApply().

  The number of quadrature points must agree across all points.
  When using @ref CEED_BASIS_NONE, the number of quadrature points is determined by the element size of `rstr`.

  @param[in,out] op         `CeedOperator` on which to provide the field
  @param[in]     field_name Name of the field (to be matched with the name used by `CeedQFunction`)
  @param[in]     rstr       `CeedElemRestriction`
  @param[in]     basis      `CeedBasis` in which the field resides or @ref CEED_BASIS_NONE if collocated with quadrature points
  @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`

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorSetField(CeedOperator op, const char *field_name, CeedElemRestriction rstr, CeedBasis basis, CeedVector vec) {
  bool               is_input = true, is_at_points, is_composite, is_immutable;
  CeedInt            num_elem = 0, num_qpts = 0, num_input_fields, num_output_fields;
  CeedQFunction      qf;
  CeedQFunctionField qf_field, *qf_input_fields, *qf_output_fields;
  CeedOperatorField *op_field;

  CeedCall(CeedOperatorIsAtPoints(op, &is_at_points));
  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCall(CeedOperatorIsImmutable(op, &is_immutable));
  CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE, "Cannot add field to composite operator.");
  CeedCheck(!is_immutable, CeedOperatorReturnCeed(op), CEED_ERROR_MAJOR, "Operator cannot be changed after set as immutable");
  CeedCheck(rstr, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE, "CeedElemRestriction rstr for field \"%s\" must be non-NULL.", field_name);
  CeedCheck(basis, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE, "CeedBasis basis for field \"%s\" must be non-NULL.", field_name);
  CeedCheck(vec, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE, "CeedVector vec for field \"%s\" must be non-NULL.", field_name);

  CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem));
  CeedCheck(rstr == CEED_ELEMRESTRICTION_NONE || !op->has_restriction || num_elem == op->num_elem, CeedOperatorReturnCeed(op), CEED_ERROR_DIMENSION,
            "CeedElemRestriction with %" CeedInt_FMT " elements incompatible with prior %" CeedInt_FMT " elements", num_elem, op->num_elem);
  {
    CeedRestrictionType rstr_type;

    CeedCall(CeedElemRestrictionGetType(rstr, &rstr_type));
    if (rstr_type == CEED_RESTRICTION_POINTS) {
      CeedCheck(is_at_points, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED,
                "CeedElemRestriction AtPoints not supported for standard operator fields");
      CeedCheck(basis == CEED_BASIS_NONE, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED,
                "CeedElemRestriction AtPoints must be used with CEED_BASIS_NONE");
      if (!op->first_points_rstr) {
        CeedCall(CeedElemRestrictionReferenceCopy(rstr, &op->first_points_rstr));
      } else {
        bool are_compatible;

        CeedCall(CeedElemRestrictionAtPointsAreCompatible(op->first_points_rstr, rstr, &are_compatible));
        CeedCheck(are_compatible, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE,
                  "CeedElemRestriction must have compatible offsets with previously set CeedElemRestriction");
      }
    }
  }

  if (basis == CEED_BASIS_NONE) CeedCall(CeedElemRestrictionGetElementSize(rstr, &num_qpts));
  else CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts));
  CeedCheck(op->num_qpts == 0 || num_qpts == op->num_qpts, CeedOperatorReturnCeed(op), CEED_ERROR_DIMENSION,
            "%s must correspond to the same number of quadrature points as previously added CeedBases. Found %" CeedInt_FMT
            " quadrature points but expected %" CeedInt_FMT " quadrature points.",
            basis == CEED_BASIS_NONE ? "CeedElemRestriction" : "CeedBasis", num_qpts, op->num_qpts);

  CeedCall(CeedOperatorGetQFunction(op, &qf));
  CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &qf_input_fields, &num_output_fields, &qf_output_fields));
  CeedCall(CeedQFunctionDestroy(&qf));
  for (CeedInt i = 0; i < num_input_fields; i++) {
    const char *qf_field_name;

    CeedCall(CeedQFunctionFieldGetName(qf_input_fields[i], &qf_field_name));
    if (!strcmp(field_name, qf_field_name)) {
      qf_field = qf_input_fields[i];
      op_field = &op->input_fields[i];
      goto found;
    }
  }
  is_input = false;
  for (CeedInt i = 0; i < num_output_fields; i++) {
    const char *qf_field_name;

    CeedCall(CeedQFunctionFieldGetName(qf_output_fields[i], &qf_field_name));
    if (!strcmp(field_name, qf_field_name)) {
      qf_field = qf_output_fields[i];
      op_field = &op->output_fields[i];
      goto found;
    }
  }
  // LCOV_EXCL_START
  return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE, "CeedQFunction has no knowledge of field '%s'", field_name);
  // LCOV_EXCL_STOP
found:
  CeedCall(CeedOperatorCheckField(CeedOperatorReturnCeed(op), qf_field, rstr, basis));
  CeedCall(CeedCalloc(1, op_field));

  if (vec == CEED_VECTOR_ACTIVE) {
    CeedSize l_size;

    CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &l_size));
    if (is_input) {
      if (op->input_size == -1) op->input_size = l_size;
      CeedCheck(l_size == op->input_size, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE,
                "LVector size %" CeedSize_FMT " does not match previous size %" CeedSize_FMT "", l_size, op->input_size);
    } else {
      if (op->output_size == -1) op->output_size = l_size;
      CeedCheck(l_size == op->output_size, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE,
                "LVector size %" CeedSize_FMT " does not match previous size %" CeedSize_FMT "", l_size, op->output_size);
    }
  }

  CeedCall(CeedVectorReferenceCopy(vec, &(*op_field)->vec));
  CeedCall(CeedElemRestrictionReferenceCopy(rstr, &(*op_field)->elem_rstr));
  if (rstr != CEED_ELEMRESTRICTION_NONE && !op->has_restriction) {
    op->num_elem        = num_elem;
    op->has_restriction = true;  // Restriction set, but num_elem may be 0
  }
  CeedCall(CeedBasisReferenceCopy(basis, &(*op_field)->basis));
  if (op->num_qpts == 0 && !is_at_points) op->num_qpts = num_qpts;  // no consistent number of qpts for OperatorAtPoints
  op->num_fields += 1;
  CeedCall(CeedStringAllocCopy(field_name, (char **)&(*op_field)->field_name));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the `CeedOperator` Field of a `CeedOperator`.

  Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable.

  @param[in]  op                `CeedOperator`
  @param[out] num_input_fields  Variable to store number of input fields
  @param[out] input_fields      Variable to store input fields
  @param[out] num_output_fields Variable to store number of output fields
  @param[out] output_fields     Variable to store output fields

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedOperatorGetFields(CeedOperator op, CeedInt *num_input_fields, CeedOperatorField **input_fields, CeedInt *num_output_fields,
                          CeedOperatorField **output_fields) {
  bool          is_composite;
  CeedQFunction qf;

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Not defined for composite operator");
  CeedCall(CeedOperatorCheckReady(op));

  CeedCall(CeedOperatorGetQFunction(op, &qf));
  CeedCall(CeedQFunctionGetFields(qf, num_input_fields, NULL, num_output_fields, NULL));
  CeedCall(CeedQFunctionDestroy(&qf));
  if (input_fields) *input_fields = op->input_fields;
  if (output_fields) *output_fields = op->output_fields;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Set the arbitrary points in each element for a `CeedOperator` at points.

  Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable.

  @param[in,out] op           `CeedOperator` at points
  @param[in]     rstr_points  `CeedElemRestriction` for the coordinates of each point by element
  @param[in]     point_coords `CeedVector` holding coordinates of each point

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedOperatorAtPointsSetPoints(CeedOperator op, CeedElemRestriction rstr_points, CeedVector point_coords) {
  bool is_at_points, is_immutable;

  CeedCall(CeedOperatorIsAtPoints(op, &is_at_points));
  CeedCall(CeedOperatorIsImmutable(op, &is_immutable));
  CeedCheck(is_at_points, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Only defined for operator at points");
  CeedCheck(!is_immutable, CeedOperatorReturnCeed(op), CEED_ERROR_MAJOR, "Operator cannot be changed after set as immutable");

  if (!op->first_points_rstr) {
    CeedCall(CeedElemRestrictionReferenceCopy(rstr_points, &op->first_points_rstr));
  } else {
    bool are_compatible;

    CeedCall(CeedElemRestrictionAtPointsAreCompatible(op->first_points_rstr, rstr_points, &are_compatible));
    CeedCheck(are_compatible, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE,
              "CeedElemRestriction must have compatible offsets with previously set field CeedElemRestriction");
  }

  CeedCall(CeedElemRestrictionReferenceCopy(rstr_points, &op->rstr_points));
  CeedCall(CeedVectorReferenceCopy(point_coords, &op->point_coords));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get a boolean value indicating if the `CeedOperator` was created with `CeedOperatorCreateAtPoints`

  @param[in]  op           `CeedOperator`
  @param[out] is_at_points Variable to store at points status

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorIsAtPoints(CeedOperator op, bool *is_at_points) {
  *is_at_points = op->is_at_points;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the arbitrary points in each element for a `CeedOperator` at points.

  Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable.

  @param[in]  op           `CeedOperator` at points
  @param[out] rstr_points  Variable to hold `CeedElemRestriction` for the coordinates of each point by element
  @param[out] point_coords Variable to hold `CeedVector` holding coordinates of each point

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedOperatorAtPointsGetPoints(CeedOperator op, CeedElemRestriction *rstr_points, CeedVector *point_coords) {
  bool is_at_points;

  CeedCall(CeedOperatorIsAtPoints(op, &is_at_points));
  CeedCheck(is_at_points, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Only defined for operator at points");
  CeedCall(CeedOperatorCheckReady(op));

  if (rstr_points) {
    *rstr_points = NULL;
    CeedCall(CeedElemRestrictionReferenceCopy(op->rstr_points, rstr_points));
  }
  if (point_coords) {
    *point_coords = NULL;
    CeedCall(CeedVectorReferenceCopy(op->point_coords, point_coords));
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get a `CeedOperator` Field of a `CeedOperator` from its name.

  `op_field` is set to `NULL` if the field is not found.

  Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable.

  @param[in]  op         `CeedOperator`
  @param[in]  field_name Name of desired `CeedOperator` Field
  @param[out] op_field   `CeedOperator` Field corresponding to the name

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedOperatorGetFieldByName(CeedOperator op, const char *field_name, CeedOperatorField *op_field) {
  const char        *name;
  CeedInt            num_input_fields, num_output_fields;
  CeedOperatorField *input_fields, *output_fields;

  *op_field = NULL;
  CeedCall(CeedOperatorGetFields(op, &num_input_fields, &input_fields, &num_output_fields, &output_fields));
  for (CeedInt i = 0; i < num_input_fields; i++) {
    CeedCall(CeedOperatorFieldGetName(input_fields[i], &name));
    if (!strcmp(name, field_name)) {
      *op_field = input_fields[i];
      return CEED_ERROR_SUCCESS;
    }
  }
  for (CeedInt i = 0; i < num_output_fields; i++) {
    CeedCall(CeedOperatorFieldGetName(output_fields[i], &name));
    if (!strcmp(name, field_name)) {
      *op_field = output_fields[i];
      return CEED_ERROR_SUCCESS;
    }
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the name of a `CeedOperator` Field

  @param[in]  op_field   `CeedOperator` Field
  @param[out] field_name Variable to store the field name

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedOperatorFieldGetName(CeedOperatorField op_field, const char **field_name) {
  *field_name = op_field->field_name;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the `CeedElemRestriction` of a `CeedOperator` Field.

  Note: Caller is responsible for destroying the `rstr` with @ref CeedElemRestrictionDestroy().

  @param[in]  op_field `CeedOperator` Field
  @param[out] rstr     Variable to store `CeedElemRestriction`

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedOperatorFieldGetElemRestriction(CeedOperatorField op_field, CeedElemRestriction *rstr) {
  *rstr = NULL;
  CeedCall(CeedElemRestrictionReferenceCopy(op_field->elem_rstr, rstr));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the `CeedBasis` of a `CeedOperator` Field.

  Note: Caller is responsible for destroying the `basis` with @ref CeedBasisDestroy().

  @param[in]  op_field `CeedOperator` Field
  @param[out] basis    Variable to store `CeedBasis`

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedOperatorFieldGetBasis(CeedOperatorField op_field, CeedBasis *basis) {
  *basis = NULL;
  CeedCall(CeedBasisReferenceCopy(op_field->basis, basis));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the `CeedVector` of a `CeedOperator` Field.

  Note: Caller is responsible for destroying the `vec` with @ref CeedVectorDestroy().

  @param[in]  op_field `CeedOperator` Field
  @param[out] vec      Variable to store `CeedVector`

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedOperatorFieldGetVector(CeedOperatorField op_field, CeedVector *vec) {
  *vec = NULL;
  CeedCall(CeedVectorReferenceCopy(op_field->vec, vec));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the data of a `CeedOperator` Field.

  Any arguments set as `NULL` are ignored..

  Note: Caller is responsible for destroying the `rstr`, `basis`, and `vec`.

  @param[in]  op_field   `CeedOperator` Field
  @param[out] field_name Variable to store the field name
  @param[out] rstr       Variable to store `CeedElemRestriction`
  @param[out] basis      Variable to store `CeedBasis`
  @param[out] vec        Variable to store `CeedVector`

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedOperatorFieldGetData(CeedOperatorField op_field, const char **field_name, CeedElemRestriction *rstr, CeedBasis *basis, CeedVector *vec) {
  if (field_name) CeedCall(CeedOperatorFieldGetName(op_field, field_name));
  if (rstr) CeedCall(CeedOperatorFieldGetElemRestriction(op_field, rstr));
  if (basis) CeedCall(CeedOperatorFieldGetBasis(op_field, basis));
  if (vec) CeedCall(CeedOperatorFieldGetVector(op_field, vec));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Add a sub-operator to a composite `CeedOperator`

  @param[in,out] composite_op Composite `CeedOperator`
  @param[in]     sub_op       Sub-operator `CeedOperator`

  @return An error code: 0 - success, otherwise - failure

  @ref User
 */
int CeedOperatorCompositeAddSub(CeedOperator composite_op, CeedOperator sub_op) {
  bool is_immutable;

  CeedCheck(composite_op->is_composite, CeedOperatorReturnCeed(composite_op), CEED_ERROR_MINOR, "CeedOperator is not a composite operator");
  CeedCheck(composite_op->num_suboperators < CEED_COMPOSITE_MAX, CeedOperatorReturnCeed(composite_op), CEED_ERROR_UNSUPPORTED,
            "Cannot add additional sub-operators");
  CeedCall(CeedOperatorIsImmutable(composite_op, &is_immutable));
  CeedCheck(!is_immutable, CeedOperatorReturnCeed(composite_op), CEED_ERROR_MAJOR, "Operator cannot be changed after set as immutable");

  {
    CeedSize input_size, output_size;

    CeedCall(CeedOperatorGetActiveVectorLengths(sub_op, &input_size, &output_size));
    if (composite_op->input_size == -1) composite_op->input_size = input_size;
    if (composite_op->output_size == -1) composite_op->output_size = output_size;
    // Note, a size of -1 means no active vector restriction set, so no incompatibility
    CeedCheck((input_size == -1 || input_size == composite_op->input_size) && (output_size == -1 || output_size == composite_op->output_size),
              CeedOperatorReturnCeed(composite_op), CEED_ERROR_MAJOR,
              "Sub-operators must have compatible dimensions; composite operator of shape (%" CeedSize_FMT ", %" CeedSize_FMT
              ") not compatible with sub-operator of "
              "shape (%" CeedSize_FMT ", %" CeedSize_FMT ")",
              composite_op->input_size, composite_op->output_size, input_size, output_size);
  }

  composite_op->sub_operators[composite_op->num_suboperators] = sub_op;
  CeedCall(CeedOperatorReference(sub_op));
  composite_op->num_suboperators++;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the number of sub-operators associated with a `CeedOperator`

  @param[in]  op               `CeedOperator`
  @param[out] num_suboperators Variable to store number of sub-operators

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedOperatorCompositeGetNumSub(CeedOperator op, CeedInt *num_suboperators) {
  bool is_composite;

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCheck(is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Only defined for a composite operator");
  *num_suboperators = op->num_suboperators;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the list of sub-operators associated with a `CeedOperator`

  @param[in]  op             `CeedOperator`
  @param[out] sub_operators  Variable to store list of sub-operators

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedOperatorCompositeGetSubList(CeedOperator op, CeedOperator **sub_operators) {
  bool is_composite;

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCheck(is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Only defined for a composite operator");
  *sub_operators = op->sub_operators;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get a sub `CeedOperator` of a composite `CeedOperator` from its name.

  `sub_op` is set to `NULL` if the sub operator is not found.

  Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable.

  @param[in]  op      Composite `CeedOperator`
  @param[in]  op_name Name of desired sub `CeedOperator`
  @param[out] sub_op  Sub `CeedOperator` corresponding to the name

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedOperatorCompositeGetSubByName(CeedOperator op, const char *op_name, CeedOperator *sub_op) {
  bool          is_composite;
  CeedInt       num_sub_ops;
  CeedOperator *sub_ops;

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCheck(is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Only defined for a composite operator");
  *sub_op = NULL;
  CeedCall(CeedOperatorCompositeGetNumSub(op, &num_sub_ops));
  CeedCall(CeedOperatorCompositeGetSubList(op, &sub_ops));
  for (CeedInt i = 0; i < num_sub_ops; i++) {
    if (sub_ops[i]->name && !strcmp(op_name, sub_ops[i]->name)) {
      *sub_op = sub_ops[i];
      return CEED_ERROR_SUCCESS;
    }
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Check if a `CeedOperator` is ready to be used.

  @param[in] op `CeedOperator` to check

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorCheckReady(CeedOperator op) {
  bool          is_at_points, is_composite;
  CeedQFunction qf = NULL;

  if (op->is_interface_setup) return CEED_ERROR_SUCCESS;

  CeedCall(CeedOperatorIsAtPoints(op, &is_at_points));
  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  if (!is_composite) CeedCall(CeedOperatorGetQFunction(op, &qf));
  if (is_composite) {
    CeedInt num_suboperators;

    CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators));
    if (!num_suboperators) {
      // Empty operator setup
      op->input_size  = 0;
      op->output_size = 0;
    } else {
      CeedOperator *sub_operators;

      CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators));
      for (CeedInt i = 0; i < num_suboperators; i++) {
        CeedCall(CeedOperatorCheckReady(sub_operators[i]));
      }
      // Sub-operators could be modified after adding to composite operator
      // Need to verify no lvec incompatibility from any changes
      CeedSize input_size, output_size;
      CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size));
    }
  } else {
    CeedInt num_input_fields, num_output_fields;

    CeedCheck(op->num_fields > 0, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE, "No operator fields set");
    CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, NULL, &num_output_fields, NULL));
    CeedCheck(op->num_fields == num_input_fields + num_output_fields, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE,
              "Not all operator fields set");
    CeedCheck(op->has_restriction, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE, "At least one restriction required");
    CeedCheck(op->num_qpts > 0 || is_at_points, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE,
              "At least one non-collocated CeedBasis is required or the number of quadrature points must be set");
  }

  // Flag as immutable and ready
  op->is_interface_setup = true;
  if (qf && qf != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionSetImmutable(qf));
  CeedCall(CeedQFunctionDestroy(&qf));
  if (op->dqf && op->dqf != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionSetImmutable(op->dqf));
  if (op->dqfT && op->dqfT != CEED_QFUNCTION_NONE) CeedCall(CeedQFunctionSetImmutable(op->dqfT));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get vector lengths for the active input and/or output `CeedVector` of a `CeedOperator`.

  Note: Lengths of `-1` indicate that the CeedOperator does not have an active input and/or output.

  @param[in]  op          `CeedOperator`
  @param[out] input_size  Variable to store active input vector length, or `NULL`
  @param[out] output_size Variable to store active output vector length, or `NULL`

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorGetActiveVectorLengths(CeedOperator op, CeedSize *input_size, CeedSize *output_size) {
  bool is_composite;

  if (input_size) *input_size = op->input_size;
  if (output_size) *output_size = op->output_size;

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  if (is_composite && (op->input_size == -1 || op->output_size == -1)) {
    CeedInt       num_suboperators;
    CeedOperator *sub_operators;

    CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators));
    CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators));
    for (CeedInt i = 0; i < num_suboperators; i++) {
      CeedSize sub_input_size, sub_output_size;

      CeedCall(CeedOperatorGetActiveVectorLengths(sub_operators[i], &sub_input_size, &sub_output_size));
      if (op->input_size == -1) op->input_size = sub_input_size;
      if (op->output_size == -1) op->output_size = sub_output_size;
      // Note, a size of -1 means no active vector restriction set, so no incompatibility
      CeedCheck((sub_input_size == -1 || sub_input_size == op->input_size) && (sub_output_size == -1 || sub_output_size == op->output_size),
                CeedOperatorReturnCeed(op), CEED_ERROR_MAJOR,
                "Sub-operators must have compatible dimensions; composite operator of shape (%" CeedSize_FMT ", %" CeedSize_FMT
                ") not compatible with sub-operator of "
                "shape (%" CeedSize_FMT ", %" CeedSize_FMT ")",
                op->input_size, op->output_size, input_size, output_size);
    }
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Set reuse of `CeedQFunction` data in `CeedOperatorLinearAssemble*()` functions.

  When `reuse_assembly_data = false` (default), the `CeedQFunction` associated with this `CeedOperator` is re-assembled every time a `CeedOperatorLinearAssemble*()` function is called.
  When `reuse_assembly_data = true`, the `CeedQFunction` associated with this `CeedOperator` is reused between calls to @ref CeedOperatorSetQFunctionAssemblyDataUpdateNeeded().

  @param[in] op                  `CeedOperator`
  @param[in] reuse_assembly_data Boolean flag setting assembly data reuse

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedOperatorSetQFunctionAssemblyReuse(CeedOperator op, bool reuse_assembly_data) {
  bool is_composite;

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  if (is_composite) {
    for (CeedInt i = 0; i < op->num_suboperators; i++) {
      CeedCall(CeedOperatorSetQFunctionAssemblyReuse(op->sub_operators[i], reuse_assembly_data));
    }
  } else {
    CeedQFunctionAssemblyData data;

    CeedCall(CeedOperatorGetQFunctionAssemblyData(op, &data));
    CeedCall(CeedQFunctionAssemblyDataSetReuse(data, reuse_assembly_data));
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Mark `CeedQFunction` data as updated and the `CeedQFunction` as requiring re-assembly.

  @param[in] op                `CeedOperator`
  @param[in] needs_data_update Boolean flag setting assembly data reuse

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(CeedOperator op, bool needs_data_update) {
  bool is_composite;

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  if (is_composite) {
    CeedInt       num_suboperators;
    CeedOperator *sub_operators;

    CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators));
    CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators));
    for (CeedInt i = 0; i < num_suboperators; i++) {
      CeedCall(CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(sub_operators[i], needs_data_update));
    }
  } else {
    CeedQFunctionAssemblyData data;

    CeedCall(CeedOperatorGetQFunctionAssemblyData(op, &data));
    CeedCall(CeedQFunctionAssemblyDataSetUpdateNeeded(data, needs_data_update));
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Set name of `CeedOperator` for @ref CeedOperatorView() output

  @param[in,out] op   `CeedOperator`
  @param[in]     name Name to set, or NULL to remove previously set name

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorSetName(CeedOperator op, const char *name) {
  char  *name_copy;
  size_t name_len = name ? strlen(name) : 0;

  CeedCall(CeedFree(&op->name));
  if (name_len > 0) {
    CeedCall(CeedCalloc(name_len + 1, &name_copy));
    memcpy(name_copy, name, name_len);
    op->name = name_copy;
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get name of `CeedOperator`

  @param[in]     op   `CeedOperator`
  @param[in,out] name Address of variable to hold currently set name

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorGetName(CeedOperator op, const char **name) {
  if (op->name) {
    *name = op->name;
  } else if (!op->is_composite) {
    CeedQFunction qf;

    CeedCall(CeedOperatorGetQFunction(op, &qf));
    if (qf) CeedCall(CeedQFunctionGetName(qf, name));
    CeedCall(CeedQFunctionDestroy(&qf));
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Core logic for viewing a `CeedOperator`

  @param[in] op     `CeedOperator` to view brief summary
  @param[in] stream  Stream to write; typically `stdout` or a file
  @param[in] is_full Whether to write full operator view or terse

  @return Error code: 0 - success, otherwise - failure

  @ref Developer
**/
static int CeedOperatorView_Core(CeedOperator op, FILE *stream, bool is_full) {
  bool        has_name, is_composite, is_at_points;
  const char *name = NULL;

  CeedCall(CeedOperatorGetName(op, &name));
  has_name = name ? strlen(name) : false;
  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCall(CeedOperatorIsAtPoints(op, &is_at_points));
  if (is_composite) {
    CeedInt       num_suboperators;
    CeedOperator *sub_operators;

    CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators));
    CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators));
    fprintf(stream, "Composite CeedOperator%s%s\n", has_name ? " - " : "", has_name ? name : "");

    for (CeedInt i = 0; i < num_suboperators; i++) {
      has_name = sub_operators[i]->name;
      fprintf(stream, "  SubOperator%s %" CeedInt_FMT "%s%s%s\n", is_at_points ? " AtPoints" : "", i, has_name ? " - " : "",
              has_name ? sub_operators[i]->name : "", is_full ? ":" : "");
      if (is_full) CeedCall(CeedOperatorSingleView(sub_operators[i], 1, stream));
    }
  } else {
    fprintf(stream, "CeedOperator%s%s%s\n", is_at_points ? " AtPoints" : "", has_name ? " - " : "", has_name ? name : "");
    if (is_full) CeedCall(CeedOperatorSingleView(op, 0, stream));
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief View a `CeedOperator`

  @param[in] op     `CeedOperator` to view
  @param[in] stream Stream to write; typically `stdout` or a file

  @return Error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorView(CeedOperator op, FILE *stream) {
  CeedCall(CeedOperatorView_Core(op, stream, true));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief View a brief summary `CeedOperator`

  @param[in] op     `CeedOperator` to view brief summary
  @param[in] stream Stream to write; typically `stdout` or a file

  @return Error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorViewTerse(CeedOperator op, FILE *stream) {
  CeedCall(CeedOperatorView_Core(op, stream, false));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the `Ceed` associated with a `CeedOperator`

  @param[in]  op   `CeedOperator`
  @param[out] ceed Variable to store `Ceed`

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedOperatorGetCeed(CeedOperator op, Ceed *ceed) {
  *ceed = NULL;
  CeedCall(CeedReferenceCopy(CeedOperatorReturnCeed(op), ceed));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Return the `Ceed` associated with a `CeedOperator`

  @param[in]  op `CeedOperator`

  @return `Ceed` associated with the `op`

  @ref Advanced
**/
Ceed CeedOperatorReturnCeed(CeedOperator op) { return op->ceed; }

/**
  @brief Get the number of elements associated with a `CeedOperator`

  @param[in]  op       `CeedOperator`
  @param[out] num_elem Variable to store number of elements

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedOperatorGetNumElements(CeedOperator op, CeedInt *num_elem) {
  bool is_composite;

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Not defined for composite operator");
  *num_elem = op->num_elem;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the number of quadrature points associated with a `CeedOperator`

  @param[in]  op       `CeedOperator`
  @param[out] num_qpts Variable to store vector number of quadrature points

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedOperatorGetNumQuadraturePoints(CeedOperator op, CeedInt *num_qpts) {
  bool is_composite;

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_MINOR, "Not defined for composite operator");
  *num_qpts = op->num_qpts;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Estimate number of FLOPs required to apply `CeedOperator` on the active `CeedVector`

  @param[in]  op    `CeedOperator` to estimate FLOPs for
  @param[out] flops Address of variable to hold FLOPs estimate

  @ref Backend
**/
int CeedOperatorGetFlopsEstimate(CeedOperator op, CeedSize *flops) {
  bool is_composite;

  CeedCall(CeedOperatorCheckReady(op));

  *flops = 0;
  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  if (is_composite) {
    CeedInt num_suboperators;

    CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators));
    CeedOperator *sub_operators;
    CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators));

    // FLOPs for each suboperator
    for (CeedInt i = 0; i < num_suboperators; i++) {
      CeedSize suboperator_flops;

      CeedCall(CeedOperatorGetFlopsEstimate(sub_operators[i], &suboperator_flops));
      *flops += suboperator_flops;
    }
  } else {
    bool                is_at_points;
    CeedInt             num_input_fields, num_output_fields, num_elem = 0, num_points = 0;
    CeedQFunction       qf;
    CeedQFunctionField *qf_input_fields, *qf_output_fields;
    CeedOperatorField  *op_input_fields, *op_output_fields;

    CeedCall(CeedOperatorGetNumElements(op, &num_elem));
    if (num_elem == 0) return CEED_ERROR_SUCCESS;
    CeedCall(CeedOperatorIsAtPoints(op, &is_at_points));
    if (is_at_points) {
      CeedMemType         mem_type;
      CeedElemRestriction rstr_points = NULL;

      CeedCall(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL));
      CeedCall(CeedGetPreferredMemType(CeedOperatorReturnCeed(op), &mem_type));
      if (mem_type == CEED_MEM_DEVICE) {
        // Device backends pad out to the same number of points per element
        CeedCall(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &num_points));
      } else {
        num_points = 0;
        for (CeedInt i = 0; i < num_elem; i++) {
          CeedInt points_in_elem = 0;

          CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr_points, i, &points_in_elem));
          num_points += points_in_elem;
        }
        num_points = num_points / num_elem + (num_points % num_elem > 0);
      }
      CeedCall(CeedElemRestrictionDestroy(&rstr_points));
    }
    CeedCall(CeedOperatorGetQFunction(op, &qf));
    CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &qf_input_fields, &num_output_fields, &qf_output_fields));
    CeedCall(CeedQFunctionDestroy(&qf));
    CeedCall(CeedOperatorGetFields(op, NULL, &op_input_fields, NULL, &op_output_fields));

    // Input FLOPs
    for (CeedInt i = 0; i < num_input_fields; i++) {
      CeedVector vec;

      CeedCall(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
      if (vec == CEED_VECTOR_ACTIVE) {
        CeedEvalMode        eval_mode;
        CeedSize            rstr_flops, basis_flops;
        CeedElemRestriction rstr;
        CeedBasis           basis;

        CeedCall(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &rstr));
        CeedCall(CeedElemRestrictionGetFlopsEstimate(rstr, CEED_NOTRANSPOSE, &rstr_flops));
        CeedCall(CeedElemRestrictionDestroy(&rstr));
        *flops += rstr_flops;
        CeedCall(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
        CeedCall(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
        CeedCall(CeedBasisGetFlopsEstimate(basis, CEED_NOTRANSPOSE, eval_mode, is_at_points, num_points, &basis_flops));
        CeedCall(CeedBasisDestroy(&basis));
        *flops += basis_flops * num_elem;
      }
      CeedCall(CeedVectorDestroy(&vec));
    }
    // QF FLOPs
    {
      CeedInt       num_qpts;
      CeedSize      qf_flops;
      CeedQFunction qf;

      if (is_at_points) num_qpts = num_points;
      else CeedCall(CeedOperatorGetNumQuadraturePoints(op, &num_qpts));
      CeedCall(CeedOperatorGetQFunction(op, &qf));
      CeedCall(CeedQFunctionGetFlopsEstimate(qf, &qf_flops));
      CeedCall(CeedQFunctionDestroy(&qf));
      CeedCheck(qf_flops > -1, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE,
                "Must set CeedQFunction FLOPs estimate with CeedQFunctionSetUserFlopsEstimate");
      *flops += num_elem * num_qpts * qf_flops;
    }

    // Output FLOPs
    for (CeedInt i = 0; i < num_output_fields; i++) {
      CeedVector vec;

      CeedCall(CeedOperatorFieldGetVector(op_output_fields[i], &vec));
      if (vec == CEED_VECTOR_ACTIVE) {
        CeedEvalMode        eval_mode;
        CeedSize            rstr_flops, basis_flops;
        CeedElemRestriction rstr;
        CeedBasis           basis;

        CeedCall(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &rstr));
        CeedCall(CeedElemRestrictionGetFlopsEstimate(rstr, CEED_TRANSPOSE, &rstr_flops));
        CeedCall(CeedElemRestrictionDestroy(&rstr));
        *flops += rstr_flops;
        CeedCall(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
        CeedCall(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
        CeedCall(CeedBasisGetFlopsEstimate(basis, CEED_TRANSPOSE, eval_mode, is_at_points, num_points, &basis_flops));
        CeedCall(CeedBasisDestroy(&basis));
        *flops += basis_flops * num_elem;
      }
      CeedCall(CeedVectorDestroy(&vec));
    }
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get `CeedQFunction` global context for a `CeedOperator`.

  The caller is responsible for destroying `ctx` returned from this function via @ref CeedQFunctionContextDestroy().

  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`.
        This `CeedQFunctionContext` will be destroyed if `ctx` is the only reference to this `CeedQFunctionContext`.

  @param[in]  op  `CeedOperator`
  @param[out] ctx Variable to store `CeedQFunctionContext`

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedOperatorGetContext(CeedOperator op, CeedQFunctionContext *ctx) {
  bool                 is_composite;
  CeedQFunction        qf;
  CeedQFunctionContext qf_ctx;

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE, "Cannot retrieve CeedQFunctionContext for composite operator");
  CeedCall(CeedOperatorGetQFunction(op, &qf));
  CeedCall(CeedQFunctionGetInnerContext(qf, &qf_ctx));
  CeedCall(CeedQFunctionDestroy(&qf));
  *ctx = NULL;
  if (qf_ctx) CeedCall(CeedQFunctionContextReferenceCopy(qf_ctx, ctx));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get label for a registered `CeedQFunctionContext` field, or `NULL` if no field has been registered with this `field_name`.

  Fields are registered via `CeedQFunctionContextRegister*()` functions (eg. @ref CeedQFunctionContextRegisterDouble()).

  @param[in]  op          `CeedOperator`
  @param[in]  field_name  Name of field to retrieve label
  @param[out] field_label Variable to field label

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorGetContextFieldLabel(CeedOperator op, const char *field_name, CeedContextFieldLabel *field_label) {
  bool is_composite, field_found = false;

  CeedCall(CeedOperatorIsComposite(op, &is_composite));

  if (is_composite) {
    // Composite operator
    // -- Check if composite label already created
    for (CeedInt i = 0; i < op->num_context_labels; i++) {
      if (!strcmp(op->context_labels[i]->name, field_name)) {
        *field_label = op->context_labels[i];
        return CEED_ERROR_SUCCESS;
      }
    }

    // -- Create composite label if needed
    CeedInt               num_sub;
    CeedOperator         *sub_operators;
    CeedContextFieldLabel new_field_label;

    CeedCall(CeedCalloc(1, &new_field_label));
    CeedCall(CeedOperatorCompositeGetNumSub(op, &num_sub));
    CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators));
    CeedCall(CeedCalloc(num_sub, &new_field_label->sub_labels));
    new_field_label->num_sub_labels = num_sub;

    for (CeedInt i = 0; i < num_sub; i++) {
      if (sub_operators[i]->qf->ctx) {
        CeedContextFieldLabel new_field_label_i;

        CeedCall(CeedQFunctionContextGetFieldLabel(sub_operators[i]->qf->ctx, field_name, &new_field_label_i));
        if (new_field_label_i) {
          field_found                    = true;
          new_field_label->sub_labels[i] = new_field_label_i;
          new_field_label->name          = new_field_label_i->name;
          new_field_label->description   = new_field_label_i->description;
          if (new_field_label->type && new_field_label->type != new_field_label_i->type) {
            // LCOV_EXCL_START
            CeedCall(CeedFree(&new_field_label));
            return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE, "Incompatible field types on sub-operator contexts. %s != %s",
                             CeedContextFieldTypes[new_field_label->type], CeedContextFieldTypes[new_field_label_i->type]);
            // LCOV_EXCL_STOP
          } else {
            new_field_label->type = new_field_label_i->type;
          }
          if (new_field_label->num_values != 0 && new_field_label->num_values != new_field_label_i->num_values) {
            // LCOV_EXCL_START
            CeedCall(CeedFree(&new_field_label));
            return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE,
                             "Incompatible field number of values on sub-operator contexts. %zu != %zu", new_field_label->num_values,
                             new_field_label_i->num_values);
            // LCOV_EXCL_STOP
          } else {
            new_field_label->num_values = new_field_label_i->num_values;
          }
        }
      }
    }
    // -- Cleanup if field was found
    if (field_found) {
      *field_label = new_field_label;
    } else {
      // LCOV_EXCL_START
      CeedCall(CeedFree(&new_field_label->sub_labels));
      CeedCall(CeedFree(&new_field_label));
      *field_label = NULL;
      // LCOV_EXCL_STOP
    }
  } else {
    CeedQFunction        qf;
    CeedQFunctionContext ctx;

    // Single, non-composite operator
    CeedCall(CeedOperatorGetQFunction(op, &qf));
    CeedCall(CeedQFunctionGetInnerContext(qf, &ctx));
    CeedCall(CeedQFunctionDestroy(&qf));
    if (ctx) {
      CeedCall(CeedQFunctionContextGetFieldLabel(ctx, field_name, field_label));
    } else {
      *field_label = NULL;
    }
  }

  // Set label in operator
  if (*field_label) {
    (*field_label)->from_op = true;

    // Move new composite label to operator
    if (op->num_context_labels == 0) {
      CeedCall(CeedCalloc(1, &op->context_labels));
      op->max_context_labels = 1;
    } else if (op->num_context_labels == op->max_context_labels) {
      CeedCall(CeedRealloc(2 * op->num_context_labels, &op->context_labels));
      op->max_context_labels *= 2;
    }
    op->context_labels[op->num_context_labels] = *field_label;
    op->num_context_labels++;
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Set `CeedQFunctionContext` field holding double precision values.

  For composite operators, the values are set in all sub-operator `CeedQFunctionContext` that have a matching `field_name`.

  @param[in,out] op          `CeedOperator`
  @param[in]     field_label Label of field to set
  @param[in]     values      Values to set

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorSetContextDouble(CeedOperator op, CeedContextFieldLabel field_label, double *values) {
  return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, values);
}

/**
  @brief Get `CeedQFunctionContext` field holding double precision values, read-only.

  For composite operators, the values correspond to the first sub-operator `CeedQFunctionContext` that has a matching `field_name`.

  @param[in]  op          `CeedOperator`
  @param[in]  field_label Label of field to get
  @param[out] num_values  Number of values in the field label
  @param[out] values      Pointer to context values

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorGetContextDoubleRead(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const double **values) {
  return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, num_values, values);
}

/**
  @brief Restore `CeedQFunctionContext` field holding double precision values, read-only.

  @param[in]  op          `CeedOperator`
  @param[in]  field_label Label of field to restore
  @param[out] values      Pointer to context values

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorRestoreContextDoubleRead(CeedOperator op, CeedContextFieldLabel field_label, const double **values) {
  return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_DOUBLE, values);
}

/**
  @brief Set `CeedQFunctionContext` field holding `int32` values.

  For composite operators, the values are set in all sub-operator `CeedQFunctionContext` that have a matching `field_name`.

  @param[in,out] op          `CeedOperator`
  @param[in]     field_label Label of field to set
  @param[in]     values      Values to set

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorSetContextInt32(CeedOperator op, CeedContextFieldLabel field_label, int32_t *values) {
  return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_INT32, values);
}

/**
  @brief Get `CeedQFunctionContext` field holding `int32` values, read-only.

  For composite operators, the values correspond to the first sub-operator `CeedQFunctionContext` that has a matching `field_name`.

  @param[in]  op          `CeedOperator`
  @param[in]  field_label Label of field to get
  @param[out] num_values  Number of `int32` values in `values`
  @param[out] values      Pointer to context values

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorGetContextInt32Read(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const int32_t **values) {
  return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_INT32, num_values, values);
}

/**
  @brief Restore `CeedQFunctionContext` field holding `int32` values, read-only.

  @param[in]  op          `CeedOperator`
  @param[in]  field_label Label of field to get
  @param[out] values      Pointer to context values

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorRestoreContextInt32Read(CeedOperator op, CeedContextFieldLabel field_label, const int32_t **values) {
  return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_INT32, values);
}

/**
  @brief Set `CeedQFunctionContext` field holding boolean values.

  For composite operators, the values are set in all sub-operator `CeedQFunctionContext` that have a matching `field_name`.

  @param[in,out] op          `CeedOperator`
  @param[in]     field_label Label of field to set
  @param[in]     values      Values to set

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorSetContextBoolean(CeedOperator op, CeedContextFieldLabel field_label, bool *values) {
  return CeedOperatorContextSetGeneric(op, field_label, CEED_CONTEXT_FIELD_BOOL, values);
}

/**
  @brief Get `CeedQFunctionContext` field holding boolean values, read-only.

  For composite operators, the values correspond to the first sub-operator `CeedQFunctionContext` that has a matching `field_name`.

  @param[in]  op          `CeedOperator`
  @param[in]  field_label Label of field to get
  @param[out] num_values  Number of boolean values in `values`
  @param[out] values      Pointer to context values

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorGetContextBooleanRead(CeedOperator op, CeedContextFieldLabel field_label, size_t *num_values, const bool **values) {
  return CeedOperatorContextGetGenericRead(op, field_label, CEED_CONTEXT_FIELD_BOOL, num_values, values);
}

/**
  @brief Restore `CeedQFunctionContext` field holding boolean values, read-only.

  @param[in]  op          `CeedOperator`
  @param[in]  field_label Label of field to get
  @param[out] values      Pointer to context values

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorRestoreContextBooleanRead(CeedOperator op, CeedContextFieldLabel field_label, const bool **values) {
  return CeedOperatorContextRestoreGenericRead(op, field_label, CEED_CONTEXT_FIELD_BOOL, values);
}

/**
  @brief Apply `CeedOperator` to a `CeedVector`.

  This computes the action of the operator on the specified (active) input, yielding its (active) output.
  All inputs and outputs must be specified using @ref CeedOperatorSetField().

  @note Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable.

  @param[in]  op      `CeedOperator` to apply
  @param[in]  in      `CeedVector` containing input state or @ref CEED_VECTOR_NONE if there are no active inputs
  @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
  @param[in]  request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorApply(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) {
  bool is_composite;

  CeedCall(CeedOperatorCheckReady(op));

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  if (is_composite && op->ApplyComposite) {
    // Composite Operator
    CeedCall(op->ApplyComposite(op, in, out, request));
  } else if (!is_composite && op->Apply) {
    // Standard Operator
    CeedCall(op->Apply(op, in, out, request));
  } else {
    // Standard or composite, default to zeroing out and calling ApplyAddActive
    // Zero active output
    if (out != CEED_VECTOR_NONE) CeedCall(CeedVectorSetValue(out, 0.0));

    // ApplyAddActive
    CeedCall(CeedOperatorApplyAddActive(op, in, out, request));
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Apply `CeedOperator` to a `CeedVector` and add result to output `CeedVector`.

  This computes the action of the operator on the specified (active) input, yielding its (active) output.
  All inputs and outputs must be specified using @ref CeedOperatorSetField().

  @note Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable.
  @warning This function adds into ALL outputs, including passive outputs. To only add into the active output, use `CeedOperatorApplyAddActive()`.
  @see `CeedOperatorApplyAddActive()`

  @param[in]  op      `CeedOperator` to apply
  @param[in]  in      `CeedVector` containing input state or @ref CEED_VECTOR_NONE if there are no active inputs
  @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
  @param[in]  request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorApplyAdd(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) {
  bool is_composite;

  CeedCall(CeedOperatorCheckReady(op));

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  if (is_composite) {
    // Composite Operator
    if (op->ApplyAddComposite) {
      CeedCall(op->ApplyAddComposite(op, in, out, request));
    } else {
      CeedInt       num_suboperators;
      CeedOperator *sub_operators;

      CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators));
      CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators));
      for (CeedInt i = 0; i < num_suboperators; i++) {
        CeedCall(CeedOperatorApplyAdd(sub_operators[i], in, out, request));
      }
    }
  } else if (op->num_elem > 0) {
    // Standard Operator
    CeedCall(op->ApplyAdd(op, in, out, request));
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Apply `CeedOperator` to a `CeedVector` and add result to output `CeedVector`. Only sums into active outputs, overwrites passive outputs.

  This computes the action of the operator on the specified (active) input, yielding its (active) output.
  All inputs and outputs must be specified using @ref CeedOperatorSetField().

  @note Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable.

  @param[in]  op      `CeedOperator` to apply
  @param[in]  in      `CeedVector` containing input state or @ref CEED_VECTOR_NONE if there are no active inputs
  @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
  @param[in]  request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorApplyAddActive(CeedOperator op, CeedVector in, CeedVector out, CeedRequest *request) {
  bool is_composite;

  CeedCall(CeedOperatorCheckReady(op));

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  if (is_composite) {
    // Composite Operator
    CeedInt       num_suboperators;
    CeedOperator *sub_operators;

    CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators));
    CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators));

    // Zero all output vectors
    for (CeedInt i = 0; i < num_suboperators; i++) {
      CeedInt            num_output_fields;
      CeedOperatorField *output_fields;

      CeedCall(CeedOperatorGetFields(sub_operators[i], NULL, NULL, &num_output_fields, &output_fields));
      for (CeedInt j = 0; j < num_output_fields; j++) {
        CeedVector vec;

        CeedCall(CeedOperatorFieldGetVector(output_fields[j], &vec));
        if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) CeedCall(CeedVectorSetValue(vec, 0.0));
        CeedCall(CeedVectorDestroy(&vec));
      }
    }
    // ApplyAdd
    CeedCall(CeedOperatorApplyAdd(op, in, out, request));
  } else {
    // Standard Operator
    CeedInt            num_output_fields;
    CeedOperatorField *output_fields;

    CeedCall(CeedOperatorGetFields(op, NULL, NULL, &num_output_fields, &output_fields));
    // Zero all output vectors
    for (CeedInt i = 0; i < num_output_fields; i++) {
      CeedVector vec;

      CeedCall(CeedOperatorFieldGetVector(output_fields[i], &vec));
      if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) CeedCall(CeedVectorSetValue(vec, 0.0));
      CeedCall(CeedVectorDestroy(&vec));
    }
    // ApplyAdd
    CeedCall(CeedOperatorApplyAdd(op, in, out, request));
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Destroy temporary assembly data associated with a `CeedOperator`

  @param[in,out] op `CeedOperator` whose assembly data to destroy

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorAssemblyDataStrip(CeedOperator op) {
  bool is_composite;

  CeedCall(CeedQFunctionAssemblyDataDestroy(&op->qf_assembled));
  CeedCall(CeedOperatorAssemblyDataDestroy(&op->op_assembled));
  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  if (is_composite) {
    CeedInt       num_suboperators;
    CeedOperator *sub_operators;

    CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators));
    CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators));
    for (CeedInt i = 0; i < num_suboperators; i++) {
      CeedCall(CeedQFunctionAssemblyDataDestroy(&sub_operators[i]->qf_assembled));
      CeedCall(CeedOperatorAssemblyDataDestroy(&sub_operators[i]->op_assembled));
    }
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Destroy a `CeedOperator`

  @param[in,out] op `CeedOperator` to destroy

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedOperatorDestroy(CeedOperator *op) {
  if (!*op || --(*op)->ref_count > 0) {
    *op = NULL;
    return CEED_ERROR_SUCCESS;
  }
  // Backend destroy
  if ((*op)->Destroy) {
    CeedCall((*op)->Destroy(*op));
  }
  // Free fields
  for (CeedInt i = 0; i < (*op)->num_fields; i++) {
    if ((*op)->input_fields[i]) {
      if ((*op)->input_fields[i]->elem_rstr != CEED_ELEMRESTRICTION_NONE) {
        CeedCall(CeedElemRestrictionDestroy(&(*op)->input_fields[i]->elem_rstr));
      }
      if ((*op)->input_fields[i]->basis != CEED_BASIS_NONE) {
        CeedCall(CeedBasisDestroy(&(*op)->input_fields[i]->basis));
      }
      if ((*op)->input_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->input_fields[i]->vec != CEED_VECTOR_NONE) {
        CeedCall(CeedVectorDestroy(&(*op)->input_fields[i]->vec));
      }
      CeedCall(CeedFree(&(*op)->input_fields[i]->field_name));
      CeedCall(CeedFree(&(*op)->input_fields[i]));
    }
  }
  for (CeedInt i = 0; i < (*op)->num_fields; i++) {
    if ((*op)->output_fields[i]) {
      CeedCall(CeedElemRestrictionDestroy(&(*op)->output_fields[i]->elem_rstr));
      if ((*op)->output_fields[i]->basis != CEED_BASIS_NONE) {
        CeedCall(CeedBasisDestroy(&(*op)->output_fields[i]->basis));
      }
      if ((*op)->output_fields[i]->vec != CEED_VECTOR_ACTIVE && (*op)->output_fields[i]->vec != CEED_VECTOR_NONE) {
        CeedCall(CeedVectorDestroy(&(*op)->output_fields[i]->vec));
      }
      CeedCall(CeedFree(&(*op)->output_fields[i]->field_name));
      CeedCall(CeedFree(&(*op)->output_fields[i]));
    }
  }
  CeedCall(CeedFree(&(*op)->input_fields));
  CeedCall(CeedFree(&(*op)->output_fields));
  // Destroy AtPoints data
  CeedCall(CeedVectorDestroy(&(*op)->point_coords));
  CeedCall(CeedElemRestrictionDestroy(&(*op)->rstr_points));
  CeedCall(CeedElemRestrictionDestroy(&(*op)->first_points_rstr));
  // Destroy assembly data (must happen before destroying sub_operators)
  CeedCall(CeedOperatorAssemblyDataStrip(*op));
  // Destroy sub_operators
  for (CeedInt i = 0; i < (*op)->num_suboperators; i++) {
    if ((*op)->sub_operators[i]) {
      CeedCall(CeedOperatorDestroy(&(*op)->sub_operators[i]));
    }
  }
  CeedCall(CeedFree(&(*op)->sub_operators));
  CeedCall(CeedQFunctionDestroy(&(*op)->qf));
  CeedCall(CeedQFunctionDestroy(&(*op)->dqf));
  CeedCall(CeedQFunctionDestroy(&(*op)->dqfT));
  // Destroy any composite labels
  if ((*op)->is_composite) {
    for (CeedInt i = 0; i < (*op)->num_context_labels; i++) {
      CeedCall(CeedFree(&(*op)->context_labels[i]->sub_labels));
      CeedCall(CeedFree(&(*op)->context_labels[i]));
    }
  }
  CeedCall(CeedFree(&(*op)->context_labels));

  // Destroy fallback
  CeedCall(CeedOperatorDestroy(&(*op)->op_fallback));

  CeedCall(CeedFree(&(*op)->name));
  CeedCall(CeedDestroy(&(*op)->ceed));
  CeedCall(CeedFree(op));
  return CEED_ERROR_SUCCESS;
}

/// @}