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