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:   PetscViewerPushFormat(viewer, format);
 15:   PetscObjectView(obj, viewer);
 16:   PetscViewerPopFormat(viewer);
 17:   return (0);
 18: }

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

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

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

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

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

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

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

 87:   PetscHSetIClear(ht);

 89:   DMPlexGetChart(dm, &pStart, &pEnd);

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

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

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

120:     PetscHSetIGetSize(ht, &htsize);
121:     PetscMalloc1(htsize, &htpoints);
122:     PetscHSetIGetElems(ht, &index, htpoints);

124:     for (i = 0; i < htsize; ++i) {
125:       PetscInt hpoint = htpoints[i];
126:       PetscInt si;

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

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

147:   return 0;
148: }

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

159:   PetscHSetIClear(ht);
160:   DMPlexGetChart(dm, &pStart, &pEnd);
161:   ISGetLocalSize(patchis, &n);
162:   ISGetIndices(patchis, &patchdata);
163:   for (i = 0; i < n; ++i) {
164:     const PetscInt ownedpoint = patchdata[i];

167:     PetscHSetIAdd(ht, ownedpoint);
168:   }
169:   ISRestoreIndices(patchis, &patchdata);
170:   return 0;
171: }

173: static PetscErrorCode PCPatchCreateDefaultSF_Private(PC pc, PetscInt n, const PetscSF *sf, const PetscInt *bs)
174: {
175:   PC_PATCH *patch = (PC_PATCH *)pc->data;
176:   PetscInt  i;

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

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

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

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

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

241:     offsets[0] = 0;
242:     for (i = 1; i < n; ++i) {
243:       PetscInt nroots;

245:       PetscSFGetGraph(sf[i - 1], &nroots, NULL, NULL, NULL);
246:       offsets[i] = offsets[i - 1] + nroots * bs[i - 1];
247:     }
248:     /* Offsets are the offsets on the current process of the global dof numbering for the subspaces. */
249:     MPI_Type_contiguous(n, MPIU_INT, &contig);
250:     MPI_Type_commit(&contig);

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

266:       PetscSFGetGraph(sf[i], &nroots, &nleaves, &local, &remote);
267:       for (j = 0; j < nleaves; ++j) {
268:         PetscInt rank = remote[j].rank;
269:         PetscInt idx, rootOffset, k;

271:         PetscHMapIGet(rankToIndex, rank, &idx);
273:         /* Offset on given rank for ith subspace */
274:         rootOffset = remoteOffsets[n * idx + i];
275:         for (k = 0; k < bs[i]; ++k) {
276:           ilocal[index]        = (local ? local[j] : j) * bs[i] + k + leafOffset;
277:           iremote[index].rank  = remote[j].rank;
278:           iremote[index].index = remote[j].index * bs[i] + k + rootOffset;
279:           ++index;
280:         }
281:       }
282:       leafOffset += nleaves * bs[i];
283:     }
284:     PetscHMapIDestroy(&rankToIndex);
285:     PetscFree(remoteOffsets);
286:     PetscSFCreate(PetscObjectComm((PetscObject)pc), &patch->sectionSF);
287:     PetscSFSetGraph(patch->sectionSF, allRoots, allLeaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
288:   }
289:   return 0;
290: }

292: PetscErrorCode PCPatchSetDenseInverse(PC pc, PetscBool flg)
293: {
294:   PC_PATCH *patch = (PC_PATCH *)pc->data;
295:   patch->denseinverse = flg;
296:   return 0;
297: }

299: PetscErrorCode PCPatchGetDenseInverse(PC pc, PetscBool *flg)
300: {
301:   PC_PATCH *patch = (PC_PATCH *)pc->data;
302:   *flg = patch->denseinverse;
303:   return 0;
304: }

306: /* TODO: Docs */
307: PetscErrorCode PCPatchSetIgnoreDim(PC pc, PetscInt dim)
308: {
309:   PC_PATCH *patch = (PC_PATCH *)pc->data;
310:   patch->ignoredim = dim;
311:   return 0;
312: }

314: /* TODO: Docs */
315: PetscErrorCode PCPatchGetIgnoreDim(PC pc, PetscInt *dim)
316: {
317:   PC_PATCH *patch = (PC_PATCH *)pc->data;
318:   *dim = patch->ignoredim;
319:   return 0;
320: }

322: /* TODO: Docs */
323: PetscErrorCode PCPatchSetSaveOperators(PC pc, PetscBool flg)
324: {
325:   PC_PATCH *patch = (PC_PATCH *)pc->data;
326:   patch->save_operators = flg;
327:   return 0;
328: }

330: /* TODO: Docs */
331: PetscErrorCode PCPatchGetSaveOperators(PC pc, PetscBool *flg)
332: {
333:   PC_PATCH *patch = (PC_PATCH *)pc->data;
334:   *flg = patch->save_operators;
335:   return 0;
336: }

338: /* TODO: Docs */
339: PetscErrorCode PCPatchSetPrecomputeElementTensors(PC pc, PetscBool flg)
340: {
341:   PC_PATCH *patch = (PC_PATCH *)pc->data;
342:   patch->precomputeElementTensors = flg;
343:   return 0;
344: }

346: /* TODO: Docs */
347: PetscErrorCode PCPatchGetPrecomputeElementTensors(PC pc, PetscBool *flg)
348: {
349:   PC_PATCH *patch = (PC_PATCH *)pc->data;
350:   *flg = patch->precomputeElementTensors;
351:   return 0;
352: }

354: /* TODO: Docs */
355: PetscErrorCode PCPatchSetPartitionOfUnity(PC pc, PetscBool flg)
356: {
357:   PC_PATCH *patch = (PC_PATCH *)pc->data;
358:   patch->partition_of_unity = flg;
359:   return 0;
360: }

362: /* TODO: Docs */
363: PetscErrorCode PCPatchGetPartitionOfUnity(PC pc, PetscBool *flg)
364: {
365:   PC_PATCH *patch = (PC_PATCH *)pc->data;
366:   *flg = patch->partition_of_unity;
367:   return 0;
368: }

370: /* TODO: Docs */
371: PetscErrorCode PCPatchSetLocalComposition(PC pc, PCCompositeType type)
372: {
373:   PC_PATCH *patch = (PC_PATCH *)pc->data;
375:   patch->local_composition_type = type;
376:   return 0;
377: }

379: /* TODO: Docs */
380: PetscErrorCode PCPatchGetLocalComposition(PC pc, PCCompositeType *type)
381: {
382:   PC_PATCH *patch = (PC_PATCH *)pc->data;
383:   *type = patch->local_composition_type;
384:   return 0;
385: }

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

392:   if (patch->sub_mat_type) PetscFree(patch->sub_mat_type);
393:   PetscStrallocpy(sub_mat_type, (char **)&patch->sub_mat_type);
394:   return 0;
395: }

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

405: /* TODO: Docs */
406: PetscErrorCode PCPatchSetCellNumbering(PC pc, PetscSection cellNumbering)
407: {
408:   PC_PATCH *patch = (PC_PATCH *)pc->data;

410:   patch->cellNumbering = cellNumbering;
411:   PetscObjectReference((PetscObject)cellNumbering);
412:   return 0;
413: }

415: /* TODO: Docs */
416: PetscErrorCode PCPatchGetCellNumbering(PC pc, PetscSection *cellNumbering)
417: {
418:   PC_PATCH *patch = (PC_PATCH *)pc->data;
419:   *cellNumbering = patch->cellNumbering;
420:   return 0;
421: }

423: /* TODO: Docs */
424: PetscErrorCode PCPatchSetConstructType(PC pc, PCPatchConstructType ctype, PetscErrorCode (*func)(PC, PetscInt *, IS **, IS *, void *), void *ctx)
425: {
426:   PC_PATCH *patch = (PC_PATCH *)pc->data;

428:   patch->ctype = ctype;
429:   switch (ctype) {
430:   case PC_PATCH_STAR:
431:     patch->user_patches     = PETSC_FALSE;
432:     patch->patchconstructop = PCPatchConstruct_Star;
433:     break;
434:   case PC_PATCH_VANKA:
435:     patch->user_patches     = PETSC_FALSE;
436:     patch->patchconstructop = PCPatchConstruct_Vanka;
437:     break;
438:   case PC_PATCH_PARDECOMP:
439:     patch->user_patches     = PETSC_FALSE;
440:     patch->patchconstructop = PCPatchConstruct_Pardecomp;
441:     break;
442:   case PC_PATCH_USER:
443:   case PC_PATCH_PYTHON:
444:     patch->user_patches     = PETSC_TRUE;
445:     patch->patchconstructop = PCPatchConstruct_User;
446:     if (func) {
447:       patch->userpatchconstructionop = func;
448:       patch->userpatchconstructctx   = ctx;
449:     }
450:     break;
451:   default:
452:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "Unknown patch construction type %" PetscInt_FMT, (PetscInt)patch->ctype);
453:   }
454:   return 0;
455: }

457: /* TODO: Docs */
458: PetscErrorCode PCPatchGetConstructType(PC pc, PCPatchConstructType *ctype, PetscErrorCode (**func)(PC, PetscInt *, IS **, IS *, void *), void **ctx)
459: {
460:   PC_PATCH *patch = (PC_PATCH *)pc->data;

462:   *ctype = patch->ctype;
463:   switch (patch->ctype) {
464:   case PC_PATCH_STAR:
465:   case PC_PATCH_VANKA:
466:   case PC_PATCH_PARDECOMP:
467:     break;
468:   case PC_PATCH_USER:
469:   case PC_PATCH_PYTHON:
470:     *func = patch->userpatchconstructionop;
471:     *ctx  = patch->userpatchconstructctx;
472:     break;
473:   default:
474:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "Unknown patch construction type %" PetscInt_FMT, (PetscInt)patch->ctype);
475:   }
476:   return 0;
477: }

479: /* TODO: Docs */
480: 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)
481: {
482:   PC_PATCH *patch = (PC_PATCH *)pc->data;
483:   DM        dm, plex;
484:   PetscSF  *sfs;
485:   PetscInt  cStart, cEnd, i, j;

487:   PCGetDM(pc, &dm);
488:   DMConvert(dm, DMPLEX, &plex);
489:   dm = plex;
490:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
491:   PetscMalloc1(nsubspaces, &sfs);
492:   PetscMalloc1(nsubspaces, &patch->dofSection);
493:   PetscMalloc1(nsubspaces, &patch->bs);
494:   PetscMalloc1(nsubspaces, &patch->nodesPerCell);
495:   PetscMalloc1(nsubspaces, &patch->cellNodeMap);
496:   PetscMalloc1(nsubspaces + 1, &patch->subspaceOffsets);

498:   patch->nsubspaces       = nsubspaces;
499:   patch->totalDofsPerCell = 0;
500:   for (i = 0; i < nsubspaces; ++i) {
501:     DMGetLocalSection(dms[i], &patch->dofSection[i]);
502:     PetscObjectReference((PetscObject)patch->dofSection[i]);
503:     DMGetSectionSF(dms[i], &sfs[i]);
504:     patch->bs[i]           = bs[i];
505:     patch->nodesPerCell[i] = nodesPerCell[i];
506:     patch->totalDofsPerCell += nodesPerCell[i] * bs[i];
507:     PetscMalloc1((cEnd - cStart) * nodesPerCell[i], &patch->cellNodeMap[i]);
508:     for (j = 0; j < (cEnd - cStart) * nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j];
509:     patch->subspaceOffsets[i] = subspaceOffsets[i];
510:   }
511:   PCPatchCreateDefaultSF_Private(pc, nsubspaces, sfs, patch->bs);
512:   PetscFree(sfs);

514:   patch->subspaceOffsets[nsubspaces] = subspaceOffsets[nsubspaces];
515:   ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);
516:   ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);
517:   DMDestroy(&dm);
518:   return 0;
519: }

521: /* TODO: Docs */
522: PetscErrorCode PCPatchSetDiscretisationInfoCombined(PC pc, DM dm, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes)
523: {
524:   PC_PATCH *patch = (PC_PATCH *)pc->data;
525:   PetscInt  cStart, cEnd, i, j;

527:   patch->combined = PETSC_TRUE;
528:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
529:   DMGetNumFields(dm, &patch->nsubspaces);
530:   PetscCalloc1(patch->nsubspaces, &patch->dofSection);
531:   PetscMalloc1(patch->nsubspaces, &patch->bs);
532:   PetscMalloc1(patch->nsubspaces, &patch->nodesPerCell);
533:   PetscMalloc1(patch->nsubspaces, &patch->cellNodeMap);
534:   PetscCalloc1(patch->nsubspaces + 1, &patch->subspaceOffsets);
535:   DMGetLocalSection(dm, &patch->dofSection[0]);
536:   PetscObjectReference((PetscObject)patch->dofSection[0]);
537:   PetscSectionGetStorageSize(patch->dofSection[0], &patch->subspaceOffsets[patch->nsubspaces]);
538:   patch->totalDofsPerCell = 0;
539:   for (i = 0; i < patch->nsubspaces; ++i) {
540:     patch->bs[i]           = 1;
541:     patch->nodesPerCell[i] = nodesPerCell[i];
542:     patch->totalDofsPerCell += nodesPerCell[i];
543:     PetscMalloc1((cEnd - cStart) * nodesPerCell[i], &patch->cellNodeMap[i]);
544:     for (j = 0; j < (cEnd - cStart) * nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j];
545:   }
546:   DMGetSectionSF(dm, &patch->sectionSF);
547:   PetscObjectReference((PetscObject)patch->sectionSF);
548:   ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);
549:   ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);
550:   return 0;
551: }

553: /*@C

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

557:   Logically collective

559:   Input Parameters:
560: + pc   - The PC
561: . func - The callback
562: - ctx  - The user context

564:   Calling sequence of func:
565: $   func (PC pc,PetscInt point,Vec x,Vec f,IS cellIS,PetscInt n,const PetscInt* dofsArray,const PetscInt* dofsArrayWithAll,void* ctx)

567: +  pc               - The PC
568: .  point            - The point
569: .  x                - The input solution (not used in linear problems)
570: .  f                - The patch residual vector
571: .  cellIS           - An array of the cell numbers
572: .  n                - The size of dofsArray
573: .  dofsArray        - The dofmap for the dofs to be solved for
574: .  dofsArrayWithAll - The dofmap for all dofs on the patch
575: -  ctx              - The user context

577:   Level: advanced

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

582: .seealso: `PCPatchSetComputeOperator()`, `PCPatchGetComputeOperator()`, `PCPatchSetDiscretisationInfo()`, `PCPatchSetComputeFunctionInteriorFacets()`
583: @*/
584: PetscErrorCode PCPatchSetComputeFunction(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Vec, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx)
585: {
586:   PC_PATCH *patch = (PC_PATCH *)pc->data;

588:   patch->usercomputef    = func;
589:   patch->usercomputefctx = ctx;
590:   return 0;
591: }

593: /*@C

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

597:   Logically collective

599:   Input Parameters:
600: + pc   - The PC
601: . func - The callback
602: - ctx  - The user context

604:   Calling sequence of func:
605: $   func (PC pc,PetscInt point,Vec x,Vec f,IS facetIS,PetscInt n,const PetscInt* dofsArray,const PetscInt* dofsArrayWithAll,void* ctx)

607: +  pc               - The PC
608: .  point            - The point
609: .  x                - The input solution (not used in linear problems)
610: .  f                - The patch residual vector
611: .  facetIS          - An array of the facet numbers
612: .  n                - The size of dofsArray
613: .  dofsArray        - The dofmap for the dofs to be solved for
614: .  dofsArrayWithAll - The dofmap for all dofs on the patch
615: -  ctx              - The user context

617:   Level: advanced

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

622: .seealso: `PCPatchSetComputeOperator()`, `PCPatchGetComputeOperator()`, `PCPatchSetDiscretisationInfo()`, `PCPatchSetComputeFunction()`
623: @*/
624: PetscErrorCode PCPatchSetComputeFunctionInteriorFacets(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Vec, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx)
625: {
626:   PC_PATCH *patch = (PC_PATCH *)pc->data;

628:   patch->usercomputefintfacet    = func;
629:   patch->usercomputefintfacetctx = ctx;
630:   return 0;
631: }

