Actual source code: pcpatch.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  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: PETSC_STATIC_INLINE PetscErrorCode ObjectView(PetscObject obj, PetscViewer viewer, PetscViewerFormat format)
 13: {

 16:   PetscViewerPushFormat(viewer, format);
 17:   PetscObjectView(obj, viewer);
 18:   PetscViewerPopFormat(viewer);
 19:   return(0);
 20: }

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

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

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

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

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

 71:       /* We've been told to ignore entities of this type.*/
 72:       if (shouldIgnore && newpoint >= iStart && newpoint < iEnd) continue;
 73:       PetscHSetIAdd(ht, newpoint);
 74:     }
 75:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
 76:   }
 77:   DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);
 78:   return(0);
 79: }

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

 95:   PetscHSetIClear(ht);

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

 99:   DMHasLabel(dm, "pyop2_ghost", &isFiredrake);
100:   if (isFiredrake) {
101:     DMGetLabel(dm, "pyop2_ghost", &ghost);
102:     DMLabelCreateIndex(ghost, pStart, pEnd);
103:   } else {
104:     PetscSF sf;
105:     DMGetPointSF(dm, &sf);
106:     PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
107:     nleaves = PetscMax(nleaves, 0);
108:   }

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

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

125:     PetscHSetIGetSize(ht, &htsize);
126:     PetscMalloc1(htsize, &htpoints);
127:     PetscHSetIGetElems(ht, &index, htpoints);

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

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

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

152:   return(0);
153: }

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

166:   PetscHSetIClear(ht);
167:   DMPlexGetChart(dm, &pStart, &pEnd);
168:   ISGetLocalSize(patchis, &n);
169:   ISGetIndices(patchis, &patchdata);
170:   for (i = 0; i < n; ++i) {
171:     const PetscInt ownedpoint = patchdata[i];

173:     if (ownedpoint < pStart || ownedpoint >= pEnd) {
174:       SETERRQ3(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D was not in [%D, %D)", ownedpoint, pStart, pEnd);
175:     }
176:     PetscHSetIAdd(ht, ownedpoint);
177:   }
178:   ISRestoreIndices(patchis, &patchdata);
179:   return(0);
180: }

182: static PetscErrorCode PCPatchCreateDefaultSF_Private(PC pc, PetscInt n, const PetscSF *sf, const PetscInt *bs)
183: {
184:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
185:   PetscInt       i;

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

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

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

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

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

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

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

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

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

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

307: PetscErrorCode PCPatchSetDenseInverse(PC pc, PetscBool flg)
308: {
309:   PC_PATCH *patch = (PC_PATCH *) pc->data;
311:   patch->denseinverse = flg;
312:   return(0);
313: }

315: PetscErrorCode PCPatchGetDenseInverse(PC pc, PetscBool *flg)
316: {
317:   PC_PATCH *patch = (PC_PATCH *) pc->data;
319:   *flg = patch->denseinverse;
320:   return(0);
321: }

323: /* TODO: Docs */
324: PetscErrorCode PCPatchSetIgnoreDim(PC pc, PetscInt dim)
325: {
326:   PC_PATCH *patch = (PC_PATCH *) pc->data;
328:   patch->ignoredim = dim;
329:   return(0);
330: }

332: /* TODO: Docs */
333: PetscErrorCode PCPatchGetIgnoreDim(PC pc, PetscInt *dim)
334: {
335:   PC_PATCH *patch = (PC_PATCH *) pc->data;
337:   *dim = patch->ignoredim;
338:   return(0);
339: }

341: /* TODO: Docs */
342: PetscErrorCode PCPatchSetSaveOperators(PC pc, PetscBool flg)
343: {
344:   PC_PATCH *patch = (PC_PATCH *) pc->data;
346:   patch->save_operators = flg;
347:   return(0);
348: }

350: /* TODO: Docs */
351: PetscErrorCode PCPatchGetSaveOperators(PC pc, PetscBool *flg)
352: {
353:   PC_PATCH *patch = (PC_PATCH *) pc->data;
355:   *flg = patch->save_operators;
356:   return(0);
357: }

359: /* TODO: Docs */
360: PetscErrorCode PCPatchSetPrecomputeElementTensors(PC pc, PetscBool flg)
361: {
362:   PC_PATCH *patch = (PC_PATCH *) pc->data;
364:   patch->precomputeElementTensors = flg;
365:   return(0);
366: }

368: /* TODO: Docs */
369: PetscErrorCode PCPatchGetPrecomputeElementTensors(PC pc, PetscBool *flg)
370: {
371:   PC_PATCH *patch = (PC_PATCH *) pc->data;
373:   *flg = patch->precomputeElementTensors;
374:   return(0);
375: }

377: /* TODO: Docs */
378: PetscErrorCode PCPatchSetPartitionOfUnity(PC pc, PetscBool flg)
379: {
380:   PC_PATCH *patch = (PC_PATCH *) pc->data;
382:   patch->partition_of_unity = flg;
383:   return(0);
384: }

386: /* TODO: Docs */
387: PetscErrorCode PCPatchGetPartitionOfUnity(PC pc, PetscBool *flg)
388: {
389:   PC_PATCH *patch = (PC_PATCH *) pc->data;
391:   *flg = patch->partition_of_unity;
392:   return(0);
393: }

395: /* TODO: Docs */
396: PetscErrorCode PCPatchSetLocalComposition(PC pc, PCCompositeType type)
397: {
398:   PC_PATCH *patch = (PC_PATCH *) pc->data;
400:   if (type != PC_COMPOSITE_ADDITIVE && type != PC_COMPOSITE_MULTIPLICATIVE) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Only supports additive or multiplicative as the local type");
401:   patch->local_composition_type = type;
402:   return(0);
403: }

405: /* TODO: Docs */
406: PetscErrorCode PCPatchGetLocalComposition(PC pc, PCCompositeType *type)
407: {
408:   PC_PATCH *patch = (PC_PATCH *) pc->data;
410:   *type = patch->local_composition_type;
411:   return(0);
412: }

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

421:   if (patch->sub_mat_type) {PetscFree(patch->sub_mat_type);}
422:   PetscStrallocpy(sub_mat_type, (char **) &patch->sub_mat_type);
423:   return(0);
424: }

426: /* TODO: Docs */
427: PetscErrorCode PCPatchGetSubMatType(PC pc, MatType *sub_mat_type)
428: {
429:   PC_PATCH *patch = (PC_PATCH *) pc->data;
431:   *sub_mat_type = patch->sub_mat_type;
432:   return(0);
433: }

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

442:   patch->cellNumbering = cellNumbering;
443:   PetscObjectReference((PetscObject) cellNumbering);
444:   return(0);
445: }

447: /* TODO: Docs */
448: PetscErrorCode PCPatchGetCellNumbering(PC pc, PetscSection *cellNumbering)
449: {
450:   PC_PATCH *patch = (PC_PATCH *) pc->data;
452:   *cellNumbering = patch->cellNumbering;
453:   return(0);
454: }

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

462:   patch->ctype = ctype;
463:   switch (ctype) {
464:   case PC_PATCH_STAR:
465:     patch->user_patches     = PETSC_FALSE;
466:     patch->patchconstructop = PCPatchConstruct_Star;
467:     break;
468:   case PC_PATCH_VANKA:
469:     patch->user_patches     = PETSC_FALSE;
470:     patch->patchconstructop = PCPatchConstruct_Vanka;
471:     break;
472:   case PC_PATCH_PARDECOMP:
473:     patch->user_patches     = PETSC_FALSE;
474:     patch->patchconstructop = PCPatchConstruct_Pardecomp;
475:     break;
476:   case PC_PATCH_USER:
477:   case PC_PATCH_PYTHON:
478:     patch->user_patches     = PETSC_TRUE;
479:     patch->patchconstructop = PCPatchConstruct_User;
480:     if (func) {
481:       patch->userpatchconstructionop = func;
482:       patch->userpatchconstructctx   = ctx;
483:     }
484:     break;
485:   default:
486:     SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_USER, "Unknown patch construction type %D", (PetscInt) patch->ctype);
487:   }
488:   return(0);
489: }

491: /* TODO: Docs */
492: PetscErrorCode PCPatchGetConstructType(PC pc, PCPatchConstructType *ctype, PetscErrorCode (**func)(PC, PetscInt *, IS **, IS *, void *), void **ctx)
493: {
494:   PC_PATCH *patch = (PC_PATCH *) pc->data;

497:   *ctype = patch->ctype;
498:   switch (patch->ctype) {
499:   case PC_PATCH_STAR:
500:   case PC_PATCH_VANKA:
501:   case PC_PATCH_PARDECOMP:
502:     break;
503:   case PC_PATCH_USER:
504:   case PC_PATCH_PYTHON:
505:     *func = patch->userpatchconstructionop;
506:     *ctx  = patch->userpatchconstructctx;
507:     break;
508:   default:
509:     SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_USER, "Unknown patch construction type %D", (PetscInt) patch->ctype);
510:   }
511:   return(0);
512: }

514: /* TODO: Docs */
515: PetscErrorCode PCPatchSetDiscretisationInfo(PC pc, PetscInt nsubspaces, DM *dms, PetscInt *bs, PetscInt *nodesPerCell, const PetscInt **cellNodeMap,
516:                                             const PetscInt *subspaceOffsets, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes)
517: {
518:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
519:   DM             dm, plex;
520:   PetscSF       *sfs;
521:   PetscInt       cStart, cEnd, i, j;

525:   PCGetDM(pc, &dm);
526:   DMConvert(dm, DMPLEX, &plex);
527:   dm = plex;
528:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
529:   PetscMalloc1(nsubspaces, &sfs);
530:   PetscMalloc1(nsubspaces, &patch->dofSection);
531:   PetscMalloc1(nsubspaces, &patch->bs);
532:   PetscMalloc1(nsubspaces, &patch->nodesPerCell);
533:   PetscMalloc1(nsubspaces, &patch->cellNodeMap);
534:   PetscMalloc1(nsubspaces+1, &patch->subspaceOffsets);

536:   patch->nsubspaces       = nsubspaces;
537:   patch->totalDofsPerCell = 0;
538:   for (i = 0; i < nsubspaces; ++i) {
539:     DMGetLocalSection(dms[i], &patch->dofSection[i]);
540:     PetscObjectReference((PetscObject) patch->dofSection[i]);
541:     DMGetSectionSF(dms[i], &sfs[i]);
542:     patch->bs[i]              = bs[i];
543:     patch->nodesPerCell[i]    = nodesPerCell[i];
544:     patch->totalDofsPerCell  += nodesPerCell[i]*bs[i];
545:     PetscMalloc1((cEnd-cStart)*nodesPerCell[i], &patch->cellNodeMap[i]);
546:     for (j = 0; j < (cEnd-cStart)*nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j];
547:     patch->subspaceOffsets[i] = subspaceOffsets[i];
548:   }
549:   PCPatchCreateDefaultSF_Private(pc, nsubspaces, sfs, patch->bs);
550:   PetscFree(sfs);

552:   patch->subspaceOffsets[nsubspaces] = subspaceOffsets[nsubspaces];
553:   ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);
554:   ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);
555:   DMDestroy(&dm);
556:   return(0);
557: }

559: /* TODO: Docs */
560: PetscErrorCode PCPatchSetDiscretisationInfoCombined(PC pc, DM dm, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes)
561: {
562:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
563:   PetscInt       cStart, cEnd, i, j;

567:   patch->combined = PETSC_TRUE;
568:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
569:   DMGetNumFields(dm, &patch->nsubspaces);
570:   PetscCalloc1(patch->nsubspaces, &patch->dofSection);
571:   PetscMalloc1(patch->nsubspaces, &patch->bs);
572:   PetscMalloc1(patch->nsubspaces, &patch->nodesPerCell);
573:   PetscMalloc1(patch->nsubspaces, &patch->cellNodeMap);
574:   PetscCalloc1(patch->nsubspaces+1, &patch->subspaceOffsets);
575:   DMGetLocalSection(dm, &patch->dofSection[0]);
576:   PetscObjectReference((PetscObject) patch->dofSection[0]);
577:   PetscSectionGetStorageSize(patch->dofSection[0], &patch->subspaceOffsets[patch->nsubspaces]);
578:   patch->totalDofsPerCell = 0;
579:   for (i = 0; i < patch->nsubspaces; ++i) {
580:     patch->bs[i]             = 1;
581:     patch->nodesPerCell[i]   = nodesPerCell[i];
582:     patch->totalDofsPerCell += nodesPerCell[i];
583:     PetscMalloc1((cEnd-cStart)*nodesPerCell[i], &patch->cellNodeMap[i]);
584:     for (j = 0; j < (cEnd-cStart)*nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j];
585:   }
586:   DMGetSectionSF(dm, &patch->sectionSF);
587:   PetscObjectReference((PetscObject) patch->sectionSF);
588:   ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);
589:   ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);
590:   return(0);
591: }

593: /*@C

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

597:   Logically collective on PC

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 cellIS,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: .  cellIS           - An array of the cell 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:   Notes:
620:   The entries of F (the output residual vector) have been set to zero before the call.

622: .seealso: PCPatchSetComputeOperator(), PCPatchGetComputeOperator(), PCPatchSetDiscretisationInfo(), PCPatchSetComputeFunctionInteriorFacets()
623: @*/
624: PetscErrorCode PCPatchSetComputeFunction(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;

629:   patch->usercomputef    = func;
630:   patch->usercomputefctx = ctx;
631:   return(0);
632: }

634: /*@C

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

638:   Logically collective on PC

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

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

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

658:   Level: advanced

660:   Notes:
661:   The entries of F (the output residual vector) have been set to zero before the call.

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

670:   patch->usercomputefintfacet    = func;
671:   patch->usercomputefintfacetctx = ctx;
672:   return(0);
673: }

675: /*@C

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

679:   Logically collective on PC

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

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

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

699:   Level: advanced

701:   Notes:
702:   The matrix entries have been set to zero before the call.

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

711:   patch->usercomputeop    = func;
712:   patch->usercomputeopctx = ctx;
713:   return(0);
714: }

716: /*@C

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

720:   Logically collective on PC

722:   Input Parameters:
723: + pc   - The PC
724: . func - The callback
725: - ctx  - The user context

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

730: +  pc               - The PC
731: .  point            - The point
732: .  x                - The input solution (not used in linear problems)
733: .  mat              - The patch matrix
734: .  facetIS          - An array of the facet numbers
735: .  n                - The size of dofsArray
736: .  dofsArray        - The dofmap for the dofs to be solved for
737: .  dofsArrayWithAll - The dofmap for all dofs on the patch
738: -  ctx              - The user context

740:   Level: advanced

742:   Notes:
743:   The matrix entries have been set to zero before the call.

745: .seealso: PCPatchGetComputeOperator(), PCPatchSetComputeFunction(), PCPatchSetDiscretisationInfo()
746: @*/
747: PetscErrorCode PCPatchSetComputeOperatorInteriorFacets(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Mat, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx)
748: {
749:   PC_PATCH *patch = (PC_PATCH *) pc->data;

752:   patch->usercomputeopintfacet    = func;
753:   patch->usercomputeopintfacetctx = ctx;
754:   return(0);
755: }

757: /* On entry, ht contains the topological entities whose dofs we are responsible for solving for;
758:    on exit, cht contains all the topological entities we need to compute their residuals.
759:    In full generality this should incorporate knowledge of the sparsity pattern of the matrix;
760:    here we assume a standard FE sparsity pattern.*/
761: /* TODO: Use DMPlexGetAdjacency() */
762: static PetscErrorCode PCPatchCompleteCellPatch(PC pc, PetscHSetI ht, PetscHSetI cht)
763: {
764:   DM             dm, plex;
765:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
766:   PetscHashIter  hi;
767:   PetscInt       point;
768:   PetscInt      *star = NULL, *closure = NULL;
769:   PetscInt       ignoredim, iStart = 0, iEnd = -1, starSize, closureSize, si, ci;
770:   PetscInt      *fStar = NULL, *fClosure = NULL;
771:   PetscInt       fBegin, fEnd, fsi, fci, fStarSize, fClosureSize;

775:   PCGetDM(pc, &dm);
776:   DMConvert(dm, DMPLEX, &plex);
777:   dm = plex;
778:   DMPlexGetHeightStratum(dm, 1, &fBegin, &fEnd);
779:   PCPatchGetIgnoreDim(pc, &ignoredim);
780:   if (ignoredim >= 0) {DMPlexGetDepthStratum(dm, ignoredim, &iStart, &iEnd);}
781:   PetscHSetIClear(cht);
782:   PetscHashIterBegin(ht, hi);
783:   while (!PetscHashIterAtEnd(ht, hi)) {

785:     PetscHashIterGetKey(ht, hi, point);
786:     PetscHashIterNext(ht, hi);

788:     /* Loop over all the cells that this point connects to */
789:     DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);
790:     for (si = 0; si < starSize*2; si += 2) {
791:       const PetscInt ownedpoint = star[si];
792:       /* TODO Check for point in cht before running through closure again */
793:       /* now loop over all entities in the closure of that cell */
794:       DMPlexGetTransitiveClosure(dm, ownedpoint, PETSC_TRUE, &closureSize, &closure);
795:       for (ci = 0; ci < closureSize*2; ci += 2) {
796:         const PetscInt seenpoint = closure[ci];
797:         if (ignoredim >= 0 && seenpoint >= iStart && seenpoint < iEnd) continue;
798:         PetscHSetIAdd(cht, seenpoint);
799:         /* Facet integrals couple dofs across facets, so in that case for each of
800:          * the facets we need to add all dofs on the other side of the facet to
801:          * the seen dofs. */
802:         if(patch->usercomputeopintfacet){
803:           if(fBegin <= seenpoint && seenpoint < fEnd){
804:             DMPlexGetTransitiveClosure(dm, seenpoint, PETSC_FALSE, &fStarSize, &fStar);
805:             for (fsi = 0; fsi < fStarSize*2; fsi += 2) {
806:               DMPlexGetTransitiveClosure(dm, fStar[fsi], PETSC_TRUE, &fClosureSize, &fClosure);
807:               for (fci = 0; fci < fClosureSize*2; fci += 2) {
808:                 PetscHSetIAdd(cht, fClosure[fci]);
809:               }
810:               DMPlexRestoreTransitiveClosure(dm, fStar[fsi], PETSC_TRUE, NULL, &fClosure);
811:             }
812:             DMPlexRestoreTransitiveClosure(dm, seenpoint, PETSC_FALSE, NULL, &fStar);
813:           }
814:         }
815:       }
816:       DMPlexRestoreTransitiveClosure(dm, ownedpoint, PETSC_TRUE, NULL, &closure);
817:     }
818:     DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, NULL, &star);
819:   }
820:   DMDestroy(&dm);
821:   return(0);
822: }

