xref: /libCEED/interface/ceed-operator.c (revision d965c7a71afb833cd21b66efc6c6ec81679c68a9)
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 
529b3b7035fSJeremy L Thompson   if (!qf || qf == CEED_QFUNCTION_NONE)
530b3b7035fSJeremy L Thompson     // LCOV_EXCL_START
531b3b7035fSJeremy L Thompson     return CeedError(ceed, 1, "Operator must have a valid QFunction.");
532b3b7035fSJeremy 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 /**
782*d965c7a7SJeremy L Thompson   @brief Assemble the diagonal of a square linear CeedOperator
783b7ec98d8SJeremy L Thompson 
784b7ec98d8SJeremy L Thompson   This returns a CeedVector containing the diagonal of a linear CeedOperator.
785b7ec98d8SJeremy L Thompson 
786*d965c7a7SJeremy L Thompson   Note: Currently only non-composite CeedOperators with a single field are
787*d965c7a7SJeremy L Thompson           supported.
788*d965c7a7SJeremy L Thompson 
789b7ec98d8SJeremy L Thompson   @param op             CeedOperator to assemble CeedQFunction
790b7ec98d8SJeremy L Thompson   @param[out] assembled CeedVector to store assembled CeedOperator diagonal
791b7ec98d8SJeremy L Thompson   @param request        Address of CeedRequest for non-blocking completion, else
7924cc79fe7SJed Brown                           @ref CEED_REQUEST_IMMEDIATE
793b7ec98d8SJeremy L Thompson 
794b7ec98d8SJeremy L Thompson   @return An error code: 0 - success, otherwise - failure
795b7ec98d8SJeremy L Thompson 
7967a982d89SJeremy L. Thompson   @ref User
797b7ec98d8SJeremy L Thompson **/
7987f823360Sjeremylt int CeedOperatorAssembleLinearDiagonal(CeedOperator op, CeedVector *assembled,
7997f823360Sjeremylt                                        CeedRequest *request) {
800b7ec98d8SJeremy L Thompson   int ierr;
801b7ec98d8SJeremy L Thompson   Ceed ceed = op->ceed;
802250756a7Sjeremylt   ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr);
803b7ec98d8SJeremy L Thompson 
804b7ec98d8SJeremy L Thompson   // Use backend version, if available
805b7ec98d8SJeremy L Thompson   if (op->AssembleLinearDiagonal) {
806b7ec98d8SJeremy L Thompson     ierr = op->AssembleLinearDiagonal(op, assembled, request); CeedChk(ierr);
8077172caadSjeremylt   } else {
8087172caadSjeremylt     // Fallback to reference Ceed
8097172caadSjeremylt     if (!op->opfallback) {
8107172caadSjeremylt       ierr = CeedOperatorCreateFallback(op); CeedChk(ierr);
811b7ec98d8SJeremy L Thompson     }
8127172caadSjeremylt     // Assemble
8137172caadSjeremylt     ierr = op->opfallback->AssembleLinearDiagonal(op->opfallback, assembled,
8147172caadSjeremylt            request); CeedChk(ierr);
815b7ec98d8SJeremy L Thompson   }
816b7ec98d8SJeremy L Thompson 
817b7ec98d8SJeremy L Thompson   return 0;
818b7ec98d8SJeremy L Thompson }
819b7ec98d8SJeremy L Thompson 
820b7ec98d8SJeremy L Thompson /**
821*d965c7a7SJeremy L Thompson   @brief Assemble the point block diagonal of a square linear CeedOperator
822*d965c7a7SJeremy L Thompson 
823*d965c7a7SJeremy L Thompson   This returns a CeedVector containing the point block diagonal of a linear
824*d965c7a7SJeremy L Thompson     CeedOperator.
825*d965c7a7SJeremy L Thompson 
826*d965c7a7SJeremy L Thompson   Note: Currently only non-composite CeedOperators with a single field are
827*d965c7a7SJeremy L Thompson           supported.
828*d965c7a7SJeremy L Thompson 
829*d965c7a7SJeremy L Thompson   @param op             CeedOperator to assemble CeedQFunction
830*d965c7a7SJeremy L Thompson   @param[out] assembled CeedVector to store assembled CeedOperator point block
831*d965c7a7SJeremy L Thompson                           diagonal, provided in row-major form with an
832*d965c7a7SJeremy L Thompson                           @a ncomp * @a ncomp block at each node. The dimensions
833*d965c7a7SJeremy L Thompson                           of this vector are derived from the active vector
834*d965c7a7SJeremy L Thompson                           for the CeedOperator. The array has shape
835*d965c7a7SJeremy L Thompson                           [nodes, component out, component in].
836*d965c7a7SJeremy L Thompson   @param request        Address of CeedRequest for non-blocking completion, else
837*d965c7a7SJeremy L Thompson                           CEED_REQUEST_IMMEDIATE
838*d965c7a7SJeremy L Thompson 
839*d965c7a7SJeremy L Thompson   @return An error code: 0 - success, otherwise - failure
840*d965c7a7SJeremy L Thompson 
841*d965c7a7SJeremy L Thompson   @ref User
842*d965c7a7SJeremy L Thompson **/
843*d965c7a7SJeremy L Thompson int CeedOperatorAssembleLinearPointBlockDiagonal(CeedOperator op,
844*d965c7a7SJeremy L Thompson     CeedVector *assembled,
845*d965c7a7SJeremy L Thompson     CeedRequest *request) {
846*d965c7a7SJeremy L Thompson   int ierr;
847*d965c7a7SJeremy L Thompson   Ceed ceed = op->ceed;
848*d965c7a7SJeremy L Thompson   ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr);
849*d965c7a7SJeremy L Thompson 
850*d965c7a7SJeremy L Thompson   // Use backend version, if available
851*d965c7a7SJeremy L Thompson   if (op->AssembleLinearPointBlockDiagonal) {
852*d965c7a7SJeremy L Thompson     ierr = op->AssembleLinearPointBlockDiagonal(op, assembled, request);
853*d965c7a7SJeremy L Thompson     CeedChk(ierr);
854*d965c7a7SJeremy L Thompson   } else {
855*d965c7a7SJeremy L Thompson     // Fallback to reference Ceed
856*d965c7a7SJeremy L Thompson     if (!op->opfallback) {
857*d965c7a7SJeremy L Thompson       ierr = CeedOperatorCreateFallback(op); CeedChk(ierr);
858*d965c7a7SJeremy L Thompson     }
859*d965c7a7SJeremy L Thompson     // Assemble
860*d965c7a7SJeremy L Thompson     ierr = op->opfallback->AssembleLinearPointBlockDiagonal(op->opfallback,
861*d965c7a7SJeremy L Thompson            assembled, request); CeedChk(ierr);
862*d965c7a7SJeremy L Thompson   }
863*d965c7a7SJeremy L Thompson 
864*d965c7a7SJeremy L Thompson   return 0;
865*d965c7a7SJeremy L Thompson }
866*d965c7a7SJeremy L Thompson 
867*d965c7a7SJeremy L Thompson /**
8683bd813ffSjeremylt   @brief Build a FDM based approximate inverse for each element for a
8693bd813ffSjeremylt            CeedOperator
8703bd813ffSjeremylt 
8713bd813ffSjeremylt   This returns a CeedOperator and CeedVector to apply a Fast Diagonalization
8723bd813ffSjeremylt     Method based approximate inverse. This function obtains the simultaneous
8733bd813ffSjeremylt     diagonalization for the 1D mass and Laplacian operators,
8743bd813ffSjeremylt       M = V^T V, K = V^T S V.
8753bd813ffSjeremylt     The assembled QFunction is used to modify the eigenvalues from simultaneous
8763bd813ffSjeremylt     diagonalization and obtain an approximate inverse of the form
8773bd813ffSjeremylt       V^T S^hat V. The CeedOperator must be linear and non-composite. The
8783bd813ffSjeremylt     associated CeedQFunction must therefore also be linear.
8793bd813ffSjeremylt 
8803bd813ffSjeremylt   @param op             CeedOperator to create element inverses
8813bd813ffSjeremylt   @param[out] fdminv    CeedOperator to apply the action of a FDM based inverse
8823bd813ffSjeremylt                           for each element
8833bd813ffSjeremylt   @param request        Address of CeedRequest for non-blocking completion, else
8844cc79fe7SJed Brown                           @ref CEED_REQUEST_IMMEDIATE
8853bd813ffSjeremylt 
8863bd813ffSjeremylt   @return An error code: 0 - success, otherwise - failure
8873bd813ffSjeremylt 
8887a982d89SJeremy L. Thompson   @ref Backend
8893bd813ffSjeremylt **/
8903bd813ffSjeremylt int CeedOperatorCreateFDMElementInverse(CeedOperator op, CeedOperator *fdminv,
8913bd813ffSjeremylt                                         CeedRequest *request) {
8923bd813ffSjeremylt   int ierr;
8933bd813ffSjeremylt   Ceed ceed = op->ceed;
894713f43c3Sjeremylt   ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr);
8953bd813ffSjeremylt 
896713f43c3Sjeremylt   // Use backend version, if available
897713f43c3Sjeremylt   if (op->CreateFDMElementInverse) {
898713f43c3Sjeremylt     ierr = op->CreateFDMElementInverse(op, fdminv, request); CeedChk(ierr);
899713f43c3Sjeremylt   } else {
900713f43c3Sjeremylt     // Fallback to reference Ceed
901713f43c3Sjeremylt     if (!op->opfallback) {
902713f43c3Sjeremylt       ierr = CeedOperatorCreateFallback(op); CeedChk(ierr);
9033bd813ffSjeremylt     }
904713f43c3Sjeremylt     // Assemble
905713f43c3Sjeremylt     ierr = op->opfallback->CreateFDMElementInverse(op->opfallback, fdminv,
9063bd813ffSjeremylt            request); CeedChk(ierr);
9073bd813ffSjeremylt   }
9083bd813ffSjeremylt 
9093bd813ffSjeremylt   return 0;
9103bd813ffSjeremylt }
9113bd813ffSjeremylt 
9127a982d89SJeremy L. Thompson /**
9137a982d89SJeremy L. Thompson   @brief View a CeedOperator
9147a982d89SJeremy L. Thompson 
9157a982d89SJeremy L. Thompson   @param[in] op     CeedOperator to view
9167a982d89SJeremy L. Thompson   @param[in] stream Stream to write; typically stdout/stderr or a file
9177a982d89SJeremy L. Thompson 
9187a982d89SJeremy L. Thompson   @return Error code: 0 - success, otherwise - failure
9197a982d89SJeremy L. Thompson 
9207a982d89SJeremy L. Thompson   @ref User
9217a982d89SJeremy L. Thompson **/
9227a982d89SJeremy L. Thompson int CeedOperatorView(CeedOperator op, FILE *stream) {
9237a982d89SJeremy L. Thompson   int ierr;
9247a982d89SJeremy L. Thompson 
9257a982d89SJeremy L. Thompson   if (op->composite) {
9267a982d89SJeremy L. Thompson     fprintf(stream, "Composite CeedOperator\n");
9277a982d89SJeremy L. Thompson 
9287a982d89SJeremy L. Thompson     for (CeedInt i=0; i<op->numsub; i++) {
9297a982d89SJeremy L. Thompson       fprintf(stream, "  SubOperator [%d]:\n", i);
9307a982d89SJeremy L. Thompson       ierr = CeedOperatorSingleView(op->suboperators[i], 1, stream);
9317a982d89SJeremy L. Thompson       CeedChk(ierr);
9327a982d89SJeremy L. Thompson     }
9337a982d89SJeremy L. Thompson   } else {
9347a982d89SJeremy L. Thompson     fprintf(stream, "CeedOperator\n");
9357a982d89SJeremy L. Thompson     ierr = CeedOperatorSingleView(op, 0, stream); CeedChk(ierr);
9367a982d89SJeremy L. Thompson   }
9377a982d89SJeremy L. Thompson 
9387a982d89SJeremy L. Thompson   return 0;
9397a982d89SJeremy L. Thompson }
9403bd813ffSjeremylt 
9413bd813ffSjeremylt /**
9423bd813ffSjeremylt   @brief Apply CeedOperator to a vector
943d7b241e6Sjeremylt 
944d7b241e6Sjeremylt   This computes the action of the operator on the specified (active) input,
945d7b241e6Sjeremylt   yielding its (active) output.  All inputs and outputs must be specified using
946d7b241e6Sjeremylt   CeedOperatorSetField().
947d7b241e6Sjeremylt 
948d7b241e6Sjeremylt   @param op        CeedOperator to apply
9494cc79fe7SJed Brown   @param[in] in    CeedVector containing input state or @ref CEED_VECTOR_NONE if
95034138859Sjeremylt                   there are no active inputs
951b11c1e72Sjeremylt   @param[out] out  CeedVector to store result of applying operator (must be
9524cc79fe7SJed Brown                      distinct from @a in) or @ref CEED_VECTOR_NONE if there are no
95334138859Sjeremylt                      active outputs
954d7b241e6Sjeremylt   @param request   Address of CeedRequest for non-blocking completion, else
9554cc79fe7SJed Brown                      @ref CEED_REQUEST_IMMEDIATE
956b11c1e72Sjeremylt 
957b11c1e72Sjeremylt   @return An error code: 0 - success, otherwise - failure
958dfdf5a53Sjeremylt 
9597a982d89SJeremy L. Thompson   @ref User
960b11c1e72Sjeremylt **/
961692c2638Sjeremylt int CeedOperatorApply(CeedOperator op, CeedVector in, CeedVector out,
962692c2638Sjeremylt                       CeedRequest *request) {
963d7b241e6Sjeremylt   int ierr;
964d7b241e6Sjeremylt   Ceed ceed = op->ceed;
965250756a7Sjeremylt   ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr);
966d7b241e6Sjeremylt 
967250756a7Sjeremylt   if (op->numelements)  {
968250756a7Sjeremylt     // Standard Operator
969cae8b89aSjeremylt     if (op->Apply) {
970250756a7Sjeremylt       ierr = op->Apply(op, in, out, request); CeedChk(ierr);
971cae8b89aSjeremylt     } else {
972cae8b89aSjeremylt       // Zero all output vectors
973250756a7Sjeremylt       CeedQFunction qf = op->qf;
974cae8b89aSjeremylt       for (CeedInt i=0; i<qf->numoutputfields; i++) {
975cae8b89aSjeremylt         CeedVector vec = op->outputfields[i]->vec;
976cae8b89aSjeremylt         if (vec == CEED_VECTOR_ACTIVE)
977cae8b89aSjeremylt           vec = out;
978cae8b89aSjeremylt         if (vec != CEED_VECTOR_NONE) {
979cae8b89aSjeremylt           ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr);
980cae8b89aSjeremylt         }
981cae8b89aSjeremylt       }
982250756a7Sjeremylt       // Apply
983250756a7Sjeremylt       ierr = op->ApplyAdd(op, in, out, request); CeedChk(ierr);
984250756a7Sjeremylt     }
985250756a7Sjeremylt   } else if (op->composite) {
986250756a7Sjeremylt     // Composite Operator
987250756a7Sjeremylt     if (op->ApplyComposite) {
988250756a7Sjeremylt       ierr = op->ApplyComposite(op, in, out, request); CeedChk(ierr);
989250756a7Sjeremylt     } else {
990250756a7Sjeremylt       CeedInt numsub;
991250756a7Sjeremylt       ierr = CeedOperatorGetNumSub(op, &numsub); CeedChk(ierr);
992250756a7Sjeremylt       CeedOperator *suboperators;
993250756a7Sjeremylt       ierr = CeedOperatorGetSubList(op, &suboperators); CeedChk(ierr);
994250756a7Sjeremylt 
995250756a7Sjeremylt       // Zero all output vectors
996250756a7Sjeremylt       if (out != CEED_VECTOR_NONE) {
997cae8b89aSjeremylt         ierr = CeedVectorSetValue(out, 0.0); CeedChk(ierr);
998cae8b89aSjeremylt       }
999250756a7Sjeremylt       for (CeedInt i=0; i<numsub; i++) {
1000250756a7Sjeremylt         for (CeedInt j=0; j<suboperators[i]->qf->numoutputfields; j++) {
1001250756a7Sjeremylt           CeedVector vec = suboperators[i]->outputfields[j]->vec;
1002250756a7Sjeremylt           if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) {
1003250756a7Sjeremylt             ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr);
1004250756a7Sjeremylt           }
1005250756a7Sjeremylt         }
1006250756a7Sjeremylt       }
1007250756a7Sjeremylt       // Apply
1008250756a7Sjeremylt       for (CeedInt i=0; i<op->numsub; i++) {
1009250756a7Sjeremylt         ierr = CeedOperatorApplyAdd(op->suboperators[i], in, out, request);
1010cae8b89aSjeremylt         CeedChk(ierr);
1011cae8b89aSjeremylt       }
1012cae8b89aSjeremylt     }
1013250756a7Sjeremylt   }
1014250756a7Sjeremylt 
1015cae8b89aSjeremylt   return 0;
1016cae8b89aSjeremylt }
1017cae8b89aSjeremylt 
1018cae8b89aSjeremylt /**
1019cae8b89aSjeremylt   @brief Apply CeedOperator to a vector and add result to output vector
1020cae8b89aSjeremylt 
1021cae8b89aSjeremylt   This computes the action of the operator on the specified (active) input,
1022cae8b89aSjeremylt   yielding its (active) output.  All inputs and outputs must be specified using
1023cae8b89aSjeremylt   CeedOperatorSetField().
1024cae8b89aSjeremylt 
1025cae8b89aSjeremylt   @param op        CeedOperator to apply
1026cae8b89aSjeremylt   @param[in] in    CeedVector containing input state or NULL if there are no
1027cae8b89aSjeremylt                      active inputs
1028cae8b89aSjeremylt   @param[out] out  CeedVector to sum in result of applying operator (must be
1029cae8b89aSjeremylt                      distinct from @a in) or NULL if there are no active outputs
1030cae8b89aSjeremylt   @param request   Address of CeedRequest for non-blocking completion, else
10314cc79fe7SJed Brown                      @ref CEED_REQUEST_IMMEDIATE
1032cae8b89aSjeremylt 
1033cae8b89aSjeremylt   @return An error code: 0 - success, otherwise - failure
1034cae8b89aSjeremylt 
10357a982d89SJeremy L. Thompson   @ref User
1036cae8b89aSjeremylt **/
1037cae8b89aSjeremylt int CeedOperatorApplyAdd(CeedOperator op, CeedVector in, CeedVector out,
1038cae8b89aSjeremylt                          CeedRequest *request) {
1039cae8b89aSjeremylt   int ierr;
1040cae8b89aSjeremylt   Ceed ceed = op->ceed;
1041250756a7Sjeremylt   ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr);
1042cae8b89aSjeremylt 
1043250756a7Sjeremylt   if (op->numelements)  {
1044250756a7Sjeremylt     // Standard Operator
1045250756a7Sjeremylt     ierr = op->ApplyAdd(op, in, out, request); CeedChk(ierr);
1046250756a7Sjeremylt   } else if (op->composite) {
1047250756a7Sjeremylt     // Composite Operator
1048250756a7Sjeremylt     if (op->ApplyAddComposite) {
1049250756a7Sjeremylt       ierr = op->ApplyAddComposite(op, in, out, request); CeedChk(ierr);
1050cae8b89aSjeremylt     } else {
1051250756a7Sjeremylt       CeedInt numsub;
1052250756a7Sjeremylt       ierr = CeedOperatorGetNumSub(op, &numsub); CeedChk(ierr);
1053250756a7Sjeremylt       CeedOperator *suboperators;
1054250756a7Sjeremylt       ierr = CeedOperatorGetSubList(op, &suboperators); CeedChk(ierr);
1055250756a7Sjeremylt 
1056250756a7Sjeremylt       for (CeedInt i=0; i<numsub; i++) {
1057250756a7Sjeremylt         ierr = CeedOperatorApplyAdd(suboperators[i], in, out, request);
1058cae8b89aSjeremylt         CeedChk(ierr);
10591d7d2407SJeremy L Thompson       }
1060250756a7Sjeremylt     }
1061250756a7Sjeremylt   }
1062250756a7Sjeremylt 
1063d7b241e6Sjeremylt   return 0;
1064d7b241e6Sjeremylt }
1065d7b241e6Sjeremylt 
1066d7b241e6Sjeremylt /**
1067b11c1e72Sjeremylt   @brief Destroy a CeedOperator
1068d7b241e6Sjeremylt 
1069d7b241e6Sjeremylt   @param op CeedOperator to destroy
1070b11c1e72Sjeremylt 
1071b11c1e72Sjeremylt   @return An error code: 0 - success, otherwise - failure
1072dfdf5a53Sjeremylt 
10737a982d89SJeremy L. Thompson   @ref User
1074b11c1e72Sjeremylt **/
1075d7b241e6Sjeremylt int CeedOperatorDestroy(CeedOperator *op) {
1076d7b241e6Sjeremylt   int ierr;
1077d7b241e6Sjeremylt 
1078d7b241e6Sjeremylt   if (!*op || --(*op)->refcount > 0) return 0;
1079d7b241e6Sjeremylt   if ((*op)->Destroy) {
1080d7b241e6Sjeremylt     ierr = (*op)->Destroy(*op); CeedChk(ierr);
1081d7b241e6Sjeremylt   }
1082fe2413ffSjeremylt   ierr = CeedDestroy(&(*op)->ceed); CeedChk(ierr);
1083fe2413ffSjeremylt   // Free fields
10841d102b48SJeremy L Thompson   for (int i=0; i<(*op)->nfields; i++)
108552d6035fSJeremy L Thompson     if ((*op)->inputfields[i]) {
108615910d16Sjeremylt       if ((*op)->inputfields[i]->Erestrict != CEED_ELEMRESTRICTION_NONE) {
108771352a93Sjeremylt         ierr = CeedElemRestrictionDestroy(&(*op)->inputfields[i]->Erestrict);
108871352a93Sjeremylt         CeedChk(ierr);
108915910d16Sjeremylt       }
109071352a93Sjeremylt       if ((*op)->inputfields[i]->basis != CEED_BASIS_COLLOCATED) {
109171352a93Sjeremylt         ierr = CeedBasisDestroy(&(*op)->inputfields[i]->basis); CeedChk(ierr);
109271352a93Sjeremylt       }
109371352a93Sjeremylt       if ((*op)->inputfields[i]->vec != CEED_VECTOR_ACTIVE &&
109471352a93Sjeremylt           (*op)->inputfields[i]->vec != CEED_VECTOR_NONE ) {
109571352a93Sjeremylt         ierr = CeedVectorDestroy(&(*op)->inputfields[i]->vec); CeedChk(ierr);
109671352a93Sjeremylt       }
1097fe2413ffSjeremylt       ierr = CeedFree(&(*op)->inputfields[i]); CeedChk(ierr);
1098fe2413ffSjeremylt     }
10991d102b48SJeremy L Thompson   for (int i=0; i<(*op)->nfields; i++)
1100d0eb30a5Sjeremylt     if ((*op)->outputfields[i]) {
110171352a93Sjeremylt       ierr = CeedElemRestrictionDestroy(&(*op)->outputfields[i]->Erestrict);
110271352a93Sjeremylt       CeedChk(ierr);
110371352a93Sjeremylt       if ((*op)->outputfields[i]->basis != CEED_BASIS_COLLOCATED) {
110471352a93Sjeremylt         ierr = CeedBasisDestroy(&(*op)->outputfields[i]->basis); CeedChk(ierr);
110571352a93Sjeremylt       }
110671352a93Sjeremylt       if ((*op)->outputfields[i]->vec != CEED_VECTOR_ACTIVE &&
110771352a93Sjeremylt           (*op)->outputfields[i]->vec != CEED_VECTOR_NONE ) {
110871352a93Sjeremylt         ierr = CeedVectorDestroy(&(*op)->outputfields[i]->vec); CeedChk(ierr);
110971352a93Sjeremylt       }
1110fe2413ffSjeremylt       ierr = CeedFree(&(*op)->outputfields[i]); CeedChk(ierr);
1111fe2413ffSjeremylt     }
111252d6035fSJeremy L Thompson   // Destroy suboperators
11131d102b48SJeremy L Thompson   for (int i=0; i<(*op)->numsub; i++)
111452d6035fSJeremy L Thompson     if ((*op)->suboperators[i]) {
111552d6035fSJeremy L Thompson       ierr = CeedOperatorDestroy(&(*op)->suboperators[i]); CeedChk(ierr);
111652d6035fSJeremy L Thompson     }
1117d7b241e6Sjeremylt   ierr = CeedQFunctionDestroy(&(*op)->qf); CeedChk(ierr);
1118d7b241e6Sjeremylt   ierr = CeedQFunctionDestroy(&(*op)->dqf); CeedChk(ierr);
1119d7b241e6Sjeremylt   ierr = CeedQFunctionDestroy(&(*op)->dqfT); CeedChk(ierr);
1120fe2413ffSjeremylt 
11215107b09fSJeremy L Thompson   // Destroy fallback
11225107b09fSJeremy L Thompson   if ((*op)->opfallback) {
11235107b09fSJeremy L Thompson     ierr = (*op)->qffallback->Destroy((*op)->qffallback); CeedChk(ierr);
11245107b09fSJeremy L Thompson     ierr = CeedFree(&(*op)->qffallback); CeedChk(ierr);
11255107b09fSJeremy L Thompson     ierr = (*op)->opfallback->Destroy((*op)->opfallback); CeedChk(ierr);
11265107b09fSJeremy L Thompson     ierr = CeedFree(&(*op)->opfallback); CeedChk(ierr);
11275107b09fSJeremy L Thompson   }
11285107b09fSJeremy L Thompson 
1129fe2413ffSjeremylt   ierr = CeedFree(&(*op)->inputfields); CeedChk(ierr);
1130fe2413ffSjeremylt   ierr = CeedFree(&(*op)->outputfields); CeedChk(ierr);
113152d6035fSJeremy L Thompson   ierr = CeedFree(&(*op)->suboperators); CeedChk(ierr);
1132d7b241e6Sjeremylt   ierr = CeedFree(op); CeedChk(ierr);
1133d7b241e6Sjeremylt   return 0;
1134d7b241e6Sjeremylt }
1135d7b241e6Sjeremylt 
1136d7b241e6Sjeremylt /// @}
1137