Actual source code: network.c

petsc-3.9.4 2018-09-11
Report Typos and Errors
  1:  #include <petsc/private/dmnetworkimpl.h>
  2:  #include <petscdmplex.h>
  3:  #include <petscsf.h>

  5: /*@
  6:   DMNetworkGetPlex - Gets the Plex DM associated with this network DM

  8:   Not collective

 10:   Input Parameters:
 11: + netdm - the dm object
 12: - plexmdm - the plex dm object

 14:   Level: Advanced

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

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

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

 30:   Collective on DM

 32:   Input Parameters:
 33: + dm - the dm object
 34: . Nsubnet - number of subnetworks
 35: . NsubnetCouple - number of coupling subnetworks
 36: . nV - number of local vertices for each subnetwork and coupling subnetwork
 37: . nE - number of local edges for each subnetwork and coupling subnetwork
 38: . NV - number of global vertices (or PETSC_DETERMINE) for each subnetwork and coupling subnetwork
 39: - NE - number of global edges (or PETSC_DETERMINE) for each subnetwork and coupling subnetwork

 41:    Notes
 42:    If one processor calls this with NV (NE) of PETSC_DECIDE then all processors must, otherwise the prgram will hang.

 44:    You cannot change the sizes once they have been set. nV, nE, NV and NE are arrays of length Nsubnet+NsubnetCouple.

 46:    Level: intermediate

 48: .seealso: DMNetworkCreate()
 49: @*/
 50: PetscErrorCode DMNetworkSetSizes(DM dm,PetscInt Nsubnet,PetscInt NsubnetCouple,PetscInt nV[], PetscInt nE[], PetscInt NV[], PetscInt NE[])
 51: {
 53:   DM_Network     *network = (DM_Network*) dm->data;
 54:   PetscInt       a[2],b[2],i;

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

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

 64:   network->nsubnet  = Nsubnet + NsubnetCouple;
 65:   Nsubnet          += NsubnetCouple;
 66:   network->ncsubnet = NsubnetCouple;
 67:   PetscCalloc1(Nsubnet,&network->subnet);
 68:   for(i=0; i < Nsubnet; i++) {
 69:     if (NV) {
 71:       if (NV[i] > 0 && nV[i] > NV[i]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Subnetwork %D: Local vertex size %D cannot be larger than global vertex size %D",i,nV[i],NV[i]);
 72:     }
 73:     if (NE) {
 75:       if (NE[i] > 0 && nE[i] > NE[i]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Subnetwork %D: Local edge size %D cannot be larger than global edge size %D",i,nE[i],NE[i]);
 76:     }

 78:     a[0] = nV[i]; a[1] = nE[i];
 79:     MPIU_Allreduce(a,b,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
 80:     network->subnet[i].Nvtx = b[0]; network->subnet[i].Nedge = b[1];

 82:     network->subnet[i].id = i;

 84:     network->subnet[i].nvtx = nV[i];
 85:     network->subnet[i].vStart = network->nVertices;
 86:     network->subnet[i].vEnd   = network->subnet[i].vStart + network->subnet[i].Nvtx;
 87:     network->nVertices += network->subnet[i].nvtx;
 88:     network->NVertices += network->subnet[i].Nvtx;

 90:     network->subnet[i].nedge = nE[i];
 91:     network->subnet[i].eStart = network->nEdges;
 92:     network->subnet[i].eEnd = network->subnet[i].eStart + network->subnet[i].Nedge;
 93:     network->nEdges += network->subnet[i].nedge;
 94:     network->NEdges += network->subnet[i].Nedge;
 95:   }
 96:   return(0);
 97: }

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

102:   Logically collective on DM

104:   Input Parameters:
105: + dm - the dm object
106: . edgelist - list of edges for each subnetwork
107: - edgelistCouple - list of edges for each coupling subnetwork

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

113:   Level: intermediate

115:   Example usage:
116:   Consider the following 2 separate networks and a coupling network:

118: .vb
119:  network 0: v0 -> v1 -> v2 -> v3
120:  network 1: v1 -> v2 -> v0
121:  coupling network: network 1: v2 -> network 0: v0
122: .ve

124:  The resulting input
125:    edgelist[0] = [0 1 | 1 2 | 2 3];
126:    edgelist[1] = [1 2 | 2 0]
127:    edgelistCouple[0] = [(network)1 (v)2 (network)0 (v)0].

129: .seealso: DMNetworkCreate, DMNetworkSetSizes
130: @*/
131: PetscErrorCode DMNetworkSetEdgeList(DM dm,PetscInt *edgelist[],PetscInt *edgelistCouple[])
132: {
133:   DM_Network *network = (DM_Network*) dm->data;
134:   PetscInt   i,nsubnet,ncsubnet=network->ncsubnet;

137:   nsubnet = network->nsubnet - ncsubnet;
138:   for(i=0; i < nsubnet; i++) {
139:     network->subnet[i].edgelist = edgelist[i];
140:   }
141:   if (edgelistCouple) {
142:     PetscInt j;
143:     j = 0;
144:     nsubnet = network->nsubnet;
145:     while (i < nsubnet) network->subnet[i++].edgelist = edgelistCouple[j++];
146:   }
147:   return(0);
148: }

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

153:   Collective on DM

155:   Input Parameters
156: . DM - the dmnetwork object

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

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

164:   Level: intermediate

166: .seealso: DMNetworkSetSizes, DMNetworkSetEdgeList
167: @*/
168: PetscErrorCode DMNetworkLayoutSetUp(DM dm)
169: {
171:   DM_Network     *network = (DM_Network*)dm->data;
172:   PetscInt       dim = 1; /* One dimensional network */
173:   PetscInt       numCorners=2,spacedim=2;
174:   double         *vertexcoords=NULL;
175:   PetscInt       i,j,ndata,ctr=0,nsubnet;
176:   PetscInt       k,netid,vid;
177:   PetscInt       *edgelist_couple=NULL;

180:   if (network->nVertices) {
181:     PetscCalloc1(numCorners*network->nVertices,&vertexcoords);
182:   }

184:   /* Create the edgelist for the network by concatenating edgelists of the subnetworks */
185:   nsubnet = network->nsubnet - network->ncsubnet;
186:   PetscCalloc1(2*network->nEdges,&network->edges);
187:   for (i=0; i < nsubnet; i++) {
188:     for (j = 0; j < network->subnet[i].nedge; j++) {
189:       network->edges[2*ctr] = network->subnet[i].vStart + network->subnet[i].edgelist[2*j];
190:       network->edges[2*ctr+1] = network->subnet[i].vStart + network->subnet[i].edgelist[2*j+1];
191:       ctr++;
192:     }
193:   }

195:   i       = nsubnet; /* netid of coupling subnet */
196:   nsubnet = network->nsubnet;
197:   while (i < nsubnet) {
198:     edgelist_couple = network->subnet[i].edgelist;
199:     k = 0;
200:     for (j = 0; j < network->subnet[i].nedge; j++) {
201:       netid = edgelist_couple[k]; vid = edgelist_couple[k+1];
202:       network->edges[2*ctr] = network->subnet[netid].vStart + vid; k += 2;

204:       netid = edgelist_couple[k]; vid = edgelist_couple[k+1];
205:       network->edges[2*ctr+1] = network->subnet[netid].vStart + vid; k+=2;
206:       ctr++;
207:     }
208:     i++;
209:   }

211: #if 0
212:   for(i=0; i < network->nEdges; i++) {
213:     PetscPrintf(PETSC_COMM_SELF,"[%D %D]",network->edges[2*i],network->edges[2*i+1]);
214:     printf("\n");
215:   }
216: #endif

218: #if defined(PETSC_USE_64BIT_INDICES)
219:   {
220:     int      *edges;
221:     PetscInt ii;
222:     PetscMalloc1(network->nEdges*numCorners,&edges);
223:     for (ii=0; ii<network->nEdges*numCorners; ii++) {
224:       edges[ii] = (int) network->edges[ii];
225:     }
226:     DMPlexCreateFromCellList(PetscObjectComm((PetscObject)dm),dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,edges,spacedim,vertexcoords,&network->plex);
227:     PetscFree(edges);
228:   }
229: #else
230:     DMPlexCreateFromCellList(PetscObjectComm((PetscObject)dm),dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,network->edges,spacedim,vertexcoords,&network->plex);
231:  #endif

233:   if (network->nVertices) {
234:     PetscFree(vertexcoords);
235:   }
236:   PetscFree(network->edges);

238:   DMPlexGetChart(network->plex,&network->pStart,&network->pEnd);
239:   DMPlexGetHeightStratum(network->plex,0,&network->eStart,&network->eEnd);
240:   DMPlexGetHeightStratum(network->plex,1,&network->vStart,&network->vEnd);

242:   PetscSectionCreate(PetscObjectComm((PetscObject)dm),&network->DataSection);
243:   PetscSectionCreate(PetscObjectComm((PetscObject)dm),&network->DofSection);
244:   PetscSectionSetChart(network->DataSection,network->pStart,network->pEnd);
245:   PetscSectionSetChart(network->DofSection,network->pStart,network->pEnd);

247:   /* Create vertices and edges array for the subnetworks */
248:   for (j=0; j < network->nsubnet; j++) {
249:     PetscCalloc1(network->subnet[j].nedge,&network->subnet[j].edges);
250:     PetscCalloc1(network->subnet[j].nvtx,&network->subnet[j].vertices);
251:     /* Temporarily setting nvtx and nedge to 0 so we can use them as counters in the below for loop.
252:        These get updated when the vertices and edges are added. */
253:     network->subnet[j].nvtx = network->subnet[j].nedge = 0;
254:   }

256:   network->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
257:   PetscCalloc1(network->pEnd-network->pStart,&network->header);
258:   for (i=network->eStart; i < network->eEnd; i++) {
259:     network->header[i].index = i;   /* Global edge number */
260:     for (j=0; j < network->nsubnet; j++) {
261:       if((network->subnet[j].eStart <= i) && (i < network->subnet[j].eEnd)) {
262:         network->header[i].subnetid = j; /* Subnetwork id */
263:         network->subnet[j].edges[network->subnet[j].nedge++] = i;
264:         break;
265:       }
266:     }
267:     network->header[i].ndata = 0;
268:     ndata = network->header[i].ndata;
269:     PetscSectionAddDof(network->DataSection,i,network->dataheadersize);
270:     network->header[i].offset[ndata] = 0;
271:   }

273:   for(i=network->vStart; i < network->vEnd; i++) {
274:     network->header[i].index = i - network->vStart;
275:     for (j=0; j < network->nsubnet; j++) {
276:       if ((network->subnet[j].vStart <= i-network->vStart) && (i-network->vStart < network->subnet[j].vEnd)) {
277:         network->header[i].subnetid = j;
278:         network->subnet[j].vertices[network->subnet[j].nvtx++] = i;
279:         break;
280:       }
281:     }
282:     network->header[i].ndata = 0;
283:     ndata = network->header[i].ndata;
284:     PetscSectionAddDof(network->DataSection,i,network->dataheadersize);
285:     network->header[i].offset[ndata] = 0;
286:   }

288:   PetscMalloc1(network->pEnd-network->pStart,&network->cvalue);
289:   return(0);
290: }

292: /*@C
293:   DMNetworkGetSubnetworkInfo - Returns the info for the subnetwork

295:   Input Parameters
296: + dm   - the number object
297: - id   - the ID (integer) of the subnetwork

299:   Output Parameters
300: + nv    - number of vertices (local)
301: . ne    - number of edges (local)
302: . vtx   - local vertices for this subnetwork
303: . edge  - local edges for this subnetwork

305:   Notes:
306:   Cannot call this routine before DMNetworkLayoutSetup()

308: .seealso: DMNetworkLayoutSetUp, DMNetworkCreate
309: @*/
310: PetscErrorCode DMNetworkGetSubnetworkInfo(DM netdm,PetscInt id,PetscInt *nv, PetscInt *ne,const PetscInt **vtx, const PetscInt **edge)
311: {
312:   DM_Network     *network = (DM_Network*) netdm->data;

315:   *nv = network->subnet[id].nvtx;
316:   *ne = network->subnet[id].nedge;
317:   *vtx = network->subnet[id].vertices;
318:   *edge = network->subnet[id].edges;
319:   return(0);
320: }

322: /*@C
323:   DMNetworkRegisterComponent - Registers the network component

325:   Logically collective on DM

327:   Input Parameters
328: + dm   - the network object
329: . name - the component name
330: - size - the storage size in bytes for this component data

332:    Output Parameters
333: .   key - an integer key that defines the component

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

338:    Level: intermediate

340: .seealso: DMNetworkLayoutSetUp, DMNetworkCreate
341: @*/
342: PetscErrorCode DMNetworkRegisterComponent(DM dm,const char *name,PetscInt size,PetscInt *key)
343: {
344:   PetscErrorCode        ierr;
345:   DM_Network            *network = (DM_Network*) dm->data;
346:   DMNetworkComponent    *component=&network->component[network->ncomponent];
347:   PetscBool             flg=PETSC_FALSE;
348:   PetscInt              i;

351:   for (i=0; i < network->ncomponent; i++) {
352:     PetscStrcmp(component->name,name,&flg);
353:     if (flg) {
354:       *key = i;
355:       return(0);
356:     }
357:   }
358:   if(network->ncomponent == MAX_COMPONENTS) {
359:     SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Number of components registered exceeds the max %D",MAX_COMPONENTS);
360:   }

362:   PetscStrcpy(component->name,name);
363:   component->size = size/sizeof(DMNetworkComponentGenericDataType);
364:   *key = network->ncomponent;
365:   network->ncomponent++;
366:   return(0);
367: }

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

372:   Not Collective

374:   Input Parameters:
375: + dm - The DMNetwork object

377:   Output Paramters:
378: + vStart - The first vertex point
379: - vEnd   - One beyond the last vertex point

381:   Level: intermediate

383: .seealso: DMNetworkGetEdgeRange
384: @*/
385: PetscErrorCode DMNetworkGetVertexRange(DM dm,PetscInt *vStart,PetscInt *vEnd)
386: {
387:   DM_Network     *network = (DM_Network*)dm->data;

390:   if (vStart) *vStart = network->vStart;
391:   if (vEnd) *vEnd = network->vEnd;
392:   return(0);
393: }

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

398:   Not Collective

400:   Input Parameters:
401: + dm - The DMNetwork object

403:   Output Paramters:
404: + eStart - The first edge point
405: - eEnd   - One beyond the last edge point

407:   Level: intermediate

409: .seealso: DMNetworkGetVertexRange
410: @*/
411: PetscErrorCode DMNetworkGetEdgeRange(DM dm,PetscInt *eStart,PetscInt *eEnd)
412: {
413:   DM_Network     *network = (DM_Network*)dm->data;

416:   if (eStart) *eStart = network->eStart;
417:   if (eEnd) *eEnd = network->eEnd;
418:   return(0);
419: }

421: /*@
422:   DMNetworkGetGlobalEdgeIndex - Get the user global numbering for the edge.

424:   Not Collective

426:   Input Parameters:
427: + dm - DMNetwork object
428: - p  - edge point

430:   Output Paramters:
431: . index - user global numbering for the edge

433:   Level: intermediate

435: .seealso: DMNetworkGetGlobalVertexIndex
436: @*/
437: PetscErrorCode DMNetworkGetGlobalEdgeIndex(DM dm,PetscInt p,PetscInt *index)
438: {
439:   PetscErrorCode    ierr;
440:   DM_Network        *network = (DM_Network*)dm->data;
441:   PetscInt          offsetp;
442:   DMNetworkComponentHeader header;

445:   PetscSectionGetOffset(network->DataSection,p,&offsetp);
446:   header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
447:   *index = header->index;
448:   return(0);
449: }

451: /*@
452:   DMNetworkGetGlobalVertexIndex - Get the user global numbering for the vertex.

454:   Not Collective

456:   Input Parameters:
457: + dm - DMNetwork object
458: - p  - vertex point

460:   Output Paramters:
461: . index - user global numbering for the vertex

463:   Level: intermediate

465: .seealso: DMNetworkGetGlobalEdgeIndex
466: @*/
467: PetscErrorCode DMNetworkGetGlobalVertexIndex(DM dm,PetscInt p,PetscInt *index)
468: {
469:   PetscErrorCode    ierr;
470:   DM_Network        *network = (DM_Network*)dm->data;
471:   PetscInt          offsetp;
472:   DMNetworkComponentHeader header;

475:   PetscSectionGetOffset(network->DataSection,p,&offsetp);
476:   header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
477:   *index = header->index;
478:   return(0);
479: }

481: /*
482:   DMNetworkGetComponentKeyOffset - Gets the type along with the offset for indexing the
483:                                     component value from the component data array

485:   Not Collective

487:   Input Parameters:
488: + dm      - The DMNetwork object
489: . p       - vertex/edge point
490: - compnum - component number

492:   Output Parameters:
493: + compkey - the key obtained when registering the component
494: - offset  - offset into the component data array associated with the vertex/edge point

496:   Notes:
497:   Typical usage:

499:   DMNetworkGetComponentDataArray(dm, &arr);
500:   DMNetworkGetVertex/EdgeRange(dm,&Start,&End);
501:   Loop over vertices or edges
502:     DMNetworkGetNumComponents(dm,v,&numcomps);
503:     Loop over numcomps
504:       DMNetworkGetComponentKeyOffset(dm,v,compnum,&key,&offset);
505:       compdata = (UserCompDataType)(arr+offset);

507:   Level: intermediate

509: .seealso: DMNetworkGetNumComponents, DMNetworkGetComponentDataArray,
510: */
511: PetscErrorCode DMNetworkGetComponentKeyOffset(DM dm,PetscInt p, PetscInt compnum, PetscInt *compkey, PetscInt *offset)
512: {
513:   PetscErrorCode           ierr;
514:   PetscInt                 offsetp;
515:   DMNetworkComponentHeader header;
516:   DM_Network               *network = (DM_Network*)dm->data;

519:   PetscSectionGetOffset(network->DataSection,p,&offsetp);
520:   header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
521:   if (compkey) *compkey = header->key[compnum];
522:   if (offset) *offset  = offsetp+network->dataheadersize+header->offset[compnum];
523:   return(0);
524: }

526: /*@
527:   DMNetworkGetComponent - Returns the network component and its key

529:   Not Collective

531:   Input Parameters
532: + dm - DMNetwork object
533: . p  - edge or vertex point
534: - compnum - component number

536:   Output Parameters:
537: + compkey - the key set for this computing during registration
538: - component - the component data

540:   Notes:
541:   Typical usage:

543:   DMNetworkGetVertex/EdgeRange(dm,&Start,&End);
544:   Loop over vertices or edges
545:     DMNetworkGetNumComponents(dm,v,&numcomps);
546:     Loop over numcomps
547:       DMNetworkGetComponent(dm,v,compnum,&key,&component);

549:   Level: intermediate

551: .seealso: DMNetworkGetNumComponents, DMNetworkGetVariableOffset
552: @*/
553: PetscErrorCode DMNetworkGetComponent(DM dm, PetscInt p, PetscInt compnum, PetscInt *key, void **component)
554: {
556:   DM_Network     *network = (DM_Network*)dm->data;
557:   PetscInt       offsetd = 0;

560:   DMNetworkGetComponentKeyOffset(dm,p,compnum,key,&offsetd);
561:   *component = network->componentdataarray+offsetd;
562:   return(0);
563: }

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

568:   Not Collective

570:   Input Parameters:
571: + dm           - The DMNetwork object
572: . p            - vertex/edge point
573: . componentkey - component key returned while registering the component
574: - compvalue    - pointer to the data structure for the component

576:   Level: intermediate

578: .seealso: DMNetworkGetVertexRange, DMNetworkGetEdgeRange, DMNetworkRegisterComponent
579: @*/
580: PetscErrorCode DMNetworkAddComponent(DM dm, PetscInt p,PetscInt componentkey,void* compvalue)
581: {
582:   DM_Network               *network = (DM_Network*)dm->data;
583:   DMNetworkComponent       *component = &network->component[componentkey];
584:   DMNetworkComponentHeader header = &network->header[p];
585:   DMNetworkComponentValue  cvalue = &network->cvalue[p];
586:   PetscErrorCode           ierr;

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

591:   header->size[header->ndata] = component->size;
592:   PetscSectionAddDof(network->DataSection,p,component->size);
593:   header->key[header->ndata] = componentkey;
594:   if (header->ndata != 0) header->offset[header->ndata] = header->offset[header->ndata-1] + header->size[header->ndata-1];

596:   cvalue->data[header->ndata] = (void*)compvalue;
597:   header->ndata++;
598:   return(0);
599: }

601: /*@
602:   DMNetworkGetNumComponents - Get the number of components at a vertex/edge

604:   Not Collective

606:   Input Parameters:
607: + dm - The DMNetwork object
608: . p  - vertex/edge point

610:   Output Parameters:
611: . numcomponents - Number of components at the vertex/edge

613:   Level: intermediate

615: .seealso: DMNetworkRegisterComponent, DMNetworkAddComponent
616: @*/
617: PetscErrorCode DMNetworkGetNumComponents(DM dm,PetscInt p,PetscInt *numcomponents)
618: {
620:   PetscInt       offset;
621:   DM_Network     *network = (DM_Network*)dm->data;

624:   PetscSectionGetOffset(network->DataSection,p,&offset);
625:   *numcomponents = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata;
626:   return(0);
627: }

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

632:   Not Collective

634:   Input Parameters:
635: + dm     - The DMNetwork object
636: - p      - the edge/vertex point

638:   Output Parameters:
639: . offset - the offset

641:   Level: intermediate

643: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
644: @*/
645: PetscErrorCode DMNetworkGetVariableOffset(DM dm,PetscInt p,PetscInt *offset)
646: {
648:   DM_Network     *network = (DM_Network*)dm->data;

651:   PetscSectionGetOffset(network->plex->defaultSection,p,offset);
652:   return(0);
653: }

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

658:   Not Collective

660:   Input Parameters:
661: + dm      - The DMNetwork object
662: - p       - the edge/vertex point

664:   Output Parameters:
665: . offsetg - the offset

667:   Level: intermediate

669: .seealso: DMNetworkGetVariableOffset, DMGetLocalVector
670: @*/
671: PetscErrorCode DMNetworkGetVariableGlobalOffset(DM dm,PetscInt p,PetscInt *offsetg)
672: {
674:   DM_Network     *network = (DM_Network*)dm->data;

677:   PetscSectionGetOffset(network->plex->defaultGlobalSection,p,offsetg);
678:   if (*offsetg < 0) *offsetg = -(*offsetg + 1); /* Convert to actual global offset for ghost vertex */
679:   return(0);
680: }

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

685:   Not Collective

687:   Input Parameters:
688: + dm     - The DMNetwork object
689: - p      - the edge point

691:   Output Parameters:
692: . offset - the offset

694:   Level: intermediate

696: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
697: @*/
698: PetscErrorCode DMNetworkGetEdgeOffset(DM dm,PetscInt p,PetscInt *offset)
699: {
701:   DM_Network     *network = (DM_Network*)dm->data;


705:   PetscSectionGetOffset(network->edge.DofSection,p,offset);
706:   return(0);
707: }

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

712:   Not Collective

714:   Input Parameters:
715: + dm     - The DMNetwork object
716: - p      - the vertex point

718:   Output Parameters:
719: . offset - the offset

721:   Level: intermediate

723: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
724: @*/
725: PetscErrorCode DMNetworkGetVertexOffset(DM dm,PetscInt p,PetscInt *offset)
726: {
728:   DM_Network     *network = (DM_Network*)dm->data;


732:   p -= network->vStart;

734:   PetscSectionGetOffset(network->vertex.DofSection,p,offset);
735:   return(0);
736: }
737: /*@
738:   DMNetworkAddNumVariables - Add number of variables associated with a given point.

740:   Not Collective

742:   Input Parameters:
743: + dm   - The DMNetworkObject
744: . p    - the vertex/edge point
745: - nvar - number of additional variables

747:   Level: intermediate

749: .seealso: DMNetworkSetNumVariables
750: @*/
751: PetscErrorCode DMNetworkAddNumVariables(DM dm,PetscInt p,PetscInt nvar)
752: {
754:   DM_Network     *network = (DM_Network*)dm->data;

757:   PetscSectionAddDof(network->DofSection,p,nvar);
758:   return(0);
759: }

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

764:   Not Collective

766:   Input Parameters:
767: + dm   - The DMNetworkObject
768: - p    - the vertex/edge point

770:   Output Parameters:
771: . nvar - number of variables

773:   Level: intermediate

775: .seealso: DMNetworkAddNumVariables, DMNetworkSddNumVariables
776: @*/
777: PetscErrorCode DMNetworkGetNumVariables(DM dm,PetscInt p,PetscInt *nvar)
778: {
780:   DM_Network     *network = (DM_Network*)dm->data;

783:   PetscSectionGetDof(network->DofSection,p,nvar);
784:   return(0);
785: }

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

790:   Not Collective

792:   Input Parameters:
793: + dm   - The DMNetworkObject
794: . p    - the vertex/edge point
795: - nvar - number of variables

797:   Level: intermediate

799: .seealso: DMNetworkAddNumVariables
800: @*/
801: PetscErrorCode DMNetworkSetNumVariables(DM dm,PetscInt p,PetscInt nvar)
802: {
804:   DM_Network     *network = (DM_Network*)dm->data;

807:   PetscSectionSetDof(network->DofSection,p,nvar);
808:   return(0);
809: }

811: /* Sets up the array that holds the data for all components and its associated section. This
812:    function is called during DMSetUp() */
813: PetscErrorCode DMNetworkComponentSetUp(DM dm)
814: {
815:   PetscErrorCode           ierr;
816:   DM_Network               *network = (DM_Network*)dm->data;
817:   PetscInt                 arr_size,p,offset,offsetp,ncomp,i;
818:   DMNetworkComponentHeader header;
819:   DMNetworkComponentValue  cvalue;
820:   DMNetworkComponentGenericDataType *componentdataarray;

823:   PetscSectionSetUp(network->DataSection);
824:   PetscSectionGetStorageSize(network->DataSection,&arr_size);
825:   PetscMalloc1(arr_size,&network->componentdataarray);
826:   componentdataarray = network->componentdataarray;
827:   for (p = network->pStart; p < network->pEnd; p++) {
828:     PetscSectionGetOffset(network->DataSection,p,&offsetp);
829:     /* Copy header */
830:     header = &network->header[p];
831:     PetscMemcpy(componentdataarray+offsetp,header,network->dataheadersize*sizeof(DMNetworkComponentGenericDataType));
832:     /* Copy data */
833:     cvalue = &network->cvalue[p];
834:     ncomp = header->ndata;
835:     for (i = 0; i < ncomp; i++) {
836:       offset = offsetp + network->dataheadersize + header->offset[i];
837:       PetscMemcpy(componentdataarray+offset,cvalue->data[i],header->size[i]*sizeof(DMNetworkComponentGenericDataType));
838:     }
839:   }
840:   return(0);
841: }

843: /* Sets up the section for dofs. This routine is called during DMSetUp() */
844: PetscErrorCode DMNetworkVariablesSetUp(DM dm)
845: {
847:   DM_Network     *network = (DM_Network*)dm->data;

850:   PetscSectionSetUp(network->DofSection);
851:   return(0);
852: }

854: /*@C
855:   DMNetworkGetComponentDataArray - Returns the component data array

857:   Not Collective

859:   Input Parameters:
860: . dm - The DMNetwork Object

862:   Output Parameters:
863: . componentdataarray - array that holds data for all components

865:   Level: intermediate

867: .seealso: DMNetworkGetComponentKeyOffset, DMNetworkGetNumComponents
868: @*/
869: PetscErrorCode DMNetworkGetComponentDataArray(DM dm,DMNetworkComponentGenericDataType **componentdataarray)
870: {
871:   DM_Network     *network = (DM_Network*)dm->data;

874:   *componentdataarray = network->componentdataarray;
875:   return(0);
876: }

878: /* Get a subsection from a range of points */
879: PetscErrorCode DMNetworkGetSubSection_private(PetscSection master, PetscInt pstart, PetscInt pend,PetscSection *subsection)
880: {
882:   PetscInt       i, nvar;

885:   PetscSectionCreate(PetscObjectComm((PetscObject)master), subsection);
886:   PetscSectionSetChart(*subsection, 0, pend - pstart);
887:   for (i = pstart; i < pend; i++) {
888:     PetscSectionGetDof(master,i,&nvar);
889:     PetscSectionSetDof(*subsection, i - pstart, nvar);
890:   }

892:   PetscSectionSetUp(*subsection);
893:   return(0);
894: }

896: /* Create a submap of points with a GlobalToLocal structure */
897: PetscErrorCode DMNetworkSetSubMap_private(PetscInt pstart, PetscInt pend, ISLocalToGlobalMapping *map)
898: {
900:   PetscInt       i, *subpoints;

903:   /* Create index sets to map from "points" to "subpoints" */
904:   PetscMalloc1(pend - pstart, &subpoints);
905:   for (i = pstart; i < pend; i++) {
906:     subpoints[i - pstart] = i;
907:   }
908:   ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,pend-pstart,subpoints,PETSC_COPY_VALUES,map);
909:   PetscFree(subpoints);
910:   return(0);
911: }

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

916:   Collective

918:   Input Parameters:
919: . dm   - The DMNetworkObject

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

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

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

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

929:   Level: intermediate

931: @*/
932: PetscErrorCode DMNetworkAssembleGraphStructures(DM dm)
933: {
935:   MPI_Comm       comm;
936:   PetscMPIInt    rank, size;
937:   DM_Network     *network = (DM_Network*)dm->data;

940:   PetscObjectGetComm((PetscObject)dm,&comm);
941:   MPI_Comm_rank(comm, &rank);
942:   MPI_Comm_size(comm, &size);

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

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

952:   if (size > 1) {
953:     PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);
954:     PetscSectionCreateGlobalSection(network->vertex.DofSection, network->vertex.sf, PETSC_FALSE, PETSC_FALSE, &network->vertex.GlobalDofSection);
955:   PetscSFGetSubSF(network->plex->sf, network->edge.mapping, &network->edge.sf);
956:   PetscSectionCreateGlobalSection(network->edge.DofSection, network->edge.sf, PETSC_FALSE, PETSC_FALSE, &network->edge.GlobalDofSection);
957:   } else {
958:   /* create structures for vertex */
959:   PetscSectionClone(network->vertex.DofSection,&network->vertex.GlobalDofSection);
960:   /* create structures for edge */
961:   PetscSectionClone(network->edge.DofSection,&network->edge.GlobalDofSection);
962:   }


965:   /* Add viewers */
966:   PetscObjectSetName((PetscObject)network->edge.GlobalDofSection,"Global edge dof section");
967:   PetscObjectSetName((PetscObject)network->vertex.GlobalDofSection,"Global vertex dof section");
968:   PetscSectionViewFromOptions(network->edge.GlobalDofSection, NULL, "-edge_global_section_view");
969:   PetscSectionViewFromOptions(network->vertex.GlobalDofSection, NULL, "-vertex_global_section_view");

971:   return(0);
972: }

974: /*@
975:   DMNetworkDistribute - Distributes the network and moves associated component data.

977:   Collective

979:   Input Parameter:
980: + DM - the DMNetwork object
981: - overlap - The overlap of partitions, 0 is the default

983:   Notes:
984:   Distributes the network with <overlap>-overlapping partitioning of the edges.

986:   Level: intermediate

988: .seealso: DMNetworkCreate
989: @*/
990: PetscErrorCode DMNetworkDistribute(DM *dm,PetscInt overlap)
991: {
992:   MPI_Comm       comm;
994:   PetscMPIInt    size;
995:   DM_Network     *oldDMnetwork = (DM_Network*)((*dm)->data);
996:   DM_Network     *newDMnetwork;
997:   PetscSF        pointsf;
998:   DM             newDM;
999:   PetscPartitioner part;
1000:   PetscInt         j,e,v,offset;
1001:   DMNetworkComponentHeader header;


1005:   PetscObjectGetComm((PetscObject)*dm,&comm);
1006:   MPI_Comm_size(comm, &size);
1007:   if (size == 1) return(0);

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

1013:   /* Enable runtime options for petscpartitioner */
1014:   DMPlexGetPartitioner(oldDMnetwork->plex,&part);
1015:   PetscPartitionerSetFromOptions(part);

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

1020:   /* Distribute dof section */
1021:   PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DofSection);
1022:   PetscSFDistributeSection(pointsf,oldDMnetwork->DofSection,NULL,newDMnetwork->DofSection);
1023:   PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DataSection);

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

