Actual source code: network.c

petsc-master 2019-09-20
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:   DMNetworkSetSizes - Sets the number of subnetworks, local and global vertices and edges for each subnetwork.

 28:   Collective on dm

 30:   Input Parameters:
 31: + dm - the dm object
 32: . Nsubnet - global number of subnetworks
 33: . nV - number of local vertices for each subnetwork
 34: . nE - number of local edges for each subnetwork
 35: . NsubnetCouple - global number of coupling subnetworks
 36: - nec - number of local edges for each coupling subnetwork

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

 40:    Level: intermediate

 42: .seealso: DMNetworkCreate()
 43: @*/
 44: PetscErrorCode DMNetworkSetSizes(DM dm,PetscInt Nsubnet,PetscInt nV[], PetscInt nE[],PetscInt NsubnetCouple,PetscInt nec[])
 45: {
 47:   DM_Network     *network = (DM_Network*) dm->data;
 48:   PetscInt       a[2],b[2],i;

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

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

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

 61:   network->nsubnet  = Nsubnet + NsubnetCouple;
 62:   network->ncsubnet = NsubnetCouple;
 63:   PetscCalloc1(Nsubnet+NsubnetCouple,&network->subnet);

 65:   /* ----------------------------------------------------------
 66:    p=v or e; P=V or E
 67:    subnet[0].pStart   = 0
 68:    subnet[i+1].pStart = subnet[i].pEnd = subnet[i].pStart + (nE[i] or NV[i])
 69:    ----------------------------------------------------------------------- */
 70:   for (i=0; i < Nsubnet; i++) {
 71:     /* Get global number of vertices and edges for subnet[i] */
 72:     a[0] = nV[i]; a[1] = nE[i]; /* local number of vertices (excluding ghost) and edges */
 73:     MPIU_Allreduce(a,b,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
 74:     network->subnet[i].Nvtx = b[0];
 75:     network->subnet[i].Nedge = b[1];

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

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

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

 86:     network->subnet[i].nedge  = nE[i];
 87:     network->subnet[i].eStart = network->nEdges;
 88:     network->subnet[i].eEnd   = network->subnet[i].eStart + nE[i];
 89:     network->nEdges += nE[i];
 90:     network->NEdges += network->subnet[i].Nedge;
 91:   }

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

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

102:     network->subnet[i].nedge  = nec[i-Nsubnet];
103:     network->subnet[i].eStart = network->nEdges;
104:     network->subnet[i].eEnd = network->subnet[i].eStart + nec[i-Nsubnet];
105:     network->nEdges += nec[i-Nsubnet];
106:     network->NEdges += network->subnet[i].Nedge;
107:   }
108:   return(0);
109: }

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

114:   Logically collective on dm

116:   Input Parameters:
117: + dm - the dm object
118: . edgelist - list of edges for each subnetwork
119: - edgelistCouple - list of edges for each coupling subnetwork

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

125:   Level: intermediate

127:   Example usage:
128:   Consider the following 2 separate networks and a coupling network:

130: .vb
131:  network 0: v0 -> v1 -> v2 -> v3
132:  network 1: v1 -> v2 -> v0
133:  coupling network: network 1: v2 -> network 0: v0
134: .ve

136:  The resulting input
137:    edgelist[0] = [0 1 | 1 2 | 2 3];
138:    edgelist[1] = [1 2 | 2 0]
139:    edgelistCouple[0] = [(network)1 (v)2 (network)0 (v)0].

141: .seealso: DMNetworkCreate, DMNetworkSetSizes
142: @*/
143: PetscErrorCode DMNetworkSetEdgeList(DM dm,PetscInt *edgelist[],PetscInt *edgelistCouple[])
144: {
145:   DM_Network *network = (DM_Network*) dm->data;
146:   PetscInt   i,nsubnet,ncsubnet=network->ncsubnet;

149:   nsubnet = network->nsubnet - ncsubnet;
150:   for(i=0; i < nsubnet; i++) {
151:     network->subnet[i].edgelist = edgelist[i];
152:   }
153:   if (edgelistCouple) {
154:     PetscInt j;
155:     j = 0;
156:     nsubnet = network->nsubnet;
157:     while (i < nsubnet) network->subnet[i++].edgelist = edgelistCouple[j++];
158:   }
159:   return(0);
160: }

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

165:   Collective on dm

167:   Input Parameters
168: . DM - the dmnetwork object

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

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

176:   Level: intermediate

178: .seealso: DMNetworkSetSizes, DMNetworkSetEdgeList
179: @*/
180: PetscErrorCode DMNetworkLayoutSetUp(DM dm)
181: {
183:   DM_Network     *network = (DM_Network*)dm->data;
184:   PetscInt       numCorners=2,spacedim=2,dim = 1; /* One dimensional network */
185:   PetscReal      *vertexcoords=NULL;
186:   PetscInt       i,j,ctr,nsubnet,*eowners,np,*edges,*subnetvtx,vStart;
187:   PetscInt       k,netid,vid, *vidxlTog,*edgelist_couple=NULL;
188:   const PetscInt *cone;
189:   MPI_Comm       comm;
190:   PetscMPIInt    size,rank;

193:   PetscObjectGetComm((PetscObject)dm,&comm);
194:   MPI_Comm_rank(comm,&rank);
195:   MPI_Comm_size(comm,&size);

197:   /* Create the local edgelist for the network by concatenating local input edgelists of the subnetworks */
198:   PetscCalloc2(numCorners*network->nVertices,&vertexcoords,2*network->nEdges,&edges);
199:   nsubnet = network->nsubnet - network->ncsubnet;
200:   ctr = 0;
201:   for (i=0; i < nsubnet; i++) {
202:     for (j = 0; j < network->subnet[i].nedge; j++) {
203:       edges[2*ctr]   = network->subnet[i].vStart + network->subnet[i].edgelist[2*j];
204:       edges[2*ctr+1] = network->subnet[i].vStart + network->subnet[i].edgelist[2*j+1];
205:       ctr++;
206:     }
207:   }

209:   /* Append local coupling edgelists of the subnetworks */
210:   i       = nsubnet; /* netid of coupling subnet */
211:   nsubnet = network->nsubnet;
212:   while (i < nsubnet) {
213:     edgelist_couple = network->subnet[i].edgelist;
214:     k = 0;
215:     for (j = 0; j < network->subnet[i].nedge; j++) {
216:       netid = edgelist_couple[k]; vid = edgelist_couple[k+1];
217:       edges[2*ctr] = network->subnet[netid].vStart + vid; k += 2;

219:       netid = edgelist_couple[k]; vid = edgelist_couple[k+1];
220:       edges[2*ctr+1] = network->subnet[netid].vStart + vid; k+=2;
221:       ctr++;
222:     }
223:     i++;
224:   }
225:   /*
226:   if (rank == 0) {
227:     PetscPrintf(PETSC_COMM_SELF,"[%d] edgelist:\n",rank);
228:     for(i=0; i < network->nEdges; i++) {
229:       PetscPrintf(PETSC_COMM_SELF,"[%D %D]",edges[2*i],edges[2*i+1]);
230:       printf("\n");
231:     }
232:   }
233:    */

235:   /* Create network->plex */
236: #if defined(PETSC_USE_64BIT_INDICES)
237:   {
238:     int *edges64;
239:     np = network->nEdges*numCorners;
240:     PetscMalloc1(np,&edges64);
241:     for (i=0; i<np; i++) edges64[i] = (int)edges[i];

243:     if (size == 1) {
244:       DMPlexCreateFromCellList(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges64,spacedim,(const double*)vertexcoords,&network->plex);
245:     } else {
246:       DMPlexCreateFromCellListParallel(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges64,spacedim,(const PetscReal*)vertexcoords,NULL,&network->plex);
247:     }
248:     PetscFree(edges64);
249:   }
250: #else
251:   if (size == 1) {
252:     DMPlexCreateFromCellList(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges,spacedim,(const double*)vertexcoords,&network->plex);
253:   } else {
254:     DMPlexCreateFromCellListParallel(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges,spacedim,(const PetscReal*)vertexcoords,NULL,&network->plex);
255:   }
256: #endif

258:   DMPlexGetChart(network->plex,&network->pStart,&network->pEnd);
259:   DMPlexGetHeightStratum(network->plex,0,&network->eStart,&network->eEnd);
260:   DMPlexGetHeightStratum(network->plex,1,&network->vStart,&network->vEnd);
261:   vStart = network->vStart;

263:   PetscSectionCreate(comm,&network->DataSection);
264:   PetscSectionCreate(comm,&network->DofSection);
265:   PetscSectionSetChart(network->DataSection,network->pStart,network->pEnd);
266:   PetscSectionSetChart(network->DofSection,network->pStart,network->pEnd);

268:   network->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
269:   np = network->pEnd - network->pStart;
270:   PetscCalloc2(np,&network->header,np,&network->cvalue);

272:   /* Create vidxlTog: maps local vertex index to global index */
273:   np = network->vEnd - vStart;
274:   PetscMalloc2(np,&vidxlTog,size+1,&eowners);
275:   ctr = 0;
276:   for (i=network->eStart; i<network->eEnd; i++) {
277:     DMNetworkGetConnectedVertices(dm,i,&cone);
278:     vidxlTog[cone[0] - vStart] = edges[2*ctr];
279:     vidxlTog[cone[1] - vStart] = edges[2*ctr+1];
280:     ctr++;
281:   }
282:   PetscFree2(vertexcoords,edges);

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

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


296:   /* Get edge ownership */
297:   np = network->eEnd - network->eStart;
298:   MPI_Allgather(&np,1,MPIU_INT,eowners+1,1,MPIU_INT,comm);
299:   eowners[0] = 0;
300:   for (i=2; i<=size; i++) eowners[i] += eowners[i-1];

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

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

310:       network->header[i].ndata = 0;
311:       PetscSectionAddDof(network->DataSection,i,network->dataheadersize);
312:       network->header[i].offset[0] = 0;
313:       i++;
314:     }
315:     if (i >= network->subnet[j].eEnd) j++;
316:   }

