Actual source code: network.c

petsc-3.8.2 2017-11-09
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 local and global vertices and edges.

 30:   Collective on DM

 32:   Input Parameters:
 33: + dm - the dm object
 34: . nV - number of local vertices
 35: . nE - number of local edges
 36: . NV - number of global vertices (or PETSC_DETERMINE)
 37: - NE - number of global edges (or PETSC_DETERMINE)

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

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

 44:    Level: intermediate

 46: .seealso: DMNetworkCreate()
 47: @*/
 48: PetscErrorCode DMNetworkSetSizes(DM dm, PetscInt nV, PetscInt nE, PetscInt NV, PetscInt NE)
 49: {
 51:   DM_Network     *network = (DM_Network*) dm->data;
 52:   PetscInt       a[2],b[2];

 58:   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);
 59:   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);
 60:   if ((network->nVertices >= 0 || network->NVertices >= 0) && (network->nVertices != nV || network->NVertices != 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->nVertices,network->NVertices);
 61:   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);
 62:   if (NE < 0 || NV < 0) {
 63:     a[0] = nV; a[1] = nE;
 64:     MPIU_Allreduce(a,b,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
 65:     NV = b[0]; NE = b[1];
 66:   }
 67:   network->nVertices = nV;
 68:   network->NVertices = NV;
 69:   network->nEdges = nE;
 70:   network->NEdges = NE;
 71:   return(0);
 72: }

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

 77:   Logically collective on DM

 79:   Input Parameters:
 80: . edges - list of edges

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

 86:   Level: intermediate

 88: .seealso: DMNetworkCreate, DMNetworkSetSizes
 89: @*/
 90: PetscErrorCode DMNetworkSetEdgeList(DM dm, int edgelist[])
 91: {
 92:   DM_Network *network = (DM_Network*) dm->data;

 95:   network->edges = edgelist;
 96:   return(0);
 97: }

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

102:   Collective on DM

104:   Input Parameters
105: . DM - the dmnetwork object

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

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

113:   Level: intermediate

115: .seealso: DMNetworkSetSizes, DMNetworkSetEdgeList
116: @*/
117: PetscErrorCode DMNetworkLayoutSetUp(DM dm)
118: {
120:   DM_Network     *network = (DM_Network*) dm->data;
121:   PetscInt       dim = 1; /* One dimensional network */
122:   PetscInt       numCorners=2;
123:   PetscInt       spacedim=2;
124:   double         *vertexcoords=NULL;
125:   PetscInt       i;
126:   PetscInt       ndata;

129:   if (network->nVertices) {
130:     PetscCalloc1(numCorners*network->nVertices,&vertexcoords);
131:   }
132:   DMPlexCreateFromCellList(PetscObjectComm((PetscObject)dm),dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,network->edges,spacedim,vertexcoords,&network->plex);
133:   if (network->nVertices) {
134:     PetscFree(vertexcoords);
135:   }
136:   DMPlexGetChart(network->plex,&network->pStart,&network->pEnd);
137:   DMPlexGetHeightStratum(network->plex,0,&network->eStart,&network->eEnd);
138:   DMPlexGetHeightStratum(network->plex,1,&network->vStart,&network->vEnd);

140:   PetscSectionCreate(PetscObjectComm((PetscObject)dm),&network->DataSection);
141:   PetscSectionCreate(PetscObjectComm((PetscObject)dm),&network->DofSection);
142:   PetscSectionSetChart(network->DataSection,network->pStart,network->pEnd);
143:   PetscSectionSetChart(network->DofSection,network->pStart,network->pEnd);



147:   network->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
148:   PetscCalloc1(network->pEnd-network->pStart,&network->header);
149:   for (i = network->pStart; i < network->pEnd; i++) {
150:     if (i < network->vStart) {
151:       network->header[i].index = i;
152:     } else {
153:       network->header[i].index = i - network->vStart;
154:     }
155:     network->header[i].ndata = 0;
156:     ndata = network->header[i].ndata;
157:     PetscSectionAddDof(network->DataSection,i,network->dataheadersize);
158:     network->header[i].offset[ndata] = 0;
159:   }
160:   PetscMalloc1(network->pEnd-network->pStart,&network->cvalue);
161:   return(0);
162: }

164: /*@C
165:   DMNetworkRegisterComponent - Registers the network component

167:   Logically collective on DM

169:   Input Parameters
170: + dm   - the network object
171: . name - the component name
172: - size - the storage size in bytes for this component data

174:    Output Parameters
175: .   key - an integer key that defines the component

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

180:    Level: intermediate

182: .seealso: DMNetworkLayoutSetUp, DMNetworkCreate
183: @*/
184: PetscErrorCode DMNetworkRegisterComponent(DM dm,const char *name,PetscInt size,PetscInt *key)
185: {
186:   PetscErrorCode        ierr;
187:   DM_Network            *network = (DM_Network*) dm->data;
188:   DMNetworkComponent    *component=&network->component[network->ncomponent];
189:   PetscBool             flg=PETSC_FALSE;
190:   PetscInt              i;


194:   for (i=0; i < network->ncomponent; i++) {
195:     PetscStrcmp(component->name,name,&flg);
196:     if (flg) {
197:       *key = i;
198:       return(0);
199:     }
200:   }

202:   PetscStrcpy(component->name,name);
203:   component->size = size/sizeof(DMNetworkComponentGenericDataType);
204:   *key = network->ncomponent;
205:   network->ncomponent++;
206:   return(0);
207: }

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

212:   Not Collective

214:   Input Parameters:
215: + dm - The DMNetwork object

217:   Output Paramters:
218: + vStart - The first vertex point
219: - vEnd   - One beyond the last vertex point

221:   Level: intermediate

223: .seealso: DMNetworkGetEdgeRange
224: @*/
225: PetscErrorCode DMNetworkGetVertexRange(DM dm,PetscInt *vStart,PetscInt *vEnd)
226: {
227:   DM_Network     *network = (DM_Network*)dm->data;

230:   if (vStart) *vStart = network->vStart;
231:   if (vEnd) *vEnd = network->vEnd;
232:   return(0);
233: }

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

238:   Not Collective

240:   Input Parameters:
241: + dm - The DMNetwork object

243:   Output Paramters:
244: + eStart - The first edge point
245: - eEnd   - One beyond the last edge point

247:   Level: intermediate

249: .seealso: DMNetworkGetVertexRange
250: @*/
251: PetscErrorCode DMNetworkGetEdgeRange(DM dm,PetscInt *eStart,PetscInt *eEnd)
252: {
253:   DM_Network     *network = (DM_Network*)dm->data;

256:   if (eStart) *eStart = network->eStart;
257:   if (eEnd) *eEnd = network->eEnd;
258:   return(0);
259: }