633: /*@C

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

637:   Logically collective

639:   Input Parameters:
640: + pc   - The PC
641: . func - The callback
642: - ctx  - The user context

644:   Calling sequence of func:
645: $   func (PC pc,PetscInt point,Vec x,Mat mat,IS facetIS,PetscInt n,const PetscInt* dofsArray,const PetscInt* dofsArrayWithAll,void* ctx)

647: +  pc               - The PC
648: .  point            - The point
649: .  x                - The input solution (not used in linear problems)
650: .  mat              - The patch matrix
651: .  cellIS           - An array of the cell numbers
652: .  n                - The size of dofsArray
653: .  dofsArray        - The dofmap for the dofs to be solved for
654: .  dofsArrayWithAll - The dofmap for all dofs on the patch
655: -  ctx              - The user context

657:   Level: advanced

659:   Note:
660:   The matrix entries have been set to zero before the call.

662: .seealso: `PCPatchGetComputeOperator()`, `PCPatchSetComputeFunction()`, `PCPatchSetDiscretisationInfo()`
663: @*/
664: PetscErrorCode PCPatchSetComputeOperator(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Mat, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx)
665: {
666:   PC_PATCH *patch = (PC_PATCH *)pc->data;

668:   patch->usercomputeop    = func;
669:   patch->usercomputeopctx = ctx;
670:   return 0;
671: }

673: /*@C

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

677:   Logically collective

679:   Input Parameters:
680: + pc   - The PC
681: . func - The callback
682: - ctx  - The user context

684:   Calling sequence of func:
685: $   func (PC pc,PetscInt point,Vec x,Mat mat,IS facetIS,PetscInt n,const PetscInt* dofsArray,const PetscInt* dofsArrayWithAll,void* ctx)

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

697:   Level: advanced

699:   Note:
700:   The matrix entries have been set to zero before the call.

702: .seealso: `PCPatchGetComputeOperator()`, `PCPatchSetComputeFunction()`, `PCPatchSetDiscretisationInfo()`
703: @*/
704: PetscErrorCode PCPatchSetComputeOperatorInteriorFacets(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Mat, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx)
705: {
706:   PC_PATCH *patch = (PC_PATCH *)pc->data;

708:   patch->usercomputeopintfacet    = func;
709:   patch->usercomputeopintfacetctx = ctx;
710:   return 0;
711: }

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

729:   PCGetDM(pc, &dm);
730:   DMConvert(dm, DMPLEX, &plex);
731:   dm = plex;
732:   DMPlexGetHeightStratum(dm, 1, &fBegin, &fEnd);
733:   PCPatchGetIgnoreDim(pc, &ignoredim);
734:   if (ignoredim >= 0) DMPlexGetDepthStratum(dm, ignoredim, &iStart, &iEnd);
735:   PetscHSetIClear(cht);
736:   PetscHashIterBegin(ht, hi);
737:   while (!PetscHashIterAtEnd(ht, hi)) {
738:     PetscHashIterGetKey(ht, hi, point);
739:     PetscHashIterNext(ht, hi);

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

775: static PetscErrorCode PCPatchGetGlobalDofs(PC pc, PetscSection dofSection[], PetscInt f, PetscBool combined, PetscInt p, PetscInt *dof, PetscInt *off)
776: {
777:   if (combined) {
778:     if (f < 0) {
779:       if (dof) PetscSectionGetDof(dofSection[0], p, dof);
780:       if (off) PetscSectionGetOffset(dofSection[0], p, off);
781:     } else {
782:       if (dof) PetscSectionGetFieldDof(dofSection[0], p, f, dof);
783:       if (off) PetscSectionGetFieldOffset(dofSection[0], p, f, off);
784:     }
785:   } else {
786:     if (f < 0) {
787:       PC_PATCH *patch = (PC_PATCH *)pc->data;
788:       PetscInt  fdof, g;

790:       if (dof) {
791:         *dof = 0;
792:         for (g = 0; g < patch->nsubspaces; ++g) {
793:           PetscSectionGetDof(dofSection[g], p, &fdof);
794:           *dof += fdof;
795:         }
796:       }
797:       if (off) {
798:         *off = 0;
799:         for (g = 0; g < patch->nsubspaces; ++g) {
800:           PetscSectionGetOffset(dofSection[g], p, &fdof);
801:           *off += fdof;
802:         }
803:       }
804:     } else {
805:       if (dof) PetscSectionGetDof(dofSection[f], p, dof);
806:       if (off) PetscSectionGetOffset(dofSection[f], p, off);
807:     }
808:   }
809:   return 0;
810: }

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

824:   PetscHSetIClear(dofs);
825:   for (k = 0; k < patch->nsubspaces; ++k) {
826:     PetscInt subspaceOffset = patch->subspaceOffsets[k];
827:     PetscInt bs             = patch->bs[k];
828:     PetscInt j, l;

830:     if (subspaces_to_exclude != NULL) {
831:       PetscBool should_exclude_k = PETSC_FALSE;
832:       PetscHSetIHas(*subspaces_to_exclude, k, &should_exclude_k);
833:       if (should_exclude_k) {
834:         /* only get this subspace dofs at the base entity, not any others */
835:         PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, base, &ldof, &loff);
836:         if (0 == ldof) continue;
837:         for (j = loff; j < ldof + loff; ++j) {
838:           for (l = 0; l < bs; ++l) {
839:             PetscInt dof = bs * j + l + subspaceOffset;
840:             PetscHSetIAdd(dofs, dof);
841:           }
842:         }
843:         continue; /* skip the other dofs of this subspace */
844:       }
845:     }

847:     PetscHashIterBegin(pts, hi);
848:     while (!PetscHashIterAtEnd(pts, hi)) {
849:       PetscHashIterGetKey(pts, hi, p);
850:       PetscHashIterNext(pts, hi);
851:       PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, p, &ldof, &loff);
852:       if (0 == ldof) continue;
853:       for (j = loff; j < ldof + loff; ++j) {
854:         for (l = 0; l < bs; ++l) {
855:           PetscInt dof = bs * j + l + subspaceOffset;
856:           PetscHSetIAdd(dofs, dof);
857:         }
858:       }
859:     }
860:   }
861:   return 0;
862: }

864: /* Given two hash tables A and B, compute the keys in B that are not in A, and put them in C */
865: static PetscErrorCode PCPatchComputeSetDifference_Private(PetscHSetI A, PetscHSetI B, PetscHSetI C)
866: {
867:   PetscHashIter hi;
868:   PetscInt      key;
869:   PetscBool     flg;

871:   PetscHSetIClear(C);
872:   PetscHashIterBegin(B, hi);
873:   while (!PetscHashIterAtEnd(B, hi)) {
874:     PetscHashIterGetKey(B, hi, key);
875:     PetscHashIterNext(B, hi);
876:     PetscHSetIHas(A, key, &flg);
877:     if (!flg) PetscHSetIAdd(C, key);
878:   }
879:   return 0;
880: }

882: /*
883:   PCPatchCreateCellPatches - create patches.

885:   Input Parameter:
886:   . dm - The DMPlex object defining the mesh

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

907:   /* Used to keep track of the cells in the patch. */
908:   PetscHSetICreate(&ht);
909:   PetscHSetICreate(&cht);

911:   PCGetDM(pc, &dm);
913:   DMConvert(dm, DMPLEX, &plex);
914:   dm = plex;
915:   DMPlexGetChart(dm, &pStart, &pEnd);
916:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);

918:   if (patch->user_patches) {
919:     patch->userpatchconstructionop(pc, &patch->npatch, &patch->userIS, &patch->iterationSet, patch->userpatchconstructctx);
920:     vStart = 0;
921:     vEnd   = patch->npatch;
922:   } else if (patch->ctype == PC_PATCH_PARDECOMP) {
923:     vStart = 0;
924:     vEnd   = 1;
925:   } else if (patch->codim < 0) {
926:     if (patch->dim < 0) DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
927:     else DMPlexGetDepthStratum(dm, patch->dim, &vStart, &vEnd);
928:   } else DMPlexGetHeightStratum(dm, patch->codim, &vStart, &vEnd);
929:   patch->npatch = vEnd - vStart;

931:   /* These labels mark the owned points.  We only create patches around points that this process owns. */
932:   DMHasLabel(dm, "pyop2_ghost", &isFiredrake);
933:   if (isFiredrake) {
934:     DMGetLabel(dm, "pyop2_ghost", &ghost);
935:     DMLabelCreateIndex(ghost, pStart, pEnd);
936:   } else {
937:     PetscSF sf;

939:     DMGetPointSF(dm, &sf);
940:     PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
941:     nleaves = PetscMax(nleaves, 0);
942:   }

944:   PetscSectionCreate(PETSC_COMM_SELF, &patch->cellCounts);
945:   PetscObjectSetName((PetscObject)patch->cellCounts, "Patch Cell Layout");
946:   cellCounts = patch->cellCounts;
947:   PetscSectionSetChart(cellCounts, vStart, vEnd);
948:   PetscSectionCreate(PETSC_COMM_SELF, &patch->pointCounts);
949:   PetscObjectSetName((PetscObject)patch->pointCounts, "Patch Point Layout");
950:   pointCounts = patch->pointCounts;
951:   PetscSectionSetChart(pointCounts, vStart, vEnd);
952:   PetscSectionCreate(PETSC_COMM_SELF, &patch->extFacetCounts);
953:   PetscObjectSetName((PetscObject)patch->extFacetCounts, "Patch Exterior Facet Layout");
954:   extFacetCounts = patch->extFacetCounts;
955:   PetscSectionSetChart(extFacetCounts, vStart, vEnd);
956:   PetscSectionCreate(PETSC_COMM_SELF, &patch->intFacetCounts);
957:   PetscObjectSetName((PetscObject)patch->intFacetCounts, "Patch Interior Facet Layout");
958:   intFacetCounts = patch->intFacetCounts;
959:   PetscSectionSetChart(intFacetCounts, vStart, vEnd);
960:   /* Count cells and points in the patch surrounding each entity */
961:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
962:   for (v = vStart; v < vEnd; ++v) {
963:     PetscHashIter hi;
964:     PetscInt      chtSize, loc = -1;
965:     PetscBool     flg;

967:     if (!patch->user_patches && patch->ctype != PC_PATCH_PARDECOMP) {
968:       if (ghost) DMLabelHasPoint(ghost, v, &flg);
969:       else {
970:         PetscFindInt(v, nleaves, leaves, &loc);
971:         flg = loc >= 0 ? PETSC_TRUE : PETSC_FALSE;
972:       }
973:       /* Not an owned entity, don't make a cell patch. */
974:       if (flg) continue;
975:     }

977:     patch->patchconstructop((void *)patch, dm, v, ht);
978:     PCPatchCompleteCellPatch(pc, ht, cht);
979:     PetscHSetIGetSize(cht, &chtSize);
980:     /* empty patch, continue */
981:     if (chtSize == 0) continue;

983:     /* safe because size(cht) > 0 from above */
984:     PetscHashIterBegin(cht, hi);
985:     while (!PetscHashIterAtEnd(cht, hi)) {
986:       PetscInt point, pdof;

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

1022:   PetscSectionSetUp(cellCounts);
1023:   PetscSectionGetStorageSize(cellCounts, &numCells);
1024:   PetscMalloc1(numCells, &cellsArray);
1025:   PetscSectionSetUp(pointCounts);
1026:   PetscSectionGetStorageSize(pointCounts, &numPoints);
1027:   PetscMalloc1(numPoints, &pointsArray);

1029:   PetscSectionSetUp(intFacetCounts);
1030:   PetscSectionSetUp(extFacetCounts);
1031:   PetscSectionGetStorageSize(intFacetCounts, &numIntFacets);
1032:   PetscSectionGetStorageSize(extFacetCounts, &numExtFacets);
1033:   PetscMalloc1(numIntFacets, &intFacetsArray);
1034:   PetscMalloc1(numIntFacets * 2, &intFacetsToPatchCell);
1035:   PetscMalloc1(numExtFacets, &extFacetsArray);

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

1042:     PetscSectionGetDof(pointCounts, v, &dof);
1043:     PetscSectionGetOffset(pointCounts, v, &off);
1044:     PetscSectionGetDof(cellCounts, v, &cdof);
1045:     PetscSectionGetOffset(cellCounts, v, &coff);
1046:     PetscSectionGetDof(intFacetCounts, v, &ifdof);
1047:     PetscSectionGetOffset(intFacetCounts, v, &ifoff);
1048:     PetscSectionGetDof(extFacetCounts, v, &efdof);
1049:     PetscSectionGetOffset(extFacetCounts, v, &efoff);
1050:     if (dof <= 0) continue;
1051:     patch->patchconstructop((void *)patch, dm, v, ht);
1052:     PCPatchCompleteCellPatch(pc, ht, cht);
1053:     PetscHashIterBegin(cht, hi);
1054:     while (!PetscHashIterAtEnd(cht, hi)) {
1055:       PetscInt point;

1057:       PetscHashIterGetKey(cht, hi, point);
1058:       if (fStart <= point && point < fEnd) {
1059:         const PetscInt *support;
1060:         PetscInt        supportSize, p;
1061:         PetscBool       interior = PETSC_TRUE;
1062:         DMPlexGetSupport(dm, point, &support);
1063:         DMPlexGetSupportSize(dm, point, &supportSize);
1064:         if (supportSize == 1) {
1065:           interior = PETSC_FALSE;
1066:         } else {
1067:           for (p = 0; p < supportSize; p++) {
1068:             PetscBool found;
1069:             /* FIXME: can I do this while iterating over cht? */
1070:             PetscHSetIHas(cht, support[p], &found);
1071:             if (!found) {
1072:               interior = PETSC_FALSE;
1073:               break;
1074:             }
1075:           }
1076:         }
1077:         if (interior) {
1078:           intFacetsToPatchCell[2 * (ifoff + ifn)]     = support[0];
1079:           intFacetsToPatchCell[2 * (ifoff + ifn) + 1] = support[1];
1080:           intFacetsArray[ifoff + ifn++]               = point;
1081:         } else {
1082:           extFacetsArray[efoff + efn++] = point;
1083:         }
1084:       }
1085:       PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);
1086:       if (pdof) pointsArray[off + n++] = point;
1087:       if (point >= cStart && point < cEnd) cellsArray[coff + cn++] = point;
1088:       PetscHashIterNext(cht, hi);
1089:     }

1095:     for (ifn = 0; ifn < ifdof; ifn++) {
1096:       PetscInt  cell0  = intFacetsToPatchCell[2 * (ifoff + ifn)];
1097:       PetscInt  cell1  = intFacetsToPatchCell[2 * (ifoff + ifn) + 1];
1098:       PetscBool found0 = PETSC_FALSE, found1 = PETSC_FALSE;
1099:       for (n = 0; n < cdof; n++) {
1100:         if (!found0 && cell0 == cellsArray[coff + n]) {
1101:           intFacetsToPatchCell[2 * (ifoff + ifn)] = n;
1102:           found0                                  = PETSC_TRUE;
1103:         }
1104:         if (!found1 && cell1 == cellsArray[coff + n]) {
1105:           intFacetsToPatchCell[2 * (ifoff + ifn) + 1] = n;
1106:           found1                                      = PETSC_TRUE;
1107:         }
1108:         if (found0 && found1) break;
1109:       }
1111:     }
1112:   }
1113:   PetscHSetIDestroy(&ht);
1114:   PetscHSetIDestroy(&cht);

