Actual source code: pcpatch.c

petsc-master 2019-05-22
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:   DM              plex;
 87:   PetscBool       flg;
 88:   PetscInt        starSize;
 89:   PetscInt       *star = NULL;
 90:   PetscInt        opoint, overlapi;
 91:   PetscErrorCode  ierr;

 94:   PetscHSetIClear(ht);

 96:   DMConvert(dm, DMPLEX, &plex);
 97:   DMPlexGetChart(plex, &pStart, &pEnd);

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

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

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

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

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

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

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

152:   return(0);
153: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

509:   PCGetDM(pc, &dm);
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:   return(0);
538: }

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

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

574: /*@C

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

578:   Logically collective on PC

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

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

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

598:   Level: advanced

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

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

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

615: /*@C

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

619:   Logically collective on PC

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

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

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

639:   Level: advanced

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

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

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

656: /*@C

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

660:   Logically collective on PC

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

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

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

680:   Level: advanced

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

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

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

697: /*@C

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

701:   Logically collective on PC

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

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

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

721:   Level: advanced

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

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

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

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

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

764:     PetscHashIterGetKey(ht, hi, point);
765:     PetscHashIterNext(ht, hi);

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

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

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

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

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

856:   PetscHSetIClear(dofs);
857:   for (k = 0; k < patch->nsubspaces; ++k) {
858:     PetscInt subspaceOffset = patch->subspaceOffsets[k];
859:     PetscInt bs             = patch->bs[k];
860:     PetscInt j, l;

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

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

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

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

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

943:   /* Used to keep track of the cells in the patch. */
944:   PetscHSetICreate(&ht);
945:   PetscHSetICreate(&cht);

947:   PCGetDM(pc, &dm);
948:   if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONGSTATE, "DM not yet set on patch PC\n");
949:   DMConvert(dm, DMPLEX, &plex);
950:   DMPlexGetChart(plex, &pStart, &pEnd);
951:   DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);

953:   if (patch->user_patches) {
954:     patch->userpatchconstructionop(pc, &patch->npatch, &patch->userIS, &patch->iterationSet, patch->userpatchconstructctx);
955:     vStart = 0; vEnd = patch->npatch;
956:   } else if (patch->ctype == PC_PATCH_PARDECOMP) {
957:     vStart = 0; vEnd = 1;
958:   } else if (patch->codim < 0) {
959:     if (patch->dim < 0) {DMPlexGetDepthStratum(plex,  0,            &vStart, &vEnd);}
960:     else                {DMPlexGetDepthStratum(plex,  patch->dim,   &vStart, &vEnd);}
961:   } else                {DMPlexGetHeightStratum(plex, patch->codim, &vStart, &vEnd);}
962:   patch->npatch = vEnd - vStart;

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

972:     DMGetPointSF(dm, &sf);
973:     PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
974:     nleaves = PetscMax(nleaves, 0);
975:   }

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

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

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

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

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

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

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


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

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

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

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

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

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


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

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

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

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

1290:   /* Hash tables for artificial BC construction */
1291:   PetscHSetICreate(&ownedpts);
1292:   PetscHSetICreate(&seenpts);
1293:   PetscHSetICreate(&owneddofs);
1294:   PetscHSetICreate(&seendofs);
1295:   PetscHSetICreate(&artificialbcs);

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

1308:     PetscHMapIClear(ht);
1309:     PetscHMapIClear(htWithArtificial);
1310:     PetscHMapIClear(htWithAll);
1311:     PetscSectionGetDof(cellCounts, v, &dof);
1312:     PetscSectionGetOffset(cellCounts, v, &off);
1313:     if (dof <= 0) continue;

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

1326:       PetscHSetICreate(&globalbcdofs);
1327:       PetscSynchronizedPrintf(comm, "Patch %d: owned dofs:\n", v);
1328:       PetscHashIterBegin(owneddofs, hi);
1329:       while (!PetscHashIterAtEnd(owneddofs, hi)) {
1330:         PetscInt globalDof;

1332:         PetscHashIterGetKey(owneddofs, hi, globalDof);
1333:         PetscHashIterNext(owneddofs, hi);
1334:         PetscSynchronizedPrintf(comm, "%d ", globalDof);
1335:       }
1336:       PetscSynchronizedPrintf(comm, "\n");
1337:       PetscSynchronizedPrintf(comm, "Patch %d: seen dofs:\n", v);
1338:       PetscHashIterBegin(seendofs, hi);
1339:       while (!PetscHashIterAtEnd(seendofs, hi)) {
1340:         PetscInt globalDof;
1341:         PetscBool flg;

1343:         PetscHashIterGetKey(seendofs, hi, globalDof);
1344:         PetscHashIterNext(seendofs, hi);
1345:         PetscSynchronizedPrintf(comm, "%d ", globalDof);

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

1383:       for (i = off; i < off + dof; ++i) {
1384:         /* Walk over the cells in this patch. */
1385:         const PetscInt c    = cellsArray[i];
1386:         PetscInt       cell = c;

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

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

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

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

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

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

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

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

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

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

1529:       /* Shove it in the output data structure. */
1530:       PetscSectionGetOffset(gtolCounts, v, &goff);
1531:       PetscHashIterBegin(ht, hi);
1532:       while (!PetscHashIterAtEnd(ht, hi)) {
1533:         PetscInt globalDof, localDof;

1535:         PetscHashIterGetKey(ht, hi, globalDof);
1536:         PetscHashIterGetVal(ht, hi, localDof);
1537:         if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof;
1538:         PetscHashIterNext(ht, hi);
1539:       }

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

1564:       for (p = 0; p < Np; ++p) {
1565:         const PetscInt point = pointsArray[ooff + p];
1566:         PetscInt       globalDof, localDof;

1568:         PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof);
1569:         PetscHMapIGet(ht, globalDof, &localDof);
1570:         offsArray[(ooff + p)*Nf + k] = localDof;
1571:         if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1572:           PetscHMapIGet(htWithArtificial, globalDof, &localDof);
1573:           offsArrayWithArtificial[(ooff + p)*Nf + k] = localDof;
1574:         }
1575:         if (isNonlinear) {
1576:           PetscHMapIGet(htWithAll, globalDof, &localDof);
1577:           offsArrayWithAll[(ooff + p)*Nf + k] = localDof;
1578:         }
1579:       }
1580:     }

1582:     PetscHSetIDestroy(&globalBcs);
1583:     PetscHSetIDestroy(&ownedpts);
1584:     PetscHSetIDestroy(&seenpts);
1585:     PetscHSetIDestroy(&owneddofs);
1586:     PetscHSetIDestroy(&seendofs);
1587:     PetscHSetIDestroy(&artificialbcs);

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

1599:         if (cellNumbering) {PetscSectionGetOffset(cellNumbering, c, &cell);}
1600:         for (k = 0; k < patch->nsubspaces; ++k) {
1601:           const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
1602:           PetscInt        nodesPerCell   = patch->nodesPerCell[k];
1603:           PetscInt        subspaceOffset = patch->subspaceOffsets[k];
1604:           PetscInt        bs             = patch->bs[k];

1606:           for (j = 0; j < nodesPerCell; ++j) {
1607:             for (l = 0; l < bs; ++l) {
1608:               const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset;
1609:               PetscInt       localDof;

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

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

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

1717: static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat, PetscBool withArtificial)
1718: {
1719:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
1720:   Vec            x, y;
1721:   PetscBool      flg;
1722:   PetscInt       csize, rsize;
1723:   const char    *prefix = NULL;

1727:   if(withArtificial) {
1728:     /* would be nice if we could create a rectangular matrix of size numDofsWithArtificial x numDofs here */
1729:     x = patch->patchRHSWithArtificial[point];
1730:     y = patch->patchRHSWithArtificial[point];
1731:   } else {
1732:     x = patch->patchRHS[point];
1733:     y = patch->patchUpdate[point];
1734:   }

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

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

1772:       PetscCalloc1(rsize, &dnnz);
1773:       PetscBTCreate(rsize*rsize, &bt);
1774:       for (c = 0; c < ncell; ++c) {
1775:         const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell;
1776:         for (i = 0; i < patch->totalDofsPerCell; ++i) {
1777:           const PetscInt row = idx[i];
1778:           if (row < 0) continue;
1779:           for (j = 0; j < patch->totalDofsPerCell; ++j) {
1780:             const PetscInt col = idx[j];
1781:             const PetscInt key = row*rsize + col;
1782:             if (col < 0) continue;
1783:             if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1784:           }
1785:         }
1786:       }

1788:       if (patch->usercomputeopintfacet) {
1789:         const PetscInt *intFacetsArray = NULL;
1790:         PetscInt        i, numIntFacets, intFacetOffset;
1791:         const PetscInt *facetCells = NULL;

1793:         PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1794:         PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1795:         ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1796:         ISGetIndices(patch->intFacets, &intFacetsArray);
1797:         for (i = 0; i < numIntFacets; i++) {
1798:           const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0];
1799:           const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1];
1800:           PetscInt       celli, cellj;

1802:           for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1803:             const PetscInt row = dofsArray[(offset + cell0)*patch->totalDofsPerCell + celli];
1804:             if (row < 0) continue;
1805:             for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1806:               const PetscInt col = dofsArray[(offset + cell1)*patch->totalDofsPerCell + cellj];
1807:               const PetscInt key = row*rsize + col;
1808:               if (col < 0) continue;
1809:               if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1810:             }
1811:           }

1813:           for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1814:             const PetscInt row = dofsArray[(offset + cell1)*patch->totalDofsPerCell + celli];
1815:             if (row < 0) continue;
1816:             for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1817:               const PetscInt col = dofsArray[(offset + cell0)*patch->totalDofsPerCell + cellj];
1818:               const PetscInt key = row*rsize + col;
1819:               if (col < 0) continue;
1820:               if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1821:             }
1822:           }
1823:         }
1824:       }
1825:       PetscBTDestroy(&bt);
1826:       MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL);
1827:       PetscFree(dnnz);

