Actual source code: pcpatch.c
petsc-3.14.6 2021-03-30
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);
268: PetscSFBcastEnd(rankSF, contig, offsets, remoteOffsets);
269: MPI_Type_free(&contig);
270: PetscFree(offsets);
271: PetscSFDestroy(&rankSF);
272: /* Now remoteOffsets contains the offsets on the remote
273: * processes who communicate with me. So now we can
274: * concatenate the list of SFs into a single one. */
275: index = 0;
276: for (i = 0; i < n; ++i) {
277: const PetscSFNode *remote = NULL;
278: const PetscInt *local = NULL;
279: PetscInt nroots, nleaves, j;
281: PetscSFGetGraph(sf[i], &nroots, &nleaves, &local, &remote);
282: for (j = 0; j < nleaves; ++j) {
283: PetscInt rank = remote[j].rank;
284: PetscInt idx, rootOffset, k;
286: PetscHMapIGet(rankToIndex, rank, &idx);
287: if (idx == -1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Didn't find rank, huh?");
288: /* Offset on given rank for ith subspace */
289: rootOffset = remoteOffsets[n*idx + i];
290: for (k = 0; k < bs[i]; ++k) {
291: ilocal[index] = (local ? local[j] : j)*bs[i] + k + leafOffset;
292: iremote[index].rank = remote[j].rank;
293: iremote[index].index = remote[j].index*bs[i] + k + rootOffset;
294: ++index;
295: }
296: }
297: leafOffset += nleaves * bs[i];
298: }
299: PetscHMapIDestroy(&rankToIndex);
300: PetscFree(remoteOffsets);
301: PetscSFCreate(PetscObjectComm((PetscObject)pc), &patch->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);
1091: /* Now that we know how much space we need, run through again and actually remember the cells. */
1092: for (v = vStart; v < vEnd; v++) {
1093: PetscHashIter hi;
1094: PetscInt dof, off, cdof, coff, efdof, efoff, ifdof, ifoff, pdof, n = 0, cn = 0, ifn = 0, efn = 0;
1096: PetscSectionGetDof(pointCounts, v, &dof);
1097: PetscSectionGetOffset(pointCounts, v, &off);
1098: PetscSectionGetDof(cellCounts, v, &cdof);
1099: PetscSectionGetOffset(cellCounts, v, &coff);
1100: PetscSectionGetDof(intFacetCounts, v, &ifdof);
1101: PetscSectionGetOffset(intFacetCounts, v, &ifoff);
1102: PetscSectionGetDof(extFacetCounts, v, &efdof);
1103: PetscSectionGetOffset(extFacetCounts, v, &efoff);
1104: if (dof <= 0) continue;
1105: patch->patchconstructop((void *) patch, dm, v, ht);
1106: PCPatchCompleteCellPatch(pc, ht, cht);
1107: PetscHashIterBegin(cht, hi);
1108: while (!PetscHashIterAtEnd(cht, hi)) {
1109: PetscInt point;
1111: PetscHashIterGetKey(cht, hi, point);
1112: if (fStart <= point && point < fEnd) {
1113: const PetscInt *support;
1114: PetscInt supportSize, p;
1115: PetscBool interior = PETSC_TRUE;
1116: DMPlexGetSupport(dm, point, &support);
1117: DMPlexGetSupportSize(dm, point, &supportSize);
1118: if (supportSize == 1) {
1119: interior = PETSC_FALSE;
1120: } else {
1121: for (p = 0; p < supportSize; p++) {
1122: PetscBool found;
1123: /* FIXME: can I do this while iterating over cht? */
1124: PetscHSetIHas(cht, support[p], &found);
1125: if (!found) {
1126: interior = PETSC_FALSE;
1127: break;
1128: }
1129: }
1130: }
1131: if (interior) {
1132: intFacetsToPatchCell[2*(ifoff + ifn)] = support[0];
1133: intFacetsToPatchCell[2*(ifoff + ifn) + 1] = support[1];
1134: intFacetsArray[ifoff + ifn++] = point;
1135: } else {
1136: extFacetsArray[efoff + efn++] = point;
1137: }
1138: }
1139: PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);
1140: if (pdof) {pointsArray[off + n++] = point;}
1141: if (point >= cStart && point < cEnd) {cellsArray[coff + cn++] = point;}
1142: PetscHashIterNext(cht, hi);
1143: }
1144: 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);
1145: 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);
1146: 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);
1147: 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);
1149: for (ifn = 0; ifn < ifdof; ifn++) {
1150: PetscInt cell0 = intFacetsToPatchCell[2*(ifoff + ifn)];
1151: PetscInt cell1 = intFacetsToPatchCell[2*(ifoff + ifn) + 1];
1152: PetscBool found0 = PETSC_FALSE, found1 = PETSC_FALSE;
1153: for (n = 0; n < cdof; n++) {
1154: if (!found0 && cell0 == cellsArray[coff + n]) {
1155: intFacetsToPatchCell[2*(ifoff + ifn)] = n;
1156: found0 = PETSC_TRUE;
1157: }
1158: if (!found1 && cell1 == cellsArray[coff + n]) {
1159: intFacetsToPatchCell[2*(ifoff + ifn) + 1] = n;
1160: found1 = PETSC_TRUE;
1161: }
1162: if (found0 && found1) break;
1163: }
1164: if (!(found0 && found1)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Didn't manage to find local point numbers for facet support");
1165: }
1166: }
1167: PetscHSetIDestroy(&ht);
1168: PetscHSetIDestroy(&cht);
1170: ISCreateGeneral(PETSC_COMM_SELF, numCells, cellsArray, PETSC_OWN_POINTER, &patch->cells);
1171: PetscObjectSetName((PetscObject) patch->cells, "Patch Cells");
1172: if (patch->viewCells) {
1173: ObjectView((PetscObject) patch->cellCounts, patch->viewerCells, patch->formatCells);
1174: ObjectView((PetscObject) patch->cells, patch->viewerCells, patch->formatCells);
1175: }
1176: ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray, PETSC_OWN_POINTER, &patch->intFacets);
1177: PetscObjectSetName((PetscObject) patch->intFacets, "Patch Interior Facets");
1178: ISCreateGeneral(PETSC_COMM_SELF, 2*numIntFacets, intFacetsToPatchCell, PETSC_OWN_POINTER, &patch->intFacetsToPatchCell);
1179: PetscObjectSetName((PetscObject) patch->intFacetsToPatchCell, "Patch Interior Facets local support");
1180: if (patch->viewIntFacets) {
1181: ObjectView((PetscObject) patch->intFacetCounts, patch->viewerIntFacets, patch->formatIntFacets);
1182: ObjectView((PetscObject) patch->intFacets, patch->viewerIntFacets, patch->formatIntFacets);
1183: ObjectView((PetscObject) patch->intFacetsToPatchCell, patch->viewerIntFacets, patch->formatIntFacets);
1184: }
1185: ISCreateGeneral(PETSC_COMM_SELF, numExtFacets, extFacetsArray, PETSC_OWN_POINTER, &patch->extFacets);
1186: PetscObjectSetName((PetscObject) patch->extFacets, "Patch Exterior Facets");
1187: if (patch->viewExtFacets) {
1188: ObjectView((PetscObject) patch->extFacetCounts, patch->viewerExtFacets, patch->formatExtFacets);
1189: ObjectView((PetscObject) patch->extFacets, patch->viewerExtFacets, patch->formatExtFacets);
1190: }
1191: ISCreateGeneral(PETSC_COMM_SELF, numPoints, pointsArray, PETSC_OWN_POINTER, &patch->points);
1192: PetscObjectSetName((PetscObject) patch->points, "Patch Points");
1193: if (patch->viewPoints) {
1194: ObjectView((PetscObject) patch->pointCounts, patch->viewerPoints, patch->formatPoints);
1195: ObjectView((PetscObject) patch->points, patch->viewerPoints, patch->formatPoints);
1196: }
1197: DMDestroy(&dm);
1198: return(0);
1199: }
1201: /*
1202: * PCPatchCreateCellPatchDiscretisationInfo - Build the dof maps for cell patches
1203: *
1204: * Input Parameters:
1205: * + dm - The DMPlex object defining the mesh
1206: * . cellCounts - Section with counts of cells around each vertex
1207: * . cells - IS of the cell point indices of cells in each patch
1208: * . cellNumbering - Section mapping plex cell points to Firedrake cell indices.
1209: * . nodesPerCell - number of nodes per cell.
1210: * - cellNodeMap - map from cells to node indices (nodesPerCell * numCells)
1211: *
1212: * Output Parameters:
1213: * + dofs - IS of local dof numbers of each cell in the patch, where local is a patch local numbering
1214: * . gtolCounts - Section with counts of dofs per cell patch
1215: * - gtol - IS mapping from global dofs to local dofs for each patch.
1216: */
1217: static PetscErrorCode PCPatchCreateCellPatchDiscretisationInfo(PC pc)
1218: {
1219: PC_PATCH *patch = (PC_PATCH *) pc->data;
1220: PetscSection cellCounts = patch->cellCounts;
1221: PetscSection pointCounts = patch->pointCounts;
1222: PetscSection gtolCounts, gtolCountsWithArtificial = NULL, gtolCountsWithAll = NULL;
1223: IS cells = patch->cells;
1224: IS points = patch->points;
1225: PetscSection cellNumbering = patch->cellNumbering;
1226: PetscInt Nf = patch->nsubspaces;
1227: PetscInt numCells, numPoints;
1228: PetscInt numDofs;
1229: PetscInt numGlobalDofs, numGlobalDofsWithArtificial, numGlobalDofsWithAll;
1230: PetscInt totalDofsPerCell = patch->totalDofsPerCell;
1231: PetscInt vStart, vEnd, v;
1232: const PetscInt *cellsArray, *pointsArray;
1233: PetscInt *newCellsArray = NULL;
1234: PetscInt *dofsArray = NULL;
1235: PetscInt *dofsArrayWithArtificial = NULL;
1236: PetscInt *dofsArrayWithAll = NULL;
1237: PetscInt *offsArray = NULL;
1238: PetscInt *offsArrayWithArtificial = NULL;
1239: PetscInt *offsArrayWithAll = NULL;
1240: PetscInt *asmArray = NULL;
1241: PetscInt *asmArrayWithArtificial = NULL;
1242: PetscInt *asmArrayWithAll = NULL;
1243: PetscInt *globalDofsArray = NULL;
1244: PetscInt *globalDofsArrayWithArtificial = NULL;
1245: PetscInt *globalDofsArrayWithAll = NULL;
1246: PetscInt globalIndex = 0;
1247: PetscInt key = 0;
1248: PetscInt asmKey = 0;
1249: DM dm = NULL, plex;
1250: const PetscInt *bcNodes = NULL;
1251: PetscHMapI ht;
1252: PetscHMapI htWithArtificial;
1253: PetscHMapI htWithAll;
1254: PetscHSetI globalBcs;
1255: PetscInt numBcs;
1256: PetscHSetI ownedpts, seenpts, owneddofs, seendofs, artificialbcs;
1257: PetscInt pStart, pEnd, p, i;
1258: char option[PETSC_MAX_PATH_LEN];
1259: PetscBool isNonlinear;
1260: PetscErrorCode ierr;
1264: PCGetDM(pc, &dm);
1265: DMConvert(dm, DMPLEX, &plex);
1266: dm = plex;
1267: /* dofcounts section is cellcounts section * dofPerCell */
1268: PetscSectionGetStorageSize(cellCounts, &numCells);
1269: PetscSectionGetStorageSize(patch->pointCounts, &numPoints);
1270: numDofs = numCells * totalDofsPerCell;
1271: PetscMalloc1(numDofs, &dofsArray);
1272: PetscMalloc1(numPoints*Nf, &offsArray);
1273: PetscMalloc1(numDofs, &asmArray);
1274: PetscMalloc1(numCells, &newCellsArray);
1275: PetscSectionGetChart(cellCounts, &vStart, &vEnd);
1276: PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCounts);
1277: gtolCounts = patch->gtolCounts;
1278: PetscSectionSetChart(gtolCounts, vStart, vEnd);
1279: PetscObjectSetName((PetscObject) patch->gtolCounts, "Patch Global Index Section");
1281: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1282: PetscMalloc1(numPoints*Nf, &offsArrayWithArtificial);
1283: PetscMalloc1(numDofs, &asmArrayWithArtificial);
1284: PetscMalloc1(numDofs, &dofsArrayWithArtificial);
1285: PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithArtificial);
1286: gtolCountsWithArtificial = patch->gtolCountsWithArtificial;
1287: PetscSectionSetChart(gtolCountsWithArtificial, vStart, vEnd);
1288: PetscObjectSetName((PetscObject) patch->gtolCountsWithArtificial, "Patch Global Index Section Including Artificial BCs");
1289: }
1291: isNonlinear = patch->isNonlinear;
1292: if (isNonlinear) {
1293: PetscMalloc1(numPoints*Nf, &offsArrayWithAll);
1294: PetscMalloc1(numDofs, &asmArrayWithAll);
1295: PetscMalloc1(numDofs, &dofsArrayWithAll);
1296: PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithAll);
1297: gtolCountsWithAll = patch->gtolCountsWithAll;
1298: PetscSectionSetChart(gtolCountsWithAll, vStart, vEnd);
1299: PetscObjectSetName((PetscObject) patch->gtolCountsWithAll, "Patch Global Index Section Including All BCs");
1300: }
1302: /* Outside the patch loop, get the dofs that are globally-enforced Dirichlet
1303: conditions */
1304: PetscHSetICreate(&globalBcs);
1305: ISGetIndices(patch->ghostBcNodes, &bcNodes);
1306: ISGetSize(patch->ghostBcNodes, &numBcs);
1307: for (i = 0; i < numBcs; ++i) {
1308: PetscHSetIAdd(globalBcs, bcNodes[i]); /* these are already in concatenated numbering */
1309: }
1310: ISRestoreIndices(patch->ghostBcNodes, &bcNodes);
1311: ISDestroy(&patch->ghostBcNodes); /* memory optimisation */
1313: /* Hash tables for artificial BC construction */
1314: PetscHSetICreate(&ownedpts);
1315: PetscHSetICreate(&seenpts);
1316: PetscHSetICreate(&owneddofs);
1317: PetscHSetICreate(&seendofs);
1318: PetscHSetICreate(&artificialbcs);
1320: ISGetIndices(cells, &cellsArray);
1321: ISGetIndices(points, &pointsArray);
1322: PetscHMapICreate(&ht);
1323: PetscHMapICreate(&htWithArtificial);
1324: PetscHMapICreate(&htWithAll);
1325: for (v = vStart; v < vEnd; ++v) {
1326: PetscInt localIndex = 0;
1327: PetscInt localIndexWithArtificial = 0;
1328: PetscInt localIndexWithAll = 0;
1329: PetscInt dof, off, i, j, k, l;
1331: PetscHMapIClear(ht);
1332: PetscHMapIClear(htWithArtificial);
1333: PetscHMapIClear(htWithAll);
1334: PetscSectionGetDof(cellCounts, v, &dof);
1335: PetscSectionGetOffset(cellCounts, v, &off);
1336: if (dof <= 0) continue;
1338: /* Calculate the global numbers of the artificial BC dofs here first */
1339: patch->patchconstructop((void*)patch, dm, v, ownedpts);
1340: PCPatchCompleteCellPatch(pc, ownedpts, seenpts);
1341: PCPatchGetPointDofs(pc, ownedpts, owneddofs, v, &patch->subspaces_to_exclude);
1342: PCPatchGetPointDofs(pc, seenpts, seendofs, v, NULL);
1343: PCPatchComputeSetDifference_Private(owneddofs, seendofs, artificialbcs);
1344: if (patch->viewPatches) {
1345: PetscHSetI globalbcdofs;
1346: PetscHashIter hi;
1347: MPI_Comm comm = PetscObjectComm((PetscObject)pc);
1349: PetscHSetICreate(&globalbcdofs);
1350: PetscSynchronizedPrintf(comm, "Patch %d: owned dofs:\n", v);
1351: PetscHashIterBegin(owneddofs, hi);
1352: while (!PetscHashIterAtEnd(owneddofs, hi)) {
1353: PetscInt globalDof;
1355: PetscHashIterGetKey(owneddofs, hi, globalDof);
1356: PetscHashIterNext(owneddofs, hi);
1357: PetscSynchronizedPrintf(comm, "%d ", globalDof);
1358: }
1359: PetscSynchronizedPrintf(comm, "\n");
1360: PetscSynchronizedPrintf(comm, "Patch %d: seen dofs:\n", v);
1361: PetscHashIterBegin(seendofs, hi);
1362: while (!PetscHashIterAtEnd(seendofs, hi)) {
1363: PetscInt globalDof;
1364: PetscBool flg;
1366: PetscHashIterGetKey(seendofs, hi, globalDof);
1367: PetscHashIterNext(seendofs, hi);
1368: PetscSynchronizedPrintf(comm, "%d ", globalDof);
1370: PetscHSetIHas(globalBcs, globalDof, &flg);
1371: if (flg) {PetscHSetIAdd(globalbcdofs, globalDof);}
1372: }
1373: PetscSynchronizedPrintf(comm, "\n");
1374: PetscSynchronizedPrintf(comm, "Patch %d: global BCs:\n", v);
1375: PetscHSetIGetSize(globalbcdofs, &numBcs);
1376: if (numBcs > 0) {
1377: PetscHashIterBegin(globalbcdofs, hi);
1378: while (!PetscHashIterAtEnd(globalbcdofs, hi)) {
1379: PetscInt globalDof;
1380: PetscHashIterGetKey(globalbcdofs, hi, globalDof);
1381: PetscHashIterNext(globalbcdofs, hi);
1382: PetscSynchronizedPrintf(comm, "%d ", globalDof);
1383: }
1384: }
1385: PetscSynchronizedPrintf(comm, "\n");
1386: PetscSynchronizedPrintf(comm, "Patch %d: artificial BCs:\n", v);
1387: PetscHSetIGetSize(artificialbcs, &numBcs);
1388: if (numBcs > 0) {
1389: PetscHashIterBegin(artificialbcs, hi);
1390: while (!PetscHashIterAtEnd(artificialbcs, hi)) {
1391: PetscInt globalDof;
1392: PetscHashIterGetKey(artificialbcs, hi, globalDof);
1393: PetscHashIterNext(artificialbcs, hi);
1394: PetscSynchronizedPrintf(comm, "%d ", globalDof);
1395: }
1396: }
1397: PetscSynchronizedPrintf(comm, "\n\n");
1398: PetscHSetIDestroy(&globalbcdofs);
1399: }
1400: for (k = 0; k < patch->nsubspaces; ++k) {
1401: const PetscInt *cellNodeMap = patch->cellNodeMap[k];
1402: PetscInt nodesPerCell = patch->nodesPerCell[k];
1403: PetscInt subspaceOffset = patch->subspaceOffsets[k];
1404: PetscInt bs = patch->bs[k];
1406: for (i = off; i < off + dof; ++i) {
1407: /* Walk over the cells in this patch. */
1408: const PetscInt c = cellsArray[i];
1409: PetscInt cell = c;
1411: /* TODO Change this to an IS */
1412: if (cellNumbering) {
1413: PetscSectionGetDof(cellNumbering, c, &cell);
1414: if (cell <= 0) SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_OUTOFRANGE, "Cell %D doesn't appear in cell numbering map", c);
1415: PetscSectionGetOffset(cellNumbering, c, &cell);
1416: }
1417: newCellsArray[i] = cell;
1418: for (j = 0; j < nodesPerCell; ++j) {
1419: /* For each global dof, map it into contiguous local storage. */
1420: const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + subspaceOffset;
1421: /* finally, loop over block size */
1422: for (l = 0; l < bs; ++l) {
1423: PetscInt localDof;
1424: PetscBool isGlobalBcDof, isArtificialBcDof;
1426: /* first, check if this is either a globally enforced or locally enforced BC dof */
1427: PetscHSetIHas(globalBcs, globalDof + l, &isGlobalBcDof);
1428: PetscHSetIHas(artificialbcs, globalDof + l, &isArtificialBcDof);
1430: /* if it's either, don't ever give it a local dof number */
1431: if (isGlobalBcDof || isArtificialBcDof) {
1432: dofsArray[globalIndex] = -1; /* don't use this in assembly in this patch */
1433: } else {
1434: PetscHMapIGet(ht, globalDof + l, &localDof);
1435: if (localDof == -1) {
1436: localDof = localIndex++;
1437: PetscHMapISet(ht, globalDof + l, localDof);
1438: }
1439: if (globalIndex >= numDofs) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs);
1440: /* And store. */
1441: dofsArray[globalIndex] = localDof;
1442: }
1444: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1445: if (isGlobalBcDof) {
1446: dofsArrayWithArtificial[globalIndex] = -1; /* don't use this in assembly in this patch */
1447: } else {
1448: PetscHMapIGet(htWithArtificial, globalDof + l, &localDof);
1449: if (localDof == -1) {
1450: localDof = localIndexWithArtificial++;
1451: PetscHMapISet(htWithArtificial, globalDof + l, localDof);
1452: }
1453: if (globalIndex >= numDofs) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs);
1454: /* And store.*/
1455: dofsArrayWithArtificial[globalIndex] = localDof;
1456: }
1457: }
1459: if (isNonlinear) {
1460: /* Build the dofmap for the function space with _all_ dofs,
1461: including those in any kind of boundary condition */
1462: PetscHMapIGet(htWithAll, globalDof + l, &localDof);
1463: if (localDof == -1) {
1464: localDof = localIndexWithAll++;
1465: PetscHMapISet(htWithAll, globalDof + l, localDof);
1466: }
1467: if (globalIndex >= numDofs) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs);
1468: /* And store.*/
1469: dofsArrayWithAll[globalIndex] = localDof;
1470: }
1471: globalIndex++;
1472: }
1473: }
1474: }
1475: }
1476: /*How many local dofs in this patch? */
1477: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1478: PetscHMapIGetSize(htWithArtificial, &dof);
1479: PetscSectionSetDof(gtolCountsWithArtificial, v, dof);
1480: }
1481: if (isNonlinear) {
1482: PetscHMapIGetSize(htWithAll, &dof);
1483: PetscSectionSetDof(gtolCountsWithAll, v, dof);
1484: }
1485: PetscHMapIGetSize(ht, &dof);
1486: PetscSectionSetDof(gtolCounts, v, dof);
1487: }
1489: DMDestroy(&dm);
1490: if (globalIndex != numDofs) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Expected number of dofs (%d) doesn't match found number (%d)", numDofs, globalIndex);
1491: PetscSectionSetUp(gtolCounts);
1492: PetscSectionGetStorageSize(gtolCounts, &numGlobalDofs);
1493: PetscMalloc1(numGlobalDofs, &globalDofsArray);
1495: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1496: PetscSectionSetUp(gtolCountsWithArtificial);
1497: PetscSectionGetStorageSize(gtolCountsWithArtificial, &numGlobalDofsWithArtificial);
1498: PetscMalloc1(numGlobalDofsWithArtificial, &globalDofsArrayWithArtificial);
1499: }
1500: if (isNonlinear) {
1501: PetscSectionSetUp(gtolCountsWithAll);
1502: PetscSectionGetStorageSize(gtolCountsWithAll, &numGlobalDofsWithAll);
1503: PetscMalloc1(numGlobalDofsWithAll, &globalDofsArrayWithAll);
1504: }
1505: /* Now populate the global to local map. This could be merged into the above loop if we were willing to deal with reallocs. */
1506: for (v = vStart; v < vEnd; ++v) {
1507: PetscHashIter hi;
1508: PetscInt dof, off, Np, ooff, i, j, k, l;
1510: PetscHMapIClear(ht);
1511: PetscHMapIClear(htWithArtificial);
1512: PetscHMapIClear(htWithAll);
1513: PetscSectionGetDof(cellCounts, v, &dof);
1514: PetscSectionGetOffset(cellCounts, v, &off);
1515: PetscSectionGetDof(pointCounts, v, &Np);
1516: PetscSectionGetOffset(pointCounts, v, &ooff);
1517: if (dof <= 0) continue;
1519: for (k = 0; k < patch->nsubspaces; ++k) {
1520: const PetscInt *cellNodeMap = patch->cellNodeMap[k];
1521: PetscInt nodesPerCell = patch->nodesPerCell[k];
1522: PetscInt subspaceOffset = patch->subspaceOffsets[k];
1523: PetscInt bs = patch->bs[k];
1524: PetscInt goff;
1526: for (i = off; i < off + dof; ++i) {
1527: /* Reconstruct mapping of global-to-local on this patch. */
1528: const PetscInt c = cellsArray[i];
1529: PetscInt cell = c;
1531: if (cellNumbering) {PetscSectionGetOffset(cellNumbering, c, &cell);}
1532: for (j = 0; j < nodesPerCell; ++j) {
1533: for (l = 0; l < bs; ++l) {
1534: const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset;
1535: const PetscInt localDof = dofsArray[key];
1536: if (localDof >= 0) {PetscHMapISet(ht, globalDof, localDof);}
1537: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1538: const PetscInt localDofWithArtificial = dofsArrayWithArtificial[key];
1539: if (localDofWithArtificial >= 0) {
1540: PetscHMapISet(htWithArtificial, globalDof, localDofWithArtificial);
1541: }
1542: }
1543: if (isNonlinear) {
1544: const PetscInt localDofWithAll = dofsArrayWithAll[key];
1545: if (localDofWithAll >= 0) {
1546: PetscHMapISet(htWithAll, globalDof, localDofWithAll);
1547: }
1548: }
1549: key++;
1550: }
1551: }
1552: }
1554: /* Shove it in the output data structure. */
1555: PetscSectionGetOffset(gtolCounts, v, &goff);
1556: PetscHashIterBegin(ht, hi);
1557: while (!PetscHashIterAtEnd(ht, hi)) {
1558: PetscInt globalDof, localDof;
1560: PetscHashIterGetKey(ht, hi, globalDof);
1561: PetscHashIterGetVal(ht, hi, localDof);
1562: if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof;
1563: PetscHashIterNext(ht, hi);
1564: }
1566: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1567: PetscSectionGetOffset(gtolCountsWithArtificial, v, &goff);
1568: PetscHashIterBegin(htWithArtificial, hi);
1569: while (!PetscHashIterAtEnd(htWithArtificial, hi)) {
1570: PetscInt globalDof, localDof;
1571: PetscHashIterGetKey(htWithArtificial, hi, globalDof);
1572: PetscHashIterGetVal(htWithArtificial, hi, localDof);
1573: if (globalDof >= 0) globalDofsArrayWithArtificial[goff + localDof] = globalDof;
1574: PetscHashIterNext(htWithArtificial, hi);
1575: }
1576: }
1577: if (isNonlinear) {
1578: PetscSectionGetOffset(gtolCountsWithAll, v, &goff);
1579: PetscHashIterBegin(htWithAll, hi);
1580: while (!PetscHashIterAtEnd(htWithAll, hi)) {
1581: PetscInt globalDof, localDof;
1582: PetscHashIterGetKey(htWithAll, hi, globalDof);
1583: PetscHashIterGetVal(htWithAll, hi, localDof);
1584: if (globalDof >= 0) globalDofsArrayWithAll[goff + localDof] = globalDof;
1585: PetscHashIterNext(htWithAll, hi);
1586: }
1587: }
1589: for (p = 0; p < Np; ++p) {
1590: const PetscInt point = pointsArray[ooff + p];
1591: PetscInt globalDof, localDof;
1593: PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof);
1594: PetscHMapIGet(ht, globalDof, &localDof);
1595: offsArray[(ooff + p)*Nf + k] = localDof;
1596: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1597: PetscHMapIGet(htWithArtificial, globalDof, &localDof);
1598: offsArrayWithArtificial[(ooff + p)*Nf + k] = localDof;
1599: }
1600: if (isNonlinear) {
1601: PetscHMapIGet(htWithAll, globalDof, &localDof);
1602: offsArrayWithAll[(ooff + p)*Nf + k] = localDof;
1603: }
1604: }
1605: }
1607: PetscHSetIDestroy(&globalBcs);
1608: PetscHSetIDestroy(&ownedpts);
1609: PetscHSetIDestroy(&seenpts);
1610: PetscHSetIDestroy(&owneddofs);
1611: PetscHSetIDestroy(&seendofs);
1612: PetscHSetIDestroy(&artificialbcs);
1614: /* At this point, we have a hash table ht built that maps globalDof -> localDof.
1615: We need to create the dof table laid out cellwise first, then by subspace,
1616: as the assembler assembles cell-wise and we need to stuff the different
1617: contributions of the different function spaces to the right places. So we loop
1618: over cells, then over subspaces. */
1619: if (patch->nsubspaces > 1) { /* for nsubspaces = 1, data we need is already in dofsArray */
1620: for (i = off; i < off + dof; ++i) {
1621: const PetscInt c = cellsArray[i];
1622: PetscInt cell = c;
1624: if (cellNumbering) {PetscSectionGetOffset(cellNumbering, c, &cell);}
1625: for (k = 0; k < patch->nsubspaces; ++k) {
1626: const PetscInt *cellNodeMap = patch->cellNodeMap[k];
1627: PetscInt nodesPerCell = patch->nodesPerCell[k];
1628: PetscInt subspaceOffset = patch->subspaceOffsets[k];
1629: PetscInt bs = patch->bs[k];
1631: for (j = 0; j < nodesPerCell; ++j) {
1632: for (l = 0; l < bs; ++l) {
1633: const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset;
1634: PetscInt localDof;
1636: PetscHMapIGet(ht, globalDof, &localDof);
1637: /* If it's not in the hash table, i.e. is a BC dof,
1638: then the PetscHSetIMap above gives -1, which matches
1639: exactly the convention for PETSc's matrix assembly to
1640: ignore the dof. So we don't need to do anything here */
1641: asmArray[asmKey] = localDof;
1642: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1643: PetscHMapIGet(htWithArtificial, globalDof, &localDof);
1644: asmArrayWithArtificial[asmKey] = localDof;
1645: }
1646: if (isNonlinear) {
1647: PetscHMapIGet(htWithAll, globalDof, &localDof);
1648: asmArrayWithAll[asmKey] = localDof;
1649: }
1650: asmKey++;
1651: }
1652: }
1653: }
1654: }
1655: }
1656: }
1657: if (1 == patch->nsubspaces) {
1658: PetscArraycpy(asmArray, dofsArray, numDofs);
1659: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1660: PetscArraycpy(asmArrayWithArtificial, dofsArrayWithArtificial, numDofs);
1661: }
1662: if (isNonlinear) {
1663: PetscArraycpy(asmArrayWithAll, dofsArrayWithAll, numDofs);
1664: }
1665: }
1667: PetscHMapIDestroy(&ht);
1668: PetscHMapIDestroy(&htWithArtificial);
1669: PetscHMapIDestroy(&htWithAll);
1670: ISRestoreIndices(cells, &cellsArray);
1671: ISRestoreIndices(points, &pointsArray);
1672: PetscFree(dofsArray);
1673: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1674: PetscFree(dofsArrayWithArtificial);
1675: }
1676: if (isNonlinear) {
1677: PetscFree(dofsArrayWithAll);
1678: }
1679: /* Create placeholder section for map from points to patch dofs */
1680: PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection);
1681: PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces);
1682: if (patch->combined) {
1683: PetscInt numFields;
1684: PetscSectionGetNumFields(patch->dofSection[0], &numFields);
1685: 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);
1686: PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd);
1687: PetscSectionSetChart(patch->patchSection, pStart, pEnd);
1688: for (p = pStart; p < pEnd; ++p) {
1689: PetscInt dof, fdof, f;
1691: PetscSectionGetDof(patch->dofSection[0], p, &dof);
1692: PetscSectionSetDof(patch->patchSection, p, dof);
1693: for (f = 0; f < patch->nsubspaces; ++f) {
1694: PetscSectionGetFieldDof(patch->dofSection[0], p, f, &fdof);
1695: PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);
1696: }
1697: }
1698: } else {
1699: PetscInt pStartf, pEndf, f;
1700: pStart = PETSC_MAX_INT;
1701: pEnd = PETSC_MIN_INT;
1702: for (f = 0; f < patch->nsubspaces; ++f) {
1703: PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);
1704: pStart = PetscMin(pStart, pStartf);
1705: pEnd = PetscMax(pEnd, pEndf);
1706: }
1707: PetscSectionSetChart(patch->patchSection, pStart, pEnd);
1708: for (f = 0; f < patch->nsubspaces; ++f) {
1709: PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);
1710: for (p = pStartf; p < pEndf; ++p) {
1711: PetscInt fdof;
1712: PetscSectionGetDof(patch->dofSection[f], p, &fdof);
1713: PetscSectionAddDof(patch->patchSection, p, fdof);
1714: PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);
1715: }
1716: }
1717: }
1718: PetscSectionSetUp(patch->patchSection);
1719: PetscSectionSetUseFieldOffsets(patch->patchSection, PETSC_TRUE);
1720: /* Replace cell indices with firedrake-numbered ones. */
1721: ISGeneralSetIndices(cells, numCells, (const PetscInt *) newCellsArray, PETSC_OWN_POINTER);
1722: ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofs, globalDofsArray, PETSC_OWN_POINTER, &patch->gtol);
1723: PetscObjectSetName((PetscObject) patch->gtol, "Global Indices");
1724: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_g2l_view", patch->classname);
1725: PetscSectionViewFromOptions(patch->gtolCounts, (PetscObject) pc, option);
1726: ISViewFromOptions(patch->gtol, (PetscObject) pc, option);
1727: ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArray, PETSC_OWN_POINTER, &patch->dofs);
1728: ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArray, PETSC_OWN_POINTER, &patch->offs);
1729: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
1730: ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithArtificial, globalDofsArrayWithArtificial, PETSC_OWN_POINTER, &patch->gtolWithArtificial);
1731: ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithArtificial, PETSC_OWN_POINTER, &patch->dofsWithArtificial);
1732: ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArrayWithArtificial, PETSC_OWN_POINTER, &patch->offsWithArtificial);
1733: }
1734: if (isNonlinear) {
1735: ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithAll, globalDofsArrayWithAll, PETSC_OWN_POINTER, &patch->gtolWithAll);
1736: ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithAll, PETSC_OWN_POINTER, &patch->dofsWithAll);
1737: ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArrayWithAll, PETSC_OWN_POINTER, &patch->offsWithAll);
1738: }
1739: return(0);
1740: }
1742: static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat, PetscBool withArtificial)
1743: {
1744: PC_PATCH *patch = (PC_PATCH *) pc->data;
1745: PetscBool flg;
1746: PetscInt csize, rsize;
1747: const char *prefix = NULL;
1751: if (withArtificial) {
1752: /* would be nice if we could create a rectangular matrix of size numDofsWithArtificial x numDofs here */
1753: PetscInt pStart;
1754: PetscSectionGetChart(patch->gtolCountsWithArtificial, &pStart, NULL);
1755: PetscSectionGetDof(patch->gtolCountsWithArtificial, point + pStart, &rsize);
1756: csize = rsize;
1757: } else {
1758: PetscInt pStart;
1759: PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);
1760: PetscSectionGetDof(patch->gtolCounts, point + pStart, &rsize);
1761: csize = rsize;
1762: }
1764: MatCreate(PETSC_COMM_SELF, mat);
1765: PCGetOptionsPrefix(pc, &prefix);
1766: MatSetOptionsPrefix(*mat, prefix);
1767: MatAppendOptionsPrefix(*mat, "pc_patch_sub_");
1768: if (patch->sub_mat_type) {MatSetType(*mat, patch->sub_mat_type);}
1769: else if (!patch->sub_mat_type) {MatSetType(*mat, MATDENSE);}
1770: MatSetSizes(*mat, rsize, csize, rsize, csize);
1771: PetscObjectTypeCompare((PetscObject) *mat, MATDENSE, &flg);
1772: if (!flg) {PetscObjectTypeCompare((PetscObject)*mat, MATSEQDENSE, &flg);}
1773: /* Sparse patch matrices */
1774: if (!flg) {
1775: PetscBT bt;
1776: PetscInt *dnnz = NULL;
1777: const PetscInt *dofsArray = NULL;
1778: PetscInt pStart, pEnd, ncell, offset, c, i, j;
1780: if (withArtificial) {
1781: ISGetIndices(patch->dofsWithArtificial, &dofsArray);
1782: } else {
1783: ISGetIndices(patch->dofs, &dofsArray);
1784: }
1785: PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);
1786: point += pStart;
1787: if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);
1788: PetscSectionGetDof(patch->cellCounts, point, &ncell);
1789: PetscSectionGetOffset(patch->cellCounts, point, &offset);
1790: PetscLogEventBegin(PC_Patch_Prealloc, pc, 0, 0, 0);
1791: /* A PetscBT uses N^2 bits to store the sparsity pattern on a
1792: * patch. This is probably OK if the patches are not too big,
1793: * but uses too much memory. We therefore switch based on rsize. */
1794: if (rsize < 3000) { /* FIXME: I picked this switch value out of my hat */
1795: PetscScalar *zeroes;
1796: PetscInt rows;
1798: PetscCalloc1(rsize, &dnnz);
1799: PetscBTCreate(rsize*rsize, &bt);
1800: for (c = 0; c < ncell; ++c) {
1801: const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell;
1802: for (i = 0; i < patch->totalDofsPerCell; ++i) {
1803: const PetscInt row = idx[i];
1804: if (row < 0) continue;
1805: for (j = 0; j < patch->totalDofsPerCell; ++j) {
1806: const PetscInt col = idx[j];
1807: const PetscInt key = row*rsize + col;
1808: if (col < 0) continue;
1809: if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1810: }
1811: }
1812: }
1814: if (patch->usercomputeopintfacet) {
1815: const PetscInt *intFacetsArray = NULL;
1816: PetscInt i, numIntFacets, intFacetOffset;
1817: const PetscInt *facetCells = NULL;
1819: PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1820: PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1821: ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1822: ISGetIndices(patch->intFacets, &intFacetsArray);
1823: for (i = 0; i < numIntFacets; i++) {
1824: const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0];
1825: const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1];
1826: PetscInt celli, cellj;
1828: for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1829: const PetscInt row = dofsArray[(offset + cell0)*patch->totalDofsPerCell + celli];
1830: if (row < 0) continue;
1831: for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1832: const PetscInt col = dofsArray[(offset + cell1)*patch->totalDofsPerCell + cellj];
1833: const PetscInt key = row*rsize + col;
1834: if (col < 0) continue;
1835: if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1836: }
1837: }
1839: for (celli = 0; celli < patch->totalDofsPerCell; celli++) {
1840: const PetscInt row = dofsArray[(offset + cell1)*patch->totalDofsPerCell + celli];
1841: if (row < 0) continue;
1842: for (cellj = 0; cellj < patch->totalDofsPerCell; cellj++) {
1843: const PetscInt col = dofsArray[(offset + cell0)*patch->totalDofsPerCell + cellj];
1844: const PetscInt key = row*rsize + col;
1845: if (col < 0) continue;
1846: if (!PetscBTLookupSet(bt, key)) ++dnnz[row];
1847: }
1848: }
1849: }
1850: }
1851: PetscBTDestroy(&bt);
1852: MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL);
1853: PetscFree(dnnz);
1855: PetscCalloc1(patch->totalDofsPerCell*patch->totalDofsPerCell, &zeroes);
1856: for (c = 0; c < ncell; ++c) {
1857: const PetscInt *idx = &dofsArray[(offset + c)*patch->totalDofsPerCell];
1858: MatSetValues(*mat, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, zeroes, INSERT_VALUES);
1859: }
1860: MatGetLocalSize(*mat, &rows, NULL);
1861: for (i = 0; i < rows; ++i) {
1862: MatSetValues(*mat, 1, &i, 1, &i, zeroes, INSERT_VALUES);
1863: }
1865: if (patch->usercomputeopintfacet) {
1866: const PetscInt *intFacetsArray = NULL;
1867: PetscInt i, numIntFacets, intFacetOffset;
1868: const PetscInt *facetCells = NULL;
1870: PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1871: PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1872: ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1873: ISGetIndices(patch->intFacets, &intFacetsArray);
1874: for (i = 0; i < numIntFacets; i++) {
1875: const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0];
1876: const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1];
1877: const PetscInt *cell0idx = &dofsArray[(offset + cell0)*patch->totalDofsPerCell];
1878: const PetscInt *cell1idx = &dofsArray[(offset + cell1)*patch->totalDofsPerCell];
1879: MatSetValues(*mat, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, zeroes, INSERT_VALUES);
1880: MatSetValues(*mat, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, zeroes, INSERT_VALUES);
1881: }
1882: }
1884: MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
1885: MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);
1887: PetscFree(zeroes);
1889: } else { /* rsize too big, use MATPREALLOCATOR */
1890: Mat preallocator;
1891: PetscScalar* vals;
1893: PetscCalloc1(patch->totalDofsPerCell*patch->totalDofsPerCell, &vals);
1894: MatCreate(PETSC_COMM_SELF, &preallocator);
1895: MatSetType(preallocator, MATPREALLOCATOR);
1896: MatSetSizes(preallocator, rsize, rsize, rsize, rsize);
1897: MatSetUp(preallocator);
1899: for (c = 0; c < ncell; ++c) {
1900: const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell;
1901: MatSetValues(preallocator, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, vals, INSERT_VALUES);
1902: }
1904: if (patch->usercomputeopintfacet) {
1905: const PetscInt *intFacetsArray = NULL;
1906: PetscInt i, numIntFacets, intFacetOffset;
1907: const PetscInt *facetCells = NULL;
1909: PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
1910: PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
1911: ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
1912: ISGetIndices(patch->intFacets, &intFacetsArray);
1913: for (i = 0; i < numIntFacets; i++) {
1914: const PetscInt cell0 = facetCells[2*(intFacetOffset + i) + 0];
1915: const PetscInt cell1 = facetCells[2*(intFacetOffset + i) + 1];
1916: const PetscInt *cell0idx = &dofsArray[(offset + cell0)*patch->totalDofsPerCell];
1917: const PetscInt *cell1idx = &dofsArray[(offset + cell1)*patch->totalDofsPerCell];
1918: MatSetValues(preallocator, patch->totalDofsPerCell, cell0idx, patch->totalDofsPerCell, cell1idx, vals, INSERT_VALUES);
1919: MatSetValues(preallocator, patch->totalDofsPerCell, cell1idx, patch->totalDofsPerCell, cell0idx, vals, INSERT_VALUES);
1920: }
1921: }
1923: PetscFree(vals);
1924: MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1925: MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1926: MatPreallocatorPreallocate(preallocator, PETSC_TRUE, *mat);
1927: MatDestroy(&preallocator);
1928: }
1929: PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0);
1930: if (withArtificial) {
1931: ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);
1932: } else {
1933: ISRestoreIndices(patch->dofs, &dofsArray);
1934: }
1935: }
1936: MatSetUp(*mat);
1937: return(0);
1938: }
1940: 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)
1941: {
1942: PC_PATCH *patch = (PC_PATCH *) pc->data;
1943: DM dm, plex;
1944: PetscSection s;
1945: const PetscInt *parray, *oarray;
1946: PetscInt Nf = patch->nsubspaces, Np, poff, p, f;
1947: PetscErrorCode ierr;
1950: if (patch->precomputeElementTensors) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Precomputing element tensors not implemented with DMPlex compute operator");
1951: PCGetDM(pc, &dm);
1952: DMConvert(dm, DMPLEX, &plex);
1953: dm = plex;
1954: DMGetLocalSection(dm, &s);
1955: /* Set offset into patch */
1956: PetscSectionGetDof(patch->pointCounts, patchNum, &Np);
1957: PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);
1958: ISGetIndices(patch->points, &parray);
1959: ISGetIndices(patch->offs, &oarray);
1960: for (f = 0; f < Nf; ++f) {
1961: for (p = 0; p < Np; ++p) {
1962: const PetscInt point = parray[poff+p];
1963: PetscInt dof;
1965: PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);
1966: PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);
1967: if (patch->nsubspaces == 1) {PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);}
1968: else {PetscSectionSetOffset(patch->patchSection, point, -1);}
1969: }
1970: }
1971: ISRestoreIndices(patch->points, &parray);
1972: ISRestoreIndices(patch->offs, &oarray);
1973: if (patch->viewSection) {ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);}
1974: DMPlexComputeResidual_Patch_Internal(dm, patch->patchSection, cellIS, 0.0, x, NULL, F, ctx);
1975: DMDestroy(&dm);
1976: return(0);
1977: }
1979: PetscErrorCode PCPatchComputeFunction_Internal(PC pc, Vec x, Vec F, PetscInt point)
1980: {
1981: PC_PATCH *patch = (PC_PATCH *) pc->data;
1982: const PetscInt *dofsArray;
1983: const PetscInt *dofsArrayWithAll;
1984: const PetscInt *cellsArray;
1985: PetscInt ncell, offset, pStart, pEnd;
1986: PetscErrorCode ierr;
1989: PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);
1990: if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n");
1991: ISGetIndices(patch->dofs, &dofsArray);
1992: ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);
1993: ISGetIndices(patch->cells, &cellsArray);
1994: PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);
1996: point += pStart;
1997: if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);
1999: PetscSectionGetDof(patch->cellCounts, point, &ncell);
2000: PetscSectionGetOffset(patch->cellCounts, point, &offset);
2001: if (ncell <= 0) {
2002: PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2003: return(0);
2004: }
2005: VecSet(F, 0.0);
2006: PetscStackPush("PCPatch user callback");
2007: /* Cannot reuse the same IS because the geometry info is being cached in it */
2008: ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);
2009: patch->usercomputef(pc, point, x, F, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell,
2010: dofsArrayWithAll + offset*patch->totalDofsPerCell,
2011: patch->usercomputefctx);
2012: PetscStackPop;
2013: ISDestroy(&patch->cellIS);
2014: ISRestoreIndices(patch->dofs, &dofsArray);
2015: ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);
2016: ISRestoreIndices(patch->cells, &cellsArray);
2017: if (patch->viewMatrix) {
2018: char name[PETSC_MAX_PATH_LEN];
2020: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch vector for Point %D", point);
2021: PetscObjectSetName((PetscObject) F, name);
2022: ObjectView((PetscObject) F, patch->viewerMatrix, patch->formatMatrix);
2023: }
2024: PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2025: return(0);
2026: }
2028: 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)
2029: {
2030: PC_PATCH *patch = (PC_PATCH *) pc->data;
2031: DM dm, plex;
2032: PetscSection s;
2033: const PetscInt *parray, *oarray;
2034: PetscInt Nf = patch->nsubspaces, Np, poff, p, f;
2035: PetscErrorCode ierr;
2038: PCGetDM(pc, &dm);
2039: DMConvert(dm, DMPLEX, &plex);
2040: dm = plex;
2041: DMGetLocalSection(dm, &s);
2042: /* Set offset into patch */
2043: PetscSectionGetDof(patch->pointCounts, patchNum, &Np);
2044: PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);
2045: ISGetIndices(patch->points, &parray);
2046: ISGetIndices(patch->offs, &oarray);
2047: for (f = 0; f < Nf; ++f) {
2048: for (p = 0; p < Np; ++p) {
2049: const PetscInt point = parray[poff+p];
2050: PetscInt dof;
2052: PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);
2053: PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);
2054: if (patch->nsubspaces == 1) {PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);}
2055: else {PetscSectionSetOffset(patch->patchSection, point, -1);}
2056: }
2057: }
2058: ISRestoreIndices(patch->points, &parray);
2059: ISRestoreIndices(patch->offs, &oarray);
2060: if (patch->viewSection) {ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);}
2061: /* TODO Shut off MatViewFromOptions() in MatAssemblyEnd() here */
2062: DMPlexComputeJacobian_Patch_Internal(dm, patch->patchSection, patch->patchSection, cellIS, 0.0, 0.0, x, NULL, J, J, ctx);
2063: DMDestroy(&dm);
2064: return(0);
2065: }
2067: /* This function zeros mat on entry */
2068: PetscErrorCode PCPatchComputeOperator_Internal(PC pc, Vec x, Mat mat, PetscInt point, PetscBool withArtificial)
2069: {
2070: PC_PATCH *patch = (PC_PATCH *) pc->data;
2071: const PetscInt *dofsArray;
2072: const PetscInt *dofsArrayWithAll = NULL;
2073: const PetscInt *cellsArray;
2074: PetscInt ncell, offset, pStart, pEnd, numIntFacets, intFacetOffset;
2075: PetscBool isNonlinear;
2076: PetscErrorCode ierr;
2079: PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);
2080: isNonlinear = patch->isNonlinear;
2081: if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n");
2082: if (withArtificial) {
2083: ISGetIndices(patch->dofsWithArtificial, &dofsArray);
2084: } else {
2085: ISGetIndices(patch->dofs, &dofsArray);
2086: }
2087: if (isNonlinear) {
2088: ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);
2089: }
2090: ISGetIndices(patch->cells, &cellsArray);
2091: PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);
2093: point += pStart;
2094: if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);
2096: PetscSectionGetDof(patch->cellCounts, point, &ncell);
2097: PetscSectionGetOffset(patch->cellCounts, point, &offset);
2098: if (ncell <= 0) {
2099: PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2100: return(0);
2101: }
2102: MatZeroEntries(mat);
2103: if (patch->precomputeElementTensors) {
2104: PetscInt i;
2105: PetscInt ndof = patch->totalDofsPerCell;
2106: const PetscScalar *elementTensors;
2108: VecGetArrayRead(patch->cellMats, &elementTensors);
2109: for (i = 0; i < ncell; i++) {
2110: const PetscInt cell = cellsArray[i + offset];
2111: const PetscInt *idx = dofsArray + (offset + i)*ndof;
2112: const PetscScalar *v = elementTensors + patch->precomputedTensorLocations[cell]*ndof*ndof;
2113: MatSetValues(mat, ndof, idx, ndof, idx, v, ADD_VALUES);
2114: }
2115: VecRestoreArrayRead(patch->cellMats, &elementTensors);
2116: MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
2117: MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
2118: } else {
2119: PetscStackPush("PCPatch user callback");
2120: /* Cannot reuse the same IS because the geometry info is being cached in it */
2121: ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);
2122: patch->usercomputeop(pc, point, x, mat, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell, dofsArrayWithAll ? dofsArrayWithAll + offset*patch->totalDofsPerCell : NULL, patch->usercomputeopctx);
2123: }
2124: if (patch->usercomputeopintfacet) {
2125: PetscSectionGetDof(patch->intFacetCounts, point, &numIntFacets);
2126: PetscSectionGetOffset(patch->intFacetCounts, point, &intFacetOffset);
2127: if (numIntFacets > 0) {
2128: /* For each interior facet, grab the two cells (in local numbering, and concatenate dof numberings for those cells) */
2129: PetscInt *facetDofs = NULL, *facetDofsWithAll = NULL;
2130: const PetscInt *intFacetsArray = NULL;
2131: PetscInt idx = 0;
2132: PetscInt i, c, d;
2133: PetscInt fStart;
2134: DM dm, plex;
2135: IS facetIS = NULL;
2136: const PetscInt *facetCells = NULL;
2138: ISGetIndices(patch->intFacetsToPatchCell, &facetCells);
2139: ISGetIndices(patch->intFacets, &intFacetsArray);
2140: PCGetDM(pc, &dm);
2141: DMConvert(dm, DMPLEX, &plex);
2142: dm = plex;
2143: DMPlexGetHeightStratum(dm, 1, &fStart, NULL);
2144: /* FIXME: Pull this malloc out. */
2145: PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofs);
2146: if (dofsArrayWithAll) {
2147: PetscMalloc1(2 * patch->totalDofsPerCell * numIntFacets, &facetDofsWithAll);
2148: }
2149: if (patch->precomputeElementTensors) {
2150: PetscInt nFacetDof = 2*patch->totalDofsPerCell;
2151: const PetscScalar *elementTensors;
2153: VecGetArrayRead(patch->intFacetMats, &elementTensors);
2155: for (i = 0; i < numIntFacets; i++) {
2156: const PetscInt facet = intFacetsArray[i + intFacetOffset];
2157: const PetscScalar *v = elementTensors + patch->precomputedIntFacetTensorLocations[facet - fStart]*nFacetDof*nFacetDof;
2158: idx = 0;
2159: /*
2160: * 0--1
2161: * |\-|
2162: * |+\|
2163: * 2--3
2164: * [0, 2, 3, 0, 1, 3]
2165: */
2166: for (c = 0; c < 2; c++) {
2167: const PetscInt cell = facetCells[2*(intFacetOffset + i) + c];
2168: for (d = 0; d < patch->totalDofsPerCell; d++) {
2169: facetDofs[idx] = dofsArray[(offset + cell)*patch->totalDofsPerCell + d];
2170: idx++;
2171: }
2172: }
2173: MatSetValues(mat, nFacetDof, facetDofs, nFacetDof, facetDofs, v, ADD_VALUES);
2174: }
2175: VecRestoreArrayRead(patch->intFacetMats, &elementTensors);
2176: } else {
2177: /*
2178: * 0--1
2179: * |\-|
2180: * |+\|
2181: * 2--3
2182: * [0, 2, 3, 0, 1, 3]
2183: */
2184: for (i = 0; i < numIntFacets; i++) {
2185: for (c = 0; c < 2; c++) {
2186: const PetscInt cell = facetCells[2*(intFacetOffset + i) + c];
2187: for (d = 0; d < patch->totalDofsPerCell; d++) {
2188: facetDofs[idx] = dofsArray[(offset + cell)*patch->totalDofsPerCell + d];
2189: if (dofsArrayWithAll) {
2190: facetDofsWithAll[idx] = dofsArrayWithAll[(offset + cell)*patch->totalDofsPerCell + d];
2191: }
2192: idx++;
2193: }
2194: }
2195: }
2196: ISCreateGeneral(PETSC_COMM_SELF, numIntFacets, intFacetsArray + intFacetOffset, PETSC_USE_POINTER, &facetIS);
2197: patch->usercomputeopintfacet(pc, point, x, mat, facetIS, 2*numIntFacets*patch->totalDofsPerCell, facetDofs, facetDofsWithAll, patch->usercomputeopintfacetctx);
2198: ISDestroy(&facetIS);
2199: }
2200: ISRestoreIndices(patch->intFacetsToPatchCell, &facetCells);
2201: ISRestoreIndices(patch->intFacets, &intFacetsArray);
2202: PetscFree(facetDofs);
2203: PetscFree(facetDofsWithAll);
2204: DMDestroy(&dm);
2205: }
2206: }
2208: MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
2209: MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
2211: if (!(withArtificial || isNonlinear) && patch->denseinverse) {
2212: MatFactorInfo info;
2213: PetscBool flg;
2214: PetscObjectTypeCompare((PetscObject)mat, MATSEQDENSE, &flg);
2215: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Invalid Mat type for dense inverse");
2216: MatFactorInfoInitialize(&info);
2217: MatLUFactor(mat, NULL, NULL, &info);
2218: MatSeqDenseInvertFactors_Private(mat);
2219: }
2220: PetscStackPop;
2221: ISDestroy(&patch->cellIS);
2222: if (withArtificial) {
2223: ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);
2224: } else {
2225: ISRestoreIndices(patch->dofs, &dofsArray);
2226: }
2227: if (isNonlinear) {
2228: ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);
2229: }
2230: ISRestoreIndices(patch->cells, &cellsArray);
2231: if (patch->viewMatrix) {
2232: char name[PETSC_MAX_PATH_LEN];
2234: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch matrix for Point %D", point);
2235: PetscObjectSetName((PetscObject) mat, name);
2236: ObjectView((PetscObject) mat, patch->viewerMatrix, patch->formatMatrix);
2237: }
2238: PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2239: return(0);
2240: }
2242: static PetscErrorCode MatSetValues_PCPatch_Private(Mat mat, PetscInt m, const PetscInt idxm[],
2243: PetscInt n, const PetscInt idxn[], const PetscScalar *v, InsertMode addv)
2244: {
2245: Vec data;
2246: PetscScalar *array;
2247: PetscInt bs, nz, i, j, cell;
2250: MatShellGetContext(mat, &data);
2251: VecGetBlockSize(data, &bs);
2252: VecGetSize(data, &nz);
2253: VecGetArray(data, &array);
2254: if (m != n) SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Only for square insertion");
2255: cell = (PetscInt)(idxm[0]/bs); /* use the fact that this is called once per cell */
2256: for (i = 0; i < m; i++) {
2257: if (idxm[i] != idxn[i]) SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Row and column indices must match!");
2258: for (j = 0; j < n; j++) {
2259: const PetscScalar v_ = v[i*bs + j];
2260: /* Indexing is special to the data structure we have! */
2261: if (addv == INSERT_VALUES) {
2262: array[cell*bs*bs + i*bs + j] = v_;
2263: } else {
2264: array[cell*bs*bs + i*bs + j] += v_;
2265: }
2266: }
2267: }
2268: VecRestoreArray(data, &array);
2269: return(0);
2270: }
2272: static PetscErrorCode PCPatchPrecomputePatchTensors_Private(PC pc)
2273: {
2274: PC_PATCH *patch = (PC_PATCH *)pc->data;
2275: const PetscInt *cellsArray;
2276: PetscInt ncell, offset;
2277: const PetscInt *dofMapArray;
2278: PetscInt i, j;
2279: IS dofMap;
2280: IS cellIS;
2281: const PetscInt ndof = patch->totalDofsPerCell;
2282: PetscErrorCode ierr;
2283: Mat vecMat;
2284: PetscInt cStart, cEnd;
2285: DM dm, plex;
2288: ISGetSize(patch->cells, &ncell);
2289: if (!ncell) { /* No cells to assemble over -> skip */
2290: return(0);
2291: }
2293: PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);
2295: PCGetDM(pc, &dm);
2296: DMConvert(dm, DMPLEX, &plex);
2297: dm = plex;
2298: if (!patch->allCells) {
2299: PetscHSetI cells;
2300: PetscHashIter hi;
2301: PetscInt pStart, pEnd;
2302: PetscInt *allCells = NULL;
2303: PetscHSetICreate(&cells);
2304: ISGetIndices(patch->cells, &cellsArray);
2305: PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);
2306: for (i = pStart; i < pEnd; i++) {
2307: PetscSectionGetDof(patch->cellCounts, i, &ncell);
2308: PetscSectionGetOffset(patch->cellCounts, i, &offset);
2309: if (ncell <= 0) continue;
2310: for (j = 0; j < ncell; j++) {
2311: PetscHSetIAdd(cells, cellsArray[offset + j]);
2312: }
2313: }
2314: ISRestoreIndices(patch->cells, &cellsArray);
2315: PetscHSetIGetSize(cells, &ncell);
2316: PetscMalloc1(ncell, &allCells);
2317: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2318: PetscMalloc1(cEnd-cStart, &patch->precomputedTensorLocations);
2319: i = 0;
2320: PetscHashIterBegin(cells, hi);
2321: while (!PetscHashIterAtEnd(cells, hi)) {
2322: PetscHashIterGetKey(cells, hi, allCells[i]);
2323: patch->precomputedTensorLocations[allCells[i]] = i;
2324: PetscHashIterNext(cells, hi);
2325: i++;
2326: }
2327: PetscHSetIDestroy(&cells);
2328: ISCreateGeneral(PETSC_COMM_SELF, ncell, allCells, PETSC_OWN_POINTER, &patch->allCells);
2329: }
2330: ISGetSize(patch->allCells, &ncell);
2331: if (!patch->cellMats) {
2332: VecCreateSeq(PETSC_COMM_SELF, ncell*ndof*ndof, &patch->cellMats);
2333: VecSetBlockSize(patch->cellMats, ndof);
2334: }
2335: VecSet(patch->cellMats, 0);
2337: MatCreateShell(PETSC_COMM_SELF, ncell*ndof, ncell*ndof, ncell*ndof, ncell*ndof,
2338: (void*)patch->cellMats, &vecMat);
2339: MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void(*)(void))&MatSetValues_PCPatch_Private);
2340: ISGetSize(patch->allCells, &ncell);
2341: ISCreateStride(PETSC_COMM_SELF, ndof*ncell, 0, 1, &dofMap);
2342: ISGetIndices(dofMap, &dofMapArray);
2343: ISGetIndices(patch->allCells, &cellsArray);
2344: ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray, PETSC_USE_POINTER, &cellIS);
2345: PetscStackPush("PCPatch user callback");
2346: /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2347: patch->usercomputeop(pc, -1, NULL, vecMat, cellIS, ndof*ncell, dofMapArray, NULL, patch->usercomputeopctx);
2348: PetscStackPop;
2349: ISDestroy(&cellIS);
2350: MatDestroy(&vecMat);
2351: ISRestoreIndices(patch->allCells, &cellsArray);
2352: ISRestoreIndices(dofMap, &dofMapArray);
2353: ISDestroy(&dofMap);
2355: if (patch->usercomputeopintfacet) {
2356: PetscInt nIntFacets;
2357: IS intFacetsIS;
2358: const PetscInt *intFacetsArray = NULL;
2359: if (!patch->allIntFacets) {
2360: PetscHSetI facets;
2361: PetscHashIter hi;
2362: PetscInt pStart, pEnd, fStart, fEnd;
2363: PetscInt *allIntFacets = NULL;
2364: PetscHSetICreate(&facets);
2365: ISGetIndices(patch->intFacets, &intFacetsArray);
2366: PetscSectionGetChart(patch->intFacetCounts, &pStart, &pEnd);
2367: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2368: for (i = pStart; i < pEnd; i++) {
2369: PetscSectionGetDof(patch->intFacetCounts, i, &nIntFacets);
2370: PetscSectionGetOffset(patch->intFacetCounts, i, &offset);
2371: if (nIntFacets <= 0) continue;
2372: for (j = 0; j < nIntFacets; j++) {
2373: PetscHSetIAdd(facets, intFacetsArray[offset + j]);
2374: }
2375: }
2376: ISRestoreIndices(patch->intFacets, &intFacetsArray);
2377: PetscHSetIGetSize(facets, &nIntFacets);
2378: PetscMalloc1(nIntFacets, &allIntFacets);
2379: PetscMalloc1(fEnd-fStart, &patch->precomputedIntFacetTensorLocations);
2380: i = 0;
2381: PetscHashIterBegin(facets, hi);
2382: while (!PetscHashIterAtEnd(facets, hi)) {
2383: PetscHashIterGetKey(facets, hi, allIntFacets[i]);
2384: patch->precomputedIntFacetTensorLocations[allIntFacets[i] - fStart] = i;
2385: PetscHashIterNext(facets, hi);
2386: i++;
2387: }
2388: PetscHSetIDestroy(&facets);
2389: ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, allIntFacets, PETSC_OWN_POINTER, &patch->allIntFacets);
2390: }
2391: ISGetSize(patch->allIntFacets, &nIntFacets);
2392: if (!patch->intFacetMats) {
2393: VecCreateSeq(PETSC_COMM_SELF, nIntFacets*ndof*ndof*4, &patch->intFacetMats);
2394: VecSetBlockSize(patch->intFacetMats, ndof*2);
2395: }
2396: VecSet(patch->intFacetMats, 0);
2398: MatCreateShell(PETSC_COMM_SELF, nIntFacets*ndof*2, nIntFacets*ndof*2, nIntFacets*ndof*2, nIntFacets*ndof*2,
2399: (void*)patch->intFacetMats, &vecMat);
2400: MatShellSetOperation(vecMat, MATOP_SET_VALUES, (void(*)(void))&MatSetValues_PCPatch_Private);
2401: ISCreateStride(PETSC_COMM_SELF, 2*ndof*nIntFacets, 0, 1, &dofMap);
2402: ISGetIndices(dofMap, &dofMapArray);
2403: ISGetIndices(patch->allIntFacets, &intFacetsArray);
2404: ISCreateGeneral(PETSC_COMM_SELF, nIntFacets, intFacetsArray, PETSC_USE_POINTER, &intFacetsIS);
2405: PetscStackPush("PCPatch user callback (interior facets)");
2406: /* TODO: Fix for DMPlex compute op, this bypasses a lot of the machinery and just assembles every element tensor. */
2407: patch->usercomputeopintfacet(pc, -1, NULL, vecMat, intFacetsIS, 2*ndof*nIntFacets, dofMapArray, NULL, patch->usercomputeopintfacetctx);
2408: PetscStackPop;
2409: ISDestroy(&intFacetsIS);
2410: MatDestroy(&vecMat);
2411: ISRestoreIndices(patch->allIntFacets, &intFacetsArray);
2412: ISRestoreIndices(dofMap, &dofMapArray);
2413: ISDestroy(&dofMap);
2414: }
2415: DMDestroy(&dm);
2416: PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);
2418: return(0);
2419: }
2421: PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat, PatchScatterType scattertype)
2422: {
2423: PC_PATCH *patch = (PC_PATCH *) pc->data;
2424: const PetscScalar *xArray = NULL;
2425: PetscScalar *yArray = NULL;
2426: const PetscInt *gtolArray = NULL;
2427: PetscInt dof, offset, lidx;
2428: PetscErrorCode ierr;
2431: VecGetArrayRead(x, &xArray);
2432: VecGetArray(y, &yArray);
2433: if (scattertype == SCATTER_WITHARTIFICIAL) {
2434: PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);
2435: PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offset);
2436: ISGetIndices(patch->gtolWithArtificial, >olArray);
2437: } else if (scattertype == SCATTER_WITHALL) {
2438: PetscSectionGetDof(patch->gtolCountsWithAll, p, &dof);
2439: PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offset);
2440: ISGetIndices(patch->gtolWithAll, >olArray);
2441: } else {
2442: PetscSectionGetDof(patch->gtolCounts, p, &dof);
2443: PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2444: ISGetIndices(patch->gtol, >olArray);
2445: }
2446: if (mode == INSERT_VALUES && scat != SCATTER_FORWARD) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't insert if not scattering forward\n");
2447: if (mode == ADD_VALUES && scat != SCATTER_REVERSE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't add if not scattering reverse\n");
2448: for (lidx = 0; lidx < dof; ++lidx) {
2449: const PetscInt gidx = gtolArray[offset+lidx];
2451: if (mode == INSERT_VALUES) yArray[lidx] = xArray[gidx]; /* Forward */
2452: else yArray[gidx] += xArray[lidx]; /* Reverse */
2453: }
2454: if (scattertype == SCATTER_WITHARTIFICIAL) {
2455: ISRestoreIndices(patch->gtolWithArtificial, >olArray);
2456: } else if (scattertype == SCATTER_WITHALL) {
2457: ISRestoreIndices(patch->gtolWithAll, >olArray);
2458: } else {
2459: ISRestoreIndices(patch->gtol, >olArray);
2460: }
2461: VecRestoreArrayRead(x, &xArray);
2462: VecRestoreArray(y, &yArray);
2463: return(0);
2464: }
2466: static PetscErrorCode PCSetUp_PATCH_Linear(PC pc)
2467: {
2468: PC_PATCH *patch = (PC_PATCH *) pc->data;
2469: const char *prefix;
2470: PetscInt i;
2474: if (!pc->setupcalled) {
2475: if (!patch->save_operators && patch->denseinverse) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Can't have dense inverse without save operators");
2476: if (!patch->denseinverse) {
2477: PetscMalloc1(patch->npatch, &patch->solver);
2478: PCGetOptionsPrefix(pc, &prefix);
2479: for (i = 0; i < patch->npatch; ++i) {
2480: KSP ksp;
2481: PC subpc;
2483: KSPCreate(PETSC_COMM_SELF, &ksp);
2484: KSPSetErrorIfNotConverged(ksp, pc->erroriffailure);
2485: KSPSetOptionsPrefix(ksp, prefix);
2486: KSPAppendOptionsPrefix(ksp, "sub_");
2487: PetscObjectIncrementTabLevel((PetscObject) ksp, (PetscObject) pc, 1);
2488: KSPGetPC(ksp, &subpc);
2489: PetscObjectIncrementTabLevel((PetscObject) subpc, (PetscObject) pc, 1);
2490: PetscLogObjectParent((PetscObject) pc, (PetscObject) ksp);
2491: patch->solver[i] = (PetscObject) ksp;
2492: }
2493: }
2494: }
2495: if (patch->save_operators) {
2496: if (patch->precomputeElementTensors) {
2497: PCPatchPrecomputePatchTensors_Private(pc);
2498: }
2499: for (i = 0; i < patch->npatch; ++i) {
2500: PCPatchComputeOperator_Internal(pc, NULL, patch->mat[i], i, PETSC_FALSE);
2501: if (!patch->denseinverse) {
2502: KSPSetOperators((KSP) patch->solver[i], patch->mat[i], patch->mat[i]);
2503: } else if (patch->mat[i] && !patch->densesolve) {
2504: /* Setup matmult callback */
2505: MatGetOperation(patch->mat[i], MATOP_MULT, (void (**)(void))&patch->densesolve);
2506: }
2507: }
2508: }
2509: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2510: for (i = 0; i < patch->npatch; ++i) {
2511: /* Instead of padding patch->patchUpdate with zeros to get */
2512: /* patch->patchUpdateWithArtificial and then multiplying with the matrix, */
2513: /* just get rid of the columns that correspond to the dofs with */
2514: /* artificial bcs. That's of course fairly inefficient, hopefully we */
2515: /* can just assemble the rectangular matrix in the first place. */
2516: Mat matSquare;
2517: IS rowis;
2518: PetscInt dof;
2520: MatGetSize(patch->mat[i], &dof, NULL);
2521: if (dof == 0) {
2522: patch->matWithArtificial[i] = NULL;
2523: continue;
2524: }
2526: PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);
2527: PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);
2529: MatGetSize(matSquare, &dof, NULL);
2530: ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis);
2531: if (pc->setupcalled) {
2532: MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_REUSE_MATRIX, &patch->matWithArtificial[i]);
2533: } else {
2534: MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &patch->matWithArtificial[i]);
2535: }
2536: ISDestroy(&rowis);
2537: MatDestroy(&matSquare);
2538: }
2539: }
2540: return(0);
2541: }
2543: static PetscErrorCode PCSetUp_PATCH(PC pc)
2544: {
2545: PC_PATCH *patch = (PC_PATCH *) pc->data;
2546: PetscInt i;
2547: PetscBool isNonlinear;
2548: PetscInt maxDof = -1, maxDofWithArtificial = -1;
2552: if (!pc->setupcalled) {
2553: PetscInt pStart, pEnd, p;
2554: PetscInt localSize;
2556: PetscLogEventBegin(PC_Patch_CreatePatches, pc, 0, 0, 0);
2558: isNonlinear = patch->isNonlinear;
2559: if (!patch->nsubspaces) {
2560: DM dm, plex;
2561: PetscSection s;
2562: PetscInt cStart, cEnd, c, Nf, f, numGlobalBcs = 0, *globalBcs, *Nb, totNb = 0, **cellDofs;
2564: PCGetDM(pc, &dm);
2565: if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "Must set DM for PCPATCH or call PCPatchSetDiscretisationInfo()");
2566: DMConvert(dm, DMPLEX, &plex);
2567: dm = plex;
2568: DMGetLocalSection(dm, &s);
2569: PetscSectionGetNumFields(s, &Nf);
2570: PetscSectionGetChart(s, &pStart, &pEnd);
2571: for (p = pStart; p < pEnd; ++p) {
2572: PetscInt cdof;
2573: PetscSectionGetConstraintDof(s, p, &cdof);
2574: numGlobalBcs += cdof;
2575: }
2576: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2577: PetscMalloc3(Nf, &Nb, Nf, &cellDofs, numGlobalBcs, &globalBcs);
2578: for (f = 0; f < Nf; ++f) {
2579: PetscFE fe;
2580: PetscDualSpace sp;
2581: PetscInt cdoff = 0;
2583: DMGetField(dm, f, NULL, (PetscObject *) &fe);
2584: /* PetscFEGetNumComponents(fe, &Nc[f]); */
2585: PetscFEGetDualSpace(fe, &sp);
2586: PetscDualSpaceGetDimension(sp, &Nb[f]);
2587: totNb += Nb[f];
2589: PetscMalloc1((cEnd-cStart)*Nb[f], &cellDofs[f]);
2590: for (c = cStart; c < cEnd; ++c) {
2591: PetscInt *closure = NULL;
2592: PetscInt clSize = 0, cl;
2594: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
2595: for (cl = 0; cl < clSize*2; cl += 2) {
2596: const PetscInt p = closure[cl];
2597: PetscInt fdof, d, foff;
2599: PetscSectionGetFieldDof(s, p, f, &fdof);
2600: PetscSectionGetFieldOffset(s, p, f, &foff);
2601: for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d;
2602: }
2603: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
2604: }
2605: 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]);
2606: }
2607: numGlobalBcs = 0;
2608: for (p = pStart; p < pEnd; ++p) {
2609: const PetscInt *ind;
2610: PetscInt off, cdof, d;
2612: PetscSectionGetOffset(s, p, &off);
2613: PetscSectionGetConstraintDof(s, p, &cdof);
2614: PetscSectionGetConstraintIndices(s, p, &ind);
2615: for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d];
2616: }
2618: PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **) cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs);
2619: for (f = 0; f < Nf; ++f) {
2620: PetscFree(cellDofs[f]);
2621: }
2622: PetscFree3(Nb, cellDofs, globalBcs);
2623: PCPatchSetComputeFunction(pc, PCPatchComputeFunction_DMPlex_Private, NULL);
2624: PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL);
2625: DMDestroy(&dm);
2626: }
2628: localSize = patch->subspaceOffsets[patch->nsubspaces];
2629: VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localRHS);
2630: VecSetUp(patch->localRHS);
2631: VecDuplicate(patch->localRHS, &patch->localUpdate);
2632: PCPatchCreateCellPatches(pc);
2633: PCPatchCreateCellPatchDiscretisationInfo(pc);
2635: /* OK, now build the work vectors */
2636: PetscSectionGetChart(patch->gtolCounts, &pStart, &pEnd);
2638: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2639: PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithArtificial);
2640: }
2641: if (isNonlinear) {
2642: PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithAll);
2643: }
2644: for (p = pStart; p < pEnd; ++p) {
2645: PetscInt dof;
2647: PetscSectionGetDof(patch->gtolCounts, p, &dof);
2648: maxDof = PetscMax(maxDof, dof);
2649: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2650: const PetscInt *gtolArray, *gtolArrayWithArtificial = NULL;
2651: PetscInt numPatchDofs, offset;
2652: PetscInt numPatchDofsWithArtificial, offsetWithArtificial;
2653: PetscInt dofWithoutArtificialCounter = 0;
2654: PetscInt *patchWithoutArtificialToWithArtificialArray;
2656: PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);
2657: maxDofWithArtificial = PetscMax(maxDofWithArtificial, dof);
2659: /* Now build the mapping that for a dof in a patch WITHOUT dofs that have artificial bcs gives the */
2660: /* the index in the patch with all dofs */
2661: ISGetIndices(patch->gtol, >olArray);
2663: PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);
2664: if (numPatchDofs == 0) {
2665: patch->dofMappingWithoutToWithArtificial[p-pStart] = NULL;
2666: continue;
2667: }
2669: PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2670: ISGetIndices(patch->gtolWithArtificial, >olArrayWithArtificial);
2671: PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &numPatchDofsWithArtificial);
2672: PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offsetWithArtificial);
2674: PetscMalloc1(numPatchDofs, &patchWithoutArtificialToWithArtificialArray);
2675: for (i=0; i<numPatchDofsWithArtificial; i++) {
2676: if (gtolArrayWithArtificial[i+offsetWithArtificial] == gtolArray[offset+dofWithoutArtificialCounter]) {
2677: patchWithoutArtificialToWithArtificialArray[dofWithoutArtificialCounter] = i;
2678: dofWithoutArtificialCounter++;
2679: if (dofWithoutArtificialCounter == numPatchDofs)
2680: break;
2681: }
2682: }
2683: ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutArtificialToWithArtificialArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithArtificial[p-pStart]);
2684: ISRestoreIndices(patch->gtol, >olArray);
2685: ISRestoreIndices(patch->gtolWithArtificial, >olArrayWithArtificial);
2686: }
2687: if (isNonlinear) {
2688: const PetscInt *gtolArray, *gtolArrayWithAll = NULL;
2689: PetscInt numPatchDofs, offset;
2690: PetscInt numPatchDofsWithAll, offsetWithAll;
2691: PetscInt dofWithoutAllCounter = 0;
2692: PetscInt *patchWithoutAllToWithAllArray;
2694: /* Now build the mapping that for a dof in a patch WITHOUT dofs that have artificial bcs gives the */
2695: /* the index in the patch with all dofs */
2696: ISGetIndices(patch->gtol, >olArray);
2698: PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);
2699: if (numPatchDofs == 0) {
2700: patch->dofMappingWithoutToWithAll[p-pStart] = NULL;
2701: continue;
2702: }
2704: PetscSectionGetOffset(patch->gtolCounts, p, &offset);
2705: ISGetIndices(patch->gtolWithAll, >olArrayWithAll);
2706: PetscSectionGetDof(patch->gtolCountsWithAll, p, &numPatchDofsWithAll);
2707: PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offsetWithAll);
2709: PetscMalloc1(numPatchDofs, &patchWithoutAllToWithAllArray);
2711: for (i=0; i<numPatchDofsWithAll; i++) {
2712: if (gtolArrayWithAll[i+offsetWithAll] == gtolArray[offset+dofWithoutAllCounter]) {
2713: patchWithoutAllToWithAllArray[dofWithoutAllCounter] = i;
2714: dofWithoutAllCounter++;
2715: if (dofWithoutAllCounter == numPatchDofs)
2716: break;
2717: }
2718: }
2719: ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutAllToWithAllArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithAll[p-pStart]);
2720: ISRestoreIndices(patch->gtol, >olArray);
2721: ISRestoreIndices(patch->gtolWithAll, >olArrayWithAll);
2722: }
2723: }
2724: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2725: VecCreateSeq(PETSC_COMM_SELF, maxDofWithArtificial, &patch->patchRHSWithArtificial);
2726: VecSetUp(patch->patchRHSWithArtificial);
2727: }
2728: VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchRHS);
2729: VecSetUp(patch->patchRHS);
2730: VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchUpdate);
2731: VecSetUp(patch->patchUpdate);
2732: if (patch->save_operators) {
2733: PetscMalloc1(patch->npatch, &patch->mat);
2734: for (i = 0; i < patch->npatch; ++i) {
2735: PCPatchCreateMatrix_Private(pc, i, &patch->mat[i], PETSC_FALSE);
2736: }
2737: }
2738: PetscLogEventEnd(PC_Patch_CreatePatches, pc, 0, 0, 0);
2740: /* If desired, calculate weights for dof multiplicity */
2741: if (patch->partition_of_unity) {
2742: PetscScalar *input = NULL;
2743: PetscScalar *output = NULL;
2744: Vec global;
2746: VecDuplicate(patch->localRHS, &patch->dof_weights);
2747: if (patch->local_composition_type == PC_COMPOSITE_ADDITIVE) {
2748: for (i = 0; i < patch->npatch; ++i) {
2749: PetscInt dof;
2751: PetscSectionGetDof(patch->gtolCounts, i+pStart, &dof);
2752: if (dof <= 0) continue;
2753: VecSet(patch->patchRHS, 1.0);
2754: PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchRHS, patch->dof_weights, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);
2755: }
2756: } else {
2757: /* multiplicative is actually only locally multiplicative and globally additive. need the pou where the mesh decomposition overlaps */
2758: VecSet(patch->dof_weights, 1.0);
2759: }
2761: VecDuplicate(patch->dof_weights, &global);
2762: VecSet(global, 0.);
2764: VecGetArray(patch->dof_weights, &input);
2765: VecGetArray(global, &output);
2766: PetscSFReduceBegin(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM);
2767: PetscSFReduceEnd(patch->sectionSF, MPIU_SCALAR, input, output, MPI_SUM);
2768: VecRestoreArray(patch->dof_weights, &input);
2769: VecRestoreArray(global, &output);
2771: VecReciprocal(global);
2773: VecGetArray(patch->dof_weights, &output);
2774: VecGetArray(global, &input);
2775: PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, input, output);
2776: PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, input, output);
2777: VecRestoreArray(patch->dof_weights, &output);
2778: VecRestoreArray(global, &input);
2779: VecDestroy(&global);
2780: }
2781: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE && patch->save_operators) {
2782: PetscMalloc1(patch->npatch, &patch->matWithArtificial);
2783: }
2784: }
2785: (*patch->setupsolver)(pc);
2786: return(0);
2787: }
2789: static PetscErrorCode PCApply_PATCH_Linear(PC pc, PetscInt i, Vec x, Vec y)
2790: {
2791: PC_PATCH *patch = (PC_PATCH *) pc->data;
2792: KSP ksp;
2793: Mat op;
2794: PetscInt m, n;
2798: if (patch->denseinverse) {
2799: (*patch->densesolve)(patch->mat[i], x, y);
2800: return(0);
2801: }
2802: ksp = (KSP) patch->solver[i];
2803: if (!patch->save_operators) {
2804: Mat mat;
2806: PCPatchCreateMatrix_Private(pc, i, &mat, PETSC_FALSE);
2807: /* Populate operator here. */
2808: PCPatchComputeOperator_Internal(pc, NULL, mat, i, PETSC_FALSE);
2809: KSPSetOperators(ksp, mat, mat);
2810: /* Drop reference so the KSPSetOperators below will blow it away. */
2811: MatDestroy(&mat);
2812: }
2813: PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);
2814: if (!ksp->setfromoptionscalled) {
2815: KSPSetFromOptions(ksp);
2816: }
2817: /* Disgusting trick to reuse work vectors */
2818: KSPGetOperators(ksp, &op, NULL);
2819: MatGetLocalSize(op, &m, &n);
2820: x->map->n = m;
2821: y->map->n = n;
2822: x->map->N = m;
2823: y->map->N = n;
2824: KSPSolve(ksp, x, y);
2825: KSPCheckSolve(ksp, pc, y);
2826: PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);
2827: if (!patch->save_operators) {
2828: PC pc;
2829: KSPSetOperators(ksp, NULL, NULL);
2830: KSPGetPC(ksp, &pc);
2831: /* Destroy PC context too, otherwise the factored matrix hangs around. */
2832: PCReset(pc);
2833: }
2834: return(0);
2835: }
2837: static PetscErrorCode PCUpdateMultiplicative_PATCH_Linear(PC pc, PetscInt i, PetscInt pStart)
2838: {
2839: PC_PATCH *patch = (PC_PATCH *) pc->data;
2840: Mat multMat;
2841: PetscInt n, m;
2846: if (patch->save_operators) {
2847: multMat = patch->matWithArtificial[i];
2848: } else {
2849: /*Very inefficient, hopefully we can just assemble the rectangular matrix in the first place.*/
2850: Mat matSquare;
2851: PetscInt dof;
2852: IS rowis;
2853: PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);
2854: PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);
2855: MatGetSize(matSquare, &dof, NULL);
2856: ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis);
2857: MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &multMat);
2858: MatDestroy(&matSquare);
2859: ISDestroy(&rowis);
2860: }
2861: /* Disgusting trick to reuse work vectors */
2862: MatGetLocalSize(multMat, &m, &n);
2863: patch->patchUpdate->map->n = n;
2864: patch->patchRHSWithArtificial->map->n = m;
2865: patch->patchUpdate->map->N = n;
2866: patch->patchRHSWithArtificial->map->N = m;
2867: MatMult(multMat, patch->patchUpdate, patch->patchRHSWithArtificial);
2868: VecScale(patch->patchRHSWithArtificial, -1.0);
2869: PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHSWithArtificial, patch->localRHS, ADD_VALUES, SCATTER_REVERSE, SCATTER_WITHARTIFICIAL);
2870: if (!patch->save_operators) {
2871: MatDestroy(&multMat);
2872: }
2873: return(0);
2874: }
2876: static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y)
2877: {
2878: PC_PATCH *patch = (PC_PATCH *) pc->data;
2879: const PetscScalar *globalRHS = NULL;
2880: PetscScalar *localRHS = NULL;
2881: PetscScalar *globalUpdate = NULL;
2882: const PetscInt *bcNodes = NULL;
2883: PetscInt nsweep = patch->symmetrise_sweep ? 2 : 1;
2884: PetscInt start[2] = {0, 0};
2885: PetscInt end[2] = {-1, -1};
2886: const PetscInt inc[2] = {1, -1};
2887: const PetscScalar *localUpdate;
2888: const PetscInt *iterationSet;
2889: PetscInt pStart, numBcs, n, sweep, bc, j;
2890: PetscErrorCode ierr;
2893: PetscLogEventBegin(PC_Patch_Apply, pc, 0, 0, 0);
2894: PetscOptionsPushGetViewerOff(PETSC_TRUE);
2895: /* start, end, inc have 2 entries to manage a second backward sweep if we symmetrize */
2896: end[0] = patch->npatch;
2897: start[1] = patch->npatch-1;
2898: if (patch->user_patches) {
2899: ISGetLocalSize(patch->iterationSet, &end[0]);
2900: start[1] = end[0] - 1;
2901: ISGetIndices(patch->iterationSet, &iterationSet);
2902: }
2903: /* Scatter from global space into overlapped local spaces */
2904: VecGetArrayRead(x, &globalRHS);
2905: VecGetArray(patch->localRHS, &localRHS);
2906: PetscSFBcastBegin(patch->sectionSF, MPIU_SCALAR, globalRHS, localRHS);
2907: PetscSFBcastEnd(patch->sectionSF, MPIU_SCALAR, globalRHS, localRHS);
2908: VecRestoreArrayRead(x, &globalRHS);
2909: VecRestoreArray(patch->localRHS, &localRHS);
2911: VecSet(patch->localUpdate, 0.0);
2912: PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);
2913: PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);
2914: for (sweep = 0; sweep < nsweep; sweep++) {
2915: for (j = start[sweep]; j*inc[sweep] < end[sweep]*inc[sweep]; j += inc[sweep]) {
2916: PetscInt i = patch->user_patches ? iterationSet[j] : j;
2917: PetscInt start, len;
2919: PetscSectionGetDof(patch->gtolCounts, i+pStart, &len);
2920: PetscSectionGetOffset(patch->gtolCounts, i+pStart, &start);
2921: /* TODO: Squash out these guys in the setup as well. */
2922: if (len <= 0) continue;
2923: /* TODO: Do we need different scatters for X and Y? */
2924: PCPatch_ScatterLocal_Private(pc, i+pStart, patch->localRHS, patch->patchRHS, INSERT_VALUES, SCATTER_FORWARD, SCATTER_INTERIOR);
2925: (*patch->applysolver)(pc, i, patch->patchRHS, patch->patchUpdate);
2926: PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchUpdate, patch->localUpdate, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);
2927: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
2928: (*patch->updatemultiplicative)(pc, i, pStart);
2929: }
2930: }
2931: }
2932: PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);
2933: if (patch->user_patches) {ISRestoreIndices(patch->iterationSet, &iterationSet);}
2934: /* XXX: should we do this on the global vector? */
2935: if (patch->partition_of_unity) {
2936: VecPointwiseMult(patch->localUpdate, patch->localUpdate, patch->dof_weights);
2937: }
2938: /* Now patch->localUpdate contains the solution of the patch solves, so we need to combine them all. */
2939: VecSet(y, 0.0);
2940: VecGetArray(y, &globalUpdate);
2941: VecGetArrayRead(patch->localUpdate, &localUpdate);
2942: PetscSFReduceBegin(patch->sectionSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);
2943: PetscSFReduceEnd(patch->sectionSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);
2944: VecRestoreArrayRead(patch->localUpdate, &localUpdate);
2946: /* Now we need to send the global BC values through */
2947: VecGetArrayRead(x, &globalRHS);
2948: ISGetSize(patch->globalBcNodes, &numBcs);
2949: ISGetIndices(patch->globalBcNodes, &bcNodes);
2950: VecGetLocalSize(x, &n);
2951: for (bc = 0; bc < numBcs; ++bc) {
2952: const PetscInt idx = bcNodes[bc];
2953: if (idx < n) globalUpdate[idx] = globalRHS[idx];
2954: }
2956: ISRestoreIndices(patch->globalBcNodes, &bcNodes);
2957: VecRestoreArrayRead(x, &globalRHS);
2958: VecRestoreArray(y, &globalUpdate);
2960: PetscOptionsPopGetViewerOff();
2961: PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0);
2962: return(0);
2963: }
2965: static PetscErrorCode PCReset_PATCH_Linear(PC pc)
2966: {
2967: PC_PATCH *patch = (PC_PATCH *) pc->data;
2968: PetscInt i;
2972: if (patch->solver) {
2973: for (i = 0; i < patch->npatch; ++i) {KSPReset((KSP) patch->solver[i]);}
2974: }
2975: return(0);
2976: }
2978: static PetscErrorCode PCReset_PATCH(PC pc)
2979: {
2980: PC_PATCH *patch = (PC_PATCH *) pc->data;
2981: PetscInt i;
2986: PetscSFDestroy(&patch->sectionSF);
2987: PetscSectionDestroy(&patch->cellCounts);
2988: PetscSectionDestroy(&patch->pointCounts);
2989: PetscSectionDestroy(&patch->cellNumbering);
2990: PetscSectionDestroy(&patch->gtolCounts);
2991: ISDestroy(&patch->gtol);
2992: ISDestroy(&patch->cells);
2993: ISDestroy(&patch->points);
2994: ISDestroy(&patch->dofs);
2995: ISDestroy(&patch->offs);
2996: PetscSectionDestroy(&patch->patchSection);
2997: ISDestroy(&patch->ghostBcNodes);
2998: ISDestroy(&patch->globalBcNodes);
2999: PetscSectionDestroy(&patch->gtolCountsWithArtificial);
3000: ISDestroy(&patch->gtolWithArtificial);
3001: ISDestroy(&patch->dofsWithArtificial);
3002: ISDestroy(&patch->offsWithArtificial);
3003: PetscSectionDestroy(&patch->gtolCountsWithAll);
3004: ISDestroy(&patch->gtolWithAll);
3005: ISDestroy(&patch->dofsWithAll);
3006: ISDestroy(&patch->offsWithAll);
3007: VecDestroy(&patch->cellMats);
3008: VecDestroy(&patch->intFacetMats);
3009: ISDestroy(&patch->allCells);
3010: ISDestroy(&patch->intFacets);
3011: ISDestroy(&patch->extFacets);
3012: ISDestroy(&patch->intFacetsToPatchCell);
3013: PetscSectionDestroy(&patch->intFacetCounts);
3014: PetscSectionDestroy(&patch->extFacetCounts);
3016: if (patch->dofSection) for (i = 0; i < patch->nsubspaces; i++) {PetscSectionDestroy(&patch->dofSection[i]);}
3017: PetscFree(patch->dofSection);
3018: PetscFree(patch->bs);
3019: PetscFree(patch->nodesPerCell);
3020: if (patch->cellNodeMap) for (i = 0; i < patch->nsubspaces; i++) {PetscFree(patch->cellNodeMap[i]);}
3021: PetscFree(patch->cellNodeMap);
3022: PetscFree(patch->subspaceOffsets);
3024: (*patch->resetsolver)(pc);
3026: if (patch->subspaces_to_exclude) {
3027: PetscHSetIDestroy(&patch->subspaces_to_exclude);
3028: }
3030: VecDestroy(&patch->localRHS);
3031: VecDestroy(&patch->localUpdate);
3032: VecDestroy(&patch->patchRHS);
3033: VecDestroy(&patch->patchUpdate);
3034: VecDestroy(&patch->dof_weights);
3035: if (patch->patch_dof_weights) {
3036: for (i = 0; i < patch->npatch; ++i) {VecDestroy(&patch->patch_dof_weights[i]);}
3037: PetscFree(patch->patch_dof_weights);
3038: }
3039: if (patch->mat) {
3040: for (i = 0; i < patch->npatch; ++i) {MatDestroy(&patch->mat[i]);}
3041: PetscFree(patch->mat);
3042: }
3043: if (patch->matWithArtificial) {
3044: for (i = 0; i < patch->npatch; ++i) {MatDestroy(&patch->matWithArtificial[i]);}
3045: PetscFree(patch->matWithArtificial);
3046: }
3047: VecDestroy(&patch->patchRHSWithArtificial);
3048: if (patch->dofMappingWithoutToWithArtificial) {
3049: for (i = 0; i < patch->npatch; ++i) {ISDestroy(&patch->dofMappingWithoutToWithArtificial[i]);}
3050: PetscFree(patch->dofMappingWithoutToWithArtificial);
3052: }
3053: if (patch->dofMappingWithoutToWithAll) {
3054: for (i = 0; i < patch->npatch; ++i) {ISDestroy(&patch->dofMappingWithoutToWithAll[i]);}
3055: PetscFree(patch->dofMappingWithoutToWithAll);
3057: }
3058: PetscFree(patch->sub_mat_type);
3059: if (patch->userIS) {
3060: for (i = 0; i < patch->npatch; ++i) {ISDestroy(&patch->userIS[i]);}
3061: PetscFree(patch->userIS);
3062: }
3063: PetscFree(patch->precomputedTensorLocations);
3064: PetscFree(patch->precomputedIntFacetTensorLocations);
3066: patch->bs = NULL;
3067: patch->cellNodeMap = NULL;
3068: patch->nsubspaces = 0;
3069: ISDestroy(&patch->iterationSet);
3071: PetscViewerDestroy(&patch->viewerSection);
3072: return(0);
3073: }
3075: static PetscErrorCode PCDestroy_PATCH_Linear(PC pc)
3076: {
3077: PC_PATCH *patch = (PC_PATCH *) pc->data;
3078: PetscInt i;
3082: if (patch->solver) {
3083: for (i = 0; i < patch->npatch; ++i) {KSPDestroy((KSP *) &patch->solver[i]);}
3084: PetscFree(patch->solver);
3085: }
3086: return(0);
3087: }
3089: static PetscErrorCode PCDestroy_PATCH(PC pc)
3090: {
3091: PC_PATCH *patch = (PC_PATCH *) pc->data;
3095: PCReset_PATCH(pc);
3096: (*patch->destroysolver)(pc);
3097: PetscFree(pc->data);
3098: return(0);
3099: }
3101: static PetscErrorCode PCSetFromOptions_PATCH(PetscOptionItems *PetscOptionsObject, PC pc)
3102: {
3103: PC_PATCH *patch = (PC_PATCH *) pc->data;
3104: PCPatchConstructType patchConstructionType = PC_PATCH_STAR;
3105: char sub_mat_type[PETSC_MAX_PATH_LEN];
3106: char option[PETSC_MAX_PATH_LEN];
3107: const char *prefix;
3108: PetscBool flg, dimflg, codimflg;
3109: MPI_Comm comm;
3110: PetscInt *ifields, nfields, k;
3111: PetscErrorCode ierr;
3112: PCCompositeType loctype = PC_COMPOSITE_ADDITIVE;
3115: PetscObjectGetComm((PetscObject) pc, &comm);
3116: PetscObjectGetOptionsPrefix((PetscObject) pc, &prefix);
3117: PetscOptionsHead(PetscOptionsObject, "Patch solver options");
3119: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_save_operators", patch->classname);
3120: PetscOptionsBool(option, "Store all patch operators for lifetime of object?", "PCPatchSetSaveOperators", patch->save_operators, &patch->save_operators, &flg);
3122: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_precompute_element_tensors", patch->classname);
3123: PetscOptionsBool(option, "Compute each element tensor only once?", "PCPatchSetPrecomputeElementTensors", patch->precomputeElementTensors, &patch->precomputeElementTensors, &flg);
3124: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_partition_of_unity", patch->classname);
3125: PetscOptionsBool(option, "Weight contributions by dof multiplicity?", "PCPatchSetPartitionOfUnity", patch->partition_of_unity, &patch->partition_of_unity, &flg);
3127: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_local_type", patch->classname);
3128: PetscOptionsEnum(option,"Type of local solver composition (additive or multiplicative)","PCPatchSetLocalComposition",PCCompositeTypes,(PetscEnum)loctype,(PetscEnum*)&loctype,&flg);
3129: if (flg) { PCPatchSetLocalComposition(pc, loctype);}
3130: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_dense_inverse", patch->classname);
3131: PetscOptionsBool(option, "Compute inverses of patch matrices and apply directly? Ignores KSP/PC settings on patch.", "PCPatchSetDenseInverse", patch->denseinverse, &patch->denseinverse, &flg);
3132: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_dim", patch->classname);
3133: PetscOptionsInt(option, "What dimension of mesh point to construct patches by? (0 = vertices)", "PCPATCH", patch->dim, &patch->dim, &dimflg);
3134: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_codim", patch->classname);
3135: PetscOptionsInt(option, "What co-dimension of mesh point to construct patches by? (0 = cells)", "PCPATCH", patch->codim, &patch->codim, &codimflg);
3136: if (dimflg && codimflg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Can only set one of dimension or co-dimension");
3138: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_type", patch->classname);
3139: PetscOptionsEnum(option, "How should the patches be constructed?", "PCPatchSetConstructType", PCPatchConstructTypes, (PetscEnum) patchConstructionType, (PetscEnum *) &patchConstructionType, &flg);
3140: if (flg) {PCPatchSetConstructType(pc, patchConstructionType, NULL, NULL);}
3142: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_vanka_dim", patch->classname);
3143: PetscOptionsInt(option, "Topological dimension of entities for Vanka to ignore", "PCPATCH", patch->vankadim, &patch->vankadim, &flg);
3145: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_ignore_dim", patch->classname);
3146: PetscOptionsInt(option, "Topological dimension of entities for completion to ignore", "PCPATCH", patch->ignoredim, &patch->ignoredim, &flg);
3148: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_pardecomp_overlap", patch->classname);
3149: PetscOptionsInt(option, "What overlap should we use in construct type pardecomp?", "PCPATCH", patch->pardecomp_overlap, &patch->pardecomp_overlap, &flg);
3151: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_sub_mat_type", patch->classname);
3152: PetscOptionsFList(option, "Matrix type for patch solves", "PCPatchSetSubMatType", MatList, NULL, sub_mat_type, PETSC_MAX_PATH_LEN, &flg);
3153: if (flg) {PCPatchSetSubMatType(pc, sub_mat_type);}
3155: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_symmetrise_sweep", patch->classname);
3156: PetscOptionsBool(option, "Go start->end, end->start?", "PCPATCH", patch->symmetrise_sweep, &patch->symmetrise_sweep, &flg);
3158: /* If the user has set the number of subspaces, use that for the buffer size,
3159: otherwise use a large number */
3160: if (patch->nsubspaces <= 0) {
3161: nfields = 128;
3162: } else {
3163: nfields = patch->nsubspaces;
3164: }
3165: PetscMalloc1(nfields, &ifields);
3166: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exclude_subspaces", patch->classname);
3167: PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,option,ifields,&nfields,&flg);
3168: 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");
3169: if (flg) {
3170: PetscHSetIClear(patch->subspaces_to_exclude);
3171: for (k = 0; k < nfields; k++) {
3172: PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]);
3173: }
3174: }
3175: PetscFree(ifields);
3177: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_patches_view", patch->classname);
3178: PetscOptionsBool(option, "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg);
3179: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_cells_view", patch->classname);
3180: PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerCells, &patch->formatCells, &patch->viewCells);
3181: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_interior_facets_view", patch->classname);
3182: PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerIntFacets, &patch->formatIntFacets, &patch->viewIntFacets);
3183: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exterior_facets_view", patch->classname);
3184: PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerExtFacets, &patch->formatExtFacets, &patch->viewExtFacets);
3185: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_points_view", patch->classname);
3186: PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints);
3187: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_section_view", patch->classname);
3188: PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerSection, &patch->formatSection, &patch->viewSection);
3189: PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_mat_view", patch->classname);
3190: PetscOptionsGetViewer(comm,((PetscObject) pc)->options, prefix, option, &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix);
3191: PetscOptionsTail();
3192: patch->optionsSet = PETSC_TRUE;
3193: return(0);
3194: }
3196: static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc)
3197: {
3198: PC_PATCH *patch = (PC_PATCH*) pc->data;
3199: KSPConvergedReason reason;
3200: PetscInt i;
3201: PetscErrorCode ierr;
3204: if (!patch->save_operators) {
3205: /* Can't do this here because the sub KSPs don't have an operator attached yet. */
3206: return(0);
3207: }
3208: if (patch->denseinverse) {
3209: /* No solvers */
3210: return(0);
3211: }
3212: for (i = 0; i < patch->npatch; ++i) {
3213: if (!((KSP) patch->solver[i])->setfromoptionscalled) {
3214: KSPSetFromOptions((KSP) patch->solver[i]);
3215: }
3216: KSPSetUp((KSP) patch->solver[i]);
3217: KSPGetConvergedReason((KSP) patch->solver[i], &reason);
3218: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
3219: }
3220: return(0);
3221: }
3223: static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer)
3224: {
3225: PC_PATCH *patch = (PC_PATCH *) pc->data;
3226: PetscViewer sviewer;
3227: PetscBool isascii;
3228: PetscMPIInt rank;
3232: /* TODO Redo tabbing with set tbas in new style */
3233: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
3234: if (!isascii) return(0);
3235: MPI_Comm_rank(PetscObjectComm((PetscObject) pc), &rank);
3236: PetscViewerASCIIPushTab(viewer);
3237: PetscViewerASCIIPrintf(viewer, "Subspace Correction preconditioner with %d patches\n", patch->npatch);
3238: if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) {
3239: PetscViewerASCIIPrintf(viewer, "Schwarz type: multiplicative\n");
3240: } else {
3241: PetscViewerASCIIPrintf(viewer, "Schwarz type: additive\n");
3242: }
3243: if (patch->partition_of_unity) {PetscViewerASCIIPrintf(viewer, "Weighting by partition of unity\n");}
3244: else {PetscViewerASCIIPrintf(viewer, "Not weighting by partition of unity\n");}
3245: if (patch->symmetrise_sweep) {PetscViewerASCIIPrintf(viewer, "Symmetrising sweep (start->end, then end->start)\n");}
3246: else {PetscViewerASCIIPrintf(viewer, "Not symmetrising sweep\n");}
3247: if (!patch->precomputeElementTensors) {PetscViewerASCIIPrintf(viewer, "Not precomputing element tensors (overlapping cells rebuilt in every patch assembly)\n");}
3248: else {PetscViewerASCIIPrintf(viewer, "Precomputing element tensors (each cell assembled only once)\n");}
3249: if (!patch->save_operators) {PetscViewerASCIIPrintf(viewer, "Not saving patch operators (rebuilt every PCApply)\n");}
3250: else {PetscViewerASCIIPrintf(viewer, "Saving patch operators (rebuilt every PCSetUp)\n");}
3251: if (patch->patchconstructop == PCPatchConstruct_Star) {PetscViewerASCIIPrintf(viewer, "Patch construction operator: star\n");}
3252: else if (patch->patchconstructop == PCPatchConstruct_Vanka) {PetscViewerASCIIPrintf(viewer, "Patch construction operator: Vanka\n");}
3253: else if (patch->patchconstructop == PCPatchConstruct_User) {PetscViewerASCIIPrintf(viewer, "Patch construction operator: user-specified\n");}
3254: else {PetscViewerASCIIPrintf(viewer, "Patch construction operator: unknown\n");}
3256: if (patch->denseinverse) {
3257: PetscViewerASCIIPrintf(viewer, "Explicitly forming dense inverse and applying patch solver via MatMult.\n");
3258: } else {
3259: if (patch->isNonlinear) {
3260: PetscViewerASCIIPrintf(viewer, "SNES on patches (all same):\n");
3261: } else {
3262: PetscViewerASCIIPrintf(viewer, "KSP on patches (all same):\n");
3263: }
3264: if (patch->solver) {
3265: PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
3266: if (!rank) {
3267: PetscViewerASCIIPushTab(sviewer);
3268: PetscObjectView(patch->solver[0], sviewer);
3269: PetscViewerASCIIPopTab(sviewer);
3270: }
3271: PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
3272: } else {
3273: PetscViewerASCIIPushTab(viewer);
3274: PetscViewerASCIIPrintf(viewer, "Solver not yet set.\n");
3275: PetscViewerASCIIPopTab(viewer);
3276: }
3277: }
3278: PetscViewerASCIIPopTab(viewer);
3279: return(0);
3280: }
3282: /*MC
3283: PCPATCH - A PC object that encapsulates flexible definition of blocks for overlapping and non-overlapping
3284: small block additive preconditioners. Block definition is based on topology from
3285: a DM and equation numbering from a PetscSection.
3287: Options Database Keys:
3288: + -pc_patch_cells_view - Views the process local cell numbers for each patch
3289: . -pc_patch_points_view - Views the process local mesh point numbers for each patch
3290: . -pc_patch_g2l_view - Views the map between global dofs and patch local dofs for each patch
3291: . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary
3292: - -pc_patch_sub_mat_view - Views the matrix associated with each patch
3294: Level: intermediate
3296: .seealso: PCType, PCCreate(), PCSetType()
3297: M*/
3298: PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc)
3299: {
3300: PC_PATCH *patch;
3304: PetscNewLog(pc, &patch);
3306: if (patch->subspaces_to_exclude) {
3307: PetscHSetIDestroy(&patch->subspaces_to_exclude);
3308: }
3309: PetscHSetICreate(&patch->subspaces_to_exclude);
3311: patch->classname = "pc";
3312: patch->isNonlinear = PETSC_FALSE;
3314: /* Set some defaults */
3315: patch->combined = PETSC_FALSE;
3316: patch->save_operators = PETSC_TRUE;
3317: patch->local_composition_type = PC_COMPOSITE_ADDITIVE;
3318: patch->precomputeElementTensors = PETSC_FALSE;
3319: patch->partition_of_unity = PETSC_FALSE;
3320: patch->codim = -1;
3321: patch->dim = -1;
3322: patch->vankadim = -1;
3323: patch->ignoredim = -1;
3324: patch->pardecomp_overlap = 0;
3325: patch->patchconstructop = PCPatchConstruct_Star;
3326: patch->symmetrise_sweep = PETSC_FALSE;
3327: patch->npatch = 0;
3328: patch->userIS = NULL;
3329: patch->optionsSet = PETSC_FALSE;
3330: patch->iterationSet = NULL;
3331: patch->user_patches = PETSC_FALSE;
3332: PetscStrallocpy(MATDENSE, (char **) &patch->sub_mat_type);
3333: patch->viewPatches = PETSC_FALSE;
3334: patch->viewCells = PETSC_FALSE;
3335: patch->viewPoints = PETSC_FALSE;
3336: patch->viewSection = PETSC_FALSE;
3337: patch->viewMatrix = PETSC_FALSE;
3338: patch->densesolve = NULL;
3339: patch->setupsolver = PCSetUp_PATCH_Linear;
3340: patch->applysolver = PCApply_PATCH_Linear;
3341: patch->resetsolver = PCReset_PATCH_Linear;
3342: patch->destroysolver = PCDestroy_PATCH_Linear;
3343: patch->updatemultiplicative = PCUpdateMultiplicative_PATCH_Linear;
3344: patch->dofMappingWithoutToWithArtificial = NULL;
3345: patch->dofMappingWithoutToWithAll = NULL;
3347: pc->data = (void *) patch;
3348: pc->ops->apply = PCApply_PATCH;
3349: pc->ops->applytranspose = NULL; /* PCApplyTranspose_PATCH; */
3350: pc->ops->setup = PCSetUp_PATCH;
3351: pc->ops->reset = PCReset_PATCH;
3352: pc->ops->destroy = PCDestroy_PATCH;
3353: pc->ops->setfromoptions = PCSetFromOptions_PATCH;
3354: pc->ops->setuponblocks = PCSetUpOnBlocks_PATCH;
3355: pc->ops->view = PCView_PATCH;
3356: pc->ops->applyrichardson = NULL;
3358: return(0);
3359: }