1116:   ISCreateGeneral(PETSC_COMM_SELF, numCells, cellsArray, PETSC_OWN_POINTER, &patch->cells);
1117:   PetscObjectSetName((PetscObject)patch->cells, "Patch Cells");
1118:   if (patch->viewCells) {
1119:     ObjectView((PetscObject)patch->cellCounts, patch->viewerCells, patch->formatCells);
1120:     ObjectView((PetscObject)patch->cells, patch->viewerCells, patch->formatCells);
1121:   }
1122:   ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray, PETSC_OWN_POINTER, &patch->intFacets);
1123:   PetscObjectSetName((PetscObject)patch->intFacets, "Patch Interior Facets");
1124:   ISCreateGeneral(PETSC_COMM_SELF, 2 * numIntFacets, intFacetsToPatchCell, PETSC_OWN_POINTER, &patch->intFacetsToPatchCell);
1125:   PetscObjectSetName((PetscObject)patch->intFacetsToPatchCell, "Patch Interior Facets local support");
1126:   if (patch->viewIntFacets) {
1127:     ObjectView((PetscObject)patch->intFacetCounts, patch->viewerIntFacets, patch->formatIntFacets);
1128:     ObjectView((PetscObject)patch->intFacets, patch->viewerIntFacets, patch->formatIntFacets);
1129:     ObjectView((PetscObject)patch->intFacetsToPatchCell, patch->viewerIntFacets, patch->formatIntFacets);
1130:   }
1131:   ISCreateGeneral(PETSC_COMM_SELF, numExtFacets, extFacetsArray, PETSC_OWN_POINTER, &patch->extFacets);
1132:   PetscObjectSetName((PetscObject)patch->extFacets, "Patch Exterior Facets");
1133:   if (patch->viewExtFacets) {
1134:     ObjectView((PetscObject)patch->extFacetCounts, patch->viewerExtFacets, patch->formatExtFacets);
1135:     ObjectView((PetscObject)patch->extFacets, patch->viewerExtFacets, patch->formatExtFacets);
1136:   }
1137:   ISCreateGeneral(PETSC_COMM_SELF, numPoints, pointsArray, PETSC_OWN_POINTER, &patch->points);
1138:   PetscObjectSetName((PetscObject)patch->points, "Patch Points");
1139:   if (patch->viewPoints) {
1140:     ObjectView((PetscObject)patch->pointCounts, patch->viewerPoints, patch->formatPoints);
1141:     ObjectView((PetscObject)patch->points, patch->viewerPoints, patch->formatPoints);
1142:   }
1143:   DMDestroy(&dm);
1144:   return 0;
1145: }

1147: /*
1148:   PCPatchCreateCellPatchDiscretisationInfo - Build the dof maps for cell patches

1150:   Input Parameters:
1151:   + dm - The DMPlex object defining the mesh
1152:   . cellCounts - Section with counts of cells around each vertex
1153:   . cells - IS of the cell point indices of cells in each patch
1154:   . cellNumbering - Section mapping plex cell points to Firedrake cell indices.
1155:   . nodesPerCell - number of nodes per cell.
1156:   - cellNodeMap - map from cells to node indices (nodesPerCell * numCells)

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


1208:   PCGetDM(pc, &dm);
1209:   DMConvert(dm, DMPLEX, &plex);
1210:   dm = plex;
1211:   /* dofcounts section is cellcounts section * dofPerCell */
1212:   PetscSectionGetStorageSize(cellCounts, &numCells);
1213:   PetscSectionGetStorageSize(patch->pointCounts, &numPoints);
1214:   numDofs = numCells * totalDofsPerCell;
1215:   PetscMalloc1(numDofs, &dofsArray);
1216:   PetscMalloc1(numPoints * Nf, &offsArray);
1217:   PetscMalloc1(numDofs, &asmArray);
1218:   PetscMalloc1(numCells, &newCellsArray);
1219:   PetscSectionGetChart(cellCounts, &vStart, &vEnd);
1220:   PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCounts);
1221:   gtolCounts = patch->gtolCounts;
1222:   PetscSectionSetChart(gtolCounts, vStart, vEnd);
1223:   PetscObjectSetName((PetscObject)patch->gtolCounts, "Patch Global Index Section");

1225:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1226:     PetscMalloc1(numPoints * Nf, &offsArrayWithArtificial);
1227:     PetscMalloc1(numDofs, &asmArrayWithArtificial);
1228:     PetscMalloc1(numDofs, &dofsArrayWithArtificial);
1229:     PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithArtificial);
1230:     gtolCountsWithArtificial = patch->gtolCountsWithArtificial;
1231:     PetscSectionSetChart(gtolCountsWithArtificial, vStart, vEnd);
1232:     PetscObjectSetName((PetscObject)patch->gtolCountsWithArtificial, "Patch Global Index Section Including Artificial BCs");
1233:   }

1235:   isNonlinear = patch->isNonlinear;
1236:   if (isNonlinear) {
1237:     PetscMalloc1(numPoints * Nf, &offsArrayWithAll);
1238:     PetscMalloc1(numDofs, &asmArrayWithAll);
1239:     PetscMalloc1(numDofs, &dofsArrayWithAll);
1240:     PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithAll);
1241:     gtolCountsWithAll = patch->gtolCountsWithAll;
1242:     PetscSectionSetChart(gtolCountsWithAll, vStart, vEnd);
1243:     PetscObjectSetName((PetscObject)patch->gtolCountsWithAll, "Patch Global Index Section Including All BCs");
1244:   }

1246:   /* Outside the patch loop, get the dofs that are globally-enforced Dirichlet
1247:    conditions */
1248:   PetscHSetICreate(&globalBcs);
1249:   ISGetIndices(patch->ghostBcNodes, &bcNodes);
1250:   ISGetSize(patch->ghostBcNodes, &numBcs);
1251:   for (i = 0; i < numBcs; ++i) { PetscHSetIAdd(globalBcs, bcNodes[i]); /* these are already in concatenated numbering */ }
1252:   ISRestoreIndices(patch->ghostBcNodes, &bcNodes);
1253:   ISDestroy(&patch->ghostBcNodes); /* memory optimisation */

1255:   /* Hash tables for artificial BC construction */
1256:   PetscHSetICreate(&ownedpts);
1257:   PetscHSetICreate(&seenpts);
1258:   PetscHSetICreate(&owneddofs);
1259:   PetscHSetICreate(&seendofs);
1260:   PetscHSetICreate(&artificialbcs);

1262:   ISGetIndices(cells, &cellsArray);
1263:   ISGetIndices(points, &pointsArray);
1264:   PetscHMapICreate(&ht);
1265:   PetscHMapICreate(&htWithArtificial);
1266:   PetscHMapICreate(&htWithAll);
1267:   for (v = vStart; v < vEnd; ++v) {
1268:     PetscInt localIndex               = 0;
1269:     PetscInt localIndexWithArtificial = 0;
1270:     PetscInt localIndexWithAll        = 0;
1271:     PetscInt dof, off, i, j, k, l;

1273:     PetscHMapIClear(ht);
1274:     PetscHMapIClear(htWithArtificial);
1275:     PetscHMapIClear(htWithAll);
1276:     PetscSectionGetDof(cellCounts, v, &dof);
1277:     PetscSectionGetOffset(cellCounts, v, &off);
1278:     if (dof <= 0) continue;

1280:     /* Calculate the global numbers of the artificial BC dofs here first */
1281:     patch->patchconstructop((void *)patch, dm, v, ownedpts);
1282:     PCPatchCompleteCellPatch(pc, ownedpts, seenpts);
1283:     PCPatchGetPointDofs(pc, ownedpts, owneddofs, v, &patch->subspaces_to_exclude);
1284:     PCPatchGetPointDofs(pc, seenpts, seendofs, v, NULL);
1285:     PCPatchComputeSetDifference_Private(owneddofs, seendofs, artificialbcs);
1286:     if (patch->viewPatches) {
1287:       PetscHSetI    globalbcdofs;
1288:       PetscHashIter hi;
1289:       MPI_Comm      comm = PetscObjectComm((PetscObject)pc);

1291:       PetscHSetICreate(&globalbcdofs);
1292:       PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": owned dofs:\n", v);
1293:       PetscHashIterBegin(owneddofs, hi);
1294:       while (!PetscHashIterAtEnd(owneddofs, hi)) {
1295:         PetscInt globalDof;

1297:         PetscHashIterGetKey(owneddofs, hi, globalDof);
1298:         PetscHashIterNext(owneddofs, hi);
1299:         PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof);
1300:       }
1301:       PetscSynchronizedPrintf(comm, "\n");
1302:       PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": seen dofs:\n", v);
1303:       PetscHashIterBegin(seendofs, hi);
1304:       while (!PetscHashIterAtEnd(seendofs, hi)) {
1305:         PetscInt  globalDof;
1306:         PetscBool flg;

1308:         PetscHashIterGetKey(seendofs, hi, globalDof);
1309:         PetscHashIterNext(seendofs, hi);
1310:         PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof);

1312:         PetscHSetIHas(globalBcs, globalDof, &flg);
1313:         if (flg) PetscHSetIAdd(globalbcdofs, globalDof);
1314:       }
1315:       PetscSynchronizedPrintf(comm, "\n");
1316:       PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": global BCs:\n", v);
1317:       PetscHSetIGetSize(globalbcdofs, &numBcs);
1318:       if (numBcs > 0) {
1319:         PetscHashIterBegin(globalbcdofs, hi);
1320:         while (!PetscHashIterAtEnd(globalbcdofs, hi)) {
1321:           PetscInt globalDof;
1322:           PetscHashIterGetKey(globalbcdofs, hi, globalDof);
1323:           PetscHashIterNext(globalbcdofs, hi);
1324:           PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof);
1325:         }
1326:       }
1327:       PetscSynchronizedPrintf(comm, "\n");
1328:       PetscSynchronizedPrintf(comm, "Patch %" PetscInt_FMT ": artificial BCs:\n", v);
1329:       PetscHSetIGetSize(artificialbcs, &numBcs);
1330:       if (numBcs > 0) {
1331:         PetscHashIterBegin(artificialbcs, hi);
1332:         while (!PetscHashIterAtEnd(artificialbcs, hi)) {
1333:           PetscInt globalDof;
1334:           PetscHashIterGetKey(artificialbcs, hi, globalDof);
1335:           PetscHashIterNext(artificialbcs, hi);
1336:           PetscSynchronizedPrintf(comm, "%" PetscInt_FMT " ", globalDof);
1337:         }
1338:       }
1339:       PetscSynchronizedPrintf(comm, "\n\n");
1340:       PetscHSetIDestroy(&globalbcdofs);
1341:     }
1342:     for (k = 0; k < patch->nsubspaces; ++k) {
1343:       const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
1344:       PetscInt        nodesPerCell   = patch->nodesPerCell[k];
1345:       PetscInt        subspaceOffset = patch->subspaceOffsets[k];
1346:       PetscInt        bs             = patch->bs[k];

1348:       for (i = off; i < off + dof; ++i) {
1349:         /* Walk over the cells in this patch. */
1350:         const PetscInt c    = cellsArray[i];
1351:         PetscInt       cell = c;

1353:         /* TODO Change this to an IS */
1354:         if (cellNumbering) {
1355:           PetscSectionGetDof(cellNumbering, c, &cell);
1357:           PetscSectionGetOffset(cellNumbering, c, &cell);
1358:         }
1359:         newCellsArray[i] = cell;
1360:         for (j = 0; j < nodesPerCell; ++j) {
1361:           /* For each global dof, map it into contiguous local storage. */
1362:           const PetscInt globalDof = cellNodeMap[cell * nodesPerCell + j] * bs + subspaceOffset;
1363:           /* finally, loop over block size */
1364:           for (l = 0; l < bs; ++l) {
1365:             PetscInt  localDof;
1366:             PetscBool isGlobalBcDof, isArtificialBcDof;

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

1372:             /* if it's either, don't ever give it a local dof number */
1373:             if (isGlobalBcDof || isArtificialBcDof) {
1374:               dofsArray[globalIndex] = -1; /* don't use this in assembly in this patch */
1375:             } else {
1376:               PetscHMapIGet(ht, globalDof + l, &localDof);
1377:               if (localDof == -1) {
1378:                 localDof = localIndex++;
1379:                 PetscHMapISet(ht, globalDof + l, localDof);
1380:               }
1382:               /* And store. */
1383:               dofsArray[globalIndex] = localDof;
1384:             }

1386:             if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1387:               if (isGlobalBcDof) {
1388:                 dofsArrayWithArtificial[globalIndex] = -1; /* don't use this in assembly in this patch */
1389:               } else {
1390:                 PetscHMapIGet(htWithArtificial, globalDof + l, &localDof);
1391:                 if (localDof == -1) {
1392:                   localDof = localIndexWithArtificial++;
1393:                   PetscHMapISet(htWithArtificial, globalDof + l, localDof);
1394:                 }
1396:                 /* And store.*/
1397:                 dofsArrayWithArtificial[globalIndex] = localDof;
1398:               }
1399:             }

1401:             if (isNonlinear) {
1402:               /* Build the dofmap for the function space with _all_ dofs,
1403:    including those in any kind of boundary condition */
1404:               PetscHMapIGet(htWithAll, globalDof + l, &localDof);
1405:               if (localDof == -1) {
1406:                 localDof = localIndexWithAll++;
1407:                 PetscHMapISet(htWithAll, globalDof + l, localDof);
1408:               }
1410:               /* And store.*/
1411:               dofsArrayWithAll[globalIndex] = localDof;
1412:             }
1413:             globalIndex++;
1414:           }
1415:         }
1416:       }
1417:     }
1418:     /*How many local dofs in this patch? */
1419:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1420:       PetscHMapIGetSize(htWithArtificial, &dof);
1421:       PetscSectionSetDof(gtolCountsWithArtificial, v, dof);
1422:     }
1423:     if (isNonlinear) {
1424:       PetscHMapIGetSize(htWithAll, &dof);
1425:       PetscSectionSetDof(gtolCountsWithAll, v, dof);
1426:     }
1427:     PetscHMapIGetSize(ht, &dof);
1428:     PetscSectionSetDof(gtolCounts, v, dof);
1429:   }

1431:   DMDestroy(&dm);
1433:   PetscSectionSetUp(gtolCounts);
1434:   PetscSectionGetStorageSize(gtolCounts, &numGlobalDofs);
1435:   PetscMalloc1(numGlobalDofs, &globalDofsArray);

1437:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1438:     PetscSectionSetUp(gtolCountsWithArtificial);
1439:     PetscSectionGetStorageSize(gtolCountsWithArtificial, &numGlobalDofsWithArtificial);
1440:     PetscMalloc1(numGlobalDofsWithArtificial, &globalDofsArrayWithArtificial);
1441:   }
1442:   if (isNonlinear) {
1443:     PetscSectionSetUp(gtolCountsWithAll);
1444:     PetscSectionGetStorageSize(gtolCountsWithAll, &numGlobalDofsWithAll);
1445:     PetscMalloc1(numGlobalDofsWithAll, &globalDofsArrayWithAll);
1446:   }
1447:   /* Now populate the global to local map.  This could be merged into the above loop if we were willing to deal with reallocs. */
1448:   for (v = vStart; v < vEnd; ++v) {
1449:     PetscHashIter hi;
1450:     PetscInt      dof, off, Np, ooff, i, j, k, l;

1452:     PetscHMapIClear(ht);
1453:     PetscHMapIClear(htWithArtificial);
1454:     PetscHMapIClear(htWithAll);
1455:     PetscSectionGetDof(cellCounts, v, &dof);
1456:     PetscSectionGetOffset(cellCounts, v, &off);
1457:     PetscSectionGetDof(pointCounts, v, &Np);
1458:     PetscSectionGetOffset(pointCounts, v, &ooff);
1459:     if (dof <= 0) continue;

1461:     for (k = 0; k < patch->nsubspaces; ++k) {
1462:       const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
1463:       PetscInt        nodesPerCell   = patch->nodesPerCell[k];
1464:       PetscInt        subspaceOffset = patch->subspaceOffsets[k];
1465:       PetscInt        bs             = patch->bs[k];
1466:       PetscInt        goff;

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

1473:         if (cellNumbering) PetscSectionGetOffset(cellNumbering, c, &cell);
1474:         for (j = 0; j < nodesPerCell; ++j) {
1475:           for (l = 0; l < bs; ++l) {
1476:             const PetscInt globalDof = cellNodeMap[cell * nodesPerCell + j] * bs + l + subspaceOffset;
1477:             const PetscInt localDof  = dofsArray[key];
1478:             if (localDof >= 0) PetscHMapISet(ht, globalDof, localDof);
1479:             if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1480:               const PetscInt localDofWithArtificial = dofsArrayWithArtificial[key];
1481:               if (localDofWithArtificial >= 0) PetscHMapISet(htWithArtificial, globalDof, localDofWithArtificial);
1482:             }
1483:             if (isNonlinear) {
1484:               const PetscInt localDofWithAll = dofsArrayWithAll[key];
1485:               if (localDofWithAll >= 0) PetscHMapISet(htWithAll, globalDof, localDofWithAll);
1486:             }
1487:             key++;
1488:           }
1489:         }
1490:       }