1829:       PetscCalloc1(patch->totalDofsPerCell*patch->totalDofsPerCell, &zeroes);
1830:       for (c = 0; c < ncell; ++c) {
1831:         const PetscInt *idx = &dofsArray[(offset + c)*patch->totalDofsPerCell];
1832:         MatSetValues(*mat, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, zeroes, INSERT_VALUES);
1833:       }
1834:       MatGetLocalSize(*mat, &rows, NULL);
1835:       for (i = 0; i < rows; ++i) {
1836:         MatSetValues(*mat, 1, &i, 1, &i, zeroes, INSERT_VALUES);
1837:       }

1839:       if (patch->usercomputeopintfacet) {
1840:         const PetscInt *intFacetsArray = NULL;
1841:         PetscInt i, numIntFacets, intFacetOffset;
1842:         const PetscInt *facetCells = NULL;

1844:         PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1845:         PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1846:         ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1847:         ISGetIndices(patch->intFacets, &intFacetsArray);
1848:         for (i = 0; i < numIntFacets; i++) {
1849:           const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0];
1850:           const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1];
1851:           const PetscInt *cell0idx = &dofsArray[(offset + cell0)*patch->totalDofsPerCell];
1852:           const PetscInt *cell1idx = &dofsArray[(offset + cell1)*patch->totalDofsPerCell];
1853:           MatSetValues(*mat, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, zeroes, INSERT_VALUES);
1854:           MatSetValues(*mat, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, zeroes, INSERT_VALUES);
1855:         }
1856:       }

1858:       MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
1859:       MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);

1861:       PetscFree(zeroes);

1863:     } else { /* rsize too big, use MATPREALLOCATOR */
1864:       Mat preallocator;
1865:       PetscScalar* vals;

1867:       PetscCalloc1(patch->totalDofsPerCell*patch->totalDofsPerCell, &vals);
1868:       MatCreate(PETSC_COMM_SELF, &preallocator);
1869:       MatSetType(preallocator, MATPREALLOCATOR);
1870:       MatSetSizes(preallocator, rsize, rsize, rsize, rsize);
1871:       MatSetUp(preallocator);

1873:       for (c = 0; c < ncell; ++c) {
1874:         const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell;
1875:         MatSetValues(preallocator, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, vals, INSERT_VALUES);
1876:       }

1878:       if (patch->usercomputeopintfacet) {
1879:         const PetscInt *intFacetsArray = NULL;
1880:         PetscInt        i, numIntFacets, intFacetOffset;
1881:         const PetscInt *facetCells = NULL;

1883:         PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1884:         PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1885:         ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1886:         ISGetIndices(patch->intFacets, &intFacetsArray);
1887:         for (i = 0; i < numIntFacets; i++) {
1888:           const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0];
1889:           const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1];
1890:           const PetscInt *cell0idx = &dofsArray[(offset + cell0)*patch->totalDofsPerCell];
1891:           const PetscInt *cell1idx = &dofsArray[(offset + cell1)*patch->totalDofsPerCell];
1892:           MatSetValues(preallocator, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, vals, INSERT_VALUES);
1893:           MatSetValues(preallocator, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, vals, INSERT_VALUES);
1894:         }
1895:       }

1897:       PetscFree(vals);
1898:       MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1899:       MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1900:       MatPreallocatorPreallocate(preallocator, PETSC_TRUE, *mat);
1901:       MatDestroy(&preallocator);
1902:     }
1903:     PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0);
1904:     if(withArtificial) {
1905:       ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);
1906:     } else {
1907:       ISRestoreIndices(patch->dofs, &dofsArray);
1908:     }
1909:   }
1910:   MatSetUp(*mat);
1911:   return(0);
1912: }

1914: 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)
1915: {
1916:   PC_PATCH       *patch = (PC_PATCH *) pc->data;
1917:   DM              dm;
1918:   PetscSection    s;
1919:   const PetscInt *parray, *oarray;
1920:   PetscInt        Nf = patch->nsubspaces, Np, poff, p, f;
1921:   PetscErrorCode  ierr;

1924:   if (patch->precomputeElementTensors) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Precomputing element tensors not implemented with DMPlex compute operator");
1925:   PCGetDM(pc, &dm);
1926:   DMGetDefaultSection(dm, &s);
1927:   /* Set offset into patch */
1928:   PetscSectionGetDof(patch->pointCounts, patchNum, &Np);
1929:   PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);
1930:   ISGetIndices(patch->points, &parray);
1931:   ISGetIndices(patch->offs,   &oarray);
1932:   for (f = 0; f < Nf; ++f) {
1933:     for (p = 0; p < Np; ++p) {
1934:       const PetscInt point = parray[poff+p];
1935:       PetscInt       dof;

1937:       PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);
1938:       PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);
1939:       if (patch->nsubspaces == 1) {PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);}
1940:       else                        {PetscSectionSetOffset(patch->patchSection, point, -1);}
1941:     }
1942:   }
1943:   ISRestoreIndices(patch->points, &parray);
1944:   ISRestoreIndices(patch->offs,   &oarray);
1945:   if (patch->viewSection) {ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);}
1946:   DMPlexComputeResidual_Patch_Internal(pc->dm, patch->patchSection, cellIS, 0.0, x, NULL, F, ctx);
1947:   return(0);
1948: }

1950: PetscErrorCode PCPatchComputeFunction_Internal(PC pc, Vec x, Vec F, PetscInt point)
1951: {
1952:   PC_PATCH       *patch = (PC_PATCH *) pc->data;
1953:   const PetscInt *dofsArray;
1954:   const PetscInt *dofsArrayWithAll;
1955:   const PetscInt *cellsArray;
1956:   PetscInt        ncell, offset, pStart, pEnd;
1957:   PetscErrorCode  ierr;

1960:   PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);
1961:   if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n");
1962:   ISGetIndices(patch->dofs, &dofsArray);
1963:   ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);
1964:   ISGetIndices(patch->cells, &cellsArray);
1965:   PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);

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

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

1991:     PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch vector for Point %D", point);
1992:     PetscObjectSetName((PetscObject) F, name);
1993:     ObjectView((PetscObject) F, patch->viewerMatrix, patch->formatMatrix);
1994:   }
1995:   PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
1996:   return(0);
1997: }

1999: 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)
2000: {
2001:   PC_PATCH       *patch = (PC_PATCH *) pc->data;
2002:   DM              dm;
2003:   PetscSection    s;
2004:   const PetscInt *parray, *oarray;
2005:   PetscInt        Nf = patch->nsubspaces, Np, poff, p, f;
2006:   PetscErrorCode  ierr;

2009:   PCGetDM(pc, &dm);
2010:   DMGetDefaultSection(dm, &s);
2011:   /* Set offset into patch */
2012:   PetscSectionGetDof(patch->pointCounts, patchNum, &Np);
2013:   PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);
2014:   ISGetIndices(patch->points, &parray);
2015:   ISGetIndices(patch->offs,   &oarray);
2016:   for (f = 0; f < Nf; ++f) {
2017:     for (p = 0; p < Np; ++p) {
2018:       const PetscInt point = parray[poff+p];
2019:       PetscInt       dof;

2021:       PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);
2022:       PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);
2023:       if (patch->nsubspaces == 1) {PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);}
2024:       else                        {PetscSectionSetOffset(patch->patchSection, point, -1);}
2025:     }
2026:   }
2027:   ISRestoreIndices(patch->points, &parray);
2028:   ISRestoreIndices(patch->offs,   &oarray);
2029:   if (patch->viewSection) {ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);}
2030:   /* TODO Shut off MatViewFromOptions() in MatAssemblyEnd() here */
2031:   DMPlexComputeJacobian_Patch_Internal(pc->dm, patch->patchSection, patch->patchSection, cellIS, 0.0, 0.0, x, NULL, J, J, ctx);
2032:   return(0);
2033: }

2035: /* This function zeros mat on entry */
2036: PetscErrorCode PCPatchComputeOperator_Internal(PC pc, Vec x, Mat mat, PetscInt point, PetscBool withArtificial)
2037: {
2038:   PC_PATCH       *patch = (PC_PATCH *) pc->data;
2039:   const PetscInt *dofsArray;
2040:   const PetscInt *dofsArrayWithAll = NULL;
2041:   const PetscInt *cellsArray;
2042:   PetscInt        ncell, offset, pStart, pEnd, numIntFacets, intFacetOffset;
2043:   PetscBool       isNonlinear;
2044:   PetscErrorCode  ierr;

2047:   PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);
2048:   isNonlinear = patch->isNonlinear;
2049:   if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n");
2050:   if(withArtificial) {
2051:     ISGetIndices(patch->dofsWithArtificial, &dofsArray);
2052:   } else {
2053:     ISGetIndices(patch->dofs, &dofsArray);
2054:   }
2055:   if (isNonlinear) {
2056:     ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);
2057:   }
2058:   ISGetIndices(patch->cells, &cellsArray);
2059:   PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);

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

