Actual source code: pcpatch.c

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

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

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

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

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

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

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

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

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

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

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

 88:   PetscHSetIClear(ht);

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

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

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

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

118:     PetscHSetIGetSize(ht, &htsize);
119:     PetscMalloc1(htsize, &htpoints);
120:     PetscHSetIGetElems(ht, &index, htpoints);

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

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

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

145:   return 0;
146: }

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

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

164:     if (ownedpoint < pStart || ownedpoint >= pEnd) {
165:       SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D was not in [%D, %D)", ownedpoint, pStart, pEnd);
166:     }
167:     PetscHSetIAdd(ht, ownedpoint);
168:   }
169:   ISRestoreIndices(patchis, &patchdata);
170:   return 0;
171: }

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

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

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

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

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

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

244:     offsets[0] = 0;
245:     for (i = 1; i < n; ++i) {
246:       PetscInt nroots;

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

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

270:       PetscSFGetGraph(sf[i], &nroots, &nleaves, &local, &remote);
271:       for (j = 0; j < nleaves; ++j) {
272:         PetscInt rank = remote[j].rank;
273:         PetscInt idx, rootOffset, k;

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

296: PetscErrorCode PCPatchSetDenseInverse(PC pc, PetscBool flg)
297: {
298:   PC_PATCH *patch = (PC_PATCH *) pc->data;
299:   patch->denseinverse = flg;
300:   return 0;
301: }

303: PetscErrorCode PCPatchGetDenseInverse(PC pc, PetscBool *flg)
304: {
305:   PC_PATCH *patch = (PC_PATCH *) pc->data;
306:   *flg = patch->denseinverse;
307:   return 0;
308: }

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

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

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

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

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

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

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

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

374: /* TODO: Docs */
375: PetscErrorCode PCPatchSetLocalComposition(PC pc, PCCompositeType type)
376: {
377:   PC_PATCH *patch = (PC_PATCH *) pc->data;
379:   patch->local_composition_type = type;
380:   return 0;
381: }

383: /* TODO: Docs */
384: PetscErrorCode PCPatchGetLocalComposition(PC pc, PCCompositeType *type)
385: {
386:   PC_PATCH *patch = (PC_PATCH *) pc->data;
387:   *type = patch->local_composition_type;
388:   return 0;
389: }

391: /* TODO: Docs */
392: PetscErrorCode PCPatchSetSubMatType(PC pc, MatType sub_mat_type)
393: {
394:   PC_PATCH      *patch = (PC_PATCH *) pc->data;

396:   if (patch->sub_mat_type) PetscFree(patch->sub_mat_type);
397:   PetscStrallocpy(sub_mat_type, (char **) &patch->sub_mat_type);
398:   return 0;
399: }

401: /* TODO: Docs */
402: PetscErrorCode PCPatchGetSubMatType(PC pc, MatType *sub_mat_type)
403: {
404:   PC_PATCH *patch = (PC_PATCH *) pc->data;
405:   *sub_mat_type = patch->sub_mat_type;
406:   return 0;
407: }

409: /* TODO: Docs */
410: PetscErrorCode PCPatchSetCellNumbering(PC pc, PetscSection cellNumbering)
411: {
412:   PC_PATCH      *patch = (PC_PATCH *) pc->data;

414:   patch->cellNumbering = cellNumbering;
415:   PetscObjectReference((PetscObject) cellNumbering);
416:   return 0;
417: }

419: /* TODO: Docs */
420: PetscErrorCode PCPatchGetCellNumbering(PC pc, PetscSection *cellNumbering)
421: {
422:   PC_PATCH *patch = (PC_PATCH *) pc->data;
423:   *cellNumbering = patch->cellNumbering;
424:   return 0;
425: }

427: /* TODO: Docs */
428: PetscErrorCode PCPatchSetConstructType(PC pc, PCPatchConstructType ctype, PetscErrorCode (*func)(PC, PetscInt *, IS **, IS *, void *), void *ctx)
429: {
430:   PC_PATCH *patch = (PC_PATCH *) pc->data;

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

461: /* TODO: Docs */
462: PetscErrorCode PCPatchGetConstructType(PC pc, PCPatchConstructType *ctype, PetscErrorCode (**func)(PC, PetscInt *, IS **, IS *, void *), void **ctx)
463: {
464:   PC_PATCH *patch = (PC_PATCH *) pc->data;

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

483: /* TODO: Docs */
484: PetscErrorCode PCPatchSetDiscretisationInfo(PC pc, PetscInt nsubspaces, DM *dms, PetscInt *bs, PetscInt *nodesPerCell, const PetscInt **cellNodeMap,
485:                                             const PetscInt *subspaceOffsets, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes)
486: {
487:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
488:   DM             dm, plex;
489:   PetscSF       *sfs;
490:   PetscInt       cStart, cEnd, i, j;

492:   PCGetDM(pc, &dm);
493:   DMConvert(dm, DMPLEX, &plex);
494:   dm = plex;
495:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
496:   PetscMalloc1(nsubspaces, &sfs);
497:   PetscMalloc1(nsubspaces, &patch->dofSection);
498:   PetscMalloc1(nsubspaces, &patch->bs);
499:   PetscMalloc1(nsubspaces, &patch->nodesPerCell);
500:   PetscMalloc1(nsubspaces, &patch->cellNodeMap);
501:   PetscMalloc1(nsubspaces+1, &patch->subspaceOffsets);

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

519:   patch->subspaceOffsets[nsubspaces] = subspaceOffsets[nsubspaces];
520:   ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);
521:   ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);
522:   DMDestroy(&dm);
523:   return 0;
524: }

526: /* TODO: Docs */
527: PetscErrorCode PCPatchSetDiscretisationInfoCombined(PC pc, DM dm, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes)
528: {
529:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
530:   PetscInt       cStart, cEnd, i, j;

532:   patch->combined = PETSC_TRUE;
533:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
534:   DMGetNumFields(dm, &patch->nsubspaces);
535:   PetscCalloc1(patch->nsubspaces, &patch->dofSection);
536:   PetscMalloc1(patch->nsubspaces, &patch->bs);
537:   PetscMalloc1(patch->nsubspaces, &patch->nodesPerCell);
538:   PetscMalloc1(patch->nsubspaces, &patch->cellNodeMap);
539:   PetscCalloc1(patch->nsubspaces+1, &patch->subspaceOffsets);
540:   DMGetLocalSection(dm, &patch->dofSection[0]);
541:   PetscObjectReference((PetscObject) patch->dofSection[0]);
542:   PetscSectionGetStorageSize(patch->dofSection[0], &patch->subspaceOffsets[patch->nsubspaces]);
543:   patch->totalDofsPerCell = 0;
544:   for (i = 0; i < patch->nsubspaces; ++i) {
545:     patch->bs[i]             = 1;
546:     patch->nodesPerCell[i]   = nodesPerCell[i];
547:     patch->totalDofsPerCell += nodesPerCell[i];
548:     PetscMalloc1((cEnd-cStart)*nodesPerCell[i], &patch->cellNodeMap[i]);
549:     for (j = 0; j < (cEnd-cStart)*nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j];
550:   }
551:   DMGetSectionSF(dm, &patch->sectionSF);
552:   PetscObjectReference((PetscObject) patch->sectionSF);
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:   return 0;
556: }

558: /*@C

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

562:   Logically collective on PC

564:   Input Parameters:
565: + pc   - The PC
566: . func - The callback
567: - ctx  - The user context

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

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

582:   Level: advanced

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

587: .seealso: PCPatchSetComputeOperator(), PCPatchGetComputeOperator(), PCPatchSetDiscretisationInfo(), PCPatchSetComputeFunctionInteriorFacets()
588: @*/
589: PetscErrorCode PCPatchSetComputeFunction(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Vec, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx)
590: {
591:   PC_PATCH *patch = (PC_PATCH *) pc->data;

593:   patch->usercomputef    = func;
594:   patch->usercomputefctx = ctx;
595:   return 0;
596: }

598: /*@C

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

602:   Logically collective on PC

604:   Input Parameters:
605: + pc   - The PC
606: . func - The callback
607: - ctx  - The user context

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

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

622:   Level: advanced

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

627: .seealso: PCPatchSetComputeOperator(), PCPatchGetComputeOperator(), PCPatchSetDiscretisationInfo(), PCPatchSetComputeFunction()
628: @*/
629: PetscErrorCode PCPatchSetComputeFunctionInteriorFacets(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Vec, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx)
630: {
631:   PC_PATCH *patch = (PC_PATCH *) pc->data;

633:   patch->usercomputefintfacet    = func;
634:   patch->usercomputefintfacetctx = ctx;
635:   return 0;
636: }

638: /*@C

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

642:   Logically collective on PC

644:   Input Parameters:
645: + pc   - The PC
646: . func - The callback
647: - ctx  - The user context

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

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

662:   Level: advanced

664:   Notes:
665:   The matrix entries have been set to zero before the call.

667: .seealso: PCPatchGetComputeOperator(), PCPatchSetComputeFunction(), PCPatchSetDiscretisationInfo()
668: @*/
669: PetscErrorCode PCPatchSetComputeOperator(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Mat, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx)
670: {
671:   PC_PATCH *patch = (PC_PATCH *) pc->data;

673:   patch->usercomputeop    = func;
674:   patch->usercomputeopctx = ctx;
675:   return 0;
676: }

678: /*@C

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

682:   Logically collective on PC

684:   Input Parameters:
685: + pc   - The PC
686: . func - The callback
687: - ctx  - The user context

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

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

702:   Level: advanced

704:   Notes:
705:   The matrix entries have been set to zero before the call.

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

713:   patch->usercomputeopintfacet    = func;
714:   patch->usercomputeopintfacetctx = ctx;
715:   return 0;
716: }

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

734:   PCGetDM(pc, &dm);
735:   DMConvert(dm, DMPLEX, &plex);
736:   dm = plex;
737:   DMPlexGetHeightStratum(dm, 1, &fBegin, &fEnd);
738:   PCPatchGetIgnoreDim(pc, &ignoredim);
739:   if (ignoredim >= 0) DMPlexGetDepthStratum(dm, ignoredim, &iStart, &iEnd);
740:   PetscHSetIClear(cht);
741:   PetscHashIterBegin(ht, hi);
742:   while (!PetscHashIterAtEnd(ht, hi)) {

744:     PetscHashIterGetKey(ht, hi, point);
745:     PetscHashIterNext(ht, hi);

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

783: static PetscErrorCode PCPatchGetGlobalDofs(PC pc, PetscSection dofSection[], PetscInt f, PetscBool combined, PetscInt p, PetscInt *dof, PetscInt *off)
784: {
785:   if (combined) {
786:     if (f < 0) {
787:       if (dof) PetscSectionGetDof(dofSection[0], p, dof);
788:       if (off) PetscSectionGetOffset(dofSection[0], p, off);
789:     } else {
790:       if (dof) PetscSectionGetFieldDof(dofSection[0], p, f, dof);
791:       if (off) PetscSectionGetFieldOffset(dofSection[0], p, f, off);
792:     }
793:   } else {
794:     if (f < 0) {
795:       PC_PATCH *patch = (PC_PATCH *) pc->data;
796:       PetscInt  fdof, g;

798:       if (dof) {
799:         *dof = 0;
800:         for (g = 0; g < patch->nsubspaces; ++g) {
801:           PetscSectionGetDof(dofSection[g], p, &fdof);
802:           *dof += fdof;
803:         }
804:       }
805:       if (off) {
806:         *off = 0;
807:         for (g = 0; g < patch->nsubspaces; ++g) {
808:           PetscSectionGetOffset(dofSection[g], p, &fdof);
809:           *off += fdof;
810:         }
811:       }
812:     } else {
813:       if (dof) PetscSectionGetDof(dofSection[f], p, dof);
814:       if (off) PetscSectionGetOffset(dofSection[f], p, off);
815:     }
816:   }
817:   return 0;
818: }

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

832:   PetscHSetIClear(dofs);
833:   for (k = 0; k < patch->nsubspaces; ++k) {
834:     PetscInt subspaceOffset = patch->subspaceOffsets[k];
835:     PetscInt bs             = patch->bs[k];
836:     PetscInt j, l;

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

855:     PetscHashIterBegin(pts, hi);
856:     while (!PetscHashIterAtEnd(pts, hi)) {
857:       PetscHashIterGetKey(pts, hi, p);
858:       PetscHashIterNext(pts, hi);
859:       PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, p, &ldof, &loff);
860:       if (0 == ldof) continue;
861:       for (j = loff; j < ldof + loff; ++j) {
862:         for (l = 0; l < bs; ++l) {
863:           PetscInt dof = bs*j + l + subspaceOffset;
864:           PetscHSetIAdd(dofs, dof);
865:         }
866:       }
867:     }
868:   }
869:   return 0;
870: }

872: /* Given two hash tables A and B, compute the keys in B that are not in A, and put them in C */
873: static PetscErrorCode PCPatchComputeSetDifference_Private(PetscHSetI A, PetscHSetI B, PetscHSetI C)
874: {
875:   PetscHashIter  hi;
876:   PetscInt       key;
877:   PetscBool      flg;

879:   PetscHSetIClear(C);
880:   PetscHashIterBegin(B, hi);
881:   while (!PetscHashIterAtEnd(B, hi)) {
882:     PetscHashIterGetKey(B, hi, key);
883:     PetscHashIterNext(B, hi);
884:     PetscHSetIHas(A, key, &flg);
885:     if (!flg) PetscHSetIAdd(C, key);
886:   }
887:   return 0;
888: }

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

915:   /* Used to keep track of the cells in the patch. */
916:   PetscHSetICreate(&ht);
917:   PetscHSetICreate(&cht);

919:   PCGetDM(pc, &dm);
921:   DMConvert(dm, DMPLEX, &plex);
922:   dm = plex;
923:   DMPlexGetChart(dm, &pStart, &pEnd);
924:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);

926:   if (patch->user_patches) {
927:     patch->userpatchconstructionop(pc, &patch->npatch, &patch->userIS, &patch->iterationSet, patch->userpatchconstructctx);
928:     vStart = 0; vEnd = patch->npatch;
929:   } else if (patch->ctype == PC_PATCH_PARDECOMP) {
930:     vStart = 0; vEnd = 1;
931:   } else if (patch->codim < 0) {
932:     if (patch->dim < 0) DMPlexGetDepthStratum(dm,  0,            &vStart, &vEnd);
933:     else                DMPlexGetDepthStratum(dm,  patch->dim,   &vStart, &vEnd);
934:   } else                DMPlexGetHeightStratum(dm, patch->codim, &vStart, &vEnd);
935:   patch->npatch = vEnd - vStart;

937:   /* These labels mark the owned points.  We only create patches around points that this process owns. */
938:   DMHasLabel(dm, "pyop2_ghost", &isFiredrake);
939:   if (isFiredrake) {
940:     DMGetLabel(dm, "pyop2_ghost", &ghost);
941:     DMLabelCreateIndex(ghost, pStart, pEnd);
942:   } else {
943:     PetscSF sf;

945:     DMGetPointSF(dm, &sf);
946:     PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
947:     nleaves = PetscMax(nleaves, 0);
948:   }

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

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

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

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

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

1025:   PetscSectionSetUp(cellCounts);
1026:   PetscSectionGetStorageSize(cellCounts, &numCells);
1027:   PetscMalloc1(numCells, &cellsArray);
1028:   PetscSectionSetUp(pointCounts);
1029:   PetscSectionGetStorageSize(pointCounts, &numPoints);
1030:   PetscMalloc1(numPoints, &pointsArray);

1032:   PetscSectionSetUp(intFacetCounts);
1033:   PetscSectionSetUp(extFacetCounts);
1034:   PetscSectionGetStorageSize(intFacetCounts, &numIntFacets);
1035:   PetscSectionGetStorageSize(extFacetCounts, &numExtFacets);
1036:   PetscMalloc1(numIntFacets, &intFacetsArray);
1037:   PetscMalloc1(numIntFacets*2, &intFacetsToPatchCell);
1038:   PetscMalloc1(numExtFacets, &extFacetsArray);

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

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

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

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

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

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


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

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

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

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

1260:   /* Hash tables for artificial BC construction */
1261:   PetscHSetICreate(&ownedpts);
1262:   PetscHSetICreate(&seenpts);
1263:   PetscHSetICreate(&owneddofs);
1264:   PetscHSetICreate(&seendofs);
1265:   PetscHSetICreate(&artificialbcs);

1267:   ISGetIndices(cells, &cellsArray);
1268:   ISGetIndices(points, &pointsArray);
1269:   PetscHMapICreate(&ht);
1270:   PetscHMapICreate(&htWithArtificial);
1271:   PetscHMapICreate(&htWithAll);
1272:   for (v = vStart; v < vEnd; ++v) {
1273:     PetscInt localIndex = 0;
1274:     PetscInt localIndexWithArtificial = 0;
1275:     PetscInt localIndexWithAll = 0;
1276:     PetscInt dof, off, i, j, k, l;

1278:     PetscHMapIClear(ht);
1279:     PetscHMapIClear(htWithArtificial);
1280:     PetscHMapIClear(htWithAll);
1281:     PetscSectionGetDof(cellCounts, v, &dof);
1282:     PetscSectionGetOffset(cellCounts, v, &off);
1283:     if (dof <= 0) continue;

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

1296:       PetscHSetICreate(&globalbcdofs);
1297:       PetscSynchronizedPrintf(comm, "Patch %d: owned dofs:\n", v);
1298:       PetscHashIterBegin(owneddofs, hi);
1299:       while (!PetscHashIterAtEnd(owneddofs, hi)) {
1300:         PetscInt globalDof;

1302:         PetscHashIterGetKey(owneddofs, hi, globalDof);
1303:         PetscHashIterNext(owneddofs, hi);
1304:         PetscSynchronizedPrintf(comm, "%d ", globalDof);
1305:       }
1306:       PetscSynchronizedPrintf(comm, "\n");
1307:       PetscSynchronizedPrintf(comm, "Patch %d: seen dofs:\n", v);
1308:       PetscHashIterBegin(seendofs, hi);
1309:       while (!PetscHashIterAtEnd(seendofs, hi)) {
1310:         PetscInt globalDof;
1311:         PetscBool flg;

1313:         PetscHashIterGetKey(seendofs, hi, globalDof);
1314:         PetscHashIterNext(seendofs, hi);
1315:         PetscSynchronizedPrintf(comm, "%d ", globalDof);

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

1353:       for (i = off; i < off + dof; ++i) {
1354:         /* Walk over the cells in this patch. */
1355:         const PetscInt c    = cellsArray[i];
1356:         PetscInt       cell = c;

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

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

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

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

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

1436:   DMDestroy(&dm);
1438:   PetscSectionSetUp(gtolCounts);
1439:   PetscSectionGetStorageSize(gtolCounts, &numGlobalDofs);
1440:   PetscMalloc1(numGlobalDofs, &globalDofsArray);

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

1457:     PetscHMapIClear(ht);
1458:     PetscHMapIClear(htWithArtificial);
1459:     PetscHMapIClear(htWithAll);
1460:     PetscSectionGetDof(cellCounts, v, &dof);
1461:     PetscSectionGetOffset(cellCounts, v, &off);
1462:     PetscSectionGetDof(pointCounts, v, &Np);
1463:     PetscSectionGetOffset(pointCounts, v, &ooff);
1464:     if (dof <= 0) continue;

1466:     for (k = 0; k < patch->nsubspaces; ++k) {
1467:       const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
1468:       PetscInt        nodesPerCell   = patch->nodesPerCell[k];
1469:       PetscInt        subspaceOffset = patch->subspaceOffsets[k];
1470:       PetscInt        bs             = patch->bs[k];
1471:       PetscInt        goff;

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

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

1501:       /* Shove it in the output data structure. */
1502:       PetscSectionGetOffset(gtolCounts, v, &goff);
1503:       PetscHashIterBegin(ht, hi);
1504:       while (!PetscHashIterAtEnd(ht, hi)) {
1505:         PetscInt globalDof, localDof;

1507:         PetscHashIterGetKey(ht, hi, globalDof);
1508:         PetscHashIterGetVal(ht, hi, localDof);
1509:         if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof;
1510:         PetscHashIterNext(ht, hi);
1511:       }

1513:       if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1514:         PetscSectionGetOffset(gtolCountsWithArtificial, v, &goff);
1515:         PetscHashIterBegin(htWithArtificial, hi);
1516:         while (!PetscHashIterAtEnd(htWithArtificial, hi)) {
1517:           PetscInt globalDof, localDof;
1518:           PetscHashIterGetKey(htWithArtificial, hi, globalDof);
1519:           PetscHashIterGetVal(htWithArtificial, hi, localDof);
1520:           if (globalDof >= 0) globalDofsArrayWithArtificial[goff + localDof] = globalDof;
1521:           PetscHashIterNext(htWithArtificial, hi);
1522:         }
1523:       }
1524:       if (isNonlinear) {
1525:         PetscSectionGetOffset(gtolCountsWithAll, v, &goff);
1526:         PetscHashIterBegin(htWithAll, hi);
1527:         while (!PetscHashIterAtEnd(htWithAll, hi)) {
1528:           PetscInt globalDof, localDof;
1529:           PetscHashIterGetKey(htWithAll, hi, globalDof);
1530:           PetscHashIterGetVal(htWithAll, hi, localDof);
1531:           if (globalDof >= 0) globalDofsArrayWithAll[goff + localDof] = globalDof;
1532:           PetscHashIterNext(htWithAll, hi);
1533:         }
1534:       }

1536:       for (p = 0; p < Np; ++p) {
1537:         const PetscInt point = pointsArray[ooff + p];
1538:         PetscInt       globalDof, localDof;

1540:         PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof);
1541:         PetscHMapIGet(ht, globalDof, &localDof);
1542:         offsArray[(ooff + p)*Nf + k] = localDof;
1543:         if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1544:           PetscHMapIGet(htWithArtificial, globalDof, &localDof);
1545:           offsArrayWithArtificial[(ooff + p)*Nf + k] = localDof;
1546:         }
1547:         if (isNonlinear) {
1548:           PetscHMapIGet(htWithAll, globalDof, &localDof);
1549:           offsArrayWithAll[(ooff + p)*Nf + k] = localDof;
1550:         }
1551:       }
1552:     }

1554:     PetscHSetIDestroy(&globalBcs);
1555:     PetscHSetIDestroy(&ownedpts);
1556:     PetscHSetIDestroy(&seenpts);
1557:     PetscHSetIDestroy(&owneddofs);
1558:     PetscHSetIDestroy(&seendofs);
1559:     PetscHSetIDestroy(&artificialbcs);

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

1571:         if (cellNumbering) PetscSectionGetOffset(cellNumbering, c, &cell);
1572:         for (k = 0; k < patch->nsubspaces; ++k) {
1573:           const PetscInt *cellNodeMap    = patch->cellNodeMap[k];
1574:           PetscInt        nodesPerCell   = patch->nodesPerCell[k];
1575:           PetscInt        subspaceOffset = patch->subspaceOffsets[k];
1576:           PetscInt        bs             = patch->bs[k];

1578:           for (j = 0; j < nodesPerCell; ++j) {
1579:             for (l = 0; l < bs; ++l) {
1580:               const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset;
1581:               PetscInt       localDof;

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

1614:   PetscHMapIDestroy(&ht);
1615:   PetscHMapIDestroy(&htWithArtificial);
1616:   PetscHMapIDestroy(&htWithAll);
1617:   ISRestoreIndices(cells, &cellsArray);
1618:   ISRestoreIndices(points, &pointsArray);
1619:   PetscFree(dofsArray);
1620:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1621:     PetscFree(dofsArrayWithArtificial);
1622:   }
1623:   if (isNonlinear) {
1624:     PetscFree(dofsArrayWithAll);
1625:   }
1626:   /* Create placeholder section for map from points to patch dofs */
1627:   PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection);
1628:   PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces);
1629:   if (patch->combined) {
1630:     PetscInt numFields;
1631:     PetscSectionGetNumFields(patch->dofSection[0], &numFields);
1633:     PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd);
1634:     PetscSectionSetChart(patch->patchSection, pStart, pEnd);
1635:     for (p = pStart; p < pEnd; ++p) {
1636:       PetscInt dof, fdof, f;

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

1689: static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat, PetscBool withArtificial)
1690: {
1691:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
1692:   PetscBool      flg;
1693:   PetscInt       csize, rsize;
1694:   const char    *prefix = NULL;

1696:   if (withArtificial) {
1697:     /* would be nice if we could create a rectangular matrix of size numDofsWithArtificial x numDofs here */
1698:     PetscInt pStart;
1699:     PetscSectionGetChart(patch->gtolCountsWithArtificial, &pStart, NULL);
1700:     PetscSectionGetDof(patch->gtolCountsWithArtificial, point + pStart, &rsize);
1701:     csize = rsize;
1702:   } else {
1703:     PetscInt pStart;
1704:     PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);
1705:     PetscSectionGetDof(patch->gtolCounts, point + pStart, &rsize);
1706:     csize = rsize;
1707:   }

1709:   MatCreate(PETSC_COMM_SELF, mat);
1710:   PCGetOptionsPrefix(pc, &prefix);
1711:   MatSetOptionsPrefix(*mat, prefix);
1712:   MatAppendOptionsPrefix(*mat, "pc_patch_sub_");
1713:   if (patch->sub_mat_type)       MatSetType(*mat, patch->sub_mat_type);
1714:   else if (!patch->sub_mat_type) MatSetType(*mat, MATDENSE);
1715:   MatSetSizes(*mat, rsize, csize, rsize, csize);
1716:   PetscObjectTypeCompare((PetscObject) *mat, MATDENSE, &flg);
1717:   if (!flg) PetscObjectTypeCompare((PetscObject)*mat, MATSEQDENSE, &flg);
1718:   /* Sparse patch matrices */
1719:   if (!flg) {
1720:     PetscBT         bt;
1721:     PetscInt       *dnnz      = NULL;
1722:     const PetscInt *dofsArray = NULL;
1723:     PetscInt        pStart, pEnd, ncell, offset, c, i, j;

1725:     if (withArtificial) {
1726:       ISGetIndices(patch->dofsWithArtificial, &dofsArray);
1727:     } else {
1728:       ISGetIndices(patch->dofs, &dofsArray);
1729:     }
1730:     PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);
1731:     point += pStart;
1733:     PetscSectionGetDof(patch->cellCounts, point, &ncell);
1734:     PetscSectionGetOffset(patch->cellCounts, point, &offset);
1735:     PetscLogEventBegin(PC_Patch_Prealloc, pc, 0, 0, 0);
1736:     /* A PetscBT uses N^2 bits to store the sparsity pattern on a
1737:      * patch. This is probably OK if the patches are not too big,
1738:      * but uses too much memory. We therefore switch based on rsize. */
1739:     if (rsize < 3000) { /* FIXME: I picked this switch value out of my hat */
1740:       PetscScalar *zeroes;
1741:       PetscInt rows;

1743:       PetscCalloc1(rsize, &dnnz);
1744:       PetscBTCreate(rsize*rsize, &bt);
1745:       for (c = 0; c < ncell; ++c) {
1746:         const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell;
1747:         for (i = 0; i < patch->totalDofsPerCell; ++i) {
1748:           const PetscInt row = idx[i];
1749:           if (row < 0) continue;
1750:           for (j = 0; j < patch->totalDofsPerCell; ++j) {
1751:             const PetscInt col = idx[j];
1752:             const PetscInt key = row*rsize + col;
1753:             if (col < 0) continue;
1754:             if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1755:           }
1756:         }
1757:       }

1759:       if (patch->usercomputeopintfacet) {
1760:         const PetscInt *intFacetsArray = NULL;
1761:         PetscInt        i, numIntFacets, intFacetOffset;
1762:         const PetscInt *facetCells = NULL;

1764:         PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1765:         PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1766:         ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1767:         ISGetIndices(patch->intFacets, &intFacetsArray);
1768:         for (i = 0; i < numIntFacets; i++) {
1769:           const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0];
1770:           const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1];
1771:           PetscInt       celli, cellj;

1773:           for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1774:             const PetscInt row = dofsArray[(offset + cell0)*patch->totalDofsPerCell + celli];
1775:             if (row < 0) continue;
1776:             for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1777:               const PetscInt col = dofsArray[(offset + cell1)*patch->totalDofsPerCell + cellj];
1778:               const PetscInt key = row*rsize + col;
1779:               if (col < 0) continue;
1780:               if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1781:             }
1782:           }