824: static PetscErrorCode PCPatchGetGlobalDofs(PC pc, PetscSection dofSection[], PetscInt f, PetscBool combined, PetscInt p, PetscInt *dof, PetscInt *off)
825: {

829:   if (combined) {
830:     if (f < 0) {
831:       if (dof) {PetscSectionGetDof(dofSection[0], p, dof);}
832:       if (off) {PetscSectionGetOffset(dofSection[0], p, off);}
833:     } else {
834:       if (dof) {PetscSectionGetFieldDof(dofSection[0], p, f, dof);}
835:       if (off) {PetscSectionGetFieldOffset(dofSection[0], p, f, off);}
836:     }
837:   } else {
838:     if (f < 0) {
839:       PC_PATCH *patch = (PC_PATCH *) pc->data;
840:       PetscInt  fdof, g;

842:       if (dof) {
843:         *dof = 0;
844:         for (g = 0; g < patch->nsubspaces; ++g) {
845:           PetscSectionGetDof(dofSection[g], p, &fdof);
846:           *dof += fdof;
847:         }
848:       }
849:       if (off) {
850:         *off = 0;
851:         for (g = 0; g < patch->nsubspaces; ++g) {
852:           PetscSectionGetOffset(dofSection[g], p, &fdof);
853:           *off += fdof;
854:         }
855:       }
856:     } else {
857:       if (dof) {PetscSectionGetDof(dofSection[f], p, dof);}
858:       if (off) {PetscSectionGetOffset(dofSection[f], p, off);}
859:     }
860:   }
861:   return(0);
862: }

864: /* Given a hash table with a set of topological entities (pts), compute the degrees of
865:    freedom in global concatenated numbering on those entities.
866:    For Vanka smoothing, this needs to do something special: ignore dofs of the
867:    constraint subspace on entities that aren't the base entity we're building the patch
868:    around. */
869: static PetscErrorCode PCPatchGetPointDofs(PC pc, PetscHSetI pts, PetscHSetI dofs, PetscInt base, PetscHSetI* subspaces_to_exclude)
870: {
871:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
872:   PetscHashIter  hi;
873:   PetscInt       ldof, loff;
874:   PetscInt       k, p;

878:   PetscHSetIClear(dofs);
879:   for (k = 0; k < patch->nsubspaces; ++k) {
880:     PetscInt subspaceOffset = patch->subspaceOffsets[k];
881:     PetscInt bs             = patch->bs[k];
882:     PetscInt j, l;

884:     if (subspaces_to_exclude != NULL) {
885:       PetscBool should_exclude_k = PETSC_FALSE;
886:       PetscHSetIHas(*subspaces_to_exclude, k, &should_exclude_k);
887:       if (should_exclude_k) {
888:         /* only get this subspace dofs at the base entity, not any others */
889:         PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, base, &ldof, &loff);
890:         if (0 == ldof) continue;
891:         for (j = loff; j < ldof + loff; ++j) {
892:           for (l = 0; l < bs; ++l) {
893:             PetscInt dof = bs*j + l + subspaceOffset;
894:             PetscHSetIAdd(dofs, dof);
895:           }
896:         }
897:         continue; /* skip the other dofs of this subspace */
898:       }
899:     }

901:     PetscHashIterBegin(pts, hi);
902:     while (!PetscHashIterAtEnd(pts, hi)) {
903:       PetscHashIterGetKey(pts, hi, p);
904:       PetscHashIterNext(pts, hi);
905:       PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, p, &ldof, &loff);
906:       if (0 == ldof) continue;
907:       for (j = loff; j < ldof + loff; ++j) {
908:         for (l = 0; l < bs; ++l) {
909:           PetscInt dof = bs*j + l + subspaceOffset;
910:           PetscHSetIAdd(dofs, dof);
911:         }
912:       }
913:     }
914:   }
915:   return(0);
916: }

918: /* Given two hash tables A and B, compute the keys in B that are not in A, and put them in C */
919: static PetscErrorCode PCPatchComputeSetDifference_Private(PetscHSetI A, PetscHSetI B, PetscHSetI C)
920: {
921:   PetscHashIter  hi;
922:   PetscInt       key;
923:   PetscBool      flg;

927:   PetscHSetIClear(C);
928:   PetscHashIterBegin(B, hi);
929:   while (!PetscHashIterAtEnd(B, hi)) {
930:     PetscHashIterGetKey(B, hi, key);
931:     PetscHashIterNext(B, hi);
932:     PetscHSetIHas(A, key, &flg);
933:     if (!flg) {PetscHSetIAdd(C, key);}
934:   }
935:   return(0);
936: }

938: /*
939:  * PCPatchCreateCellPatches - create patches.
940:  *
941:  * Input Parameters:
942:  * + dm - The DMPlex object defining the mesh
943:  *
944:  * Output Parameters:
945:  * + cellCounts  - Section with counts of cells around each vertex
946:  * . cells       - IS of the cell point indices of cells in each patch
947:  * . pointCounts - Section with counts of cells around each vertex
948:  * - point       - IS of the cell point indices of cells in each patch
949:  */
950: static PetscErrorCode PCPatchCreateCellPatches(PC pc)
951: {
952:   PC_PATCH       *patch = (PC_PATCH *) pc->data;
953:   DMLabel         ghost = NULL;
954:   DM              dm, plex;
955:   PetscHSetI      ht, cht;
956:   PetscSection    cellCounts,  pointCounts, intFacetCounts, extFacetCounts;
957:   PetscInt       *cellsArray, *pointsArray, *intFacetsArray, *extFacetsArray, *intFacetsToPatchCell;
958:   PetscInt        numCells, numPoints, numIntFacets, numExtFacets;
959:   const PetscInt *leaves;
960:   PetscInt        nleaves, pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, v;
961:   PetscBool       isFiredrake;
962:   PetscErrorCode  ierr;

965:   /* Used to keep track of the cells in the patch. */
966:   PetscHSetICreate(&ht);
967:   PetscHSetICreate(&cht);

969:   PCGetDM(pc, &dm);
970:   if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONGSTATE, "DM not yet set on patch PC\n");
971:   DMConvert(dm, DMPLEX, &plex);
972:   dm = plex;
973:   DMPlexGetChart(dm, &pStart, &pEnd);
974:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);

976:   if (patch->user_patches) {
977:     patch->userpatchconstructionop(pc, &patch->npatch, &patch->userIS, &patch->iterationSet, patch->userpatchconstructctx);
978:     vStart = 0; vEnd = patch->npatch;
979:   } else if (patch->ctype == PC_PATCH_PARDECOMP) {
980:     vStart = 0; vEnd = 1;
981:   } else if (patch->codim < 0) {
982:     if (patch->dim < 0) {DMPlexGetDepthStratum(dm,  0,            &vStart, &vEnd);}
983:     else                {DMPlexGetDepthStratum(dm,  patch->dim,   &vStart, &vEnd);}
984:   } else                {DMPlexGetHeightStratum(dm, patch->codim, &vStart, &vEnd);}
985:   patch->npatch = vEnd - vStart;

987:   /* These labels mark the owned points.  We only create patches around points that this process owns. */
988:   DMHasLabel(dm, "pyop2_ghost", &isFiredrake);
989:   if (isFiredrake) {
990:     DMGetLabel(dm, "pyop2_ghost", &ghost);
991:     DMLabelCreateIndex(ghost, pStart, pEnd);
992:   } else {
993:     PetscSF sf;

995:     DMGetPointSF(dm, &sf);
996:     PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
997:     nleaves = PetscMax(nleaves, 0);
998:   }

1000:   PetscSectionCreate(PETSC_COMM_SELF, &patch->cellCounts);
1001:   PetscObjectSetName((PetscObject) patch->cellCounts, "Patch Cell Layout");
1002:   cellCounts = patch->cellCounts;
1003:   PetscSectionSetChart(cellCounts, vStart, vEnd);
1004:   PetscSectionCreate(PETSC_COMM_SELF, &patch->pointCounts);
1005:   PetscObjectSetName((PetscObject) patch->pointCounts, "Patch Point Layout");
1006:   pointCounts = patch->pointCounts;
1007:   PetscSectionSetChart(pointCounts, vStart, vEnd);
1008:   PetscSectionCreate(PETSC_COMM_SELF, &patch->extFacetCounts);
1009:   PetscObjectSetName((PetscObject) patch->extFacetCounts, "Patch Exterior Facet Layout");
1010:   extFacetCounts = patch->extFacetCounts;
1011:   PetscSectionSetChart(extFacetCounts, vStart, vEnd);
1012:   PetscSectionCreate(PETSC_COMM_SELF, &patch->intFacetCounts);
1013:   PetscObjectSetName((PetscObject) patch->intFacetCounts, "Patch Interior Facet Layout");
1014:   intFacetCounts = patch->intFacetCounts;
1015:   PetscSectionSetChart(intFacetCounts, vStart, vEnd);
1016:   /* Count cells and points in the patch surrounding each entity */
1017:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1018:   for (v = vStart; v < vEnd; ++v) {
1019:     PetscHashIter hi;
1020:     PetscInt       chtSize, loc = -1;
1021:     PetscBool      flg;

1023:     if (!patch->user_patches && patch->ctype != PC_PATCH_PARDECOMP) {
1024:       if (ghost) {DMLabelHasPoint(ghost, v, &flg);}
1025:       else       {PetscFindInt(v, nleaves, leaves, &loc); flg = loc >=0 ? PETSC_TRUE : PETSC_FALSE;}
1026:       /* Not an owned entity, don't make a cell patch. */
1027:       if (flg) continue;
1028:     }

1030:     patch->patchconstructop((void *) patch, dm, v, ht);
1031:     PCPatchCompleteCellPatch(pc, ht, cht);
1032:     PetscHSetIGetSize(cht, &chtSize);
1033:     /* empty patch, continue */
1034:     if (chtSize == 0) continue;

1036:     /* safe because size(cht) > 0 from above */
1037:     PetscHashIterBegin(cht, hi);
1038:     while (!PetscHashIterAtEnd(cht, hi)) {
1039:       PetscInt point, pdof;

1041:       PetscHashIterGetKey(cht, hi, point);
1042:       if (fStart <= point && point < fEnd) {
1043:         const PetscInt *support;
1044:         PetscInt supportSize, p;
1045:         PetscBool interior = PETSC_TRUE;
1046:         DMPlexGetSupport(dm, point, &support);
1047:         DMPlexGetSupportSize(dm, point, &supportSize);
1048:         if (supportSize == 1) {
1049:           interior = PETSC_FALSE;
1050:         } else {
1051:           for (p = 0; p < supportSize; p++) {
1052:             PetscBool found;
1053:             /* FIXME: can I do this while iterating over cht? */
1054:             PetscHSetIHas(cht, support[p], &found);
1055:             if (!found) {
1056:               interior = PETSC_FALSE;
1057:               break;
1058:             }
1059:           }
1060:         }
1061:         if (interior) {
1062:           PetscSectionAddDof(intFacetCounts, v, 1);
1063:         } else {
1064:           PetscSectionAddDof(extFacetCounts, v, 1);
1065:         }
1066:       }
1067:       PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);
1068:       if (pdof)                            {PetscSectionAddDof(pointCounts, v, 1);}
1069:       if (point >= cStart && point < cEnd) {PetscSectionAddDof(cellCounts, v, 1);}
1070:       PetscHashIterNext(cht, hi);
1071:     }
1072:   }
1073:   if (isFiredrake) {DMLabelDestroyIndex(ghost);}

1075:   PetscSectionSetUp(cellCounts);
1076:   PetscSectionGetStorageSize(cellCounts, &numCells);
1077:   PetscMalloc1(numCells, &cellsArray);
1078:   PetscSectionSetUp(pointCounts);
1079:   PetscSectionGetStorageSize(pointCounts, &numPoints);
1080:   PetscMalloc1(numPoints, &pointsArray);

1082:   PetscSectionSetUp(intFacetCounts);
1083:   PetscSectionSetUp(extFacetCounts);
1084:   PetscSectionGetStorageSize(intFacetCounts, &numIntFacets);
1085:   PetscSectionGetStorageSize(extFacetCounts, &numExtFacets);
1086:   PetscMalloc1(numIntFacets, &intFacetsArray);
1087:   PetscMalloc1(numIntFacets*2, &intFacetsToPatchCell);
1088:   PetscMalloc1(numExtFacets, &extFacetsArray);


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