2064:   PetscSectionGetDof(patch->cellCounts, point, &ncell);
2065:   PetscSectionGetOffset(patch->cellCounts, point, &offset);
2066:   if (ncell <= 0) {
2067:     PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2068:     return(0);
2069:   }
2070:   MatZeroEntries(mat);
2071:   if (patch->precomputeElementTensors) {
2072:     PetscInt           i;
2073:     PetscInt           ndof = patch->totalDofsPerCell;
2074:     const PetscScalar *elementTensors;

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

2106:       ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
2107:       ISGetIndices(patch->intFacets, &intFacetsArray);
2108:       PCGetDM(pc, &dm);
2109:       DMPlexGetHeightStratum(dm, 1, &fStart, NULL);
2110:       /* FIXME: Pull this malloc out. */
2111:       PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofs);
2112:       if (dofsArrayWithAll) {
2113:         PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofsWithAll);
2114:       }
2115:       if (patch->precomputeElementTensors) {
2116:         PetscInt           nFacetDof = 2*patch->totalDofsPerCell;
2117:         const PetscScalar *elementTensors;

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

2121:         for (i = 0; i < numIntFacets; i++) {
2122:           const PetscInt     facet = intFacetsArray[i + intFacetOffset];
2123:           const PetscScalar *v     = elementTensors + patch->precomputedIntFacetTensorLocations[facet - fStart]*nFacetDof*nFacetDof;
2124:           idx = 0;
2125:           /*
2126:            * 0--1
2127:            * |\-|
2128:            * |+\|
2129:            * 2--3
2130:            * [0, 2, 3, 0, 1, 3]
2131:            */
2132:           for (c = 0; c < 2; c++) {
2133:             const PetscInt cell = facetCells[2*(intFacetOffset + i) + c];
2134:             for (d = 0; d < patch->totalDofsPerCell; d++) {
2135:               facetDofs[idx] = dofsArray[(offset + cell)*patch->totalDofsPerCell + d];
2136:               idx++;
2137:             }
2138:           }
2139:           MatSetValues(mat, nFacetDof, facetDofs, nFacetDof, facetDofs, v, ADD_VALUES);
2140:         }
2141:         VecRestoreArrayRead(patch->intFacetMats, &elementTensors);
2142:       } else {
2143:         /*
2144:          * 0--1
2145:          * |\-|
2146:          * |+\|
2147:          * 2--3
2148:          * [0, 2, 3, 0, 1, 3]
2149:          */
2150:         for (i = 0; i < numIntFacets; i++) {
2151:           for (c = 0; c < 2; c++) {
2152:             const PetscInt cell = facetCells[2*(intFacetOffset + i) + c];
2153:             for (d = 0; d < patch->totalDofsPerCell; d++) {
2154:               facetDofs[idx] = dofsArray[(offset + cell)*patch->totalDofsPerCell + d];
2155:               if (dofsArrayWithAll) {
2156:                 facetDofsWithAll[idx] = dofsArrayWithAll[(offset + cell)*patch->totalDofsPerCell + d];
2157:               }
2158:               idx++;
2159:             }
2160:           }
2161:         }
2162:         ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray + intFacetOffset, PETSC_USE_POINTER, &facetIS);
2163:         patch->usercomputeopintfacet(pc, point, x, mat, facetIS, 2*numIntFacets*patch->totalDofsPerCell, facetDofs, facetDofsWithAll, patch->usercomputeopintfacetctx);
2164:         ISDestroy(&facetIS);
2165:       }
2166:       ISRestoreIndices(patch->intFacetsToPatchCell, &facetCells);
2167:       ISRestoreIndices(patch->intFacets, &intFacetsArray);
2168:       PetscFree(facetDofs);
2169:       PetscFree(facetDofsWithAll);
2170:     }
2171:   }

2173:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
2174:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);

2176:   PetscStackPop;
2177:   ISDestroy(&patch->cellIS);
2178:   if(withArtificial) {
2179:     ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);
2180:   } else {
2181:     ISRestoreIndices(patch->dofs, &dofsArray);
2182:   }
2183:   if (isNonlinear) {
2184:     ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);
2185:   }
2186:   ISRestoreIndices(patch->cells, &cellsArray);
2187:   if (patch->viewMatrix) {
2188:     char name[PETSC_MAX_PATH_LEN];

2190:     PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch matrix for Point %D", point);
2191:     PetscObjectSetName((PetscObject) mat, name);
2192:     ObjectView((PetscObject) mat, patch->viewerMatrix, patch->formatMatrix);
2193:   }
2194:   PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2195:   return(0);
2196: }

2198: static PetscErrorCode MatSetValues_PCPatch_Private(Mat mat, PetscInt m, const PetscInt idxm[],
2199:                                                    PetscInt n, const PetscInt idxn[], const PetscScalar *v, InsertMode addv)
2200: {
2201:   Vec            data;
2202:   PetscScalar   *array;
2203:   PetscInt       bs, nz, i, j, cell;

2206:   MatShellGetContext(mat, &data);
2207:   VecGetBlockSize(data, &bs);
2208:   VecGetSize(data, &nz);
2209:   VecGetArray(data, &array);
2210:   if (m != n) SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Only for square insertion");
2211:   cell = (PetscInt)(idxm[0]/bs); /* use the fact that this is called once per cell */
2212:   for (i = 0; i < m; i++) {
2213:     if (idxm[i] != idxn[i]) SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Row and column indices must match!");
2214:     for (j = 0; j < n; j++) {
2215:       const PetscScalar v_ = v[i*bs + j];
2216:       /* Indexing is special to the data structure we have! */
2217:       if (addv == INSERT_VALUES) {
2218:         array[cell*bs*bs + i*bs + j] = v_;
2219:       } else {
2220:         array[cell*bs*bs + i*bs + j] += v_;
2221:       }
2222:     }
2223:   }
2224:   VecRestoreArray(data, &array);
2225:   return(0);
2226: }

2228: static PetscErrorCode PCPatchPrecomputePatchTensors_Private(PC pc)
2229: {
2230:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
2231:   const PetscInt *cellsArray;
2232:   PetscInt        ncell, offset;
2233:   const PetscInt *dofMapArray;
2234:   PetscInt        i, j;
2235:   IS              dofMap;
2236:   IS              cellIS;
2237:   const PetscInt  ndof  = patch->totalDofsPerCell;
2238:   PetscErrorCode  ierr;
2239:   Mat             vecMat;
2240:   PetscInt        cStart, cEnd;
2241:   DM              dm, plex;


2244:   ISGetSize(patch->cells, &ncell);
2245:   if (!ncell) { /* No cells to assemble over -> skip */
2246:     return(0);
2247:   }

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

2251:   if (!patch->allCells) {
2252:     PetscHSetI      cells;
2253:     PetscHashIter   hi;
2254:     PetscInt pStart, pEnd;
2255:     PetscInt *allCells = NULL;
2256:     PetscHSetICreate(&cells);
2257:     ISGetIndices(patch->cells, &cellsArray);
2258:     PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);
2259:     for (i = pStart; i < pEnd; i++) {
2260:       PetscSectionGetDof(patch->cellCounts, i, &ncell);
2261:       PetscSectionGetOffset(patch->cellCounts, i, &offset);
2262:       if (ncell <= 0) continue;
2263:       for (j = 0; j < ncell; j++) {
2264:         PetscHSetIAdd(cells, cellsArray[offset + j]);
2265:       }
2266:     }
2267:     ISRestoreIndices(patch->cells, &cellsArray);
2268:     PetscHSetIGetSize(cells, &ncell);
2269:     PetscMalloc1(ncell, &allCells);
2270:     PCGetDM(pc, &dm);
2271:     DMConvert(dm, DMPLEX, &plex);
2272:     DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2273:     PetscMalloc1(cEnd-cStart, &patch->precomputedTensorLocations);
2274:     i = 0;
2275:     PetscHashIterBegin(cells, hi);
2276:     while (!PetscHashIterAtEnd(cells, hi)) {
2277:       PetscHashIterGetKey(cells, hi, allCells[i]);
2278:       patch->precomputedTensorLocations[allCells[i]] = i;
2279:       PetscHashIterNext(cells, hi);
2280:       i++;
2281:     }
2282:     PetscHSetIDestroy(&cells);
2283:     ISCreateGeneral(PETSC_COMM_SELF, ncell, allCells, PETSC_OWN_POINTER, &patch->allCells);
2284:   }
2285:   ISGetSize(patch->allCells, &ncell);
2286:   if (!patch->cellMats) {
2287:     VecCreateSeq(PETSC_COMM_SELF, ncell*ndof*ndof, &patch->cellMats);
2288:     VecSetBlockSize(patch->cellMats, ndof);
2289:   }
2290:   VecSet(patch->cellMats, 0);

