Actual source code: plexpreallocate.c

  1: #include <petsc/private/dmpleximpl.h>
  2: #include <petsc/private/isimpl.h>
  3: #include <petscsf.h>
  4: #include <petscds.h>

  6: /* get adjacencies due to point-to-point constraints that can't be found with DMPlexGetAdjacency() */
  7: static PetscErrorCode DMPlexComputeAnchorAdjacencies(DM dm, PetscBool useCone, PetscBool useClosure, PetscSection *anchorSectionAdj, PetscInt *anchorAdj[])
  8: {
  9:   PetscInt     pStart, pEnd;
 10:   PetscSection section, sectionGlobal, adjSec, aSec;
 11:   IS           aIS;

 13:   PetscFunctionBegin;
 14:   PetscCall(DMGetLocalSection(dm, &section));
 15:   PetscCall(DMGetGlobalSection(dm, &sectionGlobal));
 16:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), &adjSec));
 17:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
 18:   PetscCall(PetscSectionSetChart(adjSec, pStart, pEnd));

 20:   PetscCall(DMPlexGetAnchors(dm, &aSec, &aIS));
 21:   if (aSec) {
 22:     const PetscInt *anchors;
 23:     PetscInt        p, q, a, aSize, *offsets, aStart, aEnd, *inverse, iSize, *adj, adjSize;
 24:     PetscInt       *tmpAdjP = NULL, *tmpAdjQ = NULL;
 25:     PetscSection    inverseSec;

 27:     /* invert the constraint-to-anchor map */
 28:     PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)aSec), &inverseSec));
 29:     PetscCall(PetscSectionSetChart(inverseSec, pStart, pEnd));
 30:     PetscCall(ISGetLocalSize(aIS, &aSize));
 31:     PetscCall(ISGetIndices(aIS, &anchors));

 33:     for (p = 0; p < aSize; p++) {
 34:       PetscInt a = anchors[p];

 36:       PetscCall(PetscSectionAddDof(inverseSec, a, 1));
 37:     }
 38:     PetscCall(PetscSectionSetUp(inverseSec));
 39:     PetscCall(PetscSectionGetStorageSize(inverseSec, &iSize));
 40:     PetscCall(PetscMalloc1(iSize, &inverse));
 41:     PetscCall(PetscCalloc1(pEnd - pStart, &offsets));
 42:     PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd));
 43:     for (p = aStart; p < aEnd; p++) {
 44:       PetscInt dof, off;

 46:       PetscCall(PetscSectionGetDof(aSec, p, &dof));
 47:       PetscCall(PetscSectionGetOffset(aSec, p, &off));

 49:       for (q = 0; q < dof; q++) {
 50:         PetscInt iOff;

 52:         a = anchors[off + q];
 53:         PetscCall(PetscSectionGetOffset(inverseSec, a, &iOff));
 54:         inverse[iOff + offsets[a - pStart]++] = p;
 55:       }
 56:     }
 57:     PetscCall(ISRestoreIndices(aIS, &anchors));
 58:     PetscCall(PetscFree(offsets));

 60:     /* construct anchorAdj and adjSec
 61:      *
 62:      * loop over anchors:
 63:      *   construct anchor adjacency
 64:      *   loop over constrained:
 65:      *     construct constrained adjacency
 66:      *     if not in anchor adjacency, add to dofs
 67:      * setup adjSec, allocate anchorAdj
 68:      * loop over anchors:
 69:      *   construct anchor adjacency
 70:      *   loop over constrained:
 71:      *     construct constrained adjacency
 72:      *     if not in anchor adjacency
 73:      *       if not already in list, put in list
 74:      *   sort, unique, reduce dof count
 75:      * optional: compactify
 76:      */
 77:     for (p = pStart; p < pEnd; p++) {
 78:       PetscInt iDof, iOff, i, r, s, numAdjP = PETSC_DETERMINE;

 80:       PetscCall(PetscSectionGetDof(inverseSec, p, &iDof));
 81:       if (!iDof) continue;
 82:       PetscCall(PetscSectionGetOffset(inverseSec, p, &iOff));
 83:       PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, PETSC_TRUE, &numAdjP, &tmpAdjP));
 84:       for (i = 0; i < iDof; i++) {
 85:         PetscInt iNew = 0, qAdj, qAdjDof, qAdjCDof, numAdjQ = PETSC_DETERMINE;

 87:         q = inverse[iOff + i];
 88:         PetscCall(DMPlexGetAdjacency_Internal(dm, q, useCone, useClosure, PETSC_TRUE, &numAdjQ, &tmpAdjQ));
 89:         for (r = 0; r < numAdjQ; r++) {
 90:           qAdj = tmpAdjQ[r];
 91:           if ((qAdj < pStart) || (qAdj >= pEnd)) continue;
 92:           for (s = 0; s < numAdjP; s++) {
 93:             if (qAdj == tmpAdjP[s]) break;
 94:           }
 95:           if (s < numAdjP) continue;
 96:           PetscCall(PetscSectionGetDof(section, qAdj, &qAdjDof));
 97:           PetscCall(PetscSectionGetConstraintDof(section, qAdj, &qAdjCDof));
 98:           iNew += qAdjDof - qAdjCDof;
 99:         }
100:         PetscCall(PetscSectionAddDof(adjSec, p, iNew));
101:       }
102:     }

104:     PetscCall(PetscSectionSetUp(adjSec));
105:     PetscCall(PetscSectionGetStorageSize(adjSec, &adjSize));
106:     PetscCall(PetscMalloc1(adjSize, &adj));

