xref: /libCEED/rust/libceed-sys/c-src/interface/ceed-preconditioning.c (revision 681d0ea73ee05192cf73f31e6e4a886b41175395)

// Copyright (c) 2017-2024, 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 <assert.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <string.h>

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

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

/**
  @brief Duplicate a `CeedQFunction` with a reference `Ceed` to fallback for advanced `CeedOperator` functionality

  @param[in]  fallback_ceed `Ceed` on which to create fallback `CeedQFunction`
  @param[in]  qf            `CeedQFunction` to create fallback for
  @param[out] qf_fallback   Fallback `CeedQFunction`

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

  @ref Developer
**/
static int CeedQFunctionCreateFallback(Ceed fallback_ceed, CeedQFunction qf, CeedQFunction *qf_fallback) {
  char               *source_path_with_name = NULL;
  CeedInt             num_input_fields, num_output_fields;
  Ceed                ceed;
  CeedQFunctionField *input_fields, *output_fields;

  // Check if NULL qf passed in
  if (!qf) return CEED_ERROR_SUCCESS;

  CeedCall(CeedQFunctionGetCeed(qf, &ceed));
  CeedDebug256(ceed, 1, "---------- CeedOperator Fallback ----------\n");
  CeedDebug(ceed, "Creating fallback CeedQFunction\n");

  if (qf->source_path) {
    size_t path_len = strlen(qf->source_path), name_len = strlen(qf->kernel_name);
    CeedCall(CeedCalloc(path_len + name_len + 2, &source_path_with_name));
    memcpy(source_path_with_name, qf->source_path, path_len);
    memcpy(&source_path_with_name[path_len], ":", 1);
    memcpy(&source_path_with_name[path_len + 1], qf->kernel_name, name_len);
  } else {
    CeedCall(CeedCalloc(1, &source_path_with_name));
  }

  {
    CeedInt           vec_length;
    CeedQFunctionUser f;

    CeedCall(CeedQFunctionGetVectorLength(qf, &vec_length));
    CeedCall(CeedQFunctionGetUserFunction(qf, &f));
    CeedCall(CeedQFunctionCreateInterior(fallback_ceed, vec_length, f, source_path_with_name, qf_fallback));
  }
  {
    CeedQFunctionContext ctx;

    CeedCall(CeedQFunctionGetContext(qf, &ctx));
    CeedCall(CeedQFunctionSetContext(*qf_fallback, ctx));
  }
  CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &input_fields, &num_output_fields, &output_fields));
  for (CeedInt i = 0; i < num_input_fields; i++) {
    const char  *field_name;
    CeedInt      size;
    CeedEvalMode eval_mode;

    CeedCall(CeedQFunctionFieldGetData(input_fields[i], &field_name, &size, &eval_mode));
    CeedCall(CeedQFunctionAddInput(*qf_fallback, field_name, size, eval_mode));
  }
  for (CeedInt i = 0; i < num_output_fields; i++) {
    const char  *field_name;
    CeedInt      size;
    CeedEvalMode eval_mode;

    CeedCall(CeedQFunctionFieldGetData(output_fields[i], &field_name, &size, &eval_mode));
    CeedCall(CeedQFunctionAddOutput(*qf_fallback, field_name, size, eval_mode));
  }
  CeedCall(CeedFree(&source_path_with_name));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Duplicate a `CeedOperator` with a reference `Ceed` to fallback for advanced `CeedOperator` functionality

  @param[in,out] op `CeedOperator` to create fallback for

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

  @ref Developer
**/
static int CeedOperatorCreateFallback(CeedOperator op) {
  bool         is_composite;
  Ceed         ceed, ceed_fallback;
  CeedOperator op_fallback;

  // Check not already created
  if (op->op_fallback) return CEED_ERROR_SUCCESS;

  // Fallback Ceed
  CeedCall(CeedOperatorGetCeed(op, &ceed));
  CeedCall(CeedGetOperatorFallbackCeed(ceed, &ceed_fallback));
  if (!ceed_fallback) return CEED_ERROR_SUCCESS;

  CeedDebug256(ceed, 1, "---------- CeedOperator Fallback ----------\n");
  CeedDebug(ceed, "Creating fallback CeedOperator\n");

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

    CeedCall(CeedCompositeOperatorCreate(ceed_fallback, &op_fallback));
    CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators));
    CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators));
    for (CeedInt i = 0; i < num_suboperators; i++) {
      CeedOperator op_sub_fallback;

      CeedCall(CeedOperatorGetFallback(sub_operators[i], &op_sub_fallback));
      CeedCall(CeedCompositeOperatorAddSub(op_fallback, op_sub_fallback));
    }
  } else {
    CeedInt            num_input_fields, num_output_fields;
    CeedQFunction      qf_fallback = NULL, dqf_fallback = NULL, dqfT_fallback = NULL;
    CeedOperatorField *input_fields, *output_fields;

    CeedCall(CeedQFunctionCreateFallback(ceed_fallback, op->qf, &qf_fallback));
    CeedCall(CeedQFunctionCreateFallback(ceed_fallback, op->dqf, &dqf_fallback));
    CeedCall(CeedQFunctionCreateFallback(ceed_fallback, op->dqfT, &dqfT_fallback));
    CeedCall(CeedOperatorCreate(ceed_fallback, qf_fallback, dqf_fallback, dqfT_fallback, &op_fallback));
    CeedCall(CeedOperatorGetFields(op, &num_input_fields, &input_fields, &num_output_fields, &output_fields));
    for (CeedInt i = 0; i < num_input_fields; i++) {
      const char         *field_name;
      CeedVector          vec;
      CeedElemRestriction rstr;
      CeedBasis           basis;

      CeedCall(CeedOperatorFieldGetData(input_fields[i], &field_name, &rstr, &basis, &vec));
      CeedCall(CeedOperatorSetField(op_fallback, field_name, rstr, basis, vec));
      CeedCall(CeedVectorDestroy(&vec));
      CeedCall(CeedElemRestrictionDestroy(&rstr));
      CeedCall(CeedBasisDestroy(&basis));
    }
    for (CeedInt i = 0; i < num_output_fields; i++) {
      const char         *field_name;
      CeedVector          vec;
      CeedElemRestriction rstr;
      CeedBasis           basis;

      CeedCall(CeedOperatorFieldGetData(output_fields[i], &field_name, &rstr, &basis, &vec));
      CeedCall(CeedOperatorSetField(op_fallback, field_name, rstr, basis, vec));
      CeedCall(CeedVectorDestroy(&vec));
      CeedCall(CeedElemRestrictionDestroy(&rstr));
      CeedCall(CeedBasisDestroy(&basis));
    }
    {
      CeedQFunctionAssemblyData data;

      CeedCall(CeedOperatorGetQFunctionAssemblyData(op, &data));
      CeedCall(CeedQFunctionAssemblyDataReferenceCopy(data, &op_fallback->qf_assembled));
    }
    // Cleanup
    CeedCall(CeedQFunctionDestroy(&qf_fallback));
    CeedCall(CeedQFunctionDestroy(&dqf_fallback));
    CeedCall(CeedQFunctionDestroy(&dqfT_fallback));
  }
  CeedCall(CeedOperatorSetName(op_fallback, op->name));
  CeedCall(CeedOperatorCheckReady(op_fallback));
  // Note: No ref-counting here so we don't get caught in a reference loop.
  //       The op holds the only reference to op_fallback and is responsible for deleting itself and op_fallback.
  op->op_fallback                 = op_fallback;
  op_fallback->op_fallback_parent = op;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Core logic for assembling operator diagonal or point block diagonal

  @param[in]  op             `CeedOperator` to assemble diagonal or point block diagonal
  @param[in]  request        Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE
  @param[in]  is_point_block Boolean flag to assemble diagonal or point block diagonal
  @param[out] assembled      `CeedVector` to store assembled diagonal

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

  @ref Developer
**/
static inline int CeedSingleOperatorLinearAssembleAddDiagonal_Mesh(CeedOperator op, CeedRequest *request, const bool is_point_block,
                                                                   CeedVector assembled) {
  Ceed ceed;
  bool is_composite;

  CeedCall(CeedOperatorGetCeed(op, &ceed));
  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCheck(!is_composite, ceed, CEED_ERROR_UNSUPPORTED, "Composite operator not supported");

  // Assemble QFunction
  CeedInt             layout_qf[3];
  const CeedScalar   *assembled_qf_array;
  CeedVector          assembled_qf        = NULL;
  CeedElemRestriction assembled_elem_rstr = NULL;

  CeedCall(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &assembled_elem_rstr, request));
  CeedCall(CeedElemRestrictionGetELayout(assembled_elem_rstr, layout_qf));
  CeedCall(CeedElemRestrictionDestroy(&assembled_elem_rstr));
  CeedCall(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_HOST, &assembled_qf_array));

  // Get assembly data
  const CeedEvalMode     **eval_modes_in, **eval_modes_out;
  CeedInt                  num_active_bases_in, *num_eval_modes_in, num_active_bases_out, *num_eval_modes_out;
  CeedSize               **eval_mode_offsets_in, **eval_mode_offsets_out, num_output_components;
  CeedBasis               *active_bases_in, *active_bases_out;
  CeedElemRestriction     *active_elem_rstrs_in, *active_elem_rstrs_out;
  CeedOperatorAssemblyData data;

  CeedCall(CeedOperatorGetOperatorAssemblyData(op, &data));
  CeedCall(CeedOperatorAssemblyDataGetEvalModes(data, &num_active_bases_in, &num_eval_modes_in, &eval_modes_in, &eval_mode_offsets_in,
                                                &num_active_bases_out, &num_eval_modes_out, &eval_modes_out, &eval_mode_offsets_out,
                                                &num_output_components));
  CeedCall(CeedOperatorAssemblyDataGetBases(data, NULL, &active_bases_in, NULL, NULL, &active_bases_out, NULL));
  CeedCall(CeedOperatorAssemblyDataGetElemRestrictions(data, NULL, &active_elem_rstrs_in, NULL, &active_elem_rstrs_out));

  // Loop over all active bases (find matching input/output pairs)
  for (CeedInt b = 0; b < CeedIntMin(num_active_bases_in, num_active_bases_out); b++) {
    CeedInt             b_in, b_out, num_elem, num_nodes, num_qpts, num_comp;
    bool                has_eval_none = false;
    CeedScalar         *elem_diag_array, *identity = NULL;
    CeedVector          elem_diag;
    CeedElemRestriction diag_elem_rstr;

    if (num_active_bases_in <= num_active_bases_out) {
      b_in = b;
      for (b_out = 0; b_out < num_active_bases_out; b_out++) {
        if (active_bases_in[b_in] == active_bases_out[b_out]) {
          break;
        }
      }
      if (b_out == num_active_bases_out) {
        continue;
      }  // No matching output basis found
    } else {
      b_out = b;
      for (b_in = 0; b_in < num_active_bases_in; b_in++) {
        if (active_bases_in[b_in] == active_bases_out[b_out]) {
          break;
        }
      }
      if (b_in == num_active_bases_in) {
        continue;
      }  // No matching output basis found
    }
    CeedCheck(active_elem_rstrs_in[b_in] == active_elem_rstrs_out[b_out], ceed, CEED_ERROR_UNSUPPORTED,
              "Cannot assemble operator diagonal with different input and output active element restrictions");

    // Assemble point block diagonal restriction, if needed
    if (is_point_block) {
      CeedCall(CeedOperatorCreateActivePointBlockRestriction(active_elem_rstrs_in[b_in], &diag_elem_rstr));
    } else {
      CeedCall(CeedElemRestrictionCreateUnsignedCopy(active_elem_rstrs_in[b_in], &diag_elem_rstr));
    }

    // Create diagonal vector
    CeedCall(CeedElemRestrictionCreateVector(diag_elem_rstr, NULL, &elem_diag));

    // Assemble element operator diagonals
    CeedCall(CeedVectorSetValue(elem_diag, 0.0));
    CeedCall(CeedVectorGetArray(elem_diag, CEED_MEM_HOST, &elem_diag_array));
    CeedCall(CeedElemRestrictionGetNumElements(diag_elem_rstr, &num_elem));
    CeedCall(CeedBasisGetNumNodes(active_bases_in[b_in], &num_nodes));
    CeedCall(CeedBasisGetNumComponents(active_bases_in[b_in], &num_comp));
    if (active_bases_in[b_in] == CEED_BASIS_NONE) num_qpts = num_nodes;
    else CeedCall(CeedBasisGetNumQuadraturePoints(active_bases_in[b_in], &num_qpts));

    // Construct identity matrix for basis if required
    for (CeedInt i = 0; i < num_eval_modes_in[b_in]; i++) {
      has_eval_none = has_eval_none || (eval_modes_in[b_in][i] == CEED_EVAL_NONE);
    }
    for (CeedInt i = 0; i < num_eval_modes_out[b_out]; i++) {
      has_eval_none = has_eval_none || (eval_modes_out[b_out][i] == CEED_EVAL_NONE);
    }
    if (has_eval_none) {
      CeedCall(CeedCalloc(num_qpts * num_nodes, &identity));
      for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) identity[i * num_nodes + i] = 1.0;
    }

    // Compute the diagonal of B^T D B
    // Each element
    for (CeedSize e = 0; e < num_elem; e++) {
      // Each basis eval mode pair
      CeedInt      d_out              = 0, q_comp_out;
      CeedEvalMode eval_mode_out_prev = CEED_EVAL_NONE;

      for (CeedInt e_out = 0; e_out < num_eval_modes_out[b_out]; e_out++) {
        CeedInt           d_in              = 0, q_comp_in;
        const CeedScalar *B_t               = NULL;
        CeedEvalMode      eval_mode_in_prev = CEED_EVAL_NONE;

        CeedCall(CeedOperatorGetBasisPointer(active_bases_out[b_out], eval_modes_out[b_out][e_out], identity, &B_t));
        CeedCall(CeedBasisGetNumQuadratureComponents(active_bases_out[b_out], eval_modes_out[b_out][e_out], &q_comp_out));
        if (q_comp_out > 1) {
          if (e_out == 0 || eval_modes_out[b_out][e_out] != eval_mode_out_prev) d_out = 0;
          else B_t = &B_t[(++d_out) * num_qpts * num_nodes];
        }
        eval_mode_out_prev = eval_modes_out[b_out][e_out];

        for (CeedInt e_in = 0; e_in < num_eval_modes_in[b_in]; e_in++) {
          const CeedScalar *B = NULL;

          CeedCall(CeedOperatorGetBasisPointer(active_bases_in[b_in], eval_modes_in[b_in][e_in], identity, &B));
          CeedCall(CeedBasisGetNumQuadratureComponents(active_bases_in[b_in], eval_modes_in[b_in][e_in], &q_comp_in));
          if (q_comp_in > 1) {
            if (e_in == 0 || eval_modes_in[b_in][e_in] != eval_mode_in_prev) d_in = 0;
            else B = &B[(++d_in) * num_qpts * num_nodes];
          }
          eval_mode_in_prev = eval_modes_in[b_in][e_in];

          // Each component
          for (CeedInt c_out = 0; c_out < num_comp; c_out++) {
            // Each qpt/node pair
            for (CeedInt q = 0; q < num_qpts; q++) {
              if (is_point_block) {
                // Point Block Diagonal
                for (CeedInt c_in = 0; c_in < num_comp; c_in++) {
                  const CeedSize c_offset =
                      (eval_mode_offsets_in[b_in][e_in] + c_in) * num_output_components + eval_mode_offsets_out[b_out][e_out] + c_out;
                  const CeedScalar qf_value = assembled_qf_array[q * layout_qf[0] + c_offset * layout_qf[1] + e * layout_qf[2]];

                  for (CeedInt n = 0; n < num_nodes; n++) {
                    elem_diag_array[((e * num_comp + c_out) * num_comp + c_in) * num_nodes + n] +=
                        B_t[q * num_nodes + n] * qf_value * B[q * num_nodes + n];
                  }
                }
              } else {
                // Diagonal Only
                const CeedInt c_offset =
                    (eval_mode_offsets_in[b_in][e_in] + c_out) * num_output_components + eval_mode_offsets_out[b_out][e_out] + c_out;
                const CeedScalar qf_value = assembled_qf_array[q * layout_qf[0] + c_offset * layout_qf[1] + e * layout_qf[2]];

                for (CeedInt n = 0; n < num_nodes; n++) {
                  elem_diag_array[(e * num_comp + c_out) * num_nodes + n] += B_t[q * num_nodes + n] * qf_value * B[q * num_nodes + n];
                }
              }
            }
          }
        }
      }
    }
    CeedCall(CeedVectorRestoreArray(elem_diag, &elem_diag_array));

    // Assemble local operator diagonal
    CeedCall(CeedElemRestrictionApply(diag_elem_rstr, CEED_TRANSPOSE, elem_diag, assembled, request));

    // Cleanup
    CeedCall(CeedElemRestrictionDestroy(&diag_elem_rstr));
    CeedCall(CeedVectorDestroy(&elem_diag));
    CeedCall(CeedFree(&identity));
  }
  CeedCall(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array));
  CeedCall(CeedVectorDestroy(&assembled_qf));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Core logic for assembling operator diagonal or point block diagonal

  @param[in]  op             `CeedOperator` to assemble diagonal or point block diagonal
  @param[in]  request        Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE
  @param[in]  is_point_block Boolean flag to assemble diagonal or point block diagonal
  @param[out] assembled      `CeedVector` to store assembled diagonal

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

  @ref Developer