1028:   PetscSectionGetChart(newDMnetwork->DataSection,&newDMnetwork->pStart,&newDMnetwork->pEnd);
1029:   DMPlexGetHeightStratum(newDMnetwork->plex,0, &newDMnetwork->eStart,&newDMnetwork->eEnd);
1030:   DMPlexGetHeightStratum(newDMnetwork->plex,1,&newDMnetwork->vStart,&newDMnetwork->vEnd);
1031:   newDMnetwork->nEdges = newDMnetwork->eEnd - newDMnetwork->eStart;
1032:   newDMnetwork->nVertices = newDMnetwork->vEnd - newDMnetwork->vStart;
1033:   newDMnetwork->NVertices = oldDMnetwork->NVertices;
1034:   newDMnetwork->NEdges = oldDMnetwork->NEdges;

1036:   /* Set Dof section as the default section for dm */
1037:   DMSetDefaultSection(newDMnetwork->plex,newDMnetwork->DofSection);
1038:   DMGetDefaultGlobalSection(newDMnetwork->plex,&newDMnetwork->GlobalDofSection);

1040:   /* Set up subnetwork info in the newDM */
1041:   newDMnetwork->nsubnet = oldDMnetwork->nsubnet;
1042:   PetscCalloc1(newDMnetwork->nsubnet,&newDMnetwork->subnet);
1043:   /* Copy over the global number of vertices and edges in each subnetwork. Note that these are already
1044:      calculated in DMNetworkLayoutSetUp()
1045:   */
1046:   for(j=0; j < newDMnetwork->nsubnet; j++) {
1047:     newDMnetwork->subnet[j].Nvtx = oldDMnetwork->subnet[j].Nvtx;
1048:     newDMnetwork->subnet[j].Nedge = oldDMnetwork->subnet[j].Nedge;
1049:   }