108:     for (p = pStart; p < pEnd; p++) {
109:       PetscInt iDof, iOff, i, r, s, aOff, aOffOrig, aDof, numAdjP = PETSC_DETERMINE;

111:       PetscCall(PetscSectionGetDof(inverseSec, p, &iDof));
112:       if (!iDof) continue;
113:       PetscCall(PetscSectionGetOffset(inverseSec, p, &iOff));
114:       PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, PETSC_TRUE, &numAdjP, &tmpAdjP));
115:       PetscCall(PetscSectionGetDof(adjSec, p, &aDof));
116:       PetscCall(PetscSectionGetOffset(adjSec, p, &aOff));
117:       aOffOrig = aOff;
118:       for (i = 0; i < iDof; i++) {
119:         PetscInt qAdj, qAdjDof, qAdjCDof, qAdjOff, nd, numAdjQ = PETSC_DETERMINE;

121:         q = inverse[iOff + i];
122:         PetscCall(DMPlexGetAdjacency_Internal(dm, q, useCone, useClosure, PETSC_TRUE, &numAdjQ, &tmpAdjQ));
123:         for (r = 0; r < numAdjQ; r++) {
124:           qAdj = tmpAdjQ[r];
125:           if ((qAdj < pStart) || (qAdj >= pEnd)) continue;
126:           for (s = 0; s < numAdjP; s++) {
127:             if (qAdj == tmpAdjP[s]) break;
128:           }
129:           if (s < numAdjP) continue;
130:           PetscCall(PetscSectionGetDof(section, qAdj, &qAdjDof));
131:           PetscCall(PetscSectionGetConstraintDof(section, qAdj, &qAdjCDof));
132:           PetscCall(PetscSectionGetOffset(sectionGlobal, qAdj, &qAdjOff));
133:           for (nd = 0; nd < qAdjDof - qAdjCDof; ++nd) adj[aOff++] = (qAdjOff < 0 ? -(qAdjOff + 1) : qAdjOff) + nd;
134:         }
135:       }
136:       PetscCall(PetscSortRemoveDupsInt(&aDof, PetscSafePointerPlusOffset(adj, aOffOrig)));
137:       PetscCall(PetscSectionSetDof(adjSec, p, aDof));
138:     }
139:     *anchorAdj = adj;

141:     /* clean up */
142:     PetscCall(PetscSectionDestroy(&inverseSec));
143:     PetscCall(PetscFree(inverse));
144:     PetscCall(PetscFree(tmpAdjP));
145:     PetscCall(PetscFree(tmpAdjQ));
146:   } else {
147:     *anchorAdj = NULL;
148:     PetscCall(PetscSectionSetUp(adjSec));
149:   }
150:   *anchorSectionAdj = adjSec;
151:   PetscFunctionReturn(PETSC_SUCCESS);
152: }

154: // Determine if any of the local adjacencies match a leaf and root of the pointSF.
155: // When using isoperiodic boundary conditions, it is possible for periodic (leaf) and donor (root) pairs to be on the same rank.
156: // This check is done to ensure the adjacency in these cases is only counted for one of the mesh points rather than both.
157: static inline PetscErrorCode AdjancencyContainsLeafRootPair(PetscInt num_pairs, const PetscInt leaves[], const PetscInt roots[], PetscInt numAdj, const PetscInt tmpAdj[], PetscInt *num_leaves_found, PetscInt leaves_found[])
158: {
159:   PetscInt root_index = -1, leaf_, num_roots_found = 0;

161:   PetscFunctionBeginUser;
162:   *num_leaves_found = 0;
163:   for (PetscInt q = 0; q < numAdj; q++) {
164:     const PetscInt padj = tmpAdj[q];
165:     PetscCall(PetscFindInt(padj, num_pairs, roots, &root_index));
166:     if (root_index >= 0) {
167:       leaves_found[num_roots_found] = root_index; // Initially use leaves_found to store pair indices
168:       num_roots_found++;
169:       break;
170:     }
171:   }
172:   if (num_roots_found == 0) PetscFunctionReturn(PETSC_SUCCESS);

174:   for (PetscInt i = 0; i < num_roots_found; i++) {
175:     leaf_ = leaves[root_index];
176:     for (PetscInt q = 0; q < numAdj; q++) {
177:       if (tmpAdj[q] == leaf_) {
178:         leaves_found[*num_leaves_found] = leaf_;
179:         (*num_leaves_found)++;
180:         continue;
181:       }
182:     }
183:   }

185:   PetscCall(PetscIntSortSemiOrdered(*num_leaves_found, leaves_found));
186:   PetscFunctionReturn(PETSC_SUCCESS);
187: }

