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: PetscLogEvent PC_Patch_CreatePatches, PC_Patch_ComputeOp, PC_Patch_Solve, PC_Patch_Apply, PC_Patch_Prealloc;

 12: static inline PetscErrorCode ObjectView(PetscObject obj, PetscViewer viewer, PetscViewerFormat format)
 13: {
 14:   PetscCall(PetscViewerPushFormat(viewer, format));
 15:   PetscCall(PetscObjectView(obj, viewer));
 16:   PetscCall(PetscViewerPopFormat(viewer));
 17:   return PETSC_SUCCESS;
 18: }

 20: static PetscErrorCode PCPatchConstruct_Star(void *vpatch, DM dm, PetscInt point, PetscHSetI ht)
 21: {
 22:   PetscInt  starSize;
 23:   PetscInt *star = NULL, si;

 25:   PetscFunctionBegin;
 26:   PetscCall(PetscHSetIClear(ht));
 27:   /* To start with, add the point we care about */
 28:   PetscCall(PetscHSetIAdd(ht, point));
 29:   /* Loop over all the points that this point connects to */
 30:   PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star));
 31:   for (si = 0; si < starSize * 2; si += 2) PetscCall(PetscHSetIAdd(ht, star[si]));
 32:   PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star));
 33:   PetscFunctionReturn(PETSC_SUCCESS);
 34: }

 36: static PetscErrorCode PCPatchConstruct_Vanka(void *vpatch, DM dm, PetscInt point, PetscHSetI ht)
 37: {
 38:   PC_PATCH *patch = (PC_PATCH *)vpatch;
 39:   PetscInt  starSize;
 40:   PetscInt *star         = NULL;
 41:   PetscBool shouldIgnore = PETSC_FALSE;
 42:   PetscInt  cStart, cEnd, iStart, iEnd, si;

 44:   PetscFunctionBegin;
 45:   PetscCall(PetscHSetIClear(ht));
 46:   /* To start with, add the point we care about */
 47:   PetscCall(PetscHSetIAdd(ht, point));
 48:   /* Should we ignore any points of a certain dimension? */
 49:   if (patch->vankadim >= 0) {
 50:     shouldIgnore = PETSC_TRUE;
 51:     PetscCall(DMPlexGetDepthStratum(dm, patch->vankadim, &iStart, &iEnd));
 52:   }
 53:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
 54:   /* Loop over all the cells that this point connects to */
 55:   PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star));
 56:   for (si = 0; si < starSize * 2; si += 2) {
 57:     const PetscInt cell = star[si];
 58:     PetscInt       closureSize;
 59:     PetscInt      *closure = NULL, ci;

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

 67:       /* We've been told to ignore entities of this type.*/
 68:       if (shouldIgnore && newpoint >= iStart && newpoint < iEnd) continue;
 69:       PetscCall(PetscHSetIAdd(ht, newpoint));
 70:     }
 71:     PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
 72:   }
 73:   PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star));
 74:   PetscFunctionReturn(PETSC_SUCCESS);
 75: }

 77: static PetscErrorCode PCPatchConstruct_Pardecomp(void *vpatch, DM dm, PetscInt point, PetscHSetI ht)
 78: {
 79:   PC_PATCH       *patch   = (PC_PATCH *)vpatch;
 80:   DMLabel         ghost   = NULL;
 81:   const PetscInt *leaves  = NULL;
 82:   PetscInt        nleaves = 0, pStart, pEnd, loc;
 83:   PetscBool       isFiredrake;
 84:   PetscBool       flg;
 85:   PetscInt        starSize;
 86:   PetscInt       *star = NULL;
 87:   PetscInt        opoint, overlapi;

 89:   PetscFunctionBegin;
 90:   PetscCall(PetscHSetIClear(ht));

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

 94:   PetscCall(DMHasLabel(dm, "pyop2_ghost", &isFiredrake));
 95:   if (isFiredrake) {
 96:     PetscCall(DMGetLabel(dm, "pyop2_ghost", &ghost));
 97:     PetscCall(DMLabelCreateIndex(ghost, pStart, pEnd));
 98:   } else {
 99:     PetscSF sf;
100:     PetscCall(DMGetPointSF(dm, &sf));
101:     PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL));
102:     nleaves = PetscMax(nleaves, 0);
103:   }

105:   for (opoint = pStart; opoint < pEnd; ++opoint) {
106:     if (ghost) PetscCall(DMLabelHasPoint(ghost, opoint, &flg));
107:     else {
108:       PetscCall(PetscFindInt(opoint, nleaves, leaves, &loc));
109:       flg = loc >= 0 ? PETSC_TRUE : PETSC_FALSE;
110:     }
111:     /* Not an owned entity, don't make a cell patch. */
112:     if (flg) continue;
113:     PetscCall(PetscHSetIAdd(ht, opoint));
114:   }

116:   /* Now build the overlap for the patch */
117:   for (overlapi = 0; overlapi < patch->pardecomp_overlap; ++overlapi) {
118:     PetscInt  index    = 0;
119:     PetscInt *htpoints = NULL;
120:     PetscInt  htsize;
121:     PetscInt  i;

123:     PetscCall(PetscHSetIGetSize(ht, &htsize));
124:     PetscCall(PetscMalloc1(htsize, &htpoints));
125:     PetscCall(PetscHSetIGetElems(ht, &index, htpoints));

127:     for (i = 0; i < htsize; ++i) {
128:       PetscInt hpoint = htpoints[i];
129:       PetscInt si;

131:       PetscCall(DMPlexGetTransitiveClosure(dm, hpoint, PETSC_FALSE, &starSize, &star));
132:       for (si = 0; si < starSize * 2; si += 2) {
133:         const PetscInt starp = star[si];
134:         PetscInt       closureSize;
135:         PetscInt      *closure = NULL, ci;

137:         /* now loop over all entities in the closure of starp */
138:         PetscCall(DMPlexGetTransitiveClosure(dm, starp, PETSC_TRUE, &closureSize, &closure));
139:         for (ci = 0; ci < closureSize * 2; ci += 2) {
140:           const PetscInt closstarp = closure[ci];
141:           PetscCall(PetscHSetIAdd(ht, closstarp));
142:         }
143:         PetscCall(DMPlexRestoreTransitiveClosure(dm, starp, PETSC_TRUE, &closureSize, &closure));
144:       }
145:       PetscCall(DMPlexRestoreTransitiveClosure(dm, hpoint, PETSC_FALSE, &starSize, &star));
146:     }
147:     PetscCall(PetscFree(htpoints));
148:   }

150:   PetscFunctionReturn(PETSC_SUCCESS);
151: }

153: /* The user's already set the patches in patch->userIS. Build the hash tables */
154: static PetscErrorCode PCPatchConstruct_User(void *vpatch, DM dm, PetscInt point, PetscHSetI ht)
155: {
156:   PC_PATCH       *patch   = (PC_PATCH *)vpatch;
157:   IS              patchis = patch->userIS[point];
158:   PetscInt        n;
159:   const PetscInt *patchdata;
160:   PetscInt        pStart, pEnd, i;

162:   PetscFunctionBegin;
163:   PetscCall(PetscHSetIClear(ht));
164:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
165:   PetscCall(ISGetLocalSize(patchis, &n));
166:   PetscCall(ISGetIndices(patchis, &patchdata));
167:   for (i = 0; i < n; ++i) {
168:     const PetscInt ownedpoint = patchdata[i];

170:     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);
171:     PetscCall(PetscHSetIAdd(ht, ownedpoint));
172:   }
173:   PetscCall(ISRestoreIndices(patchis, &patchdata));
174:   PetscFunctionReturn(PETSC_SUCCESS);
175: }

177: static PetscErrorCode PCPatchCreateDefaultSF_Private(PC pc, PetscInt n, const PetscSF *sf, const PetscInt *bs)
178: {
179:   PC_PATCH *patch = (PC_PATCH *)pc->data;
180:   PetscInt  i;

182:   PetscFunctionBegin;
183:   if (n == 1 && bs[0] == 1) {
184:     patch->sectionSF = sf[0];
185:     PetscCall(PetscObjectReference((PetscObject)patch->sectionSF));
186:   } else {
187:     PetscInt     allRoots = 0, allLeaves = 0;
188:     PetscInt     leafOffset    = 0;
189:     PetscInt    *ilocal        = NULL;
190:     PetscSFNode *iremote       = NULL;
191:     PetscInt    *remoteOffsets = NULL;
192:     PetscInt     index         = 0;
193:     PetscHMapI   rankToIndex;
194:     PetscInt     numRanks = 0;
195:     PetscSFNode *remote   = NULL;
196:     PetscSF      rankSF;
197:     PetscInt    *ranks   = NULL;
198:     PetscInt    *offsets = NULL;
199:     MPI_Datatype contig;
200:     PetscHSetI   ranksUniq;

202:     /* First figure out how many dofs there are in the concatenated numbering.
203:        allRoots: number of owned global dofs;
204:        allLeaves: number of visible dofs (global + ghosted).
205:     */
206:     for (i = 0; i < n; ++i) {
207:       PetscInt nroots, nleaves;

209:       PetscCall(PetscSFGetGraph(sf[i], &nroots, &nleaves, NULL, NULL));
210:       allRoots += nroots * bs[i];
211:       allLeaves += nleaves * bs[i];
212:     }
213:     PetscCall(PetscMalloc1(allLeaves, &ilocal));
214:     PetscCall(PetscMalloc1(allLeaves, &iremote));
215:     /* Now build an SF that just contains process connectivity. */
216:     PetscCall(PetscHSetICreate(&ranksUniq));
217:     for (i = 0; i < n; ++i) {
218:       const PetscMPIInt *ranks = NULL;
219:       PetscInt           nranks, j;

221:       PetscCall(PetscSFSetUp(sf[i]));
222:       PetscCall(PetscSFGetRootRanks(sf[i], &nranks, &ranks, NULL, NULL, NULL));
223:       /* These are all the ranks who communicate with me. */
224:       for (j = 0; j < nranks; ++j) PetscCall(PetscHSetIAdd(ranksUniq, (PetscInt)ranks[j]));
225:     }
226:     PetscCall(PetscHSetIGetSize(ranksUniq, &numRanks));
227:     PetscCall(PetscMalloc1(numRanks, &remote));
228:     PetscCall(PetscMalloc1(numRanks, &ranks));
229:     PetscCall(PetscHSetIGetElems(ranksUniq, &index, ranks));

231:     PetscCall(PetscHMapICreate(&rankToIndex));
232:     for (i = 0; i < numRanks; ++i) {
233:       remote[i].rank  = ranks[i];
234:       remote[i].index = 0;
235:       PetscCall(PetscHMapISet(rankToIndex, ranks[i], i));
236:     }
237:     PetscCall(PetscFree(ranks));
238:     PetscCall(PetscHSetIDestroy(&ranksUniq));
239:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)pc), &rankSF));
240:     PetscCall(PetscSFSetGraph(rankSF, 1, numRanks, NULL, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER));
241:     PetscCall(PetscSFSetUp(rankSF));
242:     /* OK, use it to communicate the root offset on the remote processes for each subspace. */
243:     PetscCall(PetscMalloc1(n, &offsets));
244:     PetscCall(PetscMalloc1(n * numRanks, &remoteOffsets));