1051:   for(e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++ ) {
1052:     PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);
1053:     header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1054:     newDMnetwork->subnet[header->subnetid].nedge++;
1055:   }

1057:   for(v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++ ) {
1058:     PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);
1059:     header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1060:     newDMnetwork->subnet[header->subnetid].nvtx++;
1061:   }

1063:   /* Now create the vertices and edge arrays for the subnetworks */
1064:   for(j=0; j < newDMnetwork->nsubnet; j++) {
1065:     PetscCalloc1(newDMnetwork->subnet[j].nedge,&newDMnetwork->subnet[j].edges);
1066:     PetscCalloc1(newDMnetwork->subnet[j].nvtx,&newDMnetwork->subnet[j].vertices);
1067:     /* Temporarily setting nvtx and nedge to 0 so we can use them as counters in the below for loop.
1068:        These get updated when the vertices and edges are added. */
1069:     newDMnetwork->subnet[j].nvtx = newDMnetwork->subnet[j].nedge = 0;
1070:   }

1072:   /* Set the vertices and edges in each subnetwork */
1073:   for(e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++ ) {
1074:     PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);
1075:     header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1076:     newDMnetwork->subnet[header->subnetid].edges[newDMnetwork->subnet[header->subnetid].nedge++]  = e;
1077:   }