1784:           for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1785:             const PetscInt row = dofsArray[(offset + cell1)*patch->totalDofsPerCell + celli];
1786:             if (row < 0) continue;
1787:             for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1788:               const PetscInt col = dofsArray[(offset + cell0)*patch->totalDofsPerCell + cellj];
1789:               const PetscInt key = row*rsize + col;
1790:               if (col < 0) continue;
1791:               if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1792:             }
1793:           }
1794:         }
1795:       }
1796:       PetscBTDestroy(&bt);
1797:       MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL);
1798:       PetscFree(dnnz);

1800:       PetscCalloc1(patch->totalDofsPerCell*patch->totalDofsPerCell, &zeroes);
1801:       for (c = 0; c < ncell; ++c) {
1802:         const PetscInt *idx = &dofsArray[(offset + c)*patch->totalDofsPerCell];
1803:         MatSetValues(*mat, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, zeroes, INSERT_VALUES);
1804:       }
1805:       MatGetLocalSize(*mat, &rows, NULL);
1806:       for (i = 0; i < rows; ++i) {
1807:         MatSetValues(*mat, 1, &i, 1, &i, zeroes, INSERT_VALUES);
1808:       }

1810:       if (patch->usercomputeopintfacet) {
1811:         const PetscInt *intFacetsArray = NULL;
1812:         PetscInt i, numIntFacets, intFacetOffset;
1813:         const PetscInt *facetCells = NULL;

1815:         PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1816:         PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1817:         ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1818:         ISGetIndices(patch->intFacets, &intFacetsArray);
1819:         for (i = 0; i < numIntFacets; i++) {
1820:           const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0];
1821:           const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1];
1822:           const PetscInt *cell0idx = &dofsArray[(offset + cell0)*patch->totalDofsPerCell];
1823:           const PetscInt *cell1idx = &dofsArray[(offset + cell1)*patch->totalDofsPerCell];
1824:           MatSetValues(*mat, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, zeroes, INSERT_VALUES);
1825:           MatSetValues(*mat, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, zeroes, INSERT_VALUES);
1826:         }
1827:       }

