xref: /petsc/src/dm/impls/plex/plexpreallocate.c (revision 2f6eced2a19e978d64f27de66fbc6c26cc5c7934) 
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 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 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 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 PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
688 {
689   MPI_Comm       comm;
690   PetscDS        prob;
691   MatType        mtype;
692   PetscSF        sf, sfDof;
693   PetscSection   section;
694   PetscInt      *remoteOffsets;
695   PetscSection   sectionAdj[4] = {NULL, NULL, NULL, NULL};
696   PetscInt      *cols[4]       = {NULL, NULL, NULL, NULL};
697   PetscBool      useCone, useClosure;
698   PetscInt       Nf, f, idx, locRows;
699   PetscLayout    rLayout;
700   PetscBool      isSymBlock, isSymSeqBlock, isSymMPIBlock, debug = PETSC_FALSE;
701   PetscErrorCode ierr;
702 
703   PetscFunctionBegin;
704   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
705   PetscValidHeaderSpecific(A, MAT_CLASSID, 9);
706   if (dnz)  PetscValidPointer(dnz,5);  if (onz)  PetscValidPointer(onz,6);
707   if (dnzu) PetscValidPointer(dnzu,7); if (onzu) PetscValidPointer(onzu,8);
708   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
709   ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
710   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
711   ierr = PetscOptionsGetBool(NULL,NULL, "-dm_view_preallocation", &debug, NULL);CHKERRQ(ierr);
712   ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
713   ierr = PetscLogEventBegin(DMPLEX_Preallocate,dm,0,0,0);CHKERRQ(ierr);
714   /* Create dof SF based on point SF */
715   if (debug) {
716     PetscSection section, sectionGlobal;
717     PetscSF      sf;
718 
719     ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
720     ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
721     ierr = DMGetDefaultGlobalSection(dm, &sectionGlobal);CHKERRQ(ierr);
722     ierr = PetscPrintf(comm, "Input Section for Preallocation:\n");CHKERRQ(ierr);
723     ierr = PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
724     ierr = PetscPrintf(comm, "Input Global Section for Preallocation:\n");CHKERRQ(ierr);
725     ierr = PetscSectionView(sectionGlobal, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
726     ierr = PetscPrintf(comm, "Input SF for Preallocation:\n");CHKERRQ(ierr);
727     ierr = PetscSFView(sf, NULL);CHKERRQ(ierr);
728   }
729   ierr = PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets);CHKERRQ(ierr);
730   ierr = PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof);CHKERRQ(ierr);
731   ierr = PetscFree(remoteOffsets);CHKERRQ(ierr);
732   if (debug) {
733     ierr = PetscPrintf(comm, "Dof SF for Preallocation:\n");CHKERRQ(ierr);
734     ierr = PetscSFView(sfDof, NULL);CHKERRQ(ierr);
735   }
736   /* Create allocation vectors from adjacency graph */
737   ierr = MatGetLocalSize(A, &locRows, NULL);CHKERRQ(ierr);
738   ierr = PetscLayoutCreate(comm, &rLayout);CHKERRQ(ierr);
739   ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
740   ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
741   ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
742   /* There are 4 types of adjacency */
743   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
744   if (Nf < 1 || bs > 1) {
745     ierr = DMPlexGetAdjacencyUseCone(dm, &useCone);CHKERRQ(ierr);
746     ierr = DMPlexGetAdjacencyUseClosure(dm, &useClosure);CHKERRQ(ierr);
747     idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
748     ierr = DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]);CHKERRQ(ierr);
749     ierr = DMPlexUpdateAllocation_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);CHKERRQ(ierr);
750   } else {
751     for (f = 0; f < Nf; ++f) {
752       ierr = PetscDSGetAdjacency(prob, f, &useCone, &useClosure);CHKERRQ(ierr);
753       idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
754       if (!sectionAdj[idx]) {ierr = DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]);CHKERRQ(ierr);}
755       ierr = DMPlexUpdateAllocation_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);CHKERRQ(ierr);
756     }
757   }
758   ierr = PetscSFDestroy(&sfDof);CHKERRQ(ierr);
759   /* Set matrix pattern */
760   ierr = MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);CHKERRQ(ierr);
761   ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
762   /* Check for symmetric storage */
763   ierr = MatGetType(A, &mtype);CHKERRQ(ierr);
764   ierr = PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);CHKERRQ(ierr);
765   ierr = PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);CHKERRQ(ierr);
766   ierr = PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);CHKERRQ(ierr);
767   if (isSymBlock || isSymSeqBlock || isSymMPIBlock) {ierr = MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);CHKERRQ(ierr);}
768   /* Fill matrix with zeros */
769   if (fillMatrix) {
770     if (Nf < 1 || bs > 1) {
771       ierr = DMPlexGetAdjacencyUseCone(dm, &useCone);CHKERRQ(ierr);
772       ierr = DMPlexGetAdjacencyUseClosure(dm, &useClosure);CHKERRQ(ierr);
773       idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
774       ierr = DMPlexFillMatrix_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], A);CHKERRQ(ierr);
775     } else {
776       for (f = 0; f < Nf; ++f) {
777         ierr = PetscDSGetAdjacency(prob, f, &useCone, &useClosure);CHKERRQ(ierr);
778         idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
779         ierr = DMPlexFillMatrix_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], A);CHKERRQ(ierr);
780       }
781     }
782     ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
783     ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
784   }
785   ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
786   for (idx = 0; idx < 4; ++idx) {ierr = PetscSectionDestroy(&sectionAdj[idx]);CHKERRQ(ierr); ierr = PetscFree(cols[idx]);CHKERRQ(ierr);}
787   ierr = PetscLogEventEnd(DMPLEX_Preallocate,dm,0,0,0);CHKERRQ(ierr);
788   PetscFunctionReturn(0);
789 }
790 
791 #if 0
792 PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
793 {
794   PetscInt       *tmpClosure,*tmpAdj,*visits;
795   PetscInt        c,cStart,cEnd,pStart,pEnd;
796   PetscErrorCode  ierr;
797 
798   PetscFunctionBegin;
799   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
800   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
801   ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
802 
803   maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1));
804 
805   ierr    = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
806   npoints = pEnd - pStart;
807 
808   ierr = PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits);CHKERRQ(ierr);
809   ierr = PetscMemzero(lvisits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
810   ierr = PetscMemzero(visits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
811   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
812   for (c=cStart; c<cEnd; c++) {
813     PetscInt *support = tmpClosure;
814     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);CHKERRQ(ierr);
815     for (p=0; p<supportSize; p++) lvisits[support[p]]++;
816   }
817   ierr = PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
818   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
819   ierr = PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
820   ierr = PetscSFBcastEnd  (sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
821 
822   ierr = PetscSFGetRanks();CHKERRQ(ierr);
823 
824 
825   ierr = PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner);CHKERRQ(ierr);
826   for (c=cStart; c<cEnd; c++) {
827     ierr = PetscMemzero(cellmat,maxClosureSize*maxClosureSize*sizeof(PetscInt));CHKERRQ(ierr);
828     /*
829      Depth-first walk of transitive closure.
830      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.
831      This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
832      */
833   }
834 
835   ierr = PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);CHKERRQ(ierr);
836   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lonz,onz,MPI_SUM);CHKERRQ(ierr);
837   PetscFunctionReturn(0);
838 }
839 #endif
840