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;
191:   PetscInt  i;

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

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 (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 (i = 0; i < n; ++i) {
229:       const PetscMPIInt *ranks = NULL;
230:       PetscInt           nranks, j;

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 (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 (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 (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:     PetscCallMPI(MPI_Type_contiguous(n, MPIU_INT, &contig));
266:     PetscCallMPI(MPI_Type_commit(&contig));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

389: /* TODO: Docs */
390: PetscErrorCode PCPatchSetSubMatType(PC pc, MatType sub_mat_type)
391: {
392:   PC_PATCH *patch = (PC_PATCH *)pc->data;

394:   PetscFunctionBegin;
395:   if (patch->sub_mat_type) PetscCall(PetscFree(patch->sub_mat_type));
396:   PetscCall(PetscStrallocpy(sub_mat_type, (char **)&patch->sub_mat_type));
397:   PetscFunctionReturn(PETSC_SUCCESS);
398: }

400: /* TODO: Docs */
401: PetscErrorCode PCPatchGetSubMatType(PC pc, MatType *sub_mat_type)
402: {
403:   PC_PATCH *patch = (PC_PATCH *)pc->data;

405:   PetscFunctionBegin;
406:   *sub_mat_type = patch->sub_mat_type;
407:   PetscFunctionReturn(PETSC_SUCCESS);
408: }

410: /* TODO: Docs */
411: PetscErrorCode PCPatchSetCellNumbering(PC pc, PetscSection cellNumbering)
412: {
413:   PC_PATCH *patch = (PC_PATCH *)pc->data;

415:   PetscFunctionBegin;
416:   patch->cellNumbering = cellNumbering;
417:   PetscCall(PetscObjectReference((PetscObject)cellNumbering));
418:   PetscFunctionReturn(PETSC_SUCCESS);
419: }

421: /* TODO: Docs */
422: PetscErrorCode PCPatchGetCellNumbering(PC pc, PetscSection *cellNumbering)
423: {
424:   PC_PATCH *patch = (PC_PATCH *)pc->data;

426:   PetscFunctionBegin;
427:   *cellNumbering = patch->cellNumbering;
428:   PetscFunctionReturn(PETSC_SUCCESS);
429: }

431: /* TODO: Docs */
432: PetscErrorCode PCPatchSetConstructType(PC pc, PCPatchConstructType ctype, PetscErrorCode (*func)(PC, PetscInt *, IS **, IS *, void *), void *ctx)
433: {
434:   PC_PATCH *patch = (PC_PATCH *)pc->data;

436:   PetscFunctionBegin;
437:   patch->ctype = ctype;
438:   switch (ctype) {
439:   case PC_PATCH_STAR:
440:     patch->user_patches     = PETSC_FALSE;
441:     patch->patchconstructop = PCPatchConstruct_Star;
442:     break;
443:   case PC_PATCH_VANKA:
444:     patch->user_patches     = PETSC_FALSE;
445:     patch->patchconstructop = PCPatchConstruct_Vanka;
446:     break;
447:   case PC_PATCH_PARDECOMP:
448:     patch->user_patches     = PETSC_FALSE;
449:     patch->patchconstructop = PCPatchConstruct_Pardecomp;
450:     break;
451:   case PC_PATCH_USER:
452:   case PC_PATCH_PYTHON:
453:     patch->user_patches     = PETSC_TRUE;
454:     patch->patchconstructop = PCPatchConstruct_User;
455:     if (func) {
456:       patch->userpatchconstructionop = func;
457:       patch->userpatchconstructctx   = ctx;
458:     }
459:     break;
460:   default:
461:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "Unknown patch construction type %" PetscInt_FMT, (PetscInt)patch->ctype);
462:   }
463:   PetscFunctionReturn(PETSC_SUCCESS);
464: }

466: /* TODO: Docs */
467: PetscErrorCode PCPatchGetConstructType(PC pc, PCPatchConstructType *ctype, PetscErrorCode (**func)(PC, PetscInt *, IS **, IS *, void *), void **ctx)
468: {
469:   PC_PATCH *patch = (PC_PATCH *)pc->data;

471:   PetscFunctionBegin;
472:   *ctype = patch->ctype;
473:   switch (patch->ctype) {
474:   case PC_PATCH_STAR:
475:   case PC_PATCH_VANKA:
476:   case PC_PATCH_PARDECOMP:
477:     break;
478:   case PC_PATCH_USER:
479:   case PC_PATCH_PYTHON:
480:     *func = patch->userpatchconstructionop;
481:     *ctx  = patch->userpatchconstructctx;
482:     break;
483:   default:
484:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "Unknown patch construction type %" PetscInt_FMT, (PetscInt)patch->ctype);
485:   }
486:   PetscFunctionReturn(PETSC_SUCCESS);
487: }

489: /* TODO: Docs */
490: 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)
491: {
492:   PC_PATCH *patch = (PC_PATCH *)pc->data;
493:   DM        dm, plex;
494:   PetscSF  *sfs;
495:   PetscInt  cStart, cEnd, i, j;

497:   PetscFunctionBegin;
498:   PetscCall(PCGetDM(pc, &dm));
499:   PetscCall(DMConvert(dm, DMPLEX, &plex));
500:   dm = plex;
501:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
502:   PetscCall(PetscMalloc1(nsubspaces, &sfs));
503:   PetscCall(PetscMalloc1(nsubspaces, &patch->dofSection));
504:   PetscCall(PetscMalloc1(nsubspaces, &patch->bs));
505:   PetscCall(PetscMalloc1(nsubspaces, &patch->nodesPerCell));
506:   PetscCall(PetscMalloc1(nsubspaces, &patch->cellNodeMap));
507:   PetscCall(PetscMalloc1(nsubspaces + 1, &patch->subspaceOffsets));

509:   patch->nsubspaces       = nsubspaces;
510:   patch->totalDofsPerCell = 0;
511:   for (i = 0; i < nsubspaces; ++i) {
512:     PetscCall(DMGetLocalSection(dms[i], &patch->dofSection[i]));
513:     PetscCall(PetscObjectReference((PetscObject)patch->dofSection[i]));
514:     PetscCall(DMGetSectionSF(dms[i], &sfs[i]));
515:     patch->bs[i]           = bs[i];
516:     patch->nodesPerCell[i] = nodesPerCell[i];
517:     patch->totalDofsPerCell += nodesPerCell[i] * bs[i];
518:     PetscCall(PetscMalloc1((cEnd - cStart) * nodesPerCell[i], &patch->cellNodeMap[i]));
519:     for (j = 0; j < (cEnd - cStart) * nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j];
520:     patch->subspaceOffsets[i] = subspaceOffsets[i];
521:   }
522:   PetscCall(PCPatchCreateDefaultSF_Private(pc, nsubspaces, sfs, patch->bs));
523:   PetscCall(PetscFree(sfs));

525:   patch->subspaceOffsets[nsubspaces] = subspaceOffsets[nsubspaces];
526:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes));
527:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes));
528:   PetscCall(DMDestroy(&dm));
529:   PetscFunctionReturn(PETSC_SUCCESS);
530: }

532: /* TODO: Docs */
533: static PetscErrorCode PCPatchSetDiscretisationInfoCombined(PC pc, DM dm, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes)
534: {
535:   PC_PATCH *patch = (PC_PATCH *)pc->data;
536:   PetscInt  cStart, cEnd, i, j;

538:   PetscFunctionBegin;
539:   patch->combined = PETSC_TRUE;
540:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
541:   PetscCall(DMGetNumFields(dm, &patch->nsubspaces));
542:   PetscCall(PetscCalloc1(patch->nsubspaces, &patch->dofSection));
543:   PetscCall(PetscMalloc1(patch->nsubspaces, &patch->bs));
544:   PetscCall(PetscMalloc1(patch->nsubspaces, &patch->nodesPerCell));
545:   PetscCall(PetscMalloc1(patch->nsubspaces, &patch->cellNodeMap));
546:   PetscCall(PetscCalloc1(patch->nsubspaces + 1, &patch->subspaceOffsets));
547:   PetscCall(DMGetLocalSection(dm, &patch->dofSection[0]));
548:   PetscCall(PetscObjectReference((PetscObject)patch->dofSection[0]));
549:   PetscCall(PetscSectionGetStorageSize(patch->dofSection[0], &patch->subspaceOffsets[patch->nsubspaces]));
550:   patch->totalDofsPerCell = 0;
551:   for (i = 0; i < patch->nsubspaces; ++i) {
552:     patch->bs[i]           = 1;
553:     patch->nodesPerCell[i] = nodesPerCell[i];
554:     patch->totalDofsPerCell += nodesPerCell[i];
555:     PetscCall(PetscMalloc1((cEnd - cStart) * nodesPerCell[i], &patch->cellNodeMap[i]));
556:     for (j = 0; j < (cEnd - cStart) * nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j];
557:   }
558:   PetscCall(DMGetSectionSF(dm, &patch->sectionSF));
559:   PetscCall(PetscObjectReference((PetscObject)patch->sectionSF));
560:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes));
561:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes));
562:   PetscFunctionReturn(PETSC_SUCCESS);
563: }

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

568:   Logically Collective

570:   Input Parameters:
571: + pc   - The `PC`
572: . func - The callback function
573: - ctx  - The user context

575:   Calling sequence of `func`:
576: + pc               - The `PC`
577: . point            - The point
578: . x                - The input solution (not used in linear problems)
579: . f                - The patch residual vector
580: . cellIS           - An array of the cell numbers
581: . n                - The size of `dofsArray`
582: . dofsArray        - The dofmap for the dofs to be solved for
583: . dofsArrayWithAll - The dofmap for all dofs on the patch
584: - ctx              - The user context

586:   Level: advanced

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

591: .seealso: [](ch_ksp), `PCPatchSetComputeOperator()`, `PCPatchGetComputeOperator()`, `PCPatchSetDiscretisationInfo()`, `PCPatchSetComputeFunctionInteriorFacets()`
592: @*/
593: 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)
594: {
595:   PC_PATCH *patch = (PC_PATCH *)pc->data;

597:   PetscFunctionBegin;
598:   patch->usercomputef    = func;
599:   patch->usercomputefctx = ctx;
600:   PetscFunctionReturn(PETSC_SUCCESS);
601: }

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

606:   Logically Collective

608:   Input Parameters:
609: + pc   - The `PC`
610: . func - The callback function
611: - ctx  - The user context

613:   Calling sequence of `func`:
614: + pc               - The `PC`
615: . point            - The point
616: . x                - The input solution (not used in linear problems)
617: . f                - The patch residual vector
618: . facetIS          - An array of the facet numbers
619: . n                - The size of `dofsArray`
620: . dofsArray        - The dofmap for the dofs to be solved for
621: . dofsArrayWithAll - The dofmap for all dofs on the patch
622: - ctx              - The user context

624:   Level: advanced

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

629: .seealso: [](ch_ksp), `PCPatchSetComputeOperator()`, `PCPatchGetComputeOperator()`, `PCPatchSetDiscretisationInfo()`, `PCPatchSetComputeFunction()`
630: @*/
631: 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)
632: {
633:   PC_PATCH *patch = (PC_PATCH *)pc->data;

635:   PetscFunctionBegin;
636:   patch->usercomputefintfacet    = func;
637:   patch->usercomputefintfacetctx = ctx;
638:   PetscFunctionReturn(PETSC_SUCCESS);
639: }

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

644:   Logically Collective

646:   Input Parameters:
647: + pc   - The `PC`
648: . func - The callback function
649: - ctx  - The user context

651:   Calling sequence of `func`:
652: + pc               - The `PC`
653: . point            - The point
654: . x                - The input solution (not used in linear problems)
655: . mat              - The patch matrix
656: . facetIS          - An array of the cell numbers
657: . n                - The size of `dofsArray`
658: . dofsArray        - The dofmap for the dofs to be solved for
659: . dofsArrayWithAll - The dofmap for all dofs on the patch
660: - ctx              - The user context

662:   Level: advanced

664:   Note:
665:   The matrix entries have been set to zero before the call.

667: .seealso: [](ch_ksp), `PCPatchGetComputeOperator()`, `PCPatchSetComputeFunction()`, `PCPatchSetDiscretisationInfo()`
668: @*/
669: 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)
670: {
671:   PC_PATCH *patch = (PC_PATCH *)pc->data;

673:   PetscFunctionBegin;
674:   patch->usercomputeop    = func;
675:   patch->usercomputeopctx = ctx;
676:   PetscFunctionReturn(PETSC_SUCCESS);
677: }

679: /*@C

681:   PCPatchSetComputeOperatorInteriorFacets - Set the callback function used to compute facet integrals for patch matrices

683:   Logically Collective

685:   Input Parameters:
686: + pc   - The `PC`
687: . func - The callback function
688: - ctx  - The user context

690:   Calling sequence of `func`:
691: + pc               - The `PC`
692: . point            - The point
693: . x                - The input solution (not used in linear problems)
694: . mat              - The patch matrix
695: . facetIS          - An array of the facet numbers
696: . n                - The size of `dofsArray`
697: . dofsArray        - The dofmap for the dofs to be solved for
698: . dofsArrayWithAll - The dofmap for all dofs on the patch
699: - ctx              - The user context

701:   Level: advanced

703:   Note:
704:   The matrix entries have been set to zero before the call.

706: .seealso: [](ch_ksp), `PCPatchGetComputeOperator()`, `PCPatchSetComputeFunction()`, `PCPatchSetDiscretisationInfo()`
707: @*/
708: 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)
709: {
710:   PC_PATCH *patch = (PC_PATCH *)pc->data;

712:   PetscFunctionBegin;
713:   patch->usercomputeopintfacet    = func;
714:   patch->usercomputeopintfacetctx = ctx;
715:   PetscFunctionReturn(PETSC_SUCCESS);
716: }

