Actual source code: pcpatch.c

petsc-master 2019-07-20
Report Typos and Errors
  1:  #include <petsc/private/pcpatchimpl.h>
  2: #include <petsc/private/kspimpl.h>         /* For ksp->setfromoptionscalled */
  3: #include <petsc/private/dmpleximpl.h> /* For DMPlexComputeJacobian_Patch_Internal() */
  4:  #include <petscsf.h>
  5:  #include <petscbt.h>
  6:  #include <petscds.h>

  8: PetscLogEvent PC_Patch_CreatePatches, PC_Patch_ComputeOp, PC_Patch_Solve, PC_Patch_Scatter, PC_Patch_Apply, PC_Patch_Prealloc;

 10: PETSC_STATIC_INLINE PetscErrorCode ObjectView(PetscObject obj, PetscViewer viewer, PetscViewerFormat format)
 11: {

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

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

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

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

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

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

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

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

 93:   PetscHSetIClear(ht);

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

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

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

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

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

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

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

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

150:   return(0);
151: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

377: /* TODO: Docs */
378: PetscErrorCode PCPatchSetLocalComposition(PC pc, PCCompositeType type)
379: {
380:   PC_PATCH *patch = (PC_PATCH *) pc->data;
382:   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");
383:   patch->local_composition_type = type;
384:   return(0);
385: }

387: /* TODO: Docs */
388: PetscErrorCode PCPatchGetLocalComposition(PC pc, PCCompositeType *type)
389: {
390:   PC_PATCH *patch = (PC_PATCH *) pc->data;
392:   *type = patch->local_composition_type;
393:   return(0);
394: }

396: /* TODO: Docs */
397: PetscErrorCode PCPatchSetSubMatType(PC pc, MatType sub_mat_type)
398: {
399:   PC_PATCH      *patch = (PC_PATCH *) pc->data;

403:   if (patch->sub_mat_type) {PetscFree(patch->sub_mat_type);}
404:   PetscStrallocpy(sub_mat_type, (char **) &patch->sub_mat_type);
405:   return(0);
406: }

408: /* TODO: Docs */
409: PetscErrorCode PCPatchGetSubMatType(PC pc, MatType *sub_mat_type)
410: {
411:   PC_PATCH *patch = (PC_PATCH *) pc->data;
413:   *sub_mat_type = patch->sub_mat_type;
414:   return(0);
415: }

417: /* TODO: Docs */
418: PetscErrorCode PCPatchSetCellNumbering(PC pc, PetscSection cellNumbering)
419: {
420:   PC_PATCH      *patch = (PC_PATCH *) pc->data;

424:   patch->cellNumbering = cellNumbering;
425:   PetscObjectReference((PetscObject) cellNumbering);
426:   return(0);
427: }

429: /* TODO: Docs */
430: PetscErrorCode PCPatchGetCellNumbering(PC pc, PetscSection *cellNumbering)
431: {
432:   PC_PATCH *patch = (PC_PATCH *) pc->data;
434:   *cellNumbering = patch->cellNumbering;
435:   return(0);
436: }

438: /* TODO: Docs */
439: PetscErrorCode PCPatchSetConstructType(PC pc, PCPatchConstructType ctype, PetscErrorCode (*func)(PC, PetscInt *, IS **, IS *, void *), void *ctx)
440: {
441:   PC_PATCH *patch = (PC_PATCH *) pc->data;

444:   patch->ctype = ctype;
445:   switch (ctype) {
446:   case PC_PATCH_STAR:
447:     patch->user_patches     = PETSC_FALSE;
448:     patch->patchconstructop = PCPatchConstruct_Star;
449:     break;
450:   case PC_PATCH_VANKA:
451:     patch->user_patches     = PETSC_FALSE;
452:     patch->patchconstructop = PCPatchConstruct_Vanka;
453:     break;
454:   case PC_PATCH_PARDECOMP:
455:     patch->user_patches     = PETSC_FALSE;
456:     patch->patchconstructop = PCPatchConstruct_Pardecomp;
457:     break;
458:   case PC_PATCH_USER:
459:   case PC_PATCH_PYTHON:
460:     patch->user_patches     = PETSC_TRUE;
461:     patch->patchconstructop = PCPatchConstruct_User;
462:     if (func) {
463:       patch->userpatchconstructionop = func;
464:       patch->userpatchconstructctx   = ctx;
465:     }
466:     break;
467:   default:
468:     SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_USER, "Unknown patch construction type %D", (PetscInt) patch->ctype);
469:   }
470:   return(0);
471: }

473: /* TODO: Docs */
474: PetscErrorCode PCPatchGetConstructType(PC pc, PCPatchConstructType *ctype, PetscErrorCode (**func)(PC, PetscInt *, IS **, IS *, void *), void **ctx)
475: {
476:   PC_PATCH *patch = (PC_PATCH *) pc->data;

479:   *ctype = patch->ctype;
480:   switch (patch->ctype) {
481:   case PC_PATCH_STAR:
482:   case PC_PATCH_VANKA:
483:   case PC_PATCH_PARDECOMP:
484:     break;
485:   case PC_PATCH_USER:
486:   case PC_PATCH_PYTHON:
487:     *func = patch->userpatchconstructionop;
488:     *ctx  = patch->userpatchconstructctx;
489:     break;
490:   default:
491:     SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_USER, "Unknown patch construction type %D", (PetscInt) patch->ctype);
492:   }
493:   return(0);
494: }

496: /* TODO: Docs */
497: PetscErrorCode PCPatchSetDiscretisationInfo(PC pc, PetscInt nsubspaces, DM *dms, PetscInt *bs, PetscInt *nodesPerCell, const PetscInt **cellNodeMap,
498:                                             const PetscInt *subspaceOffsets, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes)
499: {
500:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
501:   DM             dm, plex;
502:   PetscSF       *sfs;
503:   PetscInt       cStart, cEnd, i, j;

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

518:   patch->nsubspaces       = nsubspaces;
519:   patch->totalDofsPerCell = 0;
520:   for (i = 0; i < nsubspaces; ++i) {
521:     DMGetDefaultSection(dms[i], &patch->dofSection[i]);
522:     PetscObjectReference((PetscObject) patch->dofSection[i]);
523:     DMGetDefaultSF(dms[i], &sfs[i]);
524:     patch->bs[i]              = bs[i];
525:     patch->nodesPerCell[i]    = nodesPerCell[i];
526:     patch->totalDofsPerCell  += nodesPerCell[i]*bs[i];
527:     PetscMalloc1((cEnd-cStart)*nodesPerCell[i], &patch->cellNodeMap[i]);
528:     for (j = 0; j < (cEnd-cStart)*nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j];
529:     patch->subspaceOffsets[i] = subspaceOffsets[i];
530:   }
531:   PCPatchCreateDefaultSF_Private(pc, nsubspaces, sfs, patch->bs);
532:   PetscFree(sfs);

534:   patch->subspaceOffsets[nsubspaces] = subspaceOffsets[nsubspaces];
535:   ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);
536:   ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);
537:   DMDestroy(&dm);
538:   return(0);
539: }

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

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

575: /*@C

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

579:   Logically collective on PC

581:   Input Parameters:
582: + pc   - The PC
583: . func - The callback
584: - ctx  - The user context

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

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

599:   Level: advanced

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

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

611:   patch->usercomputef    = func;
612:   patch->usercomputefctx = ctx;
613:   return(0);
614: }

616: /*@C

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

620:   Logically collective on PC

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

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

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

640:   Level: advanced

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

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

652:   patch->usercomputefintfacet    = func;
653:   patch->usercomputefintfacetctx = ctx;
654:   return(0);
655: }

657: /*@C

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

661:   Logically collective on PC

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

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

671: +  pc               - The PC
672: .  point            - The point
673: .  x                - The input solution (not used in linear problems)
674: .  mat              - The patch matrix
675: .  cellIS           - An array of the cell numbers
676: .  n                - The size of dofsArray
677: .  dofsArray        - The dofmap for the dofs to be solved for
678: .  dofsArrayWithAll - The dofmap for all dofs on the patch
679: -  ctx              - The user context

681:   Level: advanced

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

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

693:   patch->usercomputeop    = func;
694:   patch->usercomputeopctx = ctx;
695:   return(0);
696: }

698: /*@C

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

702:   Logically collective on PC

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

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

712: +  pc               - The PC
713: .  point            - The point
714: .  x                - The input solution (not used in linear problems)
715: .  mat              - The patch matrix
716: .  facetIS          - An array of the facet numbers
717: .  n                - The size of dofsArray
718: .  dofsArray        - The dofmap for the dofs to be solved for
719: .  dofsArrayWithAll - The dofmap for all dofs on the patch
720: -  ctx              - The user context

722:   Level: advanced

724:   Notes:
725:   The matrix entries have been set to zero before the call.

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

734:   patch->usercomputeopintfacet    = func;
735:   patch->usercomputeopintfacetctx = ctx;
736:   return(0);
737: }

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

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

767:     PetscHashIterGetKey(ht, hi, point);
768:     PetscHashIterNext(ht, hi);

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

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

811:   if (combined) {
812:     if (f < 0) {
813:       if (dof) {PetscSectionGetDof(dofSection[0], p, dof);}
814:       if (off) {PetscSectionGetOffset(dofSection[0], p, off);}
815:     } else {
816:       if (dof) {PetscSectionGetFieldDof(dofSection[0], p, f, dof);}
817:       if (off) {PetscSectionGetFieldOffset(dofSection[0], p, f, off);}
818:     }
819:   } else {
820:     if (f < 0) {
821:       PC_PATCH *patch = (PC_PATCH *) pc->data;
822:       PetscInt  fdof, g;

824:       if (dof) {
825:         *dof = 0;
826:         for (g = 0; g < patch->nsubspaces; ++g) {
827:           PetscSectionGetDof(dofSection[g], p, &fdof);
828:           *dof += fdof;
829:         }
830:       }
831:       if (off) {
832:         *off = 0;
833:         for (g = 0; g < patch->nsubspaces; ++g) {
834:           PetscSectionGetOffset(dofSection[g], p, &fdof);
835:           *off += fdof;
836:         }
837:       }
838:     } else {
839:       if (dof) {PetscSectionGetDof(dofSection[f], p, dof);}
840:       if (off) {PetscSectionGetOffset(dofSection[f], p, off);}
841:     }
842:   }
843:   return(0);
844: }

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

860:   PetscHSetIClear(dofs);
861:   for (k = 0; k < patch->nsubspaces; ++k) {
862:     PetscInt subspaceOffset = patch->subspaceOffsets[k];
863:     PetscInt bs             = patch->bs[k];
864:     PetscInt j, l;

866:     if (subspaces_to_exclude != NULL) {
867:       PetscBool should_exclude_k = PETSC_FALSE;
868:       PetscHSetIHas(*subspaces_to_exclude, k, &should_exclude_k);
869:       if (should_exclude_k) {
870:         /* only get this subspace dofs at the base entity, not any others */
871:         PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, base, &ldof, &loff);
872:         if (0 == ldof) continue;
873:         for (j = loff; j < ldof + loff; ++j) {
874:           for (l = 0; l < bs; ++l) {
875:             PetscInt dof = bs*j + l + subspaceOffset;
876:             PetscHSetIAdd(dofs, dof);
877:           }
878:         }
879:         continue; /* skip the other dofs of this subspace */
880:       }
881:     }

883:     PetscHashIterBegin(pts, hi);
884:     while (!PetscHashIterAtEnd(pts, hi)) {
885:       PetscHashIterGetKey(pts, hi, p);
886:       PetscHashIterNext(pts, hi);
887:       PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, p, &ldof, &loff);
888:       if (0 == ldof) continue;
889:       for (j = loff; j < ldof + loff; ++j) {
890:         for (l = 0; l < bs; ++l) {
891:           PetscInt dof = bs*j + l + subspaceOffset;
892:           PetscHSetIAdd(dofs, dof);
893:         }
894:       }
895:     }
896:   }
897:   return(0);
898: }

900: /* Given two hash tables A and B, compute the keys in B that are not in A, and put them in C */
901: static PetscErrorCode PCPatchComputeSetDifference_Private(PetscHSetI A, PetscHSetI B, PetscHSetI C)
902: {
903:   PetscHashIter  hi;
904:   PetscInt       key;
905:   PetscBool      flg;

909:   PetscHSetIClear(C);
910:   PetscHashIterBegin(B, hi);
911:   while (!PetscHashIterAtEnd(B, hi)) {
912:     PetscHashIterGetKey(B, hi, key);
913:     PetscHashIterNext(B, hi);
914:     PetscHSetIHas(A, key, &flg);
915:     if (!flg) {PetscHSetIAdd(C, key);}
916:   }
917:   return(0);
918: }

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

947:   /* Used to keep track of the cells in the patch. */
948:   PetscHSetICreate(&ht);
949:   PetscHSetICreate(&cht);

951:   PCGetDM(pc, &dm);
952:   if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONGSTATE, "DM not yet set on patch PC\n");
953:   DMConvert(dm, DMPLEX, &plex);
954:   dm = plex;
955:   DMPlexGetChart(dm, &pStart, &pEnd);
956:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);

958:   if (patch->user_patches) {
959:     patch->userpatchconstructionop(pc, &patch->npatch, &patch->userIS, &patch->iterationSet, patch->userpatchconstructctx);
960:     vStart = 0; vEnd = patch->npatch;
961:   } else if (patch->ctype == PC_PATCH_PARDECOMP) {
962:     vStart = 0; vEnd = 1;
963:   } else if (patch->codim < 0) {
964:     if (patch->dim < 0) {DMPlexGetDepthStratum(dm,  0,            &vStart, &vEnd);}
965:     else                {DMPlexGetDepthStratum(dm,  patch->dim,   &vStart, &vEnd);}
966:   } else                {DMPlexGetHeightStratum(dm, patch->codim, &vStart, &vEnd);}
967:   patch->npatch = vEnd - vStart;

