xref: /petsc/src/dm/impls/plex/plexpreallocate.c (revision 2fa40bb9206b96114faa7cb222621ec184d31cd2) 
1 #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 #include <petsc/private/isimpl.h>
3 #include <petscsf.h>
4 #include <petscds.h>
5 
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;
12   PetscErrorCode ierr;
13 
14   PetscFunctionBegin;
15   ierr = DMGetLocalSection(dm, &section);CHKERRQ(ierr);
16   ierr = DMGetGlobalSection(dm, &sectionGlobal);CHKERRQ(ierr);
17   ierr = PetscSectionCreate(PetscObjectComm((PetscObject) section), &adjSec);CHKERRQ(ierr);
18   ierr = PetscSectionGetChart(section,&pStart,&pEnd);CHKERRQ(ierr);
19   ierr = PetscSectionSetChart(adjSec,pStart,pEnd);CHKERRQ(ierr);
20 
21   ierr = DMPlexGetAnchors(dm,&aSec,&aIS);CHKERRQ(ierr);
22   if (aSec) {
23     const PetscInt *anchors;
24     PetscInt       p, q, a, aSize, *offsets, aStart, aEnd, *inverse, iSize, *adj, adjSize;
25     PetscInt       *tmpAdjP = NULL, *tmpAdjQ = NULL;
26     PetscSection   inverseSec;
27 
28     /* invert the constraint-to-anchor map */
29     ierr = PetscSectionCreate(PetscObjectComm((PetscObject)aSec),&inverseSec);CHKERRQ(ierr);
30     ierr = PetscSectionSetChart(inverseSec,pStart,pEnd);CHKERRQ(ierr);
31     ierr = ISGetLocalSize(aIS, &aSize);CHKERRQ(ierr);
32     ierr = ISGetIndices(aIS, &anchors);CHKERRQ(ierr);
33 
34     for (p = 0; p < aSize; p++) {
35       PetscInt a = anchors[p];
36 
37       ierr = PetscSectionAddDof(inverseSec,a,1);CHKERRQ(ierr);
38     }
39     ierr = PetscSectionSetUp(inverseSec);CHKERRQ(ierr);
40     ierr = PetscSectionGetStorageSize(inverseSec,&iSize);CHKERRQ(ierr);
41     ierr = PetscMalloc1(iSize,&inverse);CHKERRQ(ierr);
42     ierr = PetscCalloc1(pEnd-pStart,&offsets);CHKERRQ(ierr);
43     ierr = PetscSectionGetChart(aSec,&aStart,&aEnd);CHKERRQ(ierr);
44     for (p = aStart; p < aEnd; p++) {
45       PetscInt dof, off;
46 
47       ierr = PetscSectionGetDof(aSec, p, &dof);CHKERRQ(ierr);
48       ierr = PetscSectionGetOffset(aSec, p, &off);CHKERRQ(ierr);
49 
50       for (q = 0; q < dof; q++) {
51         PetscInt iOff;
52 
53         a = anchors[off + q];
54         ierr = PetscSectionGetOffset(inverseSec, a, &iOff);CHKERRQ(ierr);
55         inverse[iOff + offsets[a-pStart]++] = p;
56       }
57     }
58     ierr = ISRestoreIndices(aIS, &anchors);CHKERRQ(ierr);
59     ierr = PetscFree(offsets);CHKERRQ(ierr);
60 
61     /* construct anchorAdj and adjSec
62      *
63      * loop over anchors:
64      *   construct anchor adjacency
65      *   loop over constrained:
66      *     construct constrained adjacency
67      *     if not in anchor adjacency, add to dofs
68      * setup adjSec, allocate anchorAdj
69      * loop over anchors:
70      *   construct anchor adjacency
71      *   loop over constrained:
72      *     construct constrained adjacency
73      *     if not in anchor adjacency
74      *       if not already in list, put in list
75      *   sort, unique, reduce dof count
76      * optional: compactify
77      */
78     for (p = pStart; p < pEnd; p++) {
79       PetscInt iDof, iOff, i, r, s, numAdjP = PETSC_DETERMINE;
80 
81       ierr = PetscSectionGetDof(inverseSec,p,&iDof);CHKERRQ(ierr);
82       if (!iDof) continue;
83       ierr = PetscSectionGetOffset(inverseSec,p,&iOff);CHKERRQ(ierr);
84       ierr = DMPlexGetAdjacency_Internal(dm,p,useCone,useClosure,PETSC_TRUE,&numAdjP,&tmpAdjP);CHKERRQ(ierr);
85       for (i = 0; i < iDof; i++) {
86         PetscInt iNew = 0, qAdj, qAdjDof, qAdjCDof, numAdjQ = PETSC_DETERMINE;
87 
88         q = inverse[iOff + i];
89         ierr = DMPlexGetAdjacency_Internal(dm,q,useCone,useClosure,PETSC_TRUE,&numAdjQ,&tmpAdjQ);CHKERRQ(ierr);
90         for (r = 0; r < numAdjQ; r++) {
91           qAdj = tmpAdjQ[r];
92           if ((qAdj < pStart) || (qAdj >= pEnd)) continue;
93           for (s = 0; s < numAdjP; s++) {
94             if (qAdj == tmpAdjP[s]) break;
95           }
96           if (s < numAdjP) continue;
97           ierr  = PetscSectionGetDof(section,qAdj,&qAdjDof);CHKERRQ(ierr);
98           ierr  = PetscSectionGetConstraintDof(section,qAdj,&qAdjCDof);CHKERRQ(ierr);
99           iNew += qAdjDof - qAdjCDof;
100         }
101         ierr = PetscSectionAddDof(adjSec,p,iNew);CHKERRQ(ierr);
102       }
103     }
104 
105     ierr = PetscSectionSetUp(adjSec);CHKERRQ(ierr);
106     ierr = PetscSectionGetStorageSize(adjSec,&adjSize);CHKERRQ(ierr);
107     ierr = PetscMalloc1(adjSize,&adj);CHKERRQ(ierr);
108 
109     for (p = pStart; p < pEnd; p++) {
110       PetscInt iDof, iOff, i, r, s, aOff, aOffOrig, aDof, numAdjP = PETSC_DETERMINE;
111 
112       ierr = PetscSectionGetDof(inverseSec,p,&iDof);CHKERRQ(ierr);
113       if (!iDof) continue;
114       ierr = PetscSectionGetOffset(inverseSec,p,&iOff);CHKERRQ(ierr);
115       ierr = DMPlexGetAdjacency_Internal(dm,p,useCone,useClosure,PETSC_TRUE,&numAdjP,&tmpAdjP);CHKERRQ(ierr);
116       ierr = PetscSectionGetDof(adjSec,p,&aDof);CHKERRQ(ierr);
117       ierr = PetscSectionGetOffset(adjSec,p,&aOff);CHKERRQ(ierr);
118       aOffOrig = aOff;
119       for (i = 0; i < iDof; i++) {
120         PetscInt qAdj, qAdjDof, qAdjCDof, qAdjOff, nd, numAdjQ = PETSC_DETERMINE;
121 
122         q = inverse[iOff + i];
123         ierr = DMPlexGetAdjacency_Internal(dm,q,useCone,useClosure,PETSC_TRUE,&numAdjQ,&tmpAdjQ);CHKERRQ(ierr);
124         for (r = 0; r < numAdjQ; r++) {
125           qAdj = tmpAdjQ[r];
126           if ((qAdj < pStart) || (qAdj >= pEnd)) continue;
127           for (s = 0; s < numAdjP; s++) {
128             if (qAdj == tmpAdjP[s]) break;
129           }
130           if (s < numAdjP) continue;
131           ierr  = PetscSectionGetDof(section,qAdj,&qAdjDof);CHKERRQ(ierr);
132           ierr  = PetscSectionGetConstraintDof(section,qAdj,&qAdjCDof);CHKERRQ(ierr);
133           ierr  = PetscSectionGetOffset(sectionGlobal,qAdj,&qAdjOff);CHKERRQ(ierr);
134           for (nd = 0; nd < qAdjDof-qAdjCDof; ++nd) {
135             adj[aOff++] = (qAdjOff < 0 ? -(qAdjOff+1) : qAdjOff) + nd;
136           }
137         }
138       }
139       ierr = PetscSortRemoveDupsInt(&aDof,&adj[aOffOrig]);CHKERRQ(ierr);
140       ierr = PetscSectionSetDof(adjSec,p,aDof);CHKERRQ(ierr);
141     }
142     *anchorAdj = adj;
143 
144     /* clean up */
145     ierr = PetscSectionDestroy(&inverseSec);CHKERRQ(ierr);
146     ierr = PetscFree(inverse);CHKERRQ(ierr);
147     ierr = PetscFree(tmpAdjP);CHKERRQ(ierr);
148     ierr = PetscFree(tmpAdjQ);CHKERRQ(ierr);
149   }
150   else {
151     *anchorAdj = NULL;
152     ierr = PetscSectionSetUp(adjSec);CHKERRQ(ierr);
153   }
154   *anchorSectionAdj = adjSec;
155   PetscFunctionReturn(0);
156 }
157 
158 static PetscErrorCode DMPlexCreateAdjacencySection_Static(DM dm, PetscInt bs, PetscSF sfDof, PetscBool useCone, PetscBool useClosure, PetscBool useAnchors, PetscSection *sA, PetscInt **colIdx)
159 {
160   MPI_Comm           comm;
161   PetscMPIInt        size;
162   PetscBool          doCommLocal, doComm, debug = PETSC_FALSE;
163   PetscSF            sf, sfAdj;
164   PetscSection       section, sectionGlobal, leafSectionAdj, rootSectionAdj, sectionAdj, anchorSectionAdj;
165   PetscInt           nroots, nleaves, l, p, r;
166   const PetscInt    *leaves;
167   const PetscSFNode *remotes;
168   PetscInt           dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols;
169   PetscInt          *tmpAdj = NULL, *adj, *rootAdj, *anchorAdj = NULL, *cols, *remoteOffsets;
170   PetscInt           adjSize;
171   PetscErrorCode     ierr;
172 
173   PetscFunctionBegin;
174   ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
175   ierr = PetscOptionsGetBool(NULL,NULL, "-dm_view_preallocation", &debug, NULL);CHKERRQ(ierr);
176   ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr);
177   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
178   ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
179   ierr = DMGetLocalSection(dm, &section);CHKERRQ(ierr);
180   ierr = DMGetGlobalSection(dm, &sectionGlobal);CHKERRQ(ierr);
181   ierr = PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);CHKERRQ(ierr);
182   doCommLocal = (size > 1) && (nroots >= 0) ? PETSC_TRUE : PETSC_FALSE;
183   ierr = MPIU_Allreduce(&doCommLocal, &doComm, 1, MPIU_BOOL, MPI_LAND, comm);CHKERRMPI(ierr);
184   /* Create section for dof adjacency (dof ==> # adj dof) */
185   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
186   ierr = PetscSectionGetStorageSize(section, &numDof);CHKERRQ(ierr);
187   ierr = PetscSectionCreate(comm, &leafSectionAdj);CHKERRQ(ierr);
188   ierr = PetscSectionSetChart(leafSectionAdj, 0, numDof);CHKERRQ(ierr);
189   ierr = PetscSectionCreate(comm, &rootSectionAdj);CHKERRQ(ierr);
190   ierr = PetscSectionSetChart(rootSectionAdj, 0, numDof);CHKERRQ(ierr);
191   /*   Fill in the ghost dofs on the interface */
192   ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes);CHKERRQ(ierr);
193   /*
194    section        - maps points to (# dofs, local dofs)
195    sectionGlobal  - maps points to (# dofs, global dofs)
196    leafSectionAdj - maps unowned local dofs to # adj dofs
197    rootSectionAdj - maps   owned local dofs to # adj dofs
198    adj            - adj global dofs indexed by leafSectionAdj
199    rootAdj        - adj global dofs indexed by rootSectionAdj
200    sf    - describes shared points across procs
201    sfDof - describes shared dofs across procs
202    sfAdj - describes shared adjacent dofs across procs
203    ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point.
204   (0). If there are point-to-point constraints, add the adjacencies of constrained points to anchors in anchorAdj
205        (This is done in DMPlexComputeAnchorAdjacencies())
206     1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj
207        Reduce those counts to rootSectionAdj (now redundantly counting some interface points)
208     2. Visit owned points on interface, count adjacencies placing in rootSectionAdj
209        Create sfAdj connecting rootSectionAdj and leafSectionAdj
210     3. Visit unowned points on interface, write adjacencies to adj
211        Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies)
212     4. Visit owned points on interface, write adjacencies to rootAdj
213        Remove redundancy in rootAdj
214    ** The last two traversals use transitive closure
215     5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj)
216        Allocate memory addressed by sectionAdj (cols)
217     6. Visit all owned points in the subdomain, insert dof adjacencies into cols
218    ** Knowing all the column adjacencies, check ownership and sum into dnz and onz
219   */
220   ierr = DMPlexComputeAnchorAdjacencies(dm, useCone, useClosure, &anchorSectionAdj, &anchorAdj);CHKERRQ(ierr);
221   for (l = 0; l < nleaves; ++l) {
222     PetscInt dof, off, d, q, anDof;
223     PetscInt p = leaves[l], numAdj = PETSC_DETERMINE;
224 
225     if ((p < pStart) || (p >= pEnd)) continue;
226     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
227     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
228     ierr = DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);CHKERRQ(ierr);
229     for (q = 0; q < numAdj; ++q) {
230       const PetscInt padj = tmpAdj[q];
231       PetscInt ndof, ncdof;
232 
233       if ((padj < pStart) || (padj >= pEnd)) continue;
234       ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
235       ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
236       for (d = off; d < off+dof; ++d) {
237         ierr = PetscSectionAddDof(leafSectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
238       }
239     }
240     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
241     if (anDof) {
242       for (d = off; d < off+dof; ++d) {
243         ierr = PetscSectionAddDof(leafSectionAdj, d, anDof);CHKERRQ(ierr);
244       }
245     }
246   }
247   ierr = PetscSectionSetUp(leafSectionAdj);CHKERRQ(ierr);
248   if (debug) {
249     ierr = PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n");CHKERRQ(ierr);
250     ierr = PetscSectionView(leafSectionAdj, NULL);CHKERRQ(ierr);
251   }
252   /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */
253   if (doComm) {
254     ierr = PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);CHKERRQ(ierr);
255     ierr = PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);CHKERRQ(ierr);
256   }
257   if (debug) {
258     ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots:\n");CHKERRQ(ierr);
259     ierr = PetscSectionView(rootSectionAdj, NULL);CHKERRQ(ierr);
260   }
261   /* Add in local adjacency sizes for owned dofs on interface (roots) */
262   for (p = pStart; p < pEnd; ++p) {
263     PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof;
264 
265     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
266     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
267     if (!dof) continue;
268     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
269     if (adof <= 0) continue;
270     ierr = DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);CHKERRQ(ierr);
271     for (q = 0; q < numAdj; ++q) {
272       const PetscInt padj = tmpAdj[q];
273       PetscInt ndof, ncdof;
274 
275       if ((padj < pStart) || (padj >= pEnd)) continue;
276       ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
277       ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
278       for (d = off; d < off+dof; ++d) {
279         ierr = PetscSectionAddDof(rootSectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
280       }
281     }
282     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
283     if (anDof) {
284       for (d = off; d < off+dof; ++d) {
285         ierr = PetscSectionAddDof(rootSectionAdj, d, anDof);CHKERRQ(ierr);
286       }
287     }
288   }
289   ierr = PetscSectionSetUp(rootSectionAdj);CHKERRQ(ierr);
290   if (debug) {
291     ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after local additions:\n");CHKERRQ(ierr);
292     ierr = PetscSectionView(rootSectionAdj, NULL);CHKERRQ(ierr);
293   }
294   /* Create adj SF based on dof SF */
295   ierr = PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets);CHKERRQ(ierr);
296   ierr = PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj);CHKERRQ(ierr);
297   ierr = PetscFree(remoteOffsets);CHKERRQ(ierr);
298   if (debug && size > 1) {
299     ierr = PetscPrintf(comm, "Adjacency SF for Preallocation:\n");CHKERRQ(ierr);
300     ierr = PetscSFView(sfAdj, NULL);CHKERRQ(ierr);
301   }
302   /* Create leaf adjacency */
303   ierr = PetscSectionSetUp(leafSectionAdj);CHKERRQ(ierr);
304   ierr = PetscSectionGetStorageSize(leafSectionAdj, &adjSize);CHKERRQ(ierr);
305   ierr = PetscCalloc1(adjSize, &adj);CHKERRQ(ierr);
306   for (l = 0; l < nleaves; ++l) {
307     PetscInt dof, off, d, q, anDof, anOff;
308     PetscInt p = leaves[l], numAdj = PETSC_DETERMINE;
309 
310     if ((p < pStart) || (p >= pEnd)) continue;
311     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
312     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
313     ierr = DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);CHKERRQ(ierr);
314     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
315     ierr = PetscSectionGetOffset(anchorSectionAdj, p, &anOff);CHKERRQ(ierr);
316     for (d = off; d < off+dof; ++d) {
317       PetscInt aoff, i = 0;
318 
319       ierr = PetscSectionGetOffset(leafSectionAdj, d, &aoff);CHKERRQ(ierr);
320       for (q = 0; q < numAdj; ++q) {
321         const PetscInt padj = tmpAdj[q];
322         PetscInt ndof, ncdof, ngoff, nd;
323 
324         if ((padj < pStart) || (padj >= pEnd)) continue;
325         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
326         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
327         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
328         for (nd = 0; nd < ndof-ncdof; ++nd) {
329           adj[aoff+i] = (ngoff < 0 ? -(ngoff+1) : ngoff) + nd;
330           ++i;
331         }
332       }
333       for (q = 0; q < anDof; q++) {
334         adj[aoff+i] = anchorAdj[anOff+q];
335         ++i;
336       }
337     }
338   }
339   /* Debugging */
340   if (debug) {
341     IS tmp;
342     ierr = PetscPrintf(comm, "Leaf adjacency indices\n");CHKERRQ(ierr);
343     ierr = ISCreateGeneral(comm, adjSize, adj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
344     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
345     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
346   }
347   /* Gather adjacent indices to root */
348   ierr = PetscSectionGetStorageSize(rootSectionAdj, &adjSize);CHKERRQ(ierr);
349   ierr = PetscMalloc1(adjSize, &rootAdj);CHKERRQ(ierr);
350   for (r = 0; r < adjSize; ++r) rootAdj[r] = -1;
351   if (doComm) {
352     const PetscInt *indegree;
353     PetscInt       *remoteadj, radjsize = 0;
354 
355     ierr = PetscSFComputeDegreeBegin(sfAdj, &indegree);CHKERRQ(ierr);
356     ierr = PetscSFComputeDegreeEnd(sfAdj, &indegree);CHKERRQ(ierr);
357     for (p = 0; p < adjSize; ++p) radjsize += indegree[p];
358     ierr = PetscMalloc1(radjsize, &remoteadj);CHKERRQ(ierr);
359     ierr = PetscSFGatherBegin(sfAdj, MPIU_INT, adj, remoteadj);CHKERRQ(ierr);
360     ierr = PetscSFGatherEnd(sfAdj, MPIU_INT, adj, remoteadj);CHKERRQ(ierr);
361     for (p = 0, l = 0, r = 0; p < adjSize; ++p, l = PetscMax(p, l + indegree[p-1])) {
362       PetscInt s;
363       for (s = 0; s < indegree[p]; ++s, ++r) rootAdj[l+s] = remoteadj[r];
364     }
365     if (r != radjsize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %d != %d", r, radjsize);
366     if (l != adjSize)  SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %d != %d", l, adjSize);
367     ierr = PetscFree(remoteadj);CHKERRQ(ierr);
368   }
369   ierr = PetscSFDestroy(&sfAdj);CHKERRQ(ierr);
370   ierr = PetscFree(adj);CHKERRQ(ierr);
371   /* Debugging */
372   if (debug) {
373     IS tmp;
374     ierr = PetscPrintf(comm, "Root adjacency indices after gather\n");CHKERRQ(ierr);
375     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
376     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
377     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
378   }
379   /* Add in local adjacency indices for owned dofs on interface (roots) */
380   for (p = pStart; p < pEnd; ++p) {
381     PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof, anOff;
382 
383     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
384     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
385     if (!dof) continue;
386     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
387     if (adof <= 0) continue;
388     ierr = DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);CHKERRQ(ierr);
389     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
390     ierr = PetscSectionGetOffset(anchorSectionAdj, p, &anOff);CHKERRQ(ierr);
391     for (d = off; d < off+dof; ++d) {
392       PetscInt adof, aoff, i;
393 
394       ierr = PetscSectionGetDof(rootSectionAdj, d, &adof);CHKERRQ(ierr);
395       ierr = PetscSectionGetOffset(rootSectionAdj, d, &aoff);CHKERRQ(ierr);
396       i    = adof-1;
397       for (q = 0; q < anDof; q++) {
398         rootAdj[aoff+i] = anchorAdj[anOff+q];
399         --i;
400       }
401       for (q = 0; q < numAdj; ++q) {
402         const PetscInt padj = tmpAdj[q];
403         PetscInt ndof, ncdof, ngoff, nd;
404 
405         if ((padj < pStart) || (padj >= pEnd)) continue;
406         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
407         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
408         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
409         for (nd = 0; nd < ndof-ncdof; ++nd) {
410           rootAdj[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
411           --i;
412         }
413       }
414     }
415   }
416   /* Debugging */
417   if (debug) {
418     IS tmp;
419     ierr = PetscPrintf(comm, "Root adjacency indices\n");CHKERRQ(ierr);
420     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
421     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
422     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
423   }
424   /* Compress indices */
425   ierr = PetscSectionSetUp(rootSectionAdj);CHKERRQ(ierr);
426   for (p = pStart; p < pEnd; ++p) {
427     PetscInt dof, cdof, off, d;
428     PetscInt adof, aoff;
429 
430     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
431     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
432     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
433     if (!dof) continue;
434     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
435     if (adof <= 0) continue;
436     for (d = off; d < off+dof-cdof; ++d) {
437       ierr = PetscSectionGetDof(rootSectionAdj, d, &adof);CHKERRQ(ierr);
438       ierr = PetscSectionGetOffset(rootSectionAdj, d, &aoff);CHKERRQ(ierr);
439       ierr = PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]);CHKERRQ(ierr);
440       ierr = PetscSectionSetDof(rootSectionAdj, d, adof);CHKERRQ(ierr);
441     }
442   }
443   /* Debugging */
444   if (debug) {
445     IS tmp;
446     ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after compression:\n");CHKERRQ(ierr);
447     ierr = PetscSectionView(rootSectionAdj, NULL);CHKERRQ(ierr);
448     ierr = PetscPrintf(comm, "Root adjacency indices after compression\n");CHKERRQ(ierr);
449     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
450     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
451     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
452   }
453   /* Build adjacency section: Maps global indices to sets of adjacent global indices */
454   ierr = PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd);CHKERRQ(ierr);
455   ierr = PetscSectionCreate(comm, &sectionAdj);CHKERRQ(ierr);
456   ierr = PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd);CHKERRQ(ierr);
457   for (p = pStart; p < pEnd; ++p) {
458     PetscInt  numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof;
459     PetscBool found  = PETSC_TRUE;
460 
461     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
462     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
463     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
464     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
465     for (d = 0; d < dof-cdof; ++d) {
466       PetscInt ldof, rdof;
467 
468       ierr = PetscSectionGetDof(leafSectionAdj, off+d, &ldof);CHKERRQ(ierr);
469       ierr = PetscSectionGetDof(rootSectionAdj, off+d, &rdof);CHKERRQ(ierr);
470       if (ldof > 0) {
471         /* We do not own this point */
472       } else if (rdof > 0) {
473         ierr = PetscSectionSetDof(sectionAdj, goff+d, rdof);CHKERRQ(ierr);
474       } else {
475         found = PETSC_FALSE;
476       }
477     }
478     if (found) continue;
479     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
480     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
481     ierr = DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);CHKERRQ(ierr);
482     for (q = 0; q < numAdj; ++q) {
483       const PetscInt padj = tmpAdj[q];
484       PetscInt ndof, ncdof, noff;
485 
486       if ((padj < pStart) || (padj >= pEnd)) continue;
487       ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
488       ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
489       ierr = PetscSectionGetOffset(section, padj, &noff);CHKERRQ(ierr);
490       for (d = goff; d < goff+dof-cdof; ++d) {
491         ierr = PetscSectionAddDof(sectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
492       }
493     }
494     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
495     if (anDof) {
496       for (d = goff; d < goff+dof-cdof; ++d) {
497         ierr = PetscSectionAddDof(sectionAdj, d, anDof);CHKERRQ(ierr);
498       }
499     }
500   }
501   ierr = PetscSectionSetUp(sectionAdj);CHKERRQ(ierr);
502   if (debug) {
503     ierr = PetscPrintf(comm, "Adjacency Section for Preallocation:\n");CHKERRQ(ierr);
504     ierr = PetscSectionView(sectionAdj, NULL);CHKERRQ(ierr);
505   }
506   /* Get adjacent indices */
507   ierr = PetscSectionGetStorageSize(sectionAdj, &numCols);CHKERRQ(ierr);
508   ierr = PetscMalloc1(numCols, &cols);CHKERRQ(ierr);
509   for (p = pStart; p < pEnd; ++p) {
510     PetscInt  numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof, anOff;
511     PetscBool found  = PETSC_TRUE;
512 
513     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
514     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
515     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
516     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
517     for (d = 0; d < dof-cdof; ++d) {
518       PetscInt ldof, rdof;
519 
520       ierr = PetscSectionGetDof(leafSectionAdj, off+d, &ldof);CHKERRQ(ierr);
521       ierr = PetscSectionGetDof(rootSectionAdj, off+d, &rdof);CHKERRQ(ierr);
522       if (ldof > 0) {
523         /* We do not own this point */
524       } else if (rdof > 0) {
525         PetscInt aoff, roff;
526 
527         ierr = PetscSectionGetOffset(sectionAdj, goff+d, &aoff);CHKERRQ(ierr);
528         ierr = PetscSectionGetOffset(rootSectionAdj, off+d, &roff);CHKERRQ(ierr);
529         ierr = PetscArraycpy(&cols[aoff], &rootAdj[roff], rdof);CHKERRQ(ierr);
530       } else {
531         found = PETSC_FALSE;
532       }
533     }
534     if (found) continue;
535     ierr = DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);CHKERRQ(ierr);
536     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
537     ierr = PetscSectionGetOffset(anchorSectionAdj, p, &anOff);CHKERRQ(ierr);
538     for (d = goff; d < goff+dof-cdof; ++d) {
539       PetscInt adof, aoff, i = 0;
540 
541       ierr = PetscSectionGetDof(sectionAdj, d, &adof);CHKERRQ(ierr);
542       ierr = PetscSectionGetOffset(sectionAdj, d, &aoff);CHKERRQ(ierr);
543       for (q = 0; q < numAdj; ++q) {
544         const PetscInt  padj = tmpAdj[q];
545         PetscInt        ndof, ncdof, ngoff, nd;
546         const PetscInt *ncind;
547 
548         /* Adjacent points may not be in the section chart */
549         if ((padj < pStart) || (padj >= pEnd)) continue;
550         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
551         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
552         ierr = PetscSectionGetConstraintIndices(section, padj, &ncind);CHKERRQ(ierr);
553         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
554         for (nd = 0; nd < ndof-ncdof; ++nd, ++i) {
555           cols[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
556         }
557       }
558       for (q = 0; q < anDof; q++, i++) {
559         cols[aoff+i] = anchorAdj[anOff + q];
560       }
561       if (i != adof) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of entries %D != %D for dof %D (point %D)", i, adof, d, p);
562     }
563   }
564   ierr = PetscSectionDestroy(&anchorSectionAdj);CHKERRQ(ierr);
565   ierr = PetscSectionDestroy(&leafSectionAdj);CHKERRQ(ierr);
566   ierr = PetscSectionDestroy(&rootSectionAdj);CHKERRQ(ierr);
567   ierr = PetscFree(anchorAdj);CHKERRQ(ierr);
568   ierr = PetscFree(rootAdj);CHKERRQ(ierr);
569   ierr = PetscFree(tmpAdj);CHKERRQ(ierr);
570   /* Debugging */
571   if (debug) {
572     IS tmp;
573     ierr = PetscPrintf(comm, "Column indices\n");CHKERRQ(ierr);
574     ierr = ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
575     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
576     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
577   }
578 
579   *sA     = sectionAdj;
580   *colIdx = cols;
581   PetscFunctionReturn(0);
582 }
583 
584 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[])
585 {
586   PetscSection   section;
587   PetscInt       rStart, rEnd, r, pStart, pEnd, p;
588   PetscErrorCode ierr;
589 
590   PetscFunctionBegin;
591   /* This loop needs to change to a loop over points, then field dofs, which means we need to look both sections */
592   ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
593   if (rStart%bs || rEnd%bs) SETERRQ3(PetscObjectComm((PetscObject) rLayout), PETSC_ERR_ARG_WRONG, "Invalid layout [%d, %d) for matrix, must be divisible by block size %d", rStart, rEnd, bs);
594   if (f >= 0 && bs == 1) {
595     ierr = DMGetLocalSection(dm, &section);CHKERRQ(ierr);
596     ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
597     for (p = pStart; p < pEnd; ++p) {
598       PetscInt rS, rE;
599 
600       ierr = DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE);CHKERRQ(ierr);
601       for (r = rS; r < rE; ++r) {
602         PetscInt numCols, cStart, c;
603 
604         ierr = PetscSectionGetDof(sectionAdj, r, &numCols);CHKERRQ(ierr);
605         ierr = PetscSectionGetOffset(sectionAdj, r, &cStart);CHKERRQ(ierr);
606         for (c = cStart; c < cStart+numCols; ++c) {
607           if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
608             ++dnz[r-rStart];
609             if (cols[c] >= r) ++dnzu[r-rStart];
610           } else {
611             ++onz[r-rStart];
612             if (cols[c] >= r) ++onzu[r-rStart];
613           }
614         }
615       }
616     }
617   } else {
618     /* Only loop over blocks of rows */
619     for (r = rStart/bs; r < rEnd/bs; ++r) {
620       const PetscInt row = r*bs;
621       PetscInt       numCols, cStart, c;
622 
623       ierr = PetscSectionGetDof(sectionAdj, row, &numCols);CHKERRQ(ierr);
624       ierr = PetscSectionGetOffset(sectionAdj, row, &cStart);CHKERRQ(ierr);
625       for (c = cStart; c < cStart+numCols; ++c) {
626         if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
627           ++dnz[r-rStart/bs];
628           if (cols[c] >= row) ++dnzu[r-rStart/bs];
629         } else {
630           ++onz[r-rStart/bs];
631           if (cols[c] >= row) ++onzu[r-rStart/bs];
632         }
633       }
634     }
635     for (r = 0; r < (rEnd - rStart)/bs; ++r) {
636       dnz[r]  /= bs;
637       onz[r]  /= bs;
638       dnzu[r] /= bs;
639       onzu[r] /= bs;
640     }
641   }
642   PetscFunctionReturn(0);
643 }
644 
645 static PetscErrorCode DMPlexFillMatrix_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], Mat A)
646 {
647   PetscSection   section;
648   PetscScalar   *values;
649   PetscInt       rStart, rEnd, r, pStart, pEnd, p, len, maxRowLen = 0;
650   PetscErrorCode ierr;
651 
652   PetscFunctionBegin;
653   ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
654   for (r = rStart; r < rEnd; ++r) {
655     ierr      = PetscSectionGetDof(sectionAdj, r, &len);CHKERRQ(ierr);
656     maxRowLen = PetscMax(maxRowLen, len);
657   }
658   ierr = PetscCalloc1(maxRowLen, &values);CHKERRQ(ierr);
659   if (f >=0 && bs == 1) {
660     ierr = DMGetLocalSection(dm, &section);CHKERRQ(ierr);
661     ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
662     for (p = pStart; p < pEnd; ++p) {
663       PetscInt rS, rE;
664 
665       ierr = DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE);CHKERRQ(ierr);
666       for (r = rS; r < rE; ++r) {
667         PetscInt numCols, cStart;
668 
669         ierr = PetscSectionGetDof(sectionAdj, r, &numCols);CHKERRQ(ierr);
670         ierr = PetscSectionGetOffset(sectionAdj, r, &cStart);CHKERRQ(ierr);
671         ierr = MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);CHKERRQ(ierr);
672       }
673     }
674   } else {
675     for (r = rStart; r < rEnd; ++r) {
676       PetscInt numCols, cStart;
677 
678       ierr = PetscSectionGetDof(sectionAdj, r, &numCols);CHKERRQ(ierr);
679       ierr = PetscSectionGetOffset(sectionAdj, r, &cStart);CHKERRQ(ierr);
680       ierr = MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);CHKERRQ(ierr);
681     }
682   }
683   ierr = PetscFree(values);CHKERRQ(ierr);
684   PetscFunctionReturn(0);
685 }
686 
687 /*@C
688   DMPlexPreallocateOperator - Calculate the matrix nonzero pattern based upon the information in the DM,
689   the PetscDS it contains, and the default PetscSection.
690 
691   Collective
692 
693   Input Parameters:
694 + dm   - The DMPlex
695 . bs   - The matrix blocksize
696 . dnz  - An array to hold the number of nonzeros in the diagonal block
697 . onz  - An array to hold the number of nonzeros in the off-diagonal block
698 . dnzu - An array to hold the number of nonzeros in the upper triangle of the diagonal block
699 . onzu - An array to hold the number of nonzeros in the upper triangle of the off-diagonal block
700 - fillMatrix - If PETSC_TRUE, fill the matrix with zeros
701 
702   Output Parameter:
703 . A - The preallocated matrix
704 
705   Level: advanced
706 
707 .seealso: DMCreateMatrix()
708 @*/
709 PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
710 {
711   MPI_Comm       comm;
712   PetscDS        prob;
713   MatType        mtype;
714   PetscSF        sf, sfDof;
715   PetscSection   section;
716   PetscInt      *remoteOffsets;
717   PetscSection   sectionAdj[4] = {NULL, NULL, NULL, NULL};
718   PetscInt      *cols[4]       = {NULL, NULL, NULL, NULL};
719   PetscBool      useCone, useClosure;
720   PetscInt       Nf, f, idx, locRows;
721   PetscLayout    rLayout;
722   PetscBool      isSymBlock, isSymSeqBlock, isSymMPIBlock, debug = PETSC_FALSE;
723   PetscMPIInt    size;
724   PetscErrorCode ierr;
725 
726   PetscFunctionBegin;
727   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
728   PetscValidHeaderSpecific(A, MAT_CLASSID, 7);
729   if (dnz) PetscValidPointer(dnz,3);
730   if (onz) PetscValidPointer(onz,4);
731   if (dnzu) PetscValidPointer(dnzu,5);
732   if (onzu) PetscValidPointer(onzu,6);
733   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
734   ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
735   ierr = DMGetLocalSection(dm, &section);CHKERRQ(ierr);
736   ierr = PetscOptionsGetBool(NULL,NULL, "-dm_view_preallocation", &debug, NULL);CHKERRQ(ierr);
737   ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
738   ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr);
739   ierr = PetscLogEventBegin(DMPLEX_Preallocate,dm,0,0,0);CHKERRQ(ierr);
740   /* Create dof SF based on point SF */
741   if (debug) {
742     PetscSection section, sectionGlobal;
743     PetscSF      sf;
744 
745     ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
746     ierr = DMGetLocalSection(dm, &section);CHKERRQ(ierr);
747     ierr = DMGetGlobalSection(dm, &sectionGlobal);CHKERRQ(ierr);
748     ierr = PetscPrintf(comm, "Input Section for Preallocation:\n");CHKERRQ(ierr);
749     ierr = PetscSectionView(section, NULL);CHKERRQ(ierr);
750     ierr = PetscPrintf(comm, "Input Global Section for Preallocation:\n");CHKERRQ(ierr);
751     ierr = PetscSectionView(sectionGlobal, NULL);CHKERRQ(ierr);
752     if (size > 1) {
753       ierr = PetscPrintf(comm, "Input SF for Preallocation:\n");CHKERRQ(ierr);
754       ierr = PetscSFView(sf, NULL);CHKERRQ(ierr);
755     }
756   }
757   ierr = PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets);CHKERRQ(ierr);
758   ierr = PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof);CHKERRQ(ierr);
759   ierr = PetscFree(remoteOffsets);CHKERRQ(ierr);
760   if (debug && size > 1) {
761     ierr = PetscPrintf(comm, "Dof SF for Preallocation:\n");CHKERRQ(ierr);
762     ierr = PetscSFView(sfDof, NULL);CHKERRQ(ierr);
763   }
764   /* Create allocation vectors from adjacency graph */
765   ierr = MatGetLocalSize(A, &locRows, NULL);CHKERRQ(ierr);
766   ierr = PetscLayoutCreate(comm, &rLayout);CHKERRQ(ierr);
767   ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
768   ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
769   ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
770   /* There are 4 types of adjacency */
771   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
772   if (Nf < 1 || bs > 1) {
773     ierr = DMGetBasicAdjacency(dm, &useCone, &useClosure);CHKERRQ(ierr);
774     idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
775     ierr = DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]);CHKERRQ(ierr);
776     ierr = DMPlexUpdateAllocation_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);CHKERRQ(ierr);
777   } else {
778     for (f = 0; f < Nf; ++f) {
779       ierr = DMGetAdjacency(dm, f, &useCone, &useClosure);CHKERRQ(ierr);
780       idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
781       if (!sectionAdj[idx]) {ierr = DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]);CHKERRQ(ierr);}
782       ierr = DMPlexUpdateAllocation_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);CHKERRQ(ierr);
783     }
784   }
785   ierr = PetscSFDestroy(&sfDof);CHKERRQ(ierr);
786   /* Set matrix pattern */
787   ierr = MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);CHKERRQ(ierr);
788   ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
789   /* Check for symmetric storage */
790   ierr = MatGetType(A, &mtype);CHKERRQ(ierr);
791   ierr = PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);CHKERRQ(ierr);
792   ierr = PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);CHKERRQ(ierr);
793   ierr = PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);CHKERRQ(ierr);
794   if (isSymBlock || isSymSeqBlock || isSymMPIBlock) {ierr = MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);CHKERRQ(ierr);}
795   /* Fill matrix with zeros */
796   if (fillMatrix) {
797     if (Nf < 1 || bs > 1) {
798       ierr = DMGetBasicAdjacency(dm, &useCone, &useClosure);CHKERRQ(ierr);
799       idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
800       ierr = DMPlexFillMatrix_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], A);CHKERRQ(ierr);
801     } else {
802       for (f = 0; f < Nf; ++f) {
803         ierr = DMGetAdjacency(dm, f, &useCone, &useClosure);CHKERRQ(ierr);
804         idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
805         ierr = DMPlexFillMatrix_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], A);CHKERRQ(ierr);
806       }
807     }
808     ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
809     ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
810   }
811   ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
812   for (idx = 0; idx < 4; ++idx) {ierr = PetscSectionDestroy(&sectionAdj[idx]);CHKERRQ(ierr); ierr = PetscFree(cols[idx]);CHKERRQ(ierr);}
813   ierr = PetscLogEventEnd(DMPLEX_Preallocate,dm,0,0,0);CHKERRQ(ierr);
814   PetscFunctionReturn(0);
815 }
816 
817 #if 0
818 PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
819 {
820   PetscInt       *tmpClosure,*tmpAdj,*visits;
821   PetscInt        c,cStart,cEnd,pStart,pEnd;
822   PetscErrorCode  ierr;
823 
824   PetscFunctionBegin;
825   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
826   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
827   ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
828 
829   maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1));
830 
831   ierr    = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
832   npoints = pEnd - pStart;
833 
834   ierr = PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits);CHKERRQ(ierr);
835   ierr = PetscArrayzero(lvisits,pEnd-pStart);CHKERRQ(ierr);
836   ierr = PetscArrayzero(visits,pEnd-pStart);CHKERRQ(ierr);
837   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
838   for (c=cStart; c<cEnd; c++) {
839     PetscInt *support = tmpClosure;
840     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);CHKERRQ(ierr);
841     for (p=0; p<supportSize; p++) lvisits[support[p]]++;
842   }
843   ierr = PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
844   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
845   ierr = PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits,MPI_REPLACE);CHKERRQ(ierr);
846   ierr = PetscSFBcastEnd  (sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
847 
848   ierr = PetscSFGetRootRanks();CHKERRQ(ierr);
849 
850   ierr = PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner);CHKERRQ(ierr);
851   for (c=cStart; c<cEnd; c++) {
852     ierr = PetscArrayzero(cellmat,maxClosureSize*maxClosureSize);CHKERRQ(ierr);
853     /*
854      Depth-first walk of transitive closure.
855      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.
856      This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
857      */
858   }
859 
860   ierr = PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);CHKERRQ(ierr);
861   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lonz,onz,MPI_SUM);CHKERRQ(ierr);
862   PetscFunctionReturn(0);
863 }
864 #endif
865