1492:       /* Shove it in the output data structure. */
1493:       PetscSectionGetOffset(gtolCounts, v, &goff);
1494:       PetscHashIterBegin(ht, hi);
1495:       while (!PetscHashIterAtEnd(ht, hi)) {
1496:         PetscInt globalDof, localDof;

1498:         PetscHashIterGetKey(ht, hi, globalDof);
1499:         PetscHashIterGetVal(ht, hi, localDof);
1500:         if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof;
1501:         PetscHashIterNext(ht, hi);
1502:       }

1504:       if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1505:         PetscSectionGetOffset(gtolCountsWithArtificial, v, &goff);
1506:         PetscHashIterBegin(htWithArtificial, hi);
1507:         while (!PetscHashIterAtEnd(htWithArtificial, hi)) {
1508:           PetscInt globalDof, localDof;
1509:           PetscHashIterGetKey(htWithArtificial, hi, globalDof);
1510:           PetscHashIterGetVal(htWithArtificial, hi, localDof);
1511:           if (globalDof >= 0) globalDofsArrayWithArtificial[goff + localDof] = globalDof;
1512:           PetscHashIterNext(htWithArtificial, hi);
1513:         }
1514:       }
1515:       if (isNonlinear) {
1516:         PetscSectionGetOffset(gtolCountsWithAll, v, &goff);
1517:         PetscHashIterBegin(htWithAll, hi);
1518:         while (!PetscHashIterAtEnd(htWithAll, hi)) {
1519:           PetscInt globalDof, localDof;
1520:           PetscHashIterGetKey(htWithAll, hi, globalDof);
1521:           PetscHashIterGetVal(htWithAll, hi, localDof);
1522:           if (globalDof >= 0) globalDofsArrayWithAll[goff + localDof] = globalDof;
1523:           PetscHashIterNext(htWithAll, hi);
1524:         }
1525:       }

1527:       for (p = 0; p < Np; ++p) {
1528:         const PetscInt point = pointsArray[ooff + p];
1529:         PetscInt       globalDof, localDof;

1531:         PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof);
1532:         PetscHMapIGet(ht, globalDof, &localDof);
1533:         offsArray[(ooff + p) * Nf + k] = localDof;
1534:         if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1535:           PetscHMapIGet(htWithArtificial, globalDof, &localDof);
1536:           offsArrayWithArtificial[(ooff + p) * Nf + k] = localDof;
1537:         }
1538:         if (isNonlinear) {
1539:           PetscHMapIGet(htWithAll, globalDof, &localDof);
1540:           offsArrayWithAll[(ooff + p) * Nf + k] = localDof;
1541:         }
1542:       }
1543:     }

1545:     PetscHSetIDestroy(&globalBcs);
1546:     PetscHSetIDestroy(&ownedpts);
1547:     PetscHSetIDestroy(&seenpts);
1548:     PetscHSetIDestroy(&owneddofs);
1549:     PetscHSetIDestroy(&seendofs);
1550:     PetscHSetIDestroy(&artificialbcs);

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

1562:         if (cellNumbering) PetscSectionGetOffset(cellNumbering, c, &cell);
1563:         for (k = 0; k < patch->nsubspaces; ++k) {
1564:           const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
1565:           PetscInt        nodesPerCell   = patch->nodesPerCell[k];
1566:           PetscInt        subspaceOffset = patch->subspaceOffsets[k];
1567:           PetscInt        bs             = patch->bs[k];

1569:           for (j = 0; j < nodesPerCell; ++j) {
1570:             for (l = 0; l < bs; ++l) {
1571:               const PetscInt globalDof = cellNodeMap[cell * nodesPerCell + j] * bs + l + subspaceOffset;
1572:               PetscInt       localDof;

1574:               PetscHMapIGet(ht, globalDof, &localDof);
1575:               /* If it's not in the hash table, i.e. is a BC dof,
1576:    then the PetscHSetIMap above gives -1, which matches
1577:    exactly the convention for PETSc's matrix assembly to
1578:    ignore the dof. So we don't need to do anything here */
1579:               asmArray[asmKey] = localDof;
1580:               if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1581:                 PetscHMapIGet(htWithArtificial, globalDof, &localDof);
1582:                 asmArrayWithArtificial[asmKey] = localDof;
1583:               }
1584:               if (isNonlinear) {
1585:                 PetscHMapIGet(htWithAll, globalDof, &localDof);
1586:                 asmArrayWithAll[asmKey] = localDof;
1587:               }
1588:               asmKey++;
1589:             }
1590:           }
1591:         }
1592:       }
1593:     }
1594:   }
1595:   if (1 == patch->nsubspaces) {
1596:     PetscArraycpy(asmArray, dofsArray, numDofs);
1597:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscArraycpy(asmArrayWithArtificial, dofsArrayWithArtificial, numDofs);
1598:     if (isNonlinear) PetscArraycpy(asmArrayWithAll, dofsArrayWithAll, numDofs);
1599:   }

1601:   PetscHMapIDestroy(&ht);
1602:   PetscHMapIDestroy(&htWithArtificial);
1603:   PetscHMapIDestroy(&htWithAll);
1604:   ISRestoreIndices(cells, &cellsArray);
1605:   ISRestoreIndices(points, &pointsArray);
1606:   PetscFree(dofsArray);
1607:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscFree(dofsArrayWithArtificial);
1608:   if (isNonlinear) PetscFree(dofsArrayWithAll);
1609:   /* Create placeholder section for map from points to patch dofs */
1610:   PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection);
1611:   PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces);
1612:   if (patch->combined) {
1613:     PetscInt numFields;
1614:     PetscSectionGetNumFields(patch->dofSection[0], &numFields);
1616:     PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd);
1617:     PetscSectionSetChart(patch->patchSection, pStart, pEnd);
1618:     for (p = pStart; p < pEnd; ++p) {
1619:       PetscInt dof, fdof, f;

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

1672: static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat, PetscBool withArtificial)
1673: {
1674:   PC_PATCH   *patch = (PC_PATCH *)pc->data;
1675:   PetscBool   flg;
1676:   PetscInt    csize, rsize;
1677:   const char *prefix = NULL;

1679:   if (withArtificial) {
1680:     /* would be nice if we could create a rectangular matrix of size numDofsWithArtificial x numDofs here */
1681:     PetscInt pStart;
1682:     PetscSectionGetChart(patch->gtolCountsWithArtificial, &pStart, NULL);
1683:     PetscSectionGetDof(patch->gtolCountsWithArtificial, point + pStart, &rsize);
1684:     csize = rsize;
1685:   } else {
1686:     PetscInt pStart;
1687:     PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);
1688:     PetscSectionGetDof(patch->gtolCounts, point + pStart, &rsize);
1689:     csize = rsize;
1690:   }

1692:   MatCreate(PETSC_COMM_SELF, mat);
1693:   PCGetOptionsPrefix(pc, &prefix);
1694:   MatSetOptionsPrefix(*mat, prefix);
1695:   MatAppendOptionsPrefix(*mat, "pc_patch_sub_");
1696:   if (patch->sub_mat_type) MatSetType(*mat, patch->sub_mat_type);
1697:   else if (!patch->sub_mat_type) MatSetType(*mat, MATDENSE);
1698:   MatSetSizes(*mat, rsize, csize, rsize, csize);
1699:   PetscObjectTypeCompare((PetscObject)*mat, MATDENSE, &flg);
1700:   if (!flg) PetscObjectTypeCompare((PetscObject)*mat, MATSEQDENSE, &flg);
1701:   /* Sparse patch matrices */
1702:   if (!flg) {
1703:     PetscBT         bt;
1704:     PetscInt       *dnnz      = NULL;
1705:     const PetscInt *dofsArray = NULL;
1706:     PetscInt        pStart, pEnd, ncell, offset, c, i, j;

1708:     if (withArtificial) {
1709:       ISGetIndices(patch->dofsWithArtificial, &dofsArray);
1710:     } else {
1711:       ISGetIndices(patch->dofs, &dofsArray);
1712:     }
1713:     PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);
1714:     point += pStart;
1716:     PetscSectionGetDof(patch->cellCounts, point, &ncell);
1717:     PetscSectionGetOffset(patch->cellCounts, point, &offset);
1718:     PetscLogEventBegin(PC_Patch_Prealloc, pc, 0, 0, 0);
1719:     /* A PetscBT uses N^2 bits to store the sparsity pattern on a
1720:    * patch. This is probably OK if the patches are not too big,
1721:    * but uses too much memory. We therefore switch based on rsize. */
1722:     if (rsize < 3000) { /* FIXME: I picked this switch value out of my hat */
1723:       PetscScalar *zeroes;
1724:       PetscInt     rows;

1726:       PetscCalloc1(rsize, &dnnz);
1727:       PetscBTCreate(rsize * rsize, &bt);
1728:       for (c = 0; c < ncell; ++c) {
1729:         const PetscInt *idx = dofsArray + (offset + c) * patch->totalDofsPerCell;
1730:         for (i = 0; i < patch->totalDofsPerCell; ++i) {
1731:           const PetscInt row = idx[i];
1732:           if (row < 0) continue;
1733:           for (j = 0; j < patch->totalDofsPerCell; ++j) {
1734:             const PetscInt col = idx[j];
1735:             const PetscInt key = row * rsize + col;
1736:             if (col < 0) continue;
1737:             if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1738:           }
1739:         }
1740:       }

1742:       if (patch->usercomputeopintfacet) {
1743:         const PetscInt *intFacetsArray = NULL;
1744:         PetscInt        i, numIntFacets, intFacetOffset;
1745:         const PetscInt *facetCells = NULL;

1747:         PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1748:         PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1749:         ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1750:         ISGetIndices(patch->intFacets, &intFacetsArray);
1751:         for (i = 0; i < numIntFacets; i++) {
1752:           const PetscInt cell0 = facetCells[2 * (intFacetOffset + i) + 0];
1753:           const PetscInt cell1 = facetCells[2 * (intFacetOffset + i) + 1];
1754:           PetscInt       celli, cellj;

1756:           for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1757:             const PetscInt row = dofsArray[(offset + cell0) * patch->totalDofsPerCell + celli];
1758:             if (row < 0) continue;
1759:             for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1760:               const PetscInt col = dofsArray[(offset + cell1) * patch->totalDofsPerCell + cellj];
1761:               const PetscInt key = row * rsize + col;
1762:               if (col < 0) continue;
1763:               if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1764:             }
1765:           }

1767:           for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1768:             const PetscInt row = dofsArray[(offset + cell1) * patch->totalDofsPerCell + celli];
1769:             if (row < 0) continue;
1770:             for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1771:               const PetscInt col = dofsArray[(offset + cell0) * patch->totalDofsPerCell + cellj];
1772:               const PetscInt key = row * rsize + col;
1773:               if (col < 0) continue;
1774:               if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1775:             }
1776:           }
1777:         }
1778:       }
1779:       PetscBTDestroy(&bt);
1780:       MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL);
1781:       PetscFree(dnnz);

1783:       PetscCalloc1(patch->totalDofsPerCell * patch->totalDofsPerCell, &zeroes);
1784:       for (c = 0; c < ncell; ++c) {
1785:         const PetscInt *idx = &dofsArray[(offset + c) * patch->totalDofsPerCell];
1786:         MatSetValues(*mat, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, zeroes, INSERT_VALUES);
1787:       }
1788:       MatGetLocalSize(*mat, &rows, NULL);
1789:       for (i = 0; i < rows; ++i) MatSetValues(*mat, 1, &i, 1, &i, zeroes, INSERT_VALUES);

1791:       if (patch->usercomputeopintfacet) {
1792:         const PetscInt *intFacetsArray = NULL;
1793:         PetscInt        i, numIntFacets, intFacetOffset;
1794:         const PetscInt *facetCells = NULL;

1796:         PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1797:         PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1798:         ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1799:         ISGetIndices(patch->intFacets, &intFacetsArray);
1800:         for (i = 0; i < numIntFacets; i++) {
1801:           const PetscInt  cell0    = facetCells[2 * (intFacetOffset + i) + 0];
1802:           const PetscInt  cell1    = facetCells[2 * (intFacetOffset + i) + 1];
1803:           const PetscInt *cell0idx = &dofsArray[(offset + cell0) * patch->totalDofsPerCell];
1804:           const PetscInt *cell1idx = &dofsArray[(offset + cell1) * patch->totalDofsPerCell];
1805:           MatSetValues(*mat, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, zeroes, INSERT_VALUES);
1806:           MatSetValues(*mat, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, zeroes, INSERT_VALUES);
1807:         }
1808:       }

1810:       MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
1811:       MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);

1813:       PetscFree(zeroes);

1815:     } else { /* rsize too big, use MATPREALLOCATOR */
1816:       Mat          preallocator;
1817:       PetscScalar *vals;

1819:       PetscCalloc1(patch->totalDofsPerCell * patch->totalDofsPerCell, &vals);
1820:       MatCreate(PETSC_COMM_SELF, &preallocator);
1821:       MatSetType(preallocator, MATPREALLOCATOR);
1822:       MatSetSizes(preallocator, rsize, rsize, rsize, rsize);
1823:       MatSetUp(preallocator);

1825:       for (c = 0; c < ncell; ++c) {
1826:         const PetscInt *idx = dofsArray + (offset + c) * patch->totalDofsPerCell;
1827:         MatSetValues(preallocator, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, vals, INSERT_VALUES);
1828:       }

1830:       if (patch->usercomputeopintfacet) {
1831:         const PetscInt *intFacetsArray = NULL;
1832:         PetscInt        i, numIntFacets, intFacetOffset;
1833:         const PetscInt *facetCells = NULL;

1835:         PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1836:         PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1837:         ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1838:         ISGetIndices(patch->intFacets, &intFacetsArray);
1839:         for (i = 0; i < numIntFacets; i++) {
1840:           const PetscInt  cell0    = facetCells[2 * (intFacetOffset + i) + 0];
1841:           const PetscInt  cell1    = facetCells[2 * (intFacetOffset + i) + 1];
1842:           const PetscInt *cell0idx = &dofsArray[(offset + cell0) * patch->totalDofsPerCell];
1843:           const PetscInt *cell1idx = &dofsArray[(offset + cell1) * patch->totalDofsPerCell];
1844:           MatSetValues(preallocator, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, vals, INSERT_VALUES);
1845:           MatSetValues(preallocator, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, vals, INSERT_VALUES);
1846:         }
1847:       }

1849:       PetscFree(vals);
1850:       MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1851:       MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1852:       MatPreallocatorPreallocate(preallocator, PETSC_TRUE, *mat);
1853:       MatDestroy(&preallocator);
1854:     }
1855:     PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0);
1856:     if (withArtificial) {
1857:       ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);
1858:     } else {
1859:       ISRestoreIndices(patch->dofs, &dofsArray);
1860:     }
1861:   }
1862:   MatSetUp(*mat);
1863:   return 0;
1864: }