189: static PetscErrorCode DMPlexCreateAdjacencySection_Static(DM dm, PetscInt bs, PetscSF sfDof, PetscBool useCone, PetscBool useClosure, PetscBool useAnchors, PetscSection *sA, PetscInt **colIdx)
190: {
191:   MPI_Comm           comm;
192:   PetscMPIInt        myrank;
193:   PetscBool          doCommLocal, doComm, debug = PETSC_FALSE;
194:   PetscSF            sf, sfAdj;
195:   PetscSection       section, sectionGlobal, leafSectionAdj, rootSectionAdj, sectionAdj, anchorSectionAdj;
196:   PetscInt           nroots, nleaves, l, p, r;
197:   const PetscInt    *leaves;
198:   const PetscSFNode *remotes;
199:   PetscInt           dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols;
200:   PetscInt          *tmpAdj = NULL, *adj, *rootAdj, *anchorAdj = NULL, *cols, *remoteOffsets, *rootsMyRankPair = NULL, *leavesMyRankPair = NULL;
201:   PetscInt           adjSize, numMyRankPair = 0;

203:   PetscFunctionBegin;
204:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
205:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_view_preallocation", &debug, NULL));
206:   PetscCallMPI(MPI_Comm_rank(comm, &myrank));
207:   PetscCall(DMGetDimension(dm, &dim));
208:   PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &sf));
209:   PetscCall(DMGetLocalSection(dm, &section));
210:   PetscCall(DMGetGlobalSection(dm, &sectionGlobal));
211:   PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
212:   doCommLocal = nroots >= 0 ? PETSC_TRUE : PETSC_FALSE;
213:   PetscCallMPI(MPIU_Allreduce(&doCommLocal, &doComm, 1, MPI_C_BOOL, MPI_LAND, comm));
214:   /* Create section for dof adjacency (dof ==> # adj dof) */
215:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
216:   PetscCall(PetscSectionGetStorageSize(section, &numDof));
217:   PetscCall(PetscSectionCreate(comm, &leafSectionAdj));
218:   PetscCall(PetscSectionSetChart(leafSectionAdj, 0, numDof));
219:   PetscCall(PetscSectionCreate(comm, &rootSectionAdj));
220:   PetscCall(PetscSectionSetChart(rootSectionAdj, 0, numDof));
221:   PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes));

