xref: /petsc/src/dm/impls/plex/plexpreallocate.c (revision a9fb644384c31e1fad7fabb80d80dc0b100418b9) 
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, PetscSection section, PetscSection sectionGlobal, PetscBool useCone, PetscBool useClosure, PetscSection *anchorSectionAdj, PetscInt *anchorAdj[])
10 {
11   PetscInt       pStart, pEnd;
12   PetscSection   adjSec, aSec;
13   IS             aIS;
14   PetscErrorCode ierr;
15 
16   PetscFunctionBegin;
17 
18   ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&adjSec);CHKERRQ(ierr);
19   ierr = PetscSectionGetChart(section,&pStart,&pEnd);CHKERRQ(ierr);
20   ierr = PetscSectionSetChart(adjSec,pStart,pEnd);CHKERRQ(ierr);
21 
22   ierr = DMPlexGetAnchors(dm,&aSec,&aIS);CHKERRQ(ierr);
23   if (aSec) {
24     const PetscInt *anchors;
25     PetscInt       p, q, a, aSize, *offsets, aStart, aEnd, *inverse, iSize, *adj, adjSize;
26     PetscInt       *tmpAdjP = NULL, *tmpAdjQ = NULL;
27     PetscSection   inverseSec;
28 
29     /* invert the constraint-to-anchor map */
30     ierr = PetscSectionCreate(PetscObjectComm((PetscObject)aSec),&inverseSec);CHKERRQ(ierr);
31     ierr = PetscSectionSetChart(inverseSec,pStart,pEnd);CHKERRQ(ierr);
32     ierr = ISGetLocalSize(aIS, &aSize);CHKERRQ(ierr);
33     ierr = ISGetIndices(aIS, &anchors);CHKERRQ(ierr);
34 
35     for (p = 0; p < aSize; p++) {
36       PetscInt a = anchors[p];
37 
38       ierr = PetscSectionAddDof(inverseSec,a,1);CHKERRQ(ierr);
39     }
40     ierr = PetscSectionSetUp(inverseSec);CHKERRQ(ierr);
41     ierr = PetscSectionGetStorageSize(inverseSec,&iSize);CHKERRQ(ierr);
42     ierr = PetscMalloc1(iSize,&inverse);CHKERRQ(ierr);
43     ierr = PetscCalloc1(pEnd-pStart,&offsets);CHKERRQ(ierr);
44     ierr = PetscSectionGetChart(aSec,&aStart,&aEnd);CHKERRQ(ierr);
45     for (p = aStart; p < aEnd; p++) {
46       PetscInt dof, off;
47 
48       ierr = PetscSectionGetDof(aSec, p, &dof);CHKERRQ(ierr);
49       ierr = PetscSectionGetOffset(aSec, p, &off);CHKERRQ(ierr);
50 
51       for (q = 0; q < dof; q++) {
52         PetscInt iOff;
53 
54         a = anchors[off + q];
55         ierr = PetscSectionGetOffset(inverseSec, a, &iOff);CHKERRQ(ierr);
56         inverse[iOff + offsets[a-pStart]++] = p;
57       }
58     }
59     ierr = ISRestoreIndices(aIS, &anchors);CHKERRQ(ierr);
60     ierr = PetscFree(offsets);CHKERRQ(ierr);
61 
62     /* construct anchorAdj and adjSec
63      *
64      * loop over anchors:
65      *   construct anchor adjacency
66      *   loop over constrained:
67      *     construct constrained adjacency
68      *     if not in anchor adjacency, add to dofs
69      * setup adjSec, allocate anchorAdj
70      * loop over anchors:
71      *   construct anchor adjacency
72      *   loop over constrained:
73      *     construct constrained adjacency
74      *     if not in anchor adjacency
75      *       if not already in list, put in list
76      *   sort, unique, reduce dof count
77      * optional: compactify
78      */
79     for (p = pStart; p < pEnd; p++) {
80       PetscInt iDof, iOff, i, r, s, numAdjP = PETSC_DETERMINE;
81 
82       ierr = PetscSectionGetDof(inverseSec,p,&iDof);CHKERRQ(ierr);
83       if (!iDof) continue;
84       ierr = PetscSectionGetOffset(inverseSec,p,&iOff);CHKERRQ(ierr);
85       ierr = DMPlexGetAdjacency_Internal(dm,p,useCone,useClosure,PETSC_TRUE,&numAdjP,&tmpAdjP);CHKERRQ(ierr);
86       for (i = 0; i < iDof; i++) {
87         PetscInt iNew = 0, qAdj, qAdjDof, qAdjCDof, numAdjQ = PETSC_DETERMINE;
88 
89         q = inverse[iOff + i];
90         ierr = DMPlexGetAdjacency_Internal(dm,q,useCone,useClosure,PETSC_TRUE,&numAdjQ,&tmpAdjQ);CHKERRQ(ierr);
91         for (r = 0; r < numAdjQ; r++) {
92           qAdj = tmpAdjQ[r];
93           if ((qAdj < pStart) || (qAdj >= pEnd)) continue;
94           for (s = 0; s < numAdjP; s++) {
95             if (qAdj == tmpAdjP[s]) break;
96           }
97           if (s < numAdjP) continue;
98           ierr  = PetscSectionGetDof(section,qAdj,&qAdjDof);CHKERRQ(ierr);
99           ierr  = PetscSectionGetConstraintDof(section,qAdj,&qAdjCDof);CHKERRQ(ierr);
100           iNew += qAdjDof - qAdjCDof;
101         }
102         ierr = PetscSectionAddDof(adjSec,p,iNew);CHKERRQ(ierr);
103       }
104     }
105 
106     ierr = PetscSectionSetUp(adjSec);CHKERRQ(ierr);
107     ierr = PetscSectionGetStorageSize(adjSec,&adjSize);CHKERRQ(ierr);
108     ierr = PetscMalloc1(adjSize,&adj);CHKERRQ(ierr);
109 
110     for (p = pStart; p < pEnd; p++) {
111       PetscInt iDof, iOff, i, r, s, aOff, aOffOrig, aDof, numAdjP = PETSC_DETERMINE;
112 
113       ierr = PetscSectionGetDof(inverseSec,p,&iDof);CHKERRQ(ierr);
114       if (!iDof) continue;
115       ierr = PetscSectionGetOffset(inverseSec,p,&iOff);CHKERRQ(ierr);
116       ierr = DMPlexGetAdjacency_Internal(dm,p,useCone,useClosure,PETSC_TRUE,&numAdjP,&tmpAdjP);CHKERRQ(ierr);
117       ierr = PetscSectionGetDof(adjSec,p,&aDof);CHKERRQ(ierr);
118       ierr = PetscSectionGetOffset(adjSec,p,&aOff);CHKERRQ(ierr);
119       aOffOrig = aOff;
120       for (i = 0; i < iDof; i++) {
121         PetscInt qAdj, qAdjDof, qAdjCDof, qAdjOff, nd, numAdjQ = PETSC_DETERMINE;
122 
123         q = inverse[iOff + i];
124         ierr = DMPlexGetAdjacency_Internal(dm,q,useCone,useClosure,PETSC_TRUE,&numAdjQ,&tmpAdjQ);CHKERRQ(ierr);
125         for (r = 0; r < numAdjQ; r++) {
126           qAdj = tmpAdjQ[r];
127           if ((qAdj < pStart) || (qAdj >= pEnd)) continue;
128           for (s = 0; s < numAdjP; s++) {
129             if (qAdj == tmpAdjP[s]) break;
130           }
131           if (s < numAdjP) continue;
132           ierr  = PetscSectionGetDof(section,qAdj,&qAdjDof);CHKERRQ(ierr);
133           ierr  = PetscSectionGetConstraintDof(section,qAdj,&qAdjCDof);CHKERRQ(ierr);
134           ierr  = PetscSectionGetOffset(sectionGlobal,qAdj,&qAdjOff);CHKERRQ(ierr);
135           for (nd = 0; nd < qAdjDof-qAdjCDof; ++nd) {
136             adj[aOff++] = (qAdjOff < 0 ? -(qAdjOff+1) : qAdjOff) + nd;
137           }
138         }
139       }
140       ierr = PetscSortRemoveDupsInt(&aDof,&adj[aOffOrig]);CHKERRQ(ierr);
141       ierr = PetscSectionSetDof(adjSec,p,aDof);CHKERRQ(ierr);
142     }
143     *anchorAdj = adj;
144 
145     /* clean up */
146     ierr = PetscSectionDestroy(&inverseSec);CHKERRQ(ierr);
147     ierr = PetscFree(inverse);CHKERRQ(ierr);
148     ierr = PetscFree(tmpAdjP);CHKERRQ(ierr);
149     ierr = PetscFree(tmpAdjQ);CHKERRQ(ierr);
150   }
151   else {
152     *anchorAdj = NULL;
153     ierr = PetscSectionSetUp(adjSec);CHKERRQ(ierr);
154   }
155   *anchorSectionAdj = adjSec;
156   PetscFunctionReturn(0);
157 }
158 
159 #undef __FUNCT__
160 #define __FUNCT__ "DMPlexCreateAdjacencySection_Static"
161 PetscErrorCode DMPlexCreateAdjacencySection_Static(DM dm, PetscInt bs, PetscSF sfDof, PetscSection section, PetscSection sectionGlobal, PetscBool useCone, PetscBool useClosure, PetscBool useAnchors, PetscSection *sA, PetscInt **colIdx)
162 {
163   MPI_Comm           comm;
164   PetscMPIInt        size;
165   PetscBool          doCommLocal, doComm, debug = PETSC_FALSE;
166   PetscSF            sf, sfAdj;
167   PetscSection       leafSectionAdj, rootSectionAdj, sectionAdj, anchorSectionAdj;
168   PetscInt           nroots, nleaves, l, p, r;
169   const PetscInt    *leaves;
170   const PetscSFNode *remotes;
171   PetscInt           dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols;
172   PetscInt          *tmpAdj = NULL, *adj, *rootAdj, *anchorAdj = NULL, *cols, *remoteOffsets;
173   PetscInt           adjSize;
174   PetscErrorCode     ierr;
175 
176   PetscFunctionBegin;
177   ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
178   ierr = PetscOptionsGetBool(NULL, "-dm_view_preallocation", &debug, NULL);CHKERRQ(ierr);
179   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
180   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
181   ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
182   ierr = PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);CHKERRQ(ierr);
183   doCommLocal = (size > 1) && (nroots >= 0) ? PETSC_TRUE : PETSC_FALSE;
184   ierr = MPI_Allreduce(&doCommLocal, &doComm, 1, MPIU_BOOL, MPI_LAND, comm);CHKERRQ(ierr);
185   /* Create section for dof adjacency (dof ==> # adj dof) */
186   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
187   ierr = PetscSectionGetStorageSize(section, &numDof);CHKERRQ(ierr);
188   ierr = PetscSectionCreate(comm, &leafSectionAdj);CHKERRQ(ierr);
189   ierr = PetscSectionSetChart(leafSectionAdj, 0, numDof);CHKERRQ(ierr);
190   ierr = PetscSectionCreate(comm, &rootSectionAdj);CHKERRQ(ierr);
191   ierr = PetscSectionSetChart(rootSectionAdj, 0, numDof);CHKERRQ(ierr);
192   /*   Fill in the ghost dofs on the interface */
193   ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes);CHKERRQ(ierr);
194   /*
195    section        - maps points to (# dofs, local dofs)
196    sectionGlobal  - maps points to (# dofs, global dofs)
197    leafSectionAdj - maps unowned local dofs to # adj dofs
198    rootSectionAdj - maps   owned local dofs to # adj dofs
199    adj            - adj global dofs indexed by leafSectionAdj
200    rootAdj        - adj global dofs indexed by rootSectionAdj
201    sf    - describes shared points across procs
202    sfDof - describes shared dofs across procs
203    sfAdj - describes shared adjacent dofs across procs
204    ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point.
205   (0). If there are point-to-point constraints, add the adjacencies of constrained points to anchors in anchorAdj
206        (This is done in DMPlexComputeAnchorAdjacencies())
207     1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj
208        Reduce those counts to rootSectionAdj (now redundantly counting some interface points)
209     2. Visit owned points on interface, count adjacencies placing in rootSectionAdj
210        Create sfAdj connecting rootSectionAdj and leafSectionAdj
211     3. Visit unowned points on interface, write adjacencies to adj
212        Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies)
213     4. Visit owned points on interface, write adjacencies to rootAdj
214        Remove redundancy in rootAdj
215    ** The last two traversals use transitive closure
216     5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj)
217        Allocate memory addressed by sectionAdj (cols)
218     6. Visit all owned points in the subdomain, insert dof adjacencies into cols
219    ** Knowing all the column adjacencies, check ownership and sum into dnz and onz
220   */
221   ierr = DMPlexComputeAnchorAdjacencies(dm, section, sectionGlobal, useCone, useClosure, &anchorSectionAdj, &anchorAdj);CHKERRQ(ierr);
222   for (l = 0; l < nleaves; ++l) {
223     PetscInt dof, off, d, q, anDof;
224     PetscInt p = leaves[l], numAdj = PETSC_DETERMINE;
225 
226     if ((p < pStart) || (p >= pEnd)) continue;
227     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
228     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
229     ierr = DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);CHKERRQ(ierr);
230     for (q = 0; q < numAdj; ++q) {
231       const PetscInt padj = tmpAdj[q];
232       PetscInt ndof, ncdof;
233 
234       if ((padj < pStart) || (padj >= pEnd)) continue;
235       ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
236       ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
237       for (d = off; d < off+dof; ++d) {
238         ierr = PetscSectionAddDof(leafSectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
239       }
240     }
241     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
242     if (anDof) {
243       for (d = off; d < off+dof; ++d) {
244         ierr = PetscSectionAddDof(leafSectionAdj, d, anDof);CHKERRQ(ierr);
245       }
246     }
247   }
248   ierr = PetscSectionSetUp(leafSectionAdj);CHKERRQ(ierr);
249   if (debug) {
250     ierr = PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n");CHKERRQ(ierr);
251     ierr = PetscSectionView(leafSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
252   }
253   /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */
254   if (doComm) {
255     ierr = PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);CHKERRQ(ierr);
256     ierr = PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);CHKERRQ(ierr);
257   }
258   if (debug) {
259     ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots:\n");CHKERRQ(ierr);
260     ierr = PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
261   }
262   /* Add in local adjacency sizes for owned dofs on interface (roots) */
263   for (p = pStart; p < pEnd; ++p) {
264     PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof;
265 
266     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
267     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
268     if (!dof) continue;
269     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
270     if (adof <= 0) continue;
271     ierr = DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);CHKERRQ(ierr);
272     for (q = 0; q < numAdj; ++q) {
273       const PetscInt padj = tmpAdj[q];
274       PetscInt ndof, ncdof;
275 
276       if ((padj < pStart) || (padj >= pEnd)) continue;
277       ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
278       ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
279       for (d = off; d < off+dof; ++d) {
280         ierr = PetscSectionAddDof(rootSectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
281       }
282     }
283     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
284     if (anDof) {
285       for (d = off; d < off+dof; ++d) {
286         ierr = PetscSectionAddDof(rootSectionAdj, d, anDof);CHKERRQ(ierr);
287       }
288     }
289   }
290   ierr = PetscSectionSetUp(rootSectionAdj);CHKERRQ(ierr);
291   if (debug) {
292     ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after local additions:\n");CHKERRQ(ierr);
293     ierr = PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
294   }
295   /* Create adj SF based on dof SF */
296   ierr = PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets);CHKERRQ(ierr);
297   ierr = PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj);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 #undef __FUNCT__
585 #define __FUNCT__ "DMPlexUpdateAllocation_Static"
586 PetscErrorCode DMPlexUpdateAllocation_Static(PetscLayout rLayout, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[])
587 {
588   PetscInt       rStart, rEnd, r, pStart, pEnd, p;
589   PetscErrorCode ierr;
590 
591   PetscFunctionBegin;
592   /* This loop needs to change to a loop over points, then field dofs, which means we need to look both sections */
593   ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
594   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);
595   if (f >= 0 && bs == 1) {
596     ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
597     for (p = pStart; p < pEnd; ++p) {
598       PetscInt goff, loff, lfoff, fdof, fcdof, rS, rE;
599 
600       ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
601       ierr = PetscSectionGetOffset(section, p, &loff);CHKERRQ(ierr);
602       ierr = PetscSectionGetFieldOffset(section, p, f, &lfoff);CHKERRQ(ierr);
603       ierr = PetscSectionGetFieldDof(section, p, f, &fdof);CHKERRQ(ierr);
604       ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr);
605       rS   = goff + (lfoff - loff);
606       rE   = rS + fdof - fcdof;
607       for (r = rS; r < rE; ++r) {
608         PetscInt numCols, cStart, c;
609 
610         ierr = PetscSectionGetDof(sectionAdj, r, &numCols);CHKERRQ(ierr);
611         ierr = PetscSectionGetOffset(sectionAdj, r, &cStart);CHKERRQ(ierr);
612         for (c = cStart; c < cStart+numCols; ++c) {
613           if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
614             ++dnz[r-rStart];
615             if (cols[c] >= r) ++dnzu[r-rStart];
616           } else {
617             ++onz[r-rStart];
618             if (cols[c] >= r) ++onzu[r-rStart];
619           }
620         }
621       }
622     }
623   } else {
624     /* Only loop over blocks of rows */
625     for (r = rStart/bs; r < rEnd/bs; ++r) {
626       const PetscInt row = r*bs;
627       PetscInt       numCols, cStart, c;
628 
629       ierr = PetscSectionGetDof(sectionAdj, row, &numCols);CHKERRQ(ierr);
630       ierr = PetscSectionGetOffset(sectionAdj, row, &cStart);CHKERRQ(ierr);
631       for (c = cStart; c < cStart+numCols; ++c) {
632         if ((cols[c] >= rStart*bs) && (cols[c] < rEnd*bs)) {
633           ++dnz[r-rStart];
634           if (cols[c] >= row) ++dnzu[r-rStart];
635         } else {
636           ++onz[r-rStart];
637           if (cols[c] >= row) ++onzu[r-rStart];
638         }
639       }
640     }
641     for (r = 0; r < (rEnd - rStart)/bs; ++r) {
642       dnz[r]  /= bs;
643       onz[r]  /= bs;
644       dnzu[r] /= bs;
645       onzu[r] /= bs;
646     }
647   }
648   PetscFunctionReturn(0);
649 }
650 
651 #undef __FUNCT__
652 #define __FUNCT__ "DMPlexFillMatrix_Static"
653 PetscErrorCode DMPlexFillMatrix_Static(PetscLayout rLayout, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], Mat A)
654 {
655   PetscScalar   *values;
656   PetscInt       rStart, rEnd, r, pStart, pEnd, p, len, maxRowLen = 0;
657   PetscErrorCode ierr;
658 
659   PetscFunctionBegin;
660   ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
661   for (r = rStart; r < rEnd; ++r) {
662     ierr      = PetscSectionGetDof(sectionAdj, r, &len);CHKERRQ(ierr);
663     maxRowLen = PetscMax(maxRowLen, len);
664   }
665   ierr = PetscCalloc1(maxRowLen, &values);CHKERRQ(ierr);
666   if (f >=0 && bs == 1) {
667     ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
668     for (p = pStart; p < pEnd; ++p) {
669       PetscInt goff, loff, lfoff, fdof, fcdof, rS, rE;
670 
671       ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
672       ierr = PetscSectionGetOffset(section, p, &loff);CHKERRQ(ierr);
673       ierr = PetscSectionGetFieldOffset(section, p, f, &lfoff);CHKERRQ(ierr);
674       ierr = PetscSectionGetFieldDof(section, p, f, &fdof);CHKERRQ(ierr);
675       ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr);
676       rS   = goff + (lfoff - loff);
677       rE   = rS + fdof - fcdof;
678       for (r = rS; r < rE; ++r) {
679         PetscInt numCols, cStart, c;
680 
681         ierr = PetscSectionGetDof(sectionAdj, r, &numCols);CHKERRQ(ierr);
682         ierr = PetscSectionGetOffset(sectionAdj, r, &cStart);CHKERRQ(ierr);
683         ierr = MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);CHKERRQ(ierr);
684       }
685     }
686   } else {
687     for (r = rStart; r < rEnd; ++r) {
688       PetscInt numCols, cStart;
689 
690       ierr = PetscSectionGetDof(sectionAdj, r, &numCols);CHKERRQ(ierr);
691       ierr = PetscSectionGetOffset(sectionAdj, r, &cStart);CHKERRQ(ierr);
692       ierr = MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);CHKERRQ(ierr);
693     }
694   }
695   ierr = PetscFree(values);CHKERRQ(ierr);
696   PetscFunctionReturn(0);
697 }
698 
699 #undef __FUNCT__
700 #define __FUNCT__ "DMPlexPreallocateOperator"
701 PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
702 {
703   MPI_Comm       comm;
704   PetscDS        prob;
705   MatType        mtype;
706   PetscSF        sf, sfDof;
707   PetscInt      *remoteOffsets;
708   PetscSection   sectionAdj[4] = {NULL, NULL, NULL, NULL};
709   PetscInt      *cols[4]       = {NULL, NULL, NULL, NULL};
710   PetscBool      useCone, useClosure;
711   PetscInt       Nf, f, idx, locRows, r;
712   PetscLayout    rLayout;
713   PetscBool      isSymBlock, isSymSeqBlock, isSymMPIBlock, debug = PETSC_FALSE;
714   PetscErrorCode ierr;
715 
716   PetscFunctionBegin;
717   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
718   PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 3);
719   PetscValidHeaderSpecific(sectionGlobal, PETSC_SECTION_CLASSID, 4);
720   PetscValidHeaderSpecific(A, MAT_CLASSID, 9);
721   if (dnz)  PetscValidPointer(dnz,5);  if (onz)  PetscValidPointer(onz,6);
722   if (dnzu) PetscValidPointer(dnzu,7); if (onzu) PetscValidPointer(onzu,8);
723   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
724   ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
725   ierr = PetscOptionsGetBool(NULL, "-dm_view_preallocation", &debug, NULL);CHKERRQ(ierr);
726   ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
727   ierr = PetscLogEventBegin(DMPLEX_Preallocate,dm,0,0,0);CHKERRQ(ierr);
728   /* Create dof SF based on point SF */
729   if (debug) {
730     PetscSF sf;
731 
732     ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
733     ierr = PetscPrintf(comm, "Input Section for Preallocation:\n");CHKERRQ(ierr);
734     ierr = PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
735     ierr = PetscPrintf(comm, "Input Global Section for Preallocation:\n");CHKERRQ(ierr);
736     ierr = PetscSectionView(sectionGlobal, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
737     ierr = PetscPrintf(comm, "Input SF for Preallocation:\n");CHKERRQ(ierr);
738     ierr = PetscSFView(sf, NULL);CHKERRQ(ierr);
739   }
740   ierr = PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets);CHKERRQ(ierr);
741   ierr = PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof);CHKERRQ(ierr);
742   if (debug) {
743     ierr = PetscPrintf(comm, "Dof SF for Preallocation:\n");CHKERRQ(ierr);
744     ierr = PetscSFView(sfDof, NULL);CHKERRQ(ierr);
745   }
746   /* Create allocation vectors from adjacency graph */
747   ierr = MatGetLocalSize(A, &locRows, NULL);CHKERRQ(ierr);
748   ierr = PetscLayoutCreate(comm, &rLayout);CHKERRQ(ierr);
749   ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
750   ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
751   ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
752   /* There are 4 types of adjacency */
753   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
754   if (Nf < 1 || bs > 1) {
755     ierr = DMPlexGetAdjacencyUseCone(dm, &useCone);CHKERRQ(ierr);
756     ierr = DMPlexGetAdjacencyUseClosure(dm, &useClosure);CHKERRQ(ierr);
757     idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
758     ierr = DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, section, sectionGlobal, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]);CHKERRQ(ierr);
759     ierr = DMPlexUpdateAllocation_Static(rLayout, bs, section, sectionGlobal, -1, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);CHKERRQ(ierr);
760   } else {
761     for (f = 0; f < Nf; ++f) {
762       ierr = PetscDSGetAdjacency(prob, f, &useCone, &useClosure);CHKERRQ(ierr);
763       idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
764       if (!sectionAdj[idx]) {ierr = DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, section, sectionGlobal, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]);CHKERRQ(ierr);}
765       ierr = DMPlexUpdateAllocation_Static(rLayout, bs, section, sectionGlobal, f, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);CHKERRQ(ierr);
766     }
767   }
768   ierr = PetscSFDestroy(&sfDof);CHKERRQ(ierr);
769   /* Set matrix pattern */
770   ierr = MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);CHKERRQ(ierr);
771   ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
772   /* Check for symmetric storage */
773   ierr = MatGetType(A, &mtype);CHKERRQ(ierr);
774   ierr = PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);CHKERRQ(ierr);
775   ierr = PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);CHKERRQ(ierr);
776   ierr = PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);CHKERRQ(ierr);
777   if (isSymBlock || isSymSeqBlock || isSymMPIBlock) {ierr = MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);CHKERRQ(ierr);}
778   /* Fill matrix with zeros */
779   if (fillMatrix) {
780     if (Nf < 1 || bs > 1) {
781       ierr = DMPlexGetAdjacencyUseCone(dm, &useCone);CHKERRQ(ierr);
782       ierr = DMPlexGetAdjacencyUseClosure(dm, &useClosure);CHKERRQ(ierr);
783       idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
784       ierr = DMPlexFillMatrix_Static(rLayout, bs, section, sectionGlobal, -1, sectionAdj[idx], cols[idx], A);CHKERRQ(ierr);
785     } else {
786       for (f = 0; f < Nf; ++f) {
787         ierr = PetscDSGetAdjacency(prob, f, &useCone, &useClosure);CHKERRQ(ierr);
788         idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
789         ierr = DMPlexFillMatrix_Static(rLayout, bs, section, sectionGlobal, f, sectionAdj[idx], cols[idx], A);CHKERRQ(ierr);
790       }
791     }
792     ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
793     ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
794   }
795   ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
796   for (idx = 0; idx < 4; ++idx) {ierr = PetscSectionDestroy(&sectionAdj[idx]);CHKERRQ(ierr); ierr = PetscFree(cols[idx]);CHKERRQ(ierr);}
797   ierr = PetscLogEventEnd(DMPLEX_Preallocate,dm,0,0,0);CHKERRQ(ierr);
798   PetscFunctionReturn(0);
799 }
800 
801 #if 0
802 #undef __FUNCT__
803 #define __FUNCT__ "DMPlexPreallocateOperator_2"
804 PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
805 {
806   PetscInt       *tmpClosure,*tmpAdj,*visits;
807   PetscInt        c,cStart,cEnd,pStart,pEnd;
808   PetscErrorCode  ierr;
809 
810   PetscFunctionBegin;
811   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
812   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
813   ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
814 
815   maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1));
816 
817   ierr    = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
818   npoints = pEnd - pStart;
819 
820   ierr = PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits);CHKERRQ(ierr);
821   ierr = PetscMemzero(lvisits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
822   ierr = PetscMemzero(visits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
823   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
824   for (c=cStart; c<cEnd; c++) {
825     PetscInt *support = tmpClosure;
826     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);CHKERRQ(ierr);
827     for (p=0; p<supportSize; p++) lvisits[support[p]]++;
828   }
829   ierr = PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
830   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
831   ierr = PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
832   ierr = PetscSFBcastEnd  (sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
833 
834   ierr = PetscSFGetRanks();CHKERRQ(ierr);
835 
836 
837   ierr = PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner);CHKERRQ(ierr);
838   for (c=cStart; c<cEnd; c++) {
839     ierr = PetscMemzero(cellmat,maxClosureSize*maxClosureSize*sizeof(PetscInt));CHKERRQ(ierr);
840     /*
841      Depth-first walk of transitive closure.
842      At each leaf frame f of transitive closure that we see, add 1/visits[f] to each pair (p,q) not marked as done in cellmat.
843      This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
844      */
845   }
846 
847   ierr = PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);CHKERRQ(ierr);
848   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lonz,onz,MPI_SUM);CHKERRQ(ierr);
849   PetscFunctionReturn(0);
850 }
851 #endif
852