246:     offsets[0] = 0;
247:     for (i = 1; i < n; ++i) {
248:       PetscInt nroots;

250:       PetscCall(PetscSFGetGraph(sf[i - 1], &nroots, NULL, NULL, NULL));
251:       offsets[i] = offsets[i - 1] + nroots * bs[i - 1];
252:     }
253:     /* Offsets are the offsets on the current process of the global dof numbering for the subspaces. */
254:     PetscCallMPI(MPI_Type_contiguous(n, MPIU_INT, &contig));
255:     PetscCallMPI(MPI_Type_commit(&contig));

257:     PetscCall(PetscSFBcastBegin(rankSF, contig, offsets, remoteOffsets, MPI_REPLACE));
258:     PetscCall(PetscSFBcastEnd(rankSF, contig, offsets, remoteOffsets, MPI_REPLACE));
259:     PetscCallMPI(MPI_Type_free(&contig));
260:     PetscCall(PetscFree(offsets));
261:     PetscCall(PetscSFDestroy(&rankSF));
262:     /* Now remoteOffsets contains the offsets on the remote
263:       processes who communicate with me.  So now we can
264:       concatenate the list of SFs into a single one. */
265:     index = 0;
266:     for (i = 0; i < n; ++i) {
267:       const PetscSFNode *remote = NULL;
268:       const PetscInt    *local  = NULL;
269:       PetscInt           nroots, nleaves, j;

271:       PetscCall(PetscSFGetGraph(sf[i], &nroots, &nleaves, &local, &remote));
272:       for (j = 0; j < nleaves; ++j) {
273:         PetscInt rank = remote[j].rank;
274:         PetscInt idx, rootOffset, k;

276:         PetscCall(PetscHMapIGet(rankToIndex, rank, &idx));
277:         PetscCheck(idx != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Didn't find rank, huh?");
278:         /* Offset on given rank for ith subspace */
279:         rootOffset = remoteOffsets[n * idx + i];
280:         for (k = 0; k < bs[i]; ++k) {
281:           ilocal[index]        = (local ? local[j] : j) * bs[i] + k + leafOffset;
282:           iremote[index].rank  = remote[j].rank;
283:           iremote[index].index = remote[j].index * bs[i] + k + rootOffset;
284:           ++index;
285:         }
286:       }
287:       leafOffset += nleaves * bs[i];
288:     }
289:     PetscCall(PetscHMapIDestroy(&rankToIndex));
290:     PetscCall(PetscFree(remoteOffsets));
291:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)pc), &patch->sectionSF));
292:     PetscCall(PetscSFSetGraph(patch->sectionSF, allRoots, allLeaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
293:   }
294:   PetscFunctionReturn(PETSC_SUCCESS);
295: }

297: PetscErrorCode PCPatchSetDenseInverse(PC pc, PetscBool flg)
298: {
299:   PC_PATCH *patch = (PC_PATCH *)pc->data;
300:   PetscFunctionBegin;
301:   patch->denseinverse = flg;
302:   PetscFunctionReturn(PETSC_SUCCESS);
303: }

305: PetscErrorCode PCPatchGetDenseInverse(PC pc, PetscBool *flg)
306: {
307:   PC_PATCH *patch = (PC_PATCH *)pc->data;
308:   PetscFunctionBegin;
309:   *flg = patch->denseinverse;
310:   PetscFunctionReturn(PETSC_SUCCESS);
311: }

313: /* TODO: Docs */
314: PetscErrorCode PCPatchSetIgnoreDim(PC pc, PetscInt dim)
315: {
316:   PC_PATCH *patch = (PC_PATCH *)pc->data;
317:   PetscFunctionBegin;
318:   patch->ignoredim = dim;
319:   PetscFunctionReturn(PETSC_SUCCESS);
320: }

322: /* TODO: Docs */
323: PetscErrorCode PCPatchGetIgnoreDim(PC pc, PetscInt *dim)
324: {
325:   PC_PATCH *patch = (PC_PATCH *)pc->data;
326:   PetscFunctionBegin;
327:   *dim = patch->ignoredim;
328:   PetscFunctionReturn(PETSC_SUCCESS);
329: }

331: /* TODO: Docs */
332: PetscErrorCode PCPatchSetSaveOperators(PC pc, PetscBool flg)
333: {
334:   PC_PATCH *patch = (PC_PATCH *)pc->data;
335:   PetscFunctionBegin;
336:   patch->save_operators = flg;
337:   PetscFunctionReturn(PETSC_SUCCESS);
338: }

340: /* TODO: Docs */
341: PetscErrorCode PCPatchGetSaveOperators(PC pc, PetscBool *flg)
342: {
343:   PC_PATCH *patch = (PC_PATCH *)pc->data;
344:   PetscFunctionBegin;
345:   *flg = patch->save_operators;
346:   PetscFunctionReturn(PETSC_SUCCESS);
347: }

349: /* TODO: Docs */
350: PetscErrorCode PCPatchSetPrecomputeElementTensors(PC pc, PetscBool flg)
351: {
352:   PC_PATCH *patch = (PC_PATCH *)pc->data;
353:   PetscFunctionBegin;
354:   patch->precomputeElementTensors = flg;
355:   PetscFunctionReturn(PETSC_SUCCESS);
356: }

358: /* TODO: Docs */
359: PetscErrorCode PCPatchGetPrecomputeElementTensors(PC pc, PetscBool *flg)
360: {
361:   PC_PATCH *patch = (PC_PATCH *)pc->data;
362:   PetscFunctionBegin;
363:   *flg = patch->precomputeElementTensors;
364:   PetscFunctionReturn(PETSC_SUCCESS);
365: }

367: /* TODO: Docs */
368: PetscErrorCode PCPatchSetPartitionOfUnity(PC pc, PetscBool flg)
369: {
370:   PC_PATCH *patch = (PC_PATCH *)pc->data;
371:   PetscFunctionBegin;
372:   patch->partition_of_unity = flg;
373:   PetscFunctionReturn(PETSC_SUCCESS);
374: }

376: /* TODO: Docs */
377: PetscErrorCode PCPatchGetPartitionOfUnity(PC pc, PetscBool *flg)
378: {
379:   PC_PATCH *patch = (PC_PATCH *)pc->data;
380:   PetscFunctionBegin;
381:   *flg = patch->partition_of_unity;
382:   PetscFunctionReturn(PETSC_SUCCESS);
383: }