718: /* On entry, ht contains the topological entities whose dofs we are responsible for solving for;
719:    on exit, cht contains all the topological entities we need to compute their residuals.
720:    In full generality this should incorporate knowledge of the sparsity pattern of the matrix;
721:    here we assume a standard FE sparsity pattern.*/
722: /* TODO: Use DMPlexGetAdjacency() */
723: static PetscErrorCode PCPatchCompleteCellPatch(PC pc, PetscHSetI ht, PetscHSetI cht)
724: {
725:   DM            dm, plex;
726:   PC_PATCH     *patch = (PC_PATCH *)pc->data;
727:   PetscHashIter hi;
728:   PetscInt      point;
729:   PetscInt     *star = NULL, *closure = NULL;
730:   PetscInt      ignoredim, iStart = 0, iEnd = -1, starSize, closureSize, si, ci;
731:   PetscInt     *fStar = NULL, *fClosure = NULL;
732:   PetscInt      fBegin, fEnd, fsi, fci, fStarSize, fClosureSize;

734:   PetscFunctionBegin;
735:   PetscCall(PCGetDM(pc, &dm));
736:   PetscCall(DMConvert(dm, DMPLEX, &plex));
737:   dm = plex;
738:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fBegin, &fEnd));
739:   PetscCall(PCPatchGetIgnoreDim(pc, &ignoredim));
740:   if (ignoredim >= 0) PetscCall(DMPlexGetDepthStratum(dm, ignoredim, &iStart, &iEnd));
741:   PetscCall(PetscHSetIClear(cht));
742:   PetscHashIterBegin(ht, hi);
743:   while (!PetscHashIterAtEnd(ht, hi)) {
744:     PetscHashIterGetKey(ht, hi, point);
745:     PetscHashIterNext(ht, hi);

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

781: static PetscErrorCode PCPatchGetGlobalDofs(PC pc, PetscSection dofSection[], PetscInt f, PetscBool combined, PetscInt p, PetscInt *dof, PetscInt *off)
782: {
783:   PetscFunctionBegin;
784:   if (combined) {
785:     if (f < 0) {
786:       if (dof) PetscCall(PetscSectionGetDof(dofSection[0], p, dof));
787:       if (off) PetscCall(PetscSectionGetOffset(dofSection[0], p, off));
788:     } else {
789:       if (dof) PetscCall(PetscSectionGetFieldDof(dofSection[0], p, f, dof));
790:       if (off) PetscCall(PetscSectionGetFieldOffset(dofSection[0], p, f, off));
791:     }
792:   } else {
793:     if (f < 0) {
794:       PC_PATCH *patch = (PC_PATCH *)pc->data;
795:       PetscInt  fdof, g;

797:       if (dof) {
798:         *dof = 0;
799:         for (g = 0; g < patch->nsubspaces; ++g) {
800:           PetscCall(PetscSectionGetDof(dofSection[g], p, &fdof));
801:           *dof += fdof;
802:         }
803:       }
804:       if (off) {
805:         *off = 0;
806:         for (g = 0; g < patch->nsubspaces; ++g) {
807:           PetscCall(PetscSectionGetOffset(dofSection[g], p, &fdof));
808:           *off += fdof;
809:         }
810:       }
811:     } else {
812:       if (dof) PetscCall(PetscSectionGetDof(dofSection[f], p, dof));
813:       if (off) PetscCall(PetscSectionGetOffset(dofSection[f], p, off));
814:     }
815:   }
816:   PetscFunctionReturn(PETSC_SUCCESS);
817: }

819: /* Given a hash table with a set of topological entities (pts), compute the degrees of
820:    freedom in global concatenated numbering on those entities.
821:    For Vanka smoothing, this needs to do something special: ignore dofs of the
822:    constraint subspace on entities that aren't the base entity we're building the patch
823:    around. */
824: static PetscErrorCode PCPatchGetPointDofs(PC pc, PetscHSetI pts, PetscHSetI dofs, PetscInt base, PetscHSetI *subspaces_to_exclude)
825: {
826:   PC_PATCH     *patch = (PC_PATCH *)pc->data;
827:   PetscHashIter hi;
828:   PetscInt      ldof, loff;
829:   PetscInt      k, p;

831:   PetscFunctionBegin;
832:   PetscCall(PetscHSetIClear(dofs));
833:   for (k = 0; k < patch->nsubspaces; ++k) {
834:     PetscInt subspaceOffset = patch->subspaceOffsets[k];
835:     PetscInt bs             = patch->bs[k];
836:     PetscInt j, l;

838:     if (subspaces_to_exclude != NULL) {
839:       PetscBool should_exclude_k = PETSC_FALSE;
840:       PetscCall(PetscHSetIHas(*subspaces_to_exclude, k, &should_exclude_k));
841:       if (should_exclude_k) {
842:         /* only get this subspace dofs at the base entity, not any others */
843:         PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, base, &ldof, &loff));
844:         if (0 == ldof) continue;
845:         for (j = loff; j < ldof + loff; ++j) {
846:           for (l = 0; l < bs; ++l) {
847:             PetscInt dof = bs * j + l + subspaceOffset;
848:             PetscCall(PetscHSetIAdd(dofs, dof));
849:           }
850:         }
851:         continue; /* skip the other dofs of this subspace */
852:       }
853:     }

855:     PetscHashIterBegin(pts, hi);
856:     while (!PetscHashIterAtEnd(pts, hi)) {
857:       PetscHashIterGetKey(pts, hi, p);
858:       PetscHashIterNext(pts, hi);
859:       PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, p, &ldof, &loff));
860:       if (0 == ldof) continue;
861:       for (j = loff; j < ldof + loff; ++j) {
862:         for (l = 0; l < bs; ++l) {
863:           PetscInt dof = bs * j + l + subspaceOffset;
864:           PetscCall(PetscHSetIAdd(dofs, dof));
865:         }
866:       }
867:     }
868:   }
869:   PetscFunctionReturn(PETSC_SUCCESS);
870: }

872: /* Given two hash tables A and B, compute the keys in B that are not in A, and put them in C */
873: static PetscErrorCode PCPatchComputeSetDifference_Private(PetscHSetI A, PetscHSetI B, PetscHSetI C)
874: {
875:   PetscHashIter hi;
876:   PetscInt      key;
877:   PetscBool     flg;

879:   PetscFunctionBegin;
880:   PetscCall(PetscHSetIClear(C));
881:   PetscHashIterBegin(B, hi);
882:   while (!PetscHashIterAtEnd(B, hi)) {
883:     PetscHashIterGetKey(B, hi, key);
884:     PetscHashIterNext(B, hi);
885:     PetscCall(PetscHSetIHas(A, key, &flg));
886:     if (!flg) PetscCall(PetscHSetIAdd(C, key));
887:   }
888:   PetscFunctionReturn(PETSC_SUCCESS);
889: }

891: // PetscClangLinter pragma disable: -fdoc-sowing-chars
892: /*
893:   PCPatchCreateCellPatches - create patches.

895:   Input Parameter:
896:   . dm - The DMPlex object defining the mesh

898:   Output Parameters:
899:   + cellCounts  - Section with counts of cells around each vertex
900:   . cells       - IS of the cell point indices of cells in each patch
901:   . pointCounts - Section with counts of cells around each vertex
902:   - point       - IS of the cell point indices of cells in each patch
903:  */
904: static PetscErrorCode PCPatchCreateCellPatches(PC pc)
905: {
906:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
907:   DMLabel         ghost = NULL;
908:   DM              dm, plex;
909:   PetscHSetI      ht = NULL, cht = NULL;
910:   PetscSection    cellCounts, pointCounts, intFacetCounts, extFacetCounts;
911:   PetscInt       *cellsArray, *pointsArray, *intFacetsArray, *extFacetsArray, *intFacetsToPatchCell;
912:   PetscInt        numCells, numPoints, numIntFacets, numExtFacets;
913:   const PetscInt *leaves;
914:   PetscInt        nleaves, pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, v;
915:   PetscBool       isFiredrake;

917:   PetscFunctionBegin;
918:   /* Used to keep track of the cells in the patch. */
919:   PetscCall(PetscHSetICreate(&ht));
920:   PetscCall(PetscHSetICreate(&cht));

922:   PetscCall(PCGetDM(pc, &dm));
923:   PetscCheck(dm, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "DM not yet set on patch PC");
924:   PetscCall(DMConvert(dm, DMPLEX, &plex));
925:   dm = plex;
926:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
927:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));

929:   if (patch->user_patches) {
930:     PetscCall(patch->userpatchconstructionop(pc, &patch->npatch, &patch->userIS, &patch->iterationSet, patch->userpatchconstructctx));
931:     vStart = 0;
932:     vEnd   = patch->npatch;
933:   } else if (patch->ctype == PC_PATCH_PARDECOMP) {
934:     vStart = 0;
935:     vEnd   = 1;
936:   } else if (patch->codim < 0) {
937:     if (patch->dim < 0) PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
938:     else PetscCall(DMPlexGetDepthStratum(dm, patch->dim, &vStart, &vEnd));
939:   } else PetscCall(DMPlexGetHeightStratum(dm, patch->codim, &vStart, &vEnd));
940:   patch->npatch = vEnd - vStart;

942:   /* These labels mark the owned points.  We only create patches around points that this process owns. */
943:   PetscCall(DMHasLabel(dm, "pyop2_ghost", &isFiredrake));
944:   if (isFiredrake) {
945:     PetscCall(DMGetLabel(dm, "pyop2_ghost", &ghost));
946:     PetscCall(DMLabelCreateIndex(ghost, pStart, pEnd));
947:   } else {
948:     PetscSF sf;

950:     PetscCall(DMGetPointSF(dm, &sf));
951:     PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL));
952:     nleaves = PetscMax(nleaves, 0);
953:   }

955:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->cellCounts));
956:   PetscCall(PetscObjectSetName((PetscObject)patch->cellCounts, "Patch Cell Layout"));
957:   cellCounts = patch->cellCounts;
958:   PetscCall(PetscSectionSetChart(cellCounts, vStart, vEnd));
959:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->pointCounts));
960:   PetscCall(PetscObjectSetName((PetscObject)patch->pointCounts, "Patch Point Layout"));
961:   pointCounts = patch->pointCounts;
962:   PetscCall(PetscSectionSetChart(pointCounts, vStart, vEnd));
963:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->extFacetCounts));
964:   PetscCall(PetscObjectSetName((PetscObject)patch->extFacetCounts, "Patch Exterior Facet Layout"));
965:   extFacetCounts = patch->extFacetCounts;
966:   PetscCall(PetscSectionSetChart(extFacetCounts, vStart, vEnd));
967:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->intFacetCounts));
968:   PetscCall(PetscObjectSetName((PetscObject)patch->intFacetCounts, "Patch Interior Facet Layout"));
969:   intFacetCounts = patch->intFacetCounts;
970:   PetscCall(PetscSectionSetChart(intFacetCounts, vStart, vEnd));
971:   /* Count cells and points in the patch surrounding each entity */
972:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
973:   for (v = vStart; v < vEnd; ++v) {
974:     PetscHashIter hi;
975:     PetscInt      chtSize, loc = -1;
976:     PetscBool     flg;

978:     if (!patch->user_patches && patch->ctype != PC_PATCH_PARDECOMP) {
979:       if (ghost) PetscCall(DMLabelHasPoint(ghost, v, &flg));
980:       else {
981:         PetscCall(PetscFindInt(v, nleaves, leaves, &loc));
982:         flg = loc >= 0 ? PETSC_TRUE : PETSC_FALSE;
983:       }
984:       /* Not an owned entity, don't make a cell patch. */
985:       if (flg) continue;
986:     }

988:     PetscCall(patch->patchconstructop((void *)patch, dm, v, ht));
989:     PetscCall(PCPatchCompleteCellPatch(pc, ht, cht));
990:     PetscCall(PetscHSetIGetSize(cht, &chtSize));
991:     /* empty patch, continue */
992:     if (chtSize == 0) continue;

994:     /* safe because size(cht) > 0 from above */
995:     PetscHashIterBegin(cht, hi);
996:     while (!PetscHashIterAtEnd(cht, hi)) {
997:       PetscInt point, pdof;

999:       PetscHashIterGetKey(cht, hi, point);
1000:       if (fStart <= point && point < fEnd) {
1001:         const PetscInt *support;
1002:         PetscInt        supportSize, p;
1003:         PetscBool       interior = PETSC_TRUE;
1004:         PetscCall(DMPlexGetSupport(dm, point, &support));
1005:         PetscCall(DMPlexGetSupportSize(dm, point, &supportSize));
1006:         if (supportSize == 1) {
1007:           interior = PETSC_FALSE;
1008:         } else {
1009:           for (p = 0; p < supportSize; p++) {
1010:             PetscBool found;
1011:             /* FIXME: can I do this while iterating over cht? */
1012:             PetscCall(PetscHSetIHas(cht, support[p], &found));
1013:             if (!found) {
1014:               interior = PETSC_FALSE;
1015:               break;
1016:             }
1017:           }
1018:         }
1019:         if (interior) {
1020:           PetscCall(PetscSectionAddDof(intFacetCounts, v, 1));
1021:         } else {
1022:           PetscCall(PetscSectionAddDof(extFacetCounts, v, 1));
1023:         }
1024:       }
1025:       PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL));
1026:       if (pdof) PetscCall(PetscSectionAddDof(pointCounts, v, 1));
1027:       if (point >= cStart && point < cEnd) PetscCall(PetscSectionAddDof(cellCounts, v, 1));
1028:       PetscHashIterNext(cht, hi);
1029:     }
1030:   }
1031:   if (isFiredrake) PetscCall(DMLabelDestroyIndex(ghost));

1033:   PetscCall(PetscSectionSetUp(cellCounts));
1034:   PetscCall(PetscSectionGetStorageSize(cellCounts, &numCells));
1035:   PetscCall(PetscMalloc1(numCells, &cellsArray));
1036:   PetscCall(PetscSectionSetUp(pointCounts));
1037:   PetscCall(PetscSectionGetStorageSize(pointCounts, &numPoints));
1038:   PetscCall(PetscMalloc1(numPoints, &pointsArray));

1040:   PetscCall(PetscSectionSetUp(intFacetCounts));
1041:   PetscCall(PetscSectionSetUp(extFacetCounts));
1042:   PetscCall(PetscSectionGetStorageSize(intFacetCounts, &numIntFacets));
1043:   PetscCall(PetscSectionGetStorageSize(extFacetCounts, &numExtFacets));
1044:   PetscCall(PetscMalloc1(numIntFacets, &intFacetsArray));
1045:   PetscCall(PetscMalloc1(numIntFacets * 2, &intFacetsToPatchCell));
1046:   PetscCall(PetscMalloc1(numExtFacets, &extFacetsArray));

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

1053:     PetscCall(PetscSectionGetDof(pointCounts, v, &dof));
1054:     PetscCall(PetscSectionGetOffset(pointCounts, v, &off));
1055:     PetscCall(PetscSectionGetDof(cellCounts, v, &cdof));
1056:     PetscCall(PetscSectionGetOffset(cellCounts, v, &coff));
1057:     PetscCall(PetscSectionGetDof(intFacetCounts, v, &ifdof));
1058:     PetscCall(PetscSectionGetOffset(intFacetCounts, v, &ifoff));
1059:     PetscCall(PetscSectionGetDof(extFacetCounts, v, &efdof));
1060:     PetscCall(PetscSectionGetOffset(extFacetCounts, v, &efoff));
1061:     if (dof <= 0) continue;
1062:     PetscCall(patch->patchconstructop((void *)patch, dm, v, ht));
1063:     PetscCall(PCPatchCompleteCellPatch(pc, ht, cht));
1064:     PetscHashIterBegin(cht, hi);
1065:     while (!PetscHashIterAtEnd(cht, hi)) {
1066:       PetscInt point;

1068:       PetscHashIterGetKey(cht, hi, point);
1069:       if (fStart <= point && point < fEnd) {
1070:         const PetscInt *support;
1071:         PetscInt        supportSize, p;
1072:         PetscBool       interior = PETSC_TRUE;
1073:         PetscCall(DMPlexGetSupport(dm, point, &support));
1074:         PetscCall(DMPlexGetSupportSize(dm, point, &supportSize));
1075:         if (supportSize == 1) {
1076:           interior = PETSC_FALSE;
1077:         } else {
1078:           for (p = 0; p < supportSize; p++) {
1079:             PetscBool found;
1080:             /* FIXME: can I do this while iterating over cht? */
1081:             PetscCall(PetscHSetIHas(cht, support[p], &found));
1082:             if (!found) {
1083:               interior = PETSC_FALSE;
1084:               break;
1085:             }
1086:           }
1087:         }
1088:         if (interior) {
1089:           intFacetsToPatchCell[2 * (ifoff + ifn)]     = support[0];
1090:           intFacetsToPatchCell[2 * (ifoff + ifn) + 1] = support[1];
1091:           intFacetsArray[ifoff + ifn++]               = point;
1092:         } else {
1093:           extFacetsArray[efoff + efn++] = point;
1094:         }
1095:       }
1096:       PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL));
1097:       if (pdof) pointsArray[off + n++] = point;
1098:       if (point >= cStart && point < cEnd) cellsArray[coff + cn++] = point;
1099:       PetscHashIterNext(cht, hi);
1100:     }
1101:     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);
1102:     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);
1103:     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);
1104:     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);

1106:     for (ifn = 0; ifn < ifdof; ifn++) {
1107:       PetscInt  cell0  = intFacetsToPatchCell[2 * (ifoff + ifn)];
1108:       PetscInt  cell1  = intFacetsToPatchCell[2 * (ifoff + ifn) + 1];
1109:       PetscBool found0 = PETSC_FALSE, found1 = PETSC_FALSE;
1110:       for (n = 0; n < cdof; n++) {
1111:         if (!found0 && cell0 == cellsArray[coff + n]) {
1112:           intFacetsToPatchCell[2 * (ifoff + ifn)] = n;
1113:           found0                                  = PETSC_TRUE;
1114:         }
1115:         if (!found1 && cell1 == cellsArray[coff + n]) {
1116:           intFacetsToPatchCell[2 * (ifoff + ifn) + 1] = n;
1117:           found1                                      = PETSC_TRUE;
1118:         }
1119:         if (found0 && found1) break;
1120:       }
1121:       PetscCheck(found0 && found1, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Didn't manage to find local point numbers for facet support");
1122:     }
1123:   }
1124:   PetscCall(PetscHSetIDestroy(&ht));
1125:   PetscCall(PetscHSetIDestroy(&cht));