261: /*@
262:   DMNetworkGetGlobalEdgeIndex - Get the user global numbering for the edge.

264:   Not Collective

266:   Input Parameters:
267: + dm - DMNetwork object
268: - p  - edge point

270:   Output Paramters:
271: . index - user global numbering for the edge

273:   Level: intermediate

275: .seealso: DMNetworkGetGlobalVertexIndex
276: @*/
277: PetscErrorCode DMNetworkGetGlobalEdgeIndex(DM dm,PetscInt p,PetscInt *index)
278: {
279:   PetscErrorCode    ierr;
280:   DM_Network        *network = (DM_Network*)dm->data;
281:   PetscInt          offsetp;
282:   DMNetworkComponentHeader header;

285:   PetscSectionGetOffset(network->DataSection,p,&offsetp);
286:   header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
287:   *index = header->index;
288:   return(0);
289: }

291: /*@
292:   DMNetworkGetGlobalVertexIndex - Get the user global numbering for the vertex.

294:   Not Collective

296:   Input Parameters:
297: + dm - DMNetwork object
298: - p  - vertex point

300:   Output Paramters:
301: . index - user global numbering for the vertex

303:   Level: intermediate

305: .seealso: DMNetworkGetGlobalEdgeIndex
306: @*/
307: PetscErrorCode DMNetworkGetGlobalVertexIndex(DM dm,PetscInt p,PetscInt *index)
308: {
309:   PetscErrorCode    ierr;
310:   DM_Network        *network = (DM_Network*)dm->data;
311:   PetscInt          offsetp;
312:   DMNetworkComponentHeader header;

315:   PetscSectionGetOffset(network->DataSection,p,&offsetp);
316:   header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
317:   *index = header->index;
318:   return(0);
319: }

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

324:   Not Collective

326:   Input Parameters:
327: + dm           - The DMNetwork object
328: . p            - vertex/edge point
329: . componentkey - component key returned while registering the component
330: - compvalue    - pointer to the data structure for the component

332:   Level: intermediate

334: .seealso: DMNetworkGetVertexRange, DMNetworkGetEdgeRange, DMNetworkRegisterComponent
335: @*/
336: PetscErrorCode DMNetworkAddComponent(DM dm, PetscInt p,PetscInt componentkey,void* compvalue)
337: {
338:   DM_Network               *network = (DM_Network*)dm->data;
339:   DMNetworkComponent       *component = &network->component[componentkey];
340:   DMNetworkComponentHeader header = &network->header[p];
341:   DMNetworkComponentValue  cvalue = &network->cvalue[p];
342:   PetscErrorCode           ierr;

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

347:   header->size[header->ndata] = component->size;
348:   PetscSectionAddDof(network->DataSection,p,component->size);
349:   header->key[header->ndata] = componentkey;
350:   if (header->ndata != 0) header->offset[header->ndata] = header->offset[header->ndata-1] + header->size[header->ndata-1];

352:   cvalue->data[header->ndata] = (void*)compvalue;
353:   header->ndata++;
354:   return(0);
355: }

357: /*@
358:   DMNetworkGetNumComponents - Get the number of components at a vertex/edge

360:   Not Collective

362:   Input Parameters:
363: + dm - The DMNetwork object
364: . p  - vertex/edge point

366:   Output Parameters:
367: . numcomponents - Number of components at the vertex/edge

369:   Level: intermediate

371: .seealso: DMNetworkRegisterComponent, DMNetworkAddComponent
372: @*/
373: PetscErrorCode DMNetworkGetNumComponents(DM dm,PetscInt p,PetscInt *numcomponents)
374: {
376:   PetscInt       offset;
377:   DM_Network     *network = (DM_Network*)dm->data;

380:   PetscSectionGetOffset(network->DataSection,p,&offset);
381:   *numcomponents = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata;
382:   return(0);
383: }

385: /*@
386:   DMNetworkGetComponentKeyOffset - Gets the type along with the offset for indexing the
387:                                     component value from the component data array

389:   Not Collective

391:   Input Parameters:
392: + dm      - The DMNetwork object
393: . p       - vertex/edge point
394: - compnum - component number

396:   Output Parameters:
397: + compkey - the key obtained when registering the component
398: - offset  - offset into the component data array associated with the vertex/edge point

400:   Notes:
401:   Typical usage:

403:   DMNetworkGetComponentDataArray(dm, &arr);
404:   DMNetworkGetVertex/EdgeRange(dm,&Start,&End);
405:   Loop over vertices or edges
406:     DMNetworkGetNumComponents(dm,v,&numcomps);
407:     Loop over numcomps
408:       DMNetworkGetComponentKeyOffset(dm,v,compnum,&key,&offset);
409:       compdata = (UserCompDataType)(arr+offset);

411:   Level: intermediate

413: .seealso: DMNetworkGetNumComponents, DMNetworkGetComponentDataArray,
414: @*/
415: PetscErrorCode DMNetworkGetComponentKeyOffset(DM dm,PetscInt p, PetscInt compnum, PetscInt *compkey, PetscInt *offset)
416: {
417:   PetscErrorCode           ierr;
418:   PetscInt                 offsetp;
419:   DMNetworkComponentHeader header;
420:   DM_Network               *network = (DM_Network*)dm->data;

423:   PetscSectionGetOffset(network->DataSection,p,&offsetp);
424:   header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
425:   if (compkey) *compkey = header->key[compnum];
426:   if (offset) *offset  = offsetp+network->dataheadersize+header->offset[compnum];
427:   return(0);
428: }

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

433:   Not Collective

435:   Input Parameters:
436: + dm     - The DMNetwork object
437: - p      - the edge/vertex point

439:   Output Parameters:
440: . offset - the offset

442:   Level: intermediate

444: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
445: @*/
446: PetscErrorCode DMNetworkGetVariableOffset(DM dm,PetscInt p,PetscInt *offset)
447: {
449:   DM_Network     *network = (DM_Network*)dm->data;

452:   PetscSectionGetOffset(network->plex->defaultSection,p,offset);
453:   return(0);
454: }

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

459:   Not Collective

461:   Input Parameters:
462: + dm      - The DMNetwork object
463: - p       - the edge/vertex point

465:   Output Parameters:
466: . offsetg - the offset

468:   Level: intermediate