318:   /* Count network->subnet[*].nvtx */
319:   for (i=vStart; i<network->vEnd; i++) { /* local vertices, including ghosts */
320:     k = vidxlTog[i-vStart];
321:     for (j=0; j < network->nsubnet; j++) {
322:       if (network->subnet[j].vStart <= k && k < network->subnet[j].vEnd) {
323:         network->subnet[j].nvtx++;
324:         break;
325:       }
326:     }
327:   }

329:   /* Set network->subnet[*].vertices on array network->subnetvtx */
330:   subnetvtx = network->subnetvtx;
331:   for (j=0; j<network->nsubnet; j++) {
332:     network->subnet[j].vertices = subnetvtx;
333:     subnetvtx                  += network->subnet[j].nvtx;
334:     network->subnet[j].nvtx = 0;
335:   }

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

341:     k = vidxlTog[i-vStart];
342:     for (j=0; j < network->nsubnet; j++) {
343:       if (network->subnet[j].vStart <= k && k < network->subnet[j].vEnd) {
344:         network->header[i].subnetid = j;
345:         network->subnet[j].vertices[network->subnet[j].nvtx++] = i;
346:         break;
347:       }
348:     }

350:     network->header[i].ndata = 0;
351:     PetscSectionAddDof(network->DataSection,i,network->dataheadersize);
352:     network->header[i].offset[0] = 0;
353:   }

355:   PetscFree2(vidxlTog,eowners);
356:   return(0);
357: }

359: /*@C
360:   DMNetworkGetSubnetworkInfo - Returns the info for the subnetwork

362:   Input Parameters
363: + dm - the DM object
364: - id   - the ID (integer) of the subnetwork

366:   Output Parameters
367: + nv    - number of vertices (local)
368: . ne    - number of edges (local)
369: . vtx   - local vertices for this subnetwork
370: - edge  - local edges for this subnetwork

372:   Notes:
373:   Cannot call this routine before DMNetworkLayoutSetup()

375:   Level: intermediate

377: .seealso: DMNetworkLayoutSetUp, DMNetworkCreate
378: @*/
379: PetscErrorCode DMNetworkGetSubnetworkInfo(DM dm,PetscInt id,PetscInt *nv, PetscInt *ne,const PetscInt **vtx, const PetscInt **edge)
380: {
381:   DM_Network *network = (DM_Network*)dm->data;

384:   *nv   = network->subnet[id].nvtx;
385:   *ne   = network->subnet[id].nedge;
386:   *vtx  = network->subnet[id].vertices;
387:   *edge = network->subnet[id].edges;
388:   return(0);
389: }

391: /*@C
392:   DMNetworkGetSubnetworkCoupleInfo - Returns the info for the coupling subnetwork

394:   Input Parameters
395: + dm - the DM object
396: - id   - the ID (integer) of the coupling subnetwork

398:   Output Parameters
399: + ne - number of edges (local)
400: - edge  - local edges for this coupling subnetwork

402:   Notes:
403:   Cannot call this routine before DMNetworkLayoutSetup()

405:   Level: intermediate

407: .seealso: DMNetworkGetSubnetworkInfo, DMNetworkLayoutSetUp, DMNetworkCreate
408: @*/
409: PetscErrorCode DMNetworkGetSubnetworkCoupleInfo(DM dm,PetscInt id,PetscInt *ne,const PetscInt **edge)
410: {
411:   DM_Network *net = (DM_Network*)dm->data;
412:   PetscInt   id1 = id + net->nsubnet - net->ncsubnet;

415:   *ne   = net->subnet[id1].nedge;
416:   *edge = net->subnet[id1].edges;
417:   return(0);
418: }

420: /*@C
421:   DMNetworkRegisterComponent - Registers the network component

423:   Logically collective on dm

425:   Input Parameters
426: + dm   - the network object
427: . name - the component name
428: - size - the storage size in bytes for this component data

430:    Output Parameters
431: .   key - an integer key that defines the component

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

436:    Level: intermediate

438: .seealso: DMNetworkLayoutSetUp, DMNetworkCreate
439: @*/
440: PetscErrorCode DMNetworkRegisterComponent(DM dm,const char *name,size_t size,PetscInt *key)
441: {
442:   PetscErrorCode        ierr;
443:   DM_Network            *network = (DM_Network*) dm->data;
444:   DMNetworkComponent    *component=&network->component[network->ncomponent];
445:   PetscBool             flg=PETSC_FALSE;
446:   PetscInt              i;

449:   for (i=0; i < network->ncomponent; i++) {
450:     PetscStrcmp(component->name,name,&flg);
451:     if (flg) {
452:       *key = i;
453:       return(0);
454:     }
455:   }
456:   if(network->ncomponent == MAX_COMPONENTS) {
457:     SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Number of components registered exceeds the max %D",MAX_COMPONENTS);
458:   }

460:   PetscStrcpy(component->name,name);
461:   component->size = size/sizeof(DMNetworkComponentGenericDataType);
462:   *key = network->ncomponent;
463:   network->ncomponent++;
464:   return(0);
465: }

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

470:   Not Collective

472:   Input Parameters:
473: . dm - The DMNetwork object

475:   Output Paramters:
476: + vStart - The first vertex point
477: - vEnd   - One beyond the last vertex point

479:   Level: intermediate

481: .seealso: DMNetworkGetEdgeRange
482: @*/
483: PetscErrorCode DMNetworkGetVertexRange(DM dm,PetscInt *vStart,PetscInt *vEnd)
484: {
485:   DM_Network     *network = (DM_Network*)dm->data;

488:   if (vStart) *vStart = network->vStart;
489:   if (vEnd) *vEnd = network->vEnd;
490:   return(0);
491: }

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

496:   Not Collective

498:   Input Parameters:
499: . dm - The DMNetwork object

501:   Output Paramters:
502: + eStart - The first edge point
503: - eEnd   - One beyond the last edge point

505:   Level: intermediate

507: .seealso: DMNetworkGetVertexRange
508: @*/
509: PetscErrorCode DMNetworkGetEdgeRange(DM dm,PetscInt *eStart,PetscInt *eEnd)
510: {
511:   DM_Network     *network = (DM_Network*)dm->data;

514:   if (eStart) *eStart = network->eStart;
515:   if (eEnd) *eEnd = network->eEnd;
516:   return(0);
517: }

519: /*@
520:   DMNetworkGetGlobalEdgeIndex - Get the user global numbering for the edge.

522:   Not Collective

524:   Input Parameters:
525: + dm - DMNetwork object
526: - p  - edge point

528:   Output Paramters:
529: . index - user global numbering for the edge

531:   Level: intermediate

533: .seealso: DMNetworkGetGlobalVertexIndex
534: @*/
535: PetscErrorCode DMNetworkGetGlobalEdgeIndex(DM dm,PetscInt p,PetscInt *index)
536: {
537:   PetscErrorCode    ierr;
538:   DM_Network        *network = (DM_Network*)dm->data;
539:   PetscInt          offsetp;
540:   DMNetworkComponentHeader header;

543:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
544:   PetscSectionGetOffset(network->DataSection,p,&offsetp);
545:   header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
546:   *index = header->index;
547:   return(0);
548: }

550: /*@
551:   DMNetworkGetGlobalVertexIndex - Get the user global numbering for the vertex.

553:   Not Collective

555:   Input Parameters:
556: + dm - DMNetwork object
557: - p  - vertex point

559:   Output Paramters:
560: . index - user global numbering for the vertex

562:   Level: intermediate

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

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

581: /*
582:   DMNetworkGetComponentKeyOffset - Gets the type along with the offset for indexing the
583:                                     component value from the component data array

585:   Not Collective

587:   Input Parameters:
588: + dm      - The DMNetwork object
589: . p       - vertex/edge point
590: - compnum - component number

592:   Output Parameters:
593: + compkey - the key obtained when registering the component
594: - offset  - offset into the component data array associated with the vertex/edge point

596:   Notes:
597:   Typical usage:

599:   DMNetworkGetComponentDataArray(dm, &arr);
600:   DMNetworkGetVertex/EdgeRange(dm,&Start,&End);
601:   Loop over vertices or edges
602:     DMNetworkGetNumComponents(dm,v,&numcomps);
603:     Loop over numcomps
604:       DMNetworkGetComponentKeyOffset(dm,v,compnum,&key,&offset);
605:       compdata = (UserCompDataType)(arr+offset);

607:   Level: intermediate

609: .seealso: DMNetworkGetNumComponents, DMNetworkGetComponentDataArray,
610: */
611: PetscErrorCode DMNetworkGetComponentKeyOffset(DM dm,PetscInt p, PetscInt compnum, PetscInt *compkey, PetscInt *offset)
612: {
613:   PetscErrorCode           ierr;
614:   PetscInt                 offsetp;
615:   DMNetworkComponentHeader header;
616:   DM_Network               *network = (DM_Network*)dm->data;

619:   PetscSectionGetOffset(network->DataSection,p,&offsetp);
620:   header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
621:   if (compkey) *compkey = header->key[compnum];
622:   if (offset) *offset  = offsetp+network->dataheadersize+header->offset[compnum];
623:   return(0);
624: }