1866: 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)
1867: {
1868:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
1869:   DM              dm, plex;
1870:   PetscSection    s;
1871:   const PetscInt *parray, *oarray;
1872:   PetscInt        Nf = patch->nsubspaces, Np, poff, p, f;

1875:   PCGetDM(pc, &dm);
1876:   DMConvert(dm, DMPLEX, &plex);
1877:   dm = plex;
1878:   DMGetLocalSection(dm, &s);
1879:   /* Set offset into patch */
1880:   PetscSectionGetDof(patch->pointCounts, patchNum, &Np);
1881:   PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);
1882:   ISGetIndices(patch->points, &parray);
1883:   ISGetIndices(patch->offs, &oarray);
1884:   for (f = 0; f < Nf; ++f) {
1885:     for (p = 0; p < Np; ++p) {
1886:       const PetscInt point = parray[poff + p];
1887:       PetscInt       dof;

1889:       PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);
1890:       PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff + p) * Nf + f]);
1891:       if (patch->nsubspaces == 1) PetscSectionSetOffset(patch->patchSection, point, oarray[(poff + p) * Nf + f]);
1892:       else PetscSectionSetOffset(patch->patchSection, point, -1);
1893:     }
1894:   }
1895:   ISRestoreIndices(patch->points, &parray);
1896:   ISRestoreIndices(patch->offs, &oarray);
1897:   if (patch->viewSection) ObjectView((PetscObject)patch->patchSection, patch->viewerSection, patch->formatSection);
1898:   DMPlexComputeResidual_Patch_Internal(dm, patch->patchSection, cellIS, 0.0, x, NULL, F, ctx);
1899:   DMDestroy(&dm);
1900:   return 0;
1901: }

1903: PetscErrorCode PCPatchComputeFunction_Internal(PC pc, Vec x, Vec F, PetscInt point)
1904: {
1905:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
1906:   const PetscInt *dofsArray;
1907:   const PetscInt *dofsArrayWithAll;
1908:   const PetscInt *cellsArray;
1909:   PetscInt        ncell, offset, pStart, pEnd;

1911:   PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);
1913:   ISGetIndices(patch->dofs, &dofsArray);
1914:   ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);
1915:   ISGetIndices(patch->cells, &cellsArray);
1916:   PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);

1918:   point += pStart;

1921:   PetscSectionGetDof(patch->cellCounts, point, &ncell);
1922:   PetscSectionGetOffset(patch->cellCounts, point, &offset);
1923:   if (ncell <= 0) {
1924:     PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
1925:     return 0;
1926:   }
1927:   VecSet(F, 0.0);
1928:   /* Cannot reuse the same IS because the geometry info is being cached in it */
1929:   ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);
1930:   PetscCallBack("PCPatch callback", patch->usercomputef(pc, point, x, F, patch->cellIS, ncell * patch->totalDofsPerCell, dofsArray + offset * patch->totalDofsPerCell, dofsArrayWithAll + offset * patch->totalDofsPerCell, patch->usercomputefctx));
1931:   ISDestroy(&patch->cellIS);
1932:   ISRestoreIndices(patch->dofs, &dofsArray);
1933:   ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);
1934:   ISRestoreIndices(patch->cells, &cellsArray);
1935:   if (patch->viewMatrix) {
1936:     char name[PETSC_MAX_PATH_LEN];

1938:     PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "Patch vector for Point %" PetscInt_FMT, point);
1939:     PetscObjectSetName((PetscObject)F, name);
1940:     ObjectView((PetscObject)F, patch->viewerMatrix, patch->formatMatrix);
1941:   }
1942:   PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
1943:   return 0;
1944: }

1946: 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)
1947: {
1948:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
1949:   DM              dm, plex;
1950:   PetscSection    s;
1951:   const PetscInt *parray, *oarray;
1952:   PetscInt        Nf = patch->nsubspaces, Np, poff, p, f;

1954:   PCGetDM(pc, &dm);
1955:   DMConvert(dm, DMPLEX, &plex);
1956:   dm = plex;
1957:   DMGetLocalSection(dm, &s);
1958:   /* Set offset into patch */
1959:   PetscSectionGetDof(patch->pointCounts, patchNum, &Np);
1960:   PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);
1961:   ISGetIndices(patch->points, &parray);
1962:   ISGetIndices(patch->offs, &oarray);
1963:   for (f = 0; f < Nf; ++f) {
1964:     for (p = 0; p < Np; ++p) {
1965:       const PetscInt point = parray[poff + p];
1966:       PetscInt       dof;

1968:       PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);
1969:       PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff + p) * Nf + f]);
1970:       if (patch->nsubspaces == 1) PetscSectionSetOffset(patch->patchSection, point, oarray[(poff + p) * Nf + f]);
1971:       else PetscSectionSetOffset(patch->patchSection, point, -1);
1972:     }
1973:   }
1974:   ISRestoreIndices(patch->points, &parray);
1975:   ISRestoreIndices(patch->offs, &oarray);
1976:   if (patch->viewSection) ObjectView((PetscObject)patch->patchSection, patch->viewerSection, patch->formatSection);
1977:   /* TODO Shut off MatViewFromOptions() in MatAssemblyEnd() here */
1978:   DMPlexComputeJacobian_Patch_Internal(dm, patch->patchSection, patch->patchSection, cellIS, 0.0, 0.0, x, NULL, J, J, ctx);
1979:   DMDestroy(&dm);
1980:   return 0;
1981: }

1983: /* This function zeros mat on entry */
1984: PetscErrorCode PCPatchComputeOperator_Internal(PC pc, Vec x, Mat mat, PetscInt point, PetscBool withArtificial)
1985: {
1986:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
1987:   const PetscInt *dofsArray;
1988:   const PetscInt *dofsArrayWithAll = NULL;
1989:   const PetscInt *cellsArray;
1990:   PetscInt        ncell, offset, pStart, pEnd, numIntFacets, intFacetOffset;
1991:   PetscBool       isNonlinear;

1993:   PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);
1994:   isNonlinear = patch->isNonlinear;
1996:   if (withArtificial) {
1997:     ISGetIndices(patch->dofsWithArtificial, &dofsArray);
1998:   } else {
1999:     ISGetIndices(patch->dofs, &dofsArray);
2000:   }
2001:   if (isNonlinear) ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);
2002:   ISGetIndices(patch->cells, &cellsArray);
2003:   PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);

2005:   point += pStart;

2008:   PetscSectionGetDof(patch->cellCounts, point, &ncell);
2009:   PetscSectionGetOffset(patch->cellCounts, point, &offset);
2010:   if (ncell <= 0) {
2011:     PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2012:     return 0;
2013:   }
2014:   MatZeroEntries(mat);
2015:   if (patch->precomputeElementTensors) {
2016:     PetscInt           i;
2017:     PetscInt           ndof = patch->totalDofsPerCell;
2018:     const PetscScalar *elementTensors;

2020:     VecGetArrayRead(patch->cellMats, &elementTensors);
2021:     for (i = 0; i < ncell; i++) {
2022:       const PetscInt     cell = cellsArray[i + offset];
2023:       const PetscInt    *idx  = dofsArray + (offset + i) * ndof;
2024:       const PetscScalar *v    = elementTensors + patch->precomputedTensorLocations[cell] * ndof * ndof;
2025:       MatSetValues(mat, ndof, idx, ndof, idx, v, ADD_VALUES);
2026:     }
2027:     VecRestoreArrayRead(patch->cellMats, &elementTensors);
2028:     MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
2029:     MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
2030:   } else {
2031:     /* Cannot reuse the same IS because the geometry info is being cached in it */
2032:     ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);
2033:     PetscCallBack("PCPatch callback",
2034:                   patch->usercomputeop(pc, point, x, mat, patch->cellIS, ncell * patch->totalDofsPerCell, dofsArray + offset * patch->totalDofsPerCell, dofsArrayWithAll ? dofsArrayWithAll + offset * patch->totalDofsPerCell : NULL, patch->usercomputeopctx));
2035:   }
2036:   if (patch->usercomputeopintfacet) {
2037:     PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
2038:     PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
2039:     if (numIntFacets > 0) {
2040:       /* For each interior facet, grab the two cells (in local numbering, and concatenate dof numberings for those cells) */
2041:       PetscInt       *facetDofs = NULL, *facetDofsWithAll = NULL;
2042:       const PetscInt *intFacetsArray = NULL;
2043:       PetscInt        idx            = 0;
2044:       PetscInt        i, c, d;
2045:       PetscInt        fStart;
2046:       DM              dm, plex;
2047:       IS              facetIS    = NULL;
2048:       const PetscInt *facetCells = NULL;

2050:       ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
2051:       ISGetIndices(patch->intFacets, &intFacetsArray);
2052:       PCGetDM(pc, &dm);
2053:       DMConvert(dm, DMPLEX, &plex);
2054:       dm = plex;
2055:       DMPlexGetHeightStratum(dm, 1, &fStart, NULL);
2056:       /* FIXME: Pull this malloc out. */
2057:       PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofs);
2058:       if (dofsArrayWithAll) PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofsWithAll);
2059:       if (patch->precomputeElementTensors) {
2060:         PetscInt           nFacetDof = 2 * patch->totalDofsPerCell;
2061:         const PetscScalar *elementTensors;

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

2065:         for (i = 0; i < numIntFacets; i++) {
2066:           const PetscInt     facet = intFacetsArray[i + intFacetOffset];
2067:           const PetscScalar *v     = elementTensors + patch->precomputedIntFacetTensorLocations[facet - fStart] * nFacetDof * nFacetDof;
2068:           idx                      = 0;
2069:           /*
2070:      0--1
2071:      |\-|
2072:      |+\|
2073:      2--3
2074:      [0, 2, 3, 0, 1, 3]
2075:    */
2076:           for (c = 0; c < 2; c++) {
2077:             const PetscInt cell = facetCells[2 * (intFacetOffset + i) + c];
2078:             for (d = 0; d < patch->totalDofsPerCell; d++) {
2079:               facetDofs[idx] = dofsArray[(offset + cell) * patch->totalDofsPerCell + d];
2080:               idx++;
2081:             }
2082:           }
2083:           MatSetValues(mat, nFacetDof, facetDofs, nFacetDof, facetDofs, v, ADD_VALUES);
2084:         }
2085:         VecRestoreArrayRead(patch->intFacetMats, &elementTensors);
2086:       } else {
2087:         /*
2088:      0--1
2089:      |\-|
2090:      |+\|
2091:      2--3
2092:      [0, 2, 3, 0, 1, 3]
2093:    */
2094:         for (i = 0; i < numIntFacets; i++) {
2095:           for (c = 0; c < 2; c++) {
2096:             const PetscInt cell = facetCells[2 * (intFacetOffset + i) + c];
2097:             for (d = 0; d < patch->totalDofsPerCell; d++) {
2098:               facetDofs[idx] = dofsArray[(offset + cell) * patch->totalDofsPerCell + d];
2099:               if (dofsArrayWithAll) facetDofsWithAll[idx] = dofsArrayWithAll[(offset + cell) * patch->totalDofsPerCell + d];
2100:               idx++;
2101:             }
2102:           }
2103:         }
2104:         ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray + intFacetOffset, PETSC_USE_POINTER, &facetIS);
2105:         patch->usercomputeopintfacet(pc, point, x, mat, facetIS, 2 * numIntFacets * patch->totalDofsPerCell, facetDofs, facetDofsWithAll, patch->usercomputeopintfacetctx);
2106:         ISDestroy(&facetIS);
2107:       }
2108:       ISRestoreIndices(patch->intFacetsToPatchCell, &facetCells);
2109:       ISRestoreIndices(patch->intFacets, &intFacetsArray);
2110:       PetscFree(facetDofs);
2111:       PetscFree(facetDofsWithAll);
2112:       DMDestroy(&dm);
2113:     }
2114:   }

2116:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
2117:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);

2119:   if (!(withArtificial || isNonlinear) && patch->denseinverse) {
2120:     MatFactorInfo info;
2121:     PetscBool     flg;
2122:     PetscObjectTypeCompare((PetscObject)mat, MATSEQDENSE, &flg);
2124:     MatFactorInfoInitialize(&info);
2125:     MatLUFactor(mat, NULL, NULL, &info);
2126:     MatSeqDenseInvertFactors_Private(mat);
2127:   }
2128:   ISDestroy(&patch->cellIS);
2129:   if (withArtificial) {
2130:     ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);
2131:   } else {
2132:     ISRestoreIndices(patch->dofs, &dofsArray);
2133:   }
2134:   if (isNonlinear) ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);
2135:   ISRestoreIndices(patch->cells, &cellsArray);
2136:   if (patch->viewMatrix) {
2137:     char name[PETSC_MAX_PATH_LEN];

2139:     PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "Patch matrix for Point %" PetscInt_FMT, point);
2140:     PetscObjectSetName((PetscObject)mat, name);
2141:     ObjectView((PetscObject)mat, patch->viewerMatrix, patch->formatMatrix);
2142:   }
2143:   PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2144:   return 0;
2145: }

2147: static PetscErrorCode MatSetValues_PCPatch_Private(Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], const PetscScalar *v, InsertMode addv)
2148: {
2149:   Vec          data;
2150:   PetscScalar *array;
2151:   PetscInt     bs, nz, i, j, cell;

2153:   MatShellGetContext(mat, &data);
2154:   VecGetBlockSize(data, &bs);
2155:   VecGetSize(data, &nz);
2156:   VecGetArray(data, &array);
2158:   cell = (PetscInt)(idxm[0] / bs); /* use the fact that this is called once per cell */
2159:   for (i = 0; i < m; i++) {
2161:     for (j = 0; j < n; j++) {
2162:       const PetscScalar v_ = v[i * bs + j];
2163:       /* Indexing is special to the data structure we have! */
2164:       if (addv == INSERT_VALUES) {
2165:         array[cell * bs * bs + i * bs + j] = v_;
2166:       } else {
2167:         array[cell * bs * bs + i * bs + j] += v_;
2168:       }
2169:     }
2170:   }
2171:   VecRestoreArray(data, &array);
2172:   return 0;
2173: }

2175: static PetscErrorCode PCPatchPrecomputePatchTensors_Private(PC pc)
2176: {
2177:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
2178:   const PetscInt *cellsArray;
2179:   PetscInt        ncell, offset;
2180:   const PetscInt *dofMapArray;
2181:   PetscInt        i, j;
2182:   IS              dofMap;
2183:   IS              cellIS;
2184:   const PetscInt  ndof = patch->totalDofsPerCell;
2185:   Mat             vecMat;
2186:   PetscInt        cStart, cEnd;
2187:   DM              dm, plex;

2189:   ISGetSize(patch->cells, &ncell);
2190:   if (!ncell) { /* No cells to assemble over -> skip */
2191:     return 0;
2192:   }

2194:   PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);

2196:   PCGetDM(pc, &dm);
2197:   DMConvert(dm, DMPLEX, &plex);
2198:   dm = plex;
2199:   if (!patch->allCells) {
2200:     PetscHSetI    cells;
2201:     PetscHashIter hi;
2202:     PetscInt      pStart, pEnd;
2203:     PetscInt     *allCells = NULL;
2204:     PetscHSetICreate(&cells);
2205:     ISGetIndices(patch->cells, &cellsArray);
2206:     PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);
2207:     for (i = pStart; i < pEnd; i++) {
2208:       PetscSectionGetDof(patch->cellCounts, i, &ncell);
2209:       PetscSectionGetOffset(patch->cellCounts, i, &offset);
2210:       if (ncell <= 0) continue;
2211:       for (j = 0; j < ncell; j++) PetscHSetIAdd(cells, cellsArray[offset + j]);
2212:     }
2213:     ISRestoreIndices(patch->cells, &cellsArray);
2214:     PetscHSetIGetSize(cells, &ncell);
2215:     PetscMalloc1(ncell, &allCells);
2216:     DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2217:     PetscMalloc1(cEnd - cStart, &patch->precomputedTensorLocations);
2218:     i = 0;
2219:     PetscHashIterBegin(cells, hi);
2220:     while (!PetscHashIterAtEnd(cells, hi)) {
2221:       PetscHashIterGetKey(cells, hi, allCells[i]);
2222:       patch->precomputedTensorLocations[allCells[i]] = i;
2223:       PetscHashIterNext(cells, hi);
2224:       i++;
2225:     }
2226:     PetscHSetIDestroy(&cells);
2227:     ISCreateGeneral(PETSC_COMM_SELF, ncell, allCells, PETSC_OWN_POINTER, &patch->allCells);
2228:   }
2229:   ISGetSize(patch->allCells, &ncell);
2230:   if (!patch->cellMats) {
2231:     VecCreateSeq(PETSC_COMM_SELF, ncell * ndof * ndof, &patch->cellMats);
2232:     VecSetBlockSize(patch->cellMats, ndof);
2233:   }
2234:   VecSet(patch->cellMats, 0);