1096:     PetscSectionGetDof(pointCounts, v, &dof);
1097:     PetscSectionGetOffset(pointCounts, v, &off);
1098:     PetscSectionGetDof(cellCounts, v, &cdof);
1099:     PetscSectionGetOffset(cellCounts, v, &coff);
1100:     PetscSectionGetDof(intFacetCounts, v, &ifdof);
1101:     PetscSectionGetOffset(intFacetCounts, v, &ifoff);
1102:     PetscSectionGetDof(extFacetCounts, v, &efdof);
1103:     PetscSectionGetOffset(extFacetCounts, v, &efoff);
1104:     if (dof <= 0) continue;
1105:     patch->patchconstructop((void *) patch, dm, v, ht);
1106:     PCPatchCompleteCellPatch(pc, ht, cht);
1107:     PetscHashIterBegin(cht, hi);
1108:     while (!PetscHashIterAtEnd(cht, hi)) {
1109:       PetscInt point;

1111:       PetscHashIterGetKey(cht, hi, point);
1112:       if (fStart <= point && point < fEnd) {
1113:         const PetscInt *support;
1114:         PetscInt supportSize, p;
1115:         PetscBool interior = PETSC_TRUE;
1116:         DMPlexGetSupport(dm, point, &support);
1117:         DMPlexGetSupportSize(dm, point, &supportSize);
1118:         if (supportSize == 1) {
1119:           interior = PETSC_FALSE;
1120:         } else {
1121:           for (p = 0; p < supportSize; p++) {
1122:             PetscBool found;
1123:             /* FIXME: can I do this while iterating over cht? */
1124:             PetscHSetIHas(cht, support[p], &found);
1125:             if (!found) {
1126:               interior = PETSC_FALSE;
1127:               break;
1128:             }
1129:           }
1130:         }
1131:         if (interior) {
1132:           intFacetsToPatchCell[2*(ifoff + ifn)] = support[0];
1133:           intFacetsToPatchCell[2*(ifoff + ifn) + 1] = support[1];
1134:           intFacetsArray[ifoff + ifn++] = point;
1135:         } else {
1136:           extFacetsArray[efoff + efn++] = point;
1137:         }
1138:       }
1139:       PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);
1140:       if (pdof)                            {pointsArray[off + n++] = point;}
1141:       if (point >= cStart && point < cEnd) {cellsArray[coff + cn++] = point;}
1142:       PetscHashIterNext(cht, hi);
1143:     }
1144:     if (ifn != ifdof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of interior facets in patch %D is %D, but should be %D", v, ifn, ifdof);
1145:     if (efn != efdof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of exterior facets in patch %D is %D, but should be %D", v, efn, efdof);
1146:     if (cn != cdof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of cells in patch %D is %D, but should be %D", v, cn, cdof);
1147:     if (n  != dof)  SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of points in patch %D is %D, but should be %D", v, n, dof);

1149:     for (ifn = 0; ifn < ifdof; ifn++) {
1150:       PetscInt cell0 = intFacetsToPatchCell[2*(ifoff + ifn)];
1151:       PetscInt cell1 = intFacetsToPatchCell[2*(ifoff + ifn) + 1];
1152:       PetscBool found0 = PETSC_FALSE, found1 = PETSC_FALSE;
1153:       for (n = 0; n < cdof; n++) {
1154:         if (!found0 && cell0 == cellsArray[coff + n]) {
1155:           intFacetsToPatchCell[2*(ifoff + ifn)] = n;
1156:           found0 = PETSC_TRUE;
1157:         }
1158:         if (!found1 && cell1 == cellsArray[coff + n]) {
1159:           intFacetsToPatchCell[2*(ifoff + ifn) + 1] = n;
1160:           found1 = PETSC_TRUE;
1161:         }
1162:         if (found0 && found1) break;
1163:       }
1164:       if (!(found0 && found1)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Didn't manage to find local point numbers for facet support");
1165:     }
1166:   }
1167:   PetscHSetIDestroy(&ht);
1168:   PetscHSetIDestroy(&cht);

1170:   ISCreateGeneral(PETSC_COMM_SELF, numCells,  cellsArray,  PETSC_OWN_POINTER, &patch->cells);
1171:   PetscObjectSetName((PetscObject) patch->cells,  "Patch Cells");
1172:   if (patch->viewCells) {
1173:     ObjectView((PetscObject) patch->cellCounts, patch->viewerCells, patch->formatCells);
1174:     ObjectView((PetscObject) patch->cells,      patch->viewerCells, patch->formatCells);
1175:   }
1176:   ISCreateGeneral(PETSC_COMM_SELF, numIntFacets,  intFacetsArray,  PETSC_OWN_POINTER, &patch->intFacets);
1177:   PetscObjectSetName((PetscObject) patch->intFacets,  "Patch Interior Facets");
1178:   ISCreateGeneral(PETSC_COMM_SELF, 2*numIntFacets, intFacetsToPatchCell, PETSC_OWN_POINTER, &patch->intFacetsToPatchCell);
1179:   PetscObjectSetName((PetscObject) patch->intFacetsToPatchCell,  "Patch Interior Facets local support");
1180:   if (patch->viewIntFacets) {
1181:     ObjectView((PetscObject) patch->intFacetCounts,       patch->viewerIntFacets, patch->formatIntFacets);
1182:     ObjectView((PetscObject) patch->intFacets,            patch->viewerIntFacets, patch->formatIntFacets);
1183:     ObjectView((PetscObject) patch->intFacetsToPatchCell, patch->viewerIntFacets, patch->formatIntFacets);
1184:   }
1185:   ISCreateGeneral(PETSC_COMM_SELF, numExtFacets,  extFacetsArray,  PETSC_OWN_POINTER, &patch->extFacets);
1186:   PetscObjectSetName((PetscObject) patch->extFacets,  "Patch Exterior Facets");
1187:   if (patch->viewExtFacets) {
1188:     ObjectView((PetscObject) patch->extFacetCounts, patch->viewerExtFacets, patch->formatExtFacets);
1189:     ObjectView((PetscObject) patch->extFacets,      patch->viewerExtFacets, patch->formatExtFacets);
1190:   }
1191:   ISCreateGeneral(PETSC_COMM_SELF, numPoints, pointsArray, PETSC_OWN_POINTER, &patch->points);
1192:   PetscObjectSetName((PetscObject) patch->points, "Patch Points");
1193:   if (patch->viewPoints) {
1194:     ObjectView((PetscObject) patch->pointCounts, patch->viewerPoints, patch->formatPoints);
1195:     ObjectView((PetscObject) patch->points,      patch->viewerPoints, patch->formatPoints);
1196:   }
1197:   DMDestroy(&dm);
1198:   return(0);
1199: }

1201: /*
1202:  * PCPatchCreateCellPatchDiscretisationInfo - Build the dof maps for cell patches
1203:  *
1204:  * Input Parameters:
1205:  * + dm - The DMPlex object defining the mesh
1206:  * . cellCounts - Section with counts of cells around each vertex
1207:  * . cells - IS of the cell point indices of cells in each patch
1208:  * . cellNumbering - Section mapping plex cell points to Firedrake cell indices.
1209:  * . nodesPerCell - number of nodes per cell.
1210:  * - cellNodeMap - map from cells to node indices (nodesPerCell * numCells)
1211:  *
1212:  * Output Parameters:
1213:  * + dofs - IS of local dof numbers of each cell in the patch, where local is a patch local numbering
1214:  * . gtolCounts - Section with counts of dofs per cell patch
1215:  * - gtol - IS mapping from global dofs to local dofs for each patch.
1216:  */
1217: static PetscErrorCode PCPatchCreateCellPatchDiscretisationInfo(PC pc)
1218: {
1219:   PC_PATCH       *patch           = (PC_PATCH *) pc->data;
1220:   PetscSection    cellCounts      = patch->cellCounts;
1221:   PetscSection    pointCounts     = patch->pointCounts;
1222:   PetscSection    gtolCounts, gtolCountsWithArtificial = NULL, gtolCountsWithAll = NULL;
1223:   IS              cells           = patch->cells;
1224:   IS              points          = patch->points;
1225:   PetscSection    cellNumbering   = patch->cellNumbering;
1226:   PetscInt        Nf              = patch->nsubspaces;
1227:   PetscInt        numCells, numPoints;
1228:   PetscInt        numDofs;
1229:   PetscInt        numGlobalDofs, numGlobalDofsWithArtificial, numGlobalDofsWithAll;
1230:   PetscInt        totalDofsPerCell = patch->totalDofsPerCell;
1231:   PetscInt        vStart, vEnd, v;
1232:   const PetscInt *cellsArray, *pointsArray;
1233:   PetscInt       *newCellsArray   = NULL;
1234:   PetscInt       *dofsArray       = NULL;
1235:   PetscInt       *dofsArrayWithArtificial = NULL;
1236:   PetscInt       *dofsArrayWithAll = NULL;
1237:   PetscInt       *offsArray       = NULL;
1238:   PetscInt       *offsArrayWithArtificial = NULL;
1239:   PetscInt       *offsArrayWithAll = NULL;
1240:   PetscInt       *asmArray        = NULL;
1241:   PetscInt       *asmArrayWithArtificial = NULL;
1242:   PetscInt       *asmArrayWithAll = NULL;
1243:   PetscInt       *globalDofsArray = NULL;
1244:   PetscInt       *globalDofsArrayWithArtificial = NULL;
1245:   PetscInt       *globalDofsArrayWithAll = NULL;
1246:   PetscInt        globalIndex     = 0;
1247:   PetscInt        key             = 0;
1248:   PetscInt        asmKey          = 0;
1249:   DM              dm              = NULL, plex;
1250:   const PetscInt *bcNodes         = NULL;
1251:   PetscHMapI      ht;
1252:   PetscHMapI      htWithArtificial;
1253:   PetscHMapI      htWithAll;
1254:   PetscHSetI      globalBcs;
1255:   PetscInt        numBcs;
1256:   PetscHSetI      ownedpts, seenpts, owneddofs, seendofs, artificialbcs;
1257:   PetscInt        pStart, pEnd, p, i;
1258:   char            option[PETSC_MAX_PATH_LEN];
1259:   PetscBool       isNonlinear;
1260:   PetscErrorCode  ierr;


1264:   PCGetDM(pc, &dm); 
1265:   DMConvert(dm, DMPLEX, &plex);
1266:   dm = plex;
1267:   /* dofcounts section is cellcounts section * dofPerCell */
1268:   PetscSectionGetStorageSize(cellCounts, &numCells);
1269:   PetscSectionGetStorageSize(patch->pointCounts, &numPoints);
1270:   numDofs = numCells * totalDofsPerCell;
1271:   PetscMalloc1(numDofs, &dofsArray);
1272:   PetscMalloc1(numPoints*Nf, &offsArray);
1273:   PetscMalloc1(numDofs, &asmArray);
1274:   PetscMalloc1(numCells, &newCellsArray);
1275:   PetscSectionGetChart(cellCounts, &vStart, &vEnd);
1276:   PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCounts);
1277:   gtolCounts = patch->gtolCounts;
1278:   PetscSectionSetChart(gtolCounts, vStart, vEnd);
1279:   PetscObjectSetName((PetscObject) patch->gtolCounts, "Patch Global Index Section");

1281:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1282:     PetscMalloc1(numPoints*Nf, &offsArrayWithArtificial);
1283:     PetscMalloc1(numDofs, &asmArrayWithArtificial);
1284:     PetscMalloc1(numDofs, &dofsArrayWithArtificial);
1285:     PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithArtificial);
1286:     gtolCountsWithArtificial = patch->gtolCountsWithArtificial;
1287:     PetscSectionSetChart(gtolCountsWithArtificial, vStart, vEnd);
1288:     PetscObjectSetName((PetscObject) patch->gtolCountsWithArtificial, "Patch Global Index Section Including Artificial BCs");
1289:   }

1291:   isNonlinear = patch->isNonlinear;
1292:   if (isNonlinear) {
1293:     PetscMalloc1(numPoints*Nf, &offsArrayWithAll);
1294:     PetscMalloc1(numDofs, &asmArrayWithAll);
1295:     PetscMalloc1(numDofs, &dofsArrayWithAll);
1296:     PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithAll);
1297:     gtolCountsWithAll = patch->gtolCountsWithAll;
1298:     PetscSectionSetChart(gtolCountsWithAll, vStart, vEnd);
1299:     PetscObjectSetName((PetscObject) patch->gtolCountsWithAll, "Patch Global Index Section Including All BCs");
1300:   }

1302:   /* Outside the patch loop, get the dofs that are globally-enforced Dirichlet
1303:    conditions */
1304:   PetscHSetICreate(&globalBcs);
1305:   ISGetIndices(patch->ghostBcNodes, &bcNodes); 
1306:   ISGetSize(patch->ghostBcNodes, &numBcs); 
1307:   for (i = 0; i < numBcs; ++i) {
1308:     PetscHSetIAdd(globalBcs, bcNodes[i]); /* these are already in concatenated numbering */
1309:   }
1310:   ISRestoreIndices(patch->ghostBcNodes, &bcNodes); 
1311:   ISDestroy(&patch->ghostBcNodes);  /* memory optimisation */

1313:   /* Hash tables for artificial BC construction */
1314:   PetscHSetICreate(&ownedpts);
1315:   PetscHSetICreate(&seenpts);
1316:   PetscHSetICreate(&owneddofs);
1317:   PetscHSetICreate(&seendofs);
1318:   PetscHSetICreate(&artificialbcs);

1320:   ISGetIndices(cells, &cellsArray);
1321:   ISGetIndices(points, &pointsArray);
1322:   PetscHMapICreate(&ht);
1323:   PetscHMapICreate(&htWithArtificial);
1324:   PetscHMapICreate(&htWithAll);
1325:   for (v = vStart; v < vEnd; ++v) {
1326:     PetscInt localIndex = 0;
1327:     PetscInt localIndexWithArtificial = 0;
1328:     PetscInt localIndexWithAll = 0;
1329:     PetscInt dof, off, i, j, k, l;

1331:     PetscHMapIClear(ht);
1332:     PetscHMapIClear(htWithArtificial);
1333:     PetscHMapIClear(htWithAll);
1334:     PetscSectionGetDof(cellCounts, v, &dof);
1335:     PetscSectionGetOffset(cellCounts, v, &off);
1336:     if (dof <= 0) continue;

1338:     /* Calculate the global numbers of the artificial BC dofs here first */
1339:     patch->patchconstructop((void*)patch, dm, v, ownedpts); 
1340:     PCPatchCompleteCellPatch(pc, ownedpts, seenpts); 
1341:     PCPatchGetPointDofs(pc, ownedpts, owneddofs, v, &patch->subspaces_to_exclude); 
1342:     PCPatchGetPointDofs(pc, seenpts, seendofs, v, NULL); 
1343:     PCPatchComputeSetDifference_Private(owneddofs, seendofs, artificialbcs); 
1344:     if (patch->viewPatches) {
1345:       PetscHSetI globalbcdofs;
1346:       PetscHashIter hi;
1347:       MPI_Comm comm = PetscObjectComm((PetscObject)pc);

1349:       PetscHSetICreate(&globalbcdofs);
1350:       PetscSynchronizedPrintf(comm, "Patch %d: owned dofs:\n", v); 
1351:       PetscHashIterBegin(owneddofs, hi);
1352:       while (!PetscHashIterAtEnd(owneddofs, hi)) {
1353:         PetscInt globalDof;

1355:         PetscHashIterGetKey(owneddofs, hi, globalDof);
1356:         PetscHashIterNext(owneddofs, hi);
1357:         PetscSynchronizedPrintf(comm, "%d ", globalDof); 
1358:       }
1359:       PetscSynchronizedPrintf(comm, "\n"); 
1360:       PetscSynchronizedPrintf(comm, "Patch %d: seen dofs:\n", v); 
1361:       PetscHashIterBegin(seendofs, hi);
1362:       while (!PetscHashIterAtEnd(seendofs, hi)) {
1363:         PetscInt globalDof;
1364:         PetscBool flg;

1366:         PetscHashIterGetKey(seendofs, hi, globalDof);
1367:         PetscHashIterNext(seendofs, hi);
1368:         PetscSynchronizedPrintf(comm, "%d ", globalDof); 

1370:         PetscHSetIHas(globalBcs, globalDof, &flg);
1371:         if (flg) {PetscHSetIAdd(globalbcdofs, globalDof);}
1372:       }
1373:       PetscSynchronizedPrintf(comm, "\n"); 
1374:       PetscSynchronizedPrintf(comm, "Patch %d: global BCs:\n", v); 
1375:       PetscHSetIGetSize(globalbcdofs, &numBcs);
1376:       if (numBcs > 0) {
1377:         PetscHashIterBegin(globalbcdofs, hi);
1378:         while (!PetscHashIterAtEnd(globalbcdofs, hi)) {
1379:           PetscInt globalDof;
1380:           PetscHashIterGetKey(globalbcdofs, hi, globalDof);
1381:           PetscHashIterNext(globalbcdofs, hi);
1382:           PetscSynchronizedPrintf(comm, "%d ", globalDof);
1383:         }
1384:       }
1385:       PetscSynchronizedPrintf(comm, "\n");
1386:       PetscSynchronizedPrintf(comm, "Patch %d: artificial BCs:\n", v);
1387:       PetscHSetIGetSize(artificialbcs, &numBcs);
1388:       if (numBcs > 0) {
1389:         PetscHashIterBegin(artificialbcs, hi);
1390:         while (!PetscHashIterAtEnd(artificialbcs, hi)) {
1391:           PetscInt globalDof;
1392:           PetscHashIterGetKey(artificialbcs, hi, globalDof);
1393:           PetscHashIterNext(artificialbcs, hi);
1394:           PetscSynchronizedPrintf(comm, "%d ", globalDof); 
1395:         }
1396:       }
1397:       PetscSynchronizedPrintf(comm, "\n\n"); 
1398:       PetscHSetIDestroy(&globalbcdofs);
1399:     }
1400:    for (k = 0; k < patch->nsubspaces; ++k) {
1401:       const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
1402:       PetscInt        nodesPerCell   = patch->nodesPerCell[k];
1403:       PetscInt        subspaceOffset = patch->subspaceOffsets[k];
1404:       PetscInt        bs             = patch->bs[k];

1406:       for (i = off; i < off + dof; ++i) {
1407:         /* Walk over the cells in this patch. */
1408:         const PetscInt c    = cellsArray[i];
1409:         PetscInt       cell = c;

1411:         /* TODO Change this to an IS */
1412:         if (cellNumbering) {
1413:           PetscSectionGetDof(cellNumbering, c, &cell);
1414:           if (cell <= 0) SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_OUTOFRANGE, "Cell %D doesn't appear in cell numbering map", c);
1415:           PetscSectionGetOffset(cellNumbering, c, &cell);
1416:         }
1417:         newCellsArray[i] = cell;
1418:         for (j = 0; j < nodesPerCell; ++j) {
1419:           /* For each global dof, map it into contiguous local storage. */
1420:           const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + subspaceOffset;
1421:           /* finally, loop over block size */
1422:           for (l = 0; l < bs; ++l) {
1423:             PetscInt  localDof;
1424:             PetscBool isGlobalBcDof, isArtificialBcDof;

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

1430:             /* if it's either, don't ever give it a local dof number */
1431:             if (isGlobalBcDof || isArtificialBcDof) {
1432:               dofsArray[globalIndex] = -1; /* don't use this in assembly in this patch */
1433:             } else {
1434:               PetscHMapIGet(ht, globalDof + l, &localDof);
1435:               if (localDof == -1) {
1436:                 localDof = localIndex++;
1437:                 PetscHMapISet(ht, globalDof + l, localDof);
1438:               }
1439:               if ( globalIndex >= numDofs ) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs);
1440:               /* And store. */
1441:               dofsArray[globalIndex] = localDof;
1442:             }

1444:             if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1445:               if (isGlobalBcDof) {
1446:                 dofsArrayWithArtificial[globalIndex] = -1; /* don't use this in assembly in this patch */
1447:               } else {
1448:                 PetscHMapIGet(htWithArtificial, globalDof + l, &localDof);
1449:                 if (localDof == -1) {
1450:                   localDof = localIndexWithArtificial++;
1451:                   PetscHMapISet(htWithArtificial, globalDof + l, localDof);
1452:                 }
1453:                 if ( globalIndex >= numDofs ) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs);
1454:                 /* And store.*/
1455:                 dofsArrayWithArtificial[globalIndex] = localDof;
1456:               }
1457:             }

1459:             if(isNonlinear) {
1460:               /* Build the dofmap for the function space with _all_ dofs,
1461:                  including those in any kind of boundary condition */
1462:               PetscHMapIGet(htWithAll, globalDof + l, &localDof);
1463:               if (localDof == -1) {
1464:                 localDof = localIndexWithAll++;
1465:                 PetscHMapISet(htWithAll, globalDof + l, localDof);
1466:               }
1467:               if ( globalIndex >= numDofs ) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs);
1468:               /* And store.*/
1469:               dofsArrayWithAll[globalIndex] = localDof;
1470:             }
1471:             globalIndex++;
1472:           }
1473:         }
1474:       }
1475:     }
1476:      /*How many local dofs in this patch? */
1477:    if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1478:      PetscHMapIGetSize(htWithArtificial, &dof);
1479:      PetscSectionSetDof(gtolCountsWithArtificial, v, dof);
1480:    }
1481:    if (isNonlinear) {
1482:      PetscHMapIGetSize(htWithAll, &dof);
1483:      PetscSectionSetDof(gtolCountsWithAll, v, dof);
1484:    }
1485:     PetscHMapIGetSize(ht, &dof);
1486:     PetscSectionSetDof(gtolCounts, v, dof);
1487:   }