626: /*@
627:   DMNetworkGetComponent - Returns the network component and its key

629:   Not Collective

631:   Input Parameters
632: + dm - DMNetwork object
633: . p  - edge or vertex point
634: - compnum - component number

636:   Output Parameters:
637: + compkey - the key set for this computing during registration
638: - component - the component data

640:   Notes:
641:   Typical usage:

643:   DMNetworkGetVertex/EdgeRange(dm,&Start,&End);
644:   Loop over vertices or edges
645:     DMNetworkGetNumComponents(dm,v,&numcomps);
646:     Loop over numcomps
647:       DMNetworkGetComponent(dm,v,compnum,&key,&component);

649:   Level: intermediate

651: .seealso: DMNetworkGetNumComponents, DMNetworkGetVariableOffset
652: @*/
653: PetscErrorCode DMNetworkGetComponent(DM dm, PetscInt p, PetscInt compnum, PetscInt *key, void **component)
654: {
656:   DM_Network     *network = (DM_Network*)dm->data;
657:   PetscInt       offsetd = 0;

660:   DMNetworkGetComponentKeyOffset(dm,p,compnum,key,&offsetd);
661:   *component = network->componentdataarray+offsetd;
662:   return(0);
663: }

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

668:   Not Collective

670:   Input Parameters:
671: + dm           - The DMNetwork object
672: . p            - vertex/edge point
673: . componentkey - component key returned while registering the component
674: - compvalue    - pointer to the data structure for the component

676:   Level: intermediate

678: .seealso: DMNetworkGetVertexRange, DMNetworkGetEdgeRange, DMNetworkRegisterComponent
679: @*/
680: PetscErrorCode DMNetworkAddComponent(DM dm, PetscInt p,PetscInt componentkey,void* compvalue)
681: {
682:   DM_Network               *network = (DM_Network*)dm->data;
683:   DMNetworkComponent       *component = &network->component[componentkey];
684:   DMNetworkComponentHeader header = &network->header[p];
685:   DMNetworkComponentValue  cvalue = &network->cvalue[p];
686:   PetscErrorCode           ierr;

689:   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);

691:   header->size[header->ndata] = component->size;
692:   PetscSectionAddDof(network->DataSection,p,component->size);
693:   header->key[header->ndata] = componentkey;
694:   if (header->ndata != 0) header->offset[header->ndata] = header->offset[header->ndata-1] + header->size[header->ndata-1];

696:   cvalue->data[header->ndata] = (void*)compvalue;
697:   header->ndata++;
698:   return(0);
699: }

701: /*@
702:   DMNetworkGetNumComponents - Get the number of components at a vertex/edge

704:   Not Collective

706:   Input Parameters:
707: + dm - The DMNetwork object
708: - p  - vertex/edge point

710:   Output Parameters:
711: . numcomponents - Number of components at the vertex/edge

713:   Level: intermediate

715: .seealso: DMNetworkRegisterComponent, DMNetworkAddComponent
716: @*/
717: PetscErrorCode DMNetworkGetNumComponents(DM dm,PetscInt p,PetscInt *numcomponents)
718: {
720:   PetscInt       offset;
721:   DM_Network     *network = (DM_Network*)dm->data;

724:   PetscSectionGetOffset(network->DataSection,p,&offset);
725:   *numcomponents = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata;
726:   return(0);
727: }

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

732:   Not Collective

734:   Input Parameters:
735: + dm     - The DMNetwork object
736: - p      - the edge/vertex point

738:   Output Parameters:
739: . offset - the offset

741:   Level: intermediate

743: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
744: @*/
745: PetscErrorCode DMNetworkGetVariableOffset(DM dm,PetscInt p,PetscInt *offset)
746: {
748:   DM_Network     *network = (DM_Network*)dm->data;

751:   PetscSectionGetOffset(network->plex->localSection,p,offset);
752:   return(0);
753: }

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

758:   Not Collective

760:   Input Parameters:
761: + dm      - The DMNetwork object
762: - p       - the edge/vertex point

764:   Output Parameters:
765: . offsetg - the offset

767:   Level: intermediate

769: .seealso: DMNetworkGetVariableOffset, DMGetLocalVector
770: @*/
771: PetscErrorCode DMNetworkGetVariableGlobalOffset(DM dm,PetscInt p,PetscInt *offsetg)
772: {
774:   DM_Network     *network = (DM_Network*)dm->data;

777:   PetscSectionGetOffset(network->plex->globalSection,p,offsetg);
778:   if (*offsetg < 0) *offsetg = -(*offsetg + 1); /* Convert to actual global offset for ghost vertex */
779:   return(0);
780: }

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

785:   Not Collective

787:   Input Parameters:
788: + dm     - The DMNetwork object
789: - p      - the edge point

791:   Output Parameters:
792: . offset - the offset

794:   Level: intermediate

796: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
797: @*/
798: PetscErrorCode DMNetworkGetEdgeOffset(DM dm,PetscInt p,PetscInt *offset)
799: {
801:   DM_Network     *network = (DM_Network*)dm->data;


805:   PetscSectionGetOffset(network->edge.DofSection,p,offset);
806:   return(0);
807: }

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

812:   Not Collective

814:   Input Parameters:
815: + dm     - The DMNetwork object
816: - p      - the vertex point

818:   Output Parameters:
819: . offset - the offset

821:   Level: intermediate

823: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
824: @*/
825: PetscErrorCode DMNetworkGetVertexOffset(DM dm,PetscInt p,PetscInt *offset)
826: {
828:   DM_Network     *network = (DM_Network*)dm->data;


832:   p -= network->vStart;

834:   PetscSectionGetOffset(network->vertex.DofSection,p,offset);
835:   return(0);
836: }
837: /*@
838:   DMNetworkAddNumVariables - Add number of variables associated with a given point.

840:   Not Collective

842:   Input Parameters:
843: + dm   - The DMNetworkObject
844: . p    - the vertex/edge point
845: - nvar - number of additional variables

847:   Level: intermediate

849: .seealso: DMNetworkSetNumVariables
850: @*/
851: PetscErrorCode DMNetworkAddNumVariables(DM dm,PetscInt p,PetscInt nvar)
852: {
854:   DM_Network     *network = (DM_Network*)dm->data;

857:   PetscSectionAddDof(network->DofSection,p,nvar);
858:   return(0);
859: }

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

864:   Not Collective

866:   Input Parameters:
867: + dm   - The DMNetworkObject
868: - p    - the vertex/edge point

870:   Output Parameters:
871: . nvar - number of variables

873:   Level: intermediate

875: .seealso: DMNetworkAddNumVariables, DMNetworkSddNumVariables
876: @*/
877: PetscErrorCode DMNetworkGetNumVariables(DM dm,PetscInt p,PetscInt *nvar)
878: {
880:   DM_Network     *network = (DM_Network*)dm->data;

883:   PetscSectionGetDof(network->DofSection,p,nvar);
884:   return(0);
885: }

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

890:   Not Collective

892:   Input Parameters:
893: + dm   - The DMNetworkObject
894: . p    - the vertex/edge point
895: - nvar - number of variables

897:   Level: intermediate

899: .seealso: DMNetworkAddNumVariables
900: @*/
901: PetscErrorCode DMNetworkSetNumVariables(DM dm,PetscInt p,PetscInt nvar)
902: {
904:   DM_Network     *network = (DM_Network*)dm->data;

907:   PetscSectionSetDof(network->DofSection,p,nvar);
908:   return(0);
909: }

911: /* Sets up the array that holds the data for all components and its associated section. This
912:    function is called during DMSetUp() */
913: PetscErrorCode DMNetworkComponentSetUp(DM dm)
914: {
915:   PetscErrorCode           ierr;
916:   DM_Network               *network = (DM_Network*)dm->data;
917:   PetscInt                 arr_size,p,offset,offsetp,ncomp,i;
918:   DMNetworkComponentHeader header;
919:   DMNetworkComponentValue  cvalue;
920:   DMNetworkComponentGenericDataType *componentdataarray;

923:   PetscSectionSetUp(network->DataSection);
924:   PetscSectionGetStorageSize(network->DataSection,&arr_size);
925:   PetscMalloc1(arr_size,&network->componentdataarray);
926:   componentdataarray = network->componentdataarray;
927:   for (p = network->pStart; p < network->pEnd; p++) {
928:     PetscSectionGetOffset(network->DataSection,p,&offsetp);
929:     /* Copy header */
930:     header = &network->header[p];
931:     PetscMemcpy(componentdataarray+offsetp,header,network->dataheadersize*sizeof(DMNetworkComponentGenericDataType));
932:     /* Copy data */
933:     cvalue = &network->cvalue[p];
934:     ncomp = header->ndata;
935:     for (i = 0; i < ncomp; i++) {
936:       offset = offsetp + network->dataheadersize + header->offset[i];
937:       PetscMemcpy(componentdataarray+offset,cvalue->data[i],header->size[i]*sizeof(DMNetworkComponentGenericDataType));
938:     }
939:   }
940:   return(0);
941: }