2292:   MatCreateShell(PETSC_COMM_SELF, ncell*ndof, ncell*ndof, ncell*ndof, ncell*ndof,
2293:                         (void*)patch->cellMats, &vecMat);
2294:   MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void(*)(void))&MatSetValues_PCPatch_Private);
2295:   ISGetSize(patch->allCells, &ncell);
2296:   ISCreateStride(PETSC_COMM_SELF, ndof*ncell, 0, 1, &dofMap);
2297:   ISGetIndices(dofMap, &dofMapArray);
2298:   ISGetIndices(patch->allCells, &cellsArray);
2299:   ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray, PETSC_USE_POINTER, &cellIS);
2300:   PetscStackPush("PCPatch user callback");
2301:   /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2302:   patch->usercomputeop(pc, -1, NULL, vecMat, cellIS, ndof*ncell, dofMapArray, NULL, patch->usercomputeopctx);
2303:   PetscStackPop;
2304:   ISDestroy(&cellIS);
2305:   MatDestroy(&vecMat);
2306:   ISRestoreIndices(patch->allCells, &cellsArray);
2307:   ISRestoreIndices(dofMap, &dofMapArray);
2308:   ISDestroy(&dofMap);

2310:   if (patch->usercomputeopintfacet) {
2311:     PetscInt nIntFacets;
2312:     IS       intFacetsIS;
2313:     const PetscInt *intFacetsArray = NULL;
2314:     if (!patch->allIntFacets) {
2315:       PetscHSetI      facets;
2316:       PetscHashIter   hi;
2317:       PetscInt pStart, pEnd, fStart, fEnd;
2318:       PetscInt *allIntFacets = NULL;
2319:       PetscHSetICreate(&facets);
2320:       ISGetIndices(patch->intFacets, &intFacetsArray);
2321:       PetscSectionGetChart(patch->intFacetCounts, &pStart, &pEnd);
2322:       DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2323:       for (i = pStart; i < pEnd; i++) {
2324:         PetscSectionGetDof(patch->intFacetCounts, i, &nIntFacets);
2325:         PetscSectionGetOffset(patch->intFacetCounts, i, &offset);
2326:         if (nIntFacets <= 0) continue;
2327:         for (j = 0; j < nIntFacets; j++) {
2328:           PetscHSetIAdd(facets, intFacetsArray[offset + j]);
2329:         }
2330:       }
2331:       ISRestoreIndices(patch->intFacets, &intFacetsArray);
2332:       PetscHSetIGetSize(facets, &nIntFacets);
2333:       PetscMalloc1(nIntFacets, &allIntFacets);
2334:       PCGetDM(pc, &dm);
2335:       DMConvert(dm, DMPLEX, &plex);
2336:       PetscMalloc1(fEnd-fStart, &patch->precomputedIntFacetTensorLocations);
2337:       i = 0;
2338:       PetscHashIterBegin(facets, hi);
2339:       while (!PetscHashIterAtEnd(facets, hi)) {
2340:         PetscHashIterGetKey(facets, hi, allIntFacets[i]);
2341:         patch->precomputedIntFacetTensorLocations[allIntFacets[i] - fStart] = i;
2342:         PetscHashIterNext(facets, hi);
2343:         i++;
2344:       }
2345:       PetscHSetIDestroy(&facets);
2346:       ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, allIntFacets, PETSC_OWN_POINTER, &patch->allIntFacets);
2347:     }
2348:     ISGetSize(patch->allIntFacets, &nIntFacets);
2349:     if (!patch->intFacetMats) {
2350:       VecCreateSeq(PETSC_COMM_SELF, nIntFacets*ndof*ndof*4, &patch->intFacetMats);
2351:       VecSetBlockSize(patch->intFacetMats, ndof*2);
2352:     }
2353:     VecSet(patch->intFacetMats, 0);

2355:     MatCreateShell(PETSC_COMM_SELF, nIntFacets*ndof*2, nIntFacets*ndof*2, nIntFacets*ndof*2, nIntFacets*ndof*2,
2356:                           (void*)patch->intFacetMats, &vecMat);
2357:     MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void(*)(void))&MatSetValues_PCPatch_Private);
2358:     ISCreateStride(PETSC_COMM_SELF, 2*ndof*nIntFacets, 0, 1, &dofMap);
2359:     ISGetIndices(dofMap, &dofMapArray);
2360:     ISGetIndices(patch->allIntFacets, &intFacetsArray);
2361:     ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, intFacetsArray, PETSC_USE_POINTER, &intFacetsIS);
2362:     PetscStackPush("PCPatch user callback (interior facets)");
2363:     /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2364:     patch->usercomputeopintfacet(pc, -1, NULL, vecMat, intFacetsIS, 2*ndof*nIntFacets, dofMapArray, NULL, patch->usercomputeopintfacetctx);
2365:     PetscStackPop;
2366:     ISDestroy(&intFacetsIS);
2367:     MatDestroy(&vecMat);
2368:     ISRestoreIndices(patch->allIntFacets, &intFacetsArray);
2369:     ISRestoreIndices(dofMap, &dofMapArray);
2370:     ISDestroy(&dofMap);
2371:   }
2372:   PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);

2374:   return(0);
2375: }

2377: PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat, PatchScatterType scattertype)
2378: {
2379:   PC_PATCH          *patch     = (PC_PATCH *) pc->data;
2380:   const PetscScalar *xArray    = NULL;
2381:   PetscScalar       *yArray    = NULL;
2382:   const PetscInt    *gtolArray = NULL;
2383:   PetscInt           dof, offset, lidx;
2384:   PetscErrorCode     ierr;

2387:   PetscLogEventBegin(PC_Patch_Scatter, pc, 0, 0, 0);
2388:   VecGetArrayRead(x, &xArray);
2389:   VecGetArray(y, &yArray);
2390:   if (scattertype == SCATTER_WITHARTIFICIAL) {
2391:     PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);
2392:     PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offset);
2393:     ISGetIndices(patch->gtolWithArtificial, &gtolArray);
2394:   } else if (scattertype == SCATTER_WITHALL) {
2395:     PetscSectionGetDof(patch->gtolCountsWithAll, p, &dof);
2396:     PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offset);
2397:     ISGetIndices(patch->gtolWithAll, &gtolArray);
2398:   } else {
2399:     PetscSectionGetDof(patch->gtolCounts, p, &dof);
2400:     PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2401:     ISGetIndices(patch->gtol, &gtolArray);
2402:   }
2403:   if (mode == INSERT_VALUES && scat != SCATTER_FORWARD) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't insert if not scattering forward\n");
2404:   if (mode == ADD_VALUES    && scat != SCATTER_REVERSE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't add if not scattering reverse\n");
2405:   for (lidx = 0; lidx < dof; ++lidx) {
2406:     const PetscInt gidx = gtolArray[offset+lidx];

2408:     if (mode == INSERT_VALUES) yArray[lidx]  = xArray[gidx]; /* Forward */
2409:     else                       yArray[gidx] += xArray[lidx]; /* Reverse */
2410:   }
2411:   if (scattertype == SCATTER_WITHARTIFICIAL) {
2412:     ISRestoreIndices(patch->gtolWithArtificial, &gtolArray);
2413:   } else if (scattertype == SCATTER_WITHALL) {
2414:     ISRestoreIndices(patch->gtolWithAll, &gtolArray);
2415:   } else {
2416:     ISRestoreIndices(patch->gtol, &gtolArray);
2417:   }
2418:   VecRestoreArrayRead(x, &xArray);
2419:   VecRestoreArray(y, &yArray);
2420:   PetscLogEventEnd(PC_Patch_Scatter, pc, 0, 0, 0);
2421:   return(0);
2422: }

2424: static PetscErrorCode PCSetUp_PATCH_Linear(PC pc)
2425: {
2426:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2427:   const char    *prefix;
2428:   PetscInt       i;

2432:   if (!pc->setupcalled) {
2433:     PetscMalloc1(patch->npatch, &patch->solver);
2434:     PCGetOptionsPrefix(pc, &prefix);
2435:     for (i = 0; i < patch->npatch; ++i) {
2436:       KSP ksp;
2437:       PC  subpc;

2439:       KSPCreate(PETSC_COMM_SELF, &ksp);
2440:       KSPSetErrorIfNotConverged(ksp, pc->erroriffailure);
2441:       KSPSetOptionsPrefix(ksp, prefix);
2442:       KSPAppendOptionsPrefix(ksp, "sub_");
2443:       PetscObjectIncrementTabLevel((PetscObject) ksp, (PetscObject) pc, 1);
2444:       KSPGetPC(ksp, &subpc);
2445:       PetscObjectIncrementTabLevel((PetscObject) subpc, (PetscObject) pc, 1);
2446:       PetscLogObjectParent((PetscObject) pc, (PetscObject) ksp);
2447:       patch->solver[i] = (PetscObject) ksp;
2448:     }
2449:   }
2450:   if (patch->save_operators) {
2451:     if (patch->precomputeElementTensors) {
2452:       PCPatchPrecomputePatchTensors_Private(pc);
2453:     }
2454:     for (i = 0; i < patch->npatch; ++i) {
2455:       PCPatchComputeOperator_Internal(pc, NULL, patch->mat[i], i, PETSC_FALSE);
2456:       KSPSetOperators((KSP) patch->solver[i], patch->mat[i], patch->mat[i]);
2457:     }
2458:   }
2459:   if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2460:     for (i = 0; i < patch->npatch; ++i) {
2461:       /* Instead of padding patch->patchUpdate with zeros to get */
2462:       /* patch->patchUpdateWithArtificial and then multiplying with the matrix, */
2463:       /* just get rid of the columns that correspond to the dofs with */
2464:       /* artificial bcs. That's of course fairly inefficient, hopefully we */
2465:       /* can just assemble the rectangular matrix in the first place. */
2466:       Mat matSquare;
2467:       IS rowis;
2468:       PetscInt dof;

2470:       MatGetSize(patch->mat[i], &dof, NULL);
2471:       if (dof == 0) {
2472:         patch->matWithArtificial[i] = NULL;
2473:         continue;
2474:       }

2476:       PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);
2477:       PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);