1079:   for(v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++ ) {
1080:     PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);
1081:     header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1082:     newDMnetwork->subnet[header->subnetid].vertices[newDMnetwork->subnet[header->subnetid].nvtx++]  = v;
1083:   }

1085:   /* Destroy point SF */
1086:   PetscSFDestroy(&pointsf);

1088:   DMDestroy(dm);
1089:   *dm  = newDM;
1090:   return(0);
1091: }

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

1096:   Input Parameters:
1097: + masterSF - the original SF structure
1098: - map      - a ISLocalToGlobal mapping that contains the subset of points

1100:   Output Parameters:
1101: . subSF    - a subset of the masterSF for the desired subset.
1102: */

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

1106:   PetscErrorCode        ierr;
1107:   PetscInt              nroots, nleaves, *ilocal_sub;
1108:   PetscInt              i, *ilocal_map, nroots_sub, nleaves_sub = 0;
1109:   PetscInt              *local_points, *remote_points;
1110:   PetscSFNode           *iremote_sub;
1111:   const PetscInt        *ilocal;
1112:   const PetscSFNode     *iremote;

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

1117:   /* Look for leaves that pertain to the subset of points. Get the local ordering */
1118:   PetscMalloc1(nleaves,&ilocal_map);
1119:   ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nleaves,ilocal,NULL,ilocal_map);
1120:   for (i = 0; i < nleaves; i++) {
1121:     if (ilocal_map[i] != -1) nleaves_sub += 1;
1122:   }
1123:   /* Re-number ilocal with subset numbering. Need information from roots */
1124:   PetscMalloc2(nroots,&local_points,nroots,&remote_points);
1125:   for (i = 0; i < nroots; i++) local_points[i] = i;
1126:   ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nroots,local_points,NULL,local_points);
1127:   PetscSFBcastBegin(mastersf, MPIU_INT, local_points, remote_points);
1128:   PetscSFBcastEnd(mastersf, MPIU_INT, local_points, remote_points);
1129:   /* Fill up graph using local (that is, local to the subset) numbering. */
1130:   PetscMalloc1(nleaves_sub,&ilocal_sub);
1131:   PetscMalloc1(nleaves_sub,&iremote_sub);
1132:   nleaves_sub = 0;
1133:   for (i = 0; i < nleaves; i++) {
1134:     if (ilocal_map[i] != -1) {
1135:       ilocal_sub[nleaves_sub] = ilocal_map[i];
1136:       iremote_sub[nleaves_sub].rank = iremote[i].rank;
1137:       iremote_sub[nleaves_sub].index = remote_points[ilocal[i]];
1138:       nleaves_sub += 1;
1139:     }
1140:   }
1141:   PetscFree2(local_points,remote_points);
1142:   ISLocalToGlobalMappingGetSize(map,&nroots_sub);