1127:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numCells, cellsArray, PETSC_OWN_POINTER, &patch->cells));
1128:   PetscCall(PetscObjectSetName((PetscObject)patch->cells, "Patch Cells"));
1129:   if (patch->viewCells) {
1130:     PetscCall(ObjectView((PetscObject)patch->cellCounts, patch->viewerCells, patch->formatCells));
1131:     PetscCall(ObjectView((PetscObject)patch->cells, patch->viewerCells, patch->formatCells));
1132:   }
1133:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray, PETSC_OWN_POINTER, &patch->intFacets));
1134:   PetscCall(PetscObjectSetName((PetscObject)patch->intFacets, "Patch Interior Facets"));
1135:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 2 * numIntFacets, intFacetsToPatchCell, PETSC_OWN_POINTER, &patch->intFacetsToPatchCell));
1136:   PetscCall(PetscObjectSetName((PetscObject)patch->intFacetsToPatchCell, "Patch Interior Facets local support"));
1137:   if (patch->viewIntFacets) {
1138:     PetscCall(ObjectView((PetscObject)patch->intFacetCounts, patch->viewerIntFacets, patch->formatIntFacets));
1139:     PetscCall(ObjectView((PetscObject)patch->intFacets, patch->viewerIntFacets, patch->formatIntFacets));
1140:     PetscCall(ObjectView((PetscObject)patch->intFacetsToPatchCell, patch->viewerIntFacets, patch->formatIntFacets));
1141:   }
1142:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numExtFacets, extFacetsArray, PETSC_OWN_POINTER, &patch->extFacets));
1143:   PetscCall(PetscObjectSetName((PetscObject)patch->extFacets, "Patch Exterior Facets"));
1144:   if (patch->viewExtFacets) {
1145:     PetscCall(ObjectView((PetscObject)patch->extFacetCounts, patch->viewerExtFacets, patch->formatExtFacets));
1146:     PetscCall(ObjectView((PetscObject)patch->extFacets, patch->viewerExtFacets, patch->formatExtFacets));
1147:   }
1148:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPoints, pointsArray, PETSC_OWN_POINTER, &patch->points));
1149:   PetscCall(PetscObjectSetName((PetscObject)patch->points, "Patch Points"));
1150:   if (patch->viewPoints) {
1151:     PetscCall(ObjectView((PetscObject)patch->pointCounts, patch->viewerPoints, patch->formatPoints));
1152:     PetscCall(ObjectView((PetscObject)patch->points, patch->viewerPoints, patch->formatPoints));
1153:   }
1154:   PetscCall(DMDestroy(&dm));
1155:   PetscFunctionReturn(PETSC_SUCCESS);
1156: }

1158: /*
1159:   PCPatchCreateCellPatchDiscretisationInfo - Build the dof maps for cell patches

1161:   Input Parameters:
1162:   + dm - The DMPlex object defining the mesh
1163:   . cellCounts - Section with counts of cells around each vertex
1164:   . cells - IS of the cell point indices of cells in each patch
1165:   . cellNumbering - Section mapping plex cell points to Firedrake cell indices.
1166:   . nodesPerCell - number of nodes per cell.
1167:   - cellNodeMap - map from cells to node indices (nodesPerCell * numCells)

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

1218:   PetscFunctionBegin;
1219:   PetscCall(PCGetDM(pc, &dm));
1220:   PetscCall(DMConvert(dm, DMPLEX, &plex));
1221:   dm = plex;
1222:   /* dofcounts section is cellcounts section * dofPerCell */
1223:   PetscCall(PetscSectionGetStorageSize(cellCounts, &numCells));
1224:   PetscCall(PetscSectionGetStorageSize(patch->pointCounts, &numPoints));
1225:   numDofs = numCells * totalDofsPerCell;
1226:   PetscCall(PetscMalloc1(numDofs, &dofsArray));
1227:   PetscCall(PetscMalloc1(numPoints * Nf, &offsArray));
1228:   PetscCall(PetscMalloc1(numDofs, &asmArray));
1229:   PetscCall(PetscMalloc1(numCells, &newCellsArray));
1230:   PetscCall(PetscSectionGetChart(cellCounts, &vStart, &vEnd));
1231:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCounts));
1232:   gtolCounts = patch->gtolCounts;
1233:   PetscCall(PetscSectionSetChart(gtolCounts, vStart, vEnd));
1234:   PetscCall(PetscObjectSetName((PetscObject)patch->gtolCounts, "Patch Global Index Section"));

1236:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1237:     PetscCall(PetscMalloc1(numPoints * Nf, &offsArrayWithArtificial));
1238:     PetscCall(PetscMalloc1(numDofs, &asmArrayWithArtificial));
1239:     PetscCall(PetscMalloc1(numDofs, &dofsArrayWithArtificial));
1240:     PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithArtificial));
1241:     gtolCountsWithArtificial = patch->gtolCountsWithArtificial;
1242:     PetscCall(PetscSectionSetChart(gtolCountsWithArtificial, vStart, vEnd));
1243:     PetscCall(PetscObjectSetName((PetscObject)patch->gtolCountsWithArtificial, "Patch Global Index Section Including Artificial BCs"));
1244:   }

1246:   isNonlinear = patch->isNonlinear;
1247:   if (isNonlinear) {
1248:     PetscCall(PetscMalloc1(numPoints * Nf, &offsArrayWithAll));
1249:     PetscCall(PetscMalloc1(numDofs, &asmArrayWithAll));
1250:     PetscCall(PetscMalloc1(numDofs, &dofsArrayWithAll));
1251:     PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithAll));
1252:     gtolCountsWithAll = patch->gtolCountsWithAll;
1253:     PetscCall(PetscSectionSetChart(gtolCountsWithAll, vStart, vEnd));
1254:     PetscCall(PetscObjectSetName((PetscObject)patch->gtolCountsWithAll, "Patch Global Index Section Including All BCs"));
1255:   }

1257:   /* Outside the patch loop, get the dofs that are globally-enforced Dirichlet
1258:    conditions */
1259:   PetscCall(PetscHSetICreate(&globalBcs));
1260:   PetscCall(ISGetIndices(patch->ghostBcNodes, &bcNodes));
1261:   PetscCall(ISGetSize(patch->ghostBcNodes, &numBcs));
1262:   for (i = 0; i < numBcs; ++i) { PetscCall(PetscHSetIAdd(globalBcs, bcNodes[i])); /* these are already in concatenated numbering */ }
1263:   PetscCall(ISRestoreIndices(patch->ghostBcNodes, &bcNodes));
1264:   PetscCall(ISDestroy(&patch->ghostBcNodes)); /* memory optimisation */

1266:   /* Hash tables for artificial BC construction */
1267:   PetscCall(PetscHSetICreate(&ownedpts));
1268:   PetscCall(PetscHSetICreate(&seenpts));
1269:   PetscCall(PetscHSetICreate(&owneddofs));
1270:   PetscCall(PetscHSetICreate(&seendofs));
1271:   PetscCall(PetscHSetICreate(&artificialbcs));

1273:   PetscCall(ISGetIndices(cells, &cellsArray));
1274:   PetscCall(ISGetIndices(points, &pointsArray));
1275:   PetscCall(PetscHMapICreate(&ht));
1276:   PetscCall(PetscHMapICreate(&htWithArtificial));
1277:   PetscCall(PetscHMapICreate(&htWithAll));
1278:   for (v = vStart; v < vEnd; ++v) {
1279:     PetscInt localIndex               = 0;
1280:     PetscInt localIndexWithArtificial = 0;
1281:     PetscInt localIndexWithAll        = 0;
1282:     PetscInt dof, off, i, j, k, l;

1284:     PetscCall(PetscHMapIClear(ht));
1285:     PetscCall(PetscHMapIClear(htWithArtificial));
1286:     PetscCall(PetscHMapIClear(htWithAll));
1287:     PetscCall(PetscSectionGetDof(cellCounts, v, &dof));
1288:     PetscCall(PetscSectionGetOffset(cellCounts, v, &off));
1289:     if (dof <= 0) continue;

1291:     /* Calculate the global numbers of the artificial BC dofs here first */
1292:     PetscCall(patch->patchconstructop((void *)patch, dm, v, ownedpts));
1293:     PetscCall(PCPatchCompleteCellPatch(pc, ownedpts, seenpts));
1294:     PetscCall(PCPatchGetPointDofs(pc, ownedpts, owneddofs, v, &patch->subspaces_to_exclude));
1295:     PetscCall(PCPatchGetPointDofs(pc, seenpts, seendofs, v, NULL));
1296:     PetscCall(PCPatchComputeSetDifference_Private(owneddofs, seendofs, artificialbcs));
1297:     if (patch->viewPatches) {
1298:       PetscHSetI    globalbcdofs;
1299:       PetscHashIter hi;
1300:       MPI_Comm      comm = PetscObjectComm((PetscObject)pc);

1302:       PetscCall(PetscHSetICreate(&globalbcdofs));
1303:       PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": owned dofs:\n", v));
1304:       PetscHashIterBegin(owneddofs, hi);
1305:       while (!PetscHashIterAtEnd(owneddofs, hi)) {
1306:         PetscInt globalDof;

1308:         PetscHashIterGetKey(owneddofs, hi, globalDof);
1309:         PetscHashIterNext(owneddofs, hi);
1310:         PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));
1311:       }
1312:       PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1313:       PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": seen dofs:\n", v));
1314:       PetscHashIterBegin(seendofs, hi);
1315:       while (!PetscHashIterAtEnd(seendofs, hi)) {
1316:         PetscInt  globalDof;
1317:         PetscBool flg;

1319:         PetscHashIterGetKey(seendofs, hi, globalDof);
1320:         PetscHashIterNext(seendofs, hi);
1321:         PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));

1323:         PetscCall(PetscHSetIHas(globalBcs, globalDof, &flg));
1324:         if (flg) PetscCall(PetscHSetIAdd(globalbcdofs, globalDof));
1325:       }
1326:       PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1327:       PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": global BCs:\n", v));
1328:       PetscCall(PetscHSetIGetSize(globalbcdofs, &numBcs));
1329:       if (numBcs > 0) {
1330:         PetscHashIterBegin(globalbcdofs, hi);
1331:         while (!PetscHashIterAtEnd(globalbcdofs, hi)) {
1332:           PetscInt globalDof;
1333:           PetscHashIterGetKey(globalbcdofs, hi, globalDof);
1334:           PetscHashIterNext(globalbcdofs, hi);
1335:           PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));
1336:         }
1337:       }
1338:       PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1339:       PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": artificial BCs:\n", v));
1340:       PetscCall(PetscHSetIGetSize(artificialbcs, &numBcs));
1341:       if (numBcs > 0) {
1342:         PetscHashIterBegin(artificialbcs, hi);
1343:         while (!PetscHashIterAtEnd(artificialbcs, hi)) {
1344:           PetscInt globalDof;
1345:           PetscHashIterGetKey(artificialbcs, hi, globalDof);
1346:           PetscHashIterNext(artificialbcs, hi);
1347:           PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));
1348:         }
1349:       }
1350:       PetscCall(PetscSynchronizedPrintf(comm, "\n\n"));
1351:       PetscCall(PetscHSetIDestroy(&globalbcdofs));
1352:     }
1353:     for (k = 0; k < patch->nsubspaces; ++k) {
1354:       const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
1355:       PetscInt        nodesPerCell   = patch->nodesPerCell[k];
1356:       PetscInt        subspaceOffset = patch->subspaceOffsets[k];
1357:       PetscInt        bs             = patch->bs[k];

1359:       for (i = off; i < off + dof; ++i) {
1360:         /* Walk over the cells in this patch. */
1361:         const PetscInt c    = cellsArray[i];
1362:         PetscInt       cell = c;

1364:         /* TODO Change this to an IS */
1365:         if (cellNumbering) {
1366:           PetscCall(PetscSectionGetDof(cellNumbering, c, &cell));
1367:           PetscCheck(cell > 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Cell %" PetscInt_FMT " doesn't appear in cell numbering map", c);
1368:           PetscCall(PetscSectionGetOffset(cellNumbering, c, &cell));
1369:         }
1370:         newCellsArray[i] = cell;
1371:         for (j = 0; j < nodesPerCell; ++j) {
1372:           /* For each global dof, map it into contiguous local storage. */
1373:           const PetscInt globalDof = cellNodeMap[cell * nodesPerCell + j] * bs + subspaceOffset;
1374:           /* finally, loop over block size */
1375:           for (l = 0; l < bs; ++l) {
1376:             PetscInt  localDof;
1377:             PetscBool isGlobalBcDof, isArtificialBcDof;

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

1383:             /* if it's either, don't ever give it a local dof number */
1384:             if (isGlobalBcDof || isArtificialBcDof) {
1385:               dofsArray[globalIndex] = -1; /* don't use this in assembly in this patch */
1386:             } else {
1387:               PetscCall(PetscHMapIGet(ht, globalDof + l, &localDof));
1388:               if (localDof == -1) {
1389:                 localDof = localIndex++;
1390:                 PetscCall(PetscHMapISet(ht, globalDof + l, localDof));
1391:               }
1392:               PetscCheck(globalIndex < numDofs, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %" PetscInt_FMT " than expected %" PetscInt_FMT, globalIndex + 1, numDofs);
1393:               /* And store. */
1394:               dofsArray[globalIndex] = localDof;
1395:             }

1397:             if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1398:               if (isGlobalBcDof) {
1399:                 dofsArrayWithArtificial[globalIndex] = -1; /* don't use this in assembly in this patch */
1400:               } else {
1401:                 PetscCall(PetscHMapIGet(htWithArtificial, globalDof + l, &localDof));
1402:                 if (localDof == -1) {
1403:                   localDof = localIndexWithArtificial++;
1404:                   PetscCall(PetscHMapISet(htWithArtificial, globalDof + l, localDof));
1405:                 }
1406:                 PetscCheck(globalIndex < numDofs, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %" PetscInt_FMT " than expected %" PetscInt_FMT, globalIndex + 1, numDofs);
1407:                 /* And store.*/
1408:                 dofsArrayWithArtificial[globalIndex] = localDof;
1409:               }
1410:             }

1412:             if (isNonlinear) {
1413:               /* Build the dofmap for the function space with _all_ dofs,
1414:    including those in any kind of boundary condition */
1415:               PetscCall(PetscHMapIGet(htWithAll, globalDof + l, &localDof));
1416:               if (localDof == -1) {
1417:                 localDof = localIndexWithAll++;
1418:                 PetscCall(PetscHMapISet(htWithAll, globalDof + l, localDof));
1419:               }
1420:               PetscCheck(globalIndex < numDofs, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %" PetscInt_FMT " than expected %" PetscInt_FMT, globalIndex + 1, numDofs);
1421:               /* And store.*/
1422:               dofsArrayWithAll[globalIndex] = localDof;
1423:             }
1424:             globalIndex++;
1425:           }
1426:         }
1427:       }
1428:     }
1429:     /* How many local dofs in this patch? */
1430:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1431:       PetscCall(PetscHMapIGetSize(htWithArtificial, &dof));
1432:       PetscCall(PetscSectionSetDof(gtolCountsWithArtificial, v, dof));
1433:     }
1434:     if (isNonlinear) {
1435:       PetscCall(PetscHMapIGetSize(htWithAll, &dof));
1436:       PetscCall(PetscSectionSetDof(gtolCountsWithAll, v, dof));
1437:     }
1438:     PetscCall(PetscHMapIGetSize(ht, &dof));
1439:     PetscCall(PetscSectionSetDof(gtolCounts, v, dof));
1440:   }

1442:   PetscCall(DMDestroy(&dm));
1443:   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);
1444:   PetscCall(PetscSectionSetUp(gtolCounts));
1445:   PetscCall(PetscSectionGetStorageSize(gtolCounts, &numGlobalDofs));
1446:   PetscCall(PetscMalloc1(numGlobalDofs, &globalDofsArray));