2479:       MatGetSize(matSquare, &dof, NULL);
2480:       ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis);
2481:       if(pc->setupcalled) {
2482:         MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_REUSE_MATRIX, &patch->matWithArtificial[i]);
2483:       } else {
2484:         MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &patch->matWithArtificial[i]);
2485:       }
2486:       ISDestroy(&rowis);
2487:       MatDestroy(&matSquare);
2488:     }
2489:   }
2490:   return(0);
2491: }

2493: static PetscErrorCode PCSetUp_PATCH(PC pc)
2494: {
2495:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2496:   PetscInt       i;
2497:   PetscBool       isNonlinear;

2501:   if (!pc->setupcalled) {
2502:     PetscInt pStart, pEnd, p;
2503:     PetscInt localSize;

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

2507:     isNonlinear = patch->isNonlinear;
2508:     if (!patch->nsubspaces) {
2509:       DM           dm;
2510:       PetscSection s;
2511:       PetscInt     cStart, cEnd, c, Nf, f, numGlobalBcs = 0, *globalBcs, *Nb, totNb = 0, **cellDofs;

2513:       PCGetDM(pc, &dm);
2514:       if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "Must set DM for PCPATCH or call PCPatchSetDiscretisationInfo()");
2515:       DMGetDefaultSection(dm, &s);
2516:       PetscSectionGetNumFields(s, &Nf);
2517:       PetscSectionGetChart(s, &pStart, &pEnd);
2518:       for (p = pStart; p < pEnd; ++p) {
2519:         PetscInt cdof;
2520:         PetscSectionGetConstraintDof(s, p, &cdof);
2521:         numGlobalBcs += cdof;
2522:       }
2523:       DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2524:       PetscMalloc3(Nf, &Nb, Nf, &cellDofs, numGlobalBcs, &globalBcs);
2525:       for (f = 0; f < Nf; ++f) {
2526:         PetscFE        fe;
2527:         PetscDualSpace sp;
2528:         PetscInt       cdoff = 0;

2530:         DMGetField(dm, f, NULL, (PetscObject *) &fe);
2531:         /* PetscFEGetNumComponents(fe, &Nc[f]); */
2532:         PetscFEGetDualSpace(fe, &sp);
2533:         PetscDualSpaceGetDimension(sp, &Nb[f]);
2534:         totNb += Nb[f];

2536:         PetscMalloc1((cEnd-cStart)*Nb[f], &cellDofs[f]);
2537:         for (c = cStart; c < cEnd; ++c) {
2538:           PetscInt *closure = NULL;
2539:           PetscInt  clSize  = 0, cl;

2541:           DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
2542:           for (cl = 0; cl < clSize*2; cl += 2) {
2543:             const PetscInt p = closure[cl];
2544:             PetscInt       fdof, d, foff;

2546:             PetscSectionGetFieldDof(s, p, f, &fdof);
2547:             PetscSectionGetFieldOffset(s, p, f, &foff);
2548:             for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d;
2549:           }
2550:           DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
2551:         }
2552:         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]);
2553:       }
2554:       numGlobalBcs = 0;
2555:       for (p = pStart; p < pEnd; ++p) {
2556:         const PetscInt *ind;
2557:         PetscInt        off, cdof, d;

2559:         PetscSectionGetOffset(s, p, &off);
2560:         PetscSectionGetConstraintDof(s, p, &cdof);
2561:         PetscSectionGetConstraintIndices(s, p, &ind);
2562:         for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d];
2563:       }

2565:       PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **) cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs);
2566:       for (f = 0; f < Nf; ++f) {
2567:         PetscFree(cellDofs[f]);
2568:       }
2569:       PetscFree3(Nb, cellDofs, globalBcs);
2570:       PCPatchSetComputeFunction(pc, PCPatchComputeFunction_DMPlex_Private, NULL);
2571:       PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL);
2572:     }

2574:     localSize = patch->subspaceOffsets[patch->nsubspaces];
2575:     VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localRHS);
2576:     VecSetUp(patch->localRHS);
2577:     VecDuplicate(patch->localRHS, &patch->localUpdate);
2578:     PCPatchCreateCellPatches(pc);
2579:     PCPatchCreateCellPatchDiscretisationInfo(pc);

2581:     /* OK, now build the work vectors */
2582:     PetscSectionGetChart(patch->gtolCounts, &pStart, &pEnd);
2583:     PetscMalloc1(patch->npatch, &patch->patchRHS);
2584:     PetscMalloc1(patch->npatch, &patch->patchUpdate);

2586:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2587:       PetscMalloc1(patch->npatch, &patch->patchRHSWithArtificial);
2588:       PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithArtificial);
2589:     }
2590:     if (isNonlinear) {
2591:       PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithAll);
2592:     }
2593:     for (p = pStart; p < pEnd; ++p) {
2594:       PetscInt dof;

2596:       PetscSectionGetDof(patch->gtolCounts, p, &dof);
2597:       VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchRHS[p-pStart]);
2598:       VecSetUp(patch->patchRHS[p-pStart]);
2599:       VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchUpdate[p-pStart]);
2600:       VecSetUp(patch->patchUpdate[p-pStart]);
2601:       if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2602:         const PetscInt    *gtolArray, *gtolArrayWithArtificial = NULL;
2603:         PetscInt           numPatchDofs, offset;
2604:         PetscInt           numPatchDofsWithArtificial, offsetWithArtificial;
2605:         PetscInt           dofWithoutArtificialCounter = 0;
2606:         PetscInt          *patchWithoutArtificialToWithArtificialArray;

2608:         PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);
2609:         VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchRHSWithArtificial[p-pStart]);
2610:         VecSetUp(patch->patchRHSWithArtificial[p-pStart]);

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

2616:         PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);
2617:         if (numPatchDofs == 0) {
2618:           patch->dofMappingWithoutToWithArtificial[p-pStart] = NULL;
2619:           continue;
2620:         }

2622:         PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2623:         ISGetIndices(patch->gtolWithArtificial, &gtolArrayWithArtificial);
2624:         PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &numPatchDofsWithArtificial);
2625:         PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offsetWithArtificial);

2627:         PetscMalloc1(numPatchDofs, &patchWithoutArtificialToWithArtificialArray);
2628:         for (i=0; i<numPatchDofsWithArtificial; i++) {
2629:           if (gtolArrayWithArtificial[i+offsetWithArtificial] == gtolArray[offset+dofWithoutArtificialCounter]) {
2630:             patchWithoutArtificialToWithArtificialArray[dofWithoutArtificialCounter] = i;
2631:             dofWithoutArtificialCounter++;
2632:             if (dofWithoutArtificialCounter == numPatchDofs)
2633:               break;
2634:           }
2635:         }
2636:         ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutArtificialToWithArtificialArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithArtificial[p-pStart]);
2637:         ISRestoreIndices(patch->gtol, &gtolArray);
2638:         ISRestoreIndices(patch->gtolWithArtificial, &gtolArrayWithArtificial);
2639:       }
2640:       if (isNonlinear) {
2641:         const PetscInt    *gtolArray, *gtolArrayWithAll = NULL;
2642:         PetscInt           numPatchDofs, offset;
2643:         PetscInt           numPatchDofsWithAll, offsetWithAll;
2644:         PetscInt           dofWithoutAllCounter = 0;
2645:         PetscInt          *patchWithoutAllToWithAllArray;

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

2651:         PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);
2652:         if (numPatchDofs == 0) {
2653:           patch->dofMappingWithoutToWithAll[p-pStart] = NULL;
2654:           continue;
2655:         }

2657:         PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2658:         ISGetIndices(patch->gtolWithAll, &gtolArrayWithAll);
2659:         PetscSectionGetDof(patch->gtolCountsWithAll, p, &numPatchDofsWithAll);
2660:         PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offsetWithAll);

2662:         PetscMalloc1(numPatchDofs, &patchWithoutAllToWithAllArray);