943: /* Sets up the section for dofs. This routine is called during DMSetUp() */
944: PetscErrorCode DMNetworkVariablesSetUp(DM dm)
945: {
947:   DM_Network     *network = (DM_Network*)dm->data;

950:   PetscSectionSetUp(network->DofSection);
951:   return(0);
952: }

954: /*@C
955:   DMNetworkGetComponentDataArray - Returns the component data array

957:   Not Collective

959:   Input Parameters:
960: . dm - The DMNetwork Object

962:   Output Parameters:
963: . componentdataarray - array that holds data for all components

965:   Level: intermediate

967: .seealso: DMNetworkGetComponentKeyOffset, DMNetworkGetNumComponents
968: @*/
969: PetscErrorCode DMNetworkGetComponentDataArray(DM dm,DMNetworkComponentGenericDataType **componentdataarray)
970: {
971:   DM_Network     *network = (DM_Network*)dm->data;

974:   *componentdataarray = network->componentdataarray;
975:   return(0);
976: }

978: /* Get a subsection from a range of points */
979: PetscErrorCode DMNetworkGetSubSection_private(PetscSection master, PetscInt pstart, PetscInt pend,PetscSection *subsection)
980: {
982:   PetscInt       i, nvar;

985:   PetscSectionCreate(PetscObjectComm((PetscObject)master), subsection);
986:   PetscSectionSetChart(*subsection, 0, pend - pstart);
987:   for (i = pstart; i < pend; i++) {
988:     PetscSectionGetDof(master,i,&nvar);
989:     PetscSectionSetDof(*subsection, i - pstart, nvar);
990:   }

992:   PetscSectionSetUp(*subsection);
993:   return(0);
994: }

996: /* Create a submap of points with a GlobalToLocal structure */
997: PetscErrorCode DMNetworkSetSubMap_private(PetscInt pstart, PetscInt pend, ISLocalToGlobalMapping *map)
998: {
1000:   PetscInt       i, *subpoints;

1003:   /* Create index sets to map from "points" to "subpoints" */
1004:   PetscMalloc1(pend - pstart, &subpoints);
1005:   for (i = pstart; i < pend; i++) {
1006:     subpoints[i - pstart] = i;
1007:   }
1008:   ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,pend-pstart,subpoints,PETSC_COPY_VALUES,map);
1009:   PetscFree(subpoints);
1010:   return(0);
1011: }

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

1016:   Collective

1018:   Input Parameters:
1019: . dm   - The DMNetworkObject

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

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

1025:   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]).

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

1029:   Level: intermediate

1031: @*/
1032: PetscErrorCode DMNetworkAssembleGraphStructures(DM dm)
1033: {
1035:   MPI_Comm       comm;
1036:   PetscMPIInt    rank, size;
1037:   DM_Network     *network = (DM_Network*)dm->data;

1040:   PetscObjectGetComm((PetscObject)dm,&comm);
1041:   MPI_Comm_rank(comm, &rank);
1042:   MPI_Comm_size(comm, &size);

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

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

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

1055:     PetscSectionCreateGlobalSection(network->vertex.DofSection, network->vertex.sf, PETSC_FALSE, PETSC_FALSE, &network->vertex.GlobalDofSection);
1056:     PetscSFGetSubSF(network->plex->sf, network->edge.mapping, &network->edge.sf);
1057:     PetscSectionCreateGlobalSection(network->edge.DofSection, network->edge.sf, PETSC_FALSE, PETSC_FALSE, &network->edge.GlobalDofSection);
1058:   } else {
1059:     /* create structures for vertex */
1060:     PetscSectionClone(network->vertex.DofSection,&network->vertex.GlobalDofSection);
1061:     /* create structures for edge */
1062:     PetscSectionClone(network->edge.DofSection,&network->edge.GlobalDofSection);
1063:   }

1065:   /* Add viewers */
1066:   PetscObjectSetName((PetscObject)network->edge.GlobalDofSection,"Global edge dof section");
1067:   PetscObjectSetName((PetscObject)network->vertex.GlobalDofSection,"Global vertex dof section");
1068:   PetscSectionViewFromOptions(network->edge.GlobalDofSection, NULL, "-edge_global_section_view");
1069:   PetscSectionViewFromOptions(network->vertex.GlobalDofSection, NULL, "-vertex_global_section_view");
1070:   return(0);
1071: }

1073: /*@
1074:   DMNetworkDistribute - Distributes the network and moves associated component data.

1076:   Collective

1078:   Input Parameter:
1079: + DM - the DMNetwork object
1080: - overlap - The overlap of partitions, 0 is the default

1082:   Notes:
1083:   Distributes the network with <overlap>-overlapping partitioning of the edges.

1085:   Level: intermediate

1087: .seealso: DMNetworkCreate
1088: @*/
1089: PetscErrorCode DMNetworkDistribute(DM *dm,PetscInt overlap)
1090: {
1091:   MPI_Comm       comm;
1093:   PetscMPIInt    size;
1094:   DM_Network     *oldDMnetwork = (DM_Network*)((*dm)->data);
1095:   DM_Network     *newDMnetwork;
1096:   PetscSF        pointsf=NULL;
1097:   DM             newDM;
1098:   PetscInt       j,e,v,offset,*subnetvtx;
1099:   PetscPartitioner         part;
1100:   DMNetworkComponentHeader header;

1103:   PetscObjectGetComm((PetscObject)*dm,&comm);
1104:   MPI_Comm_size(comm, &size);
1105:   if (size == 1) return(0);

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

1111:   /* Enable runtime options for petscpartitioner */
1112:   DMPlexGetPartitioner(oldDMnetwork->plex,&part);
1113:   PetscPartitionerSetFromOptions(part);

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

1118:   /* Distribute dof section */
1119:   PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DofSection);
1120:   PetscSFDistributeSection(pointsf,oldDMnetwork->DofSection,NULL,newDMnetwork->DofSection);
1121:   PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DataSection);

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

1126:   PetscSectionGetChart(newDMnetwork->DataSection,&newDMnetwork->pStart,&newDMnetwork->pEnd);
1127:   DMPlexGetHeightStratum(newDMnetwork->plex,0, &newDMnetwork->eStart,&newDMnetwork->eEnd);
1128:   DMPlexGetHeightStratum(newDMnetwork->plex,1,&newDMnetwork->vStart,&newDMnetwork->vEnd);
1129:   newDMnetwork->nEdges    = newDMnetwork->eEnd - newDMnetwork->eStart;
1130:   newDMnetwork->nVertices = newDMnetwork->vEnd - newDMnetwork->vStart;
1131:   newDMnetwork->NVertices = oldDMnetwork->NVertices;
1132:   newDMnetwork->NEdges    = oldDMnetwork->NEdges;

1134:   /* Set Dof section as the section for dm */
1135:   DMSetLocalSection(newDMnetwork->plex,newDMnetwork->DofSection);
1136:   DMGetGlobalSection(newDMnetwork->plex,&newDMnetwork->GlobalDofSection);

1138:   /* Set up subnetwork info in the newDM */
1139:   newDMnetwork->nsubnet  = oldDMnetwork->nsubnet;
1140:   newDMnetwork->ncsubnet = oldDMnetwork->ncsubnet;
1141:   PetscCalloc1(newDMnetwork->nsubnet,&newDMnetwork->subnet);
1142:   /* Copy over the global number of vertices and edges in each subnetwork. Note that these are already
1143:      calculated in DMNetworkLayoutSetUp()
1144:   */
1145:   for(j=0; j < newDMnetwork->nsubnet; j++) {
1146:     newDMnetwork->subnet[j].Nvtx  = oldDMnetwork->subnet[j].Nvtx;
1147:     newDMnetwork->subnet[j].Nedge = oldDMnetwork->subnet[j].Nedge;
1148:   }

1150:   for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++ ) {
1151:     PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);
1152:     header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1153:     newDMnetwork->subnet[header->subnetid].nedge++;
1154:   }

1156:   for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++ ) {
1157:     PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);
1158:     header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1159:     newDMnetwork->subnet[header->subnetid].nvtx++;
1160:   }

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

1166:   for (j=0; j<newDMnetwork->nsubnet; j++) {
1167:     PetscCalloc1(newDMnetwork->subnet[j].nedge,&newDMnetwork->subnet[j].edges);
1168:     newDMnetwork->subnet[j].vertices = subnetvtx;
1169:     subnetvtx                       += newDMnetwork->subnet[j].nvtx;

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

1176:   /* Set the vertices and edges in each subnetwork */
1177:   for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++ ) {
1178:     PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);
1179:     header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1180:     newDMnetwork->subnet[header->subnetid].edges[newDMnetwork->subnet[header->subnetid].nedge++] = e;
1181:   }

1183:   for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++ ) {
1184:     PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);
1185:     header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1186:     newDMnetwork->subnet[header->subnetid].vertices[newDMnetwork->subnet[header->subnetid].nvtx++] = v;
1187:   }

