Actual source code: network.c

petsc-master 2020-04-04
Report Typos and Errors
  1:  #include <petsc/private/dmnetworkimpl.h>

  3: /*@
  4:   DMNetworkGetPlex - Gets the Plex DM associated with this network DM

  6:   Not collective

  8:   Input Parameters:
  9: + netdm - the dm object
 10: - plexmdm - the plex dm object

 12:   Level: Advanced

 14: .seealso: DMNetworkCreate()
 15: @*/
 16: PetscErrorCode DMNetworkGetPlex(DM netdm, DM *plexdm)
 17: {
 18:   DM_Network *network = (DM_Network*) netdm->data;

 21:   *plexdm = network->plex;
 22:   return(0);
 23: }

 25: /*@
 26:   DMNetworkGetSizes - Gets the the number of subnetworks and coupling subnetworks

 28:   Collective on dm

 30:   Input Parameters:
 31: + dm - the dm object
 32: . Nsubnet - global number of subnetworks
 33: - NsubnetCouple - global number of coupling subnetworks

 35:   Level: beginner

 37: .seealso: DMNetworkCreate()
 38: @*/
 39: PetscErrorCode DMNetworkGetSizes(DM netdm, PetscInt *Nsubnet, PetscInt *Ncsubnet)
 40: {
 41:   DM_Network *network = (DM_Network*) netdm->data;

 44:   *Nsubnet = network->nsubnet;
 45:   *Ncsubnet = network->ncsubnet;
 46:   return(0);
 47: }

 49: /*@
 50:   DMNetworkSetSizes - Sets the number of subnetworks, local and global vertices and edges for each subnetwork.

 52:   Collective on dm

 54:   Input Parameters:
 55: + dm - the dm object
 56: . Nsubnet - global number of subnetworks
 57: . nV - number of local vertices for each subnetwork
 58: . nE - number of local edges for each subnetwork
 59: . NsubnetCouple - global number of coupling subnetworks
 60: - nec - number of local edges for each coupling subnetwork

 62:    You cannot change the sizes once they have been set. nV, nE are arrays of length Nsubnet, and nec is array of length NsubnetCouple.

 64:    Level: beginner

 66: .seealso: DMNetworkCreate()
 67: @*/
 68: PetscErrorCode DMNetworkSetSizes(DM dm,PetscInt Nsubnet,PetscInt nV[], PetscInt nE[],PetscInt NsubnetCouple,PetscInt nec[])
 69: {
 71:   DM_Network     *network = (DM_Network*) dm->data;
 72:   PetscInt       a[2],b[2],i;

 76:   if (Nsubnet <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subnetworks %D cannot be less than 1",Nsubnet);
 77:   if (NsubnetCouple < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of coupling subnetworks %D cannot be less than 0",NsubnetCouple);

 81:   if (network->nsubnet != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Network sizes alread set, cannot resize the network");

 83:   if (!nV || !nE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Local vertex size or edge size must be provided");

 85:   network->nsubnet  = Nsubnet + NsubnetCouple;
 86:   network->ncsubnet = NsubnetCouple;
 87:   PetscCalloc1(Nsubnet+NsubnetCouple,&network->subnet);

 89:   /* ----------------------------------------------------------
 90:    p=v or e; P=V or E
 91:    subnet[0].pStart   = 0
 92:    subnet[i+1].pStart = subnet[i].pEnd = subnet[i].pStart + (nE[i] or NV[i])
 93:    ----------------------------------------------------------------------- */
 94:   for (i=0; i < Nsubnet; i++) {
 95:     /* Get global number of vertices and edges for subnet[i] */
 96:     a[0] = nV[i]; a[1] = nE[i]; /* local number of vertices (excluding ghost) and edges */
 97:     MPIU_Allreduce(a,b,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
 98:     network->subnet[i].Nvtx = b[0];
 99:     network->subnet[i].Nedge = b[1];

101:     network->subnet[i].nvtx   = nV[i]; /* local nvtx, without ghost */

103:     /* global subnet[].vStart and vEnd, used by DMNetworkLayoutSetUp() */
104:     network->subnet[i].vStart = network->NVertices;
105:     network->subnet[i].vEnd   = network->subnet[i].vStart + network->subnet[i].Nvtx;

107:     network->nVertices += nV[i];
108:     network->NVertices += network->subnet[i].Nvtx;

110:     network->subnet[i].nedge  = nE[i];
111:     network->subnet[i].eStart = network->nEdges;
112:     network->subnet[i].eEnd   = network->subnet[i].eStart + nE[i];
113:     network->nEdges += nE[i];
114:     network->NEdges += network->subnet[i].Nedge;
115:   }

117:   /* coupling subnetwork */
118:   for (; i < Nsubnet+NsubnetCouple; i++) {
119:     /* Get global number of coupling edges for subnet[i] */
120:     MPIU_Allreduce(nec+(i-Nsubnet),&network->subnet[i].Nedge,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));

122:     network->subnet[i].nvtx   = 0; /* We design coupling subnetwork such that it does not have its own vertices */
123:     network->subnet[i].vStart = network->nVertices;
124:     network->subnet[i].vEnd   = network->subnet[i].vStart;

126:     network->subnet[i].nedge  = nec[i-Nsubnet];
127:     network->subnet[i].eStart = network->nEdges;
128:     network->subnet[i].eEnd = network->subnet[i].eStart + nec[i-Nsubnet];
129:     network->nEdges += nec[i-Nsubnet];
130:     network->NEdges += network->subnet[i].Nedge;
131:   }
132:   return(0);
133: }

135: /*@
136:   DMNetworkSetEdgeList - Sets the list of local edges (vertex connectivity) for the network

138:   Logically collective on dm

140:   Input Parameters:
141: + dm - the dm object
142: . edgelist - list of edges for each subnetwork
143: - edgelistCouple - list of edges for each coupling subnetwork

145:   Notes:
146:   There is no copy involved in this operation, only the pointer is referenced. The edgelist should
147:   not be destroyed before the call to DMNetworkLayoutSetUp

149:   Level: beginner

151:   Example usage:
152:   Consider the following 2 separate networks and a coupling network:

154: .vb
155:  network 0: v0 -> v1 -> v2 -> v3
156:  network 1: v1 -> v2 -> v0
157:  coupling network: network 1: v2 -> network 0: v0
158: .ve

160:  The resulting input
161:    edgelist[0] = [0 1 | 1 2 | 2 3];
162:    edgelist[1] = [1 2 | 2 0]
163:    edgelistCouple[0] = [(network)1 (v)2 (network)0 (v)0].

165: .seealso: DMNetworkCreate, DMNetworkSetSizes
166: @*/
167: PetscErrorCode DMNetworkSetEdgeList(DM dm,PetscInt *edgelist[],PetscInt *edgelistCouple[])
168: {
169:   DM_Network *network = (DM_Network*) dm->data;
170:   PetscInt   i;

173:   for (i=0; i < (network->nsubnet-network->ncsubnet); i++) network->subnet[i].edgelist = edgelist[i];
174:   if (network->ncsubnet) {
175:     PetscInt j = 0;
176:     if (!edgelistCouple) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Must provide edgelist_couple");
177:     while (i < network->nsubnet) network->subnet[i++].edgelist = edgelistCouple[j++];
178:   }
179:   return(0);
180: }

182: /*@
183:   DMNetworkLayoutSetUp - Sets up the bare layout (graph) for the network

185:   Collective on dm

187:   Input Parameters:
188: . DM - the dmnetwork object

190:   Notes:
191:   This routine should be called after the network sizes and edgelists have been provided. It creates
192:   the bare layout of the network and sets up the network to begin insertion of components.

194:   All the components should be registered before calling this routine.

196:   Level: beginner

198: .seealso: DMNetworkSetSizes, DMNetworkSetEdgeList
199: @*/
200: PetscErrorCode DMNetworkLayoutSetUp(DM dm)
201: {
203:   DM_Network     *network = (DM_Network*)dm->data;
204:   PetscInt       numCorners=2,spacedim=2,dim = 1; /* One dimensional network */
205:   PetscReal      *vertexcoords=NULL;
206:   PetscInt       i,j,ctr,nsubnet,*eowners,np,*edges,*subnetvtx,vStart;
207:   PetscInt       k,netid,vid, *vidxlTog,*edgelist_couple=NULL;
208:   const PetscInt *cone;
209:   MPI_Comm       comm;
210:   PetscMPIInt    size,rank;

213:   PetscObjectGetComm((PetscObject)dm,&comm);
214:   MPI_Comm_rank(comm,&rank);
215:   MPI_Comm_size(comm,&size);

217:   /* Create the local edgelist for the network by concatenating local input edgelists of the subnetworks */
218:   PetscCalloc2(numCorners*network->nVertices,&vertexcoords,2*network->nEdges,&edges);
219:   nsubnet = network->nsubnet - network->ncsubnet;
220:   ctr = 0;
221:   for (i=0; i < nsubnet; i++) {
222:     for (j = 0; j < network->subnet[i].nedge; j++) {
223:       edges[2*ctr]   = network->subnet[i].vStart + network->subnet[i].edgelist[2*j];
224:       edges[2*ctr+1] = network->subnet[i].vStart + network->subnet[i].edgelist[2*j+1];
225:       ctr++;
226:     }
227:   }

229:   /* Append local coupling edgelists of the subnetworks */
230:   i       = nsubnet; /* netid of coupling subnet */
231:   nsubnet = network->nsubnet;
232:   while (i < nsubnet) {
233:     edgelist_couple = network->subnet[i].edgelist;

235:     k = 0;
236:     for (j = 0; j < network->subnet[i].nedge; j++) {
237:       netid = edgelist_couple[k]; vid = edgelist_couple[k+1];
238:       edges[2*ctr] = network->subnet[netid].vStart + vid; k += 2;

240:       netid = edgelist_couple[k]; vid = edgelist_couple[k+1];
241:       edges[2*ctr+1] = network->subnet[netid].vStart + vid; k+=2;
242:       ctr++;
243:     }
244:     i++;
245:   }
246:   /*
247:   if (rank == 0) {
248:     PetscPrintf(PETSC_COMM_SELF,"[%d] edgelist:\n",rank);
249:     for(i=0; i < network->nEdges; i++) {
250:       PetscPrintf(PETSC_COMM_SELF,"[%D %D]",edges[2*i],edges[2*i+1]);
251:       printf("\n");
252:     }
253:   }
254:    */

256:   /* Create network->plex */
257: #if defined(PETSC_USE_64BIT_INDICES)
258:   {
259:     int *edges64;
260:     np = network->nEdges*numCorners;
261:     PetscMalloc1(np,&edges64);
262:     for (i=0; i<np; i++) edges64[i] = (int)edges[i];

264:     if (size == 1) {
265:       DMPlexCreateFromCellList(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges64,spacedim,(const double*)vertexcoords,&network->plex);
266:     } else {
267:       DMPlexCreateFromCellListParallel(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges64,spacedim,(const PetscReal*)vertexcoords,NULL,&network->plex);
268:     }
269:     PetscFree(edges64);
270:   }
271: #else
272:   if (size == 1) {
273:     DMPlexCreateFromCellList(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges,spacedim,(const double*)vertexcoords,&network->plex);
274:   } else {
275:     DMPlexCreateFromCellListParallel(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges,spacedim,(const PetscReal*)vertexcoords,NULL,&network->plex);
276:   }
277: #endif

279:   DMPlexGetChart(network->plex,&network->pStart,&network->pEnd);
280:   DMPlexGetHeightStratum(network->plex,0,&network->eStart,&network->eEnd);
281:   DMPlexGetHeightStratum(network->plex,1,&network->vStart,&network->vEnd);
282:   vStart = network->vStart;

284:   PetscSectionCreate(comm,&network->DataSection);
285:   PetscSectionCreate(comm,&network->DofSection);
286:   PetscSectionSetChart(network->DataSection,network->pStart,network->pEnd);
287:   PetscSectionSetChart(network->DofSection,network->pStart,network->pEnd);

289:   network->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
290:   np = network->pEnd - network->pStart;
291:   PetscCalloc2(np,&network->header,np,&network->cvalue);

293:   /* Create vidxlTog: maps local vertex index to global index */
294:   np = network->vEnd - vStart;
295:   PetscMalloc2(np,&vidxlTog,size+1,&eowners);
296:   ctr = 0;
297:   for (i=network->eStart; i<network->eEnd; i++) {
298:     DMNetworkGetConnectedVertices(dm,i,&cone);
299:     vidxlTog[cone[0] - vStart] = edges[2*ctr];
300:     vidxlTog[cone[1] - vStart] = edges[2*ctr+1];
301:     ctr++;
302:   }
303:   PetscFree2(vertexcoords,edges);

305:   /* Create vertices and edges array for the subnetworks */
306:   for (j=0; j < network->nsubnet; j++) {
307:     PetscCalloc1(network->subnet[j].nedge,&network->subnet[j].edges);

309:     /* Temporarily setting nvtx and nedge to 0 so we can use them as counters in the below for loop.
310:        These get updated when the vertices and edges are added. */
311:     network->subnet[j].nvtx  = 0;
312:     network->subnet[j].nedge = 0;
313:   }
314:   PetscCalloc1(np,&network->subnetvtx);


317:   /* Get edge ownership */
318:   np = network->eEnd - network->eStart;
319:   MPI_Allgather(&np,1,MPIU_INT,eowners+1,1,MPIU_INT,comm);
320:   eowners[0] = 0;
321:   for (i=2; i<=size; i++) eowners[i] += eowners[i-1];

323:   i = 0; j = 0;
324:   while (i < np) { /* local edges, including coupling edges */
325:     network->header[i].index = i + eowners[rank];   /* Global edge index */

327:     if (j < network->nsubnet && i < network->subnet[j].eEnd) {
328:       network->header[i].subnetid = j; /* Subnetwork id */
329:       network->subnet[j].edges[network->subnet[j].nedge++] = i;

331:       network->header[i].ndata = 0;
332:       PetscSectionAddDof(network->DataSection,i,network->dataheadersize);
333:       network->header[i].offset[0] = 0;
334:       network->header[i].offsetvarrel[0] = 0;
335:       i++;
336:     }
337:     if (i >= network->subnet[j].eEnd) j++;
338:   }

340:   /* Count network->subnet[*].nvtx */
341:   for (i=vStart; i<network->vEnd; i++) { /* local vertices, including ghosts */
342:     k = vidxlTog[i-vStart];
343:     for (j=0; j < network->nsubnet; j++) {
344:       if (network->subnet[j].vStart <= k && k < network->subnet[j].vEnd) {
345:         network->subnet[j].nvtx++;
346:         break;
347:       }
348:     }
349:   }

351:   /* Set network->subnet[*].vertices on array network->subnetvtx */
352:   subnetvtx = network->subnetvtx;
353:   for (j=0; j<network->nsubnet; j++) {
354:     network->subnet[j].vertices = subnetvtx;
355:     subnetvtx                  += network->subnet[j].nvtx;
356:     network->subnet[j].nvtx = 0;
357:   }

359:   /* Set vertex array for the subnetworks */
360:   for (i=vStart; i<network->vEnd; i++) { /* local vertices, including ghosts */
361:     network->header[i].index = vidxlTog[i-vStart]; /*  Global vertex index */

363:     k = vidxlTog[i-vStart];
364:     for (j=0; j < network->nsubnet; j++) {
365:       if (network->subnet[j].vStart <= k && k < network->subnet[j].vEnd) {
366:         network->header[i].subnetid = j;
367:         network->subnet[j].vertices[network->subnet[j].nvtx++] = i;
368:         break;
369:       }
370:     }

372:     network->header[i].ndata = 0;
373:     PetscSectionAddDof(network->DataSection,i,network->dataheadersize);
374:     network->header[i].offset[0] = 0;
375:     network->header[i].offsetvarrel[0] = 0;
376:   }

378:   PetscFree2(vidxlTog,eowners);
379:   return(0);
380: }

382: /*@C
383:   DMNetworkGetSubnetworkInfo - Returns the info for the subnetwork

385:   Input Parameters:
386: + dm - the DM object
387: - id   - the ID (integer) of the subnetwork

389:   Output Parameters:
390: + nv    - number of vertices (local)
391: . ne    - number of edges (local)
392: . vtx   - local vertices for this subnetwork
393: - edge  - local edges for this subnetwork

395:   Notes:
396:   Cannot call this routine before DMNetworkLayoutSetup()

398:   Level: intermediate

400: .seealso: DMNetworkLayoutSetUp, DMNetworkCreate
401: @*/
402: PetscErrorCode DMNetworkGetSubnetworkInfo(DM dm,PetscInt id,PetscInt *nv, PetscInt *ne,const PetscInt **vtx, const PetscInt **edge)
403: {
404:   DM_Network *network = (DM_Network*)dm->data;

407:   if (id >= network->nsubnet) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subnet ID %D exceeds the num of subnets %D",id,network->nsubnet);
408:   *nv   = network->subnet[id].nvtx;
409:   *ne   = network->subnet[id].nedge;
410:   *vtx  = network->subnet[id].vertices;
411:   *edge = network->subnet[id].edges;
412:   return(0);
413: }

415: /*@C
416:   DMNetworkGetSubnetworkCoupleInfo - Returns the info for the coupling subnetwork

418:   Input Parameters:
419: + dm - the DM object
420: - id   - the ID (integer) of the coupling subnetwork

422:   Output Parameters:
423: + ne - number of edges (local)
424: - edge  - local edges for this coupling subnetwork

426:   Notes:
427:   Cannot call this routine before DMNetworkLayoutSetup()

429:   Level: intermediate

431: .seealso: DMNetworkGetSubnetworkInfo, DMNetworkLayoutSetUp, DMNetworkCreate
432: @*/
433: PetscErrorCode DMNetworkGetSubnetworkCoupleInfo(DM dm,PetscInt id,PetscInt *ne,const PetscInt **edge)
434: {
435:   DM_Network *net = (DM_Network*)dm->data;
436:   PetscInt   id1;

439:   if (net->ncsubnet) {
440:     if (id >= net->ncsubnet) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subnet ID %D exceeds the num of coupling subnets %D",id,net->ncsubnet);

442:     id1   = id + net->nsubnet - net->ncsubnet;
443:     *ne   = net->subnet[id1].nedge;
444:     *edge = net->subnet[id1].edges;
445:   } else {
446:     *ne   = 0;
447:     *edge = NULL;
448:   }
449:   return(0);
450: }

452: /*@C
453:   DMNetworkRegisterComponent - Registers the network component

455:   Logically collective on dm

457:   Input Parameters:
458: + dm   - the network object
459: . name - the component name
460: - size - the storage size in bytes for this component data

462:    Output Parameters:
463: .   key - an integer key that defines the component

465:    Notes
466:    This routine should be called by all processors before calling DMNetworkLayoutSetup().

468:    Level: beginner

470: .seealso: DMNetworkLayoutSetUp, DMNetworkCreate
471: @*/
472: PetscErrorCode DMNetworkRegisterComponent(DM dm,const char *name,size_t size,PetscInt *key)
473: {
474:   PetscErrorCode        ierr;
475:   DM_Network            *network = (DM_Network*) dm->data;
476:   DMNetworkComponent    *component=&network->component[network->ncomponent];
477:   PetscBool             flg=PETSC_FALSE;
478:   PetscInt              i;

481:   for (i=0; i < network->ncomponent; i++) {
482:     PetscStrcmp(component->name,name,&flg);
483:     if (flg) {
484:       *key = i;
485:       return(0);
486:     }
487:   }
488:   if(network->ncomponent == MAX_COMPONENTS) {
489:     SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Number of components registered exceeds the max %D",MAX_COMPONENTS);
490:   }

492:   PetscStrcpy(component->name,name);
493:   component->size = size/sizeof(DMNetworkComponentGenericDataType);
494:   *key = network->ncomponent;
495:   network->ncomponent++;
496:   return(0);
497: }

499: /*@
500:   DMNetworkGetVertexRange - Get the bounds [start, end) for the vertices.

502:   Not Collective

504:   Input Parameters:
505: . dm - The DMNetwork object

507:   Output Parameters:
508: + vStart - The first vertex point
509: - vEnd   - One beyond the last vertex point

511:   Level: beginner

513: .seealso: DMNetworkGetEdgeRange
514: @*/
515: PetscErrorCode DMNetworkGetVertexRange(DM dm,PetscInt *vStart,PetscInt *vEnd)
516: {
517:   DM_Network     *network = (DM_Network*)dm->data;

520:   if (vStart) *vStart = network->vStart;
521:   if (vEnd) *vEnd = network->vEnd;
522:   return(0);
523: }

525: /*@
526:   DMNetworkGetEdgeRange - Get the bounds [start, end) for the edges.

528:   Not Collective

530:   Input Parameters:
531: . dm - The DMNetwork object

533:   Output Parameters:
534: + eStart - The first edge point
535: - eEnd   - One beyond the last edge point

537:   Level: beginner

539: .seealso: DMNetworkGetVertexRange
540: @*/
541: PetscErrorCode DMNetworkGetEdgeRange(DM dm,PetscInt *eStart,PetscInt *eEnd)
542: {
543:   DM_Network     *network = (DM_Network*)dm->data;

546:   if (eStart) *eStart = network->eStart;
547:   if (eEnd) *eEnd = network->eEnd;
548:   return(0);
549: }

551: /*@
552:   DMNetworkGetGlobalEdgeIndex - Get the user global numbering for the edge.

554:   Not Collective

556:   Input Parameters:
557: + dm - DMNetwork object
558: - p  - edge point

560:   Output Parameters:
561: . index - user global numbering for the edge

563:   Level: intermediate

565: .seealso: DMNetworkGetGlobalVertexIndex
566: @*/
567: PetscErrorCode DMNetworkGetGlobalEdgeIndex(DM dm,PetscInt p,PetscInt *index)
568: {
569:   PetscErrorCode    ierr;
570:   DM_Network        *network = (DM_Network*)dm->data;
571:   PetscInt          offsetp;
572:   DMNetworkComponentHeader header;

575:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
576:   PetscSectionGetOffset(network->DataSection,p,&offsetp);
577:   header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
578:   *index = header->index;
579:   return(0);
580: }

582: /*@
583:   DMNetworkGetGlobalVertexIndex - Get the user global numbering for the vertex.

585:   Not Collective

587:   Input Parameters:
588: + dm - DMNetwork object
589: - p  - vertex point

591:   Output Parameters:
592: . index - user global numbering for the vertex

594:   Level: intermediate

596: .seealso: DMNetworkGetGlobalEdgeIndex
597: @*/
598: PetscErrorCode DMNetworkGetGlobalVertexIndex(DM dm,PetscInt p,PetscInt *index)
599: {
600:   PetscErrorCode    ierr;
601:   DM_Network        *network = (DM_Network*)dm->data;
602:   PetscInt          offsetp;
603:   DMNetworkComponentHeader header;

606:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
607:   PetscSectionGetOffset(network->DataSection,p,&offsetp);
608:   header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
609:   *index = header->index;
610:   return(0);
611: }

613: /*
614:   DMNetworkGetComponentKeyOffset - Gets the type along with the offset for indexing the
615:                                     component value from the component data array

617:   Not Collective

619:   Input Parameters:
620: + dm      - The DMNetwork object
621: . p       - vertex/edge point
622: - compnum - component number

624:   Output Parameters:
625: + compkey - the key obtained when registering the component
626: - offset  - offset into the component data array associated with the vertex/edge point

628:   Notes:
629:   Typical usage:

631:   DMNetworkGetComponentDataArray(dm, &arr);
632:   DMNetworkGetVertex/EdgeRange(dm,&Start,&End);
633:   Loop over vertices or edges
634:     DMNetworkGetNumComponents(dm,v,&numcomps);
635:     Loop over numcomps
636:       DMNetworkGetComponentKeyOffset(dm,v,compnum,&key,&offset);
637:       compdata = (UserCompDataType)(arr+offset);

639:   Level: intermediate

641: .seealso: DMNetworkGetNumComponents, DMNetworkGetComponentDataArray,
642: */
643: PetscErrorCode DMNetworkGetComponentKeyOffset(DM dm,PetscInt p, PetscInt compnum, PetscInt *compkey, PetscInt *offset)
644: {
645:   PetscErrorCode           ierr;
646:   PetscInt                 offsetp;
647:   DMNetworkComponentHeader header;
648:   DM_Network               *network = (DM_Network*)dm->data;

651:   PetscSectionGetOffset(network->DataSection,p,&offsetp);
652:   header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
653:   if (compkey) *compkey = header->key[compnum];
654:   if (offset) *offset  = offsetp+network->dataheadersize+header->offset[compnum];
655:   return(0);
656: }

658: /*@
659:   DMNetworkGetComponent - Returns the network component and its key

661:   Not Collective

663:   Input Parameters:
664: + dm - DMNetwork object
665: . p  - edge or vertex point
666: - compnum - component number

668:   Output Parameters:
669: + compkey - the key set for this computing during registration
670: - component - the component data

672:   Notes:
673:   Typical usage:

675:   DMNetworkGetVertex/EdgeRange(dm,&Start,&End);
676:   Loop over vertices or edges
677:     DMNetworkGetNumComponents(dm,v,&numcomps);
678:     Loop over numcomps
679:       DMNetworkGetComponent(dm,v,compnum,&key,&component);

681:   Level: beginner

683: .seealso: DMNetworkGetNumComponents, DMNetworkGetVariableOffset
684: @*/
685: PetscErrorCode DMNetworkGetComponent(DM dm, PetscInt p, PetscInt compnum, PetscInt *key, void **component)
686: {
688:   DM_Network     *network = (DM_Network*)dm->data;
689:   PetscInt       offsetd = 0;

692:   DMNetworkGetComponentKeyOffset(dm,p,compnum,key,&offsetd);
693:   *component = network->componentdataarray+offsetd;
694:   return(0);
695: }

697: /*@
698:   DMNetworkAddComponent - Adds a network component at the given point (vertex/edge)

700:   Not Collective

702:   Input Parameters:
703: + dm           - The DMNetwork object
704: . p            - vertex/edge point
705: . componentkey - component key returned while registering the component
706: - compvalue    - pointer to the data structure for the component

708:   Level: beginner

710: .seealso: DMNetworkGetVertexRange, DMNetworkGetEdgeRange, DMNetworkRegisterComponent
711: @*/
712: PetscErrorCode DMNetworkAddComponent(DM dm, PetscInt p,PetscInt componentkey,void* compvalue)
713: {
714:   DM_Network               *network = (DM_Network*)dm->data;
715:   DMNetworkComponent       *component = &network->component[componentkey];
716:   DMNetworkComponentHeader header = &network->header[p];
717:   DMNetworkComponentValue  cvalue = &network->cvalue[p];
718:   PetscErrorCode           ierr;

721:   if (header->ndata == MAX_DATA_AT_POINT) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Number of components at a point exceeds the max %D",MAX_DATA_AT_POINT);

723:   header->size[header->ndata] = component->size;
724:   PetscSectionAddDof(network->DataSection,p,component->size);
725:   header->key[header->ndata] = componentkey;
726:   if (header->ndata != 0) header->offset[header->ndata] = header->offset[header->ndata-1] + header->size[header->ndata-1];
727:   header->nvar[header->ndata] = 0;

729:   cvalue->data[header->ndata] = (void*)compvalue;
730:   header->ndata++;
731:   return(0);
732: }

734: /*@
735:   DMNetworkSetComponentNumVariables - Sets the number of variables for a component

737:   Not Collective

739:   Input Parameters:
740: + dm           - The DMNetwork object
741: . p            - vertex/edge point
742: . compnum      - component number (First component added = 0, second = 1, ...)
743: - nvar         - number of variables for the component

745:   Level: beginner

747: .seealso: DMNetworkAddComponent(), DMNetworkGetNumComponents(),DMNetworkRegisterComponent()
748: @*/
749: PetscErrorCode DMNetworkSetComponentNumVariables(DM dm, PetscInt p,PetscInt compnum,PetscInt nvar)
750: {
751:   DM_Network               *network = (DM_Network*)dm->data;
752:   DMNetworkComponentHeader header = &network->header[p];
753:   PetscErrorCode           ierr;

756:   DMNetworkAddNumVariables(dm,p,nvar);
757:   header->nvar[compnum] = nvar;
758:   if (compnum != 0) header->offsetvarrel[compnum] = header->offsetvarrel[compnum-1] + header->nvar[compnum-1];
759:   return(0);
760: }

762: /*@
763:   DMNetworkGetNumComponents - Get the number of components at a vertex/edge

765:   Not Collective

767:   Input Parameters:
768: + dm - The DMNetwork object
769: - p  - vertex/edge point

771:   Output Parameters:
772: . numcomponents - Number of components at the vertex/edge

774:   Level: beginner

776: .seealso: DMNetworkRegisterComponent, DMNetworkAddComponent
777: @*/
778: PetscErrorCode DMNetworkGetNumComponents(DM dm,PetscInt p,PetscInt *numcomponents)
779: {
781:   PetscInt       offset;
782:   DM_Network     *network = (DM_Network*)dm->data;

785:   PetscSectionGetOffset(network->DataSection,p,&offset);
786:   *numcomponents = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata;
787:   return(0);
788: }

790: /*@
791:   DMNetworkGetVariableOffset - Get the offset for accessing the variable associated with the given vertex/edge from the local vector.

793:   Not Collective

795:   Input Parameters:
796: + dm     - The DMNetwork object
797: - p      - the edge/vertex point

799:   Output Parameters:
800: . offset - the offset

802:   Level: beginner

804: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
805: @*/
806: PetscErrorCode DMNetworkGetVariableOffset(DM dm,PetscInt p,PetscInt *offset)
807: {
809:   DM_Network     *network = (DM_Network*)dm->data;

812:   PetscSectionGetOffset(network->plex->localSection,p,offset);
813:   return(0);
814: }

816: /*@
817:   DMNetworkGetVariableGlobalOffset - Get the global offset for the variable associated with the given vertex/edge from the global vector.

819:   Not Collective

821:   Input Parameters:
822: + dm      - The DMNetwork object
823: - p       - the edge/vertex point

825:   Output Parameters:
826: . offsetg - the offset

828:   Level: beginner

830: .seealso: DMNetworkGetVariableOffset, DMGetLocalVector
831: @*/
832: PetscErrorCode DMNetworkGetVariableGlobalOffset(DM dm,PetscInt p,PetscInt *offsetg)
833: {
835:   DM_Network     *network = (DM_Network*)dm->data;

838:   PetscSectionGetOffset(network->plex->globalSection,p,offsetg);
839:   if (*offsetg < 0) *offsetg = -(*offsetg + 1); /* Convert to actual global offset for ghost vertex */
840:   return(0);
841: }

843: /*@
844:   DMNetworkGetComponentVariableOffset - Get the offset for accessing the variable associated with a component for the given vertex/edge from the local vector.

846:   Not Collective

848:   Input Parameters:
849: + dm     - The DMNetwork object
850: . compnum - component number
851: - p      - the edge/vertex point

853:   Output Parameters:
854: . offset - the offset

856:   Level: intermediate

858: .seealso: DMNetworkGetVariableGlobalOffset(), DMGetLocalVector(), DMNetworkSetComponentNumVariables()
859: @*/
860: PetscErrorCode DMNetworkGetComponentVariableOffset(DM dm,PetscInt p,PetscInt compnum,PetscInt *offset)
861: {
863:   DM_Network     *network = (DM_Network*)dm->data;
864:   PetscInt       offsetp,offsetd;
865:   DMNetworkComponentHeader header;

868:   DMNetworkGetVariableOffset(dm,p,&offsetp);
869:   PetscSectionGetOffset(network->DataSection,p,&offsetd);
870:   header = (DMNetworkComponentHeader)(network->componentdataarray+offsetd);
871:   *offset = offsetp + header->offsetvarrel[compnum];
872:   return(0);
873: }

875: /*@
876:   DMNetworkGetComponentVariableGlobalOffset - Get the global offset for accessing the variable associated with a component for the given vertex/edge from the local vector.

878:   Not Collective

880:   Input Parameters:
881: + dm     - The DMNetwork object
882: . compnum - component number
883: - p      - the edge/vertex point

885:   Output Parameters:
886: . offsetg - the global offset

888:   Level: intermediate

890: .seealso: DMNetworkGetVariableGlobalOffset(), DMNetworkGetComponentVariableOffset(), DMGetLocalVector(), DMNetworkSetComponentNumVariables()
891: @*/
892: PetscErrorCode DMNetworkGetComponentVariableGlobalOffset(DM dm,PetscInt p,PetscInt compnum,PetscInt *offsetg)
893: {
895:   DM_Network     *network = (DM_Network*)dm->data;
896:   PetscInt       offsetp,offsetd;
897:   DMNetworkComponentHeader header;

900:   DMNetworkGetVariableGlobalOffset(dm,p,&offsetp);
901:   PetscSectionGetOffset(network->DataSection,p,&offsetd);
902:   header = (DMNetworkComponentHeader)(network->componentdataarray+offsetd);
903:   *offsetg = offsetp + header->offsetvarrel[compnum];
904:   return(0);
905: }

907: /*@
908:   DMNetworkGetEdgeOffset - Get the offset for accessing the variable associated with the given edge from the local subvector.

910:   Not Collective

912:   Input Parameters:
913: + dm     - The DMNetwork object
914: - p      - the edge point

916:   Output Parameters:
917: . offset - the offset

919:   Level: intermediate

921: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
922: @*/
923: PetscErrorCode DMNetworkGetEdgeOffset(DM dm,PetscInt p,PetscInt *offset)
924: {
926:   DM_Network     *network = (DM_Network*)dm->data;


930:   PetscSectionGetOffset(network->edge.DofSection,p,offset);
931:   return(0);
932: }

934: /*@
935:   DMNetworkGetVertexOffset - Get the offset for accessing the variable associated with the given vertex from the local subvector.

937:   Not Collective

939:   Input Parameters:
940: + dm     - The DMNetwork object
941: - p      - the vertex point

943:   Output Parameters:
944: . offset - the offset

946:   Level: intermediate

948: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
949: @*/
950: PetscErrorCode DMNetworkGetVertexOffset(DM dm,PetscInt p,PetscInt *offset)
951: {
953:   DM_Network     *network = (DM_Network*)dm->data;


957:   p -= network->vStart;

959:   PetscSectionGetOffset(network->vertex.DofSection,p,offset);
960:   return(0);
961: }
962: /*@
963:   DMNetworkAddNumVariables - Add number of variables associated with a given point.

965:   Not Collective

967:   Input Parameters:
968: + dm   - The DMNetworkObject
969: . p    - the vertex/edge point
970: - nvar - number of additional variables

972:   Level: beginner

974: .seealso: DMNetworkSetNumVariables
975: @*/
976: PetscErrorCode DMNetworkAddNumVariables(DM dm,PetscInt p,PetscInt nvar)
977: {
979:   DM_Network     *network = (DM_Network*)dm->data;

982:   PetscSectionAddDof(network->DofSection,p,nvar);
983:   return(0);
984: }

986: /*@
987:   DMNetworkGetNumVariables - Gets number of variables for a vertex/edge point.

989:   Not Collective

991:   Input Parameters:
992: + dm   - The DMNetworkObject
993: - p    - the vertex/edge point

995:   Output Parameters:
996: . nvar - number of variables

998:   Level: beginner

1000: .seealso: DMNetworkAddNumVariables, DMNetworkSddNumVariables
1001: @*/
1002: PetscErrorCode DMNetworkGetNumVariables(DM dm,PetscInt p,PetscInt *nvar)
1003: {
1005:   DM_Network     *network = (DM_Network*)dm->data;

1008:   PetscSectionGetDof(network->DofSection,p,nvar);
1009:   return(0);
1010: }

1012: /*@
1013:   DMNetworkSetNumVariables - Sets number of variables for a vertex/edge point.

1015:   Not Collective

1017:   Input Parameters:
1018: + dm   - The DMNetworkObject
1019: . p    - the vertex/edge point
1020: - nvar - number of variables

1022:   Level: beginner

1024: .seealso: DMNetworkAddNumVariables
1025: @*/
1026: PetscErrorCode DMNetworkSetNumVariables(DM dm,PetscInt p,PetscInt nvar)
1027: {
1029:   DM_Network     *network = (DM_Network*)dm->data;

1032:   PetscSectionSetDof(network->DofSection,p,nvar);
1033:   return(0);
1034: }

1036: /* Sets up the array that holds the data for all components and its associated section. This
1037:    function is called during DMSetUp() */
1038: PetscErrorCode DMNetworkComponentSetUp(DM dm)
1039: {
1040:   PetscErrorCode           ierr;
1041:   DM_Network               *network = (DM_Network*)dm->data;
1042:   PetscInt                 arr_size,p,offset,offsetp,ncomp,i;
1043:   DMNetworkComponentHeader header;
1044:   DMNetworkComponentValue  cvalue;
1045:   DMNetworkComponentGenericDataType *componentdataarray;

1048:   PetscSectionSetUp(network->DataSection);
1049:   PetscSectionGetStorageSize(network->DataSection,&arr_size);
1050:   PetscMalloc1(arr_size,&network->componentdataarray);
1051:   componentdataarray = network->componentdataarray;
1052:   for (p = network->pStart; p < network->pEnd; p++) {
1053:     PetscSectionGetOffset(network->DataSection,p,&offsetp);
1054:     /* Copy header */
1055:     header = &network->header[p];
1056:     PetscMemcpy(componentdataarray+offsetp,header,network->dataheadersize*sizeof(DMNetworkComponentGenericDataType));
1057:     /* Copy data */
1058:     cvalue = &network->cvalue[p];
1059:     ncomp = header->ndata;
1060:     for (i = 0; i < ncomp; i++) {
1061:       offset = offsetp + network->dataheadersize + header->offset[i];
1062:       PetscMemcpy(componentdataarray+offset,cvalue->data[i],header->size[i]*sizeof(DMNetworkComponentGenericDataType));
1063:     }
1064:   }
1065:   return(0);
1066: }

1068: /* Sets up the section for dofs. This routine is called during DMSetUp() */
1069: PetscErrorCode DMNetworkVariablesSetUp(DM dm)
1070: {
1072:   DM_Network     *network = (DM_Network*)dm->data;

1075:   PetscSectionSetUp(network->DofSection);
1076:   return(0);
1077: }

1079: /*
1080:   DMNetworkGetComponentDataArray - Returns the component data array

1082:   Not Collective

1084:   Input Parameters:
1085: . dm - The DMNetwork Object

1087:   Output Parameters:
1088: . componentdataarray - array that holds data for all components

1090:   Level: intermediate

1092: .seealso: DMNetworkGetComponentKeyOffset, DMNetworkGetNumComponents
1093: */
1094: PetscErrorCode DMNetworkGetComponentDataArray(DM dm,DMNetworkComponentGenericDataType **componentdataarray)
1095: {
1096:   DM_Network     *network = (DM_Network*)dm->data;

1099:   *componentdataarray = network->componentdataarray;
1100:   return(0);
1101: }

1103: /* Get a subsection from a range of points */
1104: PetscErrorCode DMNetworkGetSubSection_private(PetscSection master, PetscInt pstart, PetscInt pend,PetscSection *subsection)
1105: {
1107:   PetscInt       i, nvar;

1110:   PetscSectionCreate(PetscObjectComm((PetscObject)master), subsection);
1111:   PetscSectionSetChart(*subsection, 0, pend - pstart);
1112:   for (i = pstart; i < pend; i++) {
1113:     PetscSectionGetDof(master,i,&nvar);
1114:     PetscSectionSetDof(*subsection, i - pstart, nvar);
1115:   }

1117:   PetscSectionSetUp(*subsection);
1118:   return(0);
1119: }

1121: /* Create a submap of points with a GlobalToLocal structure */
1122: PetscErrorCode DMNetworkSetSubMap_private(PetscInt pstart, PetscInt pend, ISLocalToGlobalMapping *map)
1123: {
1125:   PetscInt       i, *subpoints;

1128:   /* Create index sets to map from "points" to "subpoints" */
1129:   PetscMalloc1(pend - pstart, &subpoints);
1130:   for (i = pstart; i < pend; i++) {
1131:     subpoints[i - pstart] = i;
1132:   }
1133:   ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,pend-pstart,subpoints,PETSC_COPY_VALUES,map);
1134:   PetscFree(subpoints);
1135:   return(0);
1136: }

1138: /*@
1139:   DMNetworkAssembleGraphStructures - Assembles vertex and edge data structures. Must be called after DMNetworkDistribute.

1141:   Collective

1143:   Input Parameters:
1144: . dm   - The DMNetworkObject

1146:   Note: the routine will create alternative orderings for the vertices and edges. Assume global network points are:

1148:   points = [0 1 2 3 4 5 6]

1150:   where edges = [0,1,2,3] and vertices = [4,5,6]. The new orderings will be specific to the subset (i.e vertices = [0,1,2] <- [4,5,6]).

1152:   With this new ordering a local PetscSection, global PetscSection and PetscSF will be created specific to the subset.

1154:   Level: intermediate

1156: @*/
1157: PetscErrorCode DMNetworkAssembleGraphStructures(DM dm)
1158: {
1160:   MPI_Comm       comm;
1161:   PetscMPIInt    rank, size;
1162:   DM_Network     *network = (DM_Network*)dm->data;

1165:   PetscObjectGetComm((PetscObject)dm,&comm);
1166:   MPI_Comm_rank(comm, &rank);
1167:   MPI_Comm_size(comm, &size);

1169:   /* Create maps for vertices and edges */
1170:   DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);
1171:   DMNetworkSetSubMap_private(network->eStart,network->eEnd,&network->edge.mapping);

1173:   /* Create local sub-sections */
1174:   DMNetworkGetSubSection_private(network->DofSection,network->vStart,network->vEnd,&network->vertex.DofSection);
1175:   DMNetworkGetSubSection_private(network->DofSection,network->eStart,network->eEnd,&network->edge.DofSection);

1177:   if (size > 1) {
1178:     PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);

1180:     PetscSectionCreateGlobalSection(network->vertex.DofSection, network->vertex.sf, PETSC_FALSE, PETSC_FALSE, &network->vertex.GlobalDofSection);
1181:     PetscSFGetSubSF(network->plex->sf, network->edge.mapping, &network->edge.sf);
1182:     PetscSectionCreateGlobalSection(network->edge.DofSection, network->edge.sf, PETSC_FALSE, PETSC_FALSE, &network->edge.GlobalDofSection);
1183:   } else {
1184:     /* create structures for vertex */
1185:     PetscSectionClone(network->vertex.DofSection,&network->vertex.GlobalDofSection);
1186:     /* create structures for edge */
1187:     PetscSectionClone(network->edge.DofSection,&network->edge.GlobalDofSection);
1188:   }

1190:   /* Add viewers */
1191:   PetscObjectSetName((PetscObject)network->edge.GlobalDofSection,"Global edge dof section");
1192:   PetscObjectSetName((PetscObject)network->vertex.GlobalDofSection,"Global vertex dof section");
1193:   PetscSectionViewFromOptions(network->edge.GlobalDofSection, NULL, "-edge_global_section_view");
1194:   PetscSectionViewFromOptions(network->vertex.GlobalDofSection, NULL, "-vertex_global_section_view");
1195:   return(0);
1196: }

1198: /*@
1199:   DMNetworkDistribute - Distributes the network and moves associated component data.

1201:   Collective

1203:   Input Parameter:
1204: + DM - the DMNetwork object
1205: - overlap - The overlap of partitions, 0 is the default

1207:   Notes:
1208:   Distributes the network with <overlap>-overlapping partitioning of the edges.

1210:   Level: intermediate

1212: .seealso: DMNetworkCreate
1213: @*/
1214: PetscErrorCode DMNetworkDistribute(DM *dm,PetscInt overlap)
1215: {
1216:   MPI_Comm       comm;
1218:   PetscMPIInt    size;
1219:   DM_Network     *oldDMnetwork = (DM_Network*)((*dm)->data);
1220:   DM_Network     *newDMnetwork;
1221:   PetscSF        pointsf=NULL;
1222:   DM             newDM;
1223:   PetscInt       j,e,v,offset,*subnetvtx;
1224:   PetscPartitioner         part;
1225:   DMNetworkComponentHeader header;

1228:   PetscObjectGetComm((PetscObject)*dm,&comm);
1229:   MPI_Comm_size(comm, &size);
1230:   if (size == 1) return(0);

1232:   DMNetworkCreate(PetscObjectComm((PetscObject)*dm),&newDM);
1233:   newDMnetwork = (DM_Network*)newDM->data;
1234:   newDMnetwork->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);

1236:   /* Enable runtime options for petscpartitioner */
1237:   DMPlexGetPartitioner(oldDMnetwork->plex,&part);
1238:   PetscPartitionerSetFromOptions(part);

1240:   /* Distribute plex dm and dof section */
1241:   DMPlexDistribute(oldDMnetwork->plex,overlap,&pointsf,&newDMnetwork->plex);

1243:   /* Distribute dof section */
1244:   PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DofSection);
1245:   PetscSFDistributeSection(pointsf,oldDMnetwork->DofSection,NULL,newDMnetwork->DofSection);
1246:   PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DataSection);