223:   // Store leaf-root pairs if remote.rank is the current rank
224:   if (nleaves >= 0) PetscCall(PetscMalloc2(nleaves, &rootsMyRankPair, nleaves, &leavesMyRankPair));
225:   for (PetscInt l = 0; l < nleaves; l++) {
226:     if (remotes[l].rank == myrank) {
227:       rootsMyRankPair[numMyRankPair]  = remotes[l].index;
228:       leavesMyRankPair[numMyRankPair] = leaves[l];
229:       numMyRankPair++;
230:     }
231:   }
232:   PetscCall(PetscIntSortSemiOrdered(numMyRankPair, rootsMyRankPair));
233:   PetscCall(PetscIntSortSemiOrdered(numMyRankPair, leavesMyRankPair));
234:   if (debug) {
235:     PetscCall(PetscPrintf(comm, "Roots on the same rank:\n"));
236:     PetscCall(PetscIntView(numMyRankPair, rootsMyRankPair, NULL));
237:   }
238:   /*
239:    section        - maps points to (# dofs, local dofs)
240:    sectionGlobal  - maps points to (# dofs, global dofs)
241:    leafSectionAdj - maps unowned local dofs to # adj dofs
242:    rootSectionAdj - maps   owned local dofs to # adj dofs
243:    adj            - adj global dofs indexed by leafSectionAdj
244:    rootAdj        - adj global dofs indexed by rootSectionAdj
245:    sf    - describes shared points across procs
246:    sfDof - describes shared dofs across procs
247:    sfAdj - describes shared adjacent dofs across procs
248:    ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point.
249:   (0). If there are point-to-point constraints, add the adjacencies of constrained points to anchors in anchorAdj
250:        (This is done in DMPlexComputeAnchorAdjacencies())
251:     1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj
252:        Reduce those counts to rootSectionAdj (now redundantly counting some interface points)
253:     2. Visit owned points on interface, count adjacencies placing in rootSectionAdj
254:        Create sfAdj connecting rootSectionAdj and leafSectionAdj
255:     3. Visit unowned points on interface, write adjacencies to adj
256:        Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies)
257:     4. Visit owned points on interface, write adjacencies to rootAdj
258:        Remove redundancy in rootAdj
259:    ** The last two traversals use transitive closure
260:     5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj)
261:        Allocate memory addressed by sectionAdj (cols)
262:     6. Visit all owned points in the subdomain, insert dof adjacencies into cols
263:    ** Knowing all the column adjacencies, check ownership and sum into dnz and onz
264:   */
265:   PetscCall(DMPlexComputeAnchorAdjacencies(dm, useCone, useClosure, &anchorSectionAdj, &anchorAdj));
266:   for (l = 0; l < nleaves; ++l) {
267:     PetscInt dof, off, d, q, anDof;
268:     PetscInt p = leaves[l], numAdj = PETSC_DETERMINE;

270:     if ((p < pStart) || (p >= pEnd)) continue;
271:     PetscCall(PetscSectionGetDof(section, p, &dof));
272:     PetscCall(PetscSectionGetOffset(section, p, &off));
273:     PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj));
274:     for (q = 0; q < numAdj; ++q) {
275:       const PetscInt padj = tmpAdj[q];
276:       PetscInt       ndof, ncdof;

278:       if ((padj < pStart) || (padj >= pEnd)) continue;
279:       PetscCall(PetscSectionGetDof(section, padj, &ndof));
280:       PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof));
281:       for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(leafSectionAdj, d, ndof - ncdof));
282:     }
283:     PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof));
284:     if (anDof) {
285:       for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(leafSectionAdj, d, anDof));
286:     }
287:   }
288:   PetscCall(PetscSectionSetUp(leafSectionAdj));
289:   if (debug) {
290:     PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n"));
291:     PetscCall(PetscSectionView(leafSectionAdj, NULL));
292:   }
293:   /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */
294:   if (doComm) {
295:     PetscCall(PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM));
296:     PetscCall(PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM));
297:     PetscCall(PetscSectionInvalidateMaxDof_Internal(rootSectionAdj));
298:   }
299:   if (debug) {
300:     PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots:\n"));
301:     PetscCall(PetscSectionView(rootSectionAdj, NULL));
302:   }
303:   /* Add in local adjacency sizes for owned dofs on interface (roots) */
304:   for (p = pStart; p < pEnd; ++p) {
305:     PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof;

307:     PetscCall(PetscSectionGetDof(section, p, &dof));
308:     PetscCall(PetscSectionGetOffset(section, p, &off));
309:     if (!dof) continue;
310:     PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof));
311:     if (adof <= 0) continue;
312:     PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj));
313:     for (q = 0; q < numAdj; ++q) {
314:       const PetscInt padj = tmpAdj[q];
315:       PetscInt       ndof, ncdof;

317:       if ((padj < pStart) || (padj >= pEnd)) continue;
318:       PetscCall(PetscSectionGetDof(section, padj, &ndof));
319:       PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof));
320:       for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(rootSectionAdj, d, ndof - ncdof));
321:     }
322:     PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof));
323:     if (anDof) {
324:       for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(rootSectionAdj, d, anDof));
325:     }
326:   }
327:   PetscCall(PetscSectionSetUp(rootSectionAdj));
328:   if (debug) {
329:     PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots after local additions:\n"));
330:     PetscCall(PetscSectionView(rootSectionAdj, NULL));
331:   }
332:   /* Create adj SF based on dof SF */
333:   PetscCall(PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets));
334:   PetscCall(PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj));
335:   PetscCall(PetscFree(remoteOffsets));
336:   if (debug) {
337:     PetscCall(PetscPrintf(comm, "Adjacency SF for Preallocation:\n"));
338:     PetscCall(PetscSFView(sfAdj, NULL));
339:   }
340:   /* Create leaf adjacency */
341:   PetscCall(PetscSectionSetUp(leafSectionAdj));
342:   PetscCall(PetscSectionGetStorageSize(leafSectionAdj, &adjSize));
343:   PetscCall(PetscCalloc1(adjSize, &adj));
344:   for (l = 0; l < nleaves; ++l) {
345:     PetscInt dof, off, d, q, anDof, anOff;
346:     PetscInt p = leaves[l], numAdj = PETSC_DETERMINE;

348:     if ((p < pStart) || (p >= pEnd)) continue;
349:     PetscCall(PetscSectionGetDof(section, p, &dof));
350:     PetscCall(PetscSectionGetOffset(section, p, &off));
351:     PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj));
352:     PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof));
353:     PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff));
354:     for (d = off; d < off + dof; ++d) {
355:       PetscInt aoff, i = 0;

357:       PetscCall(PetscSectionGetOffset(leafSectionAdj, d, &aoff));
358:       for (q = 0; q < numAdj; ++q) {
359:         const PetscInt padj = tmpAdj[q];
360:         PetscInt       ndof, ncdof, ngoff, nd;

362:         if ((padj < pStart) || (padj >= pEnd)) continue;
363:         PetscCall(PetscSectionGetDof(section, padj, &ndof));
364:         PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof));
365:         PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff));
366:         for (nd = 0; nd < ndof - ncdof; ++nd) {
367:           adj[aoff + i] = (ngoff < 0 ? -(ngoff + 1) : ngoff) + nd;
368:           ++i;
369:         }
370:       }
371:       for (q = 0; q < anDof; q++) {
372:         adj[aoff + i] = anchorAdj[anOff + q];
373:         ++i;
374:       }
375:     }
376:   }
377:   /* Debugging */
378:   if (debug) {
379:     PetscCall(PetscPrintf(comm, "Leaf adjacency indices\n"));
380:     PetscCall(PetscSectionArrayView(leafSectionAdj, adj, PETSC_INT, NULL));
381:   }
382:   /* Gather adjacent indices to root */
383:   PetscCall(PetscSectionGetStorageSize(rootSectionAdj, &adjSize));
384:   PetscCall(PetscMalloc1(adjSize, &rootAdj));
385:   for (r = 0; r < adjSize; ++r) rootAdj[r] = -1;
386:   if (doComm) {
387:     const PetscInt *indegree;
388:     PetscInt       *remoteadj, radjsize = 0;

390:     PetscCall(PetscSFComputeDegreeBegin(sfAdj, &indegree));
391:     PetscCall(PetscSFComputeDegreeEnd(sfAdj, &indegree));
392:     for (p = 0; p < adjSize; ++p) radjsize += indegree[p];
393:     PetscCall(PetscMalloc1(radjsize, &remoteadj));
394:     PetscCall(PetscSFGatherBegin(sfAdj, MPIU_INT, adj, remoteadj));
395:     PetscCall(PetscSFGatherEnd(sfAdj, MPIU_INT, adj, remoteadj));
396:     for (p = 0, l = 0, r = 0; p < adjSize; ++p, l = PetscMax(p, l + indegree[p - 1])) {
397:       PetscInt s;
398:       for (s = 0; s < indegree[p]; ++s, ++r) rootAdj[l + s] = remoteadj[r];
399:     }
400:     PetscCheck(r == radjsize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %" PetscInt_FMT " != %" PetscInt_FMT, r, radjsize);
401:     PetscCheck(l == adjSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %" PetscInt_FMT " != %" PetscInt_FMT, l, adjSize);
402:     PetscCall(PetscFree(remoteadj));
403:   }
404:   PetscCall(PetscSFDestroy(&sfAdj));
405:   PetscCall(PetscFree(adj));
406:   /* Debugging */
407:   if (debug) {
408:     PetscCall(PetscPrintf(comm, "Root adjacency indices after gather\n"));
409:     PetscCall(PetscSectionArrayView(rootSectionAdj, rootAdj, PETSC_INT, NULL));
410:   }
411:   /* Add in local adjacency indices for owned dofs on interface (roots) */
412:   for (p = pStart; p < pEnd; ++p) {
413:     PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof, anOff;

415:     PetscCall(PetscSectionGetDof(section, p, &dof));
416:     PetscCall(PetscSectionGetOffset(section, p, &off));
417:     if (!dof) continue;
418:     PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof));
419:     if (adof <= 0) continue;
420:     PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj));
421:     PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof));
422:     PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff));
423:     for (d = off; d < off + dof; ++d) {
424:       PetscInt adof, aoff, i;

426:       PetscCall(PetscSectionGetDof(rootSectionAdj, d, &adof));
427:       PetscCall(PetscSectionGetOffset(rootSectionAdj, d, &aoff));
428:       i = adof - 1;
429:       for (q = 0; q < anDof; q++) {
430:         rootAdj[aoff + i] = anchorAdj[anOff + q];
431:         --i;
432:       }
433:       for (q = 0; q < numAdj; ++q) {
434:         const PetscInt padj = tmpAdj[q];
435:         PetscInt       ndof, ncdof, ngoff, nd;

437:         if ((padj < pStart) || (padj >= pEnd)) continue;
438:         PetscCall(PetscSectionGetDof(section, padj, &ndof));
439:         PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof));
440:         PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff));
441:         for (nd = 0; nd < ndof - ncdof; ++nd) {
442:           rootAdj[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd;
443:           --i;
444:         }
445:       }
446:     }
447:   }
448:   /* Debugging */
449:   if (debug) {
450:     PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots before compression:\n"));
451:     PetscCall(PetscSectionArrayView(rootSectionAdj, rootAdj, PETSC_INT, NULL));
452:   }
453:   /* Compress indices */
454:   PetscCall(PetscSectionSetUp(rootSectionAdj));
455:   for (p = pStart; p < pEnd; ++p) {
456:     PetscInt dof, cdof, off, d;
457:     PetscInt adof, aoff;

459:     PetscCall(PetscSectionGetDof(section, p, &dof));
460:     PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
461:     PetscCall(PetscSectionGetOffset(section, p, &off));
462:     if (!dof) continue;
463:     PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof));
464:     if (adof <= 0) continue;
465:     for (d = off; d < off + dof - cdof; ++d) {
466:       PetscCall(PetscSectionGetDof(rootSectionAdj, d, &adof));
467:       PetscCall(PetscSectionGetOffset(rootSectionAdj, d, &aoff));
468:       PetscCall(PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]));
469:       PetscCall(PetscSectionSetDof(rootSectionAdj, d, adof));
470:     }
471:   }
472:   /* Debugging */
473:   if (debug) {
474:     PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots after compression:\n"));
475:     PetscCall(PetscSectionArrayView(rootSectionAdj, rootAdj, PETSC_INT, NULL));
476:   }
477:   /* Build adjacency section: Maps global indices to sets of adjacent global indices */
478:   PetscCall(PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd));
479:   PetscCall(PetscSectionCreate(comm, &sectionAdj));
480:   PetscCall(PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd));

