Actual source code: pcpatch.c

      1: #include <petsc/private/pcpatchimpl.h>
      2: #include <petsc/private/kspimpl.h>
      3: #include <petsc/private/vecimpl.h>
      4: #include <petsc/private/dmpleximpl.h>
      5: #include <petscsf.h>
      6: #include <petscbt.h>
      7: #include <petscds.h>
      8: #include <../src/mat/impls/dense/seq/dense.h>

     10: PetscBool  PCPatchcite       = PETSC_FALSE;
     11: const char PCPatchCitation[] = "@article{FarrellKnepleyWechsungMitchell2020,\n"
     12:                                "title   = {{PCPATCH}: software for the topological construction of multigrid relaxation methods},\n"
     13:                                "author  = {Patrick E Farrell and Matthew G Knepley and Lawrence Mitchell and Florian Wechsung},\n"
     14:                                "journal = {ACM Transaction on Mathematical Software},\n"
     15:                                "eprint  = {http://arxiv.org/abs/1912.08516},\n"
     16:                                "volume  = {47},\n"
     17:                                "number  = {3},\n"
     18:                                "pages   = {1--22},\n"
     19:                                "year    = {2021},\n"
     20:                                "petsc_uses={KSP,DMPlex}\n}\n";

     22: PetscLogEvent PC_Patch_CreatePatches, PC_Patch_ComputeOp, PC_Patch_Solve, PC_Patch_Apply, PC_Patch_Prealloc;

     24: static inline PetscErrorCode ObjectView(PetscObject obj, PetscViewer viewer, PetscViewerFormat format)
     25: {
     26:   PetscCall(PetscViewerPushFormat(viewer, format));
     27:   PetscCall(PetscObjectView(obj, viewer));
     28:   PetscCall(PetscViewerPopFormat(viewer));
     29:   return PETSC_SUCCESS;
     30: }

     32: static PetscErrorCode PCPatchConstruct_Star(void *vpatch, DM dm, PetscInt point, PetscHSetI ht)
     33: {
     34:   PetscInt  starSize;
     35:   PetscInt *star = NULL, si;

     37:   PetscFunctionBegin;
     38:   PetscCall(PetscHSetIClear(ht));
     39:   /* To start with, add the point we care about */
     40:   PetscCall(PetscHSetIAdd(ht, point));
     41:   /* Loop over all the points that this point connects to */
     42:   PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star));
     43:   for (si = 0; si < starSize * 2; si += 2) PetscCall(PetscHSetIAdd(ht, star[si]));
     44:   PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star));
     45:   PetscFunctionReturn(PETSC_SUCCESS);
     46: }

     48: static PetscErrorCode PCPatchConstruct_Vanka(void *vpatch, DM dm, PetscInt point, PetscHSetI ht)
     49: {
     50:   PC_PATCH *patch = (PC_PATCH *)vpatch;
     51:   PetscInt  starSize;
     52:   PetscInt *star         = NULL;
     53:   PetscBool shouldIgnore = PETSC_FALSE;
     54:   PetscInt  cStart, cEnd, iStart, iEnd, si;

     56:   PetscFunctionBegin;
     57:   PetscCall(PetscHSetIClear(ht));
     58:   /* To start with, add the point we care about */
     59:   PetscCall(PetscHSetIAdd(ht, point));
     60:   /* Should we ignore any points of a certain dimension? */
     61:   if (patch->vankadim >= 0) {
     62:     shouldIgnore = PETSC_TRUE;
     63:     PetscCall(DMPlexGetDepthStratum(dm, patch->vankadim, &iStart, &iEnd));
     64:   }
     65:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
     66:   /* Loop over all the cells that this point connects to */
     67:   PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star));
     68:   for (si = 0; si < starSize * 2; si += 2) {
     69:     const PetscInt cell = star[si];
     70:     PetscInt       closureSize;
     71:     PetscInt      *closure = NULL, ci;

     73:     if (cell < cStart || cell >= cEnd) continue;
     74:     /* now loop over all entities in the closure of that cell */
     75:     PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
     76:     for (ci = 0; ci < closureSize * 2; ci += 2) {
     77:       const PetscInt newpoint = closure[ci];

     79:       /* We've been told to ignore entities of this type.*/
     80:       if (shouldIgnore && newpoint >= iStart && newpoint < iEnd) continue;
     81:       PetscCall(PetscHSetIAdd(ht, newpoint));
     82:     }
     83:     PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
     84:   }
     85:   PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star));
     86:   PetscFunctionReturn(PETSC_SUCCESS);
     87: }

     89: static PetscErrorCode PCPatchConstruct_Pardecomp(void *vpatch, DM dm, PetscInt point, PetscHSetI ht)
     90: {
     91:   PC_PATCH       *patch   = (PC_PATCH *)vpatch;
     92:   DMLabel         ghost   = NULL;
     93:   const PetscInt *leaves  = NULL;
     94:   PetscInt        nleaves = 0, pStart, pEnd, loc;
     95:   PetscBool       isFiredrake;
     96:   PetscBool       flg;
     97:   PetscInt        starSize;
     98:   PetscInt       *star = NULL;
     99:   PetscInt        opoint, overlapi;

    101:   PetscFunctionBegin;
    102:   PetscCall(PetscHSetIClear(ht));

    104:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));

    106:   PetscCall(DMHasLabel(dm, "pyop2_ghost", &isFiredrake));
    107:   if (isFiredrake) {
    108:     PetscCall(DMGetLabel(dm, "pyop2_ghost", &ghost));
    109:     PetscCall(DMLabelCreateIndex(ghost, pStart, pEnd));
    110:   } else {
    111:     PetscSF sf;
    112:     PetscCall(DMGetPointSF(dm, &sf));
    113:     PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL));
    114:     nleaves = PetscMax(nleaves, 0);
    115:   }

    117:   for (opoint = pStart; opoint < pEnd; ++opoint) {
    118:     if (ghost) PetscCall(DMLabelHasPoint(ghost, opoint, &flg));
    119:     else {
    120:       PetscCall(PetscFindInt(opoint, nleaves, leaves, &loc));
    121:       flg = loc >= 0 ? PETSC_TRUE : PETSC_FALSE;
    122:     }
    123:     /* Not an owned entity, don't make a cell patch. */
    124:     if (flg) continue;
    125:     PetscCall(PetscHSetIAdd(ht, opoint));
    126:   }

    128:   /* Now build the overlap for the patch */
    129:   for (overlapi = 0; overlapi < patch->pardecomp_overlap; ++overlapi) {
    130:     PetscInt  index    = 0;
    131:     PetscInt *htpoints = NULL;
    132:     PetscInt  htsize;
    133:     PetscInt  i;

    135:     PetscCall(PetscHSetIGetSize(ht, &htsize));
    136:     PetscCall(PetscMalloc1(htsize, &htpoints));
    137:     PetscCall(PetscHSetIGetElems(ht, &index, htpoints));

    139:     for (i = 0; i < htsize; ++i) {
    140:       PetscInt hpoint = htpoints[i];
    141:       PetscInt si;

    143:       PetscCall(DMPlexGetTransitiveClosure(dm, hpoint, PETSC_FALSE, &starSize, &star));
    144:       for (si = 0; si < starSize * 2; si += 2) {
    145:         const PetscInt starp = star[si];
    146:         PetscInt       closureSize;
    147:         PetscInt      *closure = NULL, ci;

    149:         /* now loop over all entities in the closure of starp */
    150:         PetscCall(DMPlexGetTransitiveClosure(dm, starp, PETSC_TRUE, &closureSize, &closure));
    151:         for (ci = 0; ci < closureSize * 2; ci += 2) {
    152:           const PetscInt closstarp = closure[ci];
    153:           PetscCall(PetscHSetIAdd(ht, closstarp));
    154:         }
    155:         PetscCall(DMPlexRestoreTransitiveClosure(dm, starp, PETSC_TRUE, &closureSize, &closure));
    156:       }
    157:       PetscCall(DMPlexRestoreTransitiveClosure(dm, hpoint, PETSC_FALSE, &starSize, &star));
    158:     }
    159:     PetscCall(PetscFree(htpoints));
    160:   }
    161:   PetscFunctionReturn(PETSC_SUCCESS);
    162: }

    164: /* The user's already set the patches in patch->userIS. Build the hash tables */
    165: static PetscErrorCode PCPatchConstruct_User(void *vpatch, DM dm, PetscInt point, PetscHSetI ht)
    166: {
    167:   PC_PATCH       *patch   = (PC_PATCH *)vpatch;
    168:   IS              patchis = patch->userIS[point];
    169:   PetscInt        n;
    170:   const PetscInt *patchdata;
    171:   PetscInt        pStart, pEnd, i;

    173:   PetscFunctionBegin;
    174:   PetscCall(PetscHSetIClear(ht));
    175:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    176:   PetscCall(ISGetLocalSize(patchis, &n));
    177:   PetscCall(ISGetIndices(patchis, &patchdata));
    178:   for (i = 0; i < n; ++i) {
    179:     const PetscInt ownedpoint = patchdata[i];

    181:     PetscCheck(ownedpoint >= pStart && ownedpoint < pEnd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %" PetscInt_FMT " was not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", ownedpoint, pStart, pEnd);
    182:     PetscCall(PetscHSetIAdd(ht, ownedpoint));
    183:   }
    184:   PetscCall(ISRestoreIndices(patchis, &patchdata));
    185:   PetscFunctionReturn(PETSC_SUCCESS);
    186: }

    188: static PetscErrorCode PCPatchCreateDefaultSF_Private(PC pc, PetscInt n, const PetscSF *sf, const PetscInt *bs)
    189: {
    190:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    192:   PetscFunctionBegin;
    193:   if (n == 1 && bs[0] == 1) {
    194:     patch->sectionSF = sf[0];
    195:     PetscCall(PetscObjectReference((PetscObject)patch->sectionSF));
    196:   } else {
    197:     PetscInt     allRoots = 0, allLeaves = 0;
    198:     PetscInt     leafOffset    = 0;
    199:     PetscInt    *ilocal        = NULL;
    200:     PetscSFNode *iremote       = NULL;
    201:     PetscInt    *remoteOffsets = NULL;
    202:     PetscInt     index         = 0;
    203:     PetscHMapI   rankToIndex;
    204:     PetscInt     numRanks = 0;
    205:     PetscSFNode *remote   = NULL;
    206:     PetscSF      rankSF;
    207:     PetscInt    *ranks   = NULL;
    208:     PetscInt    *offsets = NULL;
    209:     MPI_Datatype contig;
    210:     PetscHSetI   ranksUniq;
    211:     PetscMPIInt  in;

    213:     /* First figure out how many dofs there are in the concatenated numbering.
    214:        allRoots: number of owned global dofs;
    215:        allLeaves: number of visible dofs (global + ghosted).
    216:     */
    217:     for (PetscInt i = 0; i < n; ++i) {
    218:       PetscInt nroots, nleaves;

    220:       PetscCall(PetscSFGetGraph(sf[i], &nroots, &nleaves, NULL, NULL));
    221:       allRoots += nroots * bs[i];
    222:       allLeaves += nleaves * bs[i];
    223:     }
    224:     PetscCall(PetscMalloc1(allLeaves, &ilocal));
    225:     PetscCall(PetscMalloc1(allLeaves, &iremote));
    226:     /* Now build an SF that just contains process connectivity. */
    227:     PetscCall(PetscHSetICreate(&ranksUniq));
    228:     for (PetscInt i = 0; i < n; ++i) {
    229:       const PetscMPIInt *ranks = NULL;
    230:       PetscMPIInt        nranks;

    232:       PetscCall(PetscSFSetUp(sf[i]));
    233:       PetscCall(PetscSFGetRootRanks(sf[i], &nranks, &ranks, NULL, NULL, NULL));
    234:       /* These are all the ranks who communicate with me. */
    235:       for (PetscMPIInt j = 0; j < nranks; ++j) PetscCall(PetscHSetIAdd(ranksUniq, (PetscInt)ranks[j]));
    236:     }
    237:     PetscCall(PetscHSetIGetSize(ranksUniq, &numRanks));
    238:     PetscCall(PetscMalloc1(numRanks, &remote));
    239:     PetscCall(PetscMalloc1(numRanks, &ranks));
    240:     PetscCall(PetscHSetIGetElems(ranksUniq, &index, ranks));

    242:     PetscCall(PetscHMapICreate(&rankToIndex));
    243:     for (PetscInt i = 0; i < numRanks; ++i) {
    244:       remote[i].rank  = ranks[i];
    245:       remote[i].index = 0;
    246:       PetscCall(PetscHMapISet(rankToIndex, ranks[i], i));
    247:     }
    248:     PetscCall(PetscFree(ranks));
    249:     PetscCall(PetscHSetIDestroy(&ranksUniq));
    250:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)pc), &rankSF));
    251:     PetscCall(PetscSFSetGraph(rankSF, 1, numRanks, NULL, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER));
    252:     PetscCall(PetscSFSetUp(rankSF));
    253:     /* OK, use it to communicate the root offset on the remote processes for each subspace. */
    254:     PetscCall(PetscMalloc1(n, &offsets));
    255:     PetscCall(PetscMalloc1(n * numRanks, &remoteOffsets));

    257:     offsets[0] = 0;
    258:     for (PetscInt i = 1; i < n; ++i) {
    259:       PetscInt nroots;

    261:       PetscCall(PetscSFGetGraph(sf[i - 1], &nroots, NULL, NULL, NULL));
    262:       offsets[i] = offsets[i - 1] + nroots * bs[i - 1];
    263:     }
    264:     /* Offsets are the offsets on the current process of the global dof numbering for the subspaces. */
    265:     PetscCall(PetscMPIIntCast(n, &in));
    266:     PetscCallMPI(MPI_Type_contiguous(in, MPIU_INT, &contig));
    267:     PetscCallMPI(MPI_Type_commit(&contig));

    269:     PetscCall(PetscSFBcastBegin(rankSF, contig, offsets, remoteOffsets, MPI_REPLACE));
    270:     PetscCall(PetscSFBcastEnd(rankSF, contig, offsets, remoteOffsets, MPI_REPLACE));
    271:     PetscCallMPI(MPI_Type_free(&contig));
    272:     PetscCall(PetscFree(offsets));
    273:     PetscCall(PetscSFDestroy(&rankSF));
    274:     /* Now remoteOffsets contains the offsets on the remote
    275:       processes who communicate with me.  So now we can
    276:       concatenate the list of SFs into a single one. */
    277:     index = 0;
    278:     for (PetscInt i = 0; i < n; ++i) {
    279:       const PetscSFNode *remote = NULL;
    280:       const PetscInt    *local  = NULL;
    281:       PetscInt           nroots, nleaves, j;

    283:       PetscCall(PetscSFGetGraph(sf[i], &nroots, &nleaves, &local, &remote));
    284:       for (j = 0; j < nleaves; ++j) {
    285:         PetscInt rank = remote[j].rank;
    286:         PetscInt idx, rootOffset, k;

    288:         PetscCall(PetscHMapIGet(rankToIndex, rank, &idx));
    289:         PetscCheck(idx != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Didn't find rank, huh?");
    290:         /* Offset on given rank for ith subspace */
    291:         rootOffset = remoteOffsets[n * idx + i];
    292:         for (k = 0; k < bs[i]; ++k) {
    293:           ilocal[index]        = (local ? local[j] : j) * bs[i] + k + leafOffset;
    294:           iremote[index].rank  = remote[j].rank;
    295:           iremote[index].index = remote[j].index * bs[i] + k + rootOffset;
    296:           ++index;
    297:         }
    298:       }
    299:       leafOffset += nleaves * bs[i];
    300:     }
    301:     PetscCall(PetscHMapIDestroy(&rankToIndex));
    302:     PetscCall(PetscFree(remoteOffsets));
    303:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)pc), &patch->sectionSF));
    304:     PetscCall(PetscSFSetGraph(patch->sectionSF, allRoots, allLeaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
    305:   }
    306:   PetscFunctionReturn(PETSC_SUCCESS);
    307: }

    309: /* TODO: Docs */
    310: static PetscErrorCode PCPatchGetIgnoreDim(PC pc, PetscInt *dim)
    311: {
    312:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    314:   PetscFunctionBegin;
    315:   *dim = patch->ignoredim;
    316:   PetscFunctionReturn(PETSC_SUCCESS);
    317: }

    319: /* TODO: Docs */
    320: PetscErrorCode PCPatchSetSaveOperators(PC pc, PetscBool flg)
    321: {
    322:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    324:   PetscFunctionBegin;
    325:   patch->save_operators = flg;
    326:   PetscFunctionReturn(PETSC_SUCCESS);
    327: }

    329: /* TODO: Docs */
    330: PetscErrorCode PCPatchGetSaveOperators(PC pc, PetscBool *flg)
    331: {
    332:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    334:   PetscFunctionBegin;
    335:   *flg = patch->save_operators;
    336:   PetscFunctionReturn(PETSC_SUCCESS);
    337: }

    339: /* TODO: Docs */
    340: PetscErrorCode PCPatchSetPrecomputeElementTensors(PC pc, PetscBool flg)
    341: {
    342:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    344:   PetscFunctionBegin;
    345:   patch->precomputeElementTensors = flg;
    346:   PetscFunctionReturn(PETSC_SUCCESS);
    347: }

    349: /* TODO: Docs */
    350: PetscErrorCode PCPatchGetPrecomputeElementTensors(PC pc, PetscBool *flg)
    351: {
    352:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    354:   PetscFunctionBegin;
    355:   *flg = patch->precomputeElementTensors;
    356:   PetscFunctionReturn(PETSC_SUCCESS);
    357: }

    359: /* TODO: Docs */
    360: PetscErrorCode PCPatchSetPartitionOfUnity(PC pc, PetscBool flg)
    361: {
    362:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    364:   PetscFunctionBegin;
    365:   patch->partition_of_unity = flg;
    366:   PetscFunctionReturn(PETSC_SUCCESS);
    367: }

    369: /* TODO: Docs */
    370: PetscErrorCode PCPatchGetPartitionOfUnity(PC pc, PetscBool *flg)
    371: {
    372:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    374:   PetscFunctionBegin;
    375:   *flg = patch->partition_of_unity;
    376:   PetscFunctionReturn(PETSC_SUCCESS);
    377: }

    379: /* TODO: Docs */
    380: static PetscErrorCode PCPatchSetLocalComposition(PC pc, PCCompositeType type)
    381: {
    382:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    384:   PetscFunctionBegin;
    385:   PetscCheck(type == PC_COMPOSITE_ADDITIVE || type == PC_COMPOSITE_MULTIPLICATIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Only supports additive or multiplicative as the local type");
    386:   patch->local_composition_type = type;
    387:   PetscFunctionReturn(PETSC_SUCCESS);
    388: }

    390: /* TODO: Docs */
    391: PetscErrorCode PCPatchGetSubKSP(PC pc, PetscInt *npatch, KSP **ksp)
    392: {
    393:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    395:   PetscFunctionBegin;
    396:   PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Need to call PCSetUp() on PC (or KSPSetUp() on the outer KSP object) before calling here");
    397:   PetscCall(PetscMalloc1(patch->npatch, ksp));
    398:   for (PetscInt i = 0; i < patch->npatch; ++i) (*ksp)[i] = (KSP)patch->solver[i];
    399:   if (npatch) *npatch = patch->npatch;
    400:   PetscFunctionReturn(PETSC_SUCCESS);
    401: }

    403: /* TODO: Docs */
    404: PetscErrorCode PCPatchSetSubMatType(PC pc, MatType sub_mat_type)
    405: {
    406:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    408:   PetscFunctionBegin;
    409:   if (patch->sub_mat_type) PetscCall(PetscFree(patch->sub_mat_type));
    410:   PetscCall(PetscStrallocpy(sub_mat_type, (char **)&patch->sub_mat_type));
    411:   PetscFunctionReturn(PETSC_SUCCESS);
    412: }

    414: /* TODO: Docs */
    415: PetscErrorCode PCPatchGetSubMatType(PC pc, MatType *sub_mat_type)
    416: {
    417:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    419:   PetscFunctionBegin;
    420:   *sub_mat_type = patch->sub_mat_type;
    421:   PetscFunctionReturn(PETSC_SUCCESS);
    422: }

    424: /* TODO: Docs */
    425: PetscErrorCode PCPatchSetCellNumbering(PC pc, PetscSection cellNumbering)
    426: {
    427:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    429:   PetscFunctionBegin;
    430:   patch->cellNumbering = cellNumbering;
    431:   PetscCall(PetscObjectReference((PetscObject)cellNumbering));
    432:   PetscFunctionReturn(PETSC_SUCCESS);
    433: }

    435: /* TODO: Docs */
    436: PetscErrorCode PCPatchGetCellNumbering(PC pc, PetscSection *cellNumbering)
    437: {
    438:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    440:   PetscFunctionBegin;
    441:   *cellNumbering = patch->cellNumbering;
    442:   PetscFunctionReturn(PETSC_SUCCESS);
    443: }

    445: /* TODO: Docs */
    446: PetscErrorCode PCPatchSetConstructType(PC pc, PCPatchConstructType ctype, PetscErrorCode (*func)(PC, PetscInt *, IS **, IS *, void *), void *ctx)
    447: {
    448:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    450:   PetscFunctionBegin;
    451:   patch->ctype = ctype;
    452:   switch (ctype) {
    453:   case PC_PATCH_STAR:
    454:     patch->user_patches     = PETSC_FALSE;
    455:     patch->patchconstructop = PCPatchConstruct_Star;
    456:     break;
    457:   case PC_PATCH_VANKA:
    458:     patch->user_patches     = PETSC_FALSE;
    459:     patch->patchconstructop = PCPatchConstruct_Vanka;
    460:     break;
    461:   case PC_PATCH_PARDECOMP:
    462:     patch->user_patches     = PETSC_FALSE;
    463:     patch->patchconstructop = PCPatchConstruct_Pardecomp;
    464:     break;
    465:   case PC_PATCH_USER:
    466:   case PC_PATCH_PYTHON:
    467:     patch->user_patches     = PETSC_TRUE;
    468:     patch->patchconstructop = PCPatchConstruct_User;
    469:     if (func) {
    470:       patch->userpatchconstructionop = func;
    471:       patch->userpatchconstructctx   = ctx;
    472:     }
    473:     break;
    474:   default:
    475:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "Unknown patch construction type %" PetscInt_FMT, (PetscInt)patch->ctype);
    476:   }
    477:   PetscFunctionReturn(PETSC_SUCCESS);
    478: }

    480: /* TODO: Docs */
    481: PetscErrorCode PCPatchGetConstructType(PC pc, PCPatchConstructType *ctype, PetscErrorCode (**func)(PC, PetscInt *, IS **, IS *, void *), void **ctx)
    482: {
    483:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    485:   PetscFunctionBegin;
    486:   *ctype = patch->ctype;
    487:   switch (patch->ctype) {
    488:   case PC_PATCH_STAR:
    489:   case PC_PATCH_VANKA:
    490:   case PC_PATCH_PARDECOMP:
    491:     break;
    492:   case PC_PATCH_USER:
    493:   case PC_PATCH_PYTHON:
    494:     *func = patch->userpatchconstructionop;
    495:     *ctx  = patch->userpatchconstructctx;
    496:     break;
    497:   default:
    498:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "Unknown patch construction type %" PetscInt_FMT, (PetscInt)patch->ctype);
    499:   }
    500:   PetscFunctionReturn(PETSC_SUCCESS);
    501: }

    503: /* TODO: Docs */
    504: PetscErrorCode PCPatchSetDiscretisationInfo(PC pc, PetscInt nsubspaces, DM *dms, PetscInt *bs, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, const PetscInt *subspaceOffsets, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes)
    505: {
    506:   PC_PATCH *patch = (PC_PATCH *)pc->data;
    507:   DM        dm, plex;
    508:   PetscSF  *sfs;
    509:   PetscInt  cStart, cEnd, i, j;

    511:   PetscFunctionBegin;
    512:   PetscCall(PCGetDM(pc, &dm));
    513:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    514:   dm = plex;
    515:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
    516:   PetscCall(PetscMalloc1(nsubspaces, &sfs));
    517:   PetscCall(PetscMalloc1(nsubspaces, &patch->dofSection));
    518:   PetscCall(PetscMalloc1(nsubspaces, &patch->bs));
    519:   PetscCall(PetscMalloc1(nsubspaces, &patch->nodesPerCell));
    520:   PetscCall(PetscMalloc1(nsubspaces, &patch->cellNodeMap));
    521:   PetscCall(PetscMalloc1(nsubspaces + 1, &patch->subspaceOffsets));

    523:   patch->nsubspaces       = nsubspaces;
    524:   patch->totalDofsPerCell = 0;
    525:   for (i = 0; i < nsubspaces; ++i) {
    526:     PetscCall(DMGetLocalSection(dms[i], &patch->dofSection[i]));
    527:     PetscCall(PetscObjectReference((PetscObject)patch->dofSection[i]));
    528:     PetscCall(DMGetSectionSF(dms[i], &sfs[i]));
    529:     patch->bs[i]           = bs[i];
    530:     patch->nodesPerCell[i] = nodesPerCell[i];
    531:     patch->totalDofsPerCell += nodesPerCell[i] * bs[i];
    532:     PetscCall(PetscMalloc1((cEnd - cStart) * nodesPerCell[i], &patch->cellNodeMap[i]));
    533:     for (j = 0; j < (cEnd - cStart) * nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j];
    534:     patch->subspaceOffsets[i] = subspaceOffsets[i];
    535:   }
    536:   PetscCall(PCPatchCreateDefaultSF_Private(pc, nsubspaces, sfs, patch->bs));
    537:   PetscCall(PetscFree(sfs));

    539:   patch->subspaceOffsets[nsubspaces] = subspaceOffsets[nsubspaces];
    540:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes));
    541:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes));
    542:   PetscCall(DMDestroy(&dm));
    543:   PetscFunctionReturn(PETSC_SUCCESS);
    544: }

    546: /* TODO: Docs */
    547: static PetscErrorCode PCPatchSetDiscretisationInfoCombined(PC pc, DM dm, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes)
    548: {
    549:   PC_PATCH *patch = (PC_PATCH *)pc->data;
    550:   PetscInt  cStart, cEnd, i, j;

    552:   PetscFunctionBegin;
    553:   patch->combined = PETSC_TRUE;
    554:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
    555:   PetscCall(DMGetNumFields(dm, &patch->nsubspaces));
    556:   PetscCall(PetscCalloc1(patch->nsubspaces, &patch->dofSection));
    557:   PetscCall(PetscMalloc1(patch->nsubspaces, &patch->bs));
    558:   PetscCall(PetscMalloc1(patch->nsubspaces, &patch->nodesPerCell));
    559:   PetscCall(PetscMalloc1(patch->nsubspaces, &patch->cellNodeMap));
    560:   PetscCall(PetscCalloc1(patch->nsubspaces + 1, &patch->subspaceOffsets));
    561:   PetscCall(DMGetLocalSection(dm, &patch->dofSection[0]));
    562:   PetscCall(PetscObjectReference((PetscObject)patch->dofSection[0]));
    563:   PetscCall(PetscSectionGetStorageSize(patch->dofSection[0], &patch->subspaceOffsets[patch->nsubspaces]));
    564:   patch->totalDofsPerCell = 0;
    565:   for (i = 0; i < patch->nsubspaces; ++i) {
    566:     patch->bs[i]           = 1;
    567:     patch->nodesPerCell[i] = nodesPerCell[i];
    568:     patch->totalDofsPerCell += nodesPerCell[i];
    569:     PetscCall(PetscMalloc1((cEnd - cStart) * nodesPerCell[i], &patch->cellNodeMap[i]));
    570:     for (j = 0; j < (cEnd - cStart) * nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j];
    571:   }
    572:   PetscCall(DMGetSectionSF(dm, &patch->sectionSF));
    573:   PetscCall(PetscObjectReference((PetscObject)patch->sectionSF));
    574:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes));
    575:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes));
    576:   PetscFunctionReturn(PETSC_SUCCESS);
    577: }

    579: /*@C
    580:   PCPatchSetComputeFunction - Set the callback function used to compute patch residuals

    582:   Logically Collective

    584:   Input Parameters:
    585: + pc   - The `PC`
    586: . func - The callback function
    587: - ctx  - The user context

    589:   Calling sequence of `func`:
    590: + pc               - The `PC`
    591: . point            - The point
    592: . x                - The input solution (not used in linear problems)
    593: . f                - The patch residual vector
    594: . cellIS           - An array of the cell numbers
    595: . n                - The size of `dofsArray`
    596: . dofsArray        - The dofmap for the dofs to be solved for
    597: . dofsArrayWithAll - The dofmap for all dofs on the patch
    598: - ctx              - The user context

    600:   Level: advanced

    602:   Note:
    603:   The entries of `f` (the output residual vector) have been set to zero before the call.

    605: .seealso: [](ch_ksp), `PCPatchSetComputeOperator()`, `PCPatchGetComputeOperator()`, `PCPatchSetDiscretisationInfo()`, `PCPatchSetComputeFunctionInteriorFacets()`
    606: @*/
    607: PetscErrorCode PCPatchSetComputeFunction(PC pc, PetscErrorCode (*func)(PC pc, PetscInt point, Vec x, Vec f, IS cellIS, PetscInt n, const PetscInt *dofsArray, const PetscInt *dofsArrayWithAll, void *ctx), void *ctx)
    608: {
    609:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    611:   PetscFunctionBegin;
    612:   patch->usercomputef    = func;
    613:   patch->usercomputefctx = ctx;
    614:   PetscFunctionReturn(PETSC_SUCCESS);
    615: }

    617: /*@C
    618:   PCPatchSetComputeFunctionInteriorFacets - Set the callback function used to compute facet integrals for patch residuals

    620:   Logically Collective

    622:   Input Parameters:
    623: + pc   - The `PC`
    624: . func - The callback function
    625: - ctx  - The user context

    627:   Calling sequence of `func`:
    628: + pc               - The `PC`
    629: . point            - The point
    630: . x                - The input solution (not used in linear problems)
    631: . f                - The patch residual vector
    632: . facetIS          - An array of the facet numbers
    633: . n                - The size of `dofsArray`
    634: . dofsArray        - The dofmap for the dofs to be solved for
    635: . dofsArrayWithAll - The dofmap for all dofs on the patch
    636: - ctx              - The user context

    638:   Level: advanced

    640:   Note:
    641:   The entries of `f` (the output residual vector) have been set to zero before the call.

    643: .seealso: [](ch_ksp), `PCPatchSetComputeOperator()`, `PCPatchGetComputeOperator()`, `PCPatchSetDiscretisationInfo()`, `PCPatchSetComputeFunction()`
    644: @*/
    645: PetscErrorCode PCPatchSetComputeFunctionInteriorFacets(PC pc, PetscErrorCode (*func)(PC pc, PetscInt point, Vec x, Vec f, IS facetIS, PetscInt n, const PetscInt *dofsArray, const PetscInt *dofsArrayWithAll, void *ctx), void *ctx)
    646: {
    647:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    649:   PetscFunctionBegin;
    650:   patch->usercomputefintfacet    = func;
    651:   patch->usercomputefintfacetctx = ctx;
    652:   PetscFunctionReturn(PETSC_SUCCESS);
    653: }

    655: /*@C
    656:   PCPatchSetComputeOperator - Set the callback function used to compute patch matrices

    658:   Logically Collective

    660:   Input Parameters:
    661: + pc   - The `PC`
    662: . func - The callback function
    663: - ctx  - The user context

    665:   Calling sequence of `func`:
    666: + pc               - The `PC`
    667: . point            - The point
    668: . x                - The input solution (not used in linear problems)
    669: . mat              - The patch matrix
    670: . facetIS          - An array of the cell numbers
    671: . n                - The size of `dofsArray`
    672: . dofsArray        - The dofmap for the dofs to be solved for
    673: . dofsArrayWithAll - The dofmap for all dofs on the patch
    674: - ctx              - The user context

    676:   Level: advanced

    678:   Note:
    679:   The matrix entries have been set to zero before the call.

    681: .seealso: [](ch_ksp), `PCPatchGetComputeOperator()`, `PCPatchSetComputeFunction()`, `PCPatchSetDiscretisationInfo()`
    682: @*/
    683: PetscErrorCode PCPatchSetComputeOperator(PC pc, PetscErrorCode (*func)(PC pc, PetscInt point, Vec x, Mat mat, IS facetIS, PetscInt n, const PetscInt *dofsArray, const PetscInt *dofsArrayWithAll, void *ctx), void *ctx)
    684: {
    685:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    687:   PetscFunctionBegin;
    688:   patch->usercomputeop    = func;
    689:   patch->usercomputeopctx = ctx;
    690:   PetscFunctionReturn(PETSC_SUCCESS);
    691: }

    693: /*@C
    694:   PCPatchSetComputeOperatorInteriorFacets - Set the callback function used to compute facet integrals for patch matrices

    696:   Logically Collective

    698:   Input Parameters:
    699: + pc   - The `PC`
    700: . func - The callback function
    701: - ctx  - The user context

    703:   Calling sequence of `func`:
    704: + pc               - The `PC`
    705: . point            - The point
    706: . x                - The input solution (not used in linear problems)
    707: . mat              - The patch matrix
    708: . facetIS          - An array of the facet numbers
    709: . n                - The size of `dofsArray`
    710: . dofsArray        - The dofmap for the dofs to be solved for
    711: . dofsArrayWithAll - The dofmap for all dofs on the patch
    712: - ctx              - The user context

    714:   Level: advanced

    716:   Note:
    717:   The matrix entries have been set to zero before the call.

    719: .seealso: [](ch_ksp), `PCPatchGetComputeOperator()`, `PCPatchSetComputeFunction()`, `PCPatchSetDiscretisationInfo()`
    720: @*/
    721: PetscErrorCode PCPatchSetComputeOperatorInteriorFacets(PC pc, PetscErrorCode (*func)(PC pc, PetscInt point, Vec x, Mat mat, IS facetIS, PetscInt n, const PetscInt *dofsArray, const PetscInt *dofsArrayWithAll, void *ctx), void *ctx)
    722: {
    723:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    725:   PetscFunctionBegin;
    726:   patch->usercomputeopintfacet    = func;
    727:   patch->usercomputeopintfacetctx = ctx;
    728:   PetscFunctionReturn(PETSC_SUCCESS);
    729: }

    731: /* On entry, ht contains the topological entities whose dofs we are responsible for solving for;
    732:    on exit, cht contains all the topological entities we need to compute their residuals.
    733:    In full generality this should incorporate knowledge of the sparsity pattern of the matrix;
    734:    here we assume a standard FE sparsity pattern.*/
    735: /* TODO: Use DMPlexGetAdjacency() */
    736: static PetscErrorCode PCPatchCompleteCellPatch(PC pc, PetscHSetI ht, PetscHSetI cht)
    737: {
    738:   DM            dm, plex;
    739:   PC_PATCH     *patch = (PC_PATCH *)pc->data;
    740:   PetscHashIter hi;
    741:   PetscInt      point;
    742:   PetscInt     *star = NULL, *closure = NULL;
    743:   PetscInt      ignoredim, iStart = 0, iEnd = -1, starSize, closureSize, si, ci;
    744:   PetscInt     *fStar = NULL, *fClosure = NULL;
    745:   PetscInt      fBegin, fEnd, fsi, fci, fStarSize, fClosureSize;

    747:   PetscFunctionBegin;
    748:   PetscCall(PCGetDM(pc, &dm));
    749:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    750:   dm = plex;
    751:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fBegin, &fEnd));
    752:   PetscCall(PCPatchGetIgnoreDim(pc, &ignoredim));
    753:   if (ignoredim >= 0) PetscCall(DMPlexGetDepthStratum(dm, ignoredim, &iStart, &iEnd));
    754:   PetscCall(PetscHSetIClear(cht));
    755:   PetscHashIterBegin(ht, hi);
    756:   while (!PetscHashIterAtEnd(ht, hi)) {
    757:     PetscHashIterGetKey(ht, hi, point);
    758:     PetscHashIterNext(ht, hi);

    760:     /* Loop over all the cells that this point connects to */
    761:     PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star));
    762:     for (si = 0; si < starSize * 2; si += 2) {
    763:       const PetscInt ownedpoint = star[si];
    764:       /* TODO Check for point in cht before running through closure again */
    765:       /* now loop over all entities in the closure of that cell */
    766:       PetscCall(DMPlexGetTransitiveClosure(dm, ownedpoint, PETSC_TRUE, &closureSize, &closure));
    767:       for (ci = 0; ci < closureSize * 2; ci += 2) {
    768:         const PetscInt seenpoint = closure[ci];
    769:         if (ignoredim >= 0 && seenpoint >= iStart && seenpoint < iEnd) continue;
    770:         PetscCall(PetscHSetIAdd(cht, seenpoint));
    771:         /* Facet integrals couple dofs across facets, so in that case for each of
    772:           the facets we need to add all dofs on the other side of the facet to
    773:           the seen dofs. */
    774:         if (patch->usercomputeopintfacet) {
    775:           if (fBegin <= seenpoint && seenpoint < fEnd) {
    776:             PetscCall(DMPlexGetTransitiveClosure(dm, seenpoint, PETSC_FALSE, &fStarSize, &fStar));
    777:             for (fsi = 0; fsi < fStarSize * 2; fsi += 2) {
    778:               PetscCall(DMPlexGetTransitiveClosure(dm, fStar[fsi], PETSC_TRUE, &fClosureSize, &fClosure));
    779:               for (fci = 0; fci < fClosureSize * 2; fci += 2) PetscCall(PetscHSetIAdd(cht, fClosure[fci]));
    780:               PetscCall(DMPlexRestoreTransitiveClosure(dm, fStar[fsi], PETSC_TRUE, NULL, &fClosure));
    781:             }
    782:             PetscCall(DMPlexRestoreTransitiveClosure(dm, seenpoint, PETSC_FALSE, NULL, &fStar));
    783:           }
    784:         }
    785:       }
    786:       PetscCall(DMPlexRestoreTransitiveClosure(dm, ownedpoint, PETSC_TRUE, NULL, &closure));
    787:     }
    788:     PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, NULL, &star));
    789:   }
    790:   PetscCall(DMDestroy(&dm));
    791:   PetscFunctionReturn(PETSC_SUCCESS);
    792: }

    794: static PetscErrorCode PCPatchGetGlobalDofs(PC pc, PetscSection dofSection[], PetscInt f, PetscBool combined, PetscInt p, PetscInt *dof, PetscInt *off)
    795: {
    796:   PetscFunctionBegin;
    797:   if (combined) {
    798:     if (f < 0) {
    799:       if (dof) PetscCall(PetscSectionGetDof(dofSection[0], p, dof));
    800:       if (off) PetscCall(PetscSectionGetOffset(dofSection[0], p, off));
    801:     } else {
    802:       if (dof) PetscCall(PetscSectionGetFieldDof(dofSection[0], p, f, dof));
    803:       if (off) PetscCall(PetscSectionGetFieldOffset(dofSection[0], p, f, off));
    804:     }
    805:   } else {
    806:     if (f < 0) {
    807:       PC_PATCH *patch = (PC_PATCH *)pc->data;
    808:       PetscInt  fdof, g;

    810:       if (dof) {
    811:         *dof = 0;
    812:         for (g = 0; g < patch->nsubspaces; ++g) {
    813:           PetscCall(PetscSectionGetDof(dofSection[g], p, &fdof));
    814:           *dof += fdof;
    815:         }
    816:       }
    817:       if (off) {
    818:         *off = 0;
    819:         for (g = 0; g < patch->nsubspaces; ++g) {
    820:           PetscCall(PetscSectionGetOffset(dofSection[g], p, &fdof));
    821:           *off += fdof;
    822:         }
    823:       }
    824:     } else {
    825:       if (dof) PetscCall(PetscSectionGetDof(dofSection[f], p, dof));
    826:       if (off) PetscCall(PetscSectionGetOffset(dofSection[f], p, off));
    827:     }
    828:   }
    829:   PetscFunctionReturn(PETSC_SUCCESS);
    830: }

    832: /* Given a hash table with a set of topological entities (pts), compute the degrees of
    833:    freedom in global concatenated numbering on those entities.
    834:    For Vanka smoothing, this needs to do something special: ignore dofs of the
    835:    constraint subspace on entities that aren't the base entity we're building the patch
    836:    around. */
    837: static PetscErrorCode PCPatchGetPointDofs(PC pc, PetscHSetI pts, PetscHSetI dofs, PetscInt base, PetscHSetI *subspaces_to_exclude)
    838: {
    839:   PC_PATCH     *patch = (PC_PATCH *)pc->data;
    840:   PetscHashIter hi;
    841:   PetscInt      ldof, loff;
    842:   PetscInt      k, p;

    844:   PetscFunctionBegin;
    845:   PetscCall(PetscHSetIClear(dofs));
    846:   for (k = 0; k < patch->nsubspaces; ++k) {
    847:     PetscInt subspaceOffset = patch->subspaceOffsets[k];
    848:     PetscInt bs             = patch->bs[k];
    849:     PetscInt j, l;

    851:     if (subspaces_to_exclude != NULL) {
    852:       PetscBool should_exclude_k = PETSC_FALSE;
    853:       PetscCall(PetscHSetIHas(*subspaces_to_exclude, k, &should_exclude_k));
    854:       if (should_exclude_k) {
    855:         /* only get this subspace dofs at the base entity, not any others */
    856:         PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, base, &ldof, &loff));
    857:         if (0 == ldof) continue;
    858:         for (j = loff; j < ldof + loff; ++j) {
    859:           for (l = 0; l < bs; ++l) {
    860:             PetscInt dof = bs * j + l + subspaceOffset;
    861:             PetscCall(PetscHSetIAdd(dofs, dof));
    862:           }
    863:         }
    864:         continue; /* skip the other dofs of this subspace */
    865:       }
    866:     }

    868:     PetscHashIterBegin(pts, hi);
    869:     while (!PetscHashIterAtEnd(pts, hi)) {
    870:       PetscHashIterGetKey(pts, hi, p);
    871:       PetscHashIterNext(pts, hi);
    872:       PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, p, &ldof, &loff));
    873:       if (0 == ldof) continue;
    874:       for (j = loff; j < ldof + loff; ++j) {
    875:         for (l = 0; l < bs; ++l) {
    876:           PetscInt dof = bs * j + l + subspaceOffset;
    877:           PetscCall(PetscHSetIAdd(dofs, dof));
    878:         }
    879:       }
    880:     }
    881:   }
    882:   PetscFunctionReturn(PETSC_SUCCESS);
    883: }

    885: /* Given two hash tables A and B, compute the keys in B that are not in A, and put them in C */
    886: static PetscErrorCode PCPatchComputeSetDifference_Private(PetscHSetI A, PetscHSetI B, PetscHSetI C)
    887: {
    888:   PetscHashIter hi;
    889:   PetscInt      key;
    890:   PetscBool     flg;

    892:   PetscFunctionBegin;
    893:   PetscCall(PetscHSetIClear(C));
    894:   PetscHashIterBegin(B, hi);
    895:   while (!PetscHashIterAtEnd(B, hi)) {
    896:     PetscHashIterGetKey(B, hi, key);
    897:     PetscHashIterNext(B, hi);
    898:     PetscCall(PetscHSetIHas(A, key, &flg));
    899:     if (!flg) PetscCall(PetscHSetIAdd(C, key));
    900:   }
    901:   PetscFunctionReturn(PETSC_SUCCESS);
    902: }

    904: // PetscClangLinter pragma disable: -fdoc-sowing-chars
    905: /*
    906:   PCPatchCreateCellPatches - create patches.

    908:   Input Parameter:
    909:   . dm - The DMPlex object defining the mesh

    911:   Output Parameters:
    912:   + cellCounts  - Section with counts of cells around each vertex
    913:   . cells       - IS of the cell point indices of cells in each patch
    914:   . pointCounts - Section with counts of cells around each vertex
    915:   - point       - IS of the cell point indices of cells in each patch
    916:  */
    917: static PetscErrorCode PCPatchCreateCellPatches(PC pc)
    918: {
    919:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
    920:   DMLabel         ghost = NULL;
    921:   DM              dm, plex;
    922:   PetscHSetI      ht = NULL, cht = NULL;
    923:   PetscSection    cellCounts, pointCounts, intFacetCounts, extFacetCounts;
    924:   PetscInt       *cellsArray, *pointsArray, *intFacetsArray, *extFacetsArray, *intFacetsToPatchCell;
    925:   PetscInt        numCells, numPoints, numIntFacets, numExtFacets;
    926:   const PetscInt *leaves;
    927:   PetscInt        nleaves, pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, v;
    928:   PetscBool       isFiredrake;

    930:   PetscFunctionBegin;
    931:   /* Used to keep track of the cells in the patch. */
    932:   PetscCall(PetscHSetICreate(&ht));
    933:   PetscCall(PetscHSetICreate(&cht));

    935:   PetscCall(PCGetDM(pc, &dm));
    936:   PetscCheck(dm, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "DM not yet set on patch PC");
    937:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    938:   dm = plex;
    939:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    940:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));

    942:   if (patch->user_patches) {
    943:     PetscCall(patch->userpatchconstructionop(pc, &patch->npatch, &patch->userIS, &patch->iterationSet, patch->userpatchconstructctx));
    944:     vStart = 0;
    945:     vEnd   = patch->npatch;
    946:   } else if (patch->ctype == PC_PATCH_PARDECOMP) {
    947:     vStart = 0;
    948:     vEnd   = 1;
    949:   } else if (patch->codim < 0) {
    950:     if (patch->dim < 0) PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
    951:     else PetscCall(DMPlexGetDepthStratum(dm, patch->dim, &vStart, &vEnd));
    952:   } else PetscCall(DMPlexGetHeightStratum(dm, patch->codim, &vStart, &vEnd));
    953:   patch->npatch = vEnd - vStart;

    955:   /* These labels mark the owned points.  We only create patches around points that this process owns. */
    956:   PetscCall(DMHasLabel(dm, "pyop2_ghost", &isFiredrake));
    957:   if (isFiredrake) {
    958:     PetscCall(DMGetLabel(dm, "pyop2_ghost", &ghost));
    959:     PetscCall(DMLabelCreateIndex(ghost, pStart, pEnd));
    960:   } else {
    961:     PetscSF sf;

    963:     PetscCall(DMGetPointSF(dm, &sf));
    964:     PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL));
    965:     nleaves = PetscMax(nleaves, 0);
    966:   }

    968:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->cellCounts));
    969:   PetscCall(PetscObjectSetName((PetscObject)patch->cellCounts, "Patch Cell Layout"));
    970:   cellCounts = patch->cellCounts;
    971:   PetscCall(PetscSectionSetChart(cellCounts, vStart, vEnd));
    972:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->pointCounts));
    973:   PetscCall(PetscObjectSetName((PetscObject)patch->pointCounts, "Patch Point Layout"));
    974:   pointCounts = patch->pointCounts;
    975:   PetscCall(PetscSectionSetChart(pointCounts, vStart, vEnd));
    976:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->extFacetCounts));
    977:   PetscCall(PetscObjectSetName((PetscObject)patch->extFacetCounts, "Patch Exterior Facet Layout"));
    978:   extFacetCounts = patch->extFacetCounts;
    979:   PetscCall(PetscSectionSetChart(extFacetCounts, vStart, vEnd));
    980:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->intFacetCounts));
    981:   PetscCall(PetscObjectSetName((PetscObject)patch->intFacetCounts, "Patch Interior Facet Layout"));
    982:   intFacetCounts = patch->intFacetCounts;
    983:   PetscCall(PetscSectionSetChart(intFacetCounts, vStart, vEnd));
    984:   /* Count cells and points in the patch surrounding each entity */
    985:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
    986:   for (v = vStart; v < vEnd; ++v) {
    987:     PetscHashIter hi;
    988:     PetscInt      chtSize, loc = -1;
    989:     PetscBool     flg;

    991:     if (!patch->user_patches && patch->ctype != PC_PATCH_PARDECOMP) {
    992:       if (ghost) PetscCall(DMLabelHasPoint(ghost, v, &flg));
    993:       else {
    994:         PetscCall(PetscFindInt(v, nleaves, leaves, &loc));
    995:         flg = loc >= 0 ? PETSC_TRUE : PETSC_FALSE;
    996:       }
    997:       /* Not an owned entity, don't make a cell patch. */
    998:       if (flg) continue;
    999:     }

    1001:     PetscCall(patch->patchconstructop((void *)patch, dm, v, ht));
    1002:     PetscCall(PCPatchCompleteCellPatch(pc, ht, cht));
    1003:     PetscCall(PetscHSetIGetSize(cht, &chtSize));
    1004:     /* empty patch, continue */
    1005:     if (chtSize == 0) continue;

    1007:     /* safe because size(cht) > 0 from above */
    1008:     PetscHashIterBegin(cht, hi);
    1009:     while (!PetscHashIterAtEnd(cht, hi)) {
    1010:       PetscInt point, pdof;

    1012:       PetscHashIterGetKey(cht, hi, point);
    1013:       if (fStart <= point && point < fEnd) {
    1014:         const PetscInt *support;
    1015:         PetscInt        supportSize, p;
    1016:         PetscBool       interior = PETSC_TRUE;
    1017:         PetscCall(DMPlexGetSupport(dm, point, &support));
    1018:         PetscCall(DMPlexGetSupportSize(dm, point, &supportSize));
    1019:         if (supportSize == 1) {
    1020:           interior = PETSC_FALSE;
    1021:         } else {
    1022:           for (p = 0; p < supportSize; p++) {
    1023:             PetscBool found;
    1024:             /* FIXME: can I do this while iterating over cht? */
    1025:             PetscCall(PetscHSetIHas(cht, support[p], &found));
    1026:             if (!found) {
    1027:               interior = PETSC_FALSE;
    1028:               break;
    1029:             }
    1030:           }
    1031:         }
    1032:         if (interior) {
    1033:           PetscCall(PetscSectionAddDof(intFacetCounts, v, 1));
    1034:         } else {
    1035:           PetscCall(PetscSectionAddDof(extFacetCounts, v, 1));
    1036:         }
    1037:       }
    1038:       PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL));
    1039:       if (pdof) PetscCall(PetscSectionAddDof(pointCounts, v, 1));
    1040:       if (point >= cStart && point < cEnd) PetscCall(PetscSectionAddDof(cellCounts, v, 1));
    1041:       PetscHashIterNext(cht, hi);
    1042:     }
    1043:   }
    1044:   if (isFiredrake) PetscCall(DMLabelDestroyIndex(ghost));

    1046:   PetscCall(PetscSectionSetUp(cellCounts));
    1047:   PetscCall(PetscSectionGetStorageSize(cellCounts, &numCells));
    1048:   PetscCall(PetscMalloc1(numCells, &cellsArray));
    1049:   PetscCall(PetscSectionSetUp(pointCounts));
    1050:   PetscCall(PetscSectionGetStorageSize(pointCounts, &numPoints));
    1051:   PetscCall(PetscMalloc1(numPoints, &pointsArray));

    1053:   PetscCall(PetscSectionSetUp(intFacetCounts));
    1054:   PetscCall(PetscSectionSetUp(extFacetCounts));
    1055:   PetscCall(PetscSectionGetStorageSize(intFacetCounts, &numIntFacets));
    1056:   PetscCall(PetscSectionGetStorageSize(extFacetCounts, &numExtFacets));
    1057:   PetscCall(PetscMalloc1(numIntFacets, &intFacetsArray));
    1058:   PetscCall(PetscMalloc1(numIntFacets * 2, &intFacetsToPatchCell));
    1059:   PetscCall(PetscMalloc1(numExtFacets, &extFacetsArray));

    1061:   /* Now that we know how much space we need, run through again and actually remember the cells. */
    1062:   for (v = vStart; v < vEnd; v++) {
    1063:     PetscHashIter hi;
    1064:     PetscInt      dof, off, cdof, coff, efdof, efoff, ifdof, ifoff, pdof, n = 0, cn = 0, ifn = 0, efn = 0;

    1066:     PetscCall(PetscSectionGetDof(pointCounts, v, &dof));
    1067:     PetscCall(PetscSectionGetOffset(pointCounts, v, &off));
    1068:     PetscCall(PetscSectionGetDof(cellCounts, v, &cdof));
    1069:     PetscCall(PetscSectionGetOffset(cellCounts, v, &coff));
    1070:     PetscCall(PetscSectionGetDof(intFacetCounts, v, &ifdof));
    1071:     PetscCall(PetscSectionGetOffset(intFacetCounts, v, &ifoff));
    1072:     PetscCall(PetscSectionGetDof(extFacetCounts, v, &efdof));
    1073:     PetscCall(PetscSectionGetOffset(extFacetCounts, v, &efoff));
    1074:     if (dof <= 0) continue;
    1075:     PetscCall(patch->patchconstructop((void *)patch, dm, v, ht));
    1076:     PetscCall(PCPatchCompleteCellPatch(pc, ht, cht));
    1077:     PetscHashIterBegin(cht, hi);
    1078:     while (!PetscHashIterAtEnd(cht, hi)) {
    1079:       PetscInt point;

    1081:       PetscHashIterGetKey(cht, hi, point);
    1082:       if (fStart <= point && point < fEnd) {
    1083:         const PetscInt *support;
    1084:         PetscInt        supportSize, p;
    1085:         PetscBool       interior = PETSC_TRUE;
    1086:         PetscCall(DMPlexGetSupport(dm, point, &support));
    1087:         PetscCall(DMPlexGetSupportSize(dm, point, &supportSize));
    1088:         if (supportSize == 1) {
    1089:           interior = PETSC_FALSE;
    1090:         } else {
    1091:           for (p = 0; p < supportSize; p++) {
    1092:             PetscBool found;
    1093:             /* FIXME: can I do this while iterating over cht? */
    1094:             PetscCall(PetscHSetIHas(cht, support[p], &found));
    1095:             if (!found) {
    1096:               interior = PETSC_FALSE;
    1097:               break;
    1098:             }
    1099:           }
    1100:         }
    1101:         if (interior) {
    1102:           intFacetsToPatchCell[2 * (ifoff + ifn)]     = support[0];
    1103:           intFacetsToPatchCell[2 * (ifoff + ifn) + 1] = support[1];
    1104:           intFacetsArray[ifoff + ifn++]               = point;
    1105:         } else {
    1106:           extFacetsArray[efoff + efn++] = point;
    1107:         }
    1108:       }
    1109:       PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL));
    1110:       if (pdof) pointsArray[off + n++] = point;
    1111:       if (point >= cStart && point < cEnd) cellsArray[coff + cn++] = point;
    1112:       PetscHashIterNext(cht, hi);
    1113:     }
    1114:     PetscCheck(ifn == ifdof, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of interior facets in patch %" PetscInt_FMT " is %" PetscInt_FMT ", but should be %" PetscInt_FMT, v, ifn, ifdof);
    1115:     PetscCheck(efn == efdof, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of exterior facets in patch %" PetscInt_FMT " is %" PetscInt_FMT ", but should be %" PetscInt_FMT, v, efn, efdof);
    1116:     PetscCheck(cn == cdof, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of cells in patch %" PetscInt_FMT " is %" PetscInt_FMT ", but should be %" PetscInt_FMT, v, cn, cdof);
    1117:     PetscCheck(n == dof, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of points in patch %" PetscInt_FMT " is %" PetscInt_FMT ", but should be %" PetscInt_FMT, v, n, dof);

    1119:     for (ifn = 0; ifn < ifdof; ifn++) {
    1120:       PetscInt  cell0  = intFacetsToPatchCell[2 * (ifoff + ifn)];
    1121:       PetscInt  cell1  = intFacetsToPatchCell[2 * (ifoff + ifn) + 1];
    1122:       PetscBool found0 = PETSC_FALSE, found1 = PETSC_FALSE;
    1123:       for (n = 0; n < cdof; n++) {
    1124:         if (!found0 && cell0 == cellsArray[coff + n]) {
    1125:           intFacetsToPatchCell[2 * (ifoff + ifn)] = n;
    1126:           found0                                  = PETSC_TRUE;
    1127:         }
    1128:         if (!found1 && cell1 == cellsArray[coff + n]) {
    1129:           intFacetsToPatchCell[2 * (ifoff + ifn) + 1] = n;
    1130:           found1                                      = PETSC_TRUE;
    1131:         }
    1132:         if (found0 && found1) break;
    1133:       }
    1134:       PetscCheck(found0 && found1, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Didn't manage to find local point numbers for facet support");
    1135:     }
    1136:   }
    1137:   PetscCall(PetscHSetIDestroy(&ht));
    1138:   PetscCall(PetscHSetIDestroy(&cht));

    1140:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numCells, cellsArray, PETSC_OWN_POINTER, &patch->cells));
    1141:   PetscCall(PetscObjectSetName((PetscObject)patch->cells, "Patch Cells"));
    1142:   if (patch->viewCells) {
    1143:     PetscCall(ObjectView((PetscObject)patch->cellCounts, patch->viewerCells, patch->formatCells));
    1144:     PetscCall(ObjectView((PetscObject)patch->cells, patch->viewerCells, patch->formatCells));
    1145:   }
    1146:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray, PETSC_OWN_POINTER, &patch->intFacets));
    1147:   PetscCall(PetscObjectSetName((PetscObject)patch->intFacets, "Patch Interior Facets"));
    1148:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 2 * numIntFacets, intFacetsToPatchCell, PETSC_OWN_POINTER, &patch->intFacetsToPatchCell));
    1149:   PetscCall(PetscObjectSetName((PetscObject)patch->intFacetsToPatchCell, "Patch Interior Facets local support"));
    1150:   if (patch->viewIntFacets) {
    1151:     PetscCall(ObjectView((PetscObject)patch->intFacetCounts, patch->viewerIntFacets, patch->formatIntFacets));
    1152:     PetscCall(ObjectView((PetscObject)patch->intFacets, patch->viewerIntFacets, patch->formatIntFacets));
    1153:     PetscCall(ObjectView((PetscObject)patch->intFacetsToPatchCell, patch->viewerIntFacets, patch->formatIntFacets));
    1154:   }
    1155:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numExtFacets, extFacetsArray, PETSC_OWN_POINTER, &patch->extFacets));
    1156:   PetscCall(PetscObjectSetName((PetscObject)patch->extFacets, "Patch Exterior Facets"));
    1157:   if (patch->viewExtFacets) {
    1158:     PetscCall(ObjectView((PetscObject)patch->extFacetCounts, patch->viewerExtFacets, patch->formatExtFacets));
    1159:     PetscCall(ObjectView((PetscObject)patch->extFacets, patch->viewerExtFacets, patch->formatExtFacets));
    1160:   }
    1161:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPoints, pointsArray, PETSC_OWN_POINTER, &patch->points));
    1162:   PetscCall(PetscObjectSetName((PetscObject)patch->points, "Patch Points"));
    1163:   if (patch->viewPoints) {
    1164:     PetscCall(ObjectView((PetscObject)patch->pointCounts, patch->viewerPoints, patch->formatPoints));
    1165:     PetscCall(ObjectView((PetscObject)patch->points, patch->viewerPoints, patch->formatPoints));
    1166:   }
    1167:   PetscCall(DMDestroy(&dm));
    1168:   PetscFunctionReturn(PETSC_SUCCESS);
    1169: }

    1171: /*
    1172:   PCPatchCreateCellPatchDiscretisationInfo - Build the dof maps for cell patches

    1174:   Input Parameters:
    1175:   + dm - The DMPlex object defining the mesh
    1176:   . cellCounts - Section with counts of cells around each vertex
    1177:   . cells - IS of the cell point indices of cells in each patch
    1178:   . cellNumbering - Section mapping plex cell points to Firedrake cell indices.
    1179:   . nodesPerCell - number of nodes per cell.
    1180:   - cellNodeMap - map from cells to node indices (nodesPerCell * numCells)

    1182:   Output Parameters:
    1183:   + dofs - IS of local dof numbers of each cell in the patch, where local is a patch local numbering
    1184:   . gtolCounts - Section with counts of dofs per cell patch
    1185:   - gtol - IS mapping from global dofs to local dofs for each patch.
    1186:  */
    1187: static PetscErrorCode PCPatchCreateCellPatchDiscretisationInfo(PC pc)
    1188: {
    1189:   PC_PATCH       *patch       = (PC_PATCH *)pc->data;
    1190:   PetscSection    cellCounts  = patch->cellCounts;
    1191:   PetscSection    pointCounts = patch->pointCounts;
    1192:   PetscSection    gtolCounts, gtolCountsWithArtificial = NULL, gtolCountsWithAll = NULL;
    1193:   IS              cells         = patch->cells;
    1194:   IS              points        = patch->points;
    1195:   PetscSection    cellNumbering = patch->cellNumbering;
    1196:   PetscInt        Nf            = patch->nsubspaces;
    1197:   PetscInt        numCells, numPoints;
    1198:   PetscInt        numDofs;
    1199:   PetscInt        numGlobalDofs, numGlobalDofsWithArtificial, numGlobalDofsWithAll;
    1200:   PetscInt        totalDofsPerCell = patch->totalDofsPerCell;
    1201:   PetscInt        vStart, vEnd, v;
    1202:   const PetscInt *cellsArray, *pointsArray;
    1203:   PetscInt       *newCellsArray                 = NULL;
    1204:   PetscInt       *dofsArray                     = NULL;
    1205:   PetscInt       *dofsArrayWithArtificial       = NULL;
    1206:   PetscInt       *dofsArrayWithAll              = NULL;
    1207:   PetscInt       *offsArray                     = NULL;
    1208:   PetscInt       *offsArrayWithArtificial       = NULL;
    1209:   PetscInt       *offsArrayWithAll              = NULL;
    1210:   PetscInt       *asmArray                      = NULL;
    1211:   PetscInt       *asmArrayWithArtificial        = NULL;
    1212:   PetscInt       *asmArrayWithAll               = NULL;
    1213:   PetscInt       *globalDofsArray               = NULL;
    1214:   PetscInt       *globalDofsArrayWithArtificial = NULL;
    1215:   PetscInt       *globalDofsArrayWithAll        = NULL;
    1216:   PetscInt        globalIndex                   = 0;
    1217:   PetscInt        key                           = 0;
    1218:   PetscInt        asmKey                        = 0;
    1219:   DM              dm                            = NULL, plex;
    1220:   const PetscInt *bcNodes                       = NULL;
    1221:   PetscHMapI      ht;
    1222:   PetscHMapI      htWithArtificial;
    1223:   PetscHMapI      htWithAll;
    1224:   PetscHSetI      globalBcs;
    1225:   PetscInt        numBcs;
    1226:   PetscHSetI      ownedpts, seenpts, owneddofs, seendofs, artificialbcs;
    1227:   PetscInt        pStart, pEnd, p, i;
    1228:   char            option[PETSC_MAX_PATH_LEN];
    1229:   PetscBool       isNonlinear;

    1231:   PetscFunctionBegin;
    1232:   PetscCall(PCGetDM(pc, &dm));
    1233:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    1234:   dm = plex;
    1235:   /* dofcounts section is cellcounts section * dofPerCell */
    1236:   PetscCall(PetscSectionGetStorageSize(cellCounts, &numCells));
    1237:   PetscCall(PetscSectionGetStorageSize(patch->pointCounts, &numPoints));
    1238:   numDofs = numCells * totalDofsPerCell;
    1239:   PetscCall(PetscMalloc1(numDofs, &dofsArray));
    1240:   PetscCall(PetscMalloc1(numPoints * Nf, &offsArray));
    1241:   PetscCall(PetscMalloc1(numDofs, &asmArray));
    1242:   PetscCall(PetscMalloc1(numCells, &newCellsArray));
    1243:   PetscCall(PetscSectionGetChart(cellCounts, &vStart, &vEnd));
    1244:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCounts));
    1245:   gtolCounts = patch->gtolCounts;
    1246:   PetscCall(PetscSectionSetChart(gtolCounts, vStart, vEnd));
    1247:   PetscCall(PetscObjectSetName((PetscObject)patch->gtolCounts, "Patch Global Index Section"));

    1249:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
    1250:     PetscCall(PetscMalloc1(numPoints * Nf, &offsArrayWithArtificial));
    1251:     PetscCall(PetscMalloc1(numDofs, &asmArrayWithArtificial));
    1252:     PetscCall(PetscMalloc1(numDofs, &dofsArrayWithArtificial));
    1253:     PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithArtificial));
    1254:     gtolCountsWithArtificial = patch->gtolCountsWithArtificial;
    1255:     PetscCall(PetscSectionSetChart(gtolCountsWithArtificial, vStart, vEnd));
    1256:     PetscCall(PetscObjectSetName((PetscObject)patch->gtolCountsWithArtificial, "Patch Global Index Section Including Artificial BCs"));
    1257:   }

    1259:   isNonlinear = patch->isNonlinear;
    1260:   if (isNonlinear) {
    1261:     PetscCall(PetscMalloc1(numPoints * Nf, &offsArrayWithAll));
    1262:     PetscCall(PetscMalloc1(numDofs, &asmArrayWithAll));
    1263:     PetscCall(PetscMalloc1(numDofs, &dofsArrayWithAll));
    1264:     PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithAll));
    1265:     gtolCountsWithAll = patch->gtolCountsWithAll;
    1266:     PetscCall(PetscSectionSetChart(gtolCountsWithAll, vStart, vEnd));
    1267:     PetscCall(PetscObjectSetName((PetscObject)patch->gtolCountsWithAll, "Patch Global Index Section Including All BCs"));
    1268:   }

    1270:   /* Outside the patch loop, get the dofs that are globally-enforced Dirichlet
    1271:    conditions */
    1272:   PetscCall(PetscHSetICreate(&globalBcs));
    1273:   PetscCall(ISGetIndices(patch->ghostBcNodes, &bcNodes));
    1274:   PetscCall(ISGetSize(patch->ghostBcNodes, &numBcs));
    1275:   for (i = 0; i < numBcs; ++i) PetscCall(PetscHSetIAdd(globalBcs, bcNodes[i])); /* these are already in concatenated numbering */
    1276:   PetscCall(ISRestoreIndices(patch->ghostBcNodes, &bcNodes));
    1277:   PetscCall(ISDestroy(&patch->ghostBcNodes)); /* memory optimisation */

    1279:   /* Hash tables for artificial BC construction */
    1280:   PetscCall(PetscHSetICreate(&ownedpts));
    1281:   PetscCall(PetscHSetICreate(&seenpts));
    1282:   PetscCall(PetscHSetICreate(&owneddofs));
    1283:   PetscCall(PetscHSetICreate(&seendofs));
    1284:   PetscCall(PetscHSetICreate(&artificialbcs));

    1286:   PetscCall(ISGetIndices(cells, &cellsArray));
    1287:   PetscCall(ISGetIndices(points, &pointsArray));
    1288:   PetscCall(PetscHMapICreate(&ht));
    1289:   PetscCall(PetscHMapICreate(&htWithArtificial));
    1290:   PetscCall(PetscHMapICreate(&htWithAll));
    1291:   for (v = vStart; v < vEnd; ++v) {
    1292:     PetscInt localIndex               = 0;
    1293:     PetscInt localIndexWithArtificial = 0;
    1294:     PetscInt localIndexWithAll        = 0;
    1295:     PetscInt dof, off, i, j, k, l;

    1297:     PetscCall(PetscHMapIClear(ht));
    1298:     PetscCall(PetscHMapIClear(htWithArtificial));
    1299:     PetscCall(PetscHMapIClear(htWithAll));
    1300:     PetscCall(PetscSectionGetDof(cellCounts, v, &dof));
    1301:     PetscCall(PetscSectionGetOffset(cellCounts, v, &off));
    1302:     if (dof <= 0) continue;

    1304:     /* Calculate the global numbers of the artificial BC dofs here first */
    1305:     PetscCall(patch->patchconstructop((void *)patch, dm, v, ownedpts));
    1306:     PetscCall(PCPatchCompleteCellPatch(pc, ownedpts, seenpts));
    1307:     PetscCall(PCPatchGetPointDofs(pc, ownedpts, owneddofs, v, &patch->subspaces_to_exclude));
    1308:     PetscCall(PCPatchGetPointDofs(pc, seenpts, seendofs, v, NULL));
    1309:     PetscCall(PCPatchComputeSetDifference_Private(owneddofs, seendofs, artificialbcs));
    1310:     if (patch->viewPatches) {
    1311:       PetscHSetI    globalbcdofs;
    1312:       PetscHashIter hi;
    1313:       MPI_Comm      comm = PetscObjectComm((PetscObject)pc);

    1315:       PetscCall(PetscHSetICreate(&globalbcdofs));
    1316:       PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": owned dofs:\n", v));
    1317:       PetscHashIterBegin(owneddofs, hi);
    1318:       while (!PetscHashIterAtEnd(owneddofs, hi)) {
    1319:         PetscInt globalDof;

    1321:         PetscHashIterGetKey(owneddofs, hi, globalDof);
    1322:         PetscHashIterNext(owneddofs, hi);
    1323:         PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));
    1324:       }
    1325:       PetscCall(PetscSynchronizedPrintf(comm, "\n"));
    1326:       PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": seen dofs:\n", v));
    1327:       PetscHashIterBegin(seendofs, hi);
    1328:       while (!PetscHashIterAtEnd(seendofs, hi)) {
    1329:         PetscInt  globalDof;
    1330:         PetscBool flg;

    1332:         PetscHashIterGetKey(seendofs, hi, globalDof);
    1333:         PetscHashIterNext(seendofs, hi);
    1334:         PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));

    1336:         PetscCall(PetscHSetIHas(globalBcs, globalDof, &flg));
    1337:         if (flg) PetscCall(PetscHSetIAdd(globalbcdofs, globalDof));
    1338:       }
    1339:       PetscCall(PetscSynchronizedPrintf(comm, "\n"));
    1340:       PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": global BCs:\n", v));
    1341:       PetscCall(PetscHSetIGetSize(globalbcdofs, &numBcs));
    1342:       if (numBcs > 0) {
    1343:         PetscHashIterBegin(globalbcdofs, hi);
    1344:         while (!PetscHashIterAtEnd(globalbcdofs, hi)) {
    1345:           PetscInt globalDof;
    1346:           PetscHashIterGetKey(globalbcdofs, hi, globalDof);
    1347:           PetscHashIterNext(globalbcdofs, hi);
    1348:           PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));
    1349:         }
    1350:       }
    1351:       PetscCall(PetscSynchronizedPrintf(comm, "\n"));
    1352:       PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": artificial BCs:\n", v));
    1353:       PetscCall(PetscHSetIGetSize(artificialbcs, &numBcs));
    1354:       if (numBcs > 0) {
    1355:         PetscHashIterBegin(artificialbcs, hi);
    1356:         while (!PetscHashIterAtEnd(artificialbcs, hi)) {
    1357:           PetscInt globalDof;
    1358:           PetscHashIterGetKey(artificialbcs, hi, globalDof);
    1359:           PetscHashIterNext(artificialbcs, hi);
    1360:           PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));
    1361:         }
    1362:       }
    1363:       PetscCall(PetscSynchronizedPrintf(comm, "\n\n"));
    1364:       PetscCall(PetscHSetIDestroy(&globalbcdofs));
    1365:     }
    1366:     for (k = 0; k < patch->nsubspaces; ++k) {
    1367:       const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
    1368:       PetscInt        nodesPerCell   = patch->nodesPerCell[k];
    1369:       PetscInt        subspaceOffset = patch->subspaceOffsets[k];
    1370:       PetscInt        bs             = patch->bs[k];

    1372:       for (i = off; i < off + dof; ++i) {
    1373:         /* Walk over the cells in this patch. */
    1374:         const PetscInt c    = cellsArray[i];
    1375:         PetscInt       cell = c;

    1377:         /* TODO Change this to an IS */
    1378:         if (cellNumbering) {
    1379:           PetscCall(PetscSectionGetDof(cellNumbering, c, &cell));
    1380:           PetscCheck(cell > 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Cell %" PetscInt_FMT " doesn't appear in cell numbering map", c);
    1381:           PetscCall(PetscSectionGetOffset(cellNumbering, c, &cell));
    1382:         }
    1383:         newCellsArray[i] = cell;
    1384:         for (j = 0; j < nodesPerCell; ++j) {
    1385:           /* For each global dof, map it into contiguous local storage. */
    1386:           const PetscInt globalDof = cellNodeMap[cell * nodesPerCell + j] * bs + subspaceOffset;
    1387:           /* finally, loop over block size */
    1388:           for (l = 0; l < bs; ++l) {
    1389:             PetscInt  localDof;
    1390:             PetscBool isGlobalBcDof, isArtificialBcDof;

    1392:             /* first, check if this is either a globally enforced or locally enforced BC dof */
    1393:             PetscCall(PetscHSetIHas(globalBcs, globalDof + l, &isGlobalBcDof));
    1394:             PetscCall(PetscHSetIHas(artificialbcs, globalDof + l, &isArtificialBcDof));

    1396:             /* if it's either, don't ever give it a local dof number */
    1397:             if (isGlobalBcDof || isArtificialBcDof) {
    1398:               dofsArray[globalIndex] = -1; /* don't use this in assembly in this patch */
    1399:             } else {
    1400:               PetscCall(PetscHMapIGet(ht, globalDof + l, &localDof));
    1401:               if (localDof == -1) {
    1402:                 localDof = localIndex++;
    1403:                 PetscCall(PetscHMapISet(ht, globalDof + l, localDof));
    1404:               }
    1405:               PetscCheck(globalIndex < numDofs, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %" PetscInt_FMT " than expected %" PetscInt_FMT, globalIndex + 1, numDofs);
    1406:               /* And store. */
    1407:               dofsArray[globalIndex] = localDof;
    1408:             }

    1410:             if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
    1411:               if (isGlobalBcDof) {
    1412:                 dofsArrayWithArtificial[globalIndex] = -1; /* don't use this in assembly in this patch */
    1413:               } else {
    1414:                 PetscCall(PetscHMapIGet(htWithArtificial, globalDof + l, &localDof));
    1415:                 if (localDof == -1) {
    1416:                   localDof = localIndexWithArtificial++;
    1417:                   PetscCall(PetscHMapISet(htWithArtificial, globalDof + l, localDof));
    1418:                 }
    1419:                 PetscCheck(globalIndex < numDofs, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %" PetscInt_FMT " than expected %" PetscInt_FMT, globalIndex + 1, numDofs);
    1420:                 /* And store.*/
    1421:                 dofsArrayWithArtificial[globalIndex] = localDof;
    1422:               }
    1423:             }

    1425:             if (isNonlinear) {
    1426:               /* Build the dofmap for the function space with _all_ dofs,
    1427:    including those in any kind of boundary condition */
    1428:               PetscCall(PetscHMapIGet(htWithAll, globalDof + l, &localDof));
    1429:               if (localDof == -1) {
    1430:                 localDof = localIndexWithAll++;
    1431:                 PetscCall(PetscHMapISet(htWithAll, globalDof + l, localDof));
    1432:               }
    1433:               PetscCheck(globalIndex < numDofs, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %" PetscInt_FMT " than expected %" PetscInt_FMT, globalIndex + 1, numDofs);
    1434:               /* And store.*/
    1435:               dofsArrayWithAll[globalIndex] = localDof;
    1436:             }
    1437:             globalIndex++;
    1438:           }
    1439:         }
    1440:       }
    1441:     }
    1442:     /* How many local dofs in this patch? */
    1443:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
    1444:       PetscCall(PetscHMapIGetSize(htWithArtificial, &dof));
    1445:       PetscCall(PetscSectionSetDof(gtolCountsWithArtificial, v, dof));
    1446:     }
    1447:     if (isNonlinear) {
    1448:       PetscCall(PetscHMapIGetSize(htWithAll, &dof));
    1449:       PetscCall(PetscSectionSetDof(gtolCountsWithAll, v, dof));
    1450:     }
    1451:     PetscCall(PetscHMapIGetSize(ht, &dof));
    1452:     PetscCall(PetscSectionSetDof(gtolCounts, v, dof));
    1453:   }

    1455:   PetscCall(DMDestroy(&dm));
    1456:   PetscCheck(globalIndex == numDofs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Expected number of dofs (%" PetscInt_FMT ") doesn't match found number (%" PetscInt_FMT ")", numDofs, globalIndex);
    1457:   PetscCall(PetscSectionSetUp(gtolCounts));
    1458:   PetscCall(PetscSectionGetStorageSize(gtolCounts, &numGlobalDofs));
    1459:   PetscCall(PetscMalloc1(numGlobalDofs, &globalDofsArray));

    1461:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
    1462:     PetscCall(PetscSectionSetUp(gtolCountsWithArtificial));
    1463:     PetscCall(PetscSectionGetStorageSize(gtolCountsWithArtificial, &numGlobalDofsWithArtificial));
    1464:     PetscCall(PetscMalloc1(numGlobalDofsWithArtificial, &globalDofsArrayWithArtificial));
    1465:   }
    1466:   if (isNonlinear) {
    1467:     PetscCall(PetscSectionSetUp(gtolCountsWithAll));
    1468:     PetscCall(PetscSectionGetStorageSize(gtolCountsWithAll, &numGlobalDofsWithAll));
    1469:     PetscCall(PetscMalloc1(numGlobalDofsWithAll, &globalDofsArrayWithAll));
    1470:   }
    1471:   /* Now populate the global to local map.  This could be merged into the above loop if we were willing to deal with reallocs. */
    1472:   for (v = vStart; v < vEnd; ++v) {
    1473:     PetscHashIter hi;
    1474:     PetscInt      dof, off, Np, ooff, i, j, k, l;

    1476:     PetscCall(PetscHMapIClear(ht));
    1477:     PetscCall(PetscHMapIClear(htWithArtificial));
    1478:     PetscCall(PetscHMapIClear(htWithAll));
    1479:     PetscCall(PetscSectionGetDof(cellCounts, v, &dof));
    1480:     PetscCall(PetscSectionGetOffset(cellCounts, v, &off));
    1481:     PetscCall(PetscSectionGetDof(pointCounts, v, &Np));
    1482:     PetscCall(PetscSectionGetOffset(pointCounts, v, &ooff));
    1483:     if (dof <= 0) continue;

    1485:     for (k = 0; k < patch->nsubspaces; ++k) {
    1486:       const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
    1487:       PetscInt        nodesPerCell   = patch->nodesPerCell[k];
    1488:       PetscInt        subspaceOffset = patch->subspaceOffsets[k];
    1489:       PetscInt        bs             = patch->bs[k];
    1490:       PetscInt        goff;

    1492:       for (i = off; i < off + dof; ++i) {
    1493:         /* Reconstruct mapping of global-to-local on this patch. */
    1494:         const PetscInt c    = cellsArray[i];
    1495:         PetscInt       cell = c;

    1497:         if (cellNumbering) PetscCall(PetscSectionGetOffset(cellNumbering, c, &cell));
    1498:         for (j = 0; j < nodesPerCell; ++j) {
    1499:           for (l = 0; l < bs; ++l) {
    1500:             const PetscInt globalDof = cellNodeMap[cell * nodesPerCell + j] * bs + l + subspaceOffset;
    1501:             const PetscInt localDof  = dofsArray[key];
    1502:             if (localDof >= 0) PetscCall(PetscHMapISet(ht, globalDof, localDof));
    1503:             if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
    1504:               const PetscInt localDofWithArtificial = dofsArrayWithArtificial[key];
    1505:               if (localDofWithArtificial >= 0) PetscCall(PetscHMapISet(htWithArtificial, globalDof, localDofWithArtificial));
    1506:             }
    1507:             if (isNonlinear) {
    1508:               const PetscInt localDofWithAll = dofsArrayWithAll[key];
    1509:               if (localDofWithAll >= 0) PetscCall(PetscHMapISet(htWithAll, globalDof, localDofWithAll));
    1510:             }
    1511:             key++;
    1512:           }
    1513:         }
    1514:       }

    1516:       /* Shove it in the output data structure. */
    1517:       PetscCall(PetscSectionGetOffset(gtolCounts, v, &goff));
    1518:       PetscHashIterBegin(ht, hi);
    1519:       while (!PetscHashIterAtEnd(ht, hi)) {
    1520:         PetscInt globalDof, localDof;

    1522:         PetscHashIterGetKey(ht, hi, globalDof);
    1523:         PetscHashIterGetVal(ht, hi, localDof);
    1524:         if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof;
    1525:         PetscHashIterNext(ht, hi);
    1526:       }

    1528:       if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
    1529:         PetscCall(PetscSectionGetOffset(gtolCountsWithArtificial, v, &goff));
    1530:         PetscHashIterBegin(htWithArtificial, hi);
    1531:         while (!PetscHashIterAtEnd(htWithArtificial, hi)) {
    1532:           PetscInt globalDof, localDof;
    1533:           PetscHashIterGetKey(htWithArtificial, hi, globalDof);
    1534:           PetscHashIterGetVal(htWithArtificial, hi, localDof);
    1535:           if (globalDof >= 0) globalDofsArrayWithArtificial[goff + localDof] = globalDof;
    1536:           PetscHashIterNext(htWithArtificial, hi);
    1537:         }
    1538:       }
    1539:       if (isNonlinear) {
    1540:         PetscCall(PetscSectionGetOffset(gtolCountsWithAll, v, &goff));
    1541:         PetscHashIterBegin(htWithAll, hi);
    1542:         while (!PetscHashIterAtEnd(htWithAll, hi)) {
    1543:           PetscInt globalDof, localDof;
    1544:           PetscHashIterGetKey(htWithAll, hi, globalDof);
    1545:           PetscHashIterGetVal(htWithAll, hi, localDof);
    1546:           if (globalDof >= 0) globalDofsArrayWithAll[goff + localDof] = globalDof;
    1547:           PetscHashIterNext(htWithAll, hi);
    1548:         }
    1549:       }

    1551:       for (p = 0; p < Np; ++p) {
    1552:         const PetscInt point = pointsArray[ooff + p];
    1553:         PetscInt       globalDof, localDof;

    1555:         PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof));
    1556:         PetscCall(PetscHMapIGet(ht, globalDof, &localDof));
    1557:         offsArray[(ooff + p) * Nf + k] = localDof;
    1558:         if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
    1559:           PetscCall(PetscHMapIGet(htWithArtificial, globalDof, &localDof));
    1560:           offsArrayWithArtificial[(ooff + p) * Nf + k] = localDof;
    1561:         }
    1562:         if (isNonlinear) {
    1563:           PetscCall(PetscHMapIGet(htWithAll, globalDof, &localDof));
    1564:           offsArrayWithAll[(ooff + p) * Nf + k] = localDof;
    1565:         }
    1566:       }
    1567:     }

    1569:     PetscCall(PetscHSetIDestroy(&globalBcs));
    1570:     PetscCall(PetscHSetIDestroy(&ownedpts));
    1571:     PetscCall(PetscHSetIDestroy(&seenpts));
    1572:     PetscCall(PetscHSetIDestroy(&owneddofs));
    1573:     PetscCall(PetscHSetIDestroy(&seendofs));
    1574:     PetscCall(PetscHSetIDestroy(&artificialbcs));

    1576:     /* At this point, we have a hash table ht built that maps globalDof -> localDof.
    1577:    We need to create the dof table laid out cellwise first, then by subspace,
    1578:    as the assembler assembles cell-wise and we need to stuff the different
    1579:    contributions of the different function spaces to the right places. So we loop
    1580:    over cells, then over subspaces. */
    1581:     if (patch->nsubspaces > 1) { /* for nsubspaces = 1, data we need is already in dofsArray */
    1582:       for (i = off; i < off + dof; ++i) {
    1583:         const PetscInt c    = cellsArray[i];
    1584:         PetscInt       cell = c;

    1586:         if (cellNumbering) PetscCall(PetscSectionGetOffset(cellNumbering, c, &cell));
    1587:         for (k = 0; k < patch->nsubspaces; ++k) {
    1588:           const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
    1589:           PetscInt        nodesPerCell   = patch->nodesPerCell[k];
    1590:           PetscInt        subspaceOffset = patch->subspaceOffsets[k];
    1591:           PetscInt        bs             = patch->bs[k];

    1593:           for (j = 0; j < nodesPerCell; ++j) {
    1594:             for (l = 0; l < bs; ++l) {
    1595:               const PetscInt globalDof = cellNodeMap[cell * nodesPerCell + j] * bs + l + subspaceOffset;
    1596:               PetscInt       localDof;

    1598:               PetscCall(PetscHMapIGet(ht, globalDof, &localDof));
    1599:               /* If it's not in the hash table, i.e. is a BC dof,
    1600:    then the PetscHSetIMap above gives -1, which matches
    1601:    exactly the convention for PETSc's matrix assembly to
    1602:    ignore the dof. So we don't need to do anything here */
    1603:               asmArray[asmKey] = localDof;
    1604:               if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
    1605:                 PetscCall(PetscHMapIGet(htWithArtificial, globalDof, &localDof));
    1606:                 asmArrayWithArtificial[asmKey] = localDof;
    1607:               }
    1608:               if (isNonlinear) {
    1609:                 PetscCall(PetscHMapIGet(htWithAll, globalDof, &localDof));
    1610:                 asmArrayWithAll[asmKey] = localDof;
    1611:               }
    1612:               asmKey++;
    1613:             }
    1614:           }
    1615:         }
    1616:       }
    1617:     }
    1618:   }
    1619:   if (1 == patch->nsubspaces) {
    1620:     PetscCall(PetscArraycpy(asmArray, dofsArray, numDofs));
    1621:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscCall(PetscArraycpy(asmArrayWithArtificial, dofsArrayWithArtificial, numDofs));
    1622:     if (isNonlinear) PetscCall(PetscArraycpy(asmArrayWithAll, dofsArrayWithAll, numDofs));
    1623:   }

    1625:   PetscCall(PetscHMapIDestroy(&ht));
    1626:   PetscCall(PetscHMapIDestroy(&htWithArtificial));
    1627:   PetscCall(PetscHMapIDestroy(&htWithAll));
    1628:   PetscCall(ISRestoreIndices(cells, &cellsArray));
    1629:   PetscCall(ISRestoreIndices(points, &pointsArray));
    1630:   PetscCall(PetscFree(dofsArray));
    1631:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscCall(PetscFree(dofsArrayWithArtificial));
    1632:   if (isNonlinear) PetscCall(PetscFree(dofsArrayWithAll));
    1633:   /* Create placeholder section for map from points to patch dofs */
    1634:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection));
    1635:   PetscCall(PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces));
    1636:   if (patch->combined) {
    1637:     PetscInt numFields;
    1638:     PetscCall(PetscSectionGetNumFields(patch->dofSection[0], &numFields));
    1639:     PetscCheck(numFields == patch->nsubspaces, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Mismatch between number of section fields %" PetscInt_FMT " and number of subspaces %" PetscInt_FMT, numFields, patch->nsubspaces);
    1640:     PetscCall(PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd));
    1641:     PetscCall(PetscSectionSetChart(patch->patchSection, pStart, pEnd));
    1642:     for (p = pStart; p < pEnd; ++p) {
    1643:       PetscInt dof, fdof, f;

    1645:       PetscCall(PetscSectionGetDof(patch->dofSection[0], p, &dof));
    1646:       PetscCall(PetscSectionSetDof(patch->patchSection, p, dof));
    1647:       for (f = 0; f < patch->nsubspaces; ++f) {
    1648:         PetscCall(PetscSectionGetFieldDof(patch->dofSection[0], p, f, &fdof));
    1649:         PetscCall(PetscSectionSetFieldDof(patch->patchSection, p, f, fdof));
    1650:       }
    1651:     }
    1652:   } else {
    1653:     PetscInt pStartf, pEndf, f;
    1654:     pStart = PETSC_INT_MAX;
    1655:     pEnd   = PETSC_INT_MIN;
    1656:     for (f = 0; f < patch->nsubspaces; ++f) {
    1657:       PetscCall(PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf));
    1658:       pStart = PetscMin(pStart, pStartf);
    1659:       pEnd   = PetscMax(pEnd, pEndf);
    1660:     }
    1661:     PetscCall(PetscSectionSetChart(patch->patchSection, pStart, pEnd));
    1662:     for (f = 0; f < patch->nsubspaces; ++f) {
    1663:       PetscCall(PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf));
    1664:       for (p = pStartf; p < pEndf; ++p) {
    1665:         PetscInt fdof;
    1666:         PetscCall(PetscSectionGetDof(patch->dofSection[f], p, &fdof));
    1667:         PetscCall(PetscSectionAddDof(patch->patchSection, p, fdof));
    1668:         PetscCall(PetscSectionSetFieldDof(patch->patchSection, p, f, fdof));
    1669:       }
    1670:     }
    1671:   }
    1672:   PetscCall(PetscSectionSetUp(patch->patchSection));
    1673:   PetscCall(PetscSectionSetUseFieldOffsets(patch->patchSection, PETSC_TRUE));
    1674:   /* Replace cell indices with firedrake-numbered ones. */
    1675:   PetscCall(ISGeneralSetIndices(cells, numCells, (const PetscInt *)newCellsArray, PETSC_OWN_POINTER));
    1676:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofs, globalDofsArray, PETSC_OWN_POINTER, &patch->gtol));
    1677:   PetscCall(PetscObjectSetName((PetscObject)patch->gtol, "Global Indices"));
    1678:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_g2l_view", patch->classname));
    1679:   PetscCall(PetscSectionViewFromOptions(patch->gtolCounts, (PetscObject)pc, option));
    1680:   PetscCall(ISViewFromOptions(patch->gtol, (PetscObject)pc, option));
    1681:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArray, PETSC_OWN_POINTER, &patch->dofs));
    1682:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPoints * Nf, offsArray, PETSC_OWN_POINTER, &patch->offs));
    1683:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
    1684:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithArtificial, globalDofsArrayWithArtificial, PETSC_OWN_POINTER, &patch->gtolWithArtificial));
    1685:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithArtificial, PETSC_OWN_POINTER, &patch->dofsWithArtificial));
    1686:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPoints * Nf, offsArrayWithArtificial, PETSC_OWN_POINTER, &patch->offsWithArtificial));
    1687:   }
    1688:   if (isNonlinear) {
    1689:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithAll, globalDofsArrayWithAll, PETSC_OWN_POINTER, &patch->gtolWithAll));
    1690:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithAll, PETSC_OWN_POINTER, &patch->dofsWithAll));
    1691:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPoints * Nf, offsArrayWithAll, PETSC_OWN_POINTER, &patch->offsWithAll));
    1692:   }
    1693:   PetscFunctionReturn(PETSC_SUCCESS);
    1694: }

    1696: static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat, PetscBool withArtificial)
    1697: {
    1698:   PC_PATCH   *patch = (PC_PATCH *)pc->data;
    1699:   PetscBool   flg;
    1700:   PetscInt    csize, rsize;
    1701:   const char *prefix = NULL;

    1703:   PetscFunctionBegin;
    1704:   if (withArtificial) {
    1705:     /* would be nice if we could create a rectangular matrix of size numDofsWithArtificial x numDofs here */
    1706:     PetscInt pStart;
    1707:     PetscCall(PetscSectionGetChart(patch->gtolCountsWithArtificial, &pStart, NULL));
    1708:     PetscCall(PetscSectionGetDof(patch->gtolCountsWithArtificial, point + pStart, &rsize));
    1709:     csize = rsize;
    1710:   } else {
    1711:     PetscInt pStart;
    1712:     PetscCall(PetscSectionGetChart(patch->gtolCounts, &pStart, NULL));
    1713:     PetscCall(PetscSectionGetDof(patch->gtolCounts, point + pStart, &rsize));
    1714:     csize = rsize;
    1715:   }

    1717:   PetscCall(MatCreate(PETSC_COMM_SELF, mat));
    1718:   PetscCall(PCGetOptionsPrefix(pc, &prefix));
    1719:   PetscCall(MatSetOptionsPrefix(*mat, prefix));
    1720:   PetscCall(MatAppendOptionsPrefix(*mat, "pc_patch_sub_"));
    1721:   if (patch->sub_mat_type) PetscCall(MatSetType(*mat, patch->sub_mat_type));
    1722:   else if (!patch->sub_mat_type) PetscCall(MatSetType(*mat, MATDENSE));
    1723:   PetscCall(MatSetSizes(*mat, rsize, csize, rsize, csize));
    1724:   PetscCall(PetscObjectTypeCompare((PetscObject)*mat, MATDENSE, &flg));
    1725:   if (!flg) PetscCall(PetscObjectTypeCompare((PetscObject)*mat, MATSEQDENSE, &flg));
    1726:   /* Sparse patch matrices */
    1727:   if (!flg) {
    1728:     PetscBT         bt;
    1729:     PetscInt       *dnnz      = NULL;
    1730:     const PetscInt *dofsArray = NULL;
    1731:     PetscInt        pStart, pEnd, ncell, offset, c, i, j;

    1733:     if (withArtificial) {
    1734:       PetscCall(ISGetIndices(patch->dofsWithArtificial, &dofsArray));
    1735:     } else {
    1736:       PetscCall(ISGetIndices(patch->dofs, &dofsArray));
    1737:     }
    1738:     PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));
    1739:     point += pStart;
    1740:     PetscCheck(point < pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %" PetscInt_FMT " not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, pStart, pEnd);
    1741:     PetscCall(PetscSectionGetDof(patch->cellCounts, point, &ncell));
    1742:     PetscCall(PetscSectionGetOffset(patch->cellCounts, point, &offset));
    1743:     PetscCall(PetscLogEventBegin(PC_Patch_Prealloc, pc, 0, 0, 0));
    1744:     /* A PetscBT uses N^2 bits to store the sparsity pattern on a
    1745:    * patch. This is probably OK if the patches are not too big,
    1746:    * but uses too much memory. We therefore switch based on rsize. */
    1747:     if (rsize < 3000) { /* FIXME: I picked this switch value out of my hat */
    1748:       PetscScalar *zeroes;
    1749:       PetscInt     rows;

    1751:       PetscCall(PetscCalloc1(rsize, &dnnz));
    1752:       PetscCall(PetscBTCreate(rsize * rsize, &bt));
    1753:       for (c = 0; c < ncell; ++c) {
    1754:         const PetscInt *idx = dofsArray + (offset + c) * patch->totalDofsPerCell;
    1755:         for (i = 0; i < patch->totalDofsPerCell; ++i) {
    1756:           const PetscInt row = idx[i];
    1757:           if (row < 0) continue;
    1758:           for (j = 0; j < patch->totalDofsPerCell; ++j) {
    1759:             const PetscInt col = idx[j];
    1760:             const PetscInt key = row * rsize + col;
    1761:             if (col < 0) continue;
    1762:             if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
    1763:           }
    1764:         }
    1765:       }

    1767:       if (patch->usercomputeopintfacet) {
    1768:         const PetscInt *intFacetsArray = NULL;
    1769:         PetscInt        i, numIntFacets, intFacetOffset;
    1770:         const PetscInt *facetCells = NULL;

    1772:         PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
    1773:         PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
    1774:         PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
    1775:         PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
    1776:         for (i = 0; i < numIntFacets; i++) {
    1777:           const PetscInt cell0 = facetCells[2 * (intFacetOffset + i) + 0];
    1778:           const PetscInt cell1 = facetCells[2 * (intFacetOffset + i) + 1];
    1779:           PetscInt       celli, cellj;

    1781:           for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
    1782:             const PetscInt row = dofsArray[(offset + cell0) * patch->totalDofsPerCell + celli];
    1783:             if (row < 0) continue;
    1784:             for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
    1785:               const PetscInt col = dofsArray[(offset + cell1) * patch->totalDofsPerCell + cellj];
    1786:               const PetscInt key = row * rsize + col;
    1787:               if (col < 0) continue;
    1788:               if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
    1789:             }
    1790:           }

    1792:           for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
    1793:             const PetscInt row = dofsArray[(offset + cell1) * patch->totalDofsPerCell + celli];
    1794:             if (row < 0) continue;
    1795:             for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
    1796:               const PetscInt col = dofsArray[(offset + cell0) * patch->totalDofsPerCell + cellj];
    1797:               const PetscInt key = row * rsize + col;
    1798:               if (col < 0) continue;
    1799:               if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
    1800:             }
    1801:           }
    1802:         }
    1803:       }
    1804:       PetscCall(PetscBTDestroy(&bt));
    1805:       PetscCall(MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL));
    1806:       PetscCall(PetscFree(dnnz));

    1808:       PetscCall(PetscCalloc1(patch->totalDofsPerCell * patch->totalDofsPerCell, &zeroes));
    1809:       for (c = 0; c < ncell; ++c) {
    1810:         const PetscInt *idx = &dofsArray[(offset + c) * patch->totalDofsPerCell];
    1811:         PetscCall(MatSetValues(*mat, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, zeroes, INSERT_VALUES));
    1812:       }
    1813:       PetscCall(MatGetLocalSize(*mat, &rows, NULL));
    1814:       for (i = 0; i < rows; ++i) PetscCall(MatSetValues(*mat, 1, &i, 1, &i, zeroes, INSERT_VALUES));

    1816:       if (patch->usercomputeopintfacet) {
    1817:         const PetscInt *intFacetsArray = NULL;
    1818:         PetscInt        i, numIntFacets, intFacetOffset;
    1819:         const PetscInt *facetCells = NULL;

    1821:         PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
    1822:         PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
    1823:         PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
    1824:         PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
    1825:         for (i = 0; i < numIntFacets; i++) {
    1826:           const PetscInt  cell0    = facetCells[2 * (intFacetOffset + i) + 0];
    1827:           const PetscInt  cell1    = facetCells[2 * (intFacetOffset + i) + 1];
    1828:           const PetscInt *cell0idx = &dofsArray[(offset + cell0) * patch->totalDofsPerCell];
    1829:           const PetscInt *cell1idx = &dofsArray[(offset + cell1) * patch->totalDofsPerCell];
    1830:           PetscCall(MatSetValues(*mat, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, zeroes, INSERT_VALUES));
    1831:           PetscCall(MatSetValues(*mat, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, zeroes, INSERT_VALUES));
    1832:         }
    1833:       }

    1835:       PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY));
    1836:       PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY));

    1838:       PetscCall(PetscFree(zeroes));

    1840:     } else { /* rsize too big, use MATPREALLOCATOR */
    1841:       Mat          preallocator;
    1842:       PetscScalar *vals;

    1844:       PetscCall(PetscCalloc1(patch->totalDofsPerCell * patch->totalDofsPerCell, &vals));
    1845:       PetscCall(MatCreate(PETSC_COMM_SELF, &preallocator));
    1846:       PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
    1847:       PetscCall(MatSetSizes(preallocator, rsize, rsize, rsize, rsize));
    1848:       PetscCall(MatSetUp(preallocator));

    1850:       for (c = 0; c < ncell; ++c) {
    1851:         const PetscInt *idx = dofsArray + (offset + c) * patch->totalDofsPerCell;
    1852:         PetscCall(MatSetValues(preallocator, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, vals, INSERT_VALUES));
    1853:       }

    1855:       if (patch->usercomputeopintfacet) {
    1856:         const PetscInt *intFacetsArray = NULL;
    1857:         PetscInt        i, numIntFacets, intFacetOffset;
    1858:         const PetscInt *facetCells = NULL;

    1860:         PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
    1861:         PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
    1862:         PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
    1863:         PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
    1864:         for (i = 0; i < numIntFacets; i++) {
    1865:           const PetscInt  cell0    = facetCells[2 * (intFacetOffset + i) + 0];
    1866:           const PetscInt  cell1    = facetCells[2 * (intFacetOffset + i) + 1];
    1867:           const PetscInt *cell0idx = &dofsArray[(offset + cell0) * patch->totalDofsPerCell];
    1868:           const PetscInt *cell1idx = &dofsArray[(offset + cell1) * patch->totalDofsPerCell];
    1869:           PetscCall(MatSetValues(preallocator, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, vals, INSERT_VALUES));
    1870:           PetscCall(MatSetValues(preallocator, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, vals, INSERT_VALUES));
    1871:         }
    1872:       }

    1874:       PetscCall(PetscFree(vals));
    1875:       PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
    1876:       PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));
    1877:       PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_TRUE, *mat));
    1878:       PetscCall(MatDestroy(&preallocator));
    1879:     }
    1880:     PetscCall(PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0));
    1881:     if (withArtificial) {
    1882:       PetscCall(ISRestoreIndices(patch->dofsWithArtificial, &dofsArray));
    1883:     } else {
    1884:       PetscCall(ISRestoreIndices(patch->dofs, &dofsArray));
    1885:     }
    1886:   }
    1887:   PetscCall(MatSetUp(*mat));
    1888:   PetscFunctionReturn(PETSC_SUCCESS);
    1889: }

    1891: static PetscErrorCode PCPatchComputeFunction_DMPlex_Private(PC pc, PetscInt patchNum, Vec x, Vec F, IS cellIS, PetscInt n, const PetscInt *l2p, const PetscInt *l2pWithAll, void *ctx)
    1892: {
    1893:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
    1894:   DM              dm, plex;
    1895:   PetscSection    s;
    1896:   const PetscInt *parray, *oarray;
    1897:   PetscInt        Nf = patch->nsubspaces, Np, poff, p, f;

    1899:   PetscFunctionBegin;
    1900:   PetscCheck(!patch->precomputeElementTensors, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Precomputing element tensors not implemented with DMPlex compute operator");
    1901:   PetscCall(PCGetDM(pc, &dm));
    1902:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    1903:   dm = plex;
    1904:   PetscCall(DMGetLocalSection(dm, &s));
    1905:   /* Set offset into patch */
    1906:   PetscCall(PetscSectionGetDof(patch->pointCounts, patchNum, &Np));
    1907:   PetscCall(PetscSectionGetOffset(patch->pointCounts, patchNum, &poff));
    1908:   PetscCall(ISGetIndices(patch->points, &parray));
    1909:   PetscCall(ISGetIndices(patch->offs, &oarray));
    1910:   for (f = 0; f < Nf; ++f) {
    1911:     for (p = 0; p < Np; ++p) {
    1912:       const PetscInt point = parray[poff + p];
    1913:       PetscInt       dof;

    1915:       PetscCall(PetscSectionGetFieldDof(patch->patchSection, point, f, &dof));
    1916:       PetscCall(PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff + p) * Nf + f]));
    1917:       if (patch->nsubspaces == 1) PetscCall(PetscSectionSetOffset(patch->patchSection, point, oarray[(poff + p) * Nf + f]));
    1918:       else PetscCall(PetscSectionSetOffset(patch->patchSection, point, -1));
    1919:     }
    1920:   }
    1921:   PetscCall(ISRestoreIndices(patch->points, &parray));
    1922:   PetscCall(ISRestoreIndices(patch->offs, &oarray));
    1923:   if (patch->viewSection) PetscCall(ObjectView((PetscObject)patch->patchSection, patch->viewerSection, patch->formatSection));
    1924:   PetscCall(DMPlexComputeResidual_Patch_Internal(dm, patch->patchSection, cellIS, 0.0, x, NULL, F, ctx));
    1925:   PetscCall(DMDestroy(&dm));
    1926:   PetscFunctionReturn(PETSC_SUCCESS);
    1927: }

    1929: PetscErrorCode PCPatchComputeFunction_Internal(PC pc, Vec x, Vec F, PetscInt point)
    1930: {
    1931:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
    1932:   const PetscInt *dofsArray;
    1933:   const PetscInt *dofsArrayWithAll;
    1934:   const PetscInt *cellsArray;
    1935:   PetscInt        ncell, offset, pStart, pEnd;

    1937:   PetscFunctionBegin;
    1938:   PetscCall(PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0));
    1939:   PetscCheck(patch->usercomputeop, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set callback");
    1940:   PetscCall(ISGetIndices(patch->dofs, &dofsArray));
    1941:   PetscCall(ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll));
    1942:   PetscCall(ISGetIndices(patch->cells, &cellsArray));
    1943:   PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));

    1945:   point += pStart;
    1946:   PetscCheck(point < pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %" PetscInt_FMT " not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, pStart, pEnd);

    1948:   PetscCall(PetscSectionGetDof(patch->cellCounts, point, &ncell));
    1949:   PetscCall(PetscSectionGetOffset(patch->cellCounts, point, &offset));
    1950:   if (ncell <= 0) {
    1951:     PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
    1952:     PetscFunctionReturn(PETSC_SUCCESS);
    1953:   }
    1954:   PetscCall(VecSet(F, 0.0));
    1955:   /* Cannot reuse the same IS because the geometry info is being cached in it */
    1956:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS));
    1957:   PetscCallBack("PCPatch callback", patch->usercomputef(pc, point, x, F, patch->cellIS, ncell * patch->totalDofsPerCell, dofsArray + offset * patch->totalDofsPerCell, dofsArrayWithAll + offset * patch->totalDofsPerCell, patch->usercomputefctx));
    1958:   PetscCall(ISDestroy(&patch->cellIS));
    1959:   PetscCall(ISRestoreIndices(patch->dofs, &dofsArray));
    1960:   PetscCall(ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll));
    1961:   PetscCall(ISRestoreIndices(patch->cells, &cellsArray));
    1962:   if (patch->viewMatrix) {
    1963:     char name[PETSC_MAX_PATH_LEN];

    1965:     PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "Patch vector for Point %" PetscInt_FMT, point));
    1966:     PetscCall(PetscObjectSetName((PetscObject)F, name));
    1967:     PetscCall(ObjectView((PetscObject)F, patch->viewerMatrix, patch->formatMatrix));
    1968:   }
    1969:   PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
    1970:   PetscFunctionReturn(PETSC_SUCCESS);
    1971: }

    1973: static PetscErrorCode PCPatchComputeOperator_DMPlex_Private(PC pc, PetscInt patchNum, Vec x, Mat J, IS cellIS, PetscInt n, const PetscInt *l2p, const PetscInt *l2pWithAll, void *ctx)
    1974: {
    1975:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
    1976:   DM              dm, plex;
    1977:   PetscSection    s;
    1978:   const PetscInt *parray, *oarray;
    1979:   PetscInt        Nf = patch->nsubspaces, Np, poff, p, f;

    1981:   PetscFunctionBegin;
    1982:   PetscCall(PCGetDM(pc, &dm));
    1983:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    1984:   dm = plex;
    1985:   PetscCall(DMGetLocalSection(dm, &s));
    1986:   /* Set offset into patch */
    1987:   PetscCall(PetscSectionGetDof(patch->pointCounts, patchNum, &Np));
    1988:   PetscCall(PetscSectionGetOffset(patch->pointCounts, patchNum, &poff));
    1989:   PetscCall(ISGetIndices(patch->points, &parray));
    1990:   PetscCall(ISGetIndices(patch->offs, &oarray));
    1991:   for (f = 0; f < Nf; ++f) {
    1992:     for (p = 0; p < Np; ++p) {
    1993:       const PetscInt point = parray[poff + p];
    1994:       PetscInt       dof;

    1996:       PetscCall(PetscSectionGetFieldDof(patch->patchSection, point, f, &dof));
    1997:       PetscCall(PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff + p) * Nf + f]));
    1998:       if (patch->nsubspaces == 1) PetscCall(PetscSectionSetOffset(patch->patchSection, point, oarray[(poff + p) * Nf + f]));
    1999:       else PetscCall(PetscSectionSetOffset(patch->patchSection, point, -1));
    2000:     }
    2001:   }
    2002:   PetscCall(ISRestoreIndices(patch->points, &parray));
    2003:   PetscCall(ISRestoreIndices(patch->offs, &oarray));
    2004:   if (patch->viewSection) PetscCall(ObjectView((PetscObject)patch->patchSection, patch->viewerSection, patch->formatSection));
    2005:   /* TODO Shut off MatViewFromOptions() in MatAssemblyEnd() here */
    2006:   PetscCall(DMPlexComputeJacobian_Patch_Internal(dm, patch->patchSection, patch->patchSection, cellIS, 0.0, 0.0, x, NULL, J, J, ctx));
    2007:   PetscCall(DMDestroy(&dm));
    2008:   PetscFunctionReturn(PETSC_SUCCESS);
    2009: }

    2011: /* This function zeros mat on entry */
    2012: PetscErrorCode PCPatchComputeOperator_Internal(PC pc, Vec x, Mat mat, PetscInt point, PetscBool withArtificial)
    2013: {
    2014:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
    2015:   const PetscInt *dofsArray;
    2016:   const PetscInt *dofsArrayWithAll = NULL;
    2017:   const PetscInt *cellsArray;
    2018:   PetscInt        ncell, offset, pStart, pEnd, numIntFacets, intFacetOffset;
    2019:   PetscBool       isNonlinear;

    2021:   PetscFunctionBegin;
    2022:   PetscCall(PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0));
    2023:   isNonlinear = patch->isNonlinear;
    2024:   PetscCheck(patch->usercomputeop, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set callback");
    2025:   if (withArtificial) {
    2026:     PetscCall(ISGetIndices(patch->dofsWithArtificial, &dofsArray));
    2027:   } else {
    2028:     PetscCall(ISGetIndices(patch->dofs, &dofsArray));
    2029:   }
    2030:   if (isNonlinear) PetscCall(ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll));
    2031:   PetscCall(ISGetIndices(patch->cells, &cellsArray));
    2032:   PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));

    2034:   point += pStart;
    2035:   PetscCheck(point < pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %" PetscInt_FMT " not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, pStart, pEnd);

    2037:   PetscCall(PetscSectionGetDof(patch->cellCounts, point, &ncell));
    2038:   PetscCall(PetscSectionGetOffset(patch->cellCounts, point, &offset));
    2039:   if (ncell <= 0) {
    2040:     PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
    2041:     PetscFunctionReturn(PETSC_SUCCESS);
    2042:   }
    2043:   PetscCall(MatZeroEntries(mat));
    2044:   if (patch->precomputeElementTensors) {
    2045:     PetscInt           i;
    2046:     PetscInt           ndof = patch->totalDofsPerCell;
    2047:     const PetscScalar *elementTensors;

    2049:     PetscCall(VecGetArrayRead(patch->cellMats, &elementTensors));
    2050:     for (i = 0; i < ncell; i++) {
    2051:       const PetscInt     cell = cellsArray[i + offset];
    2052:       const PetscInt    *idx  = dofsArray + (offset + i) * ndof;
    2053:       const PetscScalar *v    = elementTensors + patch->precomputedTensorLocations[cell] * ndof * ndof;
    2054:       PetscCall(MatSetValues(mat, ndof, idx, ndof, idx, v, ADD_VALUES));
    2055:     }
    2056:     PetscCall(VecRestoreArrayRead(patch->cellMats, &elementTensors));
    2057:     PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
    2058:     PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
    2059:   } else {
    2060:     /* Cannot reuse the same IS because the geometry info is being cached in it */
    2061:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS));
    2062:     PetscCallBack("PCPatch callback",
    2063:                   patch->usercomputeop(pc, point, x, mat, patch->cellIS, ncell * patch->totalDofsPerCell, dofsArray + offset * patch->totalDofsPerCell, PetscSafePointerPlusOffset(dofsArrayWithAll, offset * patch->totalDofsPerCell), patch->usercomputeopctx));
    2064:   }
    2065:   if (patch->usercomputeopintfacet) {
    2066:     PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
    2067:     PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
    2068:     if (numIntFacets > 0) {
    2069:       /* For each interior facet, grab the two cells (in local numbering, and concatenate dof numberings for those cells) */
    2070:       PetscInt       *facetDofs = NULL, *facetDofsWithAll = NULL;
    2071:       const PetscInt *intFacetsArray = NULL;
    2072:       PetscInt        idx            = 0;
    2073:       PetscInt        i, c, d;
    2074:       PetscInt        fStart;
    2075:       DM              dm, plex;
    2076:       IS              facetIS    = NULL;
    2077:       const PetscInt *facetCells = NULL;

    2079:       PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
    2080:       PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
    2081:       PetscCall(PCGetDM(pc, &dm));
    2082:       PetscCall(DMConvert(dm, DMPLEX, &plex));
    2083:       dm = plex;
    2084:       PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, NULL));
    2085:       /* FIXME: Pull this malloc out. */
    2086:       PetscCall(PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofs));
    2087:       if (dofsArrayWithAll) PetscCall(PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofsWithAll));
    2088:       if (patch->precomputeElementTensors) {
    2089:         PetscInt           nFacetDof = 2 * patch->totalDofsPerCell;
    2090:         const PetscScalar *elementTensors;

    2092:         PetscCall(VecGetArrayRead(patch->intFacetMats, &elementTensors));

    2094:         for (i = 0; i < numIntFacets; i++) {
    2095:           const PetscInt     facet = intFacetsArray[i + intFacetOffset];
    2096:           const PetscScalar *v     = elementTensors + patch->precomputedIntFacetTensorLocations[facet - fStart] * nFacetDof * nFacetDof;
    2097:           idx                      = 0;
    2098:           /*
    2099:      0--1
    2100:      |\-|
    2101:      |+\|
    2102:      2--3
    2103:      [0, 2, 3, 0, 1, 3]
    2104:    */
    2105:           for (c = 0; c < 2; c++) {
    2106:             const PetscInt cell = facetCells[2 * (intFacetOffset + i) + c];
    2107:             for (d = 0; d < patch->totalDofsPerCell; d++) {
    2108:               facetDofs[idx] = dofsArray[(offset + cell) * patch->totalDofsPerCell + d];
    2109:               idx++;
    2110:             }
    2111:           }
    2112:           PetscCall(MatSetValues(mat, nFacetDof, facetDofs, nFacetDof, facetDofs, v, ADD_VALUES));
    2113:         }
    2114:         PetscCall(VecRestoreArrayRead(patch->intFacetMats, &elementTensors));
    2115:       } else {
    2116:         /*
    2117:      0--1
    2118:      |\-|
    2119:      |+\|
    2120:      2--3
    2121:      [0, 2, 3, 0, 1, 3]
    2122:    */
    2123:         for (i = 0; i < numIntFacets; i++) {
    2124:           for (c = 0; c < 2; c++) {
    2125:             const PetscInt cell = facetCells[2 * (intFacetOffset + i) + c];
    2126:             for (d = 0; d < patch->totalDofsPerCell; d++) {
    2127:               facetDofs[idx] = dofsArray[(offset + cell) * patch->totalDofsPerCell + d];
    2128:               if (dofsArrayWithAll) facetDofsWithAll[idx] = dofsArrayWithAll[(offset + cell) * patch->totalDofsPerCell + d];
    2129:               idx++;
    2130:             }
    2131:           }
    2132:         }
    2133:         PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray + intFacetOffset, PETSC_USE_POINTER, &facetIS));
    2134:         PetscCall(patch->usercomputeopintfacet(pc, point, x, mat, facetIS, 2 * numIntFacets * patch->totalDofsPerCell, facetDofs, facetDofsWithAll, patch->usercomputeopintfacetctx));
    2135:         PetscCall(ISDestroy(&facetIS));
    2136:       }
    2137:       PetscCall(ISRestoreIndices(patch->intFacetsToPatchCell, &facetCells));
    2138:       PetscCall(ISRestoreIndices(patch->intFacets, &intFacetsArray));
    2139:       PetscCall(PetscFree(facetDofs));
    2140:       PetscCall(PetscFree(facetDofsWithAll));
    2141:       PetscCall(DMDestroy(&dm));
    2142:     }
    2143:   }

    2145:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
    2146:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));

    2148:   if (!(withArtificial || isNonlinear) && patch->denseinverse) {
    2149:     MatFactorInfo info;
    2150:     PetscBool     flg;
    2151:     PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATSEQDENSE, &flg));
    2152:     PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Invalid Mat type for dense inverse");
    2153:     PetscCall(MatFactorInfoInitialize(&info));
    2154:     PetscCall(MatLUFactor(mat, NULL, NULL, &info));
    2155:     PetscCall(MatSeqDenseInvertFactors_Private(mat));
    2156:   }
    2157:   PetscCall(ISDestroy(&patch->cellIS));
    2158:   if (withArtificial) {
    2159:     PetscCall(ISRestoreIndices(patch->dofsWithArtificial, &dofsArray));
    2160:   } else {
    2161:     PetscCall(ISRestoreIndices(patch->dofs, &dofsArray));
    2162:   }
    2163:   if (isNonlinear) PetscCall(ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll));
    2164:   PetscCall(ISRestoreIndices(patch->cells, &cellsArray));
    2165:   if (patch->viewMatrix) {
    2166:     char name[PETSC_MAX_PATH_LEN];

    2168:     PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "Patch matrix for Point %" PetscInt_FMT, point));
    2169:     PetscCall(PetscObjectSetName((PetscObject)mat, name));
    2170:     PetscCall(ObjectView((PetscObject)mat, patch->viewerMatrix, patch->formatMatrix));
    2171:   }
    2172:   PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
    2173:   PetscFunctionReturn(PETSC_SUCCESS);
    2174: }

    2176: static PetscErrorCode MatSetValues_PCPatch_Private(Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], const PetscScalar *v, InsertMode addv)
    2177: {
    2178:   Vec          data;
    2179:   PetscScalar *array;
    2180:   PetscInt     bs, nz, i, j, cell;

    2182:   PetscCall(MatShellGetContext(mat, &data));
    2183:   PetscCall(VecGetBlockSize(data, &bs));
    2184:   PetscCall(VecGetSize(data, &nz));
    2185:   PetscCall(VecGetArray(data, &array));
    2186:   PetscCheck(m == n, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Only for square insertion");
    2187:   cell = idxm[0] / bs; /* use the fact that this is called once per cell */
    2188:   for (i = 0; i < m; i++) {
    2189:     PetscCheck(idxm[i] == idxn[i], PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Row and column indices must match!");
    2190:     for (j = 0; j < n; j++) {
    2191:       const PetscScalar v_ = v[i * bs + j];
    2192:       /* Indexing is special to the data structure we have! */
    2193:       if (addv == INSERT_VALUES) {
    2194:         array[cell * bs * bs + i * bs + j] = v_;
    2195:       } else {
    2196:         array[cell * bs * bs + i * bs + j] += v_;
    2197:       }
    2198:     }
    2199:   }
    2200:   PetscCall(VecRestoreArray(data, &array));
    2201:   PetscFunctionReturn(PETSC_SUCCESS);
    2202: }

    2204: static PetscErrorCode PCPatchPrecomputePatchTensors_Private(PC pc)
    2205: {
    2206:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
    2207:   const PetscInt *cellsArray;
    2208:   PetscInt        ncell, offset;
    2209:   const PetscInt *dofMapArray;
    2210:   PetscInt        i, j;
    2211:   IS              dofMap;
    2212:   IS              cellIS;
    2213:   const PetscInt  ndof = patch->totalDofsPerCell;
    2214:   Mat             vecMat;
    2215:   PetscInt        cStart, cEnd;
    2216:   DM              dm, plex;

    2218:   PetscCall(ISGetSize(patch->cells, &ncell));
    2219:   if (!ncell) { /* No cells to assemble over -> skip */
    2220:     PetscFunctionReturn(PETSC_SUCCESS);
    2221:   }

    2223:   PetscCall(PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0));

    2225:   PetscCall(PCGetDM(pc, &dm));
    2226:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    2227:   dm = plex;
    2228:   if (!patch->allCells) {
    2229:     PetscHSetI    cells;
    2230:     PetscHashIter hi;
    2231:     PetscInt      pStart, pEnd;
    2232:     PetscInt     *allCells = NULL;
    2233:     PetscCall(PetscHSetICreate(&cells));
    2234:     PetscCall(ISGetIndices(patch->cells, &cellsArray));
    2235:     PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));
    2236:     for (i = pStart; i < pEnd; i++) {
    2237:       PetscCall(PetscSectionGetDof(patch->cellCounts, i, &ncell));
    2238:       PetscCall(PetscSectionGetOffset(patch->cellCounts, i, &offset));
    2239:       if (ncell <= 0) continue;
    2240:       for (j = 0; j < ncell; j++) PetscCall(PetscHSetIAdd(cells, cellsArray[offset + j]));
    2241:     }
    2242:     PetscCall(ISRestoreIndices(patch->cells, &cellsArray));
    2243:     PetscCall(PetscHSetIGetSize(cells, &ncell));
    2244:     PetscCall(PetscMalloc1(ncell, &allCells));
    2245:     PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
    2246:     PetscCall(PetscMalloc1(cEnd - cStart, &patch->precomputedTensorLocations));
    2247:     i = 0;
    2248:     PetscHashIterBegin(cells, hi);
    2249:     while (!PetscHashIterAtEnd(cells, hi)) {
    2250:       PetscHashIterGetKey(cells, hi, allCells[i]);
    2251:       patch->precomputedTensorLocations[allCells[i]] = i;
    2252:       PetscHashIterNext(cells, hi);
    2253:       i++;
    2254:     }
    2255:     PetscCall(PetscHSetIDestroy(&cells));
    2256:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncell, allCells, PETSC_OWN_POINTER, &patch->allCells));
    2257:   }
    2258:   PetscCall(ISGetSize(patch->allCells, &ncell));
    2259:   if (!patch->cellMats) {
    2260:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, ncell * ndof * ndof, &patch->cellMats));
    2261:     PetscCall(VecSetBlockSize(patch->cellMats, ndof));
    2262:   }
    2263:   PetscCall(VecSet(patch->cellMats, 0));

    2265:   PetscCall(MatCreateShell(PETSC_COMM_SELF, ncell * ndof, ncell * ndof, ncell * ndof, ncell * ndof, (void *)patch->cellMats, &vecMat));
    2266:   PetscCall(MatShellSetOperation(vecMat, MATOP_SET_VALUES, (PetscErrorCodeFn *)MatSetValues_PCPatch_Private));
    2267:   PetscCall(ISGetSize(patch->allCells, &ncell));
    2268:   PetscCall(ISCreateStride(PETSC_COMM_SELF, ndof * ncell, 0, 1, &dofMap));
    2269:   PetscCall(ISGetIndices(dofMap, &dofMapArray));
    2270:   PetscCall(ISGetIndices(patch->allCells, &cellsArray));
    2271:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray, PETSC_USE_POINTER, &cellIS));
    2272:   /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
    2273:   PetscCallBack("PCPatch callback", patch->usercomputeop(pc, -1, NULL, vecMat, cellIS, ndof * ncell, dofMapArray, NULL, patch->usercomputeopctx));
    2274:   PetscCall(ISDestroy(&cellIS));
    2275:   PetscCall(MatDestroy(&vecMat));
    2276:   PetscCall(ISRestoreIndices(patch->allCells, &cellsArray));
    2277:   PetscCall(ISRestoreIndices(dofMap, &dofMapArray));
    2278:   PetscCall(ISDestroy(&dofMap));

    2280:   if (patch->usercomputeopintfacet) {
    2281:     PetscInt        nIntFacets;
    2282:     IS              intFacetsIS;
    2283:     const PetscInt *intFacetsArray = NULL;
    2284:     if (!patch->allIntFacets) {
    2285:       PetscHSetI    facets;
    2286:       PetscHashIter hi;
    2287:       PetscInt      pStart, pEnd, fStart, fEnd;
    2288:       PetscInt     *allIntFacets = NULL;
    2289:       PetscCall(PetscHSetICreate(&facets));
    2290:       PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
    2291:       PetscCall(PetscSectionGetChart(patch->intFacetCounts, &pStart, &pEnd));
    2292:       PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
    2293:       for (i = pStart; i < pEnd; i++) {
    2294:         PetscCall(PetscSectionGetDof(patch->intFacetCounts, i, &nIntFacets));
    2295:         PetscCall(PetscSectionGetOffset(patch->intFacetCounts, i, &offset));
    2296:         if (nIntFacets <= 0) continue;
    2297:         for (j = 0; j < nIntFacets; j++) PetscCall(PetscHSetIAdd(facets, intFacetsArray[offset + j]));
    2298:       }
    2299:       PetscCall(ISRestoreIndices(patch->intFacets, &intFacetsArray));
    2300:       PetscCall(PetscHSetIGetSize(facets, &nIntFacets));
    2301:       PetscCall(PetscMalloc1(nIntFacets, &allIntFacets));
    2302:       PetscCall(PetscMalloc1(fEnd - fStart, &patch->precomputedIntFacetTensorLocations));
    2303:       i = 0;
    2304:       PetscHashIterBegin(facets, hi);
    2305:       while (!PetscHashIterAtEnd(facets, hi)) {
    2306:         PetscHashIterGetKey(facets, hi, allIntFacets[i]);
    2307:         patch->precomputedIntFacetTensorLocations[allIntFacets[i] - fStart] = i;
    2308:         PetscHashIterNext(facets, hi);
    2309:         i++;
    2310:       }
    2311:       PetscCall(PetscHSetIDestroy(&facets));
    2312:       PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, allIntFacets, PETSC_OWN_POINTER, &patch->allIntFacets));
    2313:     }
    2314:     PetscCall(ISGetSize(patch->allIntFacets, &nIntFacets));
    2315:     if (!patch->intFacetMats) {
    2316:       PetscCall(VecCreateSeq(PETSC_COMM_SELF, nIntFacets * ndof * ndof * 4, &patch->intFacetMats));
    2317:       PetscCall(VecSetBlockSize(patch->intFacetMats, ndof * 2));
    2318:     }
    2319:     PetscCall(VecSet(patch->intFacetMats, 0));

    2321:     PetscCall(MatCreateShell(PETSC_COMM_SELF, nIntFacets * ndof * 2, nIntFacets * ndof * 2, nIntFacets * ndof * 2, nIntFacets * ndof * 2, (void *)patch->intFacetMats, &vecMat));
    2322:     PetscCall(MatShellSetOperation(vecMat, MATOP_SET_VALUES, (PetscErrorCodeFn *)MatSetValues_PCPatch_Private));
    2323:     PetscCall(ISCreateStride(PETSC_COMM_SELF, 2 * ndof * nIntFacets, 0, 1, &dofMap));
    2324:     PetscCall(ISGetIndices(dofMap, &dofMapArray));
    2325:     PetscCall(ISGetIndices(patch->allIntFacets, &intFacetsArray));
    2326:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, intFacetsArray, PETSC_USE_POINTER, &intFacetsIS));
    2327:     /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
    2328:     PetscCallBack("PCPatch callback (interior facets)", patch->usercomputeopintfacet(pc, -1, NULL, vecMat, intFacetsIS, 2 * ndof * nIntFacets, dofMapArray, NULL, patch->usercomputeopintfacetctx));
    2329:     PetscCall(ISDestroy(&intFacetsIS));
    2330:     PetscCall(MatDestroy(&vecMat));
    2331:     PetscCall(ISRestoreIndices(patch->allIntFacets, &intFacetsArray));
    2332:     PetscCall(ISRestoreIndices(dofMap, &dofMapArray));
    2333:     PetscCall(ISDestroy(&dofMap));
    2334:   }
    2335:   PetscCall(DMDestroy(&dm));
    2336:   PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
    2337:   PetscFunctionReturn(PETSC_SUCCESS);
    2338: }

    2340: PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat, PatchScatterType scattertype)
    2341: {
    2342:   PC_PATCH          *patch     = (PC_PATCH *)pc->data;
    2343:   const PetscScalar *xArray    = NULL;
    2344:   PetscScalar       *yArray    = NULL;
    2345:   const PetscInt    *gtolArray = NULL;
    2346:   PetscInt           dof, offset, lidx;

    2348:   PetscFunctionBeginHot;
    2349:   PetscCall(VecGetArrayRead(x, &xArray));
    2350:   PetscCall(VecGetArray(y, &yArray));
    2351:   if (scattertype == SCATTER_WITHARTIFICIAL) {
    2352:     PetscCall(PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof));
    2353:     PetscCall(PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offset));
    2354:     PetscCall(ISGetIndices(patch->gtolWithArtificial, &gtolArray));
    2355:   } else if (scattertype == SCATTER_WITHALL) {
    2356:     PetscCall(PetscSectionGetDof(patch->gtolCountsWithAll, p, &dof));
    2357:     PetscCall(PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offset));
    2358:     PetscCall(ISGetIndices(patch->gtolWithAll, &gtolArray));
    2359:   } else {
    2360:     PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &dof));
    2361:     PetscCall(PetscSectionGetOffset(patch->gtolCounts, p, &offset));
    2362:     PetscCall(ISGetIndices(patch->gtol, &gtolArray));
    2363:   }
    2364:   PetscCheck(mode != INSERT_VALUES || scat == SCATTER_FORWARD, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't insert if not scattering forward");
    2365:   PetscCheck(mode != ADD_VALUES || scat == SCATTER_REVERSE, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't add if not scattering reverse");
    2366:   for (lidx = 0; lidx < dof; ++lidx) {
    2367:     const PetscInt gidx = gtolArray[offset + lidx];

    2369:     if (mode == INSERT_VALUES) yArray[lidx] = xArray[gidx]; /* Forward */
    2370:     else yArray[gidx] += xArray[lidx];                      /* Reverse */
    2371:   }
    2372:   if (scattertype == SCATTER_WITHARTIFICIAL) {
    2373:     PetscCall(ISRestoreIndices(patch->gtolWithArtificial, &gtolArray));
    2374:   } else if (scattertype == SCATTER_WITHALL) {
    2375:     PetscCall(ISRestoreIndices(patch->gtolWithAll, &gtolArray));
    2376:   } else {
    2377:     PetscCall(ISRestoreIndices(patch->gtol, &gtolArray));
    2378:   }
    2379:   PetscCall(VecRestoreArrayRead(x, &xArray));
    2380:   PetscCall(VecRestoreArray(y, &yArray));
    2381:   PetscFunctionReturn(PETSC_SUCCESS);
    2382: }

    2384: static PetscErrorCode PCSetUp_PATCH_Linear(PC pc)
    2385: {
    2386:   PC_PATCH   *patch = (PC_PATCH *)pc->data;
    2387:   const char *prefix;
    2388:   PetscInt    i;

    2390:   PetscFunctionBegin;
    2391:   if (!pc->setupcalled) {
    2392:     PetscCheck(patch->save_operators || !patch->denseinverse, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Can't have dense inverse without save operators");
    2393:     if (!patch->denseinverse) {
    2394:       PetscCall(PetscMalloc1(patch->npatch, &patch->solver));
    2395:       PetscCall(PCGetOptionsPrefix(pc, &prefix));
    2396:       for (i = 0; i < patch->npatch; ++i) {
    2397:         KSP ksp;
    2398:         PC  subpc;

    2400:         PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
    2401:         PetscCall(KSPSetNestLevel(ksp, pc->kspnestlevel));
    2402:         PetscCall(KSPSetErrorIfNotConverged(ksp, pc->erroriffailure));
    2403:         PetscCall(KSPSetOptionsPrefix(ksp, prefix));
    2404:         PetscCall(KSPAppendOptionsPrefix(ksp, "sub_"));
    2405:         PetscCall(PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1));
    2406:         PetscCall(KSPGetPC(ksp, &subpc));
    2407:         PetscCall(PetscObjectIncrementTabLevel((PetscObject)subpc, (PetscObject)pc, 1));
    2408:         patch->solver[i] = (PetscObject)ksp;
    2409:       }
    2410:     }
    2411:   }
    2412:   if (patch->save_operators) {
    2413:     if (patch->precomputeElementTensors) PetscCall(PCPatchPrecomputePatchTensors_Private(pc));
    2414:     for (i = 0; i < patch->npatch; ++i) {
    2415:       PetscCall(PCPatchComputeOperator_Internal(pc, NULL, patch->mat[i], i, PETSC_FALSE));
    2416:       if (!patch->denseinverse) {
    2417:         PetscCall(KSPSetOperators((KSP)patch->solver[i], patch->mat[i], patch->mat[i]));
    2418:       } else if (patch->mat[i] && !patch->densesolve) {
    2419:         /* Setup matmult callback */
    2420:         PetscCall(MatGetOperation(patch->mat[i], MATOP_MULT, (PetscErrorCodeFn **)&patch->densesolve));
    2421:       }
    2422:     }
    2423:   }
    2424:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
    2425:     for (i = 0; i < patch->npatch; ++i) {
    2426:       /* Instead of padding patch->patchUpdate with zeros to get */
    2427:       /* patch->patchUpdateWithArtificial and then multiplying with the matrix, */
    2428:       /* just get rid of the columns that correspond to the dofs with */
    2429:       /* artificial bcs. That's of course fairly inefficient, hopefully we */
    2430:       /* can just assemble the rectangular matrix in the first place. */
    2431:       Mat      matSquare;
    2432:       IS       rowis;
    2433:       PetscInt dof;

    2435:       PetscCall(MatGetSize(patch->mat[i], &dof, NULL));
    2436:       if (dof == 0) {
    2437:         patch->matWithArtificial[i] = NULL;
    2438:         continue;
    2439:       }

    2441:       PetscCall(PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE));
    2442:       PetscCall(PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE));

    2444:       PetscCall(MatGetSize(matSquare, &dof, NULL));
    2445:       PetscCall(ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis));
    2446:       if (pc->setupcalled) {
    2447:         PetscCall(MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_REUSE_MATRIX, &patch->matWithArtificial[i]));
    2448:       } else {
    2449:         PetscCall(MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &patch->matWithArtificial[i]));
    2450:       }
    2451:       PetscCall(ISDestroy(&rowis));
    2452:       PetscCall(MatDestroy(&matSquare));
    2453:     }
    2454:   }
    2455:   PetscFunctionReturn(PETSC_SUCCESS);
    2456: }

    2458: static PetscErrorCode PCSetUp_PATCH(PC pc)
    2459: {
    2460:   PC_PATCH *patch = (PC_PATCH *)pc->data;
    2461:   PetscInt  i;
    2462:   PetscBool isNonlinear;
    2463:   PetscInt  maxDof = -1, maxDofWithArtificial = -1;

    2465:   PetscFunctionBegin;
    2466:   if (!pc->setupcalled) {
    2467:     PetscInt pStart, pEnd, p;
    2468:     PetscInt localSize;

    2470:     PetscCall(PetscLogEventBegin(PC_Patch_CreatePatches, pc, 0, 0, 0));

    2472:     isNonlinear = patch->isNonlinear;
    2473:     if (!patch->nsubspaces) {
    2474:       DM           dm, plex;
    2475:       PetscSection s;
    2476:       PetscInt     cStart, cEnd, c, Nf, f, numGlobalBcs = 0, *globalBcs, *Nb, **cellDofs;

    2478:       PetscCall(PCGetDM(pc, &dm));
    2479:       PetscCheck(dm, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Must set DM for PCPATCH or call PCPatchSetDiscretisationInfo()");
    2480:       PetscCall(DMConvert(dm, DMPLEX, &plex));
    2481:       dm = plex;
    2482:       PetscCall(DMGetLocalSection(dm, &s));
    2483:       PetscCall(PetscSectionGetNumFields(s, &Nf));
    2484:       PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
    2485:       for (p = pStart; p < pEnd; ++p) {
    2486:         PetscInt cdof;
    2487:         PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
    2488:         numGlobalBcs += cdof;
    2489:       }
    2490:       PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
    2491:       PetscCall(PetscMalloc3(Nf, &Nb, Nf, &cellDofs, numGlobalBcs, &globalBcs));
    2492:       for (f = 0; f < Nf; ++f) {
    2493:         PetscFE        fe;
    2494:         PetscDualSpace sp;
    2495:         PetscInt       cdoff = 0;

    2497:         PetscCall(DMGetField(dm, f, NULL, (PetscObject *)&fe));
    2498:         /* PetscCall(PetscFEGetNumComponents(fe, &Nc[f])); */
    2499:         PetscCall(PetscFEGetDualSpace(fe, &sp));
    2500:         PetscCall(PetscDualSpaceGetDimension(sp, &Nb[f]));

    2502:         PetscCall(PetscMalloc1((cEnd - cStart) * Nb[f], &cellDofs[f]));
    2503:         for (c = cStart; c < cEnd; ++c) {
    2504:           PetscInt *closure = NULL;
    2505:           PetscInt  clSize  = 0, cl;

    2507:           PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
    2508:           for (cl = 0; cl < clSize * 2; cl += 2) {
    2509:             const PetscInt p = closure[cl];
    2510:             PetscInt       fdof, d, foff;

    2512:             PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
    2513:             PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
    2514:             for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d;
    2515:           }
    2516:           PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
    2517:         }
    2518:         PetscCheck(cdoff == (cEnd - cStart) * Nb[f], PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_SIZ, "Total number of cellDofs %" PetscInt_FMT " for field %" PetscInt_FMT " should be Nc (%" PetscInt_FMT ") * cellDof (%" PetscInt_FMT ")", cdoff, f, cEnd - cStart, Nb[f]);
    2519:       }
    2520:       numGlobalBcs = 0;
    2521:       for (p = pStart; p < pEnd; ++p) {
    2522:         const PetscInt *ind;
    2523:         PetscInt        off, cdof, d;

    2525:         PetscCall(PetscSectionGetOffset(s, p, &off));
    2526:         PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
    2527:         PetscCall(PetscSectionGetConstraintIndices(s, p, &ind));
    2528:         for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d];
    2529:       }

    2531:       PetscCall(PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **)cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs));
    2532:       for (f = 0; f < Nf; ++f) PetscCall(PetscFree(cellDofs[f]));
    2533:       PetscCall(PetscFree3(Nb, cellDofs, globalBcs));
    2534:       PetscCall(PCPatchSetComputeFunction(pc, PCPatchComputeFunction_DMPlex_Private, NULL));
    2535:       PetscCall(PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL));
    2536:       PetscCall(DMDestroy(&dm));
    2537:     }

    2539:     localSize = patch->subspaceOffsets[patch->nsubspaces];
    2540:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localRHS));
    2541:     PetscCall(VecSetUp(patch->localRHS));
    2542:     PetscCall(VecDuplicate(patch->localRHS, &patch->localUpdate));
    2543:     PetscCall(PCPatchCreateCellPatches(pc));
    2544:     PetscCall(PCPatchCreateCellPatchDiscretisationInfo(pc));

    2546:     /* OK, now build the work vectors */
    2547:     PetscCall(PetscSectionGetChart(patch->gtolCounts, &pStart, &pEnd));

    2549:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscCall(PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithArtificial));
    2550:     if (isNonlinear) PetscCall(PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithAll));
    2551:     for (p = pStart; p < pEnd; ++p) {
    2552:       PetscInt dof;

    2554:       PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &dof));
    2555:       maxDof = PetscMax(maxDof, dof);
    2556:       if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
    2557:         const PetscInt *gtolArray, *gtolArrayWithArtificial = NULL;
    2558:         PetscInt        numPatchDofs, offset;
    2559:         PetscInt        numPatchDofsWithArtificial, offsetWithArtificial;
    2560:         PetscInt        dofWithoutArtificialCounter = 0;
    2561:         PetscInt       *patchWithoutArtificialToWithArtificialArray;

    2563:         PetscCall(PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof));
    2564:         maxDofWithArtificial = PetscMax(maxDofWithArtificial, dof);

    2566:         /* Now build the mapping that for a dof in a patch WITHOUT dofs that have artificial bcs gives the */
    2567:         /* the index in the patch with all dofs */
    2568:         PetscCall(ISGetIndices(patch->gtol, &gtolArray));

    2570:         PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs));
    2571:         if (numPatchDofs == 0) {
    2572:           patch->dofMappingWithoutToWithArtificial[p - pStart] = NULL;
    2573:           continue;
    2574:         }

    2576:         PetscCall(PetscSectionGetOffset(patch->gtolCounts, p, &offset));
    2577:         PetscCall(ISGetIndices(patch->gtolWithArtificial, &gtolArrayWithArtificial));
    2578:         PetscCall(PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &numPatchDofsWithArtificial));
    2579:         PetscCall(PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offsetWithArtificial));

    2581:         PetscCall(PetscMalloc1(numPatchDofs, &patchWithoutArtificialToWithArtificialArray));
    2582:         for (i = 0; i < numPatchDofsWithArtificial; i++) {
    2583:           if (gtolArrayWithArtificial[i + offsetWithArtificial] == gtolArray[offset + dofWithoutArtificialCounter]) {
    2584:             patchWithoutArtificialToWithArtificialArray[dofWithoutArtificialCounter] = i;
    2585:             dofWithoutArtificialCounter++;
    2586:             if (dofWithoutArtificialCounter == numPatchDofs) break;
    2587:           }
    2588:         }
    2589:         PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutArtificialToWithArtificialArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithArtificial[p - pStart]));
    2590:         PetscCall(ISRestoreIndices(patch->gtol, &gtolArray));
    2591:         PetscCall(ISRestoreIndices(patch->gtolWithArtificial, &gtolArrayWithArtificial));
    2592:       }
    2593:     }
    2594:     for (p = pStart; p < pEnd; ++p) {
    2595:       if (isNonlinear) {
    2596:         const PetscInt *gtolArray, *gtolArrayWithAll = NULL;
    2597:         PetscInt        numPatchDofs, offset;
    2598:         PetscInt        numPatchDofsWithAll, offsetWithAll;
    2599:         PetscInt        dofWithoutAllCounter = 0;
    2600:         PetscInt       *patchWithoutAllToWithAllArray;

    2602:         /* Now build the mapping that for a dof in a patch WITHOUT dofs that have artificial bcs gives the */
    2603:         /* the index in the patch with all dofs */
    2604:         PetscCall(ISGetIndices(patch->gtol, &gtolArray));

    2606:         PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs));
    2607:         if (numPatchDofs == 0) {
    2608:           patch->dofMappingWithoutToWithAll[p - pStart] = NULL;
    2609:           continue;
    2610:         }

    2612:         PetscCall(PetscSectionGetOffset(patch->gtolCounts, p, &offset));
    2613:         PetscCall(ISGetIndices(patch->gtolWithAll, &gtolArrayWithAll));
    2614:         PetscCall(PetscSectionGetDof(patch->gtolCountsWithAll, p, &numPatchDofsWithAll));
    2615:         PetscCall(PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offsetWithAll));

    2617:         PetscCall(PetscMalloc1(numPatchDofs, &patchWithoutAllToWithAllArray));

    2619:         for (i = 0; i < numPatchDofsWithAll; i++) {
    2620:           if (gtolArrayWithAll[i + offsetWithAll] == gtolArray[offset + dofWithoutAllCounter]) {
    2621:             patchWithoutAllToWithAllArray[dofWithoutAllCounter] = i;
    2622:             dofWithoutAllCounter++;
    2623:             if (dofWithoutAllCounter == numPatchDofs) break;
    2624:           }
    2625:         }
    2626:         PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutAllToWithAllArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithAll[p - pStart]));
    2627:         PetscCall(ISRestoreIndices(patch->gtol, &gtolArray));
    2628:         PetscCall(ISRestoreIndices(patch->gtolWithAll, &gtolArrayWithAll));
    2629:       }
    2630:     }
    2631:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
    2632:       PetscCall(VecCreateSeq(PETSC_COMM_SELF, maxDofWithArtificial, &patch->patchRHSWithArtificial));
    2633:       PetscCall(VecSetUp(patch->patchRHSWithArtificial));
    2634:     }
    2635:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchRHS));
    2636:     PetscCall(VecSetUp(patch->patchRHS));
    2637:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchUpdate));
    2638:     PetscCall(VecSetUp(patch->patchUpdate));
    2639:     if (patch->save_operators) {
    2640:       PetscCall(PetscMalloc1(patch->npatch, &patch->mat));
    2641:       for (i = 0; i < patch->npatch; ++i) PetscCall(PCPatchCreateMatrix_Private(pc, i, &patch->mat[i], PETSC_FALSE));
    2642:     }
    2643:     PetscCall(PetscLogEventEnd(PC_Patch_CreatePatches, pc, 0, 0, 0));

    2645:     /* If desired, calculate weights for dof multiplicity */
    2646:     if (patch->partition_of_unity) {
    2647:       PetscScalar *input  = NULL;
    2648:       PetscScalar *output = NULL;
    2649:       Vec          global;

    2651:       PetscCall(VecDuplicate(patch->localRHS, &patch->dof_weights));
    2652:       if (patch->local_composition_type == PC_COMPOSITE_ADDITIVE) {
    2653:         for (i = 0; i < patch->npatch; ++i) {
    2654:           PetscInt dof;

    2656:           PetscCall(PetscSectionGetDof(patch->gtolCounts, i + pStart, &dof));
    2657:           if (dof <= 0) continue;
    2658:           PetscCall(VecSet(patch->patchRHS, 1.0));
    2659:           PetscCall(PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHS, patch->dof_weights, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR));
    2660:         }
    2661:       } else {
    2662:         /* multiplicative is actually only locally multiplicative and globally additive. need the pou where the mesh decomposition overlaps */
    2663:         PetscCall(VecSet(patch->dof_weights, 1.0));
    2664:       }

    2666:       PetscCall(VecDuplicate(patch->dof_weights, &global));
    2667:       PetscCall(VecSet(global, 0.));

    2669:       PetscCall(VecGetArray(patch->dof_weights, &input));
    2670:       PetscCall(VecGetArray(global, &output));
    2671:       PetscCall(PetscSFReduceBegin(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM));
    2672:       PetscCall(PetscSFReduceEnd(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM));
    2673:       PetscCall(VecRestoreArray(patch->dof_weights, &input));
    2674:       PetscCall(VecRestoreArray(global, &output));

    2676:       PetscCall(VecReciprocal(global));

    2678:       PetscCall(VecGetArray(patch->dof_weights, &output));
    2679:       PetscCall(VecGetArray(global, &input));
    2680:       PetscCall(PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, input, output, MPI_REPLACE));
    2681:       PetscCall(PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, input, output, MPI_REPLACE));
    2682:       PetscCall(VecRestoreArray(patch->dof_weights, &output));
    2683:       PetscCall(VecRestoreArray(global, &input));
    2684:       PetscCall(VecDestroy(&global));
    2685:     }
    2686:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE && patch->save_operators && !patch->isNonlinear) PetscCall(PetscMalloc1(patch->npatch, &patch->matWithArtificial));
    2687:   }
    2688:   PetscCall((*patch->setupsolver)(pc));
    2689:   PetscFunctionReturn(PETSC_SUCCESS);
    2690: }

    2692: static PetscErrorCode PCApply_PATCH_Linear(PC pc, PetscInt i, Vec x, Vec y)
    2693: {
    2694:   PC_PATCH *patch = (PC_PATCH *)pc->data;
    2695:   KSP       ksp;
    2696:   Mat       op;
    2697:   PetscInt  m, n;

    2699:   PetscFunctionBegin;
    2700:   if (patch->denseinverse) {
    2701:     PetscCall((*patch->densesolve)(patch->mat[i], x, y));
    2702:     PetscFunctionReturn(PETSC_SUCCESS);
    2703:   }
    2704:   ksp = (KSP)patch->solver[i];
    2705:   if (!patch->save_operators) {
    2706:     Mat mat;

    2708:     PetscCall(PCPatchCreateMatrix_Private(pc, i, &mat, PETSC_FALSE));
    2709:     /* Populate operator here. */
    2710:     PetscCall(PCPatchComputeOperator_Internal(pc, NULL, mat, i, PETSC_FALSE));
    2711:     PetscCall(KSPSetOperators(ksp, mat, mat));
    2712:     /* Drop reference so the KSPSetOperators below will blow it away. */
    2713:     PetscCall(MatDestroy(&mat));
    2714:   }
    2715:   PetscCall(PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0));
    2716:   if (!ksp->setfromoptionscalled) PetscCall(KSPSetFromOptions(ksp));
    2717:   /* Disgusting trick to reuse work vectors */
    2718:   PetscCall(KSPGetOperators(ksp, &op, NULL));
    2719:   PetscCall(MatGetLocalSize(op, &m, &n));
    2720:   x->map->n = m;
    2721:   y->map->n = n;
    2722:   x->map->N = m;
    2723:   y->map->N = n;
    2724:   PetscCall(KSPSolve(ksp, x, y));
    2725:   PetscCall(KSPCheckSolve(ksp, pc, y));
    2726:   PetscCall(PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0));
    2727:   if (!patch->save_operators) {
    2728:     PC pc;
    2729:     PetscCall(KSPSetOperators(ksp, NULL, NULL));
    2730:     PetscCall(KSPGetPC(ksp, &pc));
    2731:     /* Destroy PC context too, otherwise the factored matrix hangs around. */
    2732:     PetscCall(PCReset(pc));
    2733:   }
    2734:   PetscFunctionReturn(PETSC_SUCCESS);
    2735: }

    2737: static PetscErrorCode PCUpdateMultiplicative_PATCH_Linear(PC pc, PetscInt i, PetscInt pStart)
    2738: {
    2739:   PC_PATCH *patch = (PC_PATCH *)pc->data;
    2740:   Mat       multMat;
    2741:   PetscInt  n, m;

    2743:   PetscFunctionBegin;
    2744:   if (patch->save_operators) {
    2745:     multMat = patch->matWithArtificial[i];
    2746:   } else {
    2747:     /*Very inefficient, hopefully we can just assemble the rectangular matrix in the first place.*/
    2748:     Mat      matSquare;
    2749:     PetscInt dof;
    2750:     IS       rowis;
    2751:     PetscCall(PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE));
    2752:     PetscCall(PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE));
    2753:     PetscCall(MatGetSize(matSquare, &dof, NULL));
    2754:     PetscCall(ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis));
    2755:     PetscCall(MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &multMat));
    2756:     PetscCall(MatDestroy(&matSquare));
    2757:     PetscCall(ISDestroy(&rowis));
    2758:   }
    2759:   /* Disgusting trick to reuse work vectors */
    2760:   PetscCall(MatGetLocalSize(multMat, &m, &n));
    2761:   patch->patchUpdate->map->n            = n;
    2762:   patch->patchRHSWithArtificial->map->n = m;
    2763:   patch->patchUpdate->map->N            = n;
    2764:   patch->patchRHSWithArtificial->map->N = m;
    2765:   PetscCall(MatMult(multMat, patch->patchUpdate, patch->patchRHSWithArtificial));
    2766:   PetscCall(VecScale(patch->patchRHSWithArtificial, -1.0));
    2767:   PetscCall(PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHSWithArtificial, patch->localRHS, ADD_VALUES, SCATTER_REVERSE, SCATTER_WITHARTIFICIAL));
    2768:   if (!patch->save_operators) PetscCall(MatDestroy(&multMat));
    2769:   PetscFunctionReturn(PETSC_SUCCESS);
    2770: }

    2772: static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y)
    2773: {
    2774:   PC_PATCH          *patch        = (PC_PATCH *)pc->data;
    2775:   const PetscScalar *globalRHS    = NULL;
    2776:   PetscScalar       *localRHS     = NULL;
    2777:   PetscScalar       *globalUpdate = NULL;
    2778:   const PetscInt    *bcNodes      = NULL;
    2779:   PetscInt           nsweep       = patch->symmetrise_sweep ? 2 : 1;
    2780:   PetscInt           start[2]     = {0, 0};
    2781:   PetscInt           end[2]       = {-1, -1};
    2782:   const PetscInt     inc[2]       = {1, -1};
    2783:   const PetscScalar *localUpdate;
    2784:   const PetscInt    *iterationSet;
    2785:   PetscInt           pStart, numBcs, n, sweep, bc, j;

    2787:   PetscFunctionBegin;
    2788:   PetscCall(PetscLogEventBegin(PC_Patch_Apply, pc, 0, 0, 0));
    2789:   PetscCall(PetscOptionsPushCreateViewerOff(PETSC_TRUE));
    2790:   /* start, end, inc have 2 entries to manage a second backward sweep if we symmetrize */
    2791:   end[0]   = patch->npatch;
    2792:   start[1] = patch->npatch - 1;
    2793:   if (patch->user_patches) {
    2794:     PetscCall(ISGetLocalSize(patch->iterationSet, &end[0]));
    2795:     start[1] = end[0] - 1;
    2796:     PetscCall(ISGetIndices(patch->iterationSet, &iterationSet));
    2797:   }
    2798:   /* Scatter from global space into overlapped local spaces */
    2799:   PetscCall(VecGetArrayRead(x, &globalRHS));
    2800:   PetscCall(VecGetArray(patch->localRHS, &localRHS));
    2801:   PetscCall(PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, globalRHS, localRHS, MPI_REPLACE));
    2802:   PetscCall(PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, globalRHS, localRHS, MPI_REPLACE));
    2803:   PetscCall(VecRestoreArrayRead(x, &globalRHS));
    2804:   PetscCall(VecRestoreArray(patch->localRHS, &localRHS));

    2806:   PetscCall(VecSet(patch->localUpdate, 0.0));
    2807:   PetscCall(PetscSectionGetChart(patch->gtolCounts, &pStart, NULL));
    2808:   PetscCall(PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0));
    2809:   for (sweep = 0; sweep < nsweep; sweep++) {
    2810:     for (j = start[sweep]; j * inc[sweep] < end[sweep] * inc[sweep]; j += inc[sweep]) {
    2811:       PetscInt i = patch->user_patches ? iterationSet[j] : j;
    2812:       PetscInt start, len;

    2814:       PetscCall(PetscSectionGetDof(patch->gtolCounts, i + pStart, &len));
    2815:       PetscCall(PetscSectionGetOffset(patch->gtolCounts, i + pStart, &start));
    2816:       /* TODO: Squash out these guys in the setup as well. */
    2817:       if (len <= 0) continue;
    2818:       /* TODO: Do we need different scatters for X and Y? */
    2819:       PetscCall(PCPatch_ScatterLocal_Private(pc, i + pStart, patch->localRHS, patch->patchRHS, INSERT_VALUES, SCATTER_FORWARD, SCATTER_INTERIOR));
    2820:       PetscCall((*patch->applysolver)(pc, i, patch->patchRHS, patch->patchUpdate));
    2821:       PetscCall(PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchUpdate, patch->localUpdate, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR));
    2822:       if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscCall((*patch->updatemultiplicative)(pc, i, pStart));
    2823:     }
    2824:   }
    2825:   PetscCall(PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0));
    2826:   if (patch->user_patches) PetscCall(ISRestoreIndices(patch->iterationSet, &iterationSet));
    2827:   /* XXX: should we do this on the global vector? */
    2828:   if (patch->partition_of_unity) PetscCall(VecPointwiseMult(patch->localUpdate, patch->localUpdate, patch->dof_weights));
    2829:   /* Now patch->localUpdate contains the solution of the patch solves, so we need to combine them all. */
    2830:   PetscCall(VecSet(y, 0.0));
    2831:   PetscCall(VecGetArray(y, &globalUpdate));
    2832:   PetscCall(VecGetArrayRead(patch->localUpdate, &localUpdate));
    2833:   PetscCall(PetscSFReduceBegin(patch->sectionSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM));
    2834:   PetscCall(PetscSFReduceEnd(patch->sectionSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM));
    2835:   PetscCall(VecRestoreArrayRead(patch->localUpdate, &localUpdate));

    2837:   /* Now we need to send the global BC values through */
    2838:   PetscCall(VecGetArrayRead(x, &globalRHS));
    2839:   PetscCall(ISGetSize(patch->globalBcNodes, &numBcs));
    2840:   PetscCall(ISGetIndices(patch->globalBcNodes, &bcNodes));
    2841:   PetscCall(VecGetLocalSize(x, &n));
    2842:   for (bc = 0; bc < numBcs; ++bc) {
    2843:     const PetscInt idx = bcNodes[bc];
    2844:     if (idx < n) globalUpdate[idx] = globalRHS[idx];
    2845:   }

    2847:   PetscCall(ISRestoreIndices(patch->globalBcNodes, &bcNodes));
    2848:   PetscCall(VecRestoreArrayRead(x, &globalRHS));
    2849:   PetscCall(VecRestoreArray(y, &globalUpdate));

    2851:   PetscCall(PetscOptionsPopCreateViewerOff());
    2852:   PetscCall(PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0));
    2853:   PetscFunctionReturn(PETSC_SUCCESS);
    2854: }

    2856: static PetscErrorCode PCReset_PATCH_Linear(PC pc)
    2857: {
    2858:   PC_PATCH *patch = (PC_PATCH *)pc->data;
    2859:   PetscInt  i;

    2861:   PetscFunctionBegin;
    2862:   if (patch->solver) {
    2863:     for (i = 0; i < patch->npatch; ++i) PetscCall(KSPReset((KSP)patch->solver[i]));
    2864:   }
    2865:   PetscFunctionReturn(PETSC_SUCCESS);
    2866: }

    2868: static PetscErrorCode PCReset_PATCH(PC pc)
    2869: {
    2870:   PC_PATCH *patch = (PC_PATCH *)pc->data;
    2871:   PetscInt  i;

    2873:   PetscFunctionBegin;
    2874:   PetscCall(PetscSFDestroy(&patch->sectionSF));
    2875:   PetscCall(PetscSectionDestroy(&patch->cellCounts));
    2876:   PetscCall(PetscSectionDestroy(&patch->pointCounts));
    2877:   PetscCall(PetscSectionDestroy(&patch->cellNumbering));
    2878:   PetscCall(PetscSectionDestroy(&patch->gtolCounts));
    2879:   PetscCall(ISDestroy(&patch->gtol));
    2880:   PetscCall(ISDestroy(&patch->cells));
    2881:   PetscCall(ISDestroy(&patch->points));
    2882:   PetscCall(ISDestroy(&patch->dofs));
    2883:   PetscCall(ISDestroy(&patch->offs));
    2884:   PetscCall(PetscSectionDestroy(&patch->patchSection));
    2885:   PetscCall(ISDestroy(&patch->ghostBcNodes));
    2886:   PetscCall(ISDestroy(&patch->globalBcNodes));
    2887:   PetscCall(PetscSectionDestroy(&patch->gtolCountsWithArtificial));
    2888:   PetscCall(ISDestroy(&patch->gtolWithArtificial));
    2889:   PetscCall(ISDestroy(&patch->dofsWithArtificial));
    2890:   PetscCall(ISDestroy(&patch->offsWithArtificial));
    2891:   PetscCall(PetscSectionDestroy(&patch->gtolCountsWithAll));
    2892:   PetscCall(ISDestroy(&patch->gtolWithAll));
    2893:   PetscCall(ISDestroy(&patch->dofsWithAll));
    2894:   PetscCall(ISDestroy(&patch->offsWithAll));
    2895:   PetscCall(VecDestroy(&patch->cellMats));
    2896:   PetscCall(VecDestroy(&patch->intFacetMats));
    2897:   PetscCall(ISDestroy(&patch->allCells));
    2898:   PetscCall(ISDestroy(&patch->intFacets));
    2899:   PetscCall(ISDestroy(&patch->extFacets));
    2900:   PetscCall(ISDestroy(&patch->intFacetsToPatchCell));
    2901:   PetscCall(PetscSectionDestroy(&patch->intFacetCounts));
    2902:   PetscCall(PetscSectionDestroy(&patch->extFacetCounts));

    2904:   if (patch->dofSection)
    2905:     for (i = 0; i < patch->nsubspaces; i++) PetscCall(PetscSectionDestroy(&patch->dofSection[i]));
    2906:   PetscCall(PetscFree(patch->dofSection));
    2907:   PetscCall(PetscFree(patch->bs));
    2908:   PetscCall(PetscFree(patch->nodesPerCell));
    2909:   if (patch->cellNodeMap)
    2910:     for (i = 0; i < patch->nsubspaces; i++) PetscCall(PetscFree(patch->cellNodeMap[i]));
    2911:   PetscCall(PetscFree(patch->cellNodeMap));
    2912:   PetscCall(PetscFree(patch->subspaceOffsets));

    2914:   PetscCall((*patch->resetsolver)(pc));

    2916:   if (patch->subspaces_to_exclude) PetscCall(PetscHSetIDestroy(&patch->subspaces_to_exclude));

    2918:   PetscCall(VecDestroy(&patch->localRHS));
    2919:   PetscCall(VecDestroy(&patch->localUpdate));
    2920:   PetscCall(VecDestroy(&patch->patchRHS));
    2921:   PetscCall(VecDestroy(&patch->patchUpdate));
    2922:   PetscCall(VecDestroy(&patch->dof_weights));
    2923:   if (patch->patch_dof_weights) {
    2924:     for (i = 0; i < patch->npatch; ++i) PetscCall(VecDestroy(&patch->patch_dof_weights[i]));
    2925:     PetscCall(PetscFree(patch->patch_dof_weights));
    2926:   }
    2927:   if (patch->mat) {
    2928:     for (i = 0; i < patch->npatch; ++i) PetscCall(MatDestroy(&patch->mat[i]));
    2929:     PetscCall(PetscFree(patch->mat));
    2930:   }
    2931:   if (patch->matWithArtificial && !patch->isNonlinear) {
    2932:     for (i = 0; i < patch->npatch; ++i) PetscCall(MatDestroy(&patch->matWithArtificial[i]));
    2933:     PetscCall(PetscFree(patch->matWithArtificial));
    2934:   }
    2935:   PetscCall(VecDestroy(&patch->patchRHSWithArtificial));
    2936:   if (patch->dofMappingWithoutToWithArtificial) {
    2937:     for (i = 0; i < patch->npatch; ++i) PetscCall(ISDestroy(&patch->dofMappingWithoutToWithArtificial[i]));
    2938:     PetscCall(PetscFree(patch->dofMappingWithoutToWithArtificial));
    2939:   }
    2940:   if (patch->dofMappingWithoutToWithAll) {
    2941:     for (i = 0; i < patch->npatch; ++i) PetscCall(ISDestroy(&patch->dofMappingWithoutToWithAll[i]));
    2942:     PetscCall(PetscFree(patch->dofMappingWithoutToWithAll));
    2943:   }
    2944:   PetscCall(PetscFree(patch->sub_mat_type));
    2945:   if (patch->userIS) {
    2946:     for (i = 0; i < patch->npatch; ++i) PetscCall(ISDestroy(&patch->userIS[i]));
    2947:     PetscCall(PetscFree(patch->userIS));
    2948:   }
    2949:   PetscCall(PetscFree(patch->precomputedTensorLocations));
    2950:   PetscCall(PetscFree(patch->precomputedIntFacetTensorLocations));

    2952:   patch->bs          = NULL;
    2953:   patch->cellNodeMap = NULL;
    2954:   patch->nsubspaces  = 0;
    2955:   PetscCall(ISDestroy(&patch->iterationSet));

    2957:   PetscCall(PetscViewerDestroy(&patch->viewerCells));
    2958:   PetscCall(PetscViewerDestroy(&patch->viewerIntFacets));
    2959:   PetscCall(PetscViewerDestroy(&patch->viewerPoints));
    2960:   PetscCall(PetscViewerDestroy(&patch->viewerSection));
    2961:   PetscCall(PetscViewerDestroy(&patch->viewerMatrix));
    2962:   PetscFunctionReturn(PETSC_SUCCESS);
    2963: }

    2965: static PetscErrorCode PCDestroy_PATCH_Linear(PC pc)
    2966: {
    2967:   PC_PATCH *patch = (PC_PATCH *)pc->data;
    2968:   PetscInt  i;

    2970:   PetscFunctionBegin;
    2971:   if (patch->solver) {
    2972:     for (i = 0; i < patch->npatch; ++i) PetscCall(KSPDestroy((KSP *)&patch->solver[i]));
    2973:     PetscCall(PetscFree(patch->solver));
    2974:   }
    2975:   PetscFunctionReturn(PETSC_SUCCESS);
    2976: }

    2978: static PetscErrorCode PCDestroy_PATCH(PC pc)
    2979: {
    2980:   PC_PATCH *patch = (PC_PATCH *)pc->data;

    2982:   PetscFunctionBegin;
    2983:   PetscCall(PCReset_PATCH(pc));
    2984:   PetscCall((*patch->destroysolver)(pc));
    2985:   PetscCall(PetscFree(pc->data));
    2986:   PetscFunctionReturn(PETSC_SUCCESS);
    2987: }

    2989: static PetscErrorCode PCSetFromOptions_PATCH(PC pc, PetscOptionItems PetscOptionsObject)
    2990: {
    2991:   PC_PATCH            *patch                 = (PC_PATCH *)pc->data;
    2992:   PCPatchConstructType patchConstructionType = PC_PATCH_STAR;
    2993:   char                 sub_mat_type[PETSC_MAX_PATH_LEN];
    2994:   char                 option[PETSC_MAX_PATH_LEN];
    2995:   const char          *prefix;
    2996:   PetscBool            flg, dimflg, codimflg;
    2997:   MPI_Comm             comm;
    2998:   PetscInt            *ifields, nfields, k;
    2999:   PCCompositeType      loctype = PC_COMPOSITE_ADDITIVE;

    3001:   PetscFunctionBegin;
    3002:   PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
    3003:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
    3004:   PetscOptionsHeadBegin(PetscOptionsObject, "Patch solver options");

    3006:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_save_operators", patch->classname));
    3007:   PetscCall(PetscOptionsBool(option, "Store all patch operators for lifetime of object?", "PCPatchSetSaveOperators", patch->save_operators, &patch->save_operators, &flg));

    3009:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_precompute_element_tensors", patch->classname));
    3010:   PetscCall(PetscOptionsBool(option, "Compute each element tensor only once?", "PCPatchSetPrecomputeElementTensors", patch->precomputeElementTensors, &patch->precomputeElementTensors, &flg));
    3011:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_partition_of_unity", patch->classname));
    3012:   PetscCall(PetscOptionsBool(option, "Weight contributions by dof multiplicity?", "PCPatchSetPartitionOfUnity", patch->partition_of_unity, &patch->partition_of_unity, &flg));

    3014:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_local_type", patch->classname));
    3015:   PetscCall(PetscOptionsEnum(option, "Type of local solver composition (additive or multiplicative)", "PCPatchSetLocalComposition", PCCompositeTypes, (PetscEnum)loctype, (PetscEnum *)&loctype, &flg));
    3016:   if (flg) PetscCall(PCPatchSetLocalComposition(pc, loctype));
    3017:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_dense_inverse", patch->classname));
    3018:   PetscCall(PetscOptionsBool(option, "Compute inverses of patch matrices and apply directly? Ignores KSP/PC settings on patch.", "PCPatchSetDenseInverse", patch->denseinverse, &patch->denseinverse, &flg));
    3019:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_dim", patch->classname));
    3020:   PetscCall(PetscOptionsInt(option, "What dimension of mesh point to construct patches by? (0 = vertices)", "PCPATCH", patch->dim, &patch->dim, &dimflg));
    3021:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_codim", patch->classname));
    3022:   PetscCall(PetscOptionsInt(option, "What co-dimension of mesh point to construct patches by? (0 = cells)", "PCPATCH", patch->codim, &patch->codim, &codimflg));
    3023:   PetscCheck(!dimflg || !codimflg, comm, PETSC_ERR_ARG_WRONG, "Can only set one of dimension or co-dimension");

    3025:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_type", patch->classname));
    3026:   PetscCall(PetscOptionsEnum(option, "How should the patches be constructed?", "PCPatchSetConstructType", PCPatchConstructTypes, (PetscEnum)patchConstructionType, (PetscEnum *)&patchConstructionType, &flg));
    3027:   if (flg) PetscCall(PCPatchSetConstructType(pc, patchConstructionType, NULL, NULL));

    3029:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_vanka_dim", patch->classname));
    3030:   PetscCall(PetscOptionsInt(option, "Topological dimension of entities for Vanka to ignore", "PCPATCH", patch->vankadim, &patch->vankadim, &flg));

    3032:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_ignore_dim", patch->classname));
    3033:   PetscCall(PetscOptionsInt(option, "Topological dimension of entities for completion to ignore", "PCPATCH", patch->ignoredim, &patch->ignoredim, &flg));

    3035:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_pardecomp_overlap", patch->classname));
    3036:   PetscCall(PetscOptionsInt(option, "What overlap should we use in construct type pardecomp?", "PCPATCH", patch->pardecomp_overlap, &patch->pardecomp_overlap, &flg));

    3038:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_sub_mat_type", patch->classname));
    3039:   PetscCall(PetscOptionsFList(option, "Matrix type for patch solves", "PCPatchSetSubMatType", MatList, NULL, sub_mat_type, PETSC_MAX_PATH_LEN, &flg));
    3040:   if (flg) PetscCall(PCPatchSetSubMatType(pc, sub_mat_type));

    3042:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_symmetrise_sweep", patch->classname));
    3043:   PetscCall(PetscOptionsBool(option, "Go start->end, end->start?", "PCPATCH", patch->symmetrise_sweep, &patch->symmetrise_sweep, &flg));

    3045:   /* If the user has set the number of subspaces, use that for the buffer size,
    3046:    otherwise use a large number */
    3047:   if (patch->nsubspaces <= 0) {
    3048:     nfields = 128;
    3049:   } else {
    3050:     nfields = patch->nsubspaces;
    3051:   }
    3052:   PetscCall(PetscMalloc1(nfields, &ifields));
    3053:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exclude_subspaces", patch->classname));
    3054:   PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, option, ifields, &nfields, &flg));
    3055:   PetscCheck(!flg || !(patchConstructionType == PC_PATCH_USER), comm, PETSC_ERR_ARG_INCOMP, "We cannot support excluding a subspace with user patches because we do not index patches with a mesh point");
    3056:   if (flg) {
    3057:     PetscCall(PetscHSetIClear(patch->subspaces_to_exclude));
    3058:     for (k = 0; k < nfields; k++) PetscCall(PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]));
    3059:   }
    3060:   PetscCall(PetscFree(ifields));

    3062:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_patches_view", patch->classname));
    3063:   PetscCall(PetscOptionsBool(option, "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg));
    3064:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_cells_view", patch->classname));
    3065:   PetscCall(PetscOptionsCreateViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerCells, &patch->formatCells, &patch->viewCells));
    3066:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_interior_facets_view", patch->classname));
    3067:   PetscCall(PetscOptionsCreateViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerIntFacets, &patch->formatIntFacets, &patch->viewIntFacets));
    3068:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exterior_facets_view", patch->classname));
    3069:   PetscCall(PetscOptionsCreateViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerExtFacets, &patch->formatExtFacets, &patch->viewExtFacets));
    3070:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_points_view", patch->classname));
    3071:   PetscCall(PetscOptionsCreateViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints));
    3072:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_section_view", patch->classname));
    3073:   PetscCall(PetscOptionsCreateViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerSection, &patch->formatSection, &patch->viewSection));
    3074:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_mat_view", patch->classname));
    3075:   PetscCall(PetscOptionsCreateViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix));
    3076:   PetscOptionsHeadEnd();
    3077:   patch->optionsSet = PETSC_TRUE;
    3078:   PetscFunctionReturn(PETSC_SUCCESS);
    3079: }

    3081: static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc)
    3082: {
    3083:   PC_PATCH          *patch = (PC_PATCH *)pc->data;
    3084:   KSPConvergedReason reason;
    3085:   PetscInt           i;

    3087:   PetscFunctionBegin;
    3088:   if (!patch->save_operators) {
    3089:     /* Can't do this here because the sub KSPs don't have an operator attached yet. */
    3090:     PetscFunctionReturn(PETSC_SUCCESS);
    3091:   }
    3092:   if (patch->denseinverse) {
    3093:     /* No solvers */
    3094:     PetscFunctionReturn(PETSC_SUCCESS);
    3095:   }
    3096:   for (i = 0; i < patch->npatch; ++i) {
    3097:     if (!((KSP)patch->solver[i])->setfromoptionscalled) PetscCall(KSPSetFromOptions((KSP)patch->solver[i]));
    3098:     PetscCall(KSPSetUp((KSP)patch->solver[i]));
    3099:     PetscCall(KSPGetConvergedReason((KSP)patch->solver[i], &reason));
    3100:     if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
    3101:   }
    3102:   PetscFunctionReturn(PETSC_SUCCESS);
    3103: }

    3105: static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer)
    3106: {
    3107:   PC_PATCH   *patch = (PC_PATCH *)pc->data;
    3108:   PetscViewer sviewer;
    3109:   PetscBool   isascii;
    3110:   PetscMPIInt rank;

    3112:   PetscFunctionBegin;
    3113:   /* TODO Redo tabbing with set tbas in new style */
    3114:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    3115:   if (!isascii) PetscFunctionReturn(PETSC_SUCCESS);
    3116:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
    3117:   PetscCall(PetscViewerASCIIPushTab(viewer));
    3118:   PetscCall(PetscViewerASCIIPrintf(viewer, "Subspace Correction preconditioner with %" PetscInt_FMT " patches\n", patch->npatch));
    3119:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
    3120:     PetscCall(PetscViewerASCIIPrintf(viewer, "Schwarz type: multiplicative\n"));
    3121:   } else {
    3122:     PetscCall(PetscViewerASCIIPrintf(viewer, "Schwarz type: additive\n"));
    3123:   }
    3124:   if (patch->partition_of_unity) PetscCall(PetscViewerASCIIPrintf(viewer, "Weighting by partition of unity\n"));
    3125:   else PetscCall(PetscViewerASCIIPrintf(viewer, "Not weighting by partition of unity\n"));
    3126:   if (patch->symmetrise_sweep) PetscCall(PetscViewerASCIIPrintf(viewer, "Symmetrising sweep (start->end, then end->start)\n"));
    3127:   else PetscCall(PetscViewerASCIIPrintf(viewer, "Not symmetrising sweep\n"));
    3128:   if (!patch->precomputeElementTensors) PetscCall(PetscViewerASCIIPrintf(viewer, "Not precomputing element tensors (overlapping cells rebuilt in every patch assembly)\n"));
    3129:   else PetscCall(PetscViewerASCIIPrintf(viewer, "Precomputing element tensors (each cell assembled only once)\n"));
    3130:   if (!patch->save_operators) PetscCall(PetscViewerASCIIPrintf(viewer, "Not saving patch operators (rebuilt every PCApply)\n"));
    3131:   else PetscCall(PetscViewerASCIIPrintf(viewer, "Saving patch operators (rebuilt every PCSetUp)\n"));
    3132:   if (patch->patchconstructop == PCPatchConstruct_Star) PetscCall(PetscViewerASCIIPrintf(viewer, "Patch construction operator: star\n"));
    3133:   else if (patch->patchconstructop == PCPatchConstruct_Vanka) PetscCall(PetscViewerASCIIPrintf(viewer, "Patch construction operator: Vanka\n"));
    3134:   else if (patch->patchconstructop == PCPatchConstruct_User) PetscCall(PetscViewerASCIIPrintf(viewer, "Patch construction operator: user-specified\n"));
    3135:   else PetscCall(PetscViewerASCIIPrintf(viewer, "Patch construction operator: unknown\n"));

    3137:   if (patch->denseinverse) {
    3138:     PetscCall(PetscViewerASCIIPrintf(viewer, "Explicitly forming dense inverse and applying patch solver via MatMult.\n"));
    3139:   } else {
    3140:     if (patch->isNonlinear) {
    3141:       PetscCall(PetscViewerASCIIPrintf(viewer, "SNES on patches (all same):\n"));
    3142:     } else {
    3143:       PetscCall(PetscViewerASCIIPrintf(viewer, "KSP on patches (all same):\n"));
    3144:     }
    3145:     if (patch->solver) {
    3146:       PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
    3147:       if (rank == 0) {
    3148:         PetscCall(PetscViewerASCIIPushTab(sviewer));
    3149:         PetscCall(PetscObjectView(patch->solver[0], sviewer));
    3150:         PetscCall(PetscViewerASCIIPopTab(sviewer));
    3151:       }
    3152:       PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
    3153:     } else {
    3154:       PetscCall(PetscViewerASCIIPushTab(viewer));
    3155:       PetscCall(PetscViewerASCIIPrintf(viewer, "Solver not yet set.\n"));
    3156:       PetscCall(PetscViewerASCIIPopTab(viewer));
    3157:     }
    3158:   }
    3159:   PetscCall(PetscViewerASCIIPopTab(viewer));
    3160:   PetscFunctionReturn(PETSC_SUCCESS);
    3161: }

    3163: /*MC
    3164:    PCPATCH - A `PC` object that encapsulates flexible definition of blocks for overlapping and non-overlapping
    3165:    small block additive preconditioners. Block definition is based on topology from
    3166:    a `DM` and equation numbering from a `PetscSection`.

    3168:    Options Database Keys:
    3169: + -pc_patch_cells_view   - Views the process local cell numbers for each patch
    3170: . -pc_patch_points_view  - Views the process local mesh point numbers for each patch
    3171: . -pc_patch_g2l_view     - Views the map between global dofs and patch local dofs for each patch
    3172: . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary
    3173: - -pc_patch_sub_mat_view - Views the matrix associated with each patch

    3175:    Level: intermediate

    3177: .seealso: [](ch_ksp), `PCType`, `PCCreate()`, `PCSetType()`, `PCASM`, `PCJACOBI`, `PCPBJACOBI`, `PCVPBJACOBI`, `SNESPATCH`
    3178: M*/
    3179: PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc)
    3180: {
    3181:   PC_PATCH *patch;

    3183:   PetscFunctionBegin;
    3184:   PetscCall(PetscCitationsRegister(PCPatchCitation, &PCPatchcite));
    3185:   PetscCall(PetscNew(&patch));

    3187:   if (patch->subspaces_to_exclude) PetscCall(PetscHSetIDestroy(&patch->subspaces_to_exclude));
    3188:   PetscCall(PetscHSetICreate(&patch->subspaces_to_exclude));

    3190:   patch->classname   = "pc";
    3191:   patch->isNonlinear = PETSC_FALSE;

    3193:   /* Set some defaults */
    3194:   patch->combined                 = PETSC_FALSE;
    3195:   patch->save_operators           = PETSC_TRUE;
    3196:   patch->local_composition_type   = PC_COMPOSITE_ADDITIVE;
    3197:   patch->precomputeElementTensors = PETSC_FALSE;
    3198:   patch->partition_of_unity       = PETSC_FALSE;
    3199:   patch->codim                    = -1;
    3200:   patch->dim                      = -1;
    3201:   patch->vankadim                 = -1;
    3202:   patch->ignoredim                = -1;
    3203:   patch->pardecomp_overlap        = 0;
    3204:   patch->patchconstructop         = PCPatchConstruct_Star;
    3205:   patch->symmetrise_sweep         = PETSC_FALSE;
    3206:   patch->npatch                   = 0;
    3207:   patch->userIS                   = NULL;
    3208:   patch->optionsSet               = PETSC_FALSE;
    3209:   patch->iterationSet             = NULL;
    3210:   patch->user_patches             = PETSC_FALSE;
    3211:   PetscCall(PetscStrallocpy(MATDENSE, (char **)&patch->sub_mat_type));
    3212:   patch->viewPatches                       = PETSC_FALSE;
    3213:   patch->viewCells                         = PETSC_FALSE;
    3214:   patch->viewPoints                        = PETSC_FALSE;
    3215:   patch->viewSection                       = PETSC_FALSE;
    3216:   patch->viewMatrix                        = PETSC_FALSE;
    3217:   patch->densesolve                        = NULL;
    3218:   patch->setupsolver                       = PCSetUp_PATCH_Linear;
    3219:   patch->applysolver                       = PCApply_PATCH_Linear;
    3220:   patch->resetsolver                       = PCReset_PATCH_Linear;
    3221:   patch->destroysolver                     = PCDestroy_PATCH_Linear;
    3222:   patch->updatemultiplicative              = PCUpdateMultiplicative_PATCH_Linear;
    3223:   patch->dofMappingWithoutToWithArtificial = NULL;
    3224:   patch->dofMappingWithoutToWithAll        = NULL;

    3226:   pc->data                 = (void *)patch;
    3227:   pc->ops->apply           = PCApply_PATCH;
    3228:   pc->ops->applytranspose  = NULL; /* PCApplyTranspose_PATCH; */
    3229:   pc->ops->setup           = PCSetUp_PATCH;
    3230:   pc->ops->reset           = PCReset_PATCH;
    3231:   pc->ops->destroy         = PCDestroy_PATCH;
    3232:   pc->ops->setfromoptions  = PCSetFromOptions_PATCH;
    3233:   pc->ops->setuponblocks   = PCSetUpOnBlocks_PATCH;
    3234:   pc->ops->view            = PCView_PATCH;
    3235:   pc->ops->applyrichardson = NULL;
    3236:   PetscFunctionReturn(PETSC_SUCCESS);
    3237: }