1829:       MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
1830:       MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);

1832:       PetscFree(zeroes);

1834:     } else { /* rsize too big, use MATPREALLOCATOR */
1835:       Mat preallocator;
1836:       PetscScalar* vals;

1838:       PetscCalloc1(patch->totalDofsPerCell*patch->totalDofsPerCell, &vals);
1839:       MatCreate(PETSC_COMM_SELF, &preallocator);
1840:       MatSetType(preallocator, MATPREALLOCATOR);
1841:       MatSetSizes(preallocator, rsize, rsize, rsize, rsize);
1842:       MatSetUp(preallocator);

1844:       for (c = 0; c < ncell; ++c) {
1845:         const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell;
1846:         MatSetValues(preallocator, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, vals, INSERT_VALUES);
1847:       }

1849:       if (patch->usercomputeopintfacet) {
1850:         const PetscInt *intFacetsArray = NULL;
1851:         PetscInt        i, numIntFacets, intFacetOffset;
1852:         const PetscInt *facetCells = NULL;

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

1868:       PetscFree(vals);
1869:       MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1870:       MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1871:       MatPreallocatorPreallocate(preallocator, PETSC_TRUE, *mat);
1872:       MatDestroy(&preallocator);
1873:     }
1874:     PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0);
1875:     if (withArtificial) {
1876:       ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);
1877:     } else {
1878:       ISRestoreIndices(patch->dofs, &dofsArray);
1879:     }
1880:   }
1881:   MatSetUp(*mat);
1882:   return 0;
1883: }