1144:   /* Create new subSF */
1145:   PetscSFCreate(PETSC_COMM_WORLD,subSF);
1146:   PetscSFSetFromOptions(*subSF);
1147:   PetscSFSetGraph(*subSF,nroots_sub,nleaves_sub,ilocal_sub,PETSC_OWN_POINTER,iremote_sub,PETSC_COPY_VALUES);
1148:   PetscFree(ilocal_map);
1149:   PetscFree(iremote_sub);
1150:   return(0);
1151: }

1153: /*@C
1154:   DMNetworkGetSupportingEdges - Return the supporting edges for this vertex point

1156:   Not Collective

1158:   Input Parameters:
1159: + dm - The DMNetwork object
1160: - p  - the vertex point

1162:   Output Paramters:
1163: + nedges - number of edges connected to this vertex point
1164: - edges  - List of edge points

1166:   Level: intermediate

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

1172: .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices
1173: @*/
1174: PetscErrorCode DMNetworkGetSupportingEdges(DM dm,PetscInt vertex,PetscInt *nedges,const PetscInt *edges[])
1175: {
1177:   DM_Network     *network = (DM_Network*)dm->data;

1180:   DMPlexGetSupportSize(network->plex,vertex,nedges);
1181:   DMPlexGetSupport(network->plex,vertex,edges);
1182:   return(0);
1183: }

1185: /*@C
1186:   DMNetworkGetConnectedVertices - Return the connected vertices for this edge point

1188:   Not Collective

1190:   Input Parameters:
1191: + dm - The DMNetwork object
1192: - p  - the edge point

1194:   Output Paramters:
1195: . vertices  - vertices connected to this edge

1197:   Level: intermediate

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

1203: .seealso: DMNetworkCreate, DMNetworkGetSupportingEdges
1204: @*/
1205: PetscErrorCode DMNetworkGetConnectedVertices(DM dm,PetscInt edge,const PetscInt *vertices[])
1206: {
1208:   DM_Network     *network = (DM_Network*)dm->data;

1211:   DMPlexGetCone(network->plex,edge,vertices);
1212:   return(0);
1213: }