2664:         for (i=0; i<numPatchDofsWithAll; i++) {
2665:           if (gtolArrayWithAll[i+offsetWithAll] == gtolArray[offset+dofWithoutAllCounter]) {
2666:             patchWithoutAllToWithAllArray[dofWithoutAllCounter] = i;
2667:             dofWithoutAllCounter++;
2668:             if (dofWithoutAllCounter == numPatchDofs)
2669:               break;
2670:           }
2671:         }
2672:         ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutAllToWithAllArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithAll[p-pStart]);
2673:         ISRestoreIndices(patch->gtol, &gtolArray);
2674:         ISRestoreIndices(patch->gtolWithAll, &gtolArrayWithAll);
2675:       }
2676:     }
2677:     if (patch->save_operators) {
2678:       PetscMalloc1(patch->npatch, &patch->mat);
2679:       for (i = 0; i < patch->npatch; ++i) {
2680:         PCPatchCreateMatrix_Private(pc, i, &patch->mat[i], PETSC_FALSE);
2681:       }
2682:     }
2683:     PetscLogEventEnd(PC_Patch_CreatePatches, pc, 0, 0, 0);

2685:     /* If desired, calculate weights for dof multiplicity */
2686:     if (patch->partition_of_unity) {
2687:       PetscScalar *input = NULL;
2688:       PetscScalar *output = NULL;
2689:       Vec global;

2691:       VecDuplicate(patch->localRHS, &patch->dof_weights);
2692:       if(patch->local_composition_type == PC_COMPOSITE_ADDITIVE) {
2693:         for (i = 0; i < patch->npatch; ++i) {
2694:           PetscInt dof;

2696:           PetscSectionGetDof(patch->gtolCounts, i+pStart, &dof);
2697:           if (dof <= 0) continue;
2698:           VecSet(patch->patchRHS[i], 1.0);
2699:           PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchRHS[i], patch->dof_weights, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);
2700:         }
2701:       } else {
2702:         /* multiplicative is actually only locally multiplicative and globally additive. need the pou where the mesh decomposition overlaps */
2703:         VecSet(patch->dof_weights, 1.0);
2704:       }

2706:       VecDuplicate(patch->dof_weights, &global);
2707:       VecSet(global, 0.);

2709:       VecGetArray(patch->dof_weights, &input);
2710:       VecGetArray(global, &output);
2711:       PetscSFReduceBegin(patch->defaultSF, MPIU_SCALAR, input, output, MPI_SUM);
2712:       PetscSFReduceEnd(patch->defaultSF, MPIU_SCALAR, input, output, MPI_SUM);
2713:       VecRestoreArray(patch->dof_weights, &input);
2714:       VecRestoreArray(global, &output);

2716:       VecReciprocal(global);

2718:       VecGetArray(patch->dof_weights, &output);
2719:       VecGetArray(global, &input);
2720:       PetscSFBcastBegin(patch->defaultSF, MPIU_SCALAR, input, output);
2721:       PetscSFBcastEnd(patch->defaultSF, MPIU_SCALAR, input, output);
2722:       VecRestoreArray(patch->dof_weights, &output);
2723:       VecRestoreArray(global, &input);
2724:       VecDestroy(&global);
2725:     }
2726:     if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE && patch->save_operators) {
2727:       PetscMalloc1(patch->npatch, &patch->matWithArtificial);
2728:     }
2729:   }
2730:   (*patch->setupsolver)(pc);
2731:   return(0);
2732: }

2734: static PetscErrorCode PCApply_PATCH_Linear(PC pc, PetscInt i, Vec x, Vec y)
2735: {
2736:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2737:   KSP            ksp   = (KSP) patch->solver[i];

2741:   if (!patch->save_operators) {
2742:     Mat mat;

2744:     PCPatchCreateMatrix_Private(pc, i, &mat, PETSC_FALSE);
2745:     /* Populate operator here. */
2746:     PCPatchComputeOperator_Internal(pc, NULL, mat, i, PETSC_FALSE);
2747:     KSPSetOperators(ksp, mat, mat);
2748:     /* Drop reference so the KSPSetOperators below will blow it away. */
2749:     MatDestroy(&mat);
2750:   }
2751:   PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);
2752:   if (!ksp->setfromoptionscalled) {
2753:     KSPSetFromOptions(ksp);
2754:   }
2755:   KSPSolve(ksp, x, y);
2756:   KSPCheckSolve(ksp, pc, y);
2757:   PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);
2758:   if (!patch->save_operators) {
2759:     PC pc;
2760:     KSPSetOperators(ksp, NULL, NULL);
2761:     KSPGetPC(ksp, &pc);
2762:     /* Destroy PC context too, otherwise the factored matrix hangs around. */
2763:     PCReset(pc);
2764:   }
2765:   return(0);
2766: }

2768: static PetscErrorCode PCUpdateMultiplicative_PATCH_Linear(PC pc, PetscInt i, PetscInt pStart)
2769: {
2770:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2771:   Mat multMat;


2776:   if (patch->save_operators) {
2777:     multMat = patch->matWithArtificial[i];
2778:   } else {
2779:     /*Very inefficient, hopefully we can just assemble the rectangular matrix in the first place.*/
2780:     Mat matSquare;
2781:     PetscInt dof;
2782:     IS rowis;
2783:     PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);
2784:     PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);
2785:     MatGetSize(matSquare, &dof, NULL);
2786:     ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis);
2787:     MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &multMat);
2788:     MatDestroy(&matSquare);
2789:     ISDestroy(&rowis);
2790:   }
2791:   MatMult(multMat, patch->patchUpdate[i], patch->patchRHSWithArtificial[i]);
2792:   VecScale(patch->patchRHSWithArtificial[i], -1.0);
2793:   PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHSWithArtificial[i], patch->localRHS, ADD_VALUES, SCATTER_REVERSE, SCATTER_WITHARTIFICIAL);
2794:   if (!patch->save_operators) {
2795:     MatDestroy(&multMat);
2796:   }
2797:   return(0);
2798: }

2800: static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y)
2801: {
2802:   PC_PATCH          *patch    = (PC_PATCH *) pc->data;
2803:   const PetscScalar *globalRHS  = NULL;
2804:   PetscScalar       *localRHS   = NULL;
2805:   PetscScalar       *globalUpdate  = NULL;
2806:   const PetscInt    *bcNodes  = NULL;
2807:   PetscInt           nsweep   = patch->symmetrise_sweep ? 2 : 1;
2808:   PetscInt           start[2] = {0, 0};
2809:   PetscInt           end[2]   = {-1, -1};
2810:   const PetscInt     inc[2]   = {1, -1};
2811:   const PetscScalar *localUpdate;
2812:   const PetscInt    *iterationSet;
2813:   PetscInt           pStart, numBcs, n, sweep, bc, j;
2814:   PetscErrorCode     ierr;

2817:   PetscLogEventBegin(PC_Patch_Apply, pc, 0, 0, 0);
2818:   PetscOptionsPushGetViewerOff(PETSC_TRUE);
2819:   /* start, end, inc have 2 entries to manage a second backward sweep if we symmetrize */
2820:   end[0]   = patch->npatch;
2821:   start[1] = patch->npatch-1;
2822:   if (patch->user_patches) {
2823:     ISGetLocalSize(patch->iterationSet, &end[0]);
2824:     start[1] = end[0] - 1;
2825:     ISGetIndices(patch->iterationSet, &iterationSet);
2826:   }
2827:   /* Scatter from global space into overlapped local spaces */
2828:   VecGetArrayRead(x, &globalRHS);
2829:   VecGetArray(patch->localRHS, &localRHS);
2830:   PetscSFBcastBegin(patch->defaultSF, MPIU_SCALAR, globalRHS, localRHS);
2831:   PetscSFBcastEnd(patch->defaultSF, MPIU_SCALAR, globalRHS, localRHS);
2832:   VecRestoreArrayRead(x, &globalRHS);
2833:   VecRestoreArray(patch->localRHS, &localRHS);

2835:   VecSet(patch->localUpdate, 0.0);
2836:   PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);
2837:   for (sweep = 0; sweep < nsweep; sweep++) {
2838:     for (j = start[sweep]; j*inc[sweep] < end[sweep]*inc[sweep]; j += inc[sweep]) {
2839:       PetscInt i       = patch->user_patches ? iterationSet[j] : j;
2840:       PetscInt start, len;

2842:       PetscSectionGetDof(patch->gtolCounts, i+pStart, &len);
2843:       PetscSectionGetOffset(patch->gtolCounts, i+pStart, &start);
2844:       /* TODO: Squash out these guys in the setup as well. */
2845:       if (len <= 0) continue;
2846:       /* TODO: Do we need different scatters for X and Y? */
2847:       PCPatch_ScatterLocal_Private(pc, i+pStart, patch->localRHS, patch->patchRHS[i], INSERT_VALUES, SCATTER_FORWARD, SCATTER_INTERIOR);
2848:       (*patch->applysolver)(pc, i, patch->patchRHS[i], patch->patchUpdate[i]);
2849:       PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchUpdate[i], patch->localUpdate, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);
2850:       if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2851:         (*patch->updatemultiplicative)(pc, i, pStart);
2852:       }
2853:     }
2854:   }
2855:   if (patch->user_patches) {ISRestoreIndices(patch->iterationSet, &iterationSet);}
2856:   /* XXX: should we do this on the global vector? */
2857:   if (patch->partition_of_unity) {
2858:     VecPointwiseMult(patch->localUpdate, patch->localUpdate, patch->dof_weights);
2859:   }
2860:   /* Now patch->localUpdate contains the solution of the patch solves, so we need to combine them all. */
2861:   VecSet(y, 0.0);
2862:   VecGetArray(y, &globalUpdate);
2863:   VecGetArrayRead(patch->localUpdate, &localUpdate);
2864:   PetscSFReduceBegin(patch->defaultSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);
2865:   PetscSFReduceEnd(patch->defaultSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);
2866:   VecRestoreArrayRead(patch->localUpdate, &localUpdate);