**/
static inline int CeedSingleOperatorLinearAssembleAddDiagonal(CeedOperator op, CeedRequest *request, const bool is_point_block,
                                                              CeedVector assembled) {
  Ceed ceed;
  bool is_at_points;

  CeedCall(CeedOperatorGetCeed(op, &ceed));
  CeedCall(CeedOperatorIsAtPoints(op, &is_at_points));
  CeedCheck(!is_at_points, ceed, CEED_ERROR_UNSUPPORTED, "AtPoints operator not supported");
  CeedCall(CeedSingleOperatorLinearAssembleAddDiagonal_Mesh(op, request, is_point_block, assembled));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Core logic for assembling composite operator diagonal

  @param[in]  op             `CeedOperator` to assemble point block diagonal
  @param[in]  request        Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE
  @param[in]  is_point_block Boolean flag to assemble diagonal or point block diagonal
  @param[out] assembled      `CeedVector` to store assembled diagonal

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

  @ref Developer
**/
static inline int CeedCompositeOperatorLinearAssembleAddDiagonal(CeedOperator op, CeedRequest *request, const bool is_point_block,
                                                                 CeedVector assembled) {
  CeedInt       num_sub;
  CeedOperator *suboperators;

  CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub));
  CeedCall(CeedCompositeOperatorGetSubList(op, &suboperators));
  for (CeedInt i = 0; i < num_sub; i++) {
    if (is_point_block) {
      CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(suboperators[i], assembled, request));
    } else {
      CeedCall(CeedOperatorLinearAssembleAddDiagonal(suboperators[i], assembled, request));
    }
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Build nonzero pattern for non-composite CeedOperator`.

  Users should generally use @ref CeedOperatorLinearAssembleSymbolic().

  @param[in]  op     `CeedOperator` to assemble nonzero pattern
  @param[in]  offset Offset for number of entries
  @param[out] rows   Row number for each entry
  @param[out] cols   Column number for each entry

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

  @ref Developer
**/
static int CeedSingleOperatorAssembleSymbolic(CeedOperator op, CeedInt offset, CeedInt *rows, CeedInt *cols) {
  Ceed                ceed;
  bool                is_composite;
  CeedSize            num_nodes_in, num_nodes_out, local_num_entries, count = 0;
  CeedInt             num_elem_in, elem_size_in, num_comp_in, layout_er_in[3];
  CeedInt             num_elem_out, elem_size_out, num_comp_out, layout_er_out[3];
  CeedScalar         *array;
  const CeedScalar   *elem_dof_a_in, *elem_dof_a_out;
  CeedVector          index_vec_in, index_vec_out, elem_dof_in, elem_dof_out;
  CeedElemRestriction elem_rstr_in, elem_rstr_out, index_elem_rstr_in, index_elem_rstr_out;

  CeedCall(CeedOperatorGetCeed(op, &ceed));
  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCheck(!is_composite, ceed, CEED_ERROR_UNSUPPORTED, "Composite operator not supported");

  CeedCall(CeedOperatorGetActiveVectorLengths(op, &num_nodes_in, &num_nodes_out));
  CeedCall(CeedOperatorGetActiveElemRestrictions(op, &elem_rstr_in, &elem_rstr_out));
  CeedCall(CeedElemRestrictionGetNumElements(elem_rstr_in, &num_elem_in));
  CeedCall(CeedElemRestrictionGetElementSize(elem_rstr_in, &elem_size_in));
  CeedCall(CeedElemRestrictionGetNumComponents(elem_rstr_in, &num_comp_in));
  CeedCall(CeedElemRestrictionGetELayout(elem_rstr_in, layout_er_in));

  // Determine elem_dof relation for input
  CeedCall(CeedVectorCreate(ceed, num_nodes_in, &index_vec_in));
  CeedCall(CeedVectorGetArrayWrite(index_vec_in, CEED_MEM_HOST, &array));
  for (CeedSize i = 0; i < num_nodes_in; i++) array[i] = i;
  CeedCall(CeedVectorRestoreArray(index_vec_in, &array));
  CeedCall(CeedVectorCreate(ceed, num_elem_in * elem_size_in * num_comp_in, &elem_dof_in));
  CeedCall(CeedVectorSetValue(elem_dof_in, 0.0));
  CeedCall(CeedElemRestrictionCreateUnorientedCopy(elem_rstr_in, &index_elem_rstr_in));
  CeedCall(CeedElemRestrictionApply(index_elem_rstr_in, CEED_NOTRANSPOSE, index_vec_in, elem_dof_in, CEED_REQUEST_IMMEDIATE));
  CeedCall(CeedVectorGetArrayRead(elem_dof_in, CEED_MEM_HOST, &elem_dof_a_in));
  CeedCall(CeedVectorDestroy(&index_vec_in));
  CeedCall(CeedElemRestrictionDestroy(&index_elem_rstr_in));

  if (elem_rstr_in != elem_rstr_out) {
    CeedCall(CeedElemRestrictionGetNumElements(elem_rstr_out, &num_elem_out));
    CeedCheck(num_elem_in == num_elem_out, ceed, CEED_ERROR_UNSUPPORTED,
              "Active input and output operator restrictions must have the same number of elements."
              " Input has %" CeedInt_FMT " elements; output has %" CeedInt_FMT "elements.",
              num_elem_in, num_elem_out);
    CeedCall(CeedElemRestrictionGetElementSize(elem_rstr_out, &elem_size_out));
    CeedCall(CeedElemRestrictionGetNumComponents(elem_rstr_out, &num_comp_out));
    CeedCall(CeedElemRestrictionGetELayout(elem_rstr_out, layout_er_out));

    // Determine elem_dof relation for output
    CeedCall(CeedVectorCreate(ceed, num_nodes_out, &index_vec_out));
    CeedCall(CeedVectorGetArrayWrite(index_vec_out, CEED_MEM_HOST, &array));
    for (CeedSize i = 0; i < num_nodes_out; i++) array[i] = i;
    CeedCall(CeedVectorRestoreArray(index_vec_out, &array));
    CeedCall(CeedVectorCreate(ceed, num_elem_out * elem_size_out * num_comp_out, &elem_dof_out));
    CeedCall(CeedVectorSetValue(elem_dof_out, 0.0));
    CeedCall(CeedElemRestrictionCreateUnorientedCopy(elem_rstr_out, &index_elem_rstr_out));
    CeedCall(CeedElemRestrictionApply(index_elem_rstr_out, CEED_NOTRANSPOSE, index_vec_out, elem_dof_out, CEED_REQUEST_IMMEDIATE));
    CeedCall(CeedVectorGetArrayRead(elem_dof_out, CEED_MEM_HOST, &elem_dof_a_out));
    CeedCall(CeedVectorDestroy(&index_vec_out));
    CeedCall(CeedElemRestrictionDestroy(&index_elem_rstr_out));
  } else {
    num_elem_out     = num_elem_in;
    elem_size_out    = elem_size_in;
    num_comp_out     = num_comp_in;
    layout_er_out[0] = layout_er_in[0];
    layout_er_out[1] = layout_er_in[1];
    layout_er_out[2] = layout_er_in[2];
    elem_dof_a_out   = elem_dof_a_in;
  }
  local_num_entries = (CeedSize)elem_size_out * num_comp_out * elem_size_in * num_comp_in * num_elem_in;

  // Determine i, j locations for element matrices
  for (CeedInt e = 0; e < num_elem_in; e++) {
    for (CeedInt comp_in = 0; comp_in < num_comp_in; comp_in++) {
      for (CeedInt comp_out = 0; comp_out < num_comp_out; comp_out++) {
        for (CeedInt i = 0; i < elem_size_out; i++) {
          for (CeedInt j = 0; j < elem_size_in; j++) {
            const CeedInt elem_dof_index_row = i * layout_er_out[0] + comp_out * layout_er_out[1] + e * layout_er_out[2];
            const CeedInt elem_dof_index_col = j * layout_er_in[0] + comp_in * layout_er_in[1] + e * layout_er_in[2];
            const CeedInt row                = elem_dof_a_out[elem_dof_index_row];
            const CeedInt col                = elem_dof_a_in[elem_dof_index_col];

            rows[offset + count] = row;
            cols[offset + count] = col;
            count++;
          }
        }
      }
    }
  }
  CeedCheck(count == local_num_entries, ceed, CEED_ERROR_MAJOR, "Error computing assembled entries");
  CeedCall(CeedVectorRestoreArrayRead(elem_dof_in, &elem_dof_a_in));
  CeedCall(CeedVectorDestroy(&elem_dof_in));
  if (elem_rstr_in != elem_rstr_out) {
    CeedCall(CeedVectorRestoreArrayRead(elem_dof_out, &elem_dof_a_out));
    CeedCall(CeedVectorDestroy(&elem_dof_out));
  }
  CeedCall(CeedElemRestrictionDestroy(&elem_rstr_in));
  CeedCall(CeedElemRestrictionDestroy(&elem_rstr_out));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Assemble nonzero entries for non-composite `CeedOperator`.

  Users should generally use @ref CeedOperatorLinearAssemble().

  @param[in]  op     `CeedOperator` to assemble
  @param[in]  offset Offset for number of entries
  @param[out] values Values to assemble into matrix

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

  @ref Developer
**/
static int CeedSingleOperatorAssemble(CeedOperator op, CeedInt offset, CeedVector values) {
  Ceed ceed;
  bool is_composite;

  CeedCall(CeedOperatorGetCeed(op, &ceed));
  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCheck(!is_composite, ceed, CEED_ERROR_UNSUPPORTED, "Composite operator not supported");

  // Early exit for empty operator
  {
    CeedInt num_elem = 0;

    CeedCall(CeedOperatorGetNumElements(op, &num_elem));
    if (num_elem == 0) return CEED_ERROR_SUCCESS;
  }

  if (op->LinearAssembleSingle) {
    // Backend version
    CeedCall(op->LinearAssembleSingle(op, offset, values));
    return CEED_ERROR_SUCCESS;
  } else {
    // Operator fallback
    CeedOperator op_fallback;

    CeedCall(CeedOperatorGetFallback(op, &op_fallback));
    if (op_fallback) {
      CeedCall(CeedSingleOperatorAssemble(op_fallback, offset, values));
      return CEED_ERROR_SUCCESS;
    }
  }

  // Assemble QFunction
  CeedInt             layout_qf[3];
  const CeedScalar   *assembled_qf_array;
  CeedVector          assembled_qf        = NULL;
  CeedElemRestriction assembled_elem_rstr = NULL;

  CeedCall(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &assembled_elem_rstr, CEED_REQUEST_IMMEDIATE));
  CeedCall(CeedElemRestrictionGetELayout(assembled_elem_rstr, layout_qf));
  CeedCall(CeedElemRestrictionDestroy(&assembled_elem_rstr));
  CeedCall(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_HOST, &assembled_qf_array));

  // Get assembly data
  CeedInt                  num_elem_in, elem_size_in, num_comp_in, num_qpts_in;
  CeedInt                  num_elem_out, elem_size_out, num_comp_out, num_qpts_out;
  CeedSize                 local_num_entries, count = 0;
  const CeedEvalMode     **eval_modes_in, **eval_modes_out;
  CeedInt                  num_active_bases_in, *num_eval_modes_in, num_active_bases_out, *num_eval_modes_out;
  CeedBasis               *active_bases_in, *active_bases_out, basis_in, basis_out;
  const CeedScalar       **B_mats_in, **B_mats_out, *B_mat_in, *B_mat_out;
  CeedElemRestriction      elem_rstr_in, elem_rstr_out;
  CeedRestrictionType      elem_rstr_type_in, elem_rstr_type_out;
  const bool              *elem_rstr_orients_in = NULL, *elem_rstr_orients_out = NULL;
  const CeedInt8          *elem_rstr_curl_orients_in = NULL, *elem_rstr_curl_orients_out = NULL;
  CeedOperatorAssemblyData data;

  CeedCall(CeedOperatorGetOperatorAssemblyData(op, &data));
  CeedCall(CeedOperatorAssemblyDataGetEvalModes(data, &num_active_bases_in, &num_eval_modes_in, &eval_modes_in, NULL, &num_active_bases_out,
                                                &num_eval_modes_out, &eval_modes_out, NULL, NULL));

  CeedCheck(num_active_bases_in == 1 && num_active_bases_out == 1, ceed, CEED_ERROR_UNSUPPORTED,
            "Cannot assemble operator with multiple active bases");
  CeedCheck(num_eval_modes_in[0] > 0 && num_eval_modes_out[0] > 0, ceed, CEED_ERROR_UNSUPPORTED, "Cannot assemble operator without inputs/outputs");

  CeedCall(CeedOperatorAssemblyDataGetBases(data, NULL, &active_bases_in, &B_mats_in, NULL, &active_bases_out, &B_mats_out));
  CeedCall(CeedOperatorGetActiveElemRestrictions(op, &elem_rstr_in, &elem_rstr_out));
  basis_in  = active_bases_in[0];
  basis_out = active_bases_out[0];
  B_mat_in  = B_mats_in[0];
  B_mat_out = B_mats_out[0];

  CeedCall(CeedElemRestrictionGetNumElements(elem_rstr_in, &num_elem_in));
  CeedCall(CeedElemRestrictionGetElementSize(elem_rstr_in, &elem_size_in));
  CeedCall(CeedElemRestrictionGetNumComponents(elem_rstr_in, &num_comp_in));
  if (basis_in == CEED_BASIS_NONE) num_qpts_in = elem_size_in;
  else CeedCall(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts_in));

  CeedCall(CeedElemRestrictionGetType(elem_rstr_in, &elem_rstr_type_in));
  if (elem_rstr_type_in == CEED_RESTRICTION_ORIENTED) {
    CeedCall(CeedElemRestrictionGetOrientations(elem_rstr_in, CEED_MEM_HOST, &elem_rstr_orients_in));
  } else if (elem_rstr_type_in == CEED_RESTRICTION_CURL_ORIENTED) {
    CeedCall(CeedElemRestrictionGetCurlOrientations(elem_rstr_in, CEED_MEM_HOST, &elem_rstr_curl_orients_in));
  }

  if (elem_rstr_in != elem_rstr_out) {
    CeedCall(CeedElemRestrictionGetNumElements(elem_rstr_out, &num_elem_out));
    CeedCheck(num_elem_in == num_elem_out, ceed, CEED_ERROR_UNSUPPORTED,
              "Active input and output operator restrictions must have the same number of elements."
              " Input has %" CeedInt_FMT " elements; output has %" CeedInt_FMT "elements.",
              num_elem_in, num_elem_out);
    CeedCall(CeedElemRestrictionGetElementSize(elem_rstr_out, &elem_size_out));
    CeedCall(CeedElemRestrictionGetNumComponents(elem_rstr_out, &num_comp_out));
    if (basis_out == CEED_BASIS_NONE) num_qpts_out = elem_size_out;
    else CeedCall(CeedBasisGetNumQuadraturePoints(basis_out, &num_qpts_out));
    CeedCheck(num_qpts_in == num_qpts_out, ceed, CEED_ERROR_UNSUPPORTED,
              "Active input and output bases must have the same number of quadrature points."
              " Input has %" CeedInt_FMT " points; output has %" CeedInt_FMT "points.",
              num_qpts_in, num_qpts_out);

    CeedCall(CeedElemRestrictionGetType(elem_rstr_out, &elem_rstr_type_out));
    if (elem_rstr_type_out == CEED_RESTRICTION_ORIENTED) {
      CeedCall(CeedElemRestrictionGetOrientations(elem_rstr_out, CEED_MEM_HOST, &elem_rstr_orients_out));
    } else if (elem_rstr_type_out == CEED_RESTRICTION_CURL_ORIENTED) {
      CeedCall(CeedElemRestrictionGetCurlOrientations(elem_rstr_out, CEED_MEM_HOST, &elem_rstr_curl_orients_out));
    }
  } else {
    num_elem_out  = num_elem_in;
    elem_size_out = elem_size_in;
    num_comp_out  = num_comp_in;
    num_qpts_out  = num_qpts_in;

    elem_rstr_orients_out      = elem_rstr_orients_in;
    elem_rstr_curl_orients_out = elem_rstr_curl_orients_in;
  }
  local_num_entries = (CeedSize)elem_size_out * num_comp_out * elem_size_in * num_comp_in * num_elem_in;

  // Loop over elements and put in data structure
  // We store B_mat_in, B_mat_out, BTD, elem_mat in row-major order
  CeedTensorContract contract;
  CeedScalar        *vals, *BTD_mat = NULL, *elem_mat = NULL, *elem_mat_b = NULL;

  CeedCall(CeedBasisGetTensorContract(basis_in, &contract));
  CeedCall(CeedCalloc(elem_size_out * num_qpts_in * num_eval_modes_in[0], &BTD_mat));
  CeedCall(CeedCalloc(elem_size_out * elem_size_in, &elem_mat));
  if (elem_rstr_curl_orients_in || elem_rstr_curl_orients_out) CeedCall(CeedCalloc(elem_size_out * elem_size_in, &elem_mat_b));

  CeedCall(CeedVectorGetArray(values, CEED_MEM_HOST, &vals));
  for (CeedSize e = 0; e < num_elem_in; e++) {
    for (CeedInt comp_in = 0; comp_in < num_comp_in; comp_in++) {
      for (CeedInt comp_out = 0; comp_out < num_comp_out; comp_out++) {
        // Compute B^T*D
        for (CeedSize n = 0; n < elem_size_out; n++) {
          for (CeedSize q = 0; q < num_qpts_in; q++) {
            for (CeedInt e_in = 0; e_in < num_eval_modes_in[0]; e_in++) {
              const CeedSize btd_index = n * (num_qpts_in * num_eval_modes_in[0]) + q * num_eval_modes_in[0] + e_in;
              CeedScalar     sum       = 0.0;

              for (CeedInt e_out = 0; e_out < num_eval_modes_out[0]; e_out++) {
                const CeedSize b_out_index     = (q * num_eval_modes_out[0] + e_out) * elem_size_out + n;
                const CeedSize eval_mode_index = ((e_in * num_comp_in + comp_in) * num_eval_modes_out[0] + e_out) * num_comp_out + comp_out;
                const CeedSize qf_index        = q * layout_qf[0] + eval_mode_index * layout_qf[1] + e * layout_qf[2];

                sum += B_mat_out[b_out_index] * assembled_qf_array[qf_index];
              }
              BTD_mat[btd_index] = sum;
            }
          }
        }

        // Form element matrix itself (for each block component)
        if (contract) {
          CeedCall(CeedTensorContractApply(contract, 1, num_qpts_in * num_eval_modes_in[0], elem_size_in, elem_size_out, BTD_mat, CEED_NOTRANSPOSE,
                                           false, B_mat_in, elem_mat));
        } else {
          CeedCall(CeedMatrixMatrixMultiply(ceed, BTD_mat, B_mat_in, elem_mat, elem_size_out, elem_size_in, num_qpts_in * num_eval_modes_in[0]));
        }

        // Transform the element matrix if required
        if (elem_rstr_orients_out) {
          const bool *elem_orients = &elem_rstr_orients_out[e * elem_size_out];

          for (CeedInt i = 0; i < elem_size_out; i++) {
            const double orient = elem_orients[i] ? -1.0 : 1.0;

            for (CeedInt j = 0; j < elem_size_in; j++) {
              elem_mat[i * elem_size_in + j] *= orient;
            }
          }
        } else if (elem_rstr_curl_orients_out) {
          const CeedInt8 *elem_curl_orients = &elem_rstr_curl_orients_out[e * 3 * elem_size_out];

          // T^T*(B^T*D*B)
          memcpy(elem_mat_b, elem_mat, elem_size_out * elem_size_in * sizeof(CeedScalar));
          for (CeedInt i = 0; i < elem_size_out; i++) {
            for (CeedInt j = 0; j < elem_size_in; j++) {
              elem_mat[i * elem_size_in + j] = elem_mat_b[i * elem_size_in + j] * elem_curl_orients[3 * i + 1] +
                                               (i > 0 ? elem_mat_b[(i - 1) * elem_size_in + j] * elem_curl_orients[3 * i - 1] : 0.0) +
                                               (i < elem_size_out - 1 ? elem_mat_b[(i + 1) * elem_size_in + j] * elem_curl_orients[3 * i + 3] : 0.0);
            }
          }
        }
        if (elem_rstr_orients_in) {
          const bool *elem_orients = &elem_rstr_orients_in[e * elem_size_in];

          for (CeedInt i = 0; i < elem_size_out; i++) {
            for (CeedInt j = 0; j < elem_size_in; j++) {
              elem_mat[i * elem_size_in + j] *= elem_orients[j] ? -1.0 : 1.0;
            }
          }
        } else if (elem_rstr_curl_orients_in) {
          const CeedInt8 *elem_curl_orients = &elem_rstr_curl_orients_in[e * 3 * elem_size_in];

          // (B^T*D*B)*T
          memcpy(elem_mat_b, elem_mat, elem_size_out * elem_size_in * sizeof(CeedScalar));
          for (CeedInt i = 0; i < elem_size_out; i++) {
            for (CeedInt j = 0; j < elem_size_in; j++) {
              elem_mat[i * elem_size_in + j] = elem_mat_b[i * elem_size_in + j] * elem_curl_orients[3 * j + 1] +
                                               (j > 0 ? elem_mat_b[i * elem_size_in + j - 1] * elem_curl_orients[3 * j - 1] : 0.0) +
                                               (j < elem_size_in - 1 ? elem_mat_b[i * elem_size_in + j + 1] * elem_curl_orients[3 * j + 3] : 0.0);
            }
          }
        }

        // Put element matrix in coordinate data structure
        for (CeedInt i = 0; i < elem_size_out; i++) {
          for (CeedInt j = 0; j < elem_size_in; j++) {
            vals[offset + count] = elem_mat[i * elem_size_in + j];
            count++;
          }
        }
      }
    }
  }
  CeedCheck(count == local_num_entries, ceed, CEED_ERROR_MAJOR, "Error computing entries");
  CeedCall(CeedVectorRestoreArray(values, &vals));

  // Cleanup
  CeedCall(CeedFree(&BTD_mat));
  CeedCall(CeedFree(&elem_mat));
  CeedCall(CeedFree(&elem_mat_b));
  if (elem_rstr_type_in == CEED_RESTRICTION_ORIENTED) {
    CeedCall(CeedElemRestrictionRestoreOrientations(elem_rstr_in, &elem_rstr_orients_in));
  } else if (elem_rstr_type_in == CEED_RESTRICTION_CURL_ORIENTED) {
    CeedCall(CeedElemRestrictionRestoreCurlOrientations(elem_rstr_in, &elem_rstr_curl_orients_in));
  }
  if (elem_rstr_in != elem_rstr_out) {
    if (elem_rstr_type_out == CEED_RESTRICTION_ORIENTED) {
      CeedCall(CeedElemRestrictionRestoreOrientations(elem_rstr_out, &elem_rstr_orients_out));
    } else if (elem_rstr_type_out == CEED_RESTRICTION_CURL_ORIENTED) {
      CeedCall(CeedElemRestrictionRestoreCurlOrientations(elem_rstr_out, &elem_rstr_curl_orients_out));
    }
  }
  CeedCall(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array));
  CeedCall(CeedVectorDestroy(&assembled_qf));
  CeedCall(CeedElemRestrictionDestroy(&elem_rstr_in));
  CeedCall(CeedElemRestrictionDestroy(&elem_rstr_out));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Count number of entries for assembled `CeedOperator`

  @param[in]  op          `CeedOperator` to assemble
  @param[out] num_entries Number of entries in assembled representation

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

  @ref Utility
**/
static int CeedSingleOperatorAssemblyCountEntries(CeedOperator op, CeedSize *num_entries) {
  bool                is_composite;
  CeedInt             num_elem_in, elem_size_in, num_comp_in, num_elem_out, elem_size_out, num_comp_out;
  Ceed                ceed;
  CeedElemRestriction rstr_in, rstr_out;

  CeedCall(CeedOperatorGetCeed(op, &ceed));
  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCheck(!is_composite, ceed, CEED_ERROR_UNSUPPORTED, "Composite operator not supported");

  CeedCall(CeedOperatorGetActiveElemRestrictions(op, &rstr_in, &rstr_out));
  CeedCall(CeedElemRestrictionGetNumElements(rstr_in, &num_elem_in));
  CeedCall(CeedElemRestrictionGetElementSize(rstr_in, &elem_size_in));
  CeedCall(CeedElemRestrictionGetNumComponents(rstr_in, &num_comp_in));
  if (rstr_in != rstr_out) {
    CeedCall(CeedElemRestrictionGetNumElements(rstr_out, &num_elem_out));
    CeedCheck(num_elem_in == num_elem_out, ceed, CEED_ERROR_UNSUPPORTED,
              "Active input and output operator restrictions must have the same number of elements."
              " Input has %" CeedInt_FMT " elements; output has %" CeedInt_FMT "elements.",
              num_elem_in, num_elem_out);
    CeedCall(CeedElemRestrictionGetElementSize(rstr_out, &elem_size_out));
    CeedCall(CeedElemRestrictionGetNumComponents(rstr_out, &num_comp_out));
  } else {
    num_elem_out  = num_elem_in;
    elem_size_out = elem_size_in;
    num_comp_out  = num_comp_in;
  }
  CeedCall(CeedElemRestrictionDestroy(&rstr_in));
  CeedCall(CeedElemRestrictionDestroy(&rstr_out));
  *num_entries = (CeedSize)elem_size_in * num_comp_in * elem_size_out * num_comp_out * num_elem_in;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Common code for creating a multigrid coarse `CeedOperator` and level transfer `CeedOperator` for a `CeedOperator`

  @param[in]  op_fine      Fine grid `CeedOperator`
  @param[in]  p_mult_fine  L-vector multiplicity in parallel gather/scatter, or `NULL` if not creating prolongation/restriction `CeedOperator`
  @param[in]  rstr_coarse  Coarse grid `CeedElemRestriction`
  @param[in]  basis_coarse Coarse grid active vector `CeedBasis`
  @param[in]  basis_c_to_f `CeedBasis` for coarse to fine interpolation, or `NULL` if not creating prolongation/restriction operators
  @param[out] op_coarse    Coarse grid `CeedOperator`
  @param[out] op_prolong   Coarse to fine `CeedOperator`, or `NULL`
  @param[out] op_restrict  Fine to coarse `CeedOperator`, or `NULL`

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

  @ref Developer
**/
static int CeedSingleOperatorMultigridLevel(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse,
                                            CeedBasis basis_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong, CeedOperator *op_restrict) {
  bool                is_composite;
  Ceed                ceed;
  CeedInt             num_comp, num_input_fields, num_output_fields;
  CeedVector          mult_vec         = NULL;
  CeedElemRestriction rstr_p_mult_fine = NULL, rstr_fine = NULL;
  CeedOperatorField  *input_fields, *output_fields;

  CeedCall(CeedOperatorGetCeed(op_fine, &ceed));

  // Check for composite operator
  CeedCall(CeedOperatorIsComposite(op_fine, &is_composite));
  CeedCheck(!is_composite, ceed, CEED_ERROR_UNSUPPORTED, "Automatic multigrid setup for composite operators not supported");

  // Coarse Grid
  CeedCall(CeedOperatorCreate(ceed, op_fine->qf, op_fine->dqf, op_fine->dqfT, op_coarse));
  CeedCall(CeedOperatorGetFields(op_fine, &num_input_fields, &input_fields, &num_output_fields, &output_fields));
  // -- Clone input fields
  for (CeedInt i = 0; i < num_input_fields; i++) {
    const char         *field_name;
    CeedVector          vec;
    CeedElemRestriction rstr  = NULL;
    CeedBasis           basis = NULL;

    CeedCall(CeedOperatorFieldGetName(input_fields[i], &field_name));
    CeedCall(CeedOperatorFieldGetVector(input_fields[i], &vec));
    if (vec == CEED_VECTOR_ACTIVE) {
      CeedCall(CeedElemRestrictionReferenceCopy(rstr_coarse, &rstr));
      CeedCall(CeedBasisReferenceCopy(basis_coarse, &basis));
      if (!rstr_fine) CeedCall(CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr_fine));
    } else {
      CeedCall(CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr));
      CeedCall(CeedOperatorFieldGetBasis(input_fields[i], &basis));
    }
    CeedCall(CeedOperatorSetField(*op_coarse, field_name, rstr, basis, vec));
    CeedCall(CeedVectorDestroy(&vec));
    CeedCall(CeedElemRestrictionDestroy(&rstr));
    CeedCall(CeedBasisDestroy(&basis));
  }
  // -- Clone output fields
  for (CeedInt i = 0; i < num_output_fields; i++) {
    const char         *field_name;
    CeedVector          vec;
    CeedElemRestriction rstr  = NULL;
    CeedBasis           basis = NULL;

    CeedCall(CeedOperatorFieldGetName(output_fields[i], &field_name));
    CeedCall(CeedOperatorFieldGetVector(output_fields[i], &vec));
    if (vec == CEED_VECTOR_ACTIVE) {
      CeedCall(CeedElemRestrictionReferenceCopy(rstr_coarse, &rstr));
      CeedCall(CeedBasisReferenceCopy(basis_coarse, &basis));
      if (!rstr_fine) CeedCall(CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr_fine));
    } else {
      CeedCall(CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr));
      CeedCall(CeedOperatorFieldGetBasis(output_fields[i], &basis));
    }
    CeedCall(CeedOperatorSetField(*op_coarse, field_name, rstr, basis, vec));
    CeedCall(CeedVectorDestroy(&vec));
    CeedCall(CeedElemRestrictionDestroy(&rstr));
    CeedCall(CeedBasisDestroy(&basis));
  }
  // -- Clone QFunctionAssemblyData
  {
    CeedQFunctionAssemblyData fine_data;

    CeedCall(CeedOperatorGetQFunctionAssemblyData(op_fine, &fine_data));
    CeedCall(CeedQFunctionAssemblyDataReferenceCopy(fine_data, &(*op_coarse)->qf_assembled));
  }

  // Multiplicity vector
  if (op_restrict || op_prolong) {
    CeedVector          mult_e_vec;
    CeedRestrictionType rstr_type;

    CeedCall(CeedElemRestrictionGetType(rstr_fine, &rstr_type));
    CeedCheck(rstr_type != CEED_RESTRICTION_CURL_ORIENTED, ceed, CEED_ERROR_UNSUPPORTED,
              "Element restrictions created with CeedElemRestrictionCreateCurlOriented are not supported");
    CeedCheck(p_mult_fine, ceed, CEED_ERROR_INCOMPATIBLE, "Prolongation or restriction operator creation requires fine grid multiplicity vector");
    CeedCall(CeedElemRestrictionCreateUnsignedCopy(rstr_fine, &rstr_p_mult_fine));
    CeedCall(CeedElemRestrictionCreateVector(rstr_fine, &mult_vec, &mult_e_vec));
    CeedCall(CeedVectorSetValue(mult_e_vec, 0.0));
    CeedCall(CeedElemRestrictionApply(rstr_p_mult_fine, CEED_NOTRANSPOSE, p_mult_fine, mult_e_vec, CEED_REQUEST_IMMEDIATE));
    CeedCall(CeedVectorSetValue(mult_vec, 0.0));
    CeedCall(CeedElemRestrictionApply(rstr_p_mult_fine, CEED_TRANSPOSE, mult_e_vec, mult_vec, CEED_REQUEST_IMMEDIATE));
    CeedCall(CeedVectorDestroy(&mult_e_vec));
    CeedCall(CeedVectorReciprocal(mult_vec));
  }

  // Clone name
  bool   has_name = op_fine->name;
  size_t name_len = op_fine->name ? strlen(op_fine->name) : 0;
  CeedCall(CeedOperatorSetName(*op_coarse, op_fine->name));

  // Check that coarse to fine basis is provided if prolong/restrict operators are requested
  CeedCheck(basis_c_to_f || (!op_restrict && !op_prolong), ceed, CEED_ERROR_INCOMPATIBLE,
            "Prolongation or restriction operator creation requires coarse-to-fine basis");

  // Restriction/Prolongation Operators
  CeedCall(CeedBasisGetNumComponents(basis_coarse, &num_comp));

  // Restriction
  if (op_restrict) {
    CeedInt             *num_comp_r_data;
    CeedQFunctionContext ctx_r;
    CeedQFunction        qf_restrict;

    CeedCall(CeedQFunctionCreateInteriorByName(ceed, "Scale", &qf_restrict));
    CeedCall(CeedCalloc(1, &num_comp_r_data));
    num_comp_r_data[0] = num_comp;
    CeedCall(CeedQFunctionContextCreate(ceed, &ctx_r));
    CeedCall(CeedQFunctionContextSetData(ctx_r, CEED_MEM_HOST, CEED_OWN_POINTER, sizeof(*num_comp_r_data), num_comp_r_data));
    CeedCall(CeedQFunctionSetContext(qf_restrict, ctx_r));
    CeedCall(CeedQFunctionContextDestroy(&ctx_r));
    CeedCall(CeedQFunctionAddInput(qf_restrict, "input", num_comp, CEED_EVAL_NONE));
    CeedCall(CeedQFunctionAddInput(qf_restrict, "scale", num_comp, CEED_EVAL_NONE));
    CeedCall(CeedQFunctionAddOutput(qf_restrict, "output", num_comp, CEED_EVAL_INTERP));
    CeedCall(CeedQFunctionSetUserFlopsEstimate(qf_restrict, num_comp));

    CeedCall(CeedOperatorCreate(ceed, qf_restrict, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, op_restrict));
    CeedCall(CeedOperatorSetField(*op_restrict, "input", rstr_fine, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
    CeedCall(CeedOperatorSetField(*op_restrict, "scale", rstr_p_mult_fine, CEED_BASIS_NONE, mult_vec));
    CeedCall(CeedOperatorSetField(*op_restrict, "output", rstr_coarse, basis_c_to_f, CEED_VECTOR_ACTIVE));

    // Set name
    char *restriction_name;

    CeedCall(CeedCalloc(17 + name_len, &restriction_name));
    sprintf(restriction_name, "restriction%s%s", has_name ? " for " : "", has_name ? op_fine->name : "");
    CeedCall(CeedOperatorSetName(*op_restrict, restriction_name));
    CeedCall(CeedFree(&restriction_name));

    // Check
    CeedCall(CeedOperatorCheckReady(*op_restrict));

    // Cleanup
    CeedCall(CeedQFunctionDestroy(&qf_restrict));
  }

  // Prolongation
  if (op_prolong) {
    CeedInt             *num_comp_p_data;
    CeedQFunctionContext ctx_p;
    CeedQFunction        qf_prolong;

    CeedCall(CeedQFunctionCreateInteriorByName(ceed, "Scale", &qf_prolong));
    CeedCall(CeedCalloc(1, &num_comp_p_data));
    num_comp_p_data[0] = num_comp;
    CeedCall(CeedQFunctionContextCreate(ceed, &ctx_p));
    CeedCall(CeedQFunctionContextSetData(ctx_p, CEED_MEM_HOST, CEED_OWN_POINTER, sizeof(*num_comp_p_data), num_comp_p_data));
    CeedCall(CeedQFunctionSetContext(qf_prolong, ctx_p));
    CeedCall(CeedQFunctionContextDestroy(&ctx_p));
    CeedCall(CeedQFunctionAddInput(qf_prolong, "input", num_comp, CEED_EVAL_INTERP));
    CeedCall(CeedQFunctionAddInput(qf_prolong, "scale", num_comp, CEED_EVAL_NONE));
    CeedCall(CeedQFunctionAddOutput(qf_prolong, "output", num_comp, CEED_EVAL_NONE));
    CeedCall(CeedQFunctionSetUserFlopsEstimate(qf_prolong, num_comp));

    CeedCall(CeedOperatorCreate(ceed, qf_prolong, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, op_prolong));
    CeedCall(CeedOperatorSetField(*op_prolong, "input", rstr_coarse, basis_c_to_f, CEED_VECTOR_ACTIVE));
    CeedCall(CeedOperatorSetField(*op_prolong, "scale", rstr_p_mult_fine, CEED_BASIS_NONE, mult_vec));
    CeedCall(CeedOperatorSetField(*op_prolong, "output", rstr_fine, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));

    // Set name
    char *prolongation_name;

    CeedCall(CeedCalloc(18 + name_len, &prolongation_name));
    sprintf(prolongation_name, "prolongation%s%s", has_name ? " for " : "", has_name ? op_fine->name : "");
    CeedCall(CeedOperatorSetName(*op_prolong, prolongation_name));
    CeedCall(CeedFree(&prolongation_name));

    // Check
    CeedCall(CeedOperatorCheckReady(*op_prolong));

    // Cleanup
    CeedCall(CeedQFunctionDestroy(&qf_prolong));
  }

  // Check
  CeedCall(CeedOperatorCheckReady(*op_coarse));

  // Cleanup
  CeedCall(CeedVectorDestroy(&mult_vec));
  CeedCall(CeedElemRestrictionDestroy(&rstr_fine));
  CeedCall(CeedElemRestrictionDestroy(&rstr_p_mult_fine));
  CeedCall(CeedBasisDestroy(&basis_c_to_f));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Build 1D mass matrix and Laplacian with perturbation

  @param[in]  interp_1d   Interpolation matrix in one dimension
  @param[in]  grad_1d     Gradient matrix in one dimension
  @param[in]  q_weight_1d Quadrature weights in one dimension
  @param[in]  P_1d        Number of basis nodes in one dimension
  @param[in]  Q_1d        Number of quadrature points in one dimension
  @param[in]  dim         Dimension of basis
  @param[out] mass        Assembled mass matrix in one dimension
  @param[out] laplace     Assembled perturbed Laplacian in one dimension

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

  @ref Developer
**/
CeedPragmaOptimizeOff
static int CeedBuildMassLaplace(const CeedScalar *interp_1d, const CeedScalar *grad_1d, const CeedScalar *q_weight_1d, CeedInt P_1d, CeedInt Q_1d,
                                CeedInt dim, CeedScalar *mass, CeedScalar *laplace) {
  for (CeedInt i = 0; i < P_1d; i++) {
    for (CeedInt j = 0; j < P_1d; j++) {
      CeedScalar sum = 0.0;
      for (CeedInt k = 0; k < Q_1d; k++) sum += interp_1d[k * P_1d + i] * q_weight_1d[k] * interp_1d[k * P_1d + j];
      mass[i + j * P_1d] = sum;
    }
  }
  // -- Laplacian
  for (CeedInt i = 0; i < P_1d; i++) {
    for (CeedInt j = 0; j < P_1d; j++) {
      CeedScalar sum = 0.0;

      for (CeedInt k = 0; k < Q_1d; k++) sum += grad_1d[k * P_1d + i] * q_weight_1d[k] * grad_1d[k * P_1d + j];
      laplace[i + j * P_1d] = sum;
    }
  }
  CeedScalar perturbation = dim > 2 ? 1e-6 : 1e-4;
  for (CeedInt i = 0; i < P_1d; i++) laplace[i + P_1d * i] += perturbation;
  return CEED_ERROR_SUCCESS;
}
CeedPragmaOptimizeOn

/// @}

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

/**
  @brief Select correct basis matrix pointer based on @ref CeedEvalMode

  @param[in]  basis     `CeedBasis` from which to get the basis matrix
  @param[in]  eval_mode Current basis evaluation mode
  @param[in]  identity  Pointer to identity matrix
  @param[out] basis_ptr `CeedBasis` pointer to set

  @ref Backend
**/
int CeedOperatorGetBasisPointer(CeedBasis basis, CeedEvalMode eval_mode, const CeedScalar *identity, const CeedScalar **basis_ptr) {
  switch (eval_mode) {
    case CEED_EVAL_NONE:
      *basis_ptr = identity;
      break;
    case CEED_EVAL_INTERP:
      CeedCall(CeedBasisGetInterp(basis, basis_ptr));
      break;
    case CEED_EVAL_GRAD:
      CeedCall(CeedBasisGetGrad(basis, basis_ptr));
      break;
    case CEED_EVAL_DIV:
      CeedCall(CeedBasisGetDiv(basis, basis_ptr));
      break;
    case CEED_EVAL_CURL:
      CeedCall(CeedBasisGetCurl(basis, basis_ptr));
      break;
    case CEED_EVAL_WEIGHT:
      break;  // Caught by QF Assembly
  }
  assert(*basis_ptr != NULL);
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Create point block restriction for active `CeedOperatorField`

  @param[in]  rstr             Original `CeedElemRestriction` for active field
  @param[out] point_block_rstr Address of the variable where the newly created `CeedElemRestriction` will be stored

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

  @ref Backend
**/
int CeedOperatorCreateActivePointBlockRestriction(CeedElemRestriction rstr, CeedElemRestriction *point_block_rstr) {
  Ceed           ceed;
  CeedInt        num_elem, num_comp, shift, elem_size, comp_stride, *point_block_offsets;
  CeedSize       l_size;
  const CeedInt *offsets;

  CeedCall(CeedElemRestrictionGetCeed(rstr, &ceed));
  CeedCall(CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets));

  // Expand offsets
  CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem));
  CeedCall(CeedElemRestrictionGetNumComponents(rstr, &num_comp));
  CeedCall(CeedElemRestrictionGetElementSize(rstr, &elem_size));
  CeedCall(CeedElemRestrictionGetCompStride(rstr, &comp_stride));
  CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &l_size));
  shift = num_comp;
  if (comp_stride != 1) shift *= num_comp;
  CeedCall(CeedCalloc(num_elem * elem_size, &point_block_offsets));
  for (CeedInt i = 0; i < num_elem * elem_size; i++) {
    point_block_offsets[i] = offsets[i] * shift;
  }

  // Create new restriction
  CeedCall(CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp * num_comp, 1, l_size * num_comp, CEED_MEM_HOST, CEED_OWN_POINTER,
                                     point_block_offsets, point_block_rstr));

  // Cleanup
  CeedCall(CeedElemRestrictionRestoreOffsets(rstr, &offsets));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get `CeedQFunctionAssemblyData`

  @param[in]  op   `CeedOperator` to assemble
  @param[out] data `CeedQFunctionAssemblyData`

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

  @ref Backend
**/
int CeedOperatorGetQFunctionAssemblyData(CeedOperator op, CeedQFunctionAssemblyData *data) {
  if (!op->qf_assembled) {
    CeedQFunctionAssemblyData data;

    CeedCall(CeedQFunctionAssemblyDataCreate(op->ceed, &data));
    op->qf_assembled = data;
  }
  *data = op->qf_assembled;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Create object holding `CeedQFunction` assembly data for `CeedOperator`

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

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

  @ref Backend
**/
int CeedQFunctionAssemblyDataCreate(Ceed ceed, CeedQFunctionAssemblyData *data) {
  CeedCall(CeedCalloc(1, data));
  (*data)->ref_count = 1;
  (*data)->ceed      = ceed;
  CeedCall(CeedReference(ceed));
  return CEED_ERROR_SUCCESS;
}

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

  @param[in,out] data `CeedQFunctionAssemblyData` to increment the reference counter

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

  @ref Backend
**/
int CeedQFunctionAssemblyDataReference(CeedQFunctionAssemblyData data) {
  data->ref_count++;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Set re-use of `CeedQFunctionAssemblyData`

  @param[in,out] data       `CeedQFunctionAssemblyData` to mark for reuse
  @param[in]     reuse_data Boolean flag indicating data re-use

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

  @ref Backend
**/
int CeedQFunctionAssemblyDataSetReuse(CeedQFunctionAssemblyData data, bool reuse_data) {
  data->reuse_data        = reuse_data;
  data->needs_data_update = true;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Mark `CeedQFunctionAssemblyData` as stale

  @param[in,out] data              `CeedQFunctionAssemblyData` to mark as stale
  @param[in]     needs_data_update Boolean flag indicating if update is needed or completed

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

  @ref Backend
**/
int CeedQFunctionAssemblyDataSetUpdateNeeded(CeedQFunctionAssemblyData data, bool needs_data_update) {
  data->needs_data_update = needs_data_update;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Determine if `CeedQFunctionAssemblyData` needs update

  @param[in]  data             `CeedQFunctionAssemblyData` to mark as stale
  @param[out] is_update_needed Boolean flag indicating if re-assembly is required

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

  @ref Backend
**/
int CeedQFunctionAssemblyDataIsUpdateNeeded(CeedQFunctionAssemblyData data, bool *is_update_needed) {
  *is_update_needed = !data->reuse_data || data->needs_data_update;
  return CEED_ERROR_SUCCESS;
}

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

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

  Note: If the value of ` *data_copy` passed to this function is non-`NULL` , then it is assumed that ` *data_copy` is a pointer to a `CeedQFunctionAssemblyData`.
        This `CeedQFunctionAssemblyData` will be destroyed if ` *data_copy` is the only reference to this `CeedQFunctionAssemblyData`.

  @param[in]     data      `CeedQFunctionAssemblyData` to copy reference to
  @param[in,out] data_copy Variable to store copied reference

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

  @ref Backend
**/
int CeedQFunctionAssemblyDataReferenceCopy(CeedQFunctionAssemblyData data, CeedQFunctionAssemblyData *data_copy) {
  CeedCall(CeedQFunctionAssemblyDataReference(data));
  CeedCall(CeedQFunctionAssemblyDataDestroy(data_copy));
  *data_copy = data;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get setup status for internal objects for `CeedQFunctionAssemblyData`

  @param[in]  data     `CeedQFunctionAssemblyData` to retrieve status
  @param[out] is_setup Boolean flag for setup status

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

  @ref Backend
**/
int CeedQFunctionAssemblyDataIsSetup(CeedQFunctionAssemblyData data, bool *is_setup) {
  *is_setup = data->is_setup;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Set internal objects for `CeedQFunctionAssemblyData`

  @param[in,out] data `CeedQFunctionAssemblyData` to set objects
  @param[in]     vec  `CeedVector` to store assembled `CeedQFunction` at quadrature points
  @param[in]     rstr `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction`

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

  @ref Backend
**/
int CeedQFunctionAssemblyDataSetObjects(CeedQFunctionAssemblyData data, CeedVector vec, CeedElemRestriction rstr) {
  CeedCall(CeedVectorReferenceCopy(vec, &data->vec));
  CeedCall(CeedElemRestrictionReferenceCopy(rstr, &data->rstr));

  data->is_setup = true;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get internal objects for `CeedQFunctionAssemblyData`

  @param[in,out] data `CeedQFunctionAssemblyData` to set objects
  @param[out]    vec  `CeedVector` to store assembled `CeedQFunction` at quadrature points
  @param[out]    rstr `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction`

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

  @ref Backend
**/
int CeedQFunctionAssemblyDataGetObjects(CeedQFunctionAssemblyData data, CeedVector *vec, CeedElemRestriction *rstr) {
  CeedCheck(data->is_setup, data->ceed, CEED_ERROR_INCOMPLETE, "Internal objects not set; must call CeedQFunctionAssemblyDataSetObjects first.");

  CeedCall(CeedVectorReferenceCopy(data->vec, vec));
  CeedCall(CeedElemRestrictionReferenceCopy(data->rstr, rstr));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Destroy `CeedQFunctionAssemblyData`

  @param[in,out] data  `CeedQFunctionAssemblyData` to destroy

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

  @ref Backend
**/
int CeedQFunctionAssemblyDataDestroy(CeedQFunctionAssemblyData *data) {
  if (!*data || --(*data)->ref_count > 0) {
    *data = NULL;
    return CEED_ERROR_SUCCESS;
  }
  CeedCall(CeedDestroy(&(*data)->ceed));
  CeedCall(CeedVectorDestroy(&(*data)->vec));
  CeedCall(CeedElemRestrictionDestroy(&(*data)->rstr));

  CeedCall(CeedFree(data));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get `CeedOperatorAssemblyData`

  @param[in]  op   `CeedOperator` to assemble
  @param[out] data `CeedOperatorAssemblyData`

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

  @ref Backend
**/
int CeedOperatorGetOperatorAssemblyData(CeedOperator op, CeedOperatorAssemblyData *data) {
  if (!op->op_assembled) {
    CeedOperatorAssemblyData data;

    CeedCall(CeedOperatorAssemblyDataCreate(op->ceed, op, &data));
    op->op_assembled = data;
  }
  *data = op->op_assembled;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Create object holding `CeedOperator` assembly data.

  The `CeedOperatorAssemblyData` holds an array with references to every active `CeedBasis` used in the `CeedOperator`.
  An array with references to the corresponding active `CeedElemRestriction` is also stored.
  For each active `CeedBasis, the `CeedOperatorAssemblyData` holds an array of all input and output @ref CeedEvalMode for this `CeedBasis`.
  The `CeedOperatorAssemblyData` holds an array of offsets for indexing into the assembled `CeedQFunction` arrays to the row representing each @ref CeedEvalMode.
  The number of input columns across all active bases for the assembled `CeedQFunction` is also stored.
  Lastly, the `CeedOperatorAssembly` data holds assembled matrices representing the full action of the `CeedBasis` for all @ref CeedEvalMode.

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

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

  @ref Backend
