xref: /petsc/src/dm/impls/plex/plexpreallocate.c (revision 2db13446558346134dc3b699c59fe1d864b3b5ed) ! 
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, "-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 = MPI_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   if (debug) {
302     ierr = PetscPrintf(comm, "Adjacency SF for Preallocation:\n");CHKERRQ(ierr);
303     ierr = PetscSFView(sfAdj, NULL);CHKERRQ(ierr);
304   }
305   /* Create leaf adjacency */
306   ierr = PetscSectionSetUp(leafSectionAdj);CHKERRQ(ierr);
307   ierr = PetscSectionGetStorageSize(leafSectionAdj, &adjSize);CHKERRQ(ierr);
308   ierr = PetscCalloc1(adjSize, &adj);CHKERRQ(ierr);
309   for (l = 0; l < nleaves; ++l) {
310     PetscInt dof, off, d, q, anDof, anOff;
311     PetscInt p = leaves[l], numAdj = PETSC_DETERMINE;
312 
313     if ((p < pStart) || (p >= pEnd)) continue;
314     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
315     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
316     ierr = DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);CHKERRQ(ierr);
317     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
318     ierr = PetscSectionGetOffset(anchorSectionAdj, p, &anOff);CHKERRQ(ierr);
319     for (d = off; d < off+dof; ++d) {
320       PetscInt aoff, i = 0;
321 
322       ierr = PetscSectionGetOffset(leafSectionAdj, d, &aoff);CHKERRQ(ierr);
323       for (q = 0; q < numAdj; ++q) {
324         const PetscInt padj = tmpAdj[q];
325         PetscInt ndof, ncdof, ngoff, nd;
326 
327         if ((padj < pStart) || (padj >= pEnd)) continue;
328         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
329         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
330         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
331         for (nd = 0; nd < ndof-ncdof; ++nd) {
332           adj[aoff+i] = (ngoff < 0 ? -(ngoff+1) : ngoff) + nd;
333           ++i;
334         }
335       }
336       for (q = 0; q < anDof; q++) {
337         adj[aoff+i] = anchorAdj[anOff+q];
338         ++i;
339       }
340     }
341   }
342   /* Debugging */
343   if (debug) {
344     IS tmp;
345     ierr = PetscPrintf(comm, "Leaf adjacency indices\n");CHKERRQ(ierr);
346     ierr = ISCreateGeneral(comm, adjSize, adj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
347     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
348     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
349   }
350   /* Gather adjacent indices to root */
351   ierr = PetscSectionGetStorageSize(rootSectionAdj, &adjSize);CHKERRQ(ierr);
352   ierr = PetscMalloc1(adjSize, &rootAdj);CHKERRQ(ierr);
353   for (r = 0; r < adjSize; ++r) rootAdj[r] = -1;
354   if (doComm) {
355     const PetscInt *indegree;
356     PetscInt       *remoteadj, radjsize = 0;
357 
358     ierr = PetscSFComputeDegreeBegin(sfAdj, &indegree);CHKERRQ(ierr);
359     ierr = PetscSFComputeDegreeEnd(sfAdj, &indegree);CHKERRQ(ierr);
360     for (p = 0; p < adjSize; ++p) radjsize += indegree[p];
361     ierr = PetscMalloc1(radjsize, &remoteadj);CHKERRQ(ierr);
362     ierr = PetscSFGatherBegin(sfAdj, MPIU_INT, adj, remoteadj);CHKERRQ(ierr);
363     ierr = PetscSFGatherEnd(sfAdj, MPIU_INT, adj, remoteadj);CHKERRQ(ierr);
364     for (p = 0, l = 0, r = 0; p < adjSize; ++p, l = PetscMax(p, l + indegree[p-1])) {
365       PetscInt s;
366       for (s = 0; s < indegree[p]; ++s, ++r) rootAdj[l+s] = remoteadj[r];
367     }
368     if (r != radjsize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %d != %d", r, radjsize);
369     if (l != adjSize)  SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %d != %d", l, adjSize);
370     ierr = PetscFree(remoteadj);CHKERRQ(ierr);
371   }
372   ierr = PetscSFDestroy(&sfAdj);CHKERRQ(ierr);
373   ierr = PetscFree(adj);CHKERRQ(ierr);
374   /* Debugging */
375   if (debug) {
376     IS tmp;
377     ierr = PetscPrintf(comm, "Root adjacency indices after gather\n");CHKERRQ(ierr);
378     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
379     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
380     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
381   }
382   /* Add in local adjacency indices for owned dofs on interface (roots) */
383   for (p = pStart; p < pEnd; ++p) {
384     PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof, anOff;
385 
386     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
387     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
388     if (!dof) continue;
389     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
390     if (adof <= 0) continue;
391     ierr = DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);CHKERRQ(ierr);
392     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
393     ierr = PetscSectionGetOffset(anchorSectionAdj, p, &anOff);CHKERRQ(ierr);
394     for (d = off; d < off+dof; ++d) {
395       PetscInt adof, aoff, i;
396 
397       ierr = PetscSectionGetDof(rootSectionAdj, d, &adof);CHKERRQ(ierr);
398       ierr = PetscSectionGetOffset(rootSectionAdj, d, &aoff);CHKERRQ(ierr);
399       i    = adof-1;
400       for (q = 0; q < anDof; q++) {
401         rootAdj[aoff+i] = anchorAdj[anOff+q];
402         --i;
403       }
404       for (q = 0; q < numAdj; ++q) {
405         const PetscInt padj = tmpAdj[q];
406         PetscInt ndof, ncdof, ngoff, nd;
407 
408         if ((padj < pStart) || (padj >= pEnd)) continue;
409         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
410         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
411         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
412         for (nd = 0; nd < ndof-ncdof; ++nd) {
413           rootAdj[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
414           --i;
415         }
416       }
417     }
418   }
419   /* Debugging */
420   if (debug) {
421     IS tmp;
422     ierr = PetscPrintf(comm, "Root adjacency indices\n");CHKERRQ(ierr);
423     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
424     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
425     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
426   }
427   /* Compress indices */
428   ierr = PetscSectionSetUp(rootSectionAdj);CHKERRQ(ierr);
429   for (p = pStart; p < pEnd; ++p) {
430     PetscInt dof, cdof, off, d;
431     PetscInt adof, aoff;
432 
433     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
434     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
435     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
436     if (!dof) continue;
437     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
438     if (adof <= 0) continue;
439     for (d = off; d < off+dof-cdof; ++d) {
440       ierr = PetscSectionGetDof(rootSectionAdj, d, &adof);CHKERRQ(ierr);
441       ierr = PetscSectionGetOffset(rootSectionAdj, d, &aoff);CHKERRQ(ierr);
442       ierr = PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]);CHKERRQ(ierr);
443       ierr = PetscSectionSetDof(rootSectionAdj, d, adof);CHKERRQ(ierr);
444     }
445   }
446   /* Debugging */
447   if (debug) {
448     IS tmp;
449     ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after compression:\n");CHKERRQ(ierr);
450     ierr = PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
451     ierr = PetscPrintf(comm, "Root adjacency indices after compression\n");CHKERRQ(ierr);
452     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
453     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
454     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
455   }
456   /* Build adjacency section: Maps global indices to sets of adjacent global indices */
457   ierr = PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd);CHKERRQ(ierr);
458   ierr = PetscSectionCreate(comm, &sectionAdj);CHKERRQ(ierr);
459   ierr = PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd);CHKERRQ(ierr);
460   for (p = pStart; p < pEnd; ++p) {
461     PetscInt  numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof;
462     PetscBool found  = PETSC_TRUE;
463 
464     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
465     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
466     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
467     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
468     for (d = 0; d < dof-cdof; ++d) {
469       PetscInt ldof, rdof;
470 
471       ierr = PetscSectionGetDof(leafSectionAdj, off+d, &ldof);CHKERRQ(ierr);
472       ierr = PetscSectionGetDof(rootSectionAdj, off+d, &rdof);CHKERRQ(ierr);
473       if (ldof > 0) {
474         /* We do not own this point */
475       } else if (rdof > 0) {
476         ierr = PetscSectionSetDof(sectionAdj, goff+d, rdof);CHKERRQ(ierr);
477       } else {
478         found = PETSC_FALSE;
479       }
480     }
481     if (found) continue;
482     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
483     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
484     ierr = DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);CHKERRQ(ierr);
485     for (q = 0; q < numAdj; ++q) {
486       const PetscInt padj = tmpAdj[q];
487       PetscInt ndof, ncdof, noff;
488 
489       if ((padj < pStart) || (padj >= pEnd)) continue;
490       ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
491       ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
492       ierr = PetscSectionGetOffset(section, padj, &noff);CHKERRQ(ierr);
493       for (d = goff; d < goff+dof-cdof; ++d) {
494         ierr = PetscSectionAddDof(sectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
495       }
496     }
497     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
498     if (anDof) {
499       for (d = goff; d < goff+dof-cdof; ++d) {
500         ierr = PetscSectionAddDof(sectionAdj, d, anDof);CHKERRQ(ierr);
501       }
502     }
503   }
504   ierr = PetscSectionSetUp(sectionAdj);CHKERRQ(ierr);
505   if (debug) {
506     ierr = PetscPrintf(comm, "Adjacency Section for Preallocation:\n");CHKERRQ(ierr);
507     ierr = PetscSectionView(sectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
508   }
509   /* Get adjacent indices */
510   ierr = PetscSectionGetStorageSize(sectionAdj, &numCols);CHKERRQ(ierr);
511   ierr = PetscMalloc1(numCols, &cols);CHKERRQ(ierr);
512   for (p = pStart; p < pEnd; ++p) {
513     PetscInt  numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof, anOff;
514     PetscBool found  = PETSC_TRUE;
515 
516     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
517     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
518     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
519     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
520     for (d = 0; d < dof-cdof; ++d) {
521       PetscInt ldof, rdof;
522 
523       ierr = PetscSectionGetDof(leafSectionAdj, off+d, &ldof);CHKERRQ(ierr);
524       ierr = PetscSectionGetDof(rootSectionAdj, off+d, &rdof);CHKERRQ(ierr);
525       if (ldof > 0) {
526         /* We do not own this point */
527       } else if (rdof > 0) {
528         PetscInt aoff, roff;
529 
530         ierr = PetscSectionGetOffset(sectionAdj, goff+d, &aoff);CHKERRQ(ierr);
531         ierr = PetscSectionGetOffset(rootSectionAdj, off+d, &roff);CHKERRQ(ierr);
532         ierr = PetscMemcpy(&cols[aoff], &rootAdj[roff], rdof * sizeof(PetscInt));CHKERRQ(ierr);
533       } else {
534         found = PETSC_FALSE;
535       }
536     }
537     if (found) continue;
538     ierr = DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);CHKERRQ(ierr);
539     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
540     ierr = PetscSectionGetOffset(anchorSectionAdj, p, &anOff);CHKERRQ(ierr);
541     for (d = goff; d < goff+dof-cdof; ++d) {
542       PetscInt adof, aoff, i = 0;
543 
544       ierr = PetscSectionGetDof(sectionAdj, d, &adof);CHKERRQ(ierr);
545       ierr = PetscSectionGetOffset(sectionAdj, d, &aoff);CHKERRQ(ierr);
546       for (q = 0; q < numAdj; ++q) {
547         const PetscInt  padj = tmpAdj[q];
548         PetscInt        ndof, ncdof, ngoff, nd;
549         const PetscInt *ncind;
550 
551         /* Adjacent points may not be in the section chart */
552         if ((padj < pStart) || (padj >= pEnd)) continue;
553         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
554         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
555         ierr = PetscSectionGetConstraintIndices(section, padj, &ncind);CHKERRQ(ierr);
556         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
557         for (nd = 0; nd < ndof-ncdof; ++nd, ++i) {
558           cols[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
559         }
560       }
561       for (q = 0; q < anDof; q++, i++) {
562         cols[aoff+i] = anchorAdj[anOff + q];
563       }
564       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);
565     }
566   }
567   ierr = PetscSectionDestroy(&anchorSectionAdj);CHKERRQ(ierr);
568   ierr = PetscSectionDestroy(&leafSectionAdj);CHKERRQ(ierr);
569   ierr = PetscSectionDestroy(&rootSectionAdj);CHKERRQ(ierr);
570   ierr = PetscFree(anchorAdj);CHKERRQ(ierr);
571   ierr = PetscFree(rootAdj);CHKERRQ(ierr);
572   ierr = PetscFree(tmpAdj);CHKERRQ(ierr);
573   /* Debugging */
574   if (debug) {
575     IS tmp;
576     ierr = PetscPrintf(comm, "Column indices\n");CHKERRQ(ierr);
577     ierr = ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
578     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
579     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
580   }
581 
582   *sA     = sectionAdj;
583   *colIdx = cols;
584   PetscFunctionReturn(0);
585 }
586 
587 #undef __FUNCT__
588 #define __FUNCT__ "DMPlexUpdateAllocation_Static"
589 PetscErrorCode DMPlexUpdateAllocation_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[])
590 {
591   PetscSection   section;
592   PetscInt       rStart, rEnd, r, pStart, pEnd, p;
593   PetscErrorCode ierr;
594 
595   PetscFunctionBegin;
596   /* This loop needs to change to a loop over points, then field dofs, which means we need to look both sections */
597   ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
598   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);
599   if (f >= 0 && bs == 1) {
600     ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
601     ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
602     for (p = pStart; p < pEnd; ++p) {
603       PetscInt goff, loff, lfoff, fdof, fcdof, rS, rE;
604 
605       ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
606       ierr = PetscSectionGetOffset(section, p, &loff);CHKERRQ(ierr);
607       ierr = PetscSectionGetFieldOffset(section, p, f, &lfoff);CHKERRQ(ierr);
608       ierr = PetscSectionGetFieldDof(section, p, f, &fdof);CHKERRQ(ierr);
609       ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr);
610       rS   = goff + (lfoff - loff);
611       rE   = rS + fdof - fcdof;
612       for (r = rS; r < rE; ++r) {
613         PetscInt numCols, cStart, c;
614 
615         ierr = PetscSectionGetDof(sectionAdj, r, &numCols);CHKERRQ(ierr);
616         ierr = PetscSectionGetOffset(sectionAdj, r, &cStart);CHKERRQ(ierr);
617         for (c = cStart; c < cStart+numCols; ++c) {
618           if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
619             ++dnz[r-rStart];
620             if (cols[c] >= r) ++dnzu[r-rStart];
621           } else {
622             ++onz[r-rStart];
623             if (cols[c] >= r) ++onzu[r-rStart];
624           }
625         }
626       }
627     }
628   } else {
629     /* Only loop over blocks of rows */
630     for (r = rStart/bs; r < rEnd/bs; ++r) {
631       const PetscInt row = r*bs;
632       PetscInt       numCols, cStart, c;
633 
634       ierr = PetscSectionGetDof(sectionAdj, row, &numCols);CHKERRQ(ierr);
635       ierr = PetscSectionGetOffset(sectionAdj, row, &cStart);CHKERRQ(ierr);
636       for (c = cStart; c < cStart+numCols; ++c) {
637         if ((cols[c] >= rStart*bs) && (cols[c] < rEnd*bs)) {
638           ++dnz[r-rStart];
639           if (cols[c] >= row) ++dnzu[r-rStart];
640         } else {
641           ++onz[r-rStart];
642           if (cols[c] >= row) ++onzu[r-rStart];
643         }
644       }
645     }
646     for (r = 0; r < (rEnd - rStart)/bs; ++r) {
647       dnz[r]  /= bs;
648       onz[r]  /= bs;
649       dnzu[r] /= bs;
650       onzu[r] /= bs;
651     }
652   }
653   PetscFunctionReturn(0);
654 }
655 
656 #undef __FUNCT__
657 #define __FUNCT__ "DMPlexFillMatrix_Static"
658 PetscErrorCode DMPlexFillMatrix_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], Mat A)
659 {
660   PetscSection   section;
661   PetscScalar   *values;
662   PetscInt       rStart, rEnd, r, pStart, pEnd, p, len, maxRowLen = 0;
663   PetscErrorCode ierr;
664 
665   PetscFunctionBegin;
666   ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
667   for (r = rStart; r < rEnd; ++r) {
668     ierr      = PetscSectionGetDof(sectionAdj, r, &len);CHKERRQ(ierr);
669     maxRowLen = PetscMax(maxRowLen, len);
670   }
671   ierr = PetscCalloc1(maxRowLen, &values);CHKERRQ(ierr);
672   if (f >=0 && bs == 1) {
673     ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
674     ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
675     for (p = pStart; p < pEnd; ++p) {
676       PetscInt goff, loff, lfoff, fdof, fcdof, rS, rE;
677 
678       ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
679       ierr = PetscSectionGetOffset(section, p, &loff);CHKERRQ(ierr);
680       ierr = PetscSectionGetFieldOffset(section, p, f, &lfoff);CHKERRQ(ierr);
681       ierr = PetscSectionGetFieldDof(section, p, f, &fdof);CHKERRQ(ierr);
682       ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr);
683       rS   = goff + (lfoff - loff);
684       rE   = rS + fdof - fcdof;
685       for (r = rS; r < rE; ++r) {
686         PetscInt numCols, cStart;
687 
688         ierr = PetscSectionGetDof(sectionAdj, r, &numCols);CHKERRQ(ierr);
689         ierr = PetscSectionGetOffset(sectionAdj, r, &cStart);CHKERRQ(ierr);
690         ierr = MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);CHKERRQ(ierr);
691       }
692     }
693   } else {
694     for (r = rStart; r < rEnd; ++r) {
695       PetscInt numCols, cStart;
696 
697       ierr = PetscSectionGetDof(sectionAdj, r, &numCols);CHKERRQ(ierr);
698       ierr = PetscSectionGetOffset(sectionAdj, r, &cStart);CHKERRQ(ierr);
699       ierr = MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);CHKERRQ(ierr);
700     }
701   }
702   ierr = PetscFree(values);CHKERRQ(ierr);
703   PetscFunctionReturn(0);
704 }
705 
706 #undef __FUNCT__
707 #define __FUNCT__ "DMPlexPreallocateOperator"
708 PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
709 {
710   MPI_Comm       comm;
711   PetscDS        prob;
712   MatType        mtype;
713   PetscSF        sf, sfDof;
714   PetscSection   section;
715   PetscInt      *remoteOffsets;
716   PetscSection   sectionAdj[4] = {NULL, NULL, NULL, NULL};
717   PetscInt      *cols[4]       = {NULL, NULL, NULL, NULL};
718   PetscBool      useCone, useClosure;
719   PetscInt       Nf, f, idx, locRows;
720   PetscLayout    rLayout;
721   PetscBool      isSymBlock, isSymSeqBlock, isSymMPIBlock, debug = PETSC_FALSE;
722   PetscErrorCode ierr;
723 
724   PetscFunctionBegin;
725   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
726   PetscValidHeaderSpecific(A, MAT_CLASSID, 9);
727   if (dnz)  PetscValidPointer(dnz,5);  if (onz)  PetscValidPointer(onz,6);
728   if (dnzu) PetscValidPointer(dnzu,7); if (onzu) PetscValidPointer(onzu,8);
729   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
730   ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
731   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
732   ierr = PetscOptionsGetBool(NULL, "-dm_view_preallocation", &debug, NULL);CHKERRQ(ierr);
733   ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
734   ierr = PetscLogEventBegin(DMPLEX_Preallocate,dm,0,0,0);CHKERRQ(ierr);
735   /* Create dof SF based on point SF */
736   if (debug) {
737     PetscSection section, sectionGlobal;
738     PetscSF      sf;
739 
740     ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
741     ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
742     ierr = DMGetDefaultGlobalSection(dm, &sectionGlobal);CHKERRQ(ierr);
743     ierr = PetscPrintf(comm, "Input Section for Preallocation:\n");CHKERRQ(ierr);
744     ierr = PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
745     ierr = PetscPrintf(comm, "Input Global Section for Preallocation:\n");CHKERRQ(ierr);
746     ierr = PetscSectionView(sectionGlobal, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
747     ierr = PetscPrintf(comm, "Input SF for Preallocation:\n");CHKERRQ(ierr);
748     ierr = PetscSFView(sf, NULL);CHKERRQ(ierr);
749   }
750   ierr = PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets);CHKERRQ(ierr);
751   ierr = PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof);CHKERRQ(ierr);
752   if (debug) {
753     ierr = PetscPrintf(comm, "Dof SF for Preallocation:\n");CHKERRQ(ierr);
754     ierr = PetscSFView(sfDof, NULL);CHKERRQ(ierr);
755   }
756   /* Create allocation vectors from adjacency graph */
757   ierr = MatGetLocalSize(A, &locRows, NULL);CHKERRQ(ierr);
758   ierr = PetscLayoutCreate(comm, &rLayout);CHKERRQ(ierr);
759   ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
760   ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
761   ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
762   /* There are 4 types of adjacency */
763   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
764   if (Nf < 1 || bs > 1) {
765     ierr = DMPlexGetAdjacencyUseCone(dm, &useCone);CHKERRQ(ierr);
766     ierr = DMPlexGetAdjacencyUseClosure(dm, &useClosure);CHKERRQ(ierr);
767     idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
768     ierr = DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]);CHKERRQ(ierr);
769     ierr = DMPlexUpdateAllocation_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);CHKERRQ(ierr);
770   } else {
771     for (f = 0; f < Nf; ++f) {
772       ierr = PetscDSGetAdjacency(prob, f, &useCone, &useClosure);CHKERRQ(ierr);
773       idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
774       if (!sectionAdj[idx]) {ierr = DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]);CHKERRQ(ierr);}
775       ierr = DMPlexUpdateAllocation_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);CHKERRQ(ierr);
776     }
777   }
778   ierr = PetscSFDestroy(&sfDof);CHKERRQ(ierr);
779   /* Set matrix pattern */
780   ierr = MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);CHKERRQ(ierr);
781   ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
782   /* Check for symmetric storage */
783   ierr = MatGetType(A, &mtype);CHKERRQ(ierr);
784   ierr = PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);CHKERRQ(ierr);
785   ierr = PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);CHKERRQ(ierr);
786   ierr = PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);CHKERRQ(ierr);
787   if (isSymBlock || isSymSeqBlock || isSymMPIBlock) {ierr = MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);CHKERRQ(ierr);}
788   /* Fill matrix with zeros */
789   if (fillMatrix) {
790     if (Nf < 1 || bs > 1) {
791       ierr = DMPlexGetAdjacencyUseCone(dm, &useCone);CHKERRQ(ierr);
792       ierr = DMPlexGetAdjacencyUseClosure(dm, &useClosure);CHKERRQ(ierr);
793       idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
794       ierr = DMPlexFillMatrix_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], A);CHKERRQ(ierr);
795     } else {
796       for (f = 0; f < Nf; ++f) {
797         ierr = PetscDSGetAdjacency(prob, f, &useCone, &useClosure);CHKERRQ(ierr);
798         idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
799         ierr = DMPlexFillMatrix_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], A);CHKERRQ(ierr);
800       }
801     }
802     ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
803     ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
804   }
805   ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
806   for (idx = 0; idx < 4; ++idx) {ierr = PetscSectionDestroy(&sectionAdj[idx]);CHKERRQ(ierr); ierr = PetscFree(cols[idx]);CHKERRQ(ierr);}
807   ierr = PetscLogEventEnd(DMPLEX_Preallocate,dm,0,0,0);CHKERRQ(ierr);
808   PetscFunctionReturn(0);
809 }
810 
811 #if 0
812 #undef __FUNCT__
813 #define __FUNCT__ "DMPlexPreallocateOperator_2"
814 PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
815 {
816   PetscInt       *tmpClosure,*tmpAdj,*visits;
817   PetscInt        c,cStart,cEnd,pStart,pEnd;
818   PetscErrorCode  ierr;
819 
820   PetscFunctionBegin;
821   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
822   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
823   ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
824 
825   maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1));
826 
827   ierr    = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
828   npoints = pEnd - pStart;
829 
830   ierr = PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits);CHKERRQ(ierr);
831   ierr = PetscMemzero(lvisits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
832   ierr = PetscMemzero(visits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
833   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
834   for (c=cStart; c<cEnd; c++) {
835     PetscInt *support = tmpClosure;
836     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);CHKERRQ(ierr);
837     for (p=0; p<supportSize; p++) lvisits[support[p]]++;
838   }
839   ierr = PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
840   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
841   ierr = PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
842   ierr = PetscSFBcastEnd  (sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
843 
844   ierr = PetscSFGetRanks();CHKERRQ(ierr);
845 
846 
847   ierr = PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner);CHKERRQ(ierr);
848   for (c=cStart; c<cEnd; c++) {
849     ierr = PetscMemzero(cellmat,maxClosureSize*maxClosureSize*sizeof(PetscInt));CHKERRQ(ierr);
850     /*
851      Depth-first walk of transitive closure.
852      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.
853      This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
854      */
855   }
856 
857   ierr = PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);CHKERRQ(ierr);
858   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lonz,onz,MPI_SUM);CHKERRQ(ierr);
859   PetscFunctionReturn(0);
860 }
861 #endif
862