1248:   /* Distribute data and associated section */
1249:   DMPlexDistributeData(newDMnetwork->plex,pointsf,oldDMnetwork->DataSection,MPIU_INT,(void*)oldDMnetwork->componentdataarray,newDMnetwork->DataSection,(void**)&newDMnetwork->componentdataarray);

1251:   PetscSectionGetChart(newDMnetwork->DataSection,&newDMnetwork->pStart,&newDMnetwork->pEnd);
1252:   DMPlexGetHeightStratum(newDMnetwork->plex,0, &newDMnetwork->eStart,&newDMnetwork->eEnd);
1253:   DMPlexGetHeightStratum(newDMnetwork->plex,1,&newDMnetwork->vStart,&newDMnetwork->vEnd);
1254:   newDMnetwork->nEdges    = newDMnetwork->eEnd - newDMnetwork->eStart;
1255:   newDMnetwork->nVertices = newDMnetwork->vEnd - newDMnetwork->vStart;
1256:   newDMnetwork->NVertices = oldDMnetwork->NVertices;
1257:   newDMnetwork->NEdges    = oldDMnetwork->NEdges;

1259:   /* Set Dof section as the section for dm */
1260:   DMSetLocalSection(newDMnetwork->plex,newDMnetwork->DofSection);
1261:   DMGetGlobalSection(newDMnetwork->plex,&newDMnetwork->GlobalDofSection);

