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, §ion);CHKERRQ(ierr);
16 ierr = DMGetGlobalSection(dm, §ionGlobal);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, §ion);CHKERRQ(ierr);
180 ierr = DMGetGlobalSection(dm, §ionGlobal);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, §ionAdj);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, §ion);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, §ion);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, §ion);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, §ion);CHKERRQ(ierr);
745 ierr = DMGetGlobalSection(dm, §ionGlobal);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, §ionAdj[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, §ionAdj[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(§ionAdj[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