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 291fd364f38SJeremy L Thompson @param[out] issetupdone 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 298fd364f38SJeremy L Thompson int CeedOperatorIsSetupDone(CeedOperator op, bool *issetupdone) { 299fd364f38SJeremy L Thompson *issetupdone = 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 /** 325c04a41a7SJeremy L Thompson @brief Get a boolean value indicating if the CeedOperator is composite 326c04a41a7SJeremy L Thompson 327c04a41a7SJeremy L Thompson @param op CeedOperator 328fd364f38SJeremy L Thompson @param[out] iscomposite Variable to store composite status 329c04a41a7SJeremy L Thompson 330c04a41a7SJeremy L Thompson @return An error code: 0 - success, otherwise - failure 331c04a41a7SJeremy L Thompson 332c04a41a7SJeremy L Thompson @ref Backend 333c04a41a7SJeremy L Thompson **/ 334c04a41a7SJeremy L Thompson 335fd364f38SJeremy L Thompson int CeedOperatorIsComposite(CeedOperator op, bool *iscomposite) { 336fd364f38SJeremy L Thompson *iscomposite = op->composite; 337c04a41a7SJeremy L Thompson return 0; 338c04a41a7SJeremy L Thompson } 339c04a41a7SJeremy L Thompson 340c04a41a7SJeremy L Thompson /** 3417a982d89SJeremy L. Thompson @brief Get the number of suboperators associated with a CeedOperator 3427a982d89SJeremy L. Thompson 3437a982d89SJeremy L. Thompson @param op CeedOperator 3447a982d89SJeremy L. Thompson @param[out] numsub Variable to store number of suboperators 3457a982d89SJeremy L. Thompson 3467a982d89SJeremy L. Thompson @return An error code: 0 - success, otherwise - failure 3477a982d89SJeremy L. Thompson 3487a982d89SJeremy L. Thompson @ref Backend 3497a982d89SJeremy L. Thompson **/ 3507a982d89SJeremy L. Thompson 3517a982d89SJeremy L. Thompson int CeedOperatorGetNumSub(CeedOperator op, CeedInt *numsub) { 3527a982d89SJeremy L. Thompson if (!op->composite) 3537a982d89SJeremy L. Thompson // LCOV_EXCL_START 3547a982d89SJeremy L. Thompson return CeedError(op->ceed, 1, "Not a composite operator"); 3557a982d89SJeremy L. Thompson // LCOV_EXCL_STOP 3567a982d89SJeremy L. Thompson 3577a982d89SJeremy L. Thompson *numsub = op->numsub; 3587a982d89SJeremy L. Thompson return 0; 3597a982d89SJeremy L. Thompson } 3607a982d89SJeremy L. Thompson 3617a982d89SJeremy L. Thompson /** 3627a982d89SJeremy L. Thompson @brief Get the list of suboperators associated with a CeedOperator 3637a982d89SJeremy L. Thompson 3647a982d89SJeremy L. Thompson @param op CeedOperator 3657a982d89SJeremy L. Thompson @param[out] suboperators Variable to store list of suboperators 3667a982d89SJeremy L. Thompson 3677a982d89SJeremy L. Thompson @return An error code: 0 - success, otherwise - failure 3687a982d89SJeremy L. Thompson 3697a982d89SJeremy L. Thompson @ref Backend 3707a982d89SJeremy L. Thompson **/ 3717a982d89SJeremy L. Thompson 3727a982d89SJeremy L. Thompson int CeedOperatorGetSubList(CeedOperator op, CeedOperator **suboperators) { 3737a982d89SJeremy L. Thompson if (!op->composite) 3747a982d89SJeremy L. Thompson // LCOV_EXCL_START 3757a982d89SJeremy L. Thompson return CeedError(op->ceed, 1, "Not a composite operator"); 3767a982d89SJeremy L. Thompson // LCOV_EXCL_STOP 3777a982d89SJeremy L. Thompson 3787a982d89SJeremy L. Thompson *suboperators = op->suboperators; 3797a982d89SJeremy L. Thompson return 0; 3807a982d89SJeremy L. Thompson } 3817a982d89SJeremy L. Thompson 3827a982d89SJeremy L. Thompson /** 3837a982d89SJeremy L. Thompson @brief Get the backend data of a CeedOperator 3847a982d89SJeremy L. Thompson 3857a982d89SJeremy L. Thompson @param op CeedOperator 3867a982d89SJeremy L. Thompson @param[out] data Variable to store data 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 CeedOperatorGetData(CeedOperator op, void **data) { 3947a982d89SJeremy L. Thompson *data = op->data; 3957a982d89SJeremy L. Thompson return 0; 3967a982d89SJeremy L. Thompson } 3977a982d89SJeremy L. Thompson 3987a982d89SJeremy L. Thompson /** 3997a982d89SJeremy L. Thompson @brief Set the backend data of a CeedOperator 4007a982d89SJeremy L. Thompson 4017a982d89SJeremy L. Thompson @param[out] op CeedOperator 4027a982d89SJeremy L. Thompson @param data Data to set 4037a982d89SJeremy L. Thompson 4047a982d89SJeremy L. Thompson @return An error code: 0 - success, otherwise - failure 4057a982d89SJeremy L. Thompson 4067a982d89SJeremy L. Thompson @ref Backend 4077a982d89SJeremy L. Thompson **/ 4087a982d89SJeremy L. Thompson 4097a982d89SJeremy L. Thompson int CeedOperatorSetData(CeedOperator op, void **data) { 4107a982d89SJeremy L. Thompson op->data = *data; 4117a982d89SJeremy L. Thompson return 0; 4127a982d89SJeremy L. Thompson } 4137a982d89SJeremy L. Thompson 4147a982d89SJeremy L. Thompson /** 4157a982d89SJeremy L. Thompson @brief Set the setup flag of a CeedOperator to True 4167a982d89SJeremy L. Thompson 4177a982d89SJeremy L. Thompson @param op CeedOperator 4187a982d89SJeremy L. Thompson 4197a982d89SJeremy L. Thompson @return An error code: 0 - success, otherwise - failure 4207a982d89SJeremy L. Thompson 4217a982d89SJeremy L. Thompson @ref Backend 4227a982d89SJeremy L. Thompson **/ 4237a982d89SJeremy L. Thompson 4247a982d89SJeremy L. Thompson int CeedOperatorSetSetupDone(CeedOperator op) { 4257a982d89SJeremy L. Thompson op->setupdone = 1; 4267a982d89SJeremy L. Thompson return 0; 4277a982d89SJeremy L. Thompson } 4287a982d89SJeremy L. Thompson 4297a982d89SJeremy L. Thompson /** 4307a982d89SJeremy L. Thompson @brief Get the CeedOperatorFields of a CeedOperator 4317a982d89SJeremy L. Thompson 4327a982d89SJeremy L. Thompson @param op CeedOperator 4337a982d89SJeremy L. Thompson @param[out] inputfields Variable to store inputfields 4347a982d89SJeremy L. Thompson @param[out] outputfields Variable to store outputfields 4357a982d89SJeremy L. Thompson 4367a982d89SJeremy L. Thompson @return An error code: 0 - success, otherwise - failure 4377a982d89SJeremy L. Thompson 4387a982d89SJeremy L. Thompson @ref Backend 4397a982d89SJeremy L. Thompson **/ 4407a982d89SJeremy L. Thompson 4417a982d89SJeremy L. Thompson int CeedOperatorGetFields(CeedOperator op, CeedOperatorField **inputfields, 4427a982d89SJeremy L. Thompson CeedOperatorField **outputfields) { 4437a982d89SJeremy L. Thompson if (op->composite) 4447a982d89SJeremy L. Thompson // LCOV_EXCL_START 4457a982d89SJeremy L. Thompson return CeedError(op->ceed, 1, "Not defined for composite operator"); 4467a982d89SJeremy L. Thompson // LCOV_EXCL_STOP 4477a982d89SJeremy L. Thompson 4487a982d89SJeremy L. Thompson if (inputfields) *inputfields = op->inputfields; 4497a982d89SJeremy L. Thompson if (outputfields) *outputfields = op->outputfields; 4507a982d89SJeremy L. Thompson return 0; 4517a982d89SJeremy L. Thompson } 4527a982d89SJeremy L. Thompson 4537a982d89SJeremy L. Thompson /** 4547a982d89SJeremy L. Thompson @brief Get the CeedElemRestriction of a CeedOperatorField 4557a982d89SJeremy L. Thompson 4567a982d89SJeremy L. Thompson @param opfield CeedOperatorField 4577a982d89SJeremy L. Thompson @param[out] rstr Variable to store CeedElemRestriction 4587a982d89SJeremy L. Thompson 4597a982d89SJeremy L. Thompson @return An error code: 0 - success, otherwise - failure 4607a982d89SJeremy L. Thompson 4617a982d89SJeremy L. Thompson @ref Backend 4627a982d89SJeremy L. Thompson **/ 4637a982d89SJeremy L. Thompson 4647a982d89SJeremy L. Thompson int CeedOperatorFieldGetElemRestriction(CeedOperatorField opfield, 4657a982d89SJeremy L. Thompson CeedElemRestriction *rstr) { 4667a982d89SJeremy L. Thompson *rstr = opfield->Erestrict; 4677a982d89SJeremy L. Thompson return 0; 4687a982d89SJeremy L. Thompson } 4697a982d89SJeremy L. Thompson 4707a982d89SJeremy L. Thompson /** 4717a982d89SJeremy L. Thompson @brief Get the CeedBasis of a CeedOperatorField 4727a982d89SJeremy L. Thompson 4737a982d89SJeremy L. Thompson @param opfield CeedOperatorField 4747a982d89SJeremy L. Thompson @param[out] basis Variable to store CeedBasis 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 CeedOperatorFieldGetBasis(CeedOperatorField opfield, CeedBasis *basis) { 4827a982d89SJeremy L. Thompson *basis = opfield->basis; 4837a982d89SJeremy L. Thompson return 0; 4847a982d89SJeremy L. Thompson } 4857a982d89SJeremy L. Thompson 4867a982d89SJeremy L. Thompson /** 4877a982d89SJeremy L. Thompson @brief Get the CeedVector of a CeedOperatorField 4887a982d89SJeremy L. Thompson 4897a982d89SJeremy L. Thompson @param opfield CeedOperatorField 4907a982d89SJeremy L. Thompson @param[out] vec Variable to store CeedVector 4917a982d89SJeremy L. Thompson 4927a982d89SJeremy L. Thompson @return An error code: 0 - success, otherwise - failure 4937a982d89SJeremy L. Thompson 4947a982d89SJeremy L. Thompson @ref Backend 4957a982d89SJeremy L. Thompson **/ 4967a982d89SJeremy L. Thompson 4977a982d89SJeremy L. Thompson int CeedOperatorFieldGetVector(CeedOperatorField opfield, CeedVector *vec) { 4987a982d89SJeremy L. Thompson *vec = opfield->vec; 4997a982d89SJeremy L. Thompson return 0; 5007a982d89SJeremy L. Thompson } 5017a982d89SJeremy L. Thompson 5027a982d89SJeremy L. Thompson /// @} 5037a982d89SJeremy L. Thompson 5047a982d89SJeremy L. Thompson /// ---------------------------------------------------------------------------- 5057a982d89SJeremy L. Thompson /// CeedOperator Public API 5067a982d89SJeremy L. Thompson /// ---------------------------------------------------------------------------- 5077a982d89SJeremy L. Thompson /// @addtogroup CeedOperatorUser 508dfdf5a53Sjeremylt /// @{ 509d7b241e6Sjeremylt 510d7b241e6Sjeremylt /** 5110219ea01SJeremy L Thompson @brief Create a CeedOperator and associate a CeedQFunction. A CeedBasis and 5120219ea01SJeremy L Thompson CeedElemRestriction can be associated with CeedQFunction fields with 5130219ea01SJeremy L Thompson \ref CeedOperatorSetField. 514d7b241e6Sjeremylt 515b11c1e72Sjeremylt @param ceed A Ceed object where the CeedOperator will be created 516d7b241e6Sjeremylt @param qf QFunction defining the action of the operator at quadrature points 51734138859Sjeremylt @param dqf QFunction defining the action of the Jacobian of @a qf (or 5184cc79fe7SJed Brown @ref CEED_QFUNCTION_NONE) 519d7b241e6Sjeremylt @param dqfT QFunction defining the action of the transpose of the Jacobian 5204cc79fe7SJed Brown of @a qf (or @ref CEED_QFUNCTION_NONE) 521b11c1e72Sjeremylt @param[out] op Address of the variable where the newly created 522b11c1e72Sjeremylt CeedOperator will be stored 523b11c1e72Sjeremylt 524b11c1e72Sjeremylt @return An error code: 0 - success, otherwise - failure 525dfdf5a53Sjeremylt 5267a982d89SJeremy L. Thompson @ref User 527d7b241e6Sjeremylt */ 528d7b241e6Sjeremylt int CeedOperatorCreate(Ceed ceed, CeedQFunction qf, CeedQFunction dqf, 529d7b241e6Sjeremylt CeedQFunction dqfT, CeedOperator *op) { 530d7b241e6Sjeremylt int ierr; 531d7b241e6Sjeremylt 5325fe0d4faSjeremylt if (!ceed->OperatorCreate) { 5335fe0d4faSjeremylt Ceed delegate; 534aefd8378Sjeremylt ierr = CeedGetObjectDelegate(ceed, &delegate, "Operator"); CeedChk(ierr); 5355fe0d4faSjeremylt 5365fe0d4faSjeremylt if (!delegate) 537c042f62fSJeremy L Thompson // LCOV_EXCL_START 5385fe0d4faSjeremylt return CeedError(ceed, 1, "Backend does not support OperatorCreate"); 539c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 5405fe0d4faSjeremylt 5415fe0d4faSjeremylt ierr = CeedOperatorCreate(delegate, qf, dqf, dqfT, op); CeedChk(ierr); 5425fe0d4faSjeremylt return 0; 5435fe0d4faSjeremylt } 5445fe0d4faSjeremylt 545b3b7035fSJeremy L Thompson if (!qf || qf == CEED_QFUNCTION_NONE) 546b3b7035fSJeremy L Thompson // LCOV_EXCL_START 547b3b7035fSJeremy L Thompson return CeedError(ceed, 1, "Operator must have a valid QFunction."); 548b3b7035fSJeremy L Thompson // LCOV_EXCL_STOP 549d7b241e6Sjeremylt ierr = CeedCalloc(1, op); CeedChk(ierr); 550d7b241e6Sjeremylt (*op)->ceed = ceed; 551d7b241e6Sjeremylt ceed->refcount++; 552d7b241e6Sjeremylt (*op)->refcount = 1; 553d7b241e6Sjeremylt (*op)->qf = qf; 554d7b241e6Sjeremylt qf->refcount++; 555442e7f0bSjeremylt if (dqf && dqf != CEED_QFUNCTION_NONE) { 556d7b241e6Sjeremylt (*op)->dqf = dqf; 557442e7f0bSjeremylt dqf->refcount++; 558442e7f0bSjeremylt } 559442e7f0bSjeremylt if (dqfT && dqfT != CEED_QFUNCTION_NONE) { 560d7b241e6Sjeremylt (*op)->dqfT = dqfT; 561442e7f0bSjeremylt dqfT->refcount++; 562442e7f0bSjeremylt } 563fe2413ffSjeremylt ierr = CeedCalloc(16, &(*op)->inputfields); CeedChk(ierr); 564fe2413ffSjeremylt ierr = CeedCalloc(16, &(*op)->outputfields); CeedChk(ierr); 565d7b241e6Sjeremylt ierr = ceed->OperatorCreate(*op); CeedChk(ierr); 566d7b241e6Sjeremylt return 0; 567d7b241e6Sjeremylt } 568d7b241e6Sjeremylt 569d7b241e6Sjeremylt /** 57052d6035fSJeremy L Thompson @brief Create an operator that composes the action of several operators 57152d6035fSJeremy L Thompson 57252d6035fSJeremy L Thompson @param ceed A Ceed object where the CeedOperator will be created 57352d6035fSJeremy L Thompson @param[out] op Address of the variable where the newly created 57452d6035fSJeremy L Thompson Composite CeedOperator will be stored 57552d6035fSJeremy L Thompson 57652d6035fSJeremy L Thompson @return An error code: 0 - success, otherwise - failure 57752d6035fSJeremy L Thompson 5787a982d89SJeremy L. Thompson @ref User 57952d6035fSJeremy L Thompson */ 58052d6035fSJeremy L Thompson int CeedCompositeOperatorCreate(Ceed ceed, CeedOperator *op) { 58152d6035fSJeremy L Thompson int ierr; 58252d6035fSJeremy L Thompson 58352d6035fSJeremy L Thompson if (!ceed->CompositeOperatorCreate) { 58452d6035fSJeremy L Thompson Ceed delegate; 585aefd8378Sjeremylt ierr = CeedGetObjectDelegate(ceed, &delegate, "Operator"); CeedChk(ierr); 58652d6035fSJeremy L Thompson 587250756a7Sjeremylt if (delegate) { 58852d6035fSJeremy L Thompson ierr = CeedCompositeOperatorCreate(delegate, op); CeedChk(ierr); 58952d6035fSJeremy L Thompson return 0; 59052d6035fSJeremy L Thompson } 591250756a7Sjeremylt } 59252d6035fSJeremy L Thompson 59352d6035fSJeremy L Thompson ierr = CeedCalloc(1, op); CeedChk(ierr); 59452d6035fSJeremy L Thompson (*op)->ceed = ceed; 59552d6035fSJeremy L Thompson ceed->refcount++; 59652d6035fSJeremy L Thompson (*op)->composite = true; 59752d6035fSJeremy L Thompson ierr = CeedCalloc(16, &(*op)->suboperators); CeedChk(ierr); 598250756a7Sjeremylt 599250756a7Sjeremylt if (ceed->CompositeOperatorCreate) { 60052d6035fSJeremy L Thompson ierr = ceed->CompositeOperatorCreate(*op); CeedChk(ierr); 601250756a7Sjeremylt } 60252d6035fSJeremy L Thompson return 0; 60352d6035fSJeremy L Thompson } 60452d6035fSJeremy L Thompson 60552d6035fSJeremy L Thompson /** 606b11c1e72Sjeremylt @brief Provide a field to a CeedOperator for use by its CeedQFunction 607d7b241e6Sjeremylt 608d7b241e6Sjeremylt This function is used to specify both active and passive fields to a 609d7b241e6Sjeremylt CeedOperator. For passive fields, a vector @arg v must be provided. Passive 610d7b241e6Sjeremylt fields can inputs or outputs (updated in-place when operator is applied). 611d7b241e6Sjeremylt 612d7b241e6Sjeremylt Active fields must be specified using this function, but their data (in a 613d7b241e6Sjeremylt CeedVector) is passed in CeedOperatorApply(). There can be at most one active 614d7b241e6Sjeremylt input and at most one active output. 615d7b241e6Sjeremylt 6168c91a0c9SJeremy L Thompson @param op CeedOperator on which to provide the field 6178795c945Sjeremylt @param fieldname Name of the field (to be matched with the name used by 6188795c945Sjeremylt CeedQFunction) 619b11c1e72Sjeremylt @param r CeedElemRestriction 6204cc79fe7SJed Brown @param b CeedBasis in which the field resides or @ref CEED_BASIS_COLLOCATED 621b11c1e72Sjeremylt if collocated with quadrature points 6224cc79fe7SJed Brown @param v CeedVector to be used by CeedOperator or @ref CEED_VECTOR_ACTIVE 6234cc79fe7SJed Brown if field is active or @ref CEED_VECTOR_NONE if using 6244cc79fe7SJed Brown @ref CEED_EVAL_WEIGHT in the QFunction 625b11c1e72Sjeremylt 626b11c1e72Sjeremylt @return An error code: 0 - success, otherwise - failure 627dfdf5a53Sjeremylt 6287a982d89SJeremy L. Thompson @ref User 629b11c1e72Sjeremylt **/ 630d7b241e6Sjeremylt int CeedOperatorSetField(CeedOperator op, const char *fieldname, 631a8d32208Sjeremylt CeedElemRestriction r, CeedBasis b, CeedVector v) { 632d7b241e6Sjeremylt int ierr; 63352d6035fSJeremy L Thompson if (op->composite) 634c042f62fSJeremy L Thompson // LCOV_EXCL_START 63552d6035fSJeremy L Thompson return CeedError(op->ceed, 1, "Cannot add field to composite operator."); 636c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 6378b067b84SJed Brown if (!r) 638c042f62fSJeremy L Thompson // LCOV_EXCL_START 639c042f62fSJeremy L Thompson return CeedError(op->ceed, 1, 640c042f62fSJeremy L Thompson "ElemRestriction r for field \"%s\" must be non-NULL.", 641c042f62fSJeremy L Thompson fieldname); 642c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 6438b067b84SJed Brown if (!b) 644c042f62fSJeremy L Thompson // LCOV_EXCL_START 645c042f62fSJeremy L Thompson return CeedError(op->ceed, 1, "Basis b for field \"%s\" must be non-NULL.", 646c042f62fSJeremy L Thompson fieldname); 647c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 6488b067b84SJed Brown if (!v) 649c042f62fSJeremy L Thompson // LCOV_EXCL_START 650c042f62fSJeremy L Thompson return CeedError(op->ceed, 1, "Vector v for field \"%s\" must be non-NULL.", 651c042f62fSJeremy L Thompson fieldname); 652c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 65352d6035fSJeremy L Thompson 654d7b241e6Sjeremylt CeedInt numelements; 655d7b241e6Sjeremylt ierr = CeedElemRestrictionGetNumElements(r, &numelements); CeedChk(ierr); 65615910d16Sjeremylt if (r != CEED_ELEMRESTRICTION_NONE && op->hasrestriction && 65715910d16Sjeremylt op->numelements != numelements) 658c042f62fSJeremy L Thompson // LCOV_EXCL_START 659d7b241e6Sjeremylt return CeedError(op->ceed, 1, 6601d102b48SJeremy L Thompson "ElemRestriction with %d elements incompatible with prior " 6611d102b48SJeremy L Thompson "%d elements", numelements, op->numelements); 662c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 66315910d16Sjeremylt if (r != CEED_ELEMRESTRICTION_NONE) { 664d7b241e6Sjeremylt op->numelements = numelements; 6652cb0afc5Sjeremylt op->hasrestriction = true; // Restriction set, but numelements may be 0 66615910d16Sjeremylt } 667d7b241e6Sjeremylt 668783c99b3SValeria Barra if (b != CEED_BASIS_COLLOCATED) { 669d7b241e6Sjeremylt CeedInt numqpoints; 670d7b241e6Sjeremylt ierr = CeedBasisGetNumQuadraturePoints(b, &numqpoints); CeedChk(ierr); 671d7b241e6Sjeremylt if (op->numqpoints && op->numqpoints != numqpoints) 672c042f62fSJeremy L Thompson // LCOV_EXCL_START 6731d102b48SJeremy L Thompson return CeedError(op->ceed, 1, "Basis with %d quadrature points " 6741d102b48SJeremy L Thompson "incompatible with prior %d points", numqpoints, 6751d102b48SJeremy L Thompson op->numqpoints); 676c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 677d7b241e6Sjeremylt op->numqpoints = numqpoints; 678d7b241e6Sjeremylt } 67915910d16Sjeremylt CeedQFunctionField qfield; 680d1bcdac9Sjeremylt CeedOperatorField *ofield; 681d7b241e6Sjeremylt for (CeedInt i=0; i<op->qf->numinputfields; i++) { 682fe2413ffSjeremylt if (!strcmp(fieldname, (*op->qf->inputfields[i]).fieldname)) { 68315910d16Sjeremylt qfield = op->qf->inputfields[i]; 684d7b241e6Sjeremylt ofield = &op->inputfields[i]; 685d7b241e6Sjeremylt goto found; 686d7b241e6Sjeremylt } 687d7b241e6Sjeremylt } 688d7b241e6Sjeremylt for (CeedInt i=0; i<op->qf->numoutputfields; i++) { 689fe2413ffSjeremylt if (!strcmp(fieldname, (*op->qf->outputfields[i]).fieldname)) { 69015910d16Sjeremylt qfield = op->qf->inputfields[i]; 691d7b241e6Sjeremylt ofield = &op->outputfields[i]; 692d7b241e6Sjeremylt goto found; 693d7b241e6Sjeremylt } 694d7b241e6Sjeremylt } 695c042f62fSJeremy L Thompson // LCOV_EXCL_START 696d7b241e6Sjeremylt return CeedError(op->ceed, 1, "QFunction has no knowledge of field '%s'", 697d7b241e6Sjeremylt fieldname); 698c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 699d7b241e6Sjeremylt found: 70015910d16Sjeremylt if (r == CEED_ELEMRESTRICTION_NONE && qfield->emode != CEED_EVAL_WEIGHT) 701b0d62198Sjeremylt // LCOV_EXCL_START 70215910d16Sjeremylt return CeedError(op->ceed, 1, "CEED_ELEMRESTRICTION_NONE can only be used " 70315910d16Sjeremylt "for a field with eval mode CEED_EVAL_WEIGHT"); 704b0d62198Sjeremylt // LCOV_EXCL_STOP 705fe2413ffSjeremylt ierr = CeedCalloc(1, ofield); CeedChk(ierr); 706fe2413ffSjeremylt (*ofield)->Erestrict = r; 70771352a93Sjeremylt r->refcount += 1; 708fe2413ffSjeremylt (*ofield)->basis = b; 70971352a93Sjeremylt if (b != CEED_BASIS_COLLOCATED) 71071352a93Sjeremylt b->refcount += 1; 711fe2413ffSjeremylt (*ofield)->vec = v; 71271352a93Sjeremylt if (v != CEED_VECTOR_ACTIVE && v != CEED_VECTOR_NONE) 71371352a93Sjeremylt v->refcount += 1; 714d7b241e6Sjeremylt op->nfields += 1; 715d7b241e6Sjeremylt return 0; 716d7b241e6Sjeremylt } 717d7b241e6Sjeremylt 718d7b241e6Sjeremylt /** 71952d6035fSJeremy L Thompson @brief Add a sub-operator to a composite CeedOperator 720288c0443SJeremy L Thompson 72134138859Sjeremylt @param[out] compositeop Composite CeedOperator 72234138859Sjeremylt @param subop Sub-operator CeedOperator 72352d6035fSJeremy L Thompson 72452d6035fSJeremy L Thompson @return An error code: 0 - success, otherwise - failure 72552d6035fSJeremy L Thompson 7267a982d89SJeremy L. Thompson @ref User 72752d6035fSJeremy L Thompson */ 728692c2638Sjeremylt int CeedCompositeOperatorAddSub(CeedOperator compositeop, CeedOperator subop) { 72952d6035fSJeremy L Thompson if (!compositeop->composite) 730c042f62fSJeremy L Thompson // LCOV_EXCL_START 7311d102b48SJeremy L Thompson return CeedError(compositeop->ceed, 1, "CeedOperator is not a composite " 7321d102b48SJeremy L Thompson "operator"); 733c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 73452d6035fSJeremy L Thompson 73552d6035fSJeremy L Thompson if (compositeop->numsub == CEED_COMPOSITE_MAX) 736c042f62fSJeremy L Thompson // LCOV_EXCL_START 7371d102b48SJeremy L Thompson return CeedError(compositeop->ceed, 1, "Cannot add additional suboperators"); 738c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 73952d6035fSJeremy L Thompson 74052d6035fSJeremy L Thompson compositeop->suboperators[compositeop->numsub] = subop; 74152d6035fSJeremy L Thompson subop->refcount++; 74252d6035fSJeremy L Thompson compositeop->numsub++; 74352d6035fSJeremy L Thompson return 0; 74452d6035fSJeremy L Thompson } 74552d6035fSJeremy L Thompson 74652d6035fSJeremy L Thompson /** 7471d102b48SJeremy L Thompson @brief Assemble a linear CeedQFunction associated with a CeedOperator 7481d102b48SJeremy L Thompson 7491d102b48SJeremy L Thompson This returns a CeedVector containing a matrix at each quadrature point 7501d102b48SJeremy L Thompson providing the action of the CeedQFunction associated with the CeedOperator. 7511d102b48SJeremy L Thompson The vector 'assembled' is of shape 7521d102b48SJeremy L Thompson [num_elements, num_input_fields, num_output_fields, num_quad_points] 7531d102b48SJeremy L Thompson and contains column-major matrices representing the action of the 7541d102b48SJeremy L Thompson CeedQFunction for a corresponding quadrature point on an element. Inputs and 7551d102b48SJeremy L Thompson outputs are in the order provided by the user when adding CeedOperator fields. 7564cc79fe7SJed Brown For example, a CeedQFunction with inputs 'u' and 'gradu' and outputs 'gradv' and 7571d102b48SJeremy L Thompson 'v', provided in that order, would result in an assembled QFunction that 7581d102b48SJeremy L Thompson consists of (1 + dim) x (dim + 1) matrices at each quadrature point acting 7591d102b48SJeremy L Thompson on the input [u, du_0, du_1] and producing the output [dv_0, dv_1, v]. 7601d102b48SJeremy L Thompson 7611d102b48SJeremy L Thompson @param op CeedOperator to assemble CeedQFunction 7621d102b48SJeremy L Thompson @param[out] assembled CeedVector to store assembled CeedQFunction at 7631d102b48SJeremy L Thompson quadrature points 7641d102b48SJeremy L Thompson @param[out] rstr CeedElemRestriction for CeedVector containing assembled 7651d102b48SJeremy L Thompson CeedQFunction 7661d102b48SJeremy L Thompson @param request Address of CeedRequest for non-blocking completion, else 7674cc79fe7SJed Brown @ref CEED_REQUEST_IMMEDIATE 7681d102b48SJeremy L Thompson 7691d102b48SJeremy L Thompson @return An error code: 0 - success, otherwise - failure 7701d102b48SJeremy L Thompson 7717a982d89SJeremy L. Thompson @ref User 7721d102b48SJeremy L Thompson **/ 77380ac2e43SJeremy L Thompson int CeedOperatorLinearAssembleQFunction(CeedOperator op, CeedVector *assembled, 7747f823360Sjeremylt CeedElemRestriction *rstr, 7757f823360Sjeremylt CeedRequest *request) { 7761d102b48SJeremy L Thompson int ierr; 7771d102b48SJeremy L Thompson Ceed ceed = op->ceed; 778250756a7Sjeremylt ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 7791d102b48SJeremy L Thompson 7807172caadSjeremylt // Backend version 78180ac2e43SJeremy L Thompson if (op->LinearAssembleQFunction) { 78280ac2e43SJeremy L Thompson ierr = op->LinearAssembleQFunction(op, assembled, rstr, request); 7831d102b48SJeremy L Thompson CeedChk(ierr); 7845107b09fSJeremy L Thompson } else { 7855107b09fSJeremy L Thompson // Fallback to reference Ceed 7865107b09fSJeremy L Thompson if (!op->opfallback) { 7875107b09fSJeremy L Thompson ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 7885107b09fSJeremy L Thompson } 7895107b09fSJeremy L Thompson // Assemble 79080ac2e43SJeremy L Thompson ierr = op->opfallback->LinearAssembleQFunction(op->opfallback, assembled, 7915107b09fSJeremy L Thompson rstr, request); CeedChk(ierr); 7925107b09fSJeremy L Thompson } 793250756a7Sjeremylt 7941d102b48SJeremy L Thompson return 0; 7951d102b48SJeremy L Thompson } 7961d102b48SJeremy L Thompson 7971d102b48SJeremy L Thompson /** 798d965c7a7SJeremy L Thompson @brief Assemble the diagonal of a square linear CeedOperator 799b7ec98d8SJeremy L Thompson 800*9e9210b8SJeremy L Thompson This overwrites a CeedVector with the diagonal of a linear CeedOperator. 801b7ec98d8SJeremy L Thompson 802c04a41a7SJeremy L Thompson Note: Currently only non-composite CeedOperators with a single field and 803c04a41a7SJeremy L Thompson composite CeedOperators with single field sub-operators are supported. 804d965c7a7SJeremy L Thompson 805b7ec98d8SJeremy L Thompson @param op CeedOperator to assemble CeedQFunction 806b7ec98d8SJeremy L Thompson @param[out] assembled CeedVector to store assembled CeedOperator diagonal 807b7ec98d8SJeremy L Thompson @param request Address of CeedRequest for non-blocking completion, else 8084cc79fe7SJed Brown @ref CEED_REQUEST_IMMEDIATE 809b7ec98d8SJeremy L Thompson 810b7ec98d8SJeremy L Thompson @return An error code: 0 - success, otherwise - failure 811b7ec98d8SJeremy L Thompson 8127a982d89SJeremy L. Thompson @ref User 813b7ec98d8SJeremy L Thompson **/ 8142bba3ffaSJeremy L Thompson int CeedOperatorLinearAssembleDiagonal(CeedOperator op, CeedVector assembled, 8157f823360Sjeremylt CeedRequest *request) { 816b7ec98d8SJeremy L Thompson int ierr; 817b7ec98d8SJeremy L Thompson Ceed ceed = op->ceed; 818250756a7Sjeremylt ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 819b7ec98d8SJeremy L Thompson 820b7ec98d8SJeremy L Thompson // Use backend version, if available 82180ac2e43SJeremy L Thompson if (op->LinearAssembleDiagonal) { 82280ac2e43SJeremy L Thompson ierr = op->LinearAssembleDiagonal(op, assembled, request); CeedChk(ierr); 823*9e9210b8SJeremy L Thompson } else if (op->LinearAssembleAddDiagonal) { 824*9e9210b8SJeremy L Thompson ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 825*9e9210b8SJeremy L Thompson return CeedOperatorLinearAssembleAddDiagonal(op, assembled, request); 826*9e9210b8SJeremy L Thompson } else { 827*9e9210b8SJeremy L Thompson // Fallback to reference Ceed 828*9e9210b8SJeremy L Thompson if (!op->opfallback) { 829*9e9210b8SJeremy L Thompson ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 830*9e9210b8SJeremy L Thompson } 831*9e9210b8SJeremy L Thompson // Assemble 832*9e9210b8SJeremy L Thompson if (op->opfallback->LinearAssembleDiagonal) { 833*9e9210b8SJeremy L Thompson ierr = op->opfallback->LinearAssembleDiagonal(op->opfallback, assembled, 834*9e9210b8SJeremy L Thompson request); CeedChk(ierr); 835*9e9210b8SJeremy L Thompson } else { 836*9e9210b8SJeremy L Thompson ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 837*9e9210b8SJeremy L Thompson return CeedOperatorLinearAssembleAddDiagonal(op, assembled, request); 838*9e9210b8SJeremy L Thompson } 839*9e9210b8SJeremy L Thompson } 840*9e9210b8SJeremy L Thompson 841*9e9210b8SJeremy L Thompson return 0; 842*9e9210b8SJeremy L Thompson } 843*9e9210b8SJeremy L Thompson 844*9e9210b8SJeremy L Thompson /** 845*9e9210b8SJeremy L Thompson @brief Assemble the diagonal of a square linear CeedOperator 846*9e9210b8SJeremy L Thompson 847*9e9210b8SJeremy L Thompson This sums into a CeedVector the diagonal of a linear CeedOperator. 848*9e9210b8SJeremy L Thompson 849*9e9210b8SJeremy L Thompson Note: Currently only non-composite CeedOperators with a single field and 850*9e9210b8SJeremy L Thompson composite CeedOperators with single field sub-operators are supported. 851*9e9210b8SJeremy L Thompson 852*9e9210b8SJeremy L Thompson @param op CeedOperator to assemble CeedQFunction 853*9e9210b8SJeremy L Thompson @param[out] assembled CeedVector to store assembled CeedOperator diagonal 854*9e9210b8SJeremy L Thompson @param request Address of CeedRequest for non-blocking completion, else 855*9e9210b8SJeremy L Thompson @ref CEED_REQUEST_IMMEDIATE 856*9e9210b8SJeremy L Thompson 857*9e9210b8SJeremy L Thompson @return An error code: 0 - success, otherwise - failure 858*9e9210b8SJeremy L Thompson 859*9e9210b8SJeremy L Thompson @ref User 860*9e9210b8SJeremy L Thompson **/ 861*9e9210b8SJeremy L Thompson int CeedOperatorLinearAssembleAddDiagonal(CeedOperator op, CeedVector assembled, 862*9e9210b8SJeremy L Thompson CeedRequest *request) { 863*9e9210b8SJeremy L Thompson int ierr; 864*9e9210b8SJeremy L Thompson Ceed ceed = op->ceed; 865*9e9210b8SJeremy L Thompson ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 866*9e9210b8SJeremy L Thompson 867*9e9210b8SJeremy L Thompson // Use backend version, if available 868*9e9210b8SJeremy L Thompson if (op->LinearAssembleAddDiagonal) { 869*9e9210b8SJeremy L Thompson ierr = op->LinearAssembleAddDiagonal(op, assembled, request); CeedChk(ierr); 8707172caadSjeremylt } else { 8717172caadSjeremylt // Fallback to reference Ceed 8727172caadSjeremylt if (!op->opfallback) { 8737172caadSjeremylt ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 874b7ec98d8SJeremy L Thompson } 8757172caadSjeremylt // Assemble 8762bba3ffaSJeremy L Thompson ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 877*9e9210b8SJeremy L Thompson ierr = op->opfallback->LinearAssembleAddDiagonal(op->opfallback, assembled, 8787172caadSjeremylt request); CeedChk(ierr); 879b7ec98d8SJeremy L Thompson } 880b7ec98d8SJeremy L Thompson 881b7ec98d8SJeremy L Thompson return 0; 882b7ec98d8SJeremy L Thompson } 883b7ec98d8SJeremy L Thompson 884b7ec98d8SJeremy L Thompson /** 885d965c7a7SJeremy L Thompson @brief Assemble the point block diagonal of a square linear CeedOperator 886d965c7a7SJeremy L Thompson 887*9e9210b8SJeremy L Thompson This overwrites a CeedVector with the point block diagonal of a linear 888d965c7a7SJeremy L Thompson CeedOperator. 889d965c7a7SJeremy L Thompson 890c04a41a7SJeremy L Thompson Note: Currently only non-composite CeedOperators with a single field and 891c04a41a7SJeremy L Thompson composite CeedOperators with single field sub-operators are supported. 892d965c7a7SJeremy L Thompson 893d965c7a7SJeremy L Thompson @param op CeedOperator to assemble CeedQFunction 894d965c7a7SJeremy L Thompson @param[out] assembled CeedVector to store assembled CeedOperator point block 895d965c7a7SJeremy L Thompson diagonal, provided in row-major form with an 896d965c7a7SJeremy L Thompson @a ncomp * @a ncomp block at each node. The dimensions 897d965c7a7SJeremy L Thompson of this vector are derived from the active vector 898d965c7a7SJeremy L Thompson for the CeedOperator. The array has shape 899d965c7a7SJeremy L Thompson [nodes, component out, component in]. 900d965c7a7SJeremy L Thompson @param request Address of CeedRequest for non-blocking completion, else 901d965c7a7SJeremy L Thompson CEED_REQUEST_IMMEDIATE 902d965c7a7SJeremy L Thompson 903d965c7a7SJeremy L Thompson @return An error code: 0 - success, otherwise - failure 904d965c7a7SJeremy L Thompson 905d965c7a7SJeremy L Thompson @ref User 906d965c7a7SJeremy L Thompson **/ 90780ac2e43SJeremy L Thompson int CeedOperatorLinearAssemblePointBlockDiagonal(CeedOperator op, 9082bba3ffaSJeremy L Thompson CeedVector assembled, CeedRequest *request) { 909d965c7a7SJeremy L Thompson int ierr; 910d965c7a7SJeremy L Thompson Ceed ceed = op->ceed; 911d965c7a7SJeremy L Thompson ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 912d965c7a7SJeremy L Thompson 913d965c7a7SJeremy L Thompson // Use backend version, if available 91480ac2e43SJeremy L Thompson if (op->LinearAssemblePointBlockDiagonal) { 91580ac2e43SJeremy L Thompson ierr = op->LinearAssemblePointBlockDiagonal(op, assembled, request); 916d965c7a7SJeremy L Thompson CeedChk(ierr); 917*9e9210b8SJeremy L Thompson } else if (op->LinearAssembleAddPointBlockDiagonal) { 918*9e9210b8SJeremy L Thompson ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 919*9e9210b8SJeremy L Thompson return CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, 920*9e9210b8SJeremy L Thompson request); 921d965c7a7SJeremy L Thompson } else { 922d965c7a7SJeremy L Thompson // Fallback to reference Ceed 923d965c7a7SJeremy L Thompson if (!op->opfallback) { 924d965c7a7SJeremy L Thompson ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 925d965c7a7SJeremy L Thompson } 926d965c7a7SJeremy L Thompson // Assemble 927*9e9210b8SJeremy L Thompson if (op->opfallback->LinearAssemblePointBlockDiagonal) { 92880ac2e43SJeremy L Thompson ierr = op->opfallback->LinearAssemblePointBlockDiagonal(op->opfallback, 929d965c7a7SJeremy L Thompson assembled, request); CeedChk(ierr); 930*9e9210b8SJeremy L Thompson } else { 931*9e9210b8SJeremy L Thompson ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 932*9e9210b8SJeremy L Thompson return CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, 933*9e9210b8SJeremy L Thompson request); 934*9e9210b8SJeremy L Thompson } 935*9e9210b8SJeremy L Thompson } 936*9e9210b8SJeremy L Thompson 937*9e9210b8SJeremy L Thompson return 0; 938*9e9210b8SJeremy L Thompson } 939*9e9210b8SJeremy L Thompson 940*9e9210b8SJeremy L Thompson /** 941*9e9210b8SJeremy L Thompson @brief Assemble the point block diagonal of a square linear CeedOperator 942*9e9210b8SJeremy L Thompson 943*9e9210b8SJeremy L Thompson This sums into a CeedVector with the point block diagonal of a linear 944*9e9210b8SJeremy L Thompson CeedOperator. 945*9e9210b8SJeremy L Thompson 946*9e9210b8SJeremy L Thompson Note: Currently only non-composite CeedOperators with a single field and 947*9e9210b8SJeremy L Thompson composite CeedOperators with single field sub-operators are supported. 948*9e9210b8SJeremy L Thompson 949*9e9210b8SJeremy L Thompson @param op CeedOperator to assemble CeedQFunction 950*9e9210b8SJeremy L Thompson @param[out] assembled CeedVector to store assembled CeedOperator point block 951*9e9210b8SJeremy L Thompson diagonal, provided in row-major form with an 952*9e9210b8SJeremy L Thompson @a ncomp * @a ncomp block at each node. The dimensions 953*9e9210b8SJeremy L Thompson of this vector are derived from the active vector 954*9e9210b8SJeremy L Thompson for the CeedOperator. The array has shape 955*9e9210b8SJeremy L Thompson [nodes, component out, component in]. 956*9e9210b8SJeremy L Thompson @param request Address of CeedRequest for non-blocking completion, else 957*9e9210b8SJeremy L Thompson CEED_REQUEST_IMMEDIATE 958*9e9210b8SJeremy L Thompson 959*9e9210b8SJeremy L Thompson @return An error code: 0 - success, otherwise - failure 960*9e9210b8SJeremy L Thompson 961*9e9210b8SJeremy L Thompson @ref User 962*9e9210b8SJeremy L Thompson **/ 963*9e9210b8SJeremy L Thompson int CeedOperatorLinearAssembleAddPointBlockDiagonal(CeedOperator op, 964*9e9210b8SJeremy L Thompson CeedVector assembled, CeedRequest *request) { 965*9e9210b8SJeremy L Thompson int ierr; 966*9e9210b8SJeremy L Thompson Ceed ceed = op->ceed; 967*9e9210b8SJeremy L Thompson ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 968*9e9210b8SJeremy L Thompson 969*9e9210b8SJeremy L Thompson // Use backend version, if available 970*9e9210b8SJeremy L Thompson if (op->LinearAssembleAddPointBlockDiagonal) { 971*9e9210b8SJeremy L Thompson ierr = op->LinearAssembleAddPointBlockDiagonal(op, assembled, request); 972*9e9210b8SJeremy L Thompson CeedChk(ierr); 973*9e9210b8SJeremy L Thompson } else { 974*9e9210b8SJeremy L Thompson // Fallback to reference Ceed 975*9e9210b8SJeremy L Thompson if (!op->opfallback) { 976*9e9210b8SJeremy L Thompson ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 977*9e9210b8SJeremy L Thompson } 978*9e9210b8SJeremy L Thompson // Assemble 979*9e9210b8SJeremy L Thompson ierr = op->opfallback->LinearAssembleAddPointBlockDiagonal(op->opfallback, 980*9e9210b8SJeremy L Thompson assembled, request); CeedChk(ierr); 981d965c7a7SJeremy L Thompson } 982d965c7a7SJeremy L Thompson 983d965c7a7SJeremy L Thompson return 0; 984d965c7a7SJeremy L Thompson } 985d965c7a7SJeremy L Thompson 986d965c7a7SJeremy L Thompson /** 9873bd813ffSjeremylt @brief Build a FDM based approximate inverse for each element for a 9883bd813ffSjeremylt CeedOperator 9893bd813ffSjeremylt 9903bd813ffSjeremylt This returns a CeedOperator and CeedVector to apply a Fast Diagonalization 9913bd813ffSjeremylt Method based approximate inverse. This function obtains the simultaneous 9923bd813ffSjeremylt diagonalization for the 1D mass and Laplacian operators, 9933bd813ffSjeremylt M = V^T V, K = V^T S V. 9943bd813ffSjeremylt The assembled QFunction is used to modify the eigenvalues from simultaneous 9953bd813ffSjeremylt diagonalization and obtain an approximate inverse of the form 9963bd813ffSjeremylt V^T S^hat V. The CeedOperator must be linear and non-composite. The 9973bd813ffSjeremylt associated CeedQFunction must therefore also be linear. 9983bd813ffSjeremylt 9993bd813ffSjeremylt @param op CeedOperator to create element inverses 10003bd813ffSjeremylt @param[out] fdminv CeedOperator to apply the action of a FDM based inverse 10013bd813ffSjeremylt for each element 10023bd813ffSjeremylt @param request Address of CeedRequest for non-blocking completion, else 10034cc79fe7SJed Brown @ref CEED_REQUEST_IMMEDIATE 10043bd813ffSjeremylt 10053bd813ffSjeremylt @return An error code: 0 - success, otherwise - failure 10063bd813ffSjeremylt 10077a982d89SJeremy L. Thompson @ref Backend 10083bd813ffSjeremylt **/ 10093bd813ffSjeremylt int CeedOperatorCreateFDMElementInverse(CeedOperator op, CeedOperator *fdminv, 10103bd813ffSjeremylt CeedRequest *request) { 10113bd813ffSjeremylt int ierr; 10123bd813ffSjeremylt Ceed ceed = op->ceed; 1013713f43c3Sjeremylt ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 10143bd813ffSjeremylt 1015713f43c3Sjeremylt // Use backend version, if available 1016713f43c3Sjeremylt if (op->CreateFDMElementInverse) { 1017713f43c3Sjeremylt ierr = op->CreateFDMElementInverse(op, fdminv, request); CeedChk(ierr); 1018713f43c3Sjeremylt } else { 1019713f43c3Sjeremylt // Fallback to reference Ceed 1020713f43c3Sjeremylt if (!op->opfallback) { 1021713f43c3Sjeremylt ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 10223bd813ffSjeremylt } 1023713f43c3Sjeremylt // Assemble 1024713f43c3Sjeremylt ierr = op->opfallback->CreateFDMElementInverse(op->opfallback, fdminv, 10253bd813ffSjeremylt request); CeedChk(ierr); 10263bd813ffSjeremylt } 10273bd813ffSjeremylt 10283bd813ffSjeremylt return 0; 10293bd813ffSjeremylt } 10303bd813ffSjeremylt 10317a982d89SJeremy L. Thompson /** 10327a982d89SJeremy L. Thompson @brief View a CeedOperator 10337a982d89SJeremy L. Thompson 10347a982d89SJeremy L. Thompson @param[in] op CeedOperator to view 10357a982d89SJeremy L. Thompson @param[in] stream Stream to write; typically stdout/stderr or a file 10367a982d89SJeremy L. Thompson 10377a982d89SJeremy L. Thompson @return Error code: 0 - success, otherwise - failure 10387a982d89SJeremy L. Thompson 10397a982d89SJeremy L. Thompson @ref User 10407a982d89SJeremy L. Thompson **/ 10417a982d89SJeremy L. Thompson int CeedOperatorView(CeedOperator op, FILE *stream) { 10427a982d89SJeremy L. Thompson int ierr; 10437a982d89SJeremy L. Thompson 10447a982d89SJeremy L. Thompson if (op->composite) { 10457a982d89SJeremy L. Thompson fprintf(stream, "Composite CeedOperator\n"); 10467a982d89SJeremy L. Thompson 10477a982d89SJeremy L. Thompson for (CeedInt i=0; i<op->numsub; i++) { 10487a982d89SJeremy L. Thompson fprintf(stream, " SubOperator [%d]:\n", i); 10497a982d89SJeremy L. Thompson ierr = CeedOperatorSingleView(op->suboperators[i], 1, stream); 10507a982d89SJeremy L. Thompson CeedChk(ierr); 10517a982d89SJeremy L. Thompson } 10527a982d89SJeremy L. Thompson } else { 10537a982d89SJeremy L. Thompson fprintf(stream, "CeedOperator\n"); 10547a982d89SJeremy L. Thompson ierr = CeedOperatorSingleView(op, 0, stream); CeedChk(ierr); 10557a982d89SJeremy L. Thompson } 10567a982d89SJeremy L. Thompson 10577a982d89SJeremy L. Thompson return 0; 10587a982d89SJeremy L. Thompson } 10593bd813ffSjeremylt 10603bd813ffSjeremylt /** 10613bd813ffSjeremylt @brief Apply CeedOperator to a vector 1062d7b241e6Sjeremylt 1063d7b241e6Sjeremylt This computes the action of the operator on the specified (active) input, 1064d7b241e6Sjeremylt yielding its (active) output. All inputs and outputs must be specified using 1065d7b241e6Sjeremylt CeedOperatorSetField(). 1066d7b241e6Sjeremylt 1067d7b241e6Sjeremylt @param op CeedOperator to apply 10684cc79fe7SJed Brown @param[in] in CeedVector containing input state or @ref CEED_VECTOR_NONE if 106934138859Sjeremylt there are no active inputs 1070b11c1e72Sjeremylt @param[out] out CeedVector to store result of applying operator (must be 10714cc79fe7SJed Brown distinct from @a in) or @ref CEED_VECTOR_NONE if there are no 107234138859Sjeremylt active outputs 1073d7b241e6Sjeremylt @param request Address of CeedRequest for non-blocking completion, else 10744cc79fe7SJed Brown @ref CEED_REQUEST_IMMEDIATE 1075b11c1e72Sjeremylt 1076b11c1e72Sjeremylt @return An error code: 0 - success, otherwise - failure 1077dfdf5a53Sjeremylt 10787a982d89SJeremy L. Thompson @ref User 1079b11c1e72Sjeremylt **/ 1080692c2638Sjeremylt int CeedOperatorApply(CeedOperator op, CeedVector in, CeedVector out, 1081692c2638Sjeremylt CeedRequest *request) { 1082d7b241e6Sjeremylt int ierr; 1083d7b241e6Sjeremylt Ceed ceed = op->ceed; 1084250756a7Sjeremylt ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 1085d7b241e6Sjeremylt 1086250756a7Sjeremylt if (op->numelements) { 1087250756a7Sjeremylt // Standard Operator 1088cae8b89aSjeremylt if (op->Apply) { 1089250756a7Sjeremylt ierr = op->Apply(op, in, out, request); CeedChk(ierr); 1090cae8b89aSjeremylt } else { 1091cae8b89aSjeremylt // Zero all output vectors 1092250756a7Sjeremylt CeedQFunction qf = op->qf; 1093cae8b89aSjeremylt for (CeedInt i=0; i<qf->numoutputfields; i++) { 1094cae8b89aSjeremylt CeedVector vec = op->outputfields[i]->vec; 1095cae8b89aSjeremylt if (vec == CEED_VECTOR_ACTIVE) 1096cae8b89aSjeremylt vec = out; 1097cae8b89aSjeremylt if (vec != CEED_VECTOR_NONE) { 1098cae8b89aSjeremylt ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr); 1099cae8b89aSjeremylt } 1100cae8b89aSjeremylt } 1101250756a7Sjeremylt // Apply 1102250756a7Sjeremylt ierr = op->ApplyAdd(op, in, out, request); CeedChk(ierr); 1103250756a7Sjeremylt } 1104250756a7Sjeremylt } else if (op->composite) { 1105250756a7Sjeremylt // Composite Operator 1106250756a7Sjeremylt if (op->ApplyComposite) { 1107250756a7Sjeremylt ierr = op->ApplyComposite(op, in, out, request); CeedChk(ierr); 1108250756a7Sjeremylt } else { 1109250756a7Sjeremylt CeedInt numsub; 1110250756a7Sjeremylt ierr = CeedOperatorGetNumSub(op, &numsub); CeedChk(ierr); 1111250756a7Sjeremylt CeedOperator *suboperators; 1112250756a7Sjeremylt ierr = CeedOperatorGetSubList(op, &suboperators); CeedChk(ierr); 1113250756a7Sjeremylt 1114250756a7Sjeremylt // Zero all output vectors 1115250756a7Sjeremylt if (out != CEED_VECTOR_NONE) { 1116cae8b89aSjeremylt ierr = CeedVectorSetValue(out, 0.0); CeedChk(ierr); 1117cae8b89aSjeremylt } 1118250756a7Sjeremylt for (CeedInt i=0; i<numsub; i++) { 1119250756a7Sjeremylt for (CeedInt j=0; j<suboperators[i]->qf->numoutputfields; j++) { 1120250756a7Sjeremylt CeedVector vec = suboperators[i]->outputfields[j]->vec; 1121250756a7Sjeremylt if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) { 1122250756a7Sjeremylt ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr); 1123250756a7Sjeremylt } 1124250756a7Sjeremylt } 1125250756a7Sjeremylt } 1126250756a7Sjeremylt // Apply 1127250756a7Sjeremylt for (CeedInt i=0; i<op->numsub; i++) { 1128250756a7Sjeremylt ierr = CeedOperatorApplyAdd(op->suboperators[i], in, out, request); 1129cae8b89aSjeremylt CeedChk(ierr); 1130cae8b89aSjeremylt } 1131cae8b89aSjeremylt } 1132250756a7Sjeremylt } 1133250756a7Sjeremylt 1134cae8b89aSjeremylt return 0; 1135cae8b89aSjeremylt } 1136cae8b89aSjeremylt 1137cae8b89aSjeremylt /** 1138cae8b89aSjeremylt @brief Apply CeedOperator to a vector and add result to output vector 1139cae8b89aSjeremylt 1140cae8b89aSjeremylt This computes the action of the operator on the specified (active) input, 1141cae8b89aSjeremylt yielding its (active) output. All inputs and outputs must be specified using 1142cae8b89aSjeremylt CeedOperatorSetField(). 1143cae8b89aSjeremylt 1144cae8b89aSjeremylt @param op CeedOperator to apply 1145cae8b89aSjeremylt @param[in] in CeedVector containing input state or NULL if there are no 1146cae8b89aSjeremylt active inputs 1147cae8b89aSjeremylt @param[out] out CeedVector to sum in result of applying operator (must be 1148cae8b89aSjeremylt distinct from @a in) or NULL if there are no active outputs 1149cae8b89aSjeremylt @param request Address of CeedRequest for non-blocking completion, else 11504cc79fe7SJed Brown @ref CEED_REQUEST_IMMEDIATE 1151cae8b89aSjeremylt 1152cae8b89aSjeremylt @return An error code: 0 - success, otherwise - failure 1153cae8b89aSjeremylt 11547a982d89SJeremy L. Thompson @ref User 1155cae8b89aSjeremylt **/ 1156cae8b89aSjeremylt int CeedOperatorApplyAdd(CeedOperator op, CeedVector in, CeedVector out, 1157cae8b89aSjeremylt CeedRequest *request) { 1158cae8b89aSjeremylt int ierr; 1159cae8b89aSjeremylt Ceed ceed = op->ceed; 1160250756a7Sjeremylt ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 1161cae8b89aSjeremylt 1162250756a7Sjeremylt if (op->numelements) { 1163250756a7Sjeremylt // Standard Operator 1164250756a7Sjeremylt ierr = op->ApplyAdd(op, in, out, request); CeedChk(ierr); 1165250756a7Sjeremylt } else if (op->composite) { 1166250756a7Sjeremylt // Composite Operator 1167250756a7Sjeremylt if (op->ApplyAddComposite) { 1168250756a7Sjeremylt ierr = op->ApplyAddComposite(op, in, out, request); CeedChk(ierr); 1169cae8b89aSjeremylt } else { 1170250756a7Sjeremylt CeedInt numsub; 1171250756a7Sjeremylt ierr = CeedOperatorGetNumSub(op, &numsub); CeedChk(ierr); 1172250756a7Sjeremylt CeedOperator *suboperators; 1173250756a7Sjeremylt ierr = CeedOperatorGetSubList(op, &suboperators); CeedChk(ierr); 1174250756a7Sjeremylt 1175250756a7Sjeremylt for (CeedInt i=0; i<numsub; i++) { 1176250756a7Sjeremylt ierr = CeedOperatorApplyAdd(suboperators[i], in, out, request); 1177cae8b89aSjeremylt CeedChk(ierr); 11781d7d2407SJeremy L Thompson } 1179250756a7Sjeremylt } 1180250756a7Sjeremylt } 1181250756a7Sjeremylt 1182d7b241e6Sjeremylt return 0; 1183d7b241e6Sjeremylt } 1184d7b241e6Sjeremylt 1185d7b241e6Sjeremylt /** 1186b11c1e72Sjeremylt @brief Destroy a CeedOperator 1187d7b241e6Sjeremylt 1188d7b241e6Sjeremylt @param op CeedOperator to destroy 1189b11c1e72Sjeremylt 1190b11c1e72Sjeremylt @return An error code: 0 - success, otherwise - failure 1191dfdf5a53Sjeremylt 11927a982d89SJeremy L. Thompson @ref User 1193b11c1e72Sjeremylt **/ 1194d7b241e6Sjeremylt int CeedOperatorDestroy(CeedOperator *op) { 1195d7b241e6Sjeremylt int ierr; 1196d7b241e6Sjeremylt 1197d7b241e6Sjeremylt if (!*op || --(*op)->refcount > 0) return 0; 1198d7b241e6Sjeremylt if ((*op)->Destroy) { 1199d7b241e6Sjeremylt ierr = (*op)->Destroy(*op); CeedChk(ierr); 1200d7b241e6Sjeremylt } 1201fe2413ffSjeremylt ierr = CeedDestroy(&(*op)->ceed); CeedChk(ierr); 1202fe2413ffSjeremylt // Free fields 12031d102b48SJeremy L Thompson for (int i=0; i<(*op)->nfields; i++) 120452d6035fSJeremy L Thompson if ((*op)->inputfields[i]) { 120515910d16Sjeremylt if ((*op)->inputfields[i]->Erestrict != CEED_ELEMRESTRICTION_NONE) { 120671352a93Sjeremylt ierr = CeedElemRestrictionDestroy(&(*op)->inputfields[i]->Erestrict); 120771352a93Sjeremylt CeedChk(ierr); 120815910d16Sjeremylt } 120971352a93Sjeremylt if ((*op)->inputfields[i]->basis != CEED_BASIS_COLLOCATED) { 121071352a93Sjeremylt ierr = CeedBasisDestroy(&(*op)->inputfields[i]->basis); CeedChk(ierr); 121171352a93Sjeremylt } 121271352a93Sjeremylt if ((*op)->inputfields[i]->vec != CEED_VECTOR_ACTIVE && 121371352a93Sjeremylt (*op)->inputfields[i]->vec != CEED_VECTOR_NONE ) { 121471352a93Sjeremylt ierr = CeedVectorDestroy(&(*op)->inputfields[i]->vec); CeedChk(ierr); 121571352a93Sjeremylt } 1216fe2413ffSjeremylt ierr = CeedFree(&(*op)->inputfields[i]); CeedChk(ierr); 1217fe2413ffSjeremylt } 12181d102b48SJeremy L Thompson for (int i=0; i<(*op)->nfields; i++) 1219d0eb30a5Sjeremylt if ((*op)->outputfields[i]) { 122071352a93Sjeremylt ierr = CeedElemRestrictionDestroy(&(*op)->outputfields[i]->Erestrict); 122171352a93Sjeremylt CeedChk(ierr); 122271352a93Sjeremylt if ((*op)->outputfields[i]->basis != CEED_BASIS_COLLOCATED) { 122371352a93Sjeremylt ierr = CeedBasisDestroy(&(*op)->outputfields[i]->basis); CeedChk(ierr); 122471352a93Sjeremylt } 122571352a93Sjeremylt if ((*op)->outputfields[i]->vec != CEED_VECTOR_ACTIVE && 122671352a93Sjeremylt (*op)->outputfields[i]->vec != CEED_VECTOR_NONE ) { 122771352a93Sjeremylt ierr = CeedVectorDestroy(&(*op)->outputfields[i]->vec); CeedChk(ierr); 122871352a93Sjeremylt } 1229fe2413ffSjeremylt ierr = CeedFree(&(*op)->outputfields[i]); CeedChk(ierr); 1230fe2413ffSjeremylt } 123152d6035fSJeremy L Thompson // Destroy suboperators 12321d102b48SJeremy L Thompson for (int i=0; i<(*op)->numsub; i++) 123352d6035fSJeremy L Thompson if ((*op)->suboperators[i]) { 123452d6035fSJeremy L Thompson ierr = CeedOperatorDestroy(&(*op)->suboperators[i]); CeedChk(ierr); 123552d6035fSJeremy L Thompson } 1236d7b241e6Sjeremylt ierr = CeedQFunctionDestroy(&(*op)->qf); CeedChk(ierr); 1237d7b241e6Sjeremylt ierr = CeedQFunctionDestroy(&(*op)->dqf); CeedChk(ierr); 1238d7b241e6Sjeremylt ierr = CeedQFunctionDestroy(&(*op)->dqfT); CeedChk(ierr); 1239fe2413ffSjeremylt 12405107b09fSJeremy L Thompson // Destroy fallback 12415107b09fSJeremy L Thompson if ((*op)->opfallback) { 12425107b09fSJeremy L Thompson ierr = (*op)->qffallback->Destroy((*op)->qffallback); CeedChk(ierr); 12435107b09fSJeremy L Thompson ierr = CeedFree(&(*op)->qffallback); CeedChk(ierr); 12445107b09fSJeremy L Thompson ierr = (*op)->opfallback->Destroy((*op)->opfallback); CeedChk(ierr); 12455107b09fSJeremy L Thompson ierr = CeedFree(&(*op)->opfallback); CeedChk(ierr); 12465107b09fSJeremy L Thompson } 12475107b09fSJeremy L Thompson 1248fe2413ffSjeremylt ierr = CeedFree(&(*op)->inputfields); CeedChk(ierr); 1249fe2413ffSjeremylt ierr = CeedFree(&(*op)->outputfields); CeedChk(ierr); 125052d6035fSJeremy L Thompson ierr = CeedFree(&(*op)->suboperators); CeedChk(ierr); 1251d7b241e6Sjeremylt ierr = CeedFree(op); CeedChk(ierr); 1252d7b241e6Sjeremylt return 0; 1253d7b241e6Sjeremylt } 1254d7b241e6Sjeremylt 1255d7b241e6Sjeremylt /// @} 1256