1263:   /* Set up subnetwork info in the newDM */
1264:   newDMnetwork->nsubnet  = oldDMnetwork->nsubnet;
1265:   newDMnetwork->ncsubnet = oldDMnetwork->ncsubnet;
1266:   PetscCalloc1(newDMnetwork->nsubnet,&newDMnetwork->subnet);
1267:   /* Copy over the global number of vertices and edges in each subnetwork. Note that these are already
1268:      calculated in DMNetworkLayoutSetUp()
1269:   */
1270:   for(j=0; j < newDMnetwork->nsubnet; j++) {
1271:     newDMnetwork->subnet[j].Nvtx  = oldDMnetwork->subnet[j].Nvtx;
1272:     newDMnetwork->subnet[j].Nedge = oldDMnetwork->subnet[j].Nedge;
1273:   }

1275:   for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++ ) {
1276:     PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);
1277:     header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1278:     newDMnetwork->subnet[header->subnetid].nedge++;
1279:   }

1281:   for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++ ) {
1282:     PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);
1283:     header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1284:     newDMnetwork->subnet[header->subnetid].nvtx++;
1285:   }

1287:   /* Now create the vertices and edge arrays for the subnetworks */
1288:   PetscCalloc1(newDMnetwork->vEnd-newDMnetwork->vStart,&newDMnetwork->subnetvtx);
1289:   subnetvtx = newDMnetwork->subnetvtx;