969:   /* These labels mark the owned points.  We only create patches around points that this process owns. */
970:   DMHasLabel(dm, "pyop2_ghost", &isFiredrake);
971:   if (isFiredrake) {
972:     DMGetLabel(dm, "pyop2_ghost", &ghost);
973:     DMLabelCreateIndex(ghost, pStart, pEnd);
974:   } else {
975:     PetscSF sf;

977:     DMGetPointSF(dm, &sf);
978:     PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
979:     nleaves = PetscMax(nleaves, 0);
980:   }

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

1005:     if (!patch->user_patches && patch->ctype != PC_PATCH_PARDECOMP) {
1006:       if (ghost) {DMLabelHasPoint(ghost, v, &flg);}
1007:       else       {PetscFindInt(v, nleaves, leaves, &loc); flg = loc >=0 ? PETSC_TRUE : PETSC_FALSE;}
1008:       /* Not an owned entity, don't make a cell patch. */
1009:       if (flg) continue;
1010:     }

1012:     patch->patchconstructop((void *) patch, dm, v, ht);
1013:     PCPatchCompleteCellPatch(pc, ht, cht);
1014:     PetscHSetIGetSize(cht, &chtSize);
1015:     /* empty patch, continue */
1016:     if (chtSize == 0) continue;

1018:     /* safe because size(cht) > 0 from above */
1019:     PetscHashIterBegin(cht, hi);
1020:     while (!PetscHashIterAtEnd(cht, hi)) {
1021:       PetscInt point, pdof;

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

1057:   PetscSectionSetUp(cellCounts);
1058:   PetscSectionGetStorageSize(cellCounts, &numCells);
1059:   PetscMalloc1(numCells, &cellsArray);
1060:   PetscSectionSetUp(pointCounts);
1061:   PetscSectionGetStorageSize(pointCounts, &numPoints);
1062:   PetscMalloc1(numPoints, &pointsArray);

1064:   PetscSectionSetUp(intFacetCounts);
1065:   PetscSectionSetUp(extFacetCounts);
1066:   PetscSectionGetStorageSize(intFacetCounts, &numIntFacets);
1067:   PetscSectionGetStorageSize(extFacetCounts, &numExtFacets);
1068:   PetscMalloc1(numIntFacets, &intFacetsArray);
1069:   PetscMalloc1(numIntFacets*2, &intFacetsToPatchCell);
1070:   PetscMalloc1(numExtFacets, &extFacetsArray);


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

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

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

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

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

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


1246:   PCGetDM(pc, &dm);
1247:   DMConvert(dm, DMPLEX, &plex);
1248:   dm = plex;
1249:   /* dofcounts section is cellcounts section * dofPerCell */
1250:   PetscSectionGetStorageSize(cellCounts, &numCells);
1251:   PetscSectionGetStorageSize(patch->pointCounts, &numPoints);
1252:   numDofs = numCells * totalDofsPerCell;
1253:   PetscMalloc1(numDofs, &dofsArray);
1254:   PetscMalloc1(numPoints*Nf, &offsArray);
1255:   PetscMalloc1(numDofs, &asmArray);
1256:   PetscMalloc1(numCells, &newCellsArray);
1257:   PetscSectionGetChart(cellCounts, &vStart, &vEnd);
1258:   PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCounts);
1259:   gtolCounts = patch->gtolCounts;
1260:   PetscSectionSetChart(gtolCounts, vStart, vEnd);
1261:   PetscObjectSetName((PetscObject) patch->gtolCounts, "Patch Global Index Section");

1263:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1264:     PetscMalloc1(numPoints*Nf, &offsArrayWithArtificial);
1265:     PetscMalloc1(numDofs, &asmArrayWithArtificial);
1266:     PetscMalloc1(numDofs, &dofsArrayWithArtificial);
1267:     PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithArtificial);
1268:     gtolCountsWithArtificial = patch->gtolCountsWithArtificial;
1269:     PetscSectionSetChart(gtolCountsWithArtificial, vStart, vEnd);
1270:     PetscObjectSetName((PetscObject) patch->gtolCountsWithArtificial, "Patch Global Index Section Including Artificial BCs");
1271:   }

1273:   isNonlinear = patch->isNonlinear;
1274:   if (isNonlinear) {
1275:     PetscMalloc1(numPoints*Nf, &offsArrayWithAll);
1276:     PetscMalloc1(numDofs, &asmArrayWithAll);
1277:     PetscMalloc1(numDofs, &dofsArrayWithAll);
1278:     PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithAll);
1279:     gtolCountsWithAll = patch->gtolCountsWithAll;
1280:     PetscSectionSetChart(gtolCountsWithAll, vStart, vEnd);
1281:     PetscObjectSetName((PetscObject) patch->gtolCountsWithAll, "Patch Global Index Section Including All BCs");
1282:   }

1284:   /* Outside the patch loop, get the dofs that are globally-enforced Dirichlet
1285:    conditions */
1286:   PetscHSetICreate(&globalBcs);
1287:   ISGetIndices(patch->ghostBcNodes, &bcNodes);
1288:   ISGetSize(patch->ghostBcNodes, &numBcs);
1289:   for (i = 0; i < numBcs; ++i) {
1290:     PetscHSetIAdd(globalBcs, bcNodes[i]); /* these are already in concatenated numbering */
1291:   }
1292:   ISRestoreIndices(patch->ghostBcNodes, &bcNodes);
1293:   ISDestroy(&patch->ghostBcNodes);  /* memory optimisation */

1295:   /* Hash tables for artificial BC construction */
1296:   PetscHSetICreate(&ownedpts);
1297:   PetscHSetICreate(&seenpts);
1298:   PetscHSetICreate(&owneddofs);
1299:   PetscHSetICreate(&seendofs);
1300:   PetscHSetICreate(&artificialbcs);

1302:   ISGetIndices(cells, &cellsArray);
1303:   ISGetIndices(points, &pointsArray);
1304:   PetscHMapICreate(&ht);
1305:   PetscHMapICreate(&htWithArtificial);
1306:   PetscHMapICreate(&htWithAll);
1307:   for (v = vStart; v < vEnd; ++v) {
1308:     PetscInt localIndex = 0;
1309:     PetscInt localIndexWithArtificial = 0;
1310:     PetscInt localIndexWithAll = 0;
1311:     PetscInt dof, off, i, j, k, l;

1313:     PetscHMapIClear(ht);
1314:     PetscHMapIClear(htWithArtificial);
1315:     PetscHMapIClear(htWithAll);
1316:     PetscSectionGetDof(cellCounts, v, &dof);
1317:     PetscSectionGetOffset(cellCounts, v, &off);
1318:     if (dof <= 0) continue;

1320:     /* Calculate the global numbers of the artificial BC dofs here first */
1321:     patch->patchconstructop((void*)patch, dm, v, ownedpts);
1322:     PCPatchCompleteCellPatch(pc, ownedpts, seenpts);
1323:     PCPatchGetPointDofs(pc, ownedpts, owneddofs, v, &patch->subspaces_to_exclude);
1324:     PCPatchGetPointDofs(pc, seenpts, seendofs, v, NULL);
1325:     PCPatchComputeSetDifference_Private(owneddofs, seendofs, artificialbcs);
1326:     if (patch->viewPatches) {
1327:       PetscHSetI globalbcdofs;
1328:       PetscHashIter hi;
1329:       MPI_Comm comm = PetscObjectComm((PetscObject)pc);

1331:       PetscHSetICreate(&globalbcdofs);
1332:       PetscSynchronizedPrintf(comm, "Patch %d: owned dofs:\n", v);
1333:       PetscHashIterBegin(owneddofs, hi);
1334:       while (!PetscHashIterAtEnd(owneddofs, hi)) {
1335:         PetscInt globalDof;

1337:         PetscHashIterGetKey(owneddofs, hi, globalDof);
1338:         PetscHashIterNext(owneddofs, hi);
1339:         PetscSynchronizedPrintf(comm, "%d ", globalDof);
1340:       }
1341:       PetscSynchronizedPrintf(comm, "\n");
1342:       PetscSynchronizedPrintf(comm, "Patch %d: seen dofs:\n", v);
1343:       PetscHashIterBegin(seendofs, hi);
1344:       while (!PetscHashIterAtEnd(seendofs, hi)) {
1345:         PetscInt globalDof;
1346:         PetscBool flg;

1348:         PetscHashIterGetKey(seendofs, hi, globalDof);
1349:         PetscHashIterNext(seendofs, hi);
1350:         PetscSynchronizedPrintf(comm, "%d ", globalDof);

1352:         PetscHSetIHas(globalBcs, globalDof, &flg);
1353:         if (flg) {PetscHSetIAdd(globalbcdofs, globalDof);}
1354:       }
1355:       PetscSynchronizedPrintf(comm, "\n");
1356:       PetscSynchronizedPrintf(comm, "Patch %d: global BCs:\n", v);
1357:       PetscHSetIGetSize(globalbcdofs, &numBcs);
1358:       if (numBcs > 0) {
1359:         PetscHashIterBegin(globalbcdofs, hi);
1360:         while (!PetscHashIterAtEnd(globalbcdofs, hi)) {
1361:           PetscInt globalDof;
1362:           PetscHashIterGetKey(globalbcdofs, hi, globalDof);
1363:           PetscHashIterNext(globalbcdofs, hi);
1364:           PetscSynchronizedPrintf(comm, "%d ", globalDof);
1365:         }
1366:       }
1367:       PetscSynchronizedPrintf(comm, "\n");
1368:       PetscSynchronizedPrintf(comm, "Patch %d: artificial BCs:\n", v);
1369:       PetscHSetIGetSize(artificialbcs, &numBcs);
1370:       if (numBcs > 0) {
1371:         PetscHashIterBegin(artificialbcs, hi);
1372:         while (!PetscHashIterAtEnd(artificialbcs, hi)) {
1373:           PetscInt globalDof;
1374:           PetscHashIterGetKey(artificialbcs, hi, globalDof);
1375:           PetscHashIterNext(artificialbcs, hi);
1376:           PetscSynchronizedPrintf(comm, "%d ", globalDof);
1377:         }
1378:       }
1379:       PetscSynchronizedPrintf(comm, "\n\n");
1380:       PetscHSetIDestroy(&globalbcdofs);
1381:     }
1382:    for (k = 0; k < patch->nsubspaces; ++k) {
1383:       const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
1384:       PetscInt        nodesPerCell   = patch->nodesPerCell[k];
1385:       PetscInt        subspaceOffset = patch->subspaceOffsets[k];
1386:       PetscInt        bs             = patch->bs[k];

1388:       for (i = off; i < off + dof; ++i) {
1389:         /* Walk over the cells in this patch. */
1390:         const PetscInt c    = cellsArray[i];
1391:         PetscInt       cell = c;

1393:         /* TODO Change this to an IS */
1394:         if (cellNumbering) {
1395:           PetscSectionGetDof(cellNumbering, c, &cell);
1396:           if (cell <= 0) SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_OUTOFRANGE, "Cell %D doesn't appear in cell numbering map", c);
1397:           PetscSectionGetOffset(cellNumbering, c, &cell);
1398:         }
1399:         newCellsArray[i] = cell;
1400:         for (j = 0; j < nodesPerCell; ++j) {
1401:           /* For each global dof, map it into contiguous local storage. */
1402:           const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + subspaceOffset;
1403:           /* finally, loop over block size */
1404:           for (l = 0; l < bs; ++l) {
1405:             PetscInt  localDof;
1406:             PetscBool isGlobalBcDof, isArtificialBcDof;

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

1412:             /* if it's either, don't ever give it a local dof number */
1413:             if (isGlobalBcDof || isArtificialBcDof) {
1414:               dofsArray[globalIndex] = -1; /* don't use this in assembly in this patch */
1415:             } else {
1416:               PetscHMapIGet(ht, globalDof + l, &localDof);
1417:               if (localDof == -1) {
1418:                 localDof = localIndex++;
1419:                 PetscHMapISet(ht, globalDof + l, localDof);
1420:               }
1421:               if ( globalIndex >= numDofs ) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs);
1422:               /* And store. */
1423:               dofsArray[globalIndex] = localDof;
1424:             }

1426:             if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1427:               if (isGlobalBcDof) {
1428:                 dofsArrayWithArtificial[globalIndex] = -1; /* don't use this in assembly in this patch */
1429:               } else {
1430:                 PetscHMapIGet(htWithArtificial, globalDof + l, &localDof);
1431:                 if (localDof == -1) {
1432:                   localDof = localIndexWithArtificial++;
1433:                   PetscHMapISet(htWithArtificial, globalDof + l, localDof);
1434:                 }
1435:                 if ( globalIndex >= numDofs ) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs);
1436:                 /* And store.*/
1437:                 dofsArrayWithArtificial[globalIndex] = localDof;
1438:               }
1439:             }

1441:             if(isNonlinear) {
1442:               /* Build the dofmap for the function space with _all_ dofs,
1443:                  including those in any kind of boundary condition */
1444:               PetscHMapIGet(htWithAll, globalDof + l, &localDof);
1445:               if (localDof == -1) {
1446:                 localDof = localIndexWithAll++;
1447:                 PetscHMapISet(htWithAll, globalDof + l, localDof);
1448:               }
1449:               if ( globalIndex >= numDofs ) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs);
1450:               /* And store.*/
1451:               dofsArrayWithAll[globalIndex] = localDof;
1452:             }
1453:             globalIndex++;
1454:           }
1455:         }
1456:       }
1457:     }
1458:      /*How many local dofs in this patch? */
1459:    if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1460:      PetscHMapIGetSize(htWithArtificial, &dof);
1461:      PetscSectionSetDof(gtolCountsWithArtificial, v, dof);
1462:    }
1463:    if (isNonlinear) {
1464:      PetscHMapIGetSize(htWithAll, &dof);
1465:      PetscSectionSetDof(gtolCountsWithAll, v, dof);
1466:    }
1467:     PetscHMapIGetSize(ht, &dof);
1468:     PetscSectionSetDof(gtolCounts, v, dof);
1469:   }