470: .seealso: DMNetworkGetVariableOffset, DMGetLocalVector
471: @*/
472: PetscErrorCode DMNetworkGetVariableGlobalOffset(DM dm,PetscInt p,PetscInt *offsetg)
473: {
475:   DM_Network     *network = (DM_Network*)dm->data;

478:   PetscSectionGetOffset(network->plex->defaultGlobalSection,p,offsetg);
479:   if (*offsetg < 0) *offsetg = -(*offsetg + 1); /* Convert to actual global offset for ghost vertex */
480:   return(0);
481: }

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

486:   Not Collective

488:   Input Parameters:
489: + dm     - The DMNetwork object
490: - p      - the edge point

492:   Output Parameters:
493: . offset - the offset

495:   Level: intermediate

497: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
498: @*/
499: PetscErrorCode DMNetworkGetEdgeOffset(DM dm,PetscInt p,PetscInt *offset)
500: {
502:   DM_Network     *network = (DM_Network*)dm->data;


506:   PetscSectionGetOffset(network->edge.DofSection,p,offset);
507:   return(0);
508: }

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

513:   Not Collective

515:   Input Parameters:
516: + dm     - The DMNetwork object
517: - p      - the vertex point

519:   Output Parameters:
520: . offset - the offset

522:   Level: intermediate

524: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
525: @*/
526: PetscErrorCode DMNetworkGetVertexOffset(DM dm,PetscInt p,PetscInt *offset)
527: {
529:   DM_Network     *network = (DM_Network*)dm->data;


533:   p -= network->vStart;

535:   PetscSectionGetOffset(network->vertex.DofSection,p,offset);
536:   return(0);
537: }
538: /*@
539:   DMNetworkAddNumVariables - Add number of variables associated with a given point.

541:   Not Collective

543:   Input Parameters:
544: + dm   - The DMNetworkObject
545: . p    - the vertex/edge point
546: - nvar - number of additional variables

548:   Level: intermediate

550: .seealso: DMNetworkSetNumVariables
551: @*/
552: PetscErrorCode DMNetworkAddNumVariables(DM dm,PetscInt p,PetscInt nvar)
553: {
555:   DM_Network     *network = (DM_Network*)dm->data;

558:   PetscSectionAddDof(network->DofSection,p,nvar);
559:   return(0);
560: }

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

565:   Not Collective

567:   Input Parameters:
568: + dm   - The DMNetworkObject
569: - p    - the vertex/edge point

571:   Output Parameters:
572: . nvar - number of variables

574:   Level: intermediate

576: .seealso: DMNetworkAddNumVariables, DMNetworkSddNumVariables
577: @*/
578: PetscErrorCode DMNetworkGetNumVariables(DM dm,PetscInt p,PetscInt *nvar)
579: {
581:   DM_Network     *network = (DM_Network*)dm->data;

584:   PetscSectionGetDof(network->DofSection,p,nvar);
585:   return(0);
586: }

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

591:   Not Collective

593:   Input Parameters:
594: + dm   - The DMNetworkObject
595: . p    - the vertex/edge point
596: - nvar - number of variables

598:   Level: intermediate

600: .seealso: DMNetworkAddNumVariables
601: @*/
602: PetscErrorCode DMNetworkSetNumVariables(DM dm,PetscInt p,PetscInt nvar)
603: {
605:   DM_Network     *network = (DM_Network*)dm->data;

608:   PetscSectionSetDof(network->DofSection,p,nvar);
609:   return(0);
610: }

612: /* Sets up the array that holds the data for all components and its associated section. This
613:    function is called during DMSetUp() */
614: PetscErrorCode DMNetworkComponentSetUp(DM dm)
615: {
616:   PetscErrorCode              ierr;
617:   DM_Network     *network = (DM_Network*)dm->data;
618:   PetscInt                    arr_size;
619:   PetscInt                    p,offset,offsetp;
620:   DMNetworkComponentHeader header;
621:   DMNetworkComponentValue  cvalue;
622:   DMNetworkComponentGenericDataType      *componentdataarray;
623:   PetscInt ncomp, i;

626:   PetscSectionSetUp(network->DataSection);
627:   PetscSectionGetStorageSize(network->DataSection,&arr_size);
628:   PetscMalloc1(arr_size,&network->componentdataarray);
629:   componentdataarray = network->componentdataarray;
630:   for (p = network->pStart; p < network->pEnd; p++) {
631:     PetscSectionGetOffset(network->DataSection,p,&offsetp);
632:     /* Copy header */
633:     header = &network->header[p];
634:     PetscMemcpy(componentdataarray+offsetp,header,network->dataheadersize*sizeof(DMNetworkComponentGenericDataType));
635:     /* Copy data */
636:     cvalue = &network->cvalue[p];
637:     ncomp = header->ndata;
638:     for (i = 0; i < ncomp; i++) {
639:       offset = offsetp + network->dataheadersize + header->offset[i];
640:       PetscMemcpy(componentdataarray+offset,cvalue->data[i],header->size[i]*sizeof(DMNetworkComponentGenericDataType));
641:     }
642:   }
643:   return(0);
644: }

646: /* Sets up the section for dofs. This routine is called during DMSetUp() */
647: PetscErrorCode DMNetworkVariablesSetUp(DM dm)
648: {
650:   DM_Network     *network = (DM_Network*)dm->data;

653:   PetscSectionSetUp(network->DofSection);
654:   return(0);
655: }

657: /*@C
658:   DMNetworkGetComponentDataArray - Returns the component data array

660:   Not Collective

662:   Input Parameters:
663: . dm - The DMNetwork Object

665:   Output Parameters:
666: . componentdataarray - array that holds data for all components

668:   Level: intermediate

670: .seealso: DMNetworkGetComponentKeyOffset, DMNetworkGetNumComponents
671: @*/
672: PetscErrorCode DMNetworkGetComponentDataArray(DM dm,DMNetworkComponentGenericDataType **componentdataarray)
673: {
674:   DM_Network     *network = (DM_Network*)dm->data;

677:   *componentdataarray = network->componentdataarray;
678:   return(0);
679: }

681: /* Get a subsection from a range of points */
682: PetscErrorCode DMNetworkGetSubSection_private(PetscSection master, PetscInt pstart, PetscInt pend,PetscSection *subsection)
683: {
685:   PetscInt       i, nvar;

688:   PetscSectionCreate(PetscObjectComm((PetscObject)master), subsection);
689:   PetscSectionSetChart(*subsection, 0, pend - pstart);
690:   for (i = pstart; i < pend; i++) {
691:     PetscSectionGetDof(master,i,&nvar);
692:     PetscSectionSetDof(*subsection, i - pstart, nvar);
693:   }

695:   PetscSectionSetUp(*subsection);
696:   return(0);
697: }