1885: 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)
1886: {
1887:   PC_PATCH       *patch = (PC_PATCH *) pc->data;
1888:   DM              dm, plex;
1889:   PetscSection    s;
1890:   const PetscInt *parray, *oarray;
1891:   PetscInt        Nf = patch->nsubspaces, Np, poff, p, f;

1894:   PCGetDM(pc, &dm);
1895:   DMConvert(dm, DMPLEX, &plex);
1896:   dm = plex;
1897:   DMGetLocalSection(dm, &s);
1898:   /* Set offset into patch */
1899:   PetscSectionGetDof(patch->pointCounts, patchNum, &Np);
1900:   PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);
1901:   ISGetIndices(patch->points, &parray);
1902:   ISGetIndices(patch->offs,   &oarray);
1903:   for (f = 0; f < Nf; ++f) {
1904:     for (p = 0; p < Np; ++p) {
1905:       const PetscInt point = parray[poff+p];
1906:       PetscInt       dof;

1908:       PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);
1909:       PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);
1910:       if (patch->nsubspaces == 1) PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);
1911:       else                        PetscSectionSetOffset(patch->patchSection, point, -1);
1912:     }
1913:   }
1914:   ISRestoreIndices(patch->points, &parray);
1915:   ISRestoreIndices(patch->offs,   &oarray);
1916:   if (patch->viewSection) ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);
1917:   DMPlexComputeResidual_Patch_Internal(dm, patch->patchSection, cellIS, 0.0, x, NULL, F, ctx);
1918:   DMDestroy(&dm);
1919:   return 0;
1920: }

1922: PetscErrorCode PCPatchComputeFunction_Internal(PC pc, Vec x, Vec F, PetscInt point)
1923: {
1924:   PC_PATCH       *patch = (PC_PATCH *) pc->data;
1925:   const PetscInt *dofsArray;
1926:   const PetscInt *dofsArrayWithAll;
1927:   const PetscInt *cellsArray;
1928:   PetscInt        ncell, offset, pStart, pEnd;
1929:   PetscErrorCode  ierr;

1931:   PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);
1933:   ISGetIndices(patch->dofs, &dofsArray);
1934:   ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);
1935:   ISGetIndices(patch->cells, &cellsArray);
1936:   PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);

1938:   point += pStart;

1941:   PetscSectionGetDof(patch->cellCounts, point, &ncell);
1942:   PetscSectionGetOffset(patch->cellCounts, point, &offset);
1943:   if (ncell <= 0) {
1944:     PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
1945:     return 0;
1946:   }
1947:   VecSet(F, 0.0);
1948:   PetscStackPush("PCPatch user callback");
1949:   /* Cannot reuse the same IS because the geometry info is being cached in it */
1950:   ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);
1951:   patch->usercomputef(pc, point, x, F, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell,
1952:                                                                                             dofsArrayWithAll + offset*patch->totalDofsPerCell,
1953:                                                                                             patch->usercomputefctx);
1954:   PetscStackPop;
1955:   ISDestroy(&patch->cellIS);
1956:   ISRestoreIndices(patch->dofs, &dofsArray);
1957:   ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);
1958:   ISRestoreIndices(patch->cells, &cellsArray);
1959:   if (patch->viewMatrix) {
1960:     char name[PETSC_MAX_PATH_LEN];

1962:     PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch vector for Point %D", point);
1963:     PetscObjectSetName((PetscObject) F, name);
1964:     ObjectView((PetscObject) F, patch->viewerMatrix, patch->formatMatrix);
1965:   }
1966:   PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
1967:   return 0;
1968: }

1970: 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)
1971: {
1972:   PC_PATCH       *patch = (PC_PATCH *) pc->data;
1973:   DM              dm, plex;
1974:   PetscSection    s;
1975:   const PetscInt *parray, *oarray;
1976:   PetscInt        Nf = patch->nsubspaces, Np, poff, p, f;

1978:   PCGetDM(pc, &dm);
1979:   DMConvert(dm, DMPLEX, &plex);
1980:   dm = plex;
1981:   DMGetLocalSection(dm, &s);
1982:   /* Set offset into patch */
1983:   PetscSectionGetDof(patch->pointCounts, patchNum, &Np);
1984:   PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);
1985:   ISGetIndices(patch->points, &parray);
1986:   ISGetIndices(patch->offs,   &oarray);
1987:   for (f = 0; f < Nf; ++f) {
1988:     for (p = 0; p < Np; ++p) {
1989:       const PetscInt point = parray[poff+p];
1990:       PetscInt       dof;

1992:       PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);
1993:       PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);
1994:       if (patch->nsubspaces == 1) PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);
1995:       else                        PetscSectionSetOffset(patch->patchSection, point, -1);
1996:     }
1997:   }
1998:   ISRestoreIndices(patch->points, &parray);
1999:   ISRestoreIndices(patch->offs,   &oarray);
2000:   if (patch->viewSection) ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);
2001:   /* TODO Shut off MatViewFromOptions() in MatAssemblyEnd() here */
2002:   DMPlexComputeJacobian_Patch_Internal(dm, patch->patchSection, patch->patchSection, cellIS, 0.0, 0.0, x, NULL, J, J, ctx);
2003:   DMDestroy(&dm);
2004:   return 0;
2005: }

2007: /* This function zeros mat on entry */
2008: PetscErrorCode PCPatchComputeOperator_Internal(PC pc, Vec x, Mat mat, PetscInt point, PetscBool withArtificial)
2009: {
2010:   PC_PATCH       *patch = (PC_PATCH *) pc->data;
2011:   const PetscInt *dofsArray;
2012:   const PetscInt *dofsArrayWithAll = NULL;
2013:   const PetscInt *cellsArray;
2014:   PetscInt        ncell, offset, pStart, pEnd, numIntFacets, intFacetOffset;
2015:   PetscBool       isNonlinear;

2017:   PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);
2018:   isNonlinear = patch->isNonlinear;
2020:   if (withArtificial) {
2021:     ISGetIndices(patch->dofsWithArtificial, &dofsArray);
2022:   } else {
2023:     ISGetIndices(patch->dofs, &dofsArray);
2024:   }
2025:   if (isNonlinear) {
2026:     ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);
2027:   }
2028:   ISGetIndices(patch->cells, &cellsArray);
2029:   PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);

2031:   point += pStart;