1448:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1449:     PetscCall(PetscSectionSetUp(gtolCountsWithArtificial));
1450:     PetscCall(PetscSectionGetStorageSize(gtolCountsWithArtificial, &numGlobalDofsWithArtificial));
1451:     PetscCall(PetscMalloc1(numGlobalDofsWithArtificial, &globalDofsArrayWithArtificial));
1452:   }
1453:   if (isNonlinear) {
1454:     PetscCall(PetscSectionSetUp(gtolCountsWithAll));
1455:     PetscCall(PetscSectionGetStorageSize(gtolCountsWithAll, &numGlobalDofsWithAll));
1456:     PetscCall(PetscMalloc1(numGlobalDofsWithAll, &globalDofsArrayWithAll));
1457:   }
1458:   /* Now populate the global to local map.  This could be merged into the above loop if we were willing to deal with reallocs. */
1459:   for (v = vStart; v < vEnd; ++v) {
1460:     PetscHashIter hi;
1461:     PetscInt      dof, off, Np, ooff, i, j, k, l;

1463:     PetscCall(PetscHMapIClear(ht));
1464:     PetscCall(PetscHMapIClear(htWithArtificial));
1465:     PetscCall(PetscHMapIClear(htWithAll));
1466:     PetscCall(PetscSectionGetDof(cellCounts, v, &dof));
1467:     PetscCall(PetscSectionGetOffset(cellCounts, v, &off));
1468:     PetscCall(PetscSectionGetDof(pointCounts, v, &Np));
1469:     PetscCall(PetscSectionGetOffset(pointCounts, v, &ooff));
1470:     if (dof <= 0) continue;

1472:     for (k = 0; k < patch->nsubspaces; ++k) {
1473:       const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
1474:       PetscInt        nodesPerCell   = patch->nodesPerCell[k];
1475:       PetscInt        subspaceOffset = patch->subspaceOffsets[k];
1476:       PetscInt        bs             = patch->bs[k];
1477:       PetscInt        goff;

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

1484:         if (cellNumbering) PetscCall(PetscSectionGetOffset(cellNumbering, c, &cell));
1485:         for (j = 0; j < nodesPerCell; ++j) {
1486:           for (l = 0; l < bs; ++l) {
1487:             const PetscInt globalDof = cellNodeMap[cell * nodesPerCell + j] * bs + l + subspaceOffset;
1488:             const PetscInt localDof  = dofsArray[key];
1489:             if (localDof >= 0) PetscCall(PetscHMapISet(ht, globalDof, localDof));
1490:             if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1491:               const PetscInt localDofWithArtificial = dofsArrayWithArtificial[key];
1492:               if (localDofWithArtificial >= 0) PetscCall(PetscHMapISet(htWithArtificial, globalDof, localDofWithArtificial));
1493:             }
1494:             if (isNonlinear) {
1495:               const PetscInt localDofWithAll = dofsArrayWithAll[key];
1496:               if (localDofWithAll >= 0) PetscCall(PetscHMapISet(htWithAll, globalDof, localDofWithAll));
1497:             }
1498:             key++;
1499:           }
1500:         }
1501:       }

1503:       /* Shove it in the output data structure. */
1504:       PetscCall(PetscSectionGetOffset(gtolCounts, v, &goff));
1505:       PetscHashIterBegin(ht, hi);
1506:       while (!PetscHashIterAtEnd(ht, hi)) {
1507:         PetscInt globalDof, localDof;

1509:         PetscHashIterGetKey(ht, hi, globalDof);
1510:         PetscHashIterGetVal(ht, hi, localDof);
1511:         if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof;
1512:         PetscHashIterNext(ht, hi);
1513:       }

1515:       if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1516:         PetscCall(PetscSectionGetOffset(gtolCountsWithArtificial, v, &goff));
1517:         PetscHashIterBegin(htWithArtificial, hi);
1518:         while (!PetscHashIterAtEnd(htWithArtificial, hi)) {
1519:           PetscInt globalDof, localDof;
1520:           PetscHashIterGetKey(htWithArtificial, hi, globalDof);
1521:           PetscHashIterGetVal(htWithArtificial, hi, localDof);
1522:           if (globalDof >= 0) globalDofsArrayWithArtificial[goff + localDof] = globalDof;
1523:           PetscHashIterNext(htWithArtificial, hi);
1524:         }
1525:       }
1526:       if (isNonlinear) {
1527:         PetscCall(PetscSectionGetOffset(gtolCountsWithAll, v, &goff));
1528:         PetscHashIterBegin(htWithAll, hi);
1529:         while (!PetscHashIterAtEnd(htWithAll, hi)) {
1530:           PetscInt globalDof, localDof;
1531:           PetscHashIterGetKey(htWithAll, hi, globalDof);
1532:           PetscHashIterGetVal(htWithAll, hi, localDof);
1533:           if (globalDof >= 0) globalDofsArrayWithAll[goff + localDof] = globalDof;
1534:           PetscHashIterNext(htWithAll, hi);
1535:         }
1536:       }

1538:       for (p = 0; p < Np; ++p) {
1539:         const PetscInt point = pointsArray[ooff + p];
1540:         PetscInt       globalDof, localDof;

1542:         PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof));
1543:         PetscCall(PetscHMapIGet(ht, globalDof, &localDof));
1544:         offsArray[(ooff + p) * Nf + k] = localDof;
1545:         if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1546:           PetscCall(PetscHMapIGet(htWithArtificial, globalDof, &localDof));
1547:           offsArrayWithArtificial[(ooff + p) * Nf + k] = localDof;
1548:         }
1549:         if (isNonlinear) {
1550:           PetscCall(PetscHMapIGet(htWithAll, globalDof, &localDof));
1551:           offsArrayWithAll[(ooff + p) * Nf + k] = localDof;
1552:         }
1553:       }
1554:     }

1556:     PetscCall(PetscHSetIDestroy(&globalBcs));
1557:     PetscCall(PetscHSetIDestroy(&ownedpts));
1558:     PetscCall(PetscHSetIDestroy(&seenpts));
1559:     PetscCall(PetscHSetIDestroy(&owneddofs));
1560:     PetscCall(PetscHSetIDestroy(&seendofs));
1561:     PetscCall(PetscHSetIDestroy(&artificialbcs));

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

1573:         if (cellNumbering) PetscCall(PetscSectionGetOffset(cellNumbering, c, &cell));
1574:         for (k = 0; k < patch->nsubspaces; ++k) {
1575:           const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
1576:           PetscInt        nodesPerCell   = patch->nodesPerCell[k];
1577:           PetscInt        subspaceOffset = patch->subspaceOffsets[k];
1578:           PetscInt        bs             = patch->bs[k];

1580:           for (j = 0; j < nodesPerCell; ++j) {
1581:             for (l = 0; l < bs; ++l) {
1582:               const PetscInt globalDof = cellNodeMap[cell * nodesPerCell + j] * bs + l + subspaceOffset;
1583:               PetscInt       localDof;

1585:               PetscCall(PetscHMapIGet(ht, globalDof, &localDof));
1586:               /* If it's not in the hash table, i.e. is a BC dof,
1587:    then the PetscHSetIMap above gives -1, which matches
1588:    exactly the convention for PETSc's matrix assembly to
1589:    ignore the dof. So we don't need to do anything here */
1590:               asmArray[asmKey] = localDof;
1591:               if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1592:                 PetscCall(PetscHMapIGet(htWithArtificial, globalDof, &localDof));
1593:                 asmArrayWithArtificial[asmKey] = localDof;
1594:               }
1595:               if (isNonlinear) {
1596:                 PetscCall(PetscHMapIGet(htWithAll, globalDof, &localDof));
1597:                 asmArrayWithAll[asmKey] = localDof;
1598:               }
1599:               asmKey++;
1600:             }
1601:           }
1602:         }
1603:       }
1604:     }
1605:   }
1606:   if (1 == patch->nsubspaces) {
1607:     PetscCall(PetscArraycpy(asmArray, dofsArray, numDofs));
1608:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscCall(PetscArraycpy(asmArrayWithArtificial, dofsArrayWithArtificial, numDofs));
1609:     if (isNonlinear) PetscCall(PetscArraycpy(asmArrayWithAll, dofsArrayWithAll, numDofs));
1610:   }

1612:   PetscCall(PetscHMapIDestroy(&ht));
1613:   PetscCall(PetscHMapIDestroy(&htWithArtificial));
1614:   PetscCall(PetscHMapIDestroy(&htWithAll));
1615:   PetscCall(ISRestoreIndices(cells, &cellsArray));
1616:   PetscCall(ISRestoreIndices(points, &pointsArray));
1617:   PetscCall(PetscFree(dofsArray));
1618:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscCall(PetscFree(dofsArrayWithArtificial));
1619:   if (isNonlinear) PetscCall(PetscFree(dofsArrayWithAll));
1620:   /* Create placeholder section for map from points to patch dofs */
1621:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection));
1622:   PetscCall(PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces));
1623:   if (patch->combined) {
1624:     PetscInt numFields;
1625:     PetscCall(PetscSectionGetNumFields(patch->dofSection[0], &numFields));
1626:     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);
1627:     PetscCall(PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd));
1628:     PetscCall(PetscSectionSetChart(patch->patchSection, pStart, pEnd));
1629:     for (p = pStart; p < pEnd; ++p) {
1630:       PetscInt dof, fdof, f;

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

1683: static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat, PetscBool withArtificial)
1684: {
1685:   PC_PATCH   *patch = (PC_PATCH *)pc->data;
1686:   PetscBool   flg;
1687:   PetscInt    csize, rsize;
1688:   const char *prefix = NULL;

1690:   PetscFunctionBegin;
1691:   if (withArtificial) {
1692:     /* would be nice if we could create a rectangular matrix of size numDofsWithArtificial x numDofs here */
1693:     PetscInt pStart;
1694:     PetscCall(PetscSectionGetChart(patch->gtolCountsWithArtificial, &pStart, NULL));
1695:     PetscCall(PetscSectionGetDof(patch->gtolCountsWithArtificial, point + pStart, &rsize));
1696:     csize = rsize;
1697:   } else {
1698:     PetscInt pStart;
1699:     PetscCall(PetscSectionGetChart(patch->gtolCounts, &pStart, NULL));
1700:     PetscCall(PetscSectionGetDof(patch->gtolCounts, point + pStart, &rsize));
1701:     csize = rsize;
1702:   }

1704:   PetscCall(MatCreate(PETSC_COMM_SELF, mat));
1705:   PetscCall(PCGetOptionsPrefix(pc, &prefix));
1706:   PetscCall(MatSetOptionsPrefix(*mat, prefix));
1707:   PetscCall(MatAppendOptionsPrefix(*mat, "pc_patch_sub_"));
1708:   if (patch->sub_mat_type) PetscCall(MatSetType(*mat, patch->sub_mat_type));
1709:   else if (!patch->sub_mat_type) PetscCall(MatSetType(*mat, MATDENSE));
1710:   PetscCall(MatSetSizes(*mat, rsize, csize, rsize, csize));
1711:   PetscCall(PetscObjectTypeCompare((PetscObject)*mat, MATDENSE, &flg));
1712:   if (!flg) PetscCall(PetscObjectTypeCompare((PetscObject)*mat, MATSEQDENSE, &flg));
1713:   /* Sparse patch matrices */
1714:   if (!flg) {
1715:     PetscBT         bt;
1716:     PetscInt       *dnnz      = NULL;
1717:     const PetscInt *dofsArray = NULL;
1718:     PetscInt        pStart, pEnd, ncell, offset, c, i, j;

1720:     if (withArtificial) {
1721:       PetscCall(ISGetIndices(patch->dofsWithArtificial, &dofsArray));
1722:     } else {
1723:       PetscCall(ISGetIndices(patch->dofs, &dofsArray));
1724:     }
1725:     PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));
1726:     point += pStart;
1727:     PetscCheck(point < pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %" PetscInt_FMT " not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, pStart, pEnd);
1728:     PetscCall(PetscSectionGetDof(patch->cellCounts, point, &ncell));
1729:     PetscCall(PetscSectionGetOffset(patch->cellCounts, point, &offset));
1730:     PetscCall(PetscLogEventBegin(PC_Patch_Prealloc, pc, 0, 0, 0));
1731:     /* A PetscBT uses N^2 bits to store the sparsity pattern on a
1732:    * patch. This is probably OK if the patches are not too big,
1733:    * but uses too much memory. We therefore switch based on rsize. */
1734:     if (rsize < 3000) { /* FIXME: I picked this switch value out of my hat */
1735:       PetscScalar *zeroes;
1736:       PetscInt     rows;

1738:       PetscCall(PetscCalloc1(rsize, &dnnz));
1739:       PetscCall(PetscBTCreate(rsize * rsize, &bt));
1740:       for (c = 0; c < ncell; ++c) {
1741:         const PetscInt *idx = dofsArray + (offset + c) * patch->totalDofsPerCell;
1742:         for (i = 0; i < patch->totalDofsPerCell; ++i) {
1743:           const PetscInt row = idx[i];
1744:           if (row < 0) continue;
1745:           for (j = 0; j < patch->totalDofsPerCell; ++j) {
1746:             const PetscInt col = idx[j];
1747:             const PetscInt key = row * rsize + col;
1748:             if (col < 0) continue;
1749:             if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1750:           }
1751:         }
1752:       }

1754:       if (patch->usercomputeopintfacet) {
1755:         const PetscInt *intFacetsArray = NULL;
1756:         PetscInt        i, numIntFacets, intFacetOffset;
1757:         const PetscInt *facetCells = NULL;

1759:         PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
1760:         PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
1761:         PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
1762:         PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
1763:         for (i = 0; i < numIntFacets; i++) {
1764:           const PetscInt cell0 = facetCells[2 * (intFacetOffset + i) + 0];
1765:           const PetscInt cell1 = facetCells[2 * (intFacetOffset + i) + 1];
1766:           PetscInt       celli, cellj;

1768:           for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1769:             const PetscInt row = dofsArray[(offset + cell0) * patch->totalDofsPerCell + celli];
1770:             if (row < 0) continue;
1771:             for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1772:               const PetscInt col = dofsArray[(offset + cell1) * patch->totalDofsPerCell + cellj];
1773:               const PetscInt key = row * rsize + col;
1774:               if (col < 0) continue;
1775:               if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1776:             }
1777:           }

1779:           for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1780:             const PetscInt row = dofsArray[(offset + cell1) * patch->totalDofsPerCell + celli];
1781:             if (row < 0) continue;
1782:             for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1783:               const PetscInt col = dofsArray[(offset + cell0) * patch->totalDofsPerCell + cellj];
1784:               const PetscInt key = row * rsize + col;
1785:               if (col < 0) continue;
1786:               if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1787:             }
1788:           }
1789:         }
1790:       }
1791:       PetscCall(PetscBTDestroy(&bt));
1792:       PetscCall(MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL));
1793:       PetscCall(PetscFree(dnnz));

1795:       PetscCall(PetscCalloc1(patch->totalDofsPerCell * patch->totalDofsPerCell, &zeroes));
1796:       for (c = 0; c < ncell; ++c) {
1797:         const PetscInt *idx = &dofsArray[(offset + c) * patch->totalDofsPerCell];
1798:         PetscCall(MatSetValues(*mat, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, zeroes, INSERT_VALUES));
1799:       }
1800:       PetscCall(MatGetLocalSize(*mat, &rows, NULL));
1801:       for (i = 0; i < rows; ++i) PetscCall(MatSetValues(*mat, 1, &i, 1, &i, zeroes, INSERT_VALUES));