385: /* TODO: Docs */
386: PetscErrorCode PCPatchSetLocalComposition(PC pc, PCCompositeType type)
387: {
388:   PC_PATCH *patch = (PC_PATCH *)pc->data;
389:   PetscFunctionBegin;
390:   PetscCheck(type == PC_COMPOSITE_ADDITIVE || type == PC_COMPOSITE_MULTIPLICATIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Only supports additive or multiplicative as the local type");
391:   patch->local_composition_type = type;
392:   PetscFunctionReturn(PETSC_SUCCESS);
393: }

395: /* TODO: Docs */
396: PetscErrorCode PCPatchGetLocalComposition(PC pc, PCCompositeType *type)
397: {
398:   PC_PATCH *patch = (PC_PATCH *)pc->data;
399:   PetscFunctionBegin;
400:   *type = patch->local_composition_type;
401:   PetscFunctionReturn(PETSC_SUCCESS);
402: }

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

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

415: /* TODO: Docs */
416: PetscErrorCode PCPatchGetSubMatType(PC pc, MatType *sub_mat_type)
417: {
418:   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;
439:   PetscFunctionBegin;
440:   *cellNumbering = patch->cellNumbering;
441:   PetscFunctionReturn(PETSC_SUCCESS);
442: }

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

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

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

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

502: /* TODO: Docs */
503: 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)
504: {
505:   PC_PATCH *patch = (PC_PATCH *)pc->data;
506:   DM        dm, plex;
507:   PetscSF  *sfs;
508:   PetscInt  cStart, cEnd, i, j;

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

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

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

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

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

578: /*@C

580:   PCPatchSetComputeFunction - Set the callback used to compute patch residuals

582:   Logically Collective

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

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

601:   Level: advanced

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

606: .seealso: `PCPatchSetComputeOperator()`, `PCPatchGetComputeOperator()`, `PCPatchSetDiscretisationInfo()`, `PCPatchSetComputeFunctionInteriorFacets()`
607: @*/
608: PetscErrorCode PCPatchSetComputeFunction(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Vec, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx)
609: {
610:   PC_PATCH *patch = (PC_PATCH *)pc->data;

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

618: /*@C

620:   PCPatchSetComputeFunctionInteriorFacets - Set the callback used to compute facet integrals for patch residuals

622:   Logically Collective

624:   Input Parameters:
625: + pc   - The `PC`
626: . func - The callback
627: - ctx  - The user context

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

641:   Level: advanced

643:   Note:
644:   The entries of F (the output residual vector) have been set to zero before the call.

646: .seealso: `PCPatchSetComputeOperator()`, `PCPatchGetComputeOperator()`, `PCPatchSetDiscretisationInfo()`, `PCPatchSetComputeFunction()`
647: @*/
648: PetscErrorCode PCPatchSetComputeFunctionInteriorFacets(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Vec, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx)
649: {
650:   PC_PATCH *patch = (PC_PATCH *)pc->data;

652:   PetscFunctionBegin;
653:   patch->usercomputefintfacet    = func;
654:   patch->usercomputefintfacetctx = ctx;
655:   PetscFunctionReturn(PETSC_SUCCESS);
656: }

658: /*@C

660:   PCPatchSetComputeOperator - Set the callback used to compute patch matrices

662:   Logically Collective

664:   Input Parameters:
665: + pc   - The `PC`
666: . func - The callback
667: - ctx  - The user context

669:   Calling sequence of `func`:
670: $ PetscErrorCode func(PC pc, PetscInt point, Vec x, Mat mat, IS facetIS, PetscInt n, const PetscInt* dofsArray, const PetscInt* dofsArrayWithAll, void* ctx)
671: +  pc               - The `PC`
672: .  point            - The point
673: .  x                - The input solution (not used in linear problems)
674: .  mat              - The patch matrix
675: .  cellIS           - An array of the cell numbers
676: .  n                - The size of `dofsArray`
677: .  dofsArray        - The dofmap for the dofs to be solved for
678: .  dofsArrayWithAll - The dofmap for all dofs on the patch
679: -  ctx              - The user context

681:   Level: advanced

683:   Note:
684:   The matrix entries have been set to zero before the call.

686: .seealso: `PCPatchGetComputeOperator()`, `PCPatchSetComputeFunction()`, `PCPatchSetDiscretisationInfo()`
687: @*/
688: PetscErrorCode PCPatchSetComputeOperator(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Mat, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx)
689: {
690:   PC_PATCH *patch = (PC_PATCH *)pc->data;

692:   PetscFunctionBegin;
693:   patch->usercomputeop    = func;
694:   patch->usercomputeopctx = ctx;
695:   PetscFunctionReturn(PETSC_SUCCESS);
696: }

698: /*@C

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

702:   Logically Collective

704:   Input Parameters:
705: + pc   - The `PC`
706: . func - The callback
707: - ctx  - The user context

709:   Calling sequence of `func`:
710: $  PetscErrorCode func(PC pc, PetscInt point, Vec x, Mat mat, IS facetIS, PetscInt n, const PetscInt* dofsArray, const PetscInt* dofsArrayWithAll, void* ctx)
711: +  pc               - The `PC`
712: .  point            - The point
713: .  x                - The input solution (not used in linear problems)
714: .  mat              - The patch matrix
715: .  facetIS          - An array of the facet numbers
716: .  n                - The size of `dofsArray`
717: .  dofsArray        - The dofmap for the dofs to be solved for
718: .  dofsArrayWithAll - The dofmap for all dofs on the patch
719: -  ctx              - The user context

721:   Level: advanced

723:   Note:
724:   The matrix entries have been set to zero before the call.

726: .seealso: `PCPatchGetComputeOperator()`, `PCPatchSetComputeFunction()`, `PCPatchSetDiscretisationInfo()`
727: @*/
728: PetscErrorCode PCPatchSetComputeOperatorInteriorFacets(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Mat, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx)
729: {
730:   PC_PATCH *patch = (PC_PATCH *)pc->data;

732:   PetscFunctionBegin;
733:   patch->usercomputeopintfacet    = func;
734:   patch->usercomputeopintfacetctx = ctx;
735:   PetscFunctionReturn(PETSC_SUCCESS);
736: }

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

754:   PetscFunctionBegin;
755:   PetscCall(PCGetDM(pc, &dm));
756:   PetscCall(DMConvert(dm, DMPLEX, &plex));
757:   dm = plex;
758:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fBegin, &fEnd));
759:   PetscCall(PCPatchGetIgnoreDim(pc, &ignoredim));
760:   if (ignoredim >= 0) PetscCall(DMPlexGetDepthStratum(dm, ignoredim, &iStart, &iEnd));
761:   PetscCall(PetscHSetIClear(cht));
762:   PetscHashIterBegin(ht, hi);
763:   while (!PetscHashIterAtEnd(ht, hi)) {
764:     PetscHashIterGetKey(ht, hi, point);
765:     PetscHashIterNext(ht, hi);

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

801: static PetscErrorCode PCPatchGetGlobalDofs(PC pc, PetscSection dofSection[], PetscInt f, PetscBool combined, PetscInt p, PetscInt *dof, PetscInt *off)
802: {
803:   PetscFunctionBegin;
804:   if (combined) {
805:     if (f < 0) {
806:       if (dof) PetscCall(PetscSectionGetDof(dofSection[0], p, dof));
807:       if (off) PetscCall(PetscSectionGetOffset(dofSection[0], p, off));
808:     } else {
809:       if (dof) PetscCall(PetscSectionGetFieldDof(dofSection[0], p, f, dof));
810:       if (off) PetscCall(PetscSectionGetFieldOffset(dofSection[0], p, f, off));
811:     }
812:   } else {
813:     if (f < 0) {
814:       PC_PATCH *patch = (PC_PATCH *)pc->data;
815:       PetscInt  fdof, g;

817:       if (dof) {
818:         *dof = 0;
819:         for (g = 0; g < patch->nsubspaces; ++g) {
820:           PetscCall(PetscSectionGetDof(dofSection[g], p, &fdof));
821:           *dof += fdof;
822:         }
823:       }
824:       if (off) {
825:         *off = 0;
826:         for (g = 0; g < patch->nsubspaces; ++g) {
827:           PetscCall(PetscSectionGetOffset(dofSection[g], p, &fdof));
828:           *off += fdof;
829:         }
830:       }
831:     } else {
832:       if (dof) PetscCall(PetscSectionGetDof(dofSection[f], p, dof));
833:       if (off) PetscCall(PetscSectionGetOffset(dofSection[f], p, off));
834:     }
835:   }
836:   PetscFunctionReturn(PETSC_SUCCESS);
837: }

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

851:   PetscFunctionBegin;
852:   PetscCall(PetscHSetIClear(dofs));
853:   for (k = 0; k < patch->nsubspaces; ++k) {
854:     PetscInt subspaceOffset = patch->subspaceOffsets[k];
855:     PetscInt bs             = patch->bs[k];
856:     PetscInt j, l;

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

875:     PetscHashIterBegin(pts, hi);
876:     while (!PetscHashIterAtEnd(pts, hi)) {
877:       PetscHashIterGetKey(pts, hi, p);
878:       PetscHashIterNext(pts, hi);
879:       PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, p, &ldof, &loff));
880:       if (0 == ldof) continue;
881:       for (j = loff; j < ldof + loff; ++j) {
882:         for (l = 0; l < bs; ++l) {
883:           PetscInt dof = bs * j + l + subspaceOffset;
884:           PetscCall(PetscHSetIAdd(dofs, dof));
885:         }
886:       }
887:     }
888:   }
889:   PetscFunctionReturn(PETSC_SUCCESS);
890: }

892: /* Given two hash tables A and B, compute the keys in B that are not in A, and put them in C */
893: static PetscErrorCode PCPatchComputeSetDifference_Private(PetscHSetI A, PetscHSetI B, PetscHSetI C)
894: {
895:   PetscHashIter hi;
896:   PetscInt      key;
897:   PetscBool     flg;

899:   PetscFunctionBegin;
900:   PetscCall(PetscHSetIClear(C));
901:   PetscHashIterBegin(B, hi);
902:   while (!PetscHashIterAtEnd(B, hi)) {
903:     PetscHashIterGetKey(B, hi, key);
904:     PetscHashIterNext(B, hi);
905:     PetscCall(PetscHSetIHas(A, key, &flg));
906:     if (!flg) PetscCall(PetscHSetIAdd(C, key));
907:   }
908:   PetscFunctionReturn(PETSC_SUCCESS);
909: }

911: /*
912:   PCPatchCreateCellPatches - create patches.

914:   Input Parameter:
915:   . dm - The DMPlex object defining the mesh

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

936:   PetscFunctionBegin;
937:   /* Used to keep track of the cells in the patch. */
938:   PetscCall(PetscHSetICreate(&ht));
939:   PetscCall(PetscHSetICreate(&cht));

941:   PetscCall(PCGetDM(pc, &dm));
942:   PetscCheck(dm, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "DM not yet set on patch PC");
943:   PetscCall(DMConvert(dm, DMPLEX, &plex));
944:   dm = plex;
945:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
946:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));

948:   if (patch->user_patches) {
949:     PetscCall(patch->userpatchconstructionop(pc, &patch->npatch, &patch->userIS, &patch->iterationSet, patch->userpatchconstructctx));
950:     vStart = 0;
951:     vEnd   = patch->npatch;
952:   } else if (patch->ctype == PC_PATCH_PARDECOMP) {
953:     vStart = 0;
954:     vEnd   = 1;
955:   } else if (patch->codim < 0) {
956:     if (patch->dim < 0) PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
957:     else PetscCall(DMPlexGetDepthStratum(dm, patch->dim, &vStart, &vEnd));
958:   } else PetscCall(DMPlexGetHeightStratum(dm, patch->codim, &vStart, &vEnd));
959:   patch->npatch = vEnd - vStart;

961:   /* These labels mark the owned points.  We only create patches around points that this process owns. */
962:   PetscCall(DMHasLabel(dm, "pyop2_ghost", &isFiredrake));
963:   if (isFiredrake) {
964:     PetscCall(DMGetLabel(dm, "pyop2_ghost", &ghost));
965:     PetscCall(DMLabelCreateIndex(ghost, pStart, pEnd));
966:   } else {
967:     PetscSF sf;

969:     PetscCall(DMGetPointSF(dm, &sf));
970:     PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL));
971:     nleaves = PetscMax(nleaves, 0);
972:   }

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

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

1007:     PetscCall(patch->patchconstructop((void *)patch, dm, v, ht));
1008:     PetscCall(PCPatchCompleteCellPatch(pc, ht, cht));
1009:     PetscCall(PetscHSetIGetSize(cht, &chtSize));
1010:     /* empty patch, continue */
1011:     if (chtSize == 0) continue;

