1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 #include <ceed-impl.h> 18 #include <string.h> 19 20 /** 21 @file 22 Implementation of public CeedOperator interfaces 23 24 @defgroup CeedOperator CeedOperator: composed FE-type operations on vectors 25 @{ 26 */ 27 28 /** 29 @brief Create an operator from element restriction, basis, and QFunction 30 31 @param ceed A Ceed object where the CeedOperator will be created 32 @param qf QFunction defining the action of the operator at quadrature points 33 @param dqf QFunction defining the action of the Jacobian of @a qf (or NULL) 34 @param dqfT QFunction defining the action of the transpose of the Jacobian 35 of @a qf (or NULL) 36 @param[out] op Address of the variable where the newly created 37 CeedOperator will be stored 38 39 @return An error code: 0 - success, otherwise - failure 40 */ 41 int CeedOperatorCreate(Ceed ceed, CeedQFunction qf, CeedQFunction dqf, 42 CeedQFunction dqfT, CeedOperator *op) { 43 int ierr; 44 45 if (!ceed->OperatorCreate) return CeedError(ceed, 1, 46 "Backend does not support OperatorCreate"); 47 ierr = CeedCalloc(1,op); CeedChk(ierr); 48 (*op)->ceed = ceed; 49 ceed->refcount++; 50 (*op)->refcount = 1; 51 (*op)->qf = qf; 52 qf->refcount++; 53 (*op)->dqf = dqf; 54 if (dqf) dqf->refcount++; 55 (*op)->dqfT = dqfT; 56 if (dqfT) dqfT->refcount++; 57 ierr = ceed->OperatorCreate(*op); CeedChk(ierr); 58 return 0; 59 } 60 61 /** 62 @brief Provide a field to a CeedOperator for use by its CeedQFunction 63 64 This function is used to specify both active and passive fields to a 65 CeedOperator. For passive fields, a vector @arg v must be provided. Passive 66 fields can inputs or outputs (updated in-place when operator is applied). 67 68 Active fields must be specified using this function, but their data (in a 69 CeedVector) is passed in CeedOperatorApply(). There can be at most one active 70 input and at most one active output. 71 72 @param op Ceedoperator on which to provide the field 73 @param fieldname Name of the field (to be matched with the name used by CeedQFunction) 74 @param r CeedElemRestriction 75 @param b CeedBasis in which the field resides or CEED_BASIS_COLOCATED 76 if collocated with quadrature points 77 @param v CeedVector to be used by CeedOperator or CEED_VECTOR_ACTIVE 78 if field is active or CEED_VECTOR_NONE if using 79 CEED_EVAL_WEIGHT in the qfunction 80 81 @return An error code: 0 - success, otherwise - failure 82 **/ 83 int CeedOperatorSetField(CeedOperator op, const char *fieldname, 84 CeedElemRestriction r, CeedBasis b, 85 CeedVector v) { 86 int ierr; 87 CeedInt numelements; 88 ierr = CeedElemRestrictionGetNumElements(r, &numelements); CeedChk(ierr); 89 if (op->numelements && op->numelements != numelements) 90 return CeedError(op->ceed, 1, 91 "ElemRestriction with %d elements incompatible with prior %d elements", 92 numelements, op->numelements); 93 op->numelements = numelements; 94 95 if (b != CEED_BASIS_COLOCATED) { 96 CeedInt numqpoints; 97 ierr = CeedBasisGetNumQuadraturePoints(b, &numqpoints); CeedChk(ierr); 98 if (op->numqpoints && op->numqpoints != numqpoints) 99 return CeedError(op->ceed, 1, 100 "Basis with %d quadrature points incompatible with prior %d points", 101 numqpoints, op->numqpoints); 102 op->numqpoints = numqpoints; 103 } 104 struct CeedOperatorField *ofield; 105 for (CeedInt i=0; i<op->qf->numinputfields; i++) { 106 if (!strcmp(fieldname, op->qf->inputfields[i].fieldname)) { 107 ofield = &op->inputfields[i]; 108 goto found; 109 } 110 } 111 for (CeedInt i=0; i<op->qf->numoutputfields; i++) { 112 if (!strcmp(fieldname, op->qf->outputfields[i].fieldname)) { 113 ofield = &op->outputfields[i]; 114 goto found; 115 } 116 } 117 return CeedError(op->ceed, 1, "QFunction has no knowledge of field '%s'", 118 fieldname); 119 found: 120 ofield->Erestrict = r; 121 ofield->basis = b; 122 ofield->vec = v; 123 op->nfields += 1; 124 return 0; 125 } 126 127 /** 128 @brief Apply CeedOperator to a vector 129 130 This computes the action of the operator on the specified (active) input, 131 yielding its (active) output. All inputs and outputs must be specified using 132 CeedOperatorSetField(). 133 134 @param op CeedOperator to apply 135 @param[in] in CeedVector containing input state or NULL if there are no 136 active inputs 137 @param[out] out CeedVector to store result of applying operator (must be 138 distinct from @a in) or NULL if there are no active outputs 139 @param request Address of CeedRequest for non-blocking completion, else 140 CEED_REQUEST_IMMEDIATE 141 142 @return An error code: 0 - success, otherwise - failure 143 **/ 144 int CeedOperatorApply(CeedOperator op, CeedVector in, 145 CeedVector out, CeedRequest *request) { 146 int ierr; 147 Ceed ceed = op->ceed; 148 CeedQFunction qf = op->qf; 149 150 if (op->nfields == 0) return CeedError(ceed, 1, "No operator fields set"); 151 if (op->nfields < qf->numinputfields + qf->numoutputfields) return CeedError( 152 ceed, 1, "Not all operator fields set"); 153 if (op->numelements == 0) return CeedError(ceed, 1, 154 "At least one restriction required"); 155 if (op->numqpoints == 0) return CeedError(ceed, 1, 156 "At least one non-colocated basis required"); 157 ierr = op->Apply(op, in, out, request); CeedChk(ierr); 158 return 0; 159 } 160 161 /** 162 @brief Destroy a CeedOperator 163 164 @param op CeedOperator to destroy 165 166 @return An error code: 0 - success, otherwise - failure 167 **/ 168 int CeedOperatorDestroy(CeedOperator *op) { 169 int ierr; 170 171 if (!*op || --(*op)->refcount > 0) return 0; 172 if ((*op)->Destroy) { 173 ierr = (*op)->Destroy(*op); CeedChk(ierr); 174 } 175 ierr = CeedQFunctionDestroy(&(*op)->qf); CeedChk(ierr); 176 ierr = CeedQFunctionDestroy(&(*op)->dqf); CeedChk(ierr); 177 ierr = CeedQFunctionDestroy(&(*op)->dqfT); CeedChk(ierr); 178 ierr = CeedDestroy(&(*op)->ceed); CeedChk(ierr); 179 ierr = CeedFree(op); CeedChk(ierr); 180 return 0; 181 } 182 183 /// @} 184