699: /* Create a submap of points with a GlobalToLocal structure */
700: PetscErrorCode DMNetworkSetSubMap_private(PetscInt pstart, PetscInt pend, ISLocalToGlobalMapping *map)
701: {
703:   PetscInt       i, *subpoints;

706:   /* Create index sets to map from "points" to "subpoints" */
707:   PetscMalloc1(pend - pstart, &subpoints);
708:   for (i = pstart; i < pend; i++) {
709:     subpoints[i - pstart] = i;
710:   }
711:   ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,pend-pstart,subpoints,PETSC_COPY_VALUES,map);
712:   PetscFree(subpoints);
713:   return(0);
714: }

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

719:   Collective

721:   Input Parameters:
722: . dm   - The DMNetworkObject

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

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

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

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

732:   Level: intermediate

734: @*/
735: PetscErrorCode DMNetworkAssembleGraphStructures(DM dm)
736: {
738:   MPI_Comm       comm;
739:   PetscMPIInt    rank, size;
740:   DM_Network     *network = (DM_Network*)dm->data;

743:   PetscObjectGetComm((PetscObject)dm,&comm);
744:   MPI_Comm_rank(comm, &rank);
745:   MPI_Comm_size(comm, &size);

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

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

755:   if (size > 1) {
756:     PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);
757:     PetscSectionCreateGlobalSection(network->vertex.DofSection, network->vertex.sf, PETSC_FALSE, PETSC_FALSE, &network->vertex.GlobalDofSection);
758:   PetscSFGetSubSF(network->plex->sf, network->edge.mapping, &network->edge.sf);
759:   PetscSectionCreateGlobalSection(network->edge.DofSection, network->edge.sf, PETSC_FALSE, PETSC_FALSE, &network->edge.GlobalDofSection);
760:   } else {
761:   /* create structures for vertex */
762:   PetscSectionClone(network->vertex.DofSection,&network->vertex.GlobalDofSection);
763:   /* create structures for edge */
764:   PetscSectionClone(network->edge.DofSection,&network->edge.GlobalDofSection);
765:   }


768:   /* Add viewers */
769:   PetscObjectSetName((PetscObject)network->edge.GlobalDofSection,"Global edge dof section");
770:   PetscObjectSetName((PetscObject)network->vertex.GlobalDofSection,"Global vertex dof section");
771:   PetscSectionViewFromOptions(network->edge.GlobalDofSection, NULL, "-edge_global_section_view");
772:   PetscSectionViewFromOptions(network->vertex.GlobalDofSection, NULL, "-vertex_global_section_view");

774:   return(0);
775: }

777: /*@
778:   DMNetworkDistribute - Distributes the network and moves associated component data.

780:   Collective

782:   Input Parameter:
783: + DM - the DMNetwork object
784: - overlap - The overlap of partitions, 0 is the default

786:   Notes:
787:   Distributes the network with <overlap>-overlapping partitioning of the edges.

789:   Level: intermediate

791: .seealso: DMNetworkCreate
792: @*/
793: PetscErrorCode DMNetworkDistribute(DM *dm,PetscInt overlap)
794: {
795:   MPI_Comm       comm;
797:   PetscMPIInt    size;
798:   DM_Network     *oldDMnetwork = (DM_Network*)((*dm)->data);
799:   DM_Network     *newDMnetwork;
800:   PetscSF        pointsf;
801:   DM             newDM;
802:   PetscPartitioner part;


806:   PetscObjectGetComm((PetscObject)*dm,&comm);
807:   MPI_Comm_size(comm, &size);
808:   if (size == 1) return(0);

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

814:   /* Enable runtime options for petscpartitioner */
815:   DMPlexGetPartitioner(oldDMnetwork->plex,&part);
816:   PetscPartitionerSetFromOptions(part);

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

821:   /* Distribute dof section */
822:   PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DofSection);
823:   PetscSFDistributeSection(pointsf,oldDMnetwork->DofSection,NULL,newDMnetwork->DofSection);
824:   PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DataSection);

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

829:   PetscSectionGetChart(newDMnetwork->DataSection,&newDMnetwork->pStart,&newDMnetwork->pEnd);
830:   DMPlexGetHeightStratum(newDMnetwork->plex,0, &newDMnetwork->eStart,&newDMnetwork->eEnd);
831:   DMPlexGetHeightStratum(newDMnetwork->plex,1,&newDMnetwork->vStart,&newDMnetwork->vEnd);
832:   newDMnetwork->nEdges = newDMnetwork->eEnd - newDMnetwork->eStart;
833:   newDMnetwork->nVertices = newDMnetwork->vEnd - newDMnetwork->vStart;
834:   newDMnetwork->NVertices = oldDMnetwork->NVertices;
835:   newDMnetwork->NEdges = oldDMnetwork->NEdges;

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

841:   /* Destroy point SF */
842:   PetscSFDestroy(&pointsf);

844:   DMDestroy(dm);
845:   *dm  = newDM;
846:   return(0);
847: }

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

852:   Input Parameters:
853: + masterSF - the original SF structure
854: - map      - a ISLocalToGlobal mapping that contains the subset of points

856:   Output Parameters:
857: . subSF    - a subset of the masterSF for the desired subset.
858: */

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

862:   PetscErrorCode        ierr;
863:   PetscInt              nroots, nleaves, *ilocal_sub;
864:   PetscInt              i, *ilocal_map, nroots_sub, nleaves_sub = 0;
865:   PetscInt              *local_points, *remote_points;
866:   PetscSFNode           *iremote_sub;
867:   const PetscInt        *ilocal;
868:   const PetscSFNode     *iremote;

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

873:   /* Look for leaves that pertain to the subset of points. Get the local ordering */
874:   PetscMalloc1(nleaves,&ilocal_map);
875:   ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nleaves,ilocal,NULL,ilocal_map);
876:   for (i = 0; i < nleaves; i++) {
877:     if (ilocal_map[i] != -1) nleaves_sub += 1;
878:   }
879:   /* Re-number ilocal with subset numbering. Need information from roots */
880:   PetscMalloc2(nroots,&local_points,nroots,&remote_points);
881:   for (i = 0; i < nroots; i++) local_points[i] = i;
882:   ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nroots,local_points,NULL,local_points);
883:   PetscSFBcastBegin(mastersf, MPIU_INT, local_points, remote_points);
884:   PetscSFBcastEnd(mastersf, MPIU_INT, local_points, remote_points);
885:   /* Fill up graph using local (that is, local to the subset) numbering. */
886:   PetscMalloc1(nleaves_sub,&ilocal_sub);
887:   PetscMalloc1(nleaves_sub,&iremote_sub);
888:   nleaves_sub = 0;
889:   for (i = 0; i < nleaves; i++) {
890:     if (ilocal_map[i] != -1) {
891:       ilocal_sub[nleaves_sub] = ilocal_map[i];
892:       iremote_sub[nleaves_sub].rank = iremote[i].rank;
893:       iremote_sub[nleaves_sub].index = remote_points[ilocal[i]];
894:       nleaves_sub += 1;
895:     }
896:   }
897:   PetscFree2(local_points,remote_points);
898:   ISLocalToGlobalMappingGetSize(map,&nroots_sub);

