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