1013:     /* safe because size(cht) > 0 from above */
1014:     PetscHashIterBegin(cht, hi);
1015:     while (!PetscHashIterAtEnd(cht, hi)) {
1016:       PetscInt point, pdof;

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

1052:   PetscCall(PetscSectionSetUp(cellCounts));
1053:   PetscCall(PetscSectionGetStorageSize(cellCounts, &numCells));
1054:   PetscCall(PetscMalloc1(numCells, &cellsArray));
1055:   PetscCall(PetscSectionSetUp(pointCounts));
1056:   PetscCall(PetscSectionGetStorageSize(pointCounts, &numPoints));
1057:   PetscCall(PetscMalloc1(numPoints, &pointsArray));

1059:   PetscCall(PetscSectionSetUp(intFacetCounts));
1060:   PetscCall(PetscSectionSetUp(extFacetCounts));
1061:   PetscCall(PetscSectionGetStorageSize(intFacetCounts, &numIntFacets));
1062:   PetscCall(PetscSectionGetStorageSize(extFacetCounts, &numExtFacets));
1063:   PetscCall(PetscMalloc1(numIntFacets, &intFacetsArray));
1064:   PetscCall(PetscMalloc1(numIntFacets * 2, &intFacetsToPatchCell));
1065:   PetscCall(PetscMalloc1(numExtFacets, &extFacetsArray));

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

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

1087:       PetscHashIterGetKey(cht, hi, point);
1088:       if (fStart <= point && point < fEnd) {
1089:         const PetscInt *support;
1090:         PetscInt        supportSize, p;
1091:         PetscBool       interior = PETSC_TRUE;
1092:         PetscCall(DMPlexGetSupport(dm, point, &support));
1093:         PetscCall(DMPlexGetSupportSize(dm, point, &supportSize));
1094:         if (supportSize == 1) {
1095:           interior = PETSC_FALSE;
1096:         } else {
1097:           for (p = 0; p < supportSize; p++) {
1098:             PetscBool found;
1099:             /* FIXME: can I do this while iterating over cht? */
1100:             PetscCall(PetscHSetIHas(cht, support[p], &found));
1101:             if (!found) {
1102:               interior = PETSC_FALSE;
1103:               break;
1104:             }
1105:           }
1106:         }
1107:         if (interior) {
1108:           intFacetsToPatchCell[2 * (ifoff + ifn)]     = support[0];
1109:           intFacetsToPatchCell[2 * (ifoff + ifn) + 1] = support[1];
1110:           intFacetsArray[ifoff + ifn++]               = point;
1111:         } else {
1112:           extFacetsArray[efoff + efn++] = point;
1113:         }
1114:       }
1115:       PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL));
1116:       if (pdof) pointsArray[off + n++] = point;
1117:       if (point >= cStart && point < cEnd) cellsArray[coff + cn++] = point;
1118:       PetscHashIterNext(cht, hi);
1119:     }
1120:     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);
1121:     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);
1122:     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);
1123:     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);

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

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

1177: /*
1178:   PCPatchCreateCellPatchDiscretisationInfo - Build the dof maps for cell patches

1180:   Input Parameters:
1181:   + dm - The DMPlex object defining the mesh
1182:   . cellCounts - Section with counts of cells around each vertex
1183:   . cells - IS of the cell point indices of cells in each patch
1184:   . cellNumbering - Section mapping plex cell points to Firedrake cell indices.
1185:   . nodesPerCell - number of nodes per cell.
1186:   - cellNodeMap - map from cells to node indices (nodesPerCell * numCells)

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

1237:   PetscFunctionBegin;

1239:   PetscCall(PCGetDM(pc, &dm));
1240:   PetscCall(DMConvert(dm, DMPLEX, &plex));
1241:   dm = plex;
1242:   /* dofcounts section is cellcounts section * dofPerCell */
1243:   PetscCall(PetscSectionGetStorageSize(cellCounts, &numCells));
1244:   PetscCall(PetscSectionGetStorageSize(patch->pointCounts, &numPoints));
1245:   numDofs = numCells * totalDofsPerCell;
1246:   PetscCall(PetscMalloc1(numDofs, &dofsArray));
1247:   PetscCall(PetscMalloc1(numPoints * Nf, &offsArray));
1248:   PetscCall(PetscMalloc1(numDofs, &asmArray));
1249:   PetscCall(PetscMalloc1(numCells, &newCellsArray));
1250:   PetscCall(PetscSectionGetChart(cellCounts, &vStart, &vEnd));
1251:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCounts));
1252:   gtolCounts = patch->gtolCounts;
1253:   PetscCall(PetscSectionSetChart(gtolCounts, vStart, vEnd));
1254:   PetscCall(PetscObjectSetName((PetscObject)patch->gtolCounts, "Patch Global Index Section"));

1256:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1257:     PetscCall(PetscMalloc1(numPoints * Nf, &offsArrayWithArtificial));
1258:     PetscCall(PetscMalloc1(numDofs, &asmArrayWithArtificial));
1259:     PetscCall(PetscMalloc1(numDofs, &dofsArrayWithArtificial));
1260:     PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithArtificial));
1261:     gtolCountsWithArtificial = patch->gtolCountsWithArtificial;
1262:     PetscCall(PetscSectionSetChart(gtolCountsWithArtificial, vStart, vEnd));
1263:     PetscCall(PetscObjectSetName((PetscObject)patch->gtolCountsWithArtificial, "Patch Global Index Section Including Artificial BCs"));
1264:   }

1266:   isNonlinear = patch->isNonlinear;
1267:   if (isNonlinear) {
1268:     PetscCall(PetscMalloc1(numPoints * Nf, &offsArrayWithAll));
1269:     PetscCall(PetscMalloc1(numDofs, &asmArrayWithAll));
1270:     PetscCall(PetscMalloc1(numDofs, &dofsArrayWithAll));
1271:     PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithAll));
1272:     gtolCountsWithAll = patch->gtolCountsWithAll;
1273:     PetscCall(PetscSectionSetChart(gtolCountsWithAll, vStart, vEnd));
1274:     PetscCall(PetscObjectSetName((PetscObject)patch->gtolCountsWithAll, "Patch Global Index Section Including All BCs"));
1275:   }

1277:   /* Outside the patch loop, get the dofs that are globally-enforced Dirichlet
1278:    conditions */
1279:   PetscCall(PetscHSetICreate(&globalBcs));
1280:   PetscCall(ISGetIndices(patch->ghostBcNodes, &bcNodes));
1281:   PetscCall(ISGetSize(patch->ghostBcNodes, &numBcs));
1282:   for (i = 0; i < numBcs; ++i) { PetscCall(PetscHSetIAdd(globalBcs, bcNodes[i])); /* these are already in concatenated numbering */ }
1283:   PetscCall(ISRestoreIndices(patch->ghostBcNodes, &bcNodes));
1284:   PetscCall(ISDestroy(&patch->ghostBcNodes)); /* memory optimisation */

1286:   /* Hash tables for artificial BC construction */
1287:   PetscCall(PetscHSetICreate(&ownedpts));
1288:   PetscCall(PetscHSetICreate(&seenpts));
1289:   PetscCall(PetscHSetICreate(&owneddofs));
1290:   PetscCall(PetscHSetICreate(&seendofs));
1291:   PetscCall(PetscHSetICreate(&artificialbcs));

1293:   PetscCall(ISGetIndices(cells, &cellsArray));
1294:   PetscCall(ISGetIndices(points, &pointsArray));
1295:   PetscCall(PetscHMapICreate(&ht));
1296:   PetscCall(PetscHMapICreate(&htWithArtificial));
1297:   PetscCall(PetscHMapICreate(&htWithAll));
1298:   for (v = vStart; v < vEnd; ++v) {
1299:     PetscInt localIndex               = 0;
1300:     PetscInt localIndexWithArtificial = 0;
1301:     PetscInt localIndexWithAll        = 0;
1302:     PetscInt dof, off, i, j, k, l;

1304:     PetscCall(PetscHMapIClear(ht));
1305:     PetscCall(PetscHMapIClear(htWithArtificial));
1306:     PetscCall(PetscHMapIClear(htWithAll));
1307:     PetscCall(PetscSectionGetDof(cellCounts, v, &dof));
1308:     PetscCall(PetscSectionGetOffset(cellCounts, v, &off));
1309:     if (dof <= 0) continue;

1311:     /* Calculate the global numbers of the artificial BC dofs here first */
1312:     PetscCall(patch->patchconstructop((void *)patch, dm, v, ownedpts));
1313:     PetscCall(PCPatchCompleteCellPatch(pc, ownedpts, seenpts));
1314:     PetscCall(PCPatchGetPointDofs(pc, ownedpts, owneddofs, v, &patch->subspaces_to_exclude));
1315:     PetscCall(PCPatchGetPointDofs(pc, seenpts, seendofs, v, NULL));
1316:     PetscCall(PCPatchComputeSetDifference_Private(owneddofs, seendofs, artificialbcs));
1317:     if (patch->viewPatches) {
1318:       PetscHSetI    globalbcdofs;
1319:       PetscHashIter hi;
1320:       MPI_Comm      comm = PetscObjectComm((PetscObject)pc);

1322:       PetscCall(PetscHSetICreate(&globalbcdofs));
1323:       PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": owned dofs:\n", v));
1324:       PetscHashIterBegin(owneddofs, hi);
1325:       while (!PetscHashIterAtEnd(owneddofs, hi)) {
1326:         PetscInt globalDof;

1328:         PetscHashIterGetKey(owneddofs, hi, globalDof);
1329:         PetscHashIterNext(owneddofs, hi);
1330:         PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));
1331:       }
1332:       PetscCall(PetscSynchronizedPrintf(comm, "\n"));
1333:       PetscCall(PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": seen dofs:\n", v));
1334:       PetscHashIterBegin(seendofs, hi);
1335:       while (!PetscHashIterAtEnd(seendofs, hi)) {
1336:         PetscInt  globalDof;
1337:         PetscBool flg;

1339:         PetscHashIterGetKey(seendofs, hi, globalDof);
1340:         PetscHashIterNext(seendofs, hi);
1341:         PetscCall(PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof));

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

1379:       for (i = off; i < off + dof; ++i) {
1380:         /* Walk over the cells in this patch. */
1381:         const PetscInt c    = cellsArray[i];
1382:         PetscInt       cell = c;

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

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

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

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

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

1462:   PetscCall(DMDestroy(&dm));
1463:   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);
1464:   PetscCall(PetscSectionSetUp(gtolCounts));
1465:   PetscCall(PetscSectionGetStorageSize(gtolCounts, &numGlobalDofs));
1466:   PetscCall(PetscMalloc1(numGlobalDofs, &globalDofsArray));

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