1291:   for (j=0; j<newDMnetwork->nsubnet; j++) {
1292:     PetscCalloc1(newDMnetwork->subnet[j].nedge,&newDMnetwork->subnet[j].edges);
1293:     newDMnetwork->subnet[j].vertices = subnetvtx;
1294:     subnetvtx                       += newDMnetwork->subnet[j].nvtx;

1296:     /* Temporarily setting nvtx and nedge to 0 so we can use them as counters in the below for loop.
1297:        These get updated when the vertices and edges are added. */
1298:     newDMnetwork->subnet[j].nvtx = newDMnetwork->subnet[j].nedge = 0;
1299:   }

1301:   /* Set the vertices and edges in each subnetwork */
1302:   for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++ ) {
1303:     PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);
1304:     header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1305:     newDMnetwork->subnet[header->subnetid].edges[newDMnetwork->subnet[header->subnetid].nedge++] = e;
1306:   }

1308:   for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++ ) {
1309:     PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);
1310:     header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1311:     newDMnetwork->subnet[header->subnetid].vertices[newDMnetwork->subnet[header->subnetid].nvtx++] = v;
1312:   }

1314:   newDM->setupcalled = (*dm)->setupcalled;
1315:   newDMnetwork->distributecalled = PETSC_TRUE;

1317:   /* Destroy point SF */
1318:   PetscSFDestroy(&pointsf);