**/
int CeedOperatorAssemblyDataCreate(Ceed ceed, CeedOperator op, CeedOperatorAssemblyData *data) {
  CeedInt             num_active_bases_in = 0, num_active_bases_out = 0, offset = 0;
  CeedInt             num_input_fields, *num_eval_modes_in = NULL, num_output_fields, *num_eval_modes_out = NULL;
  CeedSize          **eval_mode_offsets_in = NULL, **eval_mode_offsets_out = NULL;
  CeedEvalMode      **eval_modes_in = NULL, **eval_modes_out = NULL;
  CeedQFunctionField *qf_fields;
  CeedQFunction       qf;
  CeedOperatorField  *op_fields;
  bool                is_composite;

  CeedCall(CeedOperatorIsComposite(op, &is_composite));
  CeedCheck(!is_composite, ceed, CEED_ERROR_INCOMPATIBLE, "Can only create CeedOperator assembly data for non-composite operators.");

  // Allocate
  CeedCall(CeedCalloc(1, data));
  (*data)->ceed = ceed;
  CeedCall(CeedReference(ceed));

  // Build OperatorAssembly data
  CeedCall(CeedOperatorGetQFunction(op, &qf));

  // Determine active input basis
  CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &qf_fields, NULL, NULL));
  CeedCall(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL));
  for (CeedInt i = 0; i < num_input_fields; i++) {
    CeedVector vec;

    CeedCall(CeedOperatorFieldGetVector(op_fields[i], &vec));
    if (vec == CEED_VECTOR_ACTIVE) {
      CeedInt      index = -1, num_comp, q_comp;
      CeedEvalMode eval_mode;
      CeedBasis    basis_in = NULL;

      CeedCall(CeedOperatorFieldGetBasis(op_fields[i], &basis_in));
      CeedCall(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode));
      CeedCall(CeedBasisGetNumComponents(basis_in, &num_comp));
      CeedCall(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp));
      for (CeedInt i = 0; i < num_active_bases_in; i++) {
        if ((*data)->active_bases_in[i] == basis_in) index = i;
      }
      if (index == -1) {
        CeedElemRestriction elem_rstr_in;

        index = num_active_bases_in;
        CeedCall(CeedRealloc(num_active_bases_in + 1, &(*data)->active_bases_in));
        (*data)->active_bases_in[num_active_bases_in] = NULL;
        CeedCall(CeedBasisReferenceCopy(basis_in, &(*data)->active_bases_in[num_active_bases_in]));
        CeedCall(CeedRealloc(num_active_bases_in + 1, &(*data)->active_elem_rstrs_in));
        (*data)->active_elem_rstrs_in[num_active_bases_in] = NULL;
        CeedCall(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr_in));
        CeedCall(CeedElemRestrictionReferenceCopy(elem_rstr_in, &(*data)->active_elem_rstrs_in[num_active_bases_in]));
        CeedCall(CeedElemRestrictionDestroy(&elem_rstr_in));
        CeedCall(CeedRealloc(num_active_bases_in + 1, &num_eval_modes_in));
        num_eval_modes_in[index] = 0;
        CeedCall(CeedRealloc(num_active_bases_in + 1, &eval_modes_in));
        eval_modes_in[index] = NULL;
        CeedCall(CeedRealloc(num_active_bases_in + 1, &eval_mode_offsets_in));
        eval_mode_offsets_in[index] = NULL;
        CeedCall(CeedRealloc(num_active_bases_in + 1, &(*data)->assembled_bases_in));
        (*data)->assembled_bases_in[index] = NULL;
        num_active_bases_in++;
      }
      if (eval_mode != CEED_EVAL_WEIGHT) {
        // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly
        CeedCall(CeedRealloc(num_eval_modes_in[index] + q_comp, &eval_modes_in[index]));
        CeedCall(CeedRealloc(num_eval_modes_in[index] + q_comp, &eval_mode_offsets_in[index]));
        for (CeedInt d = 0; d < q_comp; d++) {
          eval_modes_in[index][num_eval_modes_in[index] + d]        = eval_mode;
          eval_mode_offsets_in[index][num_eval_modes_in[index] + d] = offset;
          offset += num_comp;
        }
        num_eval_modes_in[index] += q_comp;
      }
      CeedCall(CeedBasisDestroy(&basis_in));
    }
    CeedCall(CeedVectorDestroy(&vec));
  }

  // Determine active output basis
  CeedCall(CeedQFunctionGetFields(qf, NULL, NULL, &num_output_fields, &qf_fields));
  CeedCall(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields));
  offset = 0;
  for (CeedInt i = 0; i < num_output_fields; i++) {
    CeedVector vec;

    CeedCall(CeedOperatorFieldGetVector(op_fields[i], &vec));
    if (vec == CEED_VECTOR_ACTIVE) {
      CeedInt      index = -1, num_comp, q_comp;
      CeedEvalMode eval_mode;
      CeedBasis    basis_out = NULL;

      CeedCall(CeedOperatorFieldGetBasis(op_fields[i], &basis_out));
      CeedCall(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode));
      CeedCall(CeedBasisGetNumComponents(basis_out, &num_comp));
      CeedCall(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp));
      for (CeedInt i = 0; i < num_active_bases_out; i++) {
        if ((*data)->active_bases_out[i] == basis_out) index = i;
      }
      if (index == -1) {
        CeedElemRestriction elem_rstr_out;

        index = num_active_bases_out;
        CeedCall(CeedRealloc(num_active_bases_out + 1, &(*data)->active_bases_out));
        (*data)->active_bases_out[num_active_bases_out] = NULL;
        CeedCall(CeedBasisReferenceCopy(basis_out, &(*data)->active_bases_out[num_active_bases_out]));
        CeedCall(CeedRealloc(num_active_bases_out + 1, &(*data)->active_elem_rstrs_out));
        (*data)->active_elem_rstrs_out[num_active_bases_out] = NULL;
        CeedCall(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr_out));
        CeedCall(CeedElemRestrictionReferenceCopy(elem_rstr_out, &(*data)->active_elem_rstrs_out[num_active_bases_out]));
        CeedCall(CeedElemRestrictionDestroy(&elem_rstr_out));
        CeedCall(CeedRealloc(num_active_bases_out + 1, &num_eval_modes_out));
        num_eval_modes_out[index] = 0;
        CeedCall(CeedRealloc(num_active_bases_out + 1, &eval_modes_out));
        eval_modes_out[index] = NULL;
        CeedCall(CeedRealloc(num_active_bases_out + 1, &eval_mode_offsets_out));
        eval_mode_offsets_out[index] = NULL;
        CeedCall(CeedRealloc(num_active_bases_out + 1, &(*data)->assembled_bases_out));
        (*data)->assembled_bases_out[index] = NULL;
        num_active_bases_out++;
      }
      if (eval_mode != CEED_EVAL_WEIGHT) {
        // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly
        CeedCall(CeedRealloc(num_eval_modes_out[index] + q_comp, &eval_modes_out[index]));
        CeedCall(CeedRealloc(num_eval_modes_out[index] + q_comp, &eval_mode_offsets_out[index]));
        for (CeedInt d = 0; d < q_comp; d++) {
          eval_modes_out[index][num_eval_modes_out[index] + d]        = eval_mode;
          eval_mode_offsets_out[index][num_eval_modes_out[index] + d] = offset;
          offset += num_comp;
        }
        num_eval_modes_out[index] += q_comp;
      }
      CeedCall(CeedBasisDestroy(&basis_out));
    }
    CeedCall(CeedVectorDestroy(&vec));
  }
  (*data)->num_active_bases_in   = num_active_bases_in;
  (*data)->num_eval_modes_in     = num_eval_modes_in;
  (*data)->eval_modes_in         = eval_modes_in;
  (*data)->eval_mode_offsets_in  = eval_mode_offsets_in;
  (*data)->num_active_bases_out  = num_active_bases_out;
  (*data)->num_eval_modes_out    = num_eval_modes_out;
  (*data)->eval_modes_out        = eval_modes_out;
  (*data)->eval_mode_offsets_out = eval_mode_offsets_out;
  (*data)->num_output_components = offset;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get `CeedOperator` @ref CeedEvalMode for assembly.

  Note: See @ref CeedOperatorAssemblyDataCreate() for a full description of the data stored in this object.

  @param[in]  data                  `CeedOperatorAssemblyData`
  @param[out] num_active_bases_in   Total number of active bases for input
  @param[out] num_eval_modes_in     Pointer to hold array of numbers of input @ref CeedEvalMode, or `NULL`.
                                      `eval_modes_in[0]` holds an array of eval modes for the first active `CeedBasis`.
  @param[out] eval_modes_in         Pointer to hold arrays of input @ref CeedEvalMode, or `NULL`
  @param[out] eval_mode_offsets_in  Pointer to hold arrays of input offsets at each quadrature point
  @param[out] num_active_bases_out  Total number of active bases for output
  @param[out] num_eval_modes_out    Pointer to hold array of numbers of output @ref CeedEvalMode, or `NULL`
  @param[out] eval_modes_out        Pointer to hold arrays of output @ref CeedEvalMode, or `NULL`
  @param[out] eval_mode_offsets_out Pointer to hold arrays of output offsets at each quadrature point
  @param[out] num_output_components The number of columns in the assembled `CeedQFunction` matrix for each quadrature point, including contributions of all active bases

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

  @ref Backend