1189:   newDM->setupcalled = (*dm)->setupcalled;
1190:   newDMnetwork->distributecalled = PETSC_TRUE;

1192:   /* Destroy point SF */
1193:   PetscSFDestroy(&pointsf);

1195:   DMDestroy(dm);
1196:   *dm  = newDM;
1197:   return(0);
1198: }

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

1203:   Input Parameters:
1204: + masterSF - the original SF structure
1205: - map      - a ISLocalToGlobal mapping that contains the subset of points

1207:   Output Parameters:
1208: . subSF    - a subset of the masterSF for the desired subset.
1209: */

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

1213:   PetscErrorCode        ierr;
1214:   PetscInt              nroots, nleaves, *ilocal_sub;
1215:   PetscInt              i, *ilocal_map, nroots_sub, nleaves_sub = 0;
1216:   PetscInt              *local_points, *remote_points;
1217:   PetscSFNode           *iremote_sub;
1218:   const PetscInt        *ilocal;
1219:   const PetscSFNode     *iremote;

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

1224:   /* Look for leaves that pertain to the subset of points. Get the local ordering */
1225:   PetscMalloc1(nleaves,&ilocal_map);
1226:   ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nleaves,ilocal,NULL,ilocal_map);
1227:   for (i = 0; i < nleaves; i++) {
1228:     if (ilocal_map[i] != -1) nleaves_sub += 1;
1229:   }
1230:   /* Re-number ilocal with subset numbering. Need information from roots */
1231:   PetscMalloc2(nroots,&local_points,nroots,&remote_points);
1232:   for (i = 0; i < nroots; i++) local_points[i] = i;
1233:   ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nroots,local_points,NULL,local_points);
1234:   PetscSFBcastBegin(mastersf, MPIU_INT, local_points, remote_points);
1235:   PetscSFBcastEnd(mastersf, MPIU_INT, local_points, remote_points);
1236:   /* Fill up graph using local (that is, local to the subset) numbering. */
1237:   PetscMalloc1(nleaves_sub,&ilocal_sub);
1238:   PetscMalloc1(nleaves_sub,&iremote_sub);
1239:   nleaves_sub = 0;
1240:   for (i = 0; i < nleaves; i++) {
1241:     if (ilocal_map[i] != -1) {
1242:       ilocal_sub[nleaves_sub] = ilocal_map[i];
1243:       iremote_sub[nleaves_sub].rank = iremote[i].rank;
1244:       iremote_sub[nleaves_sub].index = remote_points[ilocal[i]];
1245:       nleaves_sub += 1;
1246:     }
1247:   }
1248:   PetscFree2(local_points,remote_points);
1249:   ISLocalToGlobalMappingGetSize(map,&nroots_sub);

1251:   /* Create new subSF */
1252:   PetscSFCreate(PETSC_COMM_WORLD,subSF);
1253:   PetscSFSetFromOptions(*subSF);
1254:   PetscSFSetGraph(*subSF,nroots_sub,nleaves_sub,ilocal_sub,PETSC_OWN_POINTER,iremote_sub,PETSC_COPY_VALUES);
1255:   PetscFree(ilocal_map);
1256:   PetscFree(iremote_sub);
1257:   return(0);
1258: }

1260: /*@C
1261:   DMNetworkGetSupportingEdges - Return the supporting edges for this vertex point

1263:   Not Collective

1265:   Input Parameters:
1266: + dm - The DMNetwork object
1267: - p  - the vertex point

1269:   Output Paramters:
1270: + nedges - number of edges connected to this vertex point
1271: - edges  - List of edge points

1273:   Level: intermediate

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

1279: .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices
1280: @*/
1281: PetscErrorCode DMNetworkGetSupportingEdges(DM dm,PetscInt vertex,PetscInt *nedges,const PetscInt *edges[])
1282: {
1284:   DM_Network     *network = (DM_Network*)dm->data;

1287:   DMPlexGetSupportSize(network->plex,vertex,nedges);
1288:   DMPlexGetSupport(network->plex,vertex,edges);
1289:   return(0);
1290: }

1292: /*@C
1293:   DMNetworkGetConnectedVertices - Return the connected vertices for this edge point

1295:   Not Collective

1297:   Input Parameters:
1298: + dm - The DMNetwork object
1299: - p  - the edge point

1301:   Output Paramters:
1302: . vertices  - vertices connected to this edge

1304:   Level: intermediate

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

1310: .seealso: DMNetworkCreate, DMNetworkGetSupportingEdges
1311: @*/
1312: PetscErrorCode DMNetworkGetConnectedVertices(DM dm,PetscInt edge,const PetscInt *vertices[])
1313: {
1315:   DM_Network     *network = (DM_Network*)dm->data;

1318:   DMPlexGetCone(network->plex,edge,vertices);
1319:   return(0);
1320: }

1322: /*@
1323:   DMNetworkIsGhostVertex - Returns TRUE if the vertex is a ghost vertex

1325:   Not Collective

1327:   Input Parameters:
1328: + dm - The DMNetwork object
1329: - p  - the vertex point

1331:   Output Parameter:
1332: . isghost - TRUE if the vertex is a ghost point

1334:   Level: intermediate

1336: .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices, DMNetworkGetVertexRange
1337: @*/
1338: PetscErrorCode DMNetworkIsGhostVertex(DM dm,PetscInt p,PetscBool *isghost)
1339: {
1341:   DM_Network     *network = (DM_Network*)dm->data;
1342:   PetscInt       offsetg;
1343:   PetscSection   sectiong;

1346:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
1347:   *isghost = PETSC_FALSE;
1348:   DMGetGlobalSection(network->plex,&sectiong);
1349:   PetscSectionGetOffset(sectiong,p,&offsetg);
1350:   if (offsetg < 0) *isghost = PETSC_TRUE;
1351:   return(0);
1352: }

1354: PetscErrorCode DMSetUp_Network(DM dm)
1355: {
1357:   DM_Network     *network=(DM_Network*)dm->data;

1360:   DMNetworkComponentSetUp(dm);
1361:   DMNetworkVariablesSetUp(dm);

1363:   DMSetLocalSection(network->plex,network->DofSection);
1364:   DMGetGlobalSection(network->plex,&network->GlobalDofSection);

1366:   dm->setupcalled = PETSC_TRUE;
1367:   DMViewFromOptions(dm,NULL,"-dm_view");
1368:   return(0);
1369: }

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

1375:     Collective

1377:     Input Parameters:
1378: +   dm - The DMNetwork object
1379: .   eflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for edges
1380: -   vflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for vertices

1382:     Level: intermediate

1384: @*/
1385: PetscErrorCode DMNetworkHasJacobian(DM dm,PetscBool eflg,PetscBool vflg)
1386: {
1387:   DM_Network     *network=(DM_Network*)dm->data;
1389:   PetscInt       nVertices = network->nVertices;

1392:   network->userEdgeJacobian   = eflg;
1393:   network->userVertexJacobian = vflg;

1395:   if (eflg && !network->Je) {
1396:     PetscCalloc1(3*network->nEdges,&network->Je);
1397:   }

1399:   if (vflg && !network->Jv && nVertices) {
1400:     PetscInt       i,*vptr,nedges,vStart=network->vStart;
1401:     PetscInt       nedges_total;
1402:     const PetscInt *edges;

1404:     /* count nvertex_total */
1405:     nedges_total = 0;
1406:     PetscMalloc1(nVertices+1,&vptr);

1408:     vptr[0] = 0;
1409:     for (i=0; i<nVertices; i++) {
1410:       DMNetworkGetSupportingEdges(dm,i+vStart,&nedges,&edges);
1411:       nedges_total += nedges;
1412:       vptr[i+1] = vptr[i] + 2*nedges + 1;
1413:     }

1415:     PetscCalloc1(2*nedges_total+nVertices,&network->Jv);
1416:     network->Jvptr = vptr;
1417:   }
1418:   return(0);
1419: }

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

1424:     Not Collective

1426:     Input Parameters:
1427: +   dm - The DMNetwork object
1428: .   p  - the edge point
1429: -   J - array (size = 3) of Jacobian submatrices for this edge point:
1430:         J[0]: this edge
1431:         J[1] and J[2]: connected vertices, obtained by calling DMNetworkGetConnectedVertices()

1433:     Level: intermediate

1435: .seealso: DMNetworkVertexSetMatrix
1436: @*/
1437: PetscErrorCode DMNetworkEdgeSetMatrix(DM dm,PetscInt p,Mat J[])
1438: {
1439:   DM_Network     *network=(DM_Network*)dm->data;

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

1444:   if (J) {
1445:     network->Je[3*p]   = J[0];
1446:     network->Je[3*p+1] = J[1];
1447:     network->Je[3*p+2] = J[2];
1448:   }
1449:   return(0);
1450: }

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

1455:     Not Collective

1457:     Input Parameters:
1458: +   dm - The DMNetwork object
1459: .   p  - the vertex point
1460: -   J - array of Jacobian (size = 2*(num of supporting edges) + 1) submatrices for this vertex point:
1461:         J[0]:       this vertex
1462:         J[1+2*i]:   i-th supporting edge
1463:         J[1+2*i+1]: i-th connected vertex

1465:     Level: intermediate