1471:   DMDestroy(&dm);
1472:   if (globalIndex != numDofs) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Expected number of dofs (%d) doesn't match found number (%d)", numDofs, globalIndex);
1473:   PetscSectionSetUp(gtolCounts);
1474:   PetscSectionGetStorageSize(gtolCounts, &numGlobalDofs);
1475:   PetscMalloc1(numGlobalDofs, &globalDofsArray);

1477:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1478:     PetscSectionSetUp(gtolCountsWithArtificial);
1479:     PetscSectionGetStorageSize(gtolCountsWithArtificial, &numGlobalDofsWithArtificial);
1480:     PetscMalloc1(numGlobalDofsWithArtificial, &globalDofsArrayWithArtificial);
1481:   }
1482:   if (isNonlinear) {
1483:     PetscSectionSetUp(gtolCountsWithAll);
1484:     PetscSectionGetStorageSize(gtolCountsWithAll, &numGlobalDofsWithAll);
1485:     PetscMalloc1(numGlobalDofsWithAll, &globalDofsArrayWithAll);
1486:   }
1487:   /* Now populate the global to local map.  This could be merged into the above loop if we were willing to deal with reallocs. */
1488:   for (v = vStart; v < vEnd; ++v) {
1489:     PetscHashIter hi;
1490:     PetscInt      dof, off, Np, ooff, i, j, k, l;

1492:     PetscHMapIClear(ht);
1493:     PetscHMapIClear(htWithArtificial);
1494:     PetscHMapIClear(htWithAll);
1495:     PetscSectionGetDof(cellCounts, v, &dof);
1496:     PetscSectionGetOffset(cellCounts, v, &off);
1497:     PetscSectionGetDof(pointCounts, v, &Np);
1498:     PetscSectionGetOffset(pointCounts, v, &ooff);
1499:     if (dof <= 0) continue;

1501:     for (k = 0; k < patch->nsubspaces; ++k) {
1502:       const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
1503:       PetscInt        nodesPerCell   = patch->nodesPerCell[k];
1504:       PetscInt        subspaceOffset = patch->subspaceOffsets[k];
1505:       PetscInt        bs             = patch->bs[k];
1506:       PetscInt        goff;

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

1513:         if (cellNumbering) {PetscSectionGetOffset(cellNumbering, c, &cell);}
1514:         for (j = 0; j < nodesPerCell; ++j) {
1515:           for (l = 0; l < bs; ++l) {
1516:             const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset;
1517:             const PetscInt localDof  = dofsArray[key];
1518:             if (localDof >= 0) {PetscHMapISet(ht, globalDof, localDof);}
1519:             if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1520:               const PetscInt localDofWithArtificial = dofsArrayWithArtificial[key];
1521:               if (localDofWithArtificial >= 0) {
1522:                 PetscHMapISet(htWithArtificial, globalDof, localDofWithArtificial);
1523:               }
1524:             }
1525:             if (isNonlinear) {
1526:               const PetscInt localDofWithAll = dofsArrayWithAll[key];
1527:               if (localDofWithAll >= 0) {
1528:                 PetscHMapISet(htWithAll, globalDof, localDofWithAll);
1529:               }
1530:             }
1531:             key++;
1532:           }
1533:         }
1534:       }

1536:       /* Shove it in the output data structure. */
1537:       PetscSectionGetOffset(gtolCounts, v, &goff);
1538:       PetscHashIterBegin(ht, hi);
1539:       while (!PetscHashIterAtEnd(ht, hi)) {
1540:         PetscInt globalDof, localDof;

1542:         PetscHashIterGetKey(ht, hi, globalDof);
1543:         PetscHashIterGetVal(ht, hi, localDof);
1544:         if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof;
1545:         PetscHashIterNext(ht, hi);
1546:       }

1548:       if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1549:         PetscSectionGetOffset(gtolCountsWithArtificial, v, &goff);
1550:         PetscHashIterBegin(htWithArtificial, hi);
1551:         while (!PetscHashIterAtEnd(htWithArtificial, hi)) {
1552:           PetscInt globalDof, localDof;
1553:           PetscHashIterGetKey(htWithArtificial, hi, globalDof);
1554:           PetscHashIterGetVal(htWithArtificial, hi, localDof);
1555:           if (globalDof >= 0) globalDofsArrayWithArtificial[goff + localDof] = globalDof;
1556:           PetscHashIterNext(htWithArtificial, hi);
1557:         }
1558:       }
1559:       if (isNonlinear) {
1560:         PetscSectionGetOffset(gtolCountsWithAll, v, &goff);
1561:         PetscHashIterBegin(htWithAll, hi);
1562:         while (!PetscHashIterAtEnd(htWithAll, hi)) {
1563:           PetscInt globalDof, localDof;
1564:           PetscHashIterGetKey(htWithAll, hi, globalDof);
1565:           PetscHashIterGetVal(htWithAll, hi, localDof);
1566:           if (globalDof >= 0) globalDofsArrayWithAll[goff + localDof] = globalDof;
1567:           PetscHashIterNext(htWithAll, hi);
1568:         }
1569:       }

1571:       for (p = 0; p < Np; ++p) {
1572:         const PetscInt point = pointsArray[ooff + p];
1573:         PetscInt       globalDof, localDof;

1575:         PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof);
1576:         PetscHMapIGet(ht, globalDof, &localDof);
1577:         offsArray[(ooff + p)*Nf + k] = localDof;
1578:         if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1579:           PetscHMapIGet(htWithArtificial, globalDof, &localDof);
1580:           offsArrayWithArtificial[(ooff + p)*Nf + k] = localDof;
1581:         }
1582:         if (isNonlinear) {
1583:           PetscHMapIGet(htWithAll, globalDof, &localDof);
1584:           offsArrayWithAll[(ooff + p)*Nf + k] = localDof;
1585:         }
1586:       }
1587:     }

1589:     PetscHSetIDestroy(&globalBcs);
1590:     PetscHSetIDestroy(&ownedpts);
1591:     PetscHSetIDestroy(&seenpts);
1592:     PetscHSetIDestroy(&owneddofs);
1593:     PetscHSetIDestroy(&seendofs);
1594:     PetscHSetIDestroy(&artificialbcs);

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

1606:         if (cellNumbering) {PetscSectionGetOffset(cellNumbering, c, &cell);}
1607:         for (k = 0; k < patch->nsubspaces; ++k) {
1608:           const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
1609:           PetscInt        nodesPerCell   = patch->nodesPerCell[k];
1610:           PetscInt        subspaceOffset = patch->subspaceOffsets[k];
1611:           PetscInt        bs             = patch->bs[k];

1613:           for (j = 0; j < nodesPerCell; ++j) {
1614:             for (l = 0; l < bs; ++l) {
1615:               const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset;
1616:               PetscInt       localDof;

1618:               PetscHMapIGet(ht, globalDof, &localDof);
1619:               /* If it's not in the hash table, i.e. is a BC dof,
1620:                then the PetscHSetIMap above gives -1, which matches
1621:                exactly the convention for PETSc's matrix assembly to
1622:                ignore the dof. So we don't need to do anything here */
1623:               asmArray[asmKey] = localDof;
1624:               if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1625:                 PetscHMapIGet(htWithArtificial, globalDof, &localDof);
1626:                 asmArrayWithArtificial[asmKey] = localDof;
1627:               }
1628:               if (isNonlinear) {
1629:                 PetscHMapIGet(htWithAll, globalDof, &localDof);
1630:                 asmArrayWithAll[asmKey] = localDof;
1631:               }
1632:               asmKey++;
1633:             }
1634:           }
1635:         }
1636:       }
1637:     }
1638:   }
1639:   if (1 == patch->nsubspaces) {
1640:     PetscArraycpy(asmArray, dofsArray, numDofs);
1641:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1642:       PetscArraycpy(asmArrayWithArtificial, dofsArrayWithArtificial, numDofs);
1643:     }
1644:     if (isNonlinear) {
1645:       PetscArraycpy(asmArrayWithAll, dofsArrayWithAll, numDofs);
1646:     }
1647:   }

1649:   PetscHMapIDestroy(&ht);
1650:   PetscHMapIDestroy(&htWithArtificial);
1651:   PetscHMapIDestroy(&htWithAll);
1652:   ISRestoreIndices(cells, &cellsArray);
1653:   ISRestoreIndices(points, &pointsArray);
1654:   PetscFree(dofsArray);
1655:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1656:     PetscFree(dofsArrayWithArtificial);
1657:   }
1658:   if (isNonlinear) {
1659:     PetscFree(dofsArrayWithAll);
1660:   }
1661:   /* Create placeholder section for map from points to patch dofs */
1662:   PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection);
1663:   PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces);
1664:   if (patch->combined) {
1665:     PetscInt numFields;
1666:     PetscSectionGetNumFields(patch->dofSection[0], &numFields);
1667:     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);
1668:     PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd);
1669:     PetscSectionSetChart(patch->patchSection, pStart, pEnd);
1670:     for (p = pStart; p < pEnd; ++p) {
1671:       PetscInt dof, fdof, f;

1673:       PetscSectionGetDof(patch->dofSection[0], p, &dof);
1674:       PetscSectionSetDof(patch->patchSection, p, dof);
1675:       for (f = 0; f < patch->nsubspaces; ++f) {
1676:         PetscSectionGetFieldDof(patch->dofSection[0], p, f, &fdof);
1677:         PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);
1678:       }
1679:     }
1680:   } else {
1681:     PetscInt pStartf, pEndf, f;
1682:     pStart = PETSC_MAX_INT;
1683:     pEnd = PETSC_MIN_INT;
1684:     for (f = 0; f < patch->nsubspaces; ++f) {
1685:       PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);
1686:       pStart = PetscMin(pStart, pStartf);
1687:       pEnd = PetscMax(pEnd, pEndf);
1688:     }
1689:     PetscSectionSetChart(patch->patchSection, pStart, pEnd);
1690:     for (f = 0; f < patch->nsubspaces; ++f) {
1691:       PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);
1692:       for (p = pStartf; p < pEndf; ++p) {
1693:         PetscInt fdof;
1694:         PetscSectionGetDof(patch->dofSection[f], p, &fdof);
1695:         PetscSectionAddDof(patch->patchSection, p, fdof);
1696:         PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);
1697:       }
1698:     }
1699:   }
1700:   PetscSectionSetUp(patch->patchSection);
1701:   PetscSectionSetUseFieldOffsets(patch->patchSection, PETSC_TRUE);
1702:   /* Replace cell indices with firedrake-numbered ones. */
1703:   ISGeneralSetIndices(cells, numCells, (const PetscInt *) newCellsArray, PETSC_OWN_POINTER);
1704:   ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofs, globalDofsArray, PETSC_OWN_POINTER, &patch->gtol);
1705:   PetscObjectSetName((PetscObject) patch->gtol, "Global Indices");
1706:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_g2l_view", patch->classname);
1707:   PetscSectionViewFromOptions(patch->gtolCounts, (PetscObject) pc, option);
1708:   ISViewFromOptions(patch->gtol, (PetscObject) pc, option);
1709:   ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArray, PETSC_OWN_POINTER, &patch->dofs);
1710:   ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArray, PETSC_OWN_POINTER, &patch->offs);
1711:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1712:     ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithArtificial, globalDofsArrayWithArtificial, PETSC_OWN_POINTER, &patch->gtolWithArtificial);
1713:     ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithArtificial, PETSC_OWN_POINTER, &patch->dofsWithArtificial);
1714:     ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArrayWithArtificial, PETSC_OWN_POINTER, &patch->offsWithArtificial);
1715:   }
1716:   if (isNonlinear) {
1717:     ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithAll, globalDofsArrayWithAll, PETSC_OWN_POINTER, &patch->gtolWithAll);
1718:     ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithAll, PETSC_OWN_POINTER, &patch->dofsWithAll);
1719:     ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArrayWithAll, PETSC_OWN_POINTER, &patch->offsWithAll);
1720:   }
1721:   return(0);
1722: }

