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