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