2868:   /* Now we need to send the global BC values through */
2869:   VecGetArrayRead(x, &globalRHS);
2870:   ISGetSize(patch->globalBcNodes, &numBcs);
2871:   ISGetIndices(patch->globalBcNodes, &bcNodes);
2872:   VecGetLocalSize(x, &n);
2873:   for (bc = 0; bc < numBcs; ++bc) {
2874:     const PetscInt idx = bcNodes[bc];
2875:     if (idx < n) globalUpdate[idx] = globalRHS[idx];
2876:   }

2878:   ISRestoreIndices(patch->globalBcNodes, &bcNodes);
2879:   VecRestoreArrayRead(x, &globalRHS);
2880:   VecRestoreArray(y, &globalUpdate);

2882:   PetscOptionsPopGetViewerOff();
2883:   PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0);
2884:   return(0);
2885: }

2887: static PetscErrorCode PCReset_PATCH_Linear(PC pc)
2888: {
2889:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2890:   PetscInt       i;

2894:   if (patch->solver) {
2895:     for (i = 0; i < patch->npatch; ++i) {KSPReset((KSP) patch->solver[i]);}
2896:   }
2897:   return(0);
2898: }

2900: static PetscErrorCode PCReset_PATCH(PC pc)
2901: {
2902:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2903:   PetscInt       i;


2908:   PetscSFDestroy(&patch->defaultSF);
2909:   PetscSectionDestroy(&patch->cellCounts);
2910:   PetscSectionDestroy(&patch->pointCounts);
2911:   PetscSectionDestroy(&patch->cellNumbering);
2912:   PetscSectionDestroy(&patch->gtolCounts);
2913:   ISDestroy(&patch->gtol);
2914:   ISDestroy(&patch->cells);
2915:   ISDestroy(&patch->points);
2916:   ISDestroy(&patch->dofs);
2917:   ISDestroy(&patch->offs);
2918:   PetscSectionDestroy(&patch->patchSection);
2919:   ISDestroy(&patch->ghostBcNodes);
2920:   ISDestroy(&patch->globalBcNodes);
2921:   PetscSectionDestroy(&patch->gtolCountsWithArtificial);
2922:   ISDestroy(&patch->gtolWithArtificial);
2923:   ISDestroy(&patch->dofsWithArtificial);
2924:   ISDestroy(&patch->offsWithArtificial);
2925:   PetscSectionDestroy(&patch->gtolCountsWithAll);
2926:   ISDestroy(&patch->gtolWithAll);
2927:   ISDestroy(&patch->dofsWithAll);
2928:   ISDestroy(&patch->offsWithAll);
2929:   VecDestroy(&patch->cellMats);
2930:   VecDestroy(&patch->intFacetMats);
2931:   ISDestroy(&patch->allCells);
2932:   ISDestroy(&patch->intFacets);
2933:   ISDestroy(&patch->extFacets);
2934:   ISDestroy(&patch->intFacetsToPatchCell);
2935:   PetscSectionDestroy(&patch->intFacetCounts);
2936:   PetscSectionDestroy(&patch->extFacetCounts);

2938:   if (patch->dofSection) for (i = 0; i < patch->nsubspaces; i++) {PetscSectionDestroy(&patch->dofSection[i]);}
2939:   PetscFree(patch->dofSection);
2940:   PetscFree(patch->bs);
2941:   PetscFree(patch->nodesPerCell);
2942:   if (patch->cellNodeMap) for (i = 0; i < patch->nsubspaces; i++) {PetscFree(patch->cellNodeMap[i]);}
2943:   PetscFree(patch->cellNodeMap);
2944:   PetscFree(patch->subspaceOffsets);

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

2948:   if (patch->subspaces_to_exclude) {
2949:     PetscHSetIDestroy(&patch->subspaces_to_exclude);
2950:   }

2952:   VecDestroy(&patch->localRHS);
2953:   VecDestroy(&patch->localUpdate);
2954:   if (patch->patchRHS) {
2955:     for (i = 0; i < patch->npatch; ++i) {VecDestroy(&patch->patchRHS[i]);}
2956:     PetscFree(patch->patchRHS);
2957:   }
2958:   if (patch->patchUpdate) {
2959:     for (i = 0; i < patch->npatch; ++i) {VecDestroy(&patch->patchUpdate[i]);}
2960:     PetscFree(patch->patchUpdate);
2961:   }
2962:   VecDestroy(&patch->dof_weights);
2963:   if (patch->patch_dof_weights) {
2964:     for (i = 0; i < patch->npatch; ++i) {VecDestroy(&patch->patch_dof_weights[i]);}
2965:     PetscFree(patch->patch_dof_weights);
2966:   }
2967:   if (patch->mat) {
2968:     for (i = 0; i < patch->npatch; ++i) {MatDestroy(&patch->mat[i]);}
2969:     PetscFree(patch->mat);
2970:   }
2971:   if (patch->matWithArtificial) {
2972:     for (i = 0; i < patch->npatch; ++i) {MatDestroy(&patch->matWithArtificial[i]);}
2973:     PetscFree(patch->matWithArtificial);
2974:   }
2975:   if (patch->patchRHSWithArtificial) {
2976:     for (i = 0; i < patch->npatch; ++i) {VecDestroy(&patch->patchRHSWithArtificial[i]);}
2977:     PetscFree(patch->patchRHSWithArtificial);
2978:   }
2979:   if(patch->dofMappingWithoutToWithArtificial) {
2980:     for (i = 0; i < patch->npatch; ++i) {ISDestroy(&patch->dofMappingWithoutToWithArtificial[i]);}
2981:     PetscFree(patch->dofMappingWithoutToWithArtificial);

2983:   }
2984:   if(patch->dofMappingWithoutToWithAll) {
2985:     for (i = 0; i < patch->npatch; ++i) {ISDestroy(&patch->dofMappingWithoutToWithAll[i]);}
2986:     PetscFree(patch->dofMappingWithoutToWithAll);

2988:   }
2989:   PetscFree(patch->sub_mat_type);
2990:   if (patch->userIS) {
2991:     for (i = 0; i < patch->npatch; ++i) {ISDestroy(&patch->userIS[i]);}
2992:     PetscFree(patch->userIS);
2993:   }
2994:   PetscFree(patch->precomputedTensorLocations);
2995:   PetscFree(patch->precomputedIntFacetTensorLocations);
2996: 
2997:   patch->bs          = 0;
2998:   patch->cellNodeMap = NULL;
2999:   patch->nsubspaces  = 0;
3000:   ISDestroy(&patch->iterationSet);

3002:   PetscViewerDestroy(&patch->viewerSection);
3003:   return(0);
3004: }

3006: static PetscErrorCode PCDestroy_PATCH_Linear(PC pc)
3007: {
3008:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
3009:   PetscInt       i;

3013:   if (patch->solver) {
3014:     for (i = 0; i < patch->npatch; ++i) {KSPDestroy((KSP *) &patch->solver[i]);}
3015:     PetscFree(patch->solver);
3016:   }
3017:   return(0);
3018: }

3020: static PetscErrorCode PCDestroy_PATCH(PC pc)
3021: {
3022:   PC_PATCH      *patch = (PC_PATCH *) pc->data;

3026:   PCReset_PATCH(pc);
3027:   (*patch->destroysolver)(pc);
3028:   PetscFree(pc->data);
3029:   return(0);
3030: }

3032: static PetscErrorCode PCSetFromOptions_PATCH(PetscOptionItems *PetscOptionsObject, PC pc)
3033: {
3034:   PC_PATCH            *patch = (PC_PATCH *) pc->data;
3035:   PCPatchConstructType patchConstructionType = PC_PATCH_STAR;
3036:   char                 sub_mat_type[PETSC_MAX_PATH_LEN];
3037:   char                 option[PETSC_MAX_PATH_LEN];
3038:   const char          *prefix;
3039:   PetscBool            flg, dimflg, codimflg;
3040:   MPI_Comm             comm;
3041:   PetscInt            *ifields, nfields, k;
3042:   PetscErrorCode       ierr;
3043:   PCCompositeType loctype = PC_COMPOSITE_ADDITIVE;

3046:   PetscObjectGetComm((PetscObject) pc, &comm);
3047:   PetscObjectGetOptionsPrefix((PetscObject) pc, &prefix);
3048:   PetscOptionsHead(PetscOptionsObject, "Patch solver options");

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

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

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

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

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

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

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

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

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

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

3088:   /* If the user has set the number of subspaces, use that for the buffer size,
3089:      otherwise use a large number */
3090:   if (patch->nsubspaces <= 0) {
3091:     nfields = 128;
3092:   } else {
3093:     nfields = patch->nsubspaces;
3094:   }
3095:   PetscMalloc1(nfields, &ifields);
3096:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exclude_subspaces", patch->classname);
3097:   PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,option,ifields,&nfields,&flg);
3098:   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");
3099:   if (flg) {
3100:     PetscHSetIClear(patch->subspaces_to_exclude);
3101:     for (k = 0; k < nfields; k++) {
3102:       PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]);
3103:     }
3104:   }
3105:   PetscFree(ifields);