1215: /*@
1216:   DMNetworkIsGhostVertex - Returns TRUE if the vertex is a ghost vertex

1218:   Not Collective

1220:   Input Parameters:
1221: + dm - The DMNetwork object
1222: . p  - the vertex point

1224:   Output Parameter:
1225: . isghost - TRUE if the vertex is a ghost point

1227:   Level: intermediate

1229: .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices, DMNetworkGetVertexRange
1230: @*/
1231: PetscErrorCode DMNetworkIsGhostVertex(DM dm,PetscInt p,PetscBool *isghost)
1232: {
1234:   DM_Network     *network = (DM_Network*)dm->data;
1235:   PetscInt       offsetg;
1236:   PetscSection   sectiong;

1239:   *isghost = PETSC_FALSE;
1240:   DMGetDefaultGlobalSection(network->plex,&sectiong);
1241:   PetscSectionGetOffset(sectiong,p,&offsetg);
1242:   if (offsetg < 0) *isghost = PETSC_TRUE;
1243:   return(0);
1244: }

1246: PetscErrorCode DMSetUp_Network(DM dm)
1247: {
1249:   DM_Network     *network=(DM_Network*)dm->data;

1252:   DMNetworkComponentSetUp(dm);
1253:   DMNetworkVariablesSetUp(dm);

1255:   DMSetDefaultSection(network->plex,network->DofSection);
1256:   DMGetDefaultGlobalSection(network->plex,&network->GlobalDofSection);
1257:   return(0);
1258: }

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

1264:     Collective

1266:     Input Parameters:
1267: +   dm - The DMNetwork object
1268: .   eflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for edges
1269: -   vflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for vertices

1271:     Level: intermediate

1273: @*/
1274: PetscErrorCode DMNetworkHasJacobian(DM dm,PetscBool eflg,PetscBool vflg)
1275: {
1276:   DM_Network     *network=(DM_Network*)dm->data;

1280:   network->userEdgeJacobian   = eflg;
1281:   network->userVertexJacobian = vflg;

1283:   if (eflg && !network->Je) {
1284:     PetscCalloc1(3*network->nEdges,&network->Je);
1285:   }

1287:   if (vflg && !network->Jv) {
1288:     PetscInt       i,*vptr,nedges,vStart=network->vStart;
1289:     PetscInt       nVertices = network->nVertices,nedges_total;
1290:     const PetscInt *edges;

1292:     /* count nvertex_total */
1293:     nedges_total = 0;
1294:     PetscMalloc1(nVertices+1,&vptr);

1296:     vptr[0] = 0;
1297:     for (i=0; i<nVertices; i++) {
1298:       DMNetworkGetSupportingEdges(dm,i+vStart,&nedges,&edges);
1299:       nedges_total += nedges;
1300:       vptr[i+1] = vptr[i] + 2*nedges + 1;
1301:     }

1303:     PetscCalloc1(2*nedges_total+nVertices,&network->Jv);
1304:     network->Jvptr = vptr;
1305:   }
1306:   return(0);
1307: }

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

1312:     Not Collective

1314:     Input Parameters:
1315: +   dm - The DMNetwork object
1316: .   p  - the edge point
1317: -   J - array (size = 3) of Jacobian submatrices for this edge point:
1318:         J[0]: this edge
1319:         J[1] and J[2]: connected vertices, obtained by calling DMNetworkGetConnectedVertices()

1321:     Level: intermediate

1323: .seealso: DMNetworkVertexSetMatrix
1324: @*/
1325: PetscErrorCode DMNetworkEdgeSetMatrix(DM dm,PetscInt p,Mat J[])
1326: {
1327:   DM_Network     *network=(DM_Network*)dm->data;

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

1332:   if (J) {
1333:     network->Je[3*p]   = J[0];
1334:     network->Je[3*p+1] = J[1];
1335:     network->Je[3*p+2] = J[2];
1336:   }
1337:   return(0);
1338: }

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

1343:     Not Collective

1345:     Input Parameters:
1346: +   dm - The DMNetwork object
1347: .   p  - the vertex point
1348: -   J - array of Jacobian (size = 2*(num of supporting edges) + 1) submatrices for this vertex point:
1349:         J[0]:       this vertex
1350:         J[1+2*i]:   i-th supporting edge
1351:         J[1+2*i+1]: i-th connected vertex

1353:     Level: intermediate

1355: .seealso: DMNetworkEdgeSetMatrix
1356: @*/
1357: PetscErrorCode DMNetworkVertexSetMatrix(DM dm,PetscInt p,Mat J[])
1358: {
1360:   DM_Network     *network=(DM_Network*)dm->data;
1361:   PetscInt       i,*vptr,nedges,vStart=network->vStart;
1362:   const PetscInt *edges;

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

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

1371:     /* Set Jacobian for each supporting edge and connected vertex */
1372:     DMNetworkGetSupportingEdges(dm,p,&nedges,&edges);
1373:     for (i=1; i<=2*nedges; i++) network->Jv[vptr[p-vStart]+i] = J[i];
1374:   }
1375:   return(0);
1376: }

1378: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1379: {
1381:   PetscInt       j;
1382:   PetscScalar    val=(PetscScalar)ncols;

1385:   if (!ghost) {
1386:     for (j=0; j<nrows; j++) {
1387:       VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
1388:     }
1389:   } else {
1390:     for (j=0; j<nrows; j++) {
1391:       VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
1392:     }
1393:   }
1394:   return(0);
1395: }

1397: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1398: {
1400:   PetscInt       j,ncols_u;
1401:   PetscScalar    val;

1404:   if (!ghost) {
1405:     for (j=0; j<nrows; j++) {
1406:       MatGetRow(Ju,j,&ncols_u,NULL,NULL);
1407:       val = (PetscScalar)ncols_u;
1408:       VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
1409:       MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
1410:     }
1411:   } else {
1412:     for (j=0; j<nrows; j++) {
1413:       MatGetRow(Ju,j,&ncols_u,NULL,NULL);
1414:       val = (PetscScalar)ncols_u;
1415:       VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
1416:       MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
1417:     }
1418:   }
1419:   return(0);
1420: }

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

1427:   if (Ju) {
1428:     MatSetPreallocationUserblock_private(Ju,nrows,rows,ncols,ghost,vdnz,vonz);
1429:   } else {
1430:     MatSetPreallocationDenseblock_private(nrows,rows,ncols,ghost,vdnz,vonz);
1431:   }
1432:   return(0);
1433: }

1435: PETSC_STATIC_INLINE PetscErrorCode MatSetDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1436: {
1438:   PetscInt       j,*cols;
1439:   PetscScalar    *zeros;

1442:   PetscCalloc2(ncols,&cols,nrows*ncols,&zeros);
1443:   for (j=0; j<ncols; j++) cols[j] = j+ cstart;
1444:   MatSetValues(*J,nrows,rows,ncols,cols,zeros,INSERT_VALUES);
1445:   PetscFree2(cols,zeros);
1446:   return(0);
1447: }

1449: PETSC_STATIC_INLINE PetscErrorCode MatSetUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1450: {
1452:   PetscInt       j,M,N,row,col,ncols_u;
1453:   const PetscInt *cols;
1454:   PetscScalar    zero=0.0;

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

1460:   for (row=0; row<nrows; row++) {
1461:     MatGetRow(Ju,row,&ncols_u,&cols,NULL);
1462:     for (j=0; j<ncols_u; j++) {
1463:       col = cols[j] + cstart;
1464:       MatSetValues(*J,1,&rows[row],1,&col,&zero,INSERT_VALUES);
1465:     }
1466:     MatRestoreRow(Ju,row,&ncols_u,&cols,NULL);
1467:   }
1468:   return(0);
1469: }

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

1476:   if (Ju) {
1477:     MatSetUserblock_private(Ju,nrows,rows,ncols,cstart,J);
1478:   } else {
1479:     MatSetDenseblock_private(nrows,rows,ncols,cstart,J);
1480:   }
1481:   return(0);
1482: }