1489:   DMDestroy(&dm);
1490:   if (globalIndex != numDofs) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Expected number of dofs (%d) doesn't match found number (%d)", numDofs, globalIndex);
1491:   PetscSectionSetUp(gtolCounts);
1492:   PetscSectionGetStorageSize(gtolCounts, &numGlobalDofs);
1493:   PetscMalloc1(numGlobalDofs, &globalDofsArray);

1495:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1496:     PetscSectionSetUp(gtolCountsWithArtificial);
1497:     PetscSectionGetStorageSize(gtolCountsWithArtificial, &numGlobalDofsWithArtificial);
1498:     PetscMalloc1(numGlobalDofsWithArtificial, &globalDofsArrayWithArtificial);
1499:   }
1500:   if (isNonlinear) {
1501:     PetscSectionSetUp(gtolCountsWithAll);
1502:     PetscSectionGetStorageSize(gtolCountsWithAll, &numGlobalDofsWithAll);
1503:     PetscMalloc1(numGlobalDofsWithAll, &globalDofsArrayWithAll);
1504:   }
1505:   /* Now populate the global to local map.  This could be merged into the above loop if we were willing to deal with reallocs. */
1506:   for (v = vStart; v < vEnd; ++v) {
1507:     PetscHashIter hi;
1508:     PetscInt      dof, off, Np, ooff, i, j, k, l;

1510:     PetscHMapIClear(ht);
1511:     PetscHMapIClear(htWithArtificial);
1512:     PetscHMapIClear(htWithAll);
1513:     PetscSectionGetDof(cellCounts, v, &dof);
1514:     PetscSectionGetOffset(cellCounts, v, &off);
1515:     PetscSectionGetDof(pointCounts, v, &Np);
1516:     PetscSectionGetOffset(pointCounts, v, &ooff);
1517:     if (dof <= 0) continue;

1519:     for (k = 0; k < patch->nsubspaces; ++k) {
1520:       const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
1521:       PetscInt        nodesPerCell   = patch->nodesPerCell[k];
1522:       PetscInt        subspaceOffset = patch->subspaceOffsets[k];
1523:       PetscInt        bs             = patch->bs[k];
1524:       PetscInt        goff;

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

1531:         if (cellNumbering) {PetscSectionGetOffset(cellNumbering, c, &cell);}
1532:         for (j = 0; j < nodesPerCell; ++j) {
1533:           for (l = 0; l < bs; ++l) {
1534:             const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset;
1535:             const PetscInt localDof  = dofsArray[key];
1536:             if (localDof >= 0) {PetscHMapISet(ht, globalDof, localDof);}
1537:             if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1538:               const PetscInt localDofWithArtificial = dofsArrayWithArtificial[key];
1539:               if (localDofWithArtificial >= 0) {
1540:                 PetscHMapISet(htWithArtificial, globalDof, localDofWithArtificial);
1541:               }
1542:             }
1543:             if (isNonlinear) {
1544:               const PetscInt localDofWithAll = dofsArrayWithAll[key];
1545:               if (localDofWithAll >= 0) {
1546:                 PetscHMapISet(htWithAll, globalDof, localDofWithAll);
1547:               }
1548:             }
1549:             key++;
1550:           }
1551:         }
1552:       }

1554:       /* Shove it in the output data structure. */
1555:       PetscSectionGetOffset(gtolCounts, v, &goff);
1556:       PetscHashIterBegin(ht, hi);
1557:       while (!PetscHashIterAtEnd(ht, hi)) {
1558:         PetscInt globalDof, localDof;

1560:         PetscHashIterGetKey(ht, hi, globalDof);
1561:         PetscHashIterGetVal(ht, hi, localDof);
1562:         if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof;
1563:         PetscHashIterNext(ht, hi);
1564:       }

1566:       if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1567:         PetscSectionGetOffset(gtolCountsWithArtificial, v, &goff);
1568:         PetscHashIterBegin(htWithArtificial, hi);
1569:         while (!PetscHashIterAtEnd(htWithArtificial, hi)) {
1570:           PetscInt globalDof, localDof;
1571:           PetscHashIterGetKey(htWithArtificial, hi, globalDof);
1572:           PetscHashIterGetVal(htWithArtificial, hi, localDof);
1573:           if (globalDof >= 0) globalDofsArrayWithArtificial[goff + localDof] = globalDof;
1574:           PetscHashIterNext(htWithArtificial, hi);
1575:         }
1576:       }
1577:       if (isNonlinear) {
1578:         PetscSectionGetOffset(gtolCountsWithAll, v, &goff);
1579:         PetscHashIterBegin(htWithAll, hi);
1580:         while (!PetscHashIterAtEnd(htWithAll, hi)) {
1581:           PetscInt globalDof, localDof;
1582:           PetscHashIterGetKey(htWithAll, hi, globalDof);
1583:           PetscHashIterGetVal(htWithAll, hi, localDof);
1584:           if (globalDof >= 0) globalDofsArrayWithAll[goff + localDof] = globalDof;
1585:           PetscHashIterNext(htWithAll, hi);
1586:         }
1587:       }

1589:       for (p = 0; p < Np; ++p) {
1590:         const PetscInt point = pointsArray[ooff + p];
1591:         PetscInt       globalDof, localDof;

1593:         PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof);
1594:         PetscHMapIGet(ht, globalDof, &localDof);
1595:         offsArray[(ooff + p)*Nf + k] = localDof;
1596:         if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1597:           PetscHMapIGet(htWithArtificial, globalDof, &localDof);
1598:           offsArrayWithArtificial[(ooff + p)*Nf + k] = localDof;
1599:         }
1600:         if (isNonlinear) {
1601:           PetscHMapIGet(htWithAll, globalDof, &localDof);
1602:           offsArrayWithAll[(ooff + p)*Nf + k] = localDof;
1603:         }
1604:       }
1605:     }

1607:     PetscHSetIDestroy(&globalBcs);
1608:     PetscHSetIDestroy(&ownedpts);
1609:     PetscHSetIDestroy(&seenpts);
1610:     PetscHSetIDestroy(&owneddofs);
1611:     PetscHSetIDestroy(&seendofs);
1612:     PetscHSetIDestroy(&artificialbcs);

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

1624:         if (cellNumbering) {PetscSectionGetOffset(cellNumbering, c, &cell);}
1625:         for (k = 0; k < patch->nsubspaces; ++k) {
1626:           const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
1627:           PetscInt        nodesPerCell   = patch->nodesPerCell[k];
1628:           PetscInt        subspaceOffset = patch->subspaceOffsets[k];
1629:           PetscInt        bs             = patch->bs[k];

1631:           for (j = 0; j < nodesPerCell; ++j) {
1632:             for (l = 0; l < bs; ++l) {
1633:               const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset;
1634:               PetscInt       localDof;

1636:               PetscHMapIGet(ht, globalDof, &localDof);
1637:               /* If it's not in the hash table, i.e. is a BC dof,
1638:                then the PetscHSetIMap above gives -1, which matches
1639:                exactly the convention for PETSc's matrix assembly to
1640:                ignore the dof. So we don't need to do anything here */
1641:               asmArray[asmKey] = localDof;
1642:               if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1643:                 PetscHMapIGet(htWithArtificial, globalDof, &localDof);
1644:                 asmArrayWithArtificial[asmKey] = localDof;
1645:               }
1646:               if (isNonlinear) {
1647:                 PetscHMapIGet(htWithAll, globalDof, &localDof);
1648:                 asmArrayWithAll[asmKey] = localDof;
1649:               }
1650:               asmKey++;
1651:             }
1652:           }
1653:         }
1654:       }
1655:     }
1656:   }
1657:   if (1 == patch->nsubspaces) {
1658:     PetscArraycpy(asmArray, dofsArray, numDofs);
1659:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1660:       PetscArraycpy(asmArrayWithArtificial, dofsArrayWithArtificial, numDofs);
1661:     }
1662:     if (isNonlinear) {
1663:       PetscArraycpy(asmArrayWithAll, dofsArrayWithAll, numDofs);
1664:     }
1665:   }

1667:   PetscHMapIDestroy(&ht);
1668:   PetscHMapIDestroy(&htWithArtificial);
1669:   PetscHMapIDestroy(&htWithAll);
1670:   ISRestoreIndices(cells, &cellsArray);
1671:   ISRestoreIndices(points, &pointsArray);
1672:   PetscFree(dofsArray);
1673:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1674:     PetscFree(dofsArrayWithArtificial);
1675:   }
1676:   if (isNonlinear) {
1677:     PetscFree(dofsArrayWithAll);
1678:   }
1679:   /* Create placeholder section for map from points to patch dofs */
1680:   PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection);
1681:   PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces);
1682:   if (patch->combined) {
1683:     PetscInt numFields;
1684:     PetscSectionGetNumFields(patch->dofSection[0], &numFields);
1685:     if (numFields != patch->nsubspaces) SETERRQ2(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Mismatch between number of section fields %D and number of subspaces %D", numFields, patch->nsubspaces);
1686:     PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd);
1687:     PetscSectionSetChart(patch->patchSection, pStart, pEnd);
1688:     for (p = pStart; p < pEnd; ++p) {
1689:       PetscInt dof, fdof, f;

1691:       PetscSectionGetDof(patch->dofSection[0], p, &dof);
1692:       PetscSectionSetDof(patch->patchSection, p, dof);
1693:       for (f = 0; f < patch->nsubspaces; ++f) {
1694:         PetscSectionGetFieldDof(patch->dofSection[0], p, f, &fdof);
1695:         PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);
1696:       }
1697:     }
1698:   } else {
1699:     PetscInt pStartf, pEndf, f;
1700:     pStart = PETSC_MAX_INT;
1701:     pEnd = PETSC_MIN_INT;
1702:     for (f = 0; f < patch->nsubspaces; ++f) {
1703:       PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);
1704:       pStart = PetscMin(pStart, pStartf);
1705:       pEnd = PetscMax(pEnd, pEndf);
1706:     }
1707:     PetscSectionSetChart(patch->patchSection, pStart, pEnd);
1708:     for (f = 0; f < patch->nsubspaces; ++f) {
1709:       PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);
1710:       for (p = pStartf; p < pEndf; ++p) {
1711:         PetscInt fdof;
1712:         PetscSectionGetDof(patch->dofSection[f], p, &fdof);
1713:         PetscSectionAddDof(patch->patchSection, p, fdof);
1714:         PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);
1715:       }
1716:     }
1717:   }
1718:   PetscSectionSetUp(patch->patchSection);
1719:   PetscSectionSetUseFieldOffsets(patch->patchSection, PETSC_TRUE);
1720:   /* Replace cell indices with firedrake-numbered ones. */
1721:   ISGeneralSetIndices(cells, numCells, (const PetscInt *) newCellsArray, PETSC_OWN_POINTER);
1722:   ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofs, globalDofsArray, PETSC_OWN_POINTER, &patch->gtol);
1723:   PetscObjectSetName((PetscObject) patch->gtol, "Global Indices");
1724:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_g2l_view", patch->classname);
1725:   PetscSectionViewFromOptions(patch->gtolCounts, (PetscObject) pc, option);
1726:   ISViewFromOptions(patch->gtol, (PetscObject) pc, option);
1727:   ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArray, PETSC_OWN_POINTER, &patch->dofs);
1728:   ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArray, PETSC_OWN_POINTER, &patch->offs);
1729:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1730:     ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithArtificial, globalDofsArrayWithArtificial, PETSC_OWN_POINTER, &patch->gtolWithArtificial);
1731:     ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithArtificial, PETSC_OWN_POINTER, &patch->dofsWithArtificial);
1732:     ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArrayWithArtificial, PETSC_OWN_POINTER, &patch->offsWithArtificial);
1733:   }
1734:   if (isNonlinear) {
1735:     ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithAll, globalDofsArrayWithAll, PETSC_OWN_POINTER, &patch->gtolWithAll);
1736:     ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithAll, PETSC_OWN_POINTER, &patch->dofsWithAll);
1737:     ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArrayWithAll, PETSC_OWN_POINTER, &patch->offsWithAll);
1738:   }
1739:   return(0);
1740: }

1742: static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat, PetscBool withArtificial)
1743: {
1744:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
1745:   PetscBool      flg;
1746:   PetscInt       csize, rsize;
1747:   const char    *prefix = NULL;

1751:   if(withArtificial) {
1752:     /* would be nice if we could create a rectangular matrix of size numDofsWithArtificial x numDofs here */
1753:     PetscInt pStart;
1754:     PetscSectionGetChart(patch->gtolCountsWithArtificial, &pStart, NULL);
1755:     PetscSectionGetDof(patch->gtolCountsWithArtificial, point + pStart, &rsize);
1756:     csize = rsize;
1757:   } else {
1758:     PetscInt pStart;
1759:     PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);
1760:     PetscSectionGetDof(patch->gtolCounts, point + pStart, &rsize);
1761:     csize = rsize;
1762:   }

1764:   MatCreate(PETSC_COMM_SELF, mat);
1765:   PCGetOptionsPrefix(pc, &prefix);
1766:   MatSetOptionsPrefix(*mat, prefix);
1767:   MatAppendOptionsPrefix(*mat, "pc_patch_sub_");
1768:   if (patch->sub_mat_type)       {MatSetType(*mat, patch->sub_mat_type);}
1769:   else if (!patch->sub_mat_type) {MatSetType(*mat, MATDENSE);}
1770:   MatSetSizes(*mat, rsize, csize, rsize, csize);
1771:   PetscObjectTypeCompare((PetscObject) *mat, MATDENSE, &flg);
1772:   if (!flg) {PetscObjectTypeCompare((PetscObject)*mat, MATSEQDENSE, &flg);}
1773:   /* Sparse patch matrices */
1774:   if (!flg) {
1775:     PetscBT         bt;
1776:     PetscInt       *dnnz      = NULL;
1777:     const PetscInt *dofsArray = NULL;
1778:     PetscInt        pStart, pEnd, ncell, offset, c, i, j;

1780:     if(withArtificial) {
1781:       ISGetIndices(patch->dofsWithArtificial, &dofsArray);
1782:     } else {
1783:       ISGetIndices(patch->dofs, &dofsArray);
1784:     }
1785:     PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);
1786:     point += pStart;
1787:     if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);
1788:     PetscSectionGetDof(patch->cellCounts, point, &ncell);
1789:     PetscSectionGetOffset(patch->cellCounts, point, &offset);
1790:     PetscLogEventBegin(PC_Patch_Prealloc, pc, 0, 0, 0);
1791:     /* A PetscBT uses N^2 bits to store the sparsity pattern on a
1792:      * patch. This is probably OK if the patches are not too big,
1793:      * but uses too much memory. We therefore switch based on rsize. */
1794:     if (rsize < 3000) { /* FIXME: I picked this switch value out of my hat */
1795:       PetscScalar *zeroes;
1796:       PetscInt rows;

1798:       PetscCalloc1(rsize, &dnnz);
1799:       PetscBTCreate(rsize*rsize, &bt);
1800:       for (c = 0; c < ncell; ++c) {
1801:         const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell;
1802:         for (i = 0; i < patch->totalDofsPerCell; ++i) {
1803:           const PetscInt row = idx[i];
1804:           if (row < 0) continue;
1805:           for (j = 0; j < patch->totalDofsPerCell; ++j) {
1806:             const PetscInt col = idx[j];
1807:             const PetscInt key = row*rsize + col;
1808:             if (col < 0) continue;
1809:             if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1810:           }
1811:         }
1812:       }

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

1819:         PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1820:         PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1821:         ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1822:         ISGetIndices(patch->intFacets, &intFacetsArray);
1823:         for (i = 0; i < numIntFacets; i++) {
1824:           const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0];
1825:           const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1];
1826:           PetscInt       celli, cellj;

