Actual source code: pcpatch.c

  1: #include <petsc/private/pcpatchimpl.h>
  2: #include <petsc/private/kspimpl.h>
  3: #include <petsc/private/vecimpl.h>
  4: #include <petsc/private/dmpleximpl.h>
  5: #include <petscsf.h>
  6: #include <petscbt.h>
  7: #include <petscds.h>
  8: #include <../src/mat/impls/dense/seq/dense.h>

 10: PetscLogEvent PC_Patch_CreatePatches, PC_Patch_ComputeOp, PC_Patch_Solve, PC_Patch_Apply, PC_Patch_Prealloc;

 12: PETSC_STATIC_INLINE PetscErrorCode ObjectView(PetscObject obj, PetscViewer viewer, PetscViewerFormat format)
 13: {

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

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

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

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

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

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

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

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

 95:   PetscHSetIClear(ht);

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

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

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

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

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

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

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

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

152:   return(0);
153: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

593: /*@C

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

597:   Logically collective on PC

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

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

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

617:   Level: advanced

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

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

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

634: /*@C

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

638:   Logically collective on PC

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

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

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

658:   Level: advanced

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

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

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

675: /*@C

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

679:   Logically collective on PC

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

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

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

699:   Level: advanced

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

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

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

716: /*@C

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

720:   Logically collective on PC

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

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

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

740:   Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1110:       PetscHashIterGetKey(cht, hi, point);
1111:       if (fStart <= point && point < fEnd) {
1112:         const PetscInt *support;
1113:         PetscInt supportSize, p;
1114:         PetscBool interior = PETSC_TRUE;
1115:         DMPlexGetSupport(dm, point, &support);
1116:         DMPlexGetSupportSize(dm, point, &supportSize);
1117:         if (supportSize == 1) {
1118:           interior = PETSC_FALSE;
1119:         } else {
1120:           for (p = 0; p < supportSize; p++) {
1121:             PetscBool found;
1122:             /* FIXME: can I do this while iterating over cht? */
1123:             PetscHSetIHas(cht, support[p], &found);
1124:             if (!found) {
1125:               interior = PETSC_FALSE;
1126:               break;
1127:             }
1128:           }
1129:         }
1130:         if (interior) {
1131:           intFacetsToPatchCell[2*(ifoff + ifn)] = support[0];
1132:           intFacetsToPatchCell[2*(ifoff + ifn) + 1] = support[1];
1133:           intFacetsArray[ifoff + ifn++] = point;
1134:         } else {
1135:           extFacetsArray[efoff + efn++] = point;
1136:         }
1137:       }
1138:       PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);
1139:       if (pdof)                            {pointsArray[off + n++] = point;}
1140:       if (point >= cStart && point < cEnd) {cellsArray[coff + cn++] = point;}
1141:       PetscHashIterNext(cht, hi);
1142:     }
1143:     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);
1144:     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);
1145:     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);
1146:     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);

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1666:   PetscHMapIDestroy(&ht);
1667:   PetscHMapIDestroy(&htWithArtificial);
1668:   PetscHMapIDestroy(&htWithAll);
1669:   ISRestoreIndices(cells, &cellsArray);
1670:   ISRestoreIndices(points, &pointsArray);
1671:   PetscFree(dofsArray);
1672:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1673:     PetscFree(dofsArrayWithArtificial);
1674:   }
1675:   if (isNonlinear) {
1676:     PetscFree(dofsArrayWithAll);
1677:   }
1678:   /* Create placeholder section for map from points to patch dofs */
1679:   PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection);
1680:   PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces);
1681:   if (patch->combined) {
1682:     PetscInt numFields;
1683:     PetscSectionGetNumFields(patch->dofSection[0], &numFields);
1684:     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);
1685:     PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd);
1686:     PetscSectionSetChart(patch->patchSection, pStart, pEnd);
1687:     for (p = pStart; p < pEnd; ++p) {
1688:       PetscInt dof, fdof, f;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1886:       PetscFree(zeroes);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2416:   return(0);
2417: }

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

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

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

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

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

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

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

2524:       PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);
2525:       PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);

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

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

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

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

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

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

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

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

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

2597:             PetscSectionGetFieldDof(s, p, f, &fdof);
2598:             PetscSectionGetFieldOffset(s, p, f, &foff);
2599:             for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d;
2600:           }
2601:           DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
2602:         }
2603:         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]);
2604:       }
2605:       numGlobalBcs = 0;
2606:       for (p = pStart; p < pEnd; ++p) {
2607:         const PetscInt *ind;
2608:         PetscInt        off, cdof, d;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2707:         PetscMalloc1(numPatchDofs, &patchWithoutAllToWithAllArray);

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

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

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

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

2759:       VecDuplicate(patch->dof_weights, &global);
2760:       VecSet(global, 0.);

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

2769:       VecReciprocal(global);

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

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

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

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

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


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

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

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

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

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

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

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

2958:   PetscOptionsPopGetViewerOff();
2959:   PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0);
2960:   return(0);
2961: }

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

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

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


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

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

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

3024:   if (patch->subspaces_to_exclude) {
3025:     PetscHSetIDestroy(&patch->subspaces_to_exclude);
3026:   }

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

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

3055:   }
3056:   PetscFree(patch->sub_mat_type);
3057:   if (patch->userIS) {
3058:     for (i = 0; i < patch->npatch; ++i) {ISDestroy(&patch->userIS[i]);}
3059:     PetscFree(patch->userIS);
3060:   }
3061:   PetscFree(patch->precomputedTensorLocations);
3062:   PetscFree(patch->precomputedIntFacetTensorLocations);

3064:   patch->bs          = NULL;
3065:   patch->cellNodeMap = NULL;
3066:   patch->nsubspaces  = 0;
3067:   ISDestroy(&patch->iterationSet);

3069:   PetscViewerDestroy(&patch->viewerSection);
3070:   return(0);
3071: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3156:   /* If the user has set the number of subspaces, use that for the buffer size,
3157:      otherwise use a large number */
3158:   if (patch->nsubspaces <= 0) {
3159:     nfields = 128;
3160:   } else {
3161:     nfields = patch->nsubspaces;
3162:   }
3163:   PetscMalloc1(nfields, &ifields);
3164:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exclude_subspaces", patch->classname);
3165:   PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,option,ifields,&nfields,&flg);
3166:   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");
3167:   if (flg) {
3168:     PetscHSetIClear(patch->subspaces_to_exclude);
3169:     for (k = 0; k < nfields; k++) {
3170:       PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]);
3171:     }
3172:   }
3173:   PetscFree(ifields);

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

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

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

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

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

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

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

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

3292:   Level: intermediate

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

3302:   PetscNewLog(pc, &patch);

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

3309:   patch->classname = "pc";
3310:   patch->isNonlinear = PETSC_FALSE;

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

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

3356:   return(0);
3357: }