1467: .seealso: DMNetworkEdgeSetMatrix
1468: @*/
1469: PetscErrorCode DMNetworkVertexSetMatrix(DM dm,PetscInt p,Mat J[])
1470: {
1472:   DM_Network     *network=(DM_Network*)dm->data;
1473:   PetscInt       i,*vptr,nedges,vStart=network->vStart;
1474:   const PetscInt *edges;

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

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

1483:     /* Set Jacobian for each supporting edge and connected vertex */
1484:     DMNetworkGetSupportingEdges(dm,p,&nedges,&edges);
1485:     for (i=1; i<=2*nedges; i++) network->Jv[vptr[p-vStart]+i] = J[i];
1486:   }
1487:   return(0);
1488: }

1490: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1491: {
1493:   PetscInt       j;
1494:   PetscScalar    val=(PetscScalar)ncols;

1497:   if (!ghost) {
1498:     for (j=0; j<nrows; j++) {
1499:       VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
1500:     }
1501:   } else {
1502:     for (j=0; j<nrows; j++) {
1503:       VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
1504:     }
1505:   }
1506:   return(0);
1507: }

1509: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1510: {
1512:   PetscInt       j,ncols_u;
1513:   PetscScalar    val;

1516:   if (!ghost) {
1517:     for (j=0; j<nrows; j++) {
1518:       MatGetRow(Ju,j,&ncols_u,NULL,NULL);
1519:       val = (PetscScalar)ncols_u;
1520:       VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
1521:       MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
1522:     }
1523:   } else {
1524:     for (j=0; j<nrows; j++) {
1525:       MatGetRow(Ju,j,&ncols_u,NULL,NULL);
1526:       val = (PetscScalar)ncols_u;
1527:       VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
1528:       MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
1529:     }
1530:   }
1531:   return(0);
1532: }

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

1539:   if (Ju) {
1540:     MatSetPreallocationUserblock_private(Ju,nrows,rows,ncols,ghost,vdnz,vonz);
1541:   } else {
1542:     MatSetPreallocationDenseblock_private(nrows,rows,ncols,ghost,vdnz,vonz);
1543:   }
1544:   return(0);
1545: }

1547: PETSC_STATIC_INLINE PetscErrorCode MatSetDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1548: {
1550:   PetscInt       j,*cols;
1551:   PetscScalar    *zeros;

1554:   PetscCalloc2(ncols,&cols,nrows*ncols,&zeros);
1555:   for (j=0; j<ncols; j++) cols[j] = j+ cstart;
1556:   MatSetValues(*J,nrows,rows,ncols,cols,zeros,INSERT_VALUES);
1557:   PetscFree2(cols,zeros);
1558:   return(0);
1559: }

1561: PETSC_STATIC_INLINE PetscErrorCode MatSetUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1562: {
1564:   PetscInt       j,M,N,row,col,ncols_u;
1565:   const PetscInt *cols;
1566:   PetscScalar    zero=0.0;

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

1572:   for (row=0; row<nrows; row++) {
1573:     MatGetRow(Ju,row,&ncols_u,&cols,NULL);
1574:     for (j=0; j<ncols_u; j++) {
1575:       col = cols[j] + cstart;
1576:       MatSetValues(*J,1,&rows[row],1,&col,&zero,INSERT_VALUES);
1577:     }
1578:     MatRestoreRow(Ju,row,&ncols_u,&cols,NULL);
1579:   }
1580:   return(0);
1581: }

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

1588:   if (Ju) {
1589:     MatSetUserblock_private(Ju,nrows,rows,ncols,cstart,J);
1590:   } else {
1591:     MatSetDenseblock_private(nrows,rows,ncols,cstart,J);
1592:   }
1593:   return(0);
1594: }

1596: /* 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.
1597: */
1598: PetscErrorCode CreateSubGlobalToLocalMapping_private(PetscSection globalsec, PetscSection localsec, ISLocalToGlobalMapping *ltog)
1599: {
1601:   PetscInt       i,size,dof;
1602:   PetscInt       *glob2loc;

1605:   PetscSectionGetStorageSize(localsec,&size);
1606:   PetscMalloc1(size,&glob2loc);

1608:   for (i = 0; i < size; i++) {
1609:     PetscSectionGetOffset(globalsec,i,&dof);
1610:     dof = (dof >= 0) ? dof : -(dof + 1);
1611:     glob2loc[i] = dof;
1612:   }

1614:   ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,size,glob2loc,PETSC_OWN_POINTER,ltog);
1615: #if 0
1616:   PetscIntView(size,glob2loc,PETSC_VIEWER_STDOUT_WORLD);
1617: #endif
1618:   return(0);
1619: }

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

1623: PetscErrorCode DMCreateMatrix_Network_Nest(DM dm,Mat *J)
1624: {
1626:   DM_Network     *network = (DM_Network*)dm->data;
1627:   PetscMPIInt    rank, size;
1628:   PetscInt       eDof,vDof;
1629:   Mat            j11,j12,j21,j22,bA[2][2];
1630:   MPI_Comm       comm;
1631:   ISLocalToGlobalMapping eISMap,vISMap;

1634:   PetscObjectGetComm((PetscObject)dm,&comm);
1635:   MPI_Comm_rank(comm,&rank);
1636:   MPI_Comm_size(comm,&size);

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

1641:   MatCreate(comm, &j11);
1642:   MatSetSizes(j11, eDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
1643:   MatSetType(j11, MATMPIAIJ);

1645:   MatCreate(comm, &j12);
1646:   MatSetSizes(j12, eDof, vDof, PETSC_DETERMINE ,PETSC_DETERMINE);
1647:   MatSetType(j12, MATMPIAIJ);

1649:   MatCreate(comm, &j21);
1650:   MatSetSizes(j21, vDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
1651:   MatSetType(j21, MATMPIAIJ);

1653:   MatCreate(comm, &j22);
1654:   MatSetSizes(j22, vDof, vDof, PETSC_DETERMINE, PETSC_DETERMINE);
1655:   MatSetType(j22, MATMPIAIJ);

1657:   bA[0][0] = j11;
1658:   bA[0][1] = j12;
1659:   bA[1][0] = j21;
1660:   bA[1][1] = j22;

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

1665:   MatSetLocalToGlobalMapping(j11,eISMap,eISMap);
1666:   MatSetLocalToGlobalMapping(j12,eISMap,vISMap);
1667:   MatSetLocalToGlobalMapping(j21,vISMap,eISMap);
1668:   MatSetLocalToGlobalMapping(j22,vISMap,vISMap);

1670:   MatSetUp(j11);
1671:   MatSetUp(j12);
1672:   MatSetUp(j21);
1673:   MatSetUp(j22);

1675:   MatCreateNest(comm,2,NULL,2,NULL,&bA[0][0],J);
1676:   MatSetUp(*J);
1677:   MatNestSetVecType(*J,VECNEST);
1678:   MatDestroy(&j11);
1679:   MatDestroy(&j12);
1680:   MatDestroy(&j21);
1681:   MatDestroy(&j22);

1683:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1684:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1685:   MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);

1687:   /* Free structures */
1688:   ISLocalToGlobalMappingDestroy(&eISMap);
1689:   ISLocalToGlobalMappingDestroy(&vISMap);
1690:   return(0);
1691: }

1693: PetscErrorCode DMCreateMatrix_Network(DM dm,Mat *J)
1694: {
1696:   DM_Network     *network = (DM_Network*)dm->data;
1697:   PetscInt       eStart,eEnd,vStart,vEnd,rstart,nrows,*rows,localSize;
1698:   PetscInt       cstart,ncols,j,e,v;
1699:   PetscBool      ghost,ghost_vc,ghost2,isNest;
1700:   Mat            Juser;
1701:   PetscSection   sectionGlobal;
1702:   PetscInt       nedges,*vptr=NULL,vc,*rows_v; /* suppress maybe-uninitialized warning */
1703:   const PetscInt *edges,*cone;
1704:   MPI_Comm       comm;
1705:   MatType        mtype;
1706:   Vec            vd_nz,vo_nz;
1707:   PetscInt       *dnnz,*onnz;
1708:   PetscScalar    *vdnz,*vonz;

1711:   mtype = dm->mattype;
1712:   PetscStrcmp(mtype,MATNEST,&isNest);
1713:   if (isNest) {
1714:     DMCreateMatrix_Network_Nest(dm,J);
1715:     return(0);
1716:   }

1718:   if (!network->userEdgeJacobian && !network->userVertexJacobian) {
1719:     /* user does not provide Jacobian blocks */
1720:     DMCreateMatrix_Plex(network->plex,J);
1721:     MatSetDM(*J,dm);
1722:     return(0);
1723:   }

1725:   MatCreate(PetscObjectComm((PetscObject)dm),J);
1726:   DMGetGlobalSection(network->plex,&sectionGlobal);
1727:   PetscSectionGetConstrainedStorageSize(sectionGlobal,&localSize);
1728:   MatSetSizes(*J,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);

1730:   MatSetType(*J,MATAIJ);
1731:   MatSetFromOptions(*J);

1733:   /* (1) Set matrix preallocation */
1734:   /*------------------------------*/
1735:   PetscObjectGetComm((PetscObject)dm,&comm);
1736:   VecCreate(comm,&vd_nz);
1737:   VecSetSizes(vd_nz,localSize,PETSC_DECIDE);
1738:   VecSetFromOptions(vd_nz);
1739:   VecSet(vd_nz,0.0);
1740:   VecDuplicate(vd_nz,&vo_nz);

1742:   /* Set preallocation for edges */
1743:   /*-----------------------------*/
1744:   DMNetworkGetEdgeRange(dm,&eStart,&eEnd);

1746:   PetscMalloc1(localSize,&rows);
1747:   for (e=eStart; e<eEnd; e++) {
1748:     /* Get row indices */
1749:     DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
1750:     DMNetworkGetNumVariables(dm,e,&nrows);
1751:     if (nrows) {
1752:       for (j=0; j<nrows; j++) rows[j] = j + rstart;

1754:       /* Set preallocation for conntected vertices */
1755:       DMNetworkGetConnectedVertices(dm,e,&cone);
1756:       for (v=0; v<2; v++) {
1757:         DMNetworkGetNumVariables(dm,cone[v],&ncols);

1759:         if (network->Je) {
1760:           Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
1761:         } else Juser = NULL;
1762:         DMNetworkIsGhostVertex(dm,cone[v],&ghost);
1763:         MatSetPreallocationblock_private(Juser,nrows,rows,ncols,ghost,vd_nz,vo_nz);
1764:       }

1766:       /* Set preallocation for edge self */
1767:       cstart = rstart;
1768:       if (network->Je) {
1769:         Juser = network->Je[3*e]; /* Jacobian(e,e) */
1770:       } else Juser = NULL;
1771:       MatSetPreallocationblock_private(Juser,nrows,rows,nrows,PETSC_FALSE,vd_nz,vo_nz);
1772:     }
1773:   }

1775:   /* Set preallocation for vertices */
1776:   /*--------------------------------*/
1777:   DMNetworkGetVertexRange(dm,&vStart,&vEnd);
1778:   if (vEnd - vStart) vptr = network->Jvptr;

1780:   for (v=vStart; v<vEnd; v++) {
1781:     /* Get row indices */
1782:     DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
1783:     DMNetworkGetNumVariables(dm,v,&nrows);
1784:     if (!nrows) continue;

1786:     DMNetworkIsGhostVertex(dm,v,&ghost);
1787:     if (ghost) {
1788:       PetscMalloc1(nrows,&rows_v);
1789:     } else {
1790:       rows_v = rows;
1791:     }

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

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

1798:     for (e=0; e<nedges; e++) {
1799:       /* Supporting edges */
1800:       DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
1801:       DMNetworkGetNumVariables(dm,edges[e],&ncols);

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

1808:       /* Connected vertices */
1809:       DMNetworkGetConnectedVertices(dm,edges[e],&cone);
1810:       vc = (v == cone[0]) ? cone[1]:cone[0];
1811:       DMNetworkIsGhostVertex(dm,vc,&ghost_vc);

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

1815:       if (network->Jv) {
1816:         Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
1817:       } else Juser = NULL;
1818:       if (ghost_vc||ghost) {
1819:         ghost2 = PETSC_TRUE;
1820:       } else {
1821:         ghost2 = PETSC_FALSE;
1822:       }
1823:       MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost2,vd_nz,vo_nz);
1824:     }

1826:     /* Set preallocation for vertex self */
1827:     DMNetworkIsGhostVertex(dm,v,&ghost);
1828:     if (!ghost) {
1829:       DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
1830:       if (network->Jv) {
1831:         Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
1832:       } else Juser = NULL;
1833:       MatSetPreallocationblock_private(Juser,nrows,rows_v,nrows,PETSC_FALSE,vd_nz,vo_nz);
1834:     }
1835:     if (ghost) {
1836:       PetscFree(rows_v);
1837:     }
1838:   }

1840:   VecAssemblyBegin(vd_nz);
1841:   VecAssemblyBegin(vo_nz);

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

1845:   VecAssemblyEnd(vd_nz);
1846:   VecAssemblyEnd(vo_nz);

1848:   VecGetArray(vd_nz,&vdnz);
1849:   VecGetArray(vo_nz,&vonz);
1850:   for (j=0; j<localSize; j++) {
1851:     dnnz[j] = (PetscInt)PetscRealPart(vdnz[j]);
1852:     onnz[j] = (PetscInt)PetscRealPart(vonz[j]);
1853:   }
1854:   VecRestoreArray(vd_nz,&vdnz);
1855:   VecRestoreArray(vo_nz,&vonz);
1856:   VecDestroy(&vd_nz);
1857:   VecDestroy(&vo_nz);

1859:   MatSeqAIJSetPreallocation(*J,0,dnnz);
1860:   MatMPIAIJSetPreallocation(*J,0,dnnz,0,onnz);
1861:   MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);