1828:           for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1829:             const PetscInt row = dofsArray[(offset + cell0)*patch->totalDofsPerCell + celli];
1830:             if (row < 0) continue;
1831:             for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1832:               const PetscInt col = dofsArray[(offset + cell1)*patch->totalDofsPerCell + cellj];
1833:               const PetscInt key = row*rsize + col;
1834:               if (col < 0) continue;
1835:               if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1836:             }
1837:           }

1839:           for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1840:             const PetscInt row = dofsArray[(offset + cell1)*patch->totalDofsPerCell + celli];
1841:             if (row < 0) continue;
1842:             for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1843:               const PetscInt col = dofsArray[(offset + cell0)*patch->totalDofsPerCell + cellj];
1844:               const PetscInt key = row*rsize + col;
1845:               if (col < 0) continue;
1846:               if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1847:             }
1848:           }
1849:         }
1850:       }
1851:       PetscBTDestroy(&bt);
1852:       MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL);
1853:       PetscFree(dnnz);

1855:       PetscCalloc1(patch->totalDofsPerCell*patch->totalDofsPerCell, &zeroes);
1856:       for (c = 0; c < ncell; ++c) {
1857:         const PetscInt *idx = &dofsArray[(offset + c)*patch->totalDofsPerCell];
1858:         MatSetValues(*mat, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, zeroes, INSERT_VALUES);
1859:       }
1860:       MatGetLocalSize(*mat, &rows, NULL);
1861:       for (i = 0; i < rows; ++i) {
1862:         MatSetValues(*mat, 1, &i, 1, &i, zeroes, INSERT_VALUES);
1863:       }

1865:       if (patch->usercomputeopintfacet) {
1866:         const PetscInt *intFacetsArray = NULL;
1867:         PetscInt i, numIntFacets, intFacetOffset;
1868:         const PetscInt *facetCells = NULL;

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

1884:       MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
1885:       MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);

1887:       PetscFree(zeroes);

1889:     } else { /* rsize too big, use MATPREALLOCATOR */
1890:       Mat preallocator;
1891:       PetscScalar* vals;

1893:       PetscCalloc1(patch->totalDofsPerCell*patch->totalDofsPerCell, &vals);
1894:       MatCreate(PETSC_COMM_SELF, &preallocator);
1895:       MatSetType(preallocator, MATPREALLOCATOR);
1896:       MatSetSizes(preallocator, rsize, rsize, rsize, rsize);
1897:       MatSetUp(preallocator);

1899:       for (c = 0; c < ncell; ++c) {
1900:         const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell;
1901:         MatSetValues(preallocator, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, vals, INSERT_VALUES);
1902:       }

1904:       if (patch->usercomputeopintfacet) {
1905:         const PetscInt *intFacetsArray = NULL;
1906:         PetscInt        i, numIntFacets, intFacetOffset;
1907:         const PetscInt *facetCells = NULL;

1909:         PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1910:         PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1911:         ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1912:         ISGetIndices(patch->intFacets, &intFacetsArray);
1913:         for (i = 0; i < numIntFacets; i++) {
1914:           const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0];
1915:           const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1];
1916:           const PetscInt *cell0idx = &dofsArray[(offset + cell0)*patch->totalDofsPerCell];
1917:           const PetscInt *cell1idx = &dofsArray[(offset + cell1)*patch->totalDofsPerCell];
1918:           MatSetValues(preallocator, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, vals, INSERT_VALUES);
1919:           MatSetValues(preallocator, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, vals, INSERT_VALUES);
1920:         }
1921:       }

1923:       PetscFree(vals);
1924:       MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1925:       MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1926:       MatPreallocatorPreallocate(preallocator, PETSC_TRUE, *mat);
1927:       MatDestroy(&preallocator);
1928:     }
1929:     PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0);
1930:     if(withArtificial) {
1931:       ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);
1932:     } else {
1933:       ISRestoreIndices(patch->dofs, &dofsArray);
1934:     }
1935:   }
1936:   MatSetUp(*mat);
1937:   return(0);
1938: }

1940: 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)
1941: {
1942:   PC_PATCH       *patch = (PC_PATCH *) pc->data;
1943:   DM              dm, plex;
1944:   PetscSection    s;
1945:   const PetscInt *parray, *oarray;
1946:   PetscInt        Nf = patch->nsubspaces, Np, poff, p, f;
1947:   PetscErrorCode  ierr;

1950:   if (patch->precomputeElementTensors) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Precomputing element tensors not implemented with DMPlex compute operator");
1951:   PCGetDM(pc, &dm);
1952:   DMConvert(dm, DMPLEX, &plex);
1953:   dm = plex;
1954:   DMGetLocalSection(dm, &s);
1955:   /* Set offset into patch */
1956:   PetscSectionGetDof(patch->pointCounts, patchNum, &Np);
1957:   PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);
1958:   ISGetIndices(patch->points, &parray);
1959:   ISGetIndices(patch->offs,   &oarray);
1960:   for (f = 0; f < Nf; ++f) {
1961:     for (p = 0; p < Np; ++p) {
1962:       const PetscInt point = parray[poff+p];
1963:       PetscInt       dof;

1965:       PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);
1966:       PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);
1967:       if (patch->nsubspaces == 1) {PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);}
1968:       else                        {PetscSectionSetOffset(patch->patchSection, point, -1);}
1969:     }
1970:   }
1971:   ISRestoreIndices(patch->points, &parray);
1972:   ISRestoreIndices(patch->offs,   &oarray);
1973:   if (patch->viewSection) {ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);}
1974:   DMPlexComputeResidual_Patch_Internal(dm, patch->patchSection, cellIS, 0.0, x, NULL, F, ctx);
1975:   DMDestroy(&dm);
1976:   return(0);
1977: }

1979: PetscErrorCode PCPatchComputeFunction_Internal(PC pc, Vec x, Vec F, PetscInt point)
1980: {
1981:   PC_PATCH       *patch = (PC_PATCH *) pc->data;
1982:   const PetscInt *dofsArray;
1983:   const PetscInt *dofsArrayWithAll;
1984:   const PetscInt *cellsArray;
1985:   PetscInt        ncell, offset, pStart, pEnd;
1986:   PetscErrorCode  ierr;

1989:   PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);
1990:   if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n");
1991:   ISGetIndices(patch->dofs, &dofsArray);
1992:   ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);
1993:   ISGetIndices(patch->cells, &cellsArray);
1994:   PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);

1996:   point += pStart;
1997:   if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);

1999:   PetscSectionGetDof(patch->cellCounts, point, &ncell);
2000:   PetscSectionGetOffset(patch->cellCounts, point, &offset);
2001:   if (ncell <= 0) {
2002:     PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2003:     return(0);
2004:   }
2005:   VecSet(F, 0.0);
2006:   PetscStackPush("PCPatch user callback");
2007:   /* Cannot reuse the same IS because the geometry info is being cached in it */
2008:   ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);
2009:   patch->usercomputef(pc, point, x, F, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell,
2010:                                                                                             dofsArrayWithAll + offset*patch->totalDofsPerCell,
2011:                                                                                             patch->usercomputefctx);
2012:   PetscStackPop;
2013:   ISDestroy(&patch->cellIS);
2014:   ISRestoreIndices(patch->dofs, &dofsArray);
2015:   ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);
2016:   ISRestoreIndices(patch->cells, &cellsArray);
2017:   if (patch->viewMatrix) {
2018:     char name[PETSC_MAX_PATH_LEN];

2020:     PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch vector for Point %D", point);
2021:     PetscObjectSetName((PetscObject) F, name);
2022:     ObjectView((PetscObject) F, patch->viewerMatrix, patch->formatMatrix);
2023:   }
2024:   PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2025:   return(0);
2026: }

2028: 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)
2029: {
2030:   PC_PATCH       *patch = (PC_PATCH *) pc->data;
2031:   DM              dm, plex;
2032:   PetscSection    s;
2033:   const PetscInt *parray, *oarray;
2034:   PetscInt        Nf = patch->nsubspaces, Np, poff, p, f;
2035:   PetscErrorCode  ierr;

2038:   PCGetDM(pc, &dm);
2039:   DMConvert(dm, DMPLEX, &plex);
2040:   dm = plex;
2041:   DMGetLocalSection(dm, &s);
2042:   /* Set offset into patch */
2043:   PetscSectionGetDof(patch->pointCounts, patchNum, &Np);
2044:   PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);
2045:   ISGetIndices(patch->points, &parray);
2046:   ISGetIndices(patch->offs,   &oarray);
2047:   for (f = 0; f < Nf; ++f) {
2048:     for (p = 0; p < Np; ++p) {
2049:       const PetscInt point = parray[poff+p];
2050:       PetscInt       dof;

2052:       PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);
2053:       PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);
2054:       if (patch->nsubspaces == 1) {PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);}
2055:       else                        {PetscSectionSetOffset(patch->patchSection, point, -1);}
2056:     }
2057:   }
2058:   ISRestoreIndices(patch->points, &parray);
2059:   ISRestoreIndices(patch->offs,   &oarray);
2060:   if (patch->viewSection) {ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);}
2061:   /* TODO Shut off MatViewFromOptions() in MatAssemblyEnd() here */
2062:   DMPlexComputeJacobian_Patch_Internal(dm, patch->patchSection, patch->patchSection, cellIS, 0.0, 0.0, x, NULL, J, J, ctx);
2063:   DMDestroy(&dm);
2064:   return(0);
2065: }

2067: /* This function zeros mat on entry */
2068: PetscErrorCode PCPatchComputeOperator_Internal(PC pc, Vec x, Mat mat, PetscInt point, PetscBool withArtificial)
2069: {
2070:   PC_PATCH       *patch = (PC_PATCH *) pc->data;
2071:   const PetscInt *dofsArray;
2072:   const PetscInt *dofsArrayWithAll = NULL;
2073:   const PetscInt *cellsArray;
2074:   PetscInt        ncell, offset, pStart, pEnd, numIntFacets, intFacetOffset;
2075:   PetscBool       isNonlinear;
2076:   PetscErrorCode  ierr;

2079:   PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);
2080:   isNonlinear = patch->isNonlinear;
2081:   if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n");
2082:   if(withArtificial) {
2083:     ISGetIndices(patch->dofsWithArtificial, &dofsArray);
2084:   } else {
2085:     ISGetIndices(patch->dofs, &dofsArray);
2086:   }
2087:   if (isNonlinear) {
2088:     ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);
2089:   }
2090:   ISGetIndices(patch->cells, &cellsArray);
2091:   PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);

2093:   point += pStart;
2094:   if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);

2096:   PetscSectionGetDof(patch->cellCounts, point, &ncell);
2097:   PetscSectionGetOffset(patch->cellCounts, point, &offset);
2098:   if (ncell <= 0) {
2099:     PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2100:     return(0);
2101:   }
2102:   MatZeroEntries(mat);
2103:   if (patch->precomputeElementTensors) {
2104:     PetscInt           i;
2105:     PetscInt           ndof = patch->totalDofsPerCell;
2106:     const PetscScalar *elementTensors;

2108:     VecGetArrayRead(patch->cellMats, &elementTensors);
2109:     for (i = 0; i < ncell; i++) {
2110:       const PetscInt     cell = cellsArray[i + offset];
2111:       const PetscInt    *idx  = dofsArray + (offset + i)*ndof;
2112:       const PetscScalar *v    = elementTensors + patch->precomputedTensorLocations[cell]*ndof*ndof;
2113:       MatSetValues(mat, ndof, idx, ndof, idx, v, ADD_VALUES);
2114:     }
2115:     VecRestoreArrayRead(patch->cellMats, &elementTensors);
2116:     MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
2117:     MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
2118:   } else {
2119:     PetscStackPush("PCPatch user callback");
2120:     /* Cannot reuse the same IS because the geometry info is being cached in it */
2121:     ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);
2122:     patch->usercomputeop(pc, point, x, mat, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell, dofsArrayWithAll ? dofsArrayWithAll + offset*patch->totalDofsPerCell : NULL, patch->usercomputeopctx);
2123:   }
2124:   if (patch->usercomputeopintfacet) {
2125:     PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
2126:     PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
2127:     if (numIntFacets > 0) {
2128:       /* For each interior facet, grab the two cells (in local numbering, and concatenate dof numberings for those cells) */
2129:       PetscInt       *facetDofs = NULL, *facetDofsWithAll = NULL;
2130:       const PetscInt *intFacetsArray = NULL;
2131:       PetscInt        idx = 0;
2132:       PetscInt        i, c, d;
2133:       PetscInt        fStart;
2134:       DM              dm, plex;
2135:       IS              facetIS = NULL;
2136:       const PetscInt *facetCells = NULL;

2138:       ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
2139:       ISGetIndices(patch->intFacets, &intFacetsArray);
2140:       PCGetDM(pc, &dm);
2141:       DMConvert(dm, DMPLEX, &plex);
2142:       dm = plex;
2143:       DMPlexGetHeightStratum(dm, 1, &fStart, NULL);
2144:       /* FIXME: Pull this malloc out. */
2145:       PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofs);
2146:       if (dofsArrayWithAll) {
2147:         PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofsWithAll);
2148:       }
2149:       if (patch->precomputeElementTensors) {
2150:         PetscInt           nFacetDof = 2*patch->totalDofsPerCell;
2151:         const PetscScalar *elementTensors;

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

2155:         for (i = 0; i < numIntFacets; i++) {
2156:           const PetscInt     facet = intFacetsArray[i + intFacetOffset];
2157:           const PetscScalar *v     = elementTensors + patch->precomputedIntFacetTensorLocations[facet - fStart]*nFacetDof*nFacetDof;
2158:           idx = 0;
2159:           /*
2160:            * 0--1
2161:            * |\-|
2162:            * |+\|
2163:            * 2--3
2164:            * [0, 2, 3, 0, 1, 3]
2165:            */
2166:           for (c = 0; c < 2; c++) {
2167:             const PetscInt cell = facetCells[2*(intFacetOffset + i) + c];
2168:             for (d = 0; d < patch->totalDofsPerCell; d++) {
2169:               facetDofs[idx] = dofsArray[(offset + cell)*patch->totalDofsPerCell + d];
2170:               idx++;
2171:             }
2172:           }
2173:           MatSetValues(mat, nFacetDof, facetDofs, nFacetDof, facetDofs, v, ADD_VALUES);
2174:         }
2175:         VecRestoreArrayRead(patch->intFacetMats, &elementTensors);
2176:       } else {
2177:         /*
2178:          * 0--1
2179:          * |\-|
2180:          * |+\|
2181:          * 2--3
2182:          * [0, 2, 3, 0, 1, 3]
2183:          */
2184:         for (i = 0; i < numIntFacets; i++) {
2185:           for (c = 0; c < 2; c++) {
2186:             const PetscInt cell = facetCells[2*(intFacetOffset + i) + c];
2187:             for (d = 0; d < patch->totalDofsPerCell; d++) {
2188:               facetDofs[idx] = dofsArray[(offset + cell)*patch->totalDofsPerCell + d];
2189:               if (dofsArrayWithAll) {
2190:                 facetDofsWithAll[idx] = dofsArrayWithAll[(offset + cell)*patch->totalDofsPerCell + d];
2191:               }
2192:               idx++;
2193:             }
2194:           }
2195:         }
2196:         ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray + intFacetOffset, PETSC_USE_POINTER, &facetIS);
2197:         patch->usercomputeopintfacet(pc, point, x, mat, facetIS, 2*numIntFacets*patch->totalDofsPerCell, facetDofs, facetDofsWithAll, patch->usercomputeopintfacetctx);
2198:         ISDestroy(&facetIS);
2199:       }
2200:       ISRestoreIndices(patch->intFacetsToPatchCell, &facetCells);
2201:       ISRestoreIndices(patch->intFacets, &intFacetsArray);
2202:       PetscFree(facetDofs);
2203:       PetscFree(facetDofsWithAll);
2204:       DMDestroy(&dm);
2205:     }
2206:   }

2208:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
2209:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);

2211:   if (!(withArtificial || isNonlinear) && patch->denseinverse) {
2212:     MatFactorInfo info;
2213:     PetscBool     flg;
2214:     PetscObjectTypeCompare((PetscObject)mat, MATSEQDENSE, &flg);
2215:     if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Invalid Mat type for dense inverse");
2216:     MatFactorInfoInitialize(&info);
2217:     MatLUFactor(mat, NULL, NULL, &info);
2218:     MatSeqDenseInvertFactors_Private(mat);
2219:   }
2220:   PetscStackPop;
2221:   ISDestroy(&patch->cellIS);
2222:   if(withArtificial) {
2223:     ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);
2224:   } else {
2225:     ISRestoreIndices(patch->dofs, &dofsArray);
2226:   }
2227:   if (isNonlinear) {
2228:     ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);
2229:   }
2230:   ISRestoreIndices(patch->cells, &cellsArray);
2231:   if (patch->viewMatrix) {
2232:     char name[PETSC_MAX_PATH_LEN];

2234:     PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch matrix for Point %D", point);
2235:     PetscObjectSetName((PetscObject) mat, name);
2236:     ObjectView((PetscObject) mat, patch->viewerMatrix, patch->formatMatrix);
2237:   }
2238:   PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2239:   return(0);
2240: }