1484: /* 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.
1485: */
1486: PetscErrorCode CreateSubGlobalToLocalMapping_private(PetscSection globalsec, PetscSection localsec, ISLocalToGlobalMapping *ltog)
1487: {
1489:   PetscInt       i, size, dof;
1490:   PetscInt       *glob2loc;

1493:   PetscSectionGetStorageSize(localsec,&size);
1494:   PetscMalloc1(size,&glob2loc);

1496:   for (i = 0; i < size; i++) {
1497:     PetscSectionGetOffset(globalsec,i,&dof);
1498:     dof = (dof >= 0) ? dof : -(dof + 1);
1499:     glob2loc[i] = dof;
1500:   }

1502:   ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,size,glob2loc,PETSC_OWN_POINTER,ltog);
1503: #if 0
1504:   PetscIntView(size,glob2loc,PETSC_VIEWER_STDOUT_WORLD);
1505: #endif
1506:   return(0);
1507: }

1509:  #include <petsc/private/matimpl.h>
1510: PetscErrorCode DMCreateMatrix_Network(DM dm,Mat *J)
1511: {
1513:   PetscMPIInt    rank, size;
1514:   DM_Network     *network = (DM_Network*)dm->data;
1515:   PetscInt       eStart,eEnd,vStart,vEnd,rstart,nrows,*rows,localSize;
1516:   PetscInt       cstart,ncols,j,e,v;
1517:   PetscBool      ghost,ghost_vc,ghost2,isNest;
1518:   Mat            Juser;
1519:   PetscSection   sectionGlobal;
1520:   PetscInt       nedges,*vptr=NULL,vc,*rows_v; /* suppress maybe-uninitialized warning */
1521:   const PetscInt *edges,*cone;
1522:   MPI_Comm       comm;
1523:   MatType        mtype;
1524:   Vec            vd_nz,vo_nz;
1525:   PetscInt       *dnnz,*onnz;
1526:   PetscScalar    *vdnz,*vonz;

1529:   mtype = dm->mattype;
1530:   PetscStrcmp(mtype, MATNEST, &isNest);

1532:   if (isNest) {
1533:     /* DMCreateMatrix_Network_Nest(); */
1534:     PetscInt   eDof, vDof;
1535:     Mat        j11, j12, j21, j22, bA[2][2];
1536:     ISLocalToGlobalMapping eISMap, vISMap;

1538:     PetscObjectGetComm((PetscObject)dm,&comm);
1539:     MPI_Comm_rank(comm,&rank);
1540:     MPI_Comm_size(comm,&size);

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

1545:     MatCreate(comm, &j11);
1546:     MatSetSizes(j11, eDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
1547:     MatSetType(j11, MATMPIAIJ);

1549:     MatCreate(comm, &j12);
1550:     MatSetSizes(j12, eDof, vDof, PETSC_DETERMINE ,PETSC_DETERMINE);
1551:     MatSetType(j12, MATMPIAIJ);

1553:     MatCreate(comm, &j21);
1554:     MatSetSizes(j21, vDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
1555:     MatSetType(j21, MATMPIAIJ);

1557:     MatCreate(comm, &j22);
1558:     MatSetSizes(j22, vDof, vDof, PETSC_DETERMINE, PETSC_DETERMINE);
1559:     MatSetType(j22, MATMPIAIJ);

1561:     bA[0][0] = j11;
1562:     bA[0][1] = j12;
1563:     bA[1][0] = j21;
1564:     bA[1][1] = j22;

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

1569:     MatSetLocalToGlobalMapping(j11,eISMap,eISMap);
1570:     MatSetLocalToGlobalMapping(j12,eISMap,vISMap);
1571:     MatSetLocalToGlobalMapping(j21,vISMap,eISMap);
1572:     MatSetLocalToGlobalMapping(j22,vISMap,vISMap);

1574:     MatSetUp(j11);
1575:     MatSetUp(j12);
1576:     MatSetUp(j21);
1577:     MatSetUp(j22);

1579:     MatCreateNest(comm,2,NULL,2,NULL,&bA[0][0],J);
1580:     MatSetUp(*J);
1581:     MatNestSetVecType(*J,VECNEST);
1582:     MatDestroy(&j11);
1583:     MatDestroy(&j12);
1584:     MatDestroy(&j21);
1585:     MatDestroy(&j22);

1587:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1588:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);

1590:     /* Free structures */
1591:     ISLocalToGlobalMappingDestroy(&eISMap);
1592:     ISLocalToGlobalMappingDestroy(&vISMap);

1594:     return(0);
1595:   } else if (!network->userEdgeJacobian && !network->userVertexJacobian) {
1596:     /* user does not provide Jacobian blocks */
1597:     DMCreateMatrix(network->plex,J);
1598:     MatSetDM(*J,dm);
1599:     return(0);
1600:   }

1602:   MatCreate(PetscObjectComm((PetscObject)dm),J);
1603:   DMGetDefaultGlobalSection(network->plex,&sectionGlobal);
1604:   PetscSectionGetConstrainedStorageSize(sectionGlobal,&localSize);
1605:   MatSetSizes(*J,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);

1607:   MatSetType(*J,MATAIJ);
1608:   MatSetFromOptions(*J);

1610:   /* (1) Set matrix preallocation */
1611:   /*------------------------------*/
1612:   PetscObjectGetComm((PetscObject)dm,&comm);
1613:   VecCreate(comm,&vd_nz);
1614:   VecSetSizes(vd_nz,localSize,PETSC_DECIDE);
1615:   VecSetFromOptions(vd_nz);
1616:   VecSet(vd_nz,0.0);
1617:   VecDuplicate(vd_nz,&vo_nz);

1619:   /* Set preallocation for edges */
1620:   /*-----------------------------*/
1621:   DMNetworkGetEdgeRange(dm,&eStart,&eEnd);

1623:   PetscMalloc1(localSize,&rows);
1624:   for (e=eStart; e<eEnd; e++) {
1625:     /* Get row indices */
1626:     DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
1627:     DMNetworkGetNumVariables(dm,e,&nrows);
1628:     if (nrows) {
1629:       for (j=0; j<nrows; j++) rows[j] = j + rstart;

1631:       /* Set preallocation for conntected vertices */
1632:       DMNetworkGetConnectedVertices(dm,e,&cone);
1633:       for (v=0; v<2; v++) {
1634:         DMNetworkGetNumVariables(dm,cone[v],&ncols);

1636:         if (network->Je) {
1637:           Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
1638:         } else Juser = NULL;
1639:         DMNetworkIsGhostVertex(dm,cone[v],&ghost);
1640:         MatSetPreallocationblock_private(Juser,nrows,rows,ncols,ghost,vd_nz,vo_nz);
1641:       }

1643:       /* Set preallocation for edge self */
1644:       cstart = rstart;
1645:       if (network->Je) {
1646:         Juser = network->Je[3*e]; /* Jacobian(e,e) */
1647:       } else Juser = NULL;
1648:       MatSetPreallocationblock_private(Juser,nrows,rows,nrows,PETSC_FALSE,vd_nz,vo_nz);
1649:     }
1650:   }

1652:   /* Set preallocation for vertices */
1653:   /*--------------------------------*/
1654:   DMNetworkGetVertexRange(dm,&vStart,&vEnd);
1655:   if (vEnd - vStart) vptr = network->Jvptr;

1657:   for (v=vStart; v<vEnd; v++) {
1658:     /* Get row indices */
1659:     DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
1660:     DMNetworkGetNumVariables(dm,v,&nrows);
1661:     if (!nrows) continue;

1663:     DMNetworkIsGhostVertex(dm,v,&ghost);
1664:     if (ghost) {
1665:       PetscMalloc1(nrows,&rows_v);
1666:     } else {
1667:       rows_v = rows;
1668:     }

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

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

1675:     for (e=0; e<nedges; e++) {
1676:       /* Supporting edges */
1677:       DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
1678:       DMNetworkGetNumVariables(dm,edges[e],&ncols);

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

1685:       /* Connected vertices */
1686:       DMNetworkGetConnectedVertices(dm,edges[e],&cone);
1687:       vc = (v == cone[0]) ? cone[1]:cone[0];
1688:       DMNetworkIsGhostVertex(dm,vc,&ghost_vc);

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

1692:       if (network->Jv) {
1693:         Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
1694:       } else Juser = NULL;
1695:       if (ghost_vc||ghost) {
1696:         ghost2 = PETSC_TRUE;
1697:       } else {
1698:         ghost2 = PETSC_FALSE;
1699:       }
1700:       MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost2,vd_nz,vo_nz);
1701:     }

1703:     /* Set preallocation for vertex self */
1704:     DMNetworkIsGhostVertex(dm,v,&ghost);
1705:     if (!ghost) {
1706:       DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
1707:       if (network->Jv) {
1708:         Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
1709:       } else Juser = NULL;
1710:       MatSetPreallocationblock_private(Juser,nrows,rows_v,nrows,PETSC_FALSE,vd_nz,vo_nz);
1711:     }
1712:     if (ghost) {
1713:       PetscFree(rows_v);
1714:     }
1715:   }

1717:   VecAssemblyBegin(vd_nz);
1718:   VecAssemblyBegin(vo_nz);

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

1722:   VecAssemblyEnd(vd_nz);
1723:   VecAssemblyEnd(vo_nz);

1725:   VecGetArray(vd_nz,&vdnz);
1726:   VecGetArray(vo_nz,&vonz);
1727:   for (j=0; j<localSize; j++) {
1728:     dnnz[j] = (PetscInt)PetscRealPart(vdnz[j]);
1729:     onnz[j] = (PetscInt)PetscRealPart(vonz[j]);
1730:   }
1731:   VecRestoreArray(vd_nz,&vdnz);
1732:   VecRestoreArray(vo_nz,&vonz);
1733:   VecDestroy(&vd_nz);
1734:   VecDestroy(&vo_nz);

1736:   MatSeqAIJSetPreallocation(*J,0,dnnz);
1737:   MatMPIAIJSetPreallocation(*J,0,dnnz,0,onnz);
1738:   MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);