2236:   MatCreateShell(PETSC_COMM_SELF, ncell * ndof, ncell * ndof, ncell * ndof, ncell * ndof, (void *)patch->cellMats, &vecMat);
2237:   MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void (*)(void)) & MatSetValues_PCPatch_Private);
2238:   ISGetSize(patch->allCells, &ncell);
2239:   ISCreateStride(PETSC_COMM_SELF, ndof * ncell, 0, 1, &dofMap);
2240:   ISGetIndices(dofMap, &dofMapArray);
2241:   ISGetIndices(patch->allCells, &cellsArray);
2242:   ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray, PETSC_USE_POINTER, &cellIS);
2243:   /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2244:   PetscCallBack("PCPatch callback", patch->usercomputeop(pc, -1, NULL, vecMat, cellIS, ndof * ncell, dofMapArray, NULL, patch->usercomputeopctx));
2245:   ISDestroy(&cellIS);
2246:   MatDestroy(&vecMat);
2247:   ISRestoreIndices(patch->allCells, &cellsArray);
2248:   ISRestoreIndices(dofMap, &dofMapArray);
2249:   ISDestroy(&dofMap);

2251:   if (patch->usercomputeopintfacet) {
2252:     PetscInt        nIntFacets;
2253:     IS              intFacetsIS;
2254:     const PetscInt *intFacetsArray = NULL;
2255:     if (!patch->allIntFacets) {
2256:       PetscHSetI    facets;
2257:       PetscHashIter hi;
2258:       PetscInt      pStart, pEnd, fStart, fEnd;
2259:       PetscInt     *allIntFacets = NULL;
2260:       PetscHSetICreate(&facets);
2261:       ISGetIndices(patch->intFacets, &intFacetsArray);
2262:       PetscSectionGetChart(patch->intFacetCounts, &pStart, &pEnd);
2263:       DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2264:       for (i = pStart; i < pEnd; i++) {
2265:         PetscSectionGetDof(patch->intFacetCounts, i, &nIntFacets);
2266:         PetscSectionGetOffset(patch->intFacetCounts, i, &offset);
2267:         if (nIntFacets <= 0) continue;
2268:         for (j = 0; j < nIntFacets; j++) PetscHSetIAdd(facets, intFacetsArray[offset + j]);
2269:       }
2270:       ISRestoreIndices(patch->intFacets, &intFacetsArray);
2271:       PetscHSetIGetSize(facets, &nIntFacets);
2272:       PetscMalloc1(nIntFacets, &allIntFacets);
2273:       PetscMalloc1(fEnd - fStart, &patch->precomputedIntFacetTensorLocations);
2274:       i = 0;
2275:       PetscHashIterBegin(facets, hi);
2276:       while (!PetscHashIterAtEnd(facets, hi)) {
2277:         PetscHashIterGetKey(facets, hi, allIntFacets[i]);
2278:         patch->precomputedIntFacetTensorLocations[allIntFacets[i] - fStart] = i;
2279:         PetscHashIterNext(facets, hi);
2280:         i++;
2281:       }
2282:       PetscHSetIDestroy(&facets);
2283:       ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, allIntFacets, PETSC_OWN_POINTER, &patch->allIntFacets);
2284:     }
2285:     ISGetSize(patch->allIntFacets, &nIntFacets);
2286:     if (!patch->intFacetMats) {
2287:       VecCreateSeq(PETSC_COMM_SELF, nIntFacets * ndof * ndof * 4, &patch->intFacetMats);
2288:       VecSetBlockSize(patch->intFacetMats, ndof * 2);
2289:     }
2290:     VecSet(patch->intFacetMats, 0);

2292:     MatCreateShell(PETSC_COMM_SELF, nIntFacets * ndof * 2, nIntFacets * ndof * 2, nIntFacets * ndof * 2, nIntFacets * ndof * 2, (void *)patch->intFacetMats, &vecMat);
2293:     MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void (*)(void)) & MatSetValues_PCPatch_Private);
2294:     ISCreateStride(PETSC_COMM_SELF, 2 * ndof * nIntFacets, 0, 1, &dofMap);
2295:     ISGetIndices(dofMap, &dofMapArray);
2296:     ISGetIndices(patch->allIntFacets, &intFacetsArray);
2297:     ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, intFacetsArray, PETSC_USE_POINTER, &intFacetsIS);
2298:     /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2299:     PetscCallBack("PCPatch callback (interior facets)", patch->usercomputeopintfacet(pc, -1, NULL, vecMat, intFacetsIS, 2 * ndof * nIntFacets, dofMapArray, NULL, patch->usercomputeopintfacetctx));
2300:     ISDestroy(&intFacetsIS);
2301:     MatDestroy(&vecMat);
2302:     ISRestoreIndices(patch->allIntFacets, &intFacetsArray);
2303:     ISRestoreIndices(dofMap, &dofMapArray);
2304:     ISDestroy(&dofMap);
2305:   }
2306:   DMDestroy(&dm);
2307:   PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);

2309:   return 0;
2310: }

2312: PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat, PatchScatterType scattertype)
2313: {
2314:   PC_PATCH          *patch     = (PC_PATCH *)pc->data;
2315:   const PetscScalar *xArray    = NULL;
2316:   PetscScalar       *yArray    = NULL;
2317:   const PetscInt    *gtolArray = NULL;
2318:   PetscInt           dof, offset, lidx;

2321:   VecGetArrayRead(x, &xArray);
2322:   VecGetArray(y, &yArray);
2323:   if (scattertype == SCATTER_WITHARTIFICIAL) {
2324:     PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);
2325:     PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offset);
2326:     ISGetIndices(patch->gtolWithArtificial, &gtolArray);
2327:   } else if (scattertype == SCATTER_WITHALL) {
2328:     PetscSectionGetDof(patch->gtolCountsWithAll, p, &dof);
2329:     PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offset);
2330:     ISGetIndices(patch->gtolWithAll, &gtolArray);
2331:   } else {
2332:     PetscSectionGetDof(patch->gtolCounts, p, &dof);
2333:     PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2334:     ISGetIndices(patch->gtol, &gtolArray);
2335:   }
2338:   for (lidx = 0; lidx < dof; ++lidx) {
2339:     const PetscInt gidx = gtolArray[offset + lidx];

2341:     if (mode == INSERT_VALUES) yArray[lidx] = xArray[gidx]; /* Forward */
2342:     else yArray[gidx] += xArray[lidx];                      /* Reverse */
2343:   }
2344:   if (scattertype == SCATTER_WITHARTIFICIAL) {
2345:     ISRestoreIndices(patch->gtolWithArtificial, &gtolArray);
2346:   } else if (scattertype == SCATTER_WITHALL) {
2347:     ISRestoreIndices(patch->gtolWithAll, &gtolArray);
2348:   } else {
2349:     ISRestoreIndices(patch->gtol, &gtolArray);
2350:   }
2351:   VecRestoreArrayRead(x, &xArray);
2352:   VecRestoreArray(y, &yArray);
2353:   return 0;
2354: }

2356: static PetscErrorCode PCSetUp_PATCH_Linear(PC pc)
2357: {
2358:   PC_PATCH   *patch = (PC_PATCH *)pc->data;
2359:   const char *prefix;
2360:   PetscInt    i;

2362:   if (!pc->setupcalled) {
2364:     if (!patch->denseinverse) {
2365:       PetscMalloc1(patch->npatch, &patch->solver);
2366:       PCGetOptionsPrefix(pc, &prefix);
2367:       for (i = 0; i < patch->npatch; ++i) {
2368:         KSP ksp;
2369:         PC  subpc;

2371:         KSPCreate(PETSC_COMM_SELF, &ksp);
2372:         KSPSetErrorIfNotConverged(ksp, pc->erroriffailure);
2373:         KSPSetOptionsPrefix(ksp, prefix);
2374:         KSPAppendOptionsPrefix(ksp, "sub_");
2375:         PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1);
2376:         KSPGetPC(ksp, &subpc);
2377:         PetscObjectIncrementTabLevel((PetscObject)subpc, (PetscObject)pc, 1);
2378:         patch->solver[i] = (PetscObject)ksp;
2379:       }
2380:     }
2381:   }
2382:   if (patch->save_operators) {
2383:     if (patch->precomputeElementTensors) PCPatchPrecomputePatchTensors_Private(pc);
2384:     for (i = 0; i < patch->npatch; ++i) {
2385:       PCPatchComputeOperator_Internal(pc, NULL, patch->mat[i], i, PETSC_FALSE);
2386:       if (!patch->denseinverse) {
2387:         KSPSetOperators((KSP)patch->solver[i], patch->mat[i], patch->mat[i]);
2388:       } else if (patch->mat[i] && !patch->densesolve) {
2389:         /* Setup matmult callback */
2390:         MatGetOperation(patch->mat[i], MATOP_MULT, (void (**)(void)) & patch->densesolve);
2391:       }
2392:     }
2393:   }
2394:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2395:     for (i = 0; i < patch->npatch; ++i) {
2396:       /* Instead of padding patch->patchUpdate with zeros to get */
2397:       /* patch->patchUpdateWithArtificial and then multiplying with the matrix, */
2398:       /* just get rid of the columns that correspond to the dofs with */
2399:       /* artificial bcs. That's of course fairly inefficient, hopefully we */
2400:       /* can just assemble the rectangular matrix in the first place. */
2401:       Mat      matSquare;
2402:       IS       rowis;
2403:       PetscInt dof;

2405:       MatGetSize(patch->mat[i], &dof, NULL);
2406:       if (dof == 0) {
2407:         patch->matWithArtificial[i] = NULL;
2408:         continue;
2409:       }

2411:       PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);
2412:       PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);

2414:       MatGetSize(matSquare, &dof, NULL);
2415:       ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis);
2416:       if (pc->setupcalled) {
2417:         MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_REUSE_MATRIX, &patch->matWithArtificial[i]);
2418:       } else {
2419:         MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &patch->matWithArtificial[i]);
2420:       }
2421:       ISDestroy(&rowis);
2422:       MatDestroy(&matSquare);
2423:     }
2424:   }
2425:   return 0;
2426: }

2428: static PetscErrorCode PCSetUp_PATCH(PC pc)
2429: {
2430:   PC_PATCH *patch = (PC_PATCH *)pc->data;
2431:   PetscInt  i;
2432:   PetscBool isNonlinear;
2433:   PetscInt  maxDof = -1, maxDofWithArtificial = -1;

2435:   if (!pc->setupcalled) {
2436:     PetscInt pStart, pEnd, p;
2437:     PetscInt localSize;

2439:     PetscLogEventBegin(PC_Patch_CreatePatches, pc, 0, 0, 0);

2441:     isNonlinear = patch->isNonlinear;
2442:     if (!patch->nsubspaces) {
2443:       DM           dm, plex;
2444:       PetscSection s;
2445:       PetscInt     cStart, cEnd, c, Nf, f, numGlobalBcs = 0, *globalBcs, *Nb, **cellDofs;

2447:       PCGetDM(pc, &dm);
2449:       DMConvert(dm, DMPLEX, &plex);
2450:       dm = plex;
2451:       DMGetLocalSection(dm, &s);
2452:       PetscSectionGetNumFields(s, &Nf);
2453:       PetscSectionGetChart(s, &pStart, &pEnd);
2454:       for (p = pStart; p < pEnd; ++p) {
2455:         PetscInt cdof;
2456:         PetscSectionGetConstraintDof(s, p, &cdof);
2457:         numGlobalBcs += cdof;
2458:       }
2459:       DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2460:       PetscMalloc3(Nf, &Nb, Nf, &cellDofs, numGlobalBcs, &globalBcs);
2461:       for (f = 0; f < Nf; ++f) {
2462:         PetscFE        fe;
2463:         PetscDualSpace sp;
2464:         PetscInt       cdoff = 0;

2466:         DMGetField(dm, f, NULL, (PetscObject *)&fe);
2467:         /* PetscFEGetNumComponents(fe, &Nc[f]); */
2468:         PetscFEGetDualSpace(fe, &sp);
2469:         PetscDualSpaceGetDimension(sp, &Nb[f]);

2471:         PetscMalloc1((cEnd - cStart) * Nb[f], &cellDofs[f]);
2472:         for (c = cStart; c < cEnd; ++c) {
2473:           PetscInt *closure = NULL;
2474:           PetscInt  clSize  = 0, cl;

2476:           DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
2477:           for (cl = 0; cl < clSize * 2; cl += 2) {
2478:             const PetscInt p = closure[cl];
2479:             PetscInt       fdof, d, foff;

2481:             PetscSectionGetFieldDof(s, p, f, &fdof);
2482:             PetscSectionGetFieldOffset(s, p, f, &foff);
2483:             for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d;
2484:           }
2485:           DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
2486:         }
2488:       }
2489:       numGlobalBcs = 0;
2490:       for (p = pStart; p < pEnd; ++p) {
2491:         const PetscInt *ind;
2492:         PetscInt        off, cdof, d;

2494:         PetscSectionGetOffset(s, p, &off);
2495:         PetscSectionGetConstraintDof(s, p, &cdof);
2496:         PetscSectionGetConstraintIndices(s, p, &ind);
2497:         for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d];
2498:       }

2500:       PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **)cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs);
2501:       for (f = 0; f < Nf; ++f) PetscFree(cellDofs[f]);
2502:       PetscFree3(Nb, cellDofs, globalBcs);
2503:       PCPatchSetComputeFunction(pc, PCPatchComputeFunction_DMPlex_Private, NULL);
2504:       PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL);
2505:       DMDestroy(&dm);
2506:     }

2508:     localSize = patch->subspaceOffsets[patch->nsubspaces];
2509:     VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localRHS);
2510:     VecSetUp(patch->localRHS);
2511:     VecDuplicate(patch->localRHS, &patch->localUpdate);
2512:     PCPatchCreateCellPatches(pc);
2513:     PCPatchCreateCellPatchDiscretisationInfo(pc);

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

2518:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithArtificial);
2519:     if (isNonlinear) PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithAll);
2520:     for (p = pStart; p < pEnd; ++p) {
2521:       PetscInt dof;

2523:       PetscSectionGetDof(patch->gtolCounts, p, &dof);
2524:       maxDof = PetscMax(maxDof, dof);
2525:       if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2526:         const PetscInt *gtolArray, *gtolArrayWithArtificial = NULL;
2527:         PetscInt        numPatchDofs, offset;
2528:         PetscInt        numPatchDofsWithArtificial, offsetWithArtificial;
2529:         PetscInt        dofWithoutArtificialCounter = 0;
2530:         PetscInt       *patchWithoutArtificialToWithArtificialArray;

2532:         PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);
2533:         maxDofWithArtificial = PetscMax(maxDofWithArtificial, dof);

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

2539:         PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);
2540:         if (numPatchDofs == 0) {
2541:           patch->dofMappingWithoutToWithArtificial[p - pStart] = NULL;
2542:           continue;
2543:         }

2545:         PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2546:         ISGetIndices(patch->gtolWithArtificial, &gtolArrayWithArtificial);
2547:         PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &numPatchDofsWithArtificial);
2548:         PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offsetWithArtificial);

