Actual source code: pcpatch.c
1: #include <petsc/private/pcpatchimpl.h>
2: #include <petsc/private/kspimpl.h>
3: #include <petsc/private/vecimpl.h>
4: #include <petsc/private/dmpleximpl.h>
5: #include <petscsf.h>
6: #include <petscbt.h>
7: #include <petscds.h>
8: #include <../src/mat/impls/dense/seq/dense.h>
10: PetscLogEvent PC_Patch_CreatePatches, PC_Patch_ComputeOp, PC_Patch_Solve, PC_Patch_Apply, PC_Patch_Prealloc;
12: PETSC_STATIC_INLINE PetscErrorCode ObjectView(PetscObject obj, PetscViewer viewer, PetscViewerFormat format)
13: {
16: PetscViewerPushFormat(viewer, format);
17: PetscObjectView(obj, viewer);
18: PetscViewerPopFormat(viewer);
19: return(0);
20: }
22: static PetscErrorCode PCPatchConstruct_Star(void *vpatch, DM dm, PetscInt point, PetscHSetI ht)
23: {
24: PetscInt starSize;
25: PetscInt *star = NULL, si;
29: PetscHSetIClear(ht);
30: /* To start with, add the point we care about */
31: PetscHSetIAdd(ht, point);
32: /* Loop over all the points that this point connects to */
33: DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);
34: for (si = 0; si < starSize*2; si += 2) {PetscHSetIAdd(ht, star[si]);}
35: DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);
36: return(0);
37: }
39: static PetscErrorCode PCPatchConstruct_Vanka(void *vpatch, DM dm, PetscInt point, PetscHSetI ht)
40: {
41: PC_PATCH *patch = (PC_PATCH *) vpatch;
42: PetscInt starSize;
43: PetscInt *star = NULL;
44: PetscBool shouldIgnore = PETSC_FALSE;
45: PetscInt cStart, cEnd, iStart, iEnd, si;
49: PetscHSetIClear(ht);
50: /* To start with, add the point we care about */
51: PetscHSetIAdd(ht, point);
52: /* Should we ignore any points of a certain dimension? */
53: if (patch->vankadim >= 0) {
54: shouldIgnore = PETSC_TRUE;
55: DMPlexGetDepthStratum(dm, patch->vankadim, &iStart, &iEnd);
56: }
57: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
58: /* Loop over all the cells that this point connects to */
59: DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);
60: for (si = 0; si < starSize*2; si += 2) {
61: const PetscInt cell = star[si];
62: PetscInt closureSize;
63: PetscInt *closure = NULL, ci;
65: if (cell < cStart || cell >= cEnd) continue;
66: /* now loop over all entities in the closure of that cell */
67: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
68: for (ci = 0; ci < closureSize*2; ci += 2) {
69: const PetscInt newpoint = closure[ci];
71: /* We've been told to ignore entities of this type.*/
72: if (shouldIgnore && newpoint >= iStart && newpoint < iEnd) continue;
73: PetscHSetIAdd(ht, newpoint);
74: }
75: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
76: }
77: DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);
78: return(0);
79: }
81: static PetscErrorCode PCPatchConstruct_Pardecomp(void *vpatch, DM dm, PetscInt point, PetscHSetI ht)
82: {
83: PC_PATCH *patch = (PC_PATCH *) vpatch;
84: DMLabel ghost = NULL;
85: const PetscInt *leaves;
86: PetscInt nleaves, pStart, pEnd, loc;
87: PetscBool isFiredrake;
88: PetscBool flg;
89: PetscInt starSize;
90: PetscInt *star = NULL;
91: PetscInt opoint, overlapi;
92: PetscErrorCode ierr;
95: PetscHSetIClear(ht);
97: DMPlexGetChart(dm, &pStart, &pEnd);
99: DMHasLabel(dm, "pyop2_ghost", &isFiredrake);
100: if (isFiredrake) {
101: DMGetLabel(dm, "pyop2_ghost", &ghost);
102: DMLabelCreateIndex(ghost, pStart, pEnd);
103: } else {
104: PetscSF sf;
105: DMGetPointSF(dm, &sf);
106: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
107: nleaves = PetscMax(nleaves, 0);
108: }
110: for (opoint = pStart; opoint < pEnd; ++opoint) {
111: if (ghost) {DMLabelHasPoint(ghost, opoint, &flg);}
112: else {PetscFindInt(opoint, nleaves, leaves, &loc); flg = loc >=0 ? PETSC_TRUE : PETSC_FALSE;}
113: /* Not an owned entity, don't make a cell patch. */
114: if (flg) continue;
115: PetscHSetIAdd(ht, opoint);
116: }
118: /* Now build the overlap for the patch */
119: for (overlapi = 0; overlapi < patch->pardecomp_overlap; ++overlapi) {
120: PetscInt index = 0;
121: PetscInt *htpoints = NULL;
122: PetscInt htsize;
123: PetscInt i;
125: PetscHSetIGetSize(ht, &htsize);
126: PetscMalloc1(htsize, &htpoints);
127: PetscHSetIGetElems(ht, &index, htpoints);
129: for (i = 0; i < htsize; ++i) {
130: PetscInt hpoint = htpoints[i];
131: PetscInt si;
133: DMPlexGetTransitiveClosure(dm, hpoint, PETSC_FALSE, &starSize, &star);
134: for (si = 0; si < starSize*2; si += 2) {
135: const PetscInt starp = star[si];
136: PetscInt closureSize;
137: PetscInt *closure = NULL, ci;
139: /* now loop over all entities in the closure of starp */
140: DMPlexGetTransitiveClosure(dm, starp, PETSC_TRUE, &closureSize, &closure);
141: for (ci = 0; ci < closureSize*2; ci += 2) {
142: const PetscInt closstarp = closure[ci];
143: PetscHSetIAdd(ht, closstarp);
144: }
145: DMPlexRestoreTransitiveClosure(dm, starp, PETSC_TRUE, &closureSize, &closure);
146: }
147: DMPlexRestoreTransitiveClosure(dm, hpoint, PETSC_FALSE, &starSize, &star);
148: }
149: PetscFree(htpoints);
150: }
152: return(0);
153: }
155: /* The user's already set the patches in patch->userIS. Build the hash tables */
156: static PetscErrorCode PCPatchConstruct_User(void *vpatch, DM dm, PetscInt point, PetscHSetI ht)
157: {
158: PC_PATCH *patch = (PC_PATCH *) vpatch;
159: IS patchis = patch->userIS[point];
160: PetscInt n;
161: const PetscInt *patchdata;
162: PetscInt pStart, pEnd, i;
163: PetscErrorCode ierr;
166: PetscHSetIClear(ht);
167: DMPlexGetChart(dm, &pStart, &pEnd);
168: ISGetLocalSize(patchis, &n);
169: ISGetIndices(patchis, &patchdata);
170: for (i = 0; i < n; ++i) {
171: const PetscInt ownedpoint = patchdata[i];
173: if (ownedpoint < pStart || ownedpoint >= pEnd) {
174: SETERRQ3(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D was not in [%D, %D)", ownedpoint, pStart, pEnd);
175: }
176: PetscHSetIAdd(ht, ownedpoint);
177: }
178: ISRestoreIndices(patchis, &patchdata);
179: return(0);
180: }
182: static PetscErrorCode PCPatchCreateDefaultSF_Private(PC pc, PetscInt n, const PetscSF *sf, const PetscInt *bs)
183: {
184: PC_PATCH *patch = (PC_PATCH *) pc->data;
185: PetscInt i;
189: if (n == 1 && bs[0] == 1) {
190: patch->sectionSF = sf[0];
191: PetscObjectReference((PetscObject) patch->sectionSF);
192: } else {
193: PetscInt allRoots = 0, allLeaves = 0;
194: PetscInt leafOffset = 0;
195: PetscInt *ilocal = NULL;
196: PetscSFNode *iremote = NULL;
197: PetscInt *remoteOffsets = NULL;
198: PetscInt index = 0;
199: PetscHMapI rankToIndex;
200: PetscInt numRanks = 0;
201: PetscSFNode *remote = NULL;
202: PetscSF rankSF;
203: PetscInt *ranks = NULL;
204: PetscInt *offsets = NULL;
205: MPI_Datatype contig;
206: PetscHSetI ranksUniq;
208: /* First figure out how many dofs there are in the concatenated numbering.
209: * allRoots: number of owned global dofs;
210: * allLeaves: number of visible dofs (global + ghosted).
211: */
212: for (i = 0; i < n; ++i) {
213: PetscInt nroots, nleaves;
215: PetscSFGetGraph(sf[i], &nroots, &nleaves, NULL, NULL);
216: allRoots += nroots * bs[i];
217: allLeaves += nleaves * bs[i];
218: }
219: PetscMalloc1(allLeaves, &ilocal);
220: PetscMalloc1(allLeaves, &iremote);
221: /* Now build an SF that just contains process connectivity. */
222: PetscHSetICreate(&ranksUniq);
223: for (i = 0; i < n; ++i) {
224: const PetscMPIInt *ranks = NULL;
225: PetscInt nranks, j;
227: PetscSFSetUp(sf[i]);
228: PetscSFGetRootRanks(sf[i], &nranks, &ranks, NULL, NULL, NULL);
229: /* These are all the ranks who communicate with me. */
230: for (j = 0; j < nranks; ++j) {
231: PetscHSetIAdd(ranksUniq, (PetscInt) ranks[j]);
232: }
233: }
234: PetscHSetIGetSize(ranksUniq, &numRanks);
235: PetscMalloc1(numRanks, &remote);
236: PetscMalloc1(numRanks, &ranks);
237: PetscHSetIGetElems(ranksUniq, &index, ranks);
239: PetscHMapICreate(&rankToIndex);
240: for (i = 0; i < numRanks; ++i) {
241: remote[i].rank = ranks[i];
242: remote[i].index = 0;
243: PetscHMapISet(rankToIndex, ranks[i], i);
244: }
245: PetscFree(ranks);
246: PetscHSetIDestroy(&ranksUniq);
247: PetscSFCreate(PetscObjectComm((PetscObject) pc), &rankSF);
248: PetscSFSetGraph(rankSF, 1, numRanks, NULL, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
249: PetscSFSetUp(rankSF);
250: /* OK, use it to communicate the root offset on the remote
251: * processes for each subspace. */
252: PetscMalloc1(n, &offsets);
253: PetscMalloc1(n*numRanks, &remoteOffsets);
255: offsets[0] = 0;
256: for (i = 1; i < n; ++i) {
257: PetscInt nroots;
259: PetscSFGetGraph(sf[i-1], &nroots, NULL, NULL, NULL);
260: offsets[i] = offsets[i-1] + nroots*bs[i-1];
261: }
262: /* Offsets are the offsets on the current process of the
263: * global dof numbering for the subspaces. */
264: MPI_Type_contiguous(n, MPIU_INT, &contig);
265: MPI_Type_commit(&contig);
267: PetscSFBcastBegin(rankSF, contig, offsets, remoteOffsets,MPI_REPLACE);
268: PetscSFBcastEnd(rankSF, contig, offsets, remoteOffsets,MPI_REPLACE);
269: MPI_Type_free(&contig);
270: PetscFree(offsets);
271: PetscSFDestroy(&rankSF);
272: /* Now remoteOffsets contains the offsets on the remote
273: * processes who communicate with me. So now we can
274: * concatenate the list of SFs into a single one. */
275: index = 0;
276: for (i = 0; i < n; ++i) {
277: const PetscSFNode *remote = NULL;
278: const PetscInt *local = NULL;
279: PetscInt nroots, nleaves, j;
281: PetscSFGetGraph(sf[i], &nroots, &nleaves, &local, &remote);
282: for (j = 0; j < nleaves; ++j) {
283: PetscInt rank = remote[j].rank;
284: PetscInt idx, rootOffset, k;
286: PetscHMapIGet(rankToIndex, rank, &idx);
287: if (idx == -1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Didn't find rank, huh?");
288: /* Offset on given rank for ith subspace */
289: rootOffset = remoteOffsets[n*idx + i];
290: for (k = 0; k < bs[i]; ++k) {
291: ilocal[index] = (local ? local[j] : j)*bs[i] + k + leafOffset;
292: iremote[index].rank = remote[j].rank;
293: iremote[index].index = remote[j].index*bs[i] + k + rootOffset;
294: ++index;
295: }
296: }
297: leafOffset += nleaves * bs[i];
298: }
299: PetscHMapIDestroy(&rankToIndex);
300: PetscFree(remoteOffsets);
301: PetscSFCreate(PetscObjectComm((PetscObject)pc), &patch->sectionSF);
302: PetscSFSetGraph(patch->sectionSF, allRoots, allLeaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
303: }
304: return(0);
305: }
307: PetscErrorCode PCPatchSetDenseInverse(PC pc, PetscBool flg)
308: {
309: PC_PATCH *patch = (PC_PATCH *) pc->data;
311: patch->denseinverse = flg;
312: return(0);
313: }
315: PetscErrorCode PCPatchGetDenseInverse(PC pc, PetscBool *flg)
316: {
317: PC_PATCH *patch = (PC_PATCH *) pc->data;
319: *flg = patch->denseinverse;
320: return(0);
321: }
323: /* TODO: Docs */
324: PetscErrorCode PCPatchSetIgnoreDim(PC pc, PetscInt dim)
325: {
326: PC_PATCH *patch = (PC_PATCH *) pc->data;
328: patch->ignoredim = dim;
329: return(0);
330: }
332: /* TODO: Docs */
333: PetscErrorCode PCPatchGetIgnoreDim(PC pc, PetscInt *dim)
334: {
335: PC_PATCH *patch = (PC_PATCH *) pc->data;
337: *dim = patch->ignoredim;
338: return(0);
339: }
341: /* TODO: Docs */
342: PetscErrorCode PCPatchSetSaveOperators(PC pc, PetscBool flg)
343: {
344: PC_PATCH *patch = (PC_PATCH *) pc->data;
346: patch->save_operators = flg;
347: return(0);
348: }
350: /* TODO: Docs */
351: PetscErrorCode PCPatchGetSaveOperators(PC pc, PetscBool *flg)
352: {
353: PC_PATCH *patch = (PC_PATCH *) pc->data;
355: *flg = patch->save_operators;
356: return(0);
357: }
359: /* TODO: Docs */
360: PetscErrorCode PCPatchSetPrecomputeElementTensors(PC pc, PetscBool flg)
361: {
362: PC_PATCH *patch = (PC_PATCH *) pc->data;
364: patch->precomputeElementTensors = flg;
365: return(0);
366: }
368: /* TODO: Docs */
369: PetscErrorCode PCPatchGetPrecomputeElementTensors(PC pc, PetscBool *flg)
370: {
371: PC_PATCH *patch = (PC_PATCH *) pc->data;
373: *flg = patch->precomputeElementTensors;
374: return(0);
375: }
377: /* TODO: Docs */
378: PetscErrorCode PCPatchSetPartitionOfUnity(PC pc, PetscBool flg)
379: {
380: PC_PATCH *patch = (PC_PATCH *) pc->data;
382: patch->partition_of_unity = flg;
383: return(0);
384: }
386: /* TODO: Docs */
387: PetscErrorCode PCPatchGetPartitionOfUnity(PC pc, PetscBool *flg)
388: {
389: PC_PATCH *patch = (PC_PATCH *) pc->data;
391: *flg = patch->partition_of_unity;
392: return(0);
393: }
395: /* TODO: Docs */
396: PetscErrorCode PCPatchSetLocalComposition(PC pc, PCCompositeType type)
397: {
398: PC_PATCH *patch = (PC_PATCH *) pc->data;
400: if (type != PC_COMPOSITE_ADDITIVE && type != PC_COMPOSITE_MULTIPLICATIVE) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Only supports additive or multiplicative as the local type");
401: patch->local_composition_type = type;
402: return(0);
403: }
405: /* TODO: Docs */
406: PetscErrorCode PCPatchGetLocalComposition(PC pc, PCCompositeType *type)
407: {
408: PC_PATCH *patch = (PC_PATCH *) pc->data;
410: *type = patch->local_composition_type;
411: return(0);
412: }
414: /* TODO: Docs */
415: PetscErrorCode PCPatchSetSubMatType(PC pc, MatType sub_mat_type)
416: {
417: PC_PATCH *patch = (PC_PATCH *) pc->data;
421: if (patch->sub_mat_type) {PetscFree(patch->sub_mat_type);}
422: PetscStrallocpy(sub_mat_type, (char **) &patch->sub_mat_type);
423: return(0);
424: }
426: /* TODO: Docs */
427: PetscErrorCode PCPatchGetSubMatType(PC pc, MatType *sub_mat_type)
428: {
429: PC_PATCH *patch = (PC_PATCH *) pc->data;
431: *sub_mat_type = patch->sub_mat_type;
432: return(0);
433: }
435: /* TODO: Docs */
436: PetscErrorCode PCPatchSetCellNumbering(PC pc, PetscSection cellNumbering)
437: {
438: PC_PATCH *patch = (PC_PATCH *) pc->data;
442: patch->cellNumbering = cellNumbering;
443: PetscObjectReference((PetscObject) cellNumbering);
444: return(0);
445: }
447: /* TODO: Docs */
448: PetscErrorCode PCPatchGetCellNumbering(PC pc, PetscSection *cellNumbering)
449: {
450: PC_PATCH *patch = (PC_PATCH *) pc->data;
452: *cellNumbering = patch->cellNumbering;
453: return(0);
454: }
456: /* TODO: Docs */
457: PetscErrorCode PCPatchSetConstructType(PC pc, PCPatchConstructType ctype, PetscErrorCode (*func)(PC, PetscInt *, IS **, IS *, void *), void *ctx)
458: {
459: PC_PATCH *patch = (PC_PATCH *) pc->data;
462: patch->ctype = ctype;
463: switch (ctype) {
464: case PC_PATCH_STAR:
465: patch->user_patches = PETSC_FALSE;
466: patch->patchconstructop = PCPatchConstruct_Star;
467: break;
468: case PC_PATCH_VANKA:
469: patch->user_patches = PETSC_FALSE;
470: patch->patchconstructop = PCPatchConstruct_Vanka;
471: break;
472: case PC_PATCH_PARDECOMP:
473: patch->user_patches = PETSC_FALSE;
474: patch->patchconstructop = PCPatchConstruct_Pardecomp;
475: break;
476: case PC_PATCH_USER:
477: case PC_PATCH_PYTHON:
478: patch->user_patches = PETSC_TRUE;
479: patch->patchconstructop = PCPatchConstruct_User;
480: if (func) {
481: patch->userpatchconstructionop = func;
482: patch->userpatchconstructctx = ctx;
483: }
484: break;
485: default:
486: SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_USER, "Unknown patch construction type %D", (PetscInt) patch->ctype);
487: }
488: return(0);
489: }
491: /* TODO: Docs */
492: PetscErrorCode PCPatchGetConstructType(PC pc, PCPatchConstructType *ctype, PetscErrorCode (**func)(PC, PetscInt *, IS **, IS *, void *), void **ctx)
493: {
494: PC_PATCH *patch = (PC_PATCH *) pc->data;
497: *ctype = patch->ctype;
498: switch (patch->ctype) {
499: case PC_PATCH_STAR:
500: case PC_PATCH_VANKA:
501: case PC_PATCH_PARDECOMP:
502: break;
503: case PC_PATCH_USER:
504: case PC_PATCH_PYTHON:
505: *func = patch->userpatchconstructionop;
506: *ctx = patch->userpatchconstructctx;
507: break;
508: default:
509: SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_USER, "Unknown patch construction type %D", (PetscInt) patch->ctype);
510: }
511: return(0);
512: }
514: /* TODO: Docs */
515: PetscErrorCode PCPatchSetDiscretisationInfo(PC pc, PetscInt nsubspaces, DM *dms, PetscInt *bs, PetscInt *nodesPerCell, const PetscInt **cellNodeMap,
516: const PetscInt *subspaceOffsets, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes)
517: {
518: PC_PATCH *patch = (PC_PATCH *) pc->data;
519: DM dm, plex;
520: PetscSF *sfs;
521: PetscInt cStart, cEnd, i, j;
525: PCGetDM(pc, &dm);
526: DMConvert(dm, DMPLEX, &plex);
527: dm = plex;
528: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
529: PetscMalloc1(nsubspaces, &sfs);
530: PetscMalloc1(nsubspaces, &patch->dofSection);
531: PetscMalloc1(nsubspaces, &patch->bs);
532: PetscMalloc1(nsubspaces, &patch->nodesPerCell);
533: PetscMalloc1(nsubspaces, &patch->cellNodeMap);
534: PetscMalloc1(nsubspaces+1, &patch->subspaceOffsets);
536: patch->nsubspaces = nsubspaces;
537: patch->totalDofsPerCell = 0;
538: for (i = 0; i < nsubspaces; ++i) {
539: DMGetLocalSection(dms[i], &patch->dofSection[i]);
540: PetscObjectReference((PetscObject) patch->dofSection[i]);
541: DMGetSectionSF(dms[i], &sfs[i]);
542: patch->bs[i] = bs[i];
543: patch->nodesPerCell[i] = nodesPerCell[i];
544: patch->totalDofsPerCell += nodesPerCell[i]*bs[i];
545: PetscMalloc1((cEnd-cStart)*nodesPerCell[i], &patch->cellNodeMap[i]);
546: for (j = 0; j < (cEnd-cStart)*nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j];
547: patch->subspaceOffsets[i] = subspaceOffsets[i];
548: }
549: PCPatchCreateDefaultSF_Private(pc, nsubspaces, sfs, patch->bs);
550: PetscFree(sfs);
552: patch->subspaceOffsets[nsubspaces] = subspaceOffsets[nsubspaces];
553: ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);
554: ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);
555: DMDestroy(&dm);
556: return(0);
557: }
559: /* TODO: Docs */
560: PetscErrorCode PCPatchSetDiscretisationInfoCombined(PC pc, DM dm, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes)
561: {
562: PC_PATCH *patch = (PC_PATCH *) pc->data;
563: PetscInt cStart, cEnd, i, j;
567: patch->combined = PETSC_TRUE;
568: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
569: DMGetNumFields(dm, &patch->nsubspaces);
570: PetscCalloc1(patch->nsubspaces, &patch->dofSection);
571: PetscMalloc1(patch->nsubspaces, &patch->bs);
572: PetscMalloc1(patch->nsubspaces, &patch->nodesPerCell);
573: PetscMalloc1(patch->nsubspaces, &patch->cellNodeMap);
574: PetscCalloc1(patch->nsubspaces+1, &patch->subspaceOffsets);
575: DMGetLocalSection(dm, &patch->dofSection[0]);
576: PetscObjectReference((PetscObject) patch->dofSection[0]);
577: PetscSectionGetStorageSize(patch->dofSection[0], &patch->subspaceOffsets[patch->nsubspaces]);
578: patch->totalDofsPerCell = 0;
579: for (i = 0; i < patch->nsubspaces; ++i) {
580: patch->bs[i] = 1;
581: patch->nodesPerCell[i] = nodesPerCell[i];
582: patch->totalDofsPerCell += nodesPerCell[i];
583: PetscMalloc1((cEnd-cStart)*nodesPerCell[i], &patch->cellNodeMap[i]);
584: for (j = 0; j < (cEnd-cStart)*nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j];
585: }
586: DMGetSectionSF(dm, &patch->sectionSF);
587: PetscObjectReference((PetscObject) patch->sectionSF);
588: ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);
589: ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);
590: return(0);
591: }
593: /*@C
595: PCPatchSetComputeFunction - Set the callback used to compute patch residuals
597: Logically collective on PC
599: Input Parameters:
600: + pc - The PC
601: . func - The callback
602: - ctx - The user context
604: Calling sequence of func:
605: $ func (PC pc,PetscInt point,Vec x,Vec f,IS cellIS,PetscInt n,const PetscInt* dofsArray,const PetscInt* dofsArrayWithAll,void* ctx)
607: + pc - The PC
608: . point - The point
609: . x - The input solution (not used in linear problems)
610: . f - The patch residual vector
611: . cellIS - An array of the cell numbers
612: . n - The size of dofsArray
613: . dofsArray - The dofmap for the dofs to be solved for
614: . dofsArrayWithAll - The dofmap for all dofs on the patch
615: - ctx - The user context
617: Level: advanced
619: Notes:
620: The entries of F (the output residual vector) have been set to zero before the call.
622: .seealso: PCPatchSetComputeOperator(), PCPatchGetComputeOperator(), PCPatchSetDiscretisationInfo(), PCPatchSetComputeFunctionInteriorFacets()
623: @*/
624: PetscErrorCode PCPatchSetComputeFunction(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Vec, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx)
625: {
626: PC_PATCH *patch = (PC_PATCH *) pc->data;
629: patch->usercomputef = func;
630: patch->usercomputefctx = ctx;
631: return(0);
632: }
634: /*@C
636: PCPatchSetComputeFunctionInteriorFacets - Set the callback used to compute facet integrals for patch residuals
638: Logically collective on PC
640: Input Parameters:
641: + pc - The PC
642: . func - The callback
643: - ctx - The user context
645: Calling sequence of func:
646: $ func (PC pc,PetscInt point,Vec x,Vec f,IS facetIS,PetscInt n,const PetscInt* dofsArray,const PetscInt* dofsArrayWithAll,void* ctx)
648: + pc - The PC
649: . point - The point
650: . x - The input solution (not used in linear problems)
651: . f - The patch residual vector
652: . facetIS - An array of the facet numbers
653: . n - The size of dofsArray
654: . dofsArray - The dofmap for the dofs to be solved for
655: . dofsArrayWithAll - The dofmap for all dofs on the patch
656: - ctx - The user context
658: Level: advanced
660: Notes:
661: The entries of F (the output residual vector) have been set to zero before the call.
663: .seealso: PCPatchSetComputeOperator(), PCPatchGetComputeOperator(), PCPatchSetDiscretisationInfo(), PCPatchSetComputeFunction()
664: @*/
665: PetscErrorCode PCPatchSetComputeFunctionInteriorFacets(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Vec, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx)
666: {
667: PC_PATCH *patch = (PC_PATCH *) pc->data;
670: patch->usercomputefintfacet = func;
671: patch->usercomputefintfacetctx = ctx;
672: return(0);
673: }
675: /*@C
677: PCPatchSetComputeOperator - Set the callback used to compute patch matrices
679: Logically collective on PC
681: Input Parameters:
682: + pc - The PC
683: . func - The callback
684: - ctx - The user context
686: Calling sequence of func:
687: $ func (PC pc,PetscInt point,Vec x,Mat mat,IS facetIS,PetscInt n,const PetscInt* dofsArray,const PetscInt* dofsArrayWithAll,void* ctx)
689: + pc - The PC
690: . point - The point
691: . x - The input solution (not used in linear problems)
692: . mat - The patch matrix
693: . cellIS - An array of the cell numbers
694: . n - The size of dofsArray
695: . dofsArray - The dofmap for the dofs to be solved for
696: . dofsArrayWithAll - The dofmap for all dofs on the patch
697: - ctx - The user context
699: Level: advanced
701: Notes:
702: The matrix entries have been set to zero before the call.
704: .seealso: PCPatchGetComputeOperator(), PCPatchSetComputeFunction(), PCPatchSetDiscretisationInfo()
705: @*/
706: PetscErrorCode PCPatchSetComputeOperator(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Mat, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx)
707: {
708: PC_PATCH *patch = (PC_PATCH *) pc->data;
711: patch->usercomputeop = func;
712: patch->usercomputeopctx = ctx;
713: return(0);
714: }
716: /*@C
718: PCPatchSetComputeOperatorInteriorFacets - Set the callback used to compute facet integrals for patch matrices
720: Logically collective on PC
722: Input Parameters:
723: + pc - The PC
724: . func - The callback
725: - ctx - The user context
727: Calling sequence of func:
728: $ func (PC pc,PetscInt point,Vec x,Mat mat,IS facetIS,PetscInt n,const PetscInt* dofsArray,const PetscInt* dofsArrayWithAll,void* ctx)
730: + pc - The PC
731: . point - The point
732: . x - The input solution (not used in linear problems)
733: . mat - The patch matrix
734: . facetIS - An array of the facet numbers
735: . n - The size of dofsArray
736: . dofsArray - The dofmap for the dofs to be solved for
737: . dofsArrayWithAll - The dofmap for all dofs on the patch
738: - ctx - The user context
740: Level: advanced
742: Notes:
743: The matrix entries have been set to zero before the call.
745: .seealso: PCPatchGetComputeOperator(), PCPatchSetComputeFunction(), PCPatchSetDiscretisationInfo()
746: @*/
747: PetscErrorCode PCPatchSetComputeOperatorInteriorFacets(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Mat, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx)
748: {
749: PC_PATCH *patch = (PC_PATCH *) pc->data;
752: patch->usercomputeopintfacet = func;
753: patch->usercomputeopintfacetctx = ctx;
754: return(0);
755: }
757: /* On entry, ht contains the topological entities whose dofs we are responsible for solving for;
758: on exit, cht contains all the topological entities we need to compute their residuals.
759: In full generality this should incorporate knowledge of the sparsity pattern of the matrix;
760: here we assume a standard FE sparsity pattern.*/
761: /* TODO: Use DMPlexGetAdjacency() */
762: static PetscErrorCode PCPatchCompleteCellPatch(PC pc, PetscHSetI ht, PetscHSetI cht)
763: {
764: DM dm, plex;
765: PC_PATCH *patch = (PC_PATCH *) pc->data;
766: PetscHashIter hi;
767: PetscInt point;
768: PetscInt *star = NULL, *closure = NULL;
769: PetscInt ignoredim, iStart = 0, iEnd = -1, starSize, closureSize, si, ci;
770: PetscInt *fStar = NULL, *fClosure = NULL;
771: PetscInt fBegin, fEnd, fsi, fci, fStarSize, fClosureSize;
775: PCGetDM(pc, &dm);
776: DMConvert(dm, DMPLEX, &plex);
777: dm = plex;
778: DMPlexGetHeightStratum(dm, 1, &fBegin, &fEnd);
779: PCPatchGetIgnoreDim(pc, &ignoredim);
780: if (ignoredim >= 0) {DMPlexGetDepthStratum(dm, ignoredim, &iStart, &iEnd);}
781: PetscHSetIClear(cht);
782: PetscHashIterBegin(ht, hi);
783: while (!PetscHashIterAtEnd(ht, hi)) {
785: PetscHashIterGetKey(ht, hi, point);
786: PetscHashIterNext(ht, hi);
788: /* Loop over all the cells that this point connects to */
789: DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);
790: for (si = 0; si < starSize*2; si += 2) {
791: const PetscInt ownedpoint = star[si];
792: /* TODO Check for point in cht before running through closure again */
793: /* now loop over all entities in the closure of that cell */
794: DMPlexGetTransitiveClosure(dm, ownedpoint, PETSC_TRUE, &closureSize, &closure);
795: for (ci = 0; ci < closureSize*2; ci += 2) {
796: const PetscInt seenpoint = closure[ci];
797: if (ignoredim >= 0 && seenpoint >= iStart && seenpoint < iEnd) continue;
798: PetscHSetIAdd(cht, seenpoint);
799: /* Facet integrals couple dofs across facets, so in that case for each of
800: * the facets we need to add all dofs on the other side of the facet to
801: * the seen dofs. */
802: if (patch->usercomputeopintfacet) {
803: if (fBegin <= seenpoint && seenpoint < fEnd) {
804: DMPlexGetTransitiveClosure(dm, seenpoint, PETSC_FALSE, &fStarSize, &fStar);
805: for (fsi = 0; fsi < fStarSize*2; fsi += 2) {
806: DMPlexGetTransitiveClosure(dm, fStar[fsi], PETSC_TRUE, &fClosureSize, &fClosure);
807: for (fci = 0; fci < fClosureSize*2; fci += 2) {
808: PetscHSetIAdd(cht, fClosure[fci]);
809: }
810: DMPlexRestoreTransitiveClosure(dm, fStar[fsi], PETSC_TRUE, NULL, &fClosure);
811: }
812: DMPlexRestoreTransitiveClosure(dm, seenpoint, PETSC_FALSE, NULL, &fStar);
813: }
814: }
815: }
816: DMPlexRestoreTransitiveClosure(dm, ownedpoint, PETSC_TRUE, NULL, &closure);
817: }
818: DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, NULL, &star);
819: }
820: DMDestroy(&dm);
821: return(0);
822: }
824: static PetscErrorCode PCPatchGetGlobalDofs(PC pc, PetscSection dofSection[], PetscInt f, PetscBool combined, PetscInt p, PetscInt *dof, PetscInt *off)
825: {
829: if (combined) {
830: if (f < 0) {
831: if (dof) {PetscSectionGetDof(dofSection[0], p, dof);}
832: if (off) {PetscSectionGetOffset(dofSection[0], p, off);}
833: } else {
834: if (dof) {PetscSectionGetFieldDof(dofSection[0], p, f, dof);}
835: if (off) {PetscSectionGetFieldOffset(dofSection[0], p, f, off);}
836: }
837: } else {
838: if (f < 0) {
839: PC_PATCH *patch = (PC_PATCH *) pc->data;
840: PetscInt fdof, g;
842: if (dof) {
843: *dof = 0;
844: for (g = 0; g < patch->nsubspaces; ++g) {
845: PetscSectionGetDof(dofSection[g], p, &fdof);
846: *dof += fdof;
847: }
848: }
849: if (off) {
850: *off = 0;
851: for (g = 0; g < patch->nsubspaces; ++g) {
852: PetscSectionGetOffset(dofSection[g], p, &fdof);
853: *off += fdof;
854: }
855: }
856: } else {
857: if (dof) {PetscSectionGetDof(dofSection[f], p, dof);}
858: if (off) {PetscSectionGetOffset(dofSection[f], p, off);}
859: }
860: }
861: return(0);
862: }
864: /* Given a hash table with a set of topological entities (pts), compute the degrees of
865: freedom in global concatenated numbering on those entities.
866: For Vanka smoothing, this needs to do something special: ignore dofs of the
867: constraint subspace on entities that aren't the base entity we're building the patch
868: around. */
869: static PetscErrorCode PCPatchGetPointDofs(PC pc, PetscHSetI pts, PetscHSetI dofs, PetscInt base, PetscHSetI* subspaces_to_exclude)
870: {
871: PC_PATCH *patch = (PC_PATCH *) pc->data;
872: PetscHashIter hi;
873: PetscInt ldof, loff;
874: PetscInt k, p;
878: PetscHSetIClear(dofs);
879: for (k = 0; k < patch->nsubspaces; ++k) {
880: PetscInt subspaceOffset = patch->subspaceOffsets[k];
881: PetscInt bs = patch->bs[k];
882: PetscInt j, l;
884: if (subspaces_to_exclude != NULL) {
885: PetscBool should_exclude_k = PETSC_FALSE;
886: PetscHSetIHas(*subspaces_to_exclude, k, &should_exclude_k);
887: if (should_exclude_k) {
888: /* only get this subspace dofs at the base entity, not any others */
889: PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, base, &ldof, &loff);
890: if (0 == ldof) continue;
891: for (j = loff; j < ldof + loff; ++j) {
892: for (l = 0; l < bs; ++l) {
893: PetscInt dof = bs*j + l + subspaceOffset;
894: PetscHSetIAdd(dofs, dof);
895: }
896: }
897: continue; /* skip the other dofs of this subspace */
898: }
899: }
901: PetscHashIterBegin(pts, hi);
902: while (!PetscHashIterAtEnd(pts, hi)) {
903: PetscHashIterGetKey(pts, hi, p);
904: PetscHashIterNext(pts, hi);
905: PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, p, &ldof, &loff);
906: if (0 == ldof) continue;
907: for (j = loff; j < ldof + loff; ++j) {
908: for (l = 0; l < bs; ++l) {
909: PetscInt dof = bs*j + l + subspaceOffset;
910: PetscHSetIAdd(dofs, dof);
911: }
912: }
913: }
914: }
915: return(0);
916: }
918: /* Given two hash tables A and B, compute the keys in B that are not in A, and put them in C */
919: static PetscErrorCode PCPatchComputeSetDifference_Private(PetscHSetI A, PetscHSetI B, PetscHSetI C)
920: {
921: PetscHashIter hi;
922: PetscInt key;
923: PetscBool flg;
927: PetscHSetIClear(C);
928: PetscHashIterBegin(B, hi);
929: while (!PetscHashIterAtEnd(B, hi)) {
930: PetscHashIterGetKey(B, hi, key);
931: PetscHashIterNext(B, hi);
932: PetscHSetIHas(A, key, &flg);
933: if (!flg) {PetscHSetIAdd(C, key);}
934: }
935: return(0);
936: }
938: /*
939: * PCPatchCreateCellPatches - create patches.
940: *
941: * Input Parameters:
942: * + dm - The DMPlex object defining the mesh
943: *
944: * Output Parameters:
945: * + cellCounts - Section with counts of cells around each vertex
946: * . cells - IS of the cell point indices of cells in each patch
947: * . pointCounts - Section with counts of cells around each vertex
948: * - point - IS of the cell point indices of cells in each patch
949: */
950: static PetscErrorCode PCPatchCreateCellPatches(PC pc)
951: {
952: PC_PATCH *patch = (PC_PATCH *) pc->data;
953: DMLabel ghost = NULL;
954: DM dm, plex;
955: PetscHSetI ht, cht;
956: PetscSection cellCounts, pointCounts, intFacetCounts, extFacetCounts;
957: PetscInt *cellsArray, *pointsArray, *intFacetsArray, *extFacetsArray, *intFacetsToPatchCell;
958: PetscInt numCells, numPoints, numIntFacets, numExtFacets;
959: const PetscInt *leaves;
960: PetscInt nleaves, pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, v;
961: PetscBool isFiredrake;
962: PetscErrorCode ierr;
965: /* Used to keep track of the cells in the patch. */
966: PetscHSetICreate(&ht);
967: PetscHSetICreate(&cht);
969: PCGetDM(pc, &dm);
970: if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONGSTATE, "DM not yet set on patch PC\n");
971: DMConvert(dm, DMPLEX, &plex);
972: dm = plex;
973: DMPlexGetChart(dm, &pStart, &pEnd);
974: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
976: if (patch->user_patches) {
977: patch->userpatchconstructionop(pc, &patch->npatch, &patch->userIS, &patch->iterationSet, patch->userpatchconstructctx);
978: vStart = 0; vEnd = patch->npatch;
979: } else if (patch->ctype == PC_PATCH_PARDECOMP) {
980: vStart = 0; vEnd = 1;
981: } else if (patch->codim < 0) {
982: if (patch->dim < 0) {DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);}
983: else {DMPlexGetDepthStratum(dm, patch->dim, &vStart, &vEnd);}
984: } else {DMPlexGetHeightStratum(dm, patch->codim, &vStart, &vEnd);}
985: patch->npatch = vEnd - vStart;
987: /* These labels mark the owned points. We only create patches around points that this process owns. */
988: DMHasLabel(dm, "pyop2_ghost", &isFiredrake);
989: if (isFiredrake) {
990: DMGetLabel(dm, "pyop2_ghost", &ghost);
991: DMLabelCreateIndex(ghost, pStart, pEnd);
992: } else {
993: PetscSF sf;
995: DMGetPointSF(dm, &sf);
996: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
997: nleaves = PetscMax(nleaves, 0);
998: }
1000: PetscSectionCreate(PETSC_COMM_SELF, &patch->cellCounts);
1001: PetscObjectSetName((PetscObject) patch->cellCounts, "Patch Cell Layout");
1002: cellCounts = patch->cellCounts;
1003: PetscSectionSetChart(cellCounts, vStart, vEnd);
1004: PetscSectionCreate(PETSC_COMM_SELF, &patch->pointCounts);
1005: PetscObjectSetName((PetscObject) patch->pointCounts, "Patch Point Layout");
1006: pointCounts = patch->pointCounts;
1007: PetscSectionSetChart(pointCounts, vStart, vEnd);
1008: PetscSectionCreate(PETSC_COMM_SELF, &patch->extFacetCounts);
1009: PetscObjectSetName((PetscObject) patch->extFacetCounts, "Patch Exterior Facet Layout");
1010: extFacetCounts = patch->extFacetCounts;
1011: PetscSectionSetChart(extFacetCounts, vStart, vEnd);
1012: PetscSectionCreate(PETSC_COMM_SELF, &patch->intFacetCounts);
1013: PetscObjectSetName((PetscObject) patch->intFacetCounts, "Patch Interior Facet Layout");
1014: intFacetCounts = patch->intFacetCounts;
1015: PetscSectionSetChart(intFacetCounts, vStart, vEnd);
1016: /* Count cells and points in the patch surrounding each entity */
1017: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1018: for (v = vStart; v < vEnd; ++v) {
1019: PetscHashIter hi;
1020: PetscInt chtSize, loc = -1;
1021: PetscBool flg;
1023: if (!patch->user_patches && patch->ctype != PC_PATCH_PARDECOMP) {
1024: if (ghost) {DMLabelHasPoint(ghost, v, &flg);}
1025: else {PetscFindInt(v, nleaves, leaves, &loc); flg = loc >=0 ? PETSC_TRUE : PETSC_FALSE;}
1026: /* Not an owned entity, don't make a cell patch. */
1027: if (flg) continue;
1028: }
1030: patch->patchconstructop((void *) patch, dm, v, ht);
1031: PCPatchCompleteCellPatch(pc, ht, cht);
1032: PetscHSetIGetSize(cht, &chtSize);
1033: /* empty patch, continue */
1034: if (chtSize == 0) continue;
1036: /* safe because size(cht) > 0 from above */
1037: PetscHashIterBegin(cht, hi);
1038: while (!PetscHashIterAtEnd(cht, hi)) {
1039: PetscInt point, pdof;
1041: PetscHashIterGetKey(cht, hi, point);
1042: if (fStart <= point && point < fEnd) {
1043: const PetscInt *support;
1044: PetscInt supportSize, p;
1045: PetscBool interior = PETSC_TRUE;
1046: DMPlexGetSupport(dm, point, &support);
1047: DMPlexGetSupportSize(dm, point, &supportSize);
1048: if (supportSize == 1) {
1049: interior = PETSC_FALSE;
1050: } else {
1051: for (p = 0; p < supportSize; p++) {
1052: PetscBool found;
1053: /* FIXME: can I do this while iterating over cht? */
1054: PetscHSetIHas(cht, support[p], &found);
1055: if (!found) {
1056: interior = PETSC_FALSE;
1057: break;
1058: }
1059: }
1060: }
1061: if (interior) {
1062: PetscSectionAddDof(intFacetCounts, v, 1);
1063: } else {
1064: PetscSectionAddDof(extFacetCounts, v, 1);
1065: }
1066: }
1067: PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);
1068: if (pdof) {PetscSectionAddDof(pointCounts, v, 1);}
1069: if (point >= cStart && point < cEnd) {PetscSectionAddDof(cellCounts, v, 1);}
1070: PetscHashIterNext(cht, hi);
1071: }
1072: }
1073: if (isFiredrake) {DMLabelDestroyIndex(ghost);}
1075: PetscSectionSetUp(cellCounts);
1076: PetscSectionGetStorageSize(cellCounts, &numCells);
1077: PetscMalloc1(numCells, &cellsArray);
1078: PetscSectionSetUp(pointCounts);
1079: PetscSectionGetStorageSize(pointCounts, &numPoints);
1080: PetscMalloc1(numPoints, &pointsArray);
1082: PetscSectionSetUp(intFacetCounts);
1083: PetscSectionSetUp(extFacetCounts);
1084: PetscSectionGetStorageSize(intFacetCounts, &numIntFacets);
1085: PetscSectionGetStorageSize(extFacetCounts, &numExtFacets);
1086: PetscMalloc1(numIntFacets, &intFacetsArray);
1087: PetscMalloc1(numIntFacets*2, &intFacetsToPatchCell);
1088: PetscMalloc1(numExtFacets, &extFacetsArray);
1090: /* Now that we know how much space we need, run through again and actually remember the cells. */
1091: for (v = vStart; v < vEnd; v++) {
1092: PetscHashIter hi;
1093: PetscInt dof, off, cdof, coff, efdof, efoff, ifdof, ifoff, pdof, n = 0, cn = 0, ifn = 0, efn = 0;
1095: PetscSectionGetDof(pointCounts, v, &dof);
1096: PetscSectionGetOffset(pointCounts, v, &off);
1097: PetscSectionGetDof(cellCounts, v, &cdof);
1098: PetscSectionGetOffset(cellCounts, v, &coff);
1099: PetscSectionGetDof(intFacetCounts, v, &ifdof);
1100: PetscSectionGetOffset(intFacetCounts, v, &ifoff);
1101: PetscSectionGetDof(extFacetCounts, v, &efdof);
1102: PetscSectionGetOffset(extFacetCounts, v, &efoff);
1103: if (dof <= 0) continue;
1104: patch->patchconstructop((void *) patch, dm, v, ht);
1105: PCPatchCompleteCellPatch(pc, ht, cht);
1106: PetscHashIterBegin(cht, hi);
1107: while (!PetscHashIterAtEnd(cht, hi)) {
1108: PetscInt point;
1110: PetscHashIterGetKey(cht, hi, point);
1111: if (fStart <= point && point < fEnd) {
1112: const PetscInt *support;
1113: PetscInt supportSize, p;
1114: PetscBool interior = PETSC_TRUE;
1115: DMPlexGetSupport(dm, point, &support);
1116: DMPlexGetSupportSize(dm, point, &supportSize);
1117: if (supportSize == 1) {
1118: interior = PETSC_FALSE;
1119: } else {
1120: for (p = 0; p < supportSize; p++) {
1121: PetscBool found;
1122: /* FIXME: can I do this while iterating over cht? */
1123: PetscHSetIHas(cht, support[p], &found);
1124: if (!found) {
1125: interior = PETSC_FALSE;
1126: break;
1127: }
1128: }
1129: }
1130: if (interior) {
1131: intFacetsToPatchCell[2*(ifoff + ifn)] = support[0];
1132: intFacetsToPatchCell[2*(ifoff + ifn) + 1] = support[1];
1133: intFacetsArray[ifoff + ifn++] = point;
1134: } else {
1135: extFacetsArray[efoff + efn++] = point;
1136: }
1137: }
1138: PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);
1139: if (pdof) {pointsArray[off + n++] = point;}
1140: if (point >= cStart && point < cEnd) {cellsArray[coff + cn++] = point;}
1141: PetscHashIterNext(cht, hi);
1142: }
1143: if (ifn != ifdof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of interior facets in patch %D is %D, but should be %D", v, ifn, ifdof);
1144: if (efn != efdof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of exterior facets in patch %D is %D, but should be %D", v, efn, efdof);
1145: if (cn != cdof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of cells in patch %D is %D, but should be %D", v, cn, cdof);
1146: if (n != dof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of points in patch %D is %D, but should be %D", v, n, dof);
1148: for (ifn = 0; ifn < ifdof; ifn++) {
1149: PetscInt cell0 = intFacetsToPatchCell[2*(ifoff + ifn)];
1150: PetscInt cell1 = intFacetsToPatchCell[2*(ifoff + ifn) + 1];
1151: PetscBool found0 = PETSC_FALSE, found1 = PETSC_FALSE;
1152: for (n = 0; n < cdof; n++) {
1153: if (!found0 && cell0 == cellsArray[coff + n]) {
1154: intFacetsToPatchCell[2*(ifoff + ifn)] = n;
1155: found0 = PETSC_TRUE;
1156: }
1157: if (!found1 && cell1 == cellsArray[coff + n]) {
1158: intFacetsToPatchCell[2*(ifoff + ifn) + 1] = n;
1159: found1 = PETSC_TRUE;
1160: }
1161: if (found0 && found1) break;
1162: }
1163: if (!(found0 && found1)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Didn't manage to find local point numbers for facet support");
1164: }
1165: }
1166: PetscHSetIDestroy(&ht);
1167: PetscHSetIDestroy(&cht);
1169: ISCreateGeneral(PETSC_COMM_SELF, numCells, cellsArray, PETSC_OWN_POINTER, &patch->cells);
1170: PetscObjectSetName((PetscObject) patch->cells, "Patch Cells");
1171: if (patch->viewCells) {
1172: ObjectView((PetscObject) patch->cellCounts, patch->viewerCells, patch->formatCells);
1173: ObjectView((PetscObject) patch->cells, patch->viewerCells, patch->formatCells);
1174: }
1175: ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray, PETSC_OWN_POINTER, &patch->intFacets);
1176: PetscObjectSetName((PetscObject) patch->intFacets, "Patch Interior Facets");
1177: ISCreateGeneral(PETSC_COMM_SELF, 2*numIntFacets, intFacetsToPatchCell, PETSC_OWN_POINTER, &patch->intFacetsToPatchCell);
1178: PetscObjectSetName((PetscObject) patch->intFacetsToPatchCell, "Patch Interior Facets local support");
1179: if (patch->viewIntFacets) {
1180: ObjectView((PetscObject) patch->intFacetCounts, patch->viewerIntFacets, patch->formatIntFacets);
1181: ObjectView((PetscObject) patch->intFacets, patch->viewerIntFacets, patch->formatIntFacets);
1182: ObjectView((PetscObject) patch->intFacetsToPatchCell, patch->viewerIntFacets, patch->formatIntFacets);
1183: }
1184: ISCreateGeneral(PETSC_COMM_SELF, numExtFacets, extFacetsArray, PETSC_OWN_POINTER, &patch->extFacets);
1185: PetscObjectSetName((PetscObject) patch->extFacets, "Patch Exterior Facets");
1186: if (patch->viewExtFacets) {
1187: ObjectView((PetscObject) patch->extFacetCounts, patch->viewerExtFacets, patch->formatExtFacets);
1188: ObjectView((PetscObject) patch->extFacets, patch->viewerExtFacets, patch->formatExtFacets);
1189: }
1190: ISCreateGeneral(PETSC_COMM_SELF, numPoints, pointsArray, PETSC_OWN_POINTER, &patch->points);
1191: PetscObjectSetName((PetscObject) patch->points, "Patch Points");
1192: if (patch->viewPoints) {
1193: ObjectView((PetscObject) patch->pointCounts, patch->viewerPoints, patch->formatPoints);
1194: ObjectView((PetscObject) patch->points, patch->viewerPoints, patch->formatPoints);
1195: }
1196: DMDestroy(&dm);
1197: return(0);
1198: }
1200: /*
1201: * PCPatchCreateCellPatchDiscretisationInfo - Build the dof maps for cell patches
1202: *
1203: * Input Parameters:
1204: * + dm - The DMPlex object defining the mesh
1205: * . cellCounts - Section with counts of cells around each vertex
1206: * . cells - IS of the cell point indices of cells in each patch
1207: * . cellNumbering - Section mapping plex cell points to Firedrake cell indices.
1208: * . nodesPerCell - number of nodes per cell.
1209: * - cellNodeMap - map from cells to node indices (nodesPerCell * numCells)
1210: *
1211: * Output Parameters:
1212: * + dofs - IS of local dof numbers of each cell in the patch, where local is a patch local numbering
1213: * . gtolCounts - Section with counts of dofs per cell patch
1214: * - gtol - IS mapping from global dofs to local dofs for each patch.
1215: */
1216: static PetscErrorCode PCPatchCreateCellPatchDiscretisationInfo(PC pc)
1217: {
1218: PC_PATCH *patch = (PC_PATCH *) pc->data;
1219: PetscSection cellCounts = patch->cellCounts;
1220: PetscSection pointCounts = patch->pointCounts;
1221: PetscSection gtolCounts, gtolCountsWithArtificial = NULL, gtolCountsWithAll = NULL;
1222: IS cells = patch->cells;
1223: IS points = patch->points;
1224: PetscSection cellNumbering = patch->cellNumbering;
1225: PetscInt Nf = patch->nsubspaces;
1226: PetscInt numCells, numPoints;
1227: PetscInt numDofs;
1228: PetscInt numGlobalDofs, numGlobalDofsWithArtificial, numGlobalDofsWithAll;
1229: PetscInt totalDofsPerCell = patch->totalDofsPerCell;
1230: PetscInt vStart, vEnd, v;
1231: const PetscInt *cellsArray, *pointsArray;
1232: PetscInt *newCellsArray = NULL;
1233: PetscInt *dofsArray = NULL;
1234: PetscInt *dofsArrayWithArtificial = NULL;
1235: PetscInt *dofsArrayWithAll = NULL;
1236: PetscInt *offsArray = NULL;
1237: PetscInt *offsArrayWithArtificial = NULL;
1238: PetscInt *offsArrayWithAll = NULL;
1239: PetscInt *asmArray = NULL;
1240: PetscInt *asmArrayWithArtificial = NULL;
1241: PetscInt *asmArrayWithAll = NULL;
1242: PetscInt *globalDofsArray = NULL;
1243: PetscInt *globalDofsArrayWithArtificial = NULL;
1244: PetscInt *globalDofsArrayWithAll = NULL;
1245: PetscInt globalIndex = 0;
1246: PetscInt key = 0;
1247: PetscInt asmKey = 0;
1248: DM dm = NULL, plex;
1249: const PetscInt *bcNodes = NULL;
1250: PetscHMapI ht;
1251: PetscHMapI htWithArtificial;
1252: PetscHMapI htWithAll;
1253: PetscHSetI globalBcs;
1254: PetscInt numBcs;
1255: PetscHSetI ownedpts, seenpts, owneddofs, seendofs, artificialbcs;
1256: PetscInt pStart, pEnd, p, i;
1257: char option[PETSC_MAX_PATH_LEN];
1258: PetscBool isNonlinear;
1259: PetscErrorCode ierr;
1263: PCGetDM(pc, &dm);
1264: DMConvert(dm, DMPLEX, &plex);
1265: dm = plex;
1266: /* dofcounts section is cellcounts section * dofPerCell */
1267: PetscSectionGetStorageSize(cellCounts, &numCells);
1268: PetscSectionGetStorageSize(patch->pointCounts, &numPoints);
1269: numDofs = numCells * totalDofsPerCell;
1270: PetscMalloc1(numDofs, &dofsArray);
1271: PetscMalloc1(numPoints*Nf, &offsArray);
1272: PetscMalloc1(numDofs, &asmArray);
1273: PetscMalloc1(numCells, &newCellsArray);
1274: PetscSectionGetChart(cellCounts, &vStart, &vEnd);
1275: PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCounts);
1276: gtolCounts = patch->gtolCounts;
1277: PetscSectionSetChart(gtolCounts, vStart, vEnd);
1278: PetscObjectSetName((PetscObject) patch->gtolCounts, "Patch Global Index Section");
1280: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1281: PetscMalloc1(numPoints*Nf, &offsArrayWithArtificial);
1282: PetscMalloc1(numDofs, &asmArrayWithArtificial);
1283: PetscMalloc1(numDofs, &dofsArrayWithArtificial);
1284: PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithArtificial);
1285: gtolCountsWithArtificial = patch->gtolCountsWithArtificial;
1286: PetscSectionSetChart(gtolCountsWithArtificial, vStart, vEnd);
1287: PetscObjectSetName((PetscObject) patch->gtolCountsWithArtificial, "Patch Global Index Section Including Artificial BCs");
1288: }
1290: isNonlinear = patch->isNonlinear;
1291: if (isNonlinear) {
1292: PetscMalloc1(numPoints*Nf, &offsArrayWithAll);
1293: PetscMalloc1(numDofs, &asmArrayWithAll);
1294: PetscMalloc1(numDofs, &dofsArrayWithAll);
1295: PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithAll);
1296: gtolCountsWithAll = patch->gtolCountsWithAll;
1297: PetscSectionSetChart(gtolCountsWithAll, vStart, vEnd);
1298: PetscObjectSetName((PetscObject) patch->gtolCountsWithAll, "Patch Global Index Section Including All BCs");
1299: }
1301: /* Outside the patch loop, get the dofs that are globally-enforced Dirichlet
1302: conditions */
1303: PetscHSetICreate(&globalBcs);
1304: ISGetIndices(patch->ghostBcNodes, &bcNodes);
1305: ISGetSize(patch->ghostBcNodes, &numBcs);
1306: for (i = 0; i < numBcs; ++i) {
1307: PetscHSetIAdd(globalBcs, bcNodes[i]); /* these are already in concatenated numbering */
1308: }
1309: ISRestoreIndices(patch->ghostBcNodes, &bcNodes);
1310: ISDestroy(&patch->ghostBcNodes); /* memory optimisation */
1312: /* Hash tables for artificial BC construction */
1313: PetscHSetICreate(&ownedpts);
1314: PetscHSetICreate(&seenpts);
1315: PetscHSetICreate(&owneddofs);
1316: PetscHSetICreate(&seendofs);
1317: PetscHSetICreate(&artificialbcs);
1319: ISGetIndices(cells, &cellsArray);
1320: ISGetIndices(points, &pointsArray);
1321: PetscHMapICreate(&ht);
1322: PetscHMapICreate(&htWithArtificial);
1323: PetscHMapICreate(&htWithAll);
1324: for (v = vStart; v < vEnd; ++v) {
1325: PetscInt localIndex = 0;
1326: PetscInt localIndexWithArtificial = 0;
1327: PetscInt localIndexWithAll = 0;
1328: PetscInt dof, off, i, j, k, l;
1330: PetscHMapIClear(ht);
1331: PetscHMapIClear(htWithArtificial);
1332: PetscHMapIClear(htWithAll);
1333: PetscSectionGetDof(cellCounts, v, &dof);
1334: PetscSectionGetOffset(cellCounts, v, &off);
1335: if (dof <= 0) continue;
1337: /* Calculate the global numbers of the artificial BC dofs here first */
1338: patch->patchconstructop((void*)patch, dm, v, ownedpts);
1339: PCPatchCompleteCellPatch(pc, ownedpts, seenpts);
1340: PCPatchGetPointDofs(pc, ownedpts, owneddofs, v, &patch->subspaces_to_exclude);
1341: PCPatchGetPointDofs(pc, seenpts, seendofs, v, NULL);
1342: PCPatchComputeSetDifference_Private(owneddofs, seendofs, artificialbcs);
1343: if (patch->viewPatches) {
1344: PetscHSetI globalbcdofs;
1345: PetscHashIter hi;
1346: MPI_Comm comm = PetscObjectComm((PetscObject)pc);
1348: PetscHSetICreate(&globalbcdofs);
1349: PetscSynchronizedPrintf(comm, "Patch %d: owned dofs:\n", v);
1350: PetscHashIterBegin(owneddofs, hi);
1351: while (!PetscHashIterAtEnd(owneddofs, hi)) {
1352: PetscInt globalDof;
1354: PetscHashIterGetKey(owneddofs, hi, globalDof);
1355: PetscHashIterNext(owneddofs, hi);
1356: PetscSynchronizedPrintf(comm, "%d ", globalDof);
1357: }
1358: PetscSynchronizedPrintf(comm, "\n");
1359: PetscSynchronizedPrintf(comm, "Patch %d: seen dofs:\n", v);
1360: PetscHashIterBegin(seendofs, hi);
1361: while (!PetscHashIterAtEnd(seendofs, hi)) {
1362: PetscInt globalDof;
1363: PetscBool flg;
1365: PetscHashIterGetKey(seendofs, hi, globalDof);
1366: PetscHashIterNext(seendofs, hi);
1367: PetscSynchronizedPrintf(comm, "%d ", globalDof);
1369: PetscHSetIHas(globalBcs, globalDof, &flg);
1370: if (flg) {PetscHSetIAdd(globalbcdofs, globalDof);}
1371: }
1372: PetscSynchronizedPrintf(comm, "\n");
1373: PetscSynchronizedPrintf(comm, "Patch %d: global BCs:\n", v);
1374: PetscHSetIGetSize(globalbcdofs, &numBcs);
1375: if (numBcs > 0) {
1376: PetscHashIterBegin(globalbcdofs, hi);
1377: while (!PetscHashIterAtEnd(globalbcdofs, hi)) {
1378: PetscInt globalDof;
1379: PetscHashIterGetKey(globalbcdofs, hi, globalDof);
1380: PetscHashIterNext(globalbcdofs, hi);
1381: PetscSynchronizedPrintf(comm, "%d ", globalDof);
1382: }
1383: }
1384: PetscSynchronizedPrintf(comm, "\n");
1385: PetscSynchronizedPrintf(comm, "Patch %d: artificial BCs:\n", v);
1386: PetscHSetIGetSize(artificialbcs, &numBcs);
1387: if (numBcs > 0) {
1388: PetscHashIterBegin(artificialbcs, hi);
1389: while (!PetscHashIterAtEnd(artificialbcs, hi)) {
1390: PetscInt globalDof;
1391: PetscHashIterGetKey(artificialbcs, hi, globalDof);
1392: PetscHashIterNext(artificialbcs, hi);
1393: PetscSynchronizedPrintf(comm, "%d ", globalDof);
1394: }
1395: }
1396: PetscSynchronizedPrintf(comm, "\n\n");
1397: PetscHSetIDestroy(&globalbcdofs);
1398: }
1399: for (k = 0; k < patch->nsubspaces; ++k) {
1400: const PetscInt *cellNodeMap = patch->cellNodeMap[k];
1401: PetscInt nodesPerCell = patch->nodesPerCell[k];
1402: PetscInt subspaceOffset = patch->subspaceOffsets[k];
1403: PetscInt bs = patch->bs[k];
1405: for (i = off; i < off + dof; ++i) {
1406: /* Walk over the cells in this patch. */
1407: const PetscInt c = cellsArray[i];
1408: PetscInt cell = c;
1410: /* TODO Change this to an IS */
1411: if (cellNumbering) {
1412: PetscSectionGetDof(cellNumbering, c, &cell);
1413: if (cell <= 0) SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_OUTOFRANGE, "Cell %D doesn't appear in cell numbering map", c);
1414: PetscSectionGetOffset(cellNumbering, c, &cell);
1415: }
1416: newCellsArray[i] = cell;
1417: for (j = 0; j < nodesPerCell; ++j) {
1418: /* For each global dof, map it into contiguous local storage. */
1419: const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + subspaceOffset;
1420: /* finally, loop over block size */
1421: for (l = 0; l < bs; ++l) {
1422: PetscInt localDof;
1423: PetscBool isGlobalBcDof, isArtificialBcDof;
1425: /* first, check if this is either a globally enforced or locally enforced BC dof */
1426: PetscHSetIHas(globalBcs, globalDof + l, &isGlobalBcDof);
1427: PetscHSetIHas(artificialbcs, globalDof + l, &isArtificialBcDof);
1429: /* if it's either, don't ever give it a local dof number */
1430: if (isGlobalBcDof || isArtificialBcDof) {
1431: dofsArray[globalIndex] = -1; /* don't use this in assembly in this patch */
1432: } else {
1433: PetscHMapIGet(ht, globalDof + l, &localDof);
1434: if (localDof == -1) {
1435: localDof = localIndex++;
1436: PetscHMapISet(ht, globalDof + l, localDof);
1437: }
1438: if (globalIndex >= numDofs) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs);
1439: /* And store. */
1440: dofsArray[globalIndex] = localDof;
1441: }
1443: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1444: if (isGlobalBcDof) {
1445: dofsArrayWithArtificial[globalIndex] = -1; /* don't use this in assembly in this patch */
1446: } else {
1447: PetscHMapIGet(htWithArtificial, globalDof + l, &localDof);
1448: if (localDof == -1) {
1449: localDof = localIndexWithArtificial++;
1450: PetscHMapISet(htWithArtificial, globalDof + l, localDof);
1451: }
1452: if (globalIndex >= numDofs) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs);
1453: /* And store.*/
1454: dofsArrayWithArtificial[globalIndex] = localDof;
1455: }
1456: }
1458: if (isNonlinear) {
1459: /* Build the dofmap for the function space with _all_ dofs,
1460: including those in any kind of boundary condition */
1461: PetscHMapIGet(htWithAll, globalDof + l, &localDof);
1462: if (localDof == -1) {
1463: localDof = localIndexWithAll++;
1464: PetscHMapISet(htWithAll, globalDof + l, localDof);
1465: }
1466: if (globalIndex >= numDofs) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs);
1467: /* And store.*/
1468: dofsArrayWithAll[globalIndex] = localDof;
1469: }
1470: globalIndex++;
1471: }
1472: }
1473: }
1474: }
1475: /*How many local dofs in this patch? */
1476: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1477: PetscHMapIGetSize(htWithArtificial, &dof);
1478: PetscSectionSetDof(gtolCountsWithArtificial, v, dof);
1479: }
1480: if (isNonlinear) {
1481: PetscHMapIGetSize(htWithAll, &dof);
1482: PetscSectionSetDof(gtolCountsWithAll, v, dof);
1483: }
1484: PetscHMapIGetSize(ht, &dof);
1485: PetscSectionSetDof(gtolCounts, v, dof);
1486: }
1488: DMDestroy(&dm);
1489: if (globalIndex != numDofs) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Expected number of dofs (%d) doesn't match found number (%d)", numDofs, globalIndex);
1490: PetscSectionSetUp(gtolCounts);
1491: PetscSectionGetStorageSize(gtolCounts, &numGlobalDofs);
1492: PetscMalloc1(numGlobalDofs, &globalDofsArray);
1494: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1495: PetscSectionSetUp(gtolCountsWithArtificial);
1496: PetscSectionGetStorageSize(gtolCountsWithArtificial, &numGlobalDofsWithArtificial);
1497: PetscMalloc1(numGlobalDofsWithArtificial, &globalDofsArrayWithArtificial);
1498: }
1499: if (isNonlinear) {
1500: PetscSectionSetUp(gtolCountsWithAll);
1501: PetscSectionGetStorageSize(gtolCountsWithAll, &numGlobalDofsWithAll);
1502: PetscMalloc1(numGlobalDofsWithAll, &globalDofsArrayWithAll);
1503: }
1504: /* Now populate the global to local map. This could be merged into the above loop if we were willing to deal with reallocs. */
1505: for (v = vStart; v < vEnd; ++v) {
1506: PetscHashIter hi;
1507: PetscInt dof, off, Np, ooff, i, j, k, l;
1509: PetscHMapIClear(ht);
1510: PetscHMapIClear(htWithArtificial);
1511: PetscHMapIClear(htWithAll);
1512: PetscSectionGetDof(cellCounts, v, &dof);
1513: PetscSectionGetOffset(cellCounts, v, &off);
1514: PetscSectionGetDof(pointCounts, v, &Np);
1515: PetscSectionGetOffset(pointCounts, v, &ooff);
1516: if (dof <= 0) continue;
1518: for (k = 0; k < patch->nsubspaces; ++k) {
1519: const PetscInt *cellNodeMap = patch->cellNodeMap[k];
1520: PetscInt nodesPerCell = patch->nodesPerCell[k];
1521: PetscInt subspaceOffset = patch->subspaceOffsets[k];
1522: PetscInt bs = patch->bs[k];
1523: PetscInt goff;
1525: for (i = off; i < off + dof; ++i) {
1526: /* Reconstruct mapping of global-to-local on this patch. */
1527: const PetscInt c = cellsArray[i];
1528: PetscInt cell = c;
1530: if (cellNumbering) {PetscSectionGetOffset(cellNumbering, c, &cell);}
1531: for (j = 0; j < nodesPerCell; ++j) {
1532: for (l = 0; l < bs; ++l) {
1533: const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset;
1534: const PetscInt localDof = dofsArray[key];
1535: if (localDof >= 0) {PetscHMapISet(ht, globalDof, localDof);}
1536: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1537: const PetscInt localDofWithArtificial = dofsArrayWithArtificial[key];
1538: if (localDofWithArtificial >= 0) {
1539: PetscHMapISet(htWithArtificial, globalDof, localDofWithArtificial);
1540: }
1541: }
1542: if (isNonlinear) {
1543: const PetscInt localDofWithAll = dofsArrayWithAll[key];
1544: if (localDofWithAll >= 0) {
1545: PetscHMapISet(htWithAll, globalDof, localDofWithAll);
1546: }
1547: }
1548: key++;
1549: }
1550: }
1551: }
1553: /* Shove it in the output data structure. */
1554: PetscSectionGetOffset(gtolCounts, v, &goff);
1555: PetscHashIterBegin(ht, hi);
1556: while (!PetscHashIterAtEnd(ht, hi)) {
1557: PetscInt globalDof, localDof;
1559: PetscHashIterGetKey(ht, hi, globalDof);
1560: PetscHashIterGetVal(ht, hi, localDof);
1561: if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof;
1562: PetscHashIterNext(ht, hi);
1563: }
1565: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1566: PetscSectionGetOffset(gtolCountsWithArtificial, v, &goff);
1567: PetscHashIterBegin(htWithArtificial, hi);
1568: while (!PetscHashIterAtEnd(htWithArtificial, hi)) {
1569: PetscInt globalDof, localDof;
1570: PetscHashIterGetKey(htWithArtificial, hi, globalDof);
1571: PetscHashIterGetVal(htWithArtificial, hi, localDof);
1572: if (globalDof >= 0) globalDofsArrayWithArtificial[goff + localDof] = globalDof;
1573: PetscHashIterNext(htWithArtificial, hi);
1574: }
1575: }
1576: if (isNonlinear) {
1577: PetscSectionGetOffset(gtolCountsWithAll, v, &goff);
1578: PetscHashIterBegin(htWithAll, hi);
1579: while (!PetscHashIterAtEnd(htWithAll, hi)) {
1580: PetscInt globalDof, localDof;
1581: PetscHashIterGetKey(htWithAll, hi, globalDof);
1582: PetscHashIterGetVal(htWithAll, hi, localDof);
1583: if (globalDof >= 0) globalDofsArrayWithAll[goff + localDof] = globalDof;
1584: PetscHashIterNext(htWithAll, hi);
1585: }
1586: }
1588: for (p = 0; p < Np; ++p) {
1589: const PetscInt point = pointsArray[ooff + p];
1590: PetscInt globalDof, localDof;
1592: PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof);
1593: PetscHMapIGet(ht, globalDof, &localDof);
1594: offsArray[(ooff + p)*Nf + k] = localDof;
1595: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1596: PetscHMapIGet(htWithArtificial, globalDof, &localDof);
1597: offsArrayWithArtificial[(ooff + p)*Nf + k] = localDof;
1598: }
1599: if (isNonlinear) {
1600: PetscHMapIGet(htWithAll, globalDof, &localDof);
1601: offsArrayWithAll[(ooff + p)*Nf + k] = localDof;
1602: }
1603: }
1604: }
1606: PetscHSetIDestroy(&globalBcs);
1607: PetscHSetIDestroy(&ownedpts);
1608: PetscHSetIDestroy(&seenpts);
1609: PetscHSetIDestroy(&owneddofs);
1610: PetscHSetIDestroy(&seendofs);
1611: PetscHSetIDestroy(&artificialbcs);
1613: /* At this point, we have a hash table ht built that maps globalDof -> localDof.
1614: We need to create the dof table laid out cellwise first, then by subspace,
1615: as the assembler assembles cell-wise and we need to stuff the different
1616: contributions of the different function spaces to the right places. So we loop
1617: over cells, then over subspaces. */
1618: if (patch->nsubspaces > 1) { /* for nsubspaces = 1, data we need is already in dofsArray */
1619: for (i = off; i < off + dof; ++i) {
1620: const PetscInt c = cellsArray[i];
1621: PetscInt cell = c;
1623: if (cellNumbering) {PetscSectionGetOffset(cellNumbering, c, &cell);}
1624: for (k = 0; k < patch->nsubspaces; ++k) {
1625: const PetscInt *cellNodeMap = patch->cellNodeMap[k];
1626: PetscInt nodesPerCell = patch->nodesPerCell[k];
1627: PetscInt subspaceOffset = patch->subspaceOffsets[k];
1628: PetscInt bs = patch->bs[k];
1630: for (j = 0; j < nodesPerCell; ++j) {
1631: for (l = 0; l < bs; ++l) {
1632: const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset;
1633: PetscInt localDof;
1635: PetscHMapIGet(ht, globalDof, &localDof);
1636: /* If it's not in the hash table, i.e. is a BC dof,
1637: then the PetscHSetIMap above gives -1, which matches
1638: exactly the convention for PETSc's matrix assembly to
1639: ignore the dof. So we don't need to do anything here */
1640: asmArray[asmKey] = localDof;
1641: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1642: PetscHMapIGet(htWithArtificial, globalDof, &localDof);
1643: asmArrayWithArtificial[asmKey] = localDof;
1644: }
1645: if (isNonlinear) {
1646: PetscHMapIGet(htWithAll, globalDof, &localDof);
1647: asmArrayWithAll[asmKey] = localDof;
1648: }
1649: asmKey++;
1650: }
1651: }
1652: }
1653: }
1654: }
1655: }
1656: if (1 == patch->nsubspaces) {
1657: PetscArraycpy(asmArray, dofsArray, numDofs);
1658: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1659: PetscArraycpy(asmArrayWithArtificial, dofsArrayWithArtificial, numDofs);
1660: }
1661: if (isNonlinear) {
1662: PetscArraycpy(asmArrayWithAll, dofsArrayWithAll, numDofs);
1663: }
1664: }
1666: PetscHMapIDestroy(&ht);
1667: PetscHMapIDestroy(&htWithArtificial);
1668: PetscHMapIDestroy(&htWithAll);
1669: ISRestoreIndices(cells, &cellsArray);
1670: ISRestoreIndices(points, &pointsArray);
1671: PetscFree(dofsArray);
1672: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1673: PetscFree(dofsArrayWithArtificial);
1674: }
1675: if (isNonlinear) {
1676: PetscFree(dofsArrayWithAll);
1677: }
1678: /* Create placeholder section for map from points to patch dofs */
1679: PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection);
1680: PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces);
1681: if (patch->combined) {
1682: PetscInt numFields;
1683: PetscSectionGetNumFields(patch->dofSection[0], &numFields);
1684: if (numFields != patch->nsubspaces) SETERRQ2(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Mismatch between number of section fields %D and number of subspaces %D", numFields, patch->nsubspaces);
1685: PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd);
1686: PetscSectionSetChart(patch->patchSection, pStart, pEnd);
1687: for (p = pStart; p < pEnd; ++p) {
1688: PetscInt dof, fdof, f;
1690: PetscSectionGetDof(patch->dofSection[0], p, &dof);
1691: PetscSectionSetDof(patch->patchSection, p, dof);
1692: for (f = 0; f < patch->nsubspaces; ++f) {
1693: PetscSectionGetFieldDof(patch->dofSection[0], p, f, &fdof);
1694: PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);
1695: }
1696: }
1697: } else {
1698: PetscInt pStartf, pEndf, f;
1699: pStart = PETSC_MAX_INT;
1700: pEnd = PETSC_MIN_INT;
1701: for (f = 0; f < patch->nsubspaces; ++f) {
1702: PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);
1703: pStart = PetscMin(pStart, pStartf);
1704: pEnd = PetscMax(pEnd, pEndf);
1705: }
1706: PetscSectionSetChart(patch->patchSection, pStart, pEnd);
1707: for (f = 0; f < patch->nsubspaces; ++f) {
1708: PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);
1709: for (p = pStartf; p < pEndf; ++p) {
1710: PetscInt fdof;
1711: PetscSectionGetDof(patch->dofSection[f], p, &fdof);
1712: PetscSectionAddDof(patch->patchSection, p, fdof);
1713: PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);
1714: }
1715: }
1716: }
1717: PetscSectionSetUp(patch->patchSection);
1718: PetscSectionSetUseFieldOffsets(patch->patchSection, PETSC_TRUE);
1719: /* Replace cell indices with firedrake-numbered ones. */
1720: ISGeneralSetIndices(cells, numCells, (const PetscInt *) newCellsArray, PETSC_OWN_POINTER);
1721: ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofs, globalDofsArray, PETSC_OWN_POINTER, &patch->gtol);
1722: PetscObjectSetName((PetscObject) patch->gtol, "Global Indices");
1723: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_g2l_view", patch->classname);
1724: PetscSectionViewFromOptions(patch->gtolCounts, (PetscObject) pc, option);
1725: ISViewFromOptions(patch->gtol, (PetscObject) pc, option);
1726: ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArray, PETSC_OWN_POINTER, &patch->dofs);
1727: ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArray, PETSC_OWN_POINTER, &patch->offs);
1728: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1729: ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithArtificial, globalDofsArrayWithArtificial, PETSC_OWN_POINTER, &patch->gtolWithArtificial);
1730: ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithArtificial, PETSC_OWN_POINTER, &patch->dofsWithArtificial);
1731: ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArrayWithArtificial, PETSC_OWN_POINTER, &patch->offsWithArtificial);
1732: }
1733: if (isNonlinear) {
1734: ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithAll, globalDofsArrayWithAll, PETSC_OWN_POINTER, &patch->gtolWithAll);
1735: ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithAll, PETSC_OWN_POINTER, &patch->dofsWithAll);
1736: ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArrayWithAll, PETSC_OWN_POINTER, &patch->offsWithAll);
1737: }
1738: return(0);
1739: }
1741: static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat, PetscBool withArtificial)
1742: {
1743: PC_PATCH *patch = (PC_PATCH *) pc->data;
1744: PetscBool flg;
1745: PetscInt csize, rsize;
1746: const char *prefix = NULL;
1750: if (withArtificial) {
1751: /* would be nice if we could create a rectangular matrix of size numDofsWithArtificial x numDofs here */
1752: PetscInt pStart;
1753: PetscSectionGetChart(patch->gtolCountsWithArtificial, &pStart, NULL);
1754: PetscSectionGetDof(patch->gtolCountsWithArtificial, point + pStart, &rsize);
1755: csize = rsize;
1756: } else {
1757: PetscInt pStart;
1758: PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);
1759: PetscSectionGetDof(patch->gtolCounts, point + pStart, &rsize);
1760: csize = rsize;
1761: }
1763: MatCreate(PETSC_COMM_SELF, mat);
1764: PCGetOptionsPrefix(pc, &prefix);
1765: MatSetOptionsPrefix(*mat, prefix);
1766: MatAppendOptionsPrefix(*mat, "pc_patch_sub_");
1767: if (patch->sub_mat_type) {MatSetType(*mat, patch->sub_mat_type);}
1768: else if (!patch->sub_mat_type) {MatSetType(*mat, MATDENSE);}
1769: MatSetSizes(*mat, rsize, csize, rsize, csize);
1770: PetscObjectTypeCompare((PetscObject) *mat, MATDENSE, &flg);
1771: if (!flg) {PetscObjectTypeCompare((PetscObject)*mat, MATSEQDENSE, &flg);}
1772: /* Sparse patch matrices */
1773: if (!flg) {
1774: PetscBT bt;
1775: PetscInt *dnnz = NULL;
1776: const PetscInt *dofsArray = NULL;
1777: PetscInt pStart, pEnd, ncell, offset, c, i, j;
1779: if (withArtificial) {
1780: ISGetIndices(patch->dofsWithArtificial, &dofsArray);
1781: } else {
1782: ISGetIndices(patch->dofs, &dofsArray);
1783: }
1784: PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);
1785: point += pStart;
1786: if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);
1787: PetscSectionGetDof(patch->cellCounts, point, &ncell);
1788: PetscSectionGetOffset(patch->cellCounts, point, &offset);
1789: PetscLogEventBegin(PC_Patch_Prealloc, pc, 0, 0, 0);
1790: /* A PetscBT uses N^2 bits to store the sparsity pattern on a
1791: * patch. This is probably OK if the patches are not too big,
1792: * but uses too much memory. We therefore switch based on rsize. */
1793: if (rsize < 3000) { /* FIXME: I picked this switch value out of my hat */
1794: PetscScalar *zeroes;
1795: PetscInt rows;
1797: PetscCalloc1(rsize, &dnnz);
1798: PetscBTCreate(rsize*rsize, &bt);
1799: for (c = 0; c < ncell; ++c) {
1800: const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell;
1801: for (i = 0; i < patch->totalDofsPerCell; ++i) {
1802: const PetscInt row = idx[i];
1803: if (row < 0) continue;
1804: for (j = 0; j < patch->totalDofsPerCell; ++j) {
1805: const PetscInt col = idx[j];
1806: const PetscInt key = row*rsize + col;
1807: if (col < 0) continue;
1808: if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1809: }
1810: }
1811: }
1813: if (patch->usercomputeopintfacet) {
1814: const PetscInt *intFacetsArray = NULL;
1815: PetscInt i, numIntFacets, intFacetOffset;
1816: const PetscInt *facetCells = NULL;
1818: PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1819: PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1820: ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1821: ISGetIndices(patch->intFacets, &intFacetsArray);
1822: for (i = 0; i < numIntFacets; i++) {
1823: const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0];
1824: const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1];
1825: PetscInt celli, cellj;
1827: for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1828: const PetscInt row = dofsArray[(offset + cell0)*patch->totalDofsPerCell + celli];
1829: if (row < 0) continue;
1830: for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1831: const PetscInt col = dofsArray[(offset + cell1)*patch->totalDofsPerCell + cellj];
1832: const PetscInt key = row*rsize + col;
1833: if (col < 0) continue;
1834: if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1835: }
1836: }
1838: for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1839: const PetscInt row = dofsArray[(offset + cell1)*patch->totalDofsPerCell + celli];
1840: if (row < 0) continue;
1841: for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1842: const PetscInt col = dofsArray[(offset + cell0)*patch->totalDofsPerCell + cellj];
1843: const PetscInt key = row*rsize + col;
1844: if (col < 0) continue;
1845: if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1846: }
1847: }
1848: }
1849: }
1850: PetscBTDestroy(&bt);
1851: MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL);
1852: PetscFree(dnnz);
1854: PetscCalloc1(patch->totalDofsPerCell*patch->totalDofsPerCell, &zeroes);
1855: for (c = 0; c < ncell; ++c) {
1856: const PetscInt *idx = &dofsArray[(offset + c)*patch->totalDofsPerCell];
1857: MatSetValues(*mat, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, zeroes, INSERT_VALUES);
1858: }
1859: MatGetLocalSize(*mat, &rows, NULL);
1860: for (i = 0; i < rows; ++i) {
1861: MatSetValues(*mat, 1, &i, 1, &i, zeroes, INSERT_VALUES);
1862: }
1864: if (patch->usercomputeopintfacet) {
1865: const PetscInt *intFacetsArray = NULL;
1866: PetscInt i, numIntFacets, intFacetOffset;
1867: const PetscInt *facetCells = NULL;
1869: PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1870: PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1871: ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1872: ISGetIndices(patch->intFacets, &intFacetsArray);
1873: for (i = 0; i < numIntFacets; i++) {
1874: const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0];
1875: const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1];
1876: const PetscInt *cell0idx = &dofsArray[(offset + cell0)*patch->totalDofsPerCell];
1877: const PetscInt *cell1idx = &dofsArray[(offset + cell1)*patch->totalDofsPerCell];
1878: MatSetValues(*mat, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, zeroes, INSERT_VALUES);
1879: MatSetValues(*mat, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, zeroes, INSERT_VALUES);
1880: }
1881: }
1883: MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
1884: MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);
1886: PetscFree(zeroes);
1888: } else { /* rsize too big, use MATPREALLOCATOR */
1889: Mat preallocator;
1890: PetscScalar* vals;
1892: PetscCalloc1(patch->totalDofsPerCell*patch->totalDofsPerCell, &vals);
1893: MatCreate(PETSC_COMM_SELF, &preallocator);
1894: MatSetType(preallocator, MATPREALLOCATOR);
1895: MatSetSizes(preallocator, rsize, rsize, rsize, rsize);
1896: MatSetUp(preallocator);
1898: for (c = 0; c < ncell; ++c) {
1899: const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell;
1900: MatSetValues(preallocator, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, vals, INSERT_VALUES);
1901: }
1903: if (patch->usercomputeopintfacet) {
1904: const PetscInt *intFacetsArray = NULL;
1905: PetscInt i, numIntFacets, intFacetOffset;
1906: const PetscInt *facetCells = NULL;
1908: PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1909: PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1910: ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1911: ISGetIndices(patch->intFacets, &intFacetsArray);
1912: for (i = 0; i < numIntFacets; i++) {
1913: const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0];
1914: const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1];
1915: const PetscInt *cell0idx = &dofsArray[(offset + cell0)*patch->totalDofsPerCell];
1916: const PetscInt *cell1idx = &dofsArray[(offset + cell1)*patch->totalDofsPerCell];
1917: MatSetValues(preallocator, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, vals, INSERT_VALUES);
1918: MatSetValues(preallocator, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, vals, INSERT_VALUES);
1919: }
1920: }
1922: PetscFree(vals);
1923: MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1924: MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1925: MatPreallocatorPreallocate(preallocator, PETSC_TRUE, *mat);
1926: MatDestroy(&preallocator);
1927: }
1928: PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0);
1929: if (withArtificial) {
1930: ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);
1931: } else {
1932: ISRestoreIndices(patch->dofs, &dofsArray);
1933: }
1934: }
1935: MatSetUp(*mat);
1936: return(0);
1937: }
1939: static PetscErrorCode PCPatchComputeFunction_DMPlex_Private(PC pc, PetscInt patchNum, Vec x, Vec F, IS cellIS, PetscInt n, const PetscInt *l2p, const PetscInt *l2pWithAll, void *ctx)
1940: {
1941: PC_PATCH *patch = (PC_PATCH *) pc->data;
1942: DM dm, plex;
1943: PetscSection s;
1944: const PetscInt *parray, *oarray;
1945: PetscInt Nf = patch->nsubspaces, Np, poff, p, f;
1946: PetscErrorCode ierr;
1949: if (patch->precomputeElementTensors) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Precomputing element tensors not implemented with DMPlex compute operator");
1950: PCGetDM(pc, &dm);
1951: DMConvert(dm, DMPLEX, &plex);
1952: dm = plex;
1953: DMGetLocalSection(dm, &s);
1954: /* Set offset into patch */
1955: PetscSectionGetDof(patch->pointCounts, patchNum, &Np);
1956: PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);
1957: ISGetIndices(patch->points, &parray);
1958: ISGetIndices(patch->offs, &oarray);
1959: for (f = 0; f < Nf; ++f) {
1960: for (p = 0; p < Np; ++p) {
1961: const PetscInt point = parray[poff+p];
1962: PetscInt dof;
1964: PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);
1965: PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);
1966: if (patch->nsubspaces == 1) {PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);}
1967: else {PetscSectionSetOffset(patch->patchSection, point, -1);}
1968: }
1969: }
1970: ISRestoreIndices(patch->points, &parray);
1971: ISRestoreIndices(patch->offs, &oarray);
1972: if (patch->viewSection) {ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);}
1973: DMPlexComputeResidual_Patch_Internal(dm, patch->patchSection, cellIS, 0.0, x, NULL, F, ctx);
1974: DMDestroy(&dm);
1975: return(0);
1976: }
1978: PetscErrorCode PCPatchComputeFunction_Internal(PC pc, Vec x, Vec F, PetscInt point)
1979: {
1980: PC_PATCH *patch = (PC_PATCH *) pc->data;
1981: const PetscInt *dofsArray;
1982: const PetscInt *dofsArrayWithAll;
1983: const PetscInt *cellsArray;
1984: PetscInt ncell, offset, pStart, pEnd;
1985: PetscErrorCode ierr;
1988: PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);
1989: if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n");
1990: ISGetIndices(patch->dofs, &dofsArray);
1991: ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);
1992: ISGetIndices(patch->cells, &cellsArray);
1993: PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);
1995: point += pStart;
1996: if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);
1998: PetscSectionGetDof(patch->cellCounts, point, &ncell);
1999: PetscSectionGetOffset(patch->cellCounts, point, &offset);
2000: if (ncell <= 0) {
2001: PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2002: return(0);
2003: }
2004: VecSet(F, 0.0);
2005: PetscStackPush("PCPatch user callback");
2006: /* Cannot reuse the same IS because the geometry info is being cached in it */
2007: ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);
2008: patch->usercomputef(pc, point, x, F, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell,
2009: dofsArrayWithAll + offset*patch->totalDofsPerCell,
2010: patch->usercomputefctx);
2011: PetscStackPop;
2012: ISDestroy(&patch->cellIS);
2013: ISRestoreIndices(patch->dofs, &dofsArray);
2014: ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);
2015: ISRestoreIndices(patch->cells, &cellsArray);
2016: if (patch->viewMatrix) {
2017: char name[PETSC_MAX_PATH_LEN];
2019: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch vector for Point %D", point);
2020: PetscObjectSetName((PetscObject) F, name);
2021: ObjectView((PetscObject) F, patch->viewerMatrix, patch->formatMatrix);
2022: }
2023: PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2024: return(0);
2025: }
2027: static PetscErrorCode PCPatchComputeOperator_DMPlex_Private(PC pc, PetscInt patchNum, Vec x, Mat J, IS cellIS, PetscInt n, const PetscInt *l2p, const PetscInt *l2pWithAll, void *ctx)
2028: {
2029: PC_PATCH *patch = (PC_PATCH *) pc->data;
2030: DM dm, plex;
2031: PetscSection s;
2032: const PetscInt *parray, *oarray;
2033: PetscInt Nf = patch->nsubspaces, Np, poff, p, f;
2034: PetscErrorCode ierr;
2037: PCGetDM(pc, &dm);
2038: DMConvert(dm, DMPLEX, &plex);
2039: dm = plex;
2040: DMGetLocalSection(dm, &s);
2041: /* Set offset into patch */
2042: PetscSectionGetDof(patch->pointCounts, patchNum, &Np);
2043: PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);
2044: ISGetIndices(patch->points, &parray);
2045: ISGetIndices(patch->offs, &oarray);
2046: for (f = 0; f < Nf; ++f) {
2047: for (p = 0; p < Np; ++p) {
2048: const PetscInt point = parray[poff+p];
2049: PetscInt dof;
2051: PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);
2052: PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);
2053: if (patch->nsubspaces == 1) {PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);}
2054: else {PetscSectionSetOffset(patch->patchSection, point, -1);}
2055: }
2056: }
2057: ISRestoreIndices(patch->points, &parray);
2058: ISRestoreIndices(patch->offs, &oarray);
2059: if (patch->viewSection) {ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);}
2060: /* TODO Shut off MatViewFromOptions() in MatAssemblyEnd() here */
2061: DMPlexComputeJacobian_Patch_Internal(dm, patch->patchSection, patch->patchSection, cellIS, 0.0, 0.0, x, NULL, J, J, ctx);
2062: DMDestroy(&dm);
2063: return(0);
2064: }
2066: /* This function zeros mat on entry */
2067: PetscErrorCode PCPatchComputeOperator_Internal(PC pc, Vec x, Mat mat, PetscInt point, PetscBool withArtificial)
2068: {
2069: PC_PATCH *patch = (PC_PATCH *) pc->data;
2070: const PetscInt *dofsArray;
2071: const PetscInt *dofsArrayWithAll = NULL;
2072: const PetscInt *cellsArray;
2073: PetscInt ncell, offset, pStart, pEnd, numIntFacets, intFacetOffset;
2074: PetscBool isNonlinear;
2075: PetscErrorCode ierr;
2078: PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);
2079: isNonlinear = patch->isNonlinear;
2080: if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n");
2081: if (withArtificial) {
2082: ISGetIndices(patch->dofsWithArtificial, &dofsArray);
2083: } else {
2084: ISGetIndices(patch->dofs, &dofsArray);
2085: }
2086: if (isNonlinear) {
2087: ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);
2088: }
2089: ISGetIndices(patch->cells, &cellsArray);
2090: PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);
2092: point += pStart;
2093: if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);
2095: PetscSectionGetDof(patch->cellCounts, point, &ncell);
2096: PetscSectionGetOffset(patch->cellCounts, point, &offset);
2097: if (ncell <= 0) {
2098: PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2099: return(0);
2100: }
2101: MatZeroEntries(mat);
2102: if (patch->precomputeElementTensors) {
2103: PetscInt i;
2104: PetscInt ndof = patch->totalDofsPerCell;
2105: const PetscScalar *elementTensors;
2107: VecGetArrayRead(patch->cellMats, &elementTensors);
2108: for (i = 0; i < ncell; i++) {
2109: const PetscInt cell = cellsArray[i + offset];
2110: const PetscInt *idx = dofsArray + (offset + i)*ndof;
2111: const PetscScalar *v = elementTensors + patch->precomputedTensorLocations[cell]*ndof*ndof;
2112: MatSetValues(mat, ndof, idx, ndof, idx, v, ADD_VALUES);
2113: }
2114: VecRestoreArrayRead(patch->cellMats, &elementTensors);
2115: MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
2116: MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
2117: } else {
2118: PetscStackPush("PCPatch user callback");
2119: /* Cannot reuse the same IS because the geometry info is being cached in it */
2120: ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);
2121: patch->usercomputeop(pc, point, x, mat, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell, dofsArrayWithAll ? dofsArrayWithAll + offset*patch->totalDofsPerCell : NULL, patch->usercomputeopctx);
2122: }
2123: if (patch->usercomputeopintfacet) {
2124: PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
2125: PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
2126: if (numIntFacets > 0) {
2127: /* For each interior facet, grab the two cells (in local numbering, and concatenate dof numberings for those cells) */
2128: PetscInt *facetDofs = NULL, *facetDofsWithAll = NULL;
2129: const PetscInt *intFacetsArray = NULL;
2130: PetscInt idx = 0;
2131: PetscInt i, c, d;
2132: PetscInt fStart;
2133: DM dm, plex;
2134: IS facetIS = NULL;
2135: const PetscInt *facetCells = NULL;
2137: ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
2138: ISGetIndices(patch->intFacets, &intFacetsArray);
2139: PCGetDM(pc, &dm);
2140: DMConvert(dm, DMPLEX, &plex);
2141: dm = plex;
2142: DMPlexGetHeightStratum(dm, 1, &fStart, NULL);
2143: /* FIXME: Pull this malloc out. */
2144: PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofs);
2145: if (dofsArrayWithAll) {
2146: PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofsWithAll);
2147: }
2148: if (patch->precomputeElementTensors) {
2149: PetscInt nFacetDof = 2*patch->totalDofsPerCell;
2150: const PetscScalar *elementTensors;
2152: VecGetArrayRead(patch->intFacetMats, &elementTensors);
2154: for (i = 0; i < numIntFacets; i++) {
2155: const PetscInt facet = intFacetsArray[i + intFacetOffset];
2156: const PetscScalar *v = elementTensors + patch->precomputedIntFacetTensorLocations[facet - fStart]*nFacetDof*nFacetDof;
2157: idx = 0;
2158: /*
2159: * 0--1
2160: * |\-|
2161: * |+\|
2162: * 2--3
2163: * [0, 2, 3, 0, 1, 3]
2164: */
2165: for (c = 0; c < 2; c++) {
2166: const PetscInt cell = facetCells[2*(intFacetOffset + i) + c];
2167: for (d = 0; d < patch->totalDofsPerCell; d++) {
2168: facetDofs[idx] = dofsArray[(offset + cell)*patch->totalDofsPerCell + d];
2169: idx++;
2170: }
2171: }
2172: MatSetValues(mat, nFacetDof, facetDofs, nFacetDof, facetDofs, v, ADD_VALUES);
2173: }
2174: VecRestoreArrayRead(patch->intFacetMats, &elementTensors);
2175: } else {
2176: /*
2177: * 0--1
2178: * |\-|
2179: * |+\|
2180: * 2--3
2181: * [0, 2, 3, 0, 1, 3]
2182: */
2183: for (i = 0; i < numIntFacets; i++) {
2184: for (c = 0; c < 2; c++) {
2185: const PetscInt cell = facetCells[2*(intFacetOffset + i) + c];
2186: for (d = 0; d < patch->totalDofsPerCell; d++) {
2187: facetDofs[idx] = dofsArray[(offset + cell)*patch->totalDofsPerCell + d];
2188: if (dofsArrayWithAll) {
2189: facetDofsWithAll[idx] = dofsArrayWithAll[(offset + cell)*patch->totalDofsPerCell + d];
2190: }
2191: idx++;
2192: }
2193: }
2194: }
2195: ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray + intFacetOffset, PETSC_USE_POINTER, &facetIS);
2196: patch->usercomputeopintfacet(pc, point, x, mat, facetIS, 2*numIntFacets*patch->totalDofsPerCell, facetDofs, facetDofsWithAll, patch->usercomputeopintfacetctx);
2197: ISDestroy(&facetIS);
2198: }
2199: ISRestoreIndices(patch->intFacetsToPatchCell, &facetCells);
2200: ISRestoreIndices(patch->intFacets, &intFacetsArray);
2201: PetscFree(facetDofs);
2202: PetscFree(facetDofsWithAll);
2203: DMDestroy(&dm);
2204: }
2205: }
2207: MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
2208: MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
2210: if (!(withArtificial || isNonlinear) && patch->denseinverse) {
2211: MatFactorInfo info;
2212: PetscBool flg;
2213: PetscObjectTypeCompare((PetscObject)mat, MATSEQDENSE, &flg);
2214: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Invalid Mat type for dense inverse");
2215: MatFactorInfoInitialize(&info);
2216: MatLUFactor(mat, NULL, NULL, &info);
2217: MatSeqDenseInvertFactors_Private(mat);
2218: }
2219: PetscStackPop;
2220: ISDestroy(&patch->cellIS);
2221: if (withArtificial) {
2222: ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);
2223: } else {
2224: ISRestoreIndices(patch->dofs, &dofsArray);
2225: }
2226: if (isNonlinear) {
2227: ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);
2228: }
2229: ISRestoreIndices(patch->cells, &cellsArray);
2230: if (patch->viewMatrix) {
2231: char name[PETSC_MAX_PATH_LEN];
2233: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch matrix for Point %D", point);
2234: PetscObjectSetName((PetscObject) mat, name);
2235: ObjectView((PetscObject) mat, patch->viewerMatrix, patch->formatMatrix);
2236: }
2237: PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2238: return(0);
2239: }
2241: static PetscErrorCode MatSetValues_PCPatch_Private(Mat mat, PetscInt m, const PetscInt idxm[],
2242: PetscInt n, const PetscInt idxn[], const PetscScalar *v, InsertMode addv)
2243: {
2244: Vec data;
2245: PetscScalar *array;
2246: PetscInt bs, nz, i, j, cell;
2249: MatShellGetContext(mat, &data);
2250: VecGetBlockSize(data, &bs);
2251: VecGetSize(data, &nz);
2252: VecGetArray(data, &array);
2253: if (m != n) SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Only for square insertion");
2254: cell = (PetscInt)(idxm[0]/bs); /* use the fact that this is called once per cell */
2255: for (i = 0; i < m; i++) {
2256: if (idxm[i] != idxn[i]) SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Row and column indices must match!");
2257: for (j = 0; j < n; j++) {
2258: const PetscScalar v_ = v[i*bs + j];
2259: /* Indexing is special to the data structure we have! */
2260: if (addv == INSERT_VALUES) {
2261: array[cell*bs*bs + i*bs + j] = v_;
2262: } else {
2263: array[cell*bs*bs + i*bs + j] += v_;
2264: }
2265: }
2266: }
2267: VecRestoreArray(data, &array);
2268: return(0);
2269: }
2271: static PetscErrorCode PCPatchPrecomputePatchTensors_Private(PC pc)
2272: {
2273: PC_PATCH *patch = (PC_PATCH *)pc->data;
2274: const PetscInt *cellsArray;
2275: PetscInt ncell, offset;
2276: const PetscInt *dofMapArray;
2277: PetscInt i, j;
2278: IS dofMap;
2279: IS cellIS;
2280: const PetscInt ndof = patch->totalDofsPerCell;
2281: PetscErrorCode ierr;
2282: Mat vecMat;
2283: PetscInt cStart, cEnd;
2284: DM dm, plex;
2286: ISGetSize(patch->cells, &ncell);
2287: if (!ncell) { /* No cells to assemble over -> skip */
2288: return(0);
2289: }
2291: PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);
2293: PCGetDM(pc, &dm);
2294: DMConvert(dm, DMPLEX, &plex);
2295: dm = plex;
2296: if (!patch->allCells) {
2297: PetscHSetI cells;
2298: PetscHashIter hi;
2299: PetscInt pStart, pEnd;
2300: PetscInt *allCells = NULL;
2301: PetscHSetICreate(&cells);
2302: ISGetIndices(patch->cells, &cellsArray);
2303: PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);
2304: for (i = pStart; i < pEnd; i++) {
2305: PetscSectionGetDof(patch->cellCounts, i, &ncell);
2306: PetscSectionGetOffset(patch->cellCounts, i, &offset);
2307: if (ncell <= 0) continue;
2308: for (j = 0; j < ncell; j++) {
2309: PetscHSetIAdd(cells, cellsArray[offset + j]);
2310: }
2311: }
2312: ISRestoreIndices(patch->cells, &cellsArray);
2313: PetscHSetIGetSize(cells, &ncell);
2314: PetscMalloc1(ncell, &allCells);
2315: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2316: PetscMalloc1(cEnd-cStart, &patch->precomputedTensorLocations);
2317: i = 0;
2318: PetscHashIterBegin(cells, hi);
2319: while (!PetscHashIterAtEnd(cells, hi)) {
2320: PetscHashIterGetKey(cells, hi, allCells[i]);
2321: patch->precomputedTensorLocations[allCells[i]] = i;
2322: PetscHashIterNext(cells, hi);
2323: i++;
2324: }
2325: PetscHSetIDestroy(&cells);
2326: ISCreateGeneral(PETSC_COMM_SELF, ncell, allCells, PETSC_OWN_POINTER, &patch->allCells);
2327: }
2328: ISGetSize(patch->allCells, &ncell);
2329: if (!patch->cellMats) {
2330: VecCreateSeq(PETSC_COMM_SELF, ncell*ndof*ndof, &patch->cellMats);
2331: VecSetBlockSize(patch->cellMats, ndof);
2332: }
2333: VecSet(patch->cellMats, 0);
2335: MatCreateShell(PETSC_COMM_SELF, ncell*ndof, ncell*ndof, ncell*ndof, ncell*ndof,
2336: (void*)patch->cellMats, &vecMat);
2337: MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void(*)(void))&MatSetValues_PCPatch_Private);
2338: ISGetSize(patch->allCells, &ncell);
2339: ISCreateStride(PETSC_COMM_SELF, ndof*ncell, 0, 1, &dofMap);
2340: ISGetIndices(dofMap, &dofMapArray);
2341: ISGetIndices(patch->allCells, &cellsArray);
2342: ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray, PETSC_USE_POINTER, &cellIS);
2343: PetscStackPush("PCPatch user callback");
2344: /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2345: patch->usercomputeop(pc, -1, NULL, vecMat, cellIS, ndof*ncell, dofMapArray, NULL, patch->usercomputeopctx);
2346: PetscStackPop;
2347: ISDestroy(&cellIS);
2348: MatDestroy(&vecMat);
2349: ISRestoreIndices(patch->allCells, &cellsArray);
2350: ISRestoreIndices(dofMap, &dofMapArray);
2351: ISDestroy(&dofMap);
2353: if (patch->usercomputeopintfacet) {
2354: PetscInt nIntFacets;
2355: IS intFacetsIS;
2356: const PetscInt *intFacetsArray = NULL;
2357: if (!patch->allIntFacets) {
2358: PetscHSetI facets;
2359: PetscHashIter hi;
2360: PetscInt pStart, pEnd, fStart, fEnd;
2361: PetscInt *allIntFacets = NULL;
2362: PetscHSetICreate(&facets);
2363: ISGetIndices(patch->intFacets, &intFacetsArray);
2364: PetscSectionGetChart(patch->intFacetCounts, &pStart, &pEnd);
2365: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2366: for (i = pStart; i < pEnd; i++) {
2367: PetscSectionGetDof(patch->intFacetCounts, i, &nIntFacets);
2368: PetscSectionGetOffset(patch->intFacetCounts, i, &offset);
2369: if (nIntFacets <= 0) continue;
2370: for (j = 0; j < nIntFacets; j++) {
2371: PetscHSetIAdd(facets, intFacetsArray[offset + j]);
2372: }
2373: }
2374: ISRestoreIndices(patch->intFacets, &intFacetsArray);
2375: PetscHSetIGetSize(facets, &nIntFacets);
2376: PetscMalloc1(nIntFacets, &allIntFacets);
2377: PetscMalloc1(fEnd-fStart, &patch->precomputedIntFacetTensorLocations);
2378: i = 0;
2379: PetscHashIterBegin(facets, hi);
2380: while (!PetscHashIterAtEnd(facets, hi)) {
2381: PetscHashIterGetKey(facets, hi, allIntFacets[i]);
2382: patch->precomputedIntFacetTensorLocations[allIntFacets[i] - fStart] = i;
2383: PetscHashIterNext(facets, hi);
2384: i++;
2385: }
2386: PetscHSetIDestroy(&facets);
2387: ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, allIntFacets, PETSC_OWN_POINTER, &patch->allIntFacets);
2388: }
2389: ISGetSize(patch->allIntFacets, &nIntFacets);
2390: if (!patch->intFacetMats) {
2391: VecCreateSeq(PETSC_COMM_SELF, nIntFacets*ndof*ndof*4, &patch->intFacetMats);
2392: VecSetBlockSize(patch->intFacetMats, ndof*2);
2393: }
2394: VecSet(patch->intFacetMats, 0);
2396: MatCreateShell(PETSC_COMM_SELF, nIntFacets*ndof*2, nIntFacets*ndof*2, nIntFacets*ndof*2, nIntFacets*ndof*2,
2397: (void*)patch->intFacetMats, &vecMat);
2398: MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void(*)(void))&MatSetValues_PCPatch_Private);
2399: ISCreateStride(PETSC_COMM_SELF, 2*ndof*nIntFacets, 0, 1, &dofMap);
2400: ISGetIndices(dofMap, &dofMapArray);
2401: ISGetIndices(patch->allIntFacets, &intFacetsArray);
2402: ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, intFacetsArray, PETSC_USE_POINTER, &intFacetsIS);
2403: PetscStackPush("PCPatch user callback (interior facets)");
2404: /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2405: patch->usercomputeopintfacet(pc, -1, NULL, vecMat, intFacetsIS, 2*ndof*nIntFacets, dofMapArray, NULL, patch->usercomputeopintfacetctx);
2406: PetscStackPop;
2407: ISDestroy(&intFacetsIS);
2408: MatDestroy(&vecMat);
2409: ISRestoreIndices(patch->allIntFacets, &intFacetsArray);
2410: ISRestoreIndices(dofMap, &dofMapArray);
2411: ISDestroy(&dofMap);
2412: }
2413: DMDestroy(&dm);
2414: PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2416: return(0);
2417: }
2419: PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat, PatchScatterType scattertype)
2420: {
2421: PC_PATCH *patch = (PC_PATCH *) pc->data;
2422: const PetscScalar *xArray = NULL;
2423: PetscScalar *yArray = NULL;
2424: const PetscInt *gtolArray = NULL;
2425: PetscInt dof, offset, lidx;
2426: PetscErrorCode ierr;
2429: VecGetArrayRead(x, &xArray);
2430: VecGetArray(y, &yArray);
2431: if (scattertype == SCATTER_WITHARTIFICIAL) {
2432: PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);
2433: PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offset);
2434: ISGetIndices(patch->gtolWithArtificial, >olArray);
2435: } else if (scattertype == SCATTER_WITHALL) {
2436: PetscSectionGetDof(patch->gtolCountsWithAll, p, &dof);
2437: PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offset);
2438: ISGetIndices(patch->gtolWithAll, >olArray);
2439: } else {
2440: PetscSectionGetDof(patch->gtolCounts, p, &dof);
2441: PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2442: ISGetIndices(patch->gtol, >olArray);
2443: }
2444: if (mode == INSERT_VALUES && scat != SCATTER_FORWARD) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't insert if not scattering forward\n");
2445: if (mode == ADD_VALUES && scat != SCATTER_REVERSE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't add if not scattering reverse\n");
2446: for (lidx = 0; lidx < dof; ++lidx) {
2447: const PetscInt gidx = gtolArray[offset+lidx];
2449: if (mode == INSERT_VALUES) yArray[lidx] = xArray[gidx]; /* Forward */
2450: else yArray[gidx] += xArray[lidx]; /* Reverse */
2451: }
2452: if (scattertype == SCATTER_WITHARTIFICIAL) {
2453: ISRestoreIndices(patch->gtolWithArtificial, >olArray);
2454: } else if (scattertype == SCATTER_WITHALL) {
2455: ISRestoreIndices(patch->gtolWithAll, >olArray);
2456: } else {
2457: ISRestoreIndices(patch->gtol, >olArray);
2458: }
2459: VecRestoreArrayRead(x, &xArray);
2460: VecRestoreArray(y, &yArray);
2461: return(0);
2462: }
2464: static PetscErrorCode PCSetUp_PATCH_Linear(PC pc)
2465: {
2466: PC_PATCH *patch = (PC_PATCH *) pc->data;
2467: const char *prefix;
2468: PetscInt i;
2472: if (!pc->setupcalled) {
2473: if (!patch->save_operators && patch->denseinverse) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Can't have dense inverse without save operators");
2474: if (!patch->denseinverse) {
2475: PetscMalloc1(patch->npatch, &patch->solver);
2476: PCGetOptionsPrefix(pc, &prefix);
2477: for (i = 0; i < patch->npatch; ++i) {
2478: KSP ksp;
2479: PC subpc;
2481: KSPCreate(PETSC_COMM_SELF, &ksp);
2482: KSPSetErrorIfNotConverged(ksp, pc->erroriffailure);
2483: KSPSetOptionsPrefix(ksp, prefix);
2484: KSPAppendOptionsPrefix(ksp, "sub_");
2485: PetscObjectIncrementTabLevel((PetscObject) ksp, (PetscObject) pc, 1);
2486: KSPGetPC(ksp, &subpc);
2487: PetscObjectIncrementTabLevel((PetscObject) subpc, (PetscObject) pc, 1);
2488: PetscLogObjectParent((PetscObject) pc, (PetscObject) ksp);
2489: patch->solver[i] = (PetscObject) ksp;
2490: }
2491: }
2492: }
2493: if (patch->save_operators) {
2494: if (patch->precomputeElementTensors) {
2495: PCPatchPrecomputePatchTensors_Private(pc);
2496: }
2497: for (i = 0; i < patch->npatch; ++i) {
2498: PCPatchComputeOperator_Internal(pc, NULL, patch->mat[i], i, PETSC_FALSE);
2499: if (!patch->denseinverse) {
2500: KSPSetOperators((KSP) patch->solver[i], patch->mat[i], patch->mat[i]);
2501: } else if (patch->mat[i] && !patch->densesolve) {
2502: /* Setup matmult callback */
2503: MatGetOperation(patch->mat[i], MATOP_MULT, (void (**)(void))&patch->densesolve);
2504: }
2505: }
2506: }
2507: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2508: for (i = 0; i < patch->npatch; ++i) {
2509: /* Instead of padding patch->patchUpdate with zeros to get */
2510: /* patch->patchUpdateWithArtificial and then multiplying with the matrix, */
2511: /* just get rid of the columns that correspond to the dofs with */
2512: /* artificial bcs. That's of course fairly inefficient, hopefully we */
2513: /* can just assemble the rectangular matrix in the first place. */
2514: Mat matSquare;
2515: IS rowis;
2516: PetscInt dof;
2518: MatGetSize(patch->mat[i], &dof, NULL);
2519: if (dof == 0) {
2520: patch->matWithArtificial[i] = NULL;
2521: continue;
2522: }
2524: PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);
2525: PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);
2527: MatGetSize(matSquare, &dof, NULL);
2528: ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis);
2529: if (pc->setupcalled) {
2530: MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_REUSE_MATRIX, &patch->matWithArtificial[i]);
2531: } else {
2532: MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &patch->matWithArtificial[i]);
2533: }
2534: ISDestroy(&rowis);
2535: MatDestroy(&matSquare);
2536: }
2537: }
2538: return(0);
2539: }
2541: static PetscErrorCode PCSetUp_PATCH(PC pc)
2542: {
2543: PC_PATCH *patch = (PC_PATCH *) pc->data;
2544: PetscInt i;
2545: PetscBool isNonlinear;
2546: PetscInt maxDof = -1, maxDofWithArtificial = -1;
2550: if (!pc->setupcalled) {
2551: PetscInt pStart, pEnd, p;
2552: PetscInt localSize;
2554: PetscLogEventBegin(PC_Patch_CreatePatches, pc, 0, 0, 0);
2556: isNonlinear = patch->isNonlinear;
2557: if (!patch->nsubspaces) {
2558: DM dm, plex;
2559: PetscSection s;
2560: PetscInt cStart, cEnd, c, Nf, f, numGlobalBcs = 0, *globalBcs, *Nb, **cellDofs;
2562: PCGetDM(pc, &dm);
2563: if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "Must set DM for PCPATCH or call PCPatchSetDiscretisationInfo()");
2564: DMConvert(dm, DMPLEX, &plex);
2565: dm = plex;
2566: DMGetLocalSection(dm, &s);
2567: PetscSectionGetNumFields(s, &Nf);
2568: PetscSectionGetChart(s, &pStart, &pEnd);
2569: for (p = pStart; p < pEnd; ++p) {
2570: PetscInt cdof;
2571: PetscSectionGetConstraintDof(s, p, &cdof);
2572: numGlobalBcs += cdof;
2573: }
2574: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2575: PetscMalloc3(Nf, &Nb, Nf, &cellDofs, numGlobalBcs, &globalBcs);
2576: for (f = 0; f < Nf; ++f) {
2577: PetscFE fe;
2578: PetscDualSpace sp;
2579: PetscInt cdoff = 0;
2581: DMGetField(dm, f, NULL, (PetscObject *) &fe);
2582: /* PetscFEGetNumComponents(fe, &Nc[f]); */
2583: PetscFEGetDualSpace(fe, &sp);
2584: PetscDualSpaceGetDimension(sp, &Nb[f]);
2586: PetscMalloc1((cEnd-cStart)*Nb[f], &cellDofs[f]);
2587: for (c = cStart; c < cEnd; ++c) {
2588: PetscInt *closure = NULL;
2589: PetscInt clSize = 0, cl;
2591: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
2592: for (cl = 0; cl < clSize*2; cl += 2) {
2593: const PetscInt p = closure[cl];
2594: PetscInt fdof, d, foff;
2596: PetscSectionGetFieldDof(s, p, f, &fdof);
2597: PetscSectionGetFieldOffset(s, p, f, &foff);
2598: for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d;
2599: }
2600: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
2601: }
2602: if (cdoff != (cEnd-cStart)*Nb[f]) SETERRQ4(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_SIZ, "Total number of cellDofs %D for field %D should be Nc (%D) * cellDof (%D)", cdoff, f, cEnd-cStart, Nb[f]);
2603: }
2604: numGlobalBcs = 0;
2605: for (p = pStart; p < pEnd; ++p) {
2606: const PetscInt *ind;
2607: PetscInt off, cdof, d;
2609: PetscSectionGetOffset(s, p, &off);
2610: PetscSectionGetConstraintDof(s, p, &cdof);
2611: PetscSectionGetConstraintIndices(s, p, &ind);
2612: for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d];
2613: }
2615: PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **) cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs);
2616: for (f = 0; f < Nf; ++f) {
2617: PetscFree(cellDofs[f]);
2618: }
2619: PetscFree3(Nb, cellDofs, globalBcs);
2620: PCPatchSetComputeFunction(pc, PCPatchComputeFunction_DMPlex_Private, NULL);
2621: PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL);
2622: DMDestroy(&dm);
2623: }
2625: localSize = patch->subspaceOffsets[patch->nsubspaces];
2626: VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localRHS);
2627: VecSetUp(patch->localRHS);
2628: VecDuplicate(patch->localRHS, &patch->localUpdate);
2629: PCPatchCreateCellPatches(pc);
2630: PCPatchCreateCellPatchDiscretisationInfo(pc);
2632: /* OK, now build the work vectors */
2633: PetscSectionGetChart(patch->gtolCounts, &pStart, &pEnd);
2635: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2636: PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithArtificial);
2637: }
2638: if (isNonlinear) {
2639: PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithAll);
2640: }
2641: for (p = pStart; p < pEnd; ++p) {
2642: PetscInt dof;
2644: PetscSectionGetDof(patch->gtolCounts, p, &dof);
2645: maxDof = PetscMax(maxDof, dof);
2646: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2647: const PetscInt *gtolArray, *gtolArrayWithArtificial = NULL;
2648: PetscInt numPatchDofs, offset;
2649: PetscInt numPatchDofsWithArtificial, offsetWithArtificial;
2650: PetscInt dofWithoutArtificialCounter = 0;
2651: PetscInt *patchWithoutArtificialToWithArtificialArray;
2653: PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);
2654: maxDofWithArtificial = PetscMax(maxDofWithArtificial, dof);
2656: /* Now build the mapping that for a dof in a patch WITHOUT dofs that have artificial bcs gives the */
2657: /* the index in the patch with all dofs */
2658: ISGetIndices(patch->gtol, >olArray);
2660: PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);
2661: if (numPatchDofs == 0) {
2662: patch->dofMappingWithoutToWithArtificial[p-pStart] = NULL;
2663: continue;
2664: }
2666: PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2667: ISGetIndices(patch->gtolWithArtificial, >olArrayWithArtificial);
2668: PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &numPatchDofsWithArtificial);
2669: PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offsetWithArtificial);
2671: PetscMalloc1(numPatchDofs, &patchWithoutArtificialToWithArtificialArray);
2672: for (i=0; i<numPatchDofsWithArtificial; i++) {
2673: if (gtolArrayWithArtificial[i+offsetWithArtificial] == gtolArray[offset+dofWithoutArtificialCounter]) {
2674: patchWithoutArtificialToWithArtificialArray[dofWithoutArtificialCounter] = i;
2675: dofWithoutArtificialCounter++;
2676: if (dofWithoutArtificialCounter == numPatchDofs)
2677: break;
2678: }
2679: }
2680: ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutArtificialToWithArtificialArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithArtificial[p-pStart]);
2681: ISRestoreIndices(patch->gtol, >olArray);
2682: ISRestoreIndices(patch->gtolWithArtificial, >olArrayWithArtificial);
2683: }
2684: if (isNonlinear) {
2685: const PetscInt *gtolArray, *gtolArrayWithAll = NULL;
2686: PetscInt numPatchDofs, offset;
2687: PetscInt numPatchDofsWithAll, offsetWithAll;
2688: PetscInt dofWithoutAllCounter = 0;
2689: PetscInt *patchWithoutAllToWithAllArray;
2691: /* Now build the mapping that for a dof in a patch WITHOUT dofs that have artificial bcs gives the */
2692: /* the index in the patch with all dofs */
2693: ISGetIndices(patch->gtol, >olArray);
2695: PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);
2696: if (numPatchDofs == 0) {
2697: patch->dofMappingWithoutToWithAll[p-pStart] = NULL;
2698: continue;
2699: }
2701: PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2702: ISGetIndices(patch->gtolWithAll, >olArrayWithAll);
2703: PetscSectionGetDof(patch->gtolCountsWithAll, p, &numPatchDofsWithAll);
2704: PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offsetWithAll);
2706: PetscMalloc1(numPatchDofs, &patchWithoutAllToWithAllArray);
2708: for (i=0; i<numPatchDofsWithAll; i++) {
2709: if (gtolArrayWithAll[i+offsetWithAll] == gtolArray[offset+dofWithoutAllCounter]) {
2710: patchWithoutAllToWithAllArray[dofWithoutAllCounter] = i;
2711: dofWithoutAllCounter++;
2712: if (dofWithoutAllCounter == numPatchDofs)
2713: break;
2714: }
2715: }
2716: ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutAllToWithAllArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithAll[p-pStart]);
2717: ISRestoreIndices(patch->gtol, >olArray);
2718: ISRestoreIndices(patch->gtolWithAll, >olArrayWithAll);
2719: }
2720: }
2721: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2722: VecCreateSeq(PETSC_COMM_SELF, maxDofWithArtificial, &patch->patchRHSWithArtificial);
2723: VecSetUp(patch->patchRHSWithArtificial);
2724: }
2725: VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchRHS);
2726: VecSetUp(patch->patchRHS);
2727: VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchUpdate);
2728: VecSetUp(patch->patchUpdate);
2729: if (patch->save_operators) {
2730: PetscMalloc1(patch->npatch, &patch->mat);
2731: for (i = 0; i < patch->npatch; ++i) {
2732: PCPatchCreateMatrix_Private(pc, i, &patch->mat[i], PETSC_FALSE);
2733: }
2734: }
2735: PetscLogEventEnd(PC_Patch_CreatePatches, pc, 0, 0, 0);
2737: /* If desired, calculate weights for dof multiplicity */
2738: if (patch->partition_of_unity) {
2739: PetscScalar *input = NULL;
2740: PetscScalar *output = NULL;
2741: Vec global;
2743: VecDuplicate(patch->localRHS, &patch->dof_weights);
2744: if (patch->local_composition_type == PC_COMPOSITE_ADDITIVE) {
2745: for (i = 0; i < patch->npatch; ++i) {
2746: PetscInt dof;
2748: PetscSectionGetDof(patch->gtolCounts, i+pStart, &dof);
2749: if (dof <= 0) continue;
2750: VecSet(patch->patchRHS, 1.0);
2751: PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchRHS, patch->dof_weights, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);
2752: }
2753: } else {
2754: /* multiplicative is actually only locally multiplicative and globally additive. need the pou where the mesh decomposition overlaps */
2755: VecSet(patch->dof_weights, 1.0);
2756: }
2758: VecDuplicate(patch->dof_weights, &global);
2759: VecSet(global, 0.);
2761: VecGetArray(patch->dof_weights, &input);
2762: VecGetArray(global, &output);
2763: PetscSFReduceBegin(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM);
2764: PetscSFReduceEnd(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM);
2765: VecRestoreArray(patch->dof_weights, &input);
2766: VecRestoreArray(global, &output);
2768: VecReciprocal(global);
2770: VecGetArray(patch->dof_weights, &output);
2771: VecGetArray(global, &input);
2772: PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, input, output,MPI_REPLACE);
2773: PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, input, output,MPI_REPLACE);
2774: VecRestoreArray(patch->dof_weights, &output);
2775: VecRestoreArray(global, &input);
2776: VecDestroy(&global);
2777: }
2778: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE && patch->save_operators) {
2779: PetscMalloc1(patch->npatch, &patch->matWithArtificial);
2780: }
2781: }
2782: (*patch->setupsolver)(pc);
2783: return(0);
2784: }
2786: static PetscErrorCode PCApply_PATCH_Linear(PC pc, PetscInt i, Vec x, Vec y)
2787: {
2788: PC_PATCH *patch = (PC_PATCH *) pc->data;
2789: KSP ksp;
2790: Mat op;
2791: PetscInt m, n;
2795: if (patch->denseinverse) {
2796: (*patch->densesolve)(patch->mat[i], x, y);
2797: return(0);
2798: }
2799: ksp = (KSP) patch->solver[i];
2800: if (!patch->save_operators) {
2801: Mat mat;
2803: PCPatchCreateMatrix_Private(pc, i, &mat, PETSC_FALSE);
2804: /* Populate operator here. */
2805: PCPatchComputeOperator_Internal(pc, NULL, mat, i, PETSC_FALSE);
2806: KSPSetOperators(ksp, mat, mat);
2807: /* Drop reference so the KSPSetOperators below will blow it away. */
2808: MatDestroy(&mat);
2809: }
2810: PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);
2811: if (!ksp->setfromoptionscalled) {
2812: KSPSetFromOptions(ksp);
2813: }
2814: /* Disgusting trick to reuse work vectors */
2815: KSPGetOperators(ksp, &op, NULL);
2816: MatGetLocalSize(op, &m, &n);
2817: x->map->n = m;
2818: y->map->n = n;
2819: x->map->N = m;
2820: y->map->N = n;
2821: KSPSolve(ksp, x, y);
2822: KSPCheckSolve(ksp, pc, y);
2823: PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);
2824: if (!patch->save_operators) {
2825: PC pc;
2826: KSPSetOperators(ksp, NULL, NULL);
2827: KSPGetPC(ksp, &pc);
2828: /* Destroy PC context too, otherwise the factored matrix hangs around. */
2829: PCReset(pc);
2830: }
2831: return(0);
2832: }
2834: static PetscErrorCode PCUpdateMultiplicative_PATCH_Linear(PC pc, PetscInt i, PetscInt pStart)
2835: {
2836: PC_PATCH *patch = (PC_PATCH *) pc->data;
2837: Mat multMat;
2838: PetscInt n, m;
2843: if (patch->save_operators) {
2844: multMat = patch->matWithArtificial[i];
2845: } else {
2846: /*Very inefficient, hopefully we can just assemble the rectangular matrix in the first place.*/
2847: Mat matSquare;
2848: PetscInt dof;
2849: IS rowis;
2850: PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);
2851: PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);
2852: MatGetSize(matSquare, &dof, NULL);
2853: ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis);
2854: MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &multMat);
2855: MatDestroy(&matSquare);
2856: ISDestroy(&rowis);
2857: }
2858: /* Disgusting trick to reuse work vectors */
2859: MatGetLocalSize(multMat, &m, &n);
2860: patch->patchUpdate->map->n = n;
2861: patch->patchRHSWithArtificial->map->n = m;
2862: patch->patchUpdate->map->N = n;
2863: patch->patchRHSWithArtificial->map->N = m;
2864: MatMult(multMat, patch->patchUpdate, patch->patchRHSWithArtificial);
2865: VecScale(patch->patchRHSWithArtificial, -1.0);
2866: PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHSWithArtificial, patch->localRHS, ADD_VALUES, SCATTER_REVERSE, SCATTER_WITHARTIFICIAL);
2867: if (!patch->save_operators) {
2868: MatDestroy(&multMat);
2869: }
2870: return(0);
2871: }
2873: static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y)
2874: {
2875: PC_PATCH *patch = (PC_PATCH *) pc->data;
2876: const PetscScalar *globalRHS = NULL;
2877: PetscScalar *localRHS = NULL;
2878: PetscScalar *globalUpdate = NULL;
2879: const PetscInt *bcNodes = NULL;
2880: PetscInt nsweep = patch->symmetrise_sweep ? 2 : 1;
2881: PetscInt start[2] = {0, 0};
2882: PetscInt end[2] = {-1, -1};
2883: const PetscInt inc[2] = {1, -1};
2884: const PetscScalar *localUpdate;
2885: const PetscInt *iterationSet;
2886: PetscInt pStart, numBcs, n, sweep, bc, j;
2887: PetscErrorCode ierr;
2890: PetscLogEventBegin(PC_Patch_Apply, pc, 0, 0, 0);
2891: PetscOptionsPushGetViewerOff(PETSC_TRUE);
2892: /* start, end, inc have 2 entries to manage a second backward sweep if we symmetrize */
2893: end[0] = patch->npatch;
2894: start[1] = patch->npatch-1;
2895: if (patch->user_patches) {
2896: ISGetLocalSize(patch->iterationSet, &end[0]);
2897: start[1] = end[0] - 1;
2898: ISGetIndices(patch->iterationSet, &iterationSet);
2899: }
2900: /* Scatter from global space into overlapped local spaces */
2901: VecGetArrayRead(x, &globalRHS);
2902: VecGetArray(patch->localRHS, &localRHS);
2903: PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, globalRHS, localRHS,MPI_REPLACE);
2904: PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, globalRHS, localRHS,MPI_REPLACE);
2905: VecRestoreArrayRead(x, &globalRHS);
2906: VecRestoreArray(patch->localRHS, &localRHS);
2908: VecSet(patch->localUpdate, 0.0);
2909: PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);
2910: PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);
2911: for (sweep = 0; sweep < nsweep; sweep++) {
2912: for (j = start[sweep]; j*inc[sweep] < end[sweep]*inc[sweep]; j += inc[sweep]) {
2913: PetscInt i = patch->user_patches ? iterationSet[j] : j;
2914: PetscInt start, len;
2916: PetscSectionGetDof(patch->gtolCounts, i+pStart, &len);
2917: PetscSectionGetOffset(patch->gtolCounts, i+pStart, &start);
2918: /* TODO: Squash out these guys in the setup as well. */
2919: if (len <= 0) continue;
2920: /* TODO: Do we need different scatters for X and Y? */
2921: PCPatch_ScatterLocal_Private(pc, i+pStart, patch->localRHS, patch->patchRHS, INSERT_VALUES, SCATTER_FORWARD, SCATTER_INTERIOR);
2922: (*patch->applysolver)(pc, i, patch->patchRHS, patch->patchUpdate);
2923: PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchUpdate, patch->localUpdate, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);
2924: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2925: (*patch->updatemultiplicative)(pc, i, pStart);
2926: }
2927: }
2928: }
2929: PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);
2930: if (patch->user_patches) {ISRestoreIndices(patch->iterationSet, &iterationSet);}
2931: /* XXX: should we do this on the global vector? */
2932: if (patch->partition_of_unity) {
2933: VecPointwiseMult(patch->localUpdate, patch->localUpdate, patch->dof_weights);
2934: }
2935: /* Now patch->localUpdate contains the solution of the patch solves, so we need to combine them all. */
2936: VecSet(y, 0.0);
2937: VecGetArray(y, &globalUpdate);
2938: VecGetArrayRead(patch->localUpdate, &localUpdate);
2939: PetscSFReduceBegin(patch->sectionSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);
2940: PetscSFReduceEnd(patch->sectionSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);
2941: VecRestoreArrayRead(patch->localUpdate, &localUpdate);
2943: /* Now we need to send the global BC values through */
2944: VecGetArrayRead(x, &globalRHS);
2945: ISGetSize(patch->globalBcNodes, &numBcs);
2946: ISGetIndices(patch->globalBcNodes, &bcNodes);
2947: VecGetLocalSize(x, &n);
2948: for (bc = 0; bc < numBcs; ++bc) {
2949: const PetscInt idx = bcNodes[bc];
2950: if (idx < n) globalUpdate[idx] = globalRHS[idx];
2951: }
2953: ISRestoreIndices(patch->globalBcNodes, &bcNodes);
2954: VecRestoreArrayRead(x, &globalRHS);
2955: VecRestoreArray(y, &globalUpdate);
2957: PetscOptionsPopGetViewerOff();
2958: PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0);
2959: return(0);
2960: }
2962: static PetscErrorCode PCReset_PATCH_Linear(PC pc)
2963: {
2964: PC_PATCH *patch = (PC_PATCH *) pc->data;
2965: PetscInt i;
2969: if (patch->solver) {
2970: for (i = 0; i < patch->npatch; ++i) {KSPReset((KSP) patch->solver[i]);}
2971: }
2972: return(0);
2973: }
2975: static PetscErrorCode PCReset_PATCH(PC pc)
2976: {
2977: PC_PATCH *patch = (PC_PATCH *) pc->data;
2978: PetscInt i;
2983: PetscSFDestroy(&patch->sectionSF);
2984: PetscSectionDestroy(&patch->cellCounts);
2985: PetscSectionDestroy(&patch->pointCounts);
2986: PetscSectionDestroy(&patch->cellNumbering);
2987: PetscSectionDestroy(&patch->gtolCounts);
2988: ISDestroy(&patch->gtol);
2989: ISDestroy(&patch->cells);
2990: ISDestroy(&patch->points);
2991: ISDestroy(&patch->dofs);
2992: ISDestroy(&patch->offs);
2993: PetscSectionDestroy(&patch->patchSection);
2994: ISDestroy(&patch->ghostBcNodes);
2995: ISDestroy(&patch->globalBcNodes);
2996: PetscSectionDestroy(&patch->gtolCountsWithArtificial);
2997: ISDestroy(&patch->gtolWithArtificial);
2998: ISDestroy(&patch->dofsWithArtificial);
2999: ISDestroy(&patch->offsWithArtificial);
3000: PetscSectionDestroy(&patch->gtolCountsWithAll);
3001: ISDestroy(&patch->gtolWithAll);
3002: ISDestroy(&patch->dofsWithAll);
3003: ISDestroy(&patch->offsWithAll);
3004: VecDestroy(&patch->cellMats);
3005: VecDestroy(&patch->intFacetMats);
3006: ISDestroy(&patch->allCells);
3007: ISDestroy(&patch->intFacets);
3008: ISDestroy(&patch->extFacets);
3009: ISDestroy(&patch->intFacetsToPatchCell);
3010: PetscSectionDestroy(&patch->intFacetCounts);
3011: PetscSectionDestroy(&patch->extFacetCounts);
3013: if (patch->dofSection) for (i = 0; i < patch->nsubspaces; i++) {PetscSectionDestroy(&patch->dofSection[i]);}
3014: PetscFree(patch->dofSection);
3015: PetscFree(patch->bs);
3016: PetscFree(patch->nodesPerCell);
3017: if (patch->cellNodeMap) for (i = 0; i < patch->nsubspaces; i++) {PetscFree(patch->cellNodeMap[i]);}
3018: PetscFree(patch->cellNodeMap);
3019: PetscFree(patch->subspaceOffsets);
3021: (*patch->resetsolver)(pc);
3023: if (patch->subspaces_to_exclude) {
3024: PetscHSetIDestroy(&patch->subspaces_to_exclude);
3025: }
3027: VecDestroy(&patch->localRHS);
3028: VecDestroy(&patch->localUpdate);
3029: VecDestroy(&patch->patchRHS);
3030: VecDestroy(&patch->patchUpdate);
3031: VecDestroy(&patch->dof_weights);
3032: if (patch->patch_dof_weights) {
3033: for (i = 0; i < patch->npatch; ++i) {VecDestroy(&patch->patch_dof_weights[i]);}
3034: PetscFree(patch->patch_dof_weights);
3035: }
3036: if (patch->mat) {
3037: for (i = 0; i < patch->npatch; ++i) {MatDestroy(&patch->mat[i]);}
3038: PetscFree(patch->mat);
3039: }
3040: if (patch->matWithArtificial) {
3041: for (i = 0; i < patch->npatch; ++i) {MatDestroy(&patch->matWithArtificial[i]);}
3042: PetscFree(patch->matWithArtificial);
3043: }
3044: VecDestroy(&patch->patchRHSWithArtificial);
3045: if (patch->dofMappingWithoutToWithArtificial) {
3046: for (i = 0; i < patch->npatch; ++i) {ISDestroy(&patch->dofMappingWithoutToWithArtificial[i]);}
3047: PetscFree(patch->dofMappingWithoutToWithArtificial);
3049: }
3050: if (patch->dofMappingWithoutToWithAll) {
3051: for (i = 0; i < patch->npatch; ++i) {ISDestroy(&patch->dofMappingWithoutToWithAll[i]);}
3052: PetscFree(patch->dofMappingWithoutToWithAll);
3054: }
3055: PetscFree(patch->sub_mat_type);
3056: if (patch->userIS) {
3057: for (i = 0; i < patch->npatch; ++i) {ISDestroy(&patch->userIS[i]);}
3058: PetscFree(patch->userIS);
3059: }
3060: PetscFree(patch->precomputedTensorLocations);
3061: PetscFree(patch->precomputedIntFacetTensorLocations);
3063: patch->bs = NULL;
3064: patch->cellNodeMap = NULL;
3065: patch->nsubspaces = 0;
3066: ISDestroy(&patch->iterationSet);
3068: PetscViewerDestroy(&patch->viewerSection);
3069: return(0);
3070: }
3072: static PetscErrorCode PCDestroy_PATCH_Linear(PC pc)
3073: {
3074: PC_PATCH *patch = (PC_PATCH *) pc->data;
3075: PetscInt i;
3079: if (patch->solver) {
3080: for (i = 0; i < patch->npatch; ++i) {KSPDestroy((KSP *) &patch->solver[i]);}
3081: PetscFree(patch->solver);
3082: }
3083: return(0);
3084: }
3086: static PetscErrorCode PCDestroy_PATCH(PC pc)
3087: {
3088: PC_PATCH *patch = (PC_PATCH *) pc->data;
3092: PCReset_PATCH(pc);
3093: (*patch->destroysolver)(pc);
3094: PetscFree(pc->data);
3095: return(0);
3096: }
3098: static PetscErrorCode PCSetFromOptions_PATCH(PetscOptionItems *PetscOptionsObject, PC pc)
3099: {
3100: PC_PATCH *patch = (PC_PATCH *) pc->data;
3101: PCPatchConstructType patchConstructionType = PC_PATCH_STAR;
3102: char sub_mat_type[PETSC_MAX_PATH_LEN];
3103: char option[PETSC_MAX_PATH_LEN];
3104: const char *prefix;
3105: PetscBool flg, dimflg, codimflg;
3106: MPI_Comm comm;
3107: PetscInt *ifields, nfields, k;
3108: PetscErrorCode ierr;
3109: PCCompositeType loctype = PC_COMPOSITE_ADDITIVE;
3112: PetscObjectGetComm((PetscObject) pc, &comm);
3113: PetscObjectGetOptionsPrefix((PetscObject) pc, &prefix);
3114: PetscOptionsHead(PetscOptionsObject, "Patch solver options");
3116: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_save_operators", patch->classname);
3117: PetscOptionsBool(option, "Store all patch operators for lifetime of object?", "PCPatchSetSaveOperators", patch->save_operators, &patch->save_operators, &flg);
3119: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_precompute_element_tensors", patch->classname);
3120: PetscOptionsBool(option, "Compute each element tensor only once?", "PCPatchSetPrecomputeElementTensors", patch->precomputeElementTensors, &patch->precomputeElementTensors, &flg);
3121: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_partition_of_unity", patch->classname);
3122: PetscOptionsBool(option, "Weight contributions by dof multiplicity?", "PCPatchSetPartitionOfUnity", patch->partition_of_unity, &patch->partition_of_unity, &flg);
3124: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_local_type", patch->classname);
3125: PetscOptionsEnum(option,"Type of local solver composition (additive or multiplicative)","PCPatchSetLocalComposition",PCCompositeTypes,(PetscEnum)loctype,(PetscEnum*)&loctype,&flg);
3126: if (flg) { PCPatchSetLocalComposition(pc, loctype);}
3127: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_dense_inverse", patch->classname);
3128: PetscOptionsBool(option, "Compute inverses of patch matrices and apply directly? Ignores KSP/PC settings on patch.", "PCPatchSetDenseInverse", patch->denseinverse, &patch->denseinverse, &flg);
3129: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_dim", patch->classname);
3130: PetscOptionsInt(option, "What dimension of mesh point to construct patches by? (0 = vertices)", "PCPATCH", patch->dim, &patch->dim, &dimflg);
3131: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_codim", patch->classname);
3132: PetscOptionsInt(option, "What co-dimension of mesh point to construct patches by? (0 = cells)", "PCPATCH", patch->codim, &patch->codim, &codimflg);
3133: if (dimflg && codimflg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Can only set one of dimension or co-dimension");
3135: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_type", patch->classname);
3136: PetscOptionsEnum(option, "How should the patches be constructed?", "PCPatchSetConstructType", PCPatchConstructTypes, (PetscEnum) patchConstructionType, (PetscEnum *) &patchConstructionType, &flg);
3137: if (flg) {PCPatchSetConstructType(pc, patchConstructionType, NULL, NULL);}
3139: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_vanka_dim", patch->classname);
3140: PetscOptionsInt(option, "Topological dimension of entities for Vanka to ignore", "PCPATCH", patch->vankadim, &patch->vankadim, &flg);
3142: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_ignore_dim", patch->classname);
3143: PetscOptionsInt(option, "Topological dimension of entities for completion to ignore", "PCPATCH", patch->ignoredim, &patch->ignoredim, &flg);
3145: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_pardecomp_overlap", patch->classname);
3146: PetscOptionsInt(option, "What overlap should we use in construct type pardecomp?", "PCPATCH", patch->pardecomp_overlap, &patch->pardecomp_overlap, &flg);
3148: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_sub_mat_type", patch->classname);
3149: PetscOptionsFList(option, "Matrix type for patch solves", "PCPatchSetSubMatType", MatList, NULL, sub_mat_type, PETSC_MAX_PATH_LEN, &flg);
3150: if (flg) {PCPatchSetSubMatType(pc, sub_mat_type);}
3152: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_symmetrise_sweep", patch->classname);
3153: PetscOptionsBool(option, "Go start->end, end->start?", "PCPATCH", patch->symmetrise_sweep, &patch->symmetrise_sweep, &flg);
3155: /* If the user has set the number of subspaces, use that for the buffer size,
3156: otherwise use a large number */
3157: if (patch->nsubspaces <= 0) {
3158: nfields = 128;
3159: } else {
3160: nfields = patch->nsubspaces;
3161: }
3162: PetscMalloc1(nfields, &ifields);
3163: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exclude_subspaces", patch->classname);
3164: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,option,ifields,&nfields,&flg);
3165: if (flg && (patchConstructionType == PC_PATCH_USER)) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "We cannot support excluding a subspace with user patches because we do not index patches with a mesh point");
3166: if (flg) {
3167: PetscHSetIClear(patch->subspaces_to_exclude);
3168: for (k = 0; k < nfields; k++) {
3169: PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]);
3170: }
3171: }
3172: PetscFree(ifields);
3174: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_patches_view", patch->classname);
3175: PetscOptionsBool(option, "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg);
3176: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_cells_view", patch->classname);
3177: PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerCells, &patch->formatCells, &patch->viewCells);
3178: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_interior_facets_view", patch->classname);
3179: PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerIntFacets, &patch->formatIntFacets, &patch->viewIntFacets);
3180: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exterior_facets_view", patch->classname);
3181: PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerExtFacets, &patch->formatExtFacets, &patch->viewExtFacets);
3182: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_points_view", patch->classname);
3183: PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints);
3184: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_section_view", patch->classname);
3185: PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerSection, &patch->formatSection, &patch->viewSection);
3186: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_mat_view", patch->classname);
3187: PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix);
3188: PetscOptionsTail();
3189: patch->optionsSet = PETSC_TRUE;
3190: return(0);
3191: }
3193: static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc)
3194: {
3195: PC_PATCH *patch = (PC_PATCH*) pc->data;
3196: KSPConvergedReason reason;
3197: PetscInt i;
3198: PetscErrorCode ierr;
3201: if (!patch->save_operators) {
3202: /* Can't do this here because the sub KSPs don't have an operator attached yet. */
3203: return(0);
3204: }
3205: if (patch->denseinverse) {
3206: /* No solvers */
3207: return(0);
3208: }
3209: for (i = 0; i < patch->npatch; ++i) {
3210: if (!((KSP) patch->solver[i])->setfromoptionscalled) {
3211: KSPSetFromOptions((KSP) patch->solver[i]);
3212: }
3213: KSPSetUp((KSP) patch->solver[i]);
3214: KSPGetConvergedReason((KSP) patch->solver[i], &reason);
3215: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
3216: }
3217: return(0);
3218: }
3220: static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer)
3221: {
3222: PC_PATCH *patch = (PC_PATCH *) pc->data;
3223: PetscViewer sviewer;
3224: PetscBool isascii;
3225: PetscMPIInt rank;
3229: /* TODO Redo tabbing with set tbas in new style */
3230: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
3231: if (!isascii) return(0);
3232: MPI_Comm_rank(PetscObjectComm((PetscObject) pc), &rank);
3233: PetscViewerASCIIPushTab(viewer);
3234: PetscViewerASCIIPrintf(viewer, "Subspace Correction preconditioner with %d patches\n", patch->npatch);
3235: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
3236: PetscViewerASCIIPrintf(viewer, "Schwarz type: multiplicative\n");
3237: } else {
3238: PetscViewerASCIIPrintf(viewer, "Schwarz type: additive\n");
3239: }
3240: if (patch->partition_of_unity) {PetscViewerASCIIPrintf(viewer, "Weighting by partition of unity\n");}
3241: else {PetscViewerASCIIPrintf(viewer, "Not weighting by partition of unity\n");}
3242: if (patch->symmetrise_sweep) {PetscViewerASCIIPrintf(viewer, "Symmetrising sweep (start->end, then end->start)\n");}
3243: else {PetscViewerASCIIPrintf(viewer, "Not symmetrising sweep\n");}
3244: if (!patch->precomputeElementTensors) {PetscViewerASCIIPrintf(viewer, "Not precomputing element tensors (overlapping cells rebuilt in every patch assembly)\n");}
3245: else {PetscViewerASCIIPrintf(viewer, "Precomputing element tensors (each cell assembled only once)\n");}
3246: if (!patch->save_operators) {PetscViewerASCIIPrintf(viewer, "Not saving patch operators (rebuilt every PCApply)\n");}
3247: else {PetscViewerASCIIPrintf(viewer, "Saving patch operators (rebuilt every PCSetUp)\n");}
3248: if (patch->patchconstructop == PCPatchConstruct_Star) {PetscViewerASCIIPrintf(viewer, "Patch construction operator: star\n");}
3249: else if (patch->patchconstructop == PCPatchConstruct_Vanka) {PetscViewerASCIIPrintf(viewer, "Patch construction operator: Vanka\n");}
3250: else if (patch->patchconstructop == PCPatchConstruct_User) {PetscViewerASCIIPrintf(viewer, "Patch construction operator: user-specified\n");}
3251: else {PetscViewerASCIIPrintf(viewer, "Patch construction operator: unknown\n");}
3253: if (patch->denseinverse) {
3254: PetscViewerASCIIPrintf(viewer, "Explicitly forming dense inverse and applying patch solver via MatMult.\n");
3255: } else {
3256: if (patch->isNonlinear) {
3257: PetscViewerASCIIPrintf(viewer, "SNES on patches (all same):\n");
3258: } else {
3259: PetscViewerASCIIPrintf(viewer, "KSP on patches (all same):\n");
3260: }
3261: if (patch->solver) {
3262: PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
3263: if (rank == 0) {
3264: PetscViewerASCIIPushTab(sviewer);
3265: PetscObjectView(patch->solver[0], sviewer);
3266: PetscViewerASCIIPopTab(sviewer);
3267: }
3268: PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
3269: } else {
3270: PetscViewerASCIIPushTab(viewer);
3271: PetscViewerASCIIPrintf(viewer, "Solver not yet set.\n");
3272: PetscViewerASCIIPopTab(viewer);
3273: }
3274: }
3275: PetscViewerASCIIPopTab(viewer);
3276: return(0);
3277: }
3279: /*MC
3280: PCPATCH - A PC object that encapsulates flexible definition of blocks for overlapping and non-overlapping
3281: small block additive preconditioners. Block definition is based on topology from
3282: a DM and equation numbering from a PetscSection.
3284: Options Database Keys:
3285: + -pc_patch_cells_view - Views the process local cell numbers for each patch
3286: . -pc_patch_points_view - Views the process local mesh point numbers for each patch
3287: . -pc_patch_g2l_view - Views the map between global dofs and patch local dofs for each patch
3288: . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary
3289: - -pc_patch_sub_mat_view - Views the matrix associated with each patch
3291: Level: intermediate
3293: .seealso: PCType, PCCreate(), PCSetType()
3294: M*/
3295: PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc)
3296: {
3297: PC_PATCH *patch;
3301: PetscNewLog(pc, &patch);
3303: if (patch->subspaces_to_exclude) {
3304: PetscHSetIDestroy(&patch->subspaces_to_exclude);
3305: }
3306: PetscHSetICreate(&patch->subspaces_to_exclude);
3308: patch->classname = "pc";
3309: patch->isNonlinear = PETSC_FALSE;
3311: /* Set some defaults */
3312: patch->combined = PETSC_FALSE;
3313: patch->save_operators = PETSC_TRUE;
3314: patch->local_composition_type = PC_COMPOSITE_ADDITIVE;
3315: patch->precomputeElementTensors = PETSC_FALSE;
3316: patch->partition_of_unity = PETSC_FALSE;
3317: patch->codim = -1;
3318: patch->dim = -1;
3319: patch->vankadim = -1;
3320: patch->ignoredim = -1;
3321: patch->pardecomp_overlap = 0;
3322: patch->patchconstructop = PCPatchConstruct_Star;
3323: patch->symmetrise_sweep = PETSC_FALSE;
3324: patch->npatch = 0;
3325: patch->userIS = NULL;
3326: patch->optionsSet = PETSC_FALSE;
3327: patch->iterationSet = NULL;
3328: patch->user_patches = PETSC_FALSE;
3329: PetscStrallocpy(MATDENSE, (char **) &patch->sub_mat_type);
3330: patch->viewPatches = PETSC_FALSE;
3331: patch->viewCells = PETSC_FALSE;
3332: patch->viewPoints = PETSC_FALSE;
3333: patch->viewSection = PETSC_FALSE;
3334: patch->viewMatrix = PETSC_FALSE;
3335: patch->densesolve = NULL;
3336: patch->setupsolver = PCSetUp_PATCH_Linear;
3337: patch->applysolver = PCApply_PATCH_Linear;
3338: patch->resetsolver = PCReset_PATCH_Linear;
3339: patch->destroysolver = PCDestroy_PATCH_Linear;
3340: patch->updatemultiplicative = PCUpdateMultiplicative_PATCH_Linear;
3341: patch->dofMappingWithoutToWithArtificial = NULL;
3342: patch->dofMappingWithoutToWithAll = NULL;
3344: pc->data = (void *) patch;
3345: pc->ops->apply = PCApply_PATCH;
3346: pc->ops->applytranspose = NULL; /* PCApplyTranspose_PATCH; */
3347: pc->ops->setup = PCSetUp_PATCH;
3348: pc->ops->reset = PCReset_PATCH;
3349: pc->ops->destroy = PCDestroy_PATCH;
3350: pc->ops->setfromoptions = PCSetFromOptions_PATCH;
3351: pc->ops->setuponblocks = PCSetUpOnBlocks_PATCH;
3352: pc->ops->view = PCView_PATCH;
3353: pc->ops->applyrichardson = NULL;
3355: return(0);
3356: }