Actual source code: network.c

petsc-master 2019-05-18
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]; network->subnet[i].Nedge = b[1];

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

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

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

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

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

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

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

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

113:   Logically collective on DM

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

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

124:   Level: intermediate

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

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

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

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

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

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

164:   Collective on DM

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

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

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

175:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

374:   Level: intermediate

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

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

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

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

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

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

404:   Level: intermediate

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

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

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

422:   Logically collective on DM

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

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

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

435:    Level: intermediate

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

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

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

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

469:   Not Collective

471:   Input Parameters:
472: + dm - The DMNetwork object

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

478:   Level: intermediate

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

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

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

495:   Not Collective

497:   Input Parameters:
498: + dm - The DMNetwork object

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

504:   Level: intermediate

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

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

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

521:   Not Collective

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

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

530:   Level: intermediate

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

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

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

552:   Not Collective

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

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

561:   Level: intermediate

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

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

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

584:   Not Collective

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

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

595:   Notes:
596:   Typical usage:

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

606:   Level: intermediate

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

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

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

628:   Not Collective

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

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

639:   Notes:
640:   Typical usage:

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

648:   Level: intermediate

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

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

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

667:   Not Collective

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

675:   Level: intermediate

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

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

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

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

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

703:   Not Collective

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

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

712:   Level: intermediate

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

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

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

731:   Not Collective

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

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

740:   Level: intermediate

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

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

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

757:   Not Collective

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

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

766:   Level: intermediate

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

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

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

784:   Not Collective

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

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

793:   Level: intermediate

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


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

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

811:   Not Collective

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

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

820:   Level: intermediate

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


831:   p -= network->vStart;

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

839:   Not Collective

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

846:   Level: intermediate

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

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

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

863:   Not Collective

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

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

872:   Level: intermediate

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

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

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

889:   Not Collective

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

896:   Level: intermediate

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

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

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

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

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

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

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

956:   Not Collective

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

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

964:   Level: intermediate

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

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

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

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

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

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

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

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

1015:   Collective

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

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

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

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

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

1028:   Level: intermediate

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

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

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

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

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


1064:   /* Add viewers */
1065:   PetscObjectSetName((PetscObject)network->edge.GlobalDofSection,"Global edge dof section");
1066:   PetscObjectSetName((PetscObject)network->vertex.GlobalDofSection,"Global vertex dof section");
1067:   PetscSectionViewFromOptions(network->edge.GlobalDofSection, NULL, "-edge_global_section_view");
1068:   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 default section for dm */
1135:   DMSetSection(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;

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

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

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

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

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

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

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

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

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

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

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

1262:   Not Collective

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

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

1272:   Level: intermediate

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

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

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

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

1294:   Not Collective

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

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

1303:   Level: intermediate

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

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

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

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

1324:   Not Collective

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

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

1333:   Level: intermediate

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

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

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

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

1362:   DMSetSection(network->plex,network->DofSection);
1363:   DMGetGlobalSection(network->plex,&network->GlobalDofSection);

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

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

1374:     Collective

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

1381:     Level: intermediate

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

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

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

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

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

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

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

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

1423:     Not Collective

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

1432:     Level: intermediate

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

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

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

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

1454:     Not Collective

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

1464:     Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1862:   PetscFree2(dnnz,onnz);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1990:   for(j=0; j<network->nsubnet; j++) {
1991:     PetscFree(network->subnet[j].edges);
1992:   }
1993:   PetscFree(network->subnetvtx);

1995:   PetscFree(network->subnet);
1996:   PetscFree(network->componentdataarray);
1997:   PetscFree2(network->header,network->cvalue);
1998:   PetscFree(network);
1999:   return(0);
2000: }

2002: PetscErrorCode DMView_Network(DM dm,PetscViewer viewer)
2003: {
2005:   DM_Network     *network = (DM_Network*)dm->data;
2006:   PetscBool      iascii;
2007:   PetscMPIInt    rank;
2008:   PetscInt       p,nsubnet;

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

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

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

2059: PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l)
2060: {
2062:   DM_Network     *network = (DM_Network*)dm->data;

2065:   DMGlobalToLocalBegin(network->plex,g,mode,l);
2066:   return(0);
2067: }

2069: PetscErrorCode DMGlobalToLocalEnd_Network(DM dm, Vec g, InsertMode mode, Vec l)
2070: {
2072:   DM_Network     *network = (DM_Network*)dm->data;

2075:   DMGlobalToLocalEnd(network->plex,g,mode,l);
2076:   return(0);
2077: }

2079: PetscErrorCode DMLocalToGlobalBegin_Network(DM dm, Vec l, InsertMode mode, Vec g)
2080: {
2082:   DM_Network     *network = (DM_Network*)dm->data;

2085:   DMLocalToGlobalBegin(network->plex,l,mode,g);
2086:   return(0);
2087: }

2089: PetscErrorCode DMLocalToGlobalEnd_Network(DM dm, Vec l, InsertMode mode, Vec g)
2090: {
2092:   DM_Network     *network = (DM_Network*)dm->data;

2095:   DMLocalToGlobalEnd(network->plex,l,mode,g);
2096:   return(0);
2097: }