900:   /* Create new subSF */
901:   PetscSFCreate(PETSC_COMM_WORLD,subSF);
902:   PetscSFSetFromOptions(*subSF);
903:   PetscSFSetGraph(*subSF,nroots_sub,nleaves_sub,ilocal_sub,PETSC_OWN_POINTER,iremote_sub,PETSC_COPY_VALUES);
904:   PetscFree(ilocal_map);
905:   PetscFree(iremote_sub);
906:   return(0);
907: }

909: /*@C
910:   DMNetworkGetSupportingEdges - Return the supporting edges for this vertex point

912:   Not Collective

914:   Input Parameters:
915: + dm - The DMNetwork object
916: - p  - the vertex point

918:   Output Paramters:
919: + nedges - number of edges connected to this vertex point
920: - edges  - List of edge points

922:   Level: intermediate

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

928: .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices
929: @*/
930: PetscErrorCode DMNetworkGetSupportingEdges(DM dm,PetscInt vertex,PetscInt *nedges,const PetscInt *edges[])
931: {
933:   DM_Network     *network = (DM_Network*)dm->data;

936:   DMPlexGetSupportSize(network->plex,vertex,nedges);
937:   DMPlexGetSupport(network->plex,vertex,edges);
938:   return(0);
939: }

941: /*@C
942:   DMNetworkGetConnectedVertices - Return the connected vertices for this edge point

944:   Not Collective

946:   Input Parameters:
947: + dm - The DMNetwork object
948: - p  - the edge point

950:   Output Paramters:
951: . vertices  - vertices connected to this edge

953:   Level: intermediate

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

959: .seealso: DMNetworkCreate, DMNetworkGetSupportingEdges
960: @*/
961: PetscErrorCode DMNetworkGetConnectedVertices(DM dm,PetscInt edge,const PetscInt *vertices[])
962: {
964:   DM_Network     *network = (DM_Network*)dm->data;

967:   DMPlexGetCone(network->plex,edge,vertices);
968:   return(0);
969: }

971: /*@
972:   DMNetworkIsGhostVertex - Returns TRUE if the vertex is a ghost vertex

974:   Not Collective

976:   Input Parameters:
977: + dm - The DMNetwork object
978: . p  - the vertex point

980:   Output Parameter:
981: . isghost - TRUE if the vertex is a ghost point

983:   Level: intermediate

985: .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices, DMNetworkGetVertexRange
986: @*/
987: PetscErrorCode DMNetworkIsGhostVertex(DM dm,PetscInt p,PetscBool *isghost)
988: {
990:   DM_Network     *network = (DM_Network*)dm->data;
991:   PetscInt       offsetg;
992:   PetscSection   sectiong;

995:   *isghost = PETSC_FALSE;
996:   DMGetDefaultGlobalSection(network->plex,&sectiong);
997:   PetscSectionGetOffset(sectiong,p,&offsetg);
998:   if (offsetg < 0) *isghost = PETSC_TRUE;
999:   return(0);
1000: }

1002: PetscErrorCode DMSetUp_Network(DM dm)
1003: {
1005:   DM_Network     *network=(DM_Network*)dm->data;

1008:   DMNetworkComponentSetUp(dm);
1009:   DMNetworkVariablesSetUp(dm);

1011:   DMSetDefaultSection(network->plex,network->DofSection);
1012:   DMGetDefaultGlobalSection(network->plex,&network->GlobalDofSection);
1013:   return(0);
1014: }

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

1020:     Collective

1022:     Input Parameters:
1023: +   dm - The DMNetwork object
1024: .   eflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for edges
1025: -   vflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for vertices

1027:     Level: intermediate

1029: @*/
1030: PetscErrorCode DMNetworkHasJacobian(DM dm,PetscBool eflg,PetscBool vflg)
1031: {
1032:   DM_Network     *network=(DM_Network*)dm->data;

1035:   network->userEdgeJacobian   = eflg;
1036:   network->userVertexJacobian = vflg;
1037:   return(0);
1038: }

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

1043:     Not Collective

1045:     Input Parameters:
1046: +   dm - The DMNetwork object
1047: .   p  - the edge point
1048: -   J - array (size = 3) of Jacobian submatrices for this edge point:
1049:         J[0]: this edge
1050:         J[1] and J[2]: connected vertices, obtained by calling DMNetworkGetConnectedVertices()

1052:     Level: intermediate

1054: .seealso: DMNetworkVertexSetMatrix
1055: @*/
1056: PetscErrorCode DMNetworkEdgeSetMatrix(DM dm,PetscInt p,Mat J[])
1057: {
1059:   DM_Network     *network=(DM_Network*)dm->data;

1062:   if (!network->userEdgeJacobian) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkEdgeSetMatrix");
1063:   if (!network->Je) {
1064:     PetscCalloc1(3*network->nEdges,&network->Je);
1065:   }
1066:   network->Je[3*p]   = J[0];
1067:   network->Je[3*p+1] = J[1];
1068:   network->Je[3*p+2] = J[2];
1069:   return(0);
1070: }

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

1075:     Not Collective

1077:     Input Parameters:
1078: +   dm - The DMNetwork object
1079: .   p  - the vertex point
1080: -   J - array of Jacobian (size = 2*(num of supporting edges) + 1) submatrices for this vertex point:
1081:         J[0]:       this vertex
1082:         J[1+2*i]:   i-th supporting edge
1083:         J[1+2*i+1]: i-th connected vertex

1085:     Level: intermediate