2550:         PetscMalloc1(numPatchDofs, &patchWithoutArtificialToWithArtificialArray);
2551:         for (i = 0; i < numPatchDofsWithArtificial; i++) {
2552:           if (gtolArrayWithArtificial[i + offsetWithArtificial] == gtolArray[offset + dofWithoutArtificialCounter]) {
2553:             patchWithoutArtificialToWithArtificialArray[dofWithoutArtificialCounter] = i;
2554:             dofWithoutArtificialCounter++;
2555:             if (dofWithoutArtificialCounter == numPatchDofs) break;
2556:           }
2557:         }
2558:         ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutArtificialToWithArtificialArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithArtificial[p - pStart]);
2559:         ISRestoreIndices(patch->gtol, &gtolArray);
2560:         ISRestoreIndices(patch->gtolWithArtificial, &gtolArrayWithArtificial);
2561:       }
2562:     }
2563:     for (p = pStart; p < pEnd; ++p) {
2564:       if (isNonlinear) {
2565:         const PetscInt *gtolArray, *gtolArrayWithAll = NULL;
2566:         PetscInt        numPatchDofs, offset;
2567:         PetscInt        numPatchDofsWithAll, offsetWithAll;
2568:         PetscInt        dofWithoutAllCounter = 0;
2569:         PetscInt       *patchWithoutAllToWithAllArray;

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

2575:         PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);
2576:         if (numPatchDofs == 0) {
2577:           patch->dofMappingWithoutToWithAll[p - pStart] = NULL;
2578:           continue;
2579:         }

2581:         PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2582:         ISGetIndices(patch->gtolWithAll, &gtolArrayWithAll);
2583:         PetscSectionGetDof(patch->gtolCountsWithAll, p, &numPatchDofsWithAll);
2584:         PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offsetWithAll);

2586:         PetscMalloc1(numPatchDofs, &patchWithoutAllToWithAllArray);

2588:         for (i = 0; i < numPatchDofsWithAll; i++) {
2589:           if (gtolArrayWithAll[i + offsetWithAll] == gtolArray[offset + dofWithoutAllCounter]) {
2590:             patchWithoutAllToWithAllArray[dofWithoutAllCounter] = i;
2591:             dofWithoutAllCounter++;
2592:             if (dofWithoutAllCounter == numPatchDofs) break;
2593:           }
2594:         }
2595:         ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutAllToWithAllArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithAll[p - pStart]);
2596:         ISRestoreIndices(patch->gtol, &gtolArray);
2597:         ISRestoreIndices(patch->gtolWithAll, &gtolArrayWithAll);
2598:       }
2599:     }
2600:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2601:       VecCreateSeq(PETSC_COMM_SELF, maxDofWithArtificial, &patch->patchRHSWithArtificial);
2602:       VecSetUp(patch->patchRHSWithArtificial);
2603:     }
2604:     VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchRHS);
2605:     VecSetUp(patch->patchRHS);
2606:     VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchUpdate);
2607:     VecSetUp(patch->patchUpdate);
2608:     if (patch->save_operators) {
2609:       PetscMalloc1(patch->npatch, &patch->mat);
2610:       for (i = 0; i < patch->npatch; ++i) PCPatchCreateMatrix_Private(pc, i, &patch->mat[i], PETSC_FALSE);
2611:     }
2612:     PetscLogEventEnd(PC_Patch_CreatePatches, pc, 0, 0, 0);

2614:     /* If desired, calculate weights for dof multiplicity */
2615:     if (patch->partition_of_unity) {
2616:       PetscScalar *input  = NULL;
2617:       PetscScalar *output = NULL;
2618:       Vec          global;

2620:       VecDuplicate(patch->localRHS, &patch->dof_weights);
2621:       if (patch->local_composition_type == PC_COMPOSITE_ADDITIVE) {
2622:         for (i = 0; i < patch->npatch; ++i) {
2623:           PetscInt dof;

2625:           PetscSectionGetDof(patch->gtolCounts, i + pStart, &dof);
2626:           if (dof <= 0) continue;
2627:           VecSet(patch->patchRHS, 1.0);
2628:           PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHS, patch->dof_weights, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);
2629:         }
2630:       } else {
2631:         /* multiplicative is actually only locally multiplicative and globally additive. need the pou where the mesh decomposition overlaps */
2632:         VecSet(patch->dof_weights, 1.0);
2633:       }

2635:       VecDuplicate(patch->dof_weights, &global);
2636:       VecSet(global, 0.);

2638:       VecGetArray(patch->dof_weights, &input);
2639:       VecGetArray(global, &output);
2640:       PetscSFReduceBegin(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM);
2641:       PetscSFReduceEnd(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM);
2642:       VecRestoreArray(patch->dof_weights, &input);
2643:       VecRestoreArray(global, &output);

2645:       VecReciprocal(global);

2647:       VecGetArray(patch->dof_weights, &output);
2648:       VecGetArray(global, &input);
2649:       PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, input, output, MPI_REPLACE);
2650:       PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, input, output, MPI_REPLACE);
2651:       VecRestoreArray(patch->dof_weights, &output);
2652:       VecRestoreArray(global, &input);
2653:       VecDestroy(&global);
2654:     }
2655:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE && patch->save_operators && !patch->isNonlinear) PetscMalloc1(patch->npatch, &patch->matWithArtificial);
2656:   }
2657:   (*patch->setupsolver)(pc);
2658:   return 0;
2659: }

2661: static PetscErrorCode PCApply_PATCH_Linear(PC pc, PetscInt i, Vec x, Vec y)
2662: {
2663:   PC_PATCH *patch = (PC_PATCH *)pc->data;
2664:   KSP       ksp;
2665:   Mat       op;
2666:   PetscInt  m, n;

2668:   if (patch->denseinverse) {
2669:     (*patch->densesolve)(patch->mat[i], x, y);
2670:     return 0;
2671:   }
2672:   ksp = (KSP)patch->solver[i];
2673:   if (!patch->save_operators) {
2674:     Mat mat;

2676:     PCPatchCreateMatrix_Private(pc, i, &mat, PETSC_FALSE);
2677:     /* Populate operator here. */
2678:     PCPatchComputeOperator_Internal(pc, NULL, mat, i, PETSC_FALSE);
2679:     KSPSetOperators(ksp, mat, mat);
2680:     /* Drop reference so the KSPSetOperators below will blow it away. */
2681:     MatDestroy(&mat);
2682:   }
2683:   PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);
2684:   if (!ksp->setfromoptionscalled) KSPSetFromOptions(ksp);
2685:   /* Disgusting trick to reuse work vectors */
2686:   KSPGetOperators(ksp, &op, NULL);
2687:   MatGetLocalSize(op, &m, &n);
2688:   x->map->n = m;
2689:   y->map->n = n;
2690:   x->map->N = m;
2691:   y->map->N = n;
2692:   KSPSolve(ksp, x, y);
2693:   KSPCheckSolve(ksp, pc, y);
2694:   PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);
2695:   if (!patch->save_operators) {
2696:     PC pc;
2697:     KSPSetOperators(ksp, NULL, NULL);
2698:     KSPGetPC(ksp, &pc);
2699:     /* Destroy PC context too, otherwise the factored matrix hangs around. */
2700:     PCReset(pc);
2701:   }
2702:   return 0;
2703: }

2705: static PetscErrorCode PCUpdateMultiplicative_PATCH_Linear(PC pc, PetscInt i, PetscInt pStart)
2706: {
2707:   PC_PATCH *patch = (PC_PATCH *)pc->data;
2708:   Mat       multMat;
2709:   PetscInt  n, m;


2712:   if (patch->save_operators) {
2713:     multMat = patch->matWithArtificial[i];
2714:   } else {
2715:     /*Very inefficient, hopefully we can just assemble the rectangular matrix in the first place.*/
2716:     Mat      matSquare;
2717:     PetscInt dof;
2718:     IS       rowis;
2719:     PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);
2720:     PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);
2721:     MatGetSize(matSquare, &dof, NULL);
2722:     ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis);
2723:     MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &multMat);
2724:     MatDestroy(&matSquare);
2725:     ISDestroy(&rowis);
2726:   }
2727:   /* Disgusting trick to reuse work vectors */
2728:   MatGetLocalSize(multMat, &m, &n);
2729:   patch->patchUpdate->map->n            = n;
2730:   patch->patchRHSWithArtificial->map->n = m;
2731:   patch->patchUpdate->map->N            = n;
2732:   patch->patchRHSWithArtificial->map->N = m;
2733:   MatMult(multMat, patch->patchUpdate, patch->patchRHSWithArtificial);
2734:   VecScale(patch->patchRHSWithArtificial, -1.0);
2735:   PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHSWithArtificial, patch->localRHS, ADD_VALUES, SCATTER_REVERSE, SCATTER_WITHARTIFICIAL);
2736:   if (!patch->save_operators) MatDestroy(&multMat);
2737:   return 0;
2738: }

2740: static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y)
2741: {
2742:   PC_PATCH          *patch        = (PC_PATCH *)pc->data;
2743:   const PetscScalar *globalRHS    = NULL;
2744:   PetscScalar       *localRHS     = NULL;
2745:   PetscScalar       *globalUpdate = NULL;
2746:   const PetscInt    *bcNodes      = NULL;
2747:   PetscInt           nsweep       = patch->symmetrise_sweep ? 2 : 1;
2748:   PetscInt           start[2]     = {0, 0};
2749:   PetscInt           end[2]       = {-1, -1};
2750:   const PetscInt     inc[2]       = {1, -1};
2751:   const PetscScalar *localUpdate;
2752:   const PetscInt    *iterationSet;
2753:   PetscInt           pStart, numBcs, n, sweep, bc, j;

2755:   PetscLogEventBegin(PC_Patch_Apply, pc, 0, 0, 0);
2756:   PetscOptionsPushGetViewerOff(PETSC_TRUE);
2757:   /* start, end, inc have 2 entries to manage a second backward sweep if we symmetrize */
2758:   end[0]   = patch->npatch;
2759:   start[1] = patch->npatch - 1;
2760:   if (patch->user_patches) {
2761:     ISGetLocalSize(patch->iterationSet, &end[0]);
2762:     start[1] = end[0] - 1;
2763:     ISGetIndices(patch->iterationSet, &iterationSet);
2764:   }
2765:   /* Scatter from global space into overlapped local spaces */
2766:   VecGetArrayRead(x, &globalRHS);
2767:   VecGetArray(patch->localRHS, &localRHS);
2768:   PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, globalRHS, localRHS, MPI_REPLACE);
2769:   PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, globalRHS, localRHS, MPI_REPLACE);
2770:   VecRestoreArrayRead(x, &globalRHS);
2771:   VecRestoreArray(patch->localRHS, &localRHS);

2773:   VecSet(patch->localUpdate, 0.0);
2774:   PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);
2775:   PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);
2776:   for (sweep = 0; sweep < nsweep; sweep++) {
2777:     for (j = start[sweep]; j * inc[sweep] < end[sweep] * inc[sweep]; j += inc[sweep]) {
2778:       PetscInt i = patch->user_patches ? iterationSet[j] : j;
2779:       PetscInt start, len;

2781:       PetscSectionGetDof(patch->gtolCounts, i + pStart, &len);
2782:       PetscSectionGetOffset(patch->gtolCounts, i + pStart, &start);
2783:       /* TODO: Squash out these guys in the setup as well. */
2784:       if (len <= 0) continue;
2785:       /* TODO: Do we need different scatters for X and Y? */
2786:       PCPatch_ScatterLocal_Private(pc, i + pStart, patch->localRHS, patch->patchRHS, INSERT_VALUES, SCATTER_FORWARD, SCATTER_INTERIOR);
2787:       (*patch->applysolver)(pc, i, patch->patchRHS, patch->patchUpdate);
2788:       PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchUpdate, patch->localUpdate, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);
2789:       if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) (*patch->updatemultiplicative)(pc, i, pStart);
2790:     }
2791:   }
2792:   PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);
2793:   if (patch->user_patches) ISRestoreIndices(patch->iterationSet, &iterationSet);
2794:   /* XXX: should we do this on the global vector? */
2795:   if (patch->partition_of_unity) VecPointwiseMult(patch->localUpdate, patch->localUpdate, patch->dof_weights);
2796:   /* Now patch->localUpdate contains the solution of the patch solves, so we need to combine them all. */
2797:   VecSet(y, 0.0);
2798:   VecGetArray(y, &globalUpdate);
2799:   VecGetArrayRead(patch->localUpdate, &localUpdate);
2800:   PetscSFReduceBegin(patch->sectionSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);
2801:   PetscSFReduceEnd(patch->sectionSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);
2802:   VecRestoreArrayRead(patch->localUpdate, &localUpdate);

2804:   /* Now we need to send the global BC values through */
2805:   VecGetArrayRead(x, &globalRHS);
2806:   ISGetSize(patch->globalBcNodes, &numBcs);
2807:   ISGetIndices(patch->globalBcNodes, &bcNodes);
2808:   VecGetLocalSize(x, &n);
2809:   for (bc = 0; bc < numBcs; ++bc) {
2810:     const PetscInt idx = bcNodes[bc];
2811:     if (idx < n) globalUpdate[idx] = globalRHS[idx];
2812:   }

2814:   ISRestoreIndices(patch->globalBcNodes, &bcNodes);
2815:   VecRestoreArrayRead(x, &globalRHS);
2816:   VecRestoreArray(y, &globalUpdate);

2818:   PetscOptionsPopGetViewerOff();
2819:   PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0);
2820:   return 0;
2821: }

2823: static PetscErrorCode PCReset_PATCH_Linear(PC pc)
2824: {
2825:   PC_PATCH *patch = (PC_PATCH *)pc->data;
2826:   PetscInt  i;

2828:   if (patch->solver) {
2829:     for (i = 0; i < patch->npatch; ++i) KSPReset((KSP)patch->solver[i]);
2830:   }
2831:   return 0;
2832: }

