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