1483:     PetscCall(PetscHMapIClear(ht));
1484:     PetscCall(PetscHMapIClear(htWithArtificial));
1485:     PetscCall(PetscHMapIClear(htWithAll));
1486:     PetscCall(PetscSectionGetDof(cellCounts, v, &dof));
1487:     PetscCall(PetscSectionGetOffset(cellCounts, v, &off));
1488:     PetscCall(PetscSectionGetDof(pointCounts, v, &Np));
1489:     PetscCall(PetscSectionGetOffset(pointCounts, v, &ooff));
1490:     if (dof <= 0) continue;

1492:     for (k = 0; k < patch->nsubspaces; ++k) {
1493:       const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
1494:       PetscInt        nodesPerCell   = patch->nodesPerCell[k];
1495:       PetscInt        subspaceOffset = patch->subspaceOffsets[k];
1496:       PetscInt        bs             = patch->bs[k];
1497:       PetscInt        goff;

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

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

1523:       /* Shove it in the output data structure. */
1524:       PetscCall(PetscSectionGetOffset(gtolCounts, v, &goff));
1525:       PetscHashIterBegin(ht, hi);
1526:       while (!PetscHashIterAtEnd(ht, hi)) {
1527:         PetscInt globalDof, localDof;

1529:         PetscHashIterGetKey(ht, hi, globalDof);
1530:         PetscHashIterGetVal(ht, hi, localDof);
1531:         if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof;
1532:         PetscHashIterNext(ht, hi);
1533:       }

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

1558:       for (p = 0; p < Np; ++p) {
1559:         const PetscInt point = pointsArray[ooff + p];
1560:         PetscInt       globalDof, localDof;

1562:         PetscCall(PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof));
1563:         PetscCall(PetscHMapIGet(ht, globalDof, &localDof));
1564:         offsArray[(ooff + p) * Nf + k] = localDof;
1565:         if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1566:           PetscCall(PetscHMapIGet(htWithArtificial, globalDof, &localDof));
1567:           offsArrayWithArtificial[(ooff + p) * Nf + k] = localDof;
1568:         }
1569:         if (isNonlinear) {
1570:           PetscCall(PetscHMapIGet(htWithAll, globalDof, &localDof));
1571:           offsArrayWithAll[(ooff + p) * Nf + k] = localDof;
1572:         }
1573:       }
1574:     }

1576:     PetscCall(PetscHSetIDestroy(&globalBcs));
1577:     PetscCall(PetscHSetIDestroy(&ownedpts));
1578:     PetscCall(PetscHSetIDestroy(&seenpts));
1579:     PetscCall(PetscHSetIDestroy(&owneddofs));
1580:     PetscCall(PetscHSetIDestroy(&seendofs));
1581:     PetscCall(PetscHSetIDestroy(&artificialbcs));

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

1593:         if (cellNumbering) PetscCall(PetscSectionGetOffset(cellNumbering, c, &cell));
1594:         for (k = 0; k < patch->nsubspaces; ++k) {
1595:           const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
1596:           PetscInt        nodesPerCell   = patch->nodesPerCell[k];
1597:           PetscInt        subspaceOffset = patch->subspaceOffsets[k];
1598:           PetscInt        bs             = patch->bs[k];

1600:           for (j = 0; j < nodesPerCell; ++j) {
1601:             for (l = 0; l < bs; ++l) {
1602:               const PetscInt globalDof = cellNodeMap[cell * nodesPerCell + j] * bs + l + subspaceOffset;
1603:               PetscInt       localDof;

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

1632:   PetscCall(PetscHMapIDestroy(&ht));
1633:   PetscCall(PetscHMapIDestroy(&htWithArtificial));
1634:   PetscCall(PetscHMapIDestroy(&htWithAll));
1635:   PetscCall(ISRestoreIndices(cells, &cellsArray));
1636:   PetscCall(ISRestoreIndices(points, &pointsArray));
1637:   PetscCall(PetscFree(dofsArray));
1638:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscCall(PetscFree(dofsArrayWithArtificial));
1639:   if (isNonlinear) PetscCall(PetscFree(dofsArrayWithAll));
1640:   /* Create placeholder section for map from points to patch dofs */
1641:   PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection));
1642:   PetscCall(PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces));
1643:   if (patch->combined) {
1644:     PetscInt numFields;
1645:     PetscCall(PetscSectionGetNumFields(patch->dofSection[0], &numFields));
1646:     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);
1647:     PetscCall(PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd));
1648:     PetscCall(PetscSectionSetChart(patch->patchSection, pStart, pEnd));
1649:     for (p = pStart; p < pEnd; ++p) {
1650:       PetscInt dof, fdof, f;

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

1703: static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat, PetscBool withArtificial)
1704: {
1705:   PC_PATCH   *patch = (PC_PATCH *)pc->data;
1706:   PetscBool   flg;
1707:   PetscInt    csize, rsize;
1708:   const char *prefix = NULL;

1710:   PetscFunctionBegin;
1711:   if (withArtificial) {
1712:     /* would be nice if we could create a rectangular matrix of size numDofsWithArtificial x numDofs here */
1713:     PetscInt pStart;
1714:     PetscCall(PetscSectionGetChart(patch->gtolCountsWithArtificial, &pStart, NULL));
1715:     PetscCall(PetscSectionGetDof(patch->gtolCountsWithArtificial, point + pStart, &rsize));
1716:     csize = rsize;
1717:   } else {
1718:     PetscInt pStart;
1719:     PetscCall(PetscSectionGetChart(patch->gtolCounts, &pStart, NULL));
1720:     PetscCall(PetscSectionGetDof(patch->gtolCounts, point + pStart, &rsize));
1721:     csize = rsize;
1722:   }

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

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

1758:       PetscCall(PetscCalloc1(rsize, &dnnz));
1759:       PetscCall(PetscBTCreate(rsize * rsize, &bt));
1760:       for (c = 0; c < ncell; ++c) {
1761:         const PetscInt *idx = dofsArray + (offset + c) * patch->totalDofsPerCell;
1762:         for (i = 0; i < patch->totalDofsPerCell; ++i) {
1763:           const PetscInt row = idx[i];
1764:           if (row < 0) continue;
1765:           for (j = 0; j < patch->totalDofsPerCell; ++j) {
1766:             const PetscInt col = idx[j];
1767:             const PetscInt key = row * rsize + col;
1768:             if (col < 0) continue;
1769:             if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1770:           }
1771:         }
1772:       }

1774:       if (patch->usercomputeopintfacet) {
1775:         const PetscInt *intFacetsArray = NULL;
1776:         PetscInt        i, numIntFacets, intFacetOffset;
1777:         const PetscInt *facetCells = NULL;

1779:         PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
1780:         PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
1781:         PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
1782:         PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
1783:         for (i = 0; i < numIntFacets; i++) {
1784:           const PetscInt cell0 = facetCells[2 * (intFacetOffset + i) + 0];
1785:           const PetscInt cell1 = facetCells[2 * (intFacetOffset + i) + 1];
1786:           PetscInt       celli, cellj;

1788:           for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1789:             const PetscInt row = dofsArray[(offset + cell0) * patch->totalDofsPerCell + celli];
1790:             if (row < 0) continue;
1791:             for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1792:               const PetscInt col = dofsArray[(offset + cell1) * patch->totalDofsPerCell + cellj];
1793:               const PetscInt key = row * rsize + col;
1794:               if (col < 0) continue;
1795:               if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1796:             }
1797:           }

1799:           for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1800:             const PetscInt row = dofsArray[(offset + cell1) * patch->totalDofsPerCell + celli];
1801:             if (row < 0) continue;
1802:             for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1803:               const PetscInt col = dofsArray[(offset + cell0) * patch->totalDofsPerCell + cellj];
1804:               const PetscInt key = row * rsize + col;
1805:               if (col < 0) continue;
1806:               if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1807:             }
1808:           }
1809:         }
1810:       }
1811:       PetscCall(PetscBTDestroy(&bt));
1812:       PetscCall(MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL));
1813:       PetscCall(PetscFree(dnnz));

1815:       PetscCall(PetscCalloc1(patch->totalDofsPerCell * patch->totalDofsPerCell, &zeroes));
1816:       for (c = 0; c < ncell; ++c) {
1817:         const PetscInt *idx = &dofsArray[(offset + c) * patch->totalDofsPerCell];
1818:         PetscCall(MatSetValues(*mat, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, zeroes, INSERT_VALUES));
1819:       }
1820:       PetscCall(MatGetLocalSize(*mat, &rows, NULL));
1821:       for (i = 0; i < rows; ++i) PetscCall(MatSetValues(*mat, 1, &i, 1, &i, zeroes, INSERT_VALUES));

1823:       if (patch->usercomputeopintfacet) {
1824:         const PetscInt *intFacetsArray = NULL;
1825:         PetscInt        i, numIntFacets, intFacetOffset;
1826:         const PetscInt *facetCells = NULL;

1828:         PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
1829:         PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
1830:         PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
1831:         PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
1832:         for (i = 0; i < numIntFacets; i++) {
1833:           const PetscInt  cell0    = facetCells[2 * (intFacetOffset + i) + 0];
1834:           const PetscInt  cell1    = facetCells[2 * (intFacetOffset + i) + 1];
1835:           const PetscInt *cell0idx = &dofsArray[(offset + cell0) * patch->totalDofsPerCell];
1836:           const PetscInt *cell1idx = &dofsArray[(offset + cell1) * patch->totalDofsPerCell];
1837:           PetscCall(MatSetValues(*mat, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, zeroes, INSERT_VALUES));
1838:           PetscCall(MatSetValues(*mat, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, zeroes, INSERT_VALUES));
1839:         }
1840:       }

1842:       PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY));
1843:       PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY));

1845:       PetscCall(PetscFree(zeroes));