482:   PetscInt *exclude_leaves, num_exclude_leaves, max_adjacency_size = 0;
483:   PetscCall(DMPlexGetMaxAdjacencySize_Internal(dm, useAnchors, &max_adjacency_size));
484:   PetscCall(PetscMalloc1(max_adjacency_size, &exclude_leaves));

486:   for (p = pStart; p < pEnd; ++p) {
487:     PetscInt  numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof;
488:     PetscBool found  = PETSC_TRUE;

490:     PetscCall(PetscSectionGetDof(section, p, &dof));
491:     PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
492:     PetscCall(PetscSectionGetOffset(section, p, &off));
493:     PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff));
494:     for (d = 0; d < dof - cdof; ++d) {
495:       PetscInt ldof, rdof;

497:       PetscCall(PetscSectionGetDof(leafSectionAdj, off + d, &ldof));
498:       PetscCall(PetscSectionGetDof(rootSectionAdj, off + d, &rdof));
499:       if (ldof > 0) {
500:         /* We do not own this point */
501:       } else if (rdof > 0) {
502:         PetscCall(PetscSectionSetDof(sectionAdj, goff + d, rdof));
503:       } else {
504:         found = PETSC_FALSE;
505:       }
506:     }
507:     if (found) continue;
508:     PetscCall(PetscSectionGetDof(section, p, &dof));
509:     PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff));
510:     PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj));
511:     PetscCall(AdjancencyContainsLeafRootPair(numMyRankPair, leavesMyRankPair, rootsMyRankPair, numAdj, tmpAdj, &num_exclude_leaves, exclude_leaves));
512:     for (q = 0; q < numAdj; ++q) {
513:       const PetscInt padj = tmpAdj[q];
514:       PetscInt       ndof, ncdof, noff, count;

516:       if ((padj < pStart) || (padj >= pEnd)) continue;
517:       PetscCall(PetscSectionGetDof(section, padj, &ndof));
518:       PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof));
519:       PetscCall(PetscSectionGetOffset(section, padj, &noff));
520:       // If leaf-root pair are both on this rank, only count root
521:       PetscCall(PetscFindInt(padj, num_exclude_leaves, exclude_leaves, &count));
522:       if (count >= 0) continue;
523:       for (d = goff; d < goff + dof - cdof; ++d) PetscCall(PetscSectionAddDof(sectionAdj, d, ndof - ncdof));
524:     }
525:     PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof));
526:     if (anDof) {
527:       for (d = goff; d < goff + dof - cdof; ++d) PetscCall(PetscSectionAddDof(sectionAdj, d, anDof));
528:     }
529:   }
530:   PetscCall(PetscSectionSetUp(sectionAdj));
531:   if (debug) {
532:     PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation:\n"));
533:     PetscCall(PetscSectionView(sectionAdj, NULL));
534:   }
535:   /* Get adjacent indices */
536:   PetscCall(PetscSectionGetStorageSize(sectionAdj, &numCols));
537:   PetscCall(PetscMalloc1(numCols, &cols));
538:   for (p = pStart; p < pEnd; ++p) {
539:     PetscInt  numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof, anOff;
540:     PetscBool found  = PETSC_TRUE;

542:     PetscCall(PetscSectionGetDof(section, p, &dof));
543:     PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
544:     PetscCall(PetscSectionGetOffset(section, p, &off));
545:     PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff));
546:     for (d = 0; d < dof - cdof; ++d) {
547:       PetscInt ldof, rdof;

549:       PetscCall(PetscSectionGetDof(leafSectionAdj, off + d, &ldof));
550:       PetscCall(PetscSectionGetDof(rootSectionAdj, off + d, &rdof));
551:       if (ldof > 0) {
552:         /* We do not own this point */
553:       } else if (rdof > 0) {
554:         PetscInt aoff, roff;

556:         PetscCall(PetscSectionGetOffset(sectionAdj, goff + d, &aoff));
557:         PetscCall(PetscSectionGetOffset(rootSectionAdj, off + d, &roff));
558:         PetscCall(PetscArraycpy(&cols[aoff], &rootAdj[roff], rdof));
559:       } else {
560:         found = PETSC_FALSE;
561:       }
562:     }
563:     if (found) continue;
564:     PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj));
565:     PetscCall(AdjancencyContainsLeafRootPair(numMyRankPair, leavesMyRankPair, rootsMyRankPair, numAdj, tmpAdj, &num_exclude_leaves, exclude_leaves));
566:     PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof));
567:     PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff));
568:     for (d = goff; d < goff + dof - cdof; ++d) {
569:       PetscInt adof, aoff, i = 0;

571:       PetscCall(PetscSectionGetDof(sectionAdj, d, &adof));
572:       PetscCall(PetscSectionGetOffset(sectionAdj, d, &aoff));
573:       for (q = 0; q < numAdj; ++q) {
574:         const PetscInt padj = tmpAdj[q], *ncind;
575:         PetscInt       ndof, ncdof, ngoff, nd, count;

577:         /* Adjacent points may not be in the section chart */
578:         if ((padj < pStart) || (padj >= pEnd)) continue;
579:         PetscCall(PetscSectionGetDof(section, padj, &ndof));
580:         PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof));
581:         PetscCall(PetscSectionGetConstraintIndices(section, padj, &ncind));
582:         PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff));
583:         // If leaf-root pair are both on this rank, only count root
584:         PetscCall(PetscFindInt(padj, num_exclude_leaves, exclude_leaves, &count));
585:         if (count >= 0) continue;
586:         for (nd = 0; nd < ndof - ncdof; ++nd, ++i) cols[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd;
587:       }
588:       for (q = 0; q < anDof; q++, i++) cols[aoff + i] = anchorAdj[anOff + q];
589:       PetscCheck(i == adof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of entries %" PetscInt_FMT " != %" PetscInt_FMT " for dof %" PetscInt_FMT " (point %" PetscInt_FMT ")", i, adof, d, p);
590:     }
591:   }
592:   PetscCall(PetscSectionDestroy(&anchorSectionAdj));
593:   PetscCall(PetscSectionDestroy(&leafSectionAdj));
594:   PetscCall(PetscSectionDestroy(&rootSectionAdj));
595:   PetscCall(PetscFree(exclude_leaves));
596:   PetscCall(PetscFree2(rootsMyRankPair, leavesMyRankPair));
597:   PetscCall(PetscFree(anchorAdj));
598:   PetscCall(PetscFree(rootAdj));
599:   PetscCall(PetscFree(tmpAdj));
600:   /* Debugging */
601:   if (debug) {
602:     PetscCall(PetscPrintf(comm, "Column indices\n"));
603:     PetscCall(PetscSectionArrayView(sectionAdj, cols, PETSC_INT, NULL));
604:   }