2242: static PetscErrorCode MatSetValues_PCPatch_Private(Mat mat, PetscInt m, const PetscInt idxm[],
2243:                                                    PetscInt n, const PetscInt idxn[], const PetscScalar *v, InsertMode addv)
2244: {
2245:   Vec            data;
2246:   PetscScalar   *array;
2247:   PetscInt       bs, nz, i, j, cell;

2250:   MatShellGetContext(mat, &data);
2251:   VecGetBlockSize(data, &bs);
2252:   VecGetSize(data, &nz);
2253:   VecGetArray(data, &array);
2254:   if (m != n) SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Only for square insertion");
2255:   cell = (PetscInt)(idxm[0]/bs); /* use the fact that this is called once per cell */
2256:   for (i = 0; i < m; i++) {
2257:     if (idxm[i] != idxn[i]) SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Row and column indices must match!");
2258:     for (j = 0; j < n; j++) {
2259:       const PetscScalar v_ = v[i*bs + j];
2260:       /* Indexing is special to the data structure we have! */
2261:       if (addv == INSERT_VALUES) {
2262:         array[cell*bs*bs + i*bs + j] = v_;
2263:       } else {
2264:         array[cell*bs*bs + i*bs + j] += v_;
2265:       }
2266:     }
2267:   }
2268:   VecRestoreArray(data, &array);
2269:   return(0);
2270: }

2272: static PetscErrorCode PCPatchPrecomputePatchTensors_Private(PC pc)
2273: {
2274:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
2275:   const PetscInt *cellsArray;
2276:   PetscInt        ncell, offset;
2277:   const PetscInt *dofMapArray;
2278:   PetscInt        i, j;
2279:   IS              dofMap;
2280:   IS              cellIS;
2281:   const PetscInt  ndof  = patch->totalDofsPerCell;
2282:   PetscErrorCode  ierr;
2283:   Mat             vecMat;
2284:   PetscInt        cStart, cEnd;
2285:   DM              dm, plex;


2288:   ISGetSize(patch->cells, &ncell);
2289:   if (!ncell) { /* No cells to assemble over -> skip */
2290:     return(0);
2291:   }

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

2295:   PCGetDM(pc, &dm);
2296:   DMConvert(dm, DMPLEX, &plex);
2297:   dm = plex;
2298:   if (!patch->allCells) {
2299:     PetscHSetI      cells;
2300:     PetscHashIter   hi;
2301:     PetscInt pStart, pEnd;
2302:     PetscInt *allCells = NULL;
2303:     PetscHSetICreate(&cells);
2304:     ISGetIndices(patch->cells, &cellsArray);
2305:     PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);
2306:     for (i = pStart; i < pEnd; i++) {
2307:       PetscSectionGetDof(patch->cellCounts, i, &ncell);
2308:       PetscSectionGetOffset(patch->cellCounts, i, &offset);
2309:       if (ncell <= 0) continue;
2310:       for (j = 0; j < ncell; j++) {
2311:         PetscHSetIAdd(cells, cellsArray[offset + j]);
2312:       }
2313:     }
2314:     ISRestoreIndices(patch->cells, &cellsArray);
2315:     PetscHSetIGetSize(cells, &ncell);
2316:     PetscMalloc1(ncell, &allCells);
2317:     DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2318:     PetscMalloc1(cEnd-cStart, &patch->precomputedTensorLocations);
2319:     i = 0;
2320:     PetscHashIterBegin(cells, hi);
2321:     while (!PetscHashIterAtEnd(cells, hi)) {
2322:       PetscHashIterGetKey(cells, hi, allCells[i]);
2323:       patch->precomputedTensorLocations[allCells[i]] = i;
2324:       PetscHashIterNext(cells, hi);
2325:       i++;
2326:     }
2327:     PetscHSetIDestroy(&cells);
2328:     ISCreateGeneral(PETSC_COMM_SELF, ncell, allCells, PETSC_OWN_POINTER, &patch->allCells);
2329:   }
2330:   ISGetSize(patch->allCells, &ncell);
2331:   if (!patch->cellMats) {
2332:     VecCreateSeq(PETSC_COMM_SELF, ncell*ndof*ndof, &patch->cellMats);
2333:     VecSetBlockSize(patch->cellMats, ndof);
2334:   }
2335:   VecSet(patch->cellMats, 0);

2337:   MatCreateShell(PETSC_COMM_SELF, ncell*ndof, ncell*ndof, ncell*ndof, ncell*ndof,
2338:                         (void*)patch->cellMats, &vecMat);
2339:   MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void(*)(void))&MatSetValues_PCPatch_Private);
2340:   ISGetSize(patch->allCells, &ncell);
2341:   ISCreateStride(PETSC_COMM_SELF, ndof*ncell, 0, 1, &dofMap);
2342:   ISGetIndices(dofMap, &dofMapArray);
2343:   ISGetIndices(patch->allCells, &cellsArray);
2344:   ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray, PETSC_USE_POINTER, &cellIS);
2345:   PetscStackPush("PCPatch user callback");
2346:   /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2347:   patch->usercomputeop(pc, -1, NULL, vecMat, cellIS, ndof*ncell, dofMapArray, NULL, patch->usercomputeopctx);
2348:   PetscStackPop;
2349:   ISDestroy(&cellIS);
2350:   MatDestroy(&vecMat);
2351:   ISRestoreIndices(patch->allCells, &cellsArray);
2352:   ISRestoreIndices(dofMap, &dofMapArray);
2353:   ISDestroy(&dofMap);

2355:   if (patch->usercomputeopintfacet) {
2356:     PetscInt nIntFacets;
2357:     IS       intFacetsIS;
2358:     const PetscInt *intFacetsArray = NULL;
2359:     if (!patch->allIntFacets) {
2360:       PetscHSetI      facets;
2361:       PetscHashIter   hi;
2362:       PetscInt pStart, pEnd, fStart, fEnd;
2363:       PetscInt *allIntFacets = NULL;
2364:       PetscHSetICreate(&facets);
2365:       ISGetIndices(patch->intFacets, &intFacetsArray);
2366:       PetscSectionGetChart(patch->intFacetCounts, &pStart, &pEnd);
2367:       DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2368:       for (i = pStart; i < pEnd; i++) {
2369:         PetscSectionGetDof(patch->intFacetCounts, i, &nIntFacets);
2370:         PetscSectionGetOffset(patch->intFacetCounts, i, &offset);
2371:         if (nIntFacets <= 0) continue;
2372:         for (j = 0; j < nIntFacets; j++) {
2373:           PetscHSetIAdd(facets, intFacetsArray[offset + j]);
2374:         }
2375:       }
2376:       ISRestoreIndices(patch->intFacets, &intFacetsArray);
2377:       PetscHSetIGetSize(facets, &nIntFacets);
2378:       PetscMalloc1(nIntFacets, &allIntFacets);
2379:       PetscMalloc1(fEnd-fStart, &patch->precomputedIntFacetTensorLocations);
2380:       i = 0;
2381:       PetscHashIterBegin(facets, hi);
2382:       while (!PetscHashIterAtEnd(facets, hi)) {
2383:         PetscHashIterGetKey(facets, hi, allIntFacets[i]);
2384:         patch->precomputedIntFacetTensorLocations[allIntFacets[i] - fStart] = i;
2385:         PetscHashIterNext(facets, hi);
2386:         i++;
2387:       }
2388:       PetscHSetIDestroy(&facets);
2389:       ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, allIntFacets, PETSC_OWN_POINTER, &patch->allIntFacets);
2390:     }
2391:     ISGetSize(patch->allIntFacets, &nIntFacets);
2392:     if (!patch->intFacetMats) {
2393:       VecCreateSeq(PETSC_COMM_SELF, nIntFacets*ndof*ndof*4, &patch->intFacetMats);
2394:       VecSetBlockSize(patch->intFacetMats, ndof*2);
2395:     }
2396:     VecSet(patch->intFacetMats, 0);

2398:     MatCreateShell(PETSC_COMM_SELF, nIntFacets*ndof*2, nIntFacets*ndof*2, nIntFacets*ndof*2, nIntFacets*ndof*2,
2399:                           (void*)patch->intFacetMats, &vecMat);
2400:     MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void(*)(void))&MatSetValues_PCPatch_Private);
2401:     ISCreateStride(PETSC_COMM_SELF, 2*ndof*nIntFacets, 0, 1, &dofMap);
2402:     ISGetIndices(dofMap, &dofMapArray);
2403:     ISGetIndices(patch->allIntFacets, &intFacetsArray);
2404:     ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, intFacetsArray, PETSC_USE_POINTER, &intFacetsIS);
2405:     PetscStackPush("PCPatch user callback (interior facets)");
2406:     /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2407:     patch->usercomputeopintfacet(pc, -1, NULL, vecMat, intFacetsIS, 2*ndof*nIntFacets, dofMapArray, NULL, patch->usercomputeopintfacetctx);
2408:     PetscStackPop;
2409:     ISDestroy(&intFacetsIS);
2410:     MatDestroy(&vecMat);
2411:     ISRestoreIndices(patch->allIntFacets, &intFacetsArray);
2412:     ISRestoreIndices(dofMap, &dofMapArray);
2413:     ISDestroy(&dofMap);
2414:   }
2415:   DMDestroy(&dm);
2416:   PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);

2418:   return(0);
2419: }

2421: PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat, PatchScatterType scattertype)
2422: {
2423:   PC_PATCH          *patch     = (PC_PATCH *) pc->data;
2424:   const PetscScalar *xArray    = NULL;
2425:   PetscScalar       *yArray    = NULL;
2426:   const PetscInt    *gtolArray = NULL;
2427:   PetscInt           dof, offset, lidx;
2428:   PetscErrorCode     ierr;

2431:   VecGetArrayRead(x, &xArray);
2432:   VecGetArray(y, &yArray);
2433:   if (scattertype == SCATTER_WITHARTIFICIAL) {
2434:     PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);
2435:     PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offset);
2436:     ISGetIndices(patch->gtolWithArtificial, &gtolArray);
2437:   } else if (scattertype == SCATTER_WITHALL) {
2438:     PetscSectionGetDof(patch->gtolCountsWithAll, p, &dof);
2439:     PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offset);
2440:     ISGetIndices(patch->gtolWithAll, &gtolArray);
2441:   } else {
2442:     PetscSectionGetDof(patch->gtolCounts, p, &dof);
2443:     PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2444:     ISGetIndices(patch->gtol, &gtolArray);
2445:   }
2446:   if (mode == INSERT_VALUES && scat != SCATTER_FORWARD) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't insert if not scattering forward\n");
2447:   if (mode == ADD_VALUES    && scat != SCATTER_REVERSE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't add if not scattering reverse\n");
2448:   for (lidx = 0; lidx < dof; ++lidx) {
2449:     const PetscInt gidx = gtolArray[offset+lidx];

2451:     if (mode == INSERT_VALUES) yArray[lidx]  = xArray[gidx]; /* Forward */
2452:     else                       yArray[gidx] += xArray[lidx]; /* Reverse */
2453:   }
2454:   if (scattertype == SCATTER_WITHARTIFICIAL) {
2455:     ISRestoreIndices(patch->gtolWithArtificial, &gtolArray);
2456:   } else if (scattertype == SCATTER_WITHALL) {
2457:     ISRestoreIndices(patch->gtolWithAll, &gtolArray);
2458:   } else {
2459:     ISRestoreIndices(patch->gtol, &gtolArray);
2460:   }
2461:   VecRestoreArrayRead(x, &xArray);
2462:   VecRestoreArray(y, &yArray);
2463:   return(0);
2464: }

2466: static PetscErrorCode PCSetUp_PATCH_Linear(PC pc)
2467: {
2468:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2469:   const char    *prefix;
2470:   PetscInt       i;

2474:   if (!pc->setupcalled) {
2475:     if (!patch->save_operators && patch->denseinverse) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Can't have dense inverse without save operators");
2476:     if (!patch->denseinverse) {
2477:       PetscMalloc1(patch->npatch, &patch->solver);
2478:       PCGetOptionsPrefix(pc, &prefix);
2479:       for (i = 0; i < patch->npatch; ++i) {
2480:         KSP ksp;
2481:         PC  subpc;

2483:         KSPCreate(PETSC_COMM_SELF, &ksp);
2484:         KSPSetErrorIfNotConverged(ksp, pc->erroriffailure);
2485:         KSPSetOptionsPrefix(ksp, prefix);
2486:         KSPAppendOptionsPrefix(ksp, "sub_");
2487:         PetscObjectIncrementTabLevel((PetscObject) ksp, (PetscObject) pc, 1);
2488:         KSPGetPC(ksp, &subpc);
2489:         PetscObjectIncrementTabLevel((PetscObject) subpc, (PetscObject) pc, 1);
2490:         PetscLogObjectParent((PetscObject) pc, (PetscObject) ksp);
2491:         patch->solver[i] = (PetscObject) ksp;
2492:       }
2493:     }
2494:   }
2495:   if (patch->save_operators) {
2496:     if (patch->precomputeElementTensors) {
2497:       PCPatchPrecomputePatchTensors_Private(pc);
2498:     }
2499:     for (i = 0; i < patch->npatch; ++i) {
2500:       PCPatchComputeOperator_Internal(pc, NULL, patch->mat[i], i, PETSC_FALSE);
2501:       if (!patch->denseinverse) {
2502:         KSPSetOperators((KSP) patch->solver[i], patch->mat[i], patch->mat[i]);
2503:       } else if (patch->mat[i] && !patch->densesolve) {
2504:         /* Setup matmult callback */
2505:         MatGetOperation(patch->mat[i], MATOP_MULT, (void (**)(void))&patch->densesolve);
2506:       }
2507:     }
2508:   }
2509:   if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2510:     for (i = 0; i < patch->npatch; ++i) {
2511:       /* Instead of padding patch->patchUpdate with zeros to get */
2512:       /* patch->patchUpdateWithArtificial and then multiplying with the matrix, */
2513:       /* just get rid of the columns that correspond to the dofs with */
2514:       /* artificial bcs. That's of course fairly inefficient, hopefully we */
2515:       /* can just assemble the rectangular matrix in the first place. */
2516:       Mat matSquare;
2517:       IS rowis;
2518:       PetscInt dof;

2520:       MatGetSize(patch->mat[i], &dof, NULL);
2521:       if (dof == 0) {
2522:         patch->matWithArtificial[i] = NULL;
2523:         continue;
2524:       }

2526:       PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);
2527:       PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);

2529:       MatGetSize(matSquare, &dof, NULL);
2530:       ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis); 
2531:       if (pc->setupcalled) {
2532:         MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_REUSE_MATRIX, &patch->matWithArtificial[i]); 
2533:       } else {
2534:         MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &patch->matWithArtificial[i]); 
2535:       }
2536:       ISDestroy(&rowis); 
2537:       MatDestroy(&matSquare);
2538:     }
2539:   }
2540:   return(0);
2541: }

2543: static PetscErrorCode PCSetUp_PATCH(PC pc)
2544: {
2545:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2546:   PetscInt       i;
2547:   PetscBool      isNonlinear;
2548:   PetscInt       maxDof = -1, maxDofWithArtificial = -1;

2552:   if (!pc->setupcalled) {
2553:     PetscInt pStart, pEnd, p;
2554:     PetscInt localSize;

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

2558:     isNonlinear = patch->isNonlinear;
2559:     if (!patch->nsubspaces) {
2560:       DM           dm, plex;
2561:       PetscSection s;
2562:       PetscInt     cStart, cEnd, c, Nf, f, numGlobalBcs = 0, *globalBcs, *Nb, totNb = 0, **cellDofs;

2564:       PCGetDM(pc, &dm);
2565:       if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "Must set DM for PCPATCH or call PCPatchSetDiscretisationInfo()");
2566:       DMConvert(dm, DMPLEX, &plex);
2567:       dm = plex;
2568:       DMGetLocalSection(dm, &s);
2569:       PetscSectionGetNumFields(s, &Nf);
2570:       PetscSectionGetChart(s, &pStart, &pEnd);
2571:       for (p = pStart; p < pEnd; ++p) {
2572:         PetscInt cdof;
2573:         PetscSectionGetConstraintDof(s, p, &cdof);
2574:         numGlobalBcs += cdof;
2575:       }
2576:       DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2577:       PetscMalloc3(Nf, &Nb, Nf, &cellDofs, numGlobalBcs, &globalBcs);
2578:       for (f = 0; f < Nf; ++f) {
2579:         PetscFE        fe;
2580:         PetscDualSpace sp;
2581:         PetscInt       cdoff = 0;

2583:         DMGetField(dm, f, NULL, (PetscObject *) &fe);
2584:         /* PetscFEGetNumComponents(fe, &Nc[f]); */
2585:         PetscFEGetDualSpace(fe, &sp);
2586:         PetscDualSpaceGetDimension(sp, &Nb[f]);
2587:         totNb += Nb[f];

2589:         PetscMalloc1((cEnd-cStart)*Nb[f], &cellDofs[f]);
2590:         for (c = cStart; c < cEnd; ++c) {
2591:           PetscInt *closure = NULL;
2592:           PetscInt  clSize  = 0, cl;

2594:           DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
2595:           for (cl = 0; cl < clSize*2; cl += 2) {
2596:             const PetscInt p = closure[cl];
2597:             PetscInt       fdof, d, foff;

2599:             PetscSectionGetFieldDof(s, p, f, &fdof);
2600:             PetscSectionGetFieldOffset(s, p, f, &foff);
2601:             for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d;
2602:           }
2603:           DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
2604:         }
2605:         if (cdoff != (cEnd-cStart)*Nb[f]) SETERRQ4(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_SIZ, "Total number of cellDofs %D for field %D should be Nc (%D) * cellDof (%D)", cdoff, f, cEnd-cStart, Nb[f]);
2606:       }
2607:       numGlobalBcs = 0;
2608:       for (p = pStart; p < pEnd; ++p) {
2609:         const PetscInt *ind;
2610:         PetscInt        off, cdof, d;

2612:         PetscSectionGetOffset(s, p, &off);
2613:         PetscSectionGetConstraintDof(s, p, &cdof);
2614:         PetscSectionGetConstraintIndices(s, p, &ind);
2615:         for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d];
2616:       }

2618:       PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **) cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs);
2619:       for (f = 0; f < Nf; ++f) {
2620:         PetscFree(cellDofs[f]);
2621:       }
2622:       PetscFree3(Nb, cellDofs, globalBcs);
2623:       PCPatchSetComputeFunction(pc, PCPatchComputeFunction_DMPlex_Private, NULL);
2624:       PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL);
2625:       DMDestroy(&dm);
2626:     }