1863:   PetscFree2(dnnz,onnz);

1865:   /* (2) Set matrix entries for edges */
1866:   /*----------------------------------*/
1867:   for (e=eStart; e<eEnd; e++) {
1868:     /* Get row indices */
1869:     DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
1870:     DMNetworkGetNumVariables(dm,e,&nrows);
1871:     if (nrows) {
1872:       for (j=0; j<nrows; j++) rows[j] = j + rstart;

1874:       /* Set matrix entries for conntected vertices */
1875:       DMNetworkGetConnectedVertices(dm,e,&cone);
1876:       for (v=0; v<2; v++) {
1877:         DMNetworkGetVariableGlobalOffset(dm,cone[v],&cstart);
1878:         DMNetworkGetNumVariables(dm,cone[v],&ncols);

1880:         if (network->Je) {
1881:           Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
1882:         } else Juser = NULL;
1883:         MatSetblock_private(Juser,nrows,rows,ncols,cstart,J);
1884:       }

1886:       /* Set matrix entries for edge self */
1887:       cstart = rstart;
1888:       if (network->Je) {
1889:         Juser = network->Je[3*e]; /* Jacobian(e,e) */
1890:       } else Juser = NULL;
1891:       MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);
1892:     }
1893:   }

1895:   /* Set matrix entries for vertices */
1896:   /*---------------------------------*/
1897:   for (v=vStart; v<vEnd; v++) {
1898:     /* Get row indices */
1899:     DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
1900:     DMNetworkGetNumVariables(dm,v,&nrows);
1901:     if (!nrows) continue;

1903:     DMNetworkIsGhostVertex(dm,v,&ghost);
1904:     if (ghost) {
1905:       PetscMalloc1(nrows,&rows_v);
1906:     } else {
1907:       rows_v = rows;
1908:     }
1909:     for (j=0; j<nrows; j++) rows_v[j] = j + rstart;

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

1914:     for (e=0; e<nedges; e++) {
1915:       /* Supporting edges */
1916:       DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
1917:       DMNetworkGetNumVariables(dm,edges[e],&ncols);

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

1924:       /* Connected vertices */
1925:       DMNetworkGetConnectedVertices(dm,edges[e],&cone);
1926:       vc = (v == cone[0]) ? cone[1]:cone[0];

1928:       DMNetworkGetVariableGlobalOffset(dm,vc,&cstart);
1929:       DMNetworkGetNumVariables(dm,vc,&ncols);

1931:       if (network->Jv) {
1932:         Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
1933:       } else Juser = NULL;
1934:       MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);
1935:     }

1937:     /* Set matrix entries for vertex self */
1938:     if (!ghost) {
1939:       DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
1940:       if (network->Jv) {
1941:         Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
1942:       } else Juser = NULL;
1943:       MatSetblock_private(Juser,nrows,rows_v,nrows,cstart,J);
1944:     }
1945:     if (ghost) {
1946:       PetscFree(rows_v);
1947:     }
1948:   }
1949:   PetscFree(rows);

1951:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1952:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);

1954:   MatSetDM(*J,dm);
1955:   return(0);
1956: }

1958: PetscErrorCode DMDestroy_Network(DM dm)
1959: {
1961:   DM_Network     *network = (DM_Network*)dm->data;
1962:   PetscInt       j;

1965:   if (--network->refct > 0) return(0);
1966:   if (network->Je) {
1967:     PetscFree(network->Je);
1968:   }
1969:   if (network->Jv) {
1970:     PetscFree(network->Jvptr);
1971:     PetscFree(network->Jv);
1972:   }

1974:   ISLocalToGlobalMappingDestroy(&network->vertex.mapping);
1975:   PetscSectionDestroy(&network->vertex.DofSection);
1976:   PetscSectionDestroy(&network->vertex.GlobalDofSection);
1977:   if (network->vltog) {
1978:     PetscFree(network->vltog);
1979:   }
1980:   if (network->vertex.sf) {
1981:     PetscSFDestroy(&network->vertex.sf);
1982:   }
1983:   /* edge */
1984:   ISLocalToGlobalMappingDestroy(&network->edge.mapping);
1985:   PetscSectionDestroy(&network->edge.DofSection);
1986:   PetscSectionDestroy(&network->edge.GlobalDofSection);
1987:   if (network->edge.sf) {
1988:     PetscSFDestroy(&network->edge.sf);
1989:   }
1990:   DMDestroy(&network->plex);
1991:   PetscSectionDestroy(&network->DataSection);
1992:   PetscSectionDestroy(&network->DofSection);

1994:   for(j=0; j<network->nsubnet; j++) {
1995:     PetscFree(network->subnet[j].edges);
1996:   }
1997:   PetscFree(network->subnetvtx);

1999:   PetscFree(network->subnet);
2000:   PetscFree(network->componentdataarray);
2001:   PetscFree2(network->header,network->cvalue);
2002:   PetscFree(network);
2003:   return(0);
2004: }

