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/ceed.h> 18 #include <ceed/backend.h> 19 #include <ceed-impl.h> 20 #include <math.h> 21 #include <stdbool.h> 22 #include <stdio.h> 23 #include <string.h> 24 25 /// @file 26 /// Implementation of CeedOperator interfaces 27 28 /// ---------------------------------------------------------------------------- 29 /// CeedOperator Library Internal Functions 30 /// ---------------------------------------------------------------------------- 31 /// @addtogroup CeedOperatorDeveloper 32 /// @{ 33 34 /** 35 @brief Duplicate a CeedOperator with a reference Ceed to fallback for advanced 36 CeedOperator functionality 37 38 @param op CeedOperator to create fallback for 39 40 @return An error code: 0 - success, otherwise - failure 41 42 @ref Developer 43 **/ 44 int CeedOperatorCreateFallback(CeedOperator op) { 45 int ierr; 46 47 // Fallback Ceed 48 const char *resource, *fallbackresource; 49 ierr = CeedGetResource(op->ceed, &resource); CeedChk(ierr); 50 ierr = CeedGetOperatorFallbackResource(op->ceed, &fallbackresource); 51 CeedChk(ierr); 52 if (!strcmp(resource, fallbackresource)) 53 // LCOV_EXCL_START 54 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, 55 "Backend %s cannot create an operator" 56 "fallback to resource %s", resource, fallbackresource); 57 // LCOV_EXCL_STOP 58 59 // Fallback Ceed 60 Ceed ceedref; 61 if (!op->ceed->opfallbackceed) { 62 ierr = CeedInit(fallbackresource, &ceedref); CeedChk(ierr); 63 ceedref->opfallbackparent = op->ceed; 64 ceedref->Error = op->ceed->Error; 65 op->ceed->opfallbackceed = ceedref; 66 } 67 ceedref = op->ceed->opfallbackceed; 68 69 // Clone Op 70 CeedOperator opref; 71 ierr = CeedCalloc(1, &opref); CeedChk(ierr); 72 memcpy(opref, op, sizeof(*opref)); 73 opref->data = NULL; 74 opref->interfacesetup = false; 75 opref->backendsetup = false; 76 opref->ceed = ceedref; 77 ierr = ceedref->OperatorCreate(opref); CeedChk(ierr); 78 op->opfallback = opref; 79 80 // Clone QF 81 CeedQFunction qfref; 82 ierr = CeedCalloc(1, &qfref); CeedChk(ierr); 83 memcpy(qfref, (op->qf), sizeof(*qfref)); 84 qfref->data = NULL; 85 qfref->ceed = ceedref; 86 ierr = ceedref->QFunctionCreate(qfref); CeedChk(ierr); 87 opref->qf = qfref; 88 op->qffallback = qfref; 89 return CEED_ERROR_SUCCESS; 90 } 91 92 /** 93 @brief Check if a CeedOperator Field matches the QFunction Field 94 95 @param[in] ceed Ceed object for error handling 96 @param[in] qfield QFunction Field matching Operator Field 97 @param[in] r Operator Field ElemRestriction 98 @param[in] b Operator Field Basis 99 100 @return An error code: 0 - success, otherwise - failure 101 102 @ref Developer 103 **/ 104 static int CeedOperatorCheckField(Ceed ceed, CeedQFunctionField qfield, 105 CeedElemRestriction r, CeedBasis b) { 106 int ierr; 107 CeedEvalMode emode = qfield->emode; 108 CeedInt dim = 1, ncomp = 1, restr_ncomp = 1, size = qfield->size; 109 // Restriction 110 if (r != CEED_ELEMRESTRICTION_NONE) { 111 if (emode == CEED_EVAL_WEIGHT) { 112 // LCOV_EXCL_START 113 return CeedError(ceed, CEED_ERROR_INCOMPATIBLE, 114 "CEED_ELEMRESTRICTION_NONE should be used " 115 "for a field with eval mode CEED_EVAL_WEIGHT"); 116 // LCOV_EXCL_STOP 117 } 118 ierr = CeedElemRestrictionGetNumComponents(r, &restr_ncomp); 119 } 120 if ((r == CEED_ELEMRESTRICTION_NONE) != (emode == CEED_EVAL_WEIGHT)) { 121 // LCOV_EXCL_START 122 return CeedError(ceed, CEED_ERROR_INCOMPATIBLE, 123 "CEED_ELEMRESTRICTION_NONE and CEED_EVAL_WEIGHT " 124 "must be used together."); 125 // LCOV_EXCL_STOP 126 } 127 // Basis 128 if (b != CEED_BASIS_COLLOCATED) { 129 if (emode == CEED_EVAL_NONE) 130 return CeedError(ceed, CEED_ERROR_INCOMPATIBLE, 131 "Field '%s' configured with CEED_EVAL_NONE must be used with CEED_BASIS_COLLOCATED", 132 qfield->fieldname); 133 ierr = CeedBasisGetDimension(b, &dim); CeedChk(ierr); 134 ierr = CeedBasisGetNumComponents(b, &ncomp); CeedChk(ierr); 135 if (r != CEED_ELEMRESTRICTION_NONE && restr_ncomp != ncomp) { 136 // LCOV_EXCL_START 137 return CeedError(ceed, CEED_ERROR_DIMENSION, 138 "Field '%s' of size %d and EvalMode %s: ElemRestriction has %d components, but Basis has %d components", 139 qfield->fieldname, qfield->size, CeedEvalModes[qfield->emode], restr_ncomp, 140 ncomp); 141 // LCOV_EXCL_STOP 142 } 143 } 144 // Field size 145 switch(emode) { 146 case CEED_EVAL_NONE: 147 if (size != restr_ncomp) 148 // LCOV_EXCL_START 149 return CeedError(ceed, CEED_ERROR_DIMENSION, 150 "Field '%s' of size %d and EvalMode %s: ElemRestriction has %d components", 151 qfield->fieldname, qfield->size, CeedEvalModes[qfield->emode], restr_ncomp); 152 // LCOV_EXCL_STOP 153 break; 154 case CEED_EVAL_INTERP: 155 if (size != ncomp) 156 // LCOV_EXCL_START 157 return CeedError(ceed, CEED_ERROR_DIMENSION, 158 "Field '%s' of size %d and EvalMode %s: ElemRestriction/Basis has %d components", 159 qfield->fieldname, qfield->size, CeedEvalModes[qfield->emode], ncomp); 160 // LCOV_EXCL_STOP 161 break; 162 case CEED_EVAL_GRAD: 163 if (size != ncomp * dim) 164 // LCOV_EXCL_START 165 return CeedError(ceed, CEED_ERROR_DIMENSION, 166 "Field '%s' of size %d and EvalMode %s in %d dimensions: ElemRestriction/Basis has %d components", 167 qfield->fieldname, qfield->size, CeedEvalModes[qfield->emode], dim, ncomp); 168 // LCOV_EXCL_STOP 169 break; 170 case CEED_EVAL_WEIGHT: 171 // No addional checks required 172 break; 173 case CEED_EVAL_DIV: 174 // Not implemented 175 break; 176 case CEED_EVAL_CURL: 177 // Not implemented 178 break; 179 } 180 return CEED_ERROR_SUCCESS; 181 } 182 183 /** 184 @brief Check if a CeedOperator is ready to be used. 185 186 @param[in] op CeedOperator to check 187 188 @return An error code: 0 - success, otherwise - failure 189 190 @ref Developer 191 **/ 192 static int CeedOperatorCheckReady(CeedOperator op) { 193 int ierr; 194 Ceed ceed; 195 ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr); 196 197 if (op->interfacesetup) 198 return CEED_ERROR_SUCCESS; 199 200 CeedQFunction qf = op->qf; 201 if (op->composite) { 202 if (!op->numsub) 203 // LCOV_EXCL_START 204 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "No suboperators set"); 205 // LCOV_EXCL_STOP 206 } else { 207 if (op->nfields == 0) 208 // LCOV_EXCL_START 209 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "No operator fields set"); 210 // LCOV_EXCL_STOP 211 if (op->nfields < qf->numinputfields + qf->numoutputfields) 212 // LCOV_EXCL_START 213 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "Not all operator fields set"); 214 // LCOV_EXCL_STOP 215 if (!op->hasrestriction) 216 // LCOV_EXCL_START 217 return CeedError(ceed, CEED_ERROR_INCOMPLETE, 218 "At least one restriction required"); 219 // LCOV_EXCL_STOP 220 if (op->numqpoints == 0) 221 // LCOV_EXCL_START 222 return CeedError(ceed, CEED_ERROR_INCOMPLETE, 223 "At least one non-collocated basis required"); 224 // LCOV_EXCL_STOP 225 } 226 227 // Flag as immutable and ready 228 op->interfacesetup = true; 229 if (op->qf && op->qf != CEED_QFUNCTION_NONE) 230 // LCOV_EXCL_START 231 op->qf->operatorsset++; 232 // LCOV_EXCL_STOP 233 if (op->dqf && op->dqf != CEED_QFUNCTION_NONE) 234 // LCOV_EXCL_START 235 op->dqf->operatorsset++; 236 // LCOV_EXCL_STOP 237 if (op->dqfT && op->dqfT != CEED_QFUNCTION_NONE) 238 // LCOV_EXCL_START 239 op->dqfT->operatorsset++; 240 // LCOV_EXCL_STOP 241 return CEED_ERROR_SUCCESS; 242 } 243 244 /** 245 @brief View a field of a CeedOperator 246 247 @param[in] field Operator field to view 248 @param[in] qffield QFunction field (carries field name) 249 @param[in] fieldnumber Number of field being viewed 250 @param[in] sub true indicates sub-operator, which increases indentation; false for top-level operator 251 @param[in] in true for an input field; false for output field 252 @param[in] stream Stream to view to, e.g., stdout 253 254 @return An error code: 0 - success, otherwise - failure 255 256 @ref Utility 257 **/ 258 static int CeedOperatorFieldView(CeedOperatorField field, 259 CeedQFunctionField qffield, 260 CeedInt fieldnumber, bool sub, bool in, 261 FILE *stream) { 262 const char *pre = sub ? " " : ""; 263 const char *inout = in ? "Input" : "Output"; 264 265 fprintf(stream, "%s %s Field [%d]:\n" 266 "%s Name: \"%s\"\n", 267 pre, inout, fieldnumber, pre, qffield->fieldname); 268 269 if (field->basis == CEED_BASIS_COLLOCATED) 270 fprintf(stream, "%s Collocated basis\n", pre); 271 272 if (field->vec == CEED_VECTOR_ACTIVE) 273 fprintf(stream, "%s Active vector\n", pre); 274 else if (field->vec == CEED_VECTOR_NONE) 275 fprintf(stream, "%s No vector\n", pre); 276 return CEED_ERROR_SUCCESS; 277 } 278 279 /** 280 @brief View a single CeedOperator 281 282 @param[in] op CeedOperator to view 283 @param[in] sub Boolean flag for sub-operator 284 @param[in] stream Stream to write; typically stdout/stderr or a file 285 286 @return Error code: 0 - success, otherwise - failure 287 288 @ref Utility 289 **/ 290 int CeedOperatorSingleView(CeedOperator op, bool sub, FILE *stream) { 291 int ierr; 292 const char *pre = sub ? " " : ""; 293 294 CeedInt totalfields; 295 ierr = CeedOperatorGetNumArgs(op, &totalfields); CeedChk(ierr); 296 297 fprintf(stream, "%s %d Field%s\n", pre, totalfields, totalfields>1 ? "s" : ""); 298 299 fprintf(stream, "%s %d Input Field%s:\n", pre, op->qf->numinputfields, 300 op->qf->numinputfields>1 ? "s" : ""); 301 for (CeedInt i=0; i<op->qf->numinputfields; i++) { 302 ierr = CeedOperatorFieldView(op->inputfields[i], op->qf->inputfields[i], 303 i, sub, 1, stream); CeedChk(ierr); 304 } 305 306 fprintf(stream, "%s %d Output Field%s:\n", pre, op->qf->numoutputfields, 307 op->qf->numoutputfields>1 ? "s" : ""); 308 for (CeedInt i=0; i<op->qf->numoutputfields; i++) { 309 ierr = CeedOperatorFieldView(op->outputfields[i], op->qf->outputfields[i], 310 i, sub, 0, stream); CeedChk(ierr); 311 } 312 return CEED_ERROR_SUCCESS; 313 } 314 315 /** 316 @brief Find the active vector ElemRestriction for a CeedOperator 317 318 @param[in] op CeedOperator to find active basis for 319 @param[out] activerstr ElemRestriction for active input vector 320 321 @return An error code: 0 - success, otherwise - failure 322 323 @ref Utility 324 **/ 325 static int CeedOperatorGetActiveElemRestriction(CeedOperator op, 326 CeedElemRestriction *activerstr) { 327 *activerstr = NULL; 328 for (int i = 0; i < op->qf->numinputfields; i++) 329 if (op->inputfields[i]->vec == CEED_VECTOR_ACTIVE) { 330 *activerstr = op->inputfields[i]->Erestrict; 331 break; 332 } 333 334 if (!*activerstr) { 335 // LCOV_EXCL_START 336 int ierr; 337 Ceed ceed; 338 ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr); 339 return CeedError(ceed, CEED_ERROR_INCOMPLETE, 340 "No active ElemRestriction found!"); 341 // LCOV_EXCL_STOP 342 } 343 return CEED_ERROR_SUCCESS; 344 } 345 346 /** 347 @brief Find the active vector basis for a CeedOperator 348 349 @param[in] op CeedOperator to find active basis for 350 @param[out] activeBasis Basis for active input vector 351 352 @return An error code: 0 - success, otherwise - failure 353 354 @ ref Developer 355 **/ 356 static int CeedOperatorGetActiveBasis(CeedOperator op, 357 CeedBasis *activeBasis) { 358 *activeBasis = NULL; 359 for (int i = 0; i < op->qf->numinputfields; i++) 360 if (op->inputfields[i]->vec == CEED_VECTOR_ACTIVE) { 361 *activeBasis = op->inputfields[i]->basis; 362 break; 363 } 364 365 if (!*activeBasis) { 366 // LCOV_EXCL_START 367 int ierr; 368 Ceed ceed; 369 ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr); 370 return CeedError(ceed, CEED_ERROR_MINOR, 371 "No active basis found for automatic multigrid setup"); 372 // LCOV_EXCL_STOP 373 } 374 return CEED_ERROR_SUCCESS; 375 } 376 377 /** 378 @brief Common code for creating a multigrid coarse operator and level 379 transfer operators for a CeedOperator 380 381 @param[in] opFine Fine grid operator 382 @param[in] PMultFine L-vector multiplicity in parallel gather/scatter 383 @param[in] rstrCoarse Coarse grid restriction 384 @param[in] basisCoarse Coarse grid active vector basis 385 @param[in] basisCtoF Basis for coarse to fine interpolation 386 @param[out] opCoarse Coarse grid operator 387 @param[out] opProlong Coarse to fine operator 388 @param[out] opRestrict Fine to coarse operator 389 390 @return An error code: 0 - success, otherwise - failure 391 392 @ref Developer 393 **/ 394 static int CeedOperatorMultigridLevel_Core(CeedOperator opFine, 395 CeedVector PMultFine, CeedElemRestriction rstrCoarse, CeedBasis basisCoarse, 396 CeedBasis basisCtoF, CeedOperator *opCoarse, CeedOperator *opProlong, 397 CeedOperator *opRestrict) { 398 int ierr; 399 Ceed ceed; 400 ierr = CeedOperatorGetCeed(opFine, &ceed); CeedChk(ierr); 401 402 // Check for composite operator 403 bool isComposite; 404 ierr = CeedOperatorIsComposite(opFine, &isComposite); CeedChk(ierr); 405 if (isComposite) 406 // LCOV_EXCL_START 407 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 408 "Automatic multigrid setup for composite operators not supported"); 409 // LCOV_EXCL_STOP 410 411 // Coarse Grid 412 ierr = CeedOperatorCreate(ceed, opFine->qf, opFine->dqf, opFine->dqfT, 413 opCoarse); CeedChk(ierr); 414 CeedElemRestriction rstrFine = NULL; 415 // -- Clone input fields 416 for (int i = 0; i < opFine->qf->numinputfields; i++) { 417 if (opFine->inputfields[i]->vec == CEED_VECTOR_ACTIVE) { 418 rstrFine = opFine->inputfields[i]->Erestrict; 419 ierr = CeedOperatorSetField(*opCoarse, opFine->inputfields[i]->fieldname, 420 rstrCoarse, basisCoarse, CEED_VECTOR_ACTIVE); 421 CeedChk(ierr); 422 } else { 423 ierr = CeedOperatorSetField(*opCoarse, opFine->inputfields[i]->fieldname, 424 opFine->inputfields[i]->Erestrict, 425 opFine->inputfields[i]->basis, 426 opFine->inputfields[i]->vec); CeedChk(ierr); 427 } 428 } 429 // -- Clone output fields 430 for (int i = 0; i < opFine->qf->numoutputfields; i++) { 431 if (opFine->outputfields[i]->vec == CEED_VECTOR_ACTIVE) { 432 ierr = CeedOperatorSetField(*opCoarse, opFine->outputfields[i]->fieldname, 433 rstrCoarse, basisCoarse, CEED_VECTOR_ACTIVE); 434 CeedChk(ierr); 435 } else { 436 ierr = CeedOperatorSetField(*opCoarse, opFine->outputfields[i]->fieldname, 437 opFine->outputfields[i]->Erestrict, 438 opFine->outputfields[i]->basis, 439 opFine->outputfields[i]->vec); CeedChk(ierr); 440 } 441 } 442 443 // Multiplicity vector 444 CeedVector multVec, multE; 445 ierr = CeedElemRestrictionCreateVector(rstrFine, &multVec, &multE); 446 CeedChk(ierr); 447 ierr = CeedVectorSetValue(multE, 0.0); CeedChk(ierr); 448 ierr = CeedElemRestrictionApply(rstrFine, CEED_NOTRANSPOSE, PMultFine, multE, 449 CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 450 ierr = CeedVectorSetValue(multVec, 0.0); CeedChk(ierr); 451 ierr = CeedElemRestrictionApply(rstrFine, CEED_TRANSPOSE, multE, multVec, 452 CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 453 ierr = CeedVectorDestroy(&multE); CeedChk(ierr); 454 ierr = CeedVectorReciprocal(multVec); CeedChk(ierr); 455 456 // Restriction 457 CeedInt ncomp; 458 ierr = CeedBasisGetNumComponents(basisCoarse, &ncomp); CeedChk(ierr); 459 CeedQFunction qfRestrict; 460 ierr = CeedQFunctionCreateInteriorByName(ceed, "Scale", &qfRestrict); 461 CeedChk(ierr); 462 CeedInt *ncompRData; 463 ierr = CeedCalloc(1, &ncompRData); CeedChk(ierr); 464 ncompRData[0] = ncomp; 465 CeedQFunctionContext ctxR; 466 ierr = CeedQFunctionContextCreate(ceed, &ctxR); CeedChk(ierr); 467 ierr = CeedQFunctionContextSetData(ctxR, CEED_MEM_HOST, CEED_OWN_POINTER, 468 sizeof(*ncompRData), ncompRData); 469 CeedChk(ierr); 470 ierr = CeedQFunctionSetContext(qfRestrict, ctxR); CeedChk(ierr); 471 ierr = CeedQFunctionContextDestroy(&ctxR); CeedChk(ierr); 472 ierr = CeedQFunctionAddInput(qfRestrict, "input", ncomp, CEED_EVAL_NONE); 473 CeedChk(ierr); 474 ierr = CeedQFunctionAddInput(qfRestrict, "scale", ncomp, CEED_EVAL_NONE); 475 CeedChk(ierr); 476 ierr = CeedQFunctionAddOutput(qfRestrict, "output", ncomp, CEED_EVAL_INTERP); 477 CeedChk(ierr); 478 479 ierr = CeedOperatorCreate(ceed, qfRestrict, CEED_QFUNCTION_NONE, 480 CEED_QFUNCTION_NONE, opRestrict); 481 CeedChk(ierr); 482 ierr = CeedOperatorSetField(*opRestrict, "input", rstrFine, 483 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 484 CeedChk(ierr); 485 ierr = CeedOperatorSetField(*opRestrict, "scale", rstrFine, 486 CEED_BASIS_COLLOCATED, multVec); 487 CeedChk(ierr); 488 ierr = CeedOperatorSetField(*opRestrict, "output", rstrCoarse, basisCtoF, 489 CEED_VECTOR_ACTIVE); CeedChk(ierr); 490 491 // Prolongation 492 CeedQFunction qfProlong; 493 ierr = CeedQFunctionCreateInteriorByName(ceed, "Scale", &qfProlong); 494 CeedChk(ierr); 495 CeedInt *ncompPData; 496 ierr = CeedCalloc(1, &ncompPData); CeedChk(ierr); 497 ncompPData[0] = ncomp; 498 CeedQFunctionContext ctxP; 499 ierr = CeedQFunctionContextCreate(ceed, &ctxP); CeedChk(ierr); 500 ierr = CeedQFunctionContextSetData(ctxP, CEED_MEM_HOST, CEED_OWN_POINTER, 501 sizeof(*ncompPData), ncompPData); 502 CeedChk(ierr); 503 ierr = CeedQFunctionSetContext(qfProlong, ctxP); CeedChk(ierr); 504 ierr = CeedQFunctionContextDestroy(&ctxP); CeedChk(ierr); 505 ierr = CeedQFunctionAddInput(qfProlong, "input", ncomp, CEED_EVAL_INTERP); 506 CeedChk(ierr); 507 ierr = CeedQFunctionAddInput(qfProlong, "scale", ncomp, CEED_EVAL_NONE); 508 CeedChk(ierr); 509 ierr = CeedQFunctionAddOutput(qfProlong, "output", ncomp, CEED_EVAL_NONE); 510 CeedChk(ierr); 511 512 ierr = CeedOperatorCreate(ceed, qfProlong, CEED_QFUNCTION_NONE, 513 CEED_QFUNCTION_NONE, opProlong); 514 CeedChk(ierr); 515 ierr = CeedOperatorSetField(*opProlong, "input", rstrCoarse, basisCtoF, 516 CEED_VECTOR_ACTIVE); CeedChk(ierr); 517 ierr = CeedOperatorSetField(*opProlong, "scale", rstrFine, 518 CEED_BASIS_COLLOCATED, multVec); 519 CeedChk(ierr); 520 ierr = CeedOperatorSetField(*opProlong, "output", rstrFine, 521 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 522 CeedChk(ierr); 523 524 // Cleanup 525 ierr = CeedVectorDestroy(&multVec); CeedChk(ierr); 526 ierr = CeedBasisDestroy(&basisCtoF); CeedChk(ierr); 527 ierr = CeedQFunctionDestroy(&qfRestrict); CeedChk(ierr); 528 ierr = CeedQFunctionDestroy(&qfProlong); CeedChk(ierr); 529 return CEED_ERROR_SUCCESS; 530 } 531 532 /// @} 533 534 /// ---------------------------------------------------------------------------- 535 /// CeedOperator Backend API 536 /// ---------------------------------------------------------------------------- 537 /// @addtogroup CeedOperatorBackend 538 /// @{ 539 540 /** 541 @brief Get the Ceed associated with a CeedOperator 542 543 @param op CeedOperator 544 @param[out] ceed Variable to store Ceed 545 546 @return An error code: 0 - success, otherwise - failure 547 548 @ref Backend 549 **/ 550 551 int CeedOperatorGetCeed(CeedOperator op, Ceed *ceed) { 552 *ceed = op->ceed; 553 return CEED_ERROR_SUCCESS; 554 } 555 556 /** 557 @brief Get the number of elements associated with a CeedOperator 558 559 @param op CeedOperator 560 @param[out] numelem Variable to store number of elements 561 562 @return An error code: 0 - success, otherwise - failure 563 564 @ref Backend 565 **/ 566 567 int CeedOperatorGetNumElements(CeedOperator op, CeedInt *numelem) { 568 if (op->composite) 569 // LCOV_EXCL_START 570 return CeedError(op->ceed, CEED_ERROR_MINOR, 571 "Not defined for composite operator"); 572 // LCOV_EXCL_STOP 573 574 *numelem = op->numelements; 575 return CEED_ERROR_SUCCESS; 576 } 577 578 /** 579 @brief Get the number of quadrature points associated with a CeedOperator 580 581 @param op CeedOperator 582 @param[out] numqpts Variable to store vector number of quadrature points 583 584 @return An error code: 0 - success, otherwise - failure 585 586 @ref Backend 587 **/ 588 589 int CeedOperatorGetNumQuadraturePoints(CeedOperator op, CeedInt *numqpts) { 590 if (op->composite) 591 // LCOV_EXCL_START 592 return CeedError(op->ceed, CEED_ERROR_MINOR, 593 "Not defined for composite operator"); 594 // LCOV_EXCL_STOP 595 596 *numqpts = op->numqpoints; 597 return CEED_ERROR_SUCCESS; 598 } 599 600 /** 601 @brief Get the number of arguments associated with a CeedOperator 602 603 @param op CeedOperator 604 @param[out] numargs Variable to store vector number of arguments 605 606 @return An error code: 0 - success, otherwise - failure 607 608 @ref Backend 609 **/ 610 611 int CeedOperatorGetNumArgs(CeedOperator op, CeedInt *numargs) { 612 if (op->composite) 613 // LCOV_EXCL_START 614 return CeedError(op->ceed, CEED_ERROR_MINOR, 615 "Not defined for composite operators"); 616 // LCOV_EXCL_STOP 617 618 *numargs = op->nfields; 619 return CEED_ERROR_SUCCESS; 620 } 621 622 /** 623 @brief Get the setup status of a CeedOperator 624 625 @param op CeedOperator 626 @param[out] issetupdone Variable to store setup status 627 628 @return An error code: 0 - success, otherwise - failure 629 630 @ref Backend 631 **/ 632 633 int CeedOperatorIsSetupDone(CeedOperator op, bool *issetupdone) { 634 *issetupdone = op->backendsetup; 635 return CEED_ERROR_SUCCESS; 636 } 637 638 /** 639 @brief Get the QFunction associated with a CeedOperator 640 641 @param op CeedOperator 642 @param[out] qf Variable to store QFunction 643 644 @return An error code: 0 - success, otherwise - failure 645 646 @ref Backend 647 **/ 648 649 int CeedOperatorGetQFunction(CeedOperator op, CeedQFunction *qf) { 650 if (op->composite) 651 // LCOV_EXCL_START 652 return CeedError(op->ceed, CEED_ERROR_MINOR, 653 "Not defined for composite operator"); 654 // LCOV_EXCL_STOP 655 656 *qf = op->qf; 657 return CEED_ERROR_SUCCESS; 658 } 659 660 /** 661 @brief Get a boolean value indicating if the CeedOperator is composite 662 663 @param op CeedOperator 664 @param[out] iscomposite Variable to store composite status 665 666 @return An error code: 0 - success, otherwise - failure 667 668 @ref Backend 669 **/ 670 671 int CeedOperatorIsComposite(CeedOperator op, bool *iscomposite) { 672 *iscomposite = op->composite; 673 return CEED_ERROR_SUCCESS; 674 } 675 676 /** 677 @brief Get the number of suboperators associated with a CeedOperator 678 679 @param op CeedOperator 680 @param[out] numsub Variable to store number of suboperators 681 682 @return An error code: 0 - success, otherwise - failure 683 684 @ref Backend 685 **/ 686 687 int CeedOperatorGetNumSub(CeedOperator op, CeedInt *numsub) { 688 if (!op->composite) 689 // LCOV_EXCL_START 690 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not a composite operator"); 691 // LCOV_EXCL_STOP 692 693 *numsub = op->numsub; 694 return CEED_ERROR_SUCCESS; 695 } 696 697 /** 698 @brief Get the list of suboperators associated with a CeedOperator 699 700 @param op CeedOperator 701 @param[out] suboperators Variable to store list of suboperators 702 703 @return An error code: 0 - success, otherwise - failure 704 705 @ref Backend 706 **/ 707 708 int CeedOperatorGetSubList(CeedOperator op, CeedOperator **suboperators) { 709 if (!op->composite) 710 // LCOV_EXCL_START 711 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not a composite operator"); 712 // LCOV_EXCL_STOP 713 714 *suboperators = op->suboperators; 715 return CEED_ERROR_SUCCESS; 716 } 717 718 /** 719 @brief Get the backend data of a CeedOperator 720 721 @param op CeedOperator 722 @param[out] data Variable to store data 723 724 @return An error code: 0 - success, otherwise - failure 725 726 @ref Backend 727 **/ 728 729 int CeedOperatorGetData(CeedOperator op, void *data) { 730 *(void **)data = op->data; 731 return CEED_ERROR_SUCCESS; 732 } 733 734 /** 735 @brief Set the backend data of a CeedOperator 736 737 @param[out] op CeedOperator 738 @param data Data to set 739 740 @return An error code: 0 - success, otherwise - failure 741 742 @ref Backend 743 **/ 744 745 int CeedOperatorSetData(CeedOperator op, void *data) { 746 op->data = data; 747 return CEED_ERROR_SUCCESS; 748 } 749 750 /** 751 @brief Set the setup flag of a CeedOperator to True 752 753 @param op CeedOperator 754 755 @return An error code: 0 - success, otherwise - failure 756 757 @ref Backend 758 **/ 759 760 int CeedOperatorSetSetupDone(CeedOperator op) { 761 op->backendsetup = true; 762 return CEED_ERROR_SUCCESS; 763 } 764 765 /** 766 @brief Get the CeedOperatorFields of a CeedOperator 767 768 @param op CeedOperator 769 @param[out] inputfields Variable to store inputfields 770 @param[out] outputfields Variable to store outputfields 771 772 @return An error code: 0 - success, otherwise - failure 773 774 @ref Backend 775 **/ 776 777 int CeedOperatorGetFields(CeedOperator op, CeedOperatorField **inputfields, 778 CeedOperatorField **outputfields) { 779 if (op->composite) 780 // LCOV_EXCL_START 781 return CeedError(op->ceed, CEED_ERROR_MINOR, 782 "Not defined for composite operator"); 783 // LCOV_EXCL_STOP 784 785 if (inputfields) *inputfields = op->inputfields; 786 if (outputfields) *outputfields = op->outputfields; 787 return CEED_ERROR_SUCCESS; 788 } 789 790 /** 791 @brief Get the CeedElemRestriction of a CeedOperatorField 792 793 @param opfield CeedOperatorField 794 @param[out] rstr Variable to store CeedElemRestriction 795 796 @return An error code: 0 - success, otherwise - failure 797 798 @ref Backend 799 **/ 800 801 int CeedOperatorFieldGetElemRestriction(CeedOperatorField opfield, 802 CeedElemRestriction *rstr) { 803 *rstr = opfield->Erestrict; 804 return CEED_ERROR_SUCCESS; 805 } 806 807 /** 808 @brief Get the CeedBasis of a CeedOperatorField 809 810 @param opfield CeedOperatorField 811 @param[out] basis Variable to store CeedBasis 812 813 @return An error code: 0 - success, otherwise - failure 814 815 @ref Backend 816 **/ 817 818 int CeedOperatorFieldGetBasis(CeedOperatorField opfield, CeedBasis *basis) { 819 *basis = opfield->basis; 820 return CEED_ERROR_SUCCESS; 821 } 822 823 /** 824 @brief Get the CeedVector of a CeedOperatorField 825 826 @param opfield CeedOperatorField 827 @param[out] vec Variable to store CeedVector 828 829 @return An error code: 0 - success, otherwise - failure 830 831 @ref Backend 832 **/ 833 834 int CeedOperatorFieldGetVector(CeedOperatorField opfield, CeedVector *vec) { 835 *vec = opfield->vec; 836 return CEED_ERROR_SUCCESS; 837 } 838 839 /// @} 840 841 /// ---------------------------------------------------------------------------- 842 /// CeedOperator Public API 843 /// ---------------------------------------------------------------------------- 844 /// @addtogroup CeedOperatorUser 845 /// @{ 846 847 /** 848 @brief Create a CeedOperator and associate a CeedQFunction. A CeedBasis and 849 CeedElemRestriction can be associated with CeedQFunction fields with 850 \ref CeedOperatorSetField. 851 852 @param ceed A Ceed object where the CeedOperator will be created 853 @param qf QFunction defining the action of the operator at quadrature points 854 @param dqf QFunction defining the action of the Jacobian of @a qf (or 855 @ref CEED_QFUNCTION_NONE) 856 @param dqfT QFunction defining the action of the transpose of the Jacobian 857 of @a qf (or @ref CEED_QFUNCTION_NONE) 858 @param[out] op Address of the variable where the newly created 859 CeedOperator will be stored 860 861 @return An error code: 0 - success, otherwise - failure 862 863 @ref User 864 */ 865 int CeedOperatorCreate(Ceed ceed, CeedQFunction qf, CeedQFunction dqf, 866 CeedQFunction dqfT, CeedOperator *op) { 867 int ierr; 868 869 if (!ceed->OperatorCreate) { 870 Ceed delegate; 871 ierr = CeedGetObjectDelegate(ceed, &delegate, "Operator"); CeedChk(ierr); 872 873 if (!delegate) 874 // LCOV_EXCL_START 875 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 876 "Backend does not support OperatorCreate"); 877 // LCOV_EXCL_STOP 878 879 ierr = CeedOperatorCreate(delegate, qf, dqf, dqfT, op); CeedChk(ierr); 880 return CEED_ERROR_SUCCESS; 881 } 882 883 if (!qf || qf == CEED_QFUNCTION_NONE) 884 // LCOV_EXCL_START 885 return CeedError(ceed, CEED_ERROR_MINOR, 886 "Operator must have a valid QFunction."); 887 // LCOV_EXCL_STOP 888 ierr = CeedCalloc(1, op); CeedChk(ierr); 889 (*op)->ceed = ceed; 890 ceed->refcount++; 891 (*op)->refcount = 1; 892 (*op)->qf = qf; 893 qf->refcount++; 894 if (dqf && dqf != CEED_QFUNCTION_NONE) { 895 (*op)->dqf = dqf; 896 dqf->refcount++; 897 } 898 if (dqfT && dqfT != CEED_QFUNCTION_NONE) { 899 (*op)->dqfT = dqfT; 900 dqfT->refcount++; 901 } 902 ierr = CeedCalloc(16, &(*op)->inputfields); CeedChk(ierr); 903 ierr = CeedCalloc(16, &(*op)->outputfields); CeedChk(ierr); 904 ierr = ceed->OperatorCreate(*op); CeedChk(ierr); 905 return CEED_ERROR_SUCCESS; 906 } 907 908 /** 909 @brief Create an operator that composes the action of several operators 910 911 @param ceed A Ceed object where the CeedOperator will be created 912 @param[out] op Address of the variable where the newly created 913 Composite CeedOperator will be stored 914 915 @return An error code: 0 - success, otherwise - failure 916 917 @ref User 918 */ 919 int CeedCompositeOperatorCreate(Ceed ceed, CeedOperator *op) { 920 int ierr; 921 922 if (!ceed->CompositeOperatorCreate) { 923 Ceed delegate; 924 ierr = CeedGetObjectDelegate(ceed, &delegate, "Operator"); CeedChk(ierr); 925 926 if (delegate) { 927 ierr = CeedCompositeOperatorCreate(delegate, op); CeedChk(ierr); 928 return CEED_ERROR_SUCCESS; 929 } 930 } 931 932 ierr = CeedCalloc(1, op); CeedChk(ierr); 933 (*op)->ceed = ceed; 934 ceed->refcount++; 935 (*op)->composite = true; 936 ierr = CeedCalloc(16, &(*op)->suboperators); CeedChk(ierr); 937 938 if (ceed->CompositeOperatorCreate) { 939 ierr = ceed->CompositeOperatorCreate(*op); CeedChk(ierr); 940 } 941 return CEED_ERROR_SUCCESS; 942 } 943 944 /** 945 @brief Provide a field to a CeedOperator for use by its CeedQFunction 946 947 This function is used to specify both active and passive fields to a 948 CeedOperator. For passive fields, a vector @arg v must be provided. Passive 949 fields can inputs or outputs (updated in-place when operator is applied). 950 951 Active fields must be specified using this function, but their data (in a 952 CeedVector) is passed in CeedOperatorApply(). There can be at most one active 953 input and at most one active output. 954 955 @param op CeedOperator on which to provide the field 956 @param fieldname Name of the field (to be matched with the name used by 957 CeedQFunction) 958 @param r CeedElemRestriction 959 @param b CeedBasis in which the field resides or @ref CEED_BASIS_COLLOCATED 960 if collocated with quadrature points 961 @param v CeedVector to be used by CeedOperator or @ref CEED_VECTOR_ACTIVE 962 if field is active or @ref CEED_VECTOR_NONE if using 963 @ref CEED_EVAL_WEIGHT in the QFunction 964 965 @return An error code: 0 - success, otherwise - failure 966 967 @ref User 968 **/ 969 int CeedOperatorSetField(CeedOperator op, const char *fieldname, 970 CeedElemRestriction r, CeedBasis b, CeedVector v) { 971 int ierr; 972 if (op->composite) 973 // LCOV_EXCL_START 974 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, 975 "Cannot add field to composite operator."); 976 // LCOV_EXCL_STOP 977 if (!r) 978 // LCOV_EXCL_START 979 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, 980 "ElemRestriction r for field \"%s\" must be non-NULL.", 981 fieldname); 982 // LCOV_EXCL_STOP 983 if (!b) 984 // LCOV_EXCL_START 985 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, 986 "Basis b for field \"%s\" must be non-NULL.", 987 fieldname); 988 // LCOV_EXCL_STOP 989 if (!v) 990 // LCOV_EXCL_START 991 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, 992 "Vector v for field \"%s\" must be non-NULL.", 993 fieldname); 994 // LCOV_EXCL_STOP 995 996 CeedInt numelements; 997 ierr = CeedElemRestrictionGetNumElements(r, &numelements); CeedChk(ierr); 998 if (r != CEED_ELEMRESTRICTION_NONE && op->hasrestriction && 999 op->numelements != numelements) 1000 // LCOV_EXCL_START 1001 return CeedError(op->ceed, CEED_ERROR_DIMENSION, 1002 "ElemRestriction with %d elements incompatible with prior " 1003 "%d elements", numelements, op->numelements); 1004 // LCOV_EXCL_STOP 1005 1006 CeedInt numqpoints; 1007 if (b != CEED_BASIS_COLLOCATED) { 1008 ierr = CeedBasisGetNumQuadraturePoints(b, &numqpoints); CeedChk(ierr); 1009 if (op->numqpoints && op->numqpoints != numqpoints) 1010 // LCOV_EXCL_START 1011 return CeedError(op->ceed, CEED_ERROR_DIMENSION, 1012 "Basis with %d quadrature points " 1013 "incompatible with prior %d points", numqpoints, 1014 op->numqpoints); 1015 // LCOV_EXCL_STOP 1016 } 1017 CeedQFunctionField qfield; 1018 CeedOperatorField *ofield; 1019 for (CeedInt i=0; i<op->qf->numinputfields; i++) { 1020 if (!strcmp(fieldname, (*op->qf->inputfields[i]).fieldname)) { 1021 qfield = op->qf->inputfields[i]; 1022 ofield = &op->inputfields[i]; 1023 goto found; 1024 } 1025 } 1026 for (CeedInt i=0; i<op->qf->numoutputfields; i++) { 1027 if (!strcmp(fieldname, (*op->qf->outputfields[i]).fieldname)) { 1028 qfield = op->qf->outputfields[i]; 1029 ofield = &op->outputfields[i]; 1030 goto found; 1031 } 1032 } 1033 // LCOV_EXCL_START 1034 return CeedError(op->ceed, CEED_ERROR_INCOMPLETE, 1035 "QFunction has no knowledge of field '%s'", 1036 fieldname); 1037 // LCOV_EXCL_STOP 1038 found: 1039 ierr = CeedOperatorCheckField(op->ceed, qfield, r, b); CeedChk(ierr); 1040 ierr = CeedCalloc(1, ofield); CeedChk(ierr); 1041 1042 (*ofield)->vec = v; 1043 if (v != CEED_VECTOR_ACTIVE && v != CEED_VECTOR_NONE) { 1044 v->refcount += 1; 1045 } 1046 1047 (*ofield)->Erestrict = r; 1048 r->refcount += 1; 1049 if (r != CEED_ELEMRESTRICTION_NONE) { 1050 op->numelements = numelements; 1051 op->hasrestriction = true; // Restriction set, but numelements may be 0 1052 } 1053 1054 (*ofield)->basis = b; 1055 if (b != CEED_BASIS_COLLOCATED) { 1056 op->numqpoints = numqpoints; 1057 b->refcount += 1; 1058 } 1059 1060 op->nfields += 1; 1061 size_t len = strlen(fieldname); 1062 char *tmp; 1063 ierr = CeedCalloc(len+1, &tmp); CeedChk(ierr); 1064 memcpy(tmp, fieldname, len+1); 1065 (*ofield)->fieldname = tmp; 1066 return CEED_ERROR_SUCCESS; 1067 } 1068 1069 /** 1070 @brief Add a sub-operator to a composite CeedOperator 1071 1072 @param[out] compositeop Composite CeedOperator 1073 @param subop Sub-operator CeedOperator 1074 1075 @return An error code: 0 - success, otherwise - failure 1076 1077 @ref User 1078 */ 1079 int CeedCompositeOperatorAddSub(CeedOperator compositeop, CeedOperator subop) { 1080 if (!compositeop->composite) 1081 // LCOV_EXCL_START 1082 return CeedError(compositeop->ceed, CEED_ERROR_MINOR, 1083 "CeedOperator is not a composite operator"); 1084 // LCOV_EXCL_STOP 1085 1086 if (compositeop->numsub == CEED_COMPOSITE_MAX) 1087 // LCOV_EXCL_START 1088 return CeedError(compositeop->ceed, CEED_ERROR_UNSUPPORTED, 1089 "Cannot add additional suboperators"); 1090 // LCOV_EXCL_STOP 1091 1092 compositeop->suboperators[compositeop->numsub] = subop; 1093 subop->refcount++; 1094 compositeop->numsub++; 1095 return CEED_ERROR_SUCCESS; 1096 } 1097 1098 /** 1099 @brief Assemble a linear CeedQFunction associated with a CeedOperator 1100 1101 This returns a CeedVector containing a matrix at each quadrature point 1102 providing the action of the CeedQFunction associated with the CeedOperator. 1103 The vector 'assembled' is of shape 1104 [num_elements, num_input_fields, num_output_fields, num_quad_points] 1105 and contains column-major matrices representing the action of the 1106 CeedQFunction for a corresponding quadrature point on an element. Inputs and 1107 outputs are in the order provided by the user when adding CeedOperator fields. 1108 For example, a CeedQFunction with inputs 'u' and 'gradu' and outputs 'gradv' and 1109 'v', provided in that order, would result in an assembled QFunction that 1110 consists of (1 + dim) x (dim + 1) matrices at each quadrature point acting 1111 on the input [u, du_0, du_1] and producing the output [dv_0, dv_1, v]. 1112 1113 @param op CeedOperator to assemble CeedQFunction 1114 @param[out] assembled CeedVector to store assembled CeedQFunction at 1115 quadrature points 1116 @param[out] rstr CeedElemRestriction for CeedVector containing assembled 1117 CeedQFunction 1118 @param request Address of CeedRequest for non-blocking completion, else 1119 @ref CEED_REQUEST_IMMEDIATE 1120 1121 @return An error code: 0 - success, otherwise - failure 1122 1123 @ref User 1124 **/ 1125 int CeedOperatorLinearAssembleQFunction(CeedOperator op, CeedVector *assembled, 1126 CeedElemRestriction *rstr, 1127 CeedRequest *request) { 1128 int ierr; 1129 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1130 1131 // Backend version 1132 if (op->LinearAssembleQFunction) { 1133 return op->LinearAssembleQFunction(op, assembled, rstr, request); 1134 } else { 1135 // Fallback to reference Ceed 1136 if (!op->opfallback) { 1137 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1138 } 1139 // Assemble 1140 return CeedOperatorLinearAssembleQFunction(op->opfallback, assembled, 1141 rstr, request); 1142 } 1143 } 1144 1145 /** 1146 @brief Assemble the diagonal of a square linear CeedOperator 1147 1148 This overwrites a CeedVector with the diagonal of a linear CeedOperator. 1149 1150 Note: Currently only non-composite CeedOperators with a single field and 1151 composite CeedOperators with single field sub-operators are supported. 1152 1153 @param op CeedOperator to assemble CeedQFunction 1154 @param[out] assembled CeedVector to store assembled CeedOperator diagonal 1155 @param request Address of CeedRequest for non-blocking completion, else 1156 @ref CEED_REQUEST_IMMEDIATE 1157 1158 @return An error code: 0 - success, otherwise - failure 1159 1160 @ref User 1161 **/ 1162 int CeedOperatorLinearAssembleDiagonal(CeedOperator op, CeedVector assembled, 1163 CeedRequest *request) { 1164 int ierr; 1165 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1166 1167 // Use backend version, if available 1168 if (op->LinearAssembleDiagonal) { 1169 return op->LinearAssembleDiagonal(op, assembled, request); 1170 } else if (op->LinearAssembleAddDiagonal) { 1171 ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 1172 return CeedOperatorLinearAssembleAddDiagonal(op, assembled, request); 1173 } else { 1174 // Fallback to reference Ceed 1175 if (!op->opfallback) { 1176 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1177 } 1178 // Assemble 1179 return CeedOperatorLinearAssembleDiagonal(op->opfallback, assembled, 1180 request); 1181 } 1182 } 1183 1184 /** 1185 @brief Assemble the diagonal of a square linear CeedOperator 1186 1187 This sums into a CeedVector the diagonal of a linear CeedOperator. 1188 1189 Note: Currently only non-composite CeedOperators with a single field and 1190 composite CeedOperators with single field sub-operators are supported. 1191 1192 @param op CeedOperator to assemble CeedQFunction 1193 @param[out] assembled CeedVector to store assembled CeedOperator diagonal 1194 @param request Address of CeedRequest for non-blocking completion, else 1195 @ref CEED_REQUEST_IMMEDIATE 1196 1197 @return An error code: 0 - success, otherwise - failure 1198 1199 @ref User 1200 **/ 1201 int CeedOperatorLinearAssembleAddDiagonal(CeedOperator op, CeedVector assembled, 1202 CeedRequest *request) { 1203 int ierr; 1204 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1205 1206 // Use backend version, if available 1207 if (op->LinearAssembleAddDiagonal) { 1208 return op->LinearAssembleAddDiagonal(op, assembled, request); 1209 } else { 1210 // Fallback to reference Ceed 1211 if (!op->opfallback) { 1212 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1213 } 1214 // Assemble 1215 return CeedOperatorLinearAssembleAddDiagonal(op->opfallback, assembled, 1216 request); 1217 } 1218 } 1219 1220 /** 1221 @brief Assemble the point block diagonal of a square linear CeedOperator 1222 1223 This overwrites a CeedVector with the point block diagonal of a linear 1224 CeedOperator. 1225 1226 Note: Currently only non-composite CeedOperators with a single field and 1227 composite CeedOperators with single field sub-operators are supported. 1228 1229 @param op CeedOperator to assemble CeedQFunction 1230 @param[out] assembled CeedVector to store assembled CeedOperator point block 1231 diagonal, provided in row-major form with an 1232 @a ncomp * @a ncomp block at each node. The dimensions 1233 of this vector are derived from the active vector 1234 for the CeedOperator. The array has shape 1235 [nodes, component out, component in]. 1236 @param request Address of CeedRequest for non-blocking completion, else 1237 CEED_REQUEST_IMMEDIATE 1238 1239 @return An error code: 0 - success, otherwise - failure 1240 1241 @ref User 1242 **/ 1243 int CeedOperatorLinearAssemblePointBlockDiagonal(CeedOperator op, 1244 CeedVector assembled, CeedRequest *request) { 1245 int ierr; 1246 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1247 1248 // Use backend version, if available 1249 if (op->LinearAssemblePointBlockDiagonal) { 1250 return op->LinearAssemblePointBlockDiagonal(op, assembled, request); 1251 } else if (op->LinearAssembleAddPointBlockDiagonal) { 1252 ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 1253 return CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, 1254 request); 1255 } else { 1256 // Fallback to reference Ceed 1257 if (!op->opfallback) { 1258 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1259 } 1260 // Assemble 1261 return CeedOperatorLinearAssemblePointBlockDiagonal(op->opfallback, 1262 assembled, request); 1263 } 1264 } 1265 1266 /** 1267 @brief Assemble the point block diagonal of a square linear CeedOperator 1268 1269 This sums into a CeedVector with the point block diagonal of a linear 1270 CeedOperator. 1271 1272 Note: Currently only non-composite CeedOperators with a single field and 1273 composite CeedOperators with single field sub-operators are supported. 1274 1275 @param op CeedOperator to assemble CeedQFunction 1276 @param[out] assembled CeedVector to store assembled CeedOperator point block 1277 diagonal, provided in row-major form with an 1278 @a ncomp * @a ncomp block at each node. The dimensions 1279 of this vector are derived from the active vector 1280 for the CeedOperator. The array has shape 1281 [nodes, component out, component in]. 1282 @param request Address of CeedRequest for non-blocking completion, else 1283 CEED_REQUEST_IMMEDIATE 1284 1285 @return An error code: 0 - success, otherwise - failure 1286 1287 @ref User 1288 **/ 1289 int CeedOperatorLinearAssembleAddPointBlockDiagonal(CeedOperator op, 1290 CeedVector assembled, CeedRequest *request) { 1291 int ierr; 1292 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1293 1294 // Use backend version, if available 1295 if (op->LinearAssembleAddPointBlockDiagonal) { 1296 return op->LinearAssembleAddPointBlockDiagonal(op, assembled, request); 1297 } else { 1298 // Fallback to reference Ceed 1299 if (!op->opfallback) { 1300 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1301 } 1302 // Assemble 1303 return CeedOperatorLinearAssembleAddPointBlockDiagonal(op->opfallback, 1304 assembled, request); 1305 } 1306 } 1307 1308 /** 1309 @brief Build nonzero pattern for non-composite operator. 1310 1311 Users should generally use CeedOperatorLinearAssembleSymbolic() 1312 1313 @ref Developer 1314 **/ 1315 int CeedSingleOperatorAssembleSymbolic(CeedOperator op, CeedInt offset, 1316 CeedInt *rows, CeedInt *cols) { 1317 int ierr; 1318 Ceed ceed = op->ceed; 1319 if (op->composite) 1320 // LCOV_EXCL_START 1321 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 1322 "Composite operator not supported"); 1323 // LCOV_EXCL_STOP 1324 1325 CeedElemRestriction rstrin; 1326 ierr = CeedOperatorGetActiveElemRestriction(op, &rstrin); CeedChk(ierr); 1327 CeedInt nelem, elemsize, nnodes, ncomp; 1328 ierr = CeedElemRestrictionGetNumElements(rstrin, &nelem); CeedChk(ierr); 1329 ierr = CeedElemRestrictionGetElementSize(rstrin, &elemsize); CeedChk(ierr); 1330 ierr = CeedElemRestrictionGetLVectorSize(rstrin, &nnodes); CeedChk(ierr); 1331 ierr = CeedElemRestrictionGetNumComponents(rstrin, &ncomp); CeedChk(ierr); 1332 CeedInt layout_er[3]; 1333 ierr = CeedElemRestrictionGetELayout(rstrin, &layout_er); CeedChk(ierr); 1334 1335 CeedInt local_nentries = elemsize*ncomp * elemsize*ncomp * nelem; 1336 1337 // Determine elem_dof relation 1338 CeedVector index_vec; 1339 ierr = CeedVectorCreate(ceed, nnodes, &index_vec); CeedChk(ierr); 1340 CeedScalar *array; 1341 ierr = CeedVectorGetArray(index_vec, CEED_MEM_HOST, &array); CeedChk(ierr); 1342 for (CeedInt i = 0; i < nnodes; ++i) { 1343 array[i] = i; 1344 } 1345 ierr = CeedVectorRestoreArray(index_vec, &array); CeedChk(ierr); 1346 CeedVector elem_dof; 1347 ierr = CeedVectorCreate(ceed, nelem * elemsize * ncomp, &elem_dof); 1348 CeedChk(ierr); 1349 ierr = CeedVectorSetValue(elem_dof, 0.0); CeedChk(ierr); 1350 CeedElemRestrictionApply(rstrin, CEED_NOTRANSPOSE, index_vec, 1351 elem_dof, CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 1352 const CeedScalar *elem_dof_a; 1353 ierr = CeedVectorGetArrayRead(elem_dof, CEED_MEM_HOST, &elem_dof_a); 1354 CeedChk(ierr); 1355 ierr = CeedVectorDestroy(&index_vec); CeedChk(ierr); 1356 1357 // Determine i, j locations for element matrices 1358 CeedInt count = 0; 1359 for (int e = 0; e < nelem; ++e) { 1360 for (int compin = 0; compin < ncomp; ++compin) { 1361 for (int compout = 0; compout < ncomp; ++compout) { 1362 for (int i = 0; i < elemsize; ++i) { 1363 for (int j = 0; j < elemsize; ++j) { 1364 const CeedInt edof_index_row = (i)*layout_er[0] + 1365 (compout)*layout_er[1] + e*layout_er[2]; 1366 const CeedInt edof_index_col = (j)*layout_er[0] + 1367 (compin)*layout_er[1] + e*layout_er[2]; 1368 1369 const CeedInt row = elem_dof_a[edof_index_row]; 1370 const CeedInt col = elem_dof_a[edof_index_col]; 1371 1372 rows[offset + count] = row; 1373 cols[offset + count] = col; 1374 count++; 1375 } 1376 } 1377 } 1378 } 1379 } 1380 if (count != local_nentries) 1381 // LCOV_EXCL_START 1382 return CeedError(ceed, CEED_ERROR_MAJOR, "Error computing assembled entries"); 1383 // LCOV_EXCL_STOP 1384 ierr = CeedVectorRestoreArrayRead(elem_dof, &elem_dof_a); CeedChk(ierr); 1385 ierr = CeedVectorDestroy(&elem_dof); CeedChk(ierr); 1386 1387 return CEED_ERROR_SUCCESS; 1388 } 1389 1390 /** 1391 @brief Assemble nonzero entries for non-composite operator 1392 1393 Users should generally use CeedOperatorLinearAssemble() 1394 1395 @ref Developer 1396 **/ 1397 int CeedSingleOperatorAssemble(CeedOperator op, CeedInt offset, 1398 CeedVector values) { 1399 int ierr; 1400 Ceed ceed = op->ceed;; 1401 if (op->composite) 1402 // LCOV_EXCL_START 1403 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 1404 "Composite operator not supported"); 1405 // LCOV_EXCL_STOP 1406 1407 // Assemble QFunction 1408 CeedQFunction qf; 1409 ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr); 1410 CeedInt numinputfields, numoutputfields; 1411 ierr= CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields); 1412 CeedChk(ierr); 1413 CeedVector assembledqf; 1414 CeedElemRestriction rstr_q; 1415 ierr = CeedOperatorLinearAssembleQFunction( 1416 op, &assembledqf, &rstr_q, CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 1417 1418 CeedInt qflength; 1419 ierr = CeedVectorGetLength(assembledqf, &qflength); CeedChk(ierr); 1420 1421 CeedOperatorField *input_fields; 1422 CeedOperatorField *output_fields; 1423 ierr = CeedOperatorGetFields(op, &input_fields, &output_fields); CeedChk(ierr); 1424 1425 // Determine active input basis 1426 CeedQFunctionField *qffields; 1427 ierr = CeedQFunctionGetFields(qf, &qffields, NULL); CeedChk(ierr); 1428 CeedInt numemodein = 0, dim = 1; 1429 CeedEvalMode *emodein = NULL; 1430 CeedBasis basisin = NULL; 1431 CeedElemRestriction rstrin = NULL; 1432 for (CeedInt i=0; i<numinputfields; i++) { 1433 CeedVector vec; 1434 ierr = CeedOperatorFieldGetVector(input_fields[i], &vec); CeedChk(ierr); 1435 if (vec == CEED_VECTOR_ACTIVE) { 1436 ierr = CeedOperatorFieldGetBasis(input_fields[i], &basisin); 1437 CeedChk(ierr); 1438 ierr = CeedBasisGetDimension(basisin, &dim); CeedChk(ierr); 1439 ierr = CeedOperatorFieldGetElemRestriction(input_fields[i], &rstrin); 1440 CeedChk(ierr); 1441 CeedEvalMode emode; 1442 ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode); 1443 CeedChk(ierr); 1444 switch (emode) { 1445 case CEED_EVAL_NONE: 1446 case CEED_EVAL_INTERP: 1447 ierr = CeedRealloc(numemodein + 1, &emodein); CeedChk(ierr); 1448 emodein[numemodein] = emode; 1449 numemodein += 1; 1450 break; 1451 case CEED_EVAL_GRAD: 1452 ierr = CeedRealloc(numemodein + dim, &emodein); CeedChk(ierr); 1453 for (CeedInt d=0; d<dim; d++) { 1454 emodein[numemodein+d] = emode; 1455 } 1456 numemodein += dim; 1457 break; 1458 case CEED_EVAL_WEIGHT: 1459 case CEED_EVAL_DIV: 1460 case CEED_EVAL_CURL: 1461 break; // Caught by QF Assembly 1462 } 1463 } 1464 } 1465 1466 // Determine active output basis 1467 ierr = CeedQFunctionGetFields(qf, NULL, &qffields); CeedChk(ierr); 1468 CeedInt numemodeout = 0; 1469 CeedEvalMode *emodeout = NULL; 1470 CeedBasis basisout = NULL; 1471 CeedElemRestriction rstrout = NULL; 1472 for (CeedInt i=0; i<numoutputfields; i++) { 1473 CeedVector vec; 1474 ierr = CeedOperatorFieldGetVector(output_fields[i], &vec); CeedChk(ierr); 1475 if (vec == CEED_VECTOR_ACTIVE) { 1476 ierr = CeedOperatorFieldGetBasis(output_fields[i], &basisout); 1477 CeedChk(ierr); 1478 ierr = CeedOperatorFieldGetElemRestriction(output_fields[i], &rstrout); 1479 CeedChk(ierr); 1480 CeedEvalMode emode; 1481 ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode); 1482 CeedChk(ierr); 1483 switch (emode) { 1484 case CEED_EVAL_NONE: 1485 case CEED_EVAL_INTERP: 1486 ierr = CeedRealloc(numemodeout + 1, &emodeout); CeedChk(ierr); 1487 emodeout[numemodeout] = emode; 1488 numemodeout += 1; 1489 break; 1490 case CEED_EVAL_GRAD: 1491 ierr = CeedRealloc(numemodeout + dim, &emodeout); CeedChk(ierr); 1492 for (CeedInt d=0; d<dim; d++) { 1493 emodeout[numemodeout+d] = emode; 1494 } 1495 numemodeout += dim; 1496 break; 1497 case CEED_EVAL_WEIGHT: 1498 case CEED_EVAL_DIV: 1499 case CEED_EVAL_CURL: 1500 break; // Caught by QF Assembly 1501 } 1502 } 1503 } 1504 1505 CeedInt nelem, elemsize, nqpts, ncomp; 1506 ierr = CeedElemRestrictionGetNumElements(rstrin, &nelem); CeedChk(ierr); 1507 ierr = CeedElemRestrictionGetElementSize(rstrin, &elemsize); CeedChk(ierr); 1508 ierr = CeedElemRestrictionGetNumComponents(rstrin, &ncomp); CeedChk(ierr); 1509 ierr = CeedBasisGetNumQuadraturePoints(basisin, &nqpts); CeedChk(ierr); 1510 1511 CeedInt local_nentries = elemsize*ncomp * elemsize*ncomp * nelem; 1512 1513 // loop over elements and put in data structure 1514 const CeedScalar *interpin, *gradin; 1515 ierr = CeedBasisGetInterp(basisin, &interpin); CeedChk(ierr); 1516 ierr = CeedBasisGetGrad(basisin, &gradin); CeedChk(ierr); 1517 1518 const CeedScalar *assembledqfarray; 1519 ierr = CeedVectorGetArrayRead(assembledqf, CEED_MEM_HOST, &assembledqfarray); 1520 CeedChk(ierr); 1521 1522 CeedInt layout_qf[3]; 1523 ierr = CeedElemRestrictionGetELayout(rstr_q, &layout_qf); CeedChk(ierr); 1524 ierr = CeedElemRestrictionDestroy(&rstr_q); CeedChk(ierr); 1525 1526 // we store Bmat_in, Bmat_out, BTD, elem_mat in row-major order 1527 CeedScalar Bmat_in[(nqpts * numemodein) * elemsize]; 1528 CeedScalar Bmat_out[(nqpts * numemodeout) * elemsize]; 1529 CeedScalar Dmat[numemodeout * numemodein * nqpts]; // logically 3-tensor 1530 CeedScalar BTD[elemsize * nqpts*numemodein]; 1531 CeedScalar elem_mat[elemsize * elemsize]; 1532 int count = 0; 1533 CeedScalar *vals; 1534 ierr = CeedVectorGetArray(values, CEED_MEM_HOST, &vals); CeedChk(ierr); 1535 for (int e = 0; e < nelem; ++e) { 1536 for (int compin = 0; compin < ncomp; ++compin) { 1537 for (int compout = 0; compout < ncomp; ++compout) { 1538 for (int ell = 0; ell < (nqpts * numemodein) * elemsize; ++ell) { 1539 Bmat_in[ell] = 0.0; 1540 } 1541 for (int ell = 0; ell < (nqpts * numemodeout) * elemsize; ++ell) { 1542 Bmat_out[ell] = 0.0; 1543 } 1544 // Store block-diagonal D matrix as collection of small dense blocks 1545 for (int ell = 0; ell < numemodein*numemodeout*nqpts; ++ell) { 1546 Dmat[ell] = 0.0; 1547 } 1548 // form element matrix itself (for each block component) 1549 for (int ell = 0; ell < elemsize*elemsize; ++ell) { 1550 elem_mat[ell] = 0.0; 1551 } 1552 for (int q = 0; q < nqpts; ++q) { 1553 for (int n = 0; n < elemsize; ++n) { 1554 CeedInt din = -1; 1555 for (int ein = 0; ein < numemodein; ++ein) { 1556 const int qq = numemodein*q; 1557 if (emodein[ein] == CEED_EVAL_INTERP) { 1558 Bmat_in[(qq+ein)*elemsize + n] += interpin[q * elemsize + n]; 1559 } else if (emodein[ein] == CEED_EVAL_GRAD) { 1560 din += 1; 1561 Bmat_in[(qq+ein)*elemsize + n] += 1562 gradin[(din*nqpts+q) * elemsize + n]; 1563 } else { 1564 // LCOV_EXCL_START 1565 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Not implemented!"); 1566 // LCOV_EXCL_STOP 1567 } 1568 } 1569 CeedInt dout = -1; 1570 for (int eout = 0; eout < numemodeout; ++eout) { 1571 const int qq = numemodeout*q; 1572 if (emodeout[eout] == CEED_EVAL_INTERP) { 1573 Bmat_out[(qq+eout)*elemsize + n] += interpin[q * elemsize + n]; 1574 } else if (emodeout[eout] == CEED_EVAL_GRAD) { 1575 dout += 1; 1576 Bmat_out[(qq+eout)*elemsize + n] += 1577 gradin[(dout*nqpts+q) * elemsize + n]; 1578 } else { 1579 // LCOV_EXCL_START 1580 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Not implemented!"); 1581 // LCOV_EXCL_STOP 1582 } 1583 } 1584 } 1585 for (int ei = 0; ei < numemodeout; ++ei) { 1586 for (int ej = 0; ej < numemodein; ++ej) { 1587 const int emode_index = ((ei*ncomp+compin)*numemodein+ej)*ncomp+compout; 1588 const int index = q*layout_qf[0] + emode_index*layout_qf[1] + e*layout_qf[2]; 1589 Dmat[(ei*numemodein+ej)*nqpts + q] += assembledqfarray[index]; 1590 } 1591 } 1592 } 1593 // Compute B^T*D 1594 for (int ell = 0; ell < elemsize*nqpts*numemodein; ++ell) { 1595 BTD[ell] = 0.0; 1596 } 1597 for (int j = 0; j<elemsize; ++j) { 1598 for (int q = 0; q<nqpts; ++q) { 1599 int qq = numemodeout*q; 1600 for (int ei = 0; ei < numemodein; ++ei) { 1601 for (int ej = 0; ej < numemodeout; ++ej) { 1602 BTD[j*(nqpts*numemodein) + (qq+ei)] += 1603 Bmat_out[(qq+ej)*elemsize + j] * Dmat[(ei*numemodein+ej)*nqpts + q]; 1604 } 1605 } 1606 } 1607 } 1608 1609 ierr = CeedMatrixMultiply(ceed, BTD, Bmat_in, elem_mat, elemsize, 1610 elemsize, nqpts*numemodein); CeedChk(ierr); 1611 1612 // put element matrix in coordinate data structure 1613 for (int i = 0; i < elemsize; ++i) { 1614 for (int j = 0; j < elemsize; ++j) { 1615 vals[offset + count] = elem_mat[i*elemsize + j]; 1616 count++; 1617 } 1618 } 1619 } 1620 } 1621 } 1622 if (count != local_nentries) 1623 // LCOV_EXCL_START 1624 return CeedError(ceed, CEED_ERROR_MAJOR, "Error computing entries"); 1625 // LCOV_EXCL_STOP 1626 ierr = CeedVectorRestoreArray(values, &vals); CeedChk(ierr); 1627 1628 ierr = CeedVectorRestoreArrayRead(assembledqf, &assembledqfarray); 1629 CeedChk(ierr); 1630 ierr = CeedVectorDestroy(&assembledqf); CeedChk(ierr); 1631 ierr = CeedFree(&emodein); CeedChk(ierr); 1632 ierr = CeedFree(&emodeout); CeedChk(ierr); 1633 1634 return CEED_ERROR_SUCCESS; 1635 } 1636 1637 /** 1638 @ref Utility 1639 **/ 1640 int CeedSingleOperatorAssemblyCountEntries(CeedOperator op, CeedInt *nentries) { 1641 int ierr; 1642 CeedElemRestriction rstr; 1643 CeedInt nelem, elemsize, ncomp; 1644 1645 if (op->composite) 1646 // LCOV_EXCL_START 1647 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, 1648 "Composite operator not supported"); 1649 // LCOV_EXCL_STOP 1650 ierr = CeedOperatorGetActiveElemRestriction(op, &rstr); CeedChk(ierr); 1651 ierr = CeedElemRestrictionGetNumElements(rstr, &nelem); CeedChk(ierr); 1652 ierr = CeedElemRestrictionGetElementSize(rstr, &elemsize); CeedChk(ierr); 1653 ierr = CeedElemRestrictionGetNumComponents(rstr, &ncomp); CeedChk(ierr); 1654 *nentries = elemsize*ncomp * elemsize*ncomp * nelem; 1655 1656 return CEED_ERROR_SUCCESS; 1657 } 1658 1659 /** 1660 @brief Fully assemble the nonzero pattern of a linear operator. 1661 1662 Expected to be used in conjunction with CeedOperatorLinearAssemble() 1663 1664 The assembly routines use coordinate format, with nentries tuples of the 1665 form (i, j, value) which indicate that value should be added to the matrix 1666 in entry (i, j). Note that the (i, j) pairs are not unique and may repeat. 1667 This function returns the number of entries and their (i, j) locations, 1668 while CeedOperatorLinearAssemble() provides the values in the same 1669 ordering. 1670 1671 This will generally be slow unless your operator is low-order. 1672 1673 @param[in] op CeedOperator to assemble 1674 @param[out] nentries Number of entries in coordinate nonzero pattern. 1675 @param[out] rows Row number for each entry. 1676 @param[out] cols Column number for each entry. 1677 1678 @ref User 1679 **/ 1680 int CeedOperatorLinearAssembleSymbolic(CeedOperator op, 1681 CeedInt *nentries, CeedInt **rows, CeedInt **cols) { 1682 int ierr; 1683 CeedInt numsub, single_entries; 1684 CeedOperator *suboperators; 1685 bool isComposite; 1686 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1687 1688 // Use backend version, if available 1689 if (op->LinearAssembleSymbolic) { 1690 return op->LinearAssembleSymbolic(op, nentries, rows, cols); 1691 } else { 1692 // Check for valid fallback resource 1693 const char *resource, *fallbackresource; 1694 ierr = CeedGetResource(op->ceed, &resource); CeedChk(ierr); 1695 ierr = CeedGetOperatorFallbackResource(op->ceed, &fallbackresource); 1696 if (strcmp(fallbackresource, "") && strcmp(resource, fallbackresource)) { 1697 // Fallback to reference Ceed 1698 if (!op->opfallback) { 1699 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1700 } 1701 // Assemble 1702 return CeedOperatorLinearAssembleSymbolic(op->opfallback, nentries, rows, cols); 1703 } 1704 } 1705 1706 // if neither backend nor fallback resource provides 1707 // LinearAssembleSymbolic, continue with interface-level implementation 1708 1709 // count entries and allocate rows, cols arrays 1710 ierr = CeedOperatorIsComposite(op, &isComposite); CeedChk(ierr); 1711 *nentries = 0; 1712 if (isComposite) { 1713 ierr = CeedOperatorGetNumSub(op, &numsub); CeedChk(ierr); 1714 ierr = CeedOperatorGetSubList(op, &suboperators); CeedChk(ierr); 1715 for (int k = 0; k < numsub; ++k) { 1716 ierr = CeedSingleOperatorAssemblyCountEntries(suboperators[k], 1717 &single_entries); CeedChk(ierr); 1718 *nentries += single_entries; 1719 } 1720 } else { 1721 ierr = CeedSingleOperatorAssemblyCountEntries(op, 1722 &single_entries); CeedChk(ierr); 1723 *nentries += single_entries; 1724 } 1725 ierr = CeedCalloc(*nentries, rows); CeedChk(ierr); 1726 ierr = CeedCalloc(*nentries, cols); CeedChk(ierr); 1727 1728 // assemble nonzero locations 1729 CeedInt offset = 0; 1730 if (isComposite) { 1731 ierr = CeedOperatorGetNumSub(op, &numsub); CeedChk(ierr); 1732 ierr = CeedOperatorGetSubList(op, &suboperators); CeedChk(ierr); 1733 for (int k = 0; k < numsub; ++k) { 1734 ierr = CeedSingleOperatorAssembleSymbolic(suboperators[k], offset, *rows, 1735 *cols); CeedChk(ierr); 1736 ierr = CeedSingleOperatorAssemblyCountEntries(suboperators[k], &single_entries); 1737 CeedChk(ierr); 1738 offset += single_entries; 1739 } 1740 } else { 1741 ierr = CeedSingleOperatorAssembleSymbolic(op, offset, *rows, *cols); 1742 CeedChk(ierr); 1743 } 1744 1745 return CEED_ERROR_SUCCESS; 1746 } 1747 1748 /** 1749 @brief Fully assemble the nonzero entries of a linear operator. 1750 1751 Expected to be used in conjunction with CeedOperatorLinearAssembleSymbolic() 1752 1753 The assembly routines use coordinate format, with nentries tuples of the 1754 form (i, j, value) which indicate that value should be added to the matrix 1755 in entry (i, j). Note that the (i, j) pairs are not unique and may repeat. 1756 This function returns the values of the nonzero entries to be added, their 1757 (i, j) locations are provided by CeedOperatorLinearAssembleSymbolic() 1758 1759 This will generally be slow unless your operator is low-order. 1760 1761 @param[in] op CeedOperator to assemble 1762 @param[out] values Values to assemble into matrix 1763 1764 @ref User 1765 **/ 1766 int CeedOperatorLinearAssemble(CeedOperator op, CeedVector values) { 1767 int ierr; 1768 CeedInt numsub, single_entries; 1769 CeedOperator *suboperators; 1770 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1771 1772 // Use backend version, if available 1773 if (op->LinearAssemble) { 1774 return op->LinearAssemble(op, values); 1775 } else { 1776 // Check for valid fallback resource 1777 const char *resource, *fallbackresource; 1778 ierr = CeedGetResource(op->ceed, &resource); CeedChk(ierr); 1779 ierr = CeedGetOperatorFallbackResource(op->ceed, &fallbackresource); 1780 if (strcmp(fallbackresource, "") && strcmp(resource, fallbackresource)) { 1781 // Fallback to reference Ceed 1782 if (!op->opfallback) { 1783 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1784 } 1785 // Assemble 1786 return CeedOperatorLinearAssemble(op->opfallback, values); 1787 } 1788 } 1789 1790 // if neither backend nor fallback resource provides 1791 // LinearAssemble, continue with interface-level implementation 1792 1793 bool isComposite; 1794 ierr = CeedOperatorIsComposite(op, &isComposite); CeedChk(ierr); 1795 1796 CeedInt offset = 0; 1797 if (isComposite) { 1798 ierr = CeedOperatorGetNumSub(op, &numsub); CeedChk(ierr); 1799 ierr = CeedOperatorGetSubList(op, &suboperators); CeedChk(ierr); 1800 for (int k = 0; k < numsub; ++k) { 1801 ierr = CeedSingleOperatorAssemble(suboperators[k], offset, values); 1802 CeedChk(ierr); 1803 ierr = CeedSingleOperatorAssemblyCountEntries(suboperators[k], &single_entries); 1804 CeedChk(ierr); 1805 offset += single_entries; 1806 } 1807 } else { 1808 ierr = CeedSingleOperatorAssemble(op, offset, values); CeedChk(ierr); 1809 } 1810 1811 return CEED_ERROR_SUCCESS; 1812 } 1813 1814 /** 1815 @brief Create a multigrid coarse operator and level transfer operators 1816 for a CeedOperator, creating the prolongation basis from the 1817 fine and coarse grid interpolation 1818 1819 @param[in] opFine Fine grid operator 1820 @param[in] PMultFine L-vector multiplicity in parallel gather/scatter 1821 @param[in] rstrCoarse Coarse grid restriction 1822 @param[in] basisCoarse Coarse grid active vector basis 1823 @param[out] opCoarse Coarse grid operator 1824 @param[out] opProlong Coarse to fine operator 1825 @param[out] opRestrict Fine to coarse operator 1826 1827 @return An error code: 0 - success, otherwise - failure 1828 1829 @ref User 1830 **/ 1831 int CeedOperatorMultigridLevelCreate(CeedOperator opFine, CeedVector PMultFine, 1832 CeedElemRestriction rstrCoarse, CeedBasis basisCoarse, 1833 CeedOperator *opCoarse, CeedOperator *opProlong, CeedOperator *opRestrict) { 1834 int ierr; 1835 Ceed ceed; 1836 ierr = CeedOperatorGetCeed(opFine, &ceed); CeedChk(ierr); 1837 1838 // Check for compatible quadrature spaces 1839 CeedBasis basisFine; 1840 ierr = CeedOperatorGetActiveBasis(opFine, &basisFine); CeedChk(ierr); 1841 CeedInt Qf, Qc; 1842 ierr = CeedBasisGetNumQuadraturePoints(basisFine, &Qf); CeedChk(ierr); 1843 ierr = CeedBasisGetNumQuadraturePoints(basisCoarse, &Qc); CeedChk(ierr); 1844 if (Qf != Qc) 1845 // LCOV_EXCL_START 1846 return CeedError(ceed, CEED_ERROR_DIMENSION, 1847 "Bases must have compatible quadrature spaces"); 1848 // LCOV_EXCL_STOP 1849 1850 // Coarse to fine basis 1851 CeedInt Pf, Pc, Q = Qf; 1852 bool isTensorF, isTensorC; 1853 ierr = CeedBasisIsTensor(basisFine, &isTensorF); CeedChk(ierr); 1854 ierr = CeedBasisIsTensor(basisCoarse, &isTensorC); CeedChk(ierr); 1855 CeedScalar *interpC, *interpF, *interpCtoF, *tau; 1856 if (isTensorF && isTensorC) { 1857 ierr = CeedBasisGetNumNodes1D(basisFine, &Pf); CeedChk(ierr); 1858 ierr = CeedBasisGetNumNodes1D(basisCoarse, &Pc); CeedChk(ierr); 1859 ierr = CeedBasisGetNumQuadraturePoints1D(basisCoarse, &Q); CeedChk(ierr); 1860 } else if (!isTensorF && !isTensorC) { 1861 ierr = CeedBasisGetNumNodes(basisFine, &Pf); CeedChk(ierr); 1862 ierr = CeedBasisGetNumNodes(basisCoarse, &Pc); CeedChk(ierr); 1863 } else { 1864 // LCOV_EXCL_START 1865 return CeedError(ceed, CEED_ERROR_MINOR, 1866 "Bases must both be tensor or non-tensor"); 1867 // LCOV_EXCL_STOP 1868 } 1869 1870 ierr = CeedMalloc(Q*Pf, &interpF); CeedChk(ierr); 1871 ierr = CeedMalloc(Q*Pc, &interpC); CeedChk(ierr); 1872 ierr = CeedCalloc(Pc*Pf, &interpCtoF); CeedChk(ierr); 1873 ierr = CeedMalloc(Q, &tau); CeedChk(ierr); 1874 if (isTensorF) { 1875 memcpy(interpF, basisFine->interp1d, Q*Pf*sizeof basisFine->interp1d[0]); 1876 memcpy(interpC, basisCoarse->interp1d, Q*Pc*sizeof basisCoarse->interp1d[0]); 1877 } else { 1878 memcpy(interpF, basisFine->interp, Q*Pf*sizeof basisFine->interp[0]); 1879 memcpy(interpC, basisCoarse->interp, Q*Pc*sizeof basisCoarse->interp[0]); 1880 } 1881 1882 // -- QR Factorization, interpF = Q R 1883 ierr = CeedQRFactorization(ceed, interpF, tau, Q, Pf); CeedChk(ierr); 1884 1885 // -- Apply Qtranspose, interpC = Qtranspose interpC 1886 CeedHouseholderApplyQ(interpC, interpF, tau, CEED_TRANSPOSE, 1887 Q, Pc, Pf, Pc, 1); 1888 1889 // -- Apply Rinv, interpCtoF = Rinv interpC 1890 for (CeedInt j=0; j<Pc; j++) { // Column j 1891 interpCtoF[j+Pc*(Pf-1)] = interpC[j+Pc*(Pf-1)]/interpF[Pf*Pf-1]; 1892 for (CeedInt i=Pf-2; i>=0; i--) { // Row i 1893 interpCtoF[j+Pc*i] = interpC[j+Pc*i]; 1894 for (CeedInt k=i+1; k<Pf; k++) 1895 interpCtoF[j+Pc*i] -= interpF[k+Pf*i]*interpCtoF[j+Pc*k]; 1896 interpCtoF[j+Pc*i] /= interpF[i+Pf*i]; 1897 } 1898 } 1899 ierr = CeedFree(&tau); CeedChk(ierr); 1900 ierr = CeedFree(&interpC); CeedChk(ierr); 1901 ierr = CeedFree(&interpF); CeedChk(ierr); 1902 1903 // Complete with interpCtoF versions of code 1904 if (isTensorF) { 1905 ierr = CeedOperatorMultigridLevelCreateTensorH1(opFine, PMultFine, 1906 rstrCoarse, basisCoarse, interpCtoF, opCoarse, opProlong, opRestrict); 1907 CeedChk(ierr); 1908 } else { 1909 ierr = CeedOperatorMultigridLevelCreateH1(opFine, PMultFine, 1910 rstrCoarse, basisCoarse, interpCtoF, opCoarse, opProlong, opRestrict); 1911 CeedChk(ierr); 1912 } 1913 1914 // Cleanup 1915 ierr = CeedFree(&interpCtoF); CeedChk(ierr); 1916 return CEED_ERROR_SUCCESS; 1917 } 1918 1919 /** 1920 @brief Create a multigrid coarse operator and level transfer operators 1921 for a CeedOperator with a tensor basis for the active basis 1922 1923 @param[in] opFine Fine grid operator 1924 @param[in] PMultFine L-vector multiplicity in parallel gather/scatter 1925 @param[in] rstrCoarse Coarse grid restriction 1926 @param[in] basisCoarse Coarse grid active vector basis 1927 @param[in] interpCtoF Matrix for coarse to fine interpolation 1928 @param[out] opCoarse Coarse grid operator 1929 @param[out] opProlong Coarse to fine operator 1930 @param[out] opRestrict Fine to coarse operator 1931 1932 @return An error code: 0 - success, otherwise - failure 1933 1934 @ref User 1935 **/ 1936 int CeedOperatorMultigridLevelCreateTensorH1(CeedOperator opFine, 1937 CeedVector PMultFine, CeedElemRestriction rstrCoarse, CeedBasis basisCoarse, 1938 const CeedScalar *interpCtoF, CeedOperator *opCoarse, 1939 CeedOperator *opProlong, CeedOperator *opRestrict) { 1940 int ierr; 1941 Ceed ceed; 1942 ierr = CeedOperatorGetCeed(opFine, &ceed); CeedChk(ierr); 1943 1944 // Check for compatible quadrature spaces 1945 CeedBasis basisFine; 1946 ierr = CeedOperatorGetActiveBasis(opFine, &basisFine); CeedChk(ierr); 1947 CeedInt Qf, Qc; 1948 ierr = CeedBasisGetNumQuadraturePoints(basisFine, &Qf); CeedChk(ierr); 1949 ierr = CeedBasisGetNumQuadraturePoints(basisCoarse, &Qc); CeedChk(ierr); 1950 if (Qf != Qc) 1951 // LCOV_EXCL_START 1952 return CeedError(ceed, CEED_ERROR_DIMENSION, 1953 "Bases must have compatible quadrature spaces"); 1954 // LCOV_EXCL_STOP 1955 1956 // Coarse to fine basis 1957 CeedInt dim, ncomp, nnodesCoarse, P1dFine, P1dCoarse; 1958 ierr = CeedBasisGetDimension(basisFine, &dim); CeedChk(ierr); 1959 ierr = CeedBasisGetNumComponents(basisFine, &ncomp); CeedChk(ierr); 1960 ierr = CeedBasisGetNumNodes1D(basisFine, &P1dFine); CeedChk(ierr); 1961 ierr = CeedElemRestrictionGetElementSize(rstrCoarse, &nnodesCoarse); 1962 CeedChk(ierr); 1963 P1dCoarse = dim == 1 ? nnodesCoarse : 1964 dim == 2 ? sqrt(nnodesCoarse) : 1965 cbrt(nnodesCoarse); 1966 CeedScalar *qref, *qweight, *grad; 1967 ierr = CeedCalloc(P1dFine, &qref); CeedChk(ierr); 1968 ierr = CeedCalloc(P1dFine, &qweight); CeedChk(ierr); 1969 ierr = CeedCalloc(P1dFine*P1dCoarse*dim, &grad); CeedChk(ierr); 1970 CeedBasis basisCtoF; 1971 ierr = CeedBasisCreateTensorH1(ceed, dim, ncomp, P1dCoarse, P1dFine, 1972 interpCtoF, grad, qref, qweight, &basisCtoF); 1973 CeedChk(ierr); 1974 ierr = CeedFree(&qref); CeedChk(ierr); 1975 ierr = CeedFree(&qweight); CeedChk(ierr); 1976 ierr = CeedFree(&grad); CeedChk(ierr); 1977 1978 // Core code 1979 ierr = CeedOperatorMultigridLevel_Core(opFine, PMultFine, rstrCoarse, 1980 basisCoarse, basisCtoF, opCoarse, 1981 opProlong, opRestrict); 1982 CeedChk(ierr); 1983 return CEED_ERROR_SUCCESS; 1984 } 1985 1986 /** 1987 @brief Create a multigrid coarse operator and level transfer operators 1988 for a CeedOperator with a non-tensor basis for the active vector 1989 1990 @param[in] opFine Fine grid operator 1991 @param[in] PMultFine L-vector multiplicity in parallel gather/scatter 1992 @param[in] rstrCoarse Coarse grid restriction 1993 @param[in] basisCoarse Coarse grid active vector basis 1994 @param[in] interpCtoF Matrix for coarse to fine interpolation 1995 @param[out] opCoarse Coarse grid operator 1996 @param[out] opProlong Coarse to fine operator 1997 @param[out] opRestrict Fine to coarse operator 1998 1999 @return An error code: 0 - success, otherwise - failure 2000 2001 @ref User 2002 **/ 2003 int CeedOperatorMultigridLevelCreateH1(CeedOperator opFine, 2004 CeedVector PMultFine, 2005 CeedElemRestriction rstrCoarse, 2006 CeedBasis basisCoarse, 2007 const CeedScalar *interpCtoF, 2008 CeedOperator *opCoarse, 2009 CeedOperator *opProlong, 2010 CeedOperator *opRestrict) { 2011 int ierr; 2012 Ceed ceed; 2013 ierr = CeedOperatorGetCeed(opFine, &ceed); CeedChk(ierr); 2014 2015 // Check for compatible quadrature spaces 2016 CeedBasis basisFine; 2017 ierr = CeedOperatorGetActiveBasis(opFine, &basisFine); CeedChk(ierr); 2018 CeedInt Qf, Qc; 2019 ierr = CeedBasisGetNumQuadraturePoints(basisFine, &Qf); CeedChk(ierr); 2020 ierr = CeedBasisGetNumQuadraturePoints(basisCoarse, &Qc); CeedChk(ierr); 2021 if (Qf != Qc) 2022 // LCOV_EXCL_START 2023 return CeedError(ceed, CEED_ERROR_DIMENSION, 2024 "Bases must have compatible quadrature spaces"); 2025 // LCOV_EXCL_STOP 2026 2027 // Coarse to fine basis 2028 CeedElemTopology topo; 2029 ierr = CeedBasisGetTopology(basisFine, &topo); CeedChk(ierr); 2030 CeedInt dim, ncomp, nnodesCoarse, nnodesFine; 2031 ierr = CeedBasisGetDimension(basisFine, &dim); CeedChk(ierr); 2032 ierr = CeedBasisGetNumComponents(basisFine, &ncomp); CeedChk(ierr); 2033 ierr = CeedBasisGetNumNodes(basisFine, &nnodesFine); CeedChk(ierr); 2034 ierr = CeedElemRestrictionGetElementSize(rstrCoarse, &nnodesCoarse); 2035 CeedChk(ierr); 2036 CeedScalar *qref, *qweight, *grad; 2037 ierr = CeedCalloc(nnodesFine, &qref); CeedChk(ierr); 2038 ierr = CeedCalloc(nnodesFine, &qweight); CeedChk(ierr); 2039 ierr = CeedCalloc(nnodesFine*nnodesCoarse*dim, &grad); CeedChk(ierr); 2040 CeedBasis basisCtoF; 2041 ierr = CeedBasisCreateH1(ceed, topo, ncomp, nnodesCoarse, nnodesFine, 2042 interpCtoF, grad, qref, qweight, &basisCtoF); 2043 CeedChk(ierr); 2044 ierr = CeedFree(&qref); CeedChk(ierr); 2045 ierr = CeedFree(&qweight); CeedChk(ierr); 2046 ierr = CeedFree(&grad); CeedChk(ierr); 2047 2048 // Core code 2049 ierr = CeedOperatorMultigridLevel_Core(opFine, PMultFine, rstrCoarse, 2050 basisCoarse, basisCtoF, opCoarse, 2051 opProlong, opRestrict); 2052 CeedChk(ierr); 2053 return CEED_ERROR_SUCCESS; 2054 } 2055 2056 /** 2057 @brief Build a FDM based approximate inverse for each element for a 2058 CeedOperator 2059 2060 This returns a CeedOperator and CeedVector to apply a Fast Diagonalization 2061 Method based approximate inverse. This function obtains the simultaneous 2062 diagonalization for the 1D mass and Laplacian operators, 2063 M = V^T V, K = V^T S V. 2064 The assembled QFunction is used to modify the eigenvalues from simultaneous 2065 diagonalization and obtain an approximate inverse of the form 2066 V^T S^hat V. The CeedOperator must be linear and non-composite. The 2067 associated CeedQFunction must therefore also be linear. 2068 2069 @param op CeedOperator to create element inverses 2070 @param[out] fdminv CeedOperator to apply the action of a FDM based inverse 2071 for each element 2072 @param request Address of CeedRequest for non-blocking completion, else 2073 @ref CEED_REQUEST_IMMEDIATE 2074 2075 @return An error code: 0 - success, otherwise - failure 2076 2077 @ref Backend 2078 **/ 2079 int CeedOperatorCreateFDMElementInverse(CeedOperator op, CeedOperator *fdminv, 2080 CeedRequest *request) { 2081 int ierr; 2082 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 2083 2084 // Use backend version, if available 2085 if (op->CreateFDMElementInverse) { 2086 ierr = op->CreateFDMElementInverse(op, fdminv, request); CeedChk(ierr); 2087 } else { 2088 // Fallback to reference Ceed 2089 if (!op->opfallback) { 2090 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 2091 } 2092 // Assemble 2093 ierr = op->opfallback->CreateFDMElementInverse(op->opfallback, fdminv, 2094 request); CeedChk(ierr); 2095 } 2096 return CEED_ERROR_SUCCESS; 2097 } 2098 2099 /** 2100 @brief View a CeedOperator 2101 2102 @param[in] op CeedOperator to view 2103 @param[in] stream Stream to write; typically stdout/stderr or a file 2104 2105 @return Error code: 0 - success, otherwise - failure 2106 2107 @ref User 2108 **/ 2109 int CeedOperatorView(CeedOperator op, FILE *stream) { 2110 int ierr; 2111 2112 if (op->composite) { 2113 fprintf(stream, "Composite CeedOperator\n"); 2114 2115 for (CeedInt i=0; i<op->numsub; i++) { 2116 fprintf(stream, " SubOperator [%d]:\n", i); 2117 ierr = CeedOperatorSingleView(op->suboperators[i], 1, stream); 2118 CeedChk(ierr); 2119 } 2120 } else { 2121 fprintf(stream, "CeedOperator\n"); 2122 ierr = CeedOperatorSingleView(op, 0, stream); CeedChk(ierr); 2123 } 2124 return CEED_ERROR_SUCCESS; 2125 } 2126 2127 /** 2128 @brief Apply CeedOperator to a vector 2129 2130 This computes the action of the operator on the specified (active) input, 2131 yielding its (active) output. All inputs and outputs must be specified using 2132 CeedOperatorSetField(). 2133 2134 @param op CeedOperator to apply 2135 @param[in] in CeedVector containing input state or @ref CEED_VECTOR_NONE if 2136 there are no active inputs 2137 @param[out] out CeedVector to store result of applying operator (must be 2138 distinct from @a in) or @ref CEED_VECTOR_NONE if there are no 2139 active outputs 2140 @param request Address of CeedRequest for non-blocking completion, else 2141 @ref CEED_REQUEST_IMMEDIATE 2142 2143 @return An error code: 0 - success, otherwise - failure 2144 2145 @ref User 2146 **/ 2147 int CeedOperatorApply(CeedOperator op, CeedVector in, CeedVector out, 2148 CeedRequest *request) { 2149 int ierr; 2150 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 2151 2152 if (op->numelements) { 2153 // Standard Operator 2154 if (op->Apply) { 2155 ierr = op->Apply(op, in, out, request); CeedChk(ierr); 2156 } else { 2157 // Zero all output vectors 2158 CeedQFunction qf = op->qf; 2159 for (CeedInt i=0; i<qf->numoutputfields; i++) { 2160 CeedVector vec = op->outputfields[i]->vec; 2161 if (vec == CEED_VECTOR_ACTIVE) 2162 vec = out; 2163 if (vec != CEED_VECTOR_NONE) { 2164 ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr); 2165 } 2166 } 2167 // Apply 2168 ierr = op->ApplyAdd(op, in, out, request); CeedChk(ierr); 2169 } 2170 } else if (op->composite) { 2171 // Composite Operator 2172 if (op->ApplyComposite) { 2173 ierr = op->ApplyComposite(op, in, out, request); CeedChk(ierr); 2174 } else { 2175 CeedInt numsub; 2176 ierr = CeedOperatorGetNumSub(op, &numsub); CeedChk(ierr); 2177 CeedOperator *suboperators; 2178 ierr = CeedOperatorGetSubList(op, &suboperators); CeedChk(ierr); 2179 2180 // Zero all output vectors 2181 if (out != CEED_VECTOR_NONE) { 2182 ierr = CeedVectorSetValue(out, 0.0); CeedChk(ierr); 2183 } 2184 for (CeedInt i=0; i<numsub; i++) { 2185 for (CeedInt j=0; j<suboperators[i]->qf->numoutputfields; j++) { 2186 CeedVector vec = suboperators[i]->outputfields[j]->vec; 2187 if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) { 2188 ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr); 2189 } 2190 } 2191 } 2192 // Apply 2193 for (CeedInt i=0; i<op->numsub; i++) { 2194 ierr = CeedOperatorApplyAdd(op->suboperators[i], in, out, request); 2195 CeedChk(ierr); 2196 } 2197 } 2198 } 2199 return CEED_ERROR_SUCCESS; 2200 } 2201 2202 /** 2203 @brief Apply CeedOperator to a vector and add result to output vector 2204 2205 This computes the action of the operator on the specified (active) input, 2206 yielding its (active) output. All inputs and outputs must be specified using 2207 CeedOperatorSetField(). 2208 2209 @param op CeedOperator to apply 2210 @param[in] in CeedVector containing input state or NULL if there are no 2211 active inputs 2212 @param[out] out CeedVector to sum in result of applying operator (must be 2213 distinct from @a in) or NULL if there are no active outputs 2214 @param request Address of CeedRequest for non-blocking completion, else 2215 @ref CEED_REQUEST_IMMEDIATE 2216 2217 @return An error code: 0 - success, otherwise - failure 2218 2219 @ref User 2220 **/ 2221 int CeedOperatorApplyAdd(CeedOperator op, CeedVector in, CeedVector out, 2222 CeedRequest *request) { 2223 int ierr; 2224 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 2225 2226 if (op->numelements) { 2227 // Standard Operator 2228 ierr = op->ApplyAdd(op, in, out, request); CeedChk(ierr); 2229 } else if (op->composite) { 2230 // Composite Operator 2231 if (op->ApplyAddComposite) { 2232 ierr = op->ApplyAddComposite(op, in, out, request); CeedChk(ierr); 2233 } else { 2234 CeedInt numsub; 2235 ierr = CeedOperatorGetNumSub(op, &numsub); CeedChk(ierr); 2236 CeedOperator *suboperators; 2237 ierr = CeedOperatorGetSubList(op, &suboperators); CeedChk(ierr); 2238 2239 for (CeedInt i=0; i<numsub; i++) { 2240 ierr = CeedOperatorApplyAdd(suboperators[i], in, out, request); 2241 CeedChk(ierr); 2242 } 2243 } 2244 } 2245 return CEED_ERROR_SUCCESS; 2246 } 2247 2248 /** 2249 @brief Destroy a CeedOperator 2250 2251 @param op CeedOperator to destroy 2252 2253 @return An error code: 0 - success, otherwise - failure 2254 2255 @ref User 2256 **/ 2257 int CeedOperatorDestroy(CeedOperator *op) { 2258 int ierr; 2259 2260 if (!*op || --(*op)->refcount > 0) return CEED_ERROR_SUCCESS; 2261 if ((*op)->Destroy) { 2262 ierr = (*op)->Destroy(*op); CeedChk(ierr); 2263 } 2264 ierr = CeedDestroy(&(*op)->ceed); CeedChk(ierr); 2265 // Free fields 2266 for (int i=0; i<(*op)->nfields; i++) 2267 if ((*op)->inputfields[i]) { 2268 if ((*op)->inputfields[i]->Erestrict != CEED_ELEMRESTRICTION_NONE) { 2269 ierr = CeedElemRestrictionDestroy(&(*op)->inputfields[i]->Erestrict); 2270 CeedChk(ierr); 2271 } 2272 if ((*op)->inputfields[i]->basis != CEED_BASIS_COLLOCATED) { 2273 ierr = CeedBasisDestroy(&(*op)->inputfields[i]->basis); CeedChk(ierr); 2274 } 2275 if ((*op)->inputfields[i]->vec != CEED_VECTOR_ACTIVE && 2276 (*op)->inputfields[i]->vec != CEED_VECTOR_NONE ) { 2277 ierr = CeedVectorDestroy(&(*op)->inputfields[i]->vec); CeedChk(ierr); 2278 } 2279 ierr = CeedFree(&(*op)->inputfields[i]->fieldname); CeedChk(ierr); 2280 ierr = CeedFree(&(*op)->inputfields[i]); CeedChk(ierr); 2281 } 2282 for (int i=0; i<(*op)->nfields; i++) 2283 if ((*op)->outputfields[i]) { 2284 ierr = CeedElemRestrictionDestroy(&(*op)->outputfields[i]->Erestrict); 2285 CeedChk(ierr); 2286 if ((*op)->outputfields[i]->basis != CEED_BASIS_COLLOCATED) { 2287 ierr = CeedBasisDestroy(&(*op)->outputfields[i]->basis); CeedChk(ierr); 2288 } 2289 if ((*op)->outputfields[i]->vec != CEED_VECTOR_ACTIVE && 2290 (*op)->outputfields[i]->vec != CEED_VECTOR_NONE ) { 2291 ierr = CeedVectorDestroy(&(*op)->outputfields[i]->vec); CeedChk(ierr); 2292 } 2293 ierr = CeedFree(&(*op)->outputfields[i]->fieldname); CeedChk(ierr); 2294 ierr = CeedFree(&(*op)->outputfields[i]); CeedChk(ierr); 2295 } 2296 // Destroy suboperators 2297 for (int i=0; i<(*op)->numsub; i++) 2298 if ((*op)->suboperators[i]) { 2299 ierr = CeedOperatorDestroy(&(*op)->suboperators[i]); CeedChk(ierr); 2300 } 2301 if ((*op)->qf) 2302 // LCOV_EXCL_START 2303 (*op)->qf->operatorsset--; 2304 // LCOV_EXCL_STOP 2305 ierr = CeedQFunctionDestroy(&(*op)->qf); CeedChk(ierr); 2306 if ((*op)->dqf && (*op)->dqf != CEED_QFUNCTION_NONE) 2307 // LCOV_EXCL_START 2308 (*op)->dqf->operatorsset--; 2309 // LCOV_EXCL_STOP 2310 ierr = CeedQFunctionDestroy(&(*op)->dqf); CeedChk(ierr); 2311 if ((*op)->dqfT && (*op)->dqfT != CEED_QFUNCTION_NONE) 2312 // LCOV_EXCL_START 2313 (*op)->dqfT->operatorsset--; 2314 // LCOV_EXCL_STOP 2315 ierr = CeedQFunctionDestroy(&(*op)->dqfT); CeedChk(ierr); 2316 2317 // Destroy fallback 2318 if ((*op)->opfallback) { 2319 ierr = (*op)->qffallback->Destroy((*op)->qffallback); CeedChk(ierr); 2320 ierr = CeedFree(&(*op)->qffallback); CeedChk(ierr); 2321 ierr = (*op)->opfallback->Destroy((*op)->opfallback); CeedChk(ierr); 2322 ierr = CeedFree(&(*op)->opfallback); CeedChk(ierr); 2323 } 2324 2325 ierr = CeedFree(&(*op)->inputfields); CeedChk(ierr); 2326 ierr = CeedFree(&(*op)->outputfields); CeedChk(ierr); 2327 ierr = CeedFree(&(*op)->suboperators); CeedChk(ierr); 2328 ierr = CeedFree(op); CeedChk(ierr); 2329 return CEED_ERROR_SUCCESS; 2330 } 2331 2332 /// @} 2333