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 Create an operator from element restriction, basis, and QFunction 30 31 @param ceed The Ceed library context on which to create the operator 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 Newly created CeedOperator 37 @return Error code, 0 on success 38 */ 39 int CeedOperatorCreate(Ceed ceed, CeedQFunction qf, CeedQFunction dqf, 40 CeedQFunction dqfT, CeedOperator *op) { 41 int ierr; 42 43 if (!ceed->OperatorCreate) return CeedError(ceed, 1, 44 "Backend does not support OperatorCreate"); 45 ierr = CeedCalloc(1,op); CeedChk(ierr); 46 (*op)->ceed = ceed; 47 ceed->refcount++; 48 (*op)->refcount = 1; 49 (*op)->qf = qf; 50 qf->refcount++; 51 (*op)->dqf = dqf; 52 if (dqf) dqf->refcount++; 53 (*op)->dqfT = dqfT; 54 if (dqfT) dqfT->refcount++; 55 ierr = ceed->OperatorCreate(*op); CeedChk(ierr); 56 return 0; 57 } 58 59 /** 60 Provide a field to a CeedOperator for use by its CeedQFunction 61 62 This function is used to specify both active and passive fields to a 63 CeedOperator. For passive fields, a vector @arg v must be provided. Passive 64 fields can inputs or outputs (updated in-place when operator is applied). 65 66 Active fields must be specified using this function, but their data (in a 67 CeedVector) is passed in CeedOperatorApply(). There can be at most one active 68 input and at most one active output. 69 70 @param op the operator on which to provide the field 71 @param fieldname name of the field (to be matched with the name used by CeedQFunction) 72 @param r element restriction 73 @param b basis in which the field resides or CEED_BASIS_COLOCATED if collocated 74 with quadrature points 75 @param v vector to be used by CeedOperator, CEED_VECTOR_ACTIVE if field is 76 active, or CEED_VECTOR_NONE if using CEED_EVAL_WEIGHT in the qfunction 77 */ 78 int CeedOperatorSetField(CeedOperator op, const char *fieldname, 79 CeedElemRestriction r, CeedBasis b, 80 CeedVector v) { 81 int ierr; 82 CeedInt numelements; 83 ierr = CeedElemRestrictionGetNumElements(r, &numelements); CeedChk(ierr); 84 if (op->numelements && op->numelements != numelements) 85 return CeedError(op->ceed, 1, 86 "ElemRestriction with %d elements incompatible with prior %d elements", 87 numelements, op->numelements); 88 op->numelements = numelements; 89 90 if (b != CEED_BASIS_COLOCATED) { 91 CeedInt numqpoints; 92 ierr = CeedBasisGetNumQuadraturePoints(b, &numqpoints); CeedChk(ierr); 93 if (op->numqpoints && op->numqpoints != numqpoints) 94 return CeedError(op->ceed, 1, 95 "Basis with %d quadrature points incompatible with prior %d points", 96 numqpoints, op->numqpoints); 97 op->numqpoints = numqpoints; 98 } 99 struct CeedOperatorField *ofield; 100 for (CeedInt i=0; i<op->qf->numinputfields; i++) { 101 if (!strcmp(fieldname, op->qf->inputfields[i].fieldname)) { 102 ofield = &op->inputfields[i]; 103 goto found; 104 } 105 } 106 for (CeedInt i=0; i<op->qf->numoutputfields; i++) { 107 if (!strcmp(fieldname, op->qf->outputfields[i].fieldname)) { 108 ofield = &op->outputfields[i]; 109 goto found; 110 } 111 } 112 return CeedError(op->ceed, 1, "QFunction has no knowledge of field '%s'", 113 fieldname); 114 found: 115 ofield->Erestrict = r; 116 ofield->basis = b; 117 ofield->vec = v; 118 op->nfields += 1; 119 return 0; 120 } 121 122 /** 123 Apply CeedOperator to a vector 124 125 This computes the action of the operator on the specified (active) input, 126 yielding its (active) output. All inputs and outputs must be specified using 127 CeedOperatorSetField(). 128 129 @param op CeedOperator to apply 130 @param in CeedVector containing input state or NULL if there are no active 131 inputs 132 @param out CeedVector to store result of applying operator (must be 133 distinct from @a in) or NULL if there are no active outputs 134 @param request Address of CeedRequest for non-blocking completion, else 135 CEED_REQUEST_IMMEDIATE 136 @return Error code, 0 on success 137 */ 138 int CeedOperatorApply(CeedOperator op, CeedVector in, 139 CeedVector out, CeedRequest *request) { 140 int ierr; 141 Ceed ceed = op->ceed; 142 CeedQFunction qf = op->qf; 143 144 if (op->nfields == 0) return CeedError(ceed, 1, "No operator fields set"); 145 if (op->nfields < qf->numinputfields + qf->numoutputfields) return CeedError( 146 ceed, 1, "Not all operator fields set"); 147 if (op->numelements == 0) return CeedError(ceed, 1, 148 "At least one restriction required"); 149 if (op->numqpoints == 0) return CeedError(ceed, 1, 150 "At least one non-colocated basis required"); 151 ierr = op->Apply(op, in, out, request); CeedChk(ierr); 152 return 0; 153 } 154 155 /** 156 Destroy a CeedOperator 157 158 @param op CeedOperator to destroy 159 @return Error code, 0 on success 160 */ 161 int CeedOperatorDestroy(CeedOperator *op) { 162 int ierr; 163 164 if (!*op || --(*op)->refcount > 0) return 0; 165 if ((*op)->Destroy) { 166 ierr = (*op)->Destroy(*op); CeedChk(ierr); 167 } 168 ierr = CeedQFunctionDestroy(&(*op)->qf); CeedChk(ierr); 169 ierr = CeedQFunctionDestroy(&(*op)->dqf); CeedChk(ierr); 170 ierr = CeedQFunctionDestroy(&(*op)->dqfT); CeedChk(ierr); 171 ierr = CeedDestroy(&(*op)->ceed); CeedChk(ierr); 172 ierr = CeedFree(op); CeedChk(ierr); 173 return 0; 174 } 175 176 /// @} 177