1847:     } else { /* rsize too big, use MATPREALLOCATOR */
1848:       Mat          preallocator;
1849:       PetscScalar *vals;

1851:       PetscCall(PetscCalloc1(patch->totalDofsPerCell * patch->totalDofsPerCell, &vals));
1852:       PetscCall(MatCreate(PETSC_COMM_SELF, &preallocator));
1853:       PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
1854:       PetscCall(MatSetSizes(preallocator, rsize, rsize, rsize, rsize));
1855:       PetscCall(MatSetUp(preallocator));

1857:       for (c = 0; c < ncell; ++c) {
1858:         const PetscInt *idx = dofsArray + (offset + c) * patch->totalDofsPerCell;
1859:         PetscCall(MatSetValues(preallocator, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, vals, INSERT_VALUES));
1860:       }

1862:       if (patch->usercomputeopintfacet) {
1863:         const PetscInt *intFacetsArray = NULL;
1864:         PetscInt        i, numIntFacets, intFacetOffset;
1865:         const PetscInt *facetCells = NULL;

1867:         PetscCall(PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets));
1868:         PetscCall(PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset));
1869:         PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
1870:         PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
1871:         for (i = 0; i < numIntFacets; i++) {
1872:           const PetscInt  cell0    = facetCells[2 * (intFacetOffset + i) + 0];
1873:           const PetscInt  cell1    = facetCells[2 * (intFacetOffset + i) + 1];
1874:           const PetscInt *cell0idx = &dofsArray[(offset + cell0) * patch->totalDofsPerCell];
1875:           const PetscInt *cell1idx = &dofsArray[(offset + cell1) * patch->totalDofsPerCell];
1876:           PetscCall(MatSetValues(preallocator, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, vals, INSERT_VALUES));
1877:           PetscCall(MatSetValues(preallocator, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, vals, INSERT_VALUES));
1878:         }
1879:       }

1881:       PetscCall(PetscFree(vals));
1882:       PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
1883:       PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));
1884:       PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_TRUE, *mat));
1885:       PetscCall(MatDestroy(&preallocator));
1886:     }
1887:     PetscCall(PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0));
1888:     if (withArtificial) {
1889:       PetscCall(ISRestoreIndices(patch->dofsWithArtificial, &dofsArray));
1890:     } else {
1891:       PetscCall(ISRestoreIndices(patch->dofs, &dofsArray));
1892:     }
1893:   }
1894:   PetscCall(MatSetUp(*mat));
1895:   PetscFunctionReturn(PETSC_SUCCESS);
1896: }

1898: 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)
1899: {
1900:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
1901:   DM              dm, plex;
1902:   PetscSection    s;
1903:   const PetscInt *parray, *oarray;
1904:   PetscInt        Nf = patch->nsubspaces, Np, poff, p, f;

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

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

1936: PetscErrorCode PCPatchComputeFunction_Internal(PC pc, Vec x, Vec F, PetscInt point)
1937: {
1938:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
1939:   const PetscInt *dofsArray;
1940:   const PetscInt *dofsArrayWithAll;
1941:   const PetscInt *cellsArray;
1942:   PetscInt        ncell, offset, pStart, pEnd;

1944:   PetscFunctionBegin;
1945:   PetscCall(PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0));
1946:   PetscCheck(patch->usercomputeop, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set callback");
1947:   PetscCall(ISGetIndices(patch->dofs, &dofsArray));
1948:   PetscCall(ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll));
1949:   PetscCall(ISGetIndices(patch->cells, &cellsArray));
1950:   PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));

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

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

1972:     PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "Patch vector for Point %" PetscInt_FMT, point));
1973:     PetscCall(PetscObjectSetName((PetscObject)F, name));
1974:     PetscCall(ObjectView((PetscObject)F, patch->viewerMatrix, patch->formatMatrix));
1975:   }
1976:   PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
1977:   PetscFunctionReturn(PETSC_SUCCESS);
1978: }

1980: 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)
1981: {
1982:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
1983:   DM              dm, plex;
1984:   PetscSection    s;
1985:   const PetscInt *parray, *oarray;
1986:   PetscInt        Nf = patch->nsubspaces, Np, poff, p, f;

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

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

2018: /* This function zeros mat on entry */
2019: PetscErrorCode PCPatchComputeOperator_Internal(PC pc, Vec x, Mat mat, PetscInt point, PetscBool withArtificial)
2020: {
2021:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
2022:   const PetscInt *dofsArray;
2023:   const PetscInt *dofsArrayWithAll = NULL;
2024:   const PetscInt *cellsArray;
2025:   PetscInt        ncell, offset, pStart, pEnd, numIntFacets, intFacetOffset;
2026:   PetscBool       isNonlinear;

2028:   PetscFunctionBegin;
2029:   PetscCall(PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0));
2030:   isNonlinear = patch->isNonlinear;
2031:   PetscCheck(patch->usercomputeop, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set callback");
2032:   if (withArtificial) {
2033:     PetscCall(ISGetIndices(patch->dofsWithArtificial, &dofsArray));
2034:   } else {
2035:     PetscCall(ISGetIndices(patch->dofs, &dofsArray));
2036:   }
2037:   if (isNonlinear) PetscCall(ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll));
2038:   PetscCall(ISGetIndices(patch->cells, &cellsArray));
2039:   PetscCall(PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd));

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

2044:   PetscCall(PetscSectionGetDof(patch->cellCounts, point, &ncell));
2045:   PetscCall(PetscSectionGetOffset(patch->cellCounts, point, &offset));
2046:   if (ncell <= 0) {
2047:     PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
2048:     PetscFunctionReturn(PETSC_SUCCESS);
2049:   }
2050:   PetscCall(MatZeroEntries(mat));
2051:   if (patch->precomputeElementTensors) {
2052:     PetscInt           i;
2053:     PetscInt           ndof = patch->totalDofsPerCell;
2054:     const PetscScalar *elementTensors;

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

2086:       PetscCall(ISGetIndices(patch->intFacetsToPatchCell, &facetCells));
2087:       PetscCall(ISGetIndices(patch->intFacets, &intFacetsArray));
2088:       PetscCall(PCGetDM(pc, &dm));
2089:       PetscCall(DMConvert(dm, DMPLEX, &plex));
2090:       dm = plex;
2091:       PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, NULL));
2092:       /* FIXME: Pull this malloc out. */
2093:       PetscCall(PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofs));
2094:       if (dofsArrayWithAll) PetscCall(PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofsWithAll));
2095:       if (patch->precomputeElementTensors) {
2096:         PetscInt           nFacetDof = 2 * patch->totalDofsPerCell;
2097:         const PetscScalar *elementTensors;

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

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

2152:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
2153:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));

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

2175:     PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "Patch matrix for Point %" PetscInt_FMT, point));
2176:     PetscCall(PetscObjectSetName((PetscObject)mat, name));
2177:     PetscCall(ObjectView((PetscObject)mat, patch->viewerMatrix, patch->formatMatrix));
2178:   }
2179:   PetscCall(PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0));
2180:   PetscFunctionReturn(PETSC_SUCCESS);
2181: }

2183: static PetscErrorCode MatSetValues_PCPatch_Private(Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], const PetscScalar *v, InsertMode addv)
2184: {
2185:   Vec          data;
2186:   PetscScalar *array;
2187:   PetscInt     bs, nz, i, j, cell;

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

2211: static PetscErrorCode PCPatchPrecomputePatchTensors_Private(PC pc)
2212: {
2213:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
2214:   const PetscInt *cellsArray;
2215:   PetscInt        ncell, offset;
2216:   const PetscInt *dofMapArray;
2217:   PetscInt        i, j;
2218:   IS              dofMap;
2219:   IS              cellIS;
2220:   const PetscInt  ndof = patch->totalDofsPerCell;
2221:   Mat             vecMat;
2222:   PetscInt        cStart, cEnd;
2223:   DM              dm, plex;

2225:   PetscCall(ISGetSize(patch->cells, &ncell));
2226:   if (!ncell) { /* No cells to assemble over -> skip */
2227:     PetscFunctionReturn(PETSC_SUCCESS);
2228:   }

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

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

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

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

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

2345:   PetscFunctionReturn(PETSC_SUCCESS);
2346: }

2348: PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat, PatchScatterType scattertype)
2349: {
2350:   PC_PATCH          *patch     = (PC_PATCH *)pc->data;
2351:   const PetscScalar *xArray    = NULL;
2352:   PetscScalar       *yArray    = NULL;
2353:   const PetscInt    *gtolArray = NULL;
2354:   PetscInt           dof, offset, lidx;

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

2377:     if (mode == INSERT_VALUES) yArray[lidx] = xArray[gidx]; /* Forward */
2378:     else yArray[gidx] += xArray[lidx];                      /* Reverse */
2379:   }
2380:   if (scattertype == SCATTER_WITHARTIFICIAL) {
2381:     PetscCall(ISRestoreIndices(patch->gtolWithArtificial, &gtolArray));
2382:   } else if (scattertype == SCATTER_WITHALL) {
2383:     PetscCall(ISRestoreIndices(patch->gtolWithAll, &gtolArray));
2384:   } else {
2385:     PetscCall(ISRestoreIndices(patch->gtol, &gtolArray));
2386:   }
2387:   PetscCall(VecRestoreArrayRead(x, &xArray));
2388:   PetscCall(VecRestoreArray(y, &yArray));
2389:   PetscFunctionReturn(PETSC_SUCCESS);
2390: }

2392: static PetscErrorCode PCSetUp_PATCH_Linear(PC pc)
2393: {
2394:   PC_PATCH   *patch = (PC_PATCH *)pc->data;
2395:   const char *prefix;
2396:   PetscInt    i;

2398:   PetscFunctionBegin;
2399:   if (!pc->setupcalled) {
2400:     PetscCheck(patch->save_operators || !patch->denseinverse, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Can't have dense inverse without save operators");
2401:     if (!patch->denseinverse) {
2402:       PetscCall(PetscMalloc1(patch->npatch, &patch->solver));
2403:       PetscCall(PCGetOptionsPrefix(pc, &prefix));
2404:       for (i = 0; i < patch->npatch; ++i) {
2405:         KSP ksp;
2406:         PC  subpc;

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

2442:       PetscCall(MatGetSize(patch->mat[i], &dof, NULL));
2443:       if (dof == 0) {
2444:         patch->matWithArtificial[i] = NULL;
2445:         continue;
2446:       }

2448:       PetscCall(PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE));
2449:       PetscCall(PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE));

2451:       PetscCall(MatGetSize(matSquare, &dof, NULL));
2452:       PetscCall(ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis));
2453:       if (pc->setupcalled) {
2454:         PetscCall(MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_REUSE_MATRIX, &patch->matWithArtificial[i]));
2455:       } else {
2456:         PetscCall(MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &patch->matWithArtificial[i]));
2457:       }
2458:       PetscCall(ISDestroy(&rowis));
2459:       PetscCall(MatDestroy(&matSquare));
2460:     }
2461:   }
2462:   PetscFunctionReturn(PETSC_SUCCESS);
2463: }

