xref: /libCEED/rust/libceed-sys/c-src/interface/ceed-operator.c (revision b3b7035fa845d84f3b306f1f89cf7c89e3d6d1a4)
1d7b241e6Sjeremylt // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
2d7b241e6Sjeremylt // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
3d7b241e6Sjeremylt // reserved. See files LICENSE and NOTICE for details.
4d7b241e6Sjeremylt //
5d7b241e6Sjeremylt // This file is part of CEED, a collection of benchmarks, miniapps, software
6d7b241e6Sjeremylt // libraries and APIs for efficient high-order finite element and spectral
7d7b241e6Sjeremylt // element discretizations for exascale applications. For more information and
8d7b241e6Sjeremylt // source code availability see http://github.com/ceed.
9d7b241e6Sjeremylt //
10d7b241e6Sjeremylt // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
11d7b241e6Sjeremylt // a collaborative effort of two U.S. Department of Energy organizations (Office
12d7b241e6Sjeremylt // of Science and the National Nuclear Security Administration) responsible for
13d7b241e6Sjeremylt // the planning and preparation of a capable exascale ecosystem, including
14d7b241e6Sjeremylt // software, applications, hardware, advanced system engineering and early
15d7b241e6Sjeremylt // testbed platforms, in support of the nation's exascale computing imperative.
16d7b241e6Sjeremylt 
17d7b241e6Sjeremylt #include <ceed-impl.h>
18d863ab9bSjeremylt #include <ceed-backend.h>
19d7b241e6Sjeremylt #include <string.h>
203bd813ffSjeremylt #include <math.h>
21d7b241e6Sjeremylt 
22dfdf5a53Sjeremylt /// @file
237a982d89SJeremy L. Thompson /// Implementation of CeedOperator interfaces
247a982d89SJeremy L. Thompson 
257a982d89SJeremy L. Thompson /// ----------------------------------------------------------------------------
267a982d89SJeremy L. Thompson /// CeedOperator Library Internal Functions
277a982d89SJeremy L. Thompson /// ----------------------------------------------------------------------------
287a982d89SJeremy L. Thompson /// @addtogroup CeedOperatorDeveloper
297a982d89SJeremy L. Thompson /// @{
307a982d89SJeremy L. Thompson 
317a982d89SJeremy L. Thompson /**
327a982d89SJeremy L. Thompson   @brief Duplicate a CeedOperator with a reference Ceed to fallback for advanced
337a982d89SJeremy L. Thompson            CeedOperator functionality
347a982d89SJeremy L. Thompson 
357a982d89SJeremy L. Thompson   @param op           CeedOperator to create fallback for
367a982d89SJeremy L. Thompson 
377a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
387a982d89SJeremy L. Thompson 
397a982d89SJeremy L. Thompson   @ref Developer
407a982d89SJeremy L. Thompson **/
417a982d89SJeremy L. Thompson int CeedOperatorCreateFallback(CeedOperator op) {
427a982d89SJeremy L. Thompson   int ierr;
437a982d89SJeremy L. Thompson 
447a982d89SJeremy L. Thompson   // Fallback Ceed
457a982d89SJeremy L. Thompson   const char *resource, *fallbackresource;
467a982d89SJeremy L. Thompson   ierr = CeedGetResource(op->ceed, &resource); CeedChk(ierr);
477a982d89SJeremy L. Thompson   ierr = CeedGetOperatorFallbackResource(op->ceed, &fallbackresource);
487a982d89SJeremy L. Thompson   CeedChk(ierr);
497a982d89SJeremy L. Thompson   if (!strcmp(resource, fallbackresource))
507a982d89SJeremy L. Thompson     // LCOV_EXCL_START
517a982d89SJeremy L. Thompson     return CeedError(op->ceed, 1, "Backend %s cannot create an operator"
527a982d89SJeremy L. Thompson                      "fallback to resource %s", resource, fallbackresource);
537a982d89SJeremy L. Thompson   // LCOV_EXCL_STOP
547a982d89SJeremy L. Thompson 
557a982d89SJeremy L. Thompson   // Fallback Ceed
567a982d89SJeremy L. Thompson   Ceed ceedref;
577a982d89SJeremy L. Thompson   if (!op->ceed->opfallbackceed) {
587a982d89SJeremy L. Thompson     ierr = CeedInit(fallbackresource, &ceedref); CeedChk(ierr);
597a982d89SJeremy L. Thompson     ceedref->opfallbackparent = op->ceed;
607a982d89SJeremy L. Thompson     op->ceed->opfallbackceed = ceedref;
617a982d89SJeremy L. Thompson   }
627a982d89SJeremy L. Thompson   ceedref = op->ceed->opfallbackceed;
637a982d89SJeremy L. Thompson 
647a982d89SJeremy L. Thompson   // Clone Op
657a982d89SJeremy L. Thompson   CeedOperator opref;
667a982d89SJeremy L. Thompson   ierr = CeedCalloc(1, &opref); CeedChk(ierr);
677a982d89SJeremy L. Thompson   memcpy(opref, op, sizeof(*opref)); CeedChk(ierr);
687a982d89SJeremy L. Thompson   opref->data = NULL;
697a982d89SJeremy L. Thompson   opref->setupdone = 0;
707a982d89SJeremy L. Thompson   opref->ceed = ceedref;
717a982d89SJeremy L. Thompson   ierr = ceedref->OperatorCreate(opref); CeedChk(ierr);
727a982d89SJeremy L. Thompson   op->opfallback = opref;
737a982d89SJeremy L. Thompson 
747a982d89SJeremy L. Thompson   // Clone QF
757a982d89SJeremy L. Thompson   CeedQFunction qfref;
767a982d89SJeremy L. Thompson   ierr = CeedCalloc(1, &qfref); CeedChk(ierr);
777a982d89SJeremy L. Thompson   memcpy(qfref, (op->qf), sizeof(*qfref)); CeedChk(ierr);
787a982d89SJeremy L. Thompson   qfref->data = NULL;
797a982d89SJeremy L. Thompson   qfref->ceed = ceedref;
807a982d89SJeremy L. Thompson   ierr = ceedref->QFunctionCreate(qfref); CeedChk(ierr);
817a982d89SJeremy L. Thompson   opref->qf = qfref;
827a982d89SJeremy L. Thompson   op->qffallback = qfref;
837a982d89SJeremy L. Thompson 
847a982d89SJeremy L. Thompson   return 0;
857a982d89SJeremy L. Thompson }
867a982d89SJeremy L. Thompson 
877a982d89SJeremy L. Thompson /**
887a982d89SJeremy L. Thompson   @brief Check if a CeedOperator is ready to be used.
897a982d89SJeremy L. Thompson 
907a982d89SJeremy L. Thompson   @param[in] ceed Ceed object for error handling
917a982d89SJeremy L. Thompson   @param[in] op   CeedOperator to check
927a982d89SJeremy L. Thompson 
937a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
947a982d89SJeremy L. Thompson 
957a982d89SJeremy L. Thompson   @ref Developer
967a982d89SJeremy L. Thompson **/
977a982d89SJeremy L. Thompson static int CeedOperatorCheckReady(Ceed ceed, CeedOperator op) {
987a982d89SJeremy L. Thompson   CeedQFunction qf = op->qf;
997a982d89SJeremy L. Thompson 
1007a982d89SJeremy L. Thompson   if (op->composite) {
1017a982d89SJeremy L. Thompson     if (!op->numsub)
1027a982d89SJeremy L. Thompson       // LCOV_EXCL_START
1037a982d89SJeremy L. Thompson       return CeedError(ceed, 1, "No suboperators set");
1047a982d89SJeremy L. Thompson     // LCOV_EXCL_STOP
1057a982d89SJeremy L. Thompson   } else {
1067a982d89SJeremy L. Thompson     if (op->nfields == 0)
1077a982d89SJeremy L. Thompson       // LCOV_EXCL_START
1087a982d89SJeremy L. Thompson       return CeedError(ceed, 1, "No operator fields set");
1097a982d89SJeremy L. Thompson     // LCOV_EXCL_STOP
1107a982d89SJeremy L. Thompson     if (op->nfields < qf->numinputfields + qf->numoutputfields)
1117a982d89SJeremy L. Thompson       // LCOV_EXCL_START
1127a982d89SJeremy L. Thompson       return CeedError(ceed, 1, "Not all operator fields set");
1137a982d89SJeremy L. Thompson     // LCOV_EXCL_STOP
1147a982d89SJeremy L. Thompson     if (!op->hasrestriction)
1157a982d89SJeremy L. Thompson       // LCOV_EXCL_START
1167a982d89SJeremy L. Thompson       return CeedError(ceed, 1,"At least one restriction required");
1177a982d89SJeremy L. Thompson     // LCOV_EXCL_STOP
1187a982d89SJeremy L. Thompson     if (op->numqpoints == 0)
1197a982d89SJeremy L. Thompson       // LCOV_EXCL_START
1207a982d89SJeremy L. Thompson       return CeedError(ceed, 1,"At least one non-collocated basis required");
1217a982d89SJeremy L. Thompson     // LCOV_EXCL_STOP
1227a982d89SJeremy L. Thompson   }
1237a982d89SJeremy L. Thompson 
1247a982d89SJeremy L. Thompson   return 0;
1257a982d89SJeremy L. Thompson }
1267a982d89SJeremy L. Thompson 
1277a982d89SJeremy L. Thompson /**
1287a982d89SJeremy L. Thompson   @brief View a field of a CeedOperator
1297a982d89SJeremy L. Thompson 
1307a982d89SJeremy L. Thompson   @param[in] field       Operator field to view
1314c4400c7SValeria Barra   @param[in] qffield     QFunction field (carries field name)
1327a982d89SJeremy L. Thompson   @param[in] fieldnumber Number of field being viewed
1334c4400c7SValeria Barra   @param[in] sub         true indicates sub-operator, which increases indentation; false for top-level operator
1344c4400c7SValeria Barra   @param[in] in          true for an input field; false for output field
1357a982d89SJeremy L. Thompson   @param[in] stream      Stream to view to, e.g., stdout
1367a982d89SJeremy L. Thompson 
1377a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
1387a982d89SJeremy L. Thompson 
1397a982d89SJeremy L. Thompson   @ref Utility
1407a982d89SJeremy L. Thompson **/
1417a982d89SJeremy L. Thompson static int CeedOperatorFieldView(CeedOperatorField field,
1427a982d89SJeremy L. Thompson                                  CeedQFunctionField qffield,
1437a982d89SJeremy L. Thompson                                  CeedInt fieldnumber, bool sub, bool in,
1447a982d89SJeremy L. Thompson                                  FILE *stream) {
1457a982d89SJeremy L. Thompson   const char *pre = sub ? "  " : "";
1467a982d89SJeremy L. Thompson   const char *inout = in ? "Input" : "Output";
1477a982d89SJeremy L. Thompson 
1487a982d89SJeremy L. Thompson   fprintf(stream, "%s    %s Field [%d]:\n"
1497a982d89SJeremy L. Thompson           "%s      Name: \"%s\"\n",
1507a982d89SJeremy L. Thompson           pre, inout, fieldnumber, pre, qffield->fieldname);
1517a982d89SJeremy L. Thompson 
1527a982d89SJeremy L. Thompson   if (field->basis == CEED_BASIS_COLLOCATED)
1537a982d89SJeremy L. Thompson     fprintf(stream, "%s      Collocated basis\n", pre);
1547a982d89SJeremy L. Thompson 
1557a982d89SJeremy L. Thompson   if (field->vec == CEED_VECTOR_ACTIVE)
1567a982d89SJeremy L. Thompson     fprintf(stream, "%s      Active vector\n", pre);
1577a982d89SJeremy L. Thompson   else if (field->vec == CEED_VECTOR_NONE)
1587a982d89SJeremy L. Thompson     fprintf(stream, "%s      No vector\n", pre);
1597a982d89SJeremy L. Thompson 
1607a982d89SJeremy L. Thompson   return 0;
1617a982d89SJeremy L. Thompson }
1627a982d89SJeremy L. Thompson 
1637a982d89SJeremy L. Thompson /**
1647a982d89SJeremy L. Thompson   @brief View a single CeedOperator
1657a982d89SJeremy L. Thompson 
1667a982d89SJeremy L. Thompson   @param[in] op     CeedOperator to view
1677a982d89SJeremy L. Thompson   @param[in] sub    Boolean flag for sub-operator
1687a982d89SJeremy L. Thompson   @param[in] stream Stream to write; typically stdout/stderr or a file
1697a982d89SJeremy L. Thompson 
1707a982d89SJeremy L. Thompson   @return Error code: 0 - success, otherwise - failure
1717a982d89SJeremy L. Thompson 
1727a982d89SJeremy L. Thompson   @ref Utility
1737a982d89SJeremy L. Thompson **/
1747a982d89SJeremy L. Thompson int CeedOperatorSingleView(CeedOperator op, bool sub, FILE *stream) {
1757a982d89SJeremy L. Thompson   int ierr;
1767a982d89SJeremy L. Thompson   const char *pre = sub ? "  " : "";
1777a982d89SJeremy L. Thompson 
1787a982d89SJeremy L. Thompson   CeedInt totalfields;
1797a982d89SJeremy L. Thompson   ierr = CeedOperatorGetNumArgs(op, &totalfields); CeedChk(ierr);
1807a982d89SJeremy L. Thompson 
1817a982d89SJeremy L. Thompson   fprintf(stream, "%s  %d Field%s\n", pre, totalfields, totalfields>1 ? "s" : "");
1827a982d89SJeremy L. Thompson 
1837a982d89SJeremy L. Thompson   fprintf(stream, "%s  %d Input Field%s:\n", pre, op->qf->numinputfields,
1847a982d89SJeremy L. Thompson           op->qf->numinputfields>1 ? "s" : "");
1857a982d89SJeremy L. Thompson   for (CeedInt i=0; i<op->qf->numinputfields; i++) {
1867a982d89SJeremy L. Thompson     ierr = CeedOperatorFieldView(op->inputfields[i], op->qf->inputfields[i],
1877a982d89SJeremy L. Thompson                                  i, sub, 1, stream); CeedChk(ierr);
1887a982d89SJeremy L. Thompson   }
1897a982d89SJeremy L. Thompson 
1907a982d89SJeremy L. Thompson   fprintf(stream, "%s  %d Output Field%s:\n", pre, op->qf->numoutputfields,
1917a982d89SJeremy L. Thompson           op->qf->numoutputfields>1 ? "s" : "");
1927a982d89SJeremy L. Thompson   for (CeedInt i=0; i<op->qf->numoutputfields; i++) {
1937a982d89SJeremy L. Thompson     ierr = CeedOperatorFieldView(op->outputfields[i], op->qf->inputfields[i],
1947a982d89SJeremy L. Thompson                                  i, sub, 0, stream); CeedChk(ierr);
1957a982d89SJeremy L. Thompson   }
1967a982d89SJeremy L. Thompson 
1977a982d89SJeremy L. Thompson   return 0;
1987a982d89SJeremy L. Thompson }
1997a982d89SJeremy L. Thompson 
2007a982d89SJeremy L. Thompson /// @}
2017a982d89SJeremy L. Thompson 
2027a982d89SJeremy L. Thompson /// ----------------------------------------------------------------------------
2037a982d89SJeremy L. Thompson /// CeedOperator Backend API
2047a982d89SJeremy L. Thompson /// ----------------------------------------------------------------------------
2057a982d89SJeremy L. Thompson /// @addtogroup CeedOperatorBackend
2067a982d89SJeremy L. Thompson /// @{
2077a982d89SJeremy L. Thompson 
2087a982d89SJeremy L. Thompson /**
2097a982d89SJeremy L. Thompson   @brief Get the Ceed associated with a CeedOperator
2107a982d89SJeremy L. Thompson 
2117a982d89SJeremy L. Thompson   @param op              CeedOperator
2127a982d89SJeremy L. Thompson   @param[out] ceed       Variable to store Ceed
2137a982d89SJeremy L. Thompson 
2147a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
2157a982d89SJeremy L. Thompson 
2167a982d89SJeremy L. Thompson   @ref Backend
2177a982d89SJeremy L. Thompson **/
2187a982d89SJeremy L. Thompson 
2197a982d89SJeremy L. Thompson int CeedOperatorGetCeed(CeedOperator op, Ceed *ceed) {
2207a982d89SJeremy L. Thompson   *ceed = op->ceed;
2217a982d89SJeremy L. Thompson   return 0;
2227a982d89SJeremy L. Thompson }
2237a982d89SJeremy L. Thompson 
2247a982d89SJeremy L. Thompson /**
2257a982d89SJeremy L. Thompson   @brief Get the number of elements associated with a CeedOperator
2267a982d89SJeremy L. Thompson 
2277a982d89SJeremy L. Thompson   @param op              CeedOperator
2287a982d89SJeremy L. Thompson   @param[out] numelem    Variable to store number of elements
2297a982d89SJeremy L. Thompson 
2307a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
2317a982d89SJeremy L. Thompson 
2327a982d89SJeremy L. Thompson   @ref Backend
2337a982d89SJeremy L. Thompson **/
2347a982d89SJeremy L. Thompson 
2357a982d89SJeremy L. Thompson int CeedOperatorGetNumElements(CeedOperator op, CeedInt *numelem) {
2367a982d89SJeremy L. Thompson   if (op->composite)
2377a982d89SJeremy L. Thompson     // LCOV_EXCL_START
2387a982d89SJeremy L. Thompson     return CeedError(op->ceed, 1, "Not defined for composite operator");
2397a982d89SJeremy L. Thompson   // LCOV_EXCL_STOP
2407a982d89SJeremy L. Thompson 
2417a982d89SJeremy L. Thompson   *numelem = op->numelements;
2427a982d89SJeremy L. Thompson   return 0;
2437a982d89SJeremy L. Thompson }
2447a982d89SJeremy L. Thompson 
2457a982d89SJeremy L. Thompson /**
2467a982d89SJeremy L. Thompson   @brief Get the number of quadrature points associated with a CeedOperator
2477a982d89SJeremy L. Thompson 
2487a982d89SJeremy L. Thompson   @param op              CeedOperator
2497a982d89SJeremy L. Thompson   @param[out] numqpts    Variable to store vector number of quadrature points
2507a982d89SJeremy L. Thompson 
2517a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
2527a982d89SJeremy L. Thompson 
2537a982d89SJeremy L. Thompson   @ref Backend
2547a982d89SJeremy L. Thompson **/
2557a982d89SJeremy L. Thompson 
2567a982d89SJeremy L. Thompson int CeedOperatorGetNumQuadraturePoints(CeedOperator op, CeedInt *numqpts) {
2577a982d89SJeremy L. Thompson   if (op->composite)
2587a982d89SJeremy L. Thompson     // LCOV_EXCL_START
2597a982d89SJeremy L. Thompson     return CeedError(op->ceed, 1, "Not defined for composite operator");
2607a982d89SJeremy L. Thompson   // LCOV_EXCL_STOP
2617a982d89SJeremy L. Thompson 
2627a982d89SJeremy L. Thompson   *numqpts = op->numqpoints;
2637a982d89SJeremy L. Thompson   return 0;
2647a982d89SJeremy L. Thompson }
2657a982d89SJeremy L. Thompson 
2667a982d89SJeremy L. Thompson /**
2677a982d89SJeremy L. Thompson   @brief Get the number of arguments associated with a CeedOperator
2687a982d89SJeremy L. Thompson 
2697a982d89SJeremy L. Thompson   @param op              CeedOperator
2707a982d89SJeremy L. Thompson   @param[out] numargs    Variable to store vector number of arguments
2717a982d89SJeremy L. Thompson 
2727a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
2737a982d89SJeremy L. Thompson 
2747a982d89SJeremy L. Thompson   @ref Backend
2757a982d89SJeremy L. Thompson **/
2767a982d89SJeremy L. Thompson 
2777a982d89SJeremy L. Thompson int CeedOperatorGetNumArgs(CeedOperator op, CeedInt *numargs) {
2787a982d89SJeremy L. Thompson   if (op->composite)
2797a982d89SJeremy L. Thompson     // LCOV_EXCL_START
2807a982d89SJeremy L. Thompson     return CeedError(op->ceed, 1, "Not defined for composite operators");
2817a982d89SJeremy L. Thompson   // LCOV_EXCL_STOP
2827a982d89SJeremy L. Thompson 
2837a982d89SJeremy L. Thompson   *numargs = op->nfields;
2847a982d89SJeremy L. Thompson   return 0;
2857a982d89SJeremy L. Thompson }
2867a982d89SJeremy L. Thompson 
2877a982d89SJeremy L. Thompson /**
2887a982d89SJeremy L. Thompson   @brief Get the setup status of a CeedOperator
2897a982d89SJeremy L. Thompson 
2907a982d89SJeremy L. Thompson   @param op             CeedOperator
2917a982d89SJeremy L. Thompson   @param[out] setupdone Variable to store setup status
2927a982d89SJeremy L. Thompson 
2937a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
2947a982d89SJeremy L. Thompson 
2957a982d89SJeremy L. Thompson   @ref Backend
2967a982d89SJeremy L. Thompson **/
2977a982d89SJeremy L. Thompson 
2987a982d89SJeremy L. Thompson int CeedOperatorGetSetupStatus(CeedOperator op, bool *setupdone) {
2997a982d89SJeremy L. Thompson   *setupdone = op->setupdone;
3007a982d89SJeremy L. Thompson   return 0;
3017a982d89SJeremy L. Thompson }
3027a982d89SJeremy L. Thompson 
3037a982d89SJeremy L. Thompson /**
3047a982d89SJeremy L. Thompson   @brief Get the QFunction associated with a CeedOperator
3057a982d89SJeremy L. Thompson 
3067a982d89SJeremy L. Thompson   @param op              CeedOperator
3077a982d89SJeremy L. Thompson   @param[out] qf         Variable to store QFunction
3087a982d89SJeremy L. Thompson 
3097a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
3107a982d89SJeremy L. Thompson 
3117a982d89SJeremy L. Thompson   @ref Backend
3127a982d89SJeremy L. Thompson **/
3137a982d89SJeremy L. Thompson 
3147a982d89SJeremy L. Thompson int CeedOperatorGetQFunction(CeedOperator op, CeedQFunction *qf) {
3157a982d89SJeremy L. Thompson   if (op->composite)
3167a982d89SJeremy L. Thompson     // LCOV_EXCL_START
3177a982d89SJeremy L. Thompson     return CeedError(op->ceed, 1, "Not defined for composite operator");
3187a982d89SJeremy L. Thompson   // LCOV_EXCL_STOP
3197a982d89SJeremy L. Thompson 
3207a982d89SJeremy L. Thompson   *qf = op->qf;
3217a982d89SJeremy L. Thompson   return 0;
3227a982d89SJeremy L. Thompson }
3237a982d89SJeremy L. Thompson 
3247a982d89SJeremy L. Thompson /**
3257a982d89SJeremy L. Thompson   @brief Get the number of suboperators associated with a CeedOperator
3267a982d89SJeremy L. Thompson 
3277a982d89SJeremy L. Thompson   @param op              CeedOperator
3287a982d89SJeremy L. Thompson   @param[out] numsub     Variable to store number of suboperators
3297a982d89SJeremy L. Thompson 
3307a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
3317a982d89SJeremy L. Thompson 
3327a982d89SJeremy L. Thompson   @ref Backend
3337a982d89SJeremy L. Thompson **/
3347a982d89SJeremy L. Thompson 
3357a982d89SJeremy L. Thompson int CeedOperatorGetNumSub(CeedOperator op, CeedInt *numsub) {
3367a982d89SJeremy L. Thompson   if (!op->composite)
3377a982d89SJeremy L. Thompson     // LCOV_EXCL_START
3387a982d89SJeremy L. Thompson     return CeedError(op->ceed, 1, "Not a composite operator");
3397a982d89SJeremy L. Thompson   // LCOV_EXCL_STOP
3407a982d89SJeremy L. Thompson 
3417a982d89SJeremy L. Thompson   *numsub = op->numsub;
3427a982d89SJeremy L. Thompson   return 0;
3437a982d89SJeremy L. Thompson }
3447a982d89SJeremy L. Thompson 
3457a982d89SJeremy L. Thompson /**
3467a982d89SJeremy L. Thompson   @brief Get the list of suboperators associated with a CeedOperator
3477a982d89SJeremy L. Thompson 
3487a982d89SJeremy L. Thompson   @param op                CeedOperator
3497a982d89SJeremy L. Thompson   @param[out] suboperators Variable to store list of suboperators
3507a982d89SJeremy L. Thompson 
3517a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
3527a982d89SJeremy L. Thompson 
3537a982d89SJeremy L. Thompson   @ref Backend
3547a982d89SJeremy L. Thompson **/
3557a982d89SJeremy L. Thompson 
3567a982d89SJeremy L. Thompson int CeedOperatorGetSubList(CeedOperator op, CeedOperator **suboperators) {
3577a982d89SJeremy L. Thompson   if (!op->composite)
3587a982d89SJeremy L. Thompson     // LCOV_EXCL_START
3597a982d89SJeremy L. Thompson     return CeedError(op->ceed, 1, "Not a composite operator");
3607a982d89SJeremy L. Thompson   // LCOV_EXCL_STOP
3617a982d89SJeremy L. Thompson 
3627a982d89SJeremy L. Thompson   *suboperators = op->suboperators;
3637a982d89SJeremy L. Thompson   return 0;
3647a982d89SJeremy L. Thompson }
3657a982d89SJeremy L. Thompson 
3667a982d89SJeremy L. Thompson /**
3677a982d89SJeremy L. Thompson   @brief Get the backend data of a CeedOperator
3687a982d89SJeremy L. Thompson 
3697a982d89SJeremy L. Thompson   @param op              CeedOperator
3707a982d89SJeremy L. Thompson   @param[out] data       Variable to store data
3717a982d89SJeremy L. Thompson 
3727a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
3737a982d89SJeremy L. Thompson 
3747a982d89SJeremy L. Thompson   @ref Backend
3757a982d89SJeremy L. Thompson **/
3767a982d89SJeremy L. Thompson 
3777a982d89SJeremy L. Thompson int CeedOperatorGetData(CeedOperator op, void **data) {
3787a982d89SJeremy L. Thompson   *data = op->data;
3797a982d89SJeremy L. Thompson   return 0;
3807a982d89SJeremy L. Thompson }
3817a982d89SJeremy L. Thompson 
3827a982d89SJeremy L. Thompson /**
3837a982d89SJeremy L. Thompson   @brief Set the backend data of a CeedOperator
3847a982d89SJeremy L. Thompson 
3857a982d89SJeremy L. Thompson   @param[out] op         CeedOperator
3867a982d89SJeremy L. Thompson   @param data            Data to set
3877a982d89SJeremy L. Thompson 
3887a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
3897a982d89SJeremy L. Thompson 
3907a982d89SJeremy L. Thompson   @ref Backend
3917a982d89SJeremy L. Thompson **/
3927a982d89SJeremy L. Thompson 
3937a982d89SJeremy L. Thompson int CeedOperatorSetData(CeedOperator op, void **data) {
3947a982d89SJeremy L. Thompson   op->data = *data;
3957a982d89SJeremy L. Thompson   return 0;
3967a982d89SJeremy L. Thompson }
3977a982d89SJeremy L. Thompson 
3987a982d89SJeremy L. Thompson /**
3997a982d89SJeremy L. Thompson   @brief Set the setup flag of a CeedOperator to True
4007a982d89SJeremy L. Thompson 
4017a982d89SJeremy L. Thompson   @param op              CeedOperator
4027a982d89SJeremy L. Thompson 
4037a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
4047a982d89SJeremy L. Thompson 
4057a982d89SJeremy L. Thompson   @ref Backend
4067a982d89SJeremy L. Thompson **/
4077a982d89SJeremy L. Thompson 
4087a982d89SJeremy L. Thompson int CeedOperatorSetSetupDone(CeedOperator op) {
4097a982d89SJeremy L. Thompson   op->setupdone = 1;
4107a982d89SJeremy L. Thompson   return 0;
4117a982d89SJeremy L. Thompson }
4127a982d89SJeremy L. Thompson 
4137a982d89SJeremy L. Thompson /**
4147a982d89SJeremy L. Thompson   @brief Get the CeedOperatorFields of a CeedOperator
4157a982d89SJeremy L. Thompson 
4167a982d89SJeremy L. Thompson   @param op                 CeedOperator
4177a982d89SJeremy L. Thompson   @param[out] inputfields   Variable to store inputfields
4187a982d89SJeremy L. Thompson   @param[out] outputfields  Variable to store outputfields
4197a982d89SJeremy L. Thompson 
4207a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
4217a982d89SJeremy L. Thompson 
4227a982d89SJeremy L. Thompson   @ref Backend
4237a982d89SJeremy L. Thompson **/
4247a982d89SJeremy L. Thompson 
4257a982d89SJeremy L. Thompson int CeedOperatorGetFields(CeedOperator op, CeedOperatorField **inputfields,
4267a982d89SJeremy L. Thompson                           CeedOperatorField **outputfields) {
4277a982d89SJeremy L. Thompson   if (op->composite)
4287a982d89SJeremy L. Thompson     // LCOV_EXCL_START
4297a982d89SJeremy L. Thompson     return CeedError(op->ceed, 1, "Not defined for composite operator");
4307a982d89SJeremy L. Thompson   // LCOV_EXCL_STOP
4317a982d89SJeremy L. Thompson 
4327a982d89SJeremy L. Thompson   if (inputfields) *inputfields = op->inputfields;
4337a982d89SJeremy L. Thompson   if (outputfields) *outputfields = op->outputfields;
4347a982d89SJeremy L. Thompson   return 0;
4357a982d89SJeremy L. Thompson }
4367a982d89SJeremy L. Thompson 
4377a982d89SJeremy L. Thompson /**
4387a982d89SJeremy L. Thompson   @brief Get the CeedElemRestriction of a CeedOperatorField
4397a982d89SJeremy L. Thompson 
4407a982d89SJeremy L. Thompson   @param opfield         CeedOperatorField
4417a982d89SJeremy L. Thompson   @param[out] rstr       Variable to store CeedElemRestriction
4427a982d89SJeremy L. Thompson 
4437a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
4447a982d89SJeremy L. Thompson 
4457a982d89SJeremy L. Thompson   @ref Backend
4467a982d89SJeremy L. Thompson **/
4477a982d89SJeremy L. Thompson 
4487a982d89SJeremy L. Thompson int CeedOperatorFieldGetElemRestriction(CeedOperatorField opfield,
4497a982d89SJeremy L. Thompson                                         CeedElemRestriction *rstr) {
4507a982d89SJeremy L. Thompson   *rstr = opfield->Erestrict;
4517a982d89SJeremy L. Thompson   return 0;
4527a982d89SJeremy L. Thompson }
4537a982d89SJeremy L. Thompson 
4547a982d89SJeremy L. Thompson /**
4557a982d89SJeremy L. Thompson   @brief Get the CeedBasis of a CeedOperatorField
4567a982d89SJeremy L. Thompson 
4577a982d89SJeremy L. Thompson   @param opfield         CeedOperatorField
4587a982d89SJeremy L. Thompson   @param[out] basis      Variable to store CeedBasis
4597a982d89SJeremy L. Thompson 
4607a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
4617a982d89SJeremy L. Thompson 
4627a982d89SJeremy L. Thompson   @ref Backend
4637a982d89SJeremy L. Thompson **/
4647a982d89SJeremy L. Thompson 
4657a982d89SJeremy L. Thompson int CeedOperatorFieldGetBasis(CeedOperatorField opfield, CeedBasis *basis) {
4667a982d89SJeremy L. Thompson   *basis = opfield->basis;
4677a982d89SJeremy L. Thompson   return 0;
4687a982d89SJeremy L. Thompson }
4697a982d89SJeremy L. Thompson 
4707a982d89SJeremy L. Thompson /**
4717a982d89SJeremy L. Thompson   @brief Get the CeedVector of a CeedOperatorField
4727a982d89SJeremy L. Thompson 
4737a982d89SJeremy L. Thompson   @param opfield         CeedOperatorField
4747a982d89SJeremy L. Thompson   @param[out] vec        Variable to store CeedVector
4757a982d89SJeremy L. Thompson 
4767a982d89SJeremy L. Thompson   @return An error code: 0 - success, otherwise - failure
4777a982d89SJeremy L. Thompson 
4787a982d89SJeremy L. Thompson   @ref Backend
4797a982d89SJeremy L. Thompson **/
4807a982d89SJeremy L. Thompson 
4817a982d89SJeremy L. Thompson int CeedOperatorFieldGetVector(CeedOperatorField opfield, CeedVector *vec) {
4827a982d89SJeremy L. Thompson   *vec = opfield->vec;
4837a982d89SJeremy L. Thompson   return 0;
4847a982d89SJeremy L. Thompson }
4857a982d89SJeremy L. Thompson 
4867a982d89SJeremy L. Thompson /// @}
4877a982d89SJeremy L. Thompson 
4887a982d89SJeremy L. Thompson /// ----------------------------------------------------------------------------
4897a982d89SJeremy L. Thompson /// CeedOperator Public API
4907a982d89SJeremy L. Thompson /// ----------------------------------------------------------------------------
4917a982d89SJeremy L. Thompson /// @addtogroup CeedOperatorUser
492dfdf5a53Sjeremylt /// @{
493d7b241e6Sjeremylt 
494d7b241e6Sjeremylt /**
4950219ea01SJeremy L Thompson   @brief Create a CeedOperator and associate a CeedQFunction. A CeedBasis and
4960219ea01SJeremy L Thompson            CeedElemRestriction can be associated with CeedQFunction fields with
4970219ea01SJeremy L Thompson            \ref CeedOperatorSetField.
498d7b241e6Sjeremylt 
499b11c1e72Sjeremylt   @param ceed    A Ceed object where the CeedOperator will be created
500d7b241e6Sjeremylt   @param qf      QFunction defining the action of the operator at quadrature points
50134138859Sjeremylt   @param dqf     QFunction defining the action of the Jacobian of @a qf (or
5024cc79fe7SJed Brown                    @ref CEED_QFUNCTION_NONE)
503d7b241e6Sjeremylt   @param dqfT    QFunction defining the action of the transpose of the Jacobian
5044cc79fe7SJed Brown                    of @a qf (or @ref CEED_QFUNCTION_NONE)
505b11c1e72Sjeremylt   @param[out] op Address of the variable where the newly created
506b11c1e72Sjeremylt                      CeedOperator will be stored
507b11c1e72Sjeremylt 
508b11c1e72Sjeremylt   @return An error code: 0 - success, otherwise - failure
509dfdf5a53Sjeremylt 
5107a982d89SJeremy L. Thompson   @ref User
511d7b241e6Sjeremylt  */
512d7b241e6Sjeremylt int CeedOperatorCreate(Ceed ceed, CeedQFunction qf, CeedQFunction dqf,
513d7b241e6Sjeremylt                        CeedQFunction dqfT, CeedOperator *op) {
514d7b241e6Sjeremylt   int ierr;
515d7b241e6Sjeremylt 
5165fe0d4faSjeremylt   if (!ceed->OperatorCreate) {
5175fe0d4faSjeremylt     Ceed delegate;
518aefd8378Sjeremylt     ierr = CeedGetObjectDelegate(ceed, &delegate, "Operator"); CeedChk(ierr);
5195fe0d4faSjeremylt 
5205fe0d4faSjeremylt     if (!delegate)
521c042f62fSJeremy L Thompson       // LCOV_EXCL_START
5225fe0d4faSjeremylt       return CeedError(ceed, 1, "Backend does not support OperatorCreate");
523c042f62fSJeremy L Thompson     // LCOV_EXCL_STOP
5245fe0d4faSjeremylt 
5255fe0d4faSjeremylt     ierr = CeedOperatorCreate(delegate, qf, dqf, dqfT, op); CeedChk(ierr);
5265fe0d4faSjeremylt     return 0;
5275fe0d4faSjeremylt   }
5285fe0d4faSjeremylt 
529*b3b7035fSJeremy L Thompson   if (!qf || qf == CEED_QFUNCTION_NONE)
530*b3b7035fSJeremy L Thompson     // LCOV_EXCL_START
531*b3b7035fSJeremy L Thompson     return CeedError(ceed, 1, "Operator must have a valid QFunction.");
532*b3b7035fSJeremy L Thompson   // LCOV_EXCL_STOP
533d7b241e6Sjeremylt   ierr = CeedCalloc(1, op); CeedChk(ierr);
534d7b241e6Sjeremylt   (*op)->ceed = ceed;
535d7b241e6Sjeremylt   ceed->refcount++;
536d7b241e6Sjeremylt   (*op)->refcount = 1;
537d7b241e6Sjeremylt   (*op)->qf = qf;
538d7b241e6Sjeremylt   qf->refcount++;
539442e7f0bSjeremylt   if (dqf && dqf != CEED_QFUNCTION_NONE) {
540d7b241e6Sjeremylt     (*op)->dqf = dqf;
541442e7f0bSjeremylt     dqf->refcount++;
542442e7f0bSjeremylt   }
543442e7f0bSjeremylt   if (dqfT && dqfT != CEED_QFUNCTION_NONE) {
544d7b241e6Sjeremylt     (*op)->dqfT = dqfT;
545442e7f0bSjeremylt     dqfT->refcount++;
546442e7f0bSjeremylt   }
547fe2413ffSjeremylt   ierr = CeedCalloc(16, &(*op)->inputfields); CeedChk(ierr);
548fe2413ffSjeremylt   ierr = CeedCalloc(16, &(*op)->outputfields); CeedChk(ierr);
549d7b241e6Sjeremylt   ierr = ceed->OperatorCreate(*op); CeedChk(ierr);
550d7b241e6Sjeremylt   return 0;
551d7b241e6Sjeremylt }
552d7b241e6Sjeremylt 
553d7b241e6Sjeremylt /**
55452d6035fSJeremy L Thompson   @brief Create an operator that composes the action of several operators
55552d6035fSJeremy L Thompson 
55652d6035fSJeremy L Thompson   @param ceed    A Ceed object where the CeedOperator will be created
55752d6035fSJeremy L Thompson   @param[out] op Address of the variable where the newly created
55852d6035fSJeremy L Thompson                      Composite CeedOperator will be stored
55952d6035fSJeremy L Thompson 
56052d6035fSJeremy L Thompson   @return An error code: 0 - success, otherwise - failure
56152d6035fSJeremy L Thompson 
5627a982d89SJeremy L. Thompson   @ref User
56352d6035fSJeremy L Thompson  */
56452d6035fSJeremy L Thompson int CeedCompositeOperatorCreate(Ceed ceed, CeedOperator *op) {
56552d6035fSJeremy L Thompson   int ierr;
56652d6035fSJeremy L Thompson 
56752d6035fSJeremy L Thompson   if (!ceed->CompositeOperatorCreate) {
56852d6035fSJeremy L Thompson     Ceed delegate;
569aefd8378Sjeremylt     ierr = CeedGetObjectDelegate(ceed, &delegate, "Operator"); CeedChk(ierr);
57052d6035fSJeremy L Thompson 
571250756a7Sjeremylt     if (delegate) {
57252d6035fSJeremy L Thompson       ierr = CeedCompositeOperatorCreate(delegate, op); CeedChk(ierr);
57352d6035fSJeremy L Thompson       return 0;
57452d6035fSJeremy L Thompson     }
575250756a7Sjeremylt   }
57652d6035fSJeremy L Thompson 
57752d6035fSJeremy L Thompson   ierr = CeedCalloc(1, op); CeedChk(ierr);
57852d6035fSJeremy L Thompson   (*op)->ceed = ceed;
57952d6035fSJeremy L Thompson   ceed->refcount++;
58052d6035fSJeremy L Thompson   (*op)->composite = true;
58152d6035fSJeremy L Thompson   ierr = CeedCalloc(16, &(*op)->suboperators); CeedChk(ierr);
582250756a7Sjeremylt 
583250756a7Sjeremylt   if (ceed->CompositeOperatorCreate) {
58452d6035fSJeremy L Thompson     ierr = ceed->CompositeOperatorCreate(*op); CeedChk(ierr);
585250756a7Sjeremylt   }
58652d6035fSJeremy L Thompson   return 0;
58752d6035fSJeremy L Thompson }
58852d6035fSJeremy L Thompson 
58952d6035fSJeremy L Thompson /**
590b11c1e72Sjeremylt   @brief Provide a field to a CeedOperator for use by its CeedQFunction
591d7b241e6Sjeremylt 
592d7b241e6Sjeremylt   This function is used to specify both active and passive fields to a
593d7b241e6Sjeremylt   CeedOperator.  For passive fields, a vector @arg v must be provided.  Passive
594d7b241e6Sjeremylt   fields can inputs or outputs (updated in-place when operator is applied).
595d7b241e6Sjeremylt 
596d7b241e6Sjeremylt   Active fields must be specified using this function, but their data (in a
597d7b241e6Sjeremylt   CeedVector) is passed in CeedOperatorApply().  There can be at most one active
598d7b241e6Sjeremylt   input and at most one active output.
599d7b241e6Sjeremylt 
6008c91a0c9SJeremy L Thompson   @param op         CeedOperator on which to provide the field
6018795c945Sjeremylt   @param fieldname  Name of the field (to be matched with the name used by
6028795c945Sjeremylt                       CeedQFunction)
603b11c1e72Sjeremylt   @param r          CeedElemRestriction
6044cc79fe7SJed Brown   @param b          CeedBasis in which the field resides or @ref CEED_BASIS_COLLOCATED
605b11c1e72Sjeremylt                       if collocated with quadrature points
6064cc79fe7SJed Brown   @param v          CeedVector to be used by CeedOperator or @ref CEED_VECTOR_ACTIVE
6074cc79fe7SJed Brown                       if field is active or @ref CEED_VECTOR_NONE if using
6084cc79fe7SJed Brown                       @ref CEED_EVAL_WEIGHT in the QFunction
609b11c1e72Sjeremylt 
610b11c1e72Sjeremylt   @return An error code: 0 - success, otherwise - failure
611dfdf5a53Sjeremylt 
6127a982d89SJeremy L. Thompson   @ref User
613b11c1e72Sjeremylt **/
614d7b241e6Sjeremylt int CeedOperatorSetField(CeedOperator op, const char *fieldname,
615a8d32208Sjeremylt                          CeedElemRestriction r, CeedBasis b, CeedVector v) {
616d7b241e6Sjeremylt   int ierr;
61752d6035fSJeremy L Thompson   if (op->composite)
618c042f62fSJeremy L Thompson     // LCOV_EXCL_START
61952d6035fSJeremy L Thompson     return CeedError(op->ceed, 1, "Cannot add field to composite operator.");
620c042f62fSJeremy L Thompson   // LCOV_EXCL_STOP
6218b067b84SJed Brown   if (!r)
622c042f62fSJeremy L Thompson     // LCOV_EXCL_START
623c042f62fSJeremy L Thompson     return CeedError(op->ceed, 1,
624c042f62fSJeremy L Thompson                      "ElemRestriction r for field \"%s\" must be non-NULL.",
625c042f62fSJeremy L Thompson                      fieldname);
626c042f62fSJeremy L Thompson   // LCOV_EXCL_STOP
6278b067b84SJed Brown   if (!b)
628c042f62fSJeremy L Thompson     // LCOV_EXCL_START
629c042f62fSJeremy L Thompson     return CeedError(op->ceed, 1, "Basis b for field \"%s\" must be non-NULL.",
630c042f62fSJeremy L Thompson                      fieldname);
631c042f62fSJeremy L Thompson   // LCOV_EXCL_STOP
6328b067b84SJed Brown   if (!v)
633c042f62fSJeremy L Thompson     // LCOV_EXCL_START
634c042f62fSJeremy L Thompson     return CeedError(op->ceed, 1, "Vector v for field \"%s\" must be non-NULL.",
635c042f62fSJeremy L Thompson                      fieldname);
636c042f62fSJeremy L Thompson   // LCOV_EXCL_STOP
63752d6035fSJeremy L Thompson 
638d7b241e6Sjeremylt   CeedInt numelements;
639d7b241e6Sjeremylt   ierr = CeedElemRestrictionGetNumElements(r, &numelements); CeedChk(ierr);
64015910d16Sjeremylt   if (r != CEED_ELEMRESTRICTION_NONE && op->hasrestriction &&
64115910d16Sjeremylt       op->numelements != numelements)
642c042f62fSJeremy L Thompson     // LCOV_EXCL_START
643d7b241e6Sjeremylt     return CeedError(op->ceed, 1,
6441d102b48SJeremy L Thompson                      "ElemRestriction with %d elements incompatible with prior "
6451d102b48SJeremy L Thompson                      "%d elements", numelements, op->numelements);
646c042f62fSJeremy L Thompson   // LCOV_EXCL_STOP
64715910d16Sjeremylt   if (r != CEED_ELEMRESTRICTION_NONE) {
648d7b241e6Sjeremylt     op->numelements = numelements;
6492cb0afc5Sjeremylt     op->hasrestriction = true; // Restriction set, but numelements may be 0
65015910d16Sjeremylt   }
651d7b241e6Sjeremylt 
652783c99b3SValeria Barra   if (b != CEED_BASIS_COLLOCATED) {
653d7b241e6Sjeremylt     CeedInt numqpoints;
654d7b241e6Sjeremylt     ierr = CeedBasisGetNumQuadraturePoints(b, &numqpoints); CeedChk(ierr);
655d7b241e6Sjeremylt     if (op->numqpoints && op->numqpoints != numqpoints)
656c042f62fSJeremy L Thompson       // LCOV_EXCL_START
6571d102b48SJeremy L Thompson       return CeedError(op->ceed, 1, "Basis with %d quadrature points "
6581d102b48SJeremy L Thompson                        "incompatible with prior %d points", numqpoints,
6591d102b48SJeremy L Thompson                        op->numqpoints);
660c042f62fSJeremy L Thompson     // LCOV_EXCL_STOP
661d7b241e6Sjeremylt     op->numqpoints = numqpoints;
662d7b241e6Sjeremylt   }
66315910d16Sjeremylt   CeedQFunctionField qfield;
664d1bcdac9Sjeremylt   CeedOperatorField *ofield;
665d7b241e6Sjeremylt   for (CeedInt i=0; i<op->qf->numinputfields; i++) {
666fe2413ffSjeremylt     if (!strcmp(fieldname, (*op->qf->inputfields[i]).fieldname)) {
66715910d16Sjeremylt       qfield = op->qf->inputfields[i];
668d7b241e6Sjeremylt       ofield = &op->inputfields[i];
669d7b241e6Sjeremylt       goto found;
670d7b241e6Sjeremylt     }
671d7b241e6Sjeremylt   }
672d7b241e6Sjeremylt   for (CeedInt i=0; i<op->qf->numoutputfields; i++) {
673fe2413ffSjeremylt     if (!strcmp(fieldname, (*op->qf->outputfields[i]).fieldname)) {
67415910d16Sjeremylt       qfield = op->qf->inputfields[i];
675d7b241e6Sjeremylt       ofield = &op->outputfields[i];
676d7b241e6Sjeremylt       goto found;
677d7b241e6Sjeremylt     }
678d7b241e6Sjeremylt   }
679c042f62fSJeremy L Thompson   // LCOV_EXCL_START
680d7b241e6Sjeremylt   return CeedError(op->ceed, 1, "QFunction has no knowledge of field '%s'",
681d7b241e6Sjeremylt                    fieldname);
682c042f62fSJeremy L Thompson   // LCOV_EXCL_STOP
683d7b241e6Sjeremylt found:
68415910d16Sjeremylt   if (r == CEED_ELEMRESTRICTION_NONE && qfield->emode != CEED_EVAL_WEIGHT)
685b0d62198Sjeremylt     // LCOV_EXCL_START
68615910d16Sjeremylt     return CeedError(op->ceed, 1, "CEED_ELEMRESTRICTION_NONE can only be used "
68715910d16Sjeremylt                      "for a field with eval mode CEED_EVAL_WEIGHT");
688b0d62198Sjeremylt   // LCOV_EXCL_STOP
689fe2413ffSjeremylt   ierr = CeedCalloc(1, ofield); CeedChk(ierr);
690fe2413ffSjeremylt   (*ofield)->Erestrict = r;
69171352a93Sjeremylt   r->refcount += 1;
692fe2413ffSjeremylt   (*ofield)->basis = b;
69371352a93Sjeremylt   if (b != CEED_BASIS_COLLOCATED)
69471352a93Sjeremylt     b->refcount += 1;
695fe2413ffSjeremylt   (*ofield)->vec = v;
69671352a93Sjeremylt   if (v != CEED_VECTOR_ACTIVE && v != CEED_VECTOR_NONE)
69771352a93Sjeremylt     v->refcount += 1;
698d7b241e6Sjeremylt   op->nfields += 1;
699d7b241e6Sjeremylt   return 0;
700d7b241e6Sjeremylt }
701d7b241e6Sjeremylt 
702d7b241e6Sjeremylt /**
70352d6035fSJeremy L Thompson   @brief Add a sub-operator to a composite CeedOperator
704288c0443SJeremy L Thompson 
70534138859Sjeremylt   @param[out] compositeop Composite CeedOperator
70634138859Sjeremylt   @param      subop       Sub-operator CeedOperator
70752d6035fSJeremy L Thompson 
70852d6035fSJeremy L Thompson   @return An error code: 0 - success, otherwise - failure
70952d6035fSJeremy L Thompson 
7107a982d89SJeremy L. Thompson   @ref User
71152d6035fSJeremy L Thompson  */
712692c2638Sjeremylt int CeedCompositeOperatorAddSub(CeedOperator compositeop, CeedOperator subop) {
71352d6035fSJeremy L Thompson   if (!compositeop->composite)
714c042f62fSJeremy L Thompson     // LCOV_EXCL_START
7151d102b48SJeremy L Thompson     return CeedError(compositeop->ceed, 1, "CeedOperator is not a composite "
7161d102b48SJeremy L Thompson                      "operator");
717c042f62fSJeremy L Thompson   // LCOV_EXCL_STOP
71852d6035fSJeremy L Thompson 
71952d6035fSJeremy L Thompson   if (compositeop->numsub == CEED_COMPOSITE_MAX)
720c042f62fSJeremy L Thompson     // LCOV_EXCL_START
7211d102b48SJeremy L Thompson     return CeedError(compositeop->ceed, 1, "Cannot add additional suboperators");
722c042f62fSJeremy L Thompson   // LCOV_EXCL_STOP
72352d6035fSJeremy L Thompson 
72452d6035fSJeremy L Thompson   compositeop->suboperators[compositeop->numsub] = subop;
72552d6035fSJeremy L Thompson   subop->refcount++;
72652d6035fSJeremy L Thompson   compositeop->numsub++;
72752d6035fSJeremy L Thompson   return 0;
72852d6035fSJeremy L Thompson }
72952d6035fSJeremy L Thompson 
73052d6035fSJeremy L Thompson /**
7311d102b48SJeremy L Thompson   @brief Assemble a linear CeedQFunction associated with a CeedOperator
7321d102b48SJeremy L Thompson 
7331d102b48SJeremy L Thompson   This returns a CeedVector containing a matrix at each quadrature point
7341d102b48SJeremy L Thompson     providing the action of the CeedQFunction associated with the CeedOperator.
7351d102b48SJeremy L Thompson     The vector 'assembled' is of shape
7361d102b48SJeremy L Thompson       [num_elements, num_input_fields, num_output_fields, num_quad_points]
7371d102b48SJeremy L Thompson     and contains column-major matrices representing the action of the
7381d102b48SJeremy L Thompson     CeedQFunction for a corresponding quadrature point on an element. Inputs and
7391d102b48SJeremy L Thompson     outputs are in the order provided by the user when adding CeedOperator fields.
7404cc79fe7SJed Brown     For example, a CeedQFunction with inputs 'u' and 'gradu' and outputs 'gradv' and
7411d102b48SJeremy L Thompson     'v', provided in that order, would result in an assembled QFunction that
7421d102b48SJeremy L Thompson     consists of (1 + dim) x (dim + 1) matrices at each quadrature point acting
7431d102b48SJeremy L Thompson     on the input [u, du_0, du_1] and producing the output [dv_0, dv_1, v].
7441d102b48SJeremy L Thompson 
7451d102b48SJeremy L Thompson   @param op             CeedOperator to assemble CeedQFunction
7461d102b48SJeremy L Thompson   @param[out] assembled CeedVector to store assembled CeedQFunction at
7471d102b48SJeremy L Thompson                           quadrature points
7481d102b48SJeremy L Thompson   @param[out] rstr      CeedElemRestriction for CeedVector containing assembled
7491d102b48SJeremy L Thompson                           CeedQFunction
7501d102b48SJeremy L Thompson   @param request        Address of CeedRequest for non-blocking completion, else
7514cc79fe7SJed Brown                           @ref CEED_REQUEST_IMMEDIATE
7521d102b48SJeremy L Thompson 
7531d102b48SJeremy L Thompson   @return An error code: 0 - success, otherwise - failure
7541d102b48SJeremy L Thompson 
7557a982d89SJeremy L. Thompson   @ref User
7561d102b48SJeremy L Thompson **/
7571d102b48SJeremy L Thompson int CeedOperatorAssembleLinearQFunction(CeedOperator op, CeedVector *assembled,
7587f823360Sjeremylt                                         CeedElemRestriction *rstr,
7597f823360Sjeremylt                                         CeedRequest *request) {
7601d102b48SJeremy L Thompson   int ierr;
7611d102b48SJeremy L Thompson   Ceed ceed = op->ceed;
762250756a7Sjeremylt   ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr);
7631d102b48SJeremy L Thompson 
7647172caadSjeremylt   // Backend version
7655107b09fSJeremy L Thompson   if (op->AssembleLinearQFunction) {
7661d102b48SJeremy L Thompson     ierr = op->AssembleLinearQFunction(op, assembled, rstr, request);
7671d102b48SJeremy L Thompson     CeedChk(ierr);
7685107b09fSJeremy L Thompson   } else {
7695107b09fSJeremy L Thompson     // Fallback to reference Ceed
7705107b09fSJeremy L Thompson     if (!op->opfallback) {
7715107b09fSJeremy L Thompson       ierr = CeedOperatorCreateFallback(op); CeedChk(ierr);
7725107b09fSJeremy L Thompson     }
7735107b09fSJeremy L Thompson     // Assemble
7745107b09fSJeremy L Thompson     ierr = op->opfallback->AssembleLinearQFunction(op->opfallback, assembled,
7755107b09fSJeremy L Thompson            rstr, request); CeedChk(ierr);
7765107b09fSJeremy L Thompson   }
777250756a7Sjeremylt 
7781d102b48SJeremy L Thompson   return 0;
7791d102b48SJeremy L Thompson }
7801d102b48SJeremy L Thompson 
7811d102b48SJeremy L Thompson /**
782b7ec98d8SJeremy L Thompson   @brief Assemble the diagonal of a square linear Operator
783b7ec98d8SJeremy L Thompson 
784b7ec98d8SJeremy L Thompson   This returns a CeedVector containing the diagonal of a linear CeedOperator.
785b7ec98d8SJeremy L Thompson 
786b7ec98d8SJeremy L Thompson   @param op             CeedOperator to assemble CeedQFunction
787b7ec98d8SJeremy L Thompson   @param[out] assembled CeedVector to store assembled CeedOperator diagonal
788b7ec98d8SJeremy L Thompson   @param request        Address of CeedRequest for non-blocking completion, else
7894cc79fe7SJed Brown                           @ref CEED_REQUEST_IMMEDIATE
790b7ec98d8SJeremy L Thompson 
791b7ec98d8SJeremy L Thompson   @return An error code: 0 - success, otherwise - failure
792b7ec98d8SJeremy L Thompson 
7937a982d89SJeremy L. Thompson   @ref User
794b7ec98d8SJeremy L Thompson **/
7957f823360Sjeremylt int CeedOperatorAssembleLinearDiagonal(CeedOperator op, CeedVector *assembled,
7967f823360Sjeremylt                                        CeedRequest *request) {
797b7ec98d8SJeremy L Thompson   int ierr;
798b7ec98d8SJeremy L Thompson   Ceed ceed = op->ceed;
799250756a7Sjeremylt   ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr);
800b7ec98d8SJeremy L Thompson 
801b7ec98d8SJeremy L Thompson   // Use backend version, if available
802b7ec98d8SJeremy L Thompson   if (op->AssembleLinearDiagonal) {
803b7ec98d8SJeremy L Thompson     ierr = op->AssembleLinearDiagonal(op, assembled, request); CeedChk(ierr);
8047172caadSjeremylt   } else {
8057172caadSjeremylt     // Fallback to reference Ceed
8067172caadSjeremylt     if (!op->opfallback) {
8077172caadSjeremylt       ierr = CeedOperatorCreateFallback(op); CeedChk(ierr);
808b7ec98d8SJeremy L Thompson     }
8097172caadSjeremylt     // Assemble
8107172caadSjeremylt     ierr = op->opfallback->AssembleLinearDiagonal(op->opfallback, assembled,
8117172caadSjeremylt            request); CeedChk(ierr);
812b7ec98d8SJeremy L Thompson   }
813b7ec98d8SJeremy L Thompson 
814b7ec98d8SJeremy L Thompson   return 0;
815b7ec98d8SJeremy L Thompson }
816b7ec98d8SJeremy L Thompson 
817b7ec98d8SJeremy L Thompson /**
8183bd813ffSjeremylt   @brief Build a FDM based approximate inverse for each element for a
8193bd813ffSjeremylt            CeedOperator
8203bd813ffSjeremylt 
8213bd813ffSjeremylt   This returns a CeedOperator and CeedVector to apply a Fast Diagonalization
8223bd813ffSjeremylt     Method based approximate inverse. This function obtains the simultaneous
8233bd813ffSjeremylt     diagonalization for the 1D mass and Laplacian operators,
8243bd813ffSjeremylt       M = V^T V, K = V^T S V.
8253bd813ffSjeremylt     The assembled QFunction is used to modify the eigenvalues from simultaneous
8263bd813ffSjeremylt     diagonalization and obtain an approximate inverse of the form
8273bd813ffSjeremylt       V^T S^hat V. The CeedOperator must be linear and non-composite. The
8283bd813ffSjeremylt     associated CeedQFunction must therefore also be linear.
8293bd813ffSjeremylt 
8303bd813ffSjeremylt   @param op             CeedOperator to create element inverses
8313bd813ffSjeremylt   @param[out] fdminv    CeedOperator to apply the action of a FDM based inverse
8323bd813ffSjeremylt                           for each element
8333bd813ffSjeremylt   @param request        Address of CeedRequest for non-blocking completion, else
8344cc79fe7SJed Brown                           @ref CEED_REQUEST_IMMEDIATE
8353bd813ffSjeremylt 
8363bd813ffSjeremylt   @return An error code: 0 - success, otherwise - failure
8373bd813ffSjeremylt 
8387a982d89SJeremy L. Thompson   @ref Backend
8393bd813ffSjeremylt **/
8403bd813ffSjeremylt int CeedOperatorCreateFDMElementInverse(CeedOperator op, CeedOperator *fdminv,
8413bd813ffSjeremylt                                         CeedRequest *request) {
8423bd813ffSjeremylt   int ierr;
8433bd813ffSjeremylt   Ceed ceed = op->ceed;
844713f43c3Sjeremylt   ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr);
8453bd813ffSjeremylt 
846713f43c3Sjeremylt   // Use backend version, if available
847713f43c3Sjeremylt   if (op->CreateFDMElementInverse) {
848713f43c3Sjeremylt     ierr = op->CreateFDMElementInverse(op, fdminv, request); CeedChk(ierr);
849713f43c3Sjeremylt   } else {
850713f43c3Sjeremylt     // Fallback to reference Ceed
851713f43c3Sjeremylt     if (!op->opfallback) {
852713f43c3Sjeremylt       ierr = CeedOperatorCreateFallback(op); CeedChk(ierr);
8533bd813ffSjeremylt     }
854713f43c3Sjeremylt     // Assemble
855713f43c3Sjeremylt     ierr = op->opfallback->CreateFDMElementInverse(op->opfallback, fdminv,
8563bd813ffSjeremylt            request); CeedChk(ierr);
8573bd813ffSjeremylt   }
8583bd813ffSjeremylt 
8593bd813ffSjeremylt   return 0;
8603bd813ffSjeremylt }
8613bd813ffSjeremylt 
8627a982d89SJeremy L. Thompson /**
8637a982d89SJeremy L. Thompson   @brief View a CeedOperator
8647a982d89SJeremy L. Thompson 
8657a982d89SJeremy L. Thompson   @param[in] op     CeedOperator to view
8667a982d89SJeremy L. Thompson   @param[in] stream Stream to write; typically stdout/stderr or a file
8677a982d89SJeremy L. Thompson 
8687a982d89SJeremy L. Thompson   @return Error code: 0 - success, otherwise - failure
8697a982d89SJeremy L. Thompson 
8707a982d89SJeremy L. Thompson   @ref User
8717a982d89SJeremy L. Thompson **/
8727a982d89SJeremy L. Thompson int CeedOperatorView(CeedOperator op, FILE *stream) {
8737a982d89SJeremy L. Thompson   int ierr;
8747a982d89SJeremy L. Thompson 
8757a982d89SJeremy L. Thompson   if (op->composite) {
8767a982d89SJeremy L. Thompson     fprintf(stream, "Composite CeedOperator\n");
8777a982d89SJeremy L. Thompson 
8787a982d89SJeremy L. Thompson     for (CeedInt i=0; i<op->numsub; i++) {
8797a982d89SJeremy L. Thompson       fprintf(stream, "  SubOperator [%d]:\n", i);
8807a982d89SJeremy L. Thompson       ierr = CeedOperatorSingleView(op->suboperators[i], 1, stream);
8817a982d89SJeremy L. Thompson       CeedChk(ierr);
8827a982d89SJeremy L. Thompson     }
8837a982d89SJeremy L. Thompson   } else {
8847a982d89SJeremy L. Thompson     fprintf(stream, "CeedOperator\n");
8857a982d89SJeremy L. Thompson     ierr = CeedOperatorSingleView(op, 0, stream); CeedChk(ierr);
8867a982d89SJeremy L. Thompson   }
8877a982d89SJeremy L. Thompson 
8887a982d89SJeremy L. Thompson   return 0;
8897a982d89SJeremy L. Thompson }
8903bd813ffSjeremylt 
8913bd813ffSjeremylt /**
8923bd813ffSjeremylt   @brief Apply CeedOperator to a vector
893d7b241e6Sjeremylt 
894d7b241e6Sjeremylt   This computes the action of the operator on the specified (active) input,
895d7b241e6Sjeremylt   yielding its (active) output.  All inputs and outputs must be specified using
896d7b241e6Sjeremylt   CeedOperatorSetField().
897d7b241e6Sjeremylt 
898d7b241e6Sjeremylt   @param op        CeedOperator to apply
8994cc79fe7SJed Brown   @param[in] in    CeedVector containing input state or @ref CEED_VECTOR_NONE if
90034138859Sjeremylt                   there are no active inputs
901b11c1e72Sjeremylt   @param[out] out  CeedVector to store result of applying operator (must be
9024cc79fe7SJed Brown                      distinct from @a in) or @ref CEED_VECTOR_NONE if there are no
90334138859Sjeremylt                      active outputs
904d7b241e6Sjeremylt   @param request   Address of CeedRequest for non-blocking completion, else
9054cc79fe7SJed Brown                      @ref CEED_REQUEST_IMMEDIATE
906b11c1e72Sjeremylt 
907b11c1e72Sjeremylt   @return An error code: 0 - success, otherwise - failure
908dfdf5a53Sjeremylt 
9097a982d89SJeremy L. Thompson   @ref User
910b11c1e72Sjeremylt **/
911692c2638Sjeremylt int CeedOperatorApply(CeedOperator op, CeedVector in, CeedVector out,
912692c2638Sjeremylt                       CeedRequest *request) {
913d7b241e6Sjeremylt   int ierr;
914d7b241e6Sjeremylt   Ceed ceed = op->ceed;
915250756a7Sjeremylt   ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr);
916d7b241e6Sjeremylt 
917250756a7Sjeremylt   if (op->numelements)  {
918250756a7Sjeremylt     // Standard Operator
919cae8b89aSjeremylt     if (op->Apply) {
920250756a7Sjeremylt       ierr = op->Apply(op, in, out, request); CeedChk(ierr);
921cae8b89aSjeremylt     } else {
922cae8b89aSjeremylt       // Zero all output vectors
923250756a7Sjeremylt       CeedQFunction qf = op->qf;
924cae8b89aSjeremylt       for (CeedInt i=0; i<qf->numoutputfields; i++) {
925cae8b89aSjeremylt         CeedVector vec = op->outputfields[i]->vec;
926cae8b89aSjeremylt         if (vec == CEED_VECTOR_ACTIVE)
927cae8b89aSjeremylt           vec = out;
928cae8b89aSjeremylt         if (vec != CEED_VECTOR_NONE) {
929cae8b89aSjeremylt           ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr);
930cae8b89aSjeremylt         }
931cae8b89aSjeremylt       }
932250756a7Sjeremylt       // Apply
933250756a7Sjeremylt       ierr = op->ApplyAdd(op, in, out, request); CeedChk(ierr);
934250756a7Sjeremylt     }
935250756a7Sjeremylt   } else if (op->composite) {
936250756a7Sjeremylt     // Composite Operator
937250756a7Sjeremylt     if (op->ApplyComposite) {
938250756a7Sjeremylt       ierr = op->ApplyComposite(op, in, out, request); CeedChk(ierr);
939250756a7Sjeremylt     } else {
940250756a7Sjeremylt       CeedInt numsub;
941250756a7Sjeremylt       ierr = CeedOperatorGetNumSub(op, &numsub); CeedChk(ierr);
942250756a7Sjeremylt       CeedOperator *suboperators;
943250756a7Sjeremylt       ierr = CeedOperatorGetSubList(op, &suboperators); CeedChk(ierr);
944250756a7Sjeremylt 
945250756a7Sjeremylt       // Zero all output vectors
946250756a7Sjeremylt       if (out != CEED_VECTOR_NONE) {
947cae8b89aSjeremylt         ierr = CeedVectorSetValue(out, 0.0); CeedChk(ierr);
948cae8b89aSjeremylt       }
949250756a7Sjeremylt       for (CeedInt i=0; i<numsub; i++) {
950250756a7Sjeremylt         for (CeedInt j=0; j<suboperators[i]->qf->numoutputfields; j++) {
951250756a7Sjeremylt           CeedVector vec = suboperators[i]->outputfields[j]->vec;
952250756a7Sjeremylt           if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) {
953250756a7Sjeremylt             ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr);
954250756a7Sjeremylt           }
955250756a7Sjeremylt         }
956250756a7Sjeremylt       }
957250756a7Sjeremylt       // Apply
958250756a7Sjeremylt       for (CeedInt i=0; i<op->numsub; i++) {
959250756a7Sjeremylt         ierr = CeedOperatorApplyAdd(op->suboperators[i], in, out, request);
960cae8b89aSjeremylt         CeedChk(ierr);
961cae8b89aSjeremylt       }
962cae8b89aSjeremylt     }
963250756a7Sjeremylt   }
964250756a7Sjeremylt 
965cae8b89aSjeremylt   return 0;
966cae8b89aSjeremylt }
967cae8b89aSjeremylt 
968cae8b89aSjeremylt /**
969cae8b89aSjeremylt   @brief Apply CeedOperator to a vector and add result to output vector
970cae8b89aSjeremylt 
971cae8b89aSjeremylt   This computes the action of the operator on the specified (active) input,
972cae8b89aSjeremylt   yielding its (active) output.  All inputs and outputs must be specified using
973cae8b89aSjeremylt   CeedOperatorSetField().
974cae8b89aSjeremylt 
975cae8b89aSjeremylt   @param op        CeedOperator to apply
976cae8b89aSjeremylt   @param[in] in    CeedVector containing input state or NULL if there are no
977cae8b89aSjeremylt                      active inputs
978cae8b89aSjeremylt   @param[out] out  CeedVector to sum in result of applying operator (must be
979cae8b89aSjeremylt                      distinct from @a in) or NULL if there are no active outputs
980cae8b89aSjeremylt   @param request   Address of CeedRequest for non-blocking completion, else
9814cc79fe7SJed Brown                      @ref CEED_REQUEST_IMMEDIATE
982cae8b89aSjeremylt 
983cae8b89aSjeremylt   @return An error code: 0 - success, otherwise - failure
984cae8b89aSjeremylt 
9857a982d89SJeremy L. Thompson   @ref User
986cae8b89aSjeremylt **/
987cae8b89aSjeremylt int CeedOperatorApplyAdd(CeedOperator op, CeedVector in, CeedVector out,
988cae8b89aSjeremylt                          CeedRequest *request) {
989cae8b89aSjeremylt   int ierr;
990cae8b89aSjeremylt   Ceed ceed = op->ceed;
991250756a7Sjeremylt   ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr);
992cae8b89aSjeremylt 
993250756a7Sjeremylt   if (op->numelements)  {
994250756a7Sjeremylt     // Standard Operator
995250756a7Sjeremylt     ierr = op->ApplyAdd(op, in, out, request); CeedChk(ierr);
996250756a7Sjeremylt   } else if (op->composite) {
997250756a7Sjeremylt     // Composite Operator
998250756a7Sjeremylt     if (op->ApplyAddComposite) {
999250756a7Sjeremylt       ierr = op->ApplyAddComposite(op, in, out, request); CeedChk(ierr);
1000cae8b89aSjeremylt     } else {
1001250756a7Sjeremylt       CeedInt numsub;
1002250756a7Sjeremylt       ierr = CeedOperatorGetNumSub(op, &numsub); CeedChk(ierr);
1003250756a7Sjeremylt       CeedOperator *suboperators;
1004250756a7Sjeremylt       ierr = CeedOperatorGetSubList(op, &suboperators); CeedChk(ierr);
1005250756a7Sjeremylt 
1006250756a7Sjeremylt       for (CeedInt i=0; i<numsub; i++) {
1007250756a7Sjeremylt         ierr = CeedOperatorApplyAdd(suboperators[i], in, out, request);
1008cae8b89aSjeremylt         CeedChk(ierr);
10091d7d2407SJeremy L Thompson       }
1010250756a7Sjeremylt     }
1011250756a7Sjeremylt   }
1012250756a7Sjeremylt 
1013d7b241e6Sjeremylt   return 0;
1014d7b241e6Sjeremylt }
1015d7b241e6Sjeremylt 
1016d7b241e6Sjeremylt /**
1017b11c1e72Sjeremylt   @brief Destroy a CeedOperator
1018d7b241e6Sjeremylt 
1019d7b241e6Sjeremylt   @param op CeedOperator to destroy
1020b11c1e72Sjeremylt 
1021b11c1e72Sjeremylt   @return An error code: 0 - success, otherwise - failure
1022dfdf5a53Sjeremylt 
10237a982d89SJeremy L. Thompson   @ref User
1024b11c1e72Sjeremylt **/
1025d7b241e6Sjeremylt int CeedOperatorDestroy(CeedOperator *op) {
1026d7b241e6Sjeremylt   int ierr;
1027d7b241e6Sjeremylt 
1028d7b241e6Sjeremylt   if (!*op || --(*op)->refcount > 0) return 0;
1029d7b241e6Sjeremylt   if ((*op)->Destroy) {
1030d7b241e6Sjeremylt     ierr = (*op)->Destroy(*op); CeedChk(ierr);
1031d7b241e6Sjeremylt   }
1032fe2413ffSjeremylt   ierr = CeedDestroy(&(*op)->ceed); CeedChk(ierr);
1033fe2413ffSjeremylt   // Free fields
10341d102b48SJeremy L Thompson   for (int i=0; i<(*op)->nfields; i++)
103552d6035fSJeremy L Thompson     if ((*op)->inputfields[i]) {
103615910d16Sjeremylt       if ((*op)->inputfields[i]->Erestrict != CEED_ELEMRESTRICTION_NONE) {
103771352a93Sjeremylt         ierr = CeedElemRestrictionDestroy(&(*op)->inputfields[i]->Erestrict);
103871352a93Sjeremylt         CeedChk(ierr);
103915910d16Sjeremylt       }
104071352a93Sjeremylt       if ((*op)->inputfields[i]->basis != CEED_BASIS_COLLOCATED) {
104171352a93Sjeremylt         ierr = CeedBasisDestroy(&(*op)->inputfields[i]->basis); CeedChk(ierr);
104271352a93Sjeremylt       }
104371352a93Sjeremylt       if ((*op)->inputfields[i]->vec != CEED_VECTOR_ACTIVE &&
104471352a93Sjeremylt           (*op)->inputfields[i]->vec != CEED_VECTOR_NONE ) {
104571352a93Sjeremylt         ierr = CeedVectorDestroy(&(*op)->inputfields[i]->vec); CeedChk(ierr);
104671352a93Sjeremylt       }
1047fe2413ffSjeremylt       ierr = CeedFree(&(*op)->inputfields[i]); CeedChk(ierr);
1048fe2413ffSjeremylt     }
10491d102b48SJeremy L Thompson   for (int i=0; i<(*op)->nfields; i++)
1050d0eb30a5Sjeremylt     if ((*op)->outputfields[i]) {
105171352a93Sjeremylt       ierr = CeedElemRestrictionDestroy(&(*op)->outputfields[i]->Erestrict);
105271352a93Sjeremylt       CeedChk(ierr);
105371352a93Sjeremylt       if ((*op)->outputfields[i]->basis != CEED_BASIS_COLLOCATED) {
105471352a93Sjeremylt         ierr = CeedBasisDestroy(&(*op)->outputfields[i]->basis); CeedChk(ierr);
105571352a93Sjeremylt       }
105671352a93Sjeremylt       if ((*op)->outputfields[i]->vec != CEED_VECTOR_ACTIVE &&
105771352a93Sjeremylt           (*op)->outputfields[i]->vec != CEED_VECTOR_NONE ) {
105871352a93Sjeremylt         ierr = CeedVectorDestroy(&(*op)->outputfields[i]->vec); CeedChk(ierr);
105971352a93Sjeremylt       }
1060fe2413ffSjeremylt       ierr = CeedFree(&(*op)->outputfields[i]); CeedChk(ierr);
1061fe2413ffSjeremylt     }
106252d6035fSJeremy L Thompson   // Destroy suboperators
10631d102b48SJeremy L Thompson   for (int i=0; i<(*op)->numsub; i++)
106452d6035fSJeremy L Thompson     if ((*op)->suboperators[i]) {
106552d6035fSJeremy L Thompson       ierr = CeedOperatorDestroy(&(*op)->suboperators[i]); CeedChk(ierr);
106652d6035fSJeremy L Thompson     }
1067d7b241e6Sjeremylt   ierr = CeedQFunctionDestroy(&(*op)->qf); CeedChk(ierr);
1068d7b241e6Sjeremylt   ierr = CeedQFunctionDestroy(&(*op)->dqf); CeedChk(ierr);
1069d7b241e6Sjeremylt   ierr = CeedQFunctionDestroy(&(*op)->dqfT); CeedChk(ierr);
1070fe2413ffSjeremylt 
10715107b09fSJeremy L Thompson   // Destroy fallback
10725107b09fSJeremy L Thompson   if ((*op)->opfallback) {
10735107b09fSJeremy L Thompson     ierr = (*op)->qffallback->Destroy((*op)->qffallback); CeedChk(ierr);
10745107b09fSJeremy L Thompson     ierr = CeedFree(&(*op)->qffallback); CeedChk(ierr);
10755107b09fSJeremy L Thompson     ierr = (*op)->opfallback->Destroy((*op)->opfallback); CeedChk(ierr);
10765107b09fSJeremy L Thompson     ierr = CeedFree(&(*op)->opfallback); CeedChk(ierr);
10775107b09fSJeremy L Thompson   }
10785107b09fSJeremy L Thompson 
1079fe2413ffSjeremylt   ierr = CeedFree(&(*op)->inputfields); CeedChk(ierr);
1080fe2413ffSjeremylt   ierr = CeedFree(&(*op)->outputfields); CeedChk(ierr);
108152d6035fSJeremy L Thompson   ierr = CeedFree(&(*op)->suboperators); CeedChk(ierr);
1082d7b241e6Sjeremylt   ierr = CeedFree(op); CeedChk(ierr);
1083d7b241e6Sjeremylt   return 0;
1084d7b241e6Sjeremylt }
1085d7b241e6Sjeremylt 
1086d7b241e6Sjeremylt /// @}
1087