**/
int CeedOperatorAssemblyDataGetEvalModes(CeedOperatorAssemblyData data, CeedInt *num_active_bases_in, CeedInt **num_eval_modes_in,
                                         const CeedEvalMode ***eval_modes_in, CeedSize ***eval_mode_offsets_in, CeedInt *num_active_bases_out,
                                         CeedInt **num_eval_modes_out, const CeedEvalMode ***eval_modes_out, CeedSize ***eval_mode_offsets_out,
                                         CeedSize *num_output_components) {
  if (num_active_bases_in) *num_active_bases_in = data->num_active_bases_in;
  if (num_eval_modes_in) *num_eval_modes_in = data->num_eval_modes_in;
  if (eval_modes_in) *eval_modes_in = (const CeedEvalMode **)data->eval_modes_in;
  if (eval_mode_offsets_in) *eval_mode_offsets_in = data->eval_mode_offsets_in;
  if (num_active_bases_out) *num_active_bases_out = data->num_active_bases_out;
  if (num_eval_modes_out) *num_eval_modes_out = data->num_eval_modes_out;
  if (eval_modes_out) *eval_modes_out = (const CeedEvalMode **)data->eval_modes_out;
  if (eval_mode_offsets_out) *eval_mode_offsets_out = data->eval_mode_offsets_out;
  if (num_output_components) *num_output_components = data->num_output_components;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get `CeedOperator` `CeedBasis` data for assembly.

  Note: See @ref CeedOperatorAssemblyDataCreate() for a full description of the data stored in this object.

  @param[in]  data                 `CeedOperatorAssemblyData`
  @param[out] num_active_bases_in  Number of active input bases, or `NULL`
  @param[out] active_bases_in      Pointer to hold active input `CeedBasis`, or `NULL`
  @param[out] assembled_bases_in   Pointer to hold assembled active input `B` , or `NULL`
  @param[out] num_active_bases_out Number of active output bases, or `NULL`
  @param[out] active_bases_out     Pointer to hold active output `CeedBasis`, or `NULL`
  @param[out] assembled_bases_out  Pointer to hold assembled active output `B` , or `NULL`

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

  @ref Backend
**/
int CeedOperatorAssemblyDataGetBases(CeedOperatorAssemblyData data, CeedInt *num_active_bases_in, CeedBasis **active_bases_in,
                                     const CeedScalar ***assembled_bases_in, CeedInt *num_active_bases_out, CeedBasis **active_bases_out,
                                     const CeedScalar ***assembled_bases_out) {
  // Assemble B_in, B_out if needed
  if (assembled_bases_in && !data->assembled_bases_in[0]) {
    CeedInt num_qpts;

    if (data->active_bases_in[0] == CEED_BASIS_NONE) CeedCall(CeedElemRestrictionGetElementSize(data->active_elem_rstrs_in[0], &num_qpts));
    else CeedCall(CeedBasisGetNumQuadraturePoints(data->active_bases_in[0], &num_qpts));
    for (CeedInt b = 0; b < data->num_active_bases_in; b++) {
      bool        has_eval_none = false;
      CeedInt     num_nodes;
      CeedScalar *B_in = NULL, *identity = NULL;

      CeedCall(CeedElemRestrictionGetElementSize(data->active_elem_rstrs_in[b], &num_nodes));
      CeedCall(CeedCalloc(num_qpts * num_nodes * data->num_eval_modes_in[b], &B_in));

      for (CeedInt i = 0; i < data->num_eval_modes_in[b]; i++) {
        has_eval_none = has_eval_none || (data->eval_modes_in[b][i] == CEED_EVAL_NONE);
      }
      if (has_eval_none) {
        CeedCall(CeedCalloc(num_qpts * num_nodes, &identity));
        for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) {
          identity[i * num_nodes + i] = 1.0;
        }
      }

      for (CeedInt q = 0; q < num_qpts; q++) {
        for (CeedInt n = 0; n < num_nodes; n++) {
          CeedInt      d_in              = 0, q_comp_in;
          CeedEvalMode eval_mode_in_prev = CEED_EVAL_NONE;

          for (CeedInt e_in = 0; e_in < data->num_eval_modes_in[b]; e_in++) {
            const CeedInt     qq = data->num_eval_modes_in[b] * q;
            const CeedScalar *B  = NULL;

            CeedCall(CeedOperatorGetBasisPointer(data->active_bases_in[b], data->eval_modes_in[b][e_in], identity, &B));
            CeedCall(CeedBasisGetNumQuadratureComponents(data->active_bases_in[b], data->eval_modes_in[b][e_in], &q_comp_in));
            if (q_comp_in > 1) {
              if (e_in == 0 || data->eval_modes_in[b][e_in] != eval_mode_in_prev) d_in = 0;
              else B = &B[(++d_in) * num_qpts * num_nodes];
            }
            eval_mode_in_prev                 = data->eval_modes_in[b][e_in];
            B_in[(qq + e_in) * num_nodes + n] = B[q * num_nodes + n];
          }
        }
      }
      if (identity) CeedCall(CeedFree(&identity));
      data->assembled_bases_in[b] = B_in;
    }
  }

  if (assembled_bases_out && !data->assembled_bases_out[0]) {
    CeedInt num_qpts;

    if (data->active_bases_out[0] == CEED_BASIS_NONE) CeedCall(CeedElemRestrictionGetElementSize(data->active_elem_rstrs_out[0], &num_qpts));
    else CeedCall(CeedBasisGetNumQuadraturePoints(data->active_bases_out[0], &num_qpts));
    for (CeedInt b = 0; b < data->num_active_bases_out; b++) {
      bool        has_eval_none = false;
      CeedInt     num_nodes;
      CeedScalar *B_out = NULL, *identity = NULL;

      CeedCall(CeedElemRestrictionGetElementSize(data->active_elem_rstrs_out[b], &num_nodes));
      CeedCall(CeedCalloc(num_qpts * num_nodes * data->num_eval_modes_out[b], &B_out));

      for (CeedInt i = 0; i < data->num_eval_modes_out[b]; i++) {
        has_eval_none = has_eval_none || (data->eval_modes_out[b][i] == CEED_EVAL_NONE);
      }
      if (has_eval_none) {
        CeedCall(CeedCalloc(num_qpts * num_nodes, &identity));
        for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) {
          identity[i * num_nodes + i] = 1.0;
        }
      }

      for (CeedInt q = 0; q < num_qpts; q++) {
        for (CeedInt n = 0; n < num_nodes; n++) {
          CeedInt      d_out              = 0, q_comp_out;
          CeedEvalMode eval_mode_out_prev = CEED_EVAL_NONE;

          for (CeedInt e_out = 0; e_out < data->num_eval_modes_out[b]; e_out++) {
            const CeedInt     qq = data->num_eval_modes_out[b] * q;
            const CeedScalar *B  = NULL;

            CeedCall(CeedOperatorGetBasisPointer(data->active_bases_out[b], data->eval_modes_out[b][e_out], identity, &B));
            CeedCall(CeedBasisGetNumQuadratureComponents(data->active_bases_out[b], data->eval_modes_out[b][e_out], &q_comp_out));
            if (q_comp_out > 1) {
              if (e_out == 0 || data->eval_modes_out[b][e_out] != eval_mode_out_prev) d_out = 0;
              else B = &B[(++d_out) * num_qpts * num_nodes];
            }
            eval_mode_out_prev                  = data->eval_modes_out[b][e_out];
            B_out[(qq + e_out) * num_nodes + n] = B[q * num_nodes + n];
          }
        }
      }
      if (identity) CeedCall(CeedFree(&identity));
      data->assembled_bases_out[b] = B_out;
    }
  }

  // Pass out assembled data
  if (num_active_bases_in) *num_active_bases_in = data->num_active_bases_in;
  if (active_bases_in) *active_bases_in = data->active_bases_in;
  if (assembled_bases_in) *assembled_bases_in = (const CeedScalar **)data->assembled_bases_in;
  if (num_active_bases_out) *num_active_bases_out = data->num_active_bases_out;
  if (active_bases_out) *active_bases_out = data->active_bases_out;
  if (assembled_bases_out) *assembled_bases_out = (const CeedScalar **)data->assembled_bases_out;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get `CeedOperator` `CeedBasis` data for assembly.

  Note: See @ref CeedOperatorAssemblyDataCreate() for a full description of the data stored in this object.

  @param[in]  data                      `CeedOperatorAssemblyData`
  @param[out] num_active_elem_rstrs_in  Number of active input element restrictions, or `NULL`
  @param[out] active_elem_rstrs_in      Pointer to hold active input `CeedElemRestriction`, or `NULL`
  @param[out] num_active_elem_rstrs_out Number of active output element restrictions, or `NULL`
  @param[out] active_elem_rstrs_out     Pointer to hold active output `CeedElemRestriction`, or `NULL`

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

  @ref Backend