2034:   PetscSectionGetDof(patch->cellCounts, point, &ncell);
2035:   PetscSectionGetOffset(patch->cellCounts, point, &offset);
2036:   if (ncell <= 0) {
2037:     PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2038:     return 0;
2039:   }
2040:   MatZeroEntries(mat);
2041:   if (patch->precomputeElementTensors) {
2042:     PetscInt           i;
2043:     PetscInt           ndof = patch->totalDofsPerCell;
2044:     const PetscScalar *elementTensors;

2046:     VecGetArrayRead(patch->cellMats, &elementTensors);
2047:     for (i = 0; i < ncell; i++) {
2048:       const PetscInt     cell = cellsArray[i + offset];
2049:       const PetscInt    *idx  = dofsArray + (offset + i)*ndof;
2050:       const PetscScalar *v    = elementTensors + patch->precomputedTensorLocations[cell]*ndof*ndof;
2051:       MatSetValues(mat, ndof, idx, ndof, idx, v, ADD_VALUES);
2052:     }
2053:     VecRestoreArrayRead(patch->cellMats, &elementTensors);
2054:     MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
2055:     MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
2056:   } else {
2057:     PetscStackPush("PCPatch user callback");
2058:     /* Cannot reuse the same IS because the geometry info is being cached in it */
2059:     ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);
2060:     patch->usercomputeop(pc, point, x, mat, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell, dofsArrayWithAll ? dofsArrayWithAll + offset*patch->totalDofsPerCell : NULL, patch->usercomputeopctx);
2061:   }
2062:   if (patch->usercomputeopintfacet) {
2063:     PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
2064:     PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
2065:     if (numIntFacets > 0) {
2066:       /* For each interior facet, grab the two cells (in local numbering, and concatenate dof numberings for those cells) */
2067:       PetscInt       *facetDofs = NULL, *facetDofsWithAll = NULL;
2068:       const PetscInt *intFacetsArray = NULL;
2069:       PetscInt        idx = 0;
2070:       PetscInt        i, c, d;
2071:       PetscInt        fStart;
2072:       DM              dm, plex;
2073:       IS              facetIS = NULL;
2074:       const PetscInt *facetCells = NULL;

2076:       ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
2077:       ISGetIndices(patch->intFacets, &intFacetsArray);
2078:       PCGetDM(pc, &dm);
2079:       DMConvert(dm, DMPLEX, &plex);
2080:       dm = plex;
2081:       DMPlexGetHeightStratum(dm, 1, &fStart, NULL);
2082:       /* FIXME: Pull this malloc out. */
2083:       PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofs);
2084:       if (dofsArrayWithAll) {
2085:         PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofsWithAll);
2086:       }
2087:       if (patch->precomputeElementTensors) {
2088:         PetscInt           nFacetDof = 2*patch->totalDofsPerCell;
2089:         const PetscScalar *elementTensors;

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

2093:         for (i = 0; i < numIntFacets; i++) {
2094:           const PetscInt     facet = intFacetsArray[i + intFacetOffset];
2095:           const PetscScalar *v     = elementTensors + patch->precomputedIntFacetTensorLocations[facet - fStart]*nFacetDof*nFacetDof;
2096:           idx = 0;
2097:           /*
2098:            * 0--1
2099:            * |\-|
2100:            * |+\|
2101:            * 2--3
2102:            * [0, 2, 3, 0, 1, 3]
2103:            */
2104:           for (c = 0; c < 2; c++) {
2105:             const PetscInt cell = facetCells[2*(intFacetOffset + i) + c];
2106:             for (d = 0; d < patch->totalDofsPerCell; d++) {
2107:               facetDofs[idx] = dofsArray[(offset + cell)*patch->totalDofsPerCell + d];
2108:               idx++;
2109:             }
2110:           }
2111:           MatSetValues(mat, nFacetDof, facetDofs, nFacetDof, facetDofs, v, ADD_VALUES);
2112:         }
2113:         VecRestoreArrayRead(patch->intFacetMats, &elementTensors);
2114:       } else {
2115:         /*
2116:          * 0--1
2117:          * |\-|
2118:          * |+\|
2119:          * 2--3
2120:          * [0, 2, 3, 0, 1, 3]
2121:          */
2122:         for (i = 0; i < numIntFacets; i++) {
2123:           for (c = 0; c < 2; c++) {
2124:             const PetscInt cell = facetCells[2*(intFacetOffset + i) + c];
2125:             for (d = 0; d < patch->totalDofsPerCell; d++) {
2126:               facetDofs[idx] = dofsArray[(offset + cell)*patch->totalDofsPerCell + d];
2127:               if (dofsArrayWithAll) {
2128:                 facetDofsWithAll[idx] = dofsArrayWithAll[(offset + cell)*patch->totalDofsPerCell + d];
2129:               }
2130:               idx++;
2131:             }
2132:           }
2133:         }
2134:         ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray + intFacetOffset, PETSC_USE_POINTER, &facetIS);
2135:         patch->usercomputeopintfacet(pc, point, x, mat, facetIS, 2*numIntFacets*patch->totalDofsPerCell, facetDofs, facetDofsWithAll, patch->usercomputeopintfacetctx);
2136:         ISDestroy(&facetIS);
2137:       }
2138:       ISRestoreIndices(patch->intFacetsToPatchCell, &facetCells);
2139:       ISRestoreIndices(patch->intFacets, &intFacetsArray);
2140:       PetscFree(facetDofs);
2141:       PetscFree(facetDofsWithAll);
2142:       DMDestroy(&dm);
2143:     }
2144:   }

2146:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
2147:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);

2149:   if (!(withArtificial || isNonlinear) && patch->denseinverse) {
2150:     MatFactorInfo info;
2151:     PetscBool     flg;
2152:     PetscObjectTypeCompare((PetscObject)mat, MATSEQDENSE, &flg);
2154:     MatFactorInfoInitialize(&info);
2155:     MatLUFactor(mat, NULL, NULL, &info);
2156:     MatSeqDenseInvertFactors_Private(mat);
2157:   }
2158:   PetscStackPop;
2159:   ISDestroy(&patch->cellIS);
2160:   if (withArtificial) {
2161:     ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);
2162:   } else {
2163:     ISRestoreIndices(patch->dofs, &dofsArray);
2164:   }
2165:   if (isNonlinear) {
2166:     ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);
2167:   }
2168:   ISRestoreIndices(patch->cells, &cellsArray);
2169:   if (patch->viewMatrix) {
2170:     char name[PETSC_MAX_PATH_LEN];

2172:     PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch matrix for Point %D", point);
2173:     PetscObjectSetName((PetscObject) mat, name);
2174:     ObjectView((PetscObject) mat, patch->viewerMatrix, patch->formatMatrix);
2175:   }
2176:   PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2177:   return 0;
2178: }

2180: static PetscErrorCode MatSetValues_PCPatch_Private(Mat mat, PetscInt m, const PetscInt idxm[],
2181:                                                    PetscInt n, const PetscInt idxn[], const PetscScalar *v, InsertMode addv)
2182: {
2183:   Vec            data;
2184:   PetscScalar   *array;
2185:   PetscInt       bs, nz, i, j, cell;

2187:   MatShellGetContext(mat, &data);
2188:   VecGetBlockSize(data, &bs);
2189:   VecGetSize(data, &nz);
2190:   VecGetArray(data, &array);
2192:   cell = (PetscInt)(idxm[0]/bs); /* use the fact that this is called once per cell */
2193:   for (i = 0; i < m; i++) {
2195:     for (j = 0; j < n; j++) {
2196:       const PetscScalar v_ = v[i*bs + j];
2197:       /* Indexing is special to the data structure we have! */
2198:       if (addv == INSERT_VALUES) {
2199:         array[cell*bs*bs + i*bs + j] = v_;
2200:       } else {
2201:         array[cell*bs*bs + i*bs + j] += v_;
2202:       }
2203:     }
2204:   }
2205:   VecRestoreArray(data, &array);
2206:   return 0;
2207: }

2209: static PetscErrorCode PCPatchPrecomputePatchTensors_Private(PC pc)
2210: {
2211:   PC_PATCH       *patch = (PC_PATCH *)pc->data;
2212:   const PetscInt *cellsArray;
2213:   PetscInt        ncell, offset;
2214:   const PetscInt *dofMapArray;
2215:   PetscInt        i, j;
2216:   IS              dofMap;
2217:   IS              cellIS;
2218:   const PetscInt  ndof  = patch->totalDofsPerCell;
2219:   PetscErrorCode  ierr;
2220:   Mat             vecMat;
2221:   PetscInt        cStart, cEnd;
2222:   DM              dm, plex;

2224:   ISGetSize(patch->cells, &ncell);
2225:   if (!ncell) { /* No cells to assemble over -> skip */
2226:     return 0;
2227:   }

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

2231:   PCGetDM(pc, &dm);
2232:   DMConvert(dm, DMPLEX, &plex);
2233:   dm = plex;
2234:   if (!patch->allCells) {
2235:     PetscHSetI      cells;
2236:     PetscHashIter   hi;
2237:     PetscInt        pStart, pEnd;
2238:     PetscInt        *allCells = NULL;
2239:     PetscHSetICreate(&cells);
2240:     ISGetIndices(patch->cells, &cellsArray);
2241:     PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);
2242:     for (i = pStart; i < pEnd; i++) {
2243:       PetscSectionGetDof(patch->cellCounts, i, &ncell);
2244:       PetscSectionGetOffset(patch->cellCounts, i, &offset);
2245:       if (ncell <= 0) continue;
2246:       for (j = 0; j < ncell; j++) {
2247:         PetscHSetIAdd(cells, cellsArray[offset + j]);
2248:       }
2249:     }
2250:     ISRestoreIndices(patch->cells, &cellsArray);
2251:     PetscHSetIGetSize(cells, &ncell);
2252:     PetscMalloc1(ncell, &allCells);
2253:     DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2254:     PetscMalloc1(cEnd-cStart, &patch->precomputedTensorLocations);
2255:     i = 0;
2256:     PetscHashIterBegin(cells, hi);
2257:     while (!PetscHashIterAtEnd(cells, hi)) {
2258:       PetscHashIterGetKey(cells, hi, allCells[i]);
2259:       patch->precomputedTensorLocations[allCells[i]] = i;
2260:       PetscHashIterNext(cells, hi);
2261:       i++;
2262:     }
2263:     PetscHSetIDestroy(&cells);
2264:     ISCreateGeneral(PETSC_COMM_SELF, ncell, allCells, PETSC_OWN_POINTER, &patch->allCells);
2265:   }
2266:   ISGetSize(patch->allCells, &ncell);
2267:   if (!patch->cellMats) {
2268:     VecCreateSeq(PETSC_COMM_SELF, ncell*ndof*ndof, &patch->cellMats);
2269:     VecSetBlockSize(patch->cellMats, ndof);
2270:   }
2271:   VecSet(patch->cellMats, 0);

2273:   MatCreateShell(PETSC_COMM_SELF, ncell*ndof, ncell*ndof, ncell*ndof, ncell*ndof,
2274:                         (void*)patch->cellMats, &vecMat);
2275:   MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void(*)(void))&MatSetValues_PCPatch_Private);
2276:   ISGetSize(patch->allCells, &ncell);
2277:   ISCreateStride(PETSC_COMM_SELF, ndof*ncell, 0, 1, &dofMap);
2278:   ISGetIndices(dofMap, &dofMapArray);
2279:   ISGetIndices(patch->allCells, &cellsArray);
2280:   ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray, PETSC_USE_POINTER, &cellIS);
2281:   PetscStackPush("PCPatch user callback");
2282:   /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2283:   patch->usercomputeop(pc, -1, NULL, vecMat, cellIS, ndof*ncell, dofMapArray, NULL, patch->usercomputeopctx);
2284:   PetscStackPop;
2285:   ISDestroy(&cellIS);
2286:   MatDestroy(&vecMat);
2287:   ISRestoreIndices(patch->allCells, &cellsArray);
2288:   ISRestoreIndices(dofMap, &dofMapArray);
2289:   ISDestroy(&dofMap);