2834: static PetscErrorCode PCReset_PATCH(PC pc)
2835: {
2836:   PC_PATCH *patch = (PC_PATCH *)pc->data;
2837:   PetscInt  i;


2840:   PetscSFDestroy(&patch->sectionSF);
2841:   PetscSectionDestroy(&patch->cellCounts);
2842:   PetscSectionDestroy(&patch->pointCounts);
2843:   PetscSectionDestroy(&patch->cellNumbering);
2844:   PetscSectionDestroy(&patch->gtolCounts);
2845:   ISDestroy(&patch->gtol);
2846:   ISDestroy(&patch->cells);
2847:   ISDestroy(&patch->points);
2848:   ISDestroy(&patch->dofs);
2849:   ISDestroy(&patch->offs);
2850:   PetscSectionDestroy(&patch->patchSection);
2851:   ISDestroy(&patch->ghostBcNodes);
2852:   ISDestroy(&patch->globalBcNodes);
2853:   PetscSectionDestroy(&patch->gtolCountsWithArtificial);
2854:   ISDestroy(&patch->gtolWithArtificial);
2855:   ISDestroy(&patch->dofsWithArtificial);
2856:   ISDestroy(&patch->offsWithArtificial);
2857:   PetscSectionDestroy(&patch->gtolCountsWithAll);
2858:   ISDestroy(&patch->gtolWithAll);
2859:   ISDestroy(&patch->dofsWithAll);
2860:   ISDestroy(&patch->offsWithAll);
2861:   VecDestroy(&patch->cellMats);
2862:   VecDestroy(&patch->intFacetMats);
2863:   ISDestroy(&patch->allCells);
2864:   ISDestroy(&patch->intFacets);
2865:   ISDestroy(&patch->extFacets);
2866:   ISDestroy(&patch->intFacetsToPatchCell);
2867:   PetscSectionDestroy(&patch->intFacetCounts);
2868:   PetscSectionDestroy(&patch->extFacetCounts);

2870:   if (patch->dofSection)
2871:     for (i = 0; i < patch->nsubspaces; i++) PetscSectionDestroy(&patch->dofSection[i]);
2872:   PetscFree(patch->dofSection);
2873:   PetscFree(patch->bs);
2874:   PetscFree(patch->nodesPerCell);
2875:   if (patch->cellNodeMap)
2876:     for (i = 0; i < patch->nsubspaces; i++) PetscFree(patch->cellNodeMap[i]);
2877:   PetscFree(patch->cellNodeMap);
2878:   PetscFree(patch->subspaceOffsets);

2880:   (*patch->resetsolver)(pc);

2882:   if (patch->subspaces_to_exclude) PetscHSetIDestroy(&patch->subspaces_to_exclude);

2884:   VecDestroy(&patch->localRHS);
2885:   VecDestroy(&patch->localUpdate);
2886:   VecDestroy(&patch->patchRHS);
2887:   VecDestroy(&patch->patchUpdate);
2888:   VecDestroy(&patch->dof_weights);
2889:   if (patch->patch_dof_weights) {
2890:     for (i = 0; i < patch->npatch; ++i) VecDestroy(&patch->patch_dof_weights[i]);
2891:     PetscFree(patch->patch_dof_weights);
2892:   }
2893:   if (patch->mat) {
2894:     for (i = 0; i < patch->npatch; ++i) MatDestroy(&patch->mat[i]);
2895:     PetscFree(patch->mat);
2896:   }
2897:   if (patch->matWithArtificial && !patch->isNonlinear) {
2898:     for (i = 0; i < patch->npatch; ++i) MatDestroy(&patch->matWithArtificial[i]);
2899:     PetscFree(patch->matWithArtificial);
2900:   }
2901:   VecDestroy(&patch->patchRHSWithArtificial);
2902:   if (patch->dofMappingWithoutToWithArtificial) {
2903:     for (i = 0; i < patch->npatch; ++i) ISDestroy(&patch->dofMappingWithoutToWithArtificial[i]);
2904:     PetscFree(patch->dofMappingWithoutToWithArtificial);
2905:   }
2906:   if (patch->dofMappingWithoutToWithAll) {
2907:     for (i = 0; i < patch->npatch; ++i) ISDestroy(&patch->dofMappingWithoutToWithAll[i]);
2908:     PetscFree(patch->dofMappingWithoutToWithAll);
2909:   }
2910:   PetscFree(patch->sub_mat_type);
2911:   if (patch->userIS) {
2912:     for (i = 0; i < patch->npatch; ++i) ISDestroy(&patch->userIS[i]);
2913:     PetscFree(patch->userIS);
2914:   }
2915:   PetscFree(patch->precomputedTensorLocations);
2916:   PetscFree(patch->precomputedIntFacetTensorLocations);

2918:   patch->bs          = NULL;
2919:   patch->cellNodeMap = NULL;
2920:   patch->nsubspaces  = 0;
2921:   ISDestroy(&patch->iterationSet);

2923:   PetscViewerDestroy(&patch->viewerSection);
2924:   return 0;
2925: }

2927: static PetscErrorCode PCDestroy_PATCH_Linear(PC pc)
2928: {
2929:   PC_PATCH *patch = (PC_PATCH *)pc->data;
2930:   PetscInt  i;

2932:   if (patch->solver) {
2933:     for (i = 0; i < patch->npatch; ++i) KSPDestroy((KSP *)&patch->solver[i]);
2934:     PetscFree(patch->solver);
2935:   }
2936:   return 0;
2937: }

2939: static PetscErrorCode PCDestroy_PATCH(PC pc)
2940: {
2941:   PC_PATCH *patch = (PC_PATCH *)pc->data;

2943:   PCReset_PATCH(pc);
2944:   (*patch->destroysolver)(pc);
2945:   PetscFree(pc->data);
2946:   return 0;
2947: }

2949: static PetscErrorCode PCSetFromOptions_PATCH(PC pc, PetscOptionItems *PetscOptionsObject)
2950: {
2951:   PC_PATCH            *patch                 = (PC_PATCH *)pc->data;
2952:   PCPatchConstructType patchConstructionType = PC_PATCH_STAR;
2953:   char                 sub_mat_type[PETSC_MAX_PATH_LEN];
2954:   char                 option[PETSC_MAX_PATH_LEN];
2955:   const char          *prefix;
2956:   PetscBool            flg, dimflg, codimflg;
2957:   MPI_Comm             comm;
2958:   PetscInt            *ifields, nfields, k;
2959:   PCCompositeType      loctype = PC_COMPOSITE_ADDITIVE;

2961:   PetscObjectGetComm((PetscObject)pc, &comm);
2962:   PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix);
2963:   PetscOptionsHeadBegin(PetscOptionsObject, "Patch solver options");

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

2968:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_precompute_element_tensors", patch->classname);
2969:   PetscOptionsBool(option, "Compute each element tensor only once?", "PCPatchSetPrecomputeElementTensors", patch->precomputeElementTensors, &patch->precomputeElementTensors, &flg);
2970:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_partition_of_unity", patch->classname);
2971:   PetscOptionsBool(option, "Weight contributions by dof multiplicity?", "PCPatchSetPartitionOfUnity", patch->partition_of_unity, &patch->partition_of_unity, &flg);

2973:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_local_type", patch->classname);
2974:   PetscOptionsEnum(option, "Type of local solver composition (additive or multiplicative)", "PCPatchSetLocalComposition", PCCompositeTypes, (PetscEnum)loctype, (PetscEnum *)&loctype, &flg);
2975:   if (flg) PCPatchSetLocalComposition(pc, loctype);
2976:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_dense_inverse", patch->classname);
2977:   PetscOptionsBool(option, "Compute inverses of patch matrices and apply directly? Ignores KSP/PC settings on patch.", "PCPatchSetDenseInverse", patch->denseinverse, &patch->denseinverse, &flg);
2978:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_dim", patch->classname);
2979:   PetscOptionsInt(option, "What dimension of mesh point to construct patches by? (0 = vertices)", "PCPATCH", patch->dim, &patch->dim, &dimflg);
2980:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_codim", patch->classname);
2981:   PetscOptionsInt(option, "What co-dimension of mesh point to construct patches by? (0 = cells)", "PCPATCH", patch->codim, &patch->codim, &codimflg);

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

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

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

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

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

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

3004:   /* If the user has set the number of subspaces, use that for the buffer size,
3005:    otherwise use a large number */
3006:   if (patch->nsubspaces <= 0) {
3007:     nfields = 128;
3008:   } else {
3009:     nfields = patch->nsubspaces;
3010:   }
3011:   PetscMalloc1(nfields, &ifields);
3012:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exclude_subspaces", patch->classname);
3013:   PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, option, ifields, &nfields, &flg);
3015:   if (flg) {
3016:     PetscHSetIClear(patch->subspaces_to_exclude);
3017:     for (k = 0; k < nfields; k++) PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]);
3018:   }
3019:   PetscFree(ifields);

3021:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_patches_view", patch->classname);
3022:   PetscOptionsBool(option, "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg);
3023:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_cells_view", patch->classname);
3024:   PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerCells, &patch->formatCells, &patch->viewCells);
3025:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_interior_facets_view", patch->classname);
3026:   PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerIntFacets, &patch->formatIntFacets, &patch->viewIntFacets);
3027:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exterior_facets_view", patch->classname);
3028:   PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerExtFacets, &patch->formatExtFacets, &patch->viewExtFacets);
3029:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_points_view", patch->classname);
3030:   PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints);
3031:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_section_view", patch->classname);
3032:   PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerSection, &patch->formatSection, &patch->viewSection);
3033:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_mat_view", patch->classname);
3034:   PetscOptionsGetViewer(comm, ((PetscObject)pc)->options, prefix, option, &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix);
3035:   PetscOptionsHeadEnd();
3036:   patch->optionsSet = PETSC_TRUE;
3037:   return 0;
3038: }

3040: static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc)
3041: {
3042:   PC_PATCH          *patch = (PC_PATCH *)pc->data;
3043:   KSPConvergedReason reason;
3044:   PetscInt           i;

3046:   if (!patch->save_operators) {
3047:     /* Can't do this here because the sub KSPs don't have an operator attached yet. */
3048:     return 0;
3049:   }
3050:   if (patch->denseinverse) {
3051:     /* No solvers */
3052:     return 0;
3053:   }
3054:   for (i = 0; i < patch->npatch; ++i) {
3055:     if (!((KSP)patch->solver[i])->setfromoptionscalled) KSPSetFromOptions((KSP)patch->solver[i]);
3056:     KSPSetUp((KSP)patch->solver[i]);
3057:     KSPGetConvergedReason((KSP)patch->solver[i], &reason);
3058:     if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
3059:   }
3060:   return 0;
3061: }

3063: static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer)
3064: {
3065:   PC_PATCH   *patch = (PC_PATCH *)pc->data;
3066:   PetscViewer sviewer;
3067:   PetscBool   isascii;
3068:   PetscMPIInt rank;

3070:   /* TODO Redo tabbing with set tbas in new style */
3071:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
3072:   if (!isascii) return 0;
3073:   MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
3074:   PetscViewerASCIIPushTab(viewer);
3075:   PetscViewerASCIIPrintf(viewer, "Subspace Correction preconditioner with %" PetscInt_FMT " patches\n", patch->npatch);
3076:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
3077:     PetscViewerASCIIPrintf(viewer, "Schwarz type: multiplicative\n");
3078:   } else {
3079:     PetscViewerASCIIPrintf(viewer, "Schwarz type: additive\n");
3080:   }
3081:   if (patch->partition_of_unity) PetscViewerASCIIPrintf(viewer, "Weighting by partition of unity\n");
3082:   else PetscViewerASCIIPrintf(viewer, "Not weighting by partition of unity\n");
3083:   if (patch->symmetrise_sweep) PetscViewerASCIIPrintf(viewer, "Symmetrising sweep (start->end, then end->start)\n");
3084:   else PetscViewerASCIIPrintf(viewer, "Not symmetrising sweep\n");
3085:   if (!patch->precomputeElementTensors) PetscViewerASCIIPrintf(viewer, "Not precomputing element tensors (overlapping cells rebuilt in every patch assembly)\n");
3086:   else PetscViewerASCIIPrintf(viewer, "Precomputing element tensors (each cell assembled only once)\n");
3087:   if (!patch->save_operators) PetscViewerASCIIPrintf(viewer, "Not saving patch operators (rebuilt every PCApply)\n");
3088:   else PetscViewerASCIIPrintf(viewer, "Saving patch operators (rebuilt every PCSetUp)\n");
3089:   if (patch->patchconstructop == PCPatchConstruct_Star) PetscViewerASCIIPrintf(viewer, "Patch construction operator: star\n");
3090:   else if (patch->patchconstructop == PCPatchConstruct_Vanka) PetscViewerASCIIPrintf(viewer, "Patch construction operator: Vanka\n");
3091:   else if (patch->patchconstructop == PCPatchConstruct_User) PetscViewerASCIIPrintf(viewer, "Patch construction operator: user-specified\n");
3092:   else PetscViewerASCIIPrintf(viewer, "Patch construction operator: unknown\n");

3094:   if (patch->denseinverse) {
3095:     PetscViewerASCIIPrintf(viewer, "Explicitly forming dense inverse and applying patch solver via MatMult.\n");
3096:   } else {
3097:     if (patch->isNonlinear) {
3098:       PetscViewerASCIIPrintf(viewer, "SNES on patches (all same):\n");
3099:     } else {
3100:       PetscViewerASCIIPrintf(viewer, "KSP on patches (all same):\n");
3101:     }
3102:     if (patch->solver) {
3103:       PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
3104:       if (rank == 0) {
3105:         PetscViewerASCIIPushTab(sviewer);
3106:         PetscObjectView(patch->solver[0], sviewer);
3107:         PetscViewerASCIIPopTab(sviewer);
3108:       }
3109:       PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
3110:     } else {
3111:       PetscViewerASCIIPushTab(viewer);
3112:       PetscViewerASCIIPrintf(viewer, "Solver not yet set.\n");
3113:       PetscViewerASCIIPopTab(viewer);
3114:     }
3115:   }
3116:   PetscViewerASCIIPopTab(viewer);
3117:   return 0;
3118: }

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

3125:    Options Database Keys:
3126: + -pc_patch_cells_view   - Views the process local cell numbers for each patch
3127: . -pc_patch_points_view  - Views the process local mesh point numbers for each patch
3128: . -pc_patch_g2l_view     - Views the map between global dofs and patch local dofs for each patch
3129: . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary
3130: - -pc_patch_sub_mat_view - Views the matrix associated with each patch

3132:    Level: intermediate

3134: .seealso: `PCType`, `PCCreate()`, `PCSetType()`, `PCASM`, `PCJACOBI`, `PCPBJACOBI`, `PCVPBJACOBI`, `SNESPATCH`
3135: M*/
3136: PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc)
3137: {
3138:   PC_PATCH *patch;

3140:   PetscNew(&patch);

3142:   if (patch->subspaces_to_exclude) PetscHSetIDestroy(&patch->subspaces_to_exclude);
3143:   PetscHSetICreate(&patch->subspaces_to_exclude);

3145:   patch->classname   = "pc";
3146:   patch->isNonlinear = PETSC_FALSE;

3148:   /* Set some defaults */
3149:   patch->combined                 = PETSC_FALSE;
3150:   patch->save_operators           = PETSC_TRUE;
3151:   patch->local_composition_type   = PC_COMPOSITE_ADDITIVE;
3152:   patch->precomputeElementTensors = PETSC_FALSE;
3153:   patch->partition_of_unity       = PETSC_FALSE;
3154:   patch->codim                    = -1;
3155:   patch->dim                      = -1;
3156:   patch->vankadim                 = -1;
3157:   patch->ignoredim                = -1;
3158:   patch->pardecomp_overlap        = 0;
3159:   patch->patchconstructop         = PCPatchConstruct_Star;
3160:   patch->symmetrise_sweep         = PETSC_FALSE;
3161:   patch->npatch                   = 0;
3162:   patch->userIS                   = NULL;
3163:   patch->optionsSet               = PETSC_FALSE;
3164:   patch->iterationSet             = NULL;
3165:   patch->user_patches             = PETSC_FALSE;
3166:   PetscStrallocpy(MATDENSE, (char **)&patch->sub_mat_type);
3167:   patch->viewPatches                       = PETSC_FALSE;
3168:   patch->viewCells                         = PETSC_FALSE;
3169:   patch->viewPoints                        = PETSC_FALSE;
3170:   patch->viewSection                       = PETSC_FALSE;
3171:   patch->viewMatrix                        = PETSC_FALSE;
3172:   patch->densesolve                        = NULL;
3173:   patch->setupsolver                       = PCSetUp_PATCH_Linear;
3174:   patch->applysolver                       = PCApply_PATCH_Linear;
3175:   patch->resetsolver                       = PCReset_PATCH_Linear;
3176:   patch->destroysolver                     = PCDestroy_PATCH_Linear;
3177:   patch->updatemultiplicative              = PCUpdateMultiplicative_PATCH_Linear;
3178:   patch->dofMappingWithoutToWithArtificial = NULL;
3179:   patch->dofMappingWithoutToWithAll        = NULL;

3181:   pc->data                 = (void *)patch;
3182:   pc->ops->apply           = PCApply_PATCH;
3183:   pc->ops->applytranspose  = NULL; /* PCApplyTranspose_PATCH; */
3184:   pc->ops->setup           = PCSetUp_PATCH;
3185:   pc->ops->reset           = PCReset_PATCH;
3186:   pc->ops->destroy         = PCDestroy_PATCH;
3187:   pc->ops->setfromoptions  = PCSetFromOptions_PATCH;
3188:   pc->ops->setuponblocks   = PCSetUpOnBlocks_PATCH;
3189:   pc->ops->view            = PCView_PATCH;
3190:   pc->ops->applyrichardson = NULL;

3192:   return 0;
3193: }