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() */
DMPlexComputeAnchorAdjacencies(DM dm,PetscBool useCone,PetscBool useClosure,PetscSection * anchorSectionAdj,PetscInt * anchorAdj[])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 = DMGetLocalSection(dm, &section);CHKERRQ(ierr);
16   ierr = DMGetGlobalSection(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 
DMPlexCreateAdjacencySection_Static(DM dm,PetscInt bs,PetscSF sfDof,PetscBool useCone,PetscBool useClosure,PetscBool useAnchors,PetscSection * sA,PetscInt ** colIdx)158 static PetscErrorCode DMPlexCreateAdjacencySection_Static(DM dm, PetscInt bs, PetscSF sfDof, PetscBool useCone, PetscBool useClosure, PetscBool useAnchors, PetscSection *sA, PetscInt **colIdx)
159 {
160   MPI_Comm           comm;
161   PetscMPIInt        size;
162   PetscBool          doCommLocal, doComm, debug = PETSC_FALSE;
163   PetscSF            sf, sfAdj;
164   PetscSection       section, sectionGlobal, leafSectionAdj, rootSectionAdj, sectionAdj, anchorSectionAdj;
165   PetscInt           nroots, nleaves, l, p, r;
166   const PetscInt    *leaves;
167   const PetscSFNode *remotes;
168   PetscInt           dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols;
169   PetscInt          *tmpAdj = NULL, *adj, *rootAdj, *anchorAdj = NULL, *cols, *remoteOffsets;
170   PetscInt           adjSize;
171   PetscErrorCode     ierr;
172 
173   PetscFunctionBegin;
174   ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
175   ierr = PetscOptionsGetBool(NULL,NULL, "-dm_view_preallocation", &debug, NULL);CHKERRQ(ierr);
176   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
177   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
178   ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
179   ierr = DMGetLocalSection(dm, &section);CHKERRQ(ierr);
180   ierr = DMGetGlobalSection(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, NULL);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, NULL);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, NULL);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 && size > 1) {
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, NULL);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, NULL);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 = PetscArraycpy(&cols[aoff], &rootAdj[roff], rdof);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 
DMPlexUpdateAllocation_Static(DM dm,PetscLayout rLayout,PetscInt bs,PetscInt f,PetscSection sectionAdj,const PetscInt cols[],PetscInt dnz[],PetscInt onz[],PetscInt dnzu[],PetscInt onzu[])584 static PetscErrorCode DMPlexUpdateAllocation_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[])
585 {
586   PetscSection   section;
587   PetscInt       rStart, rEnd, r, pStart, pEnd, p;
588   PetscErrorCode ierr;
589 
590   PetscFunctionBegin;
591   /* This loop needs to change to a loop over points, then field dofs, which means we need to look both sections */
592   ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
593   if (rStart%bs || rEnd%bs) SETERRQ3(PetscObjectComm((PetscObject) rLayout), PETSC_ERR_ARG_WRONG, "Invalid layout [%d, %d) for matrix, must be divisible by block size %d", rStart, rEnd, bs);
594   if (f >= 0 && bs == 1) {
595     ierr = DMGetLocalSection(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 = DMGetGlobalFieldOffset_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 
DMPlexFillMatrix_Static(DM dm,PetscLayout rLayout,PetscInt bs,PetscInt f,PetscSection sectionAdj,const PetscInt cols[],Mat A)645 static PetscErrorCode DMPlexFillMatrix_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], Mat A)
646 {
647   PetscSection   section;
648   PetscScalar   *values;
649   PetscInt       rStart, rEnd, r, pStart, pEnd, p, len, maxRowLen = 0;
650   PetscErrorCode ierr;
651 
652   PetscFunctionBegin;
653   ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
654   for (r = rStart; r < rEnd; ++r) {
655     ierr      = PetscSectionGetDof(sectionAdj, r, &len);CHKERRQ(ierr);
656     maxRowLen = PetscMax(maxRowLen, len);
657   }
658   ierr = PetscCalloc1(maxRowLen, &values);CHKERRQ(ierr);
659   if (f >=0 && bs == 1) {
660     ierr = DMGetLocalSection(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 = DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE);CHKERRQ(ierr);
666       for (r = rS; r < rE; ++r) {
667         PetscInt numCols, cStart;
668 
669         ierr = PetscSectionGetDof(sectionAdj, r, &numCols);CHKERRQ(ierr);
670         ierr = PetscSectionGetOffset(sectionAdj, r, &cStart);CHKERRQ(ierr);
671         ierr = MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);CHKERRQ(ierr);
672       }
673     }
674   } else {
675     for (r = rStart; r < rEnd; ++r) {
676       PetscInt numCols, cStart;
677 
678       ierr = PetscSectionGetDof(sectionAdj, r, &numCols);CHKERRQ(ierr);
679       ierr = PetscSectionGetOffset(sectionAdj, r, &cStart);CHKERRQ(ierr);
680       ierr = MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);CHKERRQ(ierr);
681     }
682   }
683   ierr = PetscFree(values);CHKERRQ(ierr);
684   PetscFunctionReturn(0);
685 }
686 
687 /*@C
688   DMPlexPreallocateOperator - Calculate the matrix nonzero pattern based upon the information in the DM,
689   the PetscDS it contains, and the default PetscSection.
690 
691   Collective
692 
693   Input Arguments:
694 + dm   - The DMPlex
695 . bs   - The matrix blocksize
696 . dnz  - An array to hold the number of nonzeros in the diagonal block
697 . onz  - An array to hold the number of nonzeros in the off-diagonal block
698 . dnzu - An array to hold the number of nonzeros in the upper triangle of the diagonal block
699 . onzu - An array to hold the number of nonzeros in the upper triangle of the off-diagonal block
700 - fillMatrix - If PETSC_TRUE, fill the matrix with zeros
701 
702   Ouput Argument:
703 . A - The preallocated matrix
704 
705   Level: advanced
706 
707 .seealso: DMCreateMatrix()
708 @*/
DMPlexPreallocateOperator(DM dm,PetscInt bs,PetscInt dnz[],PetscInt onz[],PetscInt dnzu[],PetscInt onzu[],Mat A,PetscBool fillMatrix)709 PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
710 {
711   MPI_Comm       comm;
712   PetscDS        prob;
713   MatType        mtype;
714   PetscSF        sf, sfDof;
715   PetscSection   section;
716   PetscInt      *remoteOffsets;
717   PetscSection   sectionAdj[4] = {NULL, NULL, NULL, NULL};
718   PetscInt      *cols[4]       = {NULL, NULL, NULL, NULL};
719   PetscBool      useCone, useClosure;
720   PetscInt       Nf, f, idx, locRows;
721   PetscLayout    rLayout;
722   PetscBool      isSymBlock, isSymSeqBlock, isSymMPIBlock, debug = PETSC_FALSE;
723   PetscMPIInt    size;
724   PetscErrorCode ierr;
725 
726   PetscFunctionBegin;
727   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
728   PetscValidHeaderSpecific(A, MAT_CLASSID, 9);
729   if (dnz)  PetscValidPointer(dnz,5);  if (onz)  PetscValidPointer(onz,6);
730   if (dnzu) PetscValidPointer(dnzu,7); if (onzu) PetscValidPointer(onzu,8);
731   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
732   ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
733   ierr = DMGetLocalSection(dm, &section);CHKERRQ(ierr);
734   ierr = PetscOptionsGetBool(NULL,NULL, "-dm_view_preallocation", &debug, NULL);CHKERRQ(ierr);
735   ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
736   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
737   ierr = PetscLogEventBegin(DMPLEX_Preallocate,dm,0,0,0);CHKERRQ(ierr);
738   /* Create dof SF based on point SF */
739   if (debug) {
740     PetscSection section, sectionGlobal;
741     PetscSF      sf;
742 
743     ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
744     ierr = DMGetLocalSection(dm, &section);CHKERRQ(ierr);
745     ierr = DMGetGlobalSection(dm, &sectionGlobal);CHKERRQ(ierr);
746     ierr = PetscPrintf(comm, "Input Section for Preallocation:\n");CHKERRQ(ierr);
747     ierr = PetscSectionView(section, NULL);CHKERRQ(ierr);
748     ierr = PetscPrintf(comm, "Input Global Section for Preallocation:\n");CHKERRQ(ierr);
749     ierr = PetscSectionView(sectionGlobal, NULL);CHKERRQ(ierr);
750     if (size > 1) {
751       ierr = PetscPrintf(comm, "Input SF for Preallocation:\n");CHKERRQ(ierr);
752       ierr = PetscSFView(sf, NULL);CHKERRQ(ierr);
753     }
754   }
755   ierr = PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets);CHKERRQ(ierr);
756   ierr = PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof);CHKERRQ(ierr);
757   ierr = PetscFree(remoteOffsets);CHKERRQ(ierr);
758   if (debug && size > 1) {
759     ierr = PetscPrintf(comm, "Dof SF for Preallocation:\n");CHKERRQ(ierr);
760     ierr = PetscSFView(sfDof, NULL);CHKERRQ(ierr);
761   }
762   /* Create allocation vectors from adjacency graph */
763   ierr = MatGetLocalSize(A, &locRows, NULL);CHKERRQ(ierr);
764   ierr = PetscLayoutCreate(comm, &rLayout);CHKERRQ(ierr);
765   ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
766   ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
767   ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
768   /* There are 4 types of adjacency */
769   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
770   if (Nf < 1 || bs > 1) {
771     ierr = DMGetBasicAdjacency(dm, &useCone, &useClosure);CHKERRQ(ierr);
772     idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
773     ierr = DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]);CHKERRQ(ierr);
774     ierr = DMPlexUpdateAllocation_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);CHKERRQ(ierr);
775   } else {
776     for (f = 0; f < Nf; ++f) {
777       ierr = DMGetAdjacency(dm, f, &useCone, &useClosure);CHKERRQ(ierr);
778       idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
779       if (!sectionAdj[idx]) {ierr = DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]);CHKERRQ(ierr);}
780       ierr = DMPlexUpdateAllocation_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);CHKERRQ(ierr);
781     }
782   }
783   ierr = PetscSFDestroy(&sfDof);CHKERRQ(ierr);
784   /* Set matrix pattern */
785   ierr = MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);CHKERRQ(ierr);
786   ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
787   /* Check for symmetric storage */
788   ierr = MatGetType(A, &mtype);CHKERRQ(ierr);
789   ierr = PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);CHKERRQ(ierr);
790   ierr = PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);CHKERRQ(ierr);
791   ierr = PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);CHKERRQ(ierr);
792   if (isSymBlock || isSymSeqBlock || isSymMPIBlock) {ierr = MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);CHKERRQ(ierr);}
793   /* Fill matrix with zeros */
794   if (fillMatrix) {
795     if (Nf < 1 || bs > 1) {
796       ierr = DMGetBasicAdjacency(dm, &useCone, &useClosure);CHKERRQ(ierr);
797       idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
798       ierr = DMPlexFillMatrix_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], A);CHKERRQ(ierr);
799     } else {
800       for (f = 0; f < Nf; ++f) {
801         ierr = DMGetAdjacency(dm, f, &useCone, &useClosure);CHKERRQ(ierr);
802         idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
803         ierr = DMPlexFillMatrix_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], A);CHKERRQ(ierr);
804       }
805     }
806     ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
807     ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
808   }
809   ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
810   for (idx = 0; idx < 4; ++idx) {ierr = PetscSectionDestroy(&sectionAdj[idx]);CHKERRQ(ierr); ierr = PetscFree(cols[idx]);CHKERRQ(ierr);}
811   ierr = PetscLogEventEnd(DMPLEX_Preallocate,dm,0,0,0);CHKERRQ(ierr);
812   PetscFunctionReturn(0);
813 }
814 
815 #if 0
816 PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
817 {
818   PetscInt       *tmpClosure,*tmpAdj,*visits;
819   PetscInt        c,cStart,cEnd,pStart,pEnd;
820   PetscErrorCode  ierr;
821 
822   PetscFunctionBegin;
823   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
824   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
825   ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
826 
827   maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1));
828 
829   ierr    = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
830   npoints = pEnd - pStart;
831 
832   ierr = PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits);CHKERRQ(ierr);
833   ierr = PetscArrayzero(lvisits,pEnd-pStart);CHKERRQ(ierr);
834   ierr = PetscArrayzero(visits,pEnd-pStart);CHKERRQ(ierr);
835   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
836   for (c=cStart; c<cEnd; c++) {
837     PetscInt *support = tmpClosure;
838     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);CHKERRQ(ierr);
839     for (p=0; p<supportSize; p++) lvisits[support[p]]++;
840   }
841   ierr = PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
842   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
843   ierr = PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
844   ierr = PetscSFBcastEnd  (sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
845 
846   ierr = PetscSFGetRootRanks();CHKERRQ(ierr);
847 
848 
849   ierr = PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner);CHKERRQ(ierr);
850   for (c=cStart; c<cEnd; c++) {
851     ierr = PetscArrayzero(cellmat,maxClosureSize*maxClosureSize);CHKERRQ(ierr);
852     /*
853      Depth-first walk of transitive closure.
854      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.
855      This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
856      */
857   }
858 
859   ierr = PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);CHKERRQ(ierr);
860   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lonz,onz,MPI_SUM);CHKERRQ(ierr);
861   PetscFunctionReturn(0);
862 }
863 #endif
864