**/
int CeedOperatorAssemblyDataGetElemRestrictions(CeedOperatorAssemblyData data, CeedInt *num_active_elem_rstrs_in,
                                                CeedElemRestriction **active_elem_rstrs_in, CeedInt *num_active_elem_rstrs_out,
                                                CeedElemRestriction **active_elem_rstrs_out) {
  if (num_active_elem_rstrs_in) *num_active_elem_rstrs_in = data->num_active_bases_in;
  if (active_elem_rstrs_in) *active_elem_rstrs_in = data->active_elem_rstrs_in;
  if (num_active_elem_rstrs_out) *num_active_elem_rstrs_out = data->num_active_bases_out;
  if (active_elem_rstrs_out) *active_elem_rstrs_out = data->active_elem_rstrs_out;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Destroy `CeedOperatorAssemblyData`

  @param[in,out] data `CeedOperatorAssemblyData` to destroy

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

  @ref Backend
**/
int CeedOperatorAssemblyDataDestroy(CeedOperatorAssemblyData *data) {
  if (!*data) {
    *data = NULL;
    return CEED_ERROR_SUCCESS;
  }
  CeedCall(CeedDestroy(&(*data)->ceed));
  for (CeedInt b = 0; b < (*data)->num_active_bases_in; b++) {
    CeedCall(CeedBasisDestroy(&(*data)->active_bases_in[b]));
    CeedCall(CeedElemRestrictionDestroy(&(*data)->active_elem_rstrs_in[b]));
    CeedCall(CeedFree(&(*data)->eval_modes_in[b]));
    CeedCall(CeedFree(&(*data)->eval_mode_offsets_in[b]));
    CeedCall(CeedFree(&(*data)->assembled_bases_in[b]));
  }
  for (CeedInt b = 0; b < (*data)->num_active_bases_out; b++) {
    CeedCall(CeedBasisDestroy(&(*data)->active_bases_out[b]));
    CeedCall(CeedElemRestrictionDestroy(&(*data)->active_elem_rstrs_out[b]));
    CeedCall(CeedFree(&(*data)->eval_modes_out[b]));
    CeedCall(CeedFree(&(*data)->eval_mode_offsets_out[b]));
    CeedCall(CeedFree(&(*data)->assembled_bases_out[b]));
  }
  CeedCall(CeedFree(&(*data)->active_bases_in));
  CeedCall(CeedFree(&(*data)->active_bases_out));
  CeedCall(CeedFree(&(*data)->active_elem_rstrs_in));
  CeedCall(CeedFree(&(*data)->active_elem_rstrs_out));
  CeedCall(CeedFree(&(*data)->num_eval_modes_in));
  CeedCall(CeedFree(&(*data)->num_eval_modes_out));
  CeedCall(CeedFree(&(*data)->eval_modes_in));
  CeedCall(CeedFree(&(*data)->eval_modes_out));
  CeedCall(CeedFree(&(*data)->eval_mode_offsets_in));
  CeedCall(CeedFree(&(*data)->eval_mode_offsets_out));
  CeedCall(CeedFree(&(*data)->assembled_bases_in));
  CeedCall(CeedFree(&(*data)->assembled_bases_out));

  CeedCall(CeedFree(data));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Retrieve fallback `CeedOperator` with a reference `Ceed` for advanced `CeedOperator` functionality

  @param[in]  op          `CeedOperator` to retrieve fallback for
  @param[out] op_fallback Fallback `CeedOperator`

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

  @ref Backend
**/
int CeedOperatorGetFallback(CeedOperator op, CeedOperator *op_fallback) {
  // Create if needed
  if (!op->op_fallback) CeedCall(CeedOperatorCreateFallback(op));
  if (op->op_fallback) {
    bool is_debug;
    Ceed ceed;

    CeedCall(CeedOperatorGetCeed(op, &ceed));
    CeedCall(CeedIsDebug(ceed, &is_debug));
    if (is_debug) {
      Ceed        ceed_fallback;
      const char *resource, *resource_fallback;

      CeedCall(CeedGetOperatorFallbackCeed(ceed, &ceed_fallback));
      CeedCall(CeedGetResource(ceed, &resource));
      CeedCall(CeedGetResource(ceed_fallback, &resource_fallback));

      CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "---------- CeedOperator Fallback ----------\n");
      CeedDebug(ceed, "Falling back from %s operator at address %p to %s operator at address %p\n", resource, op, resource_fallback, op->op_fallback);
    }
  }
  *op_fallback = op->op_fallback;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the parent `CeedOperator` for a fallback `CeedOperator`

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

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

  @ref Backend
**/
int CeedOperatorGetFallbackParent(CeedOperator op, CeedOperator *parent) {
  *parent = op->op_fallback_parent ? op->op_fallback_parent : NULL;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the `Ceed` context of the parent `CeedOperator` for a fallback `CeedOperator`

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

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

  @ref Backend
**/
int CeedOperatorGetFallbackParentCeed(CeedOperator op, Ceed *parent) {
  *parent = op->op_fallback_parent ? op->op_fallback_parent->ceed : op->ceed;
  return CEED_ERROR_SUCCESS;
}

/// @}

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

/**
  @brief Assemble a linear `CeedQFunction` associated with a `CeedOperator`.

  This returns a `CeedVector` containing a matrix at each quadrature point providing the action of the `CeedQFunction` associated with the `CeedOperator`.
  The vector `assembled` is of shape `[num_elements, num_input_fields, num_output_fields, num_quad_points]` and contains column-major matrices representing the action of the `CeedQFunction` for a corresponding quadrature point on an element.

  Inputs and outputs are in the order provided by the user when adding `CeedOperator` fields.
  For example, a `CeedQFunction` with inputs `u` and `gradu` and outputs `gradv` and `v` , provided in that order, would result in an assembled `CeedQFunction` that consists of `(1 + dim) x (dim + 1)` matrices at each quadrature point acting on the input ` [u, du_0, du_1]` and producing the output `[dv_0, dv_1, v]`.

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

  @param[in]  op        `CeedOperator` to assemble `CeedQFunction`
  @param[out] assembled `CeedVector` to store assembled `CeedQFunction` at quadrature points
  @param[out] rstr      `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction`
  @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 CeedOperatorLinearAssembleQFunction(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) {
  CeedCall(CeedOperatorCheckReady(op));

  if (op->LinearAssembleQFunction) {
    // Backend version
    CeedCall(op->LinearAssembleQFunction(op, assembled, rstr, request));
  } else {
    // Operator fallback
    Ceed         ceed;
    CeedOperator op_fallback;

    CeedCall(CeedOperatorGetCeed(op, &ceed));
    CeedCall(CeedOperatorGetFallback(op, &op_fallback));
    if (op_fallback) CeedCall(CeedOperatorLinearAssembleQFunction(op_fallback, assembled, rstr, request));
    else return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support CeedOperatorLinearAssembleQFunction");
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Assemble `CeedQFunction` and store result internally.

  Return copied references of stored data to the caller.
  Caller is responsible for ownership and destruction of the copied references.
  See also @ref CeedOperatorLinearAssembleQFunction().

  Note: If the value of `assembled` or `rstr` passed to this function are non-`NULL` , then it is assumed that they hold valid pointers.
        These objects will be destroyed if `*assembled` or `*rstr` is the only reference to the object.

  @param[in]  op        `CeedOperator` to assemble `CeedQFunction`
  @param[out] assembled `CeedVector` to store assembled `CeedQFunction` at quadrature points
  @param[out] rstr      `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction`
  @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 CeedOperatorLinearAssembleQFunctionBuildOrUpdate(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) {
  int (*LinearAssembleQFunctionUpdate)(CeedOperator, CeedVector, CeedElemRestriction, CeedRequest *) = NULL;
  CeedOperator op_assemble                                                                           = NULL;
  CeedOperator op_fallback_parent                                                                    = NULL;

  CeedCall(CeedOperatorCheckReady(op));

  // Determine if fallback parent or operator has implementation
  CeedCall(CeedOperatorGetFallbackParent(op, &op_fallback_parent));
  if (op_fallback_parent && op_fallback_parent->LinearAssembleQFunctionUpdate) {
    // -- Backend version for op fallback parent is faster, if it exists
    LinearAssembleQFunctionUpdate = op_fallback_parent->LinearAssembleQFunctionUpdate;
    op_assemble                   = op_fallback_parent;
  } else if (op->LinearAssembleQFunctionUpdate) {
    // -- Backend version for op
    LinearAssembleQFunctionUpdate = op->LinearAssembleQFunctionUpdate;
    op_assemble                   = op;
  }

  // Assemble QFunction
  if (LinearAssembleQFunctionUpdate) {
    // Backend or fallback parent version
    CeedQFunctionAssemblyData data;
    bool                      data_is_setup;
    CeedVector                assembled_vec  = NULL;
    CeedElemRestriction       assembled_rstr = NULL;

    CeedCall(CeedOperatorGetQFunctionAssemblyData(op, &data));
    CeedCall(CeedQFunctionAssemblyDataIsSetup(data, &data_is_setup));
    if (data_is_setup) {
      bool update_needed;

      CeedCall(CeedQFunctionAssemblyDataGetObjects(data, &assembled_vec, &assembled_rstr));
      CeedCall(CeedQFunctionAssemblyDataIsUpdateNeeded(data, &update_needed));
      if (update_needed) CeedCall(LinearAssembleQFunctionUpdate(op_assemble, assembled_vec, assembled_rstr, request));
    } else {
      CeedCall(CeedOperatorLinearAssembleQFunction(op_assemble, &assembled_vec, &assembled_rstr, request));
      CeedCall(CeedQFunctionAssemblyDataSetObjects(data, assembled_vec, assembled_rstr));
    }
    CeedCall(CeedQFunctionAssemblyDataSetUpdateNeeded(data, false));

    // Copy reference from internally held copy
    CeedCall(CeedVectorReferenceCopy(assembled_vec, assembled));
    CeedCall(CeedElemRestrictionReferenceCopy(assembled_rstr, rstr));
    CeedCall(CeedVectorDestroy(&assembled_vec));
    CeedCall(CeedElemRestrictionDestroy(&assembled_rstr));
  } else {
    // Operator fallback
    Ceed         ceed;
    CeedOperator op_fallback;

    CeedCall(CeedOperatorGetCeed(op, &ceed));
    CeedCall(CeedOperatorGetFallback(op, &op_fallback));
    if (op_fallback) CeedCall(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op_fallback, assembled, rstr, request));
    else return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support CeedOperatorLinearAssembleQFunctionUpdate");
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Assemble the diagonal of a square linear `CeedOperator`

  This overwrites a `CeedVector` with the diagonal of a linear `CeedOperator`.

  Note: Currently only non-composite `CeedOperator` with a single field and composite `CeedOperator` with single field sub-operators are supported.

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

  @param[in]  op        `CeedOperator` to assemble `CeedQFunction`
  @param[out] assembled `CeedVector` to store assembled `CeedOperator` diagonal
  @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 CeedOperatorLinearAssembleDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) {
  bool     is_composite;
  CeedSize input_size = 0, output_size = 0;
  Ceed     ceed;

  CeedCall(CeedOperatorGetCeed(op, &ceed));
  CeedCall(CeedOperatorCheckReady(op));
  CeedCall(CeedOperatorIsComposite(op, &is_composite));

  CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size));
  CeedCheck(input_size == output_size, ceed, CEED_ERROR_DIMENSION, "Operator must be square");

  // Early exit for empty operator
  if (!is_composite) {
    CeedInt num_elem = 0;

    CeedCall(CeedOperatorGetNumElements(op, &num_elem));
    if (num_elem == 0) return CEED_ERROR_SUCCESS;
  }

  if (op->LinearAssembleDiagonal) {
    // Backend version
    CeedCall(op->LinearAssembleDiagonal(op, assembled, request));
    return CEED_ERROR_SUCCESS;
  } else if (op->LinearAssembleAddDiagonal) {
    // Backend version with zeroing first
    CeedCall(CeedVectorSetValue(assembled, 0.0));
    CeedCall(op->LinearAssembleAddDiagonal(op, assembled, request));
    return CEED_ERROR_SUCCESS;
  } else {
    // Operator fallback
    CeedOperator op_fallback;

    CeedCall(CeedOperatorGetFallback(op, &op_fallback));
    if (op_fallback) {
      CeedCall(CeedOperatorLinearAssembleDiagonal(op_fallback, assembled, request));
      return CEED_ERROR_SUCCESS;
    }
  }
  // Default interface implementation
  CeedCall(CeedVectorSetValue(assembled, 0.0));
  CeedCall(CeedOperatorLinearAssembleAddDiagonal(op, assembled, request));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Assemble the diagonal of a square linear `CeedOperator`.

  This sums into a `CeedVector` the diagonal of a linear `CeedOperator`.

  Note: Currently only non-composite `CeedOperator` with a single field and composite `CeedOperator` with single field sub-operators are supported.

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

  @param[in]  op        `CeedOperator` to assemble `CeedQFunction`
  @param[out] assembled `CeedVector` to store assembled `CeedOperator` diagonal
  @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 CeedOperatorLinearAssembleAddDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) {
  bool     is_composite;
  CeedSize input_size = 0, output_size = 0;
  Ceed     ceed;

  CeedCall(CeedOperatorGetCeed(op, &ceed));
  CeedCall(CeedOperatorCheckReady(op));
  CeedCall(CeedOperatorIsComposite(op, &is_composite));

  CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size));
  CeedCheck(input_size == output_size, ceed, CEED_ERROR_DIMENSION, "Operator must be square");

  // Early exit for empty operator
  if (!is_composite) {
    CeedInt num_elem = 0;

    CeedCall(CeedOperatorGetNumElements(op, &num_elem));
    if (num_elem == 0) return CEED_ERROR_SUCCESS;
  }

  if (op->LinearAssembleAddDiagonal) {
    // Backend version
    CeedCall(op->LinearAssembleAddDiagonal(op, assembled, request));
    return CEED_ERROR_SUCCESS;
  } else {
    // Operator fallback
    CeedOperator op_fallback;

    CeedCall(CeedOperatorGetFallback(op, &op_fallback));
    if (op_fallback) {
      CeedCall(CeedOperatorLinearAssembleAddDiagonal(op_fallback, assembled, request));
      return CEED_ERROR_SUCCESS;
    }
  }
  // Default interface implementation
  if (is_composite) {
    CeedCall(CeedCompositeOperatorLinearAssembleAddDiagonal(op, request, false, assembled));
  } else {
    CeedCall(CeedSingleOperatorLinearAssembleAddDiagonal(op, request, false, assembled));
  }
  return CEED_ERROR_SUCCESS;
}

/**
   @brief Fully assemble the point-block diagonal pattern of a linear `CeedOperator`.

   Expected to be used in conjunction with @ref CeedOperatorLinearAssemblePointBlockDiagonal().

   The assembly routines use coordinate format, with `num_entries` tuples of the form `(i, j, value)` which indicate that value should be added to the matrix in entry `(i, j)`.
   Note that the `(i, j)` pairs are unique.
   This function returns the number of entries and their `(i, j)` locations, while @ref CeedOperatorLinearAssemblePointBlockDiagonal() provides the values in the same ordering.

   This will generally be slow unless your operator is low-order.

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

   @param[in]  op          `CeedOperator` to assemble
   @param[out] num_entries Number of entries in coordinate nonzero pattern
   @param[out] rows        Row number for each entry
   @param[out] cols        Column number for each entry

   @ref User
**/
int CeedOperatorLinearAssemblePointBlockDiagonalSymbolic(CeedOperator op, CeedSize *num_entries, CeedInt **rows, CeedInt **cols) {
  Ceed          ceed;
  bool          is_composite;
  CeedInt       num_active_components, num_sub_operators;
  CeedOperator *sub_operators;

  CeedCall(CeedOperatorGetCeed(op, &ceed));
  CeedCall(CeedOperatorIsComposite(op, &is_composite));

  CeedSize input_size = 0, output_size = 0;
  CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size));
  CeedCheck(input_size == output_size, ceed, CEED_ERROR_DIMENSION, "Operator must be square");

  if (is_composite) {
    CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub_operators));
    CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators));
  } else {
    sub_operators     = &op;
    num_sub_operators = 1;
  }

  // Verify operator can be assembled correctly
  {
    CeedOperatorAssemblyData data;
    CeedInt                  num_active_elem_rstrs, comp_stride;
    CeedElemRestriction     *active_elem_rstrs;

    // Get initial values to check against
    CeedCall(CeedOperatorGetOperatorAssemblyData(sub_operators[0], &data));
    CeedCall(CeedOperatorAssemblyDataGetElemRestrictions(data, &num_active_elem_rstrs, &active_elem_rstrs, NULL, NULL));
    CeedCall(CeedElemRestrictionGetCompStride(active_elem_rstrs[0], &comp_stride));
    CeedCall(CeedElemRestrictionGetNumComponents(active_elem_rstrs[0], &num_active_components));

    // Verify that all active element restrictions have same component stride and number of components
    for (CeedInt k = 0; k < num_sub_operators; k++) {
      CeedCall(CeedOperatorGetOperatorAssemblyData(sub_operators[k], &data));
      CeedCall(CeedOperatorAssemblyDataGetElemRestrictions(data, &num_active_elem_rstrs, &active_elem_rstrs, NULL, NULL));
      for (CeedInt i = 0; i < num_active_elem_rstrs; i++) {
        CeedInt comp_stride_sub, num_active_components_sub;

        CeedCall(CeedElemRestrictionGetCompStride(active_elem_rstrs[i], &comp_stride_sub));
        CeedCheck(comp_stride == comp_stride_sub, ceed, CEED_ERROR_DIMENSION,
                  "Active element restrictions must have the same component stride: %d vs %d", comp_stride, comp_stride_sub);
        CeedCall(CeedElemRestrictionGetNumComponents(active_elem_rstrs[i], &num_active_components_sub));
        CeedCheck(num_active_components == num_active_components_sub, ceed, CEED_ERROR_INCOMPATIBLE,
                  "All suboperators must have the same number of output components."
                  " Previous: %" CeedInt_FMT " Current: %" CeedInt_FMT,
                  num_active_components, num_active_components_sub);
      }
    }
  }
  *num_entries = input_size * num_active_components;
  CeedCall(CeedCalloc(*num_entries, rows));
  CeedCall(CeedCalloc(*num_entries, cols));

  for (CeedInt o = 0; o < num_sub_operators; o++) {
    CeedElemRestriction active_elem_rstr, point_block_active_elem_rstr;
    CeedInt             comp_stride, num_elem, elem_size;
    const CeedInt      *offsets, *point_block_offsets;

    CeedCall(CeedOperatorGetActiveElemRestriction(sub_operators[o], &active_elem_rstr));
    CeedCall(CeedElemRestrictionGetCompStride(active_elem_rstr, &comp_stride));
    CeedCall(CeedElemRestrictionGetNumElements(active_elem_rstr, &num_elem));
    CeedCall(CeedElemRestrictionGetElementSize(active_elem_rstr, &elem_size));
    CeedCall(CeedElemRestrictionGetOffsets(active_elem_rstr, CEED_MEM_HOST, &offsets));

    CeedCall(CeedOperatorCreateActivePointBlockRestriction(active_elem_rstr, &point_block_active_elem_rstr));
    CeedCall(CeedElemRestrictionGetOffsets(point_block_active_elem_rstr, CEED_MEM_HOST, &point_block_offsets));

    for (CeedSize i = 0; i < num_elem * elem_size; i++) {
      for (CeedInt c_out = 0; c_out < num_active_components; c_out++) {
        for (CeedInt c_in = 0; c_in < num_active_components; c_in++) {
          (*rows)[point_block_offsets[i] + c_out * num_active_components + c_in] = offsets[i] + c_out * comp_stride;
          (*cols)[point_block_offsets[i] + c_out * num_active_components + c_in] = offsets[i] + c_in * comp_stride;
        }
      }
    }

    CeedCall(CeedElemRestrictionRestoreOffsets(active_elem_rstr, &offsets));
    CeedCall(CeedElemRestrictionRestoreOffsets(point_block_active_elem_rstr, &point_block_offsets));
    CeedCall(CeedElemRestrictionDestroy(&active_elem_rstr));
    CeedCall(CeedElemRestrictionDestroy(&point_block_active_elem_rstr));
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Assemble the point block diagonal of a square linear `CeedOperator`.

  This overwrites a `CeedVector` with the point block diagonal of a linear `CeedOperator`.

  Note: Currently only non-composite `CeedOperator` with a single field and composite `CeedOperator` with single field sub-operators are supported.

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

  @param[in]  op        `CeedOperator` to assemble `CeedQFunction`
  @param[out] assembled `CeedVector` to store assembled `CeedOperator` point block diagonal, provided in row-major form with an `num_comp * num_comp` block at each node.
                          The dimensions of this vector are derived from the active vector for the `CeedOperator`.
                          The array has shape `[nodes, component out, component in]`.
  @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 CeedOperatorLinearAssemblePointBlockDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) {
  bool     is_composite;
  CeedSize input_size = 0, output_size = 0;
  Ceed     ceed;

  CeedCall(CeedOperatorGetCeed(op, &ceed));
  CeedCall(CeedOperatorCheckReady(op));
  CeedCall(CeedOperatorIsComposite(op, &is_composite));

  CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size));
  CeedCheck(input_size == output_size, ceed, CEED_ERROR_DIMENSION, "Operator must be square");

  // Early exit for empty operator
  if (!is_composite) {
    CeedInt num_elem = 0;

    CeedCall(CeedOperatorGetNumElements(op, &num_elem));
    if (num_elem == 0) return CEED_ERROR_SUCCESS;
  }

  if (op->LinearAssemblePointBlockDiagonal) {
    // Backend version
    CeedCall(op->LinearAssemblePointBlockDiagonal(op, assembled, request));
    return CEED_ERROR_SUCCESS;
  } else if (op->LinearAssembleAddPointBlockDiagonal) {
    // Backend version with zeroing first
    CeedCall(CeedVectorSetValue(assembled, 0.0));
    CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, request));
    return CEED_ERROR_SUCCESS;
  } else {
    // Operator fallback
    CeedOperator op_fallback;

    CeedCall(CeedOperatorGetFallback(op, &op_fallback));
    if (op_fallback) {
      CeedCall(CeedOperatorLinearAssemblePointBlockDiagonal(op_fallback, assembled, request));
      return CEED_ERROR_SUCCESS;
    }
  }
  // Default interface implementation
  CeedCall(CeedVectorSetValue(assembled, 0.0));
  CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, request));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Assemble the point block diagonal of a square linear `CeedOperator`.

  This sums into a `CeedVector` with the point block diagonal of a linear `CeedOperator`.

  Note: Currently only non-composite `CeedOperator` with a single field and composite `CeedOperator` with single field sub-operators are supported.

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

  @param[in]  op        `CeedOperator` to assemble `CeedQFunction`
  @param[out] assembled `CeedVector` to store assembled CeedOperator point block diagonal, provided in row-major form with an `num_comp * num_comp` block at each node.
                          The dimensions of this vector are derived from the active vector for the `CeedOperator`.
                          The array has shape `[nodes, component out, component in]`.
  @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 CeedOperatorLinearAssembleAddPointBlockDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) {
  bool     is_composite;
  CeedSize input_size = 0, output_size = 0;
  Ceed     ceed;

  CeedCall(CeedOperatorGetCeed(op, &ceed));
  CeedCall(CeedOperatorCheckReady(op));
  CeedCall(CeedOperatorIsComposite(op, &is_composite));

  CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size));
  CeedCheck(input_size == output_size, ceed, CEED_ERROR_DIMENSION, "Operator must be square");

  // Early exit for empty operator
  if (!is_composite) {
    CeedInt num_elem = 0;

    CeedCall(CeedOperatorGetNumElements(op, &num_elem));
    if (num_elem == 0) return CEED_ERROR_SUCCESS;
  }

  if (op->LinearAssembleAddPointBlockDiagonal) {
    // Backend version
    CeedCall(op->LinearAssembleAddPointBlockDiagonal(op, assembled, request));
    return CEED_ERROR_SUCCESS;
  } else {
    // Operator fallback
    CeedOperator op_fallback;

    CeedCall(CeedOperatorGetFallback(op, &op_fallback));
    if (op_fallback) {
      CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(op_fallback, assembled, request));
      return CEED_ERROR_SUCCESS;
    }
  }
  // Default interface implementation
  if (is_composite) {
    CeedCall(CeedCompositeOperatorLinearAssembleAddDiagonal(op, request, true, assembled));
  } else {
    CeedCall(CeedSingleOperatorLinearAssembleAddDiagonal(op, request, true, assembled));
  }
  return CEED_ERROR_SUCCESS;
}