1724: static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat, PetscBool withArtificial)
1725: {
1726:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
1727:   Vec            x, y;
1728:   PetscBool      flg;
1729:   PetscInt       csize, rsize;
1730:   const char    *prefix = NULL;

1734:   if(withArtificial) {
1735:     /* would be nice if we could create a rectangular matrix of size numDofsWithArtificial x numDofs here */
1736:     x = patch->patchRHSWithArtificial[point];
1737:     y = patch->patchRHSWithArtificial[point];
1738:   } else {
1739:     x = patch->patchRHS[point];
1740:     y = patch->patchUpdate[point];
1741:   }

1743:   VecGetSize(x, &csize);
1744:   VecGetSize(y, &rsize);
1745:   MatCreate(PETSC_COMM_SELF, mat);
1746:   PCGetOptionsPrefix(pc, &prefix);
1747:   MatSetOptionsPrefix(*mat, prefix);
1748:   MatAppendOptionsPrefix(*mat, "pc_patch_sub_");
1749:   if (patch->sub_mat_type)       {MatSetType(*mat, patch->sub_mat_type);}
1750:   else if (!patch->sub_mat_type) {MatSetType(*mat, MATDENSE);}
1751:   MatSetSizes(*mat, rsize, csize, rsize, csize);
1752:   PetscObjectTypeCompare((PetscObject) *mat, MATDENSE, &flg);
1753:   if (!flg) {PetscObjectTypeCompare((PetscObject)*mat, MATSEQDENSE, &flg);}
1754:   /* Sparse patch matrices */
1755:   if (!flg) {
1756:     PetscBT         bt;
1757:     PetscInt       *dnnz      = NULL;
1758:     const PetscInt *dofsArray = NULL;
1759:     PetscInt        pStart, pEnd, ncell, offset, c, i, j;

1761:     if(withArtificial) {
1762:       ISGetIndices(patch->dofsWithArtificial, &dofsArray);
1763:     } else {
1764:       ISGetIndices(patch->dofs, &dofsArray);
1765:     }
1766:     PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);
1767:     point += pStart;
1768:     if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);
1769:     PetscSectionGetDof(patch->cellCounts, point, &ncell);
1770:     PetscSectionGetOffset(patch->cellCounts, point, &offset);
1771:     PetscLogEventBegin(PC_Patch_Prealloc, pc, 0, 0, 0);
1772:     /* A PetscBT uses N^2 bits to store the sparsity pattern on a
1773:      * patch. This is probably OK if the patches are not too big,
1774:      * but uses too much memory. We therefore switch based on rsize. */
1775:     if (rsize < 3000) { /* FIXME: I picked this switch value out of my hat */
1776:       PetscScalar *zeroes;
1777:       PetscInt rows;

1779:       PetscCalloc1(rsize, &dnnz);
1780:       PetscBTCreate(rsize*rsize, &bt);
1781:       for (c = 0; c < ncell; ++c) {
1782:         const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell;
1783:         for (i = 0; i < patch->totalDofsPerCell; ++i) {
1784:           const PetscInt row = idx[i];
1785:           if (row < 0) continue;
1786:           for (j = 0; j < patch->totalDofsPerCell; ++j) {
1787:             const PetscInt col = idx[j];
1788:             const PetscInt key = row*rsize + col;
1789:             if (col < 0) continue;
1790:             if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1791:           }
1792:         }
1793:       }

1795:       if (patch->usercomputeopintfacet) {
1796:         const PetscInt *intFacetsArray = NULL;
1797:         PetscInt        i, numIntFacets, intFacetOffset;
1798:         const PetscInt *facetCells = NULL;

1800:         PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1801:         PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1802:         ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1803:         ISGetIndices(patch->intFacets, &intFacetsArray);
1804:         for (i = 0; i < numIntFacets; i++) {
1805:           const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0];
1806:           const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1];
1807:           PetscInt       celli, cellj;

1809:           for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1810:             const PetscInt row = dofsArray[(offset + cell0)*patch->totalDofsPerCell + celli];
1811:             if (row < 0) continue;
1812:             for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1813:               const PetscInt col = dofsArray[(offset + cell1)*patch->totalDofsPerCell + cellj];
1814:               const PetscInt key = row*rsize + col;
1815:               if (col < 0) continue;
1816:               if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1817:             }
1818:           }

1820:           for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1821:             const PetscInt row = dofsArray[(offset + cell1)*patch->totalDofsPerCell + celli];
1822:             if (row < 0) continue;
1823:             for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1824:               const PetscInt col = dofsArray[(offset + cell0)*patch->totalDofsPerCell + cellj];
1825:               const PetscInt key = row*rsize + col;
1826:               if (col < 0) continue;
1827:               if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1828:             }
1829:           }
1830:         }
1831:       }
1832:       PetscBTDestroy(&bt);
1833:       MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL);
1834:       PetscFree(dnnz);

1836:       PetscCalloc1(patch->totalDofsPerCell*patch->totalDofsPerCell, &zeroes);
1837:       for (c = 0; c < ncell; ++c) {
1838:         const PetscInt *idx = &dofsArray[(offset + c)*patch->totalDofsPerCell];
1839:         MatSetValues(*mat, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, zeroes, INSERT_VALUES);
1840:       }
1841:       MatGetLocalSize(*mat, &rows, NULL);
1842:       for (i = 0; i < rows; ++i) {
1843:         MatSetValues(*mat, 1, &i, 1, &i, zeroes, INSERT_VALUES);
1844:       }

1846:       if (patch->usercomputeopintfacet) {
1847:         const PetscInt *intFacetsArray = NULL;
1848:         PetscInt i, numIntFacets, intFacetOffset;
1849:         const PetscInt *facetCells = NULL;

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

1865:       MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
1866:       MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);

1868:       PetscFree(zeroes);

1870:     } else { /* rsize too big, use MATPREALLOCATOR */
1871:       Mat preallocator;
1872:       PetscScalar* vals;

1874:       PetscCalloc1(patch->totalDofsPerCell*patch->totalDofsPerCell, &vals);
1875:       MatCreate(PETSC_COMM_SELF, &preallocator);
1876:       MatSetType(preallocator, MATPREALLOCATOR);
1877:       MatSetSizes(preallocator, rsize, rsize, rsize, rsize);
1878:       MatSetUp(preallocator);

1880:       for (c = 0; c < ncell; ++c) {
1881:         const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell;
1882:         MatSetValues(preallocator, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, vals, INSERT_VALUES);
1883:       }

1885:       if (patch->usercomputeopintfacet) {
1886:         const PetscInt *intFacetsArray = NULL;
1887:         PetscInt        i, numIntFacets, intFacetOffset;
1888:         const PetscInt *facetCells = NULL;

1890:         PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1891:         PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1892:         ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1893:         ISGetIndices(patch->intFacets, &intFacetsArray);
1894:         for (i = 0; i < numIntFacets; i++) {
1895:           const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0];
1896:           const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1];
1897:           const PetscInt *cell0idx = &dofsArray[(offset + cell0)*patch->totalDofsPerCell];
1898:           const PetscInt *cell1idx = &dofsArray[(offset + cell1)*patch->totalDofsPerCell];
1899:           MatSetValues(preallocator, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, vals, INSERT_VALUES);
1900:           MatSetValues(preallocator, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, vals, INSERT_VALUES);
1901:         }
1902:       }

1904:       PetscFree(vals);
1905:       MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1906:       MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1907:       MatPreallocatorPreallocate(preallocator, PETSC_TRUE, *mat);
1908:       MatDestroy(&preallocator);
1909:     }
1910:     PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0);
1911:     if(withArtificial) {
1912:       ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);
1913:     } else {
1914:       ISRestoreIndices(patch->dofs, &dofsArray);
1915:     }
1916:   }
1917:   MatSetUp(*mat);
1918:   return(0);
1919: }

1921: 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)
1922: {
1923:   PC_PATCH       *patch = (PC_PATCH *) pc->data;
1924:   DM              dm, plex;
1925:   PetscSection    s;
1926:   const PetscInt *parray, *oarray;
1927:   PetscInt        Nf = patch->nsubspaces, Np, poff, p, f;
1928:   PetscErrorCode  ierr;

1931:   if (patch->precomputeElementTensors) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Precomputing element tensors not implemented with DMPlex compute operator");
1932:   PCGetDM(pc, &dm);
1933:   DMConvert(dm, DMPLEX, &plex);
1934:   dm = plex;
1935:   DMGetDefaultSection(dm, &s);
1936:   /* Set offset into patch */
1937:   PetscSectionGetDof(patch->pointCounts, patchNum, &Np);
1938:   PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);
1939:   ISGetIndices(patch->points, &parray);
1940:   ISGetIndices(patch->offs,   &oarray);
1941:   for (f = 0; f < Nf; ++f) {
1942:     for (p = 0; p < Np; ++p) {
1943:       const PetscInt point = parray[poff+p];
1944:       PetscInt       dof;

1946:       PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);
1947:       PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);
1948:       if (patch->nsubspaces == 1) {PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);}
1949:       else                        {PetscSectionSetOffset(patch->patchSection, point, -1);}
1950:     }
1951:   }
1952:   ISRestoreIndices(patch->points, &parray);
1953:   ISRestoreIndices(patch->offs,   &oarray);
1954:   if (patch->viewSection) {ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);}
1955:   DMPlexComputeResidual_Patch_Internal(dm, patch->patchSection, cellIS, 0.0, x, NULL, F, ctx);
1956:   DMDestroy(&dm);
1957:   return(0);
1958: }

1960: PetscErrorCode PCPatchComputeFunction_Internal(PC pc, Vec x, Vec F, PetscInt point)
1961: {
1962:   PC_PATCH       *patch = (PC_PATCH *) pc->data;
1963:   const PetscInt *dofsArray;
1964:   const PetscInt *dofsArrayWithAll;
1965:   const PetscInt *cellsArray;
1966:   PetscInt        ncell, offset, pStart, pEnd;
1967:   PetscErrorCode  ierr;

1970:   PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);
1971:   if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n");
1972:   ISGetIndices(patch->dofs, &dofsArray);
1973:   ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);
1974:   ISGetIndices(patch->cells, &cellsArray);
1975:   PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);

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

1980:   PetscSectionGetDof(patch->cellCounts, point, &ncell);
1981:   PetscSectionGetOffset(patch->cellCounts, point, &offset);
1982:   if (ncell <= 0) {
1983:     PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
1984:     return(0);
1985:   }
1986:   VecSet(F, 0.0);
1987:   PetscStackPush("PCPatch user callback");
1988:   /* Cannot reuse the same IS because the geometry info is being cached in it */
1989:   ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);
1990:   patch->usercomputef(pc, point, x, F, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell,
1991:                                                                                             dofsArrayWithAll + offset*patch->totalDofsPerCell,
1992:                                                                                             patch->usercomputefctx);
1993:   PetscStackPop;
1994:   ISDestroy(&patch->cellIS);
1995:   ISRestoreIndices(patch->dofs, &dofsArray);
1996:   ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);
1997:   ISRestoreIndices(patch->cells, &cellsArray);
1998:   if (patch->viewMatrix) {
1999:     char name[PETSC_MAX_PATH_LEN];

2001:     PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch vector for Point %D", point);
2002:     PetscObjectSetName((PetscObject) F, name);
2003:     ObjectView((PetscObject) F, patch->viewerMatrix, patch->formatMatrix);
2004:   }
2005:   PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2006:   return(0);
2007: }

2009: 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)
2010: {
2011:   PC_PATCH       *patch = (PC_PATCH *) pc->data;
2012:   DM              dm, plex;
2013:   PetscSection    s;
2014:   const PetscInt *parray, *oarray;
2015:   PetscInt        Nf = patch->nsubspaces, Np, poff, p, f;
2016:   PetscErrorCode  ierr;

2019:   PCGetDM(pc, &dm);
2020:   DMConvert(dm, DMPLEX, &plex);
2021:   dm = plex;
2022:   DMGetDefaultSection(dm, &s);
2023:   /* Set offset into patch */
2024:   PetscSectionGetDof(patch->pointCounts, patchNum, &Np);
2025:   PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);
2026:   ISGetIndices(patch->points, &parray);
2027:   ISGetIndices(patch->offs,   &oarray);
2028:   for (f = 0; f < Nf; ++f) {
2029:     for (p = 0; p < Np; ++p) {
2030:       const PetscInt point = parray[poff+p];
2031:       PetscInt       dof;

2033:       PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);
2034:       PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);
2035:       if (patch->nsubspaces == 1) {PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);}
2036:       else                        {PetscSectionSetOffset(patch->patchSection, point, -1);}
2037:     }
2038:   }
2039:   ISRestoreIndices(patch->points, &parray);
2040:   ISRestoreIndices(patch->offs,   &oarray);
2041:   if (patch->viewSection) {ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);}
2042:   /* TODO Shut off MatViewFromOptions() in MatAssemblyEnd() here */
2043:   DMPlexComputeJacobian_Patch_Internal(dm, patch->patchSection, patch->patchSection, cellIS, 0.0, 0.0, x, NULL, J, J, ctx);
2044:   DMDestroy(&dm);
2045:   return(0);
2046: }

2048: /* This function zeros mat on entry */
2049: PetscErrorCode PCPatchComputeOperator_Internal(PC pc, Vec x, Mat mat, PetscInt point, PetscBool withArtificial)
2050: {
2051:   PC_PATCH       *patch = (PC_PATCH *) pc->data;
2052:   const PetscInt *dofsArray;
2053:   const PetscInt *dofsArrayWithAll = NULL;
2054:   const PetscInt *cellsArray;
2055:   PetscInt        ncell, offset, pStart, pEnd, numIntFacets, intFacetOffset;
2056:   PetscBool       isNonlinear;
2057:   PetscErrorCode  ierr;

2060:   PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);
2061:   isNonlinear = patch->isNonlinear;
2062:   if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n");
2063:   if(withArtificial) {
2064:     ISGetIndices(patch->dofsWithArtificial, &dofsArray);
2065:   } else {
2066:     ISGetIndices(patch->dofs, &dofsArray);
2067:   }
2068:   if (isNonlinear) {
2069:     ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);
2070:   }
2071:   ISGetIndices(patch->cells, &cellsArray);
2072:   PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);

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