2465: static PetscErrorCode PCSetUp_PATCH(PC pc)
2466: {
2467:   PC_PATCH *patch = (PC_PATCH *)pc->data;
2468:   PetscInt  i;
2469:   PetscBool isNonlinear;
2470:   PetscInt  maxDof = -1, maxDofWithArtificial = -1;

2472:   PetscFunctionBegin;
2473:   if (!pc->setupcalled) {
2474:     PetscInt pStart, pEnd, p;
2475:     PetscInt localSize;

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

2479:     isNonlinear = patch->isNonlinear;
2480:     if (!patch->nsubspaces) {
2481:       DM           dm, plex;
2482:       PetscSection s;
2483:       PetscInt     cStart, cEnd, c, Nf, f, numGlobalBcs = 0, *globalBcs, *Nb, **cellDofs;

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

2504:         PetscCall(DMGetField(dm, f, NULL, (PetscObject *)&fe));
2505:         /* PetscCall(PetscFEGetNumComponents(fe, &Nc[f])); */
2506:         PetscCall(PetscFEGetDualSpace(fe, &sp));
2507:         PetscCall(PetscDualSpaceGetDimension(sp, &Nb[f]));

2509:         PetscCall(PetscMalloc1((cEnd - cStart) * Nb[f], &cellDofs[f]));
2510:         for (c = cStart; c < cEnd; ++c) {
2511:           PetscInt *closure = NULL;
2512:           PetscInt  clSize  = 0, cl;

2514:           PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
2515:           for (cl = 0; cl < clSize * 2; cl += 2) {
2516:             const PetscInt p = closure[cl];
2517:             PetscInt       fdof, d, foff;

2519:             PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
2520:             PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
2521:             for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d;
2522:           }
2523:           PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
2524:         }
2525:         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]);
2526:       }
2527:       numGlobalBcs = 0;
2528:       for (p = pStart; p < pEnd; ++p) {
2529:         const PetscInt *ind;
2530:         PetscInt        off, cdof, d;

2532:         PetscCall(PetscSectionGetOffset(s, p, &off));
2533:         PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2534:         PetscCall(PetscSectionGetConstraintIndices(s, p, &ind));
2535:         for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d];
2536:       }

2538:       PetscCall(PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **)cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs));
2539:       for (f = 0; f < Nf; ++f) PetscCall(PetscFree(cellDofs[f]));
2540:       PetscCall(PetscFree3(Nb, cellDofs, globalBcs));
2541:       PetscCall(PCPatchSetComputeFunction(pc, PCPatchComputeFunction_DMPlex_Private, NULL));
2542:       PetscCall(PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL));
2543:       PetscCall(DMDestroy(&dm));
2544:     }

2546:     localSize = patch->subspaceOffsets[patch->nsubspaces];
2547:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localRHS));
2548:     PetscCall(VecSetUp(patch->localRHS));
2549:     PetscCall(VecDuplicate(patch->localRHS, &patch->localUpdate));
2550:     PetscCall(PCPatchCreateCellPatches(pc));
2551:     PetscCall(PCPatchCreateCellPatchDiscretisationInfo(pc));

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

2556:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscCall(PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithArtificial));
2557:     if (isNonlinear) PetscCall(PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithAll));
2558:     for (p = pStart; p < pEnd; ++p) {
2559:       PetscInt dof;

2561:       PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &dof));
2562:       maxDof = PetscMax(maxDof, dof);
2563:       if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2564:         const PetscInt *gtolArray, *gtolArrayWithArtificial = NULL;
2565:         PetscInt        numPatchDofs, offset;
2566:         PetscInt        numPatchDofsWithArtificial, offsetWithArtificial;
2567:         PetscInt        dofWithoutArtificialCounter = 0;
2568:         PetscInt       *patchWithoutArtificialToWithArtificialArray;

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

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

2577:         PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs));
2578:         if (numPatchDofs == 0) {
2579:           patch->dofMappingWithoutToWithArtificial[p - pStart] = NULL;
2580:           continue;
2581:         }

2583:         PetscCall(PetscSectionGetOffset(patch->gtolCounts, p, &offset));
2584:         PetscCall(ISGetIndices(patch->gtolWithArtificial, &gtolArrayWithArtificial));
2585:         PetscCall(PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &numPatchDofsWithArtificial));
2586:         PetscCall(PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offsetWithArtificial));

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

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

2613:         PetscCall(PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs));
2614:         if (numPatchDofs == 0) {
2615:           patch->dofMappingWithoutToWithAll[p - pStart] = NULL;
2616:           continue;
2617:         }

2619:         PetscCall(PetscSectionGetOffset(patch->gtolCounts, p, &offset));
2620:         PetscCall(ISGetIndices(patch->gtolWithAll, &gtolArrayWithAll));
2621:         PetscCall(PetscSectionGetDof(patch->gtolCountsWithAll, p, &numPatchDofsWithAll));
2622:         PetscCall(PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offsetWithAll));

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

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

2652:     /* If desired, calculate weights for dof multiplicity */
2653:     if (patch->partition_of_unity) {
2654:       PetscScalar *input  = NULL;
2655:       PetscScalar *output = NULL;
2656:       Vec          global;

2658:       PetscCall(VecDuplicate(patch->localRHS, &patch->dof_weights));
2659:       if (patch->local_composition_type == PC_COMPOSITE_ADDITIVE) {
2660:         for (i = 0; i < patch->npatch; ++i) {
2661:           PetscInt dof;

2663:           PetscCall(PetscSectionGetDof(patch->gtolCounts, i + pStart, &dof));
2664:           if (dof <= 0) continue;
2665:           PetscCall(VecSet(patch->patchRHS, 1.0));
2666:           PetscCall(PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHS, patch->dof_weights, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR));
2667:         }
2668:       } else {
2669:         /* multiplicative is actually only locally multiplicative and globally additive. need the pou where the mesh decomposition overlaps */
2670:         PetscCall(VecSet(patch->dof_weights, 1.0));
2671:       }

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

2676:       PetscCall(VecGetArray(patch->dof_weights, &input));
2677:       PetscCall(VecGetArray(global, &output));
2678:       PetscCall(PetscSFReduceBegin(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM));
2679:       PetscCall(PetscSFReduceEnd(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM));
2680:       PetscCall(VecRestoreArray(patch->dof_weights, &input));
2681:       PetscCall(VecRestoreArray(global, &output));

2683:       PetscCall(VecReciprocal(global));

2685:       PetscCall(VecGetArray(patch->dof_weights, &output));
2686:       PetscCall(VecGetArray(global, &input));
2687:       PetscCall(PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, input, output, MPI_REPLACE));
2688:       PetscCall(PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, input, output, MPI_REPLACE));
2689:       PetscCall(VecRestoreArray(patch->dof_weights, &output));
2690:       PetscCall(VecRestoreArray(global, &input));
2691:       PetscCall(VecDestroy(&global));
2692:     }
2693:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE && patch->save_operators && !patch->isNonlinear) PetscCall(PetscMalloc1(patch->npatch, &patch->matWithArtificial));
2694:   }
2695:   PetscCall((*patch->setupsolver)(pc));
2696:   PetscFunctionReturn(PETSC_SUCCESS);
2697: }

