xref: /petsc/src/dm/impls/plex/plexpreallocate.c (revision efd4aadf157bf1ba2d80c2be092fcf4247860003) ! 
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 = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
16   ierr = DMGetDefaultGlobalSection(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);CHKERRQ(ierr);
177   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
178   ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
179   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
180   ierr = DMGetDefaultGlobalSection(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);CHKERRQ(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, PETSC_VIEWER_STDOUT_WORLD);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, PETSC_VIEWER_STDOUT_WORLD);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, PETSC_VIEWER_STDOUT_WORLD);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) {
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, PETSC_VIEWER_STDOUT_WORLD);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, PETSC_VIEWER_STDOUT_WORLD);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 = PetscMemcpy(&cols[aoff], &rootAdj[roff], rdof * sizeof(PetscInt));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 = DMGetDefaultSection(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 = DMPlexGetGlobalFieldOffset_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 = DMGetDefaultSection(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 = DMPlexGetGlobalFieldOffset_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 Arguments:
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   Ouput Argument:
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   PetscErrorCode ierr;
724 
725   PetscFunctionBegin;
726   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
727   PetscValidHeaderSpecific(A, MAT_CLASSID, 9);
728   if (dnz)  PetscValidPointer(dnz,5);  if (onz)  PetscValidPointer(onz,6);
729   if (dnzu) PetscValidPointer(dnzu,7); if (onzu) PetscValidPointer(onzu,8);
730   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
731   ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
732   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
733   ierr = PetscOptionsGetBool(NULL,NULL, "-dm_view_preallocation", &debug, NULL);CHKERRQ(ierr);
734   ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
735   ierr = PetscLogEventBegin(DMPLEX_Preallocate,dm,0,0,0);CHKERRQ(ierr);
736   /* Create dof SF based on point SF */
737   if (debug) {
738     PetscSection section, sectionGlobal;
739     PetscSF      sf;
740 
741     ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
742     ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
743     ierr = DMGetDefaultGlobalSection(dm, &sectionGlobal);CHKERRQ(ierr);
744     ierr = PetscPrintf(comm, "Input Section for Preallocation:\n");CHKERRQ(ierr);
745     ierr = PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
746     ierr = PetscPrintf(comm, "Input Global Section for Preallocation:\n");CHKERRQ(ierr);
747     ierr = PetscSectionView(sectionGlobal, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
748     ierr = PetscPrintf(comm, "Input SF for Preallocation:\n");CHKERRQ(ierr);
749     ierr = PetscSFView(sf, NULL);CHKERRQ(ierr);
750   }
751   ierr = PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets);CHKERRQ(ierr);
752   ierr = PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof);CHKERRQ(ierr);
753   ierr = PetscFree(remoteOffsets);CHKERRQ(ierr);
754   if (debug) {
755     ierr = PetscPrintf(comm, "Dof SF for Preallocation:\n");CHKERRQ(ierr);
756     ierr = PetscSFView(sfDof, NULL);CHKERRQ(ierr);
757   }
758   /* Create allocation vectors from adjacency graph */
759   ierr = MatGetLocalSize(A, &locRows, NULL);CHKERRQ(ierr);
760   ierr = PetscLayoutCreate(comm, &rLayout);CHKERRQ(ierr);
761   ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
762   ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
763   ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
764   /* There are 4 types of adjacency */
765   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
766   if (Nf < 1 || bs > 1) {
767     ierr = DMPlexGetAdjacencyUseCone(dm, &useCone);CHKERRQ(ierr);
768     ierr = DMPlexGetAdjacencyUseClosure(dm, &useClosure);CHKERRQ(ierr);
769     idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
770     ierr = DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]);CHKERRQ(ierr);
771     ierr = DMPlexUpdateAllocation_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);CHKERRQ(ierr);
772   } else {
773     for (f = 0; f < Nf; ++f) {
774       ierr = PetscDSGetAdjacency(prob, f, &useCone, &useClosure);CHKERRQ(ierr);
775       idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
776       if (!sectionAdj[idx]) {ierr = DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]);CHKERRQ(ierr);}
777       ierr = DMPlexUpdateAllocation_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);CHKERRQ(ierr);
778     }
779   }
780   ierr = PetscSFDestroy(&sfDof);CHKERRQ(ierr);
781   /* Set matrix pattern */
782   ierr = MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);CHKERRQ(ierr);
783   ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
784   /* Check for symmetric storage */
785   ierr = MatGetType(A, &mtype);CHKERRQ(ierr);
786   ierr = PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);CHKERRQ(ierr);
787   ierr = PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);CHKERRQ(ierr);
788   ierr = PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);CHKERRQ(ierr);
789   if (isSymBlock || isSymSeqBlock || isSymMPIBlock) {ierr = MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);CHKERRQ(ierr);}
790   /* Fill matrix with zeros */
791   if (fillMatrix) {
792     if (Nf < 1 || bs > 1) {
793       ierr = DMPlexGetAdjacencyUseCone(dm, &useCone);CHKERRQ(ierr);
794       ierr = DMPlexGetAdjacencyUseClosure(dm, &useClosure);CHKERRQ(ierr);
795       idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
796       ierr = DMPlexFillMatrix_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], A);CHKERRQ(ierr);
797     } else {
798       for (f = 0; f < Nf; ++f) {
799         ierr = PetscDSGetAdjacency(prob, f, &useCone, &useClosure);CHKERRQ(ierr);
800         idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
801         ierr = DMPlexFillMatrix_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], A);CHKERRQ(ierr);
802       }
803     }
804     ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
805     ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
806   }
807   ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
808   for (idx = 0; idx < 4; ++idx) {ierr = PetscSectionDestroy(&sectionAdj[idx]);CHKERRQ(ierr); ierr = PetscFree(cols[idx]);CHKERRQ(ierr);}
809   ierr = PetscLogEventEnd(DMPLEX_Preallocate,dm,0,0,0);CHKERRQ(ierr);
810   PetscFunctionReturn(0);
811 }
812 
813 #if 0
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