Actual source code: network.c
petsc-3.12.5 2020-03-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: DMNetworkSetSizes - Sets the number of subnetworks, local and global vertices and edges for each subnetwork.
28: Collective on dm
30: Input Parameters:
31: + dm - the dm object
32: . Nsubnet - global number of subnetworks
33: . nV - number of local vertices for each subnetwork
34: . nE - number of local edges for each subnetwork
35: . NsubnetCouple - global number of coupling subnetworks
36: - nec - number of local edges for each coupling subnetwork
38: You cannot change the sizes once they have been set. nV, nE are arrays of length Nsubnet, and nec is array of length NsubnetCouple.
40: Level: intermediate
42: .seealso: DMNetworkCreate()
43: @*/
44: PetscErrorCode DMNetworkSetSizes(DM dm,PetscInt Nsubnet,PetscInt nV[], PetscInt nE[],PetscInt NsubnetCouple,PetscInt nec[])
45: {
47: DM_Network *network = (DM_Network*) dm->data;
48: PetscInt a[2],b[2],i;
52: if (Nsubnet <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subnetworks %D cannot be less than 1",Nsubnet);
53: if (NsubnetCouple < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of coupling subnetworks %D cannot be less than 0",NsubnetCouple);
57: if (network->nsubnet != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Network sizes alread set, cannot resize the network");
59: if (!nV || !nE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Local vertex size or edge size must be provided");
61: network->nsubnet = Nsubnet + NsubnetCouple;
62: network->ncsubnet = NsubnetCouple;
63: PetscCalloc1(Nsubnet+NsubnetCouple,&network->subnet);
65: /* ----------------------------------------------------------
66: p=v or e; P=V or E
67: subnet[0].pStart = 0
68: subnet[i+1].pStart = subnet[i].pEnd = subnet[i].pStart + (nE[i] or NV[i])
69: ----------------------------------------------------------------------- */
70: for (i=0; i < Nsubnet; i++) {
71: /* Get global number of vertices and edges for subnet[i] */
72: a[0] = nV[i]; a[1] = nE[i]; /* local number of vertices (excluding ghost) and edges */
73: MPIU_Allreduce(a,b,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
74: network->subnet[i].Nvtx = b[0];
75: network->subnet[i].Nedge = b[1];
77: network->subnet[i].nvtx = nV[i]; /* local nvtx, without ghost */
79: /* global subnet[].vStart and vEnd, used by DMNetworkLayoutSetUp() */
80: network->subnet[i].vStart = network->NVertices;
81: network->subnet[i].vEnd = network->subnet[i].vStart + network->subnet[i].Nvtx;
83: network->nVertices += nV[i];
84: network->NVertices += network->subnet[i].Nvtx;
86: network->subnet[i].nedge = nE[i];
87: network->subnet[i].eStart = network->nEdges;
88: network->subnet[i].eEnd = network->subnet[i].eStart + nE[i];
89: network->nEdges += nE[i];
90: network->NEdges += network->subnet[i].Nedge;
91: }
93: /* coupling subnetwork */
94: for (; i < Nsubnet+NsubnetCouple; i++) {
95: /* Get global number of coupling edges for subnet[i] */
96: MPIU_Allreduce(nec+(i-Nsubnet),&network->subnet[i].Nedge,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
98: network->subnet[i].nvtx = 0; /* We design coupling subnetwork such that it does not have its own vertices */
99: network->subnet[i].vStart = network->nVertices;
100: network->subnet[i].vEnd = network->subnet[i].vStart;
102: network->subnet[i].nedge = nec[i-Nsubnet];
103: network->subnet[i].eStart = network->nEdges;
104: network->subnet[i].eEnd = network->subnet[i].eStart + nec[i-Nsubnet];
105: network->nEdges += nec[i-Nsubnet];
106: network->NEdges += network->subnet[i].Nedge;
107: }
108: return(0);
109: }
111: /*@
112: DMNetworkSetEdgeList - Sets the list of local edges (vertex connectivity) for the network
114: Logically collective on dm
116: Input Parameters:
117: + dm - the dm object
118: . edgelist - list of edges for each subnetwork
119: - edgelistCouple - list of edges for each coupling subnetwork
121: Notes:
122: There is no copy involved in this operation, only the pointer is referenced. The edgelist should
123: not be destroyed before the call to DMNetworkLayoutSetUp
125: Level: intermediate
127: Example usage:
128: Consider the following 2 separate networks and a coupling network:
130: .vb
131: network 0: v0 -> v1 -> v2 -> v3
132: network 1: v1 -> v2 -> v0
133: coupling network: network 1: v2 -> network 0: v0
134: .ve
136: The resulting input
137: edgelist[0] = [0 1 | 1 2 | 2 3];
138: edgelist[1] = [1 2 | 2 0]
139: edgelistCouple[0] = [(network)1 (v)2 (network)0 (v)0].
141: .seealso: DMNetworkCreate, DMNetworkSetSizes
142: @*/
143: PetscErrorCode DMNetworkSetEdgeList(DM dm,PetscInt *edgelist[],PetscInt *edgelistCouple[])
144: {
145: DM_Network *network = (DM_Network*) dm->data;
146: PetscInt i,nsubnet,ncsubnet=network->ncsubnet;
149: nsubnet = network->nsubnet - ncsubnet;
150: for(i=0; i < nsubnet; i++) {
151: network->subnet[i].edgelist = edgelist[i];
152: }
153: if (edgelistCouple) {
154: PetscInt j;
155: j = 0;
156: nsubnet = network->nsubnet;
157: while (i < nsubnet) network->subnet[i++].edgelist = edgelistCouple[j++];
158: }
159: return(0);
160: }
162: /*@
163: DMNetworkLayoutSetUp - Sets up the bare layout (graph) for the network
165: Collective on dm
167: Input Parameters
168: . DM - the dmnetwork object
170: Notes:
171: This routine should be called after the network sizes and edgelists have been provided. It creates
172: the bare layout of the network and sets up the network to begin insertion of components.
174: All the components should be registered before calling this routine.
176: Level: intermediate
178: .seealso: DMNetworkSetSizes, DMNetworkSetEdgeList
179: @*/
180: PetscErrorCode DMNetworkLayoutSetUp(DM dm)
181: {
183: DM_Network *network = (DM_Network*)dm->data;
184: PetscInt numCorners=2,spacedim=2,dim = 1; /* One dimensional network */
185: PetscReal *vertexcoords=NULL;
186: PetscInt i,j,ctr,nsubnet,*eowners,np,*edges,*subnetvtx,vStart;
187: PetscInt k,netid,vid, *vidxlTog,*edgelist_couple=NULL;
188: const PetscInt *cone;
189: MPI_Comm comm;
190: PetscMPIInt size,rank;
193: PetscObjectGetComm((PetscObject)dm,&comm);
194: MPI_Comm_rank(comm,&rank);
195: MPI_Comm_size(comm,&size);
197: /* Create the local edgelist for the network by concatenating local input edgelists of the subnetworks */
198: PetscCalloc2(numCorners*network->nVertices,&vertexcoords,2*network->nEdges,&edges);
199: nsubnet = network->nsubnet - network->ncsubnet;
200: ctr = 0;
201: for (i=0; i < nsubnet; i++) {
202: for (j = 0; j < network->subnet[i].nedge; j++) {
203: edges[2*ctr] = network->subnet[i].vStart + network->subnet[i].edgelist[2*j];
204: edges[2*ctr+1] = network->subnet[i].vStart + network->subnet[i].edgelist[2*j+1];
205: ctr++;
206: }
207: }
209: /* Append local coupling edgelists of the subnetworks */
210: i = nsubnet; /* netid of coupling subnet */
211: nsubnet = network->nsubnet;
212: while (i < nsubnet) {
213: edgelist_couple = network->subnet[i].edgelist;
214: k = 0;
215: for (j = 0; j < network->subnet[i].nedge; j++) {
216: netid = edgelist_couple[k]; vid = edgelist_couple[k+1];
217: edges[2*ctr] = network->subnet[netid].vStart + vid; k += 2;
219: netid = edgelist_couple[k]; vid = edgelist_couple[k+1];
220: edges[2*ctr+1] = network->subnet[netid].vStart + vid; k+=2;
221: ctr++;
222: }
223: i++;
224: }
225: /*
226: if (rank == 0) {
227: PetscPrintf(PETSC_COMM_SELF,"[%d] edgelist:\n",rank);
228: for(i=0; i < network->nEdges; i++) {
229: PetscPrintf(PETSC_COMM_SELF,"[%D %D]",edges[2*i],edges[2*i+1]);
230: printf("\n");
231: }
232: }
233: */
235: /* Create network->plex */
236: #if defined(PETSC_USE_64BIT_INDICES)
237: {
238: int *edges64;
239: np = network->nEdges*numCorners;
240: PetscMalloc1(np,&edges64);
241: for (i=0; i<np; i++) edges64[i] = (int)edges[i];
243: if (size == 1) {
244: DMPlexCreateFromCellList(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges64,spacedim,(const double*)vertexcoords,&network->plex);
245: } else {
246: DMPlexCreateFromCellListParallel(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges64,spacedim,(const PetscReal*)vertexcoords,NULL,&network->plex);
247: }
248: PetscFree(edges64);
249: }
250: #else
251: if (size == 1) {
252: DMPlexCreateFromCellList(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges,spacedim,(const double*)vertexcoords,&network->plex);
253: } else {
254: DMPlexCreateFromCellListParallel(comm,dim,network->nEdges,network->nVertices,numCorners,PETSC_FALSE,(const int*)edges,spacedim,(const PetscReal*)vertexcoords,NULL,&network->plex);
255: }
256: #endif
258: DMPlexGetChart(network->plex,&network->pStart,&network->pEnd);
259: DMPlexGetHeightStratum(network->plex,0,&network->eStart,&network->eEnd);
260: DMPlexGetHeightStratum(network->plex,1,&network->vStart,&network->vEnd);
261: vStart = network->vStart;
263: PetscSectionCreate(comm,&network->DataSection);
264: PetscSectionCreate(comm,&network->DofSection);
265: PetscSectionSetChart(network->DataSection,network->pStart,network->pEnd);
266: PetscSectionSetChart(network->DofSection,network->pStart,network->pEnd);
268: network->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
269: np = network->pEnd - network->pStart;
270: PetscCalloc2(np,&network->header,np,&network->cvalue);
272: /* Create vidxlTog: maps local vertex index to global index */
273: np = network->vEnd - vStart;
274: PetscMalloc2(np,&vidxlTog,size+1,&eowners);
275: ctr = 0;
276: for (i=network->eStart; i<network->eEnd; i++) {
277: DMNetworkGetConnectedVertices(dm,i,&cone);
278: vidxlTog[cone[0] - vStart] = edges[2*ctr];
279: vidxlTog[cone[1] - vStart] = edges[2*ctr+1];
280: ctr++;
281: }
282: PetscFree2(vertexcoords,edges);
284: /* Create vertices and edges array for the subnetworks */
285: for (j=0; j < network->nsubnet; j++) {
286: PetscCalloc1(network->subnet[j].nedge,&network->subnet[j].edges);
288: /* Temporarily setting nvtx and nedge to 0 so we can use them as counters in the below for loop.
289: These get updated when the vertices and edges are added. */
290: network->subnet[j].nvtx = 0;
291: network->subnet[j].nedge = 0;
292: }
293: PetscCalloc1(np,&network->subnetvtx);
296: /* Get edge ownership */
297: np = network->eEnd - network->eStart;
298: MPI_Allgather(&np,1,MPIU_INT,eowners+1,1,MPIU_INT,comm);
299: eowners[0] = 0;
300: for (i=2; i<=size; i++) eowners[i] += eowners[i-1];
302: i = 0; j = 0;
303: while (i < np) { /* local edges, including coupling edges */
304: network->header[i].index = i + eowners[rank]; /* Global edge index */
306: if (j < network->nsubnet && i < network->subnet[j].eEnd) {
307: network->header[i].subnetid = j; /* Subnetwork id */
308: network->subnet[j].edges[network->subnet[j].nedge++] = i;
310: network->header[i].ndata = 0;
311: PetscSectionAddDof(network->DataSection,i,network->dataheadersize);
312: network->header[i].offset[0] = 0;
313: i++;
314: }
315: if (i >= network->subnet[j].eEnd) j++;
316: }
318: /* Count network->subnet[*].nvtx */
319: for (i=vStart; i<network->vEnd; i++) { /* local vertices, including ghosts */
320: k = vidxlTog[i-vStart];
321: for (j=0; j < network->nsubnet; j++) {
322: if (network->subnet[j].vStart <= k && k < network->subnet[j].vEnd) {
323: network->subnet[j].nvtx++;
324: break;
325: }
326: }
327: }
329: /* Set network->subnet[*].vertices on array network->subnetvtx */
330: subnetvtx = network->subnetvtx;
331: for (j=0; j<network->nsubnet; j++) {
332: network->subnet[j].vertices = subnetvtx;
333: subnetvtx += network->subnet[j].nvtx;
334: network->subnet[j].nvtx = 0;
335: }
337: /* Set vertex array for the subnetworks */
338: for (i=vStart; i<network->vEnd; i++) { /* local vertices, including ghosts */
339: network->header[i].index = vidxlTog[i-vStart]; /* Global vertex index */
341: k = vidxlTog[i-vStart];
342: for (j=0; j < network->nsubnet; j++) {
343: if (network->subnet[j].vStart <= k && k < network->subnet[j].vEnd) {
344: network->header[i].subnetid = j;
345: network->subnet[j].vertices[network->subnet[j].nvtx++] = i;
346: break;
347: }
348: }
350: network->header[i].ndata = 0;
351: PetscSectionAddDof(network->DataSection,i,network->dataheadersize);
352: network->header[i].offset[0] = 0;
353: }
355: PetscFree2(vidxlTog,eowners);
356: return(0);
357: }
359: /*@C
360: DMNetworkGetSubnetworkInfo - Returns the info for the subnetwork
362: Input Parameters
363: + dm - the DM object
364: - id - the ID (integer) of the subnetwork
366: Output Parameters
367: + nv - number of vertices (local)
368: . ne - number of edges (local)
369: . vtx - local vertices for this subnetwork
370: - edge - local edges for this subnetwork
372: Notes:
373: Cannot call this routine before DMNetworkLayoutSetup()
375: Level: intermediate
377: .seealso: DMNetworkLayoutSetUp, DMNetworkCreate
378: @*/
379: PetscErrorCode DMNetworkGetSubnetworkInfo(DM dm,PetscInt id,PetscInt *nv, PetscInt *ne,const PetscInt **vtx, const PetscInt **edge)
380: {
381: DM_Network *network = (DM_Network*)dm->data;
384: *nv = network->subnet[id].nvtx;
385: *ne = network->subnet[id].nedge;
386: *vtx = network->subnet[id].vertices;
387: *edge = network->subnet[id].edges;
388: return(0);
389: }
391: /*@C
392: DMNetworkGetSubnetworkCoupleInfo - Returns the info for the coupling subnetwork
394: Input Parameters
395: + dm - the DM object
396: - id - the ID (integer) of the coupling subnetwork
398: Output Parameters
399: + ne - number of edges (local)
400: - edge - local edges for this coupling subnetwork
402: Notes:
403: Cannot call this routine before DMNetworkLayoutSetup()
405: Level: intermediate
407: .seealso: DMNetworkGetSubnetworkInfo, DMNetworkLayoutSetUp, DMNetworkCreate
408: @*/
409: PetscErrorCode DMNetworkGetSubnetworkCoupleInfo(DM dm,PetscInt id,PetscInt *ne,const PetscInt **edge)
410: {
411: DM_Network *net = (DM_Network*)dm->data;
412: PetscInt id1 = id + net->nsubnet - net->ncsubnet;
415: *ne = net->subnet[id1].nedge;
416: *edge = net->subnet[id1].edges;
417: return(0);
418: }
420: /*@C
421: DMNetworkRegisterComponent - Registers the network component
423: Logically collective on dm
425: Input Parameters
426: + dm - the network object
427: . name - the component name
428: - size - the storage size in bytes for this component data
430: Output Parameters
431: . key - an integer key that defines the component
433: Notes
434: This routine should be called by all processors before calling DMNetworkLayoutSetup().
436: Level: intermediate
438: .seealso: DMNetworkLayoutSetUp, DMNetworkCreate
439: @*/
440: PetscErrorCode DMNetworkRegisterComponent(DM dm,const char *name,size_t size,PetscInt *key)
441: {
442: PetscErrorCode ierr;
443: DM_Network *network = (DM_Network*) dm->data;
444: DMNetworkComponent *component=&network->component[network->ncomponent];
445: PetscBool flg=PETSC_FALSE;
446: PetscInt i;
449: for (i=0; i < network->ncomponent; i++) {
450: PetscStrcmp(component->name,name,&flg);
451: if (flg) {
452: *key = i;
453: return(0);
454: }
455: }
456: if(network->ncomponent == MAX_COMPONENTS) {
457: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Number of components registered exceeds the max %D",MAX_COMPONENTS);
458: }
460: PetscStrcpy(component->name,name);
461: component->size = size/sizeof(DMNetworkComponentGenericDataType);
462: *key = network->ncomponent;
463: network->ncomponent++;
464: return(0);
465: }
467: /*@
468: DMNetworkGetVertexRange - Get the bounds [start, end) for the vertices.
470: Not Collective
472: Input Parameters:
473: . dm - The DMNetwork object
475: Output Paramters:
476: + vStart - The first vertex point
477: - vEnd - One beyond the last vertex point
479: Level: intermediate
481: .seealso: DMNetworkGetEdgeRange
482: @*/
483: PetscErrorCode DMNetworkGetVertexRange(DM dm,PetscInt *vStart,PetscInt *vEnd)
484: {
485: DM_Network *network = (DM_Network*)dm->data;
488: if (vStart) *vStart = network->vStart;
489: if (vEnd) *vEnd = network->vEnd;
490: return(0);
491: }
493: /*@
494: DMNetworkGetEdgeRange - Get the bounds [start, end) for the edges.
496: Not Collective
498: Input Parameters:
499: . dm - The DMNetwork object
501: Output Paramters:
502: + eStart - The first edge point
503: - eEnd - One beyond the last edge point
505: Level: intermediate
507: .seealso: DMNetworkGetVertexRange
508: @*/
509: PetscErrorCode DMNetworkGetEdgeRange(DM dm,PetscInt *eStart,PetscInt *eEnd)
510: {
511: DM_Network *network = (DM_Network*)dm->data;
514: if (eStart) *eStart = network->eStart;
515: if (eEnd) *eEnd = network->eEnd;
516: return(0);
517: }
519: /*@
520: DMNetworkGetGlobalEdgeIndex - Get the user global numbering for the edge.
522: Not Collective
524: Input Parameters:
525: + dm - DMNetwork object
526: - p - edge point
528: Output Paramters:
529: . index - user global numbering for the edge
531: Level: intermediate
533: .seealso: DMNetworkGetGlobalVertexIndex
534: @*/
535: PetscErrorCode DMNetworkGetGlobalEdgeIndex(DM dm,PetscInt p,PetscInt *index)
536: {
537: PetscErrorCode ierr;
538: DM_Network *network = (DM_Network*)dm->data;
539: PetscInt offsetp;
540: DMNetworkComponentHeader header;
543: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
544: PetscSectionGetOffset(network->DataSection,p,&offsetp);
545: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
546: *index = header->index;
547: return(0);
548: }
550: /*@
551: DMNetworkGetGlobalVertexIndex - Get the user global numbering for the vertex.
553: Not Collective
555: Input Parameters:
556: + dm - DMNetwork object
557: - p - vertex point
559: Output Paramters:
560: . index - user global numbering for the vertex
562: Level: intermediate
564: .seealso: DMNetworkGetGlobalEdgeIndex
565: @*/
566: PetscErrorCode DMNetworkGetGlobalVertexIndex(DM dm,PetscInt p,PetscInt *index)
567: {
568: PetscErrorCode ierr;
569: DM_Network *network = (DM_Network*)dm->data;
570: PetscInt offsetp;
571: DMNetworkComponentHeader header;
574: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
575: PetscSectionGetOffset(network->DataSection,p,&offsetp);
576: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
577: *index = header->index;
578: return(0);
579: }
581: /*
582: DMNetworkGetComponentKeyOffset - Gets the type along with the offset for indexing the
583: component value from the component data array
585: Not Collective
587: Input Parameters:
588: + dm - The DMNetwork object
589: . p - vertex/edge point
590: - compnum - component number
592: Output Parameters:
593: + compkey - the key obtained when registering the component
594: - offset - offset into the component data array associated with the vertex/edge point
596: Notes:
597: Typical usage:
599: DMNetworkGetComponentDataArray(dm, &arr);
600: DMNetworkGetVertex/EdgeRange(dm,&Start,&End);
601: Loop over vertices or edges
602: DMNetworkGetNumComponents(dm,v,&numcomps);
603: Loop over numcomps
604: DMNetworkGetComponentKeyOffset(dm,v,compnum,&key,&offset);
605: compdata = (UserCompDataType)(arr+offset);
607: Level: intermediate
609: .seealso: DMNetworkGetNumComponents, DMNetworkGetComponentDataArray,
610: */
611: PetscErrorCode DMNetworkGetComponentKeyOffset(DM dm,PetscInt p, PetscInt compnum, PetscInt *compkey, PetscInt *offset)
612: {
613: PetscErrorCode ierr;
614: PetscInt offsetp;
615: DMNetworkComponentHeader header;
616: DM_Network *network = (DM_Network*)dm->data;
619: PetscSectionGetOffset(network->DataSection,p,&offsetp);
620: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
621: if (compkey) *compkey = header->key[compnum];
622: if (offset) *offset = offsetp+network->dataheadersize+header->offset[compnum];
623: return(0);
624: }
626: /*@
627: DMNetworkGetComponent - Returns the network component and its key
629: Not Collective
631: Input Parameters
632: + dm - DMNetwork object
633: . p - edge or vertex point
634: - compnum - component number
636: Output Parameters:
637: + compkey - the key set for this computing during registration
638: - component - the component data
640: Notes:
641: Typical usage:
643: DMNetworkGetVertex/EdgeRange(dm,&Start,&End);
644: Loop over vertices or edges
645: DMNetworkGetNumComponents(dm,v,&numcomps);
646: Loop over numcomps
647: DMNetworkGetComponent(dm,v,compnum,&key,&component);
649: Level: intermediate
651: .seealso: DMNetworkGetNumComponents, DMNetworkGetVariableOffset
652: @*/
653: PetscErrorCode DMNetworkGetComponent(DM dm, PetscInt p, PetscInt compnum, PetscInt *key, void **component)
654: {
656: DM_Network *network = (DM_Network*)dm->data;
657: PetscInt offsetd = 0;
660: DMNetworkGetComponentKeyOffset(dm,p,compnum,key,&offsetd);
661: *component = network->componentdataarray+offsetd;
662: return(0);
663: }
665: /*@
666: DMNetworkAddComponent - Adds a network component at the given point (vertex/edge)
668: Not Collective
670: Input Parameters:
671: + dm - The DMNetwork object
672: . p - vertex/edge point
673: . componentkey - component key returned while registering the component
674: - compvalue - pointer to the data structure for the component
676: Level: intermediate
678: .seealso: DMNetworkGetVertexRange, DMNetworkGetEdgeRange, DMNetworkRegisterComponent
679: @*/
680: PetscErrorCode DMNetworkAddComponent(DM dm, PetscInt p,PetscInt componentkey,void* compvalue)
681: {
682: DM_Network *network = (DM_Network*)dm->data;
683: DMNetworkComponent *component = &network->component[componentkey];
684: DMNetworkComponentHeader header = &network->header[p];
685: DMNetworkComponentValue cvalue = &network->cvalue[p];
686: PetscErrorCode ierr;
689: 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);
691: header->size[header->ndata] = component->size;
692: PetscSectionAddDof(network->DataSection,p,component->size);
693: header->key[header->ndata] = componentkey;
694: if (header->ndata != 0) header->offset[header->ndata] = header->offset[header->ndata-1] + header->size[header->ndata-1];
696: cvalue->data[header->ndata] = (void*)compvalue;
697: header->ndata++;
698: return(0);
699: }
701: /*@
702: DMNetworkGetNumComponents - Get the number of components at a vertex/edge
704: Not Collective
706: Input Parameters:
707: + dm - The DMNetwork object
708: - p - vertex/edge point
710: Output Parameters:
711: . numcomponents - Number of components at the vertex/edge
713: Level: intermediate
715: .seealso: DMNetworkRegisterComponent, DMNetworkAddComponent
716: @*/
717: PetscErrorCode DMNetworkGetNumComponents(DM dm,PetscInt p,PetscInt *numcomponents)
718: {
720: PetscInt offset;
721: DM_Network *network = (DM_Network*)dm->data;
724: PetscSectionGetOffset(network->DataSection,p,&offset);
725: *numcomponents = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata;
726: return(0);
727: }
729: /*@
730: DMNetworkGetVariableOffset - Get the offset for accessing the variable associated with the given vertex/edge from the local vector.
732: Not Collective
734: Input Parameters:
735: + dm - The DMNetwork object
736: - p - the edge/vertex point
738: Output Parameters:
739: . offset - the offset
741: Level: intermediate
743: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
744: @*/
745: PetscErrorCode DMNetworkGetVariableOffset(DM dm,PetscInt p,PetscInt *offset)
746: {
748: DM_Network *network = (DM_Network*)dm->data;
751: PetscSectionGetOffset(network->plex->localSection,p,offset);
752: return(0);
753: }
755: /*@
756: DMNetworkGetVariableGlobalOffset - Get the global offset for the variable associated with the given vertex/edge from the global vector.
758: Not Collective
760: Input Parameters:
761: + dm - The DMNetwork object
762: - p - the edge/vertex point
764: Output Parameters:
765: . offsetg - the offset
767: Level: intermediate
769: .seealso: DMNetworkGetVariableOffset, DMGetLocalVector
770: @*/
771: PetscErrorCode DMNetworkGetVariableGlobalOffset(DM dm,PetscInt p,PetscInt *offsetg)
772: {
774: DM_Network *network = (DM_Network*)dm->data;
777: PetscSectionGetOffset(network->plex->globalSection,p,offsetg);
778: if (*offsetg < 0) *offsetg = -(*offsetg + 1); /* Convert to actual global offset for ghost vertex */
779: return(0);
780: }
782: /*@
783: DMNetworkGetEdgeOffset - Get the offset for accessing the variable associated with the given edge from the local subvector.
785: Not Collective
787: Input Parameters:
788: + dm - The DMNetwork object
789: - p - the edge point
791: Output Parameters:
792: . offset - the offset
794: Level: intermediate
796: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
797: @*/
798: PetscErrorCode DMNetworkGetEdgeOffset(DM dm,PetscInt p,PetscInt *offset)
799: {
801: DM_Network *network = (DM_Network*)dm->data;
805: PetscSectionGetOffset(network->edge.DofSection,p,offset);
806: return(0);
807: }
809: /*@
810: DMNetworkGetVertexOffset - Get the offset for accessing the variable associated with the given vertex from the local subvector.
812: Not Collective
814: Input Parameters:
815: + dm - The DMNetwork object
816: - p - the vertex point
818: Output Parameters:
819: . offset - the offset
821: Level: intermediate
823: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
824: @*/
825: PetscErrorCode DMNetworkGetVertexOffset(DM dm,PetscInt p,PetscInt *offset)
826: {
828: DM_Network *network = (DM_Network*)dm->data;
832: p -= network->vStart;
834: PetscSectionGetOffset(network->vertex.DofSection,p,offset);
835: return(0);
836: }
837: /*@
838: DMNetworkAddNumVariables - Add number of variables associated with a given point.
840: Not Collective
842: Input Parameters:
843: + dm - The DMNetworkObject
844: . p - the vertex/edge point
845: - nvar - number of additional variables
847: Level: intermediate
849: .seealso: DMNetworkSetNumVariables
850: @*/
851: PetscErrorCode DMNetworkAddNumVariables(DM dm,PetscInt p,PetscInt nvar)
852: {
854: DM_Network *network = (DM_Network*)dm->data;
857: PetscSectionAddDof(network->DofSection,p,nvar);
858: return(0);
859: }
861: /*@
862: DMNetworkGetNumVariables - Gets number of variables for a vertex/edge point.
864: Not Collective
866: Input Parameters:
867: + dm - The DMNetworkObject
868: - p - the vertex/edge point
870: Output Parameters:
871: . nvar - number of variables
873: Level: intermediate
875: .seealso: DMNetworkAddNumVariables, DMNetworkSddNumVariables
876: @*/
877: PetscErrorCode DMNetworkGetNumVariables(DM dm,PetscInt p,PetscInt *nvar)
878: {
880: DM_Network *network = (DM_Network*)dm->data;
883: PetscSectionGetDof(network->DofSection,p,nvar);
884: return(0);
885: }
887: /*@
888: DMNetworkSetNumVariables - Sets number of variables for a vertex/edge point.
890: Not Collective
892: Input Parameters:
893: + dm - The DMNetworkObject
894: . p - the vertex/edge point
895: - nvar - number of variables
897: Level: intermediate
899: .seealso: DMNetworkAddNumVariables
900: @*/
901: PetscErrorCode DMNetworkSetNumVariables(DM dm,PetscInt p,PetscInt nvar)
902: {
904: DM_Network *network = (DM_Network*)dm->data;
907: PetscSectionSetDof(network->DofSection,p,nvar);
908: return(0);
909: }
911: /* Sets up the array that holds the data for all components and its associated section. This
912: function is called during DMSetUp() */
913: PetscErrorCode DMNetworkComponentSetUp(DM dm)
914: {
915: PetscErrorCode ierr;
916: DM_Network *network = (DM_Network*)dm->data;
917: PetscInt arr_size,p,offset,offsetp,ncomp,i;
918: DMNetworkComponentHeader header;
919: DMNetworkComponentValue cvalue;
920: DMNetworkComponentGenericDataType *componentdataarray;
923: PetscSectionSetUp(network->DataSection);
924: PetscSectionGetStorageSize(network->DataSection,&arr_size);
925: PetscMalloc1(arr_size,&network->componentdataarray);
926: componentdataarray = network->componentdataarray;
927: for (p = network->pStart; p < network->pEnd; p++) {
928: PetscSectionGetOffset(network->DataSection,p,&offsetp);
929: /* Copy header */
930: header = &network->header[p];
931: PetscMemcpy(componentdataarray+offsetp,header,network->dataheadersize*sizeof(DMNetworkComponentGenericDataType));
932: /* Copy data */
933: cvalue = &network->cvalue[p];
934: ncomp = header->ndata;
935: for (i = 0; i < ncomp; i++) {
936: offset = offsetp + network->dataheadersize + header->offset[i];
937: PetscMemcpy(componentdataarray+offset,cvalue->data[i],header->size[i]*sizeof(DMNetworkComponentGenericDataType));
938: }
939: }
940: return(0);
941: }
943: /* Sets up the section for dofs. This routine is called during DMSetUp() */
944: PetscErrorCode DMNetworkVariablesSetUp(DM dm)
945: {
947: DM_Network *network = (DM_Network*)dm->data;
950: PetscSectionSetUp(network->DofSection);
951: return(0);
952: }
954: /*@C
955: DMNetworkGetComponentDataArray - Returns the component data array
957: Not Collective
959: Input Parameters:
960: . dm - The DMNetwork Object
962: Output Parameters:
963: . componentdataarray - array that holds data for all components
965: Level: intermediate
967: .seealso: DMNetworkGetComponentKeyOffset, DMNetworkGetNumComponents
968: @*/
969: PetscErrorCode DMNetworkGetComponentDataArray(DM dm,DMNetworkComponentGenericDataType **componentdataarray)
970: {
971: DM_Network *network = (DM_Network*)dm->data;
974: *componentdataarray = network->componentdataarray;
975: return(0);
976: }
978: /* Get a subsection from a range of points */
979: PetscErrorCode DMNetworkGetSubSection_private(PetscSection master, PetscInt pstart, PetscInt pend,PetscSection *subsection)
980: {
982: PetscInt i, nvar;
985: PetscSectionCreate(PetscObjectComm((PetscObject)master), subsection);
986: PetscSectionSetChart(*subsection, 0, pend - pstart);
987: for (i = pstart; i < pend; i++) {
988: PetscSectionGetDof(master,i,&nvar);
989: PetscSectionSetDof(*subsection, i - pstart, nvar);
990: }
992: PetscSectionSetUp(*subsection);
993: return(0);
994: }
996: /* Create a submap of points with a GlobalToLocal structure */
997: PetscErrorCode DMNetworkSetSubMap_private(PetscInt pstart, PetscInt pend, ISLocalToGlobalMapping *map)
998: {
1000: PetscInt i, *subpoints;
1003: /* Create index sets to map from "points" to "subpoints" */
1004: PetscMalloc1(pend - pstart, &subpoints);
1005: for (i = pstart; i < pend; i++) {
1006: subpoints[i - pstart] = i;
1007: }
1008: ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,pend-pstart,subpoints,PETSC_COPY_VALUES,map);
1009: PetscFree(subpoints);
1010: return(0);
1011: }
1013: /*@
1014: DMNetworkAssembleGraphStructures - Assembles vertex and edge data structures. Must be called after DMNetworkDistribute.
1016: Collective
1018: Input Parameters:
1019: . dm - The DMNetworkObject
1021: Note: the routine will create alternative orderings for the vertices and edges. Assume global network points are:
1023: points = [0 1 2 3 4 5 6]
1025: 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]).
1027: With this new ordering a local PetscSection, global PetscSection and PetscSF will be created specific to the subset.
1029: Level: intermediate
1031: @*/
1032: PetscErrorCode DMNetworkAssembleGraphStructures(DM dm)
1033: {
1035: MPI_Comm comm;
1036: PetscMPIInt rank, size;
1037: DM_Network *network = (DM_Network*)dm->data;
1040: PetscObjectGetComm((PetscObject)dm,&comm);
1041: MPI_Comm_rank(comm, &rank);
1042: MPI_Comm_size(comm, &size);
1044: /* Create maps for vertices and edges */
1045: DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);
1046: DMNetworkSetSubMap_private(network->eStart,network->eEnd,&network->edge.mapping);
1048: /* Create local sub-sections */
1049: DMNetworkGetSubSection_private(network->DofSection,network->vStart,network->vEnd,&network->vertex.DofSection);
1050: DMNetworkGetSubSection_private(network->DofSection,network->eStart,network->eEnd,&network->edge.DofSection);
1052: if (size > 1) {
1053: PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);
1055: PetscSectionCreateGlobalSection(network->vertex.DofSection, network->vertex.sf, PETSC_FALSE, PETSC_FALSE, &network->vertex.GlobalDofSection);
1056: PetscSFGetSubSF(network->plex->sf, network->edge.mapping, &network->edge.sf);
1057: PetscSectionCreateGlobalSection(network->edge.DofSection, network->edge.sf, PETSC_FALSE, PETSC_FALSE, &network->edge.GlobalDofSection);
1058: } else {
1059: /* create structures for vertex */
1060: PetscSectionClone(network->vertex.DofSection,&network->vertex.GlobalDofSection);
1061: /* create structures for edge */
1062: PetscSectionClone(network->edge.DofSection,&network->edge.GlobalDofSection);
1063: }
1065: /* Add viewers */
1066: PetscObjectSetName((PetscObject)network->edge.GlobalDofSection,"Global edge dof section");
1067: PetscObjectSetName((PetscObject)network->vertex.GlobalDofSection,"Global vertex dof section");
1068: PetscSectionViewFromOptions(network->edge.GlobalDofSection, NULL, "-edge_global_section_view");
1069: PetscSectionViewFromOptions(network->vertex.GlobalDofSection, NULL, "-vertex_global_section_view");
1070: return(0);
1071: }
1073: /*@
1074: DMNetworkDistribute - Distributes the network and moves associated component data.
1076: Collective
1078: Input Parameter:
1079: + DM - the DMNetwork object
1080: - overlap - The overlap of partitions, 0 is the default
1082: Notes:
1083: Distributes the network with <overlap>-overlapping partitioning of the edges.
1085: Level: intermediate
1087: .seealso: DMNetworkCreate
1088: @*/
1089: PetscErrorCode DMNetworkDistribute(DM *dm,PetscInt overlap)
1090: {
1091: MPI_Comm comm;
1093: PetscMPIInt size;
1094: DM_Network *oldDMnetwork = (DM_Network*)((*dm)->data);
1095: DM_Network *newDMnetwork;
1096: PetscSF pointsf=NULL;
1097: DM newDM;
1098: PetscInt j,e,v,offset,*subnetvtx;
1099: PetscPartitioner part;
1100: DMNetworkComponentHeader header;
1103: PetscObjectGetComm((PetscObject)*dm,&comm);
1104: MPI_Comm_size(comm, &size);
1105: if (size == 1) return(0);
1107: DMNetworkCreate(PetscObjectComm((PetscObject)*dm),&newDM);
1108: newDMnetwork = (DM_Network*)newDM->data;
1109: newDMnetwork->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
1111: /* Enable runtime options for petscpartitioner */
1112: DMPlexGetPartitioner(oldDMnetwork->plex,&part);
1113: PetscPartitionerSetFromOptions(part);
1115: /* Distribute plex dm and dof section */
1116: DMPlexDistribute(oldDMnetwork->plex,overlap,&pointsf,&newDMnetwork->plex);
1118: /* Distribute dof section */
1119: PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DofSection);
1120: PetscSFDistributeSection(pointsf,oldDMnetwork->DofSection,NULL,newDMnetwork->DofSection);
1121: PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DataSection);
1123: /* Distribute data and associated section */
1124: DMPlexDistributeData(newDMnetwork->plex,pointsf,oldDMnetwork->DataSection,MPIU_INT,(void*)oldDMnetwork->componentdataarray,newDMnetwork->DataSection,(void**)&newDMnetwork->componentdataarray);
1126: PetscSectionGetChart(newDMnetwork->DataSection,&newDMnetwork->pStart,&newDMnetwork->pEnd);
1127: DMPlexGetHeightStratum(newDMnetwork->plex,0, &newDMnetwork->eStart,&newDMnetwork->eEnd);
1128: DMPlexGetHeightStratum(newDMnetwork->plex,1,&newDMnetwork->vStart,&newDMnetwork->vEnd);
1129: newDMnetwork->nEdges = newDMnetwork->eEnd - newDMnetwork->eStart;
1130: newDMnetwork->nVertices = newDMnetwork->vEnd - newDMnetwork->vStart;
1131: newDMnetwork->NVertices = oldDMnetwork->NVertices;
1132: newDMnetwork->NEdges = oldDMnetwork->NEdges;
1134: /* Set Dof section as the section for dm */
1135: DMSetLocalSection(newDMnetwork->plex,newDMnetwork->DofSection);
1136: DMGetGlobalSection(newDMnetwork->plex,&newDMnetwork->GlobalDofSection);
1138: /* Set up subnetwork info in the newDM */
1139: newDMnetwork->nsubnet = oldDMnetwork->nsubnet;
1140: newDMnetwork->ncsubnet = oldDMnetwork->ncsubnet;
1141: PetscCalloc1(newDMnetwork->nsubnet,&newDMnetwork->subnet);
1142: /* Copy over the global number of vertices and edges in each subnetwork. Note that these are already
1143: calculated in DMNetworkLayoutSetUp()
1144: */
1145: for(j=0; j < newDMnetwork->nsubnet; j++) {
1146: newDMnetwork->subnet[j].Nvtx = oldDMnetwork->subnet[j].Nvtx;
1147: newDMnetwork->subnet[j].Nedge = oldDMnetwork->subnet[j].Nedge;
1148: }
1150: for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++ ) {
1151: PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);
1152: header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1153: newDMnetwork->subnet[header->subnetid].nedge++;
1154: }
1156: for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++ ) {
1157: PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);
1158: header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1159: newDMnetwork->subnet[header->subnetid].nvtx++;
1160: }
1162: /* Now create the vertices and edge arrays for the subnetworks */
1163: PetscCalloc1(newDMnetwork->vEnd-newDMnetwork->vStart,&newDMnetwork->subnetvtx);
1164: subnetvtx = newDMnetwork->subnetvtx;
1166: for (j=0; j<newDMnetwork->nsubnet; j++) {
1167: PetscCalloc1(newDMnetwork->subnet[j].nedge,&newDMnetwork->subnet[j].edges);
1168: newDMnetwork->subnet[j].vertices = subnetvtx;
1169: subnetvtx += newDMnetwork->subnet[j].nvtx;
1171: /* Temporarily setting nvtx and nedge to 0 so we can use them as counters in the below for loop.
1172: These get updated when the vertices and edges are added. */
1173: newDMnetwork->subnet[j].nvtx = newDMnetwork->subnet[j].nedge = 0;
1174: }
1176: /* Set the vertices and edges in each subnetwork */
1177: for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++ ) {
1178: PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);
1179: header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1180: newDMnetwork->subnet[header->subnetid].edges[newDMnetwork->subnet[header->subnetid].nedge++] = e;
1181: }
1183: for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++ ) {
1184: PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);
1185: header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1186: newDMnetwork->subnet[header->subnetid].vertices[newDMnetwork->subnet[header->subnetid].nvtx++] = v;
1187: }
1189: newDM->setupcalled = (*dm)->setupcalled;
1190: newDMnetwork->distributecalled = PETSC_TRUE;
1192: /* Destroy point SF */
1193: PetscSFDestroy(&pointsf);
1195: DMDestroy(dm);
1196: *dm = newDM;
1197: return(0);
1198: }
1200: /*@C
1201: PetscSFGetSubSF - Returns an SF for a specific subset of points. Leaves are re-numbered to reflect the new ordering.
1203: Input Parameters:
1204: + masterSF - the original SF structure
1205: - map - a ISLocalToGlobal mapping that contains the subset of points
1207: Output Parameters:
1208: . subSF - a subset of the masterSF for the desired subset.
1209: */
1211: PetscErrorCode PetscSFGetSubSF(PetscSF mastersf, ISLocalToGlobalMapping map, PetscSF *subSF) {
1213: PetscErrorCode ierr;
1214: PetscInt nroots, nleaves, *ilocal_sub;
1215: PetscInt i, *ilocal_map, nroots_sub, nleaves_sub = 0;
1216: PetscInt *local_points, *remote_points;
1217: PetscSFNode *iremote_sub;
1218: const PetscInt *ilocal;
1219: const PetscSFNode *iremote;
1222: PetscSFGetGraph(mastersf,&nroots,&nleaves,&ilocal,&iremote);
1224: /* Look for leaves that pertain to the subset of points. Get the local ordering */
1225: PetscMalloc1(nleaves,&ilocal_map);
1226: ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nleaves,ilocal,NULL,ilocal_map);
1227: for (i = 0; i < nleaves; i++) {
1228: if (ilocal_map[i] != -1) nleaves_sub += 1;
1229: }
1230: /* Re-number ilocal with subset numbering. Need information from roots */
1231: PetscMalloc2(nroots,&local_points,nroots,&remote_points);
1232: for (i = 0; i < nroots; i++) local_points[i] = i;
1233: ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nroots,local_points,NULL,local_points);
1234: PetscSFBcastBegin(mastersf, MPIU_INT, local_points, remote_points);
1235: PetscSFBcastEnd(mastersf, MPIU_INT, local_points, remote_points);
1236: /* Fill up graph using local (that is, local to the subset) numbering. */
1237: PetscMalloc1(nleaves_sub,&ilocal_sub);
1238: PetscMalloc1(nleaves_sub,&iremote_sub);
1239: nleaves_sub = 0;
1240: for (i = 0; i < nleaves; i++) {
1241: if (ilocal_map[i] != -1) {
1242: ilocal_sub[nleaves_sub] = ilocal_map[i];
1243: iremote_sub[nleaves_sub].rank = iremote[i].rank;
1244: iremote_sub[nleaves_sub].index = remote_points[ilocal[i]];
1245: nleaves_sub += 1;
1246: }
1247: }
1248: PetscFree2(local_points,remote_points);
1249: ISLocalToGlobalMappingGetSize(map,&nroots_sub);
1251: /* Create new subSF */
1252: PetscSFCreate(PETSC_COMM_WORLD,subSF);
1253: PetscSFSetFromOptions(*subSF);
1254: PetscSFSetGraph(*subSF,nroots_sub,nleaves_sub,ilocal_sub,PETSC_OWN_POINTER,iremote_sub,PETSC_COPY_VALUES);
1255: PetscFree(ilocal_map);
1256: PetscFree(iremote_sub);
1257: return(0);
1258: }
1260: /*@C
1261: DMNetworkGetSupportingEdges - Return the supporting edges for this vertex point
1263: Not Collective
1265: Input Parameters:
1266: + dm - The DMNetwork object
1267: - p - the vertex point
1269: Output Paramters:
1270: + nedges - number of edges connected to this vertex point
1271: - edges - List of edge points
1273: Level: intermediate
1275: Fortran Notes:
1276: Since it returns an array, this routine is only available in Fortran 90, and you must
1277: include petsc.h90 in your code.
1279: .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices
1280: @*/
1281: PetscErrorCode DMNetworkGetSupportingEdges(DM dm,PetscInt vertex,PetscInt *nedges,const PetscInt *edges[])
1282: {
1284: DM_Network *network = (DM_Network*)dm->data;
1287: DMPlexGetSupportSize(network->plex,vertex,nedges);
1288: DMPlexGetSupport(network->plex,vertex,edges);
1289: return(0);
1290: }
1292: /*@C
1293: DMNetworkGetConnectedVertices - Return the connected vertices for this edge point
1295: Not Collective
1297: Input Parameters:
1298: + dm - The DMNetwork object
1299: - p - the edge point
1301: Output Paramters:
1302: . vertices - vertices connected to this edge
1304: Level: intermediate
1306: Fortran Notes:
1307: Since it returns an array, this routine is only available in Fortran 90, and you must
1308: include petsc.h90 in your code.
1310: .seealso: DMNetworkCreate, DMNetworkGetSupportingEdges
1311: @*/
1312: PetscErrorCode DMNetworkGetConnectedVertices(DM dm,PetscInt edge,const PetscInt *vertices[])
1313: {
1315: DM_Network *network = (DM_Network*)dm->data;
1318: DMPlexGetCone(network->plex,edge,vertices);
1319: return(0);
1320: }
1322: /*@
1323: DMNetworkIsGhostVertex - Returns TRUE if the vertex is a ghost vertex
1325: Not Collective
1327: Input Parameters:
1328: + dm - The DMNetwork object
1329: - p - the vertex point
1331: Output Parameter:
1332: . isghost - TRUE if the vertex is a ghost point
1334: Level: intermediate
1336: .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices, DMNetworkGetVertexRange
1337: @*/
1338: PetscErrorCode DMNetworkIsGhostVertex(DM dm,PetscInt p,PetscBool *isghost)
1339: {
1341: DM_Network *network = (DM_Network*)dm->data;
1342: PetscInt offsetg;
1343: PetscSection sectiong;
1346: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
1347: *isghost = PETSC_FALSE;
1348: DMGetGlobalSection(network->plex,§iong);
1349: PetscSectionGetOffset(sectiong,p,&offsetg);
1350: if (offsetg < 0) *isghost = PETSC_TRUE;
1351: return(0);
1352: }
1354: PetscErrorCode DMSetUp_Network(DM dm)
1355: {
1357: DM_Network *network=(DM_Network*)dm->data;
1360: DMNetworkComponentSetUp(dm);
1361: DMNetworkVariablesSetUp(dm);
1363: DMSetLocalSection(network->plex,network->DofSection);
1364: DMGetGlobalSection(network->plex,&network->GlobalDofSection);
1366: dm->setupcalled = PETSC_TRUE;
1367: DMViewFromOptions(dm,NULL,"-dm_view");
1368: return(0);
1369: }
1371: /*@
1372: DMNetworkHasJacobian - Sets global flag for using user's sub Jacobian matrices
1373: -- replaced by DMNetworkSetOption(network,userjacobian,PETSC_TURE)?
1375: Collective
1377: Input Parameters:
1378: + dm - The DMNetwork object
1379: . eflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for edges
1380: - vflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for vertices
1382: Level: intermediate
1384: @*/
1385: PetscErrorCode DMNetworkHasJacobian(DM dm,PetscBool eflg,PetscBool vflg)
1386: {
1387: DM_Network *network=(DM_Network*)dm->data;
1389: PetscInt nVertices = network->nVertices;
1392: network->userEdgeJacobian = eflg;
1393: network->userVertexJacobian = vflg;
1395: if (eflg && !network->Je) {
1396: PetscCalloc1(3*network->nEdges,&network->Je);
1397: }
1399: if (vflg && !network->Jv && nVertices) {
1400: PetscInt i,*vptr,nedges,vStart=network->vStart;
1401: PetscInt nedges_total;
1402: const PetscInt *edges;
1404: /* count nvertex_total */
1405: nedges_total = 0;
1406: PetscMalloc1(nVertices+1,&vptr);
1408: vptr[0] = 0;
1409: for (i=0; i<nVertices; i++) {
1410: DMNetworkGetSupportingEdges(dm,i+vStart,&nedges,&edges);
1411: nedges_total += nedges;
1412: vptr[i+1] = vptr[i] + 2*nedges + 1;
1413: }
1415: PetscCalloc1(2*nedges_total+nVertices,&network->Jv);
1416: network->Jvptr = vptr;
1417: }
1418: return(0);
1419: }
1421: /*@
1422: DMNetworkEdgeSetMatrix - Sets user-provided Jacobian matrices for this edge to the network
1424: Not Collective
1426: Input Parameters:
1427: + dm - The DMNetwork object
1428: . p - the edge point
1429: - J - array (size = 3) of Jacobian submatrices for this edge point:
1430: J[0]: this edge
1431: J[1] and J[2]: connected vertices, obtained by calling DMNetworkGetConnectedVertices()
1433: Level: intermediate
1435: .seealso: DMNetworkVertexSetMatrix
1436: @*/
1437: PetscErrorCode DMNetworkEdgeSetMatrix(DM dm,PetscInt p,Mat J[])
1438: {
1439: DM_Network *network=(DM_Network*)dm->data;
1442: if (!network->Je) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkEdgeSetMatrix");
1444: if (J) {
1445: network->Je[3*p] = J[0];
1446: network->Je[3*p+1] = J[1];
1447: network->Je[3*p+2] = J[2];
1448: }
1449: return(0);
1450: }
1452: /*@
1453: DMNetworkVertexSetMatrix - Sets user-provided Jacobian matrix for this vertex to the network
1455: Not Collective
1457: Input Parameters:
1458: + dm - The DMNetwork object
1459: . p - the vertex point
1460: - J - array of Jacobian (size = 2*(num of supporting edges) + 1) submatrices for this vertex point:
1461: J[0]: this vertex
1462: J[1+2*i]: i-th supporting edge
1463: J[1+2*i+1]: i-th connected vertex
1465: Level: intermediate
1467: .seealso: DMNetworkEdgeSetMatrix
1468: @*/
1469: PetscErrorCode DMNetworkVertexSetMatrix(DM dm,PetscInt p,Mat J[])
1470: {
1472: DM_Network *network=(DM_Network*)dm->data;
1473: PetscInt i,*vptr,nedges,vStart=network->vStart;
1474: const PetscInt *edges;
1477: if (!network->Jv) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkVertexSetMatrix");
1479: if (J) {
1480: vptr = network->Jvptr;
1481: network->Jv[vptr[p-vStart]] = J[0]; /* Set Jacobian for this vertex */
1483: /* Set Jacobian for each supporting edge and connected vertex */
1484: DMNetworkGetSupportingEdges(dm,p,&nedges,&edges);
1485: for (i=1; i<=2*nedges; i++) network->Jv[vptr[p-vStart]+i] = J[i];
1486: }
1487: return(0);
1488: }
1490: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1491: {
1493: PetscInt j;
1494: PetscScalar val=(PetscScalar)ncols;
1497: if (!ghost) {
1498: for (j=0; j<nrows; j++) {
1499: VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
1500: }
1501: } else {
1502: for (j=0; j<nrows; j++) {
1503: VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
1504: }
1505: }
1506: return(0);
1507: }
1509: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1510: {
1512: PetscInt j,ncols_u;
1513: PetscScalar val;
1516: if (!ghost) {
1517: for (j=0; j<nrows; j++) {
1518: MatGetRow(Ju,j,&ncols_u,NULL,NULL);
1519: val = (PetscScalar)ncols_u;
1520: VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
1521: MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
1522: }
1523: } else {
1524: for (j=0; j<nrows; j++) {
1525: MatGetRow(Ju,j,&ncols_u,NULL,NULL);
1526: val = (PetscScalar)ncols_u;
1527: VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
1528: MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
1529: }
1530: }
1531: return(0);
1532: }
1534: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1535: {
1539: if (Ju) {
1540: MatSetPreallocationUserblock_private(Ju,nrows,rows,ncols,ghost,vdnz,vonz);
1541: } else {
1542: MatSetPreallocationDenseblock_private(nrows,rows,ncols,ghost,vdnz,vonz);
1543: }
1544: return(0);
1545: }
1547: PETSC_STATIC_INLINE PetscErrorCode MatSetDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1548: {
1550: PetscInt j,*cols;
1551: PetscScalar *zeros;
1554: PetscCalloc2(ncols,&cols,nrows*ncols,&zeros);
1555: for (j=0; j<ncols; j++) cols[j] = j+ cstart;
1556: MatSetValues(*J,nrows,rows,ncols,cols,zeros,INSERT_VALUES);
1557: PetscFree2(cols,zeros);
1558: return(0);
1559: }
1561: PETSC_STATIC_INLINE PetscErrorCode MatSetUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1562: {
1564: PetscInt j,M,N,row,col,ncols_u;
1565: const PetscInt *cols;
1566: PetscScalar zero=0.0;
1569: MatGetSize(Ju,&M,&N);
1570: if (nrows != M || ncols != N) SETERRQ4(PetscObjectComm((PetscObject)Ju),PETSC_ERR_USER,"%D by %D must equal %D by %D",nrows,ncols,M,N);
1572: for (row=0; row<nrows; row++) {
1573: MatGetRow(Ju,row,&ncols_u,&cols,NULL);
1574: for (j=0; j<ncols_u; j++) {
1575: col = cols[j] + cstart;
1576: MatSetValues(*J,1,&rows[row],1,&col,&zero,INSERT_VALUES);
1577: }
1578: MatRestoreRow(Ju,row,&ncols_u,&cols,NULL);
1579: }
1580: return(0);
1581: }
1583: PETSC_STATIC_INLINE PetscErrorCode MatSetblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1584: {
1588: if (Ju) {
1589: MatSetUserblock_private(Ju,nrows,rows,ncols,cstart,J);
1590: } else {
1591: MatSetDenseblock_private(nrows,rows,ncols,cstart,J);
1592: }
1593: return(0);
1594: }
1596: /* 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.
1597: */
1598: PetscErrorCode CreateSubGlobalToLocalMapping_private(PetscSection globalsec, PetscSection localsec, ISLocalToGlobalMapping *ltog)
1599: {
1601: PetscInt i,size,dof;
1602: PetscInt *glob2loc;
1605: PetscSectionGetStorageSize(localsec,&size);
1606: PetscMalloc1(size,&glob2loc);
1608: for (i = 0; i < size; i++) {
1609: PetscSectionGetOffset(globalsec,i,&dof);
1610: dof = (dof >= 0) ? dof : -(dof + 1);
1611: glob2loc[i] = dof;
1612: }
1614: ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,size,glob2loc,PETSC_OWN_POINTER,ltog);
1615: #if 0
1616: PetscIntView(size,glob2loc,PETSC_VIEWER_STDOUT_WORLD);
1617: #endif
1618: return(0);
1619: }
1621: #include <petsc/private/matimpl.h>
1623: PetscErrorCode DMCreateMatrix_Network_Nest(DM dm,Mat *J)
1624: {
1626: DM_Network *network = (DM_Network*)dm->data;
1627: PetscMPIInt rank, size;
1628: PetscInt eDof,vDof;
1629: Mat j11,j12,j21,j22,bA[2][2];
1630: MPI_Comm comm;
1631: ISLocalToGlobalMapping eISMap,vISMap;
1634: PetscObjectGetComm((PetscObject)dm,&comm);
1635: MPI_Comm_rank(comm,&rank);
1636: MPI_Comm_size(comm,&size);
1638: PetscSectionGetConstrainedStorageSize(network->edge.GlobalDofSection,&eDof);
1639: PetscSectionGetConstrainedStorageSize(network->vertex.GlobalDofSection,&vDof);
1641: MatCreate(comm, &j11);
1642: MatSetSizes(j11, eDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
1643: MatSetType(j11, MATMPIAIJ);
1645: MatCreate(comm, &j12);
1646: MatSetSizes(j12, eDof, vDof, PETSC_DETERMINE ,PETSC_DETERMINE);
1647: MatSetType(j12, MATMPIAIJ);
1649: MatCreate(comm, &j21);
1650: MatSetSizes(j21, vDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
1651: MatSetType(j21, MATMPIAIJ);
1653: MatCreate(comm, &j22);
1654: MatSetSizes(j22, vDof, vDof, PETSC_DETERMINE, PETSC_DETERMINE);
1655: MatSetType(j22, MATMPIAIJ);
1657: bA[0][0] = j11;
1658: bA[0][1] = j12;
1659: bA[1][0] = j21;
1660: bA[1][1] = j22;
1662: CreateSubGlobalToLocalMapping_private(network->edge.GlobalDofSection,network->edge.DofSection,&eISMap);
1663: CreateSubGlobalToLocalMapping_private(network->vertex.GlobalDofSection,network->vertex.DofSection,&vISMap);
1665: MatSetLocalToGlobalMapping(j11,eISMap,eISMap);
1666: MatSetLocalToGlobalMapping(j12,eISMap,vISMap);
1667: MatSetLocalToGlobalMapping(j21,vISMap,eISMap);
1668: MatSetLocalToGlobalMapping(j22,vISMap,vISMap);
1670: MatSetUp(j11);
1671: MatSetUp(j12);
1672: MatSetUp(j21);
1673: MatSetUp(j22);
1675: MatCreateNest(comm,2,NULL,2,NULL,&bA[0][0],J);
1676: MatSetUp(*J);
1677: MatNestSetVecType(*J,VECNEST);
1678: MatDestroy(&j11);
1679: MatDestroy(&j12);
1680: MatDestroy(&j21);
1681: MatDestroy(&j22);
1683: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1684: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1685: MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1687: /* Free structures */
1688: ISLocalToGlobalMappingDestroy(&eISMap);
1689: ISLocalToGlobalMappingDestroy(&vISMap);
1690: return(0);
1691: }
1693: PetscErrorCode DMCreateMatrix_Network(DM dm,Mat *J)
1694: {
1696: DM_Network *network = (DM_Network*)dm->data;
1697: PetscInt eStart,eEnd,vStart,vEnd,rstart,nrows,*rows,localSize;
1698: PetscInt cstart,ncols,j,e,v;
1699: PetscBool ghost,ghost_vc,ghost2,isNest;
1700: Mat Juser;
1701: PetscSection sectionGlobal;
1702: PetscInt nedges,*vptr=NULL,vc,*rows_v; /* suppress maybe-uninitialized warning */
1703: const PetscInt *edges,*cone;
1704: MPI_Comm comm;
1705: MatType mtype;
1706: Vec vd_nz,vo_nz;
1707: PetscInt *dnnz,*onnz;
1708: PetscScalar *vdnz,*vonz;
1711: mtype = dm->mattype;
1712: PetscStrcmp(mtype,MATNEST,&isNest);
1713: if (isNest) {
1714: DMCreateMatrix_Network_Nest(dm,J);
1715: return(0);
1716: }
1718: if (!network->userEdgeJacobian && !network->userVertexJacobian) {
1719: /* user does not provide Jacobian blocks */
1720: DMCreateMatrix_Plex(network->plex,J);
1721: MatSetDM(*J,dm);
1722: return(0);
1723: }
1725: MatCreate(PetscObjectComm((PetscObject)dm),J);
1726: DMGetGlobalSection(network->plex,§ionGlobal);
1727: PetscSectionGetConstrainedStorageSize(sectionGlobal,&localSize);
1728: MatSetSizes(*J,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);
1730: MatSetType(*J,MATAIJ);
1731: MatSetFromOptions(*J);
1733: /* (1) Set matrix preallocation */
1734: /*------------------------------*/
1735: PetscObjectGetComm((PetscObject)dm,&comm);
1736: VecCreate(comm,&vd_nz);
1737: VecSetSizes(vd_nz,localSize,PETSC_DECIDE);
1738: VecSetFromOptions(vd_nz);
1739: VecSet(vd_nz,0.0);
1740: VecDuplicate(vd_nz,&vo_nz);
1742: /* Set preallocation for edges */
1743: /*-----------------------------*/
1744: DMNetworkGetEdgeRange(dm,&eStart,&eEnd);
1746: PetscMalloc1(localSize,&rows);
1747: for (e=eStart; e<eEnd; e++) {
1748: /* Get row indices */
1749: DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
1750: DMNetworkGetNumVariables(dm,e,&nrows);
1751: if (nrows) {
1752: for (j=0; j<nrows; j++) rows[j] = j + rstart;
1754: /* Set preallocation for conntected vertices */
1755: DMNetworkGetConnectedVertices(dm,e,&cone);
1756: for (v=0; v<2; v++) {
1757: DMNetworkGetNumVariables(dm,cone[v],&ncols);
1759: if (network->Je) {
1760: Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
1761: } else Juser = NULL;
1762: DMNetworkIsGhostVertex(dm,cone[v],&ghost);
1763: MatSetPreallocationblock_private(Juser,nrows,rows,ncols,ghost,vd_nz,vo_nz);
1764: }
1766: /* Set preallocation for edge self */
1767: cstart = rstart;
1768: if (network->Je) {
1769: Juser = network->Je[3*e]; /* Jacobian(e,e) */
1770: } else Juser = NULL;
1771: MatSetPreallocationblock_private(Juser,nrows,rows,nrows,PETSC_FALSE,vd_nz,vo_nz);
1772: }
1773: }
1775: /* Set preallocation for vertices */
1776: /*--------------------------------*/
1777: DMNetworkGetVertexRange(dm,&vStart,&vEnd);
1778: if (vEnd - vStart) vptr = network->Jvptr;
1780: for (v=vStart; v<vEnd; v++) {
1781: /* Get row indices */
1782: DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
1783: DMNetworkGetNumVariables(dm,v,&nrows);
1784: if (!nrows) continue;
1786: DMNetworkIsGhostVertex(dm,v,&ghost);
1787: if (ghost) {
1788: PetscMalloc1(nrows,&rows_v);
1789: } else {
1790: rows_v = rows;
1791: }
1793: for (j=0; j<nrows; j++) rows_v[j] = j + rstart;
1795: /* Get supporting edges and connected vertices */
1796: DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
1798: for (e=0; e<nedges; e++) {
1799: /* Supporting edges */
1800: DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
1801: DMNetworkGetNumVariables(dm,edges[e],&ncols);
1803: if (network->Jv) {
1804: Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */
1805: } else Juser = NULL;
1806: MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost,vd_nz,vo_nz);
1808: /* Connected vertices */
1809: DMNetworkGetConnectedVertices(dm,edges[e],&cone);
1810: vc = (v == cone[0]) ? cone[1]:cone[0];
1811: DMNetworkIsGhostVertex(dm,vc,&ghost_vc);
1813: DMNetworkGetNumVariables(dm,vc,&ncols);
1815: if (network->Jv) {
1816: Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
1817: } else Juser = NULL;
1818: if (ghost_vc||ghost) {
1819: ghost2 = PETSC_TRUE;
1820: } else {
1821: ghost2 = PETSC_FALSE;
1822: }
1823: MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost2,vd_nz,vo_nz);
1824: }
1826: /* Set preallocation for vertex self */
1827: DMNetworkIsGhostVertex(dm,v,&ghost);
1828: if (!ghost) {
1829: DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
1830: if (network->Jv) {
1831: Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
1832: } else Juser = NULL;
1833: MatSetPreallocationblock_private(Juser,nrows,rows_v,nrows,PETSC_FALSE,vd_nz,vo_nz);
1834: }
1835: if (ghost) {
1836: PetscFree(rows_v);
1837: }
1838: }
1840: VecAssemblyBegin(vd_nz);
1841: VecAssemblyBegin(vo_nz);
1843: PetscMalloc2(localSize,&dnnz,localSize,&onnz);
1845: VecAssemblyEnd(vd_nz);
1846: VecAssemblyEnd(vo_nz);
1848: VecGetArray(vd_nz,&vdnz);
1849: VecGetArray(vo_nz,&vonz);
1850: for (j=0; j<localSize; j++) {
1851: dnnz[j] = (PetscInt)PetscRealPart(vdnz[j]);
1852: onnz[j] = (PetscInt)PetscRealPart(vonz[j]);
1853: }
1854: VecRestoreArray(vd_nz,&vdnz);
1855: VecRestoreArray(vo_nz,&vonz);
1856: VecDestroy(&vd_nz);
1857: VecDestroy(&vo_nz);
1859: MatSeqAIJSetPreallocation(*J,0,dnnz);
1860: MatMPIAIJSetPreallocation(*J,0,dnnz,0,onnz);
1861: MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1863: PetscFree2(dnnz,onnz);
1865: /* (2) Set matrix entries for edges */
1866: /*----------------------------------*/
1867: for (e=eStart; e<eEnd; e++) {
1868: /* Get row indices */
1869: DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
1870: DMNetworkGetNumVariables(dm,e,&nrows);
1871: if (nrows) {
1872: for (j=0; j<nrows; j++) rows[j] = j + rstart;
1874: /* Set matrix entries for conntected vertices */
1875: DMNetworkGetConnectedVertices(dm,e,&cone);
1876: for (v=0; v<2; v++) {
1877: DMNetworkGetVariableGlobalOffset(dm,cone[v],&cstart);
1878: DMNetworkGetNumVariables(dm,cone[v],&ncols);
1880: if (network->Je) {
1881: Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
1882: } else Juser = NULL;
1883: MatSetblock_private(Juser,nrows,rows,ncols,cstart,J);
1884: }
1886: /* Set matrix entries for edge self */
1887: cstart = rstart;
1888: if (network->Je) {
1889: Juser = network->Je[3*e]; /* Jacobian(e,e) */
1890: } else Juser = NULL;
1891: MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);
1892: }
1893: }
1895: /* Set matrix entries for vertices */
1896: /*---------------------------------*/
1897: for (v=vStart; v<vEnd; v++) {
1898: /* Get row indices */
1899: DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
1900: DMNetworkGetNumVariables(dm,v,&nrows);
1901: if (!nrows) continue;
1903: DMNetworkIsGhostVertex(dm,v,&ghost);
1904: if (ghost) {
1905: PetscMalloc1(nrows,&rows_v);
1906: } else {
1907: rows_v = rows;
1908: }
1909: for (j=0; j<nrows; j++) rows_v[j] = j + rstart;
1911: /* Get supporting edges and connected vertices */
1912: DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
1914: for (e=0; e<nedges; e++) {
1915: /* Supporting edges */
1916: DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
1917: DMNetworkGetNumVariables(dm,edges[e],&ncols);
1919: if (network->Jv) {
1920: Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */
1921: } else Juser = NULL;
1922: MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);
1924: /* Connected vertices */
1925: DMNetworkGetConnectedVertices(dm,edges[e],&cone);
1926: vc = (v == cone[0]) ? cone[1]:cone[0];
1928: DMNetworkGetVariableGlobalOffset(dm,vc,&cstart);
1929: DMNetworkGetNumVariables(dm,vc,&ncols);
1931: if (network->Jv) {
1932: Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
1933: } else Juser = NULL;
1934: MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);
1935: }
1937: /* Set matrix entries for vertex self */
1938: if (!ghost) {
1939: DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
1940: if (network->Jv) {
1941: Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
1942: } else Juser = NULL;
1943: MatSetblock_private(Juser,nrows,rows_v,nrows,cstart,J);
1944: }
1945: if (ghost) {
1946: PetscFree(rows_v);
1947: }
1948: }
1949: PetscFree(rows);
1951: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1952: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1954: MatSetDM(*J,dm);
1955: return(0);
1956: }
1958: PetscErrorCode DMDestroy_Network(DM dm)
1959: {
1961: DM_Network *network = (DM_Network*)dm->data;
1962: PetscInt j;
1965: if (--network->refct > 0) return(0);
1966: if (network->Je) {
1967: PetscFree(network->Je);
1968: }
1969: if (network->Jv) {
1970: PetscFree(network->Jvptr);
1971: PetscFree(network->Jv);
1972: }
1974: ISLocalToGlobalMappingDestroy(&network->vertex.mapping);
1975: PetscSectionDestroy(&network->vertex.DofSection);
1976: PetscSectionDestroy(&network->vertex.GlobalDofSection);
1977: if (network->vltog) {
1978: PetscFree(network->vltog);
1979: }
1980: if (network->vertex.sf) {
1981: PetscSFDestroy(&network->vertex.sf);
1982: }
1983: /* edge */
1984: ISLocalToGlobalMappingDestroy(&network->edge.mapping);
1985: PetscSectionDestroy(&network->edge.DofSection);
1986: PetscSectionDestroy(&network->edge.GlobalDofSection);
1987: if (network->edge.sf) {
1988: PetscSFDestroy(&network->edge.sf);
1989: }
1990: DMDestroy(&network->plex);
1991: PetscSectionDestroy(&network->DataSection);
1992: PetscSectionDestroy(&network->DofSection);
1994: for(j=0; j<network->nsubnet; j++) {
1995: PetscFree(network->subnet[j].edges);
1996: }
1997: PetscFree(network->subnetvtx);
1999: PetscFree(network->subnet);
2000: PetscFree(network->componentdataarray);
2001: PetscFree2(network->header,network->cvalue);
2002: PetscFree(network);
2003: return(0);
2004: }
2006: PetscErrorCode DMView_Network(DM dm,PetscViewer viewer)
2007: {
2009: DM_Network *network = (DM_Network*)dm->data;
2010: PetscBool iascii;
2011: PetscMPIInt rank;
2012: PetscInt p,nsubnet;
2015: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
2016: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
2019: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2020: if (iascii) {
2021: const PetscInt *cone,*vtx,*edges;
2022: PetscInt vfrom,vto,i,j,nv,ne;
2024: nsubnet = network->nsubnet - network->ncsubnet; /* num of subnetworks */
2025: PetscViewerASCIIPushSynchronized(viewer);
2026: PetscViewerASCIISynchronizedPrintf(viewer, " [%d] nsubnet: %D; nsubnetCouple: %D; nEdges: %D; nVertices: %D\n",rank,nsubnet,network->ncsubnet,network->nEdges,network->nVertices);
2028: for (i=0; i<nsubnet; i++) {
2029: DMNetworkGetSubnetworkInfo(dm,i,&nv,&ne,&vtx,&edges);
2030: if (ne) {
2031: PetscViewerASCIISynchronizedPrintf(viewer, " Subnet %D: nEdges %D, nVertices %D\n",i,ne,nv);
2032: for (j=0; j<ne; j++) {
2033: p = edges[j];
2034: DMNetworkGetConnectedVertices(dm,p,&cone);
2035: DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);
2036: DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);
2037: DMNetworkGetGlobalEdgeIndex(dm,edges[j],&p);
2038: PetscViewerASCIISynchronizedPrintf(viewer, " edge %D: %D----> %D\n",p,vfrom,vto);
2039: }
2040: }
2041: }
2042: /* Coupling subnets */
2043: nsubnet = network->nsubnet;
2044: for (; i<nsubnet; i++) {
2045: DMNetworkGetSubnetworkInfo(dm,i,&nv,&ne,&vtx,&edges);
2046: if (ne) {
2047: PetscViewerASCIISynchronizedPrintf(viewer, " Subnet %D (couple): nEdges %D, nVertices %D\n",i,ne,nv);
2048: for (j=0; j<ne; j++) {
2049: p = edges[j];
2050: DMNetworkGetConnectedVertices(dm,p,&cone);
2051: DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);
2052: DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);
2053: PetscViewerASCIISynchronizedPrintf(viewer, " edge %D: %D----> %D\n",p,vfrom,vto);
2054: }
2055: }
2056: }
2057: PetscViewerFlush(viewer);
2058: PetscViewerASCIIPopSynchronized(viewer);
2059: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMNetwork writing", ((PetscObject)viewer)->type_name);
2060: return(0);
2061: }
2063: PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l)
2064: {
2066: DM_Network *network = (DM_Network*)dm->data;
2069: DMGlobalToLocalBegin(network->plex,g,mode,l);
2070: return(0);
2071: }
2073: PetscErrorCode DMGlobalToLocalEnd_Network(DM dm, Vec g, InsertMode mode, Vec l)
2074: {
2076: DM_Network *network = (DM_Network*)dm->data;
2079: DMGlobalToLocalEnd(network->plex,g,mode,l);
2080: return(0);
2081: }
2083: PetscErrorCode DMLocalToGlobalBegin_Network(DM dm, Vec l, InsertMode mode, Vec g)
2084: {
2086: DM_Network *network = (DM_Network*)dm->data;
2089: DMLocalToGlobalBegin(network->plex,l,mode,g);
2090: return(0);
2091: }
2093: PetscErrorCode DMLocalToGlobalEnd_Network(DM dm, Vec l, InsertMode mode, Vec g)
2094: {
2096: DM_Network *network = (DM_Network*)dm->data;
2099: DMLocalToGlobalEnd(network->plex,l,mode,g);
2100: return(0);
2101: }
2103: /*@
2104: DMNetworkGetVertexLocalToGlobalOrdering - Get vertex globle index
2106: Not collective
2108: Input Parameters
2109: + dm - the dm object
2110: - vloc - local vertex ordering, start from 0
2112: Output Parameters
2113: + vg - global vertex ordering, start from 0
2115: Level: Advanced
2117: .seealso: DMNetworkSetVertexLocalToGlobalOrdering()
2118: @*/
2119: PetscErrorCode DMNetworkGetVertexLocalToGlobalOrdering(DM dm,PetscInt vloc,PetscInt *vg)
2120: {
2121: DM_Network *network = (DM_Network*)dm->data;
2122: PetscInt *vltog = network->vltog;
2125: if (!vltog) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkSetVertexLocalToGlobalOrdering() first");
2126: *vg = vltog[vloc];
2127: return(0);
2128: }
2130: /*@
2131: DMNetworkSetVertexLocalToGlobalOrdering - Create and setup vertex local to globle map
2133: Collective
2135: Input Parameters:
2136: + dm - the dm object
2138: Level: Advanced
2140: .seealso: DMNetworkGetGlobalVertexIndex()
2141: @*/
2142: PetscErrorCode DMNetworkSetVertexLocalToGlobalOrdering(DM dm)
2143: {
2144: PetscErrorCode ierr;
2145: DM_Network *network=(DM_Network*)dm->data;
2146: MPI_Comm comm;
2147: PetscMPIInt rank,size,*displs,*recvcounts,remoterank;
2148: PetscBool ghost;
2149: PetscInt *vltog,nroots,nleaves,i,*vrange,k,N,lidx;
2150: const PetscSFNode *iremote;
2151: PetscSF vsf;
2152: Vec Vleaves,Vleaves_seq;
2153: VecScatter ctx;
2154: PetscScalar *varr,val;
2155: const PetscScalar *varr_read;
2158: PetscObjectGetComm((PetscObject)dm,&comm);
2159: MPI_Comm_size(comm,&size);
2160: MPI_Comm_rank(comm,&rank);
2162: if (size == 1) {
2163: nroots = network->vEnd - network->vStart;
2164: PetscMalloc1(nroots, &vltog);
2165: for (i=0; i<nroots; i++) vltog[i] = i;
2166: network->vltog = vltog;
2167: return(0);
2168: }
2170: if (!network->distributecalled) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkDistribute() first");
2171: if (network->vltog) {
2172: PetscFree(network->vltog);
2173: }
2175: DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);
2176: PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);
2177: vsf = network->vertex.sf;
2179: PetscMalloc3(size+1,&vrange,size+1,&displs,size,&recvcounts);
2180: PetscSFGetGraph(vsf,&nroots,&nleaves,NULL,&iremote);
2182: for (i=0; i<size; i++) { displs[i] = i; recvcounts[i] = 1;}
2184: i = nroots - nleaves; /* local number of vertices, excluding ghosts */
2185: vrange[0] = 0;
2186: MPI_Allgatherv(&i,1,MPIU_INT,vrange+1,recvcounts,displs,MPIU_INT,comm);
2187: for (i=2; i<size+1; i++) {vrange[i] += vrange[i-1];}
2189: PetscMalloc1(nroots, &vltog);
2190: network->vltog = vltog;
2192: /* Set vltog for non-ghost vertices */
2193: k = 0;
2194: for (i=0; i<nroots; i++) {
2195: DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);
2196: if (ghost) continue;
2197: vltog[i] = vrange[rank] + k++;
2198: }
2199: PetscFree3(vrange,displs,recvcounts);
2201: /* Set vltog for ghost vertices */
2202: /* (a) create parallel Vleaves and sequential Vleaves_seq to convert local iremote[*].index to global index */
2203: VecCreate(comm,&Vleaves);
2204: VecSetSizes(Vleaves,2*nleaves,PETSC_DETERMINE);
2205: VecSetFromOptions(Vleaves);
2206: VecGetArray(Vleaves,&varr);
2207: for (i=0; i<nleaves; i++) {
2208: varr[2*i] = (PetscScalar)(iremote[i].rank); /* rank of remote process */
2209: varr[2*i+1] = (PetscScalar)(iremote[i].index); /* local index in remote process */
2210: }
2211: VecRestoreArray(Vleaves,&varr);
2213: /* (b) scatter local info to remote processes via VecScatter() */
2214: VecScatterCreateToAll(Vleaves,&ctx,&Vleaves_seq);
2215: VecScatterBegin(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);
2216: VecScatterEnd(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);
2218: /* (c) convert local indices to global indices in parallel vector Vleaves */
2219: VecGetSize(Vleaves_seq,&N);
2220: VecGetArrayRead(Vleaves_seq,&varr_read);
2221: for (i=0; i<N; i+=2) {
2222: remoterank = (PetscMPIInt)PetscRealPart(varr_read[i]);
2223: if (remoterank == rank) {
2224: k = i+1; /* row number */
2225: lidx = (PetscInt)PetscRealPart(varr_read[i+1]);
2226: val = (PetscScalar)vltog[lidx]; /* global index for non-ghost vertex computed above */
2227: VecSetValues(Vleaves,1,&k,&val,INSERT_VALUES);
2228: }
2229: }
2230: VecRestoreArrayRead(Vleaves_seq,&varr_read);
2231: VecAssemblyBegin(Vleaves);
2232: VecAssemblyEnd(Vleaves);
2234: /* (d) Set vltog for ghost vertices by copying local values of Vleaves */
2235: VecGetArrayRead(Vleaves,&varr_read);
2236: k = 0;
2237: for (i=0; i<nroots; i++) {
2238: DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);
2239: if (!ghost) continue;
2240: vltog[i] = (PetscInt)PetscRealPart(varr_read[2*k+1]); k++;
2241: }
2242: VecRestoreArrayRead(Vleaves,&varr_read);
2244: VecDestroy(&Vleaves);
2245: VecDestroy(&Vleaves_seq);
2246: VecScatterDestroy(&ctx);
2247: return(0);
2248: }