1087: .seealso: DMNetworkEdgeSetMatrix
1088: @*/
1089: PetscErrorCode DMNetworkVertexSetMatrix(DM dm,PetscInt p,Mat J[])
1090: {
1092:   DM_Network     *network=(DM_Network*)dm->data;
1093:   PetscInt       i,*vptr,nedges,vStart,vEnd;
1094:   const PetscInt *edges;

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

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

1101:   if (!network->Jv) {
1102:     PetscInt nVertices = network->nVertices,nedges_total;
1103:     /* count nvertex_total */
1104:     nedges_total = 0;
1105:     DMNetworkGetVertexRange(dm,&vStart,&vEnd);
1106:     PetscMalloc1(nVertices+1,&vptr);

1108:     vptr[0] = 0;
1109:     for (i=0; i<nVertices; i++) {
1110:       DMNetworkGetSupportingEdges(dm,i+vStart,&nedges,&edges);
1111:       nedges_total += nedges;
1112:       vptr[i+1] = vptr[i] + 2*nedges + 1;
1113:     }

1115:     PetscCalloc1(2*nedges_total+nVertices,&network->Jv);
1116:     network->Jvptr = vptr;
1117:   }

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

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

1128: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1129: {
1131:   PetscInt       j;
1132:   PetscScalar    val=(PetscScalar)ncols;

1135:   if (!ghost) {
1136:     for (j=0; j<nrows; j++) {
1137:       VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
1138:     }
1139:   } else {
1140:     for (j=0; j<nrows; j++) {
1141:       VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
1142:     }
1143:   }
1144:   return(0);
1145: }

1147: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1148: {
1150:   PetscInt       j,ncols_u;
1151:   PetscScalar    val;

1154:   if (!ghost) {
1155:     for (j=0; j<nrows; j++) {
1156:       MatGetRow(Ju,j,&ncols_u,NULL,NULL);
1157:       val = (PetscScalar)ncols_u;
1158:       VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
1159:       MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
1160:     }
1161:   } else {
1162:     for (j=0; j<nrows; j++) {
1163:       MatGetRow(Ju,j,&ncols_u,NULL,NULL);
1164:       val = (PetscScalar)ncols_u;
1165:       VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
1166:       MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
1167:     }
1168:   }
1169:   return(0);
1170: }

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

1177:   if (Ju) {
1178:     MatSetPreallocationUserblock_private(Ju,nrows,rows,ncols,ghost,vdnz,vonz);
1179:   } else {
1180:     MatSetPreallocationDenseblock_private(nrows,rows,ncols,ghost,vdnz,vonz);
1181:   }
1182:   return(0);
1183: }

1185: PETSC_STATIC_INLINE PetscErrorCode MatSetDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1186: {
1188:   PetscInt       j,*cols;
1189:   PetscScalar    *zeros;

1192:   PetscCalloc2(ncols,&cols,nrows*ncols,&zeros);
1193:   for (j=0; j<ncols; j++) cols[j] = j+ cstart;
1194:   MatSetValues(*J,nrows,rows,ncols,cols,zeros,INSERT_VALUES);
1195:   PetscFree2(cols,zeros);
1196:   return(0);
1197: }

1199: PETSC_STATIC_INLINE PetscErrorCode MatSetUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1200: {
1202:   PetscInt       j,M,N,row,col,ncols_u;
1203:   const PetscInt *cols;
1204:   PetscScalar    zero=0.0;

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

1210:   for (row=0; row<nrows; row++) {
1211:     MatGetRow(Ju,row,&ncols_u,&cols,NULL);
1212:     for (j=0; j<ncols_u; j++) {
1213:       col = cols[j] + cstart;
1214:       MatSetValues(*J,1,&rows[row],1,&col,&zero,INSERT_VALUES);
1215:     }
1216:     MatRestoreRow(Ju,row,&ncols_u,&cols,NULL);
1217:   }
1218:   return(0);
1219: }

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

1226:   if (Ju) {
1227:     MatSetUserblock_private(Ju,nrows,rows,ncols,cstart,J);
1228:   } else {
1229:     MatSetDenseblock_private(nrows,rows,ncols,cstart,J);
1230:   }
1231:   return(0);
1232: }

1234: /* 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.
1235: */
1236: PetscErrorCode CreateSubGlobalToLocalMapping_private(PetscSection globalsec, PetscSection localsec, ISLocalToGlobalMapping *ltog)
1237: {
1239:   PetscInt       i, size, dof;
1240:   PetscInt       *glob2loc;

1243:   PetscSectionGetStorageSize(localsec,&size);
1244:   PetscMalloc1(size,&glob2loc);

1246:   for (i = 0; i < size; i++) {
1247:     PetscSectionGetOffset(globalsec,i,&dof);
1248:     dof = (dof >= 0) ? dof : -(dof + 1);
1249:     glob2loc[i] = dof;
1250:   }

1252:   ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,size,glob2loc,PETSC_OWN_POINTER,ltog);
1253: #if 0
1254:   PetscIntView(size,glob2loc,PETSC_VIEWER_STDOUT_WORLD);
1255: #endif
1256:   return(0);
1257: }