1320:   DMDestroy(dm);
1321:   *dm  = newDM;
1322:   return(0);
1323: }

1325: /*@C
1326:   PetscSFGetSubSF - Returns an SF for a specific subset of points. Leaves are re-numbered to reflect the new ordering.

1328:   Input Parameters:
1329: + masterSF - the original SF structure
1330: - map      - a ISLocalToGlobal mapping that contains the subset of points

1332:   Output Parameters:
1333: . subSF    - a subset of the masterSF for the desired subset.
1334: */

1336: PetscErrorCode PetscSFGetSubSF(PetscSF mastersf, ISLocalToGlobalMapping map, PetscSF *subSF) {

1338:   PetscErrorCode        ierr;
1339:   PetscInt              nroots, nleaves, *ilocal_sub;
1340:   PetscInt              i, *ilocal_map, nroots_sub, nleaves_sub = 0;
1341:   PetscInt              *local_points, *remote_points;
1342:   PetscSFNode           *iremote_sub;
1343:   const PetscInt        *ilocal;
1344:   const PetscSFNode     *iremote;

1347:   PetscSFGetGraph(mastersf,&nroots,&nleaves,&ilocal,&iremote);

1349:   /* Look for leaves that pertain to the subset of points. Get the local ordering */
1350:   PetscMalloc1(nleaves,&ilocal_map);
1351:   ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nleaves,ilocal,NULL,ilocal_map);
1352:   for (i = 0; i < nleaves; i++) {
1353:     if (ilocal_map[i] != -1) nleaves_sub += 1;
1354:   }
1355:   /* Re-number ilocal with subset numbering. Need information from roots */
1356:   PetscMalloc2(nroots,&local_points,nroots,&remote_points);
1357:   for (i = 0; i < nroots; i++) local_points[i] = i;
1358:   ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nroots,local_points,NULL,local_points);
1359:   PetscSFBcastBegin(mastersf, MPIU_INT, local_points, remote_points);
1360:   PetscSFBcastEnd(mastersf, MPIU_INT, local_points, remote_points);
1361:   /* Fill up graph using local (that is, local to the subset) numbering. */
1362:   PetscMalloc1(nleaves_sub,&ilocal_sub);
1363:   PetscMalloc1(nleaves_sub,&iremote_sub);
1364:   nleaves_sub = 0;
1365:   for (i = 0; i < nleaves; i++) {
1366:     if (ilocal_map[i] != -1) {
1367:       ilocal_sub[nleaves_sub] = ilocal_map[i];
1368:       iremote_sub[nleaves_sub].rank = iremote[i].rank;
1369:       iremote_sub[nleaves_sub].index = remote_points[ilocal[i]];
1370:       nleaves_sub += 1;
1371:     }
1372:   }
1373:   PetscFree2(local_points,remote_points);
1374:   ISLocalToGlobalMappingGetSize(map,&nroots_sub);

1376:   /* Create new subSF */
1377:   PetscSFCreate(PETSC_COMM_WORLD,subSF);
1378:   PetscSFSetFromOptions(*subSF);
1379:   PetscSFSetGraph(*subSF,nroots_sub,nleaves_sub,ilocal_sub,PETSC_OWN_POINTER,iremote_sub,PETSC_COPY_VALUES);
1380:   PetscFree(ilocal_map);
1381:   PetscFree(iremote_sub);
1382:   return(0);
1383: }

1385: /*@C
1386:   DMNetworkGetSupportingEdges - Return the supporting edges for this vertex point

1388:   Not Collective

1390:   Input Parameters:
1391: + dm - The DMNetwork object
1392: - p  - the vertex point

1394:   Output Parameters:
1395: + nedges - number of edges connected to this vertex point
1396: - edges  - List of edge points

1398:   Level: beginner

1400:   Fortran Notes:
1401:   Since it returns an array, this routine is only available in Fortran 90, and you must
1402:   include petsc.h90 in your code.

1404: .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices
1405: @*/
1406: PetscErrorCode DMNetworkGetSupportingEdges(DM dm,PetscInt vertex,PetscInt *nedges,const PetscInt *edges[])
1407: {
1409:   DM_Network     *network = (DM_Network*)dm->data;

1412:   DMPlexGetSupportSize(network->plex,vertex,nedges);
1413:   DMPlexGetSupport(network->plex,vertex,edges);
1414:   return(0);
1415: }

1417: /*@C
1418:   DMNetworkGetConnectedVertices - Return the connected vertices for this edge point

1420:   Not Collective

1422:   Input Parameters:
1423: + dm - The DMNetwork object
1424: - p  - the edge point

1426:   Output Parameters:
1427: . vertices  - vertices connected to this edge

1429:   Level: beginner

1431:   Fortran Notes:
1432:   Since it returns an array, this routine is only available in Fortran 90, and you must
1433:   include petsc.h90 in your code.

1435: .seealso: DMNetworkCreate, DMNetworkGetSupportingEdges
1436: @*/
1437: PetscErrorCode DMNetworkGetConnectedVertices(DM dm,PetscInt edge,const PetscInt *vertices[])
1438: {
1440:   DM_Network     *network = (DM_Network*)dm->data;

1443:   DMPlexGetCone(network->plex,edge,vertices);
1444:   return(0);
1445: }

1447: /*@
1448:   DMNetworkIsGhostVertex - Returns TRUE if the vertex is a ghost vertex

1450:   Not Collective

1452:   Input Parameters:
1453: + dm - The DMNetwork object
1454: - p  - the vertex point

1456:   Output Parameter:
1457: . isghost - TRUE if the vertex is a ghost point

1459:   Level: beginner

1461: .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices, DMNetworkGetVertexRange
1462: @*/
1463: PetscErrorCode DMNetworkIsGhostVertex(DM dm,PetscInt p,PetscBool *isghost)
1464: {
1466:   DM_Network     *network = (DM_Network*)dm->data;
1467:   PetscInt       offsetg;
1468:   PetscSection   sectiong;

1471:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
1472:   *isghost = PETSC_FALSE;
1473:   DMGetGlobalSection(network->plex,&sectiong);
1474:   PetscSectionGetOffset(sectiong,p,&offsetg);
1475:   if (offsetg < 0) *isghost = PETSC_TRUE;
1476:   return(0);
1477: }

1479: PetscErrorCode DMSetUp_Network(DM dm)
1480: {
1482:   DM_Network     *network=(DM_Network*)dm->data;

1485:   DMNetworkComponentSetUp(dm);
1486:   DMNetworkVariablesSetUp(dm);

1488:   DMSetLocalSection(network->plex,network->DofSection);
1489:   DMGetGlobalSection(network->plex,&network->GlobalDofSection);

1491:   dm->setupcalled = PETSC_TRUE;
1492:   DMViewFromOptions(dm,NULL,"-dm_view");
1493:   return(0);
1494: }

