Actual source code: network.c

petsc-3.7.7 2017-09-25
Report Typos and Errors
  1: #include <petsc/private/dmnetworkimpl.h>  /*I  "petscdmnetwork.h"  I*/
  2: #include <petscdmplex.h>
  3: #include <petscsf.h>

  7: /*@
  8:   DMNetworkSetSizes - Sets the local and global vertices and edges.

 10:   Collective on DM
 11:   
 12:   Input Parameters:
 13: + dm - the dm object
 14: . nV - number of local vertices
 15: . nE - number of local edges
 16: . NV - number of global vertices (or PETSC_DETERMINE)
 17: - NE - number of global edges (or PETSC_DETERMINE)

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

 22:    You cannot change the sizes once they have been set

 24:    Level: intermediate

 26: .seealso: DMNetworkCreate
 27: @*/
 28: PetscErrorCode DMNetworkSetSizes(DM dm, PetscInt nV, PetscInt nE, PetscInt NV, PetscInt NE)
 29: {
 31:   DM_Network     *network = (DM_Network*) dm->data;
 32:   PetscInt       a[2],b[2];

 38:   if (NV > 0 && nV > NV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local vertex size %D cannot be larger than global vertex size %D",nV,NV);
 39:   if (NE > 0 && nE > NE) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local edge size %D cannot be larger than global edge size %D",nE,NE);
 40:   if ((network->nNodes >= 0 || network->NNodes >= 0) && (network->nNodes != nV || network->NNodes != NV)) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot change/reset vertex sizes to %D local %D global after previously setting them to %D local %D global",nV,NV,network->nNodes,network->NNodes);
 41:   if ((network->nEdges >= 0 || network->NEdges >= 0) && (network->nEdges != nE || network->NEdges != NE)) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot change/reset edge sizes to %D local %D global after previously setting them to %D local %D global",nE,NE,network->nEdges,network->NEdges);
 42:   if (NE < 0 || NV < 0) {
 43:     a[0] = nV; a[1] = nE;
 44:     MPIU_Allreduce(a,b,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
 45:     NV = b[0]; NE = b[1];
 46:   }
 47:   network->nNodes = nV;
 48:   network->NNodes = NV;
 49:   network->nEdges = nE;
 50:   network->NEdges = NE;
 51:   return(0);
 52: }

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

 59:   Logically collective on DM

 61:   Input Parameters:
 62: . edges - list of edges

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

 68:   Level: intermediate

 70: .seealso: DMNetworkCreate, DMNetworkSetSizes
 71: @*/
 72: PetscErrorCode DMNetworkSetEdgeList(DM dm, int edgelist[])
 73: {
 74:   DM_Network *network = (DM_Network*) dm->data;
 75: 
 77:   network->edges = edgelist;
 78:   return(0);
 79: }

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

 86:   Collective on DM

 88:   Input Parameters
 89: . DM - the dmnetwork object

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

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

 97:   Level: intermediate

 99: .seealso: DMNetworkSetSizes, DMNetworkSetEdgeList
100: @*/
101: PetscErrorCode DMNetworkLayoutSetUp(DM dm)
102: {
104:   DM_Network     *network = (DM_Network*) dm->data;
105:   PetscInt       dim = 1; /* One dimensional network */
106:   PetscInt       numCorners=2;
107:   PetscInt       spacedim=2;
108:   double         *vertexcoords=NULL;
109:   PetscInt       i;
110:   PetscInt       ndata;

113:   if (network->nNodes) {
114:     PetscCalloc1(numCorners*network->nNodes,&vertexcoords);
115:   }
116:   DMPlexCreateFromCellList(PetscObjectComm((PetscObject)dm),dim,network->nEdges,network->nNodes,numCorners,PETSC_FALSE,network->edges,spacedim,vertexcoords,&network->plex);
117:   if (network->nNodes) {
118:     PetscFree(vertexcoords);
119:   }
120:   DMPlexGetChart(network->plex,&network->pStart,&network->pEnd);
121:   DMPlexGetHeightStratum(network->plex,0,&network->eStart,&network->eEnd);
122:   DMPlexGetHeightStratum(network->plex,1,&network->vStart,&network->vEnd);

124:   PetscSectionCreate(PetscObjectComm((PetscObject)dm),&network->DataSection);
125:   PetscSectionCreate(PetscObjectComm((PetscObject)dm),&network->DofSection);
126:   PetscSectionSetChart(network->DataSection,network->pStart,network->pEnd);
127:   PetscSectionSetChart(network->DofSection,network->pStart,network->pEnd);

129:   network->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
130:   PetscCalloc1(network->pEnd-network->pStart,&network->header);
131:   for (i = network->pStart; i < network->pEnd; i++) {
132:     network->header[i].ndata = 0;
133:     ndata = network->header[i].ndata;
134:     PetscSectionAddDof(network->DataSection,i,network->dataheadersize);
135:     network->header[i].offset[ndata] = 0;
136:   }
137:   PetscMalloc1(network->pEnd-network->pStart,&network->cvalue);
138:   return(0);
139: }

143: /*@
144:   DMNetworkRegisterComponent - Registers the network component

146:   Logically collective on DM

148:   Input Parameters
149: + dm   - the network object
150: . name - the component name
151: - size - the storage size in bytes for this component data

153:    Output Parameters
154: .   key - an integer key that defines the component

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

159:    Level: intermediate

161: .seealso: DMNetworkLayoutSetUp, DMNetworkCreate
162: @*/
163: PetscErrorCode DMNetworkRegisterComponent(DM dm,const char *name,PetscInt size,PetscInt *key)
164: {
165:   PetscErrorCode        ierr;
166:   DM_Network            *network = (DM_Network*) dm->data;
167:   DMNetworkComponent    *component=&network->component[network->ncomponent];
168:   PetscBool             flg=PETSC_FALSE;
169:   PetscInt              i;


173:   for (i=0; i < network->ncomponent; i++) {
174:     PetscStrcmp(component->name,name,&flg);
175:     if (flg) {
176:       *key = i;
177:       return(0);
178:     }
179:   }
180: 
181:   PetscStrcpy(component->name,name);
182:   component->size = size/sizeof(DMNetworkComponentGenericDataType);
183:   *key = network->ncomponent;
184:   network->ncomponent++;
185:   return(0);
186: }

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

193:   Not Collective

195:   Input Parameters:
196: + dm - The DMNetwork object

198:   Output Paramters:
199: + vStart - The first vertex point
200: - vEnd   - One beyond the last vertex point

202:   Level: intermediate

204: .seealso: DMNetworkGetEdgeRange
205: @*/
206: PetscErrorCode DMNetworkGetVertexRange(DM dm,PetscInt *vStart,PetscInt *vEnd)
207: {
208:   DM_Network     *network = (DM_Network*)dm->data;

211:   if (vStart) *vStart = network->vStart;
212:   if (vEnd) *vEnd = network->vEnd;
213:   return(0);
214: }

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

221:   Not Collective

223:   Input Parameters:
224: + dm - The DMNetwork object

226:   Output Paramters:
227: + eStart - The first edge point
228: - eEnd   - One beyond the last edge point

230:   Level: intermediate

232: .seealso: DMNetworkGetVertexRange
233: @*/
234: PetscErrorCode DMNetworkGetEdgeRange(DM dm,PetscInt *eStart,PetscInt *eEnd)
235: {
236:   DM_Network     *network = (DM_Network*)dm->data;

239:   if (eStart) *eStart = network->eStart;
240:   if (eEnd) *eEnd = network->eEnd;
241:   return(0);
242: }

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

249:   Not Collective

251:   Input Parameters:
252: + dm           - The DMNetwork object
253: . p            - vertex/edge point
254: . componentkey - component key returned while registering the component
255: - compvalue    - pointer to the data structure for the component

257:   Level: intermediate

259: .seealso: DMNetworkGetVertexRange, DMNetworkGetEdgeRange, DMNetworkRegisterComponent
260: @*/
261: PetscErrorCode DMNetworkAddComponent(DM dm, PetscInt p,PetscInt componentkey,void* compvalue)
262: {
263:   DM_Network               *network = (DM_Network*)dm->data;
264:   DMNetworkComponent       *component = &network->component[componentkey];
265:   DMNetworkComponentHeader header = &network->header[p];
266:   DMNetworkComponentValue  cvalue = &network->cvalue[p];
267:   PetscErrorCode           ierr;
268: 
270:   header->size[header->ndata] = component->size;
271:   PetscSectionAddDof(network->DataSection,p,component->size);
272:   header->key[header->ndata] = componentkey;
273:   if (header->ndata != 0) header->offset[header->ndata] = header->offset[header->ndata-1] + header->size[header->ndata-1];

275:   cvalue->data[header->ndata] = (void*)compvalue;
276:   header->ndata++;
277:   return(0);
278: }

282: /*@
283:   DMNetworkGetNumComponents - Get the number of components at a vertex/edge

285:   Not Collective 

287:   Input Parameters:
288: + dm - The DMNetwork object
289: . p  - vertex/edge point

291:   Output Parameters:
292: . numcomponents - Number of components at the vertex/edge

294:   Level: intermediate

296: .seealso: DMNetworkRegisterComponent, DMNetworkAddComponent
297: @*/
298: PetscErrorCode DMNetworkGetNumComponents(DM dm,PetscInt p,PetscInt *numcomponents)
299: {
301:   PetscInt       offset;
302:   DM_Network     *network = (DM_Network*)dm->data;

305:   PetscSectionGetOffset(network->DataSection,p,&offset);
306:   *numcomponents = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata;
307:   return(0);
308: }

312: /*@
313:   DMNetworkGetComponentTypeOffset - Gets the type along with the offset for indexing the 
314:                                     component value from the component data array

316:   Not Collective

318:   Input Parameters:
319: + dm      - The DMNetwork object
320: . p       - vertex/edge point
321: - compnum - component number
322:         
323:   Output Parameters:
324: + compkey - the key obtained when registering the component
325: - offset  - offset into the component data array associated with the vertex/edge point

327:   Notes:
328:   Typical usage:

330:   DMNetworkGetComponentDataArray(dm, &arr);
331:   DMNetworkGetVertex/EdgeRange(dm,&Start,&End);
332:   Loop over vertices or edges
333:     DMNetworkGetNumComponents(dm,v,&numcomps);
334:     Loop over numcomps
335:       DMNetworkGetComponentTypeOffset(dm,v,compnum,&key,&offset);
336:       compdata = (UserCompDataType)(arr+offset);
337:   
338:   Level: intermediate

340: .seealso: DMNetworkGetNumComponents, DMNetworkGetComponentDataArray, 
341: @*/
342: PetscErrorCode DMNetworkGetComponentTypeOffset(DM dm,PetscInt p, PetscInt compnum, PetscInt *compkey, PetscInt *offset)
343: {
344:   PetscErrorCode           ierr;
345:   PetscInt                 offsetp;
346:   DMNetworkComponentHeader header;
347:   DM_Network               *network = (DM_Network*)dm->data;

350:   PetscSectionGetOffset(network->DataSection,p,&offsetp);
351:   header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
352:   if (compkey) *compkey = header->key[compnum];
353:   if (offset) *offset  = offsetp+network->dataheadersize+header->offset[compnum];
354:   return(0);
355: }

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

362:   Not Collective

364:   Input Parameters:
365: + dm     - The DMNetwork object
366: - p      - the edge/vertex point

368:   Output Parameters:
369: . offset - the offset

371:   Level: intermediate

373: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
374: @*/
375: PetscErrorCode DMNetworkGetVariableOffset(DM dm,PetscInt p,PetscInt *offset)
376: {
378:   DM_Network     *network = (DM_Network*)dm->data;

381:   PetscSectionGetOffset(network->DofSection,p,offset);
382:   return(0);
383: }

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

390:   Not Collective

392:   Input Parameters:
393: + dm      - The DMNetwork object
394: - p       - the edge/vertex point

396:   Output Parameters:
397: . offsetg - the offset

399:   Level: intermediate

401: .seealso: DMNetworkGetVariableOffset, DMGetLocalVector
402: @*/
403: PetscErrorCode DMNetworkGetVariableGlobalOffset(DM dm,PetscInt p,PetscInt *offsetg)
404: {
406:   DM_Network     *network = (DM_Network*)dm->data;

409:   PetscSectionGetOffset(network->GlobalDofSection,p,offsetg);
410:   if (*offsetg < 0) *offsetg = -(*offsetg + 1); /* Convert to actual global offset for ghost node */
411:   return(0);
412: }

416: /*@ 
417:   DMNetworkAddNumVariables - Add number of variables associated with a given point.

419:   Not Collective

421:   Input Parameters:
422: + dm   - The DMNetworkObject
423: . p    - the vertex/edge point
424: - nvar - number of additional variables

426:   Level: intermediate

428: .seealso: DMNetworkSetNumVariables
429: @*/
430: PetscErrorCode DMNetworkAddNumVariables(DM dm,PetscInt p,PetscInt nvar)
431: {
433:   DM_Network     *network = (DM_Network*)dm->data;

436:   PetscSectionAddDof(network->DofSection,p,nvar);
437:   return(0);
438: }

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

445:   Not Collective

447:   Input Parameters:
448: + dm   - The DMNetworkObject
449: - p    - the vertex/edge point

451:   Output Parameters:
452: . nvar - number of variables

454:   Level: intermediate

456: .seealso: DMNetworkAddNumVariables, DMNetworkSddNumVariables
457: @*/
458: PetscErrorCode DMNetworkGetNumVariables(DM dm,PetscInt p,PetscInt *nvar)
459: {
461:   DM_Network     *network = (DM_Network*)dm->data;

464:   PetscSectionGetDof(network->DofSection,p,nvar);
465:   return(0);
466: }

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

473:   Not Collective

475:   Input Parameters:
476: + dm   - The DMNetworkObject
477: . p    - the vertex/edge point
478: - nvar - number of variables

480:   Level: intermediate

482: .seealso: DMNetworkAddNumVariables
483: @*/
484: PetscErrorCode DMNetworkSetNumVariables(DM dm,PetscInt p,PetscInt nvar)
485: {
487:   DM_Network     *network = (DM_Network*)dm->data;

490:   PetscSectionSetDof(network->DofSection,p,nvar);
491:   return(0);
492: }

494: /* Sets up the array that holds the data for all components and its associated section. This
495:    function is called during DMSetUp() */
498: PetscErrorCode DMNetworkComponentSetUp(DM dm)
499: {
500:   PetscErrorCode              ierr;
501:   DM_Network     *network = (DM_Network*)dm->data;
502:   PetscInt                    arr_size;
503:   PetscInt                    p,offset,offsetp;
504:   DMNetworkComponentHeader header;
505:   DMNetworkComponentValue  cvalue;
506:   DMNetworkComponentGenericDataType      *componentdataarray;
507:   PetscInt ncomp, i;

510:   PetscSectionSetUp(network->DataSection);
511:   PetscSectionGetStorageSize(network->DataSection,&arr_size);
512:   PetscMalloc1(arr_size,&network->componentdataarray);
513:   componentdataarray = network->componentdataarray;
514:   for (p = network->pStart; p < network->pEnd; p++) {
515:     PetscSectionGetOffset(network->DataSection,p,&offsetp);
516:     /* Copy header */
517:     header = &network->header[p];
518:     PetscMemcpy(componentdataarray+offsetp,header,network->dataheadersize*sizeof(DMNetworkComponentGenericDataType));
519:     /* Copy data */
520:     cvalue = &network->cvalue[p];
521:     ncomp = header->ndata;
522:     for (i = 0; i < ncomp; i++) {
523:       offset = offsetp + network->dataheadersize + header->offset[i];
524:       PetscMemcpy(componentdataarray+offset,cvalue->data[i],header->size[i]*sizeof(DMNetworkComponentGenericDataType));
525:     }
526:   }
527:   return(0);
528: }

530: /* Sets up the section for dofs. This routine is called during DMSetUp() */
533: PetscErrorCode DMNetworkVariablesSetUp(DM dm)
534: {
536:   DM_Network     *network = (DM_Network*)dm->data;

539:   PetscSectionSetUp(network->DofSection);
540:   return(0);
541: }

545: /*@C
546:   DMNetworkGetComponentDataArray - Returns the component data array

548:   Not Collective

550:   Input Parameters:
551: . dm - The DMNetwork Object

553:   Output Parameters:
554: . componentdataarray - array that holds data for all components

556:   Level: intermediate

558: .seealso: DMNetworkGetComponentTypeOffset, DMNetworkGetNumComponents
559: @*/
560: PetscErrorCode DMNetworkGetComponentDataArray(DM dm,DMNetworkComponentGenericDataType **componentdataarray)
561: {
562:   DM_Network     *network = (DM_Network*)dm->data;

565:   *componentdataarray = network->componentdataarray;
566:   return(0);
567: }

571: /*@
572:   DMNetworkDistribute - Distributes the network and moves associated component data.

574:   Collective

576:   Input Parameter:
577: + oldDM - the original DMNetwork object
578: - overlap - The overlap of partitions, 0 is the default

580:   Output Parameter:
581: . distDM - the distributed DMNetwork object

583:   Notes:
584:   This routine should be called only when using multiple processors.

586:   Distributes the network with <overlap>-overlapping partitioning of the edges.

588:   Level: intermediate

590: .seealso: DMNetworkCreate
591: @*/
592: PetscErrorCode DMNetworkDistribute(DM oldDM, PetscInt overlap,DM *distDM)
593: {
595:   DM_Network     *oldDMnetwork = (DM_Network*)oldDM->data;
596:   PetscSF        pointsf;
597:   DM             newDM;
598:   DM_Network     *newDMnetwork;

601:   DMNetworkCreate(PetscObjectComm((PetscObject)oldDM),&newDM);
602:   newDMnetwork = (DM_Network*)newDM->data;
603:   newDMnetwork->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
604:   /* Distribute plex dm and dof section */
605:   DMPlexDistribute(oldDMnetwork->plex,overlap,&pointsf,&newDMnetwork->plex);
606:   /* Distribute dof section */
607:   PetscSectionCreate(PetscObjectComm((PetscObject)oldDM),&newDMnetwork->DofSection);
608:   PetscSFDistributeSection(pointsf,oldDMnetwork->DofSection,NULL,newDMnetwork->DofSection);
609:   PetscSectionCreate(PetscObjectComm((PetscObject)oldDM),&newDMnetwork->DataSection);
610:   /* Distribute data and associated section */
611:   DMPlexDistributeData(newDMnetwork->plex,pointsf,oldDMnetwork->DataSection,MPIU_INT,(void*)oldDMnetwork->componentdataarray,newDMnetwork->DataSection,(void**)&newDMnetwork->componentdataarray);
612:   /* Destroy point SF */
613:   PetscSFDestroy(&pointsf);
614: 
615:   PetscSectionGetChart(newDMnetwork->DataSection,&newDMnetwork->pStart,&newDMnetwork->pEnd);
616:   DMPlexGetHeightStratum(newDMnetwork->plex,0, &newDMnetwork->eStart,&newDMnetwork->eEnd);
617:   DMPlexGetHeightStratum(newDMnetwork->plex,1,&newDMnetwork->vStart,&newDMnetwork->vEnd);
618:   newDMnetwork->nEdges = newDMnetwork->eEnd - newDMnetwork->eStart;
619:   newDMnetwork->nNodes = newDMnetwork->vEnd - newDMnetwork->vStart;
620:   newDMnetwork->NNodes = oldDMnetwork->NNodes;
621:   newDMnetwork->NEdges = oldDMnetwork->NEdges;
622:   /* Set Dof section as the default section for dm */
623:   DMSetDefaultSection(newDMnetwork->plex,newDMnetwork->DofSection);
624:   DMGetDefaultGlobalSection(newDMnetwork->plex,&newDMnetwork->GlobalDofSection);

626:   *distDM = newDM;
627:   return(0);
628: }

632: /*@C
633:   DMNetworkGetSupportingEdges - Return the supporting edges for this vertex point

635:   Not Collective

637:   Input Parameters:
638: + dm - The DMNetwork object
639: - p  - the vertex point

641:   Output Paramters:
642: + nedges - number of edges connected to this vertex point
643: - edges  - List of edge points

645:   Level: intermediate

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

651: .seealso: DMNetworkCreate, DMNetworkGetConnectedNodes
652: @*/
653: PetscErrorCode DMNetworkGetSupportingEdges(DM dm,PetscInt vertex,PetscInt *nedges,const PetscInt *edges[])
654: {
656:   DM_Network     *network = (DM_Network*)dm->data;

659:   DMPlexGetSupportSize(network->plex,vertex,nedges);
660:   DMPlexGetSupport(network->plex,vertex,edges);
661:   return(0);
662: }

666: /*@C
667:   DMNetworkGetConnectedNodes - Return the connected vertices for this edge point

669:   Not Collective

671:   Input Parameters:
672: + dm - The DMNetwork object
673: - p  - the edge point

675:   Output Paramters:
676: . vertices  - vertices connected to this edge

678:   Level: intermediate

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

684: .seealso: DMNetworkCreate, DMNetworkGetSupportingEdges
685: @*/
686: PetscErrorCode DMNetworkGetConnectedNodes(DM dm,PetscInt edge,const PetscInt *vertices[])
687: {
689:   DM_Network     *network = (DM_Network*)dm->data;

692:   DMPlexGetCone(network->plex,edge,vertices);
693:   return(0);
694: }

698: /*@
699:   DMNetworkIsGhostVertex - Returns TRUE if the vertex is a ghost vertex

701:   Not Collective

703:   Input Parameters:
704: + dm - The DMNetwork object
705: . p  - the vertex point

707:   Output Parameter:
708: . isghost - TRUE if the vertex is a ghost point 

710:   Level: intermediate

712: .seealso: DMNetworkCreate, DMNetworkGetConnectedNodes, DMNetworkGetVertexRange
713: @*/
714: PetscErrorCode DMNetworkIsGhostVertex(DM dm,PetscInt p,PetscBool *isghost)
715: {
717:   DM_Network     *network = (DM_Network*)dm->data;
718:   PetscInt       offsetg;
719:   PetscSection   sectiong;

722:   *isghost = PETSC_FALSE;
723:   DMGetDefaultGlobalSection(network->plex,&sectiong);
724:   PetscSectionGetOffset(sectiong,p,&offsetg);
725:   if (offsetg < 0) *isghost = PETSC_TRUE;
726:   return(0);
727: }

731: PetscErrorCode DMSetUp_Network(DM dm)
732: {
734:   DM_Network     *network=(DM_Network*)dm->data;

737:   DMNetworkComponentSetUp(dm);
738:   DMNetworkVariablesSetUp(dm);

740:   DMSetDefaultSection(network->plex,network->DofSection);
741:   DMGetDefaultGlobalSection(network->plex,&network->GlobalDofSection);
742:   return(0);
743: }

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

751:     Collective

753:     Input Parameters:
754: +   dm - The DMNetwork object
755: .   eflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for edges
756: -   vflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for vertices

758:     Level: intermediate

760: @*/
761: PetscErrorCode DMNetworkHasJacobian(DM dm,PetscBool eflg,PetscBool vflg)
762: {
763:   DM_Network     *network=(DM_Network*)dm->data;

766:   network->userEdgeJacobian   = eflg;
767:   network->userVertexJacobian = vflg;
768:   return(0);
769: }

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

776:     Not Collective

778:     Input Parameters:
779: +   dm - The DMNetwork object
780: .   p  - the edge point
781: -   J - array (size = 3) of Jacobian submatrices for this edge point:
782:         J[0]: this edge
783:         J[1] and J[2]: connected vertices, obtained by calling DMNetworkGetConnectedNodes()

785:     Level: intermediate

787: .seealso: DMNetworkVertexSetMatrix
788: @*/
789: PetscErrorCode DMNetworkEdgeSetMatrix(DM dm,PetscInt p,Mat J[])
790: {
792:   DM_Network     *network=(DM_Network*)dm->data;

795:   if (!network->userEdgeJacobian) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkEdgeSetMatrix");
796:   if (!network->Je) {
797:     PetscCalloc1(3*network->nEdges,&network->Je);
798:   }
799:   network->Je[3*p]   = J[0];
800:   network->Je[3*p+1] = J[1];
801:   network->Je[3*p+2] = J[2];
802:   return(0);
803: }

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

810:     Not Collective

812:     Input Parameters:
813: +   dm - The DMNetwork object
814: .   p  - the vertex point
815: -   J - array of Jacobian (size = 2*(num of supporting edges) + 1) submatrices for this vertex point: 
816:         J[0]:       this vertex
817:         J[1+2*i]:   i-th supporting edge
818:         J[1+2*i+1]: i-th connected vertex

820:     Level: intermediate

822: .seealso: DMNetworkEdgeSetMatrix
823: @*/
824: PetscErrorCode DMNetworkVertexSetMatrix(DM dm,PetscInt p,Mat J[])
825: {
827:   DM_Network     *network=(DM_Network*)dm->data;
828:   PetscInt       i,*vptr,nedges,vStart,vEnd;
829:   const PetscInt *edges;

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

834:   DMNetworkGetVertexRange(dm,&vStart,&vEnd);

836:   if (!network->Jv) {
837:     PetscInt nNodes = network->nNodes,nedges_total;
838:     /* count nvertex_total */
839:     nedges_total = 0;
840:     DMNetworkGetVertexRange(dm,&vStart,&vEnd);
841:     PetscMalloc1(nNodes+1,&vptr);
842: 
843:     vptr[0] = 0;
844:     for (i=0; i<nNodes; i++) {
845:       DMNetworkGetSupportingEdges(dm,i+vStart,&nedges,&edges);
846:       nedges_total += nedges;
847:       vptr[i+1] = vptr[i] + 2*nedges + 1;
848:     }

850:     PetscCalloc1(2*nedges_total+nNodes,&network->Jv);
851:     network->Jvptr = vptr;
852:   }

854:   vptr = network->Jvptr;
855:   network->Jv[vptr[p-vStart]] = J[0]; /* Set Jacobian for this vertex */

857:   /* Set Jacobian for each supporting edge and connected vertex */
858:   DMNetworkGetSupportingEdges(dm,p,&nedges,&edges);
859:   for (i=1; i<=2*nedges; i++) network->Jv[vptr[p-vStart]+i] = J[i];
860:   return(0);
861: }

865: PetscErrorCode MatSetDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
866: {
868:   PetscInt       j,*cols;
869:   PetscScalar    *zeros;

872:   PetscCalloc2(ncols,&cols,nrows*ncols,&zeros);
873:   for (j=0; j<ncols; j++) cols[j] = j+ cstart;
874:   MatSetValues(*J,nrows,rows,ncols,cols,zeros,INSERT_VALUES);
875:   PetscFree2(cols,zeros);
876:   return(0);
877: }

881: PetscErrorCode MatSetUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
882: {
884:   PetscInt       j,M,N,row,col,ncols_u;
885:   const PetscInt *cols;
886:   PetscScalar    zero=0.0;

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

892:   for (row=0; row<nrows; row++) {
893:     MatGetRow(Ju,row,&ncols_u,&cols,NULL);
894:     for (j=0; j<ncols_u; j++) {
895:       col = cols[j] + cstart;
896:       MatSetValues(*J,1,&rows[row],1,&col,&zero,INSERT_VALUES);
897:     }
898:     MatRestoreRow(Ju,row,&ncols_u,&cols,NULL);
899:   }
900:   return(0);
901: }

905: PetscErrorCode MatSetblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
906: {

910:   if (Ju) {
911:     MatSetUserblock_private(Ju,nrows,rows,ncols,cstart,J);
912:   } else {
913:     MatSetDenseblock_private(nrows,rows,ncols,cstart,J);
914:   }
915:   return(0);
916: }

920: PetscErrorCode DMCreateMatrix_Network(DM dm,Mat *J)
921: {
923:   DM_Network     *network = (DM_Network*) dm->data;
924:   PetscInt       eStart,eEnd,vStart,vEnd,rstart,nrows,*rows,localSize;
925:   PetscInt       cstart,ncols,j,e,v,*dnz,*onz,*dnzu,*onzu;
926:   PetscBool      ghost;
927:   Mat            Juser;
928:   PetscSection   sectionGlobal;
929:   PetscInt       nedges,*vptr=NULL,vc; /* suppress maybe-uninitialized warning */
930:   const PetscInt *edges,*cone;

933:   if (!network->userEdgeJacobian && !network->userVertexJacobian) {
934:     /* user does not provide Jacobian blocks */
935:     DMCreateMatrix(network->plex,J);
936:     MatSetDM(*J,dm);
937:     return(0);
938:   }
939: 
940:   MatCreate(PetscObjectComm((PetscObject)dm),J);
941:   DMGetDefaultGlobalSection(network->plex,&sectionGlobal);
942:   PetscSectionGetConstrainedStorageSize(sectionGlobal,&localSize);
943:   MatSetSizes(*J,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);

945:   MatSetType(*J,MATAIJ);
946:   MatSetFromOptions(*J);

948:   /* Preallocation - submatrix for an element (edge/vertex) is allocated as a dense block, 
949:    see DMCreateMatrix_Plex() */
950:   PetscCalloc4(localSize,&dnz,localSize,&onz,localSize,&dnzu,localSize,&onzu);
951:   DMPlexPreallocateOperator(network->plex,1,dnz,onz,dnzu,onzu,*J,PETSC_FALSE);
952:   PetscFree4(dnz,onz,dnzu,onzu);

954:   /* Set matrix entries for edges */
955:   /*------------------------------*/
956:   DMNetworkGetEdgeRange(dm,&eStart,&eEnd);

958:   for (e=eStart; e<eEnd; e++) {
959:     /* Get row indices */
960:     DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
961:     DMNetworkGetNumVariables(dm,e,&nrows);
962:     if (nrows) {
963:       if (!network->Je) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Uer must provide Je");
964:       PetscMalloc1(nrows,&rows);
965:       for (j=0; j<nrows; j++) rows[j] = j + rstart;

967:       /* Set matrix entries for conntected vertices */
968:       DMNetworkGetConnectedNodes(dm,e,&cone);
969:       for (v=0; v<2; v++) {
970:         DMNetworkGetVariableGlobalOffset(dm,cone[v],&cstart);
971:         DMNetworkGetNumVariables(dm,cone[v],&ncols);

973:         Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
974:         MatSetblock_private(Juser,nrows,rows,ncols,cstart,J);
975:       }

977:       /* Set matrix entries for edge self */
978:       cstart = rstart;
979:       Juser = network->Je[3*e]; /* Jacobian(e,e) */
980:       MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);

982:       PetscFree(rows);
983:     }
984:   }

986:   /* Set matrix entries for vertices */
987:   /*---------------------------------*/
988:   DMNetworkGetVertexRange(dm,&vStart,&vEnd);
989:   if (vEnd - vStart) {
990:     if (!network->Jv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Uer must provide Jv");
991:     vptr = network->Jvptr;
992:     if (!vptr) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Uer must provide vptr");
993:   }
994: 
995:   for (v=vStart; v<vEnd; v++) {
996:     /* Get row indices */
997:     DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
998:     DMNetworkGetNumVariables(dm,v,&nrows);
999:     if (!nrows) continue;

1001:     PetscMalloc1(nrows,&rows);
1002:     for (j=0; j<nrows; j++) rows[j] = j + rstart;

1004:     /* Get supporting edges and connected vertices */
1005:     DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
1006: 
1007:     for (e=0; e<nedges; e++) {
1008:       /* Supporting edges */
1009:       DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
1010:       DMNetworkGetNumVariables(dm,edges[e],&ncols);

1012:       Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */
1013:       MatSetblock_private(Juser,nrows,rows,ncols,cstart,J);

1015:       /* Connected vertices */
1016:       DMNetworkGetConnectedNodes(dm,edges[e],&cone);
1017:       vc = (v == cone[0]) ? cone[1]:cone[0];
1018: 
1019:       DMNetworkGetVariableGlobalOffset(dm,vc,&cstart);
1020:       DMNetworkGetNumVariables(dm,vc,&ncols);

1022:       Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
1023:       MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);
1024:     }

1026:     /* Set matrix entries for vertex self */
1027:     DMNetworkIsGhostVertex(dm,v,&ghost);
1028:     if (!ghost) {
1029:       DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
1030:       Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
1031:       MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);
1032:     }
1033:     PetscFree(rows);
1034:   }
1035:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1036:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);