1259: PetscErrorCode DMCreateMatrix_Network(DM dm,Mat *J)
1260: {
1262:   PetscMPIInt    rank, size;
1263:   DM_Network     *network = (DM_Network*) dm->data;
1264:   PetscInt       eStart,eEnd,vStart,vEnd,rstart,nrows,*rows,localSize;
1265:   PetscInt       cstart,ncols,j,e,v;
1266:   PetscBool      ghost,ghost_vc,ghost2,isNest;
1267:   Mat            Juser;
1268:   PetscSection   sectionGlobal;
1269:   PetscInt       nedges,*vptr=NULL,vc,*rows_v; /* suppress maybe-uninitialized warning */
1270:   const PetscInt *edges,*cone;
1271:   MPI_Comm       comm;
1272:   MatType        mtype;
1273:   Vec            vd_nz,vo_nz;
1274:   PetscInt       *dnnz,*onnz;
1275:   PetscScalar    *vdnz,*vonz;

1278:   mtype = dm->mattype;
1279:   PetscStrcmp(mtype, MATNEST, &isNest);

1281:   if (isNest) {
1282:     /* DMCreateMatrix_Network_Nest(); */
1283:     PetscInt   eDof, vDof;
1284:     Mat        j11, j12, j21, j22, bA[2][2];
1285:     ISLocalToGlobalMapping eISMap, vISMap;

1287:     PetscObjectGetComm((PetscObject)dm,&comm);
1288:     MPI_Comm_rank(comm,&rank);
1289:     MPI_Comm_size(comm,&size);

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

1294:     MatCreate(PETSC_COMM_WORLD, &j11);
1295:     MatSetSizes(j11, eDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
1296:     MatSetType(j11, MATMPIAIJ);

1298:     MatCreate(PETSC_COMM_WORLD, &j12);
1299:     MatSetSizes(j12, eDof, vDof, PETSC_DETERMINE ,PETSC_DETERMINE);
1300:     MatSetType(j12, MATMPIAIJ);

1302:     MatCreate(PETSC_COMM_WORLD, &j21);
1303:     MatSetSizes(j21, vDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
1304:     MatSetType(j21, MATMPIAIJ);

1306:     MatCreate(PETSC_COMM_WORLD, &j22);
1307:     MatSetSizes(j22, vDof, vDof, PETSC_DETERMINE, PETSC_DETERMINE);
1308:     MatSetType(j22, MATMPIAIJ);

1310:     bA[0][0] = j11;
1311:     bA[0][1] = j12;
1312:     bA[1][0] = j21;
1313:     bA[1][1] = j22;

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

1318:     MatSetLocalToGlobalMapping(j11,eISMap,eISMap);
1319:     MatSetLocalToGlobalMapping(j12,eISMap,vISMap);
1320:     MatSetLocalToGlobalMapping(j21,vISMap,eISMap);
1321:     MatSetLocalToGlobalMapping(j22,vISMap,vISMap);

1323:     MatSetUp(j11);
1324:     MatSetUp(j12);
1325:     MatSetUp(j21);
1326:     MatSetUp(j22);

1328:     MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&bA[0][0],J);
1329:     MatSetUp(*J);
1330:     MatNestSetVecType(*J,VECNEST);
1331:     MatDestroy(&j11);
1332:     MatDestroy(&j12);
1333:     MatDestroy(&j21);
1334:     MatDestroy(&j22);

1336:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1337:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);

1339:     /* Free structures */
1340:     ISLocalToGlobalMappingDestroy(&eISMap);
1341:     ISLocalToGlobalMappingDestroy(&vISMap);

1343:     return(0);
1344:   } else if (!network->userEdgeJacobian && !network->userVertexJacobian) {
1345:     /* user does not provide Jacobian blocks */
1346:     DMCreateMatrix(network->plex,J);
1347:     MatSetDM(*J,dm);
1348:     return(0);
1349:   }

1351:   MatCreate(PetscObjectComm((PetscObject)dm),J);
1352:   DMGetDefaultGlobalSection(network->plex,&sectionGlobal);
1353:   PetscSectionGetConstrainedStorageSize(sectionGlobal,&localSize);
1354:   MatSetSizes(*J,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);

1356:   MatSetType(*J,MATAIJ);
1357:   MatSetFromOptions(*J);

1359:   /* (1) Set matrix preallocation */
1360:   /*------------------------------*/
1361:   PetscObjectGetComm((PetscObject)dm,&comm);
1362:   VecCreate(comm,&vd_nz);
1363:   VecSetSizes(vd_nz,localSize,PETSC_DECIDE);
1364:   VecSetFromOptions(vd_nz);
1365:   VecSet(vd_nz,0.0);
1366:   VecDuplicate(vd_nz,&vo_nz);

1368:   /* Set preallocation for edges */
1369:   /*-----------------------------*/
1370:   DMNetworkGetEdgeRange(dm,&eStart,&eEnd);

1372:   PetscMalloc1(localSize,&rows);
1373:   for (e=eStart; e<eEnd; e++) {
1374:     /* Get row indices */
1375:     DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
1376:     DMNetworkGetNumVariables(dm,e,&nrows);
1377:     if (nrows) {
1378:       if (!network->Je) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Uer must provide Je");
1379:       for (j=0; j<nrows; j++) rows[j] = j + rstart;

1381:       /* Set preallocation for conntected vertices */
1382:       DMNetworkGetConnectedVertices(dm,e,&cone);
1383:       for (v=0; v<2; v++) {
1384:         DMNetworkGetNumVariables(dm,cone[v],&ncols);

1386:         Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
1387:         DMNetworkIsGhostVertex(dm,cone[v],&ghost);
1388:         MatSetPreallocationblock_private(Juser,nrows,rows,ncols,ghost,vd_nz,vo_nz);
1389:       }

1391:       /* Set preallocation for edge self */
1392:       cstart = rstart;
1393:       Juser = network->Je[3*e]; /* Jacobian(e,e) */
1394:       MatSetPreallocationblock_private(Juser,nrows,rows,nrows,PETSC_FALSE,vd_nz,vo_nz);
1395:     }
1396:   }

1398:   /* Set preallocation for vertices */
1399:   /*--------------------------------*/
1400:   DMNetworkGetVertexRange(dm,&vStart,&vEnd);
1401:   if (vEnd - vStart) {
1402:     if (!network->Jv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Uer must provide Jv");
1403:     vptr = network->Jvptr;
1404:     if (!vptr) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Uer must provide vptr");
1405:   }

1407:   for (v=vStart; v<vEnd; v++) {
1408:     /* Get row indices */
1409:     DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
1410:     DMNetworkGetNumVariables(dm,v,&nrows);
1411:     if (!nrows) continue;

1413:     DMNetworkIsGhostVertex(dm,v,&ghost);
1414:     if (ghost) {
1415:       PetscMalloc1(nrows,&rows_v);
1416:     } else {
1417:       rows_v = rows;
1418:     }

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

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

1425:     for (e=0; e<nedges; e++) {
1426:       /* Supporting edges */
1427:       DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
1428:       DMNetworkGetNumVariables(dm,edges[e],&ncols);

1430:       Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */
1431:       MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost,vd_nz,vo_nz);

1433:       /* Connected vertices */
1434:       DMNetworkGetConnectedVertices(dm,edges[e],&cone);
1435:       vc = (v == cone[0]) ? cone[1]:cone[0];
1436:       DMNetworkIsGhostVertex(dm,vc,&ghost_vc);

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

1440:       Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
1441:       if (ghost_vc||ghost) {
1442:         ghost2 = PETSC_TRUE;
1443:       } else {
1444:         ghost2 = PETSC_FALSE;
1445:       }
1446:       MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost2,vd_nz,vo_nz);
1447:     }

1449:     /* Set preallocation for vertex self */
1450:     DMNetworkIsGhostVertex(dm,v,&ghost);
1451:     if (!ghost) {
1452:       DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
1453:       Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
1454:       MatSetPreallocationblock_private(Juser,nrows,rows_v,nrows,PETSC_FALSE,vd_nz,vo_nz);
1455:     }
1456:     if (ghost) {
1457:       PetscFree(rows_v);
1458:     }
1459:   }

1461:   VecAssemblyBegin(vd_nz);
1462:   VecAssemblyBegin(vo_nz);

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