1496: /*@
1497:     DMNetworkHasJacobian - Sets global flag for using user's sub Jacobian matrices
1498:                             -- replaced by DMNetworkSetOption(network,userjacobian,PETSC_TURE)?

1500:     Collective

1502:     Input Parameters:
1503: +   dm - The DMNetwork object
1504: .   eflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for edges
1505: -   vflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for vertices

1507:     Level: intermediate

1509: @*/
1510: PetscErrorCode DMNetworkHasJacobian(DM dm,PetscBool eflg,PetscBool vflg)
1511: {
1512:   DM_Network     *network=(DM_Network*)dm->data;
1514:   PetscInt       nVertices = network->nVertices;

1517:   network->userEdgeJacobian   = eflg;
1518:   network->userVertexJacobian = vflg;

1520:   if (eflg && !network->Je) {
1521:     PetscCalloc1(3*network->nEdges,&network->Je);
1522:   }

1524:   if (vflg && !network->Jv && nVertices) {
1525:     PetscInt       i,*vptr,nedges,vStart=network->vStart;
1526:     PetscInt       nedges_total;
1527:     const PetscInt *edges;

1529:     /* count nvertex_total */
1530:     nedges_total = 0;
1531:     PetscMalloc1(nVertices+1,&vptr);

1533:     vptr[0] = 0;
1534:     for (i=0; i<nVertices; i++) {
1535:       DMNetworkGetSupportingEdges(dm,i+vStart,&nedges,&edges);
1536:       nedges_total += nedges;
1537:       vptr[i+1] = vptr[i] + 2*nedges + 1;
1538:     }

1540:     PetscCalloc1(2*nedges_total+nVertices,&network->Jv);
1541:     network->Jvptr = vptr;
1542:   }
1543:   return(0);
1544: }

1546: /*@
1547:     DMNetworkEdgeSetMatrix - Sets user-provided Jacobian matrices for this edge to the network

1549:     Not Collective

1551:     Input Parameters:
1552: +   dm - The DMNetwork object
1553: .   p  - the edge point
1554: -   J - array (size = 3) of Jacobian submatrices for this edge point:
1555:         J[0]: this edge
1556:         J[1] and J[2]: connected vertices, obtained by calling DMNetworkGetConnectedVertices()

1558:     Level: advanced

1560: .seealso: DMNetworkVertexSetMatrix
1561: @*/
1562: PetscErrorCode DMNetworkEdgeSetMatrix(DM dm,PetscInt p,Mat J[])
1563: {
1564:   DM_Network     *network=(DM_Network*)dm->data;

1567:   if (!network->Je) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkEdgeSetMatrix");

1569:   if (J) {
1570:     network->Je[3*p]   = J[0];
1571:     network->Je[3*p+1] = J[1];
1572:     network->Je[3*p+2] = J[2];
1573:   }
1574:   return(0);
1575: }

1577: /*@
1578:     DMNetworkVertexSetMatrix - Sets user-provided Jacobian matrix for this vertex to the network

1580:     Not Collective

1582:     Input Parameters:
1583: +   dm - The DMNetwork object
1584: .   p  - the vertex point
1585: -   J - array of Jacobian (size = 2*(num of supporting edges) + 1) submatrices for this vertex point:
1586:         J[0]:       this vertex
1587:         J[1+2*i]:   i-th supporting edge
1588:         J[1+2*i+1]: i-th connected vertex

1590:     Level: advanced

1592: .seealso: DMNetworkEdgeSetMatrix
1593: @*/
1594: PetscErrorCode DMNetworkVertexSetMatrix(DM dm,PetscInt p,Mat J[])
1595: {
1597:   DM_Network     *network=(DM_Network*)dm->data;
1598:   PetscInt       i,*vptr,nedges,vStart=network->vStart;
1599:   const PetscInt *edges;

1602:   if (!network->Jv) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkVertexSetMatrix");

1604:   if (J) {
1605:     vptr = network->Jvptr;
1606:     network->Jv[vptr[p-vStart]] = J[0]; /* Set Jacobian for this vertex */

1608:     /* Set Jacobian for each supporting edge and connected vertex */
1609:     DMNetworkGetSupportingEdges(dm,p,&nedges,&edges);
1610:     for (i=1; i<=2*nedges; i++) network->Jv[vptr[p-vStart]+i] = J[i];
1611:   }
1612:   return(0);
1613: }

1615: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1616: {
1618:   PetscInt       j;
1619:   PetscScalar    val=(PetscScalar)ncols;

1622:   if (!ghost) {
1623:     for (j=0; j<nrows; j++) {
1624:       VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
1625:     }
1626:   } else {
1627:     for (j=0; j<nrows; j++) {
1628:       VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
1629:     }
1630:   }
1631:   return(0);
1632: }

1634: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1635: {
1637:   PetscInt       j,ncols_u;
1638:   PetscScalar    val;

1641:   if (!ghost) {
1642:     for (j=0; j<nrows; j++) {
1643:       MatGetRow(Ju,j,&ncols_u,NULL,NULL);
1644:       val = (PetscScalar)ncols_u;
1645:       VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
1646:       MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
1647:     }
1648:   } else {
1649:     for (j=0; j<nrows; j++) {
1650:       MatGetRow(Ju,j,&ncols_u,NULL,NULL);
1651:       val = (PetscScalar)ncols_u;
1652:       VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
1653:       MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
1654:     }
1655:   }
1656:   return(0);
1657: }

1659: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1660: {

1664:   if (Ju) {
1665:     MatSetPreallocationUserblock_private(Ju,nrows,rows,ncols,ghost,vdnz,vonz);
1666:   } else {
1667:     MatSetPreallocationDenseblock_private(nrows,rows,ncols,ghost,vdnz,vonz);
1668:   }
1669:   return(0);
1670: }

1672: PETSC_STATIC_INLINE PetscErrorCode MatSetDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1673: {
1675:   PetscInt       j,*cols;
1676:   PetscScalar    *zeros;

1679:   PetscCalloc2(ncols,&cols,nrows*ncols,&zeros);
1680:   for (j=0; j<ncols; j++) cols[j] = j+ cstart;
1681:   MatSetValues(*J,nrows,rows,ncols,cols,zeros,INSERT_VALUES);
1682:   PetscFree2(cols,zeros);
1683:   return(0);
1684: }

1686: PETSC_STATIC_INLINE PetscErrorCode MatSetUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1687: {
1689:   PetscInt       j,M,N,row,col,ncols_u;
1690:   const PetscInt *cols;
1691:   PetscScalar    zero=0.0;

1694:   MatGetSize(Ju,&M,&N);
1695:   if (nrows != M || ncols != N) SETERRQ4(PetscObjectComm((PetscObject)Ju),PETSC_ERR_USER,"%D by %D must equal %D by %D",nrows,ncols,M,N);

1697:   for (row=0; row<nrows; row++) {
1698:     MatGetRow(Ju,row,&ncols_u,&cols,NULL);
1699:     for (j=0; j<ncols_u; j++) {
1700:       col = cols[j] + cstart;
1701:       MatSetValues(*J,1,&rows[row],1,&col,&zero,INSERT_VALUES);
1702:     }
1703:     MatRestoreRow(Ju,row,&ncols_u,&cols,NULL);
1704:   }
1705:   return(0);
1706: }

1708: PETSC_STATIC_INLINE PetscErrorCode MatSetblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1709: {

1713:   if (Ju) {
1714:     MatSetUserblock_private(Ju,nrows,rows,ncols,cstart,J);
1715:   } else {
1716:     MatSetDenseblock_private(nrows,rows,ncols,cstart,J);
1717:   }
1718:   return(0);
1719: }

1721: /* Creates a GlobalToLocal mapping with a Local and Global section. This is akin to the routine DMGetLocalToGlobalMapping but without the need of providing a dm.
1722: */
1723: PetscErrorCode CreateSubGlobalToLocalMapping_private(PetscSection globalsec, PetscSection localsec, ISLocalToGlobalMapping *ltog)
1724: {
1726:   PetscInt       i,size,dof;
1727:   PetscInt       *glob2loc;

1730:   PetscSectionGetStorageSize(localsec,&size);
1731:   PetscMalloc1(size,&glob2loc);

1733:   for (i = 0; i < size; i++) {
1734:     PetscSectionGetOffset(globalsec,i,&dof);
1735:     dof = (dof >= 0) ? dof : -(dof + 1);
1736:     glob2loc[i] = dof;
1737:   }

1739:   ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,size,glob2loc,PETSC_OWN_POINTER,ltog);
1740: #if 0
1741:   PetscIntView(size,glob2loc,PETSC_VIEWER_STDOUT_WORLD);
1742: #endif
1743:   return(0);
1744: }

1746:  #include <petsc/private/matimpl.h>