1803:       if (patch->usercomputeopintfacet) {
1804:         const PetscInt *intFacetsArray = NULL;
1805:         PetscInt        i, numIntFacets, intFacetOffset;
1806:         const PetscInt *facetCells = NULL;

1808:         PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
1809:         PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
1810:         PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
1811:         PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
1812:         for (i = 0; i < numIntFacets; i++) {
1813:           const PetscInt  cell0    = facetCells[2 * (intFacetOffset + i) + 0];
1814:           const PetscInt  cell1    = facetCells[2 * (intFacetOffset + i) + 1];
1815:           const PetscInt *cell0idx = &dofsArray[(offset + cell0) * patch->totalDofsPerCell];
1816:           const PetscInt *cell1idx = &dofsArray[(offset + cell1) * patch->totalDofsPerCell];
1817:           PetscCall(MatSetValues(*mat, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, zeroes, INSERT_VALUES));
1818:           PetscCall(MatSetValues(*mat, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, zeroes, INSERT_VALUES));
1819:         }
1820:       }

1822:       PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY));
1823:       PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY));

1825:       PetscCall(PetscFree(zeroes));

1827:     } else { /* rsize too big, use MATPREALLOCATOR */
1828:       Mat          preallocator;
1829:       PetscScalar *vals;

1831:       PetscCall(PetscCalloc1(patch->totalDofsPerCell * patch->totalDofsPerCell, &vals));
1832:       PetscCall(MatCreate(PETSC_COMM_SELF, &preallocator));
1833:       PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
1834:       PetscCall(MatSetSizes(preallocator, rsize, rsize, rsize, rsize));
1835:       PetscCall(MatSetUp(preallocator));

1837:       for (c = 0; c < ncell; ++c) {
1838:         const PetscInt *idx = dofsArray + (offset + c) * patch->totalDofsPerCell;
1839:         PetscCall(MatSetValues(preallocator, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, vals, INSERT_VALUES));
1840:       }

1842:       if (patch->usercomputeopintfacet) {
1843:         const PetscInt *intFacetsArray = NULL;
1844:         PetscInt        i, numIntFacets, intFacetOffset;
1845:         const PetscInt *facetCells = NULL;

1847:         PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
1848:         PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
1849:         PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
1850:         PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
1851:         for (i = 0; i < numIntFacets; i++) {
1852:           const PetscInt  cell0    = facetCells[2 * (intFacetOffset + i) + 0];
1853:           const PetscInt  cell1    = facetCells[2 * (intFacetOffset + i) + 1];
1854:           const PetscInt *cell0idx = &dofsArray[(offset + cell0) * patch->totalDofsPerCell];
1855:           const PetscInt *cell1idx = &dofsArray[(offset + cell1) * patch->totalDofsPerCell];
1856:           PetscCall(MatSetValues(preallocator, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, vals, INSERT_VALUES));
1857:           PetscCall(MatSetValues(preallocator, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, vals, INSERT_VALUES));
1858:         }
1859:       }

1861:       PetscCall(PetscFree(vals));
1862:       PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
1863:       PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));
1864:       PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_TRUE, *mat));
1865:       PetscCall(MatDestroy(&preallocator));
1866:     }
1867:     PetscCall(PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0));
1868:     if (withArtificial) {
1869:       PetscCall(ISRestoreIndices(patch->dofsWithArtificial, &dofsArray));
1870:     } else {
1871:       PetscCall(ISRestoreIndices(patch->dofs, &dofsArray));
1872:     }
1873:   }
1874:   PetscCall(MatSetUp(*mat));
1875:   PetscFunctionReturn(PETSC_SUCCESS);
1876: }

1878: 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)
1879: {
1880:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
1881:   DM              dm, plex;
1882:   PetscSection    s;
1883:   const PetscInt *parray, *oarray;
1884:   PetscInt        Nf = patch->nsubspaces, Np, poff, p, f;

1886:   PetscFunctionBegin;
1887:   PetscCheck(!patch->precomputeElementTensors, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Precomputing element tensors not implemented with DMPlex compute operator");
1888:   PetscCall(PCGetDM(pc, &dm));
1889:   PetscCall(DMConvert(dm, DMPLEX, &plex));
1890:   dm = plex;
1891:   PetscCall(DMGetLocalSection(dm, &s));
1892:   /* Set offset into patch */
1893:   PetscCall(PetscSectionGetDof(patch->pointCounts, patchNum, &Np));
1894:   PetscCall(PetscSectionGetOffset(patch->pointCounts, patchNum, &poff));
1895:   PetscCall(ISGetIndices(patch->points, &parray));
1896:   PetscCall(ISGetIndices(patch->offs, &oarray));
1897:   for (f = 0; f < Nf; ++f) {
1898:     for (p = 0; p < Np; ++p) {
1899:       const PetscInt point = parray[poff + p];
1900:       PetscInt       dof;

1902:       PetscCall(PetscSectionGetFieldDof(patch->patchSection, point, f, &dof));
1903:       PetscCall(PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff + p) * Nf + f]));
1904:       if (patch->nsubspaces == 1) PetscCall(PetscSectionSetOffset(patch->patchSection, point, oarray[(poff + p) * Nf + f]));
1905:       else PetscCall(PetscSectionSetOffset(patch->patchSection, point, -1));
1906:     }
1907:   }
1908:   PetscCall(ISRestoreIndices(patch->points, &parray));
1909:   PetscCall(ISRestoreIndices(patch->offs, &oarray));
1910:   if (patch->viewSection) PetscCall(ObjectView((PetscObject)patch->patchSection, patch->viewerSection, patch->formatSection));
1911:   PetscCall(DMPlexComputeResidual_Patch_Internal(dm, patch->patchSection, cellIS, 0.0, x, NULL, F, ctx));
1912:   PetscCall(DMDestroy(&dm));
1913:   PetscFunctionReturn(PETSC_SUCCESS);
1914: }

1916: PetscErrorCode PCPatchComputeFunction_Internal(PC pc, Vec x, Vec F, PetscInt point)
1917: {
1918:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
1919:   const PetscInt *dofsArray;
1920:   const PetscInt *dofsArrayWithAll;
1921:   const PetscInt *cellsArray;
1922:   PetscInt        ncell, offset, pStart, pEnd;

1924:   PetscFunctionBegin;
1925:   PetscCall(PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0));
1926:   PetscCheck(patch->usercomputeop, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set callback");
1927:   PetscCall(ISGetIndices(patch->dofs, &dofsArray));
1928:   PetscCall(ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll));
1929:   PetscCall(ISGetIndices(patch->cells, &cellsArray));
1930:   PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));

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

1935:   PetscCall(PetscSectionGetDof(patch->cellCounts, point, &ncell));
1936:   PetscCall(PetscSectionGetOffset(patch->cellCounts, point, &offset));
1937:   if (ncell <= 0) {
1938:     PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
1939:     PetscFunctionReturn(PETSC_SUCCESS);
1940:   }
1941:   PetscCall(VecSet(F, 0.0));
1942:   /* Cannot reuse the same IS because the geometry info is being cached in it */
1943:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS));
1944:   PetscCallBack("PCPatch callback", patch->usercomputef(pc, point, x, F, patch->cellIS, ncell * patch->totalDofsPerCell, dofsArray + offset * patch->totalDofsPerCell, dofsArrayWithAll + offset * patch->totalDofsPerCell, patch->usercomputefctx));
1945:   PetscCall(ISDestroy(&patch->cellIS));
1946:   PetscCall(ISRestoreIndices(patch->dofs, &dofsArray));
1947:   PetscCall(ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll));
1948:   PetscCall(ISRestoreIndices(patch->cells, &cellsArray));
1949:   if (patch->viewMatrix) {
1950:     char name[PETSC_MAX_PATH_LEN];

1952:     PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "Patch vector for Point %" PetscInt_FMT, point));
1953:     PetscCall(PetscObjectSetName((PetscObject)F, name));
1954:     PetscCall(ObjectView((PetscObject)F, patch->viewerMatrix, patch->formatMatrix));
1955:   }
1956:   PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
1957:   PetscFunctionReturn(PETSC_SUCCESS);
1958: }

1960: 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)
1961: {
1962:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
1963:   DM              dm, plex;
1964:   PetscSection    s;
1965:   const PetscInt *parray, *oarray;
1966:   PetscInt        Nf = patch->nsubspaces, Np, poff, p, f;

1968:   PetscFunctionBegin;
1969:   PetscCall(PCGetDM(pc, &dm));
1970:   PetscCall(DMConvert(dm, DMPLEX, &plex));
1971:   dm = plex;
1972:   PetscCall(DMGetLocalSection(dm, &s));
1973:   /* Set offset into patch */
1974:   PetscCall(PetscSectionGetDof(patch->pointCounts, patchNum, &Np));
1975:   PetscCall(PetscSectionGetOffset(patch->pointCounts, patchNum, &poff));
1976:   PetscCall(ISGetIndices(patch->points, &parray));
1977:   PetscCall(ISGetIndices(patch->offs, &oarray));
1978:   for (f = 0; f < Nf; ++f) {
1979:     for (p = 0; p < Np; ++p) {
1980:       const PetscInt point = parray[poff + p];
1981:       PetscInt       dof;

1983:       PetscCall(PetscSectionGetFieldDof(patch->patchSection, point, f, &dof));
1984:       PetscCall(PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff + p) * Nf + f]));
1985:       if (patch->nsubspaces == 1) PetscCall(PetscSectionSetOffset(patch->patchSection, point, oarray[(poff + p) * Nf + f]));
1986:       else PetscCall(PetscSectionSetOffset(patch->patchSection, point, -1));
1987:     }
1988:   }
1989:   PetscCall(ISRestoreIndices(patch->points, &parray));
1990:   PetscCall(ISRestoreIndices(patch->offs, &oarray));
1991:   if (patch->viewSection) PetscCall(ObjectView((PetscObject)patch->patchSection, patch->viewerSection, patch->formatSection));
1992:   /* TODO Shut off MatViewFromOptions() in MatAssemblyEnd() here */
1993:   PetscCall(DMPlexComputeJacobian_Patch_Internal(dm, patch->patchSection, patch->patchSection, cellIS, 0.0, 0.0, x, NULL, J, J, ctx));
1994:   PetscCall(DMDestroy(&dm));
1995:   PetscFunctionReturn(PETSC_SUCCESS);
1996: }

1998: /* This function zeros mat on entry */
1999: PetscErrorCode PCPatchComputeOperator_Internal(PC pc, Vec x, Mat mat, PetscInt point, PetscBool withArtificial)
2000: {
2001:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
2002:   const PetscInt *dofsArray;
2003:   const PetscInt *dofsArrayWithAll = NULL;
2004:   const PetscInt *cellsArray;
2005:   PetscInt        ncell, offset, pStart, pEnd, numIntFacets, intFacetOffset;
2006:   PetscBool       isNonlinear;

2008:   PetscFunctionBegin;
2009:   PetscCall(PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0));
2010:   isNonlinear = patch->isNonlinear;
2011:   PetscCheck(patch->usercomputeop, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set callback");
2012:   if (withArtificial) {
2013:     PetscCall(ISGetIndices(patch->dofsWithArtificial, &dofsArray));
2014:   } else {
2015:     PetscCall(ISGetIndices(patch->dofs, &dofsArray));
2016:   }
2017:   if (isNonlinear) PetscCall(ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll));
2018:   PetscCall(ISGetIndices(patch->cells, &cellsArray));
2019:   PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));

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

2024:   PetscCall(PetscSectionGetDof(patch->cellCounts, point, &ncell));
2025:   PetscCall(PetscSectionGetOffset(patch->cellCounts, point, &offset));
2026:   if (ncell <= 0) {
2027:     PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
2028:     PetscFunctionReturn(PETSC_SUCCESS);
2029:   }
2030:   PetscCall(MatZeroEntries(mat));
2031:   if (patch->precomputeElementTensors) {
2032:     PetscInt           i;
2033:     PetscInt           ndof = patch->totalDofsPerCell;
2034:     const PetscScalar *elementTensors;

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

2066:       PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
2067:       PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
2068:       PetscCall(PCGetDM(pc, &dm));
2069:       PetscCall(DMConvert(dm, DMPLEX, &plex));
2070:       dm = plex;
2071:       PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, NULL));
2072:       /* FIXME: Pull this malloc out. */
2073:       PetscCall(PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofs));
2074:       if (dofsArrayWithAll) PetscCall(PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofsWithAll));
2075:       if (patch->precomputeElementTensors) {
2076:         PetscInt           nFacetDof = 2 * patch->totalDofsPerCell;
2077:         const PetscScalar *elementTensors;

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

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

2132:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
2133:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));

2135:   if (!(withArtificial || isNonlinear) && patch->denseinverse) {
2136:     MatFactorInfo info;
2137:     PetscBool     flg;
2138:     PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATSEQDENSE, &flg));
2139:     PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Invalid Mat type for dense inverse");
2140:     PetscCall(MatFactorInfoInitialize(&info));
2141:     PetscCall(MatLUFactor(mat, NULL, NULL, &info));
2142:     PetscCall(MatSeqDenseInvertFactors_Private(mat));
2143:   }
2144:   PetscCall(ISDestroy(&patch->cellIS));
2145:   if (withArtificial) {
2146:     PetscCall(ISRestoreIndices(patch->dofsWithArtificial, &dofsArray));
2147:   } else {
2148:     PetscCall(ISRestoreIndices(patch->dofs, &dofsArray));
2149:   }
2150:   if (isNonlinear) PetscCall(ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll));
2151:   PetscCall(ISRestoreIndices(patch->cells, &cellsArray));
2152:   if (patch->viewMatrix) {
2153:     char name[PETSC_MAX_PATH_LEN];

2155:     PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "Patch matrix for Point %" PetscInt_FMT, point));
2156:     PetscCall(PetscObjectSetName((PetscObject)mat, name));
2157:     PetscCall(ObjectView((PetscObject)mat, patch->viewerMatrix, patch->formatMatrix));
2158:   }
2159:   PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
2160:   PetscFunctionReturn(PETSC_SUCCESS);
2161: }