1466:   VecAssemblyEnd(vd_nz);
1467:   VecAssemblyEnd(vo_nz);

1469:   VecGetArray(vd_nz,&vdnz);
1470:   VecGetArray(vo_nz,&vonz);
1471:   for (j=0; j<localSize; j++) {
1472:     dnnz[j] = (PetscInt)PetscRealPart(vdnz[j]);
1473:     onnz[j] = (PetscInt)PetscRealPart(vonz[j]);
1474:   }
1475:   VecRestoreArray(vd_nz,&vdnz);
1476:   VecRestoreArray(vo_nz,&vonz);
1477:   VecDestroy(&vd_nz);
1478:   VecDestroy(&vo_nz);

1480:   MatSeqAIJSetPreallocation(*J,0,dnnz);
1481:   MatMPIAIJSetPreallocation(*J,0,dnnz,0,onnz);
1482:   MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);

1484:   PetscFree2(dnnz,onnz);

1486:   /* (2) Set matrix entries for edges */
1487:   /*----------------------------------*/
1488:   for (e=eStart; e<eEnd; e++) {
1489:     /* Get row indices */
1490:     DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
1491:     DMNetworkGetNumVariables(dm,e,&nrows);
1492:     if (nrows) {
1493:       if (!network->Je) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Uer must provide Je");

1495:       for (j=0; j<nrows; j++) rows[j] = j + rstart;

1497:       /* Set matrix entries for conntected vertices */
1498:       DMNetworkGetConnectedVertices(dm,e,&cone);
1499:       for (v=0; v<2; v++) {
1500:         DMNetworkGetVariableGlobalOffset(dm,cone[v],&cstart);
1501:         DMNetworkGetNumVariables(dm,cone[v],&ncols);

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

1507:       /* Set matrix entries for edge self */
1508:       cstart = rstart;
1509:       Juser = network->Je[3*e]; /* Jacobian(e,e) */
1510:       MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);
1511:     }
1512:   }

1514:   /* Set matrix entries for vertices */
1515:   /*---------------------------------*/
1516:   for (v=vStart; v<vEnd; v++) {
1517:     /* Get row indices */
1518:     DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
1519:     DMNetworkGetNumVariables(dm,v,&nrows);
1520:     if (!nrows) continue;

1522:     DMNetworkIsGhostVertex(dm,v,&ghost);
1523:     if (ghost) {
1524:       PetscMalloc1(nrows,&rows_v);
1525:     } else {
1526:       rows_v = rows;
1527:     }
1528:     for (j=0; j<nrows; j++) rows_v[j] = j + rstart;

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

1533:     for (e=0; e<nedges; e++) {
1534:       /* Supporting edges */
1535:       DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
1536:       DMNetworkGetNumVariables(dm,edges[e],&ncols);

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

1541:       /* Connected vertices */
1542:       DMNetworkGetConnectedVertices(dm,edges[e],&cone);
1543:       vc = (v == cone[0]) ? cone[1]:cone[0];

1545:       DMNetworkGetVariableGlobalOffset(dm,vc,&cstart);
1546:       DMNetworkGetNumVariables(dm,vc,&ncols);

1548:       Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
1549:       MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);
1550:     }

1552:     /* Set matrix entries for vertex self */
1553:     if (!ghost) {
1554:       DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
1555:       Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
1556:       MatSetblock_private(Juser,nrows,rows_v,nrows,cstart,J);
1557:     }
1558:     if (ghost) {
1559:       PetscFree(rows_v);
1560:     }
1561:   }
1562:   PetscFree(rows);

1564:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1565:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);

1567:   MatSetDM(*J,dm);
1568:   return(0);
1569: }

1571: PetscErrorCode DMDestroy_Network(DM dm)
1572: {
1574:   DM_Network     *network = (DM_Network*) dm->data;

1577:   if (--network->refct > 0) return(0);
1578:   if (network->Je) {
1579:     PetscFree(network->Je);
1580:   }
1581:   if (network->Jv) {
1582:     PetscFree(network->Jvptr);
1583:     PetscFree(network->Jv);
1584:   }

1586:   ISLocalToGlobalMappingDestroy(&network->vertex.mapping);
1587:   PetscSectionDestroy(&network->vertex.DofSection);
1588:   PetscSectionDestroy(&network->vertex.GlobalDofSection);
1589:   if (network->vertex.sf) {
1590:     PetscSFDestroy(&network->vertex.sf);
1591:   }
1592:   /* edge */
1593:   ISLocalToGlobalMappingDestroy(&network->edge.mapping);
1594:   PetscSectionDestroy(&network->edge.DofSection);
1595:   PetscSectionDestroy(&network->edge.GlobalDofSection);
1596:   if (network->edge.sf) {
1597:     PetscSFDestroy(&network->edge.sf);
1598:   }
1599:   DMDestroy(&network->plex);
1600:   network->edges = NULL;
1601:   PetscSectionDestroy(&network->DataSection);
1602:   PetscSectionDestroy(&network->DofSection);

1604:   PetscFree(network->componentdataarray);
1605:   PetscFree(network->cvalue);
1606:   PetscFree(network->header);
1607:   PetscFree(network);
1608:   return(0);
1609: }

1611: PetscErrorCode DMView_Network(DM dm, PetscViewer viewer)
1612: {
1614:   DM_Network     *network = (DM_Network*) dm->data;

1617:   DMView(network->plex,viewer);
1618:   return(0);
1619: }

1621: PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l)
1622: {
1624:   DM_Network     *network = (DM_Network*) dm->data;

1627:   DMGlobalToLocalBegin(network->plex,g,mode,l);
1628:   return(0);
1629: }

1631: PetscErrorCode DMGlobalToLocalEnd_Network(DM dm, Vec g, InsertMode mode, Vec l)
1632: {
1634:   DM_Network     *network = (DM_Network*) dm->data;

1637:   DMGlobalToLocalEnd(network->plex,g,mode,l);
1638:   return(0);
1639: }

1641: PetscErrorCode DMLocalToGlobalBegin_Network(DM dm, Vec l, InsertMode mode, Vec g)
1642: {
1644:   DM_Network     *network = (DM_Network*) dm->data;

1647:   DMLocalToGlobalBegin(network->plex,l,mode,g);
1648:   return(0);
1649: }

1651: PetscErrorCode DMLocalToGlobalEnd_Network(DM dm, Vec l, InsertMode mode, Vec g)
1652: {
1654:   DM_Network     *network = (DM_Network*) dm->data;

1657:   DMLocalToGlobalEnd(network->plex,l,mode,g);
1658:   return(0);
1659: }