2699: static PetscErrorCode PCApply_PATCH_Linear(PC pc, PetscInt i, Vec x, Vec y)
2700: {
2701:   PC_PATCH *patch = (PC_PATCH *)pc->data;
2702:   KSP       ksp;
2703:   Mat       op;
2704:   PetscInt  m, n;

2706:   PetscFunctionBegin;
2707:   if (patch->denseinverse) {
2708:     PetscCall((*patch->densesolve)(patch->mat[i], x, y));
2709:     PetscFunctionReturn(PETSC_SUCCESS);
2710:   }
2711:   ksp = (KSP)patch->solver[i];
2712:   if (!patch->save_operators) {
2713:     Mat mat;

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

2744: static PetscErrorCode PCUpdateMultiplicative_PATCH_Linear(PC pc, PetscInt i, PetscInt pStart)
2745: {
2746:   PC_PATCH *patch = (PC_PATCH *)pc->data;
2747:   Mat       multMat;
2748:   PetscInt  n, m;

2750:   PetscFunctionBegin;

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

2780: static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y)
2781: {
2782:   PC_PATCH          *patch        = (PC_PATCH *)pc->data;
2783:   const PetscScalar *globalRHS    = NULL;
2784:   PetscScalar       *localRHS     = NULL;
2785:   PetscScalar       *globalUpdate = NULL;
2786:   const PetscInt    *bcNodes      = NULL;
2787:   PetscInt           nsweep       = patch->symmetrise_sweep ? 2 : 1;
2788:   PetscInt           start[2]     = {0, 0};
2789:   PetscInt           end[2]       = {-1, -1};
2790:   const PetscInt     inc[2]       = {1, -1};
2791:   const PetscScalar *localUpdate;
2792:   const PetscInt    *iterationSet;
2793:   PetscInt           pStart, numBcs, n, sweep, bc, j;

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

2814:   PetscCall(VecSet(patch->localUpdate, 0.0));
2815:   PetscCall(PetscSectionGetChart(patch->gtolCounts, &pStart, NULL));
2816:   PetscCall(PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0));
2817:   for (sweep = 0; sweep < nsweep; sweep++) {
2818:     for (j = start[sweep]; j * inc[sweep] < end[sweep] * inc[sweep]; j += inc[sweep]) {
2819:       PetscInt i = patch->user_patches ? iterationSet[j] : j;
2820:       PetscInt start, len;

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

2845:   /* Now we need to send the global BC values through */
2846:   PetscCall(VecGetArrayRead(x, &globalRHS));
2847:   PetscCall(ISGetSize(patch->globalBcNodes, &numBcs));
2848:   PetscCall(ISGetIndices(patch->globalBcNodes, &bcNodes));
2849:   PetscCall(VecGetLocalSize(x, &n));
2850:   for (bc = 0; bc < numBcs; ++bc) {
2851:     const PetscInt idx = bcNodes[bc];
2852:     if (idx < n) globalUpdate[idx] = globalRHS[idx];
2853:   }

2855:   PetscCall(ISRestoreIndices(patch->globalBcNodes, &bcNodes));
2856:   PetscCall(VecRestoreArrayRead(x, &globalRHS));
2857:   PetscCall(VecRestoreArray(y, &globalUpdate));

2859:   PetscCall(PetscOptionsPopGetViewerOff());
2860:   PetscCall(PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0));
2861:   PetscFunctionReturn(PETSC_SUCCESS);
2862: }

2864: static PetscErrorCode PCReset_PATCH_Linear(PC pc)
2865: {
2866:   PC_PATCH *patch = (PC_PATCH *)pc->data;
2867:   PetscInt  i;

2869:   PetscFunctionBegin;
2870:   if (patch->solver) {
2871:     for (i = 0; i < patch->npatch; ++i) PetscCall(KSPReset((KSP)patch->solver[i]));
2872:   }
2873:   PetscFunctionReturn(PETSC_SUCCESS);
2874: }

2876: static PetscErrorCode PCReset_PATCH(PC pc)
2877: {
2878:   PC_PATCH *patch = (PC_PATCH *)pc->data;
2879:   PetscInt  i;

2881:   PetscFunctionBegin;

2883:   PetscCall(PetscSFDestroy(&patch->sectionSF));
2884:   PetscCall(PetscSectionDestroy(&patch->cellCounts));
2885:   PetscCall(PetscSectionDestroy(&patch->pointCounts));
2886:   PetscCall(PetscSectionDestroy(&patch->cellNumbering));
2887:   PetscCall(PetscSectionDestroy(&patch->gtolCounts));
2888:   PetscCall(ISDestroy(&patch->gtol));
2889:   PetscCall(ISDestroy(&patch->cells));
2890:   PetscCall(ISDestroy(&patch->points));
2891:   PetscCall(ISDestroy(&patch->dofs));
2892:   PetscCall(ISDestroy(&patch->offs));
2893:   PetscCall(PetscSectionDestroy(&patch->patchSection));
2894:   PetscCall(ISDestroy(&patch->ghostBcNodes));
2895:   PetscCall(ISDestroy(&patch->globalBcNodes));
2896:   PetscCall(PetscSectionDestroy(&patch->gtolCountsWithArtificial));
2897:   PetscCall(ISDestroy(&patch->gtolWithArtificial));
2898:   PetscCall(ISDestroy(&patch->dofsWithArtificial));
2899:   PetscCall(ISDestroy(&patch->offsWithArtificial));
2900:   PetscCall(PetscSectionDestroy(&patch->gtolCountsWithAll));
2901:   PetscCall(ISDestroy(&patch->gtolWithAll));
2902:   PetscCall(ISDestroy(&patch->dofsWithAll));
2903:   PetscCall(ISDestroy(&patch->offsWithAll));
2904:   PetscCall(VecDestroy(&patch->cellMats));
2905:   PetscCall(VecDestroy(&patch->intFacetMats));
2906:   PetscCall(ISDestroy(&patch->allCells));
2907:   PetscCall(ISDestroy(&patch->intFacets));
2908:   PetscCall(ISDestroy(&patch->extFacets));
2909:   PetscCall(ISDestroy(&patch->intFacetsToPatchCell));
2910:   PetscCall(PetscSectionDestroy(&patch->intFacetCounts));
2911:   PetscCall(PetscSectionDestroy(&patch->extFacetCounts));

2913:   if (patch->dofSection)
2914:     for (i = 0; i < patch->nsubspaces; i++) PetscCall(PetscSectionDestroy(&patch->dofSection[i]));
2915:   PetscCall(PetscFree(patch->dofSection));
2916:   PetscCall(PetscFree(patch->bs));
2917:   PetscCall(PetscFree(patch->nodesPerCell));
2918:   if (patch->cellNodeMap)
2919:     for (i = 0; i < patch->nsubspaces; i++) PetscCall(PetscFree(patch->cellNodeMap[i]));
2920:   PetscCall(PetscFree(patch->cellNodeMap));
2921:   PetscCall(PetscFree(patch->subspaceOffsets));

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

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

2927:   PetscCall(VecDestroy(&patch->localRHS));
2928:   PetscCall(VecDestroy(&patch->localUpdate));
2929:   PetscCall(VecDestroy(&patch->patchRHS));
2930:   PetscCall(VecDestroy(&patch->patchUpdate));
2931:   PetscCall(VecDestroy(&patch->dof_weights));
2932:   if (patch->patch_dof_weights) {
2933:     for (i = 0; i < patch->npatch; ++i) PetscCall(VecDestroy(&patch->patch_dof_weights[i]));
2934:     PetscCall(PetscFree(patch->patch_dof_weights));
2935:   }
2936:   if (patch->mat) {
2937:     for (i = 0; i < patch->npatch; ++i) PetscCall(MatDestroy(&patch->mat[i]));
2938:     PetscCall(PetscFree(patch->mat));
2939:   }
2940:   if (patch->matWithArtificial && !patch->isNonlinear) {
2941:     for (i = 0; i < patch->npatch; ++i) PetscCall(MatDestroy(&patch->matWithArtificial[i]));
2942:     PetscCall(PetscFree(patch->matWithArtificial));
2943:   }
2944:   PetscCall(VecDestroy(&patch->patchRHSWithArtificial));
2945:   if (patch->dofMappingWithoutToWithArtificial) {
2946:     for (i = 0; i < patch->npatch; ++i) PetscCall(ISDestroy(&patch->dofMappingWithoutToWithArtificial[i]));
2947:     PetscCall(PetscFree(patch->dofMappingWithoutToWithArtificial));
2948:   }
2949:   if (patch->dofMappingWithoutToWithAll) {
2950:     for (i = 0; i < patch->npatch; ++i) PetscCall(ISDestroy(&patch->dofMappingWithoutToWithAll[i]));
2951:     PetscCall(PetscFree(patch->dofMappingWithoutToWithAll));
2952:   }
2953:   PetscCall(PetscFree(patch->sub_mat_type));
2954:   if (patch->userIS) {
2955:     for (i = 0; i < patch->npatch; ++i) PetscCall(ISDestroy(&patch->userIS[i]));
2956:     PetscCall(PetscFree(patch->userIS));
2957:   }
2958:   PetscCall(PetscFree(patch->precomputedTensorLocations));
2959:   PetscCall(PetscFree(patch->precomputedIntFacetTensorLocations));

2961:   patch->bs          = NULL;
2962:   patch->cellNodeMap = NULL;
2963:   patch->nsubspaces  = 0;
2964:   PetscCall(ISDestroy(&patch->iterationSet));

2966:   PetscCall(PetscViewerDestroy(&patch->viewerSection));
2967:   PetscFunctionReturn(PETSC_SUCCESS);
2968: }

2970: static PetscErrorCode PCDestroy_PATCH_Linear(PC pc)
2971: {
2972:   PC_PATCH *patch = (PC_PATCH *)pc->data;
2973:   PetscInt  i;

2975:   PetscFunctionBegin;
2976:   if (patch->solver) {
2977:     for (i = 0; i < patch->npatch; ++i) PetscCall(KSPDestroy((KSP *)&patch->solver[i]));
2978:     PetscCall(PetscFree(patch->solver));
2979:   }
2980:   PetscFunctionReturn(PETSC_SUCCESS);
2981: }

2983: static PetscErrorCode PCDestroy_PATCH(PC pc)
2984: {
2985:   PC_PATCH *patch = (PC_PATCH *)pc->data;

2987:   PetscFunctionBegin;
2988:   PetscCall(PCReset_PATCH(pc));
2989:   PetscCall((*patch->destroysolver)(pc));
2990:   PetscCall(PetscFree(pc->data));
2991:   PetscFunctionReturn(PETSC_SUCCESS);
2992: }

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

3006:   PetscFunctionBegin;
3007:   PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
3008:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
3009:   PetscOptionsHeadBegin(PetscOptionsObject, "Patch solver options");

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

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

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

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

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

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

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

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

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

3050:   /* If the user has set the number of subspaces, use that for the buffer size,
3051:    otherwise use a large number */
3052:   if (patch->nsubspaces <= 0) {
3053:     nfields = 128;
3054:   } else {
3055:     nfields = patch->nsubspaces;
3056:   }
3057:   PetscCall(PetscMalloc1(nfields, &ifields));
3058:   PetscCall(PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exclude_subspaces", patch->classname));
3059:   PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, option, ifields, &nfields, &flg));
3060:   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");
3061:   if (flg) {
3062:     PetscCall(PetscHSetIClear(patch->subspaces_to_exclude));
3063:     for (k = 0; k < nfields; k++) PetscCall(PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]));
3064:   }
3065:   PetscCall(PetscFree(ifields));

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

3086: static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc)
3087: {
3088:   PC_PATCH          *patch = (PC_PATCH *)pc->data;
3089:   KSPConvergedReason reason;
3090:   PetscInt           i;

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

3110: static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer)
3111: {
3112:   PC_PATCH   *patch = (PC_PATCH *)pc->data;
3113:   PetscViewer sviewer;
3114:   PetscBool   isascii;
3115:   PetscMPIInt rank;

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

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

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

3173:    Options Database Keys:
3174: + -pc_patch_cells_view   - Views the process local cell numbers for each patch
3175: . -pc_patch_points_view  - Views the process local mesh point numbers for each patch
3176: . -pc_patch_g2l_view     - Views the map between global dofs and patch local dofs for each patch
3177: . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary
3178: - -pc_patch_sub_mat_view - Views the matrix associated with each patch

3180:    Level: intermediate

3182: .seealso: `PCType`, `PCCreate()`, `PCSetType()`, `PCASM`, `PCJACOBI`, `PCPBJACOBI`, `PCVPBJACOBI`, `SNESPATCH`
3183: M*/
3184: PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc)
3185: {
3186:   PC_PATCH *patch;

3188:   PetscFunctionBegin;
3189:   PetscCall(PetscNew(&patch));

3191:   if (patch->subspaces_to_exclude) PetscCall(PetscHSetIDestroy(&patch->subspaces_to_exclude));
3192:   PetscCall(PetscHSetICreate(&patch->subspaces_to_exclude));

3194:   patch->classname   = "pc";
3195:   patch->isNonlinear = PETSC_FALSE;

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

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

3241:   PetscFunctionReturn(PETSC_SUCCESS);
3242: }