2163: static PetscErrorCode MatSetValues_PCPatch_Private(Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], const PetscScalar *v, InsertMode addv)
2164: {
2165:   Vec          data;
2166:   PetscScalar *array;
2167:   PetscInt     bs, nz, i, j, cell;

2169:   PetscCall(MatShellGetContext(mat, &data));
2170:   PetscCall(VecGetBlockSize(data, &bs));
2171:   PetscCall(VecGetSize(data, &nz));
2172:   PetscCall(VecGetArray(data, &array));
2173:   PetscCheck(m == n, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Only for square insertion");
2174:   cell = (PetscInt)(idxm[0] / bs); /* use the fact that this is called once per cell */
2175:   for (i = 0; i < m; i++) {
2176:     PetscCheck(idxm[i] == idxn[i], PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Row and column indices must match!");
2177:     for (j = 0; j < n; j++) {
2178:       const PetscScalar v_ = v[i * bs + j];
2179:       /* Indexing is special to the data structure we have! */
2180:       if (addv == INSERT_VALUES) {
2181:         array[cell * bs * bs + i * bs + j] = v_;
2182:       } else {
2183:         array[cell * bs * bs + i * bs + j] += v_;
2184:       }
2185:     }
2186:   }
2187:   PetscCall(VecRestoreArray(data, &array));
2188:   PetscFunctionReturn(PETSC_SUCCESS);
2189: }

2191: static PetscErrorCode PCPatchPrecomputePatchTensors_Private(PC pc)
2192: {
2193:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
2194:   const PetscInt *cellsArray;
2195:   PetscInt        ncell, offset;
2196:   const PetscInt *dofMapArray;
2197:   PetscInt        i, j;
2198:   IS              dofMap;
2199:   IS              cellIS;
2200:   const PetscInt  ndof = patch->totalDofsPerCell;
2201:   Mat             vecMat;
2202:   PetscInt        cStart, cEnd;
2203:   DM              dm, plex;

2205:   PetscCall(ISGetSize(patch->cells, &ncell));
2206:   if (!ncell) { /* No cells to assemble over -> skip */
2207:     PetscFunctionReturn(PETSC_SUCCESS);
2208:   }

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

2212:   PetscCall(PCGetDM(pc, &dm));
2213:   PetscCall(DMConvert(dm, DMPLEX, &plex));
2214:   dm = plex;
2215:   if (!patch->allCells) {
2216:     PetscHSetI    cells;
2217:     PetscHashIter hi;
2218:     PetscInt      pStart, pEnd;
2219:     PetscInt     *allCells = NULL;
2220:     PetscCall(PetscHSetICreate(&cells));
2221:     PetscCall(ISGetIndices(patch->cells, &cellsArray));
2222:     PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));
2223:     for (i = pStart; i < pEnd; i++) {
2224:       PetscCall(PetscSectionGetDof(patch->cellCounts, i, &ncell));
2225:       PetscCall(PetscSectionGetOffset(patch->cellCounts, i, &offset));
2226:       if (ncell <= 0) continue;
2227:       for (j = 0; j < ncell; j++) PetscCall(PetscHSetIAdd(cells, cellsArray[offset + j]));
2228:     }
2229:     PetscCall(ISRestoreIndices(patch->cells, &cellsArray));
2230:     PetscCall(PetscHSetIGetSize(cells, &ncell));
2231:     PetscCall(PetscMalloc1(ncell, &allCells));
2232:     PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2233:     PetscCall(PetscMalloc1(cEnd - cStart, &patch->precomputedTensorLocations));
2234:     i = 0;
2235:     PetscHashIterBegin(cells, hi);
2236:     while (!PetscHashIterAtEnd(cells, hi)) {
2237:       PetscHashIterGetKey(cells, hi, allCells[i]);
2238:       patch->precomputedTensorLocations[allCells[i]] = i;
2239:       PetscHashIterNext(cells, hi);
2240:       i++;
2241:     }
2242:     PetscCall(PetscHSetIDestroy(&cells));
2243:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncell, allCells, PETSC_OWN_POINTER, &patch->allCells));
2244:   }
2245:   PetscCall(ISGetSize(patch->allCells, &ncell));
2246:   if (!patch->cellMats) {
2247:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, ncell * ndof * ndof, &patch->cellMats));
2248:     PetscCall(VecSetBlockSize(patch->cellMats, ndof));
2249:   }
2250:   PetscCall(VecSet(patch->cellMats, 0));

2252:   PetscCall(MatCreateShell(PETSC_COMM_SELF, ncell * ndof, ncell * ndof, ncell * ndof, ncell * ndof, (void *)patch->cellMats, &vecMat));
2253:   PetscCall(MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void (*)(void)) & MatSetValues_PCPatch_Private));
2254:   PetscCall(ISGetSize(patch->allCells, &ncell));
2255:   PetscCall(ISCreateStride(PETSC_COMM_SELF, ndof * ncell, 0, 1, &dofMap));
2256:   PetscCall(ISGetIndices(dofMap, &dofMapArray));
2257:   PetscCall(ISGetIndices(patch->allCells, &cellsArray));
2258:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray, PETSC_USE_POINTER, &cellIS));
2259:   /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2260:   PetscCallBack("PCPatch callback", patch->usercomputeop(pc, -1, NULL, vecMat, cellIS, ndof * ncell, dofMapArray, NULL, patch->usercomputeopctx));
2261:   PetscCall(ISDestroy(&cellIS));
2262:   PetscCall(MatDestroy(&vecMat));
2263:   PetscCall(ISRestoreIndices(patch->allCells, &cellsArray));
2264:   PetscCall(ISRestoreIndices(dofMap, &dofMapArray));
2265:   PetscCall(ISDestroy(&dofMap));

2267:   if (patch->usercomputeopintfacet) {
2268:     PetscInt        nIntFacets;
2269:     IS              intFacetsIS;
2270:     const PetscInt *intFacetsArray = NULL;
2271:     if (!patch->allIntFacets) {
2272:       PetscHSetI    facets;
2273:       PetscHashIter hi;
2274:       PetscInt      pStart, pEnd, fStart, fEnd;
2275:       PetscInt     *allIntFacets = NULL;
2276:       PetscCall(PetscHSetICreate(&facets));
2277:       PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
2278:       PetscCall(PetscSectionGetChart(patch->intFacetCounts, &pStart, &pEnd));
2279:       PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
2280:       for (i = pStart; i < pEnd; i++) {
2281:         PetscCall(PetscSectionGetDof(patch->intFacetCounts, i, &nIntFacets));
2282:         PetscCall(PetscSectionGetOffset(patch->intFacetCounts, i, &offset));
2283:         if (nIntFacets <= 0) continue;
2284:         for (j = 0; j < nIntFacets; j++) PetscCall(PetscHSetIAdd(facets, intFacetsArray[offset + j]));
2285:       }
2286:       PetscCall(ISRestoreIndices(patch->intFacets, &intFacetsArray));
2287:       PetscCall(PetscHSetIGetSize(facets, &nIntFacets));
2288:       PetscCall(PetscMalloc1(nIntFacets, &allIntFacets));
2289:       PetscCall(PetscMalloc1(fEnd - fStart, &patch->precomputedIntFacetTensorLocations));
2290:       i = 0;
2291:       PetscHashIterBegin(facets, hi);
2292:       while (!PetscHashIterAtEnd(facets, hi)) {
2293:         PetscHashIterGetKey(facets, hi, allIntFacets[i]);
2294:         patch->precomputedIntFacetTensorLocations[allIntFacets[i] - fStart] = i;
2295:         PetscHashIterNext(facets, hi);
2296:         i++;
2297:       }
2298:       PetscCall(PetscHSetIDestroy(&facets));
2299:       PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, allIntFacets, PETSC_OWN_POINTER, &patch->allIntFacets));
2300:     }
2301:     PetscCall(ISGetSize(patch->allIntFacets, &nIntFacets));
2302:     if (!patch->intFacetMats) {
2303:       PetscCall(VecCreateSeq(PETSC_COMM_SELF, nIntFacets * ndof * ndof * 4, &patch->intFacetMats));
2304:       PetscCall(VecSetBlockSize(patch->intFacetMats, ndof * 2));
2305:     }
2306:     PetscCall(VecSet(patch->intFacetMats, 0));

2308:     PetscCall(MatCreateShell(PETSC_COMM_SELF, nIntFacets * ndof * 2, nIntFacets * ndof * 2, nIntFacets * ndof * 2, nIntFacets * ndof * 2, (void *)patch->intFacetMats, &vecMat));
2309:     PetscCall(MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void (*)(void)) & MatSetValues_PCPatch_Private));
2310:     PetscCall(ISCreateStride(PETSC_COMM_SELF, 2 * ndof * nIntFacets, 0, 1, &dofMap));
2311:     PetscCall(ISGetIndices(dofMap, &dofMapArray));
2312:     PetscCall(ISGetIndices(patch->allIntFacets, &intFacetsArray));
2313:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, intFacetsArray, PETSC_USE_POINTER, &intFacetsIS));
2314:     /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2315:     PetscCallBack("PCPatch callback (interior facets)", patch->usercomputeopintfacet(pc, -1, NULL, vecMat, intFacetsIS, 2 * ndof * nIntFacets, dofMapArray, NULL, patch->usercomputeopintfacetctx));
2316:     PetscCall(ISDestroy(&intFacetsIS));
2317:     PetscCall(MatDestroy(&vecMat));
2318:     PetscCall(ISRestoreIndices(patch->allIntFacets, &intFacetsArray));
2319:     PetscCall(ISRestoreIndices(dofMap, &dofMapArray));
2320:     PetscCall(ISDestroy(&dofMap));
2321:   }
2322:   PetscCall(DMDestroy(&dm));
2323:   PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
2324:   PetscFunctionReturn(PETSC_SUCCESS);
2325: }

2327: PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat, PatchScatterType scattertype)
2328: {
2329:   PC_PATCH          *patch     = (PC_PATCH *)pc->data;
2330:   const PetscScalar *xArray    = NULL;
2331:   PetscScalar       *yArray    = NULL;
2332:   const PetscInt    *gtolArray = NULL;
2333:   PetscInt           dof, offset, lidx;

2335:   PetscFunctionBeginHot;
2336:   PetscCall(VecGetArrayRead(x, &xArray));
2337:   PetscCall(VecGetArray(y, &yArray));
2338:   if (scattertype == SCATTER_WITHARTIFICIAL) {
2339:     PetscCall(PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof));
2340:     PetscCall(PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offset));
2341:     PetscCall(ISGetIndices(patch->gtolWithArtificial, &gtolArray));
2342:   } else if (scattertype == SCATTER_WITHALL) {
2343:     PetscCall(PetscSectionGetDof(patch->gtolCountsWithAll, p, &dof));
2344:     PetscCall(PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offset));
2345:     PetscCall(ISGetIndices(patch->gtolWithAll, &gtolArray));
2346:   } else {
2347:     PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &dof));
2348:     PetscCall(PetscSectionGetOffset(patch->gtolCounts, p, &offset));
2349:     PetscCall(ISGetIndices(patch->gtol, &gtolArray));
2350:   }
2351:   PetscCheck(mode != INSERT_VALUES || scat == SCATTER_FORWARD, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't insert if not scattering forward");
2352:   PetscCheck(mode != ADD_VALUES || scat == SCATTER_REVERSE, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't add if not scattering reverse");
2353:   for (lidx = 0; lidx < dof; ++lidx) {
2354:     const PetscInt gidx = gtolArray[offset + lidx];

2356:     if (mode == INSERT_VALUES) yArray[lidx] = xArray[gidx]; /* Forward */
2357:     else yArray[gidx] += xArray[lidx];                      /* Reverse */
2358:   }
2359:   if (scattertype == SCATTER_WITHARTIFICIAL) {
2360:     PetscCall(ISRestoreIndices(patch->gtolWithArtificial, &gtolArray));
2361:   } else if (scattertype == SCATTER_WITHALL) {
2362:     PetscCall(ISRestoreIndices(patch->gtolWithAll, &gtolArray));
2363:   } else {
2364:     PetscCall(ISRestoreIndices(patch->gtol, &gtolArray));
2365:   }
2366:   PetscCall(VecRestoreArrayRead(x, &xArray));
2367:   PetscCall(VecRestoreArray(y, &yArray));
2368:   PetscFunctionReturn(PETSC_SUCCESS);
2369: }

2371: static PetscErrorCode PCSetUp_PATCH_Linear(PC pc)
2372: {
2373:   PC_PATCH   *patch = (PC_PATCH *)pc->data;
2374:   const char *prefix;
2375:   PetscInt    i;

2377:   PetscFunctionBegin;
2378:   if (!pc->setupcalled) {
2379:     PetscCheck(patch->save_operators || !patch->denseinverse, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Can't have dense inverse without save operators");
2380:     if (!patch->denseinverse) {
2381:       PetscCall(PetscMalloc1(patch->npatch, &patch->solver));
2382:       PetscCall(PCGetOptionsPrefix(pc, &prefix));
2383:       for (i = 0; i < patch->npatch; ++i) {
2384:         KSP ksp;
2385:         PC  subpc;

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

2422:       PetscCall(MatGetSize(patch->mat[i], &dof, NULL));
2423:       if (dof == 0) {
2424:         patch->matWithArtificial[i] = NULL;
2425:         continue;
2426:       }

2428:       PetscCall(PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE));
2429:       PetscCall(PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE));

2431:       PetscCall(MatGetSize(matSquare, &dof, NULL));
2432:       PetscCall(ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis));
2433:       if (pc->setupcalled) {
2434:         PetscCall(MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_REUSE_MATRIX, &patch->matWithArtificial[i]));
2435:       } else {
2436:         PetscCall(MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &patch->matWithArtificial[i]));
2437:       }
2438:       PetscCall(ISDestroy(&rowis));
2439:       PetscCall(MatDestroy(&matSquare));
2440:     }
2441:   }
2442:   PetscFunctionReturn(PETSC_SUCCESS);
2443: }

2445: static PetscErrorCode PCSetUp_PATCH(PC pc)
2446: {
2447:   PC_PATCH *patch = (PC_PATCH *)pc->data;
2448:   PetscInt  i;
2449:   PetscBool isNonlinear;
2450:   PetscInt  maxDof = -1, maxDofWithArtificial = -1;

2452:   PetscFunctionBegin;
2453:   if (!pc->setupcalled) {
2454:     PetscInt pStart, pEnd, p;
2455:     PetscInt localSize;

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

2459:     isNonlinear = patch->isNonlinear;
2460:     if (!patch->nsubspaces) {
2461:       DM           dm, plex;
2462:       PetscSection s;
2463:       PetscInt     cStart, cEnd, c, Nf, f, numGlobalBcs = 0, *globalBcs, *Nb, **cellDofs;

2465:       PetscCall(PCGetDM(pc, &dm));
2466:       PetscCheck(dm, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Must set DM for PCPATCH or call PCPatchSetDiscretisationInfo()");
2467:       PetscCall(DMConvert(dm, DMPLEX, &plex));
2468:       dm = plex;
2469:       PetscCall(DMGetLocalSection(dm, &s));
2470:       PetscCall(PetscSectionGetNumFields(s, &Nf));
2471:       PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2472:       for (p = pStart; p < pEnd; ++p) {
2473:         PetscInt cdof;
2474:         PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2475:         numGlobalBcs += cdof;
2476:       }
2477:       PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2478:       PetscCall(PetscMalloc3(Nf, &Nb, Nf, &cellDofs, numGlobalBcs, &globalBcs));
2479:       for (f = 0; f < Nf; ++f) {
2480:         PetscFE        fe;
2481:         PetscDualSpace sp;
2482:         PetscInt       cdoff = 0;

2484:         PetscCall(DMGetField(dm, f, NULL, (PetscObject *)&fe));
2485:         /* PetscCall(PetscFEGetNumComponents(fe, &Nc[f])); */
2486:         PetscCall(PetscFEGetDualSpace(fe, &sp));
2487:         PetscCall(PetscDualSpaceGetDimension(sp, &Nb[f]));

2489:         PetscCall(PetscMalloc1((cEnd - cStart) * Nb[f], &cellDofs[f]));
2490:         for (c = cStart; c < cEnd; ++c) {
2491:           PetscInt *closure = NULL;
2492:           PetscInt  clSize  = 0, cl;

2494:           PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
2495:           for (cl = 0; cl < clSize * 2; cl += 2) {
2496:             const PetscInt p = closure[cl];
2497:             PetscInt       fdof, d, foff;

2499:             PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
2500:             PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
2501:             for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d;
2502:           }
2503:           PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
2504:         }
2505:         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]);
2506:       }
2507:       numGlobalBcs = 0;
2508:       for (p = pStart; p < pEnd; ++p) {
2509:         const PetscInt *ind;
2510:         PetscInt        off, cdof, d;

2512:         PetscCall(PetscSectionGetOffset(s, p, &off));
2513:         PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2514:         PetscCall(PetscSectionGetConstraintIndices(s, p, &ind));
2515:         for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d];
2516:       }

2518:       PetscCall(PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **)cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs));
2519:       for (f = 0; f < Nf; ++f) PetscCall(PetscFree(cellDofs[f]));
2520:       PetscCall(PetscFree3(Nb, cellDofs, globalBcs));
2521:       PetscCall(PCPatchSetComputeFunction(pc, PCPatchComputeFunction_DMPlex_Private, NULL));
2522:       PetscCall(PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL));
2523:       PetscCall(DMDestroy(&dm));
2524:     }

