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