2291:   if (patch->usercomputeopintfacet) {
2292:     PetscInt nIntFacets;
2293:     IS       intFacetsIS;
2294:     const PetscInt *intFacetsArray = NULL;
2295:     if (!patch->allIntFacets) {
2296:       PetscHSetI      facets;
2297:       PetscHashIter   hi;
2298:       PetscInt        pStart, pEnd, fStart, fEnd;
2299:       PetscInt        *allIntFacets = NULL;
2300:       PetscHSetICreate(&facets);
2301:       ISGetIndices(patch->intFacets, &intFacetsArray);
2302:       PetscSectionGetChart(patch->intFacetCounts, &pStart, &pEnd);
2303:       DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2304:       for (i = pStart; i < pEnd; i++) {
2305:         PetscSectionGetDof(patch->intFacetCounts, i, &nIntFacets);
2306:         PetscSectionGetOffset(patch->intFacetCounts, i, &offset);
2307:         if (nIntFacets <= 0) continue;
2308:         for (j = 0; j < nIntFacets; j++) {
2309:           PetscHSetIAdd(facets, intFacetsArray[offset + j]);
2310:         }
2311:       }
2312:       ISRestoreIndices(patch->intFacets, &intFacetsArray);
2313:       PetscHSetIGetSize(facets, &nIntFacets);
2314:       PetscMalloc1(nIntFacets, &allIntFacets);
2315:       PetscMalloc1(fEnd-fStart, &patch->precomputedIntFacetTensorLocations);
2316:       i = 0;
2317:       PetscHashIterBegin(facets, hi);
2318:       while (!PetscHashIterAtEnd(facets, hi)) {
2319:         PetscHashIterGetKey(facets, hi, allIntFacets[i]);
2320:         patch->precomputedIntFacetTensorLocations[allIntFacets[i] - fStart] = i;
2321:         PetscHashIterNext(facets, hi);
2322:         i++;
2323:       }
2324:       PetscHSetIDestroy(&facets);
2325:       ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, allIntFacets, PETSC_OWN_POINTER, &patch->allIntFacets);
2326:     }
2327:     ISGetSize(patch->allIntFacets, &nIntFacets);
2328:     if (!patch->intFacetMats) {
2329:       VecCreateSeq(PETSC_COMM_SELF, nIntFacets*ndof*ndof*4, &patch->intFacetMats);
2330:       VecSetBlockSize(patch->intFacetMats, ndof*2);
2331:     }
2332:     VecSet(patch->intFacetMats, 0);

2334:     MatCreateShell(PETSC_COMM_SELF, nIntFacets*ndof*2, nIntFacets*ndof*2, nIntFacets*ndof*2, nIntFacets*ndof*2,
2335:                           (void*)patch->intFacetMats, &vecMat);
2336:     MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void(*)(void))&MatSetValues_PCPatch_Private);
2337:     ISCreateStride(PETSC_COMM_SELF, 2*ndof*nIntFacets, 0, 1, &dofMap);
2338:     ISGetIndices(dofMap, &dofMapArray);
2339:     ISGetIndices(patch->allIntFacets, &intFacetsArray);
2340:     ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, intFacetsArray, PETSC_USE_POINTER, &intFacetsIS);
2341:     PetscStackPush("PCPatch user callback (interior facets)");
2342:     /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2343:     patch->usercomputeopintfacet(pc, -1, NULL, vecMat, intFacetsIS, 2*ndof*nIntFacets, dofMapArray, NULL, patch->usercomputeopintfacetctx);
2344:     PetscStackPop;
2345:     ISDestroy(&intFacetsIS);
2346:     MatDestroy(&vecMat);
2347:     ISRestoreIndices(patch->allIntFacets, &intFacetsArray);
2348:     ISRestoreIndices(dofMap, &dofMapArray);
2349:     ISDestroy(&dofMap);
2350:   }
2351:   DMDestroy(&dm);
2352:   PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);

2354:   return 0;
2355: }

2357: PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat, PatchScatterType scattertype)
2358: {
2359:   PC_PATCH          *patch     = (PC_PATCH *) pc->data;
2360:   const PetscScalar *xArray    = NULL;
2361:   PetscScalar       *yArray    = NULL;
2362:   const PetscInt    *gtolArray = NULL;
2363:   PetscInt           dof, offset, lidx;

2366:   VecGetArrayRead(x, &xArray);
2367:   VecGetArray(y, &yArray);
2368:   if (scattertype == SCATTER_WITHARTIFICIAL) {
2369:     PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);
2370:     PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offset);
2371:     ISGetIndices(patch->gtolWithArtificial, &gtolArray);
2372:   } else if (scattertype == SCATTER_WITHALL) {
2373:     PetscSectionGetDof(patch->gtolCountsWithAll, p, &dof);
2374:     PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offset);
2375:     ISGetIndices(patch->gtolWithAll, &gtolArray);
2376:   } else {
2377:     PetscSectionGetDof(patch->gtolCounts, p, &dof);
2378:     PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2379:     ISGetIndices(patch->gtol, &gtolArray);
2380:   }
2383:   for (lidx = 0; lidx < dof; ++lidx) {
2384:     const PetscInt gidx = gtolArray[offset+lidx];

2386:     if (mode == INSERT_VALUES) yArray[lidx]  = xArray[gidx]; /* Forward */
2387:     else                       yArray[gidx] += xArray[lidx]; /* Reverse */
2388:   }
2389:   if (scattertype == SCATTER_WITHARTIFICIAL) {
2390:     ISRestoreIndices(patch->gtolWithArtificial, &gtolArray);
2391:   } else if (scattertype == SCATTER_WITHALL) {
2392:     ISRestoreIndices(patch->gtolWithAll, &gtolArray);
2393:   } else {
2394:     ISRestoreIndices(patch->gtol, &gtolArray);
2395:   }
2396:   VecRestoreArrayRead(x, &xArray);
2397:   VecRestoreArray(y, &yArray);
2398:   return 0;
2399: }

2401: static PetscErrorCode PCSetUp_PATCH_Linear(PC pc)
2402: {
2403:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2404:   const char    *prefix;
2405:   PetscInt       i;

2407:   if (!pc->setupcalled) {
2409:     if (!patch->denseinverse) {
2410:       PetscMalloc1(patch->npatch, &patch->solver);
2411:       PCGetOptionsPrefix(pc, &prefix);
2412:       for (i = 0; i < patch->npatch; ++i) {
2413:         KSP ksp;
2414:         PC  subpc;

2416:         KSPCreate(PETSC_COMM_SELF, &ksp);
2417:         KSPSetErrorIfNotConverged(ksp, pc->erroriffailure);
2418:         KSPSetOptionsPrefix(ksp, prefix);
2419:         KSPAppendOptionsPrefix(ksp, "sub_");
2420:         PetscObjectIncrementTabLevel((PetscObject) ksp, (PetscObject) pc, 1);
2421:         KSPGetPC(ksp, &subpc);
2422:         PetscObjectIncrementTabLevel((PetscObject) subpc, (PetscObject) pc, 1);
2423:         PetscLogObjectParent((PetscObject) pc, (PetscObject) ksp);
2424:         patch->solver[i] = (PetscObject) ksp;
2425:       }
2426:     }
2427:   }
2428:   if (patch->save_operators) {
2429:     if (patch->precomputeElementTensors) {
2430:       PCPatchPrecomputePatchTensors_Private(pc);
2431:     }
2432:     for (i = 0; i < patch->npatch; ++i) {
2433:       PCPatchComputeOperator_Internal(pc, NULL, patch->mat[i], i, PETSC_FALSE);
2434:       if (!patch->denseinverse) {
2435:         KSPSetOperators((KSP) patch->solver[i], patch->mat[i], patch->mat[i]);
2436:       } else if (patch->mat[i] && !patch->densesolve) {
2437:         /* Setup matmult callback */
2438:         MatGetOperation(patch->mat[i], MATOP_MULT, (void (**)(void))&patch->densesolve);
2439:       }
2440:     }
2441:   }
2442:   if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2443:     for (i = 0; i < patch->npatch; ++i) {
2444:       /* Instead of padding patch->patchUpdate with zeros to get */
2445:       /* patch->patchUpdateWithArtificial and then multiplying with the matrix, */
2446:       /* just get rid of the columns that correspond to the dofs with */
2447:       /* artificial bcs. That's of course fairly inefficient, hopefully we */
2448:       /* can just assemble the rectangular matrix in the first place. */
2449:       Mat matSquare;
2450:       IS rowis;
2451:       PetscInt dof;

2453:       MatGetSize(patch->mat[i], &dof, NULL);
2454:       if (dof == 0) {
2455:         patch->matWithArtificial[i] = NULL;
2456:         continue;
2457:       }

2459:       PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);
2460:       PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);

2462:       MatGetSize(matSquare, &dof, NULL);
2463:       ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis);
2464:       if (pc->setupcalled) {
2465:         MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_REUSE_MATRIX, &patch->matWithArtificial[i]);
2466:       } else {
2467:         MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &patch->matWithArtificial[i]);
2468:       }
2469:       ISDestroy(&rowis);
2470:       MatDestroy(&matSquare);
2471:     }
2472:   }
2473:   return 0;
2474: }

2476: static PetscErrorCode PCSetUp_PATCH(PC pc)
2477: {
2478:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2479:   PetscInt       i;
2480:   PetscBool      isNonlinear;
2481:   PetscInt       maxDof = -1, maxDofWithArtificial = -1;

2483:   if (!pc->setupcalled) {
2484:     PetscInt pStart, pEnd, p;
2485:     PetscInt localSize;

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

2489:     isNonlinear = patch->isNonlinear;
2490:     if (!patch->nsubspaces) {
2491:       DM           dm, plex;
2492:       PetscSection s;
2493:       PetscInt     cStart, cEnd, c, Nf, f, numGlobalBcs = 0, *globalBcs, *Nb, **cellDofs;

2495:       PCGetDM(pc, &dm);
2497:       DMConvert(dm, DMPLEX, &plex);
2498:       dm = plex;
2499:       DMGetLocalSection(dm, &s);
2500:       PetscSectionGetNumFields(s, &Nf);
2501:       PetscSectionGetChart(s, &pStart, &pEnd);
2502:       for (p = pStart; p < pEnd; ++p) {
2503:         PetscInt cdof;
2504:         PetscSectionGetConstraintDof(s, p, &cdof);
2505:         numGlobalBcs += cdof;
2506:       }
2507:       DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2508:       PetscMalloc3(Nf, &Nb, Nf, &cellDofs, numGlobalBcs, &globalBcs);
2509:       for (f = 0; f < Nf; ++f) {
2510:         PetscFE        fe;
2511:         PetscDualSpace sp;
2512:         PetscInt       cdoff = 0;

2514:         DMGetField(dm, f, NULL, (PetscObject *) &fe);
2515:         /* PetscFEGetNumComponents(fe, &Nc[f]); */
2516:         PetscFEGetDualSpace(fe, &sp);
2517:         PetscDualSpaceGetDimension(sp, &Nb[f]);

2519:         PetscMalloc1((cEnd-cStart)*Nb[f], &cellDofs[f]);
2520:         for (c = cStart; c < cEnd; ++c) {
2521:           PetscInt *closure = NULL;
2522:           PetscInt  clSize  = 0, cl;

2524:           DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
2525:           for (cl = 0; cl < clSize*2; cl += 2) {
2526:             const PetscInt p = closure[cl];
2527:             PetscInt       fdof, d, foff;

2529:             PetscSectionGetFieldDof(s, p, f, &fdof);
2530:             PetscSectionGetFieldOffset(s, p, f, &foff);
2531:             for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d;
2532:           }
2533:           DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
2534:         }
2536:       }
2537:       numGlobalBcs = 0;
2538:       for (p = pStart; p < pEnd; ++p) {
2539:         const PetscInt *ind;
2540:         PetscInt        off, cdof, d;

2542:         PetscSectionGetOffset(s, p, &off);
2543:         PetscSectionGetConstraintDof(s, p, &cdof);
2544:         PetscSectionGetConstraintIndices(s, p, &ind);
2545:         for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d];
2546:       }

2548:       PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **) cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs);
2549:       for (f = 0; f < Nf; ++f) {
2550:         PetscFree(cellDofs[f]);
2551:       }
2552:       PetscFree3(Nb, cellDofs, globalBcs);
2553:       PCPatchSetComputeFunction(pc, PCPatchComputeFunction_DMPlex_Private, NULL);
2554:       PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL);
2555:       DMDestroy(&dm);
2556:     }