1038:   MatSetDM(*J,dm);
1039:   return(0);
1040: }

1044: PetscErrorCode DMDestroy_Network(DM dm)
1045: {
1047:   DM_Network     *network = (DM_Network*) dm->data;

1050:   if (--network->refct > 0) return(0);
1051:   if (network->Je) {
1052:     PetscFree(network->Je);
1053:   }
1054:   if (network->Jv) {
1055:     PetscFree(network->Jvptr);
1056:     PetscFree(network->Jv);
1057:   }
1058:   DMDestroy(&network->plex);
1059:   network->edges = NULL;
1060:   PetscSectionDestroy(&network->DataSection);
1061:   PetscSectionDestroy(&network->DofSection);
1062: 
1063:   PetscFree(network->componentdataarray);
1064:   PetscFree(network->cvalue);
1065:   PetscFree(network->header);
1066:   PetscFree(network);
1067:   return(0);
1068: }

1072: PetscErrorCode DMView_Network(DM dm, PetscViewer viewer)
1073: {
1075:   DM_Network     *network = (DM_Network*) dm->data;

1078:   DMView(network->plex,viewer);
1079:   return(0);
1080: }

1084: PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l)
1085: {
1087:   DM_Network     *network = (DM_Network*) dm->data;

1090:   DMGlobalToLocalBegin(network->plex,g,mode,l);
1091:   return(0);
1092: }

1096: PetscErrorCode DMGlobalToLocalEnd_Network(DM dm, Vec g, InsertMode mode, Vec l)
1097: {
1099:   DM_Network     *network = (DM_Network*) dm->data;

1102:   DMGlobalToLocalEnd(network->plex,g,mode,l);
1103:   return(0);
1104: }

1108: PetscErrorCode DMLocalToGlobalBegin_Network(DM dm, Vec l, InsertMode mode, Vec g)
1109: {
1111:   DM_Network     *network = (DM_Network*) dm->data;

1114:   DMLocalToGlobalBegin(network->plex,l,mode,g);
1115:   return(0);
1116: }

1120: PetscErrorCode DMLocalToGlobalEnd_Network(DM dm, Vec l, InsertMode mode, Vec g)
1121: {
1123:   DM_Network     *network = (DM_Network*) dm->data;

1126:   DMLocalToGlobalEnd(network->plex,l,mode,g);
1127:   return(0);
1128: }