606:   *sA     = sectionAdj;
607:   *colIdx = cols;
608:   PetscFunctionReturn(PETSC_SUCCESS);
609: }

611: static PetscErrorCode DMPlexUpdateAllocation_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[])
612: {
613:   PetscSection section;
614:   PetscInt     rStart, rEnd, r, pStart, pEnd, p;

616:   PetscFunctionBegin;
617:   /* This loop needs to change to a loop over points, then field dofs, which means we need to look both sections */
618:   PetscCall(PetscLayoutGetRange(rLayout, &rStart, &rEnd));
619:   PetscCheck((rStart % bs) == 0 && (rEnd % bs) == 0, PetscObjectComm((PetscObject)rLayout), PETSC_ERR_ARG_WRONG, "Invalid layout [%" PetscInt_FMT ", %" PetscInt_FMT ") for matrix, must be divisible by block size %" PetscInt_FMT, rStart, rEnd, bs);
620:   if (f >= 0 && bs == 1) {
621:     PetscCall(DMGetLocalSection(dm, &section));
622:     PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
623:     for (p = pStart; p < pEnd; ++p) {
624:       PetscInt rS, rE;

626:       PetscCall(DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE));
627:       for (r = rS; r < rE; ++r) {
628:         PetscInt numCols, cStart, c;

630:         PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols));
631:         PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart));
632:         for (c = cStart; c < cStart + numCols; ++c) {
633:           if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
634:             ++dnz[r - rStart];
635:             if (cols[c] >= r) ++dnzu[r - rStart];
636:           } else {
637:             ++onz[r - rStart];
638:             if (cols[c] >= r) ++onzu[r - rStart];
639:           }
640:         }
641:       }
642:     }
643:   } else {
644:     /* Only loop over blocks of rows */
645:     for (r = rStart / bs; r < rEnd / bs; ++r) {
646:       const PetscInt row = r * bs;
647:       PetscInt       numCols, cStart, c;

649:       PetscCall(PetscSectionGetDof(sectionAdj, row, &numCols));
650:       PetscCall(PetscSectionGetOffset(sectionAdj, row, &cStart));
651:       for (c = cStart; c < cStart + numCols; ++c) {
652:         if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
653:           ++dnz[r - rStart / bs];
654:           if (cols[c] >= row) ++dnzu[r - rStart / bs];
655:         } else {
656:           ++onz[r - rStart / bs];
657:           if (cols[c] >= row) ++onzu[r - rStart / bs];
658:         }
659:       }
660:     }
661:     for (r = 0; r < (rEnd - rStart) / bs; ++r) {
662:       dnz[r] /= bs;
663:       onz[r] /= bs;
664:       dnzu[r] /= bs;
665:       onzu[r] /= bs;
666:     }
667:   }
668:   PetscFunctionReturn(PETSC_SUCCESS);
669: }