2628:     localSize = patch->subspaceOffsets[patch->nsubspaces];
2629:     VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localRHS);
2630:     VecSetUp(patch->localRHS);
2631:     VecDuplicate(patch->localRHS, &patch->localUpdate);
2632:     PCPatchCreateCellPatches(pc);
2633:     PCPatchCreateCellPatchDiscretisationInfo(pc);

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

2638:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2639:       PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithArtificial);
2640:     }
2641:     if (isNonlinear) {
2642:       PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithAll);
2643:     }
2644:     for (p = pStart; p < pEnd; ++p) {
2645:       PetscInt dof;

2647:       PetscSectionGetDof(patch->gtolCounts, p, &dof);
2648:       maxDof = PetscMax(maxDof, dof);
2649:       if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2650:         const PetscInt    *gtolArray, *gtolArrayWithArtificial = NULL;
2651:         PetscInt           numPatchDofs, offset;
2652:         PetscInt           numPatchDofsWithArtificial, offsetWithArtificial;
2653:         PetscInt           dofWithoutArtificialCounter = 0;
2654:         PetscInt          *patchWithoutArtificialToWithArtificialArray;

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

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

2663:         PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);
2664:         if (numPatchDofs == 0) {
2665:           patch->dofMappingWithoutToWithArtificial[p-pStart] = NULL;
2666:           continue;
2667:         }

2669:         PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2670:         ISGetIndices(patch->gtolWithArtificial, &gtolArrayWithArtificial);
2671:         PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &numPatchDofsWithArtificial);
2672:         PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offsetWithArtificial);

2674:         PetscMalloc1(numPatchDofs, &patchWithoutArtificialToWithArtificialArray);
2675:         for (i=0; i<numPatchDofsWithArtificial; i++) {
2676:           if (gtolArrayWithArtificial[i+offsetWithArtificial] == gtolArray[offset+dofWithoutArtificialCounter]) {
2677:             patchWithoutArtificialToWithArtificialArray[dofWithoutArtificialCounter] = i;
2678:             dofWithoutArtificialCounter++;
2679:             if (dofWithoutArtificialCounter == numPatchDofs)
2680:               break;
2681:           }
2682:         }
2683:         ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutArtificialToWithArtificialArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithArtificial[p-pStart]);
2684:         ISRestoreIndices(patch->gtol, &gtolArray);
2685:         ISRestoreIndices(patch->gtolWithArtificial, &gtolArrayWithArtificial);
2686:       }
2687:       if (isNonlinear) {
2688:         const PetscInt    *gtolArray, *gtolArrayWithAll = NULL;
2689:         PetscInt           numPatchDofs, offset;
2690:         PetscInt           numPatchDofsWithAll, offsetWithAll;
2691:         PetscInt           dofWithoutAllCounter = 0;
2692:         PetscInt          *patchWithoutAllToWithAllArray;

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

2698:         PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);
2699:         if (numPatchDofs == 0) {
2700:           patch->dofMappingWithoutToWithAll[p-pStart] = NULL;
2701:           continue;
2702:         }

2704:         PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2705:         ISGetIndices(patch->gtolWithAll, &gtolArrayWithAll);
2706:         PetscSectionGetDof(patch->gtolCountsWithAll, p, &numPatchDofsWithAll);
2707:         PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offsetWithAll);

2709:         PetscMalloc1(numPatchDofs, &patchWithoutAllToWithAllArray);

2711:         for (i=0; i<numPatchDofsWithAll; i++) {
2712:           if (gtolArrayWithAll[i+offsetWithAll] == gtolArray[offset+dofWithoutAllCounter]) {
2713:             patchWithoutAllToWithAllArray[dofWithoutAllCounter] = i;
2714:             dofWithoutAllCounter++;
2715:             if (dofWithoutAllCounter == numPatchDofs)
2716:               break;
2717:           }
2718:         }
2719:         ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutAllToWithAllArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithAll[p-pStart]);
2720:         ISRestoreIndices(patch->gtol, &gtolArray);
2721:         ISRestoreIndices(patch->gtolWithAll, &gtolArrayWithAll);
2722:       }
2723:     }
2724:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2725:       VecCreateSeq(PETSC_COMM_SELF, maxDofWithArtificial, &patch->patchRHSWithArtificial);
2726:       VecSetUp(patch->patchRHSWithArtificial);
2727:     }
2728:     VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchRHS);
2729:     VecSetUp(patch->patchRHS);
2730:     VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchUpdate);
2731:     VecSetUp(patch->patchUpdate);
2732:     if (patch->save_operators) {
2733:       PetscMalloc1(patch->npatch, &patch->mat);
2734:       for (i = 0; i < patch->npatch; ++i) {
2735:         PCPatchCreateMatrix_Private(pc, i, &patch->mat[i], PETSC_FALSE);
2736:       }
2737:     }
2738:     PetscLogEventEnd(PC_Patch_CreatePatches, pc, 0, 0, 0);

2740:     /* If desired, calculate weights for dof multiplicity */
2741:     if (patch->partition_of_unity) {
2742:       PetscScalar *input = NULL;
2743:       PetscScalar *output = NULL;
2744:       Vec global;

2746:       VecDuplicate(patch->localRHS, &patch->dof_weights);
2747:       if(patch->local_composition_type == PC_COMPOSITE_ADDITIVE) {
2748:         for (i = 0; i < patch->npatch; ++i) {
2749:           PetscInt dof;

2751:           PetscSectionGetDof(patch->gtolCounts, i+pStart, &dof);
2752:           if (dof <= 0) continue;
2753:           VecSet(patch->patchRHS, 1.0);
2754:           PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchRHS, patch->dof_weights, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);
2755:         }
2756:       } else {
2757:         /* multiplicative is actually only locally multiplicative and globally additive. need the pou where the mesh decomposition overlaps */
2758:         VecSet(patch->dof_weights, 1.0);
2759:       }

2761:       VecDuplicate(patch->dof_weights, &global);
2762:       VecSet(global, 0.);

2764:       VecGetArray(patch->dof_weights, &input);
2765:       VecGetArray(global, &output);
2766:       PetscSFReduceBegin(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM);
2767:       PetscSFReduceEnd(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM);
2768:       VecRestoreArray(patch->dof_weights, &input);
2769:       VecRestoreArray(global, &output);

2771:       VecReciprocal(global);

2773:       VecGetArray(patch->dof_weights, &output);
2774:       VecGetArray(global, &input);
2775:       PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, input, output);
2776:       PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, input, output);
2777:       VecRestoreArray(patch->dof_weights, &output);
2778:       VecRestoreArray(global, &input);
2779:       VecDestroy(&global);
2780:     }
2781:     if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE && patch->save_operators) {
2782:       PetscMalloc1(patch->npatch, &patch->matWithArtificial);
2783:     }
2784:   }
2785:   (*patch->setupsolver)(pc);
2786:   return(0);
2787: }

2789: static PetscErrorCode PCApply_PATCH_Linear(PC pc, PetscInt i, Vec x, Vec y)
2790: {
2791:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2792:   KSP            ksp;
2793:   Mat            op;
2794:   PetscInt       m, n;

2798:   if (patch->denseinverse) {
2799:     (*patch->densesolve)(patch->mat[i], x, y);
2800:     return(0);
2801:   }
2802:   ksp = (KSP) patch->solver[i];
2803:   if (!patch->save_operators) {
2804:     Mat mat;

2806:     PCPatchCreateMatrix_Private(pc, i, &mat, PETSC_FALSE);
2807:     /* Populate operator here. */
2808:     PCPatchComputeOperator_Internal(pc, NULL, mat, i, PETSC_FALSE);
2809:     KSPSetOperators(ksp, mat, mat);
2810:     /* Drop reference so the KSPSetOperators below will blow it away. */
2811:     MatDestroy(&mat);
2812:   }
2813:   PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);
2814:   if (!ksp->setfromoptionscalled) {
2815:     KSPSetFromOptions(ksp);
2816:   }
2817:   /* Disgusting trick to reuse work vectors */
2818:   KSPGetOperators(ksp, &op, NULL);
2819:   MatGetLocalSize(op, &m, &n);
2820:   x->map->n = m;
2821:   y->map->n = n;
2822:   x->map->N = m;
2823:   y->map->N = n;
2824:   KSPSolve(ksp, x, y);
2825:   KSPCheckSolve(ksp, pc, y);
2826:   PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);
2827:   if (!patch->save_operators) {
2828:     PC pc;
2829:     KSPSetOperators(ksp, NULL, NULL);
2830:     KSPGetPC(ksp, &pc);
2831:     /* Destroy PC context too, otherwise the factored matrix hangs around. */
2832:     PCReset(pc);
2833:   }
2834:   return(0);
2835: }

2837: static PetscErrorCode PCUpdateMultiplicative_PATCH_Linear(PC pc, PetscInt i, PetscInt pStart)
2838: {
2839:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2840:   Mat multMat;
2841:   PetscInt n, m;


2846:   if (patch->save_operators) {
2847:     multMat = patch->matWithArtificial[i];
2848:   } else {
2849:     /*Very inefficient, hopefully we can just assemble the rectangular matrix in the first place.*/
2850:     Mat matSquare;
2851:     PetscInt dof;
2852:     IS rowis;
2853:     PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);
2854:     PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);
2855:     MatGetSize(matSquare, &dof, NULL);
2856:     ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis); 
2857:     MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &multMat); 
2858:     MatDestroy(&matSquare);
2859:     ISDestroy(&rowis); 
2860:   }
2861:   /* Disgusting trick to reuse work vectors */
2862:   MatGetLocalSize(multMat, &m, &n);
2863:   patch->patchUpdate->map->n = n;
2864:   patch->patchRHSWithArtificial->map->n = m;
2865:   patch->patchUpdate->map->N = n;
2866:   patch->patchRHSWithArtificial->map->N = m;
2867:   MatMult(multMat, patch->patchUpdate, patch->patchRHSWithArtificial); 
2868:   VecScale(patch->patchRHSWithArtificial, -1.0); 
2869:   PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHSWithArtificial, patch->localRHS, ADD_VALUES, SCATTER_REVERSE, SCATTER_WITHARTIFICIAL); 
2870:   if (!patch->save_operators) {
2871:     MatDestroy(&multMat); 
2872:   }
2873:   return(0);
2874: }

2876: static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y)
2877: {
2878:   PC_PATCH          *patch    = (PC_PATCH *) pc->data;
2879:   const PetscScalar *globalRHS  = NULL;
2880:   PetscScalar       *localRHS   = NULL;
2881:   PetscScalar       *globalUpdate  = NULL;
2882:   const PetscInt    *bcNodes  = NULL;
2883:   PetscInt           nsweep   = patch->symmetrise_sweep ? 2 : 1;
2884:   PetscInt           start[2] = {0, 0};
2885:   PetscInt           end[2]   = {-1, -1};
2886:   const PetscInt     inc[2]   = {1, -1};
2887:   const PetscScalar *localUpdate;
2888:   const PetscInt    *iterationSet;
2889:   PetscInt           pStart, numBcs, n, sweep, bc, j;
2890:   PetscErrorCode     ierr;

2893:   PetscLogEventBegin(PC_Patch_Apply, pc, 0, 0, 0);
2894:   PetscOptionsPushGetViewerOff(PETSC_TRUE);
2895:   /* start, end, inc have 2 entries to manage a second backward sweep if we symmetrize */
2896:   end[0]   = patch->npatch;
2897:   start[1] = patch->npatch-1;
2898:   if (patch->user_patches) {
2899:     ISGetLocalSize(patch->iterationSet, &end[0]);
2900:     start[1] = end[0] - 1;
2901:     ISGetIndices(patch->iterationSet, &iterationSet);
2902:   }
2903:   /* Scatter from global space into overlapped local spaces */
2904:   VecGetArrayRead(x, &globalRHS);
2905:   VecGetArray(patch->localRHS, &localRHS);
2906:   PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, globalRHS, localRHS);
2907:   PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, globalRHS, localRHS);
2908:   VecRestoreArrayRead(x, &globalRHS);
2909:   VecRestoreArray(patch->localRHS, &localRHS);

2911:   VecSet(patch->localUpdate, 0.0);
2912:   PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);
2913:   PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);
2914:   for (sweep = 0; sweep < nsweep; sweep++) {
2915:     for (j = start[sweep]; j*inc[sweep] < end[sweep]*inc[sweep]; j += inc[sweep]) {
2916:       PetscInt i       = patch->user_patches ? iterationSet[j] : j;
2917:       PetscInt start, len;

2919:       PetscSectionGetDof(patch->gtolCounts, i+pStart, &len);
2920:       PetscSectionGetOffset(patch->gtolCounts, i+pStart, &start);
2921:       /* TODO: Squash out these guys in the setup as well. */
2922:       if (len <= 0) continue;
2923:       /* TODO: Do we need different scatters for X and Y? */
2924:       PCPatch_ScatterLocal_Private(pc, i+pStart, patch->localRHS, patch->patchRHS, INSERT_VALUES, SCATTER_FORWARD, SCATTER_INTERIOR);
2925:       (*patch->applysolver)(pc, i, patch->patchRHS, patch->patchUpdate);
2926:       PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchUpdate, patch->localUpdate, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);
2927:       if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2928:         (*patch->updatemultiplicative)(pc, i, pStart);
2929:       }
2930:     }
2931:   }
2932:   PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);
2933:   if (patch->user_patches) {ISRestoreIndices(patch->iterationSet, &iterationSet);}
2934:   /* XXX: should we do this on the global vector? */
2935:   if (patch->partition_of_unity) {
2936:     VecPointwiseMult(patch->localUpdate, patch->localUpdate, patch->dof_weights);
2937:   }
2938:   /* Now patch->localUpdate contains the solution of the patch solves, so we need to combine them all. */
2939:   VecSet(y, 0.0);
2940:   VecGetArray(y, &globalUpdate);
2941:   VecGetArrayRead(patch->localUpdate, &localUpdate);
2942:   PetscSFReduceBegin(patch->sectionSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);
2943:   PetscSFReduceEnd(patch->sectionSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);
2944:   VecRestoreArrayRead(patch->localUpdate, &localUpdate);

2946:   /* Now we need to send the global BC values through */
2947:   VecGetArrayRead(x, &globalRHS);
2948:   ISGetSize(patch->globalBcNodes, &numBcs);
2949:   ISGetIndices(patch->globalBcNodes, &bcNodes);
2950:   VecGetLocalSize(x, &n);
2951:   for (bc = 0; bc < numBcs; ++bc) {
2952:     const PetscInt idx = bcNodes[bc];
2953:     if (idx < n) globalUpdate[idx] = globalRHS[idx];
2954:   }

2956:   ISRestoreIndices(patch->globalBcNodes, &bcNodes);
2957:   VecRestoreArrayRead(x, &globalRHS);
2958:   VecRestoreArray(y, &globalUpdate);

2960:   PetscOptionsPopGetViewerOff();
2961:   PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0);
2962:   return(0);
2963: }

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

2972:   if (patch->solver) {
2973:     for (i = 0; i < patch->npatch; ++i) {KSPReset((KSP) patch->solver[i]);}
2974:   }
2975:   return(0);
2976: }

