xref: /petsc/src/vec/is/sf/utils/sfutils.c (revision 2cb5e1cc91ad4e0472b8976614576d28ebef7100)
1 #include <petsc/private/sfimpl.h>       /*I  "petscsf.h"   I*/
2 #include <petsc/private/sectionimpl.h>
3 
4 /*@C
5    PetscSFSetGraphLayout - Set a parallel star forest via global indices and a PetscLayout
6 
7    Collective
8 
9    Input Arguments:
10 +  sf - star forest
11 .  layout - PetscLayout defining the global space
12 .  nleaves - number of leaf vertices on the current process, each of these references a root on any process
13 .  ilocal - locations of leaves in leafdata buffers, pass NULL for contiguous storage
14 .  localmode - copy mode for ilocal
15 -  iremote - remote locations of root vertices for each leaf on the current process
16 
17    Level: intermediate
18 
19    Developers Note: Local indices which are the identity permutation in the range [0,nleaves) are discarded as they
20    encode contiguous storage. In such case, if localmode is PETSC_OWN_POINTER, the memory is deallocated as it is not
21    needed
22 
23 .seealso: PetscSFCreate(), PetscSFView(), PetscSFSetGraph(), PetscSFGetGraph()
24 @*/
25 PetscErrorCode PetscSFSetGraphLayout(PetscSF sf,PetscLayout layout,PetscInt nleaves,const PetscInt *ilocal,PetscCopyMode localmode,const PetscInt *iremote)
26 {
27   PetscErrorCode ierr;
28   PetscInt       i,nroots;
29   PetscSFNode    *remote;
30 
31   PetscFunctionBegin;
32   ierr = PetscLayoutGetLocalSize(layout,&nroots);CHKERRQ(ierr);
33   ierr = PetscMalloc1(nleaves,&remote);CHKERRQ(ierr);
34   for (i=0; i<nleaves; i++) {
35     PetscMPIInt owner = -1;
36     ierr = PetscLayoutFindOwner(layout,iremote[i],&owner);CHKERRQ(ierr);
37     remote[i].rank  = owner;
38     remote[i].index = iremote[i] - layout->range[owner];
39   }
40   ierr = PetscSFSetGraph(sf,nroots,nleaves,ilocal,localmode,remote,PETSC_OWN_POINTER);CHKERRQ(ierr);
41   PetscFunctionReturn(0);
42 }
43 
44 /*@
45   PetscSFSetGraphSection - Sets the PetscSF graph encoding the parallel dof overlap based upon the PetscSections
46   describing the data layout.
47 
48   Input Parameters:
49 + sf - The SF
50 . localSection - PetscSection describing the local data layout
51 - globalSection - PetscSection describing the global data layout
52 
53   Level: developer
54 
55 .seealso: PetscSFSetGraph(), PetscSFSetGraphLayout()
56 @*/
57 PetscErrorCode PetscSFSetGraphSection(PetscSF sf, PetscSection localSection, PetscSection globalSection)
58 {
59   MPI_Comm       comm;
60   PetscLayout    layout;
61   const PetscInt *ranges;
62   PetscInt       *local;
63   PetscSFNode    *remote;
64   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
65   PetscMPIInt    size, rank;
66   PetscErrorCode ierr;
67 
68   PetscFunctionBegin;
69   PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 1);
70   PetscValidHeaderSpecific(localSection, PETSC_SECTION_CLASSID, 2);
71   PetscValidHeaderSpecific(globalSection, PETSC_SECTION_CLASSID, 3);
72 
73   ierr = PetscObjectGetComm((PetscObject)sf,&comm);CHKERRQ(ierr);
74   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
75   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
76   ierr = PetscSectionGetChart(globalSection, &pStart, &pEnd);CHKERRQ(ierr);
77   ierr = PetscSectionGetConstrainedStorageSize(globalSection, &nroots);CHKERRQ(ierr);
78   ierr = PetscLayoutCreate(comm, &layout);CHKERRQ(ierr);
79   ierr = PetscLayoutSetBlockSize(layout, 1);CHKERRQ(ierr);
80   ierr = PetscLayoutSetLocalSize(layout, nroots);CHKERRQ(ierr);
81   ierr = PetscLayoutSetUp(layout);CHKERRQ(ierr);
82   ierr = PetscLayoutGetRanges(layout, &ranges);CHKERRQ(ierr);
83   for (p = pStart; p < pEnd; ++p) {
84     PetscInt gdof, gcdof;
85 
86     ierr = PetscSectionGetDof(globalSection, p, &gdof);CHKERRQ(ierr);
87     ierr = PetscSectionGetConstraintDof(globalSection, p, &gcdof);CHKERRQ(ierr);
88     if (gcdof > (gdof < 0 ? -(gdof+1) : gdof)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %D has %D constraints > %D dof", p, gcdof, (gdof < 0 ? -(gdof+1) : gdof));
89     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
90   }
91   ierr = PetscMalloc1(nleaves, &local);CHKERRQ(ierr);
92   ierr = PetscMalloc1(nleaves, &remote);CHKERRQ(ierr);
93   for (p = pStart, l = 0; p < pEnd; ++p) {
94     const PetscInt *cind;
95     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;
96 
97     ierr = PetscSectionGetDof(localSection, p, &dof);CHKERRQ(ierr);
98     ierr = PetscSectionGetOffset(localSection, p, &off);CHKERRQ(ierr);
99     ierr = PetscSectionGetConstraintDof(localSection, p, &cdof);CHKERRQ(ierr);
100     ierr = PetscSectionGetConstraintIndices(localSection, p, &cind);CHKERRQ(ierr);
101     ierr = PetscSectionGetDof(globalSection, p, &gdof);CHKERRQ(ierr);
102     ierr = PetscSectionGetConstraintDof(globalSection, p, &gcdof);CHKERRQ(ierr);
103     ierr = PetscSectionGetOffset(globalSection, p, &goff);CHKERRQ(ierr);
104     if (!gdof) continue; /* Censored point */
105     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
106     if (gsize != dof-cdof) {
107       if (gsize != dof) SETERRQ4(comm, PETSC_ERR_ARG_WRONG, "Global dof %D for point %D is neither the constrained size %D, nor the unconstrained %D", gsize, p, dof-cdof, dof);
108       cdof = 0; /* Ignore constraints */
109     }
110     for (d = 0, c = 0; d < dof; ++d) {
111       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
112       local[l+d-c] = off+d;
113     }
114     if (gdof < 0) {
115       for (d = 0; d < gsize; ++d, ++l) {
116         PetscInt offset = -(goff+1) + d, r;
117 
118         ierr = PetscFindInt(offset,size+1,ranges,&r);CHKERRQ(ierr);
119         if (r < 0) r = -(r+2);
120         if ((r < 0) || (r >= size)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %D mapped to invalid process %D (%D, %D)", p, r, gdof, goff);
121         remote[l].rank  = r;
122         remote[l].index = offset - ranges[r];
123       }
124     } else {
125       for (d = 0; d < gsize; ++d, ++l) {
126         remote[l].rank  = rank;
127         remote[l].index = goff+d - ranges[rank];
128       }
129     }
130   }
131   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %D != nleaves %D", l, nleaves);
132   ierr = PetscLayoutDestroy(&layout);CHKERRQ(ierr);
133   ierr = PetscSFSetGraph(sf, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);CHKERRQ(ierr);
134   PetscFunctionReturn(0);
135 }
136 
137 static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
138 {
139   PetscErrorCode ierr;
140 
141   PetscFunctionBegin;
142   if (!s->bc) {
143     ierr = PetscSectionCreate(PETSC_COMM_SELF, &s->bc);CHKERRQ(ierr);
144     ierr = PetscSectionSetChart(s->bc, s->pStart, s->pEnd);CHKERRQ(ierr);
145   }
146   PetscFunctionReturn(0);
147 }
148 
149 /*@C
150   PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF
151 
152   Collective on sf
153 
154   Input Parameters:
155 + sf - The SF
156 - rootSection - Section defined on root space
157 
158   Output Parameters:
159 + remoteOffsets - root offsets in leaf storage, or NULL
160 - leafSection - Section defined on the leaf space
161 
162   Level: advanced
163 
164 .seealso: PetscSFCreate()
165 @*/
166 PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
167 {
168   PetscSF        embedSF;
169   const PetscInt *indices;
170   IS             selected;
171   PetscInt       numFields, nroots, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, f, c;
172   PetscBool      *sub, hasc;
173   PetscErrorCode ierr;
174 
175   PetscFunctionBegin;
176   ierr = PetscLogEventBegin(PETSCSF_DistSect,sf,0,0,0);CHKERRQ(ierr);
177   ierr = PetscSectionGetNumFields(rootSection, &numFields);CHKERRQ(ierr);
178   if (numFields) {ierr = PetscSectionSetNumFields(leafSection, numFields);CHKERRQ(ierr);}
179   ierr = PetscMalloc1(numFields+2, &sub);CHKERRQ(ierr);
180   sub[1] = rootSection->bc ? PETSC_TRUE : PETSC_FALSE;
181   for (f = 0; f < numFields; ++f) {
182     PetscSectionSym sym;
183     const char      *name   = NULL;
184     PetscInt        numComp = 0;
185 
186     sub[2 + f] = rootSection->field[f]->bc ? PETSC_TRUE : PETSC_FALSE;
187     ierr = PetscSectionGetFieldComponents(rootSection, f, &numComp);CHKERRQ(ierr);
188     ierr = PetscSectionGetFieldName(rootSection, f, &name);CHKERRQ(ierr);
189     ierr = PetscSectionGetFieldSym(rootSection, f, &sym);CHKERRQ(ierr);
190     ierr = PetscSectionSetFieldComponents(leafSection, f, numComp);CHKERRQ(ierr);
191     ierr = PetscSectionSetFieldName(leafSection, f, name);CHKERRQ(ierr);
192     ierr = PetscSectionSetFieldSym(leafSection, f, sym);CHKERRQ(ierr);
193     for (c = 0; c < rootSection->numFieldComponents[f]; ++c) {
194       ierr = PetscSectionGetComponentName(rootSection, f, c, &name);CHKERRQ(ierr);
195       ierr = PetscSectionSetComponentName(leafSection, f, c, name);CHKERRQ(ierr);
196     }
197   }
198   ierr = PetscSectionGetChart(rootSection, &rpStart, &rpEnd);CHKERRQ(ierr);
199   ierr = PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);CHKERRQ(ierr);
200   rpEnd = PetscMin(rpEnd,nroots);
201   rpEnd = PetscMax(rpStart,rpEnd);
202   /* see if we can avoid creating the embedded SF, since it can cost more than an allreduce */
203   sub[0] = (PetscBool)(nroots != rpEnd - rpStart);
204   ierr = MPIU_Allreduce(MPI_IN_PLACE, sub, 2+numFields, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)sf));CHKERRQ(ierr);
205   if (sub[0]) {
206     ierr = ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);CHKERRQ(ierr);
207     ierr = ISGetIndices(selected, &indices);CHKERRQ(ierr);
208     ierr = PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);CHKERRQ(ierr);
209     ierr = ISRestoreIndices(selected, &indices);CHKERRQ(ierr);
210     ierr = ISDestroy(&selected);CHKERRQ(ierr);
211   } else {
212     ierr = PetscObjectReference((PetscObject)sf);CHKERRQ(ierr);
213     embedSF = sf;
214   }
215   ierr = PetscSFGetLeafRange(embedSF, &lpStart, &lpEnd);CHKERRQ(ierr);
216   lpEnd++;
217 
218   ierr = PetscSectionSetChart(leafSection, lpStart, lpEnd);CHKERRQ(ierr);
219 
220   /* Constrained dof section */
221   hasc = sub[1];
222   for (f = 0; f < numFields; ++f) hasc = (PetscBool)(hasc || sub[2+f]);
223 
224   /* Could fuse these at the cost of copies and extra allocation */
225   ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);CHKERRQ(ierr);
226   ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);CHKERRQ(ierr);
227   if (sub[1]) {
228     ierr = PetscSectionCheckConstraints_Static(rootSection);CHKERRQ(ierr);
229     ierr = PetscSectionCheckConstraints_Static(leafSection);CHKERRQ(ierr);
230     ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart]);CHKERRQ(ierr);
231     ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart]);CHKERRQ(ierr);
232   }
233   for (f = 0; f < numFields; ++f) {
234     ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);CHKERRQ(ierr);
235     ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);CHKERRQ(ierr);
236     if (sub[2+f]) {
237       ierr = PetscSectionCheckConstraints_Static(rootSection->field[f]);CHKERRQ(ierr);
238       ierr = PetscSectionCheckConstraints_Static(leafSection->field[f]);CHKERRQ(ierr);
239       ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart]);CHKERRQ(ierr);
240       ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart]);CHKERRQ(ierr);
241     }
242   }
243   if (remoteOffsets) {
244     ierr = PetscMalloc1(lpEnd - lpStart, remoteOffsets);CHKERRQ(ierr);
245     ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);CHKERRQ(ierr);
246     ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);CHKERRQ(ierr);
247   }
248   ierr = PetscSectionSetUp(leafSection);CHKERRQ(ierr);
249   if (hasc) { /* need to communicate bcIndices */
250     PetscSF  bcSF;
251     PetscInt *rOffBc;
252 
253     ierr = PetscMalloc1(lpEnd - lpStart, &rOffBc);CHKERRQ(ierr);
254     if (sub[1]) {
255       ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);CHKERRQ(ierr);
256       ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);CHKERRQ(ierr);
257       ierr = PetscSFCreateSectionSF(embedSF, rootSection->bc, rOffBc, leafSection->bc, &bcSF);CHKERRQ(ierr);
258       ierr = PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices);CHKERRQ(ierr);
259       ierr = PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices);CHKERRQ(ierr);
260       ierr = PetscSFDestroy(&bcSF);CHKERRQ(ierr);
261     }
262     for (f = 0; f < numFields; ++f) {
263       if (sub[2+f]) {
264         ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);CHKERRQ(ierr);
265         ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);CHKERRQ(ierr);
266         ierr = PetscSFCreateSectionSF(embedSF, rootSection->field[f]->bc, rOffBc, leafSection->field[f]->bc, &bcSF);CHKERRQ(ierr);
267         ierr = PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices);CHKERRQ(ierr);
268         ierr = PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices);CHKERRQ(ierr);
269         ierr = PetscSFDestroy(&bcSF);CHKERRQ(ierr);
270       }
271     }
272     ierr = PetscFree(rOffBc);CHKERRQ(ierr);
273   }
274   ierr = PetscSFDestroy(&embedSF);CHKERRQ(ierr);
275   ierr = PetscFree(sub);CHKERRQ(ierr);
276   ierr = PetscLogEventEnd(PETSCSF_DistSect,sf,0,0,0);CHKERRQ(ierr);
277   PetscFunctionReturn(0);
278 }
279 
280 /*@C
281   PetscSFCreateRemoteOffsets - Create offsets for point data on remote processes
282 
283   Collective on sf
284 
285   Input Parameters:
286 + sf          - The SF
287 . rootSection - Data layout of remote points for outgoing data (this is layout for SF roots)
288 - leafSection - Data layout of local points for incoming data  (this is layout for SF leaves)
289 
290   Output Parameter:
291 . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
292 
293   Level: developer
294 
295 .seealso: PetscSFCreate()
296 @*/
297 PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
298 {
299   PetscSF         embedSF;
300   const PetscInt *indices;
301   IS              selected;
302   PetscInt        numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
303   PetscErrorCode  ierr;
304 
305   PetscFunctionBegin;
306   *remoteOffsets = NULL;
307   ierr = PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);CHKERRQ(ierr);
308   if (numRoots < 0) PetscFunctionReturn(0);
309   ierr = PetscLogEventBegin(PETSCSF_RemoteOff,sf,0,0,0);CHKERRQ(ierr);
310   ierr = PetscSectionGetChart(rootSection, &rpStart, &rpEnd);CHKERRQ(ierr);
311   ierr = PetscSectionGetChart(leafSection, &lpStart, &lpEnd);CHKERRQ(ierr);
312   ierr = ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);CHKERRQ(ierr);
313   ierr = ISGetIndices(selected, &indices);CHKERRQ(ierr);
314   ierr = PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);CHKERRQ(ierr);
315   ierr = ISRestoreIndices(selected, &indices);CHKERRQ(ierr);
316   ierr = ISDestroy(&selected);CHKERRQ(ierr);
317   ierr = PetscCalloc1(lpEnd - lpStart, remoteOffsets);CHKERRQ(ierr);
318   ierr = PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);CHKERRQ(ierr);
319   ierr = PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);CHKERRQ(ierr);
320   ierr = PetscSFDestroy(&embedSF);CHKERRQ(ierr);
321   ierr = PetscLogEventEnd(PETSCSF_RemoteOff,sf,0,0,0);CHKERRQ(ierr);
322   PetscFunctionReturn(0);
323 }
324 
325 /*@C
326   PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points
327 
328   Collective on sf
329 
330   Input Parameters:
331 + sf - The SF
332 . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
333 . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
334 - leafSection - Data layout of local points for incoming data  (this is the distributed section)
335 
336   Output Parameters:
337 - sectionSF - The new SF
338 
339   Note: Either rootSection or remoteOffsets can be specified
340 
341   Level: advanced
342 
343 .seealso: PetscSFCreate()
344 @*/
345 PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
346 {
347   MPI_Comm          comm;
348   const PetscInt    *localPoints;
349   const PetscSFNode *remotePoints;
350   PetscInt          lpStart, lpEnd;
351   PetscInt          numRoots, numSectionRoots, numPoints, numIndices = 0;
352   PetscInt          *localIndices;
353   PetscSFNode       *remoteIndices;
354   PetscInt          i, ind;
355   PetscErrorCode    ierr;
356 
357   PetscFunctionBegin;
358   PetscValidHeaderSpecific(sf,PETSCSF_CLASSID,1);
359   PetscValidPointer(rootSection,2);
360   /* Cannot check PetscValidIntPointer(remoteOffsets,3) because it can be NULL if sf does not reference any points in leafSection */
361   PetscValidPointer(leafSection,4);
362   PetscValidPointer(sectionSF,5);
363   ierr = PetscObjectGetComm((PetscObject)sf,&comm);CHKERRQ(ierr);
364   ierr = PetscSFCreate(comm, sectionSF);CHKERRQ(ierr);
365   ierr = PetscSectionGetChart(leafSection, &lpStart, &lpEnd);CHKERRQ(ierr);
366   ierr = PetscSectionGetStorageSize(rootSection, &numSectionRoots);CHKERRQ(ierr);
367   ierr = PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);CHKERRQ(ierr);
368   if (numRoots < 0) PetscFunctionReturn(0);
369   ierr = PetscLogEventBegin(PETSCSF_SectSF,sf,0,0,0);CHKERRQ(ierr);
370   for (i = 0; i < numPoints; ++i) {
371     PetscInt localPoint = localPoints ? localPoints[i] : i;
372     PetscInt dof;
373 
374     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
375       ierr = PetscSectionGetDof(leafSection, localPoint, &dof);CHKERRQ(ierr);
376       numIndices += dof;
377     }
378   }
379   ierr = PetscMalloc1(numIndices, &localIndices);CHKERRQ(ierr);
380   ierr = PetscMalloc1(numIndices, &remoteIndices);CHKERRQ(ierr);
381   /* Create new index graph */
382   for (i = 0, ind = 0; i < numPoints; ++i) {
383     PetscInt localPoint = localPoints ? localPoints[i] : i;
384     PetscInt rank       = remotePoints[i].rank;
385 
386     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
387       PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
388       PetscInt loff, dof, d;
389 
390       ierr = PetscSectionGetOffset(leafSection, localPoint, &loff);CHKERRQ(ierr);
391       ierr = PetscSectionGetDof(leafSection, localPoint, &dof);CHKERRQ(ierr);
392       for (d = 0; d < dof; ++d, ++ind) {
393         localIndices[ind]        = loff+d;
394         remoteIndices[ind].rank  = rank;
395         remoteIndices[ind].index = remoteOffset+d;
396       }
397     }
398   }
399   if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %D should be %D", ind, numIndices);
400   ierr = PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);CHKERRQ(ierr);
401   ierr = PetscSFSetUp(*sectionSF);CHKERRQ(ierr);
402   ierr = PetscLogEventEnd(PETSCSF_SectSF,sf,0,0,0);CHKERRQ(ierr);
403   PetscFunctionReturn(0);
404 }
405