671: static PetscErrorCode DMPlexFillMatrix_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], Mat A)
672: {
673:   PetscSection section;
674:   PetscScalar *values;
675:   PetscInt     rStart, rEnd, r, pStart, pEnd, p, len, maxRowLen = 0;

677:   PetscFunctionBegin;
678:   PetscCall(PetscLayoutGetRange(rLayout, &rStart, &rEnd));
679:   for (r = rStart; r < rEnd; ++r) {
680:     PetscCall(PetscSectionGetDof(sectionAdj, r, &len));
681:     maxRowLen = PetscMax(maxRowLen, len);
682:   }
683:   PetscCall(PetscCalloc1(maxRowLen, &values));
684:   if (f >= 0 && bs == 1) {
685:     PetscCall(DMGetLocalSection(dm, &section));
686:     PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
687:     for (p = pStart; p < pEnd; ++p) {
688:       PetscInt rS, rE;

690:       PetscCall(DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE));
691:       for (r = rS; r < rE; ++r) {
692:         PetscInt numCols, cStart;

694:         PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols));
695:         PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart));
696:         PetscCall(MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES));
697:       }
698:     }
699:   } else {
700:     for (r = rStart; r < rEnd; ++r) {
701:       PetscInt numCols, cStart;

703:       PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols));
704:       PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart));
705:       PetscCall(MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES));
706:     }
707:   }
708:   PetscCall(PetscFree(values));
709:   PetscFunctionReturn(PETSC_SUCCESS);
710: }

