xref: /petsc/src/dm/impls/plex/plexpreallocate.c (revision bdf89e91df64bd49f4414ba801fb07e6338da027) 
1 #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 #include <petsc-private/isimpl.h>
3 
4 #undef __FUNCT__
5 #define __FUNCT__ "DMPlexGetAdjacency_Internal"
6 static PetscErrorCode DMPlexGetAdjacency_Internal(DM dm, PetscInt p, PetscBool useClosure, const PetscInt *tmpClosure, PetscInt *adjSize, PetscInt adj[])
7 {
8   const PetscInt *star  = tmpClosure;
9   PetscInt        numAdj = 0, maxAdjSize = *adjSize, starSize, s;
10   PetscErrorCode  ierr;
11 
12   PetscFunctionBegin;
13   ierr = DMPlexGetTransitiveClosure(dm, p, useClosure, &starSize, (PetscInt**) &star);CHKERRQ(ierr);
14   for (s = 2; s < starSize*2; s += 2) {
15     const PetscInt *closure = NULL;
16     PetscInt        closureSize, c, q;
17 
18     ierr = DMPlexGetTransitiveClosure(dm, star[s], (PetscBool)!useClosure, &closureSize, (PetscInt**) &closure);CHKERRQ(ierr);
19     for (c = 0; c < closureSize*2; c += 2) {
20       for (q = 0; q < numAdj || (adj[numAdj++] = closure[c],0); ++q) {
21         if (closure[c] == adj[q]) break;
22       }
23       if (numAdj > maxAdjSize) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%D)", maxAdjSize);
24     }
25     ierr = DMPlexRestoreTransitiveClosure(dm, star[s], (PetscBool)!useClosure, &closureSize, (PetscInt**) &closure);CHKERRQ(ierr);
26   }
27   *adjSize = numAdj;
28   PetscFunctionReturn(0);
29 }
30 
31 #undef __FUNCT__
32 #define __FUNCT__ "DMPlexSetPreallocationCenterDimension"
33 PetscErrorCode DMPlexSetPreallocationCenterDimension(DM dm, PetscInt preallocCenterDim)
34 {
35   DM_Plex *mesh = (DM_Plex*) dm->data;
36 
37   PetscFunctionBegin;
38   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
39   mesh->preallocCenterDim = preallocCenterDim;
40   PetscFunctionReturn(0);
41 }
42 
43 #undef __FUNCT__
44 #define __FUNCT__ "DMPlexPreallocateOperator"
45 PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
46 {
47   DM_Plex           *mesh = (DM_Plex*) dm->data;
48   MPI_Comm           comm;
49   PetscSF            sf, sfDof, sfAdj;
50   PetscSection       leafSectionAdj, rootSectionAdj, sectionAdj;
51   PetscInt           nleaves, l, p;
52   const PetscInt    *leaves;
53   const PetscSFNode *remotes;
54   PetscInt           dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols;
55   PetscInt          *tmpClosure, *tmpAdj, *adj, *rootAdj, *cols, *remoteOffsets;
56   PetscInt           depth, maxConeSize, maxSupportSize, maxClosureSize, maxAdjSize, adjSize;
57   PetscLayout        rLayout;
58   PetscInt           locRows, rStart, rEnd, r;
59   PetscMPIInt        size;
60   PetscBool          useClosure, debug = PETSC_FALSE;
61   PetscErrorCode     ierr;
62 
63   PetscFunctionBegin;
64   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
65   ierr = PetscOptionsGetBool(NULL, "-dm_view_preallocation", &debug, NULL);CHKERRQ(ierr);
66   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
67   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
68   ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
69   /* Create dof SF based on point SF */
70   if (debug) {
71     ierr = PetscPrintf(comm, "Input Section for Preallocation:\n");CHKERRQ(ierr);
72     ierr = PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
73     ierr = PetscPrintf(comm, "Input Global Section for Preallocation:\n");CHKERRQ(ierr);
74     ierr = PetscSectionView(sectionGlobal, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
75     ierr = PetscPrintf(comm, "Input SF for Preallocation:\n");CHKERRQ(ierr);
76     ierr = PetscSFView(sf, NULL);CHKERRQ(ierr);
77   }
78   ierr = PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets);CHKERRQ(ierr);
79   ierr = PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof);CHKERRQ(ierr);
80   if (debug) {
81     ierr = PetscPrintf(comm, "Dof SF for Preallocation:\n");CHKERRQ(ierr);
82     ierr = PetscSFView(sfDof, NULL);CHKERRQ(ierr);
83   }
84   /* Create section for dof adjacency (dof ==> # adj dof) */
85   /*   FEM:   Two points p and q are adjacent if q \in closure(star(p)), preallocCenterDim = dim   */
86   /*   FVM:   Two points p and q are adjacent if q \in star(cone(p)),    preallocCenterDim = dim-1 */
87   /*   FVM++: Two points p and q are adjacent if q \in star(closure(p)), preallocCenterDim = 0     */
88   if (mesh->preallocCenterDim == dim) {
89     useClosure = PETSC_FALSE;
90   } else if (mesh->preallocCenterDim == 0) {
91     useClosure = PETSC_TRUE;
92   } else SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Do not support preallocation with center points of dimension %d", mesh->preallocCenterDim);
93 
94   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
95   ierr = PetscSectionGetStorageSize(section, &numDof);CHKERRQ(ierr);
96   ierr = PetscSectionCreate(comm, &leafSectionAdj);CHKERRQ(ierr);
97   ierr = PetscSectionSetChart(leafSectionAdj, 0, numDof);CHKERRQ(ierr);
98   ierr = PetscSectionCreate(comm, &rootSectionAdj);CHKERRQ(ierr);
99   ierr = PetscSectionSetChart(rootSectionAdj, 0, numDof);CHKERRQ(ierr);
100   /*   Fill in the ghost dofs on the interface */
101   ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes);CHKERRQ(ierr);
102   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
103   ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
104 
105   maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth),PetscPowInt(mesh->maxSupportSize,depth)) + 2;
106   maxAdjSize     = PetscPowInt(mesh->maxConeSize,depth) * PetscPowInt(mesh->maxSupportSize,depth) + 1;
107 
108   ierr = PetscMalloc2(maxClosureSize,PetscInt,&tmpClosure,maxAdjSize,PetscInt,&tmpAdj);CHKERRQ(ierr);
109 
110   /*
111    ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point.
112     1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj
113        Reduce those counts to rootSectionAdj (now redundantly counting some interface points)
114     2. Visit owned points on interface, count adjacencies placing in rootSectionAdj
115        Create sfAdj connecting rootSectionAdj and leafSectionAdj
116     3. Visit unowned points on interface, write adjacencies to adj
117        Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies)
118     4. Visit owned points on interface, write adjacencies to rootAdj
119        Remove redundancy in rootAdj
120    ** The last two traversals use transitive closure
121     5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj)
122        Allocate memory addressed by sectionAdj (cols)
123     6. Visit all owned points in the subdomain, insert dof adjacencies into cols
124    ** Knowing all the column adjacencies, check ownership and sum into dnz and onz
125   */
126 
127   for (l = 0; l < nleaves; ++l) {
128     PetscInt dof, off, d, q;
129     PetscInt p = leaves[l], numAdj = maxAdjSize;
130 
131     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
132     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
133     ierr = DMPlexGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
134     for (q = 0; q < numAdj; ++q) {
135       const PetscInt padj = tmpAdj[q];
136       PetscInt ndof, ncdof;
137 
138       if ((padj < pStart) || (padj >= pEnd)) continue;
139       ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
140       ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
141       for (d = off; d < off+dof; ++d) {
142         ierr = PetscSectionAddDof(leafSectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
143       }
144     }
145   }
146   ierr = PetscSectionSetUp(leafSectionAdj);CHKERRQ(ierr);
147   if (debug) {
148     ierr = PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n");CHKERRQ(ierr);
149     ierr = PetscSectionView(leafSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
150   }
151   /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */
152   if (size > 1) {
153     ierr = PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);CHKERRQ(ierr);
154     ierr = PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);CHKERRQ(ierr);
155   }
156   if (debug) {
157     ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots:\n");CHKERRQ(ierr);
158     ierr = PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
159   }
160   /* Add in local adjacency sizes for owned dofs on interface (roots) */
161   for (p = pStart; p < pEnd; ++p) {
162     PetscInt numAdj = maxAdjSize, adof, dof, off, d, q;
163 
164     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
165     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
166     if (!dof) continue;
167     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
168     if (adof <= 0) continue;
169     ierr = DMPlexGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
170     for (q = 0; q < numAdj; ++q) {
171       const PetscInt padj = tmpAdj[q];
172       PetscInt ndof, ncdof;
173 
174       if ((padj < pStart) || (padj >= pEnd)) continue;
175       ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
176       ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
177       for (d = off; d < off+dof; ++d) {
178         ierr = PetscSectionAddDof(rootSectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
179       }
180     }
181   }
182   ierr = PetscSectionSetUp(rootSectionAdj);CHKERRQ(ierr);
183   if (debug) {
184     ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after local additions:\n");CHKERRQ(ierr);
185     ierr = PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
186   }
187   /* Create adj SF based on dof SF */
188   ierr = PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets);CHKERRQ(ierr);
189   ierr = PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj);CHKERRQ(ierr);
190   if (debug) {
191     ierr = PetscPrintf(comm, "Adjacency SF for Preallocation:\n");CHKERRQ(ierr);
192     ierr = PetscSFView(sfAdj, NULL);CHKERRQ(ierr);
193   }
194   ierr = PetscSFDestroy(&sfDof);CHKERRQ(ierr);
195   /* Create leaf adjacency */
196   ierr = PetscSectionSetUp(leafSectionAdj);CHKERRQ(ierr);
197   ierr = PetscSectionGetStorageSize(leafSectionAdj, &adjSize);CHKERRQ(ierr);
198   ierr = PetscMalloc(adjSize * sizeof(PetscInt), &adj);CHKERRQ(ierr);
199   ierr = PetscMemzero(adj, adjSize * sizeof(PetscInt));CHKERRQ(ierr);
200   for (l = 0; l < nleaves; ++l) {
201     PetscInt dof, off, d, q;
202     PetscInt p = leaves[l], numAdj = maxAdjSize;
203 
204     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
205     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
206     ierr = DMPlexGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
207     for (d = off; d < off+dof; ++d) {
208       PetscInt aoff, i = 0;
209 
210       ierr = PetscSectionGetOffset(leafSectionAdj, d, &aoff);CHKERRQ(ierr);
211       for (q = 0; q < numAdj; ++q) {
212         const PetscInt padj = tmpAdj[q];
213         PetscInt ndof, ncdof, ngoff, nd;
214 
215         if ((padj < pStart) || (padj >= pEnd)) continue;
216         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
217         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
218         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
219         for (nd = 0; nd < ndof-ncdof; ++nd) {
220           adj[aoff+i] = (ngoff < 0 ? -(ngoff+1) : ngoff) + nd;
221           ++i;
222         }
223       }
224     }
225   }
226   /* Debugging */
227   if (debug) {
228     IS tmp;
229     ierr = PetscPrintf(comm, "Leaf adjacency indices\n");CHKERRQ(ierr);
230     ierr = ISCreateGeneral(comm, adjSize, adj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
231     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
232   }
233   /* Gather adjacenct indices to root */
234   ierr = PetscSectionGetStorageSize(rootSectionAdj, &adjSize);CHKERRQ(ierr);
235   ierr = PetscMalloc(adjSize * sizeof(PetscInt), &rootAdj);CHKERRQ(ierr);
236   for (r = 0; r < adjSize; ++r) rootAdj[r] = -1;
237   if (size > 1) {
238     ierr = PetscSFGatherBegin(sfAdj, MPIU_INT, adj, rootAdj);CHKERRQ(ierr);
239     ierr = PetscSFGatherEnd(sfAdj, MPIU_INT, adj, rootAdj);CHKERRQ(ierr);
240   }
241   ierr = PetscSFDestroy(&sfAdj);CHKERRQ(ierr);
242   ierr = PetscFree(adj);CHKERRQ(ierr);
243   /* Debugging */
244   if (debug) {
245     IS tmp;
246     ierr = PetscPrintf(comm, "Root adjacency indices after gather\n");CHKERRQ(ierr);
247     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
248     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
249   }
250   /* Add in local adjacency indices for owned dofs on interface (roots) */
251   for (p = pStart; p < pEnd; ++p) {
252     PetscInt numAdj = maxAdjSize, adof, dof, off, d, q;
253 
254     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
255     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
256     if (!dof) continue;
257     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
258     if (adof <= 0) continue;
259     ierr = DMPlexGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
260     for (d = off; d < off+dof; ++d) {
261       PetscInt adof, aoff, i;
262 
263       ierr = PetscSectionGetDof(rootSectionAdj, d, &adof);CHKERRQ(ierr);
264       ierr = PetscSectionGetOffset(rootSectionAdj, d, &aoff);CHKERRQ(ierr);
265       i    = adof-1;
266       for (q = 0; q < numAdj; ++q) {
267         const PetscInt padj = tmpAdj[q];
268         PetscInt ndof, ncdof, ngoff, nd;
269 
270         if ((padj < pStart) || (padj >= pEnd)) continue;
271         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
272         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
273         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
274         for (nd = 0; nd < ndof-ncdof; ++nd) {
275           rootAdj[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
276           --i;
277         }
278       }
279     }
280   }
281   /* Debugging */
282   if (debug) {
283     IS tmp;
284     ierr = PetscPrintf(comm, "Root adjacency indices\n");CHKERRQ(ierr);
285     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
286     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
287   }
288   /* Compress indices */
289   ierr = PetscSectionSetUp(rootSectionAdj);CHKERRQ(ierr);
290   for (p = pStart; p < pEnd; ++p) {
291     PetscInt dof, cdof, off, d;
292     PetscInt adof, aoff;
293 
294     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
295     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
296     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
297     if (!dof) continue;
298     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
299     if (adof <= 0) continue;
300     for (d = off; d < off+dof-cdof; ++d) {
301       ierr = PetscSectionGetDof(rootSectionAdj, d, &adof);CHKERRQ(ierr);
302       ierr = PetscSectionGetOffset(rootSectionAdj, d, &aoff);CHKERRQ(ierr);
303       ierr = PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]);CHKERRQ(ierr);
304       ierr = PetscSectionSetDof(rootSectionAdj, d, adof);CHKERRQ(ierr);
305     }
306   }
307   /* Debugging */
308   if (debug) {
309     IS tmp;
310     ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after compression:\n");CHKERRQ(ierr);
311     ierr = PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
312     ierr = PetscPrintf(comm, "Root adjacency indices after compression\n");CHKERRQ(ierr);
313     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
314     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
315   }
316   /* Build adjacency section: Maps global indices to sets of adjacent global indices */
317   ierr = PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd);CHKERRQ(ierr);
318   ierr = PetscSectionCreate(comm, &sectionAdj);CHKERRQ(ierr);
319   ierr = PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd);CHKERRQ(ierr);
320   for (p = pStart; p < pEnd; ++p) {
321     PetscInt  numAdj = maxAdjSize, dof, cdof, off, goff, d, q;
322     PetscBool found  = PETSC_TRUE;
323 
324     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
325     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
326     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
327     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
328     for (d = 0; d < dof-cdof; ++d) {
329       PetscInt ldof, rdof;
330 
331       ierr = PetscSectionGetDof(leafSectionAdj, off+d, &ldof);CHKERRQ(ierr);
332       ierr = PetscSectionGetDof(rootSectionAdj, off+d, &rdof);CHKERRQ(ierr);
333       if (ldof > 0) {
334         /* We do not own this point */
335       } else if (rdof > 0) {
336         ierr = PetscSectionSetDof(sectionAdj, goff+d, rdof);CHKERRQ(ierr);
337       } else {
338         found = PETSC_FALSE;
339       }
340     }
341     if (found) continue;
342     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
343     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
344     ierr = DMPlexGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
345     for (q = 0; q < numAdj; ++q) {
346       const PetscInt padj = tmpAdj[q];
347       PetscInt ndof, ncdof, noff;
348 
349       if ((padj < pStart) || (padj >= pEnd)) continue;
350       ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
351       ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
352       ierr = PetscSectionGetOffset(section, padj, &noff);CHKERRQ(ierr);
353       for (d = goff; d < goff+dof-cdof; ++d) {
354         ierr = PetscSectionAddDof(sectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
355       }
356     }
357   }
358   ierr = PetscSectionSetUp(sectionAdj);CHKERRQ(ierr);
359   if (debug) {
360     ierr = PetscPrintf(comm, "Adjacency Section for Preallocation:\n");CHKERRQ(ierr);
361     ierr = PetscSectionView(sectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
362   }
363   /* Get adjacent indices */
364   ierr = PetscSectionGetStorageSize(sectionAdj, &numCols);CHKERRQ(ierr);
365   ierr = PetscMalloc(numCols * sizeof(PetscInt), &cols);CHKERRQ(ierr);
366   for (p = pStart; p < pEnd; ++p) {
367     PetscInt  numAdj = maxAdjSize, dof, cdof, off, goff, d, q;
368     PetscBool found  = PETSC_TRUE;
369 
370     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
371     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
372     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
373     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
374     for (d = 0; d < dof-cdof; ++d) {
375       PetscInt ldof, rdof;
376 
377       ierr = PetscSectionGetDof(leafSectionAdj, off+d, &ldof);CHKERRQ(ierr);
378       ierr = PetscSectionGetDof(rootSectionAdj, off+d, &rdof);CHKERRQ(ierr);
379       if (ldof > 0) {
380         /* We do not own this point */
381       } else if (rdof > 0) {
382         PetscInt aoff, roff;
383 
384         ierr = PetscSectionGetOffset(sectionAdj, goff+d, &aoff);CHKERRQ(ierr);
385         ierr = PetscSectionGetOffset(rootSectionAdj, off+d, &roff);CHKERRQ(ierr);
386         ierr = PetscMemcpy(&cols[aoff], &rootAdj[roff], rdof * sizeof(PetscInt));CHKERRQ(ierr);
387       } else {
388         found = PETSC_FALSE;
389       }
390     }
391     if (found) continue;
392     ierr = DMPlexGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
393     for (d = goff; d < goff+dof-cdof; ++d) {
394       PetscInt adof, aoff, i = 0;
395 
396       ierr = PetscSectionGetDof(sectionAdj, d, &adof);CHKERRQ(ierr);
397       ierr = PetscSectionGetOffset(sectionAdj, d, &aoff);CHKERRQ(ierr);
398       for (q = 0; q < numAdj; ++q) {
399         const PetscInt  padj = tmpAdj[q];
400         PetscInt        ndof, ncdof, ngoff, nd;
401         const PetscInt *ncind;
402 
403         /* Adjacent points may not be in the section chart */
404         if ((padj < pStart) || (padj >= pEnd)) continue;
405         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
406         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
407         ierr = PetscSectionGetConstraintIndices(section, padj, &ncind);CHKERRQ(ierr);
408         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
409         for (nd = 0; nd < ndof-ncdof; ++nd, ++i) {
410           cols[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
411         }
412       }
413       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);
414     }
415   }
416   ierr = PetscSectionDestroy(&leafSectionAdj);CHKERRQ(ierr);
417   ierr = PetscSectionDestroy(&rootSectionAdj);CHKERRQ(ierr);
418   ierr = PetscFree(rootAdj);CHKERRQ(ierr);
419   ierr = PetscFree2(tmpClosure, tmpAdj);CHKERRQ(ierr);
420   /* Debugging */
421   if (debug) {
422     IS tmp;
423     ierr = PetscPrintf(comm, "Column indices\n");CHKERRQ(ierr);
424     ierr = ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
425     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
426   }
427   /* Create allocation vectors from adjacency graph */
428   ierr = MatGetLocalSize(A, &locRows, NULL);CHKERRQ(ierr);
429   ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)A), &rLayout);CHKERRQ(ierr);
430   ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
431   ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
432   ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
433   ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
434   ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
435   /* Only loop over blocks of rows */
436   if (rStart%bs || rEnd%bs) SETERRQ3(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid layout [%d, %d) for matrix, must be divisible by block size %d", rStart, rEnd, bs);
437   for (r = rStart/bs; r < rEnd/bs; ++r) {
438     const PetscInt row = r*bs;
439     PetscInt       numCols, cStart, c;
440 
441     ierr = PetscSectionGetDof(sectionAdj, row, &numCols);CHKERRQ(ierr);
442     ierr = PetscSectionGetOffset(sectionAdj, row, &cStart);CHKERRQ(ierr);
443     for (c = cStart; c < cStart+numCols; ++c) {
444       if ((cols[c] >= rStart*bs) && (cols[c] < rEnd*bs)) {
445         ++dnz[r-rStart];
446         if (cols[c] >= row) ++dnzu[r-rStart];
447       } else {
448         ++onz[r-rStart];
449         if (cols[c] >= row) ++onzu[r-rStart];
450       }
451     }
452   }
453   if (bs > 1) {
454     for (r = 0; r < locRows/bs; ++r) {
455       dnz[r]  /= bs;
456       onz[r]  /= bs;
457       dnzu[r] /= bs;
458       onzu[r] /= bs;
459     }
460   }
461   /* Set matrix pattern */
462   ierr = MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);CHKERRQ(ierr);
463   ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
464   /* Fill matrix with zeros */
465   if (fillMatrix) {
466     PetscScalar *values;
467     PetscInt     maxRowLen = 0;
468 
469     for (r = rStart; r < rEnd; ++r) {
470       PetscInt len;
471 
472       ierr      = PetscSectionGetDof(sectionAdj, r, &len);CHKERRQ(ierr);
473       maxRowLen = PetscMax(maxRowLen, len);
474     }
475     ierr = PetscMalloc(maxRowLen * sizeof(PetscScalar), &values);CHKERRQ(ierr);
476     ierr = PetscMemzero(values, maxRowLen * sizeof(PetscScalar));CHKERRQ(ierr);
477     for (r = rStart; r < rEnd; ++r) {
478       PetscInt numCols, cStart;
479 
480       ierr = PetscSectionGetDof(sectionAdj, r, &numCols);CHKERRQ(ierr);
481       ierr = PetscSectionGetOffset(sectionAdj, r, &cStart);CHKERRQ(ierr);
482       ierr = MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);CHKERRQ(ierr);
483     }
484     ierr = PetscFree(values);CHKERRQ(ierr);
485     ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
486     ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
487   }
488   ierr = PetscSectionDestroy(&sectionAdj);CHKERRQ(ierr);
489   ierr = PetscFree(cols);CHKERRQ(ierr);
490   PetscFunctionReturn(0);
491 }
492 
493 #if 0
494 #undef __FUNCT__
495 #define __FUNCT__ "DMPlexPreallocateOperator_2"
496 PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
497 {
498   PetscInt       *tmpClosure,*tmpAdj,*visits;
499   PetscInt        c,cStart,cEnd,pStart,pEnd;
500   PetscErrorCode  ierr;
501 
502   PetscFunctionBegin;
503   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
504   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
505   ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
506 
507   maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth),PetscPowInt(mesh->maxSupportSize,depth));
508 
509   ierr    = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
510   npoints = pEnd - pStart;
511 
512   ierr = PetscMalloc3(maxClosureSize,PetscInt,&tmpClosure,npoints,PetscInt,&lvisits,npoints,PetscInt,&visits);CHKERRQ(ierr);
513   ierr = PetscMemzero(lvisits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
514   ierr = PetscMemzero(visits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
515   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
516   for (c=cStart; c<cEnd; c++) {
517     PetscInt *support = tmpClosure;
518     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);CHKERRQ(ierr);
519     for (p=0; p<supportSize; p++) lvisits[support[p]]++;
520   }
521   ierr = PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
522   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
523   ierr = PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
524   ierr = PetscSFBcastEnd  (sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
525 
526   ierr = PetscSFGetRanks();CHKERRQ(ierr);
527 
528 
529   ierr = PetscMalloc2(maxClosureSize*maxClosureSize,PetscInt,&cellmat,npoints,PetscInt,&owner);CHKERRQ(ierr);
530   for (c=cStart; c<cEnd; c++) {
531     ierr = PetscMemzero(cellmat,maxClosureSize*maxClosureSize*sizeof(PetscInt));CHKERRQ(ierr);
532     /*
533      Depth-first walk of transitive closure.
534      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.
535      This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
536      */
537   }
538 
539   ierr = PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);CHKERRQ(ierr);
540   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lonz,onz,MPI_SUM);CHKERRQ(ierr);
541   PetscFunctionReturn(0);
542 }
543 #endif
544