2077:   PetscSectionGetDof(patch->cellCounts, point, &ncell);
2078:   PetscSectionGetOffset(patch->cellCounts, point, &offset);
2079:   if (ncell <= 0) {
2080:     PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2081:     return(0);
2082:   }
2083:   MatZeroEntries(mat);
2084:   if (patch->precomputeElementTensors) {
2085:     PetscInt           i;
2086:     PetscInt           ndof = patch->totalDofsPerCell;
2087:     const PetscScalar *elementTensors;

2089:     VecGetArrayRead(patch->cellMats, &elementTensors);
2090:     for (i = 0; i < ncell; i++) {
2091:       const PetscInt     cell = cellsArray[i + offset];
2092:       const PetscInt    *idx  = dofsArray + (offset + i)*ndof;
2093:       const PetscScalar *v    = elementTensors + patch->precomputedTensorLocations[cell]*ndof*ndof;
2094:       MatSetValues(mat, ndof, idx, ndof, idx, v, ADD_VALUES);
2095:     }
2096:     VecRestoreArrayRead(patch->cellMats, &elementTensors);
2097:     MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
2098:     MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
2099:   } else {
2100:     PetscStackPush("PCPatch user callback");
2101:     /* Cannot reuse the same IS because the geometry info is being cached in it */
2102:     ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);
2103:     patch->usercomputeop(pc, point, x, mat, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell, dofsArrayWithAll ? dofsArrayWithAll + offset*patch->totalDofsPerCell : NULL, patch->usercomputeopctx);
2104:   }
2105:   if (patch->usercomputeopintfacet) {
2106:     PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
2107:     PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
2108:     if (numIntFacets > 0) {
2109:       /* For each interior facet, grab the two cells (in local numbering, and concatenate dof numberings for those cells) */
2110:       PetscInt       *facetDofs = NULL, *facetDofsWithAll = NULL;
2111:       const PetscInt *intFacetsArray = NULL;
2112:       PetscInt        idx = 0;
2113:       PetscInt        i, c, d;
2114:       PetscInt        fStart;
2115:       DM              dm, plex;
2116:       IS              facetIS = NULL;
2117:       const PetscInt *facetCells = NULL;

2119:       ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
2120:       ISGetIndices(patch->intFacets, &intFacetsArray);
2121:       PCGetDM(pc, &dm);
2122:       DMConvert(dm, DMPLEX, &plex);
2123:       dm = plex;
2124:       DMPlexGetHeightStratum(dm, 1, &fStart, NULL);
2125:       /* FIXME: Pull this malloc out. */
2126:       PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofs);
2127:       if (dofsArrayWithAll) {
2128:         PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofsWithAll);
2129:       }
2130:       if (patch->precomputeElementTensors) {
2131:         PetscInt           nFacetDof = 2*patch->totalDofsPerCell;
2132:         const PetscScalar *elementTensors;

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

2136:         for (i = 0; i < numIntFacets; i++) {
2137:           const PetscInt     facet = intFacetsArray[i + intFacetOffset];
2138:           const PetscScalar *v     = elementTensors + patch->precomputedIntFacetTensorLocations[facet - fStart]*nFacetDof*nFacetDof;
2139:           idx = 0;
2140:           /*
2141:            * 0--1
2142:            * |\-|
2143:            * |+\|
2144:            * 2--3
2145:            * [0, 2, 3, 0, 1, 3]
2146:            */
2147:           for (c = 0; c < 2; c++) {
2148:             const PetscInt cell = facetCells[2*(intFacetOffset + i) + c];
2149:             for (d = 0; d < patch->totalDofsPerCell; d++) {
2150:               facetDofs[idx] = dofsArray[(offset + cell)*patch->totalDofsPerCell + d];
2151:               idx++;
2152:             }
2153:           }
2154:           MatSetValues(mat, nFacetDof, facetDofs, nFacetDof, facetDofs, v, ADD_VALUES);
2155:         }
2156:         VecRestoreArrayRead(patch->intFacetMats, &elementTensors);
2157:       } else {
2158:         /*
2159:          * 0--1
2160:          * |\-|
2161:          * |+\|
2162:          * 2--3
2163:          * [0, 2, 3, 0, 1, 3]
2164:          */
2165:         for (i = 0; i < numIntFacets; i++) {
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:               if (dofsArrayWithAll) {
2171:                 facetDofsWithAll[idx] = dofsArrayWithAll[(offset + cell)*patch->totalDofsPerCell + d];
2172:               }
2173:               idx++;
2174:             }
2175:           }
2176:         }
2177:         ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray + intFacetOffset, PETSC_USE_POINTER, &facetIS);
2178:         patch->usercomputeopintfacet(pc, point, x, mat, facetIS, 2*numIntFacets*patch->totalDofsPerCell, facetDofs, facetDofsWithAll, patch->usercomputeopintfacetctx);
2179:         ISDestroy(&facetIS);
2180:       }
2181:       ISRestoreIndices(patch->intFacetsToPatchCell, &facetCells);
2182:       ISRestoreIndices(patch->intFacets, &intFacetsArray);
2183:       PetscFree(facetDofs);
2184:       PetscFree(facetDofsWithAll);
2185:       DMDestroy(&dm);
2186:     }
2187:   }

2189:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
2190:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);

2192:   PetscStackPop;
2193:   ISDestroy(&patch->cellIS);
2194:   if(withArtificial) {
2195:     ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);
2196:   } else {
2197:     ISRestoreIndices(patch->dofs, &dofsArray);
2198:   }
2199:   if (isNonlinear) {
2200:     ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);
2201:   }
2202:   ISRestoreIndices(patch->cells, &cellsArray);
2203:   if (patch->viewMatrix) {
2204:     char name[PETSC_MAX_PATH_LEN];

2206:     PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch matrix for Point %D", point);
2207:     PetscObjectSetName((PetscObject) mat, name);
2208:     ObjectView((PetscObject) mat, patch->viewerMatrix, patch->formatMatrix);
2209:   }
2210:   PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2211:   return(0);
2212: }

2214: static PetscErrorCode MatSetValues_PCPatch_Private(Mat mat, PetscInt m, const PetscInt idxm[],
2215:                                                    PetscInt n, const PetscInt idxn[], const PetscScalar *v, InsertMode addv)
2216: {
2217:   Vec            data;
2218:   PetscScalar   *array;
2219:   PetscInt       bs, nz, i, j, cell;

2222:   MatShellGetContext(mat, &data);
2223:   VecGetBlockSize(data, &bs);
2224:   VecGetSize(data, &nz);
2225:   VecGetArray(data, &array);
2226:   if (m != n) SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Only for square insertion");
2227:   cell = (PetscInt)(idxm[0]/bs); /* use the fact that this is called once per cell */
2228:   for (i = 0; i < m; i++) {
2229:     if (idxm[i] != idxn[i]) SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Row and column indices must match!");
2230:     for (j = 0; j < n; j++) {
2231:       const PetscScalar v_ = v[i*bs + j];
2232:       /* Indexing is special to the data structure we have! */
2233:       if (addv == INSERT_VALUES) {
2234:         array[cell*bs*bs + i*bs + j] = v_;
2235:       } else {
2236:         array[cell*bs*bs + i*bs + j] += v_;
2237:       }
2238:     }
2239:   }
2240:   VecRestoreArray(data, &array);
2241:   return(0);
2242: }

2244: static PetscErrorCode PCPatchPrecomputePatchTensors_Private(PC pc)
2245: {
2246:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
2247:   const PetscInt *cellsArray;
2248:   PetscInt        ncell, offset;
2249:   const PetscInt *dofMapArray;
2250:   PetscInt        i, j;
2251:   IS              dofMap;
2252:   IS              cellIS;
2253:   const PetscInt  ndof  = patch->totalDofsPerCell;
2254:   PetscErrorCode  ierr;
2255:   Mat             vecMat;
2256:   PetscInt        cStart, cEnd;
2257:   DM              dm, plex;


2260:   ISGetSize(patch->cells, &ncell);
2261:   if (!ncell) { /* No cells to assemble over -> skip */
2262:     return(0);
2263:   }

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

2267:   PCGetDM(pc, &dm);
2268:   DMConvert(dm, DMPLEX, &plex);
2269:   dm = plex;
2270:   if (!patch->allCells) {
2271:     PetscHSetI      cells;
2272:     PetscHashIter   hi;
2273:     PetscInt pStart, pEnd;
2274:     PetscInt *allCells = NULL;
2275:     PetscHSetICreate(&cells);
2276:     ISGetIndices(patch->cells, &cellsArray);
2277:     PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);
2278:     for (i = pStart; i < pEnd; i++) {
2279:       PetscSectionGetDof(patch->cellCounts, i, &ncell);
2280:       PetscSectionGetOffset(patch->cellCounts, i, &offset);
2281:       if (ncell <= 0) continue;
2282:       for (j = 0; j < ncell; j++) {
2283:         PetscHSetIAdd(cells, cellsArray[offset + j]);
2284:       }
2285:     }
2286:     ISRestoreIndices(patch->cells, &cellsArray);
2287:     PetscHSetIGetSize(cells, &ncell);
2288:     PetscMalloc1(ncell, &allCells);
2289:     DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2290:     PetscMalloc1(cEnd-cStart, &patch->precomputedTensorLocations);
2291:     i = 0;
2292:     PetscHashIterBegin(cells, hi);
2293:     while (!PetscHashIterAtEnd(cells, hi)) {
2294:       PetscHashIterGetKey(cells, hi, allCells[i]);
2295:       patch->precomputedTensorLocations[allCells[i]] = i;
2296:       PetscHashIterNext(cells, hi);
2297:       i++;
2298:     }
2299:     PetscHSetIDestroy(&cells);
2300:     ISCreateGeneral(PETSC_COMM_SELF, ncell, allCells, PETSC_OWN_POINTER, &patch->allCells);
2301:   }
2302:   ISGetSize(patch->allCells, &ncell);
2303:   if (!patch->cellMats) {
2304:     VecCreateSeq(PETSC_COMM_SELF, ncell*ndof*ndof, &patch->cellMats);
2305:     VecSetBlockSize(patch->cellMats, ndof);
2306:   }
2307:   VecSet(patch->cellMats, 0);

2309:   MatCreateShell(PETSC_COMM_SELF, ncell*ndof, ncell*ndof, ncell*ndof, ncell*ndof,
2310:                         (void*)patch->cellMats, &vecMat);
2311:   MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void(*)(void))&MatSetValues_PCPatch_Private);
2312:   ISGetSize(patch->allCells, &ncell);
2313:   ISCreateStride(PETSC_COMM_SELF, ndof*ncell, 0, 1, &dofMap);
2314:   ISGetIndices(dofMap, &dofMapArray);
2315:   ISGetIndices(patch->allCells, &cellsArray);
2316:   ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray, PETSC_USE_POINTER, &cellIS);
2317:   PetscStackPush("PCPatch user callback");
2318:   /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2319:   patch->usercomputeop(pc, -1, NULL, vecMat, cellIS, ndof*ncell, dofMapArray, NULL, patch->usercomputeopctx);
2320:   PetscStackPop;
2321:   ISDestroy(&cellIS);
2322:   MatDestroy(&vecMat);
2323:   ISRestoreIndices(patch->allCells, &cellsArray);
2324:   ISRestoreIndices(dofMap, &dofMapArray);
2325:   ISDestroy(&dofMap);

2327:   if (patch->usercomputeopintfacet) {
2328:     PetscInt nIntFacets;
2329:     IS       intFacetsIS;
2330:     const PetscInt *intFacetsArray = NULL;
2331:     if (!patch->allIntFacets) {
2332:       PetscHSetI      facets;
2333:       PetscHashIter   hi;
2334:       PetscInt pStart, pEnd, fStart, fEnd;
2335:       PetscInt *allIntFacets = NULL;
2336:       PetscHSetICreate(&facets);
2337:       ISGetIndices(patch->intFacets, &intFacetsArray);
2338:       PetscSectionGetChart(patch->intFacetCounts, &pStart, &pEnd);
2339:       DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2340:       for (i = pStart; i < pEnd; i++) {
2341:         PetscSectionGetDof(patch->intFacetCounts, i, &nIntFacets);
2342:         PetscSectionGetOffset(patch->intFacetCounts, i, &offset);
2343:         if (nIntFacets <= 0) continue;
2344:         for (j = 0; j < nIntFacets; j++) {
2345:           PetscHSetIAdd(facets, intFacetsArray[offset + j]);
2346:         }
2347:       }
2348:       ISRestoreIndices(patch->intFacets, &intFacetsArray);
2349:       PetscHSetIGetSize(facets, &nIntFacets);
2350:       PetscMalloc1(nIntFacets, &allIntFacets);
2351:       PetscMalloc1(fEnd-fStart, &patch->precomputedIntFacetTensorLocations);
2352:       i = 0;
2353:       PetscHashIterBegin(facets, hi);
2354:       while (!PetscHashIterAtEnd(facets, hi)) {
2355:         PetscHashIterGetKey(facets, hi, allIntFacets[i]);
2356:         patch->precomputedIntFacetTensorLocations[allIntFacets[i] - fStart] = i;
2357:         PetscHashIterNext(facets, hi);
2358:         i++;
2359:       }
2360:       PetscHSetIDestroy(&facets);
2361:       ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, allIntFacets, PETSC_OWN_POINTER, &patch->allIntFacets);
2362:     }
2363:     ISGetSize(patch->allIntFacets, &nIntFacets);
2364:     if (!patch->intFacetMats) {
2365:       VecCreateSeq(PETSC_COMM_SELF, nIntFacets*ndof*ndof*4, &patch->intFacetMats);
2366:       VecSetBlockSize(patch->intFacetMats, ndof*2);
2367:     }
2368:     VecSet(patch->intFacetMats, 0);

2370:     MatCreateShell(PETSC_COMM_SELF, nIntFacets*ndof*2, nIntFacets*ndof*2, nIntFacets*ndof*2, nIntFacets*ndof*2,
2371:                           (void*)patch->intFacetMats, &vecMat);
2372:     MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void(*)(void))&MatSetValues_PCPatch_Private);
2373:     ISCreateStride(PETSC_COMM_SELF, 2*ndof*nIntFacets, 0, 1, &dofMap);
2374:     ISGetIndices(dofMap, &dofMapArray);
2375:     ISGetIndices(patch->allIntFacets, &intFacetsArray);
2376:     ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, intFacetsArray, PETSC_USE_POINTER, &intFacetsIS);
2377:     PetscStackPush("PCPatch user callback (interior facets)");
2378:     /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2379:     patch->usercomputeopintfacet(pc, -1, NULL, vecMat, intFacetsIS, 2*ndof*nIntFacets, dofMapArray, NULL, patch->usercomputeopintfacetctx);
2380:     PetscStackPop;
2381:     ISDestroy(&intFacetsIS);
2382:     MatDestroy(&vecMat);
2383:     ISRestoreIndices(patch->allIntFacets, &intFacetsArray);
2384:     ISRestoreIndices(dofMap, &dofMapArray);
2385:     ISDestroy(&dofMap);
2386:   }
2387:   DMDestroy(&dm);
2388:   PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);