2526:     localSize = patch->subspaceOffsets[patch->nsubspaces];
2527:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localRHS));
2528:     PetscCall(VecSetUp(patch->localRHS));
2529:     PetscCall(VecDuplicate(patch->localRHS, &patch->localUpdate));
2530:     PetscCall(PCPatchCreateCellPatches(pc));
2531:     PetscCall(PCPatchCreateCellPatchDiscretisationInfo(pc));

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

2536:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscCall(PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithArtificial));
2537:     if (isNonlinear) PetscCall(PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithAll));
2538:     for (p = pStart; p < pEnd; ++p) {
2539:       PetscInt dof;

2541:       PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &dof));
2542:       maxDof = PetscMax(maxDof, dof);
2543:       if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2544:         const PetscInt *gtolArray, *gtolArrayWithArtificial = NULL;
2545:         PetscInt        numPatchDofs, offset;
2546:         PetscInt        numPatchDofsWithArtificial, offsetWithArtificial;
2547:         PetscInt        dofWithoutArtificialCounter = 0;
2548:         PetscInt       *patchWithoutArtificialToWithArtificialArray;

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

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

2557:         PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs));
2558:         if (numPatchDofs == 0) {
2559:           patch->dofMappingWithoutToWithArtificial[p - pStart] = NULL;
2560:           continue;
2561:         }

2563:         PetscCall(PetscSectionGetOffset(patch->gtolCounts, p, &offset));
2564:         PetscCall(ISGetIndices(patch->gtolWithArtificial, &gtolArrayWithArtificial));
2565:         PetscCall(PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &numPatchDofsWithArtificial));
2566:         PetscCall(PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offsetWithArtificial));

2568:         PetscCall(PetscMalloc1(numPatchDofs, &patchWithoutArtificialToWithArtificialArray));
2569:         for (i = 0; i < numPatchDofsWithArtificial; i++) {
2570:           if (gtolArrayWithArtificial[i + offsetWithArtificial] == gtolArray[offset + dofWithoutArtificialCounter]) {
2571:             patchWithoutArtificialToWithArtificialArray[dofWithoutArtificialCounter] = i;
2572:             dofWithoutArtificialCounter++;
2573:             if (dofWithoutArtificialCounter == numPatchDofs) break;
2574:           }
2575:         }
2576:         PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutArtificialToWithArtificialArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithArtificial[p - pStart]));
2577:         PetscCall(ISRestoreIndices(patch->gtol, &gtolArray));
2578:         PetscCall(ISRestoreIndices(patch->gtolWithArtificial, &gtolArrayWithArtificial));
2579:       }
2580:     }
2581:     for (p = pStart; p < pEnd; ++p) {
2582:       if (isNonlinear) {
2583:         const PetscInt *gtolArray, *gtolArrayWithAll = NULL;
2584:         PetscInt        numPatchDofs, offset;
2585:         PetscInt        numPatchDofsWithAll, offsetWithAll;
2586:         PetscInt        dofWithoutAllCounter = 0;
2587:         PetscInt       *patchWithoutAllToWithAllArray;

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

2593:         PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs));
2594:         if (numPatchDofs == 0) {
2595:           patch->dofMappingWithoutToWithAll[p - pStart] = NULL;
2596:           continue;
2597:         }

2599:         PetscCall(PetscSectionGetOffset(patch->gtolCounts, p, &offset));
2600:         PetscCall(ISGetIndices(patch->gtolWithAll, &gtolArrayWithAll));
2601:         PetscCall(PetscSectionGetDof(patch->gtolCountsWithAll, p, &numPatchDofsWithAll));
2602:         PetscCall(PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offsetWithAll));

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

2606:         for (i = 0; i < numPatchDofsWithAll; i++) {
2607:           if (gtolArrayWithAll[i + offsetWithAll] == gtolArray[offset + dofWithoutAllCounter]) {
2608:             patchWithoutAllToWithAllArray[dofWithoutAllCounter] = i;
2609:             dofWithoutAllCounter++;
2610:             if (dofWithoutAllCounter == numPatchDofs) break;
2611:           }
2612:         }
2613:         PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutAllToWithAllArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithAll[p - pStart]));
2614:         PetscCall(ISRestoreIndices(patch->gtol, &gtolArray));
2615:         PetscCall(ISRestoreIndices(patch->gtolWithAll, &gtolArrayWithAll));
2616:       }
2617:     }
2618:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2619:       PetscCall(VecCreateSeq(PETSC_COMM_SELF, maxDofWithArtificial, &patch->patchRHSWithArtificial));
2620:       PetscCall(VecSetUp(patch->patchRHSWithArtificial));
2621:     }
2622:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchRHS));
2623:     PetscCall(VecSetUp(patch->patchRHS));
2624:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchUpdate));
2625:     PetscCall(VecSetUp(patch->patchUpdate));
2626:     if (patch->save_operators) {
2627:       PetscCall(PetscMalloc1(patch->npatch, &patch->mat));
2628:       for (i = 0; i < patch->npatch; ++i) PetscCall(PCPatchCreateMatrix_Private(pc, i, &patch->mat[i], PETSC_FALSE));
2629:     }
2630:     PetscCall(PetscLogEventEnd(PC_Patch_CreatePatches, pc, 0, 0, 0));

2632:     /* If desired, calculate weights for dof multiplicity */
2633:     if (patch->partition_of_unity) {
2634:       PetscScalar *input  = NULL;
2635:       PetscScalar *output = NULL;
2636:       Vec          global;

2638:       PetscCall(VecDuplicate(patch->localRHS, &patch->dof_weights));
2639:       if (patch->local_composition_type == PC_COMPOSITE_ADDITIVE) {
2640:         for (i = 0; i < patch->npatch; ++i) {
2641:           PetscInt dof;

2643:           PetscCall(PetscSectionGetDof(patch->gtolCounts, i + pStart, &dof));
2644:           if (dof <= 0) continue;
2645:           PetscCall(VecSet(patch->patchRHS, 1.0));
2646:           PetscCall(PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHS, patch->dof_weights, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR));
2647:         }
2648:       } else {
2649:         /* multiplicative is actually only locally multiplicative and globally additive. need the pou where the mesh decomposition overlaps */
2650:         PetscCall(VecSet(patch->dof_weights, 1.0));
2651:       }

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

2656:       PetscCall(VecGetArray(patch->dof_weights, &input));
2657:       PetscCall(VecGetArray(global, &output));
2658:       PetscCall(PetscSFReduceBegin(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM));
2659:       PetscCall(PetscSFReduceEnd(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM));
2660:       PetscCall(VecRestoreArray(patch->dof_weights, &input));
2661:       PetscCall(VecRestoreArray(global, &output));

2663:       PetscCall(VecReciprocal(global));

2665:       PetscCall(VecGetArray(patch->dof_weights, &output));
2666:       PetscCall(VecGetArray(global, &input));
2667:       PetscCall(PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, input, output, MPI_REPLACE));
2668:       PetscCall(PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, input, output, MPI_REPLACE));
2669:       PetscCall(VecRestoreArray(patch->dof_weights, &output));
2670:       PetscCall(VecRestoreArray(global, &input));
2671:       PetscCall(VecDestroy(&global));
2672:     }
2673:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE && patch->save_operators && !patch->isNonlinear) PetscCall(PetscMalloc1(patch->npatch, &patch->matWithArtificial));
2674:   }
2675:   PetscCall((*patch->setupsolver)(pc));
2676:   PetscFunctionReturn(PETSC_SUCCESS);
2677: }

2679: static PetscErrorCode PCApply_PATCH_Linear(PC pc, PetscInt i, Vec x, Vec y)
2680: {
2681:   PC_PATCH *patch = (PC_PATCH *)pc->data;
2682:   KSP       ksp;
2683:   Mat       op;
2684:   PetscInt  m, n;

2686:   PetscFunctionBegin;
2687:   if (patch->denseinverse) {
2688:     PetscCall((*patch->densesolve)(patch->mat[i], x, y));
2689:     PetscFunctionReturn(PETSC_SUCCESS);
2690:   }
2691:   ksp = (KSP)patch->solver[i];
2692:   if (!patch->save_operators) {
2693:     Mat mat;

2695:     PetscCall(PCPatchCreateMatrix_Private(pc, i, &mat, PETSC_FALSE));
2696:     /* Populate operator here. */
2697:     PetscCall(PCPatchComputeOperator_Internal(pc, NULL, mat, i, PETSC_FALSE));
2698:     PetscCall(KSPSetOperators(ksp, mat, mat));
2699:     /* Drop reference so the KSPSetOperators below will blow it away. */
2700:     PetscCall(MatDestroy(&mat));
2701:   }
2702:   PetscCall(PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0));
2703:   if (!ksp->setfromoptionscalled) PetscCall(KSPSetFromOptions(ksp));
2704:   /* Disgusting trick to reuse work vectors */
2705:   PetscCall(KSPGetOperators(ksp, &op, NULL));
2706:   PetscCall(MatGetLocalSize(op, &m, &n));
2707:   x->map->n = m;
2708:   y->map->n = n;
2709:   x->map->N = m;
2710:   y->map->N = n;
2711:   PetscCall(KSPSolve(ksp, x, y));
2712:   PetscCall(KSPCheckSolve(ksp, pc, y));
2713:   PetscCall(PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0));
2714:   if (!patch->save_operators) {
2715:     PC pc;
2716:     PetscCall(KSPSetOperators(ksp, NULL, NULL));
2717:     PetscCall(KSPGetPC(ksp, &pc));
2718:     /* Destroy PC context too, otherwise the factored matrix hangs around. */
2719:     PetscCall(PCReset(pc));
2720:   }
2721:   PetscFunctionReturn(PETSC_SUCCESS);
2722: }

2724: static PetscErrorCode PCUpdateMultiplicative_PATCH_Linear(PC pc, PetscInt i, PetscInt pStart)
2725: {
2726:   PC_PATCH *patch = (PC_PATCH *)pc->data;
2727:   Mat       multMat;
2728:   PetscInt  n, m;

2730:   PetscFunctionBegin;
2731:   if (patch->save_operators) {
2732:     multMat = patch->matWithArtificial[i];
2733:   } else {
2734:     /*Very inefficient, hopefully we can just assemble the rectangular matrix in the first place.*/
2735:     Mat      matSquare;
2736:     PetscInt dof;
2737:     IS       rowis;
2738:     PetscCall(PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE));
2739:     PetscCall(PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE));
2740:     PetscCall(MatGetSize(matSquare, &dof, NULL));
2741:     PetscCall(ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis));
2742:     PetscCall(MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &multMat));
2743:     PetscCall(MatDestroy(&matSquare));
2744:     PetscCall(ISDestroy(&rowis));
2745:   }
2746:   /* Disgusting trick to reuse work vectors */
2747:   PetscCall(MatGetLocalSize(multMat, &m, &n));
2748:   patch->patchUpdate->map->n            = n;
2749:   patch->patchRHSWithArtificial->map->n = m;
2750:   patch->patchUpdate->map->N            = n;
2751:   patch->patchRHSWithArtificial->map->N = m;
2752:   PetscCall(MatMult(multMat, patch->patchUpdate, patch->patchRHSWithArtificial));
2753:   PetscCall(VecScale(patch->patchRHSWithArtificial, -1.0));
2754:   PetscCall(PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHSWithArtificial, patch->localRHS, ADD_VALUES, SCATTER_REVERSE, SCATTER_WITHARTIFICIAL));
2755:   if (!patch->save_operators) PetscCall(MatDestroy(&multMat));
2756:   PetscFunctionReturn(PETSC_SUCCESS);
2757: }

2759: static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y)
2760: {
2761:   PC_PATCH          *patch        = (PC_PATCH *)pc->data;
2762:   const PetscScalar *globalRHS    = NULL;
2763:   PetscScalar       *localRHS     = NULL;
2764:   PetscScalar       *globalUpdate = NULL;
2765:   const PetscInt    *bcNodes      = NULL;
2766:   PetscInt           nsweep       = patch->symmetrise_sweep ? 2 : 1;
2767:   PetscInt           start[2]     = {0, 0};
2768:   PetscInt           end[2]       = {-1, -1};
2769:   const PetscInt     inc[2]       = {1, -1};
2770:   const PetscScalar *localUpdate;
2771:   const PetscInt    *iterationSet;
2772:   PetscInt           pStart, numBcs, n, sweep, bc, j;

2774:   PetscFunctionBegin;
2775:   PetscCall(PetscLogEventBegin(PC_Patch_Apply, pc, 0, 0, 0));
2776:   PetscCall(PetscOptionsPushGetViewerOff(PETSC_TRUE));
2777:   /* start, end, inc have 2 entries to manage a second backward sweep if we symmetrize */
2778:   end[0]   = patch->npatch;
2779:   start[1] = patch->npatch - 1;
2780:   if (patch->user_patches) {
2781:     PetscCall(ISGetLocalSize(patch->iterationSet, &end[0]));
2782:     start[1] = end[0] - 1;
2783:     PetscCall(ISGetIndices(patch->iterationSet, &iterationSet));
2784:   }
2785:   /* Scatter from global space into overlapped local spaces */
2786:   PetscCall(VecGetArrayRead(x, &globalRHS));
2787:   PetscCall(VecGetArray(patch->localRHS, &localRHS));
2788:   PetscCall(PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, globalRHS, localRHS, MPI_REPLACE));
2789:   PetscCall(PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, globalRHS, localRHS, MPI_REPLACE));
2790:   PetscCall(VecRestoreArrayRead(x, &globalRHS));
2791:   PetscCall(VecRestoreArray(patch->localRHS, &localRHS));

2793:   PetscCall(VecSet(patch->localUpdate, 0.0));
2794:   PetscCall(PetscSectionGetChart(patch->gtolCounts, &pStart, NULL));
2795:   PetscCall(PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0));
2796:   for (sweep = 0; sweep < nsweep; sweep++) {
2797:     for (j = start[sweep]; j * inc[sweep] < end[sweep] * inc[sweep]; j += inc[sweep]) {
2798:       PetscInt i = patch->user_patches ? iterationSet[j] : j;
2799:       PetscInt start, len;

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

2824:   /* Now we need to send the global BC values through */
2825:   PetscCall(VecGetArrayRead(x, &globalRHS));
2826:   PetscCall(ISGetSize(patch->globalBcNodes, &numBcs));
2827:   PetscCall(ISGetIndices(patch->globalBcNodes, &bcNodes));
2828:   PetscCall(VecGetLocalSize(x, &n));
2829:   for (bc = 0; bc < numBcs; ++bc) {
2830:     const PetscInt idx = bcNodes[bc];
2831:     if (idx < n) globalUpdate[idx] = globalRHS[idx];
2832:   }

2834:   PetscCall(ISRestoreIndices(patch->globalBcNodes, &bcNodes));
2835:   PetscCall(VecRestoreArrayRead(x, &globalRHS));
2836:   PetscCall(VecRestoreArray(y, &globalUpdate));

2838:   PetscCall(PetscOptionsPopGetViewerOff());
2839:   PetscCall(PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0));
2840:   PetscFunctionReturn(PETSC_SUCCESS);
2841: }

2843: static PetscErrorCode PCReset_PATCH_Linear(PC pc)
2844: {
2845:   PC_PATCH *patch = (PC_PATCH *)pc->data;
2846:   PetscInt  i;

2848:   PetscFunctionBegin;
2849:   if (patch->solver) {
2850:     for (i = 0; i < patch->npatch; ++i) PetscCall(KSPReset((KSP)patch->solver[i]));
2851:   }
2852:   PetscFunctionReturn(PETSC_SUCCESS);
2853: }

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

