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