Actual source code: network.c
petsc-3.13.6 2020-09-29
1: #include <petsc/private/dmnetworkimpl.h>
3: /*@
4: DMNetworkGetPlex - Gets the Plex DM associated with this network DM
6: Not collective
8: Input Parameters:
9: + netdm - the dm object
10: - plexmdm - the plex dm object
12: Level: Advanced
14: .seealso: DMNetworkCreate()
15: @*/
16: PetscErrorCode DMNetworkGetPlex(DM netdm, DM *plexdm)
17: {
18: DM_Network *network = (DM_Network*) netdm->data;
21: *plexdm = network->plex;
22: return(0);
23: }
25: /*@
26: DMNetworkGetSizes - Gets the the number of subnetworks and coupling subnetworks
28: Collective on dm
30: Input Parameters:
31: + dm - the dm object
32: . Nsubnet - global number of subnetworks
33: - NsubnetCouple - global number of coupling subnetworks
35: Level: beginner
37: .seealso: DMNetworkCreate()
38: @*/
39: PetscErrorCode DMNetworkGetSizes(DM netdm, PetscInt *Nsubnet, PetscInt *Ncsubnet)
40: {
41: DM_Network *network = (DM_Network*) netdm->data;
44: *Nsubnet = network->nsubnet;
45: *Ncsubnet = network->ncsubnet;
46: return(0);
47: }
49: /*@
50: DMNetworkSetSizes - Sets the number of subnetworks, local and global vertices and edges for each subnetwork.
52: Collective on dm
54: Input Parameters:
55: + dm - the dm object
56: . Nsubnet - global number of subnetworks
57: . nV - number of local vertices for each subnetwork
58: . nE - number of local edges for each subnetwork
59: . NsubnetCouple - global number of coupling subnetworks
60: - nec - number of local edges for each coupling subnetwork
62: You cannot change the sizes once they have been set. nV, nE are arrays of length Nsubnet, and nec is array of length NsubnetCouple.
64: Level: beginner
66: .seealso: DMNetworkCreate()
67: @*/
68: PetscErrorCode DMNetworkSetSizes(DM dm,PetscInt Nsubnet,PetscInt nV[], PetscInt nE[],PetscInt NsubnetCouple,PetscInt nec[])
69: {
71: DM_Network *network = (DM_Network*) dm->data;
72: PetscInt a[2],b[2],i;
76: if (Nsubnet <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subnetworks %D cannot be less than 1",Nsubnet);
77: if (NsubnetCouple < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of coupling subnetworks %D cannot be less than 0",NsubnetCouple);
81: if (network->nsubnet != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Network sizes alread set, cannot resize the network");
83: if (!nV || !nE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Local vertex size or edge size must be provided");
85: network->nsubnet = Nsubnet + NsubnetCouple;
86: network->ncsubnet = NsubnetCouple;
87: PetscCalloc1(Nsubnet+NsubnetCouple,&network->subnet);
89: /* ----------------------------------------------------------
90: p=v or e; P=V or E
91: subnet[0].pStart = 0
92: subnet[i+1].pStart = subnet[i].pEnd = subnet[i].pStart + (nE[i] or NV[i])
93: ----------------------------------------------------------------------- */
94: for (i=0; i < Nsubnet; i++) {
95: /* Get global number of vertices and edges for subnet[i] */
96: a[0] = nV[i]; a[1] = nE[i]; /* local number of vertices (excluding ghost) and edges */
97: MPIU_Allreduce(a,b,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
98: network->subnet[i].Nvtx = b[0];
99: network->subnet[i].Nedge = b[1];
101: network->subnet[i].nvtx = nV[i]; /* local nvtx, without ghost */
103: /* global subnet[].vStart and vEnd, used by DMNetworkLayoutSetUp() */
104: network->subnet[i].vStart = network->NVertices;
105: network->subnet[i].vEnd = network->subnet[i].vStart + network->subnet[i].Nvtx;
107: network->nVertices += nV[i];
108: network->NVertices += network->subnet[i].Nvtx;
110: network->subnet[i].nedge = nE[i];
111: network->subnet[i].eStart = network->nEdges;
112: network->subnet[i].eEnd = network->subnet[i].eStart + nE[i];
113: network->nEdges += nE[i];
114: network->NEdges += network->subnet[i].Nedge;
115: }
117: /* coupling subnetwork */
118: for (; i < Nsubnet+NsubnetCouple; i++) {
119: /* Get global number of coupling edges for subnet[i] */
120: MPIU_Allreduce(nec+(i-Nsubnet),&network->subnet[i].Nedge,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
122: network->subnet[i].nvtx = 0; /* We design coupling subnetwork such that it does not have its own vertices */
123: network->subnet[i].vStart = network->nVertices;
124: network->subnet[i].vEnd = network->subnet[i].vStart;
126: network->subnet[i].nedge = nec[i-Nsubnet];
127: network->subnet[i].eStart = network->nEdges;
128: network->subnet[i].eEnd = network->subnet[i].eStart + nec[i-Nsubnet];
129: network->nEdges += nec[i-Nsubnet];
130: network->NEdges += network->subnet[i].Nedge;
131: }
132: return(0);
133: }
135: /*@
136: DMNetworkSetEdgeList - Sets the list of local edges (vertex connectivity) for the network
138: Logically collective on dm
140: Input Parameters:
141: + dm - the dm object
142: . edgelist - list of edges for each subnetwork
143: - edgelistCouple - list of edges for each coupling subnetwork
145: Notes:
146: There is no copy involved in this operation, only the pointer is referenced. The edgelist should
147: not be destroyed before the call to DMNetworkLayoutSetUp
149: Level: beginner
151: Example usage:
152: Consider the following 2 separate networks and a coupling network:
154: .vb
155: network 0: v0 -> v1 -> v2 -> v3
156: network 1: v1 -> v2 -> v0
157: coupling network: network 1: v2 -> network 0: v0
158: .ve
160: The resulting input
161: edgelist[0] = [0 1 | 1 2 | 2 3];
162: edgelist[1] = [1 2 | 2 0]
163: edgelistCouple[0] = [(network)1 (v)2 (network)0 (v)0].
165: .seealso: DMNetworkCreate, DMNetworkSetSizes
166: @*/
167: PetscErrorCode DMNetworkSetEdgeList(DM dm,PetscInt *edgelist[],PetscInt *edgelistCouple[])
168: {
169: DM_Network *network = (DM_Network*) dm->data;
170: PetscInt i;
173: for (i=0; i < (network->nsubnet-network->ncsubnet); i++) network->subnet[i].edgelist = edgelist[i];
174: if (network->ncsubnet) {
175: PetscInt j = 0;
176: if (!edgelistCouple) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Must provide edgelist_couple");
177: while (i < network->nsubnet) network->subnet[i++].edgelist = edgelistCouple[j++];
178: }
179: return(0);
180: }
182: /*@
183: DMNetworkLayoutSetUp - Sets up the bare layout (graph) for the network
185: Collective on dm
187: Input Parameters:
188: . DM - the dmnetwork object
190: Notes:
191: This routine should be called after the network sizes and edgelists have been provided. It creates
192: the bare layout of the network and sets up the network to begin insertion of components.
194: All the components should be registered before calling this routine.
196: Level: beginner
198: .seealso: DMNetworkSetSizes, DMNetworkSetEdgeList
199: @*/
200: PetscErrorCode DMNetworkLayoutSetUp(DM dm)
201: {
203: DM_Network *network = (DM_Network*)dm->data;
204: PetscInt numCorners=2,spacedim=2,dim = 1; /* One dimensional network */
205: PetscReal *vertexcoords=NULL;
206: PetscInt i,j,ctr,nsubnet,*eowners,np,*edges,*subnetvtx,vStart;
207: PetscInt k,netid,vid, *vidxlTog,*edgelist_couple=NULL;
208: const PetscInt *cone;
209: MPI_Comm comm;
210: PetscMPIInt size,rank;
213: PetscObjectGetComm((PetscObject)dm,&comm);
214: MPI_Comm_rank(comm,&rank);
215: MPI_Comm_size(comm,&size);
217: /* Create the local edgelist for the network by concatenating local input edgelists of the subnetworks */
218: PetscCalloc2(numCorners*network->nVertices,&vertexcoords,2*network->nEdges,&edges);
219: nsubnet = network->nsubnet - network->ncsubnet;
220: ctr = 0;
221: for (i=0; i < nsubnet; i++) {
222: for (j = 0; j < network->subnet[i].nedge; j++) {
223: edges[2*ctr] = network->subnet[i].vStart + network->subnet[i].edgelist[2*j];
224: edges[2*ctr+1] = network->subnet[i].vStart + network->subnet[i].edgelist[2*j+1];
225: ctr++;
226: }
227: }
229: /* Append local coupling edgelists of the subnetworks */
230: i = nsubnet; /* netid of coupling subnet */
231: nsubnet = network->nsubnet;
232: while (i < nsubnet) {
233: edgelist_couple = network->subnet[i].edgelist;
235: k = 0;
236: for (j = 0; j < network->subnet[i].nedge; j++) {
237: netid = edgelist_couple[k]; vid = edgelist_couple[k+1];
238: edges[2*ctr] = network->subnet[netid].vStart + vid; k += 2;
240: netid = edgelist_couple[k]; vid = edgelist_couple[k+1];
241: edges[2*ctr+1] = network->subnet[netid].vStart + vid; k+=2;
242: ctr++;
243: }
244: i++;
245: }
246: /*
247: if (rank == 0) {
248: PetscPrintf(PETSC_COMM_SELF,"[%d] edgelist:\n",rank);
249: for(i=0; i < network->nEdges; i++) {
250: PetscPrintf(PETSC_COMM_SELF,"[%D %D]",edges[2*i],edges[2*i+1]);
251: printf("\n");
252: }
253: }
254: */
256: /* Create network->plex */
257: #if defined(PETSC_USE_64BIT_INDICES)
258: {
259: int *edges64;
260: np = network->nEdges*numCorners;
261: PetscMalloc1(np,&edges64);
262: for (i=0; i<np; i++) edges64[i] = (int)edges[i];
264: if (size == 1) {
265: DMPlexCreateFromCellList(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges64,spacedim,(const double*)vertexcoords,&network->plex);
266: } else {
267: DMPlexCreateFromCellListParallel(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges64,spacedim,(const PetscReal*)vertexcoords,NULL,&network->plex);
268: }
269: PetscFree(edges64);
270: }
271: #else
272: if (size == 1) {
273: DMPlexCreateFromCellList(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges,spacedim,(const double*)vertexcoords,&network->plex);
274: } else {
275: DMPlexCreateFromCellListParallel(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges,spacedim,(const PetscReal*)vertexcoords,NULL,&network->plex);
276: }
277: #endif
279: DMPlexGetChart(network->plex,&network->pStart,&network->pEnd);
280: DMPlexGetHeightStratum(network->plex,0,&network->eStart,&network->eEnd);
281: DMPlexGetHeightStratum(network->plex,1,&network->vStart,&network->vEnd);
282: vStart = network->vStart;
284: PetscSectionCreate(comm,&network->DataSection);
285: PetscSectionCreate(comm,&network->DofSection);
286: PetscSectionSetChart(network->DataSection,network->pStart,network->pEnd);
287: PetscSectionSetChart(network->DofSection,network->pStart,network->pEnd);
289: network->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
290: np = network->pEnd - network->pStart;
291: PetscCalloc2(np,&network->header,np,&network->cvalue);
293: /* Create vidxlTog: maps local vertex index to global index */
294: np = network->vEnd - vStart;
295: PetscMalloc2(np,&vidxlTog,size+1,&eowners);
296: ctr = 0;
297: for (i=network->eStart; i<network->eEnd; i++) {
298: DMNetworkGetConnectedVertices(dm,i,&cone);
299: vidxlTog[cone[0] - vStart] = edges[2*ctr];
300: vidxlTog[cone[1] - vStart] = edges[2*ctr+1];
301: ctr++;
302: }
303: PetscFree2(vertexcoords,edges);
305: /* Create vertices and edges array for the subnetworks */
306: for (j=0; j < network->nsubnet; j++) {
307: PetscCalloc1(network->subnet[j].nedge,&network->subnet[j].edges);
309: /* Temporarily setting nvtx and nedge to 0 so we can use them as counters in the below for loop.
310: These get updated when the vertices and edges are added. */
311: network->subnet[j].nvtx = 0;
312: network->subnet[j].nedge = 0;
313: }
314: PetscCalloc1(np,&network->subnetvtx);
317: /* Get edge ownership */
318: np = network->eEnd - network->eStart;
319: MPI_Allgather(&np,1,MPIU_INT,eowners+1,1,MPIU_INT,comm);
320: eowners[0] = 0;
321: for (i=2; i<=size; i++) eowners[i] += eowners[i-1];
323: i = 0; j = 0;
324: while (i < np) { /* local edges, including coupling edges */
325: network->header[i].index = i + eowners[rank]; /* Global edge index */
327: if (j < network->nsubnet && i < network->subnet[j].eEnd) {
328: network->header[i].subnetid = j; /* Subnetwork id */
329: network->subnet[j].edges[network->subnet[j].nedge++] = i;
331: network->header[i].ndata = 0;
332: PetscSectionAddDof(network->DataSection,i,network->dataheadersize);
333: network->header[i].offset[0] = 0;
334: network->header[i].offsetvarrel[0] = 0;
335: i++;
336: }
337: if (i >= network->subnet[j].eEnd) j++;
338: }
340: /* Count network->subnet[*].nvtx */
341: for (i=vStart; i<network->vEnd; i++) { /* local vertices, including ghosts */
342: k = vidxlTog[i-vStart];
343: for (j=0; j < network->nsubnet; j++) {
344: if (network->subnet[j].vStart <= k && k < network->subnet[j].vEnd) {
345: network->subnet[j].nvtx++;
346: break;
347: }
348: }
349: }
351: /* Set network->subnet[*].vertices on array network->subnetvtx */
352: subnetvtx = network->subnetvtx;
353: for (j=0; j<network->nsubnet; j++) {
354: network->subnet[j].vertices = subnetvtx;
355: subnetvtx += network->subnet[j].nvtx;
356: network->subnet[j].nvtx = 0;
357: }
359: /* Set vertex array for the subnetworks */
360: for (i=vStart; i<network->vEnd; i++) { /* local vertices, including ghosts */
361: network->header[i].index = vidxlTog[i-vStart]; /* Global vertex index */
363: k = vidxlTog[i-vStart];
364: for (j=0; j < network->nsubnet; j++) {
365: if (network->subnet[j].vStart <= k && k < network->subnet[j].vEnd) {
366: network->header[i].subnetid = j;
367: network->subnet[j].vertices[network->subnet[j].nvtx++] = i;
368: break;
369: }
370: }
372: network->header[i].ndata = 0;
373: PetscSectionAddDof(network->DataSection,i,network->dataheadersize);
374: network->header[i].offset[0] = 0;
375: network->header[i].offsetvarrel[0] = 0;
376: }
378: PetscFree2(vidxlTog,eowners);
379: return(0);
380: }
382: /*@C
383: DMNetworkGetSubnetworkInfo - Returns the info for the subnetwork
385: Input Parameters:
386: + dm - the DM object
387: - id - the ID (integer) of the subnetwork
389: Output Parameters:
390: + nv - number of vertices (local)
391: . ne - number of edges (local)
392: . vtx - local vertices for this subnetwork
393: - edge - local edges for this subnetwork
395: Notes:
396: Cannot call this routine before DMNetworkLayoutSetup()
398: Level: intermediate
400: .seealso: DMNetworkLayoutSetUp, DMNetworkCreate
401: @*/
402: PetscErrorCode DMNetworkGetSubnetworkInfo(DM dm,PetscInt id,PetscInt *nv, PetscInt *ne,const PetscInt **vtx, const PetscInt **edge)
403: {
404: DM_Network *network = (DM_Network*)dm->data;
407: if (id >= network->nsubnet) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subnet ID %D exceeds the num of subnets %D",id,network->nsubnet);
408: *nv = network->subnet[id].nvtx;
409: *ne = network->subnet[id].nedge;
410: *vtx = network->subnet[id].vertices;
411: *edge = network->subnet[id].edges;
412: return(0);
413: }
415: /*@C
416: DMNetworkGetSubnetworkCoupleInfo - Returns the info for the coupling subnetwork
418: Input Parameters:
419: + dm - the DM object
420: - id - the ID (integer) of the coupling subnetwork
422: Output Parameters:
423: + ne - number of edges (local)
424: - edge - local edges for this coupling subnetwork
426: Notes:
427: Cannot call this routine before DMNetworkLayoutSetup()
429: Level: intermediate
431: .seealso: DMNetworkGetSubnetworkInfo, DMNetworkLayoutSetUp, DMNetworkCreate
432: @*/
433: PetscErrorCode DMNetworkGetSubnetworkCoupleInfo(DM dm,PetscInt id,PetscInt *ne,const PetscInt **edge)
434: {
435: DM_Network *net = (DM_Network*)dm->data;
436: PetscInt id1;
439: if (net->ncsubnet) {
440: if (id >= net->ncsubnet) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subnet ID %D exceeds the num of coupling subnets %D",id,net->ncsubnet);
442: id1 = id + net->nsubnet - net->ncsubnet;
443: *ne = net->subnet[id1].nedge;
444: *edge = net->subnet[id1].edges;
445: } else {
446: *ne = 0;
447: *edge = NULL;
448: }
449: return(0);
450: }
452: /*@C
453: DMNetworkRegisterComponent - Registers the network component
455: Logically collective on dm
457: Input Parameters:
458: + dm - the network object
459: . name - the component name
460: - size - the storage size in bytes for this component data
462: Output Parameters:
463: . key - an integer key that defines the component
465: Notes
466: This routine should be called by all processors before calling DMNetworkLayoutSetup().
468: Level: beginner
470: .seealso: DMNetworkLayoutSetUp, DMNetworkCreate
471: @*/
472: PetscErrorCode DMNetworkRegisterComponent(DM dm,const char *name,size_t size,PetscInt *key)
473: {
474: PetscErrorCode ierr;
475: DM_Network *network = (DM_Network*) dm->data;
476: DMNetworkComponent *component=&network->component[network->ncomponent];
477: PetscBool flg=PETSC_FALSE;
478: PetscInt i;
481: for (i=0; i < network->ncomponent; i++) {
482: PetscStrcmp(component->name,name,&flg);
483: if (flg) {
484: *key = i;
485: return(0);
486: }
487: }
488: if(network->ncomponent == MAX_COMPONENTS) {
489: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Number of components registered exceeds the max %D",MAX_COMPONENTS);
490: }
492: PetscStrcpy(component->name,name);
493: component->size = size/sizeof(DMNetworkComponentGenericDataType);
494: *key = network->ncomponent;
495: network->ncomponent++;
496: return(0);
497: }
499: /*@
500: DMNetworkGetVertexRange - Get the bounds [start, end) for the vertices.
502: Not Collective
504: Input Parameters:
505: . dm - The DMNetwork object
507: Output Parameters:
508: + vStart - The first vertex point
509: - vEnd - One beyond the last vertex point
511: Level: beginner
513: .seealso: DMNetworkGetEdgeRange
514: @*/
515: PetscErrorCode DMNetworkGetVertexRange(DM dm,PetscInt *vStart,PetscInt *vEnd)
516: {
517: DM_Network *network = (DM_Network*)dm->data;
520: if (vStart) *vStart = network->vStart;
521: if (vEnd) *vEnd = network->vEnd;
522: return(0);
523: }
525: /*@
526: DMNetworkGetEdgeRange - Get the bounds [start, end) for the edges.
528: Not Collective
530: Input Parameters:
531: . dm - The DMNetwork object
533: Output Parameters:
534: + eStart - The first edge point
535: - eEnd - One beyond the last edge point
537: Level: beginner
539: .seealso: DMNetworkGetVertexRange
540: @*/
541: PetscErrorCode DMNetworkGetEdgeRange(DM dm,PetscInt *eStart,PetscInt *eEnd)
542: {
543: DM_Network *network = (DM_Network*)dm->data;
546: if (eStart) *eStart = network->eStart;
547: if (eEnd) *eEnd = network->eEnd;
548: return(0);
549: }
551: /*@
552: DMNetworkGetGlobalEdgeIndex - Get the user global numbering for the edge.
554: Not Collective
556: Input Parameters:
557: + dm - DMNetwork object
558: - p - edge point
560: Output Parameters:
561: . index - user global numbering for the edge
563: Level: intermediate
565: .seealso: DMNetworkGetGlobalVertexIndex
566: @*/
567: PetscErrorCode DMNetworkGetGlobalEdgeIndex(DM dm,PetscInt p,PetscInt *index)
568: {
569: PetscErrorCode ierr;
570: DM_Network *network = (DM_Network*)dm->data;
571: PetscInt offsetp;
572: DMNetworkComponentHeader header;
575: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
576: PetscSectionGetOffset(network->DataSection,p,&offsetp);
577: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
578: *index = header->index;
579: return(0);
580: }
582: /*@
583: DMNetworkGetGlobalVertexIndex - Get the user global numbering for the vertex.
585: Not Collective
587: Input Parameters:
588: + dm - DMNetwork object
589: - p - vertex point
591: Output Parameters:
592: . index - user global numbering for the vertex
594: Level: intermediate
596: .seealso: DMNetworkGetGlobalEdgeIndex
597: @*/
598: PetscErrorCode DMNetworkGetGlobalVertexIndex(DM dm,PetscInt p,PetscInt *index)
599: {
600: PetscErrorCode ierr;
601: DM_Network *network = (DM_Network*)dm->data;
602: PetscInt offsetp;
603: DMNetworkComponentHeader header;
606: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
607: PetscSectionGetOffset(network->DataSection,p,&offsetp);
608: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
609: *index = header->index;
610: return(0);
611: }
613: /*
614: DMNetworkGetComponentKeyOffset - Gets the type along with the offset for indexing the
615: component value from the component data array
617: Not Collective
619: Input Parameters:
620: + dm - The DMNetwork object
621: . p - vertex/edge point
622: - compnum - component number
624: Output Parameters:
625: + compkey - the key obtained when registering the component
626: - offset - offset into the component data array associated with the vertex/edge point
628: Notes:
629: Typical usage:
631: DMNetworkGetComponentDataArray(dm, &arr);
632: DMNetworkGetVertex/EdgeRange(dm,&Start,&End);
633: Loop over vertices or edges
634: DMNetworkGetNumComponents(dm,v,&numcomps);
635: Loop over numcomps
636: DMNetworkGetComponentKeyOffset(dm,v,compnum,&key,&offset);
637: compdata = (UserCompDataType)(arr+offset);
639: Level: intermediate
641: .seealso: DMNetworkGetNumComponents, DMNetworkGetComponentDataArray,
642: */
643: PetscErrorCode DMNetworkGetComponentKeyOffset(DM dm,PetscInt p, PetscInt compnum, PetscInt *compkey, PetscInt *offset)
644: {
645: PetscErrorCode ierr;
646: PetscInt offsetp;
647: DMNetworkComponentHeader header;
648: DM_Network *network = (DM_Network*)dm->data;
651: PetscSectionGetOffset(network->DataSection,p,&offsetp);
652: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
653: if (compkey) *compkey = header->key[compnum];
654: if (offset) *offset = offsetp+network->dataheadersize+header->offset[compnum];
655: return(0);
656: }
658: /*@
659: DMNetworkGetComponent - Returns the network component and its key
661: Not Collective
663: Input Parameters:
664: + dm - DMNetwork object
665: . p - edge or vertex point
666: - compnum - component number
668: Output Parameters:
669: + compkey - the key set for this computing during registration
670: - component - the component data
672: Notes:
673: Typical usage:
675: DMNetworkGetVertex/EdgeRange(dm,&Start,&End);
676: Loop over vertices or edges
677: DMNetworkGetNumComponents(dm,v,&numcomps);
678: Loop over numcomps
679: DMNetworkGetComponent(dm,v,compnum,&key,&component);
681: Level: beginner
683: .seealso: DMNetworkGetNumComponents, DMNetworkGetVariableOffset
684: @*/
685: PetscErrorCode DMNetworkGetComponent(DM dm, PetscInt p, PetscInt compnum, PetscInt *key, void **component)
686: {
688: DM_Network *network = (DM_Network*)dm->data;
689: PetscInt offsetd = 0;
692: DMNetworkGetComponentKeyOffset(dm,p,compnum,key,&offsetd);
693: *component = network->componentdataarray+offsetd;
694: return(0);
695: }
697: /*@
698: DMNetworkAddComponent - Adds a network component at the given point (vertex/edge)
700: Not Collective
702: Input Parameters:
703: + dm - The DMNetwork object
704: . p - vertex/edge point
705: . componentkey - component key returned while registering the component
706: - compvalue - pointer to the data structure for the component
708: Level: beginner
710: .seealso: DMNetworkGetVertexRange, DMNetworkGetEdgeRange, DMNetworkRegisterComponent
711: @*/
712: PetscErrorCode DMNetworkAddComponent(DM dm, PetscInt p,PetscInt componentkey,void* compvalue)
713: {
714: DM_Network *network = (DM_Network*)dm->data;
715: DMNetworkComponent *component = &network->component[componentkey];
716: DMNetworkComponentHeader header = &network->header[p];
717: DMNetworkComponentValue cvalue = &network->cvalue[p];
718: PetscErrorCode ierr;
721: 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);
723: header->size[header->ndata] = component->size;
724: PetscSectionAddDof(network->DataSection,p,component->size);
725: header->key[header->ndata] = componentkey;
726: if (header->ndata != 0) header->offset[header->ndata] = header->offset[header->ndata-1] + header->size[header->ndata-1];
727: header->nvar[header->ndata] = 0;
729: cvalue->data[header->ndata] = (void*)compvalue;
730: header->ndata++;
731: return(0);
732: }
734: /*@
735: DMNetworkSetComponentNumVariables - Sets the number of variables for a component
737: Not Collective
739: Input Parameters:
740: + dm - The DMNetwork object
741: . p - vertex/edge point
742: . compnum - component number (First component added = 0, second = 1, ...)
743: - nvar - number of variables for the component
745: Level: beginner
747: .seealso: DMNetworkAddComponent(), DMNetworkGetNumComponents(),DMNetworkRegisterComponent()
748: @*/
749: PetscErrorCode DMNetworkSetComponentNumVariables(DM dm, PetscInt p,PetscInt compnum,PetscInt nvar)
750: {
751: DM_Network *network = (DM_Network*)dm->data;
752: DMNetworkComponentHeader header = &network->header[p];
753: PetscErrorCode ierr;
756: DMNetworkAddNumVariables(dm,p,nvar);
757: header->nvar[compnum] = nvar;
758: if (compnum != 0) header->offsetvarrel[compnum] = header->offsetvarrel[compnum-1] + header->nvar[compnum-1];
759: return(0);
760: }
762: /*@
763: DMNetworkGetNumComponents - Get the number of components at a vertex/edge
765: Not Collective
767: Input Parameters:
768: + dm - The DMNetwork object
769: - p - vertex/edge point
771: Output Parameters:
772: . numcomponents - Number of components at the vertex/edge
774: Level: beginner
776: .seealso: DMNetworkRegisterComponent, DMNetworkAddComponent
777: @*/
778: PetscErrorCode DMNetworkGetNumComponents(DM dm,PetscInt p,PetscInt *numcomponents)
779: {
781: PetscInt offset;
782: DM_Network *network = (DM_Network*)dm->data;
785: PetscSectionGetOffset(network->DataSection,p,&offset);
786: *numcomponents = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata;
787: return(0);
788: }
790: /*@
791: DMNetworkGetVariableOffset - Get the offset for accessing the variable associated with the given vertex/edge from the local vector.
793: Not Collective
795: Input Parameters:
796: + dm - The DMNetwork object
797: - p - the edge/vertex point
799: Output Parameters:
800: . offset - the offset
802: Level: beginner
804: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
805: @*/
806: PetscErrorCode DMNetworkGetVariableOffset(DM dm,PetscInt p,PetscInt *offset)
807: {
809: DM_Network *network = (DM_Network*)dm->data;
812: PetscSectionGetOffset(network->plex->localSection,p,offset);
813: return(0);
814: }
816: /*@
817: DMNetworkGetVariableGlobalOffset - Get the global offset for the variable associated with the given vertex/edge from the global vector.
819: Not Collective
821: Input Parameters:
822: + dm - The DMNetwork object
823: - p - the edge/vertex point
825: Output Parameters:
826: . offsetg - the offset
828: Level: beginner
830: .seealso: DMNetworkGetVariableOffset, DMGetLocalVector
831: @*/
832: PetscErrorCode DMNetworkGetVariableGlobalOffset(DM dm,PetscInt p,PetscInt *offsetg)
833: {
835: DM_Network *network = (DM_Network*)dm->data;
838: PetscSectionGetOffset(network->plex->globalSection,p,offsetg);
839: if (*offsetg < 0) *offsetg = -(*offsetg + 1); /* Convert to actual global offset for ghost vertex */
840: return(0);
841: }
843: /*@
844: DMNetworkGetComponentVariableOffset - Get the offset for accessing the variable associated with a component for the given vertex/edge from the local vector.
846: Not Collective
848: Input Parameters:
849: + dm - The DMNetwork object
850: . p - the edge/vertex point
851: - compnum - component number
853: Output Parameters:
854: . offset - the offset
856: Level: intermediate
858: .seealso: DMNetworkGetVariableGlobalOffset(), DMGetLocalVector(), DMNetworkSetComponentNumVariables()
859: @*/
860: PetscErrorCode DMNetworkGetComponentVariableOffset(DM dm,PetscInt p,PetscInt compnum,PetscInt *offset)
861: {
863: DM_Network *network = (DM_Network*)dm->data;
864: PetscInt offsetp,offsetd;
865: DMNetworkComponentHeader header;
868: DMNetworkGetVariableOffset(dm,p,&offsetp);
869: PetscSectionGetOffset(network->DataSection,p,&offsetd);
870: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetd);
871: *offset = offsetp + header->offsetvarrel[compnum];
872: return(0);
873: }
875: /*@
876: DMNetworkGetComponentVariableGlobalOffset - Get the global offset for accessing the variable associated with a component for the given vertex/edge from the local vector.
878: Not Collective
880: Input Parameters:
881: + dm - The DMNetwork object
882: . p - the edge/vertex point
883: - compnum - component number
885: Output Parameters:
886: . offsetg - the global offset
888: Level: intermediate
890: .seealso: DMNetworkGetVariableGlobalOffset(), DMNetworkGetComponentVariableOffset(), DMGetLocalVector(), DMNetworkSetComponentNumVariables()
891: @*/
892: PetscErrorCode DMNetworkGetComponentVariableGlobalOffset(DM dm,PetscInt p,PetscInt compnum,PetscInt *offsetg)
893: {
895: DM_Network *network = (DM_Network*)dm->data;
896: PetscInt offsetp,offsetd;
897: DMNetworkComponentHeader header;
900: DMNetworkGetVariableGlobalOffset(dm,p,&offsetp);
901: PetscSectionGetOffset(network->DataSection,p,&offsetd);
902: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetd);
903: *offsetg = offsetp + header->offsetvarrel[compnum];
904: return(0);
905: }
907: /*@
908: DMNetworkGetEdgeOffset - Get the offset for accessing the variable associated with the given edge from the local subvector.
910: Not Collective
912: Input Parameters:
913: + dm - The DMNetwork object
914: - p - the edge point
916: Output Parameters:
917: . offset - the offset
919: Level: intermediate
921: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
922: @*/
923: PetscErrorCode DMNetworkGetEdgeOffset(DM dm,PetscInt p,PetscInt *offset)
924: {
926: DM_Network *network = (DM_Network*)dm->data;
930: PetscSectionGetOffset(network->edge.DofSection,p,offset);
931: return(0);
932: }
934: /*@
935: DMNetworkGetVertexOffset - Get the offset for accessing the variable associated with the given vertex from the local subvector.
937: Not Collective
939: Input Parameters:
940: + dm - The DMNetwork object
941: - p - the vertex point
943: Output Parameters:
944: . offset - the offset
946: Level: intermediate
948: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
949: @*/
950: PetscErrorCode DMNetworkGetVertexOffset(DM dm,PetscInt p,PetscInt *offset)
951: {
953: DM_Network *network = (DM_Network*)dm->data;
957: p -= network->vStart;
959: PetscSectionGetOffset(network->vertex.DofSection,p,offset);
960: return(0);
961: }
962: /*@
963: DMNetworkAddNumVariables - Add number of variables associated with a given point.
965: Not Collective
967: Input Parameters:
968: + dm - The DMNetworkObject
969: . p - the vertex/edge point
970: - nvar - number of additional variables
972: Level: beginner
974: .seealso: DMNetworkSetNumVariables
975: @*/
976: PetscErrorCode DMNetworkAddNumVariables(DM dm,PetscInt p,PetscInt nvar)
977: {
979: DM_Network *network = (DM_Network*)dm->data;
982: PetscSectionAddDof(network->DofSection,p,nvar);
983: return(0);
984: }
986: /*@
987: DMNetworkGetNumVariables - Gets number of variables for a vertex/edge point.
989: Not Collective
991: Input Parameters:
992: + dm - The DMNetworkObject
993: - p - the vertex/edge point
995: Output Parameters:
996: . nvar - number of variables
998: Level: beginner
1000: .seealso: DMNetworkAddNumVariables, DMNetworkSddNumVariables
1001: @*/
1002: PetscErrorCode DMNetworkGetNumVariables(DM dm,PetscInt p,PetscInt *nvar)
1003: {
1005: DM_Network *network = (DM_Network*)dm->data;
1008: PetscSectionGetDof(network->DofSection,p,nvar);
1009: return(0);
1010: }
1012: /*@
1013: DMNetworkSetNumVariables - Sets number of variables for a vertex/edge point.
1015: Not Collective
1017: Input Parameters:
1018: + dm - The DMNetworkObject
1019: . p - the vertex/edge point
1020: - nvar - number of variables
1022: Level: beginner
1024: .seealso: DMNetworkAddNumVariables
1025: @*/
1026: PetscErrorCode DMNetworkSetNumVariables(DM dm,PetscInt p,PetscInt nvar)
1027: {
1029: DM_Network *network = (DM_Network*)dm->data;
1032: PetscSectionSetDof(network->DofSection,p,nvar);
1033: return(0);
1034: }
1036: /* Sets up the array that holds the data for all components and its associated section. This
1037: function is called during DMSetUp() */
1038: PetscErrorCode DMNetworkComponentSetUp(DM dm)
1039: {
1040: PetscErrorCode ierr;
1041: DM_Network *network = (DM_Network*)dm->data;
1042: PetscInt arr_size,p,offset,offsetp,ncomp,i;
1043: DMNetworkComponentHeader header;
1044: DMNetworkComponentValue cvalue;
1045: DMNetworkComponentGenericDataType *componentdataarray;
1048: PetscSectionSetUp(network->DataSection);
1049: PetscSectionGetStorageSize(network->DataSection,&arr_size);
1050: PetscMalloc1(arr_size,&network->componentdataarray);
1051: componentdataarray = network->componentdataarray;
1052: for (p = network->pStart; p < network->pEnd; p++) {
1053: PetscSectionGetOffset(network->DataSection,p,&offsetp);
1054: /* Copy header */
1055: header = &network->header[p];
1056: PetscMemcpy(componentdataarray+offsetp,header,network->dataheadersize*sizeof(DMNetworkComponentGenericDataType));
1057: /* Copy data */
1058: cvalue = &network->cvalue[p];
1059: ncomp = header->ndata;
1060: for (i = 0; i < ncomp; i++) {
1061: offset = offsetp + network->dataheadersize + header->offset[i];
1062: PetscMemcpy(componentdataarray+offset,cvalue->data[i],header->size[i]*sizeof(DMNetworkComponentGenericDataType));
1063: }
1064: }
1065: return(0);
1066: }
1068: /* Sets up the section for dofs. This routine is called during DMSetUp() */
1069: PetscErrorCode DMNetworkVariablesSetUp(DM dm)
1070: {
1072: DM_Network *network = (DM_Network*)dm->data;
1075: PetscSectionSetUp(network->DofSection);
1076: return(0);
1077: }
1079: /*
1080: DMNetworkGetComponentDataArray - Returns the component data array
1082: Not Collective
1084: Input Parameters:
1085: . dm - The DMNetwork Object
1087: Output Parameters:
1088: . componentdataarray - array that holds data for all components
1090: Level: intermediate
1092: .seealso: DMNetworkGetComponentKeyOffset, DMNetworkGetNumComponents
1093: */
1094: PetscErrorCode DMNetworkGetComponentDataArray(DM dm,DMNetworkComponentGenericDataType **componentdataarray)
1095: {
1096: DM_Network *network = (DM_Network*)dm->data;
1099: *componentdataarray = network->componentdataarray;
1100: return(0);
1101: }
1103: /* Get a subsection from a range of points */
1104: PetscErrorCode DMNetworkGetSubSection_private(PetscSection master, PetscInt pstart, PetscInt pend,PetscSection *subsection)
1105: {
1107: PetscInt i, nvar;
1110: PetscSectionCreate(PetscObjectComm((PetscObject)master), subsection);
1111: PetscSectionSetChart(*subsection, 0, pend - pstart);
1112: for (i = pstart; i < pend; i++) {
1113: PetscSectionGetDof(master,i,&nvar);
1114: PetscSectionSetDof(*subsection, i - pstart, nvar);
1115: }
1117: PetscSectionSetUp(*subsection);
1118: return(0);
1119: }
1121: /* Create a submap of points with a GlobalToLocal structure */
1122: PetscErrorCode DMNetworkSetSubMap_private(PetscInt pstart, PetscInt pend, ISLocalToGlobalMapping *map)
1123: {
1125: PetscInt i, *subpoints;
1128: /* Create index sets to map from "points" to "subpoints" */
1129: PetscMalloc1(pend - pstart, &subpoints);
1130: for (i = pstart; i < pend; i++) {
1131: subpoints[i - pstart] = i;
1132: }
1133: ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,pend-pstart,subpoints,PETSC_COPY_VALUES,map);
1134: PetscFree(subpoints);
1135: return(0);
1136: }
1138: /*@
1139: DMNetworkAssembleGraphStructures - Assembles vertex and edge data structures. Must be called after DMNetworkDistribute.
1141: Collective
1143: Input Parameters:
1144: . dm - The DMNetworkObject
1146: Note: the routine will create alternative orderings for the vertices and edges. Assume global network points are:
1148: points = [0 1 2 3 4 5 6]
1150: where edges = [0,1,2,3] and vertices = [4,5,6]. The new orderings will be specific to the subset (i.e vertices = [0,1,2] <- [4,5,6]).
1152: With this new ordering a local PetscSection, global PetscSection and PetscSF will be created specific to the subset.
1154: Level: intermediate
1156: @*/
1157: PetscErrorCode DMNetworkAssembleGraphStructures(DM dm)
1158: {
1160: MPI_Comm comm;
1161: PetscMPIInt rank, size;
1162: DM_Network *network = (DM_Network*)dm->data;
1165: PetscObjectGetComm((PetscObject)dm,&comm);
1166: MPI_Comm_rank(comm, &rank);
1167: MPI_Comm_size(comm, &size);
1169: /* Create maps for vertices and edges */
1170: DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);
1171: DMNetworkSetSubMap_private(network->eStart,network->eEnd,&network->edge.mapping);
1173: /* Create local sub-sections */
1174: DMNetworkGetSubSection_private(network->DofSection,network->vStart,network->vEnd,&network->vertex.DofSection);
1175: DMNetworkGetSubSection_private(network->DofSection,network->eStart,network->eEnd,&network->edge.DofSection);
1177: if (size > 1) {
1178: PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);
1180: PetscSectionCreateGlobalSection(network->vertex.DofSection, network->vertex.sf, PETSC_FALSE, PETSC_FALSE, &network->vertex.GlobalDofSection);
1181: PetscSFGetSubSF(network->plex->sf, network->edge.mapping, &network->edge.sf);
1182: PetscSectionCreateGlobalSection(network->edge.DofSection, network->edge.sf, PETSC_FALSE, PETSC_FALSE, &network->edge.GlobalDofSection);
1183: } else {
1184: /* create structures for vertex */
1185: PetscSectionClone(network->vertex.DofSection,&network->vertex.GlobalDofSection);
1186: /* create structures for edge */
1187: PetscSectionClone(network->edge.DofSection,&network->edge.GlobalDofSection);
1188: }
1190: /* Add viewers */
1191: PetscObjectSetName((PetscObject)network->edge.GlobalDofSection,"Global edge dof section");
1192: PetscObjectSetName((PetscObject)network->vertex.GlobalDofSection,"Global vertex dof section");
1193: PetscSectionViewFromOptions(network->edge.GlobalDofSection, NULL, "-edge_global_section_view");
1194: PetscSectionViewFromOptions(network->vertex.GlobalDofSection, NULL, "-vertex_global_section_view");
1195: return(0);
1196: }
1198: /*@
1199: DMNetworkDistribute - Distributes the network and moves associated component data.
1201: Collective
1203: Input Parameter:
1204: + DM - the DMNetwork object
1205: - overlap - The overlap of partitions, 0 is the default
1207: Notes:
1208: Distributes the network with <overlap>-overlapping partitioning of the edges.
1210: Level: intermediate
1212: .seealso: DMNetworkCreate
1213: @*/
1214: PetscErrorCode DMNetworkDistribute(DM *dm,PetscInt overlap)
1215: {
1216: MPI_Comm comm;
1218: PetscMPIInt size;
1219: DM_Network *oldDMnetwork = (DM_Network*)((*dm)->data);
1220: DM_Network *newDMnetwork;
1221: PetscSF pointsf=NULL;
1222: DM newDM;
1223: PetscInt j,e,v,offset,*subnetvtx;
1224: PetscPartitioner part;
1225: DMNetworkComponentHeader header;
1228: PetscObjectGetComm((PetscObject)*dm,&comm);
1229: MPI_Comm_size(comm, &size);
1230: if (size == 1) return(0);
1232: DMNetworkCreate(PetscObjectComm((PetscObject)*dm),&newDM);
1233: newDMnetwork = (DM_Network*)newDM->data;
1234: newDMnetwork->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
1236: /* Enable runtime options for petscpartitioner */
1237: DMPlexGetPartitioner(oldDMnetwork->plex,&part);
1238: PetscPartitionerSetFromOptions(part);
1240: /* Distribute plex dm and dof section */
1241: DMPlexDistribute(oldDMnetwork->plex,overlap,&pointsf,&newDMnetwork->plex);
1243: /* Distribute dof section */
1244: PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DofSection);
1245: PetscSFDistributeSection(pointsf,oldDMnetwork->DofSection,NULL,newDMnetwork->DofSection);
1246: PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DataSection);
1248: /* Distribute data and associated section */
1249: DMPlexDistributeData(newDMnetwork->plex,pointsf,oldDMnetwork->DataSection,MPIU_INT,(void*)oldDMnetwork->componentdataarray,newDMnetwork->DataSection,(void**)&newDMnetwork->componentdataarray);
1251: PetscSectionGetChart(newDMnetwork->DataSection,&newDMnetwork->pStart,&newDMnetwork->pEnd);
1252: DMPlexGetHeightStratum(newDMnetwork->plex,0, &newDMnetwork->eStart,&newDMnetwork->eEnd);
1253: DMPlexGetHeightStratum(newDMnetwork->plex,1,&newDMnetwork->vStart,&newDMnetwork->vEnd);
1254: newDMnetwork->nEdges = newDMnetwork->eEnd - newDMnetwork->eStart;
1255: newDMnetwork->nVertices = newDMnetwork->vEnd - newDMnetwork->vStart;
1256: newDMnetwork->NVertices = oldDMnetwork->NVertices;
1257: newDMnetwork->NEdges = oldDMnetwork->NEdges;
1259: /* Set Dof section as the section for dm */
1260: DMSetLocalSection(newDMnetwork->plex,newDMnetwork->DofSection);
1261: DMGetGlobalSection(newDMnetwork->plex,&newDMnetwork->GlobalDofSection);
1263: /* Set up subnetwork info in the newDM */
1264: newDMnetwork->nsubnet = oldDMnetwork->nsubnet;
1265: newDMnetwork->ncsubnet = oldDMnetwork->ncsubnet;
1266: PetscCalloc1(newDMnetwork->nsubnet,&newDMnetwork->subnet);
1267: /* Copy over the global number of vertices and edges in each subnetwork. Note that these are already
1268: calculated in DMNetworkLayoutSetUp()
1269: */
1270: for(j=0; j < newDMnetwork->nsubnet; j++) {
1271: newDMnetwork->subnet[j].Nvtx = oldDMnetwork->subnet[j].Nvtx;
1272: newDMnetwork->subnet[j].Nedge = oldDMnetwork->subnet[j].Nedge;
1273: }
1275: for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++ ) {
1276: PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);
1277: header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1278: newDMnetwork->subnet[header->subnetid].nedge++;
1279: }
1281: for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++ ) {
1282: PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);
1283: header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1284: newDMnetwork->subnet[header->subnetid].nvtx++;
1285: }
1287: /* Now create the vertices and edge arrays for the subnetworks */
1288: PetscCalloc1(newDMnetwork->vEnd-newDMnetwork->vStart,&newDMnetwork->subnetvtx);
1289: subnetvtx = newDMnetwork->subnetvtx;
1291: for (j=0; j<newDMnetwork->nsubnet; j++) {
1292: PetscCalloc1(newDMnetwork->subnet[j].nedge,&newDMnetwork->subnet[j].edges);
1293: newDMnetwork->subnet[j].vertices = subnetvtx;
1294: subnetvtx += newDMnetwork->subnet[j].nvtx;
1296: /* Temporarily setting nvtx and nedge to 0 so we can use them as counters in the below for loop.
1297: These get updated when the vertices and edges are added. */
1298: newDMnetwork->subnet[j].nvtx = newDMnetwork->subnet[j].nedge = 0;
1299: }
1301: /* Set the vertices and edges in each subnetwork */
1302: for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++ ) {
1303: PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);
1304: header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1305: newDMnetwork->subnet[header->subnetid].edges[newDMnetwork->subnet[header->subnetid].nedge++] = e;
1306: }
1308: for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++ ) {
1309: PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);
1310: header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1311: newDMnetwork->subnet[header->subnetid].vertices[newDMnetwork->subnet[header->subnetid].nvtx++] = v;
1312: }
1314: newDM->setupcalled = (*dm)->setupcalled;
1315: newDMnetwork->distributecalled = PETSC_TRUE;
1317: /* Destroy point SF */
1318: PetscSFDestroy(&pointsf);
1320: DMDestroy(dm);
1321: *dm = newDM;
1322: return(0);
1323: }
1325: /*@C
1326: PetscSFGetSubSF - Returns an SF for a specific subset of points. Leaves are re-numbered to reflect the new ordering.
1328: Input Parameters:
1329: + masterSF - the original SF structure
1330: - map - a ISLocalToGlobal mapping that contains the subset of points
1332: Output Parameters:
1333: . subSF - a subset of the masterSF for the desired subset.
1334: @*/
1335: PetscErrorCode PetscSFGetSubSF(PetscSF mastersf, ISLocalToGlobalMapping map, PetscSF *subSF) {
1337: PetscErrorCode ierr;
1338: PetscInt nroots, nleaves, *ilocal_sub;
1339: PetscInt i, *ilocal_map, nroots_sub, nleaves_sub = 0;
1340: PetscInt *local_points, *remote_points;
1341: PetscSFNode *iremote_sub;
1342: const PetscInt *ilocal;
1343: const PetscSFNode *iremote;
1346: PetscSFGetGraph(mastersf,&nroots,&nleaves,&ilocal,&iremote);
1348: /* Look for leaves that pertain to the subset of points. Get the local ordering */
1349: PetscMalloc1(nleaves,&ilocal_map);
1350: ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nleaves,ilocal,NULL,ilocal_map);
1351: for (i = 0; i < nleaves; i++) {
1352: if (ilocal_map[i] != -1) nleaves_sub += 1;
1353: }
1354: /* Re-number ilocal with subset numbering. Need information from roots */
1355: PetscMalloc2(nroots,&local_points,nroots,&remote_points);
1356: for (i = 0; i < nroots; i++) local_points[i] = i;
1357: ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nroots,local_points,NULL,local_points);
1358: PetscSFBcastBegin(mastersf, MPIU_INT, local_points, remote_points);
1359: PetscSFBcastEnd(mastersf, MPIU_INT, local_points, remote_points);
1360: /* Fill up graph using local (that is, local to the subset) numbering. */
1361: PetscMalloc1(nleaves_sub,&ilocal_sub);
1362: PetscMalloc1(nleaves_sub,&iremote_sub);
1363: nleaves_sub = 0;
1364: for (i = 0; i < nleaves; i++) {
1365: if (ilocal_map[i] != -1) {
1366: ilocal_sub[nleaves_sub] = ilocal_map[i];
1367: iremote_sub[nleaves_sub].rank = iremote[i].rank;
1368: iremote_sub[nleaves_sub].index = remote_points[ilocal[i]];
1369: nleaves_sub += 1;
1370: }
1371: }
1372: PetscFree2(local_points,remote_points);
1373: ISLocalToGlobalMappingGetSize(map,&nroots_sub);
1375: /* Create new subSF */
1376: PetscSFCreate(PETSC_COMM_WORLD,subSF);
1377: PetscSFSetFromOptions(*subSF);
1378: PetscSFSetGraph(*subSF,nroots_sub,nleaves_sub,ilocal_sub,PETSC_OWN_POINTER,iremote_sub,PETSC_COPY_VALUES);
1379: PetscFree(ilocal_map);
1380: PetscFree(iremote_sub);
1381: return(0);
1382: }
1384: /*@C
1385: DMNetworkGetSupportingEdges - Return the supporting edges for this vertex point
1387: Not Collective
1389: Input Parameters:
1390: + dm - The DMNetwork object
1391: - p - the vertex point
1393: Output Parameters:
1394: + nedges - number of edges connected to this vertex point
1395: - edges - List of edge points
1397: Level: beginner
1399: Fortran Notes:
1400: Since it returns an array, this routine is only available in Fortran 90, and you must
1401: include petsc.h90 in your code.
1403: .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices
1404: @*/
1405: PetscErrorCode DMNetworkGetSupportingEdges(DM dm,PetscInt vertex,PetscInt *nedges,const PetscInt *edges[])
1406: {
1408: DM_Network *network = (DM_Network*)dm->data;
1411: DMPlexGetSupportSize(network->plex,vertex,nedges);
1412: DMPlexGetSupport(network->plex,vertex,edges);
1413: return(0);
1414: }
1416: /*@C
1417: DMNetworkGetConnectedVertices - Return the connected vertices for this edge point
1419: Not Collective
1421: Input Parameters:
1422: + dm - The DMNetwork object
1423: - p - the edge point
1425: Output Parameters:
1426: . vertices - vertices connected to this edge
1428: Level: beginner
1430: Fortran Notes:
1431: Since it returns an array, this routine is only available in Fortran 90, and you must
1432: include petsc.h90 in your code.
1434: .seealso: DMNetworkCreate, DMNetworkGetSupportingEdges
1435: @*/
1436: PetscErrorCode DMNetworkGetConnectedVertices(DM dm,PetscInt edge,const PetscInt *vertices[])
1437: {
1439: DM_Network *network = (DM_Network*)dm->data;
1442: DMPlexGetCone(network->plex,edge,vertices);
1443: return(0);
1444: }
1446: /*@
1447: DMNetworkIsGhostVertex - Returns TRUE if the vertex is a ghost vertex
1449: Not Collective
1451: Input Parameters:
1452: + dm - The DMNetwork object
1453: - p - the vertex point
1455: Output Parameter:
1456: . isghost - TRUE if the vertex is a ghost point
1458: Level: beginner
1460: .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices, DMNetworkGetVertexRange
1461: @*/
1462: PetscErrorCode DMNetworkIsGhostVertex(DM dm,PetscInt p,PetscBool *isghost)
1463: {
1465: DM_Network *network = (DM_Network*)dm->data;
1466: PetscInt offsetg;
1467: PetscSection sectiong;
1470: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
1471: *isghost = PETSC_FALSE;
1472: DMGetGlobalSection(network->plex,§iong);
1473: PetscSectionGetOffset(sectiong,p,&offsetg);
1474: if (offsetg < 0) *isghost = PETSC_TRUE;
1475: return(0);
1476: }
1478: PetscErrorCode DMSetUp_Network(DM dm)
1479: {
1481: DM_Network *network=(DM_Network*)dm->data;
1484: DMNetworkComponentSetUp(dm);
1485: DMNetworkVariablesSetUp(dm);
1487: DMSetLocalSection(network->plex,network->DofSection);
1488: DMGetGlobalSection(network->plex,&network->GlobalDofSection);
1490: dm->setupcalled = PETSC_TRUE;
1491: DMViewFromOptions(dm,NULL,"-dm_view");
1492: return(0);
1493: }
1495: /*@
1496: DMNetworkHasJacobian - Sets global flag for using user's sub Jacobian matrices
1497: -- replaced by DMNetworkSetOption(network,userjacobian,PETSC_TURE)?
1499: Collective
1501: Input Parameters:
1502: + dm - The DMNetwork object
1503: . eflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for edges
1504: - vflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for vertices
1506: Level: intermediate
1508: @*/
1509: PetscErrorCode DMNetworkHasJacobian(DM dm,PetscBool eflg,PetscBool vflg)
1510: {
1511: DM_Network *network=(DM_Network*)dm->data;
1513: PetscInt nVertices = network->nVertices;
1516: network->userEdgeJacobian = eflg;
1517: network->userVertexJacobian = vflg;
1519: if (eflg && !network->Je) {
1520: PetscCalloc1(3*network->nEdges,&network->Je);
1521: }
1523: if (vflg && !network->Jv && nVertices) {
1524: PetscInt i,*vptr,nedges,vStart=network->vStart;
1525: PetscInt nedges_total;
1526: const PetscInt *edges;
1528: /* count nvertex_total */
1529: nedges_total = 0;
1530: PetscMalloc1(nVertices+1,&vptr);
1532: vptr[0] = 0;
1533: for (i=0; i<nVertices; i++) {
1534: DMNetworkGetSupportingEdges(dm,i+vStart,&nedges,&edges);
1535: nedges_total += nedges;
1536: vptr[i+1] = vptr[i] + 2*nedges + 1;
1537: }
1539: PetscCalloc1(2*nedges_total+nVertices,&network->Jv);
1540: network->Jvptr = vptr;
1541: }
1542: return(0);
1543: }
1545: /*@
1546: DMNetworkEdgeSetMatrix - Sets user-provided Jacobian matrices for this edge to the network
1548: Not Collective
1550: Input Parameters:
1551: + dm - The DMNetwork object
1552: . p - the edge point
1553: - J - array (size = 3) of Jacobian submatrices for this edge point:
1554: J[0]: this edge
1555: J[1] and J[2]: connected vertices, obtained by calling DMNetworkGetConnectedVertices()
1557: Level: advanced
1559: .seealso: DMNetworkVertexSetMatrix
1560: @*/
1561: PetscErrorCode DMNetworkEdgeSetMatrix(DM dm,PetscInt p,Mat J[])
1562: {
1563: DM_Network *network=(DM_Network*)dm->data;
1566: if (!network->Je) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkEdgeSetMatrix");
1568: if (J) {
1569: network->Je[3*p] = J[0];
1570: network->Je[3*p+1] = J[1];
1571: network->Je[3*p+2] = J[2];
1572: }
1573: return(0);
1574: }
1576: /*@
1577: DMNetworkVertexSetMatrix - Sets user-provided Jacobian matrix for this vertex to the network
1579: Not Collective
1581: Input Parameters:
1582: + dm - The DMNetwork object
1583: . p - the vertex point
1584: - J - array of Jacobian (size = 2*(num of supporting edges) + 1) submatrices for this vertex point:
1585: J[0]: this vertex
1586: J[1+2*i]: i-th supporting edge
1587: J[1+2*i+1]: i-th connected vertex
1589: Level: advanced
1591: .seealso: DMNetworkEdgeSetMatrix
1592: @*/
1593: PetscErrorCode DMNetworkVertexSetMatrix(DM dm,PetscInt p,Mat J[])
1594: {
1596: DM_Network *network=(DM_Network*)dm->data;
1597: PetscInt i,*vptr,nedges,vStart=network->vStart;
1598: const PetscInt *edges;
1601: if (!network->Jv) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkVertexSetMatrix");
1603: if (J) {
1604: vptr = network->Jvptr;
1605: network->Jv[vptr[p-vStart]] = J[0]; /* Set Jacobian for this vertex */
1607: /* Set Jacobian for each supporting edge and connected vertex */
1608: DMNetworkGetSupportingEdges(dm,p,&nedges,&edges);
1609: for (i=1; i<=2*nedges; i++) network->Jv[vptr[p-vStart]+i] = J[i];
1610: }
1611: return(0);
1612: }
1614: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1615: {
1617: PetscInt j;
1618: PetscScalar val=(PetscScalar)ncols;
1621: if (!ghost) {
1622: for (j=0; j<nrows; j++) {
1623: VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
1624: }
1625: } else {
1626: for (j=0; j<nrows; j++) {
1627: VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
1628: }
1629: }
1630: return(0);
1631: }
1633: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1634: {
1636: PetscInt j,ncols_u;
1637: PetscScalar val;
1640: if (!ghost) {
1641: for (j=0; j<nrows; j++) {
1642: MatGetRow(Ju,j,&ncols_u,NULL,NULL);
1643: val = (PetscScalar)ncols_u;
1644: VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
1645: MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
1646: }
1647: } else {
1648: for (j=0; j<nrows; j++) {
1649: MatGetRow(Ju,j,&ncols_u,NULL,NULL);
1650: val = (PetscScalar)ncols_u;
1651: VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
1652: MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
1653: }
1654: }
1655: return(0);
1656: }
1658: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1659: {
1663: if (Ju) {
1664: MatSetPreallocationUserblock_private(Ju,nrows,rows,ncols,ghost,vdnz,vonz);
1665: } else {
1666: MatSetPreallocationDenseblock_private(nrows,rows,ncols,ghost,vdnz,vonz);
1667: }
1668: return(0);
1669: }
1671: PETSC_STATIC_INLINE PetscErrorCode MatSetDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1672: {
1674: PetscInt j,*cols;
1675: PetscScalar *zeros;
1678: PetscCalloc2(ncols,&cols,nrows*ncols,&zeros);
1679: for (j=0; j<ncols; j++) cols[j] = j+ cstart;
1680: MatSetValues(*J,nrows,rows,ncols,cols,zeros,INSERT_VALUES);
1681: PetscFree2(cols,zeros);
1682: return(0);
1683: }
1685: PETSC_STATIC_INLINE PetscErrorCode MatSetUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1686: {
1688: PetscInt j,M,N,row,col,ncols_u;
1689: const PetscInt *cols;
1690: PetscScalar zero=0.0;
1693: MatGetSize(Ju,&M,&N);
1694: if (nrows != M || ncols != N) SETERRQ4(PetscObjectComm((PetscObject)Ju),PETSC_ERR_USER,"%D by %D must equal %D by %D",nrows,ncols,M,N);
1696: for (row=0; row<nrows; row++) {
1697: MatGetRow(Ju,row,&ncols_u,&cols,NULL);
1698: for (j=0; j<ncols_u; j++) {
1699: col = cols[j] + cstart;
1700: MatSetValues(*J,1,&rows[row],1,&col,&zero,INSERT_VALUES);
1701: }
1702: MatRestoreRow(Ju,row,&ncols_u,&cols,NULL);
1703: }
1704: return(0);
1705: }
1707: PETSC_STATIC_INLINE PetscErrorCode MatSetblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1708: {
1712: if (Ju) {
1713: MatSetUserblock_private(Ju,nrows,rows,ncols,cstart,J);
1714: } else {
1715: MatSetDenseblock_private(nrows,rows,ncols,cstart,J);
1716: }
1717: return(0);
1718: }
1720: /* 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.
1721: */
1722: PetscErrorCode CreateSubGlobalToLocalMapping_private(PetscSection globalsec, PetscSection localsec, ISLocalToGlobalMapping *ltog)
1723: {
1725: PetscInt i,size,dof;
1726: PetscInt *glob2loc;
1729: PetscSectionGetStorageSize(localsec,&size);
1730: PetscMalloc1(size,&glob2loc);
1732: for (i = 0; i < size; i++) {
1733: PetscSectionGetOffset(globalsec,i,&dof);
1734: dof = (dof >= 0) ? dof : -(dof + 1);
1735: glob2loc[i] = dof;
1736: }
1738: ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,size,glob2loc,PETSC_OWN_POINTER,ltog);
1739: #if 0
1740: PetscIntView(size,glob2loc,PETSC_VIEWER_STDOUT_WORLD);
1741: #endif
1742: return(0);
1743: }
1745: #include <petsc/private/matimpl.h>
1747: PetscErrorCode DMCreateMatrix_Network_Nest(DM dm,Mat *J)
1748: {
1750: DM_Network *network = (DM_Network*)dm->data;
1751: PetscMPIInt rank, size;
1752: PetscInt eDof,vDof;
1753: Mat j11,j12,j21,j22,bA[2][2];
1754: MPI_Comm comm;
1755: ISLocalToGlobalMapping eISMap,vISMap;
1758: PetscObjectGetComm((PetscObject)dm,&comm);
1759: MPI_Comm_rank(comm,&rank);
1760: MPI_Comm_size(comm,&size);
1762: PetscSectionGetConstrainedStorageSize(network->edge.GlobalDofSection,&eDof);
1763: PetscSectionGetConstrainedStorageSize(network->vertex.GlobalDofSection,&vDof);
1765: MatCreate(comm, &j11);
1766: MatSetSizes(j11, eDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
1767: MatSetType(j11, MATMPIAIJ);
1769: MatCreate(comm, &j12);
1770: MatSetSizes(j12, eDof, vDof, PETSC_DETERMINE ,PETSC_DETERMINE);
1771: MatSetType(j12, MATMPIAIJ);
1773: MatCreate(comm, &j21);
1774: MatSetSizes(j21, vDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
1775: MatSetType(j21, MATMPIAIJ);
1777: MatCreate(comm, &j22);
1778: MatSetSizes(j22, vDof, vDof, PETSC_DETERMINE, PETSC_DETERMINE);
1779: MatSetType(j22, MATMPIAIJ);
1781: bA[0][0] = j11;
1782: bA[0][1] = j12;
1783: bA[1][0] = j21;
1784: bA[1][1] = j22;
1786: CreateSubGlobalToLocalMapping_private(network->edge.GlobalDofSection,network->edge.DofSection,&eISMap);
1787: CreateSubGlobalToLocalMapping_private(network->vertex.GlobalDofSection,network->vertex.DofSection,&vISMap);
1789: MatSetLocalToGlobalMapping(j11,eISMap,eISMap);
1790: MatSetLocalToGlobalMapping(j12,eISMap,vISMap);
1791: MatSetLocalToGlobalMapping(j21,vISMap,eISMap);
1792: MatSetLocalToGlobalMapping(j22,vISMap,vISMap);
1794: MatSetUp(j11);
1795: MatSetUp(j12);
1796: MatSetUp(j21);
1797: MatSetUp(j22);
1799: MatCreateNest(comm,2,NULL,2,NULL,&bA[0][0],J);
1800: MatSetUp(*J);
1801: MatNestSetVecType(*J,VECNEST);
1802: MatDestroy(&j11);
1803: MatDestroy(&j12);
1804: MatDestroy(&j21);
1805: MatDestroy(&j22);
1807: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1808: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1809: MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1811: /* Free structures */
1812: ISLocalToGlobalMappingDestroy(&eISMap);
1813: ISLocalToGlobalMappingDestroy(&vISMap);
1814: return(0);
1815: }
1817: PetscErrorCode DMCreateMatrix_Network(DM dm,Mat *J)
1818: {
1820: DM_Network *network = (DM_Network*)dm->data;
1821: PetscInt eStart,eEnd,vStart,vEnd,rstart,nrows,*rows,localSize;
1822: PetscInt cstart,ncols,j,e,v;
1823: PetscBool ghost,ghost_vc,ghost2,isNest;
1824: Mat Juser;
1825: PetscSection sectionGlobal;
1826: PetscInt nedges,*vptr=NULL,vc,*rows_v; /* suppress maybe-uninitialized warning */
1827: const PetscInt *edges,*cone;
1828: MPI_Comm comm;
1829: MatType mtype;
1830: Vec vd_nz,vo_nz;
1831: PetscInt *dnnz,*onnz;
1832: PetscScalar *vdnz,*vonz;
1835: mtype = dm->mattype;
1836: PetscStrcmp(mtype,MATNEST,&isNest);
1837: if (isNest) {
1838: DMCreateMatrix_Network_Nest(dm,J);
1839: MatSetDM(*J,dm);
1840: return(0);
1841: }
1843: if (!network->userEdgeJacobian && !network->userVertexJacobian) {
1844: /* user does not provide Jacobian blocks */
1845: DMCreateMatrix_Plex(network->plex,J);
1846: MatSetDM(*J,dm);
1847: return(0);
1848: }
1850: MatCreate(PetscObjectComm((PetscObject)dm),J);
1851: DMGetGlobalSection(network->plex,§ionGlobal);
1852: PetscSectionGetConstrainedStorageSize(sectionGlobal,&localSize);
1853: MatSetSizes(*J,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);
1855: MatSetType(*J,MATAIJ);
1856: MatSetFromOptions(*J);
1858: /* (1) Set matrix preallocation */
1859: /*------------------------------*/
1860: PetscObjectGetComm((PetscObject)dm,&comm);
1861: VecCreate(comm,&vd_nz);
1862: VecSetSizes(vd_nz,localSize,PETSC_DECIDE);
1863: VecSetFromOptions(vd_nz);
1864: VecSet(vd_nz,0.0);
1865: VecDuplicate(vd_nz,&vo_nz);
1867: /* Set preallocation for edges */
1868: /*-----------------------------*/
1869: DMNetworkGetEdgeRange(dm,&eStart,&eEnd);
1871: PetscMalloc1(localSize,&rows);
1872: for (e=eStart; e<eEnd; e++) {
1873: /* Get row indices */
1874: DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
1875: DMNetworkGetNumVariables(dm,e,&nrows);
1876: if (nrows) {
1877: for (j=0; j<nrows; j++) rows[j] = j + rstart;
1879: /* Set preallocation for conntected vertices */
1880: DMNetworkGetConnectedVertices(dm,e,&cone);
1881: for (v=0; v<2; v++) {
1882: DMNetworkGetNumVariables(dm,cone[v],&ncols);
1884: if (network->Je) {
1885: Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
1886: } else Juser = NULL;
1887: DMNetworkIsGhostVertex(dm,cone[v],&ghost);
1888: MatSetPreallocationblock_private(Juser,nrows,rows,ncols,ghost,vd_nz,vo_nz);
1889: }
1891: /* Set preallocation for edge self */
1892: cstart = rstart;
1893: if (network->Je) {
1894: Juser = network->Je[3*e]; /* Jacobian(e,e) */
1895: } else Juser = NULL;
1896: MatSetPreallocationblock_private(Juser,nrows,rows,nrows,PETSC_FALSE,vd_nz,vo_nz);
1897: }
1898: }
1900: /* Set preallocation for vertices */
1901: /*--------------------------------*/
1902: DMNetworkGetVertexRange(dm,&vStart,&vEnd);
1903: if (vEnd - vStart) vptr = network->Jvptr;
1905: for (v=vStart; v<vEnd; v++) {
1906: /* Get row indices */
1907: DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
1908: DMNetworkGetNumVariables(dm,v,&nrows);
1909: if (!nrows) continue;
1911: DMNetworkIsGhostVertex(dm,v,&ghost);
1912: if (ghost) {
1913: PetscMalloc1(nrows,&rows_v);
1914: } else {
1915: rows_v = rows;
1916: }
1918: for (j=0; j<nrows; j++) rows_v[j] = j + rstart;
1920: /* Get supporting edges and connected vertices */
1921: DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
1923: for (e=0; e<nedges; e++) {
1924: /* Supporting edges */
1925: DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
1926: DMNetworkGetNumVariables(dm,edges[e],&ncols);
1928: if (network->Jv) {
1929: Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */
1930: } else Juser = NULL;
1931: MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost,vd_nz,vo_nz);
1933: /* Connected vertices */
1934: DMNetworkGetConnectedVertices(dm,edges[e],&cone);
1935: vc = (v == cone[0]) ? cone[1]:cone[0];
1936: DMNetworkIsGhostVertex(dm,vc,&ghost_vc);
1938: DMNetworkGetNumVariables(dm,vc,&ncols);
1940: if (network->Jv) {
1941: Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
1942: } else Juser = NULL;
1943: if (ghost_vc||ghost) {
1944: ghost2 = PETSC_TRUE;
1945: } else {
1946: ghost2 = PETSC_FALSE;
1947: }
1948: MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost2,vd_nz,vo_nz);
1949: }
1951: /* Set preallocation for vertex self */
1952: DMNetworkIsGhostVertex(dm,v,&ghost);
1953: if (!ghost) {
1954: DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
1955: if (network->Jv) {
1956: Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
1957: } else Juser = NULL;
1958: MatSetPreallocationblock_private(Juser,nrows,rows_v,nrows,PETSC_FALSE,vd_nz,vo_nz);
1959: }
1960: if (ghost) {
1961: PetscFree(rows_v);
1962: }
1963: }
1965: VecAssemblyBegin(vd_nz);
1966: VecAssemblyBegin(vo_nz);
1968: PetscMalloc2(localSize,&dnnz,localSize,&onnz);
1970: VecAssemblyEnd(vd_nz);
1971: VecAssemblyEnd(vo_nz);
1973: VecGetArray(vd_nz,&vdnz);
1974: VecGetArray(vo_nz,&vonz);
1975: for (j=0; j<localSize; j++) {
1976: dnnz[j] = (PetscInt)PetscRealPart(vdnz[j]);
1977: onnz[j] = (PetscInt)PetscRealPart(vonz[j]);
1978: }
1979: VecRestoreArray(vd_nz,&vdnz);
1980: VecRestoreArray(vo_nz,&vonz);
1981: VecDestroy(&vd_nz);
1982: VecDestroy(&vo_nz);
1984: MatSeqAIJSetPreallocation(*J,0,dnnz);
1985: MatMPIAIJSetPreallocation(*J,0,dnnz,0,onnz);
1986: MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1988: PetscFree2(dnnz,onnz);
1990: /* (2) Set matrix entries for edges */
1991: /*----------------------------------*/
1992: for (e=eStart; e<eEnd; e++) {
1993: /* Get row indices */
1994: DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
1995: DMNetworkGetNumVariables(dm,e,&nrows);
1996: if (nrows) {
1997: for (j=0; j<nrows; j++) rows[j] = j + rstart;
1999: /* Set matrix entries for conntected vertices */
2000: DMNetworkGetConnectedVertices(dm,e,&cone);
2001: for (v=0; v<2; v++) {
2002: DMNetworkGetVariableGlobalOffset(dm,cone[v],&cstart);
2003: DMNetworkGetNumVariables(dm,cone[v],&ncols);
2005: if (network->Je) {
2006: Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
2007: } else Juser = NULL;
2008: MatSetblock_private(Juser,nrows,rows,ncols,cstart,J);
2009: }
2011: /* Set matrix entries for edge self */
2012: cstart = rstart;
2013: if (network->Je) {
2014: Juser = network->Je[3*e]; /* Jacobian(e,e) */
2015: } else Juser = NULL;
2016: MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);
2017: }
2018: }
2020: /* Set matrix entries for vertices */
2021: /*---------------------------------*/
2022: for (v=vStart; v<vEnd; v++) {
2023: /* Get row indices */
2024: DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
2025: DMNetworkGetNumVariables(dm,v,&nrows);
2026: if (!nrows) continue;
2028: DMNetworkIsGhostVertex(dm,v,&ghost);
2029: if (ghost) {
2030: PetscMalloc1(nrows,&rows_v);
2031: } else {
2032: rows_v = rows;
2033: }
2034: for (j=0; j<nrows; j++) rows_v[j] = j + rstart;
2036: /* Get supporting edges and connected vertices */
2037: DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
2039: for (e=0; e<nedges; e++) {
2040: /* Supporting edges */
2041: DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
2042: DMNetworkGetNumVariables(dm,edges[e],&ncols);
2044: if (network->Jv) {
2045: Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */
2046: } else Juser = NULL;
2047: MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);
2049: /* Connected vertices */
2050: DMNetworkGetConnectedVertices(dm,edges[e],&cone);
2051: vc = (v == cone[0]) ? cone[1]:cone[0];
2053: DMNetworkGetVariableGlobalOffset(dm,vc,&cstart);
2054: DMNetworkGetNumVariables(dm,vc,&ncols);
2056: if (network->Jv) {
2057: Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
2058: } else Juser = NULL;
2059: MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);
2060: }
2062: /* Set matrix entries for vertex self */
2063: if (!ghost) {
2064: DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
2065: if (network->Jv) {
2066: Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
2067: } else Juser = NULL;
2068: MatSetblock_private(Juser,nrows,rows_v,nrows,cstart,J);
2069: }
2070: if (ghost) {
2071: PetscFree(rows_v);
2072: }
2073: }
2074: PetscFree(rows);
2076: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
2077: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
2079: MatSetDM(*J,dm);
2080: return(0);
2081: }
2083: PetscErrorCode DMDestroy_Network(DM dm)
2084: {
2086: DM_Network *network = (DM_Network*)dm->data;
2087: PetscInt j;
2090: if (--network->refct > 0) return(0);
2091: if (network->Je) {
2092: PetscFree(network->Je);
2093: }
2094: if (network->Jv) {
2095: PetscFree(network->Jvptr);
2096: PetscFree(network->Jv);
2097: }
2099: ISLocalToGlobalMappingDestroy(&network->vertex.mapping);
2100: PetscSectionDestroy(&network->vertex.DofSection);
2101: PetscSectionDestroy(&network->vertex.GlobalDofSection);
2102: if (network->vltog) {
2103: PetscFree(network->vltog);
2104: }
2105: if (network->vertex.sf) {
2106: PetscSFDestroy(&network->vertex.sf);
2107: }
2108: /* edge */
2109: ISLocalToGlobalMappingDestroy(&network->edge.mapping);
2110: PetscSectionDestroy(&network->edge.DofSection);
2111: PetscSectionDestroy(&network->edge.GlobalDofSection);
2112: if (network->edge.sf) {
2113: PetscSFDestroy(&network->edge.sf);
2114: }
2115: DMDestroy(&network->plex);
2116: PetscSectionDestroy(&network->DataSection);
2117: PetscSectionDestroy(&network->DofSection);
2119: for(j=0; j<network->nsubnet; j++) {
2120: PetscFree(network->subnet[j].edges);
2121: }
2122: PetscFree(network->subnetvtx);
2124: PetscFree(network->subnet);
2125: PetscFree(network->componentdataarray);
2126: PetscFree2(network->header,network->cvalue);
2127: PetscFree(network);
2128: return(0);
2129: }
2131: PetscErrorCode DMView_Network(DM dm,PetscViewer viewer)
2132: {
2134: DM_Network *network = (DM_Network*)dm->data;
2135: PetscBool iascii;
2136: PetscMPIInt rank;
2137: PetscInt p,nsubnet;
2140: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
2141: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
2144: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2145: if (iascii) {
2146: const PetscInt *cone,*vtx,*edges;
2147: PetscInt vfrom,vto,i,j,nv,ne;
2149: nsubnet = network->nsubnet - network->ncsubnet; /* num of subnetworks */
2150: PetscViewerASCIIPushSynchronized(viewer);
2151: PetscViewerASCIISynchronizedPrintf(viewer, " [%d] nsubnet: %D; nsubnetCouple: %D; nEdges: %D; nVertices: %D\n",rank,nsubnet,network->ncsubnet,network->nEdges,network->nVertices);
2153: for (i=0; i<nsubnet; i++) {
2154: DMNetworkGetSubnetworkInfo(dm,i,&nv,&ne,&vtx,&edges);
2155: if (ne) {
2156: PetscViewerASCIISynchronizedPrintf(viewer, " Subnet %D: nEdges %D, nVertices %D\n",i,ne,nv);
2157: for (j=0; j<ne; j++) {
2158: p = edges[j];
2159: DMNetworkGetConnectedVertices(dm,p,&cone);
2160: DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);
2161: DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);
2162: DMNetworkGetGlobalEdgeIndex(dm,edges[j],&p);
2163: PetscViewerASCIISynchronizedPrintf(viewer, " edge %D: %D----> %D\n",p,vfrom,vto);
2164: }
2165: }
2166: }
2167: /* Coupling subnets */
2168: nsubnet = network->nsubnet;
2169: for (; i<nsubnet; i++) {
2170: DMNetworkGetSubnetworkInfo(dm,i,&nv,&ne,&vtx,&edges);
2171: if (ne) {
2172: PetscViewerASCIISynchronizedPrintf(viewer, " Subnet %D (couple): nEdges %D, nVertices %D\n",i,ne,nv);
2173: for (j=0; j<ne; j++) {
2174: p = edges[j];
2175: DMNetworkGetConnectedVertices(dm,p,&cone);
2176: DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);
2177: DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);
2178: PetscViewerASCIISynchronizedPrintf(viewer, " edge %D: %D----> %D\n",p,vfrom,vto);
2179: }
2180: }
2181: }
2182: PetscViewerFlush(viewer);
2183: PetscViewerASCIIPopSynchronized(viewer);
2184: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMNetwork writing", ((PetscObject)viewer)->type_name);
2185: return(0);
2186: }
2188: PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l)
2189: {
2191: DM_Network *network = (DM_Network*)dm->data;
2194: DMGlobalToLocalBegin(network->plex,g,mode,l);
2195: return(0);
2196: }
2198: PetscErrorCode DMGlobalToLocalEnd_Network(DM dm, Vec g, InsertMode mode, Vec l)
2199: {
2201: DM_Network *network = (DM_Network*)dm->data;
2204: DMGlobalToLocalEnd(network->plex,g,mode,l);
2205: return(0);
2206: }
2208: PetscErrorCode DMLocalToGlobalBegin_Network(DM dm, Vec l, InsertMode mode, Vec g)
2209: {
2211: DM_Network *network = (DM_Network*)dm->data;
2214: DMLocalToGlobalBegin(network->plex,l,mode,g);
2215: return(0);
2216: }
2218: PetscErrorCode DMLocalToGlobalEnd_Network(DM dm, Vec l, InsertMode mode, Vec g)
2219: {
2221: DM_Network *network = (DM_Network*)dm->data;
2224: DMLocalToGlobalEnd(network->plex,l,mode,g);
2225: return(0);
2226: }
2228: /*@
2229: DMNetworkGetVertexLocalToGlobalOrdering - Get vertex global index
2231: Not collective
2233: Input Parameters:
2234: + dm - the dm object
2235: - vloc - local vertex ordering, start from 0
2237: Output Parameters:
2238: . vg - global vertex ordering, start from 0
2240: Level: advanced
2242: .seealso: DMNetworkSetVertexLocalToGlobalOrdering()
2243: @*/
2244: PetscErrorCode DMNetworkGetVertexLocalToGlobalOrdering(DM dm,PetscInt vloc,PetscInt *vg)
2245: {
2246: DM_Network *network = (DM_Network*)dm->data;
2247: PetscInt *vltog = network->vltog;
2250: if (!vltog) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkSetVertexLocalToGlobalOrdering() first");
2251: *vg = vltog[vloc];
2252: return(0);
2253: }
2255: /*@
2256: DMNetworkSetVertexLocalToGlobalOrdering - Create and setup vertex local to global map
2258: Collective
2260: Input Parameters:
2261: . dm - the dm object
2263: Level: advanced
2265: .seealso: DMNetworkGetGlobalVertexIndex()
2266: @*/
2267: PetscErrorCode DMNetworkSetVertexLocalToGlobalOrdering(DM dm)
2268: {
2269: PetscErrorCode ierr;
2270: DM_Network *network=(DM_Network*)dm->data;
2271: MPI_Comm comm;
2272: PetscMPIInt rank,size,*displs,*recvcounts,remoterank;
2273: PetscBool ghost;
2274: PetscInt *vltog,nroots,nleaves,i,*vrange,k,N,lidx;
2275: const PetscSFNode *iremote;
2276: PetscSF vsf;
2277: Vec Vleaves,Vleaves_seq;
2278: VecScatter ctx;
2279: PetscScalar *varr,val;
2280: const PetscScalar *varr_read;
2283: PetscObjectGetComm((PetscObject)dm,&comm);
2284: MPI_Comm_size(comm,&size);
2285: MPI_Comm_rank(comm,&rank);
2287: if (size == 1) {
2288: nroots = network->vEnd - network->vStart;
2289: PetscMalloc1(nroots, &vltog);
2290: for (i=0; i<nroots; i++) vltog[i] = i;
2291: network->vltog = vltog;
2292: return(0);
2293: }
2295: if (!network->distributecalled) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkDistribute() first");
2296: if (network->vltog) {
2297: PetscFree(network->vltog);
2298: }
2300: DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);
2301: PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);
2302: vsf = network->vertex.sf;
2304: PetscMalloc3(size+1,&vrange,size+1,&displs,size,&recvcounts);
2305: PetscSFGetGraph(vsf,&nroots,&nleaves,NULL,&iremote);
2307: for (i=0; i<size; i++) { displs[i] = i; recvcounts[i] = 1;}
2309: i = nroots - nleaves; /* local number of vertices, excluding ghosts */
2310: vrange[0] = 0;
2311: MPI_Allgatherv(&i,1,MPIU_INT,vrange+1,recvcounts,displs,MPIU_INT,comm);
2312: for (i=2; i<size+1; i++) {vrange[i] += vrange[i-1];}
2314: PetscMalloc1(nroots, &vltog);
2315: network->vltog = vltog;
2317: /* Set vltog for non-ghost vertices */
2318: k = 0;
2319: for (i=0; i<nroots; i++) {
2320: DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);
2321: if (ghost) continue;
2322: vltog[i] = vrange[rank] + k++;
2323: }
2324: PetscFree3(vrange,displs,recvcounts);
2326: /* Set vltog for ghost vertices */
2327: /* (a) create parallel Vleaves and sequential Vleaves_seq to convert local iremote[*].index to global index */
2328: VecCreate(comm,&Vleaves);
2329: VecSetSizes(Vleaves,2*nleaves,PETSC_DETERMINE);
2330: VecSetFromOptions(Vleaves);
2331: VecGetArray(Vleaves,&varr);
2332: for (i=0; i<nleaves; i++) {
2333: varr[2*i] = (PetscScalar)(iremote[i].rank); /* rank of remote process */
2334: varr[2*i+1] = (PetscScalar)(iremote[i].index); /* local index in remote process */
2335: }
2336: VecRestoreArray(Vleaves,&varr);
2338: /* (b) scatter local info to remote processes via VecScatter() */
2339: VecScatterCreateToAll(Vleaves,&ctx,&Vleaves_seq);
2340: VecScatterBegin(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);
2341: VecScatterEnd(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);
2343: /* (c) convert local indices to global indices in parallel vector Vleaves */
2344: VecGetSize(Vleaves_seq,&N);
2345: VecGetArrayRead(Vleaves_seq,&varr_read);
2346: for (i=0; i<N; i+=2) {
2347: remoterank = (PetscMPIInt)PetscRealPart(varr_read[i]);
2348: if (remoterank == rank) {
2349: k = i+1; /* row number */
2350: lidx = (PetscInt)PetscRealPart(varr_read[i+1]);
2351: val = (PetscScalar)vltog[lidx]; /* global index for non-ghost vertex computed above */
2352: VecSetValues(Vleaves,1,&k,&val,INSERT_VALUES);
2353: }
2354: }
2355: VecRestoreArrayRead(Vleaves_seq,&varr_read);
2356: VecAssemblyBegin(Vleaves);
2357: VecAssemblyEnd(Vleaves);
2359: /* (d) Set vltog for ghost vertices by copying local values of Vleaves */
2360: VecGetArrayRead(Vleaves,&varr_read);
2361: k = 0;
2362: for (i=0; i<nroots; i++) {
2363: DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);
2364: if (!ghost) continue;
2365: vltog[i] = (PetscInt)PetscRealPart(varr_read[2*k+1]); k++;
2366: }
2367: VecRestoreArrayRead(Vleaves,&varr_read);
2369: VecDestroy(&Vleaves);
2370: VecDestroy(&Vleaves_seq);
2371: VecScatterDestroy(&ctx);
2372: return(0);
2373: }