3107:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_patches_view", patch->classname);
3108:   PetscOptionsBool(option, "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg);
3109:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_cells_view", patch->classname);
3110:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerCells, &patch->formatCells, &patch->viewCells);
3111:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_interior_facets_view", patch->classname);
3112:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerIntFacets, &patch->formatIntFacets, &patch->viewIntFacets);
3113:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exterior_facets_view", patch->classname);
3114:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerExtFacets, &patch->formatExtFacets, &patch->viewExtFacets);
3115:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_points_view", patch->classname);
3116:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints);
3117:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_section_view", patch->classname);
3118:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerSection, &patch->formatSection, &patch->viewSection);
3119:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_mat_view", patch->classname);
3120:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix);
3121:   PetscOptionsTail();
3122:   patch->optionsSet = PETSC_TRUE;
3123:   return(0);
3124: }

3126: static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc)
3127: {
3128:   PC_PATCH          *patch = (PC_PATCH*) pc->data;
3129:   KSPConvergedReason reason;
3130:   PetscInt           i;
3131:   PetscErrorCode     ierr;

3134:   if (!patch->save_operators) {
3135:     /* Can't do this here because the sub KSPs don't have an operator attached yet. */
3136:     return(0);
3137:   }
3138:   for (i = 0; i < patch->npatch; ++i) {
3139:     if (!((KSP) patch->solver[i])->setfromoptionscalled) {
3140:       KSPSetFromOptions((KSP) patch->solver[i]);
3141:     }
3142:     KSPSetUp((KSP) patch->solver[i]);
3143:     KSPGetConvergedReason((KSP) patch->solver[i], &reason);
3144:     if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
3145:   }
3146:   return(0);
3147: }

3149: static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer)
3150: {
3151:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
3152:   PetscViewer    sviewer;
3153:   PetscBool      isascii;
3154:   PetscMPIInt    rank;

3158:   /* TODO Redo tabbing with set tbas in new style */
3159:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
3160:   if (!isascii) return(0);
3161:   MPI_Comm_rank(PetscObjectComm((PetscObject) pc), &rank);
3162:   PetscViewerASCIIPushTab(viewer);
3163:   PetscViewerASCIIPrintf(viewer, "Subspace Correction preconditioner with %d patches\n", patch->npatch);
3164:   if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
3165:     PetscViewerASCIIPrintf(viewer, "Schwarz type: multiplicative\n");
3166:   } else {
3167:     PetscViewerASCIIPrintf(viewer, "Schwarz type: additive\n");
3168:   }
3169:   if (patch->partition_of_unity) {PetscViewerASCIIPrintf(viewer, "Weighting by partition of unity\n");}
3170:   else                           {PetscViewerASCIIPrintf(viewer, "Not weighting by partition of unity\n");}
3171:   if (patch->symmetrise_sweep) {PetscViewerASCIIPrintf(viewer, "Symmetrising sweep (start->end, then end->start)\n");}
3172:   else                         {PetscViewerASCIIPrintf(viewer, "Not symmetrising sweep\n");}
3173:   if (!patch->precomputeElementTensors) {PetscViewerASCIIPrintf(viewer, "Not precomputing element tensors (overlapping cells rebuilt in every patch assembly)\n");}
3174:   else                            {PetscViewerASCIIPrintf(viewer, "Precomputing element tensors (each cell assembled only once)\n");}
3175:   if (!patch->save_operators) {PetscViewerASCIIPrintf(viewer, "Not saving patch operators (rebuilt every PCApply)\n");}
3176:   else                        {PetscViewerASCIIPrintf(viewer, "Saving patch operators (rebuilt every PCSetUp)\n");}
3177:   if (patch->patchconstructop == PCPatchConstruct_Star)       {PetscViewerASCIIPrintf(viewer, "Patch construction operator: star\n");}
3178:   else if (patch->patchconstructop == PCPatchConstruct_Vanka) {PetscViewerASCIIPrintf(viewer, "Patch construction operator: Vanka\n");}
3179:   else if (patch->patchconstructop == PCPatchConstruct_User)  {PetscViewerASCIIPrintf(viewer, "Patch construction operator: user-specified\n");}
3180:   else                                                        {PetscViewerASCIIPrintf(viewer, "Patch construction operator: unknown\n");}

3182:   if (patch->isNonlinear) {
3183:     PetscViewerASCIIPrintf(viewer, "SNES on patches (all same):\n");
3184:   } else {
3185:     PetscViewerASCIIPrintf(viewer, "KSP on patches (all same):\n");
3186:   }
3187:   if (patch->solver) {
3188:     PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
3189:     if (!rank) {
3190:       PetscViewerASCIIPushTab(sviewer);
3191:       PetscObjectView(patch->solver[0], sviewer);
3192:       PetscViewerASCIIPopTab(sviewer);
3193:     }
3194:     PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
3195:   } else {
3196:     PetscViewerASCIIPushTab(viewer);
3197:     PetscViewerASCIIPrintf(viewer, "Solver not yet set.\n");
3198:     PetscViewerASCIIPopTab(viewer);
3199:   }
3200:   PetscViewerASCIIPopTab(viewer);
3201:   return(0);
3202: }

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

3209:   Options Database Keys:
3210: + -pc_patch_cells_view   - Views the process local cell numbers for each patch
3211: . -pc_patch_points_view  - Views the process local mesh point numbers for each patch
3212: . -pc_patch_g2l_view     - Views the map between global dofs and patch local dofs for each patch
3213: . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary
3214: - -pc_patch_sub_mat_view - Views the matrix associated with each patch

3216:   Level: intermediate

3218: .seealso: PCType, PCCreate(), PCSetType()
3219: M*/
3220: PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc)
3221: {
3222:   PC_PATCH      *patch;

3226:   PetscNewLog(pc, &patch);

3228:   if (patch->subspaces_to_exclude) {
3229:     PetscHSetIDestroy(&patch->subspaces_to_exclude);
3230:   }
3231:   PetscHSetICreate(&patch->subspaces_to_exclude);

3233:   patch->classname = "pc";
3234:   patch->isNonlinear = PETSC_FALSE;

3236:   /* Set some defaults */
3237:   patch->combined           = PETSC_FALSE;
3238:   patch->save_operators     = PETSC_TRUE;
3239:   patch->local_composition_type = PC_COMPOSITE_ADDITIVE;
3240:   patch->precomputeElementTensors = PETSC_FALSE;
3241:   patch->partition_of_unity = PETSC_FALSE;
3242:   patch->codim              = -1;
3243:   patch->dim                = -1;
3244:   patch->vankadim           = -1;
3245:   patch->ignoredim          = -1;
3246:   patch->pardecomp_overlap  = 0;
3247:   patch->patchconstructop   = PCPatchConstruct_Star;
3248:   patch->symmetrise_sweep   = PETSC_FALSE;
3249:   patch->npatch             = 0;
3250:   patch->userIS             = NULL;
3251:   patch->optionsSet         = PETSC_FALSE;
3252:   patch->iterationSet       = NULL;
3253:   patch->user_patches       = PETSC_FALSE;
3254:   PetscStrallocpy(MATDENSE, (char **) &patch->sub_mat_type);
3255:   patch->viewPatches        = PETSC_FALSE;
3256:   patch->viewCells          = PETSC_FALSE;
3257:   patch->viewPoints         = PETSC_FALSE;
3258:   patch->viewSection        = PETSC_FALSE;
3259:   patch->viewMatrix         = PETSC_FALSE;
3260:   patch->setupsolver        = PCSetUp_PATCH_Linear;
3261:   patch->applysolver        = PCApply_PATCH_Linear;
3262:   patch->resetsolver        = PCReset_PATCH_Linear;
3263:   patch->destroysolver      = PCDestroy_PATCH_Linear;
3264:   patch->updatemultiplicative = PCUpdateMultiplicative_PATCH_Linear;
3265:   patch->dofMappingWithoutToWithArtificial = NULL;
3266:   patch->dofMappingWithoutToWithAll = NULL;

3268:   pc->data                 = (void *) patch;
3269:   pc->ops->apply           = PCApply_PATCH;
3270:   pc->ops->applytranspose  = 0; /* PCApplyTranspose_PATCH; */
3271:   pc->ops->setup           = PCSetUp_PATCH;
3272:   pc->ops->reset           = PCReset_PATCH;
3273:   pc->ops->destroy         = PCDestroy_PATCH;
3274:   pc->ops->setfromoptions  = PCSetFromOptions_PATCH;
3275:   pc->ops->setuponblocks   = PCSetUpOnBlocks_PATCH;
3276:   pc->ops->view            = PCView_PATCH;
3277:   pc->ops->applyrichardson = 0;

3279:   return(0);
3280: }