2558:     localSize = patch->subspaceOffsets[patch->nsubspaces];
2559:     VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localRHS);
2560:     VecSetUp(patch->localRHS);
2561:     VecDuplicate(patch->localRHS, &patch->localUpdate);
2562:     PCPatchCreateCellPatches(pc);
2563:     PCPatchCreateCellPatchDiscretisationInfo(pc);

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

2568:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2569:       PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithArtificial);
2570:     }
2571:     if (isNonlinear) {
2572:       PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithAll);
2573:     }
2574:     for (p = pStart; p < pEnd; ++p) {
2575:       PetscInt dof;

2577:       PetscSectionGetDof(patch->gtolCounts, p, &dof);
2578:       maxDof = PetscMax(maxDof, dof);
2579:       if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2580:         const PetscInt    *gtolArray, *gtolArrayWithArtificial = NULL;
2581:         PetscInt           numPatchDofs, offset;
2582:         PetscInt           numPatchDofsWithArtificial, offsetWithArtificial;
2583:         PetscInt           dofWithoutArtificialCounter = 0;
2584:         PetscInt          *patchWithoutArtificialToWithArtificialArray;

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

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

2593:         PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);
2594:         if (numPatchDofs == 0) {
2595:           patch->dofMappingWithoutToWithArtificial[p-pStart] = NULL;
2596:           continue;
2597:         }

2599:         PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2600:         ISGetIndices(patch->gtolWithArtificial, &gtolArrayWithArtificial);
2601:         PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &numPatchDofsWithArtificial);
2602:         PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offsetWithArtificial);

2604:         PetscMalloc1(numPatchDofs, &patchWithoutArtificialToWithArtificialArray);
2605:         for (i=0; i<numPatchDofsWithArtificial; i++) {
2606:           if (gtolArrayWithArtificial[i+offsetWithArtificial] == gtolArray[offset+dofWithoutArtificialCounter]) {
2607:             patchWithoutArtificialToWithArtificialArray[dofWithoutArtificialCounter] = i;
2608:             dofWithoutArtificialCounter++;
2609:             if (dofWithoutArtificialCounter == numPatchDofs)
2610:               break;
2611:           }
2612:         }
2613:         ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutArtificialToWithArtificialArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithArtificial[p-pStart]);
2614:         ISRestoreIndices(patch->gtol, &gtolArray);
2615:         ISRestoreIndices(patch->gtolWithArtificial, &gtolArrayWithArtificial);
2616:       }
2617:     }
2618:     for (p = pStart; p < pEnd; ++p) {
2619:       if (isNonlinear) {
2620:         const PetscInt    *gtolArray, *gtolArrayWithAll = NULL;
2621:         PetscInt           numPatchDofs, offset;
2622:         PetscInt           numPatchDofsWithAll, offsetWithAll;
2623:         PetscInt           dofWithoutAllCounter = 0;
2624:         PetscInt          *patchWithoutAllToWithAllArray;

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

2630:         PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);
2631:         if (numPatchDofs == 0) {
2632:           patch->dofMappingWithoutToWithAll[p-pStart] = NULL;
2633:           continue;
2634:         }

2636:         PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2637:         ISGetIndices(patch->gtolWithAll, &gtolArrayWithAll);
2638:         PetscSectionGetDof(patch->gtolCountsWithAll, p, &numPatchDofsWithAll);
2639:         PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offsetWithAll);

2641:         PetscMalloc1(numPatchDofs, &patchWithoutAllToWithAllArray);

2643:         for (i=0; i<numPatchDofsWithAll; i++) {
2644:           if (gtolArrayWithAll[i+offsetWithAll] == gtolArray[offset+dofWithoutAllCounter]) {
2645:             patchWithoutAllToWithAllArray[dofWithoutAllCounter] = i;
2646:             dofWithoutAllCounter++;
2647:             if (dofWithoutAllCounter == numPatchDofs)
2648:               break;
2649:           }
2650:         }
2651:         ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutAllToWithAllArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithAll[p-pStart]);
2652:         ISRestoreIndices(patch->gtol, &gtolArray);
2653:         ISRestoreIndices(patch->gtolWithAll, &gtolArrayWithAll);
2654:       }
2655:     }
2656:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2657:       VecCreateSeq(PETSC_COMM_SELF, maxDofWithArtificial, &patch->patchRHSWithArtificial);
2658:       VecSetUp(patch->patchRHSWithArtificial);
2659:     }
2660:     VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchRHS);
2661:     VecSetUp(patch->patchRHS);
2662:     VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchUpdate);
2663:     VecSetUp(patch->patchUpdate);
2664:     if (patch->save_operators) {
2665:       PetscMalloc1(patch->npatch, &patch->mat);
2666:       for (i = 0; i < patch->npatch; ++i) {
2667:         PCPatchCreateMatrix_Private(pc, i, &patch->mat[i], PETSC_FALSE);
2668:       }
2669:     }
2670:     PetscLogEventEnd(PC_Patch_CreatePatches, pc, 0, 0, 0);

2672:     /* If desired, calculate weights for dof multiplicity */
2673:     if (patch->partition_of_unity) {
2674:       PetscScalar *input = NULL;
2675:       PetscScalar *output = NULL;
2676:       Vec         global;

2678:       VecDuplicate(patch->localRHS, &patch->dof_weights);
2679:       if (patch->local_composition_type == PC_COMPOSITE_ADDITIVE) {
2680:         for (i = 0; i < patch->npatch; ++i) {
2681:           PetscInt dof;

2683:           PetscSectionGetDof(patch->gtolCounts, i+pStart, &dof);
2684:           if (dof <= 0) continue;
2685:           VecSet(patch->patchRHS, 1.0);
2686:           PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchRHS, patch->dof_weights, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);
2687:         }
2688:       } else {
2689:         /* multiplicative is actually only locally multiplicative and globally additive. need the pou where the mesh decomposition overlaps */
2690:         VecSet(patch->dof_weights, 1.0);
2691:       }

2693:       VecDuplicate(patch->dof_weights, &global);
2694:       VecSet(global, 0.);

2696:       VecGetArray(patch->dof_weights, &input);
2697:       VecGetArray(global, &output);
2698:       PetscSFReduceBegin(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM);
2699:       PetscSFReduceEnd(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM);
2700:       VecRestoreArray(patch->dof_weights, &input);
2701:       VecRestoreArray(global, &output);

2703:       VecReciprocal(global);

2705:       VecGetArray(patch->dof_weights, &output);
2706:       VecGetArray(global, &input);
2707:       PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, input, output,MPI_REPLACE);
2708:       PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, input, output,MPI_REPLACE);
2709:       VecRestoreArray(patch->dof_weights, &output);
2710:       VecRestoreArray(global, &input);
2711:       VecDestroy(&global);
2712:     }
2713:     if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE && patch->save_operators && !patch->isNonlinear) {
2714:       PetscMalloc1(patch->npatch, &patch->matWithArtificial);
2715:     }
2716:   }
2717:   (*patch->setupsolver)(pc);
2718:   return 0;
2719: }

2721: static PetscErrorCode PCApply_PATCH_Linear(PC pc, PetscInt i, Vec x, Vec y)
2722: {
2723:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2724:   KSP            ksp;
2725:   Mat            op;
2726:   PetscInt       m, n;

2728:   if (patch->denseinverse) {
2729:     (*patch->densesolve)(patch->mat[i], x, y);
2730:     return 0;
2731:   }
2732:   ksp = (KSP) patch->solver[i];
2733:   if (!patch->save_operators) {
2734:     Mat mat;

2736:     PCPatchCreateMatrix_Private(pc, i, &mat, PETSC_FALSE);
2737:     /* Populate operator here. */
2738:     PCPatchComputeOperator_Internal(pc, NULL, mat, i, PETSC_FALSE);
2739:     KSPSetOperators(ksp, mat, mat);
2740:     /* Drop reference so the KSPSetOperators below will blow it away. */
2741:     MatDestroy(&mat);
2742:   }
2743:   PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);
2744:   if (!ksp->setfromoptionscalled) {
2745:     KSPSetFromOptions(ksp);
2746:   }
2747:   /* Disgusting trick to reuse work vectors */
2748:   KSPGetOperators(ksp, &op, NULL);
2749:   MatGetLocalSize(op, &m, &n);
2750:   x->map->n = m;
2751:   y->map->n = n;
2752:   x->map->N = m;
2753:   y->map->N = n;
2754:   KSPSolve(ksp, x, y);
2755:   KSPCheckSolve(ksp, pc, y);
2756:   PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);
2757:   if (!patch->save_operators) {
2758:     PC pc;
2759:     KSPSetOperators(ksp, NULL, NULL);
2760:     KSPGetPC(ksp, &pc);
2761:     /* Destroy PC context too, otherwise the factored matrix hangs around. */
2762:     PCReset(pc);
2763:   }
2764:   return 0;
2765: }

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


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

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

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

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

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

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

2882:   ISRestoreIndices(patch->globalBcNodes, &bcNodes);
2883:   VecRestoreArrayRead(x, &globalRHS);
2884:   VecRestoreArray(y, &globalUpdate);