/**
   @brief Fully assemble the nonzero pattern of a linear `CeedOperator`.

   Expected to be used in conjunction with @ref CeedOperatorLinearAssemble().

   The assembly routines use coordinate format, with `num_entries` tuples of the form `(i, j, value)` which indicate that value should be added to the matrix in entry `(i, j)`.
   Note that the `(i, j)` pairs are not unique and may repeat.
   This function returns the number of entries and their `(i, j)` locations, while @ref CeedOperatorLinearAssemble() provides the values in the same ordering.

   This will generally be slow unless your operator is low-order.

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

   @param[in]  op          `CeedOperator` to assemble
   @param[out] num_entries Number of entries in coordinate nonzero pattern
   @param[out] rows        Row number for each entry
   @param[out] cols        Column number for each entry

   @ref User
**/
int CeedOperatorLinearAssembleSymbolic(CeedOperator op, CeedSize *num_entries, CeedInt **rows, CeedInt **cols) {
  bool          is_composite;
  CeedInt       num_suboperators, offset = 0;
  CeedSize      single_entries;
  CeedOperator *sub_operators;

  CeedCall(CeedOperatorCheckReady(op));
  CeedCall(CeedOperatorIsComposite(op, &is_composite));

  if (op->LinearAssembleSymbolic) {
    // Backend version
    CeedCall(op->LinearAssembleSymbolic(op, num_entries, rows, cols));
    return CEED_ERROR_SUCCESS;
  } else {
    // Operator fallback
    CeedOperator op_fallback;

    CeedCall(CeedOperatorGetFallback(op, &op_fallback));
    if (op_fallback) {
      CeedCall(CeedOperatorLinearAssembleSymbolic(op_fallback, num_entries, rows, cols));
      return CEED_ERROR_SUCCESS;
    }
  }

  // Default interface implementation

  // Count entries and allocate rows, cols arrays
  *num_entries = 0;
  if (is_composite) {
    CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators));
    CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators));
    for (CeedInt k = 0; k < num_suboperators; ++k) {
      CeedCall(CeedSingleOperatorAssemblyCountEntries(sub_operators[k], &single_entries));
      *num_entries += single_entries;
    }
  } else {
    CeedCall(CeedSingleOperatorAssemblyCountEntries(op, &single_entries));
    *num_entries += single_entries;
  }
  CeedCall(CeedCalloc(*num_entries, rows));
  CeedCall(CeedCalloc(*num_entries, cols));

  // Assemble nonzero locations
  if (is_composite) {
    CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators));
    CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators));
    for (CeedInt k = 0; k < num_suboperators; ++k) {
      CeedCall(CeedSingleOperatorAssembleSymbolic(sub_operators[k], offset, *rows, *cols));
      CeedCall(CeedSingleOperatorAssemblyCountEntries(sub_operators[k], &single_entries));
      offset += single_entries;
    }
  } else {
    CeedCall(CeedSingleOperatorAssembleSymbolic(op, offset, *rows, *cols));
  }
  return CEED_ERROR_SUCCESS;
}

