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