2006: PetscErrorCode DMView_Network(DM dm,PetscViewer viewer)
2007: {
2009:   DM_Network     *network = (DM_Network*)dm->data;
2010:   PetscBool      iascii;
2011:   PetscMPIInt    rank;
2012:   PetscInt       p,nsubnet;

2015:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
2016:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
2019:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2020:   if (iascii) {
2021:     const PetscInt    *cone,*vtx,*edges;
2022:     PetscInt          vfrom,vto,i,j,nv,ne;

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

2028:     for (i=0; i<nsubnet; i++) {
2029:       DMNetworkGetSubnetworkInfo(dm,i,&nv,&ne,&vtx,&edges);
2030:       if (ne) {
2031:         PetscViewerASCIISynchronizedPrintf(viewer, "     Subnet %D: nEdges %D, nVertices %D\n",i,ne,nv);
2032:         for (j=0; j<ne; j++) {
2033:           p = edges[j];
2034:           DMNetworkGetConnectedVertices(dm,p,&cone);
2035:           DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);
2036:           DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);
2037:           DMNetworkGetGlobalEdgeIndex(dm,edges[j],&p);
2038:           PetscViewerASCIISynchronizedPrintf(viewer, "       edge %D: %D----> %D\n",p,vfrom,vto);
2039:         }
2040:       }
2041:     }
2042:     /* Coupling subnets */
2043:     nsubnet = network->nsubnet;
2044:     for (; i<nsubnet; i++) {
2045:       DMNetworkGetSubnetworkInfo(dm,i,&nv,&ne,&vtx,&edges);
2046:       if (ne) {
2047:         PetscViewerASCIISynchronizedPrintf(viewer, "     Subnet %D (couple): nEdges %D, nVertices %D\n",i,ne,nv);
2048:         for (j=0; j<ne; j++) {
2049:           p = edges[j];
2050:           DMNetworkGetConnectedVertices(dm,p,&cone);
2051:           DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);
2052:           DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);
2053:           PetscViewerASCIISynchronizedPrintf(viewer, "       edge %D: %D----> %D\n",p,vfrom,vto);
2054:         }
2055:       }
2056:     }
2057:     PetscViewerFlush(viewer);
2058:     PetscViewerASCIIPopSynchronized(viewer);
2059:   } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMNetwork writing", ((PetscObject)viewer)->type_name);
2060:   return(0);
2061: }

2063: PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l)
2064: {
2066:   DM_Network     *network = (DM_Network*)dm->data;

2069:   DMGlobalToLocalBegin(network->plex,g,mode,l);
2070:   return(0);
2071: }

2073: PetscErrorCode DMGlobalToLocalEnd_Network(DM dm, Vec g, InsertMode mode, Vec l)
2074: {
2076:   DM_Network     *network = (DM_Network*)dm->data;

2079:   DMGlobalToLocalEnd(network->plex,g,mode,l);
2080:   return(0);
2081: }

2083: PetscErrorCode DMLocalToGlobalBegin_Network(DM dm, Vec l, InsertMode mode, Vec g)
2084: {
2086:   DM_Network     *network = (DM_Network*)dm->data;

2089:   DMLocalToGlobalBegin(network->plex,l,mode,g);
2090:   return(0);
2091: }

2093: PetscErrorCode DMLocalToGlobalEnd_Network(DM dm, Vec l, InsertMode mode, Vec g)
2094: {
2096:   DM_Network     *network = (DM_Network*)dm->data;

2099:   DMLocalToGlobalEnd(network->plex,l,mode,g);
2100:   return(0);
2101: }

2103: /*@
2104:   DMNetworkGetVertexLocalToGlobalOrdering - Get vertex globle index

2106:   Not collective

2108:   Input Parameters
2109: + dm - the dm object
2110: - vloc - local vertex ordering, start from 0

2112:   Output Parameters
2113: +  vg  - global vertex ordering, start from 0

2115:   Level: Advanced

2117: .seealso: DMNetworkSetVertexLocalToGlobalOrdering()
2118: @*/
2119: PetscErrorCode DMNetworkGetVertexLocalToGlobalOrdering(DM dm,PetscInt vloc,PetscInt *vg)
2120: {
2121:   DM_Network  *network = (DM_Network*)dm->data;
2122:   PetscInt    *vltog = network->vltog;

2125:   if (!vltog) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkSetVertexLocalToGlobalOrdering() first");
2126:   *vg = vltog[vloc];
2127:   return(0);
2128: }

2130: /*@
2131:   DMNetworkSetVertexLocalToGlobalOrdering - Create and setup vertex local to globle map

2133:   Collective

2135:   Input Parameters:
2136: + dm - the dm object

2138:   Level: Advanced

2140: .seealso: DMNetworkGetGlobalVertexIndex()
2141: @*/
2142: PetscErrorCode DMNetworkSetVertexLocalToGlobalOrdering(DM dm)
2143: {
2144:   PetscErrorCode    ierr;
2145:   DM_Network        *network=(DM_Network*)dm->data;
2146:   MPI_Comm          comm;
2147:   PetscMPIInt       rank,size,*displs,*recvcounts,remoterank;
2148:   PetscBool         ghost;
2149:   PetscInt          *vltog,nroots,nleaves,i,*vrange,k,N,lidx;
2150:   const PetscSFNode *iremote;
2151:   PetscSF           vsf;
2152:   Vec               Vleaves,Vleaves_seq;
2153:   VecScatter        ctx;
2154:   PetscScalar       *varr,val;
2155:   const PetscScalar *varr_read;

2158:   PetscObjectGetComm((PetscObject)dm,&comm);
2159:   MPI_Comm_size(comm,&size);
2160:   MPI_Comm_rank(comm,&rank);

2162:   if (size == 1) {
2163:     nroots = network->vEnd - network->vStart;
2164:     PetscMalloc1(nroots, &vltog);
2165:     for (i=0; i<nroots; i++) vltog[i] = i;
2166:     network->vltog = vltog;
2167:     return(0);
2168:   }

2170:   if (!network->distributecalled) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkDistribute() first");
2171:   if (network->vltog) {
2172:     PetscFree(network->vltog);
2173:   }

2175:   DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);
2176:   PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);
2177:   vsf = network->vertex.sf;

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

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

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

2189:   PetscMalloc1(nroots, &vltog);
2190:   network->vltog = vltog;

2192:   /* Set vltog for non-ghost vertices */
2193:   k = 0;
2194:   for (i=0; i<nroots; i++) {
2195:     DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);
2196:     if (ghost) continue;
2197:     vltog[i] = vrange[rank] + k++;
2198:   }
2199:   PetscFree3(vrange,displs,recvcounts);

2201:   /* Set vltog for ghost vertices */
2202:   /* (a) create parallel Vleaves and sequential Vleaves_seq to convert local iremote[*].index to global index */
2203:   VecCreate(comm,&Vleaves);
2204:   VecSetSizes(Vleaves,2*nleaves,PETSC_DETERMINE);
2205:   VecSetFromOptions(Vleaves);
2206:   VecGetArray(Vleaves,&varr);
2207:   for (i=0; i<nleaves; i++) {
2208:     varr[2*i]   = (PetscScalar)(iremote[i].rank);  /* rank of remote process */
2209:     varr[2*i+1] = (PetscScalar)(iremote[i].index); /* local index in remote process */
2210:   }
2211:   VecRestoreArray(Vleaves,&varr);

2213:   /* (b) scatter local info to remote processes via VecScatter() */
2214:   VecScatterCreateToAll(Vleaves,&ctx,&Vleaves_seq);
2215:   VecScatterBegin(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);
2216:   VecScatterEnd(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);

2218:   /* (c) convert local indices to global indices in parallel vector Vleaves */
2219:   VecGetSize(Vleaves_seq,&N);
2220:   VecGetArrayRead(Vleaves_seq,&varr_read);
2221:   for (i=0; i<N; i+=2) {
2222:     remoterank = (PetscMPIInt)PetscRealPart(varr_read[i]);
2223:     if (remoterank == rank) {
2224:       k = i+1; /* row number */
2225:       lidx = (PetscInt)PetscRealPart(varr_read[i+1]);
2226:       val  = (PetscScalar)vltog[lidx]; /* global index for non-ghost vertex computed above */
2227:       VecSetValues(Vleaves,1,&k,&val,INSERT_VALUES);
2228:     }
2229:   }
2230:   VecRestoreArrayRead(Vleaves_seq,&varr_read);
2231:   VecAssemblyBegin(Vleaves);
2232:   VecAssemblyEnd(Vleaves);

2234:   /* (d) Set vltog for ghost vertices by copying local values of Vleaves */
2235:   VecGetArrayRead(Vleaves,&varr_read);
2236:   k = 0;
2237:   for (i=0; i<nroots; i++) {
2238:     DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);
2239:     if (!ghost) continue;
2240:     vltog[i] = (PetscInt)PetscRealPart(varr_read[2*k+1]); k++;
2241:   }
2242:   VecRestoreArrayRead(Vleaves,&varr_read);

2244:   VecDestroy(&Vleaves);
2245:   VecDestroy(&Vleaves_seq);
2246:   VecScatterDestroy(&ctx);
2247:   return(0);
2248: }