712: /*@
713:   DMPlexPreallocateOperator - Calculate the matrix nonzero pattern based upon the information in the `DM`,
714:   the `PetscDS` it contains, and the default `PetscSection`.

716:   Collective

718:   Input Parameters:
719: + dm         - The `DMPLEX`
720: . bs         - The matrix blocksize
721: . dnz        - An array to hold the number of nonzeros in the diagonal block
722: . onz        - An array to hold the number of nonzeros in the off-diagonal block
723: . dnzu       - An array to hold the number of nonzeros in the upper triangle of the diagonal block
724: . onzu       - An array to hold the number of nonzeros in the upper triangle of the off-diagonal block
725: - fillMatrix - If `PETSC_TRUE`, fill the matrix with zeros

727:   Output Parameter:
728: . A - The preallocated matrix

730:   Level: advanced

732: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreateMatrix()`
733: @*/
734: PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
735: {
736:   MPI_Comm     comm;
737:   MatType      mtype;
738:   PetscSF      sf, sfDof;
739:   PetscSection section;
740:   PetscInt    *remoteOffsets;
741:   PetscSection sectionAdj[4] = {NULL, NULL, NULL, NULL};
742:   PetscInt    *cols[4]       = {NULL, NULL, NULL, NULL};
743:   PetscBool    useCone, useClosure;
744:   PetscInt     Nf, f, idx, locRows;
745:   PetscLayout  rLayout;
746:   PetscBool    isSymBlock, isSymSeqBlock, isSymMPIBlock, debug = PETSC_FALSE;

748:   PetscFunctionBegin;
751:   if (dnz) PetscAssertPointer(dnz, 3);
752:   if (onz) PetscAssertPointer(onz, 4);
753:   if (dnzu) PetscAssertPointer(dnzu, 5);
754:   if (onzu) PetscAssertPointer(onzu, 6);
755:   PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &sf));
756:   PetscCall(DMGetLocalSection(dm, &section));
757:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_view_preallocation", &debug, NULL));
758:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
759:   PetscCall(PetscLogEventBegin(DMPLEX_Preallocate, dm, 0, 0, 0));
760:   /* Create dof SF based on point SF */
761:   if (debug) {
762:     PetscSection section, sectionGlobal;
763:     PetscSF      sf;

765:     PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &sf));
766:     PetscCall(DMGetLocalSection(dm, &section));
767:     PetscCall(DMGetGlobalSection(dm, &sectionGlobal));
768:     PetscCall(PetscPrintf(comm, "Input Section for Preallocation:\n"));
769:     PetscCall(PetscSectionView(section, NULL));
770:     PetscCall(PetscPrintf(comm, "Input Global Section for Preallocation:\n"));
771:     PetscCall(PetscSectionView(sectionGlobal, NULL));
772:     PetscCall(PetscPrintf(comm, "Input SF for Preallocation:\n"));
773:     PetscCall(PetscSFView(sf, NULL));
774:   }
775:   PetscCall(PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets));
776:   PetscCall(PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof));
777:   PetscCall(PetscFree(remoteOffsets));
778:   if (debug) {
779:     PetscCall(PetscPrintf(comm, "Dof SF for Preallocation:\n"));
780:     PetscCall(PetscSFView(sfDof, NULL));
781:   }
782:   /* Create allocation vectors from adjacency graph */
783:   PetscCall(MatGetLocalSize(A, &locRows, NULL));
784:   PetscCall(PetscLayoutCreate(comm, &rLayout));
785:   PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
786:   PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
787:   PetscCall(PetscLayoutSetUp(rLayout));
788:   /* There are 4 types of adjacency */
789:   PetscCall(PetscSectionGetNumFields(section, &Nf));
790:   if (Nf < 1 || bs > 1) {
791:     PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure));
792:     idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
793:     PetscCall(DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]));
794:     PetscCall(DMPlexUpdateAllocation_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu));
795:   } else {
796:     for (f = 0; f < Nf; ++f) {
797:       PetscCall(DMGetAdjacency(dm, f, &useCone, &useClosure));
798:       idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
799:       if (!sectionAdj[idx]) PetscCall(DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]));
800:       PetscCall(DMPlexUpdateAllocation_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu));
801:     }
802:   }
803:   PetscCall(PetscSFDestroy(&sfDof));
804:   /* Set matrix pattern */
805:   PetscCall(MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu));
806:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
807:   /* Check for symmetric storage */
808:   PetscCall(MatGetType(A, &mtype));
809:   PetscCall(PetscStrcmp(mtype, MATSBAIJ, &isSymBlock));
810:   PetscCall(PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock));
811:   PetscCall(PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock));
812:   if (isSymBlock || isSymSeqBlock || isSymMPIBlock) PetscCall(MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE));
813:   /* Fill matrix with zeros */
814:   if (fillMatrix) {
815:     if (Nf < 1 || bs > 1) {
816:       PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure));
817:       idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
818:       PetscCall(DMPlexFillMatrix_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], A));
819:     } else {
820:       for (f = 0; f < Nf; ++f) {
821:         PetscCall(DMGetAdjacency(dm, f, &useCone, &useClosure));
822:         idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
823:         PetscCall(DMPlexFillMatrix_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], A));
824:       }
825:     }
826:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
827:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
828:   }
829:   PetscCall(PetscLayoutDestroy(&rLayout));
830:   for (idx = 0; idx < 4; ++idx) {
831:     PetscCall(PetscSectionDestroy(&sectionAdj[idx]));
832:     PetscCall(PetscFree(cols[idx]));
833:   }
834:   PetscCall(PetscLogEventEnd(DMPLEX_Preallocate, dm, 0, 0, 0));
835:   PetscFunctionReturn(PETSC_SUCCESS);
836: }

838: #if 0
839: PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
840: {
841:   PetscInt       *tmpClosure,*tmpAdj,*visits;
842:   PetscInt        c,cStart,cEnd,pStart,pEnd;

844:   PetscFunctionBegin;
845:   PetscCall(DMGetDimension(dm, &dim));
846:   PetscCall(DMPlexGetDepth(dm, &depth));
847:   PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize));

849:   maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1));

851:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
852:   npoints = pEnd - pStart;

854:   PetscCall(PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits));
855:   PetscCall(PetscArrayzero(lvisits,pEnd-pStart));
856:   PetscCall(PetscArrayzero(visits,pEnd-pStart));
857:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
858:   for (c=cStart; c<cEnd; c++) {
859:     PetscInt *support = tmpClosure;
860:     PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support));
861:     for (p=0; p<supportSize; p++) lvisits[support[p]]++;
862:   }
863:   PetscCall(PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM));
864:   PetscCall(PetscSFReduceEnd  (sf,MPIU_INT,lvisits,visits,MPI_SUM));
865:   PetscCall(PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits,MPI_REPLACE));
866:   PetscCall(PetscSFBcastEnd  (sf,MPIU_INT,visits,lvisits));

868:   PetscCall(PetscSFGetRootRanks());

870:   PetscCall(PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner));
871:   for (c=cStart; c<cEnd; c++) {
872:     PetscCall(PetscArrayzero(cellmat,maxClosureSize*maxClosureSize));
873:     /*
874:      Depth-first walk of transitive closure.
875:      At each leaf frame f of transitive closure that we see, add 1/visits[f] to each pair (p,q) not marked as done in cellmat.
876:      This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
877:      */
878:   }

880:   PetscCall(PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM));
881:   PetscCall(PetscSFReduceEnd  (sf,MPIU_INT,lonz,onz,MPI_SUM));
882:   PetscFunctionReturn(PETSC_SUCCESS);
883: }
884: #endif