2860:   PetscFunctionBegin;
2861:   PetscCall(PetscSFDestroy(&patch->sectionSF));
2862:   PetscCall(PetscSectionDestroy(&patch->cellCounts));
2863:   PetscCall(PetscSectionDestroy(&patch->pointCounts));
2864:   PetscCall(PetscSectionDestroy(&patch->cellNumbering));
2865:   PetscCall(PetscSectionDestroy(&patch->gtolCounts));
2866:   PetscCall(ISDestroy(&patch->gtol));
2867:   PetscCall(ISDestroy(&patch->cells));
2868:   PetscCall(ISDestroy(&patch->points));
2869:   PetscCall(ISDestroy(&patch->dofs));
2870:   PetscCall(ISDestroy(&patch->offs));
2871:   PetscCall(PetscSectionDestroy(&patch->patchSection));
2872:   PetscCall(ISDestroy(&patch->ghostBcNodes));
2873:   PetscCall(ISDestroy(&patch->globalBcNodes));
2874:   PetscCall(PetscSectionDestroy(&patch->gtolCountsWithArtificial));
2875:   PetscCall(ISDestroy(&patch->gtolWithArtificial));
2876:   PetscCall(ISDestroy(&patch->dofsWithArtificial));
2877:   PetscCall(ISDestroy(&patch->offsWithArtificial));
2878:   PetscCall(PetscSectionDestroy(&patch->gtolCountsWithAll));
2879:   PetscCall(ISDestroy(&patch->gtolWithAll));
2880:   PetscCall(ISDestroy(&patch->dofsWithAll));
2881:   PetscCall(ISDestroy(&patch->offsWithAll));
2882:   PetscCall(VecDestroy(&patch->cellMats));
2883:   PetscCall(VecDestroy(&patch->intFacetMats));
2884:   PetscCall(ISDestroy(&patch->allCells));
2885:   PetscCall(ISDestroy(&patch->intFacets));
2886:   PetscCall(ISDestroy(&patch->extFacets));
2887:   PetscCall(ISDestroy(&patch->intFacetsToPatchCell));
2888:   PetscCall(PetscSectionDestroy(&patch->intFacetCounts));
2889:   PetscCall(PetscSectionDestroy(&patch->extFacetCounts));

2891:   if (patch->dofSection)
2892:     for (i = 0; i < patch->nsubspaces; i++) PetscCall(PetscSectionDestroy(&patch->dofSection[i]));
2893:   PetscCall(PetscFree(patch->dofSection));
2894:   PetscCall(PetscFree(patch->bs));
2895:   PetscCall(PetscFree(patch->nodesPerCell));
2896:   if (patch->cellNodeMap)
2897:     for (i = 0; i < patch->nsubspaces; i++) PetscCall(PetscFree(patch->cellNodeMap[i]));
2898:   PetscCall(PetscFree(patch->cellNodeMap));
2899:   PetscCall(PetscFree(patch->subspaceOffsets));

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

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

2905:   PetscCall(VecDestroy(&patch->localRHS));
2906:   PetscCall(VecDestroy(&patch->localUpdate));
2907:   PetscCall(VecDestroy(&patch->patchRHS));
2908:   PetscCall(VecDestroy(&patch->patchUpdate));
2909:   PetscCall(VecDestroy(&patch->dof_weights));
2910:   if (patch->patch_dof_weights) {
2911:     for (i = 0; i < patch->npatch; ++i) PetscCall(VecDestroy(&patch->patch_dof_weights[i]));
2912:     PetscCall(PetscFree(patch->patch_dof_weights));
2913:   }
2914:   if (patch->mat) {
2915:     for (i = 0; i < patch->npatch; ++i) PetscCall(MatDestroy(&patch->mat[i]));
2916:     PetscCall(PetscFree(patch->mat));
2917:   }
2918:   if (patch->matWithArtificial && !patch->isNonlinear) {
2919:     for (i = 0; i < patch->npatch; ++i) PetscCall(MatDestroy(&patch->matWithArtificial[i]));
2920:     PetscCall(PetscFree(patch->matWithArtificial));
2921:   }
2922:   PetscCall(VecDestroy(&patch->patchRHSWithArtificial));
2923:   if (patch->dofMappingWithoutToWithArtificial) {
2924:     for (i = 0; i < patch->npatch; ++i) PetscCall(ISDestroy(&patch->dofMappingWithoutToWithArtificial[i]));
2925:     PetscCall(PetscFree(patch->dofMappingWithoutToWithArtificial));
2926:   }
2927:   if (patch->dofMappingWithoutToWithAll) {
2928:     for (i = 0; i < patch->npatch; ++i) PetscCall(ISDestroy(&patch->dofMappingWithoutToWithAll[i]));
2929:     PetscCall(PetscFree(patch->dofMappingWithoutToWithAll));
2930:   }
2931:   PetscCall(PetscFree(patch->sub_mat_type));
2932:   if (patch->userIS) {
2933:     for (i = 0; i < patch->npatch; ++i) PetscCall(ISDestroy(&patch->userIS[i]));
2934:     PetscCall(PetscFree(patch->userIS));
2935:   }
2936:   PetscCall(PetscFree(patch->precomputedTensorLocations));
2937:   PetscCall(PetscFree(patch->precomputedIntFacetTensorLocations));

2939:   patch->bs          = NULL;
2940:   patch->cellNodeMap = NULL;
2941:   patch->nsubspaces  = 0;
2942:   PetscCall(ISDestroy(&patch->iterationSet));

2944:   PetscCall(PetscOptionsRestoreViewer(&patch->viewerCells));
2945:   PetscCall(PetscOptionsRestoreViewer(&patch->viewerIntFacets));
2946:   PetscCall(PetscOptionsRestoreViewer(&patch->viewerPoints));
2947:   PetscCall(PetscOptionsRestoreViewer(&patch->viewerSection));
2948:   PetscCall(PetscOptionsRestoreViewer(&patch->viewerMatrix));
2949:   PetscFunctionReturn(PETSC_SUCCESS);
2950: }

2952: static PetscErrorCode PCDestroy_PATCH_Linear(PC pc)
2953: {
2954:   PC_PATCH *patch = (PC_PATCH *)pc->data;
2955:   PetscInt  i;

2957:   PetscFunctionBegin;
2958:   if (patch->solver) {
2959:     for (i = 0; i < patch->npatch; ++i) PetscCall(KSPDestroy((KSP *)&patch->solver[i]));
2960:     PetscCall(PetscFree(patch->solver));
2961:   }
2962:   PetscFunctionReturn(PETSC_SUCCESS);
2963: }

2965: static PetscErrorCode PCDestroy_PATCH(PC pc)
2966: {
2967:   PC_PATCH *patch = (PC_PATCH *)pc->data;

2969:   PetscFunctionBegin;
2970:   PetscCall(PCReset_PATCH(pc));
2971:   PetscCall((*patch->destroysolver)(pc));
2972:   PetscCall(PetscFree(pc->data));
2973:   PetscFunctionReturn(PETSC_SUCCESS);
2974: }

2976: static PetscErrorCode PCSetFromOptions_PATCH(PC pc, PetscOptionItems *PetscOptionsObject)
2977: {
2978:   PC_PATCH            *patch                 = (PC_PATCH *)pc->data;
2979:   PCPatchConstructType patchConstructionType = PC_PATCH_STAR;
2980:   char                 sub_mat_type[PETSC_MAX_PATH_LEN];
2981:   char                 option[PETSC_MAX_PATH_LEN];
2982:   const char          *prefix;
2983:   PetscBool            flg, dimflg, codimflg;
2984:   MPI_Comm             comm;
2985:   PetscInt            *ifields, nfields, k;
2986:   PCCompositeType      loctype = PC_COMPOSITE_ADDITIVE;

2988:   PetscFunctionBegin;
2989:   PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
2990:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
2991:   PetscOptionsHeadBegin(PetscOptionsObject, "Patch solver options");

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

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

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

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

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

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

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

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

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

3032:   /* If the user has set the number of subspaces, use that for the buffer size,
3033:    otherwise use a large number */
3034:   if (patch->nsubspaces <= 0) {
3035:     nfields = 128;
3036:   } else {
3037:     nfields = patch->nsubspaces;
3038:   }
3039:   PetscCall(PetscMalloc1(nfields, &ifields));
3040:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exclude_subspaces", patch->classname));
3041:   PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, option, ifields, &nfields, &flg));
3042:   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");
3043:   if (flg) {
3044:     PetscCall(PetscHSetIClear(patch->subspaces_to_exclude));
3045:     for (k = 0; k < nfields; k++) PetscCall(PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]));
3046:   }
3047:   PetscCall(PetscFree(ifields));

3049:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_patches_view", patch->classname));
3050:   PetscCall(PetscOptionsBool(option, "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg));
3051:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_cells_view", patch->classname));
3052:   PetscCall(PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerCells, &patch->formatCells, &patch->viewCells));
3053:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_interior_facets_view", patch->classname));
3054:   PetscCall(PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerIntFacets, &patch->formatIntFacets, &patch->viewIntFacets));
3055:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exterior_facets_view", patch->classname));
3056:   PetscCall(PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerExtFacets, &patch->formatExtFacets, &patch->viewExtFacets));
3057:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_points_view", patch->classname));
3058:   PetscCall(PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints));
3059:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_section_view", patch->classname));
3060:   PetscCall(PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerSection, &patch->formatSection, &patch->viewSection));
3061:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_mat_view", patch->classname));
3062:   PetscCall(PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix));
3063:   PetscOptionsHeadEnd();
3064:   patch->optionsSet = PETSC_TRUE;
3065:   PetscFunctionReturn(PETSC_SUCCESS);
3066: }

3068: static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc)
3069: {
3070:   PC_PATCH          *patch = (PC_PATCH *)pc->data;
3071:   KSPConvergedReason reason;
3072:   PetscInt           i;

3074:   PetscFunctionBegin;
3075:   if (!patch->save_operators) {
3076:     /* Can't do this here because the sub KSPs don't have an operator attached yet. */
3077:     PetscFunctionReturn(PETSC_SUCCESS);
3078:   }
3079:   if (patch->denseinverse) {
3080:     /* No solvers */
3081:     PetscFunctionReturn(PETSC_SUCCESS);
3082:   }
3083:   for (i = 0; i < patch->npatch; ++i) {
3084:     if (!((KSP)patch->solver[i])->setfromoptionscalled) PetscCall(KSPSetFromOptions((KSP)patch->solver[i]));
3085:     PetscCall(KSPSetUp((KSP)patch->solver[i]));
3086:     PetscCall(KSPGetConvergedReason((KSP)patch->solver[i], &reason));
3087:     if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
3088:   }
3089:   PetscFunctionReturn(PETSC_SUCCESS);
3090: }

3092: static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer)
3093: {
3094:   PC_PATCH   *patch = (PC_PATCH *)pc->data;
3095:   PetscViewer sviewer;
3096:   PetscBool   isascii;
3097:   PetscMPIInt rank;

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

3124:   if (patch->denseinverse) {
3125:     PetscCall(PetscViewerASCIIPrintf(viewer, "Explicitly forming dense inverse and applying patch solver via MatMult.\n"));
3126:   } else {
3127:     if (patch->isNonlinear) {
3128:       PetscCall(PetscViewerASCIIPrintf(viewer, "SNES on patches (all same):\n"));
3129:     } else {
3130:       PetscCall(PetscViewerASCIIPrintf(viewer, "KSP on patches (all same):\n"));
3131:     }
3132:     if (patch->solver) {
3133:       PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
3134:       if (rank == 0) {
3135:         PetscCall(PetscViewerASCIIPushTab(sviewer));
3136:         PetscCall(PetscObjectView(patch->solver[0], sviewer));
3137:         PetscCall(PetscViewerASCIIPopTab(sviewer));
3138:       }
3139:       PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
3140:     } else {
3141:       PetscCall(PetscViewerASCIIPushTab(viewer));
3142:       PetscCall(PetscViewerASCIIPrintf(viewer, "Solver not yet set.\n"));
3143:       PetscCall(PetscViewerASCIIPopTab(viewer));
3144:     }
3145:   }
3146:   PetscCall(PetscViewerASCIIPopTab(viewer));
3147:   PetscFunctionReturn(PETSC_SUCCESS);
3148: }

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

3155:    Options Database Keys:
3156: + -pc_patch_cells_view   - Views the process local cell numbers for each patch
3157: . -pc_patch_points_view  - Views the process local mesh point numbers for each patch
3158: . -pc_patch_g2l_view     - Views the map between global dofs and patch local dofs for each patch
3159: . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary
3160: - -pc_patch_sub_mat_view - Views the matrix associated with each patch

3162:    Level: intermediate

3164: .seealso: [](ch_ksp), `PCType`, `PCCreate()`, `PCSetType()`, `PCASM`, `PCJACOBI`, `PCPBJACOBI`, `PCVPBJACOBI`, `SNESPATCH`
3165: M*/
3166: PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc)
3167: {
3168:   PC_PATCH *patch;

3170:   PetscFunctionBegin;
3171:   PetscCall(PetscCitationsRegister(PCPatchCitation, &PCPatchcite));
3172:   PetscCall(PetscNew(&patch));

3174:   if (patch->subspaces_to_exclude) PetscCall(PetscHSetIDestroy(&patch->subspaces_to_exclude));
3175:   PetscCall(PetscHSetICreate(&patch->subspaces_to_exclude));

3177:   patch->classname   = "pc";
3178:   patch->isNonlinear = PETSC_FALSE;

3180:   /* Set some defaults */
3181:   patch->combined                 = PETSC_FALSE;
3182:   patch->save_operators           = PETSC_TRUE;
3183:   patch->local_composition_type   = PC_COMPOSITE_ADDITIVE;
3184:   patch->precomputeElementTensors = PETSC_FALSE;
3185:   patch->partition_of_unity       = PETSC_FALSE;
3186:   patch->codim                    = -1;
3187:   patch->dim                      = -1;
3188:   patch->vankadim                 = -1;
3189:   patch->ignoredim                = -1;
3190:   patch->pardecomp_overlap        = 0;
3191:   patch->patchconstructop         = PCPatchConstruct_Star;
3192:   patch->symmetrise_sweep         = PETSC_FALSE;
3193:   patch->npatch                   = 0;
3194:   patch->userIS                   = NULL;
3195:   patch->optionsSet               = PETSC_FALSE;
3196:   patch->iterationSet             = NULL;
3197:   patch->user_patches             = PETSC_FALSE;
3198:   PetscCall(PetscStrallocpy(MATDENSE, (char **)&patch->sub_mat_type));
3199:   patch->viewPatches                       = PETSC_FALSE;
3200:   patch->viewCells                         = PETSC_FALSE;
3201:   patch->viewPoints                        = PETSC_FALSE;
3202:   patch->viewSection                       = PETSC_FALSE;
3203:   patch->viewMatrix                        = PETSC_FALSE;
3204:   patch->densesolve                        = NULL;
3205:   patch->setupsolver                       = PCSetUp_PATCH_Linear;
3206:   patch->applysolver                       = PCApply_PATCH_Linear;
3207:   patch->resetsolver                       = PCReset_PATCH_Linear;
3208:   patch->destroysolver                     = PCDestroy_PATCH_Linear;
3209:   patch->updatemultiplicative              = PCUpdateMultiplicative_PATCH_Linear;
3210:   patch->dofMappingWithoutToWithArtificial = NULL;
3211:   patch->dofMappingWithoutToWithAll        = NULL;

3213:   pc->data                 = (void *)patch;
3214:   pc->ops->apply           = PCApply_PATCH;
3215:   pc->ops->applytranspose  = NULL; /* PCApplyTranspose_PATCH; */
3216:   pc->ops->setup           = PCSetUp_PATCH;
3217:   pc->ops->reset           = PCReset_PATCH;
3218:   pc->ops->destroy         = PCDestroy_PATCH;
3219:   pc->ops->setfromoptions  = PCSetFromOptions_PATCH;
3220:   pc->ops->setuponblocks   = PCSetUpOnBlocks_PATCH;
3221:   pc->ops->view            = PCView_PATCH;
3222:   pc->ops->applyrichardson = NULL;
3223:   PetscFunctionReturn(PETSC_SUCCESS);
3224: }