/**
   @brief Fully assemble the nonzero entries of a linear operator.

   Expected to be used in conjunction with @ref CeedOperatorLinearAssembleSymbolic().

   The assembly routines use coordinate format, with `num_entries` tuples of the form `(i, j, value)` which indicate that value should be added to the matrix in entry `(i, j)`.
   Note that the `(i, j)` pairs are not unique and may repeat.
   This function returns the values of the nonzero entries to be added, their `(i, j)` locations are provided by @ref CeedOperatorLinearAssembleSymbolic().

   This will generally be slow unless your operator is low-order.

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

   @param[in]  op     `CeedOperator` to assemble
   @param[out] values Values to assemble into matrix

   @ref User
**/
int CeedOperatorLinearAssemble(CeedOperator op, CeedVector values) {
  bool          is_composite;
  CeedInt       num_suboperators, offset = 0;
  CeedSize      single_entries = 0;
  CeedOperator *sub_operators;

  CeedCall(CeedOperatorCheckReady(op));
  CeedCall(CeedOperatorIsComposite(op, &is_composite));

  // Early exit for empty operator
  if (!is_composite) {
    CeedInt num_elem = 0;

    CeedCall(CeedOperatorGetNumElements(op, &num_elem));
    if (num_elem == 0) return CEED_ERROR_SUCCESS;
  }

  if (op->LinearAssemble) {
    // Backend version
    CeedCall(op->LinearAssemble(op, values));
    return CEED_ERROR_SUCCESS;
  } else {
    // Operator fallback
    CeedOperator op_fallback;

    CeedCall(CeedOperatorGetFallback(op, &op_fallback));
    if (op_fallback) {
      CeedCall(CeedOperatorLinearAssemble(op_fallback, values));
      return CEED_ERROR_SUCCESS;
    }
  }

  // Default interface implementation
  CeedCall(CeedVectorSetValue(values, 0.0));
  if (is_composite) {
    CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators));
    CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators));
    for (CeedInt k = 0; k < num_suboperators; k++) {
      CeedCall(CeedSingleOperatorAssemble(sub_operators[k], offset, values));
      CeedCall(CeedSingleOperatorAssemblyCountEntries(sub_operators[k], &single_entries));
      offset += single_entries;
    }
  } else {
    CeedCall(CeedSingleOperatorAssemble(op, offset, values));
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the multiplicity of nodes across sub-operators in a composite `CeedOperator`.

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

  @param[in]  op               Composite `CeedOperator`
  @param[in]  num_skip_indices Number of sub-operators to skip
  @param[in]  skip_indices     Array of indices of sub-operators to skip
  @param[out] mult             Vector to store multiplicity (of size `l_size` )

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

  @ref User
**/
int CeedCompositeOperatorGetMultiplicity(CeedOperator op, CeedInt num_skip_indices, CeedInt *skip_indices, CeedVector mult) {
  Ceed                ceed;
  CeedInt             num_suboperators;
  CeedSize            l_vec_len;
  CeedScalar         *mult_array;
  CeedVector          ones_l_vec;
  CeedElemRestriction elem_rstr, mult_elem_rstr;
  CeedOperator       *sub_operators;

  CeedCall(CeedOperatorCheckReady(op));

  CeedCall(CeedOperatorGetCeed(op, &ceed));

  // Zero mult vector
  CeedCall(CeedVectorSetValue(mult, 0.0));

  // Get suboperators
  CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators));
  CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators));
  if (num_suboperators == 0) return CEED_ERROR_SUCCESS;

  // Work vector
  CeedCall(CeedVectorGetLength(mult, &l_vec_len));
  CeedCall(CeedVectorCreate(ceed, l_vec_len, &ones_l_vec));
  CeedCall(CeedVectorSetValue(ones_l_vec, 1.0));
  CeedCall(CeedVectorGetArray(mult, CEED_MEM_HOST, &mult_array));

  // Compute multiplicity across suboperators
  for (CeedInt i = 0; i < num_suboperators; i++) {
    const CeedScalar *sub_mult_array;
    CeedVector        sub_mult_l_vec, ones_e_vec;

    // -- Check for suboperator to skip
    for (CeedInt j = 0; j < num_skip_indices; j++) {
      if (skip_indices[j] == i) continue;
    }

    // -- Sub operator multiplicity
    CeedCall(CeedOperatorGetActiveElemRestriction(sub_operators[i], &elem_rstr));
    CeedCall(CeedElemRestrictionCreateUnorientedCopy(elem_rstr, &mult_elem_rstr));
    CeedCall(CeedElemRestrictionDestroy(&elem_rstr));
    CeedCall(CeedElemRestrictionCreateVector(mult_elem_rstr, &sub_mult_l_vec, &ones_e_vec));
    CeedCall(CeedVectorSetValue(sub_mult_l_vec, 0.0));
    CeedCall(CeedElemRestrictionApply(mult_elem_rstr, CEED_NOTRANSPOSE, ones_l_vec, ones_e_vec, CEED_REQUEST_IMMEDIATE));
    CeedCall(CeedElemRestrictionApply(mult_elem_rstr, CEED_TRANSPOSE, ones_e_vec, sub_mult_l_vec, CEED_REQUEST_IMMEDIATE));
    CeedCall(CeedVectorGetArrayRead(sub_mult_l_vec, CEED_MEM_HOST, &sub_mult_array));
    // ---- Flag every node present in the current suboperator
    for (CeedSize j = 0; j < l_vec_len; j++) {
      if (sub_mult_array[j] > 0.0) mult_array[j] += 1.0;
    }
    CeedCall(CeedVectorRestoreArrayRead(sub_mult_l_vec, &sub_mult_array));
    CeedCall(CeedVectorDestroy(&sub_mult_l_vec));
    CeedCall(CeedVectorDestroy(&ones_e_vec));
    CeedCall(CeedElemRestrictionDestroy(&mult_elem_rstr));
  }
  CeedCall(CeedVectorRestoreArray(mult, &mult_array));
  CeedCall(CeedVectorDestroy(&ones_l_vec));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Create a multigrid coarse `CeedOperator` and level transfer `CeedOperator` for a `CeedOperator`, creating the prolongation basis from the fine and coarse grid interpolation.

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

  @param[in]  op_fine      Fine grid `CeedOperator`
  @param[in]  p_mult_fine  L-vector multiplicity in parallel gather/scatter, or `NULL` if not creating prolongation/restriction `CeedOperator`
  @param[in]  rstr_coarse  Coarse grid `CeedElemRestriction`
  @param[in]  basis_coarse Coarse grid active vector `CeedBasis`
  @param[out] op_coarse    Coarse grid `CeedOperator`
  @param[out] op_prolong   Coarse to fine `CeedOperator`, or `NULL`
  @param[out] op_restrict  Fine to coarse `CeedOperator`, or `NULL`

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

  @ref User
**/
int CeedOperatorMultigridLevelCreate(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse,
                                     CeedOperator *op_coarse, CeedOperator *op_prolong, CeedOperator *op_restrict) {
  CeedBasis basis_c_to_f = NULL;

  CeedCall(CeedOperatorCheckReady(op_fine));

  // Build prolongation matrix, if required
  if (op_prolong || op_restrict) {
    CeedBasis basis_fine;

    CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine));
    CeedCall(CeedBasisCreateProjection(basis_coarse, basis_fine, &basis_c_to_f));
    CeedCall(CeedBasisDestroy(&basis_fine));
  }

  // Core code
  CeedCall(CeedSingleOperatorMultigridLevel(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, op_restrict));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Create a multigrid coarse `CeedOperator` and level transfer `CeedOperator` for a `CeedOperator` with a tensor basis for the active basis.

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

  @param[in]  op_fine       Fine grid `CeedOperator`
  @param[in]  p_mult_fine   L-vector multiplicity in parallel gather/scatter, or `NULL` if not creating prolongation/restriction `CeedOperator`
  @param[in]  rstr_coarse   Coarse grid `CeedElemRestriction`
  @param[in]  basis_coarse  Coarse grid active vector `CeedBasis`
  @param[in]  interp_c_to_f Matrix for coarse to fine interpolation, or `NULL` if not creating prolongation/restriction `CeedOperator`
  @param[out] op_coarse     Coarse grid `CeedOperator`
  @param[out] op_prolong    Coarse to fine `CeedOperator`, or `NULL`
  @param[out] op_restrict   Fine to coarse `CeedOperator`, or `NULL`

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

  @ref User