2978: static PetscErrorCode PCReset_PATCH(PC pc)
2979: {
2980:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2981:   PetscInt       i;


2986:   PetscSFDestroy(&patch->sectionSF);
2987:   PetscSectionDestroy(&patch->cellCounts);
2988:   PetscSectionDestroy(&patch->pointCounts);
2989:   PetscSectionDestroy(&patch->cellNumbering);
2990:   PetscSectionDestroy(&patch->gtolCounts);
2991:   ISDestroy(&patch->gtol);
2992:   ISDestroy(&patch->cells);
2993:   ISDestroy(&patch->points);
2994:   ISDestroy(&patch->dofs);
2995:   ISDestroy(&patch->offs);
2996:   PetscSectionDestroy(&patch->patchSection);
2997:   ISDestroy(&patch->ghostBcNodes);
2998:   ISDestroy(&patch->globalBcNodes);
2999:   PetscSectionDestroy(&patch->gtolCountsWithArtificial);
3000:   ISDestroy(&patch->gtolWithArtificial);
3001:   ISDestroy(&patch->dofsWithArtificial);
3002:   ISDestroy(&patch->offsWithArtificial);
3003:   PetscSectionDestroy(&patch->gtolCountsWithAll);
3004:   ISDestroy(&patch->gtolWithAll);
3005:   ISDestroy(&patch->dofsWithAll);
3006:   ISDestroy(&patch->offsWithAll);
3007:   VecDestroy(&patch->cellMats);
3008:   VecDestroy(&patch->intFacetMats);
3009:   ISDestroy(&patch->allCells);
3010:   ISDestroy(&patch->intFacets);
3011:   ISDestroy(&patch->extFacets);
3012:   ISDestroy(&patch->intFacetsToPatchCell);
3013:   PetscSectionDestroy(&patch->intFacetCounts);
3014:   PetscSectionDestroy(&patch->extFacetCounts);

3016:   if (patch->dofSection) for (i = 0; i < patch->nsubspaces; i++) {PetscSectionDestroy(&patch->dofSection[i]);}
3017:   PetscFree(patch->dofSection);
3018:   PetscFree(patch->bs);
3019:   PetscFree(patch->nodesPerCell);
3020:   if (patch->cellNodeMap) for (i = 0; i < patch->nsubspaces; i++) {PetscFree(patch->cellNodeMap[i]);}
3021:   PetscFree(patch->cellNodeMap);
3022:   PetscFree(patch->subspaceOffsets);

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

3026:   if (patch->subspaces_to_exclude) {
3027:     PetscHSetIDestroy(&patch->subspaces_to_exclude);
3028:   }

3030:   VecDestroy(&patch->localRHS);
3031:   VecDestroy(&patch->localUpdate);
3032:   VecDestroy(&patch->patchRHS);
3033:   VecDestroy(&patch->patchUpdate);
3034:   VecDestroy(&patch->dof_weights);
3035:   if (patch->patch_dof_weights) {
3036:     for (i = 0; i < patch->npatch; ++i) {VecDestroy(&patch->patch_dof_weights[i]);}
3037:     PetscFree(patch->patch_dof_weights);
3038:   }
3039:   if (patch->mat) {
3040:     for (i = 0; i < patch->npatch; ++i) {MatDestroy(&patch->mat[i]);}
3041:     PetscFree(patch->mat);
3042:   }
3043:   if (patch->matWithArtificial) {
3044:     for (i = 0; i < patch->npatch; ++i) {MatDestroy(&patch->matWithArtificial[i]);}
3045:     PetscFree(patch->matWithArtificial);
3046:   }
3047:   VecDestroy(&patch->patchRHSWithArtificial);
3048:   if(patch->dofMappingWithoutToWithArtificial) {
3049:     for (i = 0; i < patch->npatch; ++i) {ISDestroy(&patch->dofMappingWithoutToWithArtificial[i]);}
3050:     PetscFree(patch->dofMappingWithoutToWithArtificial);

3052:   }
3053:   if(patch->dofMappingWithoutToWithAll) {
3054:     for (i = 0; i < patch->npatch; ++i) {ISDestroy(&patch->dofMappingWithoutToWithAll[i]);}
3055:     PetscFree(patch->dofMappingWithoutToWithAll);

3057:   }
3058:   PetscFree(patch->sub_mat_type);
3059:   if (patch->userIS) {
3060:     for (i = 0; i < patch->npatch; ++i) {ISDestroy(&patch->userIS[i]);}
3061:     PetscFree(patch->userIS);
3062:   }
3063:   PetscFree(patch->precomputedTensorLocations);
3064:   PetscFree(patch->precomputedIntFacetTensorLocations);
3065:   
3066:   patch->bs          = 0;
3067:   patch->cellNodeMap = NULL;
3068:   patch->nsubspaces  = 0;
3069:   ISDestroy(&patch->iterationSet);

3071:   PetscViewerDestroy(&patch->viewerSection);
3072:   return(0);
3073: }

3075: static PetscErrorCode PCDestroy_PATCH_Linear(PC pc)
3076: {
3077:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
3078:   PetscInt       i;

3082:   if (patch->solver) {
3083:     for (i = 0; i < patch->npatch; ++i) {KSPDestroy((KSP *) &patch->solver[i]);}
3084:     PetscFree(patch->solver);
3085:   }
3086:   return(0);
3087: }

3089: static PetscErrorCode PCDestroy_PATCH(PC pc)
3090: {
3091:   PC_PATCH      *patch = (PC_PATCH *) pc->data;

3095:   PCReset_PATCH(pc);
3096:   (*patch->destroysolver)(pc);
3097:   PetscFree(pc->data);
3098:   return(0);
3099: }

3101: static PetscErrorCode PCSetFromOptions_PATCH(PetscOptionItems *PetscOptionsObject, PC pc)
3102: {
3103:   PC_PATCH            *patch = (PC_PATCH *) pc->data;
3104:   PCPatchConstructType patchConstructionType = PC_PATCH_STAR;
3105:   char                 sub_mat_type[PETSC_MAX_PATH_LEN];
3106:   char                 option[PETSC_MAX_PATH_LEN];
3107:   const char          *prefix;
3108:   PetscBool            flg, dimflg, codimflg;
3109:   MPI_Comm             comm;
3110:   PetscInt            *ifields, nfields, k;
3111:   PetscErrorCode       ierr;
3112:   PCCompositeType loctype = PC_COMPOSITE_ADDITIVE;

3115:   PetscObjectGetComm((PetscObject) pc, &comm);
3116:   PetscObjectGetOptionsPrefix((PetscObject) pc, &prefix);
3117:   PetscOptionsHead(PetscOptionsObject, "Patch solver options");

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

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

3127:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_local_type", patch->classname);
3128:   PetscOptionsEnum(option,"Type of local solver composition (additive or multiplicative)","PCPatchSetLocalComposition",PCCompositeTypes,(PetscEnum)loctype,(PetscEnum*)&loctype,&flg);
3129:   if(flg) { PCPatchSetLocalComposition(pc, loctype);}
3130:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_dense_inverse", patch->classname);
3131:   PetscOptionsBool(option, "Compute inverses of patch matrices and apply directly? Ignores KSP/PC settings on patch.", "PCPatchSetDenseInverse", patch->denseinverse, &patch->denseinverse, &flg);
3132:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_dim", patch->classname);
3133:   PetscOptionsInt(option, "What dimension of mesh point to construct patches by? (0 = vertices)", "PCPATCH", patch->dim, &patch->dim, &dimflg);
3134:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_codim", patch->classname);
3135:   PetscOptionsInt(option, "What co-dimension of mesh point to construct patches by? (0 = cells)", "PCPATCH", patch->codim, &patch->codim, &codimflg);
3136:   if (dimflg && codimflg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Can only set one of dimension or co-dimension");

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

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

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

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

3151:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_sub_mat_type", patch->classname);
3152:   PetscOptionsFList(option, "Matrix type for patch solves", "PCPatchSetSubMatType", MatList, NULL, sub_mat_type, PETSC_MAX_PATH_LEN, &flg);
3153:   if (flg) {PCPatchSetSubMatType(pc, sub_mat_type);}

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

3158:   /* If the user has set the number of subspaces, use that for the buffer size,
3159:      otherwise use a large number */
3160:   if (patch->nsubspaces <= 0) {
3161:     nfields = 128;
3162:   } else {
3163:     nfields = patch->nsubspaces;
3164:   }
3165:   PetscMalloc1(nfields, &ifields);
3166:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exclude_subspaces", patch->classname);
3167:   PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,option,ifields,&nfields,&flg);
3168:   if (flg && (patchConstructionType == PC_PATCH_USER)) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "We cannot support excluding a subspace with user patches because we do not index patches with a mesh point");
3169:   if (flg) {
3170:     PetscHSetIClear(patch->subspaces_to_exclude);
3171:     for (k = 0; k < nfields; k++) {
3172:       PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]);
3173:     }
3174:   }
3175:   PetscFree(ifields);

3177:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_patches_view", patch->classname);
3178:   PetscOptionsBool(option, "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg);
3179:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_cells_view", patch->classname);
3180:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerCells, &patch->formatCells, &patch->viewCells);
3181:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_interior_facets_view", patch->classname);
3182:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerIntFacets, &patch->formatIntFacets, &patch->viewIntFacets);
3183:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exterior_facets_view", patch->classname);
3184:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerExtFacets, &patch->formatExtFacets, &patch->viewExtFacets);
3185:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_points_view", patch->classname);
3186:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints);
3187:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_section_view", patch->classname);
3188:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerSection, &patch->formatSection, &patch->viewSection);
3189:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_mat_view", patch->classname);
3190:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix);
3191:   PetscOptionsTail();
3192:   patch->optionsSet = PETSC_TRUE;
3193:   return(0);
3194: }

3196: static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc)
3197: {
3198:   PC_PATCH          *patch = (PC_PATCH*) pc->data;
3199:   KSPConvergedReason reason;
3200:   PetscInt           i;
3201:   PetscErrorCode     ierr;

3204:   if (!patch->save_operators) {
3205:     /* Can't do this here because the sub KSPs don't have an operator attached yet. */
3206:     return(0);
3207:   }
3208:   if (patch->denseinverse) {
3209:     /* No solvers */
3210:     return(0);
3211:   }
3212:   for (i = 0; i < patch->npatch; ++i) {
3213:     if (!((KSP) patch->solver[i])->setfromoptionscalled) {
3214:       KSPSetFromOptions((KSP) patch->solver[i]);
3215:     }
3216:     KSPSetUp((KSP) patch->solver[i]);
3217:     KSPGetConvergedReason((KSP) patch->solver[i], &reason);
3218:     if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
3219:   }
3220:   return(0);
3221: }

3223: static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer)
3224: {
3225:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
3226:   PetscViewer    sviewer;
3227:   PetscBool      isascii;
3228:   PetscMPIInt    rank;

3232:   /* TODO Redo tabbing with set tbas in new style */
3233:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
3234:   if (!isascii) return(0);
3235:   MPI_Comm_rank(PetscObjectComm((PetscObject) pc), &rank);
3236:   PetscViewerASCIIPushTab(viewer);
3237:   PetscViewerASCIIPrintf(viewer, "Subspace Correction preconditioner with %d patches\n", patch->npatch);
3238:   if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
3239:     PetscViewerASCIIPrintf(viewer, "Schwarz type: multiplicative\n");
3240:   } else {
3241:     PetscViewerASCIIPrintf(viewer, "Schwarz type: additive\n");
3242:   }
3243:   if (patch->partition_of_unity) {PetscViewerASCIIPrintf(viewer, "Weighting by partition of unity\n");}
3244:   else                           {PetscViewerASCIIPrintf(viewer, "Not weighting by partition of unity\n");}
3245:   if (patch->symmetrise_sweep) {PetscViewerASCIIPrintf(viewer, "Symmetrising sweep (start->end, then end->start)\n");}
3246:   else                         {PetscViewerASCIIPrintf(viewer, "Not symmetrising sweep\n");}
3247:   if (!patch->precomputeElementTensors) {PetscViewerASCIIPrintf(viewer, "Not precomputing element tensors (overlapping cells rebuilt in every patch assembly)\n");}
3248:   else                            {PetscViewerASCIIPrintf(viewer, "Precomputing element tensors (each cell assembled only once)\n");}
3249:   if (!patch->save_operators) {PetscViewerASCIIPrintf(viewer, "Not saving patch operators (rebuilt every PCApply)\n");}
3250:   else                        {PetscViewerASCIIPrintf(viewer, "Saving patch operators (rebuilt every PCSetUp)\n");}
3251:   if (patch->patchconstructop == PCPatchConstruct_Star)       {PetscViewerASCIIPrintf(viewer, "Patch construction operator: star\n");}
3252:   else if (patch->patchconstructop == PCPatchConstruct_Vanka) {PetscViewerASCIIPrintf(viewer, "Patch construction operator: Vanka\n");}
3253:   else if (patch->patchconstructop == PCPatchConstruct_User)  {PetscViewerASCIIPrintf(viewer, "Patch construction operator: user-specified\n");}
3254:   else                                                        {PetscViewerASCIIPrintf(viewer, "Patch construction operator: unknown\n");}

3256:   if (patch->denseinverse) {
3257:     PetscViewerASCIIPrintf(viewer, "Explicitly forming dense inverse and applying patch solver via MatMult.\n");
3258:   } else {
3259:     if (patch->isNonlinear) {
3260:       PetscViewerASCIIPrintf(viewer, "SNES on patches (all same):\n");
3261:     } else {
3262:       PetscViewerASCIIPrintf(viewer, "KSP on patches (all same):\n");
3263:     }
3264:     if (patch->solver) {
3265:       PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
3266:       if (!rank) {
3267:         PetscViewerASCIIPushTab(sviewer);
3268:         PetscObjectView(patch->solver[0], sviewer);
3269:         PetscViewerASCIIPopTab(sviewer);
3270:       }
3271:       PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
3272:     } else {
3273:       PetscViewerASCIIPushTab(viewer);
3274:       PetscViewerASCIIPrintf(viewer, "Solver not yet set.\n");
3275:       PetscViewerASCIIPopTab(viewer);
3276:     }
3277:   }
3278:   PetscViewerASCIIPopTab(viewer);
3279:   return(0);
3280: }

3282: /*MC
3283:   PCPATCH - A PC object that encapsulates flexible definition of blocks for overlapping and non-overlapping
3284:             small block additive preconditioners. Block definition is based on topology from
3285:             a DM and equation numbering from a PetscSection.

3287:   Options Database Keys:
3288: + -pc_patch_cells_view   - Views the process local cell numbers for each patch
3289: . -pc_patch_points_view  - Views the process local mesh point numbers for each patch
3290: . -pc_patch_g2l_view     - Views the map between global dofs and patch local dofs for each patch
3291: . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary
3292: - -pc_patch_sub_mat_view - Views the matrix associated with each patch

3294:   Level: intermediate

3296: .seealso: PCType, PCCreate(), PCSetType()
3297: M*/
3298: PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc)
3299: {
3300:   PC_PATCH      *patch;

3304:   PetscNewLog(pc, &patch);

3306:   if (patch->subspaces_to_exclude) {
3307:     PetscHSetIDestroy(&patch->subspaces_to_exclude);
3308:   }
3309:   PetscHSetICreate(&patch->subspaces_to_exclude);

3311:   patch->classname = "pc";
3312:   patch->isNonlinear = PETSC_FALSE;

3314:   /* Set some defaults */
3315:   patch->combined           = PETSC_FALSE;
3316:   patch->save_operators     = PETSC_TRUE;
3317:   patch->local_composition_type = PC_COMPOSITE_ADDITIVE;
3318:   patch->precomputeElementTensors = PETSC_FALSE;
3319:   patch->partition_of_unity = PETSC_FALSE;
3320:   patch->codim              = -1;
3321:   patch->dim                = -1;
3322:   patch->vankadim           = -1;
3323:   patch->ignoredim          = -1;
3324:   patch->pardecomp_overlap  = 0;
3325:   patch->patchconstructop   = PCPatchConstruct_Star;
3326:   patch->symmetrise_sweep   = PETSC_FALSE;
3327:   patch->npatch             = 0;
3328:   patch->userIS             = NULL;
3329:   patch->optionsSet         = PETSC_FALSE;
3330:   patch->iterationSet       = NULL;
3331:   patch->user_patches       = PETSC_FALSE;
3332:   PetscStrallocpy(MATDENSE, (char **) &patch->sub_mat_type);
3333:   patch->viewPatches        = PETSC_FALSE;
3334:   patch->viewCells          = PETSC_FALSE;
3335:   patch->viewPoints         = PETSC_FALSE;
3336:   patch->viewSection        = PETSC_FALSE;
3337:   patch->viewMatrix         = PETSC_FALSE;
3338:   patch->densesolve         = NULL;
3339:   patch->setupsolver        = PCSetUp_PATCH_Linear;
3340:   patch->applysolver        = PCApply_PATCH_Linear;
3341:   patch->resetsolver        = PCReset_PATCH_Linear;
3342:   patch->destroysolver      = PCDestroy_PATCH_Linear;
3343:   patch->updatemultiplicative = PCUpdateMultiplicative_PATCH_Linear;
3344:   patch->dofMappingWithoutToWithArtificial = NULL;
3345:   patch->dofMappingWithoutToWithAll = NULL;

3347:   pc->data                 = (void *) patch;
3348:   pc->ops->apply           = PCApply_PATCH;
3349:   pc->ops->applytranspose  = 0; /* PCApplyTranspose_PATCH; */
3350:   pc->ops->setup           = PCSetUp_PATCH;
3351:   pc->ops->reset           = PCReset_PATCH;
3352:   pc->ops->destroy         = PCDestroy_PATCH;
3353:   pc->ops->setfromoptions  = PCSetFromOptions_PATCH;
3354:   pc->ops->setuponblocks   = PCSetUpOnBlocks_PATCH;
3355:   pc->ops->view            = PCView_PATCH;
3356:   pc->ops->applyrichardson = 0;

3358:   return(0);
3359: }