1740:   PetscFree2(dnnz,onnz);

1742:   /* (2) Set matrix entries for edges */
1743:   /*----------------------------------*/
1744:   for (e=eStart; e<eEnd; e++) {
1745:     /* Get row indices */
1746:     DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
1747:     DMNetworkGetNumVariables(dm,e,&nrows);
1748:     if (nrows) {
1749:       for (j=0; j<nrows; j++) rows[j] = j + rstart;

1751:       /* Set matrix entries for conntected vertices */
1752:       DMNetworkGetConnectedVertices(dm,e,&cone);
1753:       for (v=0; v<2; v++) {
1754:         DMNetworkGetVariableGlobalOffset(dm,cone[v],&cstart);
1755:         DMNetworkGetNumVariables(dm,cone[v],&ncols);

1757:         if (network->Je) {
1758:           Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
1759:         } else Juser = NULL;
1760:         MatSetblock_private(Juser,nrows,rows,ncols,cstart,J);
1761:       }

1763:       /* Set matrix entries for edge self */
1764:       cstart = rstart;
1765:       if (network->Je) {
1766:         Juser = network->Je[3*e]; /* Jacobian(e,e) */
1767:       } else Juser = NULL;
1768:       MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);
1769:     }
1770:   }

1772:   /* Set matrix entries for vertices */
1773:   /*---------------------------------*/
1774:   for (v=vStart; v<vEnd; v++) {
1775:     /* Get row indices */
1776:     DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
1777:     DMNetworkGetNumVariables(dm,v,&nrows);
1778:     if (!nrows) continue;

1780:     DMNetworkIsGhostVertex(dm,v,&ghost);
1781:     if (ghost) {
1782:       PetscMalloc1(nrows,&rows_v);
1783:     } else {
1784:       rows_v = rows;
1785:     }
1786:     for (j=0; j<nrows; j++) rows_v[j] = j + rstart;

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

1791:     for (e=0; e<nedges; e++) {
1792:       /* Supporting edges */
1793:       DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
1794:       DMNetworkGetNumVariables(dm,edges[e],&ncols);

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

1801:       /* Connected vertices */
1802:       DMNetworkGetConnectedVertices(dm,edges[e],&cone);
1803:       vc = (v == cone[0]) ? cone[1]:cone[0];

1805:       DMNetworkGetVariableGlobalOffset(dm,vc,&cstart);
1806:       DMNetworkGetNumVariables(dm,vc,&ncols);

1808:       if (network->Jv) {
1809:         Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
1810:       } else Juser = NULL;
1811:       MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);
1812:     }

1814:     /* Set matrix entries for vertex self */
1815:     if (!ghost) {
1816:       DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
1817:       if (network->Jv) {
1818:         Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
1819:       } else Juser = NULL;
1820:       MatSetblock_private(Juser,nrows,rows_v,nrows,cstart,J);
1821:     }
1822:     if (ghost) {
1823:       PetscFree(rows_v);
1824:     }
1825:   }
1826:   PetscFree(rows);

1828:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1829:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);

1831:   MatSetDM(*J,dm);
1832:   return(0);
1833: }

1835: PetscErrorCode DMDestroy_Network(DM dm)
1836: {
1838:   DM_Network     *network = (DM_Network*)dm->data;
1839:   PetscInt       j;

1842:   if (--network->refct > 0) return(0);
1843:   if (network->Je) {
1844:     PetscFree(network->Je);
1845:   }
1846:   if (network->Jv) {
1847:     PetscFree(network->Jvptr);
1848:     PetscFree(network->Jv);
1849:   }

1851:   ISLocalToGlobalMappingDestroy(&network->vertex.mapping);
1852:   PetscSectionDestroy(&network->vertex.DofSection);
1853:   PetscSectionDestroy(&network->vertex.GlobalDofSection);
1854:   if (network->vertex.sf) {
1855:     PetscSFDestroy(&network->vertex.sf);
1856:   }
1857:   /* edge */
1858:   ISLocalToGlobalMappingDestroy(&network->edge.mapping);
1859:   PetscSectionDestroy(&network->edge.DofSection);
1860:   PetscSectionDestroy(&network->edge.GlobalDofSection);
1861:   if (network->edge.sf) {
1862:     PetscSFDestroy(&network->edge.sf);
1863:   }
1864:   DMDestroy(&network->plex);
1865:   network->edges = NULL;
1866:   PetscSectionDestroy(&network->DataSection);
1867:   PetscSectionDestroy(&network->DofSection);

1869:   for(j=0; j < network->nsubnet; j++) {
1870:     PetscFree(network->subnet[j].edges);
1871:     PetscFree(network->subnet[j].vertices);
1872:   }
1873:   PetscFree(network->subnet);
1874:   PetscFree(network->componentdataarray);
1875:   PetscFree(network->cvalue);
1876:   PetscFree(network->header);
1877:   PetscFree(network);
1878:   return(0);
1879: }

1881: PetscErrorCode DMView_Network(DM dm,PetscViewer viewer)
1882: {
1884:   DM_Network     *network = (DM_Network*)dm->data;

1887:   DMView(network->plex,viewer);
1888:   return(0);
1889: }

1891: PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l)
1892: {
1894:   DM_Network     *network = (DM_Network*)dm->data;

1897:   DMGlobalToLocalBegin(network->plex,g,mode,l);
1898:   return(0);
1899: }

1901: PetscErrorCode DMGlobalToLocalEnd_Network(DM dm, Vec g, InsertMode mode, Vec l)
1902: {
1904:   DM_Network     *network = (DM_Network*)dm->data;

1907:   DMGlobalToLocalEnd(network->plex,g,mode,l);
1908:   return(0);
1909: }

1911: PetscErrorCode DMLocalToGlobalBegin_Network(DM dm, Vec l, InsertMode mode, Vec g)
1912: {
1914:   DM_Network     *network = (DM_Network*)dm->data;

1917:   DMLocalToGlobalBegin(network->plex,l,mode,g);
1918:   return(0);
1919: }

1921: PetscErrorCode DMLocalToGlobalEnd_Network(DM dm, Vec l, InsertMode mode, Vec g)
1922: {
1924:   DM_Network     *network = (DM_Network*)dm->data;

1927:   DMLocalToGlobalEnd(network->plex,l,mode,g);
1928:   return(0);
1929: }