1748: PetscErrorCode DMCreateMatrix_Network_Nest(DM dm,Mat *J)
1749: {
1751:   DM_Network     *network = (DM_Network*)dm->data;
1752:   PetscMPIInt    rank, size;
1753:   PetscInt       eDof,vDof;
1754:   Mat            j11,j12,j21,j22,bA[2][2];
1755:   MPI_Comm       comm;
1756:   ISLocalToGlobalMapping eISMap,vISMap;

1759:   PetscObjectGetComm((PetscObject)dm,&comm);
1760:   MPI_Comm_rank(comm,&rank);
1761:   MPI_Comm_size(comm,&size);

1763:   PetscSectionGetConstrainedStorageSize(network->edge.GlobalDofSection,&eDof);
1764:   PetscSectionGetConstrainedStorageSize(network->vertex.GlobalDofSection,&vDof);

1766:   MatCreate(comm, &j11);
1767:   MatSetSizes(j11, eDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
1768:   MatSetType(j11, MATMPIAIJ);

1770:   MatCreate(comm, &j12);
1771:   MatSetSizes(j12, eDof, vDof, PETSC_DETERMINE ,PETSC_DETERMINE);
1772:   MatSetType(j12, MATMPIAIJ);

1774:   MatCreate(comm, &j21);
1775:   MatSetSizes(j21, vDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
1776:   MatSetType(j21, MATMPIAIJ);

1778:   MatCreate(comm, &j22);
1779:   MatSetSizes(j22, vDof, vDof, PETSC_DETERMINE, PETSC_DETERMINE);
1780:   MatSetType(j22, MATMPIAIJ);

1782:   bA[0][0] = j11;
1783:   bA[0][1] = j12;
1784:   bA[1][0] = j21;
1785:   bA[1][1] = j22;

1787:   CreateSubGlobalToLocalMapping_private(network->edge.GlobalDofSection,network->edge.DofSection,&eISMap);
1788:   CreateSubGlobalToLocalMapping_private(network->vertex.GlobalDofSection,network->vertex.DofSection,&vISMap);

1790:   MatSetLocalToGlobalMapping(j11,eISMap,eISMap);
1791:   MatSetLocalToGlobalMapping(j12,eISMap,vISMap);
1792:   MatSetLocalToGlobalMapping(j21,vISMap,eISMap);
1793:   MatSetLocalToGlobalMapping(j22,vISMap,vISMap);

1795:   MatSetUp(j11);
1796:   MatSetUp(j12);
1797:   MatSetUp(j21);
1798:   MatSetUp(j22);

1800:   MatCreateNest(comm,2,NULL,2,NULL,&bA[0][0],J);
1801:   MatSetUp(*J);
1802:   MatNestSetVecType(*J,VECNEST);
1803:   MatDestroy(&j11);
1804:   MatDestroy(&j12);
1805:   MatDestroy(&j21);
1806:   MatDestroy(&j22);

1808:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1809:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1810:   MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);

1812:   /* Free structures */
1813:   ISLocalToGlobalMappingDestroy(&eISMap);
1814:   ISLocalToGlobalMappingDestroy(&vISMap);
1815:   return(0);
1816: }

1818: PetscErrorCode DMCreateMatrix_Network(DM dm,Mat *J)
1819: {
1821:   DM_Network     *network = (DM_Network*)dm->data;
1822:   PetscInt       eStart,eEnd,vStart,vEnd,rstart,nrows,*rows,localSize;
1823:   PetscInt       cstart,ncols,j,e,v;
1824:   PetscBool      ghost,ghost_vc,ghost2,isNest;
1825:   Mat            Juser;
1826:   PetscSection   sectionGlobal;
1827:   PetscInt       nedges,*vptr=NULL,vc,*rows_v; /* suppress maybe-uninitialized warning */
1828:   const PetscInt *edges,*cone;
1829:   MPI_Comm       comm;
1830:   MatType        mtype;
1831:   Vec            vd_nz,vo_nz;
1832:   PetscInt       *dnnz,*onnz;
1833:   PetscScalar    *vdnz,*vonz;

1836:   mtype = dm->mattype;
1837:   PetscStrcmp(mtype,MATNEST,&isNest);
1838:   if (isNest) {
1839:     DMCreateMatrix_Network_Nest(dm,J);
1840:     MatSetDM(*J,dm);
1841:     return(0);
1842:   }

1844:   if (!network->userEdgeJacobian && !network->userVertexJacobian) {
1845:     /* user does not provide Jacobian blocks */
1846:     DMCreateMatrix_Plex(network->plex,J);
1847:     MatSetDM(*J,dm);
1848:     return(0);
1849:   }

1851:   MatCreate(PetscObjectComm((PetscObject)dm),J);
1852:   DMGetGlobalSection(network->plex,&sectionGlobal);
1853:   PetscSectionGetConstrainedStorageSize(sectionGlobal,&localSize);
1854:   MatSetSizes(*J,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);

1856:   MatSetType(*J,MATAIJ);
1857:   MatSetFromOptions(*J);

1859:   /* (1) Set matrix preallocation */
1860:   /*------------------------------*/
1861:   PetscObjectGetComm((PetscObject)dm,&comm);
1862:   VecCreate(comm,&vd_nz);
1863:   VecSetSizes(vd_nz,localSize,PETSC_DECIDE);
1864:   VecSetFromOptions(vd_nz);
1865:   VecSet(vd_nz,0.0);
1866:   VecDuplicate(vd_nz,&vo_nz);

1868:   /* Set preallocation for edges */
1869:   /*-----------------------------*/
1870:   DMNetworkGetEdgeRange(dm,&eStart,&eEnd);

1872:   PetscMalloc1(localSize,&rows);
1873:   for (e=eStart; e<eEnd; e++) {
1874:     /* Get row indices */
1875:     DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
1876:     DMNetworkGetNumVariables(dm,e,&nrows);
1877:     if (nrows) {
1878:       for (j=0; j<nrows; j++) rows[j] = j + rstart;

1880:       /* Set preallocation for conntected vertices */
1881:       DMNetworkGetConnectedVertices(dm,e,&cone);
1882:       for (v=0; v<2; v++) {
1883:         DMNetworkGetNumVariables(dm,cone[v],&ncols);

1885:         if (network->Je) {
1886:           Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
1887:         } else Juser = NULL;
1888:         DMNetworkIsGhostVertex(dm,cone[v],&ghost);
1889:         MatSetPreallocationblock_private(Juser,nrows,rows,ncols,ghost,vd_nz,vo_nz);
1890:       }

1892:       /* Set preallocation for edge self */
1893:       cstart = rstart;
1894:       if (network->Je) {
1895:         Juser = network->Je[3*e]; /* Jacobian(e,e) */
1896:       } else Juser = NULL;
1897:       MatSetPreallocationblock_private(Juser,nrows,rows,nrows,PETSC_FALSE,vd_nz,vo_nz);
1898:     }
1899:   }

1901:   /* Set preallocation for vertices */
1902:   /*--------------------------------*/
1903:   DMNetworkGetVertexRange(dm,&vStart,&vEnd);
1904:   if (vEnd - vStart) vptr = network->Jvptr;

1906:   for (v=vStart; v<vEnd; v++) {
1907:     /* Get row indices */
1908:     DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
1909:     DMNetworkGetNumVariables(dm,v,&nrows);
1910:     if (!nrows) continue;

1912:     DMNetworkIsGhostVertex(dm,v,&ghost);
1913:     if (ghost) {
1914:       PetscMalloc1(nrows,&rows_v);
1915:     } else {
1916:       rows_v = rows;
1917:     }

1919:     for (j=0; j<nrows; j++) rows_v[j] = j + rstart;

1921:     /* Get supporting edges and connected vertices */
1922:     DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);

1924:     for (e=0; e<nedges; e++) {
1925:       /* Supporting edges */
1926:       DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
1927:       DMNetworkGetNumVariables(dm,edges[e],&ncols);

1929:       if (network->Jv) {
1930:         Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */
1931:       } else Juser = NULL;
1932:       MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost,vd_nz,vo_nz);

1934:       /* Connected vertices */
1935:       DMNetworkGetConnectedVertices(dm,edges[e],&cone);
1936:       vc = (v == cone[0]) ? cone[1]:cone[0];
1937:       DMNetworkIsGhostVertex(dm,vc,&ghost_vc);

1939:       DMNetworkGetNumVariables(dm,vc,&ncols);

1941:       if (network->Jv) {
1942:         Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
1943:       } else Juser = NULL;
1944:       if (ghost_vc||ghost) {
1945:         ghost2 = PETSC_TRUE;
1946:       } else {
1947:         ghost2 = PETSC_FALSE;
1948:       }
1949:       MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost2,vd_nz,vo_nz);
1950:     }

1952:     /* Set preallocation for vertex self */
1953:     DMNetworkIsGhostVertex(dm,v,&ghost);
1954:     if (!ghost) {
1955:       DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
1956:       if (network->Jv) {
1957:         Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
1958:       } else Juser = NULL;
1959:       MatSetPreallocationblock_private(Juser,nrows,rows_v,nrows,PETSC_FALSE,vd_nz,vo_nz);
1960:     }
1961:     if (ghost) {
1962:       PetscFree(rows_v);
1963:     }
1964:   }

1966:   VecAssemblyBegin(vd_nz);
1967:   VecAssemblyBegin(vo_nz);

1969:   PetscMalloc2(localSize,&dnnz,localSize,&onnz);

1971:   VecAssemblyEnd(vd_nz);
1972:   VecAssemblyEnd(vo_nz);

1974:   VecGetArray(vd_nz,&vdnz);
1975:   VecGetArray(vo_nz,&vonz);
1976:   for (j=0; j<localSize; j++) {
1977:     dnnz[j] = (PetscInt)PetscRealPart(vdnz[j]);
1978:     onnz[j] = (PetscInt)PetscRealPart(vonz[j]);
1979:   }
1980:   VecRestoreArray(vd_nz,&vdnz);
1981:   VecRestoreArray(vo_nz,&vonz);
1982:   VecDestroy(&vd_nz);
1983:   VecDestroy(&vo_nz);

1985:   MatSeqAIJSetPreallocation(*J,0,dnnz);
1986:   MatMPIAIJSetPreallocation(*J,0,dnnz,0,onnz);
1987:   MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);

1989:   PetscFree2(dnnz,onnz);

1991:   /* (2) Set matrix entries for edges */
1992:   /*----------------------------------*/
1993:   for (e=eStart; e<eEnd; e++) {
1994:     /* Get row indices */
1995:     DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
1996:     DMNetworkGetNumVariables(dm,e,&nrows);
1997:     if (nrows) {
1998:       for (j=0; j<nrows; j++) rows[j] = j + rstart;

2000:       /* Set matrix entries for conntected vertices */
2001:       DMNetworkGetConnectedVertices(dm,e,&cone);
2002:       for (v=0; v<2; v++) {
2003:         DMNetworkGetVariableGlobalOffset(dm,cone[v],&cstart);
2004:         DMNetworkGetNumVariables(dm,cone[v],&ncols);

2006:         if (network->Je) {
2007:           Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
2008:         } else Juser = NULL;
2009:         MatSetblock_private(Juser,nrows,rows,ncols,cstart,J);
2010:       }

2012:       /* Set matrix entries for edge self */
2013:       cstart = rstart;
2014:       if (network->Je) {
2015:         Juser = network->Je[3*e]; /* Jacobian(e,e) */
2016:       } else Juser = NULL;
2017:       MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);
2018:     }
2019:   }

2021:   /* Set matrix entries for vertices */
2022:   /*---------------------------------*/
2023:   for (v=vStart; v<vEnd; v++) {
2024:     /* Get row indices */
2025:     DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
2026:     DMNetworkGetNumVariables(dm,v,&nrows);
2027:     if (!nrows) continue;

2029:     DMNetworkIsGhostVertex(dm,v,&ghost);
2030:     if (ghost) {
2031:       PetscMalloc1(nrows,&rows_v);
2032:     } else {
2033:       rows_v = rows;
2034:     }
2035:     for (j=0; j<nrows; j++) rows_v[j] = j + rstart;

2037:     /* Get supporting edges and connected vertices */
2038:     DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);

2040:     for (e=0; e<nedges; e++) {
2041:       /* Supporting edges */
2042:       DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
2043:       DMNetworkGetNumVariables(dm,edges[e],&ncols);

2045:       if (network->Jv) {
2046:         Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */
2047:       } else Juser = NULL;
2048:       MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);

2050:       /* Connected vertices */
2051:       DMNetworkGetConnectedVertices(dm,edges[e],&cone);
2052:       vc = (v == cone[0]) ? cone[1]:cone[0];

2054:       DMNetworkGetVariableGlobalOffset(dm,vc,&cstart);
2055:       DMNetworkGetNumVariables(dm,vc,&ncols);

2057:       if (network->Jv) {
2058:         Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
2059:       } else Juser = NULL;
2060:       MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);
2061:     }

2063:     /* Set matrix entries for vertex self */
2064:     if (!ghost) {
2065:       DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
2066:       if (network->Jv) {
2067:         Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
2068:       } else Juser = NULL;
2069:       MatSetblock_private(Juser,nrows,rows_v,nrows,cstart,J);
2070:     }
2071:     if (ghost) {
2072:       PetscFree(rows_v);
2073:     }
2074:   }
2075:   PetscFree(rows);

2077:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
2078:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);

2080:   MatSetDM(*J,dm);
2081:   return(0);
2082: }

