Actual source code: network.c
petsc-3.14.6 2021-03-30
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,dim = 1; /* One dimensional network */
205: PetscInt i,j,ctr,nsubnet,*eowners,np,*edges,*subnetvtx,vStart;
206: PetscInt k,netid,vid, *vidxlTog,*edgelist_couple=NULL;
207: const PetscInt *cone;
208: MPI_Comm comm;
209: PetscMPIInt size,rank;
212: PetscObjectGetComm((PetscObject)dm,&comm);
213: MPI_Comm_rank(comm,&rank);
214: MPI_Comm_size(comm,&size);
216: /* Create the local edgelist for the network by concatenating local input edgelists of the subnetworks */
217: PetscCalloc1(2*network->nEdges,&edges);
218: nsubnet = network->nsubnet - network->ncsubnet;
219: ctr = 0;
220: for (i=0; i < nsubnet; i++) {
221: for (j = 0; j < network->subnet[i].nedge; j++) {
222: edges[2*ctr] = network->subnet[i].vStart + network->subnet[i].edgelist[2*j];
223: edges[2*ctr+1] = network->subnet[i].vStart + network->subnet[i].edgelist[2*j+1];
224: ctr++;
225: }
226: }
228: /* Append local coupling edgelists of the subnetworks */
229: i = nsubnet; /* netid of coupling subnet */
230: nsubnet = network->nsubnet;
231: while (i < nsubnet) {
232: edgelist_couple = network->subnet[i].edgelist;
234: k = 0;
235: for (j = 0; j < network->subnet[i].nedge; j++) {
236: netid = edgelist_couple[k]; vid = edgelist_couple[k+1];
237: edges[2*ctr] = network->subnet[netid].vStart + vid; k += 2;
239: netid = edgelist_couple[k]; vid = edgelist_couple[k+1];
240: edges[2*ctr+1] = network->subnet[netid].vStart + vid; k+=2;
241: ctr++;
242: }
243: i++;
244: }
245: /*
246: if (rank == 0) {
247: PetscPrintf(PETSC_COMM_SELF,"[%d] edgelist:\n",rank);
248: for (i=0; i < network->nEdges; i++) {
249: PetscPrintf(PETSC_COMM_SELF,"[%D %D]",edges[2*i],edges[2*i+1]);
250: printf("\n");
251: }
252: }
253: */
255: /* Create network->plex */
256: DMCreate(comm,&network->plex);
257: DMSetType(network->plex,DMPLEX);
258: DMSetDimension(network->plex,dim);
259: if (size == 1) {
260: DMPlexBuildFromCellList(network->plex,network->nEdges,network->nVertices,numCorners,edges);
261: } else {
262: DMPlexBuildFromCellListParallel(network->plex,network->nEdges,network->nVertices,network->NVertices,numCorners,edges,NULL);
263: }
265: DMPlexGetChart(network->plex,&network->pStart,&network->pEnd);
266: DMPlexGetHeightStratum(network->plex,0,&network->eStart,&network->eEnd);
267: DMPlexGetHeightStratum(network->plex,1,&network->vStart,&network->vEnd);
268: vStart = network->vStart;
270: PetscSectionCreate(comm,&network->DataSection);
271: PetscSectionCreate(comm,&network->DofSection);
272: PetscSectionSetChart(network->DataSection,network->pStart,network->pEnd);
273: PetscSectionSetChart(network->DofSection,network->pStart,network->pEnd);
275: network->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
276: np = network->pEnd - network->pStart;
277: PetscCalloc2(np,&network->header,np,&network->cvalue);
279: /* Create vidxlTog: maps local vertex index to global index */
280: np = network->vEnd - vStart;
281: PetscMalloc2(np,&vidxlTog,size+1,&eowners);
282: ctr = 0;
283: for (i=network->eStart; i<network->eEnd; i++) {
284: DMNetworkGetConnectedVertices(dm,i,&cone);
285: vidxlTog[cone[0] - vStart] = edges[2*ctr];
286: vidxlTog[cone[1] - vStart] = edges[2*ctr+1];
287: ctr++;
288: }
289: PetscFree(edges);
291: /* Create vertices and edges array for the subnetworks */
292: for (j=0; j < network->nsubnet; j++) {
293: PetscCalloc1(network->subnet[j].nedge,&network->subnet[j].edges);
295: /* Temporarily setting nvtx and nedge to 0 so we can use them as counters in the below for loop.
296: These get updated when the vertices and edges are added. */
297: network->subnet[j].nvtx = 0;
298: network->subnet[j].nedge = 0;
299: }
300: PetscCalloc1(np,&network->subnetvtx);
303: /* Get edge ownership */
304: np = network->eEnd - network->eStart;
305: MPI_Allgather(&np,1,MPIU_INT,eowners+1,1,MPIU_INT,comm);
306: eowners[0] = 0;
307: for (i=2; i<=size; i++) eowners[i] += eowners[i-1];
309: i = 0; j = 0;
310: while (i < np) { /* local edges, including coupling edges */
311: network->header[i].index = i + eowners[rank]; /* Global edge index */
313: if (j < network->nsubnet && i < network->subnet[j].eEnd) {
314: network->header[i].subnetid = j; /* Subnetwork id */
315: network->subnet[j].edges[network->subnet[j].nedge++] = i;
317: network->header[i].ndata = 0;
318: PetscSectionAddDof(network->DataSection,i,network->dataheadersize);
319: network->header[i].offset[0] = 0;
320: network->header[i].offsetvarrel[0] = 0;
321: i++;
322: }
323: if (i >= network->subnet[j].eEnd) j++;
324: }
326: /* Count network->subnet[*].nvtx */
327: for (i=vStart; i<network->vEnd; i++) { /* local vertices, including ghosts */
328: k = vidxlTog[i-vStart];
329: for (j=0; j < network->nsubnet; j++) {
330: if (network->subnet[j].vStart <= k && k < network->subnet[j].vEnd) {
331: network->subnet[j].nvtx++;
332: break;
333: }
334: }
335: }
337: /* Set network->subnet[*].vertices on array network->subnetvtx */
338: subnetvtx = network->subnetvtx;
339: for (j=0; j<network->nsubnet; j++) {
340: network->subnet[j].vertices = subnetvtx;
341: subnetvtx += network->subnet[j].nvtx;
342: network->subnet[j].nvtx = 0;
343: }
345: /* Set vertex array for the subnetworks */
346: for (i=vStart; i<network->vEnd; i++) { /* local vertices, including ghosts */
347: network->header[i].index = vidxlTog[i-vStart]; /* Global vertex index */
349: k = vidxlTog[i-vStart];
350: for (j=0; j < network->nsubnet; j++) {
351: if (network->subnet[j].vStart <= k && k < network->subnet[j].vEnd) {
352: network->header[i].subnetid = j;
353: network->subnet[j].vertices[network->subnet[j].nvtx++] = i;
354: break;
355: }
356: }
358: network->header[i].ndata = 0;
359: PetscSectionAddDof(network->DataSection,i,network->dataheadersize);
360: network->header[i].offset[0] = 0;
361: network->header[i].offsetvarrel[0] = 0;
362: }
364: PetscFree2(vidxlTog,eowners);
365: return(0);
366: }
368: /*@C
369: DMNetworkGetSubnetworkInfo - Returns the info for the subnetwork
371: Input Parameters:
372: + dm - the DM object
373: - id - the ID (integer) of the subnetwork
375: Output Parameters:
376: + nv - number of vertices (local)
377: . ne - number of edges (local)
378: . vtx - local vertices for this subnetwork
379: - edge - local edges for this subnetwork
381: Notes:
382: Cannot call this routine before DMNetworkLayoutSetup()
384: Level: intermediate
386: .seealso: DMNetworkLayoutSetUp, DMNetworkCreate
387: @*/
388: PetscErrorCode DMNetworkGetSubnetworkInfo(DM dm,PetscInt id,PetscInt *nv, PetscInt *ne,const PetscInt **vtx, const PetscInt **edge)
389: {
390: DM_Network *network = (DM_Network*)dm->data;
393: if (id >= network->nsubnet) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subnet ID %D exceeds the num of subnets %D",id,network->nsubnet);
394: *nv = network->subnet[id].nvtx;
395: *ne = network->subnet[id].nedge;
396: *vtx = network->subnet[id].vertices;
397: *edge = network->subnet[id].edges;
398: return(0);
399: }
401: /*@C
402: DMNetworkGetSubnetworkCoupleInfo - Returns the info for the coupling subnetwork
404: Input Parameters:
405: + dm - the DM object
406: - id - the ID (integer) of the coupling subnetwork
408: Output Parameters:
409: + ne - number of edges (local)
410: - edge - local edges for this coupling subnetwork
412: Notes:
413: Cannot call this routine before DMNetworkLayoutSetup()
415: Level: intermediate
417: .seealso: DMNetworkGetSubnetworkInfo, DMNetworkLayoutSetUp, DMNetworkCreate
418: @*/
419: PetscErrorCode DMNetworkGetSubnetworkCoupleInfo(DM dm,PetscInt id,PetscInt *ne,const PetscInt **edge)
420: {
421: DM_Network *net = (DM_Network*)dm->data;
422: PetscInt id1;
425: if (net->ncsubnet) {
426: 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);
428: id1 = id + net->nsubnet - net->ncsubnet;
429: *ne = net->subnet[id1].nedge;
430: *edge = net->subnet[id1].edges;
431: } else {
432: *ne = 0;
433: *edge = NULL;
434: }
435: return(0);
436: }
438: /*@C
439: DMNetworkRegisterComponent - Registers the network component
441: Logically collective on dm
443: Input Parameters:
444: + dm - the network object
445: . name - the component name
446: - size - the storage size in bytes for this component data
448: Output Parameters:
449: . key - an integer key that defines the component
451: Notes
452: This routine should be called by all processors before calling DMNetworkLayoutSetup().
454: Level: beginner
456: .seealso: DMNetworkLayoutSetUp, DMNetworkCreate
457: @*/
458: PetscErrorCode DMNetworkRegisterComponent(DM dm,const char *name,size_t size,PetscInt *key)
459: {
460: PetscErrorCode ierr;
461: DM_Network *network = (DM_Network*) dm->data;
462: DMNetworkComponent *component=&network->component[network->ncomponent];
463: PetscBool flg=PETSC_FALSE;
464: PetscInt i;
467: for (i=0; i < network->ncomponent; i++) {
468: PetscStrcmp(component->name,name,&flg);
469: if (flg) {
470: *key = i;
471: return(0);
472: }
473: }
474: if (network->ncomponent == MAX_COMPONENTS) {
475: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Number of components registered exceeds the max %D",MAX_COMPONENTS);
476: }
478: PetscStrcpy(component->name,name);
479: component->size = size/sizeof(DMNetworkComponentGenericDataType);
480: *key = network->ncomponent;
481: network->ncomponent++;
482: return(0);
483: }
485: /*@
486: DMNetworkGetVertexRange - Get the bounds [start, end) for the vertices.
488: Not Collective
490: Input Parameters:
491: . dm - The DMNetwork object
493: Output Parameters:
494: + vStart - The first vertex point
495: - vEnd - One beyond the last vertex point
497: Level: beginner
499: .seealso: DMNetworkGetEdgeRange
500: @*/
501: PetscErrorCode DMNetworkGetVertexRange(DM dm,PetscInt *vStart,PetscInt *vEnd)
502: {
503: DM_Network *network = (DM_Network*)dm->data;
506: if (vStart) *vStart = network->vStart;
507: if (vEnd) *vEnd = network->vEnd;
508: return(0);
509: }
511: /*@
512: DMNetworkGetEdgeRange - Get the bounds [start, end) for the edges.
514: Not Collective
516: Input Parameters:
517: . dm - The DMNetwork object
519: Output Parameters:
520: + eStart - The first edge point
521: - eEnd - One beyond the last edge point
523: Level: beginner
525: .seealso: DMNetworkGetVertexRange
526: @*/
527: PetscErrorCode DMNetworkGetEdgeRange(DM dm,PetscInt *eStart,PetscInt *eEnd)
528: {
529: DM_Network *network = (DM_Network*)dm->data;
532: if (eStart) *eStart = network->eStart;
533: if (eEnd) *eEnd = network->eEnd;
534: return(0);
535: }
537: /*@
538: DMNetworkGetGlobalEdgeIndex - Get the user global numbering for the edge.
540: Not Collective
542: Input Parameters:
543: + dm - DMNetwork object
544: - p - edge point
546: Output Parameters:
547: . index - user global numbering for the edge
549: Level: intermediate
551: .seealso: DMNetworkGetGlobalVertexIndex
552: @*/
553: PetscErrorCode DMNetworkGetGlobalEdgeIndex(DM dm,PetscInt p,PetscInt *index)
554: {
555: PetscErrorCode ierr;
556: DM_Network *network = (DM_Network*)dm->data;
557: PetscInt offsetp;
558: DMNetworkComponentHeader header;
561: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
562: PetscSectionGetOffset(network->DataSection,p,&offsetp);
563: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
564: *index = header->index;
565: return(0);
566: }
568: /*@
569: DMNetworkGetGlobalVertexIndex - Get the user global numbering for the vertex.
571: Not Collective
573: Input Parameters:
574: + dm - DMNetwork object
575: - p - vertex point
577: Output Parameters:
578: . index - user global numbering for the vertex
580: Level: intermediate
582: .seealso: DMNetworkGetGlobalEdgeIndex
583: @*/
584: PetscErrorCode DMNetworkGetGlobalVertexIndex(DM dm,PetscInt p,PetscInt *index)
585: {
586: PetscErrorCode ierr;
587: DM_Network *network = (DM_Network*)dm->data;
588: PetscInt offsetp;
589: DMNetworkComponentHeader header;
592: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
593: PetscSectionGetOffset(network->DataSection,p,&offsetp);
594: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
595: *index = header->index;
596: return(0);
597: }
599: /*
600: DMNetworkGetComponentKeyOffset - Gets the type along with the offset for indexing the
601: component value from the component data array
603: Not Collective
605: Input Parameters:
606: + dm - The DMNetwork object
607: . p - vertex/edge point
608: - compnum - component number
610: Output Parameters:
611: + compkey - the key obtained when registering the component
612: - offset - offset into the component data array associated with the vertex/edge point
614: Notes:
615: Typical usage:
617: DMNetworkGetComponentDataArray(dm, &arr);
618: DMNetworkGetVertex/EdgeRange(dm,&Start,&End);
619: Loop over vertices or edges
620: DMNetworkGetNumComponents(dm,v,&numcomps);
621: Loop over numcomps
622: DMNetworkGetComponentKeyOffset(dm,v,compnum,&key,&offset);
623: compdata = (UserCompDataType)(arr+offset);
625: Level: intermediate
627: .seealso: DMNetworkGetNumComponents, DMNetworkGetComponentDataArray,
628: */
629: PetscErrorCode DMNetworkGetComponentKeyOffset(DM dm,PetscInt p, PetscInt compnum, PetscInt *compkey, PetscInt *offset)
630: {
631: PetscErrorCode ierr;
632: PetscInt offsetp;
633: DMNetworkComponentHeader header;
634: DM_Network *network = (DM_Network*)dm->data;
637: PetscSectionGetOffset(network->DataSection,p,&offsetp);
638: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
639: if (compkey) *compkey = header->key[compnum];
640: if (offset) *offset = offsetp+network->dataheadersize+header->offset[compnum];
641: return(0);
642: }
644: /*@
645: DMNetworkGetComponent - Returns the network component and its key
647: Not Collective
649: Input Parameters:
650: + dm - DMNetwork object
651: . p - edge or vertex point
652: - compnum - component number
654: Output Parameters:
655: + compkey - the key set for this computing during registration
656: - component - the component data
658: Notes:
659: Typical usage:
661: DMNetworkGetVertex/EdgeRange(dm,&Start,&End);
662: Loop over vertices or edges
663: DMNetworkGetNumComponents(dm,v,&numcomps);
664: Loop over numcomps
665: DMNetworkGetComponent(dm,v,compnum,&key,&component);
667: Level: beginner
669: .seealso: DMNetworkGetNumComponents, DMNetworkGetVariableOffset
670: @*/
671: PetscErrorCode DMNetworkGetComponent(DM dm, PetscInt p, PetscInt compnum, PetscInt *key, void **component)
672: {
674: DM_Network *network = (DM_Network*)dm->data;
675: PetscInt offsetd = 0;
678: DMNetworkGetComponentKeyOffset(dm,p,compnum,key,&offsetd);
679: *component = network->componentdataarray+offsetd;
680: return(0);
681: }
683: /*@
684: DMNetworkAddComponent - Adds a network component at the given point (vertex/edge)
686: Not Collective
688: Input Parameters:
689: + dm - The DMNetwork object
690: . p - vertex/edge point
691: . componentkey - component key returned while registering the component
692: - compvalue - pointer to the data structure for the component
694: Level: beginner
696: .seealso: DMNetworkGetVertexRange, DMNetworkGetEdgeRange, DMNetworkRegisterComponent
697: @*/
698: PetscErrorCode DMNetworkAddComponent(DM dm, PetscInt p,PetscInt componentkey,void* compvalue)
699: {
700: DM_Network *network = (DM_Network*)dm->data;
701: DMNetworkComponent *component = &network->component[componentkey];
702: DMNetworkComponentHeader header = &network->header[p];
703: DMNetworkComponentValue cvalue = &network->cvalue[p];
704: PetscErrorCode ierr;
707: 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);
709: header->size[header->ndata] = component->size;
710: PetscSectionAddDof(network->DataSection,p,component->size);
711: header->key[header->ndata] = componentkey;
712: if (header->ndata != 0) header->offset[header->ndata] = header->offset[header->ndata-1] + header->size[header->ndata-1];
713: header->nvar[header->ndata] = 0;
715: cvalue->data[header->ndata] = (void*)compvalue;
716: header->ndata++;
717: return(0);
718: }
720: /*@
721: DMNetworkSetComponentNumVariables - Sets the number of variables for a component
723: Not Collective
725: Input Parameters:
726: + dm - The DMNetwork object
727: . p - vertex/edge point
728: . compnum - component number (First component added = 0, second = 1, ...)
729: - nvar - number of variables for the component
731: Level: beginner
733: .seealso: DMNetworkAddComponent(), DMNetworkGetNumComponents(),DMNetworkRegisterComponent()
734: @*/
735: PetscErrorCode DMNetworkSetComponentNumVariables(DM dm, PetscInt p,PetscInt compnum,PetscInt nvar)
736: {
737: DM_Network *network = (DM_Network*)dm->data;
738: DMNetworkComponentHeader header = &network->header[p];
739: PetscErrorCode ierr;
742: DMNetworkAddNumVariables(dm,p,nvar);
743: header->nvar[compnum] = nvar;
744: if (compnum != 0) header->offsetvarrel[compnum] = header->offsetvarrel[compnum-1] + header->nvar[compnum-1];
745: return(0);
746: }
748: /*@
749: DMNetworkGetNumComponents - Get the number of components at a vertex/edge
751: Not Collective
753: Input Parameters:
754: + dm - The DMNetwork object
755: - p - vertex/edge point
757: Output Parameters:
758: . numcomponents - Number of components at the vertex/edge
760: Level: beginner
762: .seealso: DMNetworkRegisterComponent, DMNetworkAddComponent
763: @*/
764: PetscErrorCode DMNetworkGetNumComponents(DM dm,PetscInt p,PetscInt *numcomponents)
765: {
767: PetscInt offset;
768: DM_Network *network = (DM_Network*)dm->data;
771: PetscSectionGetOffset(network->DataSection,p,&offset);
772: *numcomponents = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata;
773: return(0);
774: }
776: /*@
777: DMNetworkGetVariableOffset - Get the offset for accessing the variable associated with the given vertex/edge from the local vector.
779: Not Collective
781: Input Parameters:
782: + dm - The DMNetwork object
783: - p - the edge/vertex point
785: Output Parameters:
786: . offset - the offset
788: Level: beginner
790: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
791: @*/
792: PetscErrorCode DMNetworkGetVariableOffset(DM dm,PetscInt p,PetscInt *offset)
793: {
795: DM_Network *network = (DM_Network*)dm->data;
798: PetscSectionGetOffset(network->plex->localSection,p,offset);
799: return(0);
800: }
802: /*@
803: DMNetworkGetVariableGlobalOffset - Get the global offset for the variable associated with the given vertex/edge from the global vector.
805: Not Collective
807: Input Parameters:
808: + dm - The DMNetwork object
809: - p - the edge/vertex point
811: Output Parameters:
812: . offsetg - the offset
814: Level: beginner
816: .seealso: DMNetworkGetVariableOffset, DMGetLocalVector
817: @*/
818: PetscErrorCode DMNetworkGetVariableGlobalOffset(DM dm,PetscInt p,PetscInt *offsetg)
819: {
821: DM_Network *network = (DM_Network*)dm->data;
824: PetscSectionGetOffset(network->plex->globalSection,p,offsetg);
825: if (*offsetg < 0) *offsetg = -(*offsetg + 1); /* Convert to actual global offset for ghost vertex */
826: return(0);
827: }
829: /*@
830: DMNetworkGetComponentVariableOffset - Get the offset for accessing the variable associated with a component for the given vertex/edge from the local vector.
832: Not Collective
834: Input Parameters:
835: + dm - The DMNetwork object
836: . p - the edge/vertex point
837: - compnum - component number
839: Output Parameters:
840: . offset - the offset
842: Level: intermediate
844: .seealso: DMNetworkGetVariableGlobalOffset(), DMGetLocalVector(), DMNetworkSetComponentNumVariables()
845: @*/
846: PetscErrorCode DMNetworkGetComponentVariableOffset(DM dm,PetscInt p,PetscInt compnum,PetscInt *offset)
847: {
849: DM_Network *network = (DM_Network*)dm->data;
850: PetscInt offsetp,offsetd;
851: DMNetworkComponentHeader header;
854: DMNetworkGetVariableOffset(dm,p,&offsetp);
855: PetscSectionGetOffset(network->DataSection,p,&offsetd);
856: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetd);
857: *offset = offsetp + header->offsetvarrel[compnum];
858: return(0);
859: }
861: /*@
862: DMNetworkGetComponentVariableGlobalOffset - Get the global offset for accessing the variable associated with a component for the given vertex/edge from the local vector.
864: Not Collective
866: Input Parameters:
867: + dm - The DMNetwork object
868: . p - the edge/vertex point
869: - compnum - component number
871: Output Parameters:
872: . offsetg - the global offset
874: Level: intermediate
876: .seealso: DMNetworkGetVariableGlobalOffset(), DMNetworkGetComponentVariableOffset(), DMGetLocalVector(), DMNetworkSetComponentNumVariables()
877: @*/
878: PetscErrorCode DMNetworkGetComponentVariableGlobalOffset(DM dm,PetscInt p,PetscInt compnum,PetscInt *offsetg)
879: {
881: DM_Network *network = (DM_Network*)dm->data;
882: PetscInt offsetp,offsetd;
883: DMNetworkComponentHeader header;
886: DMNetworkGetVariableGlobalOffset(dm,p,&offsetp);
887: PetscSectionGetOffset(network->DataSection,p,&offsetd);
888: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetd);
889: *offsetg = offsetp + header->offsetvarrel[compnum];
890: return(0);
891: }
893: /*@
894: DMNetworkGetEdgeOffset - Get the offset for accessing the variable associated with the given edge from the local subvector.
896: Not Collective
898: Input Parameters:
899: + dm - The DMNetwork object
900: - p - the edge point
902: Output Parameters:
903: . offset - the offset
905: Level: intermediate
907: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
908: @*/
909: PetscErrorCode DMNetworkGetEdgeOffset(DM dm,PetscInt p,PetscInt *offset)
910: {
912: DM_Network *network = (DM_Network*)dm->data;
916: PetscSectionGetOffset(network->edge.DofSection,p,offset);
917: return(0);
918: }
920: /*@
921: DMNetworkGetVertexOffset - Get the offset for accessing the variable associated with the given vertex from the local subvector.
923: Not Collective
925: Input Parameters:
926: + dm - The DMNetwork object
927: - p - the vertex point
929: Output Parameters:
930: . offset - the offset
932: Level: intermediate
934: .seealso: DMNetworkGetVariableGlobalOffset, DMGetLocalVector
935: @*/
936: PetscErrorCode DMNetworkGetVertexOffset(DM dm,PetscInt p,PetscInt *offset)
937: {
939: DM_Network *network = (DM_Network*)dm->data;
943: p -= network->vStart;
945: PetscSectionGetOffset(network->vertex.DofSection,p,offset);
946: return(0);
947: }
948: /*@
949: DMNetworkAddNumVariables - Add number of variables associated with a given point.
951: Not Collective
953: Input Parameters:
954: + dm - The DMNetworkObject
955: . p - the vertex/edge point
956: - nvar - number of additional variables
958: Level: beginner
960: .seealso: DMNetworkSetNumVariables
961: @*/
962: PetscErrorCode DMNetworkAddNumVariables(DM dm,PetscInt p,PetscInt nvar)
963: {
965: DM_Network *network = (DM_Network*)dm->data;
968: PetscSectionAddDof(network->DofSection,p,nvar);
969: return(0);
970: }
972: /*@
973: DMNetworkGetNumVariables - Gets number of variables for a vertex/edge point.
975: Not Collective
977: Input Parameters:
978: + dm - The DMNetworkObject
979: - p - the vertex/edge point
981: Output Parameters:
982: . nvar - number of variables
984: Level: beginner
986: .seealso: DMNetworkAddNumVariables, DMNetworkSddNumVariables
987: @*/
988: PetscErrorCode DMNetworkGetNumVariables(DM dm,PetscInt p,PetscInt *nvar)
989: {
991: DM_Network *network = (DM_Network*)dm->data;
994: PetscSectionGetDof(network->DofSection,p,nvar);
995: return(0);
996: }
998: /*@
999: DMNetworkSetNumVariables - Sets number of variables for a vertex/edge point.
1001: Not Collective
1003: Input Parameters:
1004: + dm - The DMNetworkObject
1005: . p - the vertex/edge point
1006: - nvar - number of variables
1008: Level: beginner
1010: .seealso: DMNetworkAddNumVariables
1011: @*/
1012: PetscErrorCode DMNetworkSetNumVariables(DM dm,PetscInt p,PetscInt nvar)
1013: {
1015: DM_Network *network = (DM_Network*)dm->data;
1018: PetscSectionSetDof(network->DofSection,p,nvar);
1019: return(0);
1020: }
1022: /* Sets up the array that holds the data for all components and its associated section. This
1023: function is called during DMSetUp() */
1024: PetscErrorCode DMNetworkComponentSetUp(DM dm)
1025: {
1026: PetscErrorCode ierr;
1027: DM_Network *network = (DM_Network*)dm->data;
1028: PetscInt arr_size,p,offset,offsetp,ncomp,i;
1029: DMNetworkComponentHeader header;
1030: DMNetworkComponentValue cvalue;
1031: DMNetworkComponentGenericDataType *componentdataarray;
1034: PetscSectionSetUp(network->DataSection);
1035: PetscSectionGetStorageSize(network->DataSection,&arr_size);
1036: PetscMalloc1(arr_size,&network->componentdataarray);
1037: componentdataarray = network->componentdataarray;
1038: for (p = network->pStart; p < network->pEnd; p++) {
1039: PetscSectionGetOffset(network->DataSection,p,&offsetp);
1040: /* Copy header */
1041: header = &network->header[p];
1042: PetscMemcpy(componentdataarray+offsetp,header,network->dataheadersize*sizeof(DMNetworkComponentGenericDataType));
1043: /* Copy data */
1044: cvalue = &network->cvalue[p];
1045: ncomp = header->ndata;
1046: for (i = 0; i < ncomp; i++) {
1047: offset = offsetp + network->dataheadersize + header->offset[i];
1048: PetscMemcpy(componentdataarray+offset,cvalue->data[i],header->size[i]*sizeof(DMNetworkComponentGenericDataType));
1049: }
1050: }
1051: return(0);
1052: }
1054: /* Sets up the section for dofs. This routine is called during DMSetUp() */
1055: PetscErrorCode DMNetworkVariablesSetUp(DM dm)
1056: {
1058: DM_Network *network = (DM_Network*)dm->data;
1061: PetscSectionSetUp(network->DofSection);
1062: return(0);
1063: }
1065: /*
1066: DMNetworkGetComponentDataArray - Returns the component data array
1068: Not Collective
1070: Input Parameters:
1071: . dm - The DMNetwork Object
1073: Output Parameters:
1074: . componentdataarray - array that holds data for all components
1076: Level: intermediate
1078: .seealso: DMNetworkGetComponentKeyOffset, DMNetworkGetNumComponents
1079: */
1080: PetscErrorCode DMNetworkGetComponentDataArray(DM dm,DMNetworkComponentGenericDataType **componentdataarray)
1081: {
1082: DM_Network *network = (DM_Network*)dm->data;
1085: *componentdataarray = network->componentdataarray;
1086: return(0);
1087: }
1089: /* Get a subsection from a range of points */
1090: PetscErrorCode DMNetworkGetSubSection_private(PetscSection master, PetscInt pstart, PetscInt pend,PetscSection *subsection)
1091: {
1093: PetscInt i, nvar;
1096: PetscSectionCreate(PetscObjectComm((PetscObject)master), subsection);
1097: PetscSectionSetChart(*subsection, 0, pend - pstart);
1098: for (i = pstart; i < pend; i++) {
1099: PetscSectionGetDof(master,i,&nvar);
1100: PetscSectionSetDof(*subsection, i - pstart, nvar);
1101: }
1103: PetscSectionSetUp(*subsection);
1104: return(0);
1105: }
1107: /* Create a submap of points with a GlobalToLocal structure */
1108: PetscErrorCode DMNetworkSetSubMap_private(PetscInt pstart, PetscInt pend, ISLocalToGlobalMapping *map)
1109: {
1111: PetscInt i, *subpoints;
1114: /* Create index sets to map from "points" to "subpoints" */
1115: PetscMalloc1(pend - pstart, &subpoints);
1116: for (i = pstart; i < pend; i++) {
1117: subpoints[i - pstart] = i;
1118: }
1119: ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,pend-pstart,subpoints,PETSC_COPY_VALUES,map);
1120: PetscFree(subpoints);
1121: return(0);
1122: }
1124: /*@
1125: DMNetworkAssembleGraphStructures - Assembles vertex and edge data structures. Must be called after DMNetworkDistribute.
1127: Collective
1129: Input Parameters:
1130: . dm - The DMNetworkObject
1132: Note: the routine will create alternative orderings for the vertices and edges. Assume global network points are:
1134: points = [0 1 2 3 4 5 6]
1136: 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]).
1138: With this new ordering a local PetscSection, global PetscSection and PetscSF will be created specific to the subset.
1140: Level: intermediate
1142: @*/
1143: PetscErrorCode DMNetworkAssembleGraphStructures(DM dm)
1144: {
1146: MPI_Comm comm;
1147: PetscMPIInt rank, size;
1148: DM_Network *network = (DM_Network*)dm->data;
1151: PetscObjectGetComm((PetscObject)dm,&comm);
1152: MPI_Comm_rank(comm, &rank);
1153: MPI_Comm_size(comm, &size);
1155: /* Create maps for vertices and edges */
1156: DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);
1157: DMNetworkSetSubMap_private(network->eStart,network->eEnd,&network->edge.mapping);
1159: /* Create local sub-sections */
1160: DMNetworkGetSubSection_private(network->DofSection,network->vStart,network->vEnd,&network->vertex.DofSection);
1161: DMNetworkGetSubSection_private(network->DofSection,network->eStart,network->eEnd,&network->edge.DofSection);
1163: if (size > 1) {
1164: PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);
1166: PetscSectionCreateGlobalSection(network->vertex.DofSection, network->vertex.sf, PETSC_FALSE, PETSC_FALSE, &network->vertex.GlobalDofSection);
1167: PetscSFGetSubSF(network->plex->sf, network->edge.mapping, &network->edge.sf);
1168: PetscSectionCreateGlobalSection(network->edge.DofSection, network->edge.sf, PETSC_FALSE, PETSC_FALSE, &network->edge.GlobalDofSection);
1169: } else {
1170: /* create structures for vertex */
1171: PetscSectionClone(network->vertex.DofSection,&network->vertex.GlobalDofSection);
1172: /* create structures for edge */
1173: PetscSectionClone(network->edge.DofSection,&network->edge.GlobalDofSection);
1174: }
1176: /* Add viewers */
1177: PetscObjectSetName((PetscObject)network->edge.GlobalDofSection,"Global edge dof section");
1178: PetscObjectSetName((PetscObject)network->vertex.GlobalDofSection,"Global vertex dof section");
1179: PetscSectionViewFromOptions(network->edge.GlobalDofSection, NULL, "-edge_global_section_view");
1180: PetscSectionViewFromOptions(network->vertex.GlobalDofSection, NULL, "-vertex_global_section_view");
1181: return(0);
1182: }
1184: /*@
1185: DMNetworkDistribute - Distributes the network and moves associated component data.
1187: Collective
1189: Input Parameter:
1190: + DM - the DMNetwork object
1191: - overlap - The overlap of partitions, 0 is the default
1193: Notes:
1194: Distributes the network with <overlap>-overlapping partitioning of the edges.
1196: Level: intermediate
1198: .seealso: DMNetworkCreate
1199: @*/
1200: PetscErrorCode DMNetworkDistribute(DM *dm,PetscInt overlap)
1201: {
1202: MPI_Comm comm;
1204: PetscMPIInt size;
1205: DM_Network *oldDMnetwork = (DM_Network*)((*dm)->data);
1206: DM_Network *newDMnetwork;
1207: PetscSF pointsf=NULL;
1208: DM newDM;
1209: PetscInt j,e,v,offset,*subnetvtx;
1210: PetscPartitioner part;
1211: DMNetworkComponentHeader header;
1214: PetscObjectGetComm((PetscObject)*dm,&comm);
1215: MPI_Comm_size(comm, &size);
1216: if (size == 1) return(0);
1218: DMNetworkCreate(PetscObjectComm((PetscObject)*dm),&newDM);
1219: newDMnetwork = (DM_Network*)newDM->data;
1220: newDMnetwork->dataheadersize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(DMNetworkComponentGenericDataType);
1222: /* Enable runtime options for petscpartitioner */
1223: DMPlexGetPartitioner(oldDMnetwork->plex,&part);
1224: PetscPartitionerSetFromOptions(part);
1226: /* Distribute plex dm and dof section */
1227: DMPlexDistribute(oldDMnetwork->plex,overlap,&pointsf,&newDMnetwork->plex);
1229: /* Distribute dof section */
1230: PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DofSection);
1231: PetscSFDistributeSection(pointsf,oldDMnetwork->DofSection,NULL,newDMnetwork->DofSection);
1232: PetscSectionCreate(PetscObjectComm((PetscObject)*dm),&newDMnetwork->DataSection);
1234: /* Distribute data and associated section */
1235: DMPlexDistributeData(newDMnetwork->plex,pointsf,oldDMnetwork->DataSection,MPIU_INT,(void*)oldDMnetwork->componentdataarray,newDMnetwork->DataSection,(void**)&newDMnetwork->componentdataarray);
1237: PetscSectionGetChart(newDMnetwork->DataSection,&newDMnetwork->pStart,&newDMnetwork->pEnd);
1238: DMPlexGetHeightStratum(newDMnetwork->plex,0, &newDMnetwork->eStart,&newDMnetwork->eEnd);
1239: DMPlexGetHeightStratum(newDMnetwork->plex,1,&newDMnetwork->vStart,&newDMnetwork->vEnd);
1240: newDMnetwork->nEdges = newDMnetwork->eEnd - newDMnetwork->eStart;
1241: newDMnetwork->nVertices = newDMnetwork->vEnd - newDMnetwork->vStart;
1242: newDMnetwork->NVertices = oldDMnetwork->NVertices;
1243: newDMnetwork->NEdges = oldDMnetwork->NEdges;
1245: /* Set Dof section as the section for dm */
1246: DMSetLocalSection(newDMnetwork->plex,newDMnetwork->DofSection);
1247: DMGetGlobalSection(newDMnetwork->plex,&newDMnetwork->GlobalDofSection);
1249: /* Set up subnetwork info in the newDM */
1250: newDMnetwork->nsubnet = oldDMnetwork->nsubnet;
1251: newDMnetwork->ncsubnet = oldDMnetwork->ncsubnet;
1252: PetscCalloc1(newDMnetwork->nsubnet,&newDMnetwork->subnet);
1253: /* Copy over the global number of vertices and edges in each subnetwork. Note that these are already
1254: calculated in DMNetworkLayoutSetUp()
1255: */
1256: for (j=0; j < newDMnetwork->nsubnet; j++) {
1257: newDMnetwork->subnet[j].Nvtx = oldDMnetwork->subnet[j].Nvtx;
1258: newDMnetwork->subnet[j].Nedge = oldDMnetwork->subnet[j].Nedge;
1259: }
1261: for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++) {
1262: PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);
1263: header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1264: newDMnetwork->subnet[header->subnetid].nedge++;
1265: }
1267: for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++) {
1268: PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);
1269: header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1270: newDMnetwork->subnet[header->subnetid].nvtx++;
1271: }
1273: /* Now create the vertices and edge arrays for the subnetworks */
1274: PetscCalloc1(newDMnetwork->vEnd-newDMnetwork->vStart,&newDMnetwork->subnetvtx);
1275: subnetvtx = newDMnetwork->subnetvtx;
1277: for (j=0; j<newDMnetwork->nsubnet; j++) {
1278: PetscCalloc1(newDMnetwork->subnet[j].nedge,&newDMnetwork->subnet[j].edges);
1279: newDMnetwork->subnet[j].vertices = subnetvtx;
1280: subnetvtx += newDMnetwork->subnet[j].nvtx;
1282: /* Temporarily setting nvtx and nedge to 0 so we can use them as counters in the below for loop.
1283: These get updated when the vertices and edges are added. */
1284: newDMnetwork->subnet[j].nvtx = newDMnetwork->subnet[j].nedge = 0;
1285: }
1287: /* Set the vertices and edges in each subnetwork */
1288: for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++) {
1289: PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);
1290: header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1291: newDMnetwork->subnet[header->subnetid].edges[newDMnetwork->subnet[header->subnetid].nedge++] = e;
1292: }
1294: for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++) {
1295: PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);
1296: header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1297: newDMnetwork->subnet[header->subnetid].vertices[newDMnetwork->subnet[header->subnetid].nvtx++] = v;
1298: }
1300: newDM->setupcalled = (*dm)->setupcalled;
1301: newDMnetwork->distributecalled = PETSC_TRUE;
1303: /* Destroy point SF */
1304: PetscSFDestroy(&pointsf);
1306: DMDestroy(dm);
1307: *dm = newDM;
1308: return(0);
1309: }
1311: /*@C
1312: PetscSFGetSubSF - Returns an SF for a specific subset of points. Leaves are re-numbered to reflect the new ordering.
1314: Input Parameters:
1315: + masterSF - the original SF structure
1316: - map - a ISLocalToGlobal mapping that contains the subset of points
1318: Output Parameters:
1319: . subSF - a subset of the masterSF for the desired subset.
1321: Level: intermediate
1322: @*/
1323: PetscErrorCode PetscSFGetSubSF(PetscSF mastersf, ISLocalToGlobalMapping map, PetscSF *subSF) {
1325: PetscErrorCode ierr;
1326: PetscInt nroots, nleaves, *ilocal_sub;
1327: PetscInt i, *ilocal_map, nroots_sub, nleaves_sub = 0;
1328: PetscInt *local_points, *remote_points;
1329: PetscSFNode *iremote_sub;
1330: const PetscInt *ilocal;
1331: const PetscSFNode *iremote;
1334: PetscSFGetGraph(mastersf,&nroots,&nleaves,&ilocal,&iremote);
1336: /* Look for leaves that pertain to the subset of points. Get the local ordering */
1337: PetscMalloc1(nleaves,&ilocal_map);
1338: ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nleaves,ilocal,NULL,ilocal_map);
1339: for (i = 0; i < nleaves; i++) {
1340: if (ilocal_map[i] != -1) nleaves_sub += 1;
1341: }
1342: /* Re-number ilocal with subset numbering. Need information from roots */
1343: PetscMalloc2(nroots,&local_points,nroots,&remote_points);
1344: for (i = 0; i < nroots; i++) local_points[i] = i;
1345: ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nroots,local_points,NULL,local_points);
1346: PetscSFBcastBegin(mastersf, MPIU_INT, local_points, remote_points);
1347: PetscSFBcastEnd(mastersf, MPIU_INT, local_points, remote_points);
1348: /* Fill up graph using local (that is, local to the subset) numbering. */
1349: PetscMalloc1(nleaves_sub,&ilocal_sub);
1350: PetscMalloc1(nleaves_sub,&iremote_sub);
1351: nleaves_sub = 0;
1352: for (i = 0; i < nleaves; i++) {
1353: if (ilocal_map[i] != -1) {
1354: ilocal_sub[nleaves_sub] = ilocal_map[i];
1355: iremote_sub[nleaves_sub].rank = iremote[i].rank;
1356: iremote_sub[nleaves_sub].index = remote_points[ilocal[i]];
1357: nleaves_sub += 1;
1358: }
1359: }
1360: PetscFree2(local_points,remote_points);
1361: ISLocalToGlobalMappingGetSize(map,&nroots_sub);
1363: /* Create new subSF */
1364: PetscSFCreate(PETSC_COMM_WORLD,subSF);
1365: PetscSFSetFromOptions(*subSF);
1366: PetscSFSetGraph(*subSF,nroots_sub,nleaves_sub,ilocal_sub,PETSC_OWN_POINTER,iremote_sub,PETSC_COPY_VALUES);
1367: PetscFree(ilocal_map);
1368: PetscFree(iremote_sub);
1369: return(0);
1370: }
1372: /*@C
1373: DMNetworkGetSupportingEdges - Return the supporting edges for this vertex point
1375: Not Collective
1377: Input Parameters:
1378: + dm - The DMNetwork object
1379: - p - the vertex point
1381: Output Parameters:
1382: + nedges - number of edges connected to this vertex point
1383: - edges - List of edge points
1385: Level: beginner
1387: Fortran Notes:
1388: Since it returns an array, this routine is only available in Fortran 90, and you must
1389: include petsc.h90 in your code.
1391: .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices
1392: @*/
1393: PetscErrorCode DMNetworkGetSupportingEdges(DM dm,PetscInt vertex,PetscInt *nedges,const PetscInt *edges[])
1394: {
1396: DM_Network *network = (DM_Network*)dm->data;
1399: DMPlexGetSupportSize(network->plex,vertex,nedges);
1400: DMPlexGetSupport(network->plex,vertex,edges);
1401: return(0);
1402: }
1404: /*@C
1405: DMNetworkGetConnectedVertices - Return the connected vertices for this edge point
1407: Not Collective
1409: Input Parameters:
1410: + dm - The DMNetwork object
1411: - p - the edge point
1413: Output Parameters:
1414: . vertices - vertices connected to this edge
1416: Level: beginner
1418: Fortran Notes:
1419: Since it returns an array, this routine is only available in Fortran 90, and you must
1420: include petsc.h90 in your code.
1422: .seealso: DMNetworkCreate, DMNetworkGetSupportingEdges
1423: @*/
1424: PetscErrorCode DMNetworkGetConnectedVertices(DM dm,PetscInt edge,const PetscInt *vertices[])
1425: {
1427: DM_Network *network = (DM_Network*)dm->data;
1430: DMPlexGetCone(network->plex,edge,vertices);
1431: return(0);
1432: }
1434: /*@
1435: DMNetworkIsGhostVertex - Returns TRUE if the vertex is a ghost vertex
1437: Not Collective
1439: Input Parameters:
1440: + dm - The DMNetwork object
1441: - p - the vertex point
1443: Output Parameter:
1444: . isghost - TRUE if the vertex is a ghost point
1446: Level: beginner
1448: .seealso: DMNetworkCreate, DMNetworkGetConnectedVertices, DMNetworkGetVertexRange
1449: @*/
1450: PetscErrorCode DMNetworkIsGhostVertex(DM dm,PetscInt p,PetscBool *isghost)
1451: {
1453: DM_Network *network = (DM_Network*)dm->data;
1454: PetscInt offsetg;
1455: PetscSection sectiong;
1458: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
1459: *isghost = PETSC_FALSE;
1460: DMGetGlobalSection(network->plex,§iong);
1461: PetscSectionGetOffset(sectiong,p,&offsetg);
1462: if (offsetg < 0) *isghost = PETSC_TRUE;
1463: return(0);
1464: }
1466: PetscErrorCode DMSetUp_Network(DM dm)
1467: {
1469: DM_Network *network=(DM_Network*)dm->data;
1472: DMNetworkComponentSetUp(dm);
1473: DMNetworkVariablesSetUp(dm);
1475: DMSetLocalSection(network->plex,network->DofSection);
1476: DMGetGlobalSection(network->plex,&network->GlobalDofSection);
1478: dm->setupcalled = PETSC_TRUE;
1479: DMViewFromOptions(dm,NULL,"-dm_view");
1480: return(0);
1481: }
1483: /*@
1484: DMNetworkHasJacobian - Sets global flag for using user's sub Jacobian matrices
1485: -- replaced by DMNetworkSetOption(network,userjacobian,PETSC_TURE)?
1487: Collective
1489: Input Parameters:
1490: + dm - The DMNetwork object
1491: . eflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for edges
1492: - vflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for vertices
1494: Level: intermediate
1496: @*/
1497: PetscErrorCode DMNetworkHasJacobian(DM dm,PetscBool eflg,PetscBool vflg)
1498: {
1499: DM_Network *network=(DM_Network*)dm->data;
1501: PetscInt nVertices = network->nVertices;
1504: network->userEdgeJacobian = eflg;
1505: network->userVertexJacobian = vflg;
1507: if (eflg && !network->Je) {
1508: PetscCalloc1(3*network->nEdges,&network->Je);
1509: }
1511: if (vflg && !network->Jv && nVertices) {
1512: PetscInt i,*vptr,nedges,vStart=network->vStart;
1513: PetscInt nedges_total;
1514: const PetscInt *edges;
1516: /* count nvertex_total */
1517: nedges_total = 0;
1518: PetscMalloc1(nVertices+1,&vptr);
1520: vptr[0] = 0;
1521: for (i=0; i<nVertices; i++) {
1522: DMNetworkGetSupportingEdges(dm,i+vStart,&nedges,&edges);
1523: nedges_total += nedges;
1524: vptr[i+1] = vptr[i] + 2*nedges + 1;
1525: }
1527: PetscCalloc1(2*nedges_total+nVertices,&network->Jv);
1528: network->Jvptr = vptr;
1529: }
1530: return(0);
1531: }
1533: /*@
1534: DMNetworkEdgeSetMatrix - Sets user-provided Jacobian matrices for this edge to the network
1536: Not Collective
1538: Input Parameters:
1539: + dm - The DMNetwork object
1540: . p - the edge point
1541: - J - array (size = 3) of Jacobian submatrices for this edge point:
1542: J[0]: this edge
1543: J[1] and J[2]: connected vertices, obtained by calling DMNetworkGetConnectedVertices()
1545: Level: advanced
1547: .seealso: DMNetworkVertexSetMatrix
1548: @*/
1549: PetscErrorCode DMNetworkEdgeSetMatrix(DM dm,PetscInt p,Mat J[])
1550: {
1551: DM_Network *network=(DM_Network*)dm->data;
1554: if (!network->Je) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkEdgeSetMatrix");
1556: if (J) {
1557: network->Je[3*p] = J[0];
1558: network->Je[3*p+1] = J[1];
1559: network->Je[3*p+2] = J[2];
1560: }
1561: return(0);
1562: }
1564: /*@
1565: DMNetworkVertexSetMatrix - Sets user-provided Jacobian matrix for this vertex to the network
1567: Not Collective
1569: Input Parameters:
1570: + dm - The DMNetwork object
1571: . p - the vertex point
1572: - J - array of Jacobian (size = 2*(num of supporting edges) + 1) submatrices for this vertex point:
1573: J[0]: this vertex
1574: J[1+2*i]: i-th supporting edge
1575: J[1+2*i+1]: i-th connected vertex
1577: Level: advanced
1579: .seealso: DMNetworkEdgeSetMatrix
1580: @*/
1581: PetscErrorCode DMNetworkVertexSetMatrix(DM dm,PetscInt p,Mat J[])
1582: {
1584: DM_Network *network=(DM_Network*)dm->data;
1585: PetscInt i,*vptr,nedges,vStart=network->vStart;
1586: const PetscInt *edges;
1589: if (!network->Jv) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkVertexSetMatrix");
1591: if (J) {
1592: vptr = network->Jvptr;
1593: network->Jv[vptr[p-vStart]] = J[0]; /* Set Jacobian for this vertex */
1595: /* Set Jacobian for each supporting edge and connected vertex */
1596: DMNetworkGetSupportingEdges(dm,p,&nedges,&edges);
1597: for (i=1; i<=2*nedges; i++) network->Jv[vptr[p-vStart]+i] = J[i];
1598: }
1599: return(0);
1600: }
1602: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1603: {
1605: PetscInt j;
1606: PetscScalar val=(PetscScalar)ncols;
1609: if (!ghost) {
1610: for (j=0; j<nrows; j++) {
1611: VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
1612: }
1613: } else {
1614: for (j=0; j<nrows; j++) {
1615: VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
1616: }
1617: }
1618: return(0);
1619: }
1621: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1622: {
1624: PetscInt j,ncols_u;
1625: PetscScalar val;
1628: if (!ghost) {
1629: for (j=0; j<nrows; j++) {
1630: MatGetRow(Ju,j,&ncols_u,NULL,NULL);
1631: val = (PetscScalar)ncols_u;
1632: VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
1633: MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
1634: }
1635: } else {
1636: for (j=0; j<nrows; j++) {
1637: MatGetRow(Ju,j,&ncols_u,NULL,NULL);
1638: val = (PetscScalar)ncols_u;
1639: VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
1640: MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
1641: }
1642: }
1643: return(0);
1644: }
1646: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
1647: {
1651: if (Ju) {
1652: MatSetPreallocationUserblock_private(Ju,nrows,rows,ncols,ghost,vdnz,vonz);
1653: } else {
1654: MatSetPreallocationDenseblock_private(nrows,rows,ncols,ghost,vdnz,vonz);
1655: }
1656: return(0);
1657: }
1659: PETSC_STATIC_INLINE PetscErrorCode MatSetDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1660: {
1662: PetscInt j,*cols;
1663: PetscScalar *zeros;
1666: PetscCalloc2(ncols,&cols,nrows*ncols,&zeros);
1667: for (j=0; j<ncols; j++) cols[j] = j+ cstart;
1668: MatSetValues(*J,nrows,rows,ncols,cols,zeros,INSERT_VALUES);
1669: PetscFree2(cols,zeros);
1670: return(0);
1671: }
1673: PETSC_STATIC_INLINE PetscErrorCode MatSetUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1674: {
1676: PetscInt j,M,N,row,col,ncols_u;
1677: const PetscInt *cols;
1678: PetscScalar zero=0.0;
1681: MatGetSize(Ju,&M,&N);
1682: if (nrows != M || ncols != N) SETERRQ4(PetscObjectComm((PetscObject)Ju),PETSC_ERR_USER,"%D by %D must equal %D by %D",nrows,ncols,M,N);
1684: for (row=0; row<nrows; row++) {
1685: MatGetRow(Ju,row,&ncols_u,&cols,NULL);
1686: for (j=0; j<ncols_u; j++) {
1687: col = cols[j] + cstart;
1688: MatSetValues(*J,1,&rows[row],1,&col,&zero,INSERT_VALUES);
1689: }
1690: MatRestoreRow(Ju,row,&ncols_u,&cols,NULL);
1691: }
1692: return(0);
1693: }
1695: PETSC_STATIC_INLINE PetscErrorCode MatSetblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
1696: {
1700: if (Ju) {
1701: MatSetUserblock_private(Ju,nrows,rows,ncols,cstart,J);
1702: } else {
1703: MatSetDenseblock_private(nrows,rows,ncols,cstart,J);
1704: }
1705: return(0);
1706: }
1708: /* 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.
1709: */
1710: PetscErrorCode CreateSubGlobalToLocalMapping_private(PetscSection globalsec, PetscSection localsec, ISLocalToGlobalMapping *ltog)
1711: {
1713: PetscInt i,size,dof;
1714: PetscInt *glob2loc;
1717: PetscSectionGetStorageSize(localsec,&size);
1718: PetscMalloc1(size,&glob2loc);
1720: for (i = 0; i < size; i++) {
1721: PetscSectionGetOffset(globalsec,i,&dof);
1722: dof = (dof >= 0) ? dof : -(dof + 1);
1723: glob2loc[i] = dof;
1724: }
1726: ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,size,glob2loc,PETSC_OWN_POINTER,ltog);
1727: #if 0
1728: PetscIntView(size,glob2loc,PETSC_VIEWER_STDOUT_WORLD);
1729: #endif
1730: return(0);
1731: }
1733: #include <petsc/private/matimpl.h>
1735: PetscErrorCode DMCreateMatrix_Network_Nest(DM dm,Mat *J)
1736: {
1738: DM_Network *network = (DM_Network*)dm->data;
1739: PetscMPIInt rank, size;
1740: PetscInt eDof,vDof;
1741: Mat j11,j12,j21,j22,bA[2][2];
1742: MPI_Comm comm;
1743: ISLocalToGlobalMapping eISMap,vISMap;
1746: PetscObjectGetComm((PetscObject)dm,&comm);
1747: MPI_Comm_rank(comm,&rank);
1748: MPI_Comm_size(comm,&size);
1750: PetscSectionGetConstrainedStorageSize(network->edge.GlobalDofSection,&eDof);
1751: PetscSectionGetConstrainedStorageSize(network->vertex.GlobalDofSection,&vDof);
1753: MatCreate(comm, &j11);
1754: MatSetSizes(j11, eDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
1755: MatSetType(j11, MATMPIAIJ);
1757: MatCreate(comm, &j12);
1758: MatSetSizes(j12, eDof, vDof, PETSC_DETERMINE ,PETSC_DETERMINE);
1759: MatSetType(j12, MATMPIAIJ);
1761: MatCreate(comm, &j21);
1762: MatSetSizes(j21, vDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
1763: MatSetType(j21, MATMPIAIJ);
1765: MatCreate(comm, &j22);
1766: MatSetSizes(j22, vDof, vDof, PETSC_DETERMINE, PETSC_DETERMINE);
1767: MatSetType(j22, MATMPIAIJ);
1769: bA[0][0] = j11;
1770: bA[0][1] = j12;
1771: bA[1][0] = j21;
1772: bA[1][1] = j22;
1774: CreateSubGlobalToLocalMapping_private(network->edge.GlobalDofSection,network->edge.DofSection,&eISMap);
1775: CreateSubGlobalToLocalMapping_private(network->vertex.GlobalDofSection,network->vertex.DofSection,&vISMap);
1777: MatSetLocalToGlobalMapping(j11,eISMap,eISMap);
1778: MatSetLocalToGlobalMapping(j12,eISMap,vISMap);
1779: MatSetLocalToGlobalMapping(j21,vISMap,eISMap);
1780: MatSetLocalToGlobalMapping(j22,vISMap,vISMap);
1782: MatSetUp(j11);
1783: MatSetUp(j12);
1784: MatSetUp(j21);
1785: MatSetUp(j22);
1787: MatCreateNest(comm,2,NULL,2,NULL,&bA[0][0],J);
1788: MatSetUp(*J);
1789: MatNestSetVecType(*J,VECNEST);
1790: MatDestroy(&j11);
1791: MatDestroy(&j12);
1792: MatDestroy(&j21);
1793: MatDestroy(&j22);
1795: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1796: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1797: MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1799: /* Free structures */
1800: ISLocalToGlobalMappingDestroy(&eISMap);
1801: ISLocalToGlobalMappingDestroy(&vISMap);
1802: return(0);
1803: }
1805: PetscErrorCode DMCreateMatrix_Network(DM dm,Mat *J)
1806: {
1808: DM_Network *network = (DM_Network*)dm->data;
1809: PetscInt eStart,eEnd,vStart,vEnd,rstart,nrows,*rows,localSize;
1810: PetscInt cstart,ncols,j,e,v;
1811: PetscBool ghost,ghost_vc,ghost2,isNest;
1812: Mat Juser;
1813: PetscSection sectionGlobal;
1814: PetscInt nedges,*vptr=NULL,vc,*rows_v; /* suppress maybe-uninitialized warning */
1815: const PetscInt *edges,*cone;
1816: MPI_Comm comm;
1817: MatType mtype;
1818: Vec vd_nz,vo_nz;
1819: PetscInt *dnnz,*onnz;
1820: PetscScalar *vdnz,*vonz;
1823: mtype = dm->mattype;
1824: PetscStrcmp(mtype,MATNEST,&isNest);
1825: if (isNest) {
1826: DMCreateMatrix_Network_Nest(dm,J);
1827: MatSetDM(*J,dm);
1828: return(0);
1829: }
1831: if (!network->userEdgeJacobian && !network->userVertexJacobian) {
1832: /* user does not provide Jacobian blocks */
1833: DMCreateMatrix_Plex(network->plex,J);
1834: MatSetDM(*J,dm);
1835: return(0);
1836: }
1838: MatCreate(PetscObjectComm((PetscObject)dm),J);
1839: DMGetGlobalSection(network->plex,§ionGlobal);
1840: PetscSectionGetConstrainedStorageSize(sectionGlobal,&localSize);
1841: MatSetSizes(*J,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);
1843: MatSetType(*J,MATAIJ);
1844: MatSetFromOptions(*J);
1846: /* (1) Set matrix preallocation */
1847: /*------------------------------*/
1848: PetscObjectGetComm((PetscObject)dm,&comm);
1849: VecCreate(comm,&vd_nz);
1850: VecSetSizes(vd_nz,localSize,PETSC_DECIDE);
1851: VecSetFromOptions(vd_nz);
1852: VecSet(vd_nz,0.0);
1853: VecDuplicate(vd_nz,&vo_nz);
1855: /* Set preallocation for edges */
1856: /*-----------------------------*/
1857: DMNetworkGetEdgeRange(dm,&eStart,&eEnd);
1859: PetscMalloc1(localSize,&rows);
1860: for (e=eStart; e<eEnd; e++) {
1861: /* Get row indices */
1862: DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
1863: DMNetworkGetNumVariables(dm,e,&nrows);
1864: if (nrows) {
1865: for (j=0; j<nrows; j++) rows[j] = j + rstart;
1867: /* Set preallocation for conntected vertices */
1868: DMNetworkGetConnectedVertices(dm,e,&cone);
1869: for (v=0; v<2; v++) {
1870: DMNetworkGetNumVariables(dm,cone[v],&ncols);
1872: if (network->Je) {
1873: Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
1874: } else Juser = NULL;
1875: DMNetworkIsGhostVertex(dm,cone[v],&ghost);
1876: MatSetPreallocationblock_private(Juser,nrows,rows,ncols,ghost,vd_nz,vo_nz);
1877: }
1879: /* Set preallocation for edge self */
1880: cstart = rstart;
1881: if (network->Je) {
1882: Juser = network->Je[3*e]; /* Jacobian(e,e) */
1883: } else Juser = NULL;
1884: MatSetPreallocationblock_private(Juser,nrows,rows,nrows,PETSC_FALSE,vd_nz,vo_nz);
1885: }
1886: }
1888: /* Set preallocation for vertices */
1889: /*--------------------------------*/
1890: DMNetworkGetVertexRange(dm,&vStart,&vEnd);
1891: if (vEnd - vStart) vptr = network->Jvptr;
1893: for (v=vStart; v<vEnd; v++) {
1894: /* Get row indices */
1895: DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
1896: DMNetworkGetNumVariables(dm,v,&nrows);
1897: if (!nrows) continue;
1899: DMNetworkIsGhostVertex(dm,v,&ghost);
1900: if (ghost) {
1901: PetscMalloc1(nrows,&rows_v);
1902: } else {
1903: rows_v = rows;
1904: }
1906: for (j=0; j<nrows; j++) rows_v[j] = j + rstart;
1908: /* Get supporting edges and connected vertices */
1909: DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
1911: for (e=0; e<nedges; e++) {
1912: /* Supporting edges */
1913: DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
1914: DMNetworkGetNumVariables(dm,edges[e],&ncols);
1916: if (network->Jv) {
1917: Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */
1918: } else Juser = NULL;
1919: MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost,vd_nz,vo_nz);
1921: /* Connected vertices */
1922: DMNetworkGetConnectedVertices(dm,edges[e],&cone);
1923: vc = (v == cone[0]) ? cone[1]:cone[0];
1924: DMNetworkIsGhostVertex(dm,vc,&ghost_vc);
1926: DMNetworkGetNumVariables(dm,vc,&ncols);
1928: if (network->Jv) {
1929: Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
1930: } else Juser = NULL;
1931: if (ghost_vc||ghost) {
1932: ghost2 = PETSC_TRUE;
1933: } else {
1934: ghost2 = PETSC_FALSE;
1935: }
1936: MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost2,vd_nz,vo_nz);
1937: }
1939: /* Set preallocation for vertex self */
1940: DMNetworkIsGhostVertex(dm,v,&ghost);
1941: if (!ghost) {
1942: DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
1943: if (network->Jv) {
1944: Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
1945: } else Juser = NULL;
1946: MatSetPreallocationblock_private(Juser,nrows,rows_v,nrows,PETSC_FALSE,vd_nz,vo_nz);
1947: }
1948: if (ghost) {
1949: PetscFree(rows_v);
1950: }
1951: }
1953: VecAssemblyBegin(vd_nz);
1954: VecAssemblyBegin(vo_nz);
1956: PetscMalloc2(localSize,&dnnz,localSize,&onnz);
1958: VecAssemblyEnd(vd_nz);
1959: VecAssemblyEnd(vo_nz);
1961: VecGetArray(vd_nz,&vdnz);
1962: VecGetArray(vo_nz,&vonz);
1963: for (j=0; j<localSize; j++) {
1964: dnnz[j] = (PetscInt)PetscRealPart(vdnz[j]);
1965: onnz[j] = (PetscInt)PetscRealPart(vonz[j]);
1966: }
1967: VecRestoreArray(vd_nz,&vdnz);
1968: VecRestoreArray(vo_nz,&vonz);
1969: VecDestroy(&vd_nz);
1970: VecDestroy(&vo_nz);
1972: MatSeqAIJSetPreallocation(*J,0,dnnz);
1973: MatMPIAIJSetPreallocation(*J,0,dnnz,0,onnz);
1974: MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1976: PetscFree2(dnnz,onnz);
1978: /* (2) Set matrix entries for edges */
1979: /*----------------------------------*/
1980: for (e=eStart; e<eEnd; e++) {
1981: /* Get row indices */
1982: DMNetworkGetVariableGlobalOffset(dm,e,&rstart);
1983: DMNetworkGetNumVariables(dm,e,&nrows);
1984: if (nrows) {
1985: for (j=0; j<nrows; j++) rows[j] = j + rstart;
1987: /* Set matrix entries for conntected vertices */
1988: DMNetworkGetConnectedVertices(dm,e,&cone);
1989: for (v=0; v<2; v++) {
1990: DMNetworkGetVariableGlobalOffset(dm,cone[v],&cstart);
1991: DMNetworkGetNumVariables(dm,cone[v],&ncols);
1993: if (network->Je) {
1994: Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
1995: } else Juser = NULL;
1996: MatSetblock_private(Juser,nrows,rows,ncols,cstart,J);
1997: }
1999: /* Set matrix entries for edge self */
2000: cstart = rstart;
2001: if (network->Je) {
2002: Juser = network->Je[3*e]; /* Jacobian(e,e) */
2003: } else Juser = NULL;
2004: MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);
2005: }
2006: }
2008: /* Set matrix entries for vertices */
2009: /*---------------------------------*/
2010: for (v=vStart; v<vEnd; v++) {
2011: /* Get row indices */
2012: DMNetworkGetVariableGlobalOffset(dm,v,&rstart);
2013: DMNetworkGetNumVariables(dm,v,&nrows);
2014: if (!nrows) continue;
2016: DMNetworkIsGhostVertex(dm,v,&ghost);
2017: if (ghost) {
2018: PetscMalloc1(nrows,&rows_v);
2019: } else {
2020: rows_v = rows;
2021: }
2022: for (j=0; j<nrows; j++) rows_v[j] = j + rstart;
2024: /* Get supporting edges and connected vertices */
2025: DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
2027: for (e=0; e<nedges; e++) {
2028: /* Supporting edges */
2029: DMNetworkGetVariableGlobalOffset(dm,edges[e],&cstart);
2030: DMNetworkGetNumVariables(dm,edges[e],&ncols);
2032: if (network->Jv) {
2033: Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */
2034: } else Juser = NULL;
2035: MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);
2037: /* Connected vertices */
2038: DMNetworkGetConnectedVertices(dm,edges[e],&cone);
2039: vc = (v == cone[0]) ? cone[1]:cone[0];
2041: DMNetworkGetVariableGlobalOffset(dm,vc,&cstart);
2042: DMNetworkGetNumVariables(dm,vc,&ncols);
2044: if (network->Jv) {
2045: Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
2046: } else Juser = NULL;
2047: MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);
2048: }
2050: /* Set matrix entries for vertex self */
2051: if (!ghost) {
2052: DMNetworkGetVariableGlobalOffset(dm,v,&cstart);
2053: if (network->Jv) {
2054: Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
2055: } else Juser = NULL;
2056: MatSetblock_private(Juser,nrows,rows_v,nrows,cstart,J);
2057: }
2058: if (ghost) {
2059: PetscFree(rows_v);
2060: }
2061: }
2062: PetscFree(rows);
2064: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
2065: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
2067: MatSetDM(*J,dm);
2068: return(0);
2069: }
2071: PetscErrorCode DMDestroy_Network(DM dm)
2072: {
2074: DM_Network *network = (DM_Network*)dm->data;
2075: PetscInt j;
2078: if (--network->refct > 0) return(0);
2079: if (network->Je) {
2080: PetscFree(network->Je);
2081: }
2082: if (network->Jv) {
2083: PetscFree(network->Jvptr);
2084: PetscFree(network->Jv);
2085: }
2087: ISLocalToGlobalMappingDestroy(&network->vertex.mapping);
2088: PetscSectionDestroy(&network->vertex.DofSection);
2089: PetscSectionDestroy(&network->vertex.GlobalDofSection);
2090: if (network->vltog) {
2091: PetscFree(network->vltog);
2092: }
2093: if (network->vertex.sf) {
2094: PetscSFDestroy(&network->vertex.sf);
2095: }
2096: /* edge */
2097: ISLocalToGlobalMappingDestroy(&network->edge.mapping);
2098: PetscSectionDestroy(&network->edge.DofSection);
2099: PetscSectionDestroy(&network->edge.GlobalDofSection);
2100: if (network->edge.sf) {
2101: PetscSFDestroy(&network->edge.sf);
2102: }
2103: DMDestroy(&network->plex);
2104: PetscSectionDestroy(&network->DataSection);
2105: PetscSectionDestroy(&network->DofSection);
2107: for (j=0; j<network->nsubnet; j++) {
2108: PetscFree(network->subnet[j].edges);
2109: }
2110: PetscFree(network->subnetvtx);
2112: PetscFree(network->subnet);
2113: PetscFree(network->componentdataarray);
2114: PetscFree2(network->header,network->cvalue);
2115: PetscFree(network);
2116: return(0);
2117: }
2119: PetscErrorCode DMView_Network(DM dm,PetscViewer viewer)
2120: {
2122: DM_Network *network = (DM_Network*)dm->data;
2123: PetscBool iascii;
2124: PetscMPIInt rank;
2125: PetscInt p,nsubnet;
2128: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
2129: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
2132: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2133: if (iascii) {
2134: const PetscInt *cone,*vtx,*edges;
2135: PetscInt vfrom,vto,i,j,nv,ne;
2137: nsubnet = network->nsubnet - network->ncsubnet; /* num of subnetworks */
2138: PetscViewerASCIIPushSynchronized(viewer);
2139: PetscViewerASCIISynchronizedPrintf(viewer, " [%d] nsubnet: %D; nsubnetCouple: %D; nEdges: %D; nVertices: %D\n",rank,nsubnet,network->ncsubnet,network->nEdges,network->nVertices);
2141: for (i=0; i<nsubnet; i++) {
2142: DMNetworkGetSubnetworkInfo(dm,i,&nv,&ne,&vtx,&edges);
2143: if (ne) {
2144: PetscViewerASCIISynchronizedPrintf(viewer, " Subnet %D: nEdges %D, nVertices %D\n",i,ne,nv);
2145: for (j=0; j<ne; j++) {
2146: p = edges[j];
2147: DMNetworkGetConnectedVertices(dm,p,&cone);
2148: DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);
2149: DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);
2150: DMNetworkGetGlobalEdgeIndex(dm,edges[j],&p);
2151: PetscViewerASCIISynchronizedPrintf(viewer, " edge %D: %D----> %D\n",p,vfrom,vto);
2152: }
2153: }
2154: }
2155: /* Coupling subnets */
2156: nsubnet = network->nsubnet;
2157: for (; i<nsubnet; i++) {
2158: DMNetworkGetSubnetworkInfo(dm,i,&nv,&ne,&vtx,&edges);
2159: if (ne) {
2160: PetscViewerASCIISynchronizedPrintf(viewer, " Subnet %D (couple): nEdges %D, nVertices %D\n",i,ne,nv);
2161: for (j=0; j<ne; j++) {
2162: p = edges[j];
2163: DMNetworkGetConnectedVertices(dm,p,&cone);
2164: DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);
2165: DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);
2166: DMNetworkGetGlobalEdgeIndex(dm,edges[j],&p);
2167: PetscViewerASCIISynchronizedPrintf(viewer, " edge %D: %D----> %D\n",p,vfrom,vto);
2168: }
2169: }
2170: }
2171: PetscViewerFlush(viewer);
2172: PetscViewerASCIIPopSynchronized(viewer);
2173: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMNetwork writing", ((PetscObject)viewer)->type_name);
2174: return(0);
2175: }
2177: PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l)
2178: {
2180: DM_Network *network = (DM_Network*)dm->data;
2183: DMGlobalToLocalBegin(network->plex,g,mode,l);
2184: return(0);
2185: }
2187: PetscErrorCode DMGlobalToLocalEnd_Network(DM dm, Vec g, InsertMode mode, Vec l)
2188: {
2190: DM_Network *network = (DM_Network*)dm->data;
2193: DMGlobalToLocalEnd(network->plex,g,mode,l);
2194: return(0);
2195: }
2197: PetscErrorCode DMLocalToGlobalBegin_Network(DM dm, Vec l, InsertMode mode, Vec g)
2198: {
2200: DM_Network *network = (DM_Network*)dm->data;
2203: DMLocalToGlobalBegin(network->plex,l,mode,g);
2204: return(0);
2205: }
2207: PetscErrorCode DMLocalToGlobalEnd_Network(DM dm, Vec l, InsertMode mode, Vec g)
2208: {
2210: DM_Network *network = (DM_Network*)dm->data;
2213: DMLocalToGlobalEnd(network->plex,l,mode,g);
2214: return(0);
2215: }
2217: /*@
2218: DMNetworkGetVertexLocalToGlobalOrdering - Get vertex global index
2220: Not collective
2222: Input Parameters:
2223: + dm - the dm object
2224: - vloc - local vertex ordering, start from 0
2226: Output Parameters:
2227: . vg - global vertex ordering, start from 0
2229: Level: advanced
2231: .seealso: DMNetworkSetVertexLocalToGlobalOrdering()
2232: @*/
2233: PetscErrorCode DMNetworkGetVertexLocalToGlobalOrdering(DM dm,PetscInt vloc,PetscInt *vg)
2234: {
2235: DM_Network *network = (DM_Network*)dm->data;
2236: PetscInt *vltog = network->vltog;
2239: if (!vltog) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkSetVertexLocalToGlobalOrdering() first");
2240: *vg = vltog[vloc];
2241: return(0);
2242: }
2244: /*@
2245: DMNetworkSetVertexLocalToGlobalOrdering - Create and setup vertex local to global map
2247: Collective
2249: Input Parameters:
2250: . dm - the dm object
2252: Level: advanced
2254: .seealso: DMNetworkGetGlobalVertexIndex()
2255: @*/
2256: PetscErrorCode DMNetworkSetVertexLocalToGlobalOrdering(DM dm)
2257: {
2258: PetscErrorCode ierr;
2259: DM_Network *network=(DM_Network*)dm->data;
2260: MPI_Comm comm;
2261: PetscMPIInt rank,size,*displs,*recvcounts,remoterank;
2262: PetscBool ghost;
2263: PetscInt *vltog,nroots,nleaves,i,*vrange,k,N,lidx;
2264: const PetscSFNode *iremote;
2265: PetscSF vsf;
2266: Vec Vleaves,Vleaves_seq;
2267: VecScatter ctx;
2268: PetscScalar *varr,val;
2269: const PetscScalar *varr_read;
2272: PetscObjectGetComm((PetscObject)dm,&comm);
2273: MPI_Comm_size(comm,&size);
2274: MPI_Comm_rank(comm,&rank);
2276: if (size == 1) {
2277: nroots = network->vEnd - network->vStart;
2278: PetscMalloc1(nroots, &vltog);
2279: for (i=0; i<nroots; i++) vltog[i] = i;
2280: network->vltog = vltog;
2281: return(0);
2282: }
2284: if (!network->distributecalled) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkDistribute() first");
2285: if (network->vltog) {
2286: PetscFree(network->vltog);
2287: }
2289: DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);
2290: PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);
2291: vsf = network->vertex.sf;
2293: PetscMalloc3(size+1,&vrange,size+1,&displs,size,&recvcounts);
2294: PetscSFGetGraph(vsf,&nroots,&nleaves,NULL,&iremote);
2296: for (i=0; i<size; i++) { displs[i] = i; recvcounts[i] = 1;}
2298: i = nroots - nleaves; /* local number of vertices, excluding ghosts */
2299: vrange[0] = 0;
2300: MPI_Allgatherv(&i,1,MPIU_INT,vrange+1,recvcounts,displs,MPIU_INT,comm);
2301: for (i=2; i<size+1; i++) {vrange[i] += vrange[i-1];}
2303: PetscMalloc1(nroots, &vltog);
2304: network->vltog = vltog;
2306: /* Set vltog for non-ghost vertices */
2307: k = 0;
2308: for (i=0; i<nroots; i++) {
2309: DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);
2310: if (ghost) continue;
2311: vltog[i] = vrange[rank] + k++;
2312: }
2313: PetscFree3(vrange,displs,recvcounts);
2315: /* Set vltog for ghost vertices */
2316: /* (a) create parallel Vleaves and sequential Vleaves_seq to convert local iremote[*].index to global index */
2317: VecCreate(comm,&Vleaves);
2318: VecSetSizes(Vleaves,2*nleaves,PETSC_DETERMINE);
2319: VecSetFromOptions(Vleaves);
2320: VecGetArray(Vleaves,&varr);
2321: for (i=0; i<nleaves; i++) {
2322: varr[2*i] = (PetscScalar)(iremote[i].rank); /* rank of remote process */
2323: varr[2*i+1] = (PetscScalar)(iremote[i].index); /* local index in remote process */
2324: }
2325: VecRestoreArray(Vleaves,&varr);
2327: /* (b) scatter local info to remote processes via VecScatter() */
2328: VecScatterCreateToAll(Vleaves,&ctx,&Vleaves_seq);
2329: VecScatterBegin(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);
2330: VecScatterEnd(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);
2332: /* (c) convert local indices to global indices in parallel vector Vleaves */
2333: VecGetSize(Vleaves_seq,&N);
2334: VecGetArrayRead(Vleaves_seq,&varr_read);
2335: for (i=0; i<N; i+=2) {
2336: remoterank = (PetscMPIInt)PetscRealPart(varr_read[i]);
2337: if (remoterank == rank) {
2338: k = i+1; /* row number */
2339: lidx = (PetscInt)PetscRealPart(varr_read[i+1]);
2340: val = (PetscScalar)vltog[lidx]; /* global index for non-ghost vertex computed above */
2341: VecSetValues(Vleaves,1,&k,&val,INSERT_VALUES);
2342: }
2343: }
2344: VecRestoreArrayRead(Vleaves_seq,&varr_read);
2345: VecAssemblyBegin(Vleaves);
2346: VecAssemblyEnd(Vleaves);
2348: /* (d) Set vltog for ghost vertices by copying local values of Vleaves */
2349: VecGetArrayRead(Vleaves,&varr_read);
2350: k = 0;
2351: for (i=0; i<nroots; i++) {
2352: DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);
2353: if (!ghost) continue;
2354: vltog[i] = (PetscInt)PetscRealPart(varr_read[2*k+1]); k++;
2355: }
2356: VecRestoreArrayRead(Vleaves,&varr_read);
2358: VecDestroy(&Vleaves);
2359: VecDestroy(&Vleaves_seq);
2360: VecScatterDestroy(&ctx);
2361: return(0);
2362: }