Actual source code: pcpatch.c

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

582:   Logically Collective

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

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

600:   Level: advanced

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

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

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

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

620:   Logically Collective

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

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

638:   Level: advanced

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

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

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

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

658:   Logically Collective

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

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

676:   Level: advanced

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

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

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

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

696:   Logically Collective

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

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

714:   Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1838:       PetscCall(PetscFree(zeroes));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2676:       PetscCall(VecReciprocal(global));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3175:    Level: intermediate

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

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

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

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

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

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