2084: PetscErrorCode DMDestroy_Network(DM dm)
2085: {
2087:   DM_Network     *network = (DM_Network*)dm->data;
2088:   PetscInt       j;

2091:   if (--network->refct > 0) return(0);
2092:   if (network->Je) {
2093:     PetscFree(network->Je);
2094:   }
2095:   if (network->Jv) {
2096:     PetscFree(network->Jvptr);
2097:     PetscFree(network->Jv);
2098:   }

2100:   ISLocalToGlobalMappingDestroy(&network->vertex.mapping);
2101:   PetscSectionDestroy(&network->vertex.DofSection);
2102:   PetscSectionDestroy(&network->vertex.GlobalDofSection);
2103:   if (network->vltog) {
2104:     PetscFree(network->vltog);
2105:   }
2106:   if (network->vertex.sf) {
2107:     PetscSFDestroy(&network->vertex.sf);
2108:   }
2109:   /* edge */
2110:   ISLocalToGlobalMappingDestroy(&network->edge.mapping);
2111:   PetscSectionDestroy(&network->edge.DofSection);
2112:   PetscSectionDestroy(&network->edge.GlobalDofSection);
2113:   if (network->edge.sf) {
2114:     PetscSFDestroy(&network->edge.sf);
2115:   }
2116:   DMDestroy(&network->plex);
2117:   PetscSectionDestroy(&network->DataSection);
2118:   PetscSectionDestroy(&network->DofSection);

2120:   for(j=0; j<network->nsubnet; j++) {
2121:     PetscFree(network->subnet[j].edges);
2122:   }
2123:   PetscFree(network->subnetvtx);

2125:   PetscFree(network->subnet);
2126:   PetscFree(network->componentdataarray);
2127:   PetscFree2(network->header,network->cvalue);
2128:   PetscFree(network);
2129:   return(0);
2130: }

2132: PetscErrorCode DMView_Network(DM dm,PetscViewer viewer)
2133: {
2135:   DM_Network     *network = (DM_Network*)dm->data;
2136:   PetscBool      iascii;
2137:   PetscMPIInt    rank;
2138:   PetscInt       p,nsubnet;

2141:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
2142:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
2145:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2146:   if (iascii) {
2147:     const PetscInt    *cone,*vtx,*edges;
2148:     PetscInt          vfrom,vto,i,j,nv,ne;

2150:     nsubnet = network->nsubnet - network->ncsubnet; /* num of subnetworks */
2151:     PetscViewerASCIIPushSynchronized(viewer);
2152:     PetscViewerASCIISynchronizedPrintf(viewer, "  [%d] nsubnet: %D; nsubnetCouple: %D; nEdges: %D; nVertices: %D\n",rank,nsubnet,network->ncsubnet,network->nEdges,network->nVertices);

2154:     for (i=0; i<nsubnet; i++) {
2155:       DMNetworkGetSubnetworkInfo(dm,i,&nv,&ne,&vtx,&edges);
2156:       if (ne) {
2157:         PetscViewerASCIISynchronizedPrintf(viewer, "     Subnet %D: nEdges %D, nVertices %D\n",i,ne,nv);
2158:         for (j=0; j<ne; j++) {
2159:           p = edges[j];
2160:           DMNetworkGetConnectedVertices(dm,p,&cone);
2161:           DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);
2162:           DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);
2163:           DMNetworkGetGlobalEdgeIndex(dm,edges[j],&p);
2164:           PetscViewerASCIISynchronizedPrintf(viewer, "       edge %D: %D----> %D\n",p,vfrom,vto);
2165:         }
2166:       }
2167:     }
2168:     /* Coupling subnets */
2169:     nsubnet = network->nsubnet;
2170:     for (; i<nsubnet; i++) {
2171:       DMNetworkGetSubnetworkInfo(dm,i,&nv,&ne,&vtx,&edges);
2172:       if (ne) {
2173:         PetscViewerASCIISynchronizedPrintf(viewer, "     Subnet %D (couple): nEdges %D, nVertices %D\n",i,ne,nv);
2174:         for (j=0; j<ne; j++) {
2175:           p = edges[j];
2176:           DMNetworkGetConnectedVertices(dm,p,&cone);
2177:           DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);
2178:           DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);
2179:           PetscViewerASCIISynchronizedPrintf(viewer, "       edge %D: %D----> %D\n",p,vfrom,vto);
2180:         }
2181:       }
2182:     }
2183:     PetscViewerFlush(viewer);
2184:     PetscViewerASCIIPopSynchronized(viewer);
2185:   } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMNetwork writing", ((PetscObject)viewer)->type_name);
2186:   return(0);
2187: }

2189: PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l)
2190: {
2192:   DM_Network     *network = (DM_Network*)dm->data;

2195:   DMGlobalToLocalBegin(network->plex,g,mode,l);
2196:   return(0);
2197: }

2199: PetscErrorCode DMGlobalToLocalEnd_Network(DM dm, Vec g, InsertMode mode, Vec l)
2200: {
2202:   DM_Network     *network = (DM_Network*)dm->data;

2205:   DMGlobalToLocalEnd(network->plex,g,mode,l);
2206:   return(0);
2207: }

2209: PetscErrorCode DMLocalToGlobalBegin_Network(DM dm, Vec l, InsertMode mode, Vec g)
2210: {
2212:   DM_Network     *network = (DM_Network*)dm->data;

2215:   DMLocalToGlobalBegin(network->plex,l,mode,g);
2216:   return(0);
2217: }

2219: PetscErrorCode DMLocalToGlobalEnd_Network(DM dm, Vec l, InsertMode mode, Vec g)
2220: {
2222:   DM_Network     *network = (DM_Network*)dm->data;

2225:   DMLocalToGlobalEnd(network->plex,l,mode,g);
2226:   return(0);
2227: }

2229: /*@
2230:   DMNetworkGetVertexLocalToGlobalOrdering - Get vertex globle index

2232:   Not collective

2234:   Input Parameters:
2235: + dm - the dm object
2236: - vloc - local vertex ordering, start from 0

2238:   Output Parameters:
2239: .  vg  - global vertex ordering, start from 0

2241:   Level: advanced

2243: .seealso: DMNetworkSetVertexLocalToGlobalOrdering()
2244: @*/
2245: PetscErrorCode DMNetworkGetVertexLocalToGlobalOrdering(DM dm,PetscInt vloc,PetscInt *vg)
2246: {
2247:   DM_Network  *network = (DM_Network*)dm->data;
2248:   PetscInt    *vltog = network->vltog;

2251:   if (!vltog) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkSetVertexLocalToGlobalOrdering() first");
2252:   *vg = vltog[vloc];
2253:   return(0);
2254: }

2256: /*@
2257:   DMNetworkSetVertexLocalToGlobalOrdering - Create and setup vertex local to globle map

2259:   Collective

2261:   Input Parameters:
2262: . dm - the dm object

2264:   Level: advanced

2266: .seealso: DMNetworkGetGlobalVertexIndex()
2267: @*/
2268: PetscErrorCode DMNetworkSetVertexLocalToGlobalOrdering(DM dm)
2269: {
2270:   PetscErrorCode    ierr;
2271:   DM_Network        *network=(DM_Network*)dm->data;
2272:   MPI_Comm          comm;
2273:   PetscMPIInt       rank,size,*displs,*recvcounts,remoterank;
2274:   PetscBool         ghost;
2275:   PetscInt          *vltog,nroots,nleaves,i,*vrange,k,N,lidx;
2276:   const PetscSFNode *iremote;
2277:   PetscSF           vsf;
2278:   Vec               Vleaves,Vleaves_seq;
2279:   VecScatter        ctx;
2280:   PetscScalar       *varr,val;
2281:   const PetscScalar *varr_read;

2284:   PetscObjectGetComm((PetscObject)dm,&comm);
2285:   MPI_Comm_size(comm,&size);
2286:   MPI_Comm_rank(comm,&rank);

2288:   if (size == 1) {
2289:     nroots = network->vEnd - network->vStart;
2290:     PetscMalloc1(nroots, &vltog);
2291:     for (i=0; i<nroots; i++) vltog[i] = i;
2292:     network->vltog = vltog;
2293:     return(0);
2294:   }

2296:   if (!network->distributecalled) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkDistribute() first");
2297:   if (network->vltog) {
2298:     PetscFree(network->vltog);
2299:   }

2301:   DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);
2302:   PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);
2303:   vsf = network->vertex.sf;

2305:   PetscMalloc3(size+1,&vrange,size+1,&displs,size,&recvcounts);
2306:   PetscSFGetGraph(vsf,&nroots,&nleaves,NULL,&iremote);

2308:   for (i=0; i<size; i++) { displs[i] = i; recvcounts[i] = 1;}

2310:   i         = nroots - nleaves; /* local number of vertices, excluding ghosts */
2311:   vrange[0] = 0;
2312:   MPI_Allgatherv(&i,1,MPIU_INT,vrange+1,recvcounts,displs,MPIU_INT,comm);
2313:   for (i=2; i<size+1; i++) {vrange[i] += vrange[i-1];}

2315:   PetscMalloc1(nroots, &vltog);
2316:   network->vltog = vltog;

2318:   /* Set vltog for non-ghost vertices */
2319:   k = 0;
2320:   for (i=0; i<nroots; i++) {
2321:     DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);
2322:     if (ghost) continue;
2323:     vltog[i] = vrange[rank] + k++;
2324:   }
2325:   PetscFree3(vrange,displs,recvcounts);

2327:   /* Set vltog for ghost vertices */
2328:   /* (a) create parallel Vleaves and sequential Vleaves_seq to convert local iremote[*].index to global index */
2329:   VecCreate(comm,&Vleaves);
2330:   VecSetSizes(Vleaves,2*nleaves,PETSC_DETERMINE);
2331:   VecSetFromOptions(Vleaves);
2332:   VecGetArray(Vleaves,&varr);
2333:   for (i=0; i<nleaves; i++) {
2334:     varr[2*i]   = (PetscScalar)(iremote[i].rank);  /* rank of remote process */
2335:     varr[2*i+1] = (PetscScalar)(iremote[i].index); /* local index in remote process */
2336:   }
2337:   VecRestoreArray(Vleaves,&varr);

2339:   /* (b) scatter local info to remote processes via VecScatter() */
2340:   VecScatterCreateToAll(Vleaves,&ctx,&Vleaves_seq);
2341:   VecScatterBegin(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);
2342:   VecScatterEnd(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);

2344:   /* (c) convert local indices to global indices in parallel vector Vleaves */
2345:   VecGetSize(Vleaves_seq,&N);
2346:   VecGetArrayRead(Vleaves_seq,&varr_read);
2347:   for (i=0; i<N; i+=2) {
2348:     remoterank = (PetscMPIInt)PetscRealPart(varr_read[i]);
2349:     if (remoterank == rank) {
2350:       k = i+1; /* row number */
2351:       lidx = (PetscInt)PetscRealPart(varr_read[i+1]);
2352:       val  = (PetscScalar)vltog[lidx]; /* global index for non-ghost vertex computed above */
2353:       VecSetValues(Vleaves,1,&k,&val,INSERT_VALUES);
2354:     }
2355:   }
2356:   VecRestoreArrayRead(Vleaves_seq,&varr_read);
2357:   VecAssemblyBegin(Vleaves);
2358:   VecAssemblyEnd(Vleaves);

2360:   /* (d) Set vltog for ghost vertices by copying local values of Vleaves */
2361:   VecGetArrayRead(Vleaves,&varr_read);
2362:   k = 0;
2363:   for (i=0; i<nroots; i++) {
2364:     DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);
2365:     if (!ghost) continue;
2366:     vltog[i] = (PetscInt)PetscRealPart(varr_read[2*k+1]); k++;
2367:   }
2368:   VecRestoreArrayRead(Vleaves,&varr_read);

2370:   VecDestroy(&Vleaves);
2371:   VecDestroy(&Vleaves_seq);
2372:   VecScatterDestroy(&ctx);
2373:   return(0);
2374: }