2390:   return(0);
2391: }

2393: PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat, PatchScatterType scattertype)
2394: {
2395:   PC_PATCH          *patch     = (PC_PATCH *) pc->data;
2396:   const PetscScalar *xArray    = NULL;
2397:   PetscScalar       *yArray    = NULL;
2398:   const PetscInt    *gtolArray = NULL;
2399:   PetscInt           dof, offset, lidx;
2400:   PetscErrorCode     ierr;

2403:   PetscLogEventBegin(PC_Patch_Scatter, pc, 0, 0, 0);
2404:   VecGetArrayRead(x, &xArray);
2405:   VecGetArray(y, &yArray);
2406:   if (scattertype == SCATTER_WITHARTIFICIAL) {
2407:     PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);
2408:     PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offset);
2409:     ISGetIndices(patch->gtolWithArtificial, &gtolArray);
2410:   } else if (scattertype == SCATTER_WITHALL) {
2411:     PetscSectionGetDof(patch->gtolCountsWithAll, p, &dof);
2412:     PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offset);
2413:     ISGetIndices(patch->gtolWithAll, &gtolArray);
2414:   } else {
2415:     PetscSectionGetDof(patch->gtolCounts, p, &dof);
2416:     PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2417:     ISGetIndices(patch->gtol, &gtolArray);
2418:   }
2419:   if (mode == INSERT_VALUES && scat != SCATTER_FORWARD) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't insert if not scattering forward\n");
2420:   if (mode == ADD_VALUES    && scat != SCATTER_REVERSE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't add if not scattering reverse\n");
2421:   for (lidx = 0; lidx < dof; ++lidx) {
2422:     const PetscInt gidx = gtolArray[offset+lidx];

2424:     if (mode == INSERT_VALUES) yArray[lidx]  = xArray[gidx]; /* Forward */
2425:     else                       yArray[gidx] += xArray[lidx]; /* Reverse */
2426:   }
2427:   if (scattertype == SCATTER_WITHARTIFICIAL) {
2428:     ISRestoreIndices(patch->gtolWithArtificial, &gtolArray);
2429:   } else if (scattertype == SCATTER_WITHALL) {
2430:     ISRestoreIndices(patch->gtolWithAll, &gtolArray);
2431:   } else {
2432:     ISRestoreIndices(patch->gtol, &gtolArray);
2433:   }
2434:   VecRestoreArrayRead(x, &xArray);
2435:   VecRestoreArray(y, &yArray);
2436:   PetscLogEventEnd(PC_Patch_Scatter, pc, 0, 0, 0);
2437:   return(0);
2438: }

2440: static PetscErrorCode PCSetUp_PATCH_Linear(PC pc)
2441: {
2442:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2443:   const char    *prefix;
2444:   PetscInt       i;

2448:   if (!pc->setupcalled) {
2449:     PetscMalloc1(patch->npatch, &patch->solver);
2450:     PCGetOptionsPrefix(pc, &prefix);
2451:     for (i = 0; i < patch->npatch; ++i) {
2452:       KSP ksp;
2453:       PC  subpc;

2455:       KSPCreate(PETSC_COMM_SELF, &ksp);
2456:       KSPSetErrorIfNotConverged(ksp, pc->erroriffailure);
2457:       KSPSetOptionsPrefix(ksp, prefix);
2458:       KSPAppendOptionsPrefix(ksp, "sub_");
2459:       PetscObjectIncrementTabLevel((PetscObject) ksp, (PetscObject) pc, 1);
2460:       KSPGetPC(ksp, &subpc);
2461:       PetscObjectIncrementTabLevel((PetscObject) subpc, (PetscObject) pc, 1);
2462:       PetscLogObjectParent((PetscObject) pc, (PetscObject) ksp);
2463:       patch->solver[i] = (PetscObject) ksp;
2464:     }
2465:   }
2466:   if (patch->save_operators) {
2467:     if (patch->precomputeElementTensors) {
2468:       PCPatchPrecomputePatchTensors_Private(pc);
2469:     }
2470:     for (i = 0; i < patch->npatch; ++i) {
2471:       PCPatchComputeOperator_Internal(pc, NULL, patch->mat[i], i, PETSC_FALSE);
2472:       KSPSetOperators((KSP) patch->solver[i], patch->mat[i], patch->mat[i]);
2473:     }
2474:   }
2475:   if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2476:     for (i = 0; i < patch->npatch; ++i) {
2477:       /* Instead of padding patch->patchUpdate with zeros to get */
2478:       /* patch->patchUpdateWithArtificial and then multiplying with the matrix, */
2479:       /* just get rid of the columns that correspond to the dofs with */
2480:       /* artificial bcs. That's of course fairly inefficient, hopefully we */
2481:       /* can just assemble the rectangular matrix in the first place. */
2482:       Mat matSquare;
2483:       IS rowis;
2484:       PetscInt dof;

2486:       MatGetSize(patch->mat[i], &dof, NULL);
2487:       if (dof == 0) {
2488:         patch->matWithArtificial[i] = NULL;
2489:         continue;
2490:       }

2492:       PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);
2493:       PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);

2495:       MatGetSize(matSquare, &dof, NULL);
2496:       ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis);
2497:       if(pc->setupcalled) {
2498:         MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_REUSE_MATRIX, &patch->matWithArtificial[i]);
2499:       } else {
2500:         MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &patch->matWithArtificial[i]);
2501:       }
2502:       ISDestroy(&rowis);
2503:       MatDestroy(&matSquare);
2504:     }
2505:   }
2506:   return(0);
2507: }

2509: static PetscErrorCode PCSetUp_PATCH(PC pc)
2510: {
2511:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2512:   PetscInt       i;
2513:   PetscBool       isNonlinear;

2517:   if (!pc->setupcalled) {
2518:     PetscInt pStart, pEnd, p;
2519:     PetscInt localSize;

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

2523:     isNonlinear = patch->isNonlinear;
2524:     if (!patch->nsubspaces) {
2525:       DM           dm, plex;
2526:       PetscSection s;
2527:       PetscInt     cStart, cEnd, c, Nf, f, numGlobalBcs = 0, *globalBcs, *Nb, totNb = 0, **cellDofs;

2529:       PCGetDM(pc, &dm);
2530:       if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "Must set DM for PCPATCH or call PCPatchSetDiscretisationInfo()");
2531:       DMConvert(dm, DMPLEX, &plex);
2532:       dm = plex;
2533:       DMGetDefaultSection(dm, &s);
2534:       PetscSectionGetNumFields(s, &Nf);
2535:       PetscSectionGetChart(s, &pStart, &pEnd);
2536:       for (p = pStart; p < pEnd; ++p) {
2537:         PetscInt cdof;
2538:         PetscSectionGetConstraintDof(s, p, &cdof);
2539:         numGlobalBcs += cdof;
2540:       }
2541:       DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2542:       PetscMalloc3(Nf, &Nb, Nf, &cellDofs, numGlobalBcs, &globalBcs);
2543:       for (f = 0; f < Nf; ++f) {
2544:         PetscFE        fe;
2545:         PetscDualSpace sp;
2546:         PetscInt       cdoff = 0;

2548:         DMGetField(dm, f, NULL, (PetscObject *) &fe);
2549:         /* PetscFEGetNumComponents(fe, &Nc[f]); */
2550:         PetscFEGetDualSpace(fe, &sp);
2551:         PetscDualSpaceGetDimension(sp, &Nb[f]);
2552:         totNb += Nb[f];

2554:         PetscMalloc1((cEnd-cStart)*Nb[f], &cellDofs[f]);
2555:         for (c = cStart; c < cEnd; ++c) {
2556:           PetscInt *closure = NULL;
2557:           PetscInt  clSize  = 0, cl;

2559:           DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
2560:           for (cl = 0; cl < clSize*2; cl += 2) {
2561:             const PetscInt p = closure[cl];
2562:             PetscInt       fdof, d, foff;

2564:             PetscSectionGetFieldDof(s, p, f, &fdof);
2565:             PetscSectionGetFieldOffset(s, p, f, &foff);
2566:             for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d;
2567:           }
2568:           DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
2569:         }
2570:         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]);
2571:       }
2572:       numGlobalBcs = 0;
2573:       for (p = pStart; p < pEnd; ++p) {
2574:         const PetscInt *ind;
2575:         PetscInt        off, cdof, d;

2577:         PetscSectionGetOffset(s, p, &off);
2578:         PetscSectionGetConstraintDof(s, p, &cdof);
2579:         PetscSectionGetConstraintIndices(s, p, &ind);
2580:         for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d];
2581:       }

2583:       PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **) cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs);
2584:       for (f = 0; f < Nf; ++f) {
2585:         PetscFree(cellDofs[f]);
2586:       }
2587:       PetscFree3(Nb, cellDofs, globalBcs);
2588:       PCPatchSetComputeFunction(pc, PCPatchComputeFunction_DMPlex_Private, NULL);
2589:       PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL);
2590:       DMDestroy(&dm);
2591:     }

2593:     localSize = patch->subspaceOffsets[patch->nsubspaces];
2594:     VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localRHS);
2595:     VecSetUp(patch->localRHS);
2596:     VecDuplicate(patch->localRHS, &patch->localUpdate);
2597:     PCPatchCreateCellPatches(pc);
2598:     PCPatchCreateCellPatchDiscretisationInfo(pc);

2600:     /* OK, now build the work vectors */
2601:     PetscSectionGetChart(patch->gtolCounts, &pStart, &pEnd);
2602:     PetscMalloc1(patch->npatch, &patch->patchRHS);
2603:     PetscMalloc1(patch->npatch, &patch->patchUpdate);

2605:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2606:       PetscMalloc1(patch->npatch, &patch->patchRHSWithArtificial);
2607:       PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithArtificial);
2608:     }
2609:     if (isNonlinear) {
2610:       PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithAll);
2611:     }
2612:     for (p = pStart; p < pEnd; ++p) {
2613:       PetscInt dof;

2615:       PetscSectionGetDof(patch->gtolCounts, p, &dof);
2616:       VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchRHS[p-pStart]);
2617:       VecSetUp(patch->patchRHS[p-pStart]);
2618:       VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchUpdate[p-pStart]);
2619:       VecSetUp(patch->patchUpdate[p-pStart]);
2620:       if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2621:         const PetscInt    *gtolArray, *gtolArrayWithArtificial = NULL;
2622:         PetscInt           numPatchDofs, offset;
2623:         PetscInt           numPatchDofsWithArtificial, offsetWithArtificial;
2624:         PetscInt           dofWithoutArtificialCounter = 0;
2625:         PetscInt          *patchWithoutArtificialToWithArtificialArray;

2627:         PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);
2628:         VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchRHSWithArtificial[p-pStart]);
2629:         VecSetUp(patch->patchRHSWithArtificial[p-pStart]);

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

2635:         PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);
2636:         if (numPatchDofs == 0) {
2637:           patch->dofMappingWithoutToWithArtificial[p-pStart] = NULL;
2638:           continue;
2639:         }

2641:         PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2642:         ISGetIndices(patch->gtolWithArtificial, &gtolArrayWithArtificial);
2643:         PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &numPatchDofsWithArtificial);
2644:         PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offsetWithArtificial);

2646:         PetscMalloc1(numPatchDofs, &patchWithoutArtificialToWithArtificialArray);
2647:         for (i=0; i<numPatchDofsWithArtificial; i++) {
2648:           if (gtolArrayWithArtificial[i+offsetWithArtificial] == gtolArray[offset+dofWithoutArtificialCounter]) {
2649:             patchWithoutArtificialToWithArtificialArray[dofWithoutArtificialCounter] = i;
2650:             dofWithoutArtificialCounter++;
2651:             if (dofWithoutArtificialCounter == numPatchDofs)
2652:               break;
2653:           }
2654:         }
2655:         ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutArtificialToWithArtificialArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithArtificial[p-pStart]);
2656:         ISRestoreIndices(patch->gtol, &gtolArray);
2657:         ISRestoreIndices(patch->gtolWithArtificial, &gtolArrayWithArtificial);
2658:       }
2659:       if (isNonlinear) {
2660:         const PetscInt    *gtolArray, *gtolArrayWithAll = NULL;
2661:         PetscInt           numPatchDofs, offset;
2662:         PetscInt           numPatchDofsWithAll, offsetWithAll;
2663:         PetscInt           dofWithoutAllCounter = 0;
2664:         PetscInt          *patchWithoutAllToWithAllArray;

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

2670:         PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);
2671:         if (numPatchDofs == 0) {
2672:           patch->dofMappingWithoutToWithAll[p-pStart] = NULL;
2673:           continue;
2674:         }

2676:         PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2677:         ISGetIndices(patch->gtolWithAll, &gtolArrayWithAll);
2678:         PetscSectionGetDof(patch->gtolCountsWithAll, p, &numPatchDofsWithAll);
2679:         PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offsetWithAll);

2681:         PetscMalloc1(numPatchDofs, &patchWithoutAllToWithAllArray);

2683:         for (i=0; i<numPatchDofsWithAll; i++) {
2684:           if (gtolArrayWithAll[i+offsetWithAll] == gtolArray[offset+dofWithoutAllCounter]) {
2685:             patchWithoutAllToWithAllArray[dofWithoutAllCounter] = i;
2686:             dofWithoutAllCounter++;
2687:             if (dofWithoutAllCounter == numPatchDofs)
2688:               break;
2689:           }
2690:         }
2691:         ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutAllToWithAllArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithAll[p-pStart]);
2692:         ISRestoreIndices(patch->gtol, &gtolArray);
2693:         ISRestoreIndices(patch->gtolWithAll, &gtolArrayWithAll);
2694:       }
2695:     }
2696:     if (patch->save_operators) {
2697:       PetscMalloc1(patch->npatch, &patch->mat);
2698:       for (i = 0; i < patch->npatch; ++i) {
2699:         PCPatchCreateMatrix_Private(pc, i, &patch->mat[i], PETSC_FALSE);
2700:       }
2701:     }
2702:     PetscLogEventEnd(PC_Patch_CreatePatches, pc, 0, 0, 0);