2886:   PetscOptionsPopGetViewerOff();
2887:   PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0);
2888:   return 0;
2889: }

2891: static PetscErrorCode PCReset_PATCH_Linear(PC pc)
2892: {
2893:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
2894:   PetscInt       i;

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

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


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

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

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

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

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

2974:   }
2975:   if (patch->dofMappingWithoutToWithAll) {
2976:     for (i = 0; i < patch->npatch; ++i) ISDestroy(&patch->dofMappingWithoutToWithAll[i]);
2977:     PetscFree(patch->dofMappingWithoutToWithAll);

2979:   }
2980:   PetscFree(patch->sub_mat_type);
2981:   if (patch->userIS) {
2982:     for (i = 0; i < patch->npatch; ++i) ISDestroy(&patch->userIS[i]);
2983:     PetscFree(patch->userIS);
2984:   }
2985:   PetscFree(patch->precomputedTensorLocations);
2986:   PetscFree(patch->precomputedIntFacetTensorLocations);

2988:   patch->bs          = NULL;
2989:   patch->cellNodeMap = NULL;
2990:   patch->nsubspaces  = 0;
2991:   ISDestroy(&patch->iterationSet);

2993:   PetscViewerDestroy(&patch->viewerSection);
2994:   return 0;
2995: }

2997: static PetscErrorCode PCDestroy_PATCH_Linear(PC pc)
2998: {
2999:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
3000:   PetscInt       i;

3002:   if (patch->solver) {
3003:     for (i = 0; i < patch->npatch; ++i) KSPDestroy((KSP *) &patch->solver[i]);
3004:     PetscFree(patch->solver);
3005:   }
3006:   return 0;
3007: }

3009: static PetscErrorCode PCDestroy_PATCH(PC pc)
3010: {
3011:   PC_PATCH      *patch = (PC_PATCH *) pc->data;

3013:   PCReset_PATCH(pc);
3014:   (*patch->destroysolver)(pc);
3015:   PetscFree(pc->data);
3016:   return 0;
3017: }

3019: static PetscErrorCode PCSetFromOptions_PATCH(PetscOptionItems *PetscOptionsObject, PC pc)
3020: {
3021:   PC_PATCH            *patch = (PC_PATCH *) pc->data;
3022:   PCPatchConstructType patchConstructionType = PC_PATCH_STAR;
3023:   char                 sub_mat_type[PETSC_MAX_PATH_LEN];
3024:   char                 option[PETSC_MAX_PATH_LEN];
3025:   const char          *prefix;
3026:   PetscBool            flg, dimflg, codimflg;
3027:   MPI_Comm             comm;
3028:   PetscInt            *ifields, nfields, k;
3029:   PCCompositeType      loctype = PC_COMPOSITE_ADDITIVE;

3031:   PetscObjectGetComm((PetscObject) pc, &comm);
3032:   PetscObjectGetOptionsPrefix((PetscObject) pc, &prefix);
3033:   PetscOptionsHead(PetscOptionsObject, "Patch solver options");

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

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

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

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

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

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

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

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

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

3074:   /* If the user has set the number of subspaces, use that for the buffer size,
3075:      otherwise use a large number */
3076:   if (patch->nsubspaces <= 0) {
3077:     nfields = 128;
3078:   } else {
3079:     nfields = patch->nsubspaces;
3080:   }
3081:   PetscMalloc1(nfields, &ifields);
3082:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exclude_subspaces", patch->classname);
3083:   PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,option,ifields,&nfields,&flg);
3085:   if (flg) {
3086:     PetscHSetIClear(patch->subspaces_to_exclude);
3087:     for (k = 0; k < nfields; k++) {
3088:       PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]);
3089:     }
3090:   }
3091:   PetscFree(ifields);

3093:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_patches_view", patch->classname);
3094:   PetscOptionsBool(option, "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg);
3095:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_cells_view", patch->classname);
3096:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerCells, &patch->formatCells, &patch->viewCells);
3097:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_interior_facets_view", patch->classname);
3098:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerIntFacets, &patch->formatIntFacets, &patch->viewIntFacets);
3099:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exterior_facets_view", patch->classname);
3100:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerExtFacets, &patch->formatExtFacets, &patch->viewExtFacets);
3101:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_points_view", patch->classname);
3102:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints);
3103:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_section_view", patch->classname);
3104:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerSection, &patch->formatSection, &patch->viewSection);
3105:   PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_mat_view", patch->classname);
3106:   PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix);
3107:   PetscOptionsTail();
3108:   patch->optionsSet = PETSC_TRUE;
3109:   return 0;
3110: }

3112: static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc)
3113: {
3114:   PC_PATCH          *patch = (PC_PATCH*) pc->data;
3115:   KSPConvergedReason reason;
3116:   PetscInt           i;

3118:   if (!patch->save_operators) {
3119:     /* Can't do this here because the sub KSPs don't have an operator attached yet. */
3120:     return 0;
3121:   }
3122:   if (patch->denseinverse) {
3123:     /* No solvers */
3124:     return 0;
3125:   }
3126:   for (i = 0; i < patch->npatch; ++i) {
3127:     if (!((KSP) patch->solver[i])->setfromoptionscalled) {
3128:       KSPSetFromOptions((KSP) patch->solver[i]);
3129:     }
3130:     KSPSetUp((KSP) patch->solver[i]);
3131:     KSPGetConvergedReason((KSP) patch->solver[i], &reason);
3132:     if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
3133:   }
3134:   return 0;
3135: }

3137: static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer)
3138: {
3139:   PC_PATCH      *patch = (PC_PATCH *) pc->data;
3140:   PetscViewer    sviewer;
3141:   PetscBool      isascii;
3142:   PetscMPIInt    rank;

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

3168:   if (patch->denseinverse) {
3169:     PetscViewerASCIIPrintf(viewer, "Explicitly forming dense inverse and applying patch solver via MatMult.\n");
3170:   } else {
3171:     if (patch->isNonlinear) {
3172:       PetscViewerASCIIPrintf(viewer, "SNES on patches (all same):\n");
3173:     } else {
3174:       PetscViewerASCIIPrintf(viewer, "KSP on patches (all same):\n");
3175:     }
3176:     if (patch->solver) {
3177:       PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
3178:       if (rank == 0) {
3179:         PetscViewerASCIIPushTab(sviewer);
3180:         PetscObjectView(patch->solver[0], sviewer);
3181:         PetscViewerASCIIPopTab(sviewer);
3182:       }
3183:       PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
3184:     } else {
3185:       PetscViewerASCIIPushTab(viewer);
3186:       PetscViewerASCIIPrintf(viewer, "Solver not yet set.\n");
3187:       PetscViewerASCIIPopTab(viewer);
3188:     }
3189:   }
3190:   PetscViewerASCIIPopTab(viewer);
3191:   return 0;
3192: }

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

3199:   Options Database Keys:
3200: + -pc_patch_cells_view   - Views the process local cell numbers for each patch
3201: . -pc_patch_points_view  - Views the process local mesh point numbers for each patch
3202: . -pc_patch_g2l_view     - Views the map between global dofs and patch local dofs for each patch
3203: . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary
3204: - -pc_patch_sub_mat_view - Views the matrix associated with each patch

3206:   Level: intermediate

3208: .seealso: PCType, PCCreate(), PCSetType()
3209: M*/
3210: PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc)
3211: {
3212:   PC_PATCH      *patch;

3214:   PetscNewLog(pc, &patch);

3216:   if (patch->subspaces_to_exclude) {
3217:     PetscHSetIDestroy(&patch->subspaces_to_exclude);
3218:   }
3219:   PetscHSetICreate(&patch->subspaces_to_exclude);

3221:   patch->classname = "pc";
3222:   patch->isNonlinear = PETSC_FALSE;

3224:   /* Set some defaults */
3225:   patch->combined           = PETSC_FALSE;
3226:   patch->save_operators     = PETSC_TRUE;
3227:   patch->local_composition_type = PC_COMPOSITE_ADDITIVE;
3228:   patch->precomputeElementTensors = PETSC_FALSE;
3229:   patch->partition_of_unity = PETSC_FALSE;
3230:   patch->codim              = -1;
3231:   patch->dim                = -1;
3232:   patch->vankadim           = -1;
3233:   patch->ignoredim          = -1;
3234:   patch->pardecomp_overlap  = 0;
3235:   patch->patchconstructop   = PCPatchConstruct_Star;
3236:   patch->symmetrise_sweep   = PETSC_FALSE;
3237:   patch->npatch             = 0;
3238:   patch->userIS             = NULL;
3239:   patch->optionsSet         = PETSC_FALSE;
3240:   patch->iterationSet       = NULL;
3241:   patch->user_patches       = PETSC_FALSE;
3242:   PetscStrallocpy(MATDENSE, (char **) &patch->sub_mat_type);
3243:   patch->viewPatches        = PETSC_FALSE;
3244:   patch->viewCells          = PETSC_FALSE;
3245:   patch->viewPoints         = PETSC_FALSE;
3246:   patch->viewSection        = PETSC_FALSE;
3247:   patch->viewMatrix         = PETSC_FALSE;
3248:   patch->densesolve         = NULL;
3249:   patch->setupsolver        = PCSetUp_PATCH_Linear;
3250:   patch->applysolver        = PCApply_PATCH_Linear;
3251:   patch->resetsolver        = PCReset_PATCH_Linear;
3252:   patch->destroysolver      = PCDestroy_PATCH_Linear;
3253:   patch->updatemultiplicative = PCUpdateMultiplicative_PATCH_Linear;
3254:   patch->dofMappingWithoutToWithArtificial = NULL;
3255:   patch->dofMappingWithoutToWithAll = NULL;

3257:   pc->data                 = (void *) patch;
3258:   pc->ops->apply           = PCApply_PATCH;
3259:   pc->ops->applytranspose  = NULL; /* PCApplyTranspose_PATCH; */
3260:   pc->ops->setup           = PCSetUp_PATCH;
3261:   pc->ops->reset           = PCReset_PATCH;
3262:   pc->ops->destroy         = PCDestroy_PATCH;
3263:   pc->ops->setfromoptions  = PCSetFromOptions_PATCH;
3264:   pc->ops->setuponblocks   = PCSetUpOnBlocks_PATCH;
3265:   pc->ops->view            = PCView_PATCH;
3266:   pc->ops->applyrichardson = NULL;

3268:   return 0;
3269: }