**/
int CeedOperatorMultigridLevelCreateTensorH1(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse,
                                             const CeedScalar *interp_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong,
                                             CeedOperator *op_restrict) {
  Ceed      ceed;
  CeedInt   Q_f, Q_c;
  CeedBasis basis_fine, basis_c_to_f = NULL;

  CeedCall(CeedOperatorCheckReady(op_fine));
  CeedCall(CeedOperatorGetCeed(op_fine, &ceed));

  // Check for compatible quadrature spaces
  CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine));
  CeedCall(CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f));
  CeedCall(CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c));
  CeedCheck(Q_f == Q_c, ceed, CEED_ERROR_DIMENSION,
            "Bases must have compatible quadrature spaces."
            " Fine grid: %" CeedInt_FMT " points, Coarse grid: %" CeedInt_FMT " points",
            Q_f, Q_c);

  // Create coarse to fine basis, if required
  if (op_prolong || op_restrict) {
    CeedInt     dim, num_comp, num_nodes_c, P_1d_f, P_1d_c;
    CeedScalar *q_ref, *q_weight, *grad;

    // Check if interpolation matrix is provided
    CeedCheck(interp_c_to_f, ceed, CEED_ERROR_INCOMPATIBLE,
              "Prolongation or restriction operator creation requires coarse-to-fine interpolation matrix");
    CeedCall(CeedBasisGetDimension(basis_fine, &dim));
    CeedCall(CeedBasisGetNumComponents(basis_fine, &num_comp));
    CeedCall(CeedBasisGetNumNodes1D(basis_fine, &P_1d_f));
    CeedCall(CeedBasisDestroy(&basis_fine));
    CeedCall(CeedElemRestrictionGetElementSize(rstr_coarse, &num_nodes_c));
    P_1d_c = dim == 1 ? num_nodes_c : dim == 2 ? sqrt(num_nodes_c) : cbrt(num_nodes_c);
    CeedCall(CeedCalloc(P_1d_f, &q_ref));
    CeedCall(CeedCalloc(P_1d_f, &q_weight));
    CeedCall(CeedCalloc(P_1d_f * P_1d_c * dim, &grad));
    CeedCall(CeedBasisCreateTensorH1(ceed, dim, num_comp, P_1d_c, P_1d_f, interp_c_to_f, grad, q_ref, q_weight, &basis_c_to_f));
    CeedCall(CeedFree(&q_ref));
    CeedCall(CeedFree(&q_weight));
    CeedCall(CeedFree(&grad));
  }

  // Core code
  CeedCall(CeedSingleOperatorMultigridLevel(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, op_restrict));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Create a multigrid coarse `CeedOperator` and level transfer `CeedOperator` for a `CeedOperator` with a non-tensor basis for the active vector

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

  @param[in]  op_fine       Fine grid `CeedOperator`
  @param[in]  p_mult_fine   L-vector multiplicity in parallel gather/scatter, or `NULL` if not creating prolongation/restriction `CeedOperator`
  @param[in]  rstr_coarse   Coarse grid `CeedElemRestriction`
  @param[in]  basis_coarse  Coarse grid active vector `CeedBasis`
  @param[in]  interp_c_to_f Matrix for coarse to fine interpolation, or `NULL` if not creating prolongation/restriction `CeedOperator`
  @param[out] op_coarse     Coarse grid `CeedOperator`
  @param[out] op_prolong    Coarse to fine `CeedOperator`, or `NULL`
  @param[out] op_restrict   Fine to coarse `CeedOperator`, or `NULL`

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

  @ref User
**/
int CeedOperatorMultigridLevelCreateH1(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse,
                                       const CeedScalar *interp_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong,
                                       CeedOperator *op_restrict) {
  Ceed      ceed;
  CeedInt   Q_f, Q_c;
  CeedBasis basis_fine, basis_c_to_f = NULL;

  CeedCall(CeedOperatorCheckReady(op_fine));
  CeedCall(CeedOperatorGetCeed(op_fine, &ceed));

  // Check for compatible quadrature spaces
  CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine));
  CeedCall(CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f));
  CeedCall(CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c));
  CeedCheck(Q_f == Q_c, ceed, CEED_ERROR_DIMENSION, "Bases must have compatible quadrature spaces");

  // Coarse to fine basis
  if (op_prolong || op_restrict) {
    CeedInt          dim, num_comp, num_nodes_c, num_nodes_f;
    CeedScalar      *q_ref, *q_weight, *grad;
    CeedElemTopology topo;

    // Check if interpolation matrix is provided
    CeedCheck(interp_c_to_f, ceed, CEED_ERROR_INCOMPATIBLE,
              "Prolongation or restriction operator creation requires coarse-to-fine interpolation matrix");
    CeedCall(CeedBasisGetTopology(basis_fine, &topo));
    CeedCall(CeedBasisGetDimension(basis_fine, &dim));
    CeedCall(CeedBasisGetNumComponents(basis_fine, &num_comp));
    CeedCall(CeedBasisGetNumNodes(basis_fine, &num_nodes_f));
    CeedCall(CeedBasisDestroy(&basis_fine));
    CeedCall(CeedElemRestrictionGetElementSize(rstr_coarse, &num_nodes_c));
    CeedCall(CeedCalloc(num_nodes_f * dim, &q_ref));
    CeedCall(CeedCalloc(num_nodes_f, &q_weight));
    CeedCall(CeedCalloc(num_nodes_f * num_nodes_c * dim, &grad));
    CeedCall(CeedBasisCreateH1(ceed, topo, num_comp, num_nodes_c, num_nodes_f, interp_c_to_f, grad, q_ref, q_weight, &basis_c_to_f));
    CeedCall(CeedFree(&q_ref));
    CeedCall(CeedFree(&q_weight));
    CeedCall(CeedFree(&grad));
  }

  // Core code
  CeedCall(CeedSingleOperatorMultigridLevel(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, op_restrict));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Build a FDM based approximate inverse for each element for a `CeedOperator`.

  This returns a `CeedOperator` and `CeedVector` to apply a Fast Diagonalization Method based approximate inverse.
  This function obtains the simultaneous diagonalization for the 1D mass and Laplacian operators, \f$M = V^T V, K = V^T S V\f$.
  The assembled `CeedQFunction` is used to modify the eigenvalues from simultaneous diagonalization and obtain an approximate inverse of the form \f$V^T \hat S V\f$.
  The `CeedOperator` must be linear and non-composite.
  The associated `CeedQFunction` must therefore also be linear.

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

  @param[in]  op      `CeedOperator` to create element inverses
  @param[out] fdm_inv `CeedOperator` to apply the action of a FDM based inverse for each element
  @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 CeedOperatorCreateFDMElementInverse(CeedOperator op, CeedOperator *fdm_inv, CeedRequest *request) {
  Ceed                 ceed, ceed_parent;
  bool                 interp = false, grad = false, is_tensor_basis = true;
  CeedInt              num_input_fields, P_1d, Q_1d, num_nodes, num_qpts, dim, num_comp = 1, num_elem = 1;
  CeedScalar          *mass, *laplace, *x, *fdm_interp, *lambda, *elem_avg;
  const CeedScalar    *interp_1d, *grad_1d, *q_weight_1d;
  CeedVector           q_data;
  CeedElemRestriction  rstr  = NULL, rstr_qd_i;
  CeedBasis            basis = NULL, fdm_basis;
  CeedQFunctionContext ctx_fdm;
  CeedQFunctionField  *qf_fields;
  CeedQFunction        qf, qf_fdm;
  CeedOperatorField   *op_fields;

  CeedCall(CeedOperatorCheckReady(op));

  if (op->CreateFDMElementInverse) {
    // Backend version
    CeedCall(op->CreateFDMElementInverse(op, fdm_inv, request));
    return CEED_ERROR_SUCCESS;
  } else {
    // Operator fallback
    CeedOperator op_fallback;

    CeedCall(CeedOperatorGetFallback(op, &op_fallback));
    if (op_fallback) {
      CeedCall(CeedOperatorCreateFDMElementInverse(op_fallback, fdm_inv, request));
      return CEED_ERROR_SUCCESS;
    }
  }

  // Default interface implementation
  CeedCall(CeedOperatorGetCeed(op, &ceed));
  CeedCall(CeedOperatorGetFallbackParentCeed(op, &ceed_parent));
  CeedCall(CeedOperatorGetQFunction(op, &qf));

  // Determine active input basis
  CeedCall(CeedOperatorGetFields(op, &num_input_fields, &op_fields, NULL, NULL));
  CeedCall(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL));
  for (CeedInt i = 0; i < num_input_fields; i++) {
    CeedVector vec;

    CeedCall(CeedOperatorFieldGetVector(op_fields[i], &vec));
    if (vec == CEED_VECTOR_ACTIVE) {
      CeedEvalMode eval_mode;

      CeedCall(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode));
      interp = interp || eval_mode == CEED_EVAL_INTERP;
      grad   = grad || eval_mode == CEED_EVAL_GRAD;
      if (!basis) CeedCall(CeedOperatorFieldGetBasis(op_fields[i], &basis));
      if (!rstr) CeedCall(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr));
    }
    CeedCall(CeedVectorDestroy(&vec));
  }
  CeedCheck(basis, ceed, CEED_ERROR_BACKEND, "No active field set");
  CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d));
  CeedCall(CeedBasisGetNumNodes(basis, &num_nodes));
  CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
  CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts));
  CeedCall(CeedBasisGetDimension(basis, &dim));
  CeedCall(CeedBasisGetNumComponents(basis, &num_comp));
  CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem));

  // Build and diagonalize 1D Mass and Laplacian
  CeedCall(CeedBasisIsTensor(basis, &is_tensor_basis));
  CeedCheck(is_tensor_basis, ceed, CEED_ERROR_BACKEND, "FDMElementInverse only supported for tensor bases");
  CeedCall(CeedCalloc(P_1d * P_1d, &mass));
  CeedCall(CeedCalloc(P_1d * P_1d, &laplace));
  CeedCall(CeedCalloc(P_1d * P_1d, &x));
  CeedCall(CeedCalloc(P_1d * P_1d, &fdm_interp));
  CeedCall(CeedCalloc(P_1d, &lambda));
  // -- Build matrices
  CeedCall(CeedBasisGetInterp1D(basis, &interp_1d));
  CeedCall(CeedBasisGetGrad1D(basis, &grad_1d));
  CeedCall(CeedBasisGetQWeights(basis, &q_weight_1d));
  CeedCall(CeedBuildMassLaplace(interp_1d, grad_1d, q_weight_1d, P_1d, Q_1d, dim, mass, laplace));

  // -- Diagonalize
  CeedCall(CeedSimultaneousDiagonalization(ceed, laplace, mass, x, lambda, P_1d));
  CeedCall(CeedFree(&mass));
  CeedCall(CeedFree(&laplace));
  for (CeedInt i = 0; i < P_1d; i++) {
    for (CeedInt j = 0; j < P_1d; j++) fdm_interp[i + j * P_1d] = x[j + i * P_1d];
  }
  CeedCall(CeedFree(&x));

  {
    CeedInt             layout[3], num_modes = (interp ? 1 : 0) + (grad ? dim : 0);
    CeedScalar          max_norm = 0;
    const CeedScalar   *assembled_array, *q_weight_array;
    CeedVector          assembled = NULL, q_weight;
    CeedElemRestriction rstr_qf   = NULL;

    // Assemble QFunction
    CeedCall(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled, &rstr_qf, request));
    CeedCall(CeedElemRestrictionGetELayout(rstr_qf, layout));
    CeedCall(CeedElemRestrictionDestroy(&rstr_qf));
    CeedCall(CeedVectorNorm(assembled, CEED_NORM_MAX, &max_norm));

    // Calculate element averages
    CeedCall(CeedVectorCreate(ceed_parent, num_qpts, &q_weight));
    CeedCall(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_weight));
    CeedCall(CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array));
    CeedCall(CeedVectorGetArrayRead(q_weight, CEED_MEM_HOST, &q_weight_array));
    CeedCall(CeedCalloc(num_elem, &elem_avg));
    const CeedScalar qf_value_bound = max_norm * 100 * CEED_EPSILON;

    for (CeedInt e = 0; e < num_elem; e++) {
      CeedInt count = 0;

      for (CeedInt q = 0; q < num_qpts; q++) {
        for (CeedInt i = 0; i < num_comp * num_comp * num_modes * num_modes; i++) {
          if (fabs(assembled_array[q * layout[0] + i * layout[1] + e * layout[2]]) > qf_value_bound) {
            elem_avg[e] += assembled_array[q * layout[0] + i * layout[1] + e * layout[2]] / q_weight_array[q];
            count++;
          }
        }
      }
      if (count) {
        elem_avg[e] /= count;
      } else {
        elem_avg[e] = 1.0;
      }
    }
    CeedCall(CeedVectorRestoreArrayRead(assembled, &assembled_array));
    CeedCall(CeedVectorDestroy(&assembled));
    CeedCall(CeedVectorRestoreArrayRead(q_weight, &q_weight_array));
    CeedCall(CeedVectorDestroy(&q_weight));
  }

  // Build FDM diagonal
  {
    CeedScalar *q_data_array, *fdm_diagonal;

    CeedCall(CeedCalloc(num_comp * num_nodes, &fdm_diagonal));
    const CeedScalar fdm_diagonal_bound = num_nodes * CEED_EPSILON;
    for (CeedInt c = 0; c < num_comp; c++) {
      for (CeedInt n = 0; n < num_nodes; n++) {
        if (interp) fdm_diagonal[c * num_nodes + n] = 1.0;
        if (grad) {
          for (CeedInt d = 0; d < dim; d++) {
            CeedInt i = (n / CeedIntPow(P_1d, d)) % P_1d;
            fdm_diagonal[c * num_nodes + n] += lambda[i];
          }
        }
        if (fabs(fdm_diagonal[c * num_nodes + n]) < fdm_diagonal_bound) fdm_diagonal[c * num_nodes + n] = fdm_diagonal_bound;
      }
    }
    CeedCall(CeedVectorCreate(ceed_parent, num_elem * num_comp * num_nodes, &q_data));
    CeedCall(CeedVectorSetValue(q_data, 0.0));
    CeedCall(CeedVectorGetArrayWrite(q_data, CEED_MEM_HOST, &q_data_array));
    for (CeedInt e = 0; e < num_elem; e++) {
      for (CeedInt c = 0; c < num_comp; c++) {
        for (CeedInt n = 0; n < num_nodes; n++) {
          q_data_array[(e * num_comp + c) * num_nodes + n] = 1. / (elem_avg[e] * fdm_diagonal[c * num_nodes + n]);
        }
      }
    }
    CeedCall(CeedFree(&elem_avg));
    CeedCall(CeedFree(&fdm_diagonal));
    CeedCall(CeedVectorRestoreArray(q_data, &q_data_array));
  }

  // Setup FDM operator
  // -- Basis
  {
    CeedScalar *grad_dummy, *q_ref_dummy, *q_weight_dummy;

    CeedCall(CeedCalloc(P_1d * P_1d, &grad_dummy));
    CeedCall(CeedCalloc(P_1d, &q_ref_dummy));
    CeedCall(CeedCalloc(P_1d, &q_weight_dummy));
    CeedCall(CeedBasisCreateTensorH1(ceed_parent, dim, num_comp, P_1d, P_1d, fdm_interp, grad_dummy, q_ref_dummy, q_weight_dummy, &fdm_basis));
    CeedCall(CeedFree(&fdm_interp));
    CeedCall(CeedFree(&grad_dummy));
    CeedCall(CeedFree(&q_ref_dummy));
    CeedCall(CeedFree(&q_weight_dummy));
    CeedCall(CeedFree(&lambda));
  }

  // -- Restriction
  {
    CeedInt strides[3] = {1, num_nodes, num_nodes * num_comp};
    CeedCall(CeedElemRestrictionCreateStrided(ceed_parent, num_elem, num_nodes, num_comp,
                                              (CeedSize)num_elem * (CeedSize)num_comp * (CeedSize)num_nodes, strides, &rstr_qd_i));
  }

  // -- QFunction
  CeedCall(CeedQFunctionCreateInteriorByName(ceed_parent, "Scale", &qf_fdm));
  CeedCall(CeedQFunctionAddInput(qf_fdm, "input", num_comp, CEED_EVAL_INTERP));
  CeedCall(CeedQFunctionAddInput(qf_fdm, "scale", num_comp, CEED_EVAL_NONE));
  CeedCall(CeedQFunctionAddOutput(qf_fdm, "output", num_comp, CEED_EVAL_INTERP));
  CeedCall(CeedQFunctionSetUserFlopsEstimate(qf_fdm, num_comp));

  // -- QFunction context
  {
    CeedInt *num_comp_data;

    CeedCall(CeedCalloc(1, &num_comp_data));
    num_comp_data[0] = num_comp;
    CeedCall(CeedQFunctionContextCreate(ceed, &ctx_fdm));
    CeedCall(CeedQFunctionContextSetData(ctx_fdm, CEED_MEM_HOST, CEED_OWN_POINTER, sizeof(*num_comp_data), num_comp_data));
  }
  CeedCall(CeedQFunctionSetContext(qf_fdm, ctx_fdm));
  CeedCall(CeedQFunctionContextDestroy(&ctx_fdm));

  // -- Operator
  CeedCall(CeedOperatorCreate(ceed_parent, qf_fdm, NULL, NULL, fdm_inv));
  CeedCall(CeedOperatorSetField(*fdm_inv, "input", rstr, fdm_basis, CEED_VECTOR_ACTIVE));
  CeedCall(CeedOperatorSetField(*fdm_inv, "scale", rstr_qd_i, CEED_BASIS_NONE, q_data));
  CeedCall(CeedOperatorSetField(*fdm_inv, "output", rstr, fdm_basis, CEED_VECTOR_ACTIVE));

  // Cleanup
  CeedCall(CeedVectorDestroy(&q_data));
  CeedCall(CeedElemRestrictionDestroy(&rstr));
  CeedCall(CeedElemRestrictionDestroy(&rstr_qd_i));
  CeedCall(CeedBasisDestroy(&basis));
  CeedCall(CeedBasisDestroy(&fdm_basis));
  CeedCall(CeedQFunctionDestroy(&qf_fdm));
  return CEED_ERROR_SUCCESS;
}

/// @}