Actual source code: network.c

      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: . dm - the dm object

     11:   Output Parameters:
     12: . plexdm - the plex dm object

     14:   Level: Advanced

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

     23:   *plexdm = network->plex;
     24:   return(0);
     25: }

     27: /*@
     28:   DMNetworkGetNumSubNetworks - Gets the the number of subnetworks

     30:   Not collective

     32:   Input Parameters:
     33: . dm - the dm object

     35:   Output Parameters:
     36: + nsubnet - local number of subnetworks
     37: - Nsubnet - global number of subnetworks

     39:   Level: beginner

     41: .seealso: DMNetworkCreate(), DMNetworkSetNumSubNetworks()
     42: @*/
     43: PetscErrorCode DMNetworkGetNumSubNetworks(DM dm,PetscInt *nsubnet,PetscInt *Nsubnet)
     44: {
     45:   DM_Network *network = (DM_Network*)dm->data;

     48:   if (nsubnet) *nsubnet = network->nsubnet;
     49:   if (Nsubnet) *Nsubnet = network->Nsubnet;
     50:   return(0);
     51: }

     53: /*@
     54:   DMNetworkSetNumSubNetworks - Sets the number of subnetworks

     56:   Collective on dm

     58:   Input Parameters:
     59: + dm - the dm object
     60: . nsubnet - local number of subnetworks
     61: - Nsubnet - global number of subnetworks

     63:    Level: beginner

     65: .seealso: DMNetworkCreate(), DMNetworkGetNumSubNetworks()
     66: @*/
     67: PetscErrorCode DMNetworkSetNumSubNetworks(DM dm,PetscInt nsubnet,PetscInt Nsubnet)
     68: {
     70:   DM_Network     *network = (DM_Network*)dm->data;

     73:   if (network->Nsubnet != 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_INCOMP,"Network sizes alread set, cannot resize the network");

     79:   if (Nsubnet == PETSC_DECIDE) {
     80:     if (nsubnet < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of local subnetworks %D cannot be less than 0",nsubnet);
     81:     MPIU_Allreduce(&nsubnet,&Nsubnet,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
     82:   }
     83:   if (Nsubnet < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_INCOMP,"Number of global subnetworks %D cannot be less than 1",Nsubnet);

     85:   network->Nsubnet  = Nsubnet;
     86:   network->nsubnet  = 0;       /* initia value; will be determind by DMNetworkAddSubnetwork() */
     87:   PetscCalloc1(Nsubnet,&network->subnet);

     89:   /* num of shared vertices */
     90:   network->nsvtx = 0;
     91:   network->Nsvtx = 0;
     92:   return(0);
     93: }

     95: /*@
     96:   DMNetworkAddSubnetwork - Add a subnetwork

     98:   Collective on dm

    100:   Input Parameters:
    101: + dm - the dm object
    102: . name - name of the subnetwork
    103: . nv - number of local vertices of this subnetwork
    104: . ne - number of local edges of this subnetwork
    105: - edgelist - list of edges for this subnetwork

    107:   Output Parameters:
    108: . netnum - global index of the subnetwork

    110:   Notes:
    111:   There is no copy involved in this operation, only the pointer is referenced. The edgelist should
    112:   not be destroyed before the call to DMNetworkLayoutSetUp()

    114:   Level: beginner

    116:   Example usage:
    117:   Consider the following network:
    118: .vb
    119:  network 1: v1 -> v2 -> v0
    120: .ve

    122:  The resulting input
    123:    edgelist = [1 2 | 2 0]

    125: .seealso: DMNetworkCreate(), DMNetworkSetNumSubnetworks()
    126: @*/
    127: PetscErrorCode DMNetworkAddSubnetwork(DM dm,const char* name,PetscInt nv,PetscInt ne,PetscInt edgelist[],PetscInt *netnum)
    128: {
    130:   DM_Network     *network = (DM_Network*)dm->data;
    131:   PetscInt       i = network->nsubnet,a[2],b[2];

    134:   if (name) {
    135:     PetscStrcpy(network->subnet[i].name,name);
    136:   }

    138:   network->subnet[i].nvtx     = nv;
    139:   network->subnet[i].nedge    = ne;
    140:   network->subnet[i].edgelist = edgelist;

    142:   /* Get global number of vertices and edges for subnet[i] */
    143:   a[0] = nv; a[1] = ne;
    144:   MPIU_Allreduce(a,b,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
    145:   network->subnet[i].Nvtx  = b[0];
    146:   network->subnet[i].Nedge = b[1];

    148:   /* ----------------------------------------------------------
    149:    p=v or e;
    150:    subnet[0].pStart   = 0
    151:    subnet[i+1].pStart = subnet[i].pEnd = subnet[i].pStart + (nE[i] or NV[i])
    152:    ----------------------------------------------------------------------- */
    153:   /* GLOBAL subnet[].vStart and vEnd, used by DMNetworkLayoutSetUp() */
    154:   network->subnet[i].vStart = network->NVertices;
    155:   network->subnet[i].vEnd   = network->subnet[i].vStart + network->subnet[i].Nvtx; /* global vEnd of subnet[i] */

    157:   network->nVertices += nv;
    158:   network->NVertices += network->subnet[i].Nvtx;

    160:   /* LOCAL subnet[].eStart and eEnd, used by DMNetworkLayoutSetUp() */
    161:   network->subnet[i].eStart = network->nEdges;
    162:   network->subnet[i].eEnd   = network->subnet[i].eStart + ne;
    163:   network->nEdges += ne;
    164:   network->NEdges += network->subnet[i].Nedge;

    166:   PetscStrcpy(network->subnet[i].name,name);
    167:   if (netnum) *netnum = network->nsubnet;
    168:   network->nsubnet++;
    169:   return(0);
    170: }

    172: /*
    173:   SetUp a single svtx struct. See SVtx defined in dmnetworkimpl.h
    174:   Set gidx and type if input v=(net,idx) is a from_vertex;
    175:   Get gid, type and index in the svtx array if input v=(net,idx) is a to_vertex.

    177:   Input:  Nsvtx, svtx, net, idx, gidx
    178:   Output: gidx, svtype, svtx_idx
    179:  */
    180: static PetscErrorCode SVtxSetUp(PetscInt Nsvtx,SVtx *svtx,PetscInt net,PetscInt idx,PetscInt *gidx,SVtxType *svtype,PetscInt *svtx_idx)
    181: {
    182:   PetscInt i,j,*svto;
    183:   SVtxType vtype;

    186:   if (!Nsvtx) return(0);

    188:   vtype = SVNONE;
    189:   for (i=0; i<Nsvtx; i++) {
    190:     if (net == svtx[i].sv[0] && idx == svtx[i].sv[1]) {
    191:       /* (1) input vertex net.idx is a shared from_vertex, set its global index and output its svtype */
    192:       svtx[i].gidx = *gidx; /* set gidx */
    193:       vtype        = SVFROM;
    194:     } else { /* loop over svtx[i].n */
    195:       for (j=1; j<svtx[i].n; j++) {
    196:         svto = svtx[i].sv + 2*j;
    197:         if (net == svto[0] && idx == svto[1]) {
    198:           /* input vertex net.idx is a shared to_vertex, output its global index and its svtype */
    199:           *gidx = svtx[i].gidx; /* output gidx for to_vertex */
    200:           vtype = SVTO;
    201:         }
    202:       }
    203:     }
    204:     if (vtype != SVNONE) break;
    205:   }
    206:   if (svtype)   *svtype   = vtype;
    207:   if (svtx_idx) *svtx_idx = i;
    208:   return(0);
    209: }

    211: /*
    212:  Add a new shared vertex sv=(net,idx) to table svtas[ita]
    213: */
    214: static PetscErrorCode TableAddSVtx(PetscTable *svtas,PetscInt ita,PetscInt* tdata,PetscInt *sv_wk,PetscInt *ii,PetscInt *sedgelist,PetscInt k,DM_Network *network,PetscInt **ta2sv)
    215: {
    216:   PetscInt       net,idx,gidx,i=*ii;

    220:   net = sv_wk[2*i]   = sedgelist[k];
    221:   idx = sv_wk[2*i+1] = sedgelist[k+1];
    222:   gidx = network->subnet[net].vStart + idx;
    223:   PetscTableAdd(svtas[ita],gidx+1,tdata[ita]+1,INSERT_VALUES);
    224:   *(ta2sv[ita] + tdata[ita]) = i; /* maps tdata to index of sv_wk; sv_wk keeps (net,idx) info */
    225:   tdata[ita]++; (*ii)++;
    226:   return(0);
    227: }

    229: /*
    230:   Create an array of shared vertices. See SVtx and SVtxType in dmnetworkimpl.h

    232:   Input:  dm, Nsedgelist, sedgelist
    233:   Output: Nsvtx,svtx

    235:   Note: Output svtx is organized as
    236:         sv(net[0],idx[0]) --> sv(net[1],idx[1])
    237:                           --> sv(net[1],idx[1])
    238:                           ...
    239:                           --> sv(net[n-1],idx[n-1])
    240:         and net[0] < net[1] < ... < net[n-1]
    241:         where sv[0] has SVFROM type, sv[i], i>0, has SVTO type.
    242:  */
    243: static PetscErrorCode SVtxCreate(DM dm,PetscInt Nsedgelist,PetscInt *sedgelist,PetscInt *Nsvtx,SVtx **svtx)
    244: {
    246:   SVtx           *sedges = NULL;
    247:   PetscInt       *sv,k,j,nsv,*tdata,**ta2sv;
    248:   PetscTable     *svtas;
    249:   PetscInt       gidx,net,idx,i,nta,ita,idx_from,idx_to,n,*sv_wk;
    250:   DM_Network     *network = (DM_Network*)dm->data;
    251:   PetscTablePosition ppos;

    254:   /* (1) Crete ctables svtas */
    255:   PetscCalloc4(Nsedgelist,&svtas,Nsedgelist,&tdata,4*Nsedgelist,&sv_wk,2*Nsedgelist,&ta2sv);

    257:   j   = 0;   /* sedgelist counter */
    258:   k   = 0;   /* sedgelist vertex counter j = 4*k */
    259:   i   = 0;   /* sv_wk (vertices added to the ctables) counter */
    260:   nta = 0;   /* num of sv tables created */

    262:   /* for j=0 */
    263:   PetscTableCreate(2*Nsedgelist,network->NVertices+1,&svtas[nta]);
    264:   PetscMalloc1(2*Nsedgelist,&ta2sv[nta]);

    266:   TableAddSVtx(svtas,nta,tdata,sv_wk,&i,sedgelist,k,network,ta2sv);
    267:   TableAddSVtx(svtas,nta,tdata,sv_wk,&i,sedgelist,k+2,network,ta2sv);
    268:   nta++; k += 4;

    270:   for (j = 1; j < Nsedgelist; j++) {
    271:     for (ita = 0; ita < nta; ita++) {
    272:       /* vfrom */
    273:       net = sedgelist[k]; idx = sedgelist[k+1];
    274:       gidx = network->subnet[net].vStart + idx; /* global index of the vertex net.idx before merging shared vertices */
    275:       PetscTableFind(svtas[ita],gidx+1,&idx_from);

    277:       /* vto */
    278:       net = sedgelist[k+2]; idx = sedgelist[k+3];
    279:       gidx = network->subnet[net].vStart + idx;
    280:       PetscTableFind(svtas[ita],gidx+1,&idx_to);

    282:       if (idx_from || idx_to) { /* vfrom or vto is on table svtas[ita] */
    283:         idx_from--; idx_to--;
    284:         if (idx_from < 0) { /* vto is on svtas[ita] */
    285:           TableAddSVtx(svtas,ita,tdata,sv_wk,&i,sedgelist,k,network,ta2sv);
    286:           break;
    287:         } else if (idx_to < 0) {
    288:           TableAddSVtx(svtas,ita,tdata,sv_wk,&i,sedgelist,k+2,network,ta2sv);
    289:           break;
    290:         }
    291:       }
    292:     }

    294:     if (ita == nta) {
    295:       PetscTableCreate(2*Nsedgelist,network->NVertices+1,&svtas[nta]);
    296:       PetscMalloc1(2*Nsedgelist, &ta2sv[nta]);

    298:       TableAddSVtx(svtas,nta,tdata,sv_wk,&i,sedgelist,k,network,ta2sv);
    299:       TableAddSVtx(svtas,nta,tdata,sv_wk,&i,sedgelist,k+2,network,ta2sv);
    300:       nta++;
    301:     }
    302:     k += 4;
    303:   }

    305:   /* (2) Construct sedges from ctable
    306:      sedges: edges connect vertex sv[0]=(net[0],idx[0]) to vertices sv[k], k=1,...,n-1;
    307:      net[k], k=0, ...,n-1, are in ascending order */
    308:   PetscMalloc1(nta,&sedges);
    309:   for (nsv = 0; nsv < nta; nsv++) {
    310:     /* for a single svtx, put shared vertices in ascending order of gidx */
    311:     PetscTableGetCount(svtas[nsv],&n);
    312:     PetscCalloc1(2*n,&sv);
    313:     sedges[nsv].sv   = sv;
    314:     sedges[nsv].n    = n;
    315:     sedges[nsv].gidx = -1; /* initialization */

    317:     PetscTableGetHeadPosition(svtas[nsv],&ppos);
    318:     for (k=0; k<n; k++) { /* gidx is sorted in ascending order */
    319:       PetscTableGetNext(svtas[nsv],&ppos,&gidx,&i);
    320:       gidx--; i--;

    322:       j = ta2sv[nsv][i]; /* maps i to index of sv_wk */
    323:       sv[2*k]   = sv_wk[2*j];
    324:       sv[2*k+1] = sv_wk[2*j + 1];
    325:     }
    326:   }

    328:   for (j=0; j<nta; j++) {
    329:     PetscTableDestroy(&svtas[j]);
    330:     PetscFree(ta2sv[j]);
    331:   }
    332:   PetscFree4(svtas,tdata,sv_wk,ta2sv);

    334:   *Nsvtx = nta;
    335:   *svtx  = sedges;
    336:   return(0);
    337: }

    339: static PetscErrorCode DMNetworkLayoutSetUp_Coupling(DM dm)
    340: {
    342:   DM_Network     *network = (DM_Network*)dm->data;
    343:   PetscInt       i,j,ctr,np,*edges,*subnetvtx,vStart;
    344:   PetscInt       k,*vidxlTog,Nsv=0,Nsubnet=network->Nsubnet;
    345:   PetscInt       *sedgelist=network->sedgelist;
    346:   const PetscInt *cone;
    347:   MPI_Comm       comm;
    348:   PetscMPIInt    size,rank,*recvcounts=NULL,*displs=NULL;
    349:   PetscInt       net,idx,gidx,nmerged,e,v,vfrom,vto,*vrange,*eowners,gidx_from,net_from,sv_idx;
    350:   SVtxType       svtype = SVNONE;
    351:   SVtx           *svtx=NULL;
    352:   PetscSection   sectiong;

    355:   /* This implementation requires user input each subnet by a single processor, thus subnet[net].nvtx=subnet[net].Nvtx */
    356:   for (net=0; net<Nsubnet; net++) {
    357:     if (network->subnet[net].nvtx && network->subnet[net].nvtx != network->subnet[net].Nvtx) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"subnetwork %D local num of vertices %D != %D global num",net,network->subnet[net].nvtx,network->subnet[net].Nvtx);
    358:   }

    360:   PetscObjectGetComm((PetscObject)dm,&comm);
    361:   MPI_Comm_rank(comm,&rank);
    362:   MPI_Comm_size(comm,&size);

    364:   /* (1) Create svtx[] from sedgelist */
    365:   /* -------------------------------- */
    366:   /* Nsv: global number of SVtx; svtx: shared vertices, see SVtx in dmnetworkimpl.h */
    367:   SVtxCreate(dm,network->Nsvtx,sedgelist,&Nsv,&svtx);

    369:   /* (2) Setup svtx; Shared vto vertices are merged to their vfrom vertex with same global vetex index (gidx) */
    370:   /* -------------------------------------------------------------------------------------------------------- */
    371:   /* (2.1) compute vrage[rank]: global index of 1st local vertex in proc[rank] */
    372:   PetscMalloc3(size+1,&vrange,size,&displs,size,&recvcounts);
    373:   for (i=0; i<size; i++) {displs[i] = i; recvcounts[i] = 1;}

    375:   vrange[0] = 0;
    376:   MPI_Allgatherv(&network->nVertices,1,MPIU_INT,vrange+1,recvcounts,displs,MPIU_INT,comm);
    377:   for (i=2; i<size+1; i++) {
    378:     vrange[i] += vrange[i-1];
    379:   }

    381:   /* (2.2) Create vidxlTog: maps UN-MERGED local vertex index i to global index gidx (plex, excluding ghost vertices) */
    382:   PetscMalloc1(network->nVertices,&vidxlTog);
    383:   i = 0; gidx = 0;
    384:   nmerged = 0; /* local num of merged vertices */
    385:   network->nsvtx = 0;
    386:   for (net=0; net<Nsubnet; net++) {
    387:     for (idx=0; idx<network->subnet[net].Nvtx; idx++) {
    388:       gidx_from = gidx;
    389:       sv_idx    = -1;

    391:       SVtxSetUp(Nsv,svtx,net,idx,&gidx_from,&svtype,&sv_idx);
    392:       if (svtype == SVTO) {
    393:         if (network->subnet[net].nvtx) {/* this proc owns sv_to */
    394:           net_from = svtx[sv_idx].sv[0]; /* subnet num of its shared vertex */
    395:           if (network->subnet[net_from].nvtx == 0) {
    396:             /* this proc does not own v_from, thus a new local coupling vertex */
    397:             network->nsvtx++;
    398:           }
    399:           vidxlTog[i++] = gidx_from;
    400:           nmerged++; /* a coupling vertex -- merged */
    401:         }
    402:       } else {
    403:         if (svtype == SVFROM) {
    404:           if (network->subnet[net].nvtx) {
    405:             /* this proc owns this v_from, a new local coupling vertex */
    406:             network->nsvtx++;
    407:           }
    408:         }
    409:         if (network->subnet[net].nvtx) vidxlTog[i++] = gidx;
    410:         gidx++;
    411:       }
    412:     }
    413:   }
    414: #if defined(PETSC_USE_DEBUG)
    415:   if (i != network->nVertices) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"%D != %D nVertices",i,network->nVertices);
    416: #endif

    418:   /* (2.3) Setup svtable for querry shared vertices */
    419:   for (v=0; v<Nsv; v++) {
    420:     gidx = svtx[v].gidx;
    421:     PetscTableAdd(network->svtable,gidx+1,v+1,INSERT_VALUES);
    422:   }

    424:   /* (2.4) Shared vertices in the subnetworks are merged, update global NVertices: np = sum(local nmerged) */
    425:   MPI_Allreduce(&nmerged,&np,1,MPIU_INT,MPI_SUM,comm);
    426:   network->NVertices -= np;

    428:   PetscCalloc1(2*network->nEdges,&edges);
    429:   PetscCalloc1(network->nVertices+network->nsvtx,&network->subnetvtx);

    431:   ctr = 0;
    432:   for (net=0; net<Nsubnet; net++) {
    433:     for (j = 0; j < network->subnet[net].nedge; j++) {
    434:       /* vfrom: */
    435:       i = network->subnet[net].edgelist[2*j] + (network->subnet[net].vStart - vrange[rank]);
    436:       edges[2*ctr] = vidxlTog[i];

    438:       /* vto */
    439:       i = network->subnet[net].edgelist[2*j+1] + (network->subnet[net].vStart - vrange[rank]);
    440:       edges[2*ctr+1] = vidxlTog[i];
    441:       ctr++;
    442:     }
    443:   }
    444:   PetscFree3(vrange,displs,recvcounts);
    445:   PetscFree(vidxlTog);

    447:   /* (3) Create network->plex */
    448:   DMCreate(comm,&network->plex);
    449:   DMSetType(network->plex,DMPLEX);
    450:   DMSetDimension(network->plex,1);
    451:   if (size == 1) {
    452:     DMPlexBuildFromCellList(network->plex,network->nEdges,network->nVertices-nmerged,2,edges);
    453:   } else {
    454:     DMPlexBuildFromCellListParallel(network->plex,network->nEdges,network->nVertices-nmerged,PETSC_DECIDE,2,edges,NULL);
    455:   }

    457:   DMPlexGetChart(network->plex,&network->pStart,&network->pEnd);
    458:   DMPlexGetHeightStratum(network->plex,0,&network->eStart,&network->eEnd);
    459:   DMPlexGetHeightStratum(network->plex,1,&network->vStart,&network->vEnd);
    460:   vStart = network->vStart;

    462:   PetscSectionCreate(comm,&network->DataSection);
    463:   PetscSectionCreate(comm,&network->DofSection);
    464:   PetscSectionSetChart(network->DataSection,network->pStart,network->pEnd);
    465:   PetscSectionSetChart(network->DofSection,network->pStart,network->pEnd);

    467:   network->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
    468:   np = network->pEnd - network->pStart;
    469:   PetscCalloc2(np,&network->header,np,&network->cvalue);

    471:   /* (4) Create vidxlTog: maps MERGED plex local vertex index (including ghosts) to User's global vertex index (without merging shared vertices) */
    472:   np = network->vEnd - vStart; /* include ghost vertices */
    473:   PetscMalloc2(np,&vidxlTog,size+1,&eowners);

    475:   ctr = 0;
    476:   for (e=network->eStart; e<network->eEnd; e++) {
    477:     DMNetworkGetConnectedVertices(dm,e,&cone);
    478:     vidxlTog[cone[0] - vStart] = edges[2*ctr];
    479:     vidxlTog[cone[1] - vStart] = edges[2*ctr+1];
    480:     ctr++;
    481:   }
    482:   PetscFree(edges);

    484:   /* (5) Create vertices and edges array for the subnetworks */
    485:   subnetvtx = network->subnetvtx;
    486:   for (j=0; j < Nsubnet; j++) {
    487:     PetscCalloc1(network->subnet[j].nedge,&network->subnet[j].edges);
    488:     network->subnet[j].vertices = subnetvtx;
    489:     subnetvtx                  += network->subnet[j].nvtx;
    490:   }
    491:   network->svertices                = subnetvtx;

    493:   /* Get edge ownership */
    494:   np = network->eEnd - network->eStart; /* num of local edges */
    495:   MPI_Allgather(&np,1,MPIU_INT,eowners+1,1,MPIU_INT,comm);
    496:   eowners[0] = 0;
    497:   for (i=2; i<=size; i++) eowners[i] += eowners[i-1];

    499:   e = 0;
    500:   for (i=0; i < Nsubnet; i++) {
    501:     v = 0;
    502:     for (j = 0; j < network->subnet[i].nedge; j++) {

    504:       /* edge e */
    505:       network->header[e].index    = e + eowners[rank]; /* Global edge index */
    506:       network->header[e].subnetid = i;                 /* Subnetwork id */
    507:       network->subnet[i].edges[j] = e;
    508:       network->header[e].ndata           = 0;
    509:       network->header[e].offset[0]       = 0;
    510:       network->header[e].offsetvarrel[0] = 0;
    511:       PetscSectionAddDof(network->DataSection,e,network->dataheadersize);

    513:       /* connected vertices */
    514:       DMPlexGetCone(network->plex,e,&cone);

    516:       /* vertex cone[0] */
    517:       vfrom = network->subnet[i].edgelist[2*v];
    518:       network->header[cone[0]].index = vidxlTog[cone[0]-vStart]; /*  Global vertex index */
    519:       network->header[cone[0]].subnetid = i;                     /* Subnetwork id */
    520:       network->subnet[i].vertices[vfrom] = cone[0];              /* user's subnet[].dix = petsc's v */

    522:       /* vertex cone[1] */
    523:       vto = network->subnet[i].edgelist[2*v+1];
    524:       network->header[cone[1]].index = vidxlTog[cone[1]-vStart]; /*  Global vertex index */
    525:       network->header[cone[1]].subnetid   = i;
    526:       network->subnet[i].vertices[vto]= cone[1];

    528:       e++; v++;
    529:     }
    530:   }

    532:   /* Set vertex array for the subnetworks */
    533:   k = 0;
    534:   for (v=vStart; v<network->vEnd; v++) { /* local vertices, including ghosts */
    535:     network->header[v].ndata           = 0;
    536:     network->header[v].offset[0]       = 0;
    537:     network->header[v].offsetvarrel[0] = 0;
    538:     PetscSectionAddDof(network->DataSection,v,network->dataheadersize);

    540:     /* shared vertex */
    541:     PetscTableFind(network->svtable,vidxlTog[v-vStart]+1,&i);
    542:     if (i) network->svertices[k++] = v;
    543:   }
    544: #if defined(PETSC_USE_DEBUG)
    545:   if (k != network->nsvtx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"k %D != %D nsvtx",k,network->nsvtx);
    546: #endif

    548:   PetscFree2(vidxlTog,eowners);

    550:   network->svtx  = svtx;
    551:   network->Nsvtx = Nsv;
    552:   PetscFree(sedgelist);

    554:   /* Create a global section to be used by DMNetworkIsGhostVertex() which is a non-collective routine */
    555:   DMGetGlobalSection(network->plex,&sectiong);
    556:   return(0);
    557: }

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

    562:   Collective on dm

    564:   Input Parameters:
    565: . DM - the dmnetwork object

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

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

    573:   Level: beginner

    575: .seealso: DMNetworkSetNumSubNetworks(), DMNetworkAddSubnetwork()
    576: @*/
    577: PetscErrorCode DMNetworkLayoutSetUp(DM dm)
    578: {
    580:   DM_Network     *network = (DM_Network*)dm->data;
    581:   PetscInt       i,j,ctr,Nsubnet=network->Nsubnet,*eowners,np,*edges,*subnetvtx;
    582:   PetscInt       e,v,vfrom,vto;
    583:   const PetscInt *cone;
    584:   MPI_Comm       comm;
    585:   PetscMPIInt    size,rank;

    588:   if (network->nsubnet != network->Nsubnet) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Must call DMNetworkAddSubnetwork() %D times",network->Nsubnet);

    590:   /* Create svtable for querry shared vertices */
    591:   PetscTableCreate(network->Nsvtx,network->NVertices+1,&network->svtable);

    593:   if (network->Nsvtx) {
    594:     DMNetworkLayoutSetUp_Coupling(dm);
    595:     return(0);
    596:   }

    598:   PetscObjectGetComm((PetscObject)dm,&comm);
    599:   MPI_Comm_rank(comm,&rank);
    600:   MPI_Comm_size(comm,&size);

    602:   /* Create LOCAL edgelist for the network by concatenating local input edgelists of the subnetworks */
    603:   PetscCalloc1(2*network->nEdges,&edges);
    604:   ctr = 0;
    605:   for (i=0; i < Nsubnet; i++) {
    606:     for (j = 0; j < network->subnet[i].nedge; j++) {
    607:       edges[2*ctr]   = network->subnet[i].vStart + network->subnet[i].edgelist[2*j];
    608:       edges[2*ctr+1] = network->subnet[i].vStart + network->subnet[i].edgelist[2*j+1];
    609:       ctr++;
    610:     }
    611:   }

    613:   /* Create network->plex; One dimensional network, numCorners=2 */
    614:   DMCreate(comm,&network->plex);
    615:   DMSetType(network->plex,DMPLEX);
    616:   DMSetDimension(network->plex,1);
    617:   if (size == 1) {
    618:     DMPlexBuildFromCellList(network->plex,network->nEdges,network->nVertices,2,edges);
    619:   } else {
    620:     DMPlexBuildFromCellListParallel(network->plex,network->nEdges,network->nVertices,PETSC_DECIDE,2,edges,NULL);
    621:   }
    622:   PetscFree(edges); /* local edge list with global idx used by DMPlexBuildFromCellList() */

    624:   DMPlexGetChart(network->plex,&network->pStart,&network->pEnd);
    625:   DMPlexGetHeightStratum(network->plex,0,&network->eStart,&network->eEnd);
    626:   DMPlexGetHeightStratum(network->plex,1,&network->vStart,&network->vEnd);

    628:   PetscSectionCreate(comm,&network->DataSection);
    629:   PetscSectionCreate(comm,&network->DofSection);
    630:   PetscSectionSetChart(network->DataSection,network->pStart,network->pEnd);
    631:   PetscSectionSetChart(network->DofSection,network->pStart,network->pEnd);

    633:   network->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
    634:   np = network->pEnd - network->pStart;
    635:   PetscCalloc2(np,&network->header,np,&network->cvalue);

    637:   /* Create edge and vertex arrays for the subnetworks */
    638:   for (j=0; j < network->Nsubnet; j++) {
    639:     PetscCalloc1(network->subnet[j].nedge,&network->subnet[j].edges);
    640:   }

    642:   /* Get edge ownership */
    643:   PetscMalloc1(size+1,&eowners);
    644:   np = network->eEnd - network->eStart;
    645:   MPI_Allgather(&np,1,MPIU_INT,eowners+1,1,MPIU_INT,comm);
    646:   eowners[0] = 0;
    647:   for (i=2; i<=size; i++) eowners[i] += eowners[i-1];

    649:   /* Set network->subnet[*].vertices on array network->subnetvtx */
    650:   np = 0;
    651:   for (j=0; j<network->Nsubnet; j++) {
    652:     /* sum up subnet[j].Nvtx instead of subnet[j].nvtx, because a subnet might be owned by more than one processor;
    653:        below, subnet[i].vertices[vfrom/vto] requires vfrom/vto =0, ...,Nvtx-1
    654:      */
    655:     if (network->subnet[j].nvtx) np += network->subnet[j].Nvtx;
    656:   }

    658:   PetscCalloc1(np,&network->subnetvtx); /* Maps local vertex to local subnetwork's vertex */
    659:   subnetvtx = network->subnetvtx;
    660:   for (j=0; j<network->Nsubnet; j++) {
    661:     network->subnet[j].vertices = subnetvtx;
    662:     if (network->subnet[j].nvtx) subnetvtx += network->subnet[j].Nvtx;
    663:   }

    665:   /* Setup edge and vertex arrays for subnetworks */
    666:   e = 0;
    667:   for (i=0; i < Nsubnet; i++) {
    668:     v = 0;
    669:     for (j = 0; j < network->subnet[i].nedge; j++) {
    670:       /* edge e */
    671:       network->header[e].index    = e + eowners[rank];   /* Global edge index */
    672:       network->header[e].subnetid = i;
    673:       network->subnet[i].edges[j] = e;

    675:       network->header[e].ndata           = 0;
    676:       network->header[e].offset[0]       = 0;
    677:       network->header[e].offsetvarrel[0] = 0;
    678:       PetscSectionAddDof(network->DataSection,e,network->dataheadersize);

    680:       /* connected vertices */
    681:       DMPlexGetCone(network->plex,e,&cone);

    683:       /* vertex cone[0] */
    684:       vfrom = network->subnet[i].edgelist[2*v];     /* =subnet[i].idx, Global index! */
    685:       network->header[cone[0]].index     = vfrom + network->subnet[i].vStart; /* Global vertex index */
    686:       network->header[cone[0]].subnetid  = i;       /* Subnetwork id */
    687:       network->subnet[i].vertices[vfrom] = cone[0]; /* user's subnet[].dix = petsc's v */

    689:       /* vertex cone[1] */
    690:       vto   = network->subnet[i].edgelist[2*v+1];   /* =subnet[i].idx, Global index! */
    691:       network->header[cone[1]].index    = vto + network->subnet[i].vStart;  /* Global vertex index */
    692:       network->header[cone[1]].subnetid = i;
    693:       network->subnet[i].vertices[vto]  = cone[1];  /* user's subnet[].dix = petsc's v */

    695:       e++; v++;
    696:     }
    697:   }
    698:   PetscFree(eowners);

    700:   for (v = network->vStart; v < network->vEnd; v++) {
    701:     network->header[v].ndata           = 0;
    702:     network->header[v].offset[0]       = 0;
    703:     network->header[v].offsetvarrel[0] = 0;
    704:     PetscSectionAddDof(network->DataSection,v,network->dataheadersize);
    705:   }
    706:   return(0);
    707: }

    709: /*@C
    710:   DMNetworkGetSubnetwork - Returns the information about a requested subnetwork

    712:   Not collective

    714:   Input Parameters:
    715: + dm - the DM object
    716: - netnum - the global index of the subnetwork

    718:   Output Parameters:
    719: + nv - number of vertices (local)
    720: . ne - number of edges (local)
    721: . vtx - local vertices of the subnetwork
    722: - edge - local edges of the subnetwork

    724:   Notes:
    725:   Cannot call this routine before DMNetworkLayoutSetup()

    727:   Level: intermediate

    729: .seealso: DMNetworkCreate(), DMNetworkAddSubnetwork(), DMNetworkLayoutSetUp()
    730: @*/
    731: PetscErrorCode DMNetworkGetSubnetwork(DM dm,PetscInt netnum,PetscInt *nv,PetscInt *ne,const PetscInt **vtx,const PetscInt **edge)
    732: {
    733:   DM_Network *network = (DM_Network*)dm->data;

    736:   if (netnum >= network->Nsubnet) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subnet index %D exceeds the num of subnets %D",netnum,network->Nsubnet);
    737:   if (nv) *nv     = network->subnet[netnum].nvtx;
    738:   if (ne) *ne     = network->subnet[netnum].nedge;
    739:   if (vtx) *vtx   = network->subnet[netnum].vertices;
    740:   if (edge) *edge = network->subnet[netnum].edges;
    741:   return(0);
    742: }

    744: /*@
    745:   DMNetworkAddSharedVertices - Add shared vertices that connect two given subnetworks

    747:   Collective on dm

    749:   Input Parameters:
    750: + dm - the dm object
    751: . anetnum - first subnetwork global numbering returned by DMNetworkAddSubnetwork()
    752: . bnetnum - second subnetwork global numbering returned by DMNetworkAddSubnetwork()
    753: . nsvtx - number of vertices that are shared by the two subnetworks
    754: . asvtx - vertex index in the first subnetwork
    755: - bsvtx - vertex index in the second subnetwork

    757:   Level: beginner

    759: .seealso: DMNetworkCreate(), DMNetworkAddSubnetwork(), DMNetworkGetSharedVertices()
    760: @*/
    761: PetscErrorCode DMNetworkAddSharedVertices(DM dm,PetscInt anetnum,PetscInt bnetnum,PetscInt nsvtx,PetscInt asvtx[],PetscInt bsvtx[])
    762: {
    764:   DM_Network     *network = (DM_Network*)dm->data;
    765:   PetscInt       i,nsubnet = network->Nsubnet,*sedgelist,Nsvtx=network->Nsvtx;

    768:   if (anetnum == bnetnum) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Subnetworks must have different netnum");
    769:   if (anetnum < 0 || bnetnum < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"netnum cannot be negative");
    770:   if (!Nsvtx) {
    771:     /* allocate network->sedgelist to hold at most 2*nsubnet pairs of shared vertices */
    772:     PetscMalloc1(2*4*nsubnet,&network->sedgelist);
    773:   }

    775:   sedgelist = network->sedgelist;
    776:   for (i=0; i<nsvtx; i++) {
    777:     sedgelist[4*Nsvtx]   = anetnum; sedgelist[4*Nsvtx+1] = asvtx[i];
    778:     sedgelist[4*Nsvtx+2] = bnetnum; sedgelist[4*Nsvtx+3] = bsvtx[i];
    779:     Nsvtx++;
    780:   }
    781:   if (Nsvtx > 2*nsubnet) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"allocate more space for coupling edgelist");
    782:   network->Nsvtx = Nsvtx;
    783:   return(0);
    784: }

    786: /*@C
    787:   DMNetworkGetSharedVertices - Returns the info for the shared vertices

    789:   Not collective

    791:   Input Parameters:
    792: . dm - the DM object

    794:   Output Parameters:
    795: + nsv - number of local shared vertices
    796: - svtx - local shared vertices

    798:   Notes:
    799:   Cannot call this routine before DMNetworkLayoutSetup()

    801:   Level: intermediate

    803: .seealso: DMNetworkGetSubnetwork(), DMNetworkLayoutSetUp(), DMNetworkAddSharedVertices()
    804: @*/
    805: PetscErrorCode DMNetworkGetSharedVertices(DM dm,PetscInt *nsv,const PetscInt **svtx)
    806: {
    807:   DM_Network *net = (DM_Network*)dm->data;

    810:   if (net->Nsvtx) {
    811:     *nsv  = net->nsvtx;
    812:     *svtx = net->svertices;
    813:   } else {
    814:     *nsv  = 0;
    815:     *svtx = NULL;
    816:   }
    817:   return(0);
    818: }

    820: /*@C
    821:   DMNetworkRegisterComponent - Registers the network component

    823:   Logically collective on dm

    825:   Input Parameters:
    826: + dm - the network object
    827: . name - the component name
    828: - size - the storage size in bytes for this component data

    830:    Output Parameters:
    831: .  key - an integer key that defines the component

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

    836:    Level: beginner

    838: .seealso: DMNetworkCreate(), DMNetworkLayoutSetUp()
    839: @*/
    840: PetscErrorCode DMNetworkRegisterComponent(DM dm,const char *name,size_t size,PetscInt *key)
    841: {
    842:   PetscErrorCode        ierr;
    843:   DM_Network            *network = (DM_Network*) dm->data;
    844:   DMNetworkComponent    *component=&network->component[network->ncomponent];
    845:   PetscBool             flg=PETSC_FALSE;
    846:   PetscInt              i;

    849:   for (i=0; i < network->ncomponent; i++) {
    850:     PetscStrcmp(component->name,name,&flg);
    851:     if (flg) {
    852:       *key = i;
    853:       return(0);
    854:     }
    855:   }
    856:   if (network->ncomponent == MAX_NETCOMPONENTS) {
    857:     SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Number of components registered exceeds the max %D",MAX_NETCOMPONENTS);
    858:   }

    860:   PetscStrcpy(component->name,name);
    861:   component->size = size/sizeof(DMNetworkComponentGenericDataType);
    862:   *key = network->ncomponent;
    863:   network->ncomponent++;
    864:   return(0);
    865: }

    867: /*@
    868:   DMNetworkGetVertexRange - Get the bounds [start, end) for the local vertices

    870:   Not Collective

    872:   Input Parameters:
    873: . dm - the DMNetwork object

    875:   Output Parameters:
    876: + vStart - the first vertex point
    877: - vEnd - one beyond the last vertex point

    879:   Level: beginner

    881: .seealso: DMNetworkGetEdgeRange()
    882: @*/
    883: PetscErrorCode DMNetworkGetVertexRange(DM dm,PetscInt *vStart,PetscInt *vEnd)
    884: {
    885:   DM_Network *network = (DM_Network*)dm->data;

    888:   if (vStart) *vStart = network->vStart;
    889:   if (vEnd) *vEnd = network->vEnd;
    890:   return(0);
    891: }

    893: /*@
    894:   DMNetworkGetEdgeRange - Get the bounds [start, end) for the local edges

    896:   Not Collective

    898:   Input Parameters:
    899: . dm - the DMNetwork object

    901:   Output Parameters:
    902: + eStart - The first edge point
    903: - eEnd - One beyond the last edge point

    905:   Level: beginner

    907: .seealso: DMNetworkGetVertexRange()
    908: @*/
    909: PetscErrorCode DMNetworkGetEdgeRange(DM dm,PetscInt *eStart,PetscInt *eEnd)
    910: {
    911:   DM_Network *network = (DM_Network*)dm->data;

    914:   if (eStart) *eStart = network->eStart;
    915:   if (eEnd) *eEnd = network->eEnd;
    916:   return(0);
    917: }

    919: static PetscErrorCode DMNetworkGetIndex(DM dm,PetscInt p,PetscInt *index)
    920: {
    921:   PetscErrorCode    ierr;
    922:   DM_Network        *network = (DM_Network*)dm->data;
    923:   PetscInt          offsetp;
    924:   DMNetworkComponentHeader header;

    927:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
    928:   PetscSectionGetOffset(network->DataSection,p,&offsetp);
    929:   header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
    930:   *index = header->index;
    931:   return(0);
    932: }

    934: /*@
    935:   DMNetworkGetGlobalEdgeIndex - Get the global numbering for the edge on the network

    937:   Not Collective

    939:   Input Parameters:
    940: + dm - DMNetwork object
    941: - p - edge point

    943:   Output Parameters:
    944: . index - the global numbering for the edge

    946:   Level: intermediate

    948: .seealso: DMNetworkGetGlobalVertexIndex()
    949: @*/
    950: PetscErrorCode DMNetworkGetGlobalEdgeIndex(DM dm,PetscInt p,PetscInt *index)
    951: {

    955:   DMNetworkGetIndex(dm,p,index);
    956:   return(0);
    957: }

    959: /*@
    960:   DMNetworkGetGlobalVertexIndex - Get the global numbering for the vertex on the network

    962:   Not Collective

    964:   Input Parameters:
    965: + dm - DMNetwork object
    966: - p  - vertex point

    968:   Output Parameters:
    969: . index - the global numbering for the vertex

    971:   Level: intermediate

    973: .seealso: DMNetworkGetGlobalEdgeIndex()
    974: @*/
    975: PetscErrorCode DMNetworkGetGlobalVertexIndex(DM dm,PetscInt p,PetscInt *index)
    976: {

    980:   DMNetworkGetIndex(dm,p,index);
    981:   return(0);
    982: }

    984: /*@
    985:   DMNetworkGetNumComponents - Get the number of components at a vertex/edge

    987:   Not Collective

    989:   Input Parameters:
    990: + dm - the DMNetwork object
    991: - p - vertex/edge point

    993:   Output Parameters:
    994: . numcomponents - Number of components at the vertex/edge

    996:   Level: beginner

    998: .seealso: DMNetworkRegisterComponent(), DMNetworkAddComponent()
    999: @*/
    1000: PetscErrorCode DMNetworkGetNumComponents(DM dm,PetscInt p,PetscInt *numcomponents)
    1001: {
    1003:   PetscInt       offset;
    1004:   DM_Network     *network = (DM_Network*)dm->data;

    1007:   PetscSectionGetOffset(network->DataSection,p,&offset);
    1008:   *numcomponents = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata;
    1009:   return(0);
    1010: }

    1012: /*@
    1013:   DMNetworkGetLocalVecOffset - Get the offset for accessing the variables associated with a component at the given vertex/edge from the local vector

    1015:   Not Collective

    1017:   Input Parameters:
    1018: + dm - the DMNetwork object
    1019: . p - the edge/vertex point
    1020: - compnum - component number; use ALL_COMPONENTS if no specific component is requested

    1022:   Output Parameters:
    1023: . offset - the local offset

    1025:   Level: intermediate

    1027: .seealso: DMGetLocalVector(), DMNetworkGetComponent(), DMNetworkGetGlobalVecOffset()
    1028: @*/
    1029: PetscErrorCode DMNetworkGetLocalVecOffset(DM dm,PetscInt p,PetscInt compnum,PetscInt *offset)
    1030: {
    1032:   DM_Network     *network = (DM_Network*)dm->data;
    1033:   PetscInt       offsetp,offsetd;
    1034:   DMNetworkComponentHeader header;

    1037:   PetscSectionGetOffset(network->plex->localSection,p,&offsetp);
    1038:   if (compnum == ALL_COMPONENTS) {
    1039:     *offset = offsetp;
    1040:     return(0);
    1041:   }

    1043:   PetscSectionGetOffset(network->DataSection,p,&offsetd);
    1044:   header = (DMNetworkComponentHeader)(network->componentdataarray+offsetd);
    1045:   *offset = offsetp + header->offsetvarrel[compnum];
    1046:   return(0);
    1047: }

    1049: /*@
    1050:   DMNetworkGetGlobalVecOffset - Get the global offset for accessing the variables associated with a component for the given vertex/edge from the global vector

    1052:   Not Collective

    1054:   Input Parameters:
    1055: + dm - the DMNetwork object
    1056: . p - the edge/vertex point
    1057: - compnum - component number; use ALL_COMPONENTS if no specific component is requested

    1059:   Output Parameters:
    1060: . offsetg - the global offset

    1062:   Level: intermediate

    1064: .seealso: DMNetworkGetLocalVecOffset(), DMGetGlobalVector(), DMNetworkGetComponent()
    1065: @*/
    1066: PetscErrorCode DMNetworkGetGlobalVecOffset(DM dm,PetscInt p,PetscInt compnum,PetscInt *offsetg)
    1067: {
    1069:   DM_Network     *network = (DM_Network*)dm->data;
    1070:   PetscInt       offsetp,offsetd;
    1071:   DMNetworkComponentHeader header;

    1074:   PetscSectionGetOffset(network->plex->globalSection,p,&offsetp);
    1075:   if (offsetp < 0) offsetp = -(offsetp + 1); /* Convert to actual global offset for ghost vertex */

    1077:   if (compnum == ALL_COMPONENTS) {
    1078:     *offsetg = offsetp;
    1079:     return(0);
    1080:   }
    1081:   PetscSectionGetOffset(network->DataSection,p,&offsetd);
    1082:   header = (DMNetworkComponentHeader)(network->componentdataarray+offsetd);
    1083:   *offsetg = offsetp + header->offsetvarrel[compnum];
    1084:   return(0);
    1085: }

    1087: /*@
    1088:   DMNetworkGetEdgeOffset - Get the offset for accessing the variables associated with the given edge from the local subvector

    1090:   Not Collective

    1092:   Input Parameters:
    1093: + dm - the DMNetwork object
    1094: - p - the edge point

    1096:   Output Parameters:
    1097: . offset - the offset

    1099:   Level: intermediate

    1101: .seealso: DMNetworkGetLocalVecOffset(), DMGetLocalVector()
    1102: @*/
    1103: PetscErrorCode DMNetworkGetEdgeOffset(DM dm,PetscInt p,PetscInt *offset)
    1104: {
    1106:   DM_Network     *network = (DM_Network*)dm->data;

    1109:   PetscSectionGetOffset(network->edge.DofSection,p,offset);
    1110:   return(0);
    1111: }

    1113: /*@
    1114:   DMNetworkGetVertexOffset - Get the offset for accessing the variables associated with the given vertex from the local subvector

    1116:   Not Collective

    1118:   Input Parameters:
    1119: + dm - the DMNetwork object
    1120: - p - the vertex point

    1122:   Output Parameters:
    1123: . offset - the offset

    1125:   Level: intermediate

    1127: .seealso: DMNetworkGetEdgeOffset(), DMGetLocalVector()
    1128: @*/
    1129: PetscErrorCode DMNetworkGetVertexOffset(DM dm,PetscInt p,PetscInt *offset)
    1130: {
    1132:   DM_Network     *network = (DM_Network*)dm->data;

    1135:   p -= network->vStart;
    1136:   PetscSectionGetOffset(network->vertex.DofSection,p,offset);
    1137:   return(0);
    1138: }

    1140: /*@
    1141:   DMNetworkAddComponent - Adds a network component and number of variables at the given point (vertex/edge)

    1143:   Not Collective

    1145:   Input Parameters:
    1146: + dm - the DMNetworkObject
    1147: . p - the vertex/edge point
    1148: . componentkey - component key returned while registering the component; ignored if compvalue=NULL
    1149: . compvalue - pointer to the data structure for the component, or NULL if not required.
    1150: - nvar - number of variables for the component at the vertex/edge point

    1152:   Level: beginner

    1154: .seealso: DMNetworkGetComponent()
    1155: @*/
    1156: PetscErrorCode DMNetworkAddComponent(DM dm,PetscInt p,PetscInt componentkey,void* compvalue,PetscInt nvar)
    1157: {
    1158:   PetscErrorCode           ierr;
    1159:   DM_Network               *network = (DM_Network*)dm->data;
    1160:   DMNetworkComponent       *component = &network->component[componentkey];
    1161:   DMNetworkComponentHeader header = &network->header[p];
    1162:   DMNetworkComponentValue  cvalue = &network->cvalue[p];
    1163:   PetscBool                sharedv=PETSC_FALSE;
    1164:   PetscInt                 compnum=header->ndata;

    1167:   PetscSectionAddDof(network->DofSection,p,nvar);
    1168:   if (!compvalue) return(0);

    1170:   DMNetworkIsSharedVertex(dm,p,&sharedv);
    1171:   if (sharedv) {
    1172:     PetscBool ghost;
    1173:     DMNetworkIsGhostVertex(dm,p,&ghost);
    1174:     if (ghost) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Adding a component at a leaf(ghost) shared vertex is not supported");
    1175:   }

    1177:   if (compnum == MAX_NETCOMPONENTS) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Number of components at a point exceeds the max %D",MAX_NETCOMPONENTS);

    1179:   header->size[compnum] = component->size;
    1180:   PetscSectionAddDof(network->DataSection,p,component->size);
    1181:   header->key[compnum] = componentkey;
    1182:   if (compnum != 0) header->offset[compnum] = header->offset[compnum-1] + header->size[compnum-1];
    1183:   cvalue->data[compnum] = (void*)compvalue;

    1185:   /* variables */
    1186:   header->nvar[compnum] += nvar;
    1187:   if (compnum != 0) header->offsetvarrel[compnum] = header->offsetvarrel[compnum-1] + header->nvar[compnum-1];

    1189:   header->ndata++;
    1190:   return(0);
    1191: }

    1193: /*@
    1194:   DMNetworkGetComponent - Gets the component key, the component data, and the number of variables at a given network point

    1196:   Not Collective

    1198:   Input Parameters:
    1199: + dm - the DMNetwork object
    1200: . p - vertex/edge point
    1201: - compnum - component number; use ALL_COMPONENTS if sum up all the components

    1203:   Output Parameters:
    1204: + compkey - the key obtained when registering the component (use NULL if not required)
    1205: . component - the component data (use NULL if not required)
    1206: - nvar  - number of variables (use NULL if not required)

    1208:   Level: beginner

    1210: .seealso: DMNetworkAddComponent(), DMNetworkGetNumComponents()
    1211: @*/
    1212: PetscErrorCode DMNetworkGetComponent(DM dm,PetscInt p,PetscInt compnum,PetscInt *compkey,void **component,PetscInt *nvar)
    1213: {
    1215:   DM_Network     *network = (DM_Network*)dm->data;
    1216:   PetscInt       offset = 0;
    1217:   DMNetworkComponentHeader header;

    1220:   if (compnum == ALL_COMPONENTS) {
    1221:     PetscSectionGetDof(network->DofSection,p,nvar);
    1222:     return(0);
    1223:   }

    1225:   PetscSectionGetOffset(network->DataSection,p,&offset);
    1226:   header = (DMNetworkComponentHeader)(network->componentdataarray+offset);

    1228:   if (compnum >= 0) {
    1229:     if (compkey) *compkey = header->key[compnum];
    1230:     if (component) {
    1231:       offset += network->dataheadersize+header->offset[compnum];
    1232:       *component = network->componentdataarray+offset;
    1233:     }
    1234:   }

    1236:   if (nvar) *nvar = header->nvar[compnum];
    1237:   return(0);
    1238: }

    1240: /*
    1241:  Sets up the array that holds the data for all components and its associated section.
    1242:  It copies the data for all components in a contiguous array called componentdataarray. The component data is stored pointwise with an additional header (metadata) stored for each point. The header has metadata information such as number of components at each point, number of variables for each component, offsets for the components data, etc.
    1243: */
    1244: PetscErrorCode DMNetworkComponentSetUp(DM dm)
    1245: {
    1246:   PetscErrorCode           ierr;
    1247:   DM_Network               *network = (DM_Network*)dm->data;
    1248:   PetscInt                 arr_size,p,offset,offsetp,ncomp,i;
    1249:   MPI_Comm                 comm;
    1250:   PetscMPIInt              size,rank;
    1251:   DMNetworkComponentHeader header;
    1252:   DMNetworkComponentValue  cvalue;
    1253:   DMNetworkComponentGenericDataType *componentdataarray;

    1256:   PetscObjectGetComm((PetscObject)dm,&comm);
    1257:   MPI_Comm_size(comm,&size);
    1258:   MPI_Comm_rank(comm,&rank);
    1259: #if 0
    1260:   //------------- new
    1261:   if (size > 1) { /* Sync nvar at shared vertices for all processes */
    1262:     PetscSF        sf = network->plex->sf;
    1263:     const PetscInt *degree;
    1264:     PetscInt       i,nleaves_total,*indata,*outdata,nroots,nleaves,nsv,p,ncomp;
    1265:     const PetscInt *svtx;
    1266:     PetscBool      ghost;

    1268:     PetscSFGetGraph(sf,&nroots,&nleaves,NULL,NULL);
    1269:     PetscSFComputeDegreeBegin(sf,&degree);
    1270:     PetscSFComputeDegreeEnd(sf,&degree);
    1271:     nleaves_total=0;
    1272:     for (i=0; i<nroots; i++) nleaves_total += degree[i];
    1273:     printf("[%d] nleaves_total %d\n",rank,nleaves_total);
    1274:     MPI_Barrier(comm);

    1276:     PetscCalloc2(nleaves_total,&indata,nleaves,&outdata);

    1278:     /* Leaves copy user's ncomp to outdata */
    1279:     DMNetworkGetSharedVertices(dm,&nsv,&svtx);
    1280:     for (i=0; i<nsv; i++) {
    1281:       p = svtx[i];
    1282:       DMNetworkIsGhostVertex(dm,p,&ghost);
    1283:       if (!ghost) continue;

    1285:       header = &network->header[p];
    1286:       ncomp = header->ndata;
    1287:       printf("[%d] leaf has ncomp %d\n",rank,ncomp);
    1288:       outdata[p] = ncomp;
    1289:     }

    1291:     /* Roots gather ncomp from leaves */
    1292:     PetscSFGatherBegin(sf,MPIU_INT,outdata,indata);
    1293:     PetscSFGatherEnd(sf,MPIU_INT,outdata,indata);
    1294:     PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"## Gathered data at multi-roots from leaves\n");
    1295:     PetscIntView(nleaves_total,indata,PETSC_VIEWER_STDOUT_WORLD);

    1297:     PetscFree2(indata,outdata);
    1298:   }
    1299:   //----------------------
    1300: #endif

    1302:   PetscSectionSetUp(network->DataSection);
    1303:   PetscSectionGetStorageSize(network->DataSection,&arr_size);
    1304:   PetscMalloc1(arr_size,&network->componentdataarray);
    1305:   componentdataarray = network->componentdataarray;
    1306:   for (p = network->pStart; p < network->pEnd; p++) {
    1307:     PetscSectionGetOffset(network->DataSection,p,&offsetp);
    1308:     /* Copy header */
    1309:     header = &network->header[p];
    1310:     PetscMemcpy(componentdataarray+offsetp,header,network->dataheadersize*sizeof(DMNetworkComponentGenericDataType));
    1311:     /* Copy data */
    1312:     cvalue = &network->cvalue[p];
    1313:     ncomp = header->ndata;

    1315:     for (i = 0; i < ncomp; i++) {
    1316:       offset = offsetp + network->dataheadersize + header->offset[i];
    1317:       PetscMemcpy(componentdataarray+offset,cvalue->data[i],header->size[i]*sizeof(DMNetworkComponentGenericDataType));
    1318:     }
    1319:   }
    1320:   return(0);
    1321: }

    1323: /* Sets up the section for dofs. This routine is called during DMSetUp() */
    1324: static PetscErrorCode DMNetworkVariablesSetUp(DM dm)
    1325: {
    1327:   DM_Network     *network = (DM_Network*)dm->data;
    1328:   MPI_Comm       comm;
    1329:   PetscMPIInt    size;

    1332:   PetscObjectGetComm((PetscObject)dm,&comm);
    1333:   MPI_Comm_size(comm,&size);

    1335:   if (size > 1) { /* Sync nvar at shared vertices for all processes */
    1336:     PetscSF           sf = network->plex->sf;
    1337:     PetscInt          *local_nvar, *remote_nvar,nroots,nleaves,p=-1,i,nsv;
    1338:     const PetscInt    *svtx;
    1339:     PetscBool         ghost;
    1340:     /*
    1341:      PetscMPIInt       rank;
    1342:      const PetscInt    *ilocal;
    1343:      const PetscSFNode *iremote;
    1344:      MPI_Comm_rank(comm,&rank);
    1345:      PetscSFGetGraph(sf,&nroots,&nleaves,&ilocal,&iremote);
    1346:     */
    1347:     PetscSFGetGraph(sf,&nroots,&nleaves,NULL,NULL);
    1348:     PetscCalloc2(nroots,&local_nvar,nroots,&remote_nvar);

    1350:     /* Leaves copy user's nvar to local_nvar */
    1351:     DMNetworkGetSharedVertices(dm,&nsv,&svtx);
    1352:     for (i=0; i<nsv; i++) {
    1353:       p = svtx[i];
    1354:       DMNetworkIsGhostVertex(dm,p,&ghost);
    1355:       if (!ghost) continue;
    1356:       PetscSectionGetDof(network->DofSection,p,&local_nvar[p]);
    1357:       /* printf("[%d] Before SFReduce: leaf local_nvar[%d] = %d\n",rank,p,local_nvar[p]); */
    1358:     }

    1360:     /* Leaves add local_nvar to root remote_nvar */
    1361:     PetscSFReduceBegin(sf, MPIU_INT, local_nvar, remote_nvar, MPI_SUM);
    1362:     PetscSFReduceEnd(sf, MPIU_INT, local_nvar, remote_nvar, MPI_SUM);
    1363:     /* printf("[%d] remote_nvar[%d] = %d\n",rank,p,remote_nvar[p]); */

    1365:     /* Update roots' local_nvar */
    1366:     for (i=0; i<nsv; i++) {
    1367:       p = svtx[i];
    1368:       DMNetworkIsGhostVertex(dm,p,&ghost);
    1369:       if (ghost) continue;

    1371:       PetscSectionAddDof(network->DofSection,p,remote_nvar[p]);
    1372:       PetscSectionGetDof(network->DofSection,p,&local_nvar[p]);
    1373:       /* printf("[%d]  After SFReduce: root local_nvar[%d] = %d\n",rank,p,local_nvar[p]); */
    1374:     }

    1376:     /* Roots Bcast nvar to leaves */
    1377:     PetscSFBcastBegin(sf, MPIU_INT, local_nvar, remote_nvar,MPI_REPLACE);
    1378:     PetscSFBcastEnd(sf, MPIU_INT, local_nvar, remote_nvar,MPI_REPLACE);

    1380:     /* Leaves reset receved/remote nvar to dm */
    1381:     for (i=0; i<nsv; i++) {
    1382:       DMNetworkIsGhostVertex(dm,p,&ghost);
    1383:       if (!ghost) continue;
    1384:       p = svtx[i];
    1385:       /* printf("[%d] leaf reset nvar %d at p= %d \n",rank,remote_nvar[p],p); */
    1386:       /* DMNetworkSetNumVariables(dm,p,remote_nvar[p]) */
    1387:       PetscSectionSetDof(network->DofSection,p,remote_nvar[p]);
    1388:     }

    1390:     PetscFree2(local_nvar,remote_nvar);
    1391:   }

    1393:   PetscSectionSetUp(network->DofSection);
    1394:   return(0);
    1395: }

    1397: /* Get a subsection from a range of points */
    1398: static PetscErrorCode DMNetworkGetSubSection_private(PetscSection main,PetscInt pstart,PetscInt pend,PetscSection *subsection)
    1399: {
    1401:   PetscInt       i, nvar;

    1404:   PetscSectionCreate(PetscObjectComm((PetscObject)main), subsection);
    1405:   PetscSectionSetChart(*subsection, 0, pend - pstart);
    1406:   for (i = pstart; i < pend; i++) {
    1407:     PetscSectionGetDof(main,i,&nvar);
    1408:     PetscSectionSetDof(*subsection, i - pstart, nvar);
    1409:   }

    1411:   PetscSectionSetUp(*subsection);
    1412:   return(0);
    1413: }

    1415: /* Create a submap of points with a GlobalToLocal structure */
    1416: static PetscErrorCode DMNetworkSetSubMap_private(PetscInt pstart, PetscInt pend, ISLocalToGlobalMapping *map)
    1417: {
    1419:   PetscInt       i, *subpoints;

    1422:   /* Create index sets to map from "points" to "subpoints" */
    1423:   PetscMalloc1(pend - pstart, &subpoints);
    1424:   for (i = pstart; i < pend; i++) {
    1425:     subpoints[i - pstart] = i;
    1426:   }
    1427:   ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,pend-pstart,subpoints,PETSC_COPY_VALUES,map);
    1428:   PetscFree(subpoints);
    1429:   return(0);
    1430: }

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

    1435:   Collective on dm

    1437:   Input Parameters:
    1438: . dm - the DMNetworkObject

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

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

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

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

    1448:   Level: intermediate

    1450: @*/
    1451: PetscErrorCode DMNetworkAssembleGraphStructures(DM dm)
    1452: {
    1454:   MPI_Comm       comm;
    1455:   PetscMPIInt    rank, size;
    1456:   DM_Network     *network = (DM_Network*)dm->data;

    1459:   PetscObjectGetComm((PetscObject)dm,&comm);
    1460:   MPI_Comm_rank(comm, &rank);
    1461:   MPI_Comm_size(comm, &size);

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

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

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

    1474:     PetscSectionCreateGlobalSection(network->vertex.DofSection, network->vertex.sf, PETSC_FALSE, PETSC_FALSE, &network->vertex.GlobalDofSection);
    1475:     PetscSFGetSubSF(network->plex->sf, network->edge.mapping, &network->edge.sf);
    1476:     PetscSectionCreateGlobalSection(network->edge.DofSection, network->edge.sf, PETSC_FALSE, PETSC_FALSE, &network->edge.GlobalDofSection);
    1477:   } else {
    1478:     /* create structures for vertex */
    1479:     PetscSectionClone(network->vertex.DofSection,&network->vertex.GlobalDofSection);
    1480:     /* create structures for edge */
    1481:     PetscSectionClone(network->edge.DofSection,&network->edge.GlobalDofSection);
    1482:   }

    1484:   /* Add viewers */
    1485:   PetscObjectSetName((PetscObject)network->edge.GlobalDofSection,"Global edge dof section");
    1486:   PetscObjectSetName((PetscObject)network->vertex.GlobalDofSection,"Global vertex dof section");
    1487:   PetscSectionViewFromOptions(network->edge.GlobalDofSection, NULL, "-edge_global_section_view");
    1488:   PetscSectionViewFromOptions(network->vertex.GlobalDofSection, NULL, "-vertex_global_section_view");
    1489:   return(0);
    1490: }

    1492: /*@
    1493:   DMNetworkDistribute - Distributes the network and moves associated component data

    1495:   Collective

    1497:   Input Parameter:
    1498: + DM - the DMNetwork object
    1499: - overlap - the overlap of partitions, 0 is the default

    1501:   Notes:
    1502:   Distributes the network with <overlap>-overlapping partitioning of the edges.

    1504:   Level: intermediate

    1506: .seealso: DMNetworkCreate()
    1507: @*/
    1508: PetscErrorCode DMNetworkDistribute(DM *dm,PetscInt overlap)
    1509: {
    1510:   MPI_Comm       comm;
    1512:   PetscMPIInt    size;
    1513:   DM_Network     *oldDMnetwork = (DM_Network*)((*dm)->data);
    1514:   DM_Network     *newDMnetwork;
    1515:   PetscSF        pointsf=NULL;
    1516:   DM             newDM;
    1517:   PetscInt       j,e,v,offset,*subnetvtx,nsubnet,gidx,svtx_idx,nv;
    1518:   PetscInt       to_net,from_net,*svto;
    1519:   PetscPartitioner         part;
    1520:   DMNetworkComponentHeader header;

    1523:   PetscObjectGetComm((PetscObject)*dm,&comm);
    1524:   MPI_Comm_size(comm, &size);
    1525:   if (size == 1) return(0);

    1527:   /* This routine moves the component data to the appropriate processors. It makes use of the DataSection and the componentdataarray to move the component data to appropriate processors and returns a new DataSection and new componentdataarray. */
    1528:   DMNetworkCreate(PetscObjectComm((PetscObject)*dm),&newDM);
    1529:   newDMnetwork = (DM_Network*)newDM->data;
    1530:   newDMnetwork->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);

    1532:   /* Enable runtime options for petscpartitioner */
    1533:   DMPlexGetPartitioner(oldDMnetwork->plex,&part);
    1534:   PetscPartitionerSetFromOptions(part);

    1536:   /* Distribute plex dm */
    1537:   DMPlexDistribute(oldDMnetwork->plex,overlap,&pointsf,&newDMnetwork->plex);

    1539:   /* Distribute dof section */
    1540:   PetscSectionCreate(comm,&newDMnetwork->DofSection);
    1541:   PetscSFDistributeSection(pointsf,oldDMnetwork->DofSection,NULL,newDMnetwork->DofSection);

    1543:   /* Distribute data and associated section */
    1544:   PetscSectionCreate(comm,&newDMnetwork->DataSection);
    1545:   DMPlexDistributeData(newDMnetwork->plex,pointsf,oldDMnetwork->DataSection,MPIU_INT,(void*)oldDMnetwork->componentdataarray,newDMnetwork->DataSection,(void**)&newDMnetwork->componentdataarray);

    1547:   PetscSectionGetChart(newDMnetwork->DataSection,&newDMnetwork->pStart,&newDMnetwork->pEnd);
    1548:   DMPlexGetHeightStratum(newDMnetwork->plex,0, &newDMnetwork->eStart,&newDMnetwork->eEnd);
    1549:   DMPlexGetHeightStratum(newDMnetwork->plex,1,&newDMnetwork->vStart,&newDMnetwork->vEnd);
    1550:   newDMnetwork->nEdges    = newDMnetwork->eEnd - newDMnetwork->eStart;
    1551:   newDMnetwork->nVertices = newDMnetwork->vEnd - newDMnetwork->vStart;
    1552:   newDMnetwork->NVertices = oldDMnetwork->NVertices;
    1553:   newDMnetwork->NEdges    = oldDMnetwork->NEdges;
    1554:   newDMnetwork->svtable   = oldDMnetwork->svtable;
    1555:   oldDMnetwork->svtable   = NULL;

    1557:   /* Set Dof section as the section for dm */
    1558:   DMSetLocalSection(newDMnetwork->plex,newDMnetwork->DofSection);
    1559:   DMGetGlobalSection(newDMnetwork->plex,&newDMnetwork->GlobalDofSection);

    1561:   /* Set up subnetwork info in the newDM */
    1562:   newDMnetwork->Nsubnet = oldDMnetwork->Nsubnet;
    1563:   newDMnetwork->Nsvtx   = oldDMnetwork->Nsvtx;
    1564:   oldDMnetwork->Nsvtx   = 0;
    1565:   newDMnetwork->svtx    = oldDMnetwork->svtx; /* global vertices! */
    1566:   oldDMnetwork->svtx    = NULL;
    1567:   PetscCalloc1(newDMnetwork->Nsubnet,&newDMnetwork->subnet);

    1569:   /* Copy over the global number of vertices and edges in each subnetwork.
    1570:      Note: these are calculated in DMNetworkLayoutSetUp()
    1571:   */
    1572:   nsubnet = newDMnetwork->Nsubnet;
    1573:   for (j = 0; j < nsubnet; j++) {
    1574:     newDMnetwork->subnet[j].Nvtx  = oldDMnetwork->subnet[j].Nvtx;
    1575:     newDMnetwork->subnet[j].Nedge = oldDMnetwork->subnet[j].Nedge;
    1576:   }

    1578:   /* Get local nedges and nvtx for subnetworks */
    1579:   for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++) {
    1580:     PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);
    1581:     header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
    1582:     newDMnetwork->subnet[header->subnetid].nedge++;
    1583:   }

    1585:   for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++) {
    1586:     PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);
    1587:     header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);

    1589:     /* shared vertices: use gidx = header->index to check if v is a shared vertex */
    1590:     gidx = header->index;
    1591:     PetscTableFind(newDMnetwork->svtable,gidx+1,&svtx_idx);
    1592:     svtx_idx--;

    1594:     if (svtx_idx < 0) { /* not a shared vertex */
    1595:       newDMnetwork->subnet[header->subnetid].nvtx++;
    1596:     } else { /* a shared vertex belongs to more than one subnetworks, it is being counted by multiple subnets */
    1597:       from_net = newDMnetwork->svtx[svtx_idx].sv[0];
    1598:       newDMnetwork->subnet[from_net].nvtx++;
    1599:       for (j=1; j<newDMnetwork->svtx[svtx_idx].n; j++) {
    1600:         svto   = newDMnetwork->svtx[svtx_idx].sv + 2*j;
    1601:         to_net = svto[0];
    1602:         newDMnetwork->subnet[to_net].nvtx++;
    1603:       }
    1604:     }
    1605:   }

    1607:   /* Get total local nvtx for subnetworks */
    1608:   nv = 0;
    1609:   for (j=0; j<nsubnet; j++) nv += newDMnetwork->subnet[j].nvtx;
    1610:   nv += newDMnetwork->Nsvtx;

    1612:   /* Now create the vertices and edge arrays for the subnetworks */
    1613:   PetscCalloc1(nv,&subnetvtx);
    1614:   newDMnetwork->subnetvtx = subnetvtx;

    1616:   for (j=0; j<nsubnet; j++) {
    1617:     PetscCalloc1(newDMnetwork->subnet[j].nedge,&newDMnetwork->subnet[j].edges);
    1618:     newDMnetwork->subnet[j].vertices = subnetvtx;
    1619:     subnetvtx                       += newDMnetwork->subnet[j].nvtx;

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

    1626:   /* Set the edges and vertices in each subnetwork */
    1627:   for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++) {
    1628:     PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);
    1629:     header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
    1630:     newDMnetwork->subnet[header->subnetid].edges[newDMnetwork->subnet[header->subnetid].nedge++] = e;
    1631:   }

    1633:   nv = 0;
    1634:   for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++) {
    1635:     PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);
    1636:     header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);

    1638:     /* coupling vertices: use gidx = header->index to check if v is a coupling vertex */
    1639:     PetscTableFind(newDMnetwork->svtable,header->index+1,&svtx_idx);
    1640:     svtx_idx--;
    1641:     if (svtx_idx < 0) {
    1642:       newDMnetwork->subnet[header->subnetid].vertices[newDMnetwork->subnet[header->subnetid].nvtx++] = v;
    1643:     } else { /* a shared vertex */
    1644:       newDMnetwork->svertices[nv++] = v;

    1646:       from_net = newDMnetwork->svtx[svtx_idx].sv[0];
    1647:       newDMnetwork->subnet[from_net].vertices[newDMnetwork->subnet[from_net].nvtx++] = v;

    1649:       for (j=1; j<newDMnetwork->svtx[svtx_idx].n; j++) {
    1650:         svto   = newDMnetwork->svtx[svtx_idx].sv + 2*j;
    1651:         to_net = svto[0];
    1652:         newDMnetwork->subnet[to_net].vertices[newDMnetwork->subnet[to_net].nvtx++] = v;
    1653:       }
    1654:     }
    1655:   }
    1656:   newDMnetwork->nsvtx = nv;   /* num of local shared vertices */

    1658:   newDM->setupcalled = (*dm)->setupcalled;
    1659:   newDMnetwork->distributecalled = PETSC_TRUE;

    1661:   /* Free spaces */
    1662:   PetscSFDestroy(&pointsf);
    1663:   DMDestroy(dm);

    1665:   *dm  = newDM;
    1666:   return(0);
    1667: }

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

    1672:  Collective

    1674:   Input Parameters:
    1675: + mainSF - the original SF structure
    1676: - map - a ISLocalToGlobal mapping that contains the subset of points

    1678:   Output Parameters:
    1679: . subSF - a subset of the mainSF for the desired subset.

    1681:   Level: intermediate
    1682: @*/
    1683: PetscErrorCode PetscSFGetSubSF(PetscSF mainsf,ISLocalToGlobalMapping map,PetscSF *subSF)
    1684: {
    1685:   PetscErrorCode        ierr;
    1686:   PetscInt              nroots, nleaves, *ilocal_sub;
    1687:   PetscInt              i, *ilocal_map, nroots_sub, nleaves_sub = 0;
    1688:   PetscInt              *local_points, *remote_points;
    1689:   PetscSFNode           *iremote_sub;
    1690:   const PetscInt        *ilocal;
    1691:   const PetscSFNode     *iremote;

    1694:   PetscSFGetGraph(mainsf,&nroots,&nleaves,&ilocal,&iremote);

    1696:   /* Look for leaves that pertain to the subset of points. Get the local ordering */
    1697:   PetscMalloc1(nleaves,&ilocal_map);
    1698:   ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nleaves,ilocal,NULL,ilocal_map);
    1699:   for (i = 0; i < nleaves; i++) {
    1700:     if (ilocal_map[i] != -1) nleaves_sub += 1;
    1701:   }
    1702:   /* Re-number ilocal with subset numbering. Need information from roots */
    1703:   PetscMalloc2(nroots,&local_points,nroots,&remote_points);
    1704:   for (i = 0; i < nroots; i++) local_points[i] = i;
    1705:   ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nroots,local_points,NULL,local_points);
    1706:   PetscSFBcastBegin(mainsf, MPIU_INT, local_points, remote_points,MPI_REPLACE);
    1707:   PetscSFBcastEnd(mainsf, MPIU_INT, local_points, remote_points,MPI_REPLACE);
    1708:   /* Fill up graph using local (that is, local to the subset) numbering. */
    1709:   PetscMalloc1(nleaves_sub,&ilocal_sub);
    1710:   PetscMalloc1(nleaves_sub,&iremote_sub);
    1711:   nleaves_sub = 0;
    1712:   for (i = 0; i < nleaves; i++) {
    1713:     if (ilocal_map[i] != -1) {
    1714:       ilocal_sub[nleaves_sub] = ilocal_map[i];
    1715:       iremote_sub[nleaves_sub].rank = iremote[i].rank;
    1716:       iremote_sub[nleaves_sub].index = remote_points[ilocal[i]];
    1717:       nleaves_sub += 1;
    1718:     }
    1719:   }
    1720:   PetscFree2(local_points,remote_points);
    1721:   ISLocalToGlobalMappingGetSize(map,&nroots_sub);

    1723:   /* Create new subSF */
    1724:   PetscSFCreate(PETSC_COMM_WORLD,subSF);
    1725:   PetscSFSetFromOptions(*subSF);
    1726:   PetscSFSetGraph(*subSF,nroots_sub,nleaves_sub,ilocal_sub,PETSC_OWN_POINTER,iremote_sub,PETSC_COPY_VALUES);
    1727:   PetscFree(ilocal_map);
    1728:   PetscFree(iremote_sub);
    1729:   return(0);
    1730: }

    1732: /*@C
    1733:   DMNetworkGetSupportingEdges - Return the supporting edges for this vertex point

    1735:   Not Collective

    1737:   Input Parameters:
    1738: + dm - the DMNetwork object
    1739: - p  - the vertex point

    1741:   Output Parameters:
    1742: + nedges - number of edges connected to this vertex point
    1743: - edges  - list of edge points

    1745:   Level: beginner

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

    1751: .seealso: DMNetworkCreate(), DMNetworkGetConnectedVertices()
    1752: @*/
    1753: PetscErrorCode DMNetworkGetSupportingEdges(DM dm,PetscInt vertex,PetscInt *nedges,const PetscInt *edges[])
    1754: {
    1756:   DM_Network     *network = (DM_Network*)dm->data;

    1759:   DMPlexGetSupportSize(network->plex,vertex,nedges);
    1760:   if (edges) {DMPlexGetSupport(network->plex,vertex,edges);}
    1761:   return(0);
    1762: }

    1764: /*@C
    1765:   DMNetworkGetConnectedVertices - Return the connected vertices for this edge point

    1767:   Not Collective

    1769:   Input Parameters:
    1770: + dm - the DMNetwork object
    1771: - p - the edge point

    1773:   Output Parameters:
    1774: . vertices - vertices connected to this edge

    1776:   Level: beginner

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

    1782: .seealso: DMNetworkCreate(), DMNetworkGetSupportingEdges()
    1783: @*/
    1784: PetscErrorCode DMNetworkGetConnectedVertices(DM dm,PetscInt edge,const PetscInt *vertices[])
    1785: {
    1787:   DM_Network     *network = (DM_Network*)dm->data;

    1790:   DMPlexGetCone(network->plex,edge,vertices);
    1791:   return(0);
    1792: }

    1794: /*@
    1795:   DMNetworkIsSharedVertex - Returns TRUE if the vertex is shared by subnetworks

    1797:   Not Collective

    1799:   Input Parameters:
    1800: + dm - the DMNetwork object
    1801: - p - the vertex point

    1803:   Output Parameter:
    1804: . flag - TRUE if the vertex is shared by subnetworks

    1806:   Level: beginner

    1808: .seealso: DMNetworkAddSharedVertices(), DMNetworkIsGhostVertex()
    1809: @*/
    1810: PetscErrorCode DMNetworkIsSharedVertex(DM dm,PetscInt p,PetscBool *flag)
    1811: {
    1813:   PetscInt       i;

    1816:   *flag = PETSC_FALSE;

    1818:   if (dm->setupcalled) { /* DMNetworkGetGlobalVertexIndex() requires DMSetUp() be called */
    1819:     DM_Network     *network = (DM_Network*)dm->data;
    1820:     PetscInt       gidx;
    1821:     DMNetworkGetGlobalVertexIndex(dm,p,&gidx);
    1822:     PetscTableFind(network->svtable,gidx+1,&i);
    1823:     if (i) *flag = PETSC_TRUE;
    1824:   } else { /* would be removed? */
    1825:     PetscInt       nv;
    1826:     const PetscInt *vtx;
    1827:     DMNetworkGetSharedVertices(dm,&nv,&vtx);
    1828:     for (i=0; i<nv; i++) {
    1829:       if (p == vtx[i]) {
    1830:         *flag = PETSC_TRUE;
    1831:         break;
    1832:       }
    1833:     }
    1834:   }
    1835:   return(0);
    1836: }

    1838: /*@
    1839:   DMNetworkIsGhostVertex - Returns TRUE if the vertex is a ghost vertex

    1841:   Not Collective

    1843:   Input Parameters:
    1844: + dm - the DMNetwork object
    1845: - p - the vertex point

    1847:   Output Parameter:
    1848: . isghost - TRUE if the vertex is a ghost point

    1850:   Level: beginner

    1852: .seealso: DMNetworkGetConnectedVertices(), DMNetworkGetVertexRange(), DMNetworkIsSharedVertex()
    1853: @*/
    1854: PetscErrorCode DMNetworkIsGhostVertex(DM dm,PetscInt p,PetscBool *isghost)
    1855: {
    1857:   DM_Network     *network = (DM_Network*)dm->data;
    1858:   PetscInt       offsetg;
    1859:   PetscSection   sectiong;

    1862:   *isghost = PETSC_FALSE;
    1863:   DMGetGlobalSection(network->plex,&sectiong);
    1864:   PetscSectionGetOffset(sectiong,p,&offsetg);
    1865:   if (offsetg < 0) *isghost = PETSC_TRUE;
    1866:   return(0);
    1867: }

    1869: PetscErrorCode DMSetUp_Network(DM dm)
    1870: {
    1872:   DM_Network     *network=(DM_Network*)dm->data;

    1875:   DMNetworkComponentSetUp(dm);
    1876:   DMNetworkVariablesSetUp(dm);

    1878:   DMSetLocalSection(network->plex,network->DofSection);
    1879:   DMGetGlobalSection(network->plex,&network->GlobalDofSection);

    1881:   dm->setupcalled = PETSC_TRUE;
    1882:   DMViewFromOptions(dm,NULL,"-dm_view");
    1883:   return(0);
    1884: }

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

    1890:   Collective

    1892:   Input Parameters:
    1893: + dm - the DMNetwork object
    1894: . eflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for edges
    1895: - vflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for vertices

    1897:  Level: intermediate

    1899: @*/
    1900: PetscErrorCode DMNetworkHasJacobian(DM dm,PetscBool eflg,PetscBool vflg)
    1901: {
    1902:   DM_Network     *network=(DM_Network*)dm->data;
    1904:   PetscInt       nVertices = network->nVertices;

    1907:   network->userEdgeJacobian   = eflg;
    1908:   network->userVertexJacobian = vflg;

    1910:   if (eflg && !network->Je) {
    1911:     PetscCalloc1(3*network->nEdges,&network->Je);
    1912:   }

    1914:   if (vflg && !network->Jv && nVertices) {
    1915:     PetscInt       i,*vptr,nedges,vStart=network->vStart;
    1916:     PetscInt       nedges_total;
    1917:     const PetscInt *edges;

    1919:     /* count nvertex_total */
    1920:     nedges_total = 0;
    1921:     PetscMalloc1(nVertices+1,&vptr);

    1923:     vptr[0] = 0;
    1924:     for (i=0; i<nVertices; i++) {
    1925:       DMNetworkGetSupportingEdges(dm,i+vStart,&nedges,&edges);
    1926:       nedges_total += nedges;
    1927:       vptr[i+1] = vptr[i] + 2*nedges + 1;
    1928:     }

    1930:     PetscCalloc1(2*nedges_total+nVertices,&network->Jv);
    1931:     network->Jvptr = vptr;
    1932:   }
    1933:   return(0);
    1934: }

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

    1939:   Not Collective

    1941:   Input Parameters:
    1942: + dm - the DMNetwork object
    1943: . p - the edge point
    1944: - J - array (size = 3) of Jacobian submatrices for this edge point:
    1945:         J[0]: this edge
    1946:         J[1] and J[2]: connected vertices, obtained by calling DMNetworkGetConnectedVertices()

    1948:   Level: advanced

    1950: .seealso: DMNetworkVertexSetMatrix()
    1951: @*/
    1952: PetscErrorCode DMNetworkEdgeSetMatrix(DM dm,PetscInt p,Mat J[])
    1953: {
    1954:   DM_Network *network=(DM_Network*)dm->data;

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

    1959:   if (J) {
    1960:     network->Je[3*p]   = J[0];
    1961:     network->Je[3*p+1] = J[1];
    1962:     network->Je[3*p+2] = J[2];
    1963:   }
    1964:   return(0);
    1965: }

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

    1970:   Not Collective

    1972:   Input Parameters:
    1973: + dm - The DMNetwork object
    1974: . p - the vertex point
    1975: - J - array of Jacobian (size = 2*(num of supporting edges) + 1) submatrices for this vertex point:
    1976:         J[0]:       this vertex
    1977:         J[1+2*i]:   i-th supporting edge
    1978:         J[1+2*i+1]: i-th connected vertex

    1980:   Level: advanced

    1982: .seealso: DMNetworkEdgeSetMatrix()
    1983: @*/
    1984: PetscErrorCode DMNetworkVertexSetMatrix(DM dm,PetscInt p,Mat J[])
    1985: {
    1987:   DM_Network     *network=(DM_Network*)dm->data;
    1988:   PetscInt       i,*vptr,nedges,vStart=network->vStart;
    1989:   const PetscInt *edges;

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

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

    1998:     /* Set Jacobian for each supporting edge and connected vertex */
    1999:     DMNetworkGetSupportingEdges(dm,p,&nedges,&edges);
    2000:     for (i=1; i<=2*nedges; i++) network->Jv[vptr[p-vStart]+i] = J[i];
    2001:   }
    2002:   return(0);
    2003: }

    2005: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
    2006: {
    2008:   PetscInt       j;
    2009:   PetscScalar    val=(PetscScalar)ncols;

    2012:   if (!ghost) {
    2013:     for (j=0; j<nrows; j++) {
    2014:       VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
    2015:     }
    2016:   } else {
    2017:     for (j=0; j<nrows; j++) {
    2018:       VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
    2019:     }
    2020:   }
    2021:   return(0);
    2022: }

    2024: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
    2025: {
    2027:   PetscInt       j,ncols_u;
    2028:   PetscScalar    val;

    2031:   if (!ghost) {
    2032:     for (j=0; j<nrows; j++) {
    2033:       MatGetRow(Ju,j,&ncols_u,NULL,NULL);
    2034:       val = (PetscScalar)ncols_u;
    2035:       VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
    2036:       MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
    2037:     }
    2038:   } else {
    2039:     for (j=0; j<nrows; j++) {
    2040:       MatGetRow(Ju,j,&ncols_u,NULL,NULL);
    2041:       val = (PetscScalar)ncols_u;
    2042:       VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
    2043:       MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
    2044:     }
    2045:   }
    2046:   return(0);
    2047: }

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

    2054:   if (Ju) {
    2055:     MatSetPreallocationUserblock_private(Ju,nrows,rows,ncols,ghost,vdnz,vonz);
    2056:   } else {
    2057:     MatSetPreallocationDenseblock_private(nrows,rows,ncols,ghost,vdnz,vonz);
    2058:   }
    2059:   return(0);
    2060: }

    2062: PETSC_STATIC_INLINE PetscErrorCode MatSetDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
    2063: {
    2065:   PetscInt       j,*cols;
    2066:   PetscScalar    *zeros;

    2069:   PetscCalloc2(ncols,&cols,nrows*ncols,&zeros);
    2070:   for (j=0; j<ncols; j++) cols[j] = j+ cstart;
    2071:   MatSetValues(*J,nrows,rows,ncols,cols,zeros,INSERT_VALUES);
    2072:   PetscFree2(cols,zeros);
    2073:   return(0);
    2074: }

    2076: PETSC_STATIC_INLINE PetscErrorCode MatSetUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
    2077: {
    2079:   PetscInt       j,M,N,row,col,ncols_u;
    2080:   const PetscInt *cols;
    2081:   PetscScalar    zero=0.0;

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

    2087:   for (row=0; row<nrows; row++) {
    2088:     MatGetRow(Ju,row,&ncols_u,&cols,NULL);
    2089:     for (j=0; j<ncols_u; j++) {
    2090:       col = cols[j] + cstart;
    2091:       MatSetValues(*J,1,&rows[row],1,&col,&zero,INSERT_VALUES);
    2092:     }
    2093:     MatRestoreRow(Ju,row,&ncols_u,&cols,NULL);
    2094:   }
    2095:   return(0);
    2096: }

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

    2103:   if (Ju) {
    2104:     MatSetUserblock_private(Ju,nrows,rows,ncols,cstart,J);
    2105:   } else {
    2106:     MatSetDenseblock_private(nrows,rows,ncols,cstart,J);
    2107:   }
    2108:   return(0);
    2109: }

    2111: /* 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.
    2112: */
    2113: PetscErrorCode CreateSubGlobalToLocalMapping_private(PetscSection globalsec, PetscSection localsec, ISLocalToGlobalMapping *ltog)
    2114: {
    2116:   PetscInt       i,size,dof;
    2117:   PetscInt       *glob2loc;

    2120:   PetscSectionGetStorageSize(localsec,&size);
    2121:   PetscMalloc1(size,&glob2loc);

    2123:   for (i = 0; i < size; i++) {
    2124:     PetscSectionGetOffset(globalsec,i,&dof);
    2125:     dof = (dof >= 0) ? dof : -(dof + 1);
    2126:     glob2loc[i] = dof;
    2127:   }

    2129:   ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,size,glob2loc,PETSC_OWN_POINTER,ltog);
    2130: #if 0
    2131:   PetscIntView(size,glob2loc,PETSC_VIEWER_STDOUT_WORLD);
    2132: #endif
    2133:   return(0);
    2134: }

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

    2138: PetscErrorCode DMCreateMatrix_Network_Nest(DM dm,Mat *J)
    2139: {
    2141:   DM_Network     *network = (DM_Network*)dm->data;
    2142:   PetscMPIInt    rank, size;
    2143:   PetscInt       eDof,vDof;
    2144:   Mat            j11,j12,j21,j22,bA[2][2];
    2145:   MPI_Comm       comm;
    2146:   ISLocalToGlobalMapping eISMap,vISMap;

    2149:   PetscObjectGetComm((PetscObject)dm,&comm);
    2150:   MPI_Comm_rank(comm,&rank);
    2151:   MPI_Comm_size(comm,&size);

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

    2156:   MatCreate(comm, &j11);
    2157:   MatSetSizes(j11, eDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
    2158:   MatSetType(j11, MATMPIAIJ);

    2160:   MatCreate(comm, &j12);
    2161:   MatSetSizes(j12, eDof, vDof, PETSC_DETERMINE ,PETSC_DETERMINE);
    2162:   MatSetType(j12, MATMPIAIJ);

    2164:   MatCreate(comm, &j21);
    2165:   MatSetSizes(j21, vDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
    2166:   MatSetType(j21, MATMPIAIJ);

    2168:   MatCreate(comm, &j22);
    2169:   MatSetSizes(j22, vDof, vDof, PETSC_DETERMINE, PETSC_DETERMINE);
    2170:   MatSetType(j22, MATMPIAIJ);

    2172:   bA[0][0] = j11;
    2173:   bA[0][1] = j12;
    2174:   bA[1][0] = j21;
    2175:   bA[1][1] = j22;

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

    2180:   MatSetLocalToGlobalMapping(j11,eISMap,eISMap);
    2181:   MatSetLocalToGlobalMapping(j12,eISMap,vISMap);
    2182:   MatSetLocalToGlobalMapping(j21,vISMap,eISMap);
    2183:   MatSetLocalToGlobalMapping(j22,vISMap,vISMap);

    2185:   MatSetUp(j11);
    2186:   MatSetUp(j12);
    2187:   MatSetUp(j21);
    2188:   MatSetUp(j22);

    2190:   MatCreateNest(comm,2,NULL,2,NULL,&bA[0][0],J);
    2191:   MatSetUp(*J);
    2192:   MatNestSetVecType(*J,VECNEST);
    2193:   MatDestroy(&j11);
    2194:   MatDestroy(&j12);
    2195:   MatDestroy(&j21);
    2196:   MatDestroy(&j22);

    2198:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
    2199:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
    2200:   MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);

    2202:   /* Free structures */
    2203:   ISLocalToGlobalMappingDestroy(&eISMap);
    2204:   ISLocalToGlobalMappingDestroy(&vISMap);
    2205:   return(0);
    2206: }

    2208: PetscErrorCode DMCreateMatrix_Network(DM dm,Mat *J)
    2209: {
    2211:   DM_Network     *network = (DM_Network*)dm->data;
    2212:   PetscInt       eStart,eEnd,vStart,vEnd,rstart,nrows,*rows,localSize;
    2213:   PetscInt       cstart,ncols,j,e,v;
    2214:   PetscBool      ghost,ghost_vc,ghost2,isNest;
    2215:   Mat            Juser;
    2216:   PetscSection   sectionGlobal;
    2217:   PetscInt       nedges,*vptr=NULL,vc,*rows_v; /* suppress maybe-uninitialized warning */
    2218:   const PetscInt *edges,*cone;
    2219:   MPI_Comm       comm;
    2220:   MatType        mtype;
    2221:   Vec            vd_nz,vo_nz;
    2222:   PetscInt       *dnnz,*onnz;
    2223:   PetscScalar    *vdnz,*vonz;

    2226:   mtype = dm->mattype;
    2227:   PetscStrcmp(mtype,MATNEST,&isNest);
    2228:   if (isNest) {
    2229:     DMCreateMatrix_Network_Nest(dm,J);
    2230:     MatSetDM(*J,dm);
    2231:     return(0);
    2232:   }

    2234:   if (!network->userEdgeJacobian && !network->userVertexJacobian) {
    2235:     /* user does not provide Jacobian blocks */
    2236:     DMCreateMatrix_Plex(network->plex,J);
    2237:     MatSetDM(*J,dm);
    2238:     return(0);
    2239:   }

    2241:   MatCreate(PetscObjectComm((PetscObject)dm),J);
    2242:   DMGetGlobalSection(network->plex,&sectionGlobal);
    2243:   PetscSectionGetConstrainedStorageSize(sectionGlobal,&localSize);
    2244:   MatSetSizes(*J,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);

    2246:   MatSetType(*J,MATAIJ);
    2247:   MatSetFromOptions(*J);

    2249:   /* (1) Set matrix preallocation */
    2250:   /*------------------------------*/
    2251:   PetscObjectGetComm((PetscObject)dm,&comm);
    2252:   VecCreate(comm,&vd_nz);
    2253:   VecSetSizes(vd_nz,localSize,PETSC_DECIDE);
    2254:   VecSetFromOptions(vd_nz);
    2255:   VecSet(vd_nz,0.0);
    2256:   VecDuplicate(vd_nz,&vo_nz);

    2258:   /* Set preallocation for edges */
    2259:   /*-----------------------------*/
    2260:   DMNetworkGetEdgeRange(dm,&eStart,&eEnd);

    2262:   PetscMalloc1(localSize,&rows);
    2263:   for (e=eStart; e<eEnd; e++) {
    2264:     /* Get row indices */
    2265:     DMNetworkGetGlobalVecOffset(dm,e,ALL_COMPONENTS,&rstart);
    2266:     PetscSectionGetDof(network->DofSection,e,&nrows);
    2267:     if (nrows) {
    2268:       for (j=0; j<nrows; j++) rows[j] = j + rstart;

    2270:       /* Set preallocation for conntected vertices */
    2271:       DMNetworkGetConnectedVertices(dm,e,&cone);
    2272:       for (v=0; v<2; v++) {
    2273:         PetscSectionGetDof(network->DofSection,cone[v],&ncols);

    2275:         if (network->Je) {
    2276:           Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
    2277:         } else Juser = NULL;
    2278:         DMNetworkIsGhostVertex(dm,cone[v],&ghost);
    2279:         MatSetPreallocationblock_private(Juser,nrows,rows,ncols,ghost,vd_nz,vo_nz);
    2280:       }

    2282:       /* Set preallocation for edge self */
    2283:       cstart = rstart;
    2284:       if (network->Je) {
    2285:         Juser = network->Je[3*e]; /* Jacobian(e,e) */
    2286:       } else Juser = NULL;
    2287:       MatSetPreallocationblock_private(Juser,nrows,rows,nrows,PETSC_FALSE,vd_nz,vo_nz);
    2288:     }
    2289:   }

    2291:   /* Set preallocation for vertices */
    2292:   /*--------------------------------*/
    2293:   DMNetworkGetVertexRange(dm,&vStart,&vEnd);
    2294:   if (vEnd - vStart) vptr = network->Jvptr;

    2296:   for (v=vStart; v<vEnd; v++) {
    2297:     /* Get row indices */
    2298:     DMNetworkGetGlobalVecOffset(dm,v,ALL_COMPONENTS,&rstart);
    2299:     PetscSectionGetDof(network->DofSection,v,&nrows);
    2300:     if (!nrows) continue;

    2302:     DMNetworkIsGhostVertex(dm,v,&ghost);
    2303:     if (ghost) {
    2304:       PetscMalloc1(nrows,&rows_v);
    2305:     } else {
    2306:       rows_v = rows;
    2307:     }

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

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

    2314:     for (e=0; e<nedges; e++) {
    2315:       /* Supporting edges */
    2316:       DMNetworkGetGlobalVecOffset(dm,edges[e],ALL_COMPONENTS,&cstart);
    2317:       PetscSectionGetDof(network->DofSection,edges[e],&ncols);

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

    2324:       /* Connected vertices */
    2325:       DMNetworkGetConnectedVertices(dm,edges[e],&cone);
    2326:       vc = (v == cone[0]) ? cone[1]:cone[0];
    2327:       DMNetworkIsGhostVertex(dm,vc,&ghost_vc);

    2329:       PetscSectionGetDof(network->DofSection,vc,&ncols);

    2331:       if (network->Jv) {
    2332:         Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
    2333:       } else Juser = NULL;
    2334:       if (ghost_vc||ghost) {
    2335:         ghost2 = PETSC_TRUE;
    2336:       } else {
    2337:         ghost2 = PETSC_FALSE;
    2338:       }
    2339:       MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost2,vd_nz,vo_nz);
    2340:     }

    2342:     /* Set preallocation for vertex self */
    2343:     DMNetworkIsGhostVertex(dm,v,&ghost);
    2344:     if (!ghost) {
    2345:       DMNetworkGetGlobalVecOffset(dm,v,ALL_COMPONENTS,&cstart);
    2346:       if (network->Jv) {
    2347:         Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
    2348:       } else Juser = NULL;
    2349:       MatSetPreallocationblock_private(Juser,nrows,rows_v,nrows,PETSC_FALSE,vd_nz,vo_nz);
    2350:     }
    2351:     if (ghost) {
    2352:       PetscFree(rows_v);
    2353:     }
    2354:   }

    2356:   VecAssemblyBegin(vd_nz);
    2357:   VecAssemblyBegin(vo_nz);

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

    2361:   VecAssemblyEnd(vd_nz);
    2362:   VecAssemblyEnd(vo_nz);

    2364:   VecGetArray(vd_nz,&vdnz);
    2365:   VecGetArray(vo_nz,&vonz);
    2366:   for (j=0; j<localSize; j++) {
    2367:     dnnz[j] = (PetscInt)PetscRealPart(vdnz[j]);
    2368:     onnz[j] = (PetscInt)PetscRealPart(vonz[j]);
    2369:   }
    2370:   VecRestoreArray(vd_nz,&vdnz);
    2371:   VecRestoreArray(vo_nz,&vonz);
    2372:   VecDestroy(&vd_nz);
    2373:   VecDestroy(&vo_nz);

    2375:   MatSeqAIJSetPreallocation(*J,0,dnnz);
    2376:   MatMPIAIJSetPreallocation(*J,0,dnnz,0,onnz);
    2377:   MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);

    2379:   PetscFree2(dnnz,onnz);

    2381:   /* (2) Set matrix entries for edges */
    2382:   /*----------------------------------*/
    2383:   for (e=eStart; e<eEnd; e++) {
    2384:     /* Get row indices */
    2385:     DMNetworkGetGlobalVecOffset(dm,e,ALL_COMPONENTS,&rstart);
    2386:     PetscSectionGetDof(network->DofSection,e,&nrows);
    2387:     if (nrows) {
    2388:       for (j=0; j<nrows; j++) rows[j] = j + rstart;

    2390:       /* Set matrix entries for conntected vertices */
    2391:       DMNetworkGetConnectedVertices(dm,e,&cone);
    2392:       for (v=0; v<2; v++) {
    2393:         DMNetworkGetGlobalVecOffset(dm,cone[v],ALL_COMPONENTS,&cstart);
    2394:         PetscSectionGetDof(network->DofSection,cone[v],&ncols);

    2396:         if (network->Je) {
    2397:           Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
    2398:         } else Juser = NULL;
    2399:         MatSetblock_private(Juser,nrows,rows,ncols,cstart,J);
    2400:       }

    2402:       /* Set matrix entries for edge self */
    2403:       cstart = rstart;
    2404:       if (network->Je) {
    2405:         Juser = network->Je[3*e]; /* Jacobian(e,e) */
    2406:       } else Juser = NULL;
    2407:       MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);
    2408:     }
    2409:   }

    2411:   /* Set matrix entries for vertices */
    2412:   /*---------------------------------*/
    2413:   for (v=vStart; v<vEnd; v++) {
    2414:     /* Get row indices */
    2415:     DMNetworkGetGlobalVecOffset(dm,v,ALL_COMPONENTS,&rstart);
    2416:     PetscSectionGetDof(network->DofSection,v,&nrows);
    2417:     if (!nrows) continue;

    2419:     DMNetworkIsGhostVertex(dm,v,&ghost);
    2420:     if (ghost) {
    2421:       PetscMalloc1(nrows,&rows_v);
    2422:     } else {
    2423:       rows_v = rows;
    2424:     }
    2425:     for (j=0; j<nrows; j++) rows_v[j] = j + rstart;

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

    2430:     for (e=0; e<nedges; e++) {
    2431:       /* Supporting edges */
    2432:       DMNetworkGetGlobalVecOffset(dm,edges[e],ALL_COMPONENTS,&cstart);
    2433:       PetscSectionGetDof(network->DofSection,edges[e],&ncols);

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

    2440:       /* Connected vertices */
    2441:       DMNetworkGetConnectedVertices(dm,edges[e],&cone);
    2442:       vc = (v == cone[0]) ? cone[1]:cone[0];

    2444:       DMNetworkGetGlobalVecOffset(dm,vc,ALL_COMPONENTS,&cstart);
    2445:       PetscSectionGetDof(network->DofSection,vc,&ncols);

    2447:       if (network->Jv) {
    2448:         Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
    2449:       } else Juser = NULL;
    2450:       MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);
    2451:     }

    2453:     /* Set matrix entries for vertex self */
    2454:     if (!ghost) {
    2455:       DMNetworkGetGlobalVecOffset(dm,v,ALL_COMPONENTS,&cstart);
    2456:       if (network->Jv) {
    2457:         Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
    2458:       } else Juser = NULL;
    2459:       MatSetblock_private(Juser,nrows,rows_v,nrows,cstart,J);
    2460:     }
    2461:     if (ghost) {
    2462:       PetscFree(rows_v);
    2463:     }
    2464:   }
    2465:   PetscFree(rows);

    2467:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
    2468:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);

    2470:   MatSetDM(*J,dm);
    2471:   return(0);
    2472: }

    2474: PetscErrorCode DMDestroy_Network(DM dm)
    2475: {
    2477:   DM_Network     *network = (DM_Network*)dm->data;
    2478:   PetscInt       j;

    2481:   if (--network->refct > 0) return(0);
    2482:   if (network->Je) {
    2483:     PetscFree(network->Je);
    2484:   }
    2485:   if (network->Jv) {
    2486:     PetscFree(network->Jvptr);
    2487:     PetscFree(network->Jv);
    2488:   }

    2490:   ISLocalToGlobalMappingDestroy(&network->vertex.mapping);
    2491:   PetscSectionDestroy(&network->vertex.DofSection);
    2492:   PetscSectionDestroy(&network->vertex.GlobalDofSection);
    2493:   if (network->vltog) {
    2494:     PetscFree(network->vltog);
    2495:   }
    2496:   if (network->vertex.sf) {
    2497:     PetscSFDestroy(&network->vertex.sf);
    2498:   }
    2499:   /* edge */
    2500:   ISLocalToGlobalMappingDestroy(&network->edge.mapping);
    2501:   PetscSectionDestroy(&network->edge.DofSection);
    2502:   PetscSectionDestroy(&network->edge.GlobalDofSection);
    2503:   if (network->edge.sf) {
    2504:     PetscSFDestroy(&network->edge.sf);
    2505:   }
    2506:   DMDestroy(&network->plex);
    2507:   PetscSectionDestroy(&network->DataSection);
    2508:   PetscSectionDestroy(&network->DofSection);

    2510:   for (j=0; j<network->Nsvtx; j++) {
    2511:     PetscFree(network->svtx[j].sv);
    2512:   }
    2513:   if (network->svtx) {PetscFree(network->svtx);}

    2515:   for (j=0; j<network->Nsubnet; j++) {
    2516:     PetscFree(network->subnet[j].edges);
    2517:   }
    2518:   if (network->subnetvtx) {PetscFree(network->subnetvtx);}

    2520:   PetscTableDestroy(&network->svtable);
    2521:   PetscFree(network->subnet);
    2522:   PetscFree(network->componentdataarray);
    2523:   PetscFree2(network->header,network->cvalue);
    2524:   PetscFree(network);
    2525:   return(0);
    2526: }

    2528: PetscErrorCode DMView_Network(DM dm,PetscViewer viewer)
    2529: {
    2531:   PetscBool      iascii;
    2532:   PetscMPIInt    rank;

    2535:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
    2536:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
    2539:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
    2540:   if (iascii) {
    2541:     const PetscInt *cone,*vtx,*edges;
    2542:     PetscInt       vfrom,vto,i,j,nv,ne,ncv,p,nsubnet;
    2543:     DM_Network     *network = (DM_Network*)dm->data;

    2545:     nsubnet = network->Nsubnet; /* num of subnetworks */
    2546:     if (!rank) {
    2547:       PetscPrintf(PETSC_COMM_SELF,"  NSubnets: %D; NEdges: %D; NVertices: %D; NSharedVertices: %D.\n",nsubnet,network->NEdges,network->NVertices,network->Nsvtx);
    2548:     }

    2550:     DMNetworkGetSharedVertices(dm,&ncv,&vtx);
    2551:     PetscViewerASCIIPushSynchronized(viewer);
    2552:     PetscViewerASCIISynchronizedPrintf(viewer, "  [%d] nEdges: %D; nVertices: %D; nSharedVertices: %D\n",rank,network->nEdges,network->nVertices,ncv);

    2554:     for (i=0; i<nsubnet; i++) {
    2555:       DMNetworkGetSubnetwork(dm,i,&nv,&ne,&vtx,&edges);
    2556:       if (ne) {
    2557:         PetscViewerASCIISynchronizedPrintf(viewer, "     Subnet %D: nEdges %D, nVertices(include shared vertices) %D\n",i,ne,nv);
    2558:         for (j=0; j<ne; j++) {
    2559:           p = edges[j];
    2560:           DMNetworkGetConnectedVertices(dm,p,&cone);
    2561:           DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);
    2562:           DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);
    2563:           DMNetworkGetGlobalEdgeIndex(dm,edges[j],&p);
    2564:           PetscViewerASCIISynchronizedPrintf(viewer, "       edge %D: %D ----> %D\n",p,vfrom,vto);
    2565:         }
    2566:       }
    2567:     }

    2569:     /* Shared vertices */
    2570:     DMNetworkGetSharedVertices(dm,&ncv,&vtx);
    2571:     if (ncv) {
    2572:       SVtx       *svtx = network->svtx;
    2573:       PetscInt    gidx,svtx_idx,nvto,vfrom_net,vfrom_idx,*svto;
    2574:       PetscBool   ghost;
    2575:       PetscViewerASCIISynchronizedPrintf(viewer, "     SharedVertices:\n");
    2576:       for (i=0; i<ncv; i++) {
    2577:         DMNetworkIsGhostVertex(dm,vtx[i],&ghost);
    2578:         if (ghost) continue;

    2580:         DMNetworkGetGlobalVertexIndex(dm,vtx[i],&gidx);
    2581:         PetscTableFind(network->svtable,gidx+1,&svtx_idx);
    2582:         svtx_idx--;
    2583:         nvto = svtx[svtx_idx].n;

    2585:         vfrom_net = svtx[svtx_idx].sv[0];
    2586:         vfrom_idx = svtx[svtx_idx].sv[1];
    2587:         PetscViewerASCIISynchronizedPrintf(viewer, "       svtx %D: global index %D, subnet[%D].%D ---->\n",i,gidx,vfrom_net,vfrom_idx);
    2588:         for (j=1; j<nvto; j++) {
    2589:           svto = svtx[svtx_idx].sv + 2*j;
    2590:           PetscViewerASCIISynchronizedPrintf(viewer, "                                           ----> subnet[%D].%D\n",svto[0],svto[1]);
    2591:         }
    2592:       }
    2593:     }
    2594:     PetscViewerFlush(viewer);
    2595:     PetscViewerASCIIPopSynchronized(viewer);
    2596:   } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMNetwork writing", ((PetscObject)viewer)->type_name);
    2597:   return(0);
    2598: }

    2600: PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l)
    2601: {
    2603:   DM_Network     *network = (DM_Network*)dm->data;

    2606:   DMGlobalToLocalBegin(network->plex,g,mode,l);
    2607:   return(0);
    2608: }

    2610: PetscErrorCode DMGlobalToLocalEnd_Network(DM dm, Vec g, InsertMode mode, Vec l)
    2611: {
    2613:   DM_Network     *network = (DM_Network*)dm->data;

    2616:   DMGlobalToLocalEnd(network->plex,g,mode,l);
    2617:   return(0);
    2618: }

    2620: PetscErrorCode DMLocalToGlobalBegin_Network(DM dm, Vec l, InsertMode mode, Vec g)
    2621: {
    2623:   DM_Network     *network = (DM_Network*)dm->data;

    2626:   DMLocalToGlobalBegin(network->plex,l,mode,g);
    2627:   return(0);
    2628: }

    2630: PetscErrorCode DMLocalToGlobalEnd_Network(DM dm, Vec l, InsertMode mode, Vec g)
    2631: {
    2633:   DM_Network     *network = (DM_Network*)dm->data;

    2636:   DMLocalToGlobalEnd(network->plex,l,mode,g);
    2637:   return(0);
    2638: }

    2640: /*@
    2641:   DMNetworkGetVertexLocalToGlobalOrdering - Get vertex global index

    2643:   Not collective

    2645:   Input Parameters:
    2646: + dm - the dm object
    2647: - vloc - local vertex ordering, start from 0

    2649:   Output Parameters:
    2650: .  vg  - global vertex ordering, start from 0

    2652:   Level: advanced

    2654: .seealso: DMNetworkSetVertexLocalToGlobalOrdering()
    2655: @*/
    2656: PetscErrorCode DMNetworkGetVertexLocalToGlobalOrdering(DM dm,PetscInt vloc,PetscInt *vg)
    2657: {
    2658:   DM_Network  *network = (DM_Network*)dm->data;
    2659:   PetscInt    *vltog = network->vltog;

    2662:   if (!vltog) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkSetVertexLocalToGlobalOrdering() first");
    2663:   *vg = vltog[vloc];
    2664:   return(0);
    2665: }

    2667: /*@
    2668:   DMNetworkSetVertexLocalToGlobalOrdering - Create and setup vertex local to global map

    2670:   Collective

    2672:   Input Parameters:
    2673: . dm - the dm object

    2675:   Level: advanced

    2677: .seealso: DMNetworkGetGlobalVertexIndex()
    2678: @*/
    2679: PetscErrorCode DMNetworkSetVertexLocalToGlobalOrdering(DM dm)
    2680: {
    2681:   PetscErrorCode    ierr;
    2682:   DM_Network        *network=(DM_Network*)dm->data;
    2683:   MPI_Comm          comm;
    2684:   PetscMPIInt       rank,size,*displs=NULL,*recvcounts=NULL,remoterank;
    2685:   PetscBool         ghost;
    2686:   PetscInt          *vltog,nroots,nleaves,i,*vrange,k,N,lidx;
    2687:   const PetscSFNode *iremote;
    2688:   PetscSF           vsf;
    2689:   Vec               Vleaves,Vleaves_seq;
    2690:   VecScatter        ctx;
    2691:   PetscScalar       *varr,val;
    2692:   const PetscScalar *varr_read;

    2695:   PetscObjectGetComm((PetscObject)dm,&comm);
    2696:   MPI_Comm_size(comm,&size);
    2697:   MPI_Comm_rank(comm,&rank);

    2699:   if (size == 1) {
    2700:     nroots = network->vEnd - network->vStart;
    2701:     PetscMalloc1(nroots, &vltog);
    2702:     for (i=0; i<nroots; i++) vltog[i] = i;
    2703:     network->vltog = vltog;
    2704:     return(0);
    2705:   }

    2707:   if (!network->distributecalled) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkDistribute() first");
    2708:   if (network->vltog) {
    2709:     PetscFree(network->vltog);
    2710:   }

    2712:   DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);
    2713:   PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);
    2714:   vsf = network->vertex.sf;

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

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

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

    2726:   PetscMalloc1(nroots, &vltog);
    2727:   network->vltog = vltog;

    2729:   /* Set vltog for non-ghost vertices */
    2730:   k = 0;
    2731:   for (i=0; i<nroots; i++) {
    2732:     DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);
    2733:     if (ghost) continue;
    2734:     vltog[i] = vrange[rank] + k++;
    2735:   }
    2736:   PetscFree3(vrange,displs,recvcounts);

    2738:   /* Set vltog for ghost vertices */
    2739:   /* (a) create parallel Vleaves and sequential Vleaves_seq to convert local iremote[*].index to global index */
    2740:   VecCreate(comm,&Vleaves);
    2741:   VecSetSizes(Vleaves,2*nleaves,PETSC_DETERMINE);
    2742:   VecSetFromOptions(Vleaves);
    2743:   VecGetArray(Vleaves,&varr);
    2744:   for (i=0; i<nleaves; i++) {
    2745:     varr[2*i]   = (PetscScalar)(iremote[i].rank);  /* rank of remote process */
    2746:     varr[2*i+1] = (PetscScalar)(iremote[i].index); /* local index in remote process */
    2747:   }
    2748:   VecRestoreArray(Vleaves,&varr);

    2750:   /* (b) scatter local info to remote processes via VecScatter() */
    2751:   VecScatterCreateToAll(Vleaves,&ctx,&Vleaves_seq);
    2752:   VecScatterBegin(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);
    2753:   VecScatterEnd(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);

    2755:   /* (c) convert local indices to global indices in parallel vector Vleaves */
    2756:   VecGetSize(Vleaves_seq,&N);
    2757:   VecGetArrayRead(Vleaves_seq,&varr_read);
    2758:   for (i=0; i<N; i+=2) {
    2759:     remoterank = (PetscMPIInt)PetscRealPart(varr_read[i]);
    2760:     if (remoterank == rank) {
    2761:       k = i+1; /* row number */
    2762:       lidx = (PetscInt)PetscRealPart(varr_read[i+1]);
    2763:       val  = (PetscScalar)vltog[lidx]; /* global index for non-ghost vertex computed above */
    2764:       VecSetValues(Vleaves,1,&k,&val,INSERT_VALUES);
    2765:     }
    2766:   }
    2767:   VecRestoreArrayRead(Vleaves_seq,&varr_read);
    2768:   VecAssemblyBegin(Vleaves);
    2769:   VecAssemblyEnd(Vleaves);

    2771:   /* (d) Set vltog for ghost vertices by copying local values of Vleaves */
    2772:   VecGetArrayRead(Vleaves,&varr_read);
    2773:   k = 0;
    2774:   for (i=0; i<nroots; i++) {
    2775:     DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);
    2776:     if (!ghost) continue;
    2777:     vltog[i] = (PetscInt)PetscRealPart(varr_read[2*k+1]); k++;
    2778:   }
    2779:   VecRestoreArrayRead(Vleaves,&varr_read);

    2781:   VecDestroy(&Vleaves);
    2782:   VecDestroy(&Vleaves_seq);
    2783:   VecScatterDestroy(&ctx);
    2784:   return(0);
    2785: }