2704:     /* If desired, calculate weights for dof multiplicity */
2705:     if (patch->partition_of_unity) {
2706:       PetscScalar *input = NULL;
2707:       PetscScalar *output = NULL;
2708:       Vec global;

2710:       VecDuplicate(patch->localRHS, &patch->dof_weights);
2711:       if(patch->local_composition_type == PC_COMPOSITE_ADDITIVE) {
2712:         for (i = 0; i < patch->npatch; ++i) {
2713:           PetscInt dof;

2715:           PetscSectionGetDof(patch->gtolCounts, i+pStart, &dof);
2716:           if (dof <= 0) continue;
2717:           VecSet(patch->patchRHS[i], 1.0);
2718:           PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchRHS[i], patch->dof_weights, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);
2719:         }
2720:       } else {
2721:         /* multiplicative is actually only locally multiplicative and globally additive. need the pou where the mesh decomposition overlaps */
2722:         VecSet(patch->dof_weights, 1.0);
2723:       }

2725:       VecDuplicate(patch->dof_weights, &global);
2726:       VecSet(global, 0.);

2728:       VecGetArray(patch->dof_weights, &input);
2729:       VecGetArray(global, &output);
2730:       PetscSFReduceBegin(patch->defaultSF, MPIU_SCALAR, input, output, MPI_SUM);
2731:       PetscSFReduceEnd(patch->defaultSF, MPIU_SCALAR, input, output, MPI_SUM);
2732:       VecRestoreArray(patch->dof_weights, &input);
2733:       VecRestoreArray(global, &output);

2735:       VecReciprocal(global);

2737:       VecGetArray(patch->dof_weights, &output);
2738:       VecGetArray(global, &input);
2739:       PetscSFBcastBegin(patch->defaultSF, MPIU_SCALAR, input, output);
2740:       PetscSFBcastEnd(patch->defaultSF, MPIU_SCALAR, input, output);
2741:       VecRestoreArray(patch->dof_weights, &output);
2742:       VecRestoreArray(global, &input);
2743:       VecDestroy(&global);
2744:     }
2745:     if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE && patch->save_operators) {
2746:       PetscMalloc1(patch->npatch, &patch->matWithArtificial);
2747:     }
2748:   }
2749:   (*patch->setupsolver)(pc);
2750:   return(0);
2751: }

2753: static PetscErrorCode PCApply_PATCH_Linear(PC pc, PetscInt i, Vec x, Vec y)
2754: {
2755:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2756:   KSP            ksp   = (KSP) patch->solver[i];

2760:   if (!patch->save_operators) {
2761:     Mat mat;

2763:     PCPatchCreateMatrix_Private(pc, i, &mat, PETSC_FALSE);
2764:     /* Populate operator here. */
2765:     PCPatchComputeOperator_Internal(pc, NULL, mat, i, PETSC_FALSE);
2766:     KSPSetOperators(ksp, mat, mat);
2767:     /* Drop reference so the KSPSetOperators below will blow it away. */
2768:     MatDestroy(&mat);
2769:   }
2770:   PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);
2771:   if (!ksp->setfromoptionscalled) {
2772:     KSPSetFromOptions(ksp);
2773:   }
2774:   KSPSolve(ksp, x, y);
2775:   KSPCheckSolve(ksp, pc, y);
2776:   PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);
2777:   if (!patch->save_operators) {
2778:     PC pc;
2779:     KSPSetOperators(ksp, NULL, NULL);
2780:     KSPGetPC(ksp, &pc);
2781:     /* Destroy PC context too, otherwise the factored matrix hangs around. */
2782:     PCReset(pc);
2783:   }
2784:   return(0);
2785: }

2787: static PetscErrorCode PCUpdateMultiplicative_PATCH_Linear(PC pc, PetscInt i, PetscInt pStart)
2788: {
2789:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2790:   Mat multMat;


2795:   if (patch->save_operators) {
2796:     multMat = patch->matWithArtificial[i];
2797:   } else {
2798:     /*Very inefficient, hopefully we can just assemble the rectangular matrix in the first place.*/
2799:     Mat matSquare;
2800:     PetscInt dof;
2801:     IS rowis;
2802:     PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);
2803:     PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);
2804:     MatGetSize(matSquare, &dof, NULL);
2805:     ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis);
2806:     MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &multMat);
2807:     MatDestroy(&matSquare);
2808:     ISDestroy(&rowis);
2809:   }
2810:   MatMult(multMat, patch->patchUpdate[i], patch->patchRHSWithArtificial[i]);
2811:   VecScale(patch->patchRHSWithArtificial[i], -1.0);
2812:   PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHSWithArtificial[i], patch->localRHS, ADD_VALUES, SCATTER_REVERSE, SCATTER_WITHARTIFICIAL);
2813:   if (!patch->save_operators) {
2814:     MatDestroy(&multMat);
2815:   }
2816:   return(0);
2817: }

2819: static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y)
2820: {
2821:   PC_PATCH          *patch    = (PC_PATCH *) pc->data;
2822:   const PetscScalar *globalRHS  = NULL;
2823:   PetscScalar       *localRHS   = NULL;
2824:   PetscScalar       *globalUpdate  = NULL;
2825:   const PetscInt    *bcNodes  = NULL;
2826:   PetscInt           nsweep   = patch->symmetrise_sweep ? 2 : 1;
2827:   PetscInt           start[2] = {0, 0};
2828:   PetscInt           end[2]   = {-1, -1};
2829:   const PetscInt     inc[2]   = {1, -1};
2830:   const PetscScalar *localUpdate;
2831:   const PetscInt    *iterationSet;
2832:   PetscInt           pStart, numBcs, n, sweep, bc, j;
2833:   PetscErrorCode     ierr;

2836:   PetscLogEventBegin(PC_Patch_Apply, pc, 0, 0, 0);
2837:   PetscOptionsPushGetViewerOff(PETSC_TRUE);
2838:   /* start, end, inc have 2 entries to manage a second backward sweep if we symmetrize */
2839:   end[0]   = patch->npatch;
2840:   start[1] = patch->npatch-1;
2841:   if (patch->user_patches) {
2842:     ISGetLocalSize(patch->iterationSet, &end[0]);
2843:     start[1] = end[0] - 1;
2844:     ISGetIndices(patch->iterationSet, &iterationSet);
2845:   }
2846:   /* Scatter from global space into overlapped local spaces */
2847:   VecGetArrayRead(x, &globalRHS);
2848:   VecGetArray(patch->localRHS, &localRHS);
2849:   PetscSFBcastBegin(patch->defaultSF, MPIU_SCALAR, globalRHS, localRHS);
2850:   PetscSFBcastEnd(patch->defaultSF, MPIU_SCALAR, globalRHS, localRHS);
2851:   VecRestoreArrayRead(x, &globalRHS);
2852:   VecRestoreArray(patch->localRHS, &localRHS);

2854:   VecSet(patch->localUpdate, 0.0);
2855:   PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);
2856:   for (sweep = 0; sweep < nsweep; sweep++) {
2857:     for (j = start[sweep]; j*inc[sweep] < end[sweep]*inc[sweep]; j += inc[sweep]) {
2858:       PetscInt i       = patch->user_patches ? iterationSet[j] : j;
2859:       PetscInt start, len;

2861:       PetscSectionGetDof(patch->gtolCounts, i+pStart, &len);
2862:       PetscSectionGetOffset(patch->gtolCounts, i+pStart, &start);
2863:       /* TODO: Squash out these guys in the setup as well. */
2864:       if (len <= 0) continue;
2865:       /* TODO: Do we need different scatters for X and Y? */
2866:       PCPatch_ScatterLocal_Private(pc, i+pStart, patch->localRHS, patch->patchRHS[i], INSERT_VALUES, SCATTER_FORWARD, SCATTER_INTERIOR);
2867:       (*patch->applysolver)(pc, i, patch->patchRHS[i], patch->patchUpdate[i]);
2868:       PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchUpdate[i], patch->localUpdate, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);
2869:       if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2870:         (*patch->updatemultiplicative)(pc, i, pStart);
2871:       }
2872:     }
2873:   }
2874:   if (patch->user_patches) {ISRestoreIndices(patch->iterationSet, &iterationSet);}
2875:   /* XXX: should we do this on the global vector? */
2876:   if (patch->partition_of_unity) {
2877:     VecPointwiseMult(patch->localUpdate, patch->localUpdate, patch->dof_weights);
2878:   }
2879:   /* Now patch->localUpdate contains the solution of the patch solves, so we need to combine them all. */
2880:   VecSet(y, 0.0);
2881:   VecGetArray(y, &globalUpdate);
2882:   VecGetArrayRead(patch->localUpdate, &localUpdate);
2883:   PetscSFReduceBegin(patch->defaultSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);
2884:   PetscSFReduceEnd(patch->defaultSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);
2885:   VecRestoreArrayRead(patch->localUpdate, &localUpdate);

2887:   /* Now we need to send the global BC values through */
2888:   VecGetArrayRead(x, &globalRHS);
2889:   ISGetSize(patch->globalBcNodes, &numBcs);
2890:   ISGetIndices(patch->globalBcNodes, &bcNodes);
2891:   VecGetLocalSize(x, &n);
2892:   for (bc = 0; bc < numBcs; ++bc) {
2893:     const PetscInt idx = bcNodes[bc];
2894:     if (idx < n) globalUpdate[idx] = globalRHS[idx];
2895:   }

2897:   ISRestoreIndices(patch->globalBcNodes, &bcNodes);
2898:   VecRestoreArrayRead(x, &globalRHS);
2899:   VecRestoreArray(y, &globalUpdate);

2901:   PetscOptionsPopGetViewerOff();
2902:   PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0);
2903:   return(0);
2904: }

2906: static PetscErrorCode PCReset_PATCH_Linear(PC pc)
2907: {
2908:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2909:   PetscInt       i;

2913:   if (patch->solver) {
2914:     for (i = 0; i < patch->npatch; ++i) {KSPReset((KSP) patch->solver[i]);}
2915:   }
2916:   return(0);
2917: }

2919: static PetscErrorCode PCReset_PATCH(PC pc)
2920: {
2921:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2922:   PetscInt       i;


2927:   PetscSFDestroy(&patch->defaultSF);
2928:   PetscSectionDestroy(&patch->cellCounts);
2929:   PetscSectionDestroy(&patch->pointCounts);
2930:   PetscSectionDestroy(&patch->cellNumbering);
2931:   PetscSectionDestroy(&patch->gtolCounts);
2932:   ISDestroy(&patch->gtol);
2933:   ISDestroy(&patch->cells);
2934:   ISDestroy(&patch->points);
2935:   ISDestroy(&patch->dofs);
2936:   ISDestroy(&patch->offs);
2937:   PetscSectionDestroy(&patch->patchSection);
2938:   ISDestroy(&patch->ghostBcNodes);
2939:   ISDestroy(&patch->globalBcNodes);
2940:   PetscSectionDestroy(&patch->gtolCountsWithArtificial);
2941:   ISDestroy(&patch->gtolWithArtificial);
2942:   ISDestroy(&patch->dofsWithArtificial);
2943:   ISDestroy(&patch->offsWithArtificial);
2944:   PetscSectionDestroy(&patch->gtolCountsWithAll);
2945:   ISDestroy(&patch->gtolWithAll);
2946:   ISDestroy(&patch->dofsWithAll);
2947:   ISDestroy(&patch->offsWithAll);
2948:   VecDestroy(&patch->cellMats);
2949:   VecDestroy(&patch->intFacetMats);
2950:   ISDestroy(&patch->allCells);
2951:   ISDestroy(&patch->intFacets);
2952:   ISDestroy(&patch->extFacets);
2953:   ISDestroy(&patch->intFacetsToPatchCell);
2954:   PetscSectionDestroy(&patch->intFacetCounts);
2955:   PetscSectionDestroy(&patch->extFacetCounts);

2957:   if (patch->dofSection) for (i = 0; i < patch->nsubspaces; i++) {PetscSectionDestroy(&patch->dofSection[i]);}
2958:   PetscFree(patch->dofSection);
2959:   PetscFree(patch->bs);
2960:   PetscFree(patch->nodesPerCell);
2961:   if (patch->cellNodeMap) for (i = 0; i < patch->nsubspaces; i++) {PetscFree(patch->cellNodeMap[i]);}
2962:   PetscFree(patch->cellNodeMap);
2963:   PetscFree(patch->subspaceOffsets);

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

2967:   if (patch->subspaces_to_exclude) {
2968:     PetscHSetIDestroy(&patch->subspaces_to_exclude);
2969:   }

2971:   VecDestroy(&patch->localRHS);
2972:   VecDestroy(&patch->localUpdate);
2973:   if (patch->patchRHS) {
2974:     for (i = 0; i < patch->npatch; ++i) {VecDestroy(&patch->patchRHS[i]);}
2975:     PetscFree(patch->patchRHS);
2976:   }
2977:   if (patch->patchUpdate) {
2978:     for (i = 0; i < patch->npatch; ++i) {VecDestroy(&patch->patchUpdate[i]);}
2979:     PetscFree(patch->patchUpdate);
2980:   }
2981:   VecDestroy(&patch->dof_weights);
2982:   if (patch->patch_dof_weights) {
2983:     for (i = 0; i < patch->npatch; ++i) {VecDestroy(&patch->patch_dof_weights[i]);}
2984:     PetscFree(patch->patch_dof_weights);
2985:   }
2986:   if (patch->mat) {
2987:     for (i = 0; i < patch->npatch; ++i) {MatDestroy(&patch->mat[i]);}
2988:     PetscFree(patch->mat);
2989:   }
2990:   if (patch->matWithArtificial) {
2991:     for (i = 0; i < patch->npatch; ++i) {MatDestroy(&patch->matWithArtificial[i]);}
2992:     PetscFree(patch->matWithArtificial);
2993:   }
2994:   if (patch->patchRHSWithArtificial) {
2995:     for (i = 0; i < patch->npatch; ++i) {VecDestroy(&patch->patchRHSWithArtificial[i]);}
2996:     PetscFree(patch->patchRHSWithArtificial);
2997:   }
2998:   if(patch->dofMappingWithoutToWithArtificial) {
2999:     for (i = 0; i < patch->npatch; ++i) {ISDestroy(&patch->dofMappingWithoutToWithArtificial[i]);}
3000:     PetscFree(patch->dofMappingWithoutToWithArtificial);

3002:   }
3003:   if(patch->dofMappingWithoutToWithAll) {
3004:     for (i = 0; i < patch->npatch; ++i) {ISDestroy(&patch->dofMappingWithoutToWithAll[i]);}
3005:     PetscFree(patch->dofMappingWithoutToWithAll);

3007:   }
3008:   PetscFree(patch->sub_mat_type);
3009:   if (patch->userIS) {
3010:     for (i = 0; i < patch->npatch; ++i) {ISDestroy(&patch->userIS[i]);}
3011:     PetscFree(patch->userIS);
3012:   }
3013:   PetscFree(patch->precomputedTensorLocations);
3014:   PetscFree(patch->precomputedIntFacetTensorLocations);
3015: 
3016:   patch->bs          = 0;
3017:   patch->cellNodeMap = NULL;
3018:   patch->nsubspaces  = 0;
3019:   ISDestroy(&patch->iterationSet);

3021:   PetscViewerDestroy(&patch->viewerSection);
3022:   return(0);
3023: }

3025: static PetscErrorCode PCDestroy_PATCH_Linear(PC pc)
3026: {
3027:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
3028:   PetscInt       i;

3032:   if (patch->solver) {
3033:     for (i = 0; i < patch->npatch; ++i) {KSPDestroy((KSP *) &patch->solver[i]);}
3034:     PetscFree(patch->solver);
3035:   }
3036:   return(0);
3037: }

3039: static PetscErrorCode PCDestroy_PATCH(PC pc)
3040: {
3041:   PC_PATCH      *patch = (PC_PATCH *) pc->data;

3045:   PCReset_PATCH(pc);
3046:   (*patch->destroysolver)(pc);
3047:   PetscFree(pc->data);
3048:   return(0);
3049: }

3051: static PetscErrorCode PCSetFromOptions_PATCH(PetscOptionItems *PetscOptionsObject, PC pc)
3052: {
3053:   PC_PATCH            *patch = (PC_PATCH *) pc->data;
3054:   PCPatchConstructType patchConstructionType = PC_PATCH_STAR;
3055:   char                 sub_mat_type[PETSC_MAX_PATH_LEN];
3056:   char                 option[PETSC_MAX_PATH_LEN];
3057:   const char          *prefix;
3058:   PetscBool            flg, dimflg, codimflg;
3059:   MPI_Comm             comm;
3060:   PetscInt            *ifields, nfields, k;
3061:   PetscErrorCode       ierr;
3062:   PCCompositeType loctype = PC_COMPOSITE_ADDITIVE;

3065:   PetscObjectGetComm((PetscObject) pc, &comm);
3066:   PetscObjectGetOptionsPrefix((PetscObject) pc, &prefix);
3067:   PetscOptionsHead(PetscOptionsObject, "Patch solver options");

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

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

3077:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_local_type", patch->classname);
3078:   PetscOptionsEnum(option,"Type of local solver composition (additive or multiplicative)","PCPatchSetLocalComposition",PCCompositeTypes,(PetscEnum)loctype,(PetscEnum*)&loctype,&flg);
3079:   if(flg) { PCPatchSetLocalComposition(pc, loctype);}

3081:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_dim", patch->classname);
3082:   PetscOptionsInt(option, "What dimension of mesh point to construct patches by? (0 = vertices)", "PCPATCH", patch->dim, &patch->dim, &dimflg);
3083:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_codim", patch->classname);
3084:   PetscOptionsInt(option, "What co-dimension of mesh point to construct patches by? (0 = cells)", "PCPATCH", patch->codim, &patch->codim, &codimflg);
3085:   if (dimflg && codimflg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Can only set one of dimension or co-dimension");

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

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

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

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

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

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

3107:   /* If the user has set the number of subspaces, use that for the buffer size,
3108:      otherwise use a large number */
3109:   if (patch->nsubspaces <= 0) {
3110:     nfields = 128;
3111:   } else {
3112:     nfields = patch->nsubspaces;
3113:   }
3114:   PetscMalloc1(nfields, &ifields);
3115:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exclude_subspaces", patch->classname);
3116:   PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,option,ifields,&nfields,&flg);
3117:   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");
3118:   if (flg) {
3119:     PetscHSetIClear(patch->subspaces_to_exclude);
3120:     for (k = 0; k < nfields; k++) {
3121:       PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]);
3122:     }
3123:   }
3124:   PetscFree(ifields);

3126:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_patches_view", patch->classname);
3127:   PetscOptionsBool(option, "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg);
3128:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_cells_view", patch->classname);
3129:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerCells, &patch->formatCells, &patch->viewCells);
3130:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_interior_facets_view", patch->classname);
3131:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerIntFacets, &patch->formatIntFacets, &patch->viewIntFacets);
3132:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exterior_facets_view", patch->classname);
3133:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerExtFacets, &patch->formatExtFacets, &patch->viewExtFacets);
3134:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_points_view", patch->classname);
3135:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints);
3136:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_section_view", patch->classname);
3137:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerSection, &patch->formatSection, &patch->viewSection);
3138:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_mat_view", patch->classname);
3139:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix);
3140:   PetscOptionsTail();
3141:   patch->optionsSet = PETSC_TRUE;
3142:   return(0);
3143: }

3145: static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc)
3146: {
3147:   PC_PATCH          *patch = (PC_PATCH*) pc->data;
3148:   KSPConvergedReason reason;
3149:   PetscInt           i;
3150:   PetscErrorCode     ierr;

3153:   if (!patch->save_operators) {
3154:     /* Can't do this here because the sub KSPs don't have an operator attached yet. */
3155:     return(0);
3156:   }
3157:   for (i = 0; i < patch->npatch; ++i) {
3158:     if (!((KSP) patch->solver[i])->setfromoptionscalled) {
3159:       KSPSetFromOptions((KSP) patch->solver[i]);
3160:     }
3161:     KSPSetUp((KSP) patch->solver[i]);
3162:     KSPGetConvergedReason((KSP) patch->solver[i], &reason);
3163:     if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
3164:   }
3165:   return(0);
3166: }

3168: static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer)
3169: {
3170:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
3171:   PetscViewer    sviewer;
3172:   PetscBool      isascii;
3173:   PetscMPIInt    rank;

3177:   /* TODO Redo tabbing with set tbas in new style */
3178:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
3179:   if (!isascii) return(0);
3180:   MPI_Comm_rank(PetscObjectComm((PetscObject) pc), &rank);
3181:   PetscViewerASCIIPushTab(viewer);
3182:   PetscViewerASCIIPrintf(viewer, "Subspace Correction preconditioner with %d patches\n", patch->npatch);
3183:   if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
3184:     PetscViewerASCIIPrintf(viewer, "Schwarz type: multiplicative\n");
3185:   } else {
3186:     PetscViewerASCIIPrintf(viewer, "Schwarz type: additive\n");
3187:   }
3188:   if (patch->partition_of_unity) {PetscViewerASCIIPrintf(viewer, "Weighting by partition of unity\n");}
3189:   else                           {PetscViewerASCIIPrintf(viewer, "Not weighting by partition of unity\n");}
3190:   if (patch->symmetrise_sweep) {PetscViewerASCIIPrintf(viewer, "Symmetrising sweep (start->end, then end->start)\n");}
3191:   else                         {PetscViewerASCIIPrintf(viewer, "Not symmetrising sweep\n");}
3192:   if (!patch->precomputeElementTensors) {PetscViewerASCIIPrintf(viewer, "Not precomputing element tensors (overlapping cells rebuilt in every patch assembly)\n");}
3193:   else                            {PetscViewerASCIIPrintf(viewer, "Precomputing element tensors (each cell assembled only once)\n");}
3194:   if (!patch->save_operators) {PetscViewerASCIIPrintf(viewer, "Not saving patch operators (rebuilt every PCApply)\n");}
3195:   else                        {PetscViewerASCIIPrintf(viewer, "Saving patch operators (rebuilt every PCSetUp)\n");}
3196:   if (patch->patchconstructop == PCPatchConstruct_Star)       {PetscViewerASCIIPrintf(viewer, "Patch construction operator: star\n");}
3197:   else if (patch->patchconstructop == PCPatchConstruct_Vanka) {PetscViewerASCIIPrintf(viewer, "Patch construction operator: Vanka\n");}
3198:   else if (patch->patchconstructop == PCPatchConstruct_User)  {PetscViewerASCIIPrintf(viewer, "Patch construction operator: user-specified\n");}
3199:   else                                                        {PetscViewerASCIIPrintf(viewer, "Patch construction operator: unknown\n");}

3201:   if (patch->isNonlinear) {
3202:     PetscViewerASCIIPrintf(viewer, "SNES on patches (all same):\n");
3203:   } else {
3204:     PetscViewerASCIIPrintf(viewer, "KSP on patches (all same):\n");
3205:   }
3206:   if (patch->solver) {
3207:     PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
3208:     if (!rank) {
3209:       PetscViewerASCIIPushTab(sviewer);
3210:       PetscObjectView(patch->solver[0], sviewer);
3211:       PetscViewerASCIIPopTab(sviewer);
3212:     }
3213:     PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
3214:   } else {
3215:     PetscViewerASCIIPushTab(viewer);
3216:     PetscViewerASCIIPrintf(viewer, "Solver not yet set.\n");
3217:     PetscViewerASCIIPopTab(viewer);
3218:   }
3219:   PetscViewerASCIIPopTab(viewer);
3220:   return(0);
3221: }

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

3228:   Options Database Keys:
3229: + -pc_patch_cells_view   - Views the process local cell numbers for each patch
3230: . -pc_patch_points_view  - Views the process local mesh point numbers for each patch
3231: . -pc_patch_g2l_view     - Views the map between global dofs and patch local dofs for each patch
3232: . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary
3233: - -pc_patch_sub_mat_view - Views the matrix associated with each patch

3235:   Level: intermediate

3237: .seealso: PCType, PCCreate(), PCSetType()
3238: M*/
3239: PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc)
3240: {
3241:   PC_PATCH      *patch;

3245:   PetscNewLog(pc, &patch);

3247:   if (patch->subspaces_to_exclude) {
3248:     PetscHSetIDestroy(&patch->subspaces_to_exclude);
3249:   }
3250:   PetscHSetICreate(&patch->subspaces_to_exclude);

3252:   patch->classname = "pc";
3253:   patch->isNonlinear = PETSC_FALSE;

3255:   /* Set some defaults */
3256:   patch->combined           = PETSC_FALSE;
3257:   patch->save_operators     = PETSC_TRUE;
3258:   patch->local_composition_type = PC_COMPOSITE_ADDITIVE;
3259:   patch->precomputeElementTensors = PETSC_FALSE;
3260:   patch->partition_of_unity = PETSC_FALSE;
3261:   patch->codim              = -1;
3262:   patch->dim                = -1;
3263:   patch->vankadim           = -1;
3264:   patch->ignoredim          = -1;
3265:   patch->pardecomp_overlap  = 0;
3266:   patch->patchconstructop   = PCPatchConstruct_Star;
3267:   patch->symmetrise_sweep   = PETSC_FALSE;
3268:   patch->npatch             = 0;
3269:   patch->userIS             = NULL;
3270:   patch->optionsSet         = PETSC_FALSE;
3271:   patch->iterationSet       = NULL;
3272:   patch->user_patches       = PETSC_FALSE;
3273:   PetscStrallocpy(MATDENSE, (char **) &patch->sub_mat_type);
3274:   patch->viewPatches        = PETSC_FALSE;
3275:   patch->viewCells          = PETSC_FALSE;
3276:   patch->viewPoints         = PETSC_FALSE;
3277:   patch->viewSection        = PETSC_FALSE;
3278:   patch->viewMatrix         = PETSC_FALSE;
3279:   patch->setupsolver        = PCSetUp_PATCH_Linear;
3280:   patch->applysolver        = PCApply_PATCH_Linear;
3281:   patch->resetsolver        = PCReset_PATCH_Linear;
3282:   patch->destroysolver      = PCDestroy_PATCH_Linear;
3283:   patch->updatemultiplicative = PCUpdateMultiplicative_PATCH_Linear;
3284:   patch->dofMappingWithoutToWithArtificial = NULL;
3285:   patch->dofMappingWithoutToWithAll = NULL;

3287:   pc->data                 = (void *) patch;
3288:   pc->ops->apply           = PCApply_PATCH;
3289:   pc->ops->applytranspose  = 0; /* PCApplyTranspose_PATCH; */
3290:   pc->ops->setup           = PCSetUp_PATCH;
3291:   pc->ops->reset           = PCReset_PATCH;
3292:   pc->ops->destroy         = PCDestroy_PATCH;
3293:   pc->ops->setfromoptions  = PCSetFromOptions_PATCH;
3294:   pc->ops->setuponblocks   = PCSetUpOnBlocks_PATCH;
3295:   pc->ops->view            = PCView_PATCH;
3296:   pc->ops->applyrichardson = 0;

3298:   return(0);
3299: }