Actual source code: network.c
1: #include <petsc/private/dmnetworkimpl.h>
3: /*
4: Creates the component header and value objects for a network point
5: */
6: static PetscErrorCode SetUpNetworkHeaderComponentValue(DM dm,DMNetworkComponentHeader header,DMNetworkComponentValue cvalue)
7: {
11: /* Allocate arrays for component information */
12: PetscCalloc5(header->maxcomps,&header->size,header->maxcomps,&header->key,header->maxcomps,&header->offset,header->maxcomps,&header->nvar,header->maxcomps,&header->offsetvarrel);
13: PetscCalloc1(header->maxcomps,&cvalue->data);
15: /* The size of the header is the size of struct _p_DMNetworkComponentHeader. Since the struct contains PetscInt pointers we cannot use sizeof(struct). So, we need to explicitly calculate the size.
16: If the data header struct changes then this header size calculation needs to be updated. */
17: header->hsize = sizeof(struct _p_DMNetworkComponentHeader) + 5*header->maxcomps*sizeof(PetscInt);
18: header->hsize /= sizeof(DMNetworkComponentGenericDataType);
19: return(0);
20: }
22: /*@
23: DMNetworkGetPlex - Gets the Plex DM associated with this network DM
25: Not collective
27: Input Parameters:
28: . dm - the dm object
30: Output Parameters:
31: . plexdm - the plex dm object
33: Level: Advanced
35: .seealso: DMNetworkCreate()
36: @*/
37: PetscErrorCode DMNetworkGetPlex(DM dm,DM *plexdm)
38: {
39: DM_Network *network = (DM_Network*)dm->data;
42: *plexdm = network->plex;
43: return(0);
44: }
46: /*@
47: DMNetworkGetNumSubNetworks - Gets the the number of subnetworks
49: Not collective
51: Input Parameter:
52: . dm - the dm object
54: Output Parameters:
55: + nsubnet - local number of subnetworks
56: - Nsubnet - global number of subnetworks
58: Level: beginner
60: .seealso: DMNetworkCreate(), DMNetworkSetNumSubNetworks()
61: @*/
62: PetscErrorCode DMNetworkGetNumSubNetworks(DM dm,PetscInt *nsubnet,PetscInt *Nsubnet)
63: {
64: DM_Network *network = (DM_Network*)dm->data;
67: if (nsubnet) *nsubnet = network->nsubnet;
68: if (Nsubnet) *Nsubnet = network->Nsubnet;
69: return(0);
70: }
72: /*@
73: DMNetworkSetNumSubNetworks - Sets the number of subnetworks
75: Collective on dm
77: Input Parameters:
78: + dm - the dm object
79: . nsubnet - local number of subnetworks
80: - Nsubnet - global number of subnetworks
82: Level: beginner
84: .seealso: DMNetworkCreate(), DMNetworkGetNumSubNetworks()
85: @*/
86: PetscErrorCode DMNetworkSetNumSubNetworks(DM dm,PetscInt nsubnet,PetscInt Nsubnet)
87: {
89: DM_Network *network = (DM_Network*)dm->data;
92: if (network->Nsubnet != 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_INCOMP,"Network sizes alread set, cannot resize the network");
98: if (Nsubnet == PETSC_DECIDE) {
99: if (nsubnet < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of local subnetworks %D cannot be less than 0",nsubnet);
100: MPIU_Allreduce(&nsubnet,&Nsubnet,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
101: }
102: if (Nsubnet < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_INCOMP,"Number of global subnetworks %D cannot be less than 1",Nsubnet);
104: network->Nsubnet = Nsubnet;
105: network->nsubnet = 0; /* initia value; will be determind by DMNetworkAddSubnetwork() */
106: PetscCalloc1(Nsubnet,&network->subnet);
108: /* num of shared vertices */
109: network->nsvtx = 0;
110: network->Nsvtx = 0;
111: return(0);
112: }
114: /*@
115: DMNetworkAddSubnetwork - Add a subnetwork
117: Collective on dm
119: Input Parameters:
120: + dm - the dm object
121: . name - name of the subnetwork
122: . ne - number of local edges of this subnetwork
123: - edgelist - list of edges for this subnetwork, this is a one dimensional array with pairs of entries being the two vertices (in global numbering of the vertices) of each edge,
124: $ [first vertex of first edge, second vertex of first edge, first vertex of second edge, second vertex of second edge, etc]
126: Output Parameters:
127: . netnum - global index of the subnetwork
129: Notes:
130: There is no copy involved in this operation, only the pointer is referenced. The edgelist should
131: not be destroyed before the call to DMNetworkLayoutSetUp()
133: A network can comprise of a single subnetwork OR multiple subnetworks. For a single subnetwork, the subnetwork can be read either in serial or parallel. For a multiple subnetworks,
134: each subnetwork topology needs to be set on a unique rank and the communicator size needs to be at least equal to the number of subnetworks.
136: Level: beginner
138: Example usage:
139: Consider the following networks:
140: 1) A sigle subnetwork:
141: .vb
142: network 0:
143: rank[0]:
144: v0 --> v2; v1 --> v2
145: rank[1]:
146: v3 --> v5; v4 --> v5
147: .ve
149: The resulting input
150: network 0:
151: rank[0]:
152: ne = 2
153: edgelist = [0 2 | 1 2]
154: rank[1]:
155: ne = 2
156: edgelist = [3 5 | 4 5]
158: 2) Two subnetworks:
159: .vb
160: subnetwork 0:
161: rank[0]:
162: v0 --> v2; v2 --> v1; v1 --> v3;
163: subnetwork 1:
164: rank[1]:
165: v0 --> v3; v3 --> v2; v2 --> v1;
166: .ve
168: The resulting input
169: subnetwork 0:
170: rank[0]:
171: ne = 3
172: edgelist = [0 2 | 2 1 | 1 3]
173: rank[1]:
174: ne = 0
175: edgelist = NULL
177: subnetwork 1:
178: rank[0]:
179: ne = 0
180: edgelist = NULL
181: rank[1]:
182: edgelist = [0 3 | 3 2 | 2 1]
184: .seealso: DMNetworkCreate(), DMNetworkSetNumSubnetworks()
185: @*/
186: PetscErrorCode DMNetworkAddSubnetwork(DM dm,const char* name,PetscInt ne,PetscInt edgelist[],PetscInt *netnum)
187: {
189: DM_Network *network = (DM_Network*)dm->data;
190: PetscInt i,Nedge,j,Nvtx,nvtx;
191: PetscBT table;
194: for (i=0; i<ne; i++) {
195: if (edgelist[2*i] == edgelist[2*i+1]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Edge %D has the same vertex %D at each endpoint",i,edgelist[2*i]);
196: }
197: /* Get global total Nvtx = max(edgelist[])+1 for this subnet */
198: nvtx = -1; i = 0;
199: for (j=0; j<ne; j++) {
200: nvtx = PetscMax(nvtx, edgelist[i]); i++;
201: nvtx = PetscMax(nvtx, edgelist[i]); i++;
202: }
203: nvtx++;
204: MPIU_Allreduce(&nvtx,&Nvtx,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
206: /* Get local nvtx for this subnet */
207: PetscBTCreate(Nvtx,&table);
208: PetscBTMemzero(Nvtx,table);
209: i = 0;
210: for (j=0; j<ne; j++) {
211: PetscBTSet(table,edgelist[i]);
212: i++;
213: PetscBTSet(table,edgelist[i]);
214: i++;
215: }
216: nvtx = 0;
217: for (j=0; j<Nvtx; j++) {
218: if (PetscBTLookup(table,j)) nvtx++;
219: }
220: PetscBTDestroy(&table);
222: /* Get global total Nedge for this subnet */
223: MPIU_Allreduce(&ne,&Nedge,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
225: i = network->nsubnet;
226: if (name) {
227: PetscStrcpy(network->subnet[i].name,name);
228: }
229: network->subnet[i].nvtx = nvtx; /* include ghost vertices */
230: network->subnet[i].nedge = ne;
231: network->subnet[i].edgelist = edgelist;
232: network->subnet[i].Nvtx = Nvtx;
233: network->subnet[i].Nedge = Nedge;
235: /* ----------------------------------------------------------
236: p=v or e;
237: subnet[0].pStart = 0
238: subnet[i+1].pStart = subnet[i].pEnd = subnet[i].pStart + (nE[i] or NV[i])
239: ----------------------------------------------------------------------- */
240: /* GLOBAL subnet[].vStart and vEnd, used by DMNetworkLayoutSetUp() */
241: network->subnet[i].vStart = network->NVertices;
242: network->subnet[i].vEnd = network->subnet[i].vStart + network->subnet[i].Nvtx; /* global vEnd of subnet[i] */
244: network->nVertices += nvtx; /* include ghost vertices */
245: network->NVertices += network->subnet[i].Nvtx;
247: /* LOCAL subnet[].eStart and eEnd, used by DMNetworkLayoutSetUp() */
248: network->subnet[i].eStart = network->nEdges;
249: network->subnet[i].eEnd = network->subnet[i].eStart + ne;
250: network->nEdges += ne;
251: network->NEdges += network->subnet[i].Nedge;
253: PetscStrcpy(network->subnet[i].name,name);
254: if (netnum) *netnum = network->nsubnet;
255: network->nsubnet++;
256: return(0);
257: }
259: /*
260: SetUp a single svtx struct. See SVtx defined in dmnetworkimpl.h
261: Set gidx and type if input v=(net,idx) is a from_vertex;
262: Get gid, type and index in the svtx array if input v=(net,idx) is a to_vertex.
264: Input: Nsvtx, svtx, net, idx, gidx
265: Output: gidx, svtype, svtx_idx
266: */
267: static PetscErrorCode SVtxSetUp(PetscInt Nsvtx,SVtx *svtx,PetscInt net,PetscInt idx,PetscInt *gidx,SVtxType *svtype,PetscInt *svtx_idx)
268: {
269: PetscInt i,j,*svto;
270: SVtxType vtype;
273: if (!Nsvtx) return(0);
275: vtype = SVNONE;
276: for (i=0; i<Nsvtx; i++) {
277: if (net == svtx[i].sv[0] && idx == svtx[i].sv[1]) {
278: /* (1) input vertex net.idx is a shared from_vertex, set its global index and output its svtype */
279: svtx[i].gidx = *gidx; /* set gidx */
280: vtype = SVFROM;
281: } else { /* loop over svtx[i].n */
282: for (j=1; j<svtx[i].n; j++) {
283: svto = svtx[i].sv + 2*j;
284: if (net == svto[0] && idx == svto[1]) {
285: /* input vertex net.idx is a shared to_vertex, output its global index and its svtype */
286: *gidx = svtx[i].gidx; /* output gidx for to_vertex */
287: vtype = SVTO;
288: }
289: }
290: }
291: if (vtype != SVNONE) break;
292: }
293: if (svtype) *svtype = vtype;
294: if (svtx_idx) *svtx_idx = i;
295: return(0);
296: }
298: /*
299: Add a new shared vertex sv=(net,idx) to table svtas[ita]
300: */
301: static PetscErrorCode TableAddSVtx(PetscTable *svtas,PetscInt ita,PetscInt* tdata,PetscInt *sv_wk,PetscInt *ii,PetscInt *sedgelist,PetscInt k,DM_Network *network,PetscInt **ta2sv)
302: {
303: PetscInt net,idx,gidx,i=*ii;
307: net = sv_wk[2*i] = sedgelist[k];
308: idx = sv_wk[2*i+1] = sedgelist[k+1];
309: gidx = network->subnet[net].vStart + idx;
310: PetscTableAdd(svtas[ita],gidx+1,tdata[ita]+1,INSERT_VALUES);
311: *(ta2sv[ita] + tdata[ita]) = i; /* maps tdata to index of sv_wk; sv_wk keeps (net,idx) info */
312: tdata[ita]++; (*ii)++;
313: return(0);
314: }
316: /*
317: Create an array of shared vertices. See SVtx and SVtxType in dmnetworkimpl.h
319: Input: dm, Nsedgelist, sedgelist
320: Output: Nsvtx,svtx
322: Note: Output svtx is organized as
323: sv(net[0],idx[0]) --> sv(net[1],idx[1])
324: --> sv(net[1],idx[1])
325: ...
326: --> sv(net[n-1],idx[n-1])
327: and net[0] < net[1] < ... < net[n-1]
328: where sv[0] has SVFROM type, sv[i], i>0, has SVTO type.
329: */
330: static PetscErrorCode SVtxCreate(DM dm,PetscInt Nsedgelist,PetscInt *sedgelist,PetscInt *Nsvtx,SVtx **svtx)
331: {
333: SVtx *sedges = NULL;
334: PetscInt *sv,k,j,nsv,*tdata,**ta2sv;
335: PetscTable *svtas;
336: PetscInt gidx,net,idx,i,nta,ita,idx_from,idx_to,n,*sv_wk;
337: DM_Network *network = (DM_Network*)dm->data;
338: PetscTablePosition ppos;
341: /* (1) Crete ctables svtas */
342: PetscCalloc4(Nsedgelist,&svtas,Nsedgelist,&tdata,4*Nsedgelist,&sv_wk,2*Nsedgelist,&ta2sv);
344: j = 0; /* sedgelist counter */
345: k = 0; /* sedgelist vertex counter j = 4*k */
346: i = 0; /* sv_wk (vertices added to the ctables) counter */
347: nta = 0; /* num of sv tables created */
349: /* for j=0 */
350: PetscTableCreate(2*Nsedgelist,network->NVertices+1,&svtas[nta]);
351: PetscMalloc1(2*Nsedgelist,&ta2sv[nta]);
353: TableAddSVtx(svtas,nta,tdata,sv_wk,&i,sedgelist,k,network,ta2sv);
354: TableAddSVtx(svtas,nta,tdata,sv_wk,&i,sedgelist,k+2,network,ta2sv);
355: nta++; k += 4;
357: for (j = 1; j < Nsedgelist; j++) {
358: for (ita = 0; ita < nta; ita++) {
359: /* vfrom */
360: net = sedgelist[k]; idx = sedgelist[k+1];
361: gidx = network->subnet[net].vStart + idx; /* global index of the vertex net.idx before merging shared vertices */
362: PetscTableFind(svtas[ita],gidx+1,&idx_from);
364: /* vto */
365: net = sedgelist[k+2]; idx = sedgelist[k+3];
366: gidx = network->subnet[net].vStart + idx;
367: PetscTableFind(svtas[ita],gidx+1,&idx_to);
369: if (idx_from || idx_to) { /* vfrom or vto is on table svtas[ita] */
370: idx_from--; idx_to--;
371: if (idx_from < 0) { /* vto is on svtas[ita] */
372: TableAddSVtx(svtas,ita,tdata,sv_wk,&i,sedgelist,k,network,ta2sv);
373: break;
374: } else if (idx_to < 0) {
375: TableAddSVtx(svtas,ita,tdata,sv_wk,&i,sedgelist,k+2,network,ta2sv);
376: break;
377: }
378: }
379: }
381: if (ita == nta) {
382: PetscTableCreate(2*Nsedgelist,network->NVertices+1,&svtas[nta]);
383: PetscMalloc1(2*Nsedgelist, &ta2sv[nta]);
385: TableAddSVtx(svtas,nta,tdata,sv_wk,&i,sedgelist,k,network,ta2sv);
386: TableAddSVtx(svtas,nta,tdata,sv_wk,&i,sedgelist,k+2,network,ta2sv);
387: nta++;
388: }
389: k += 4;
390: }
392: /* (2) Construct sedges from ctable
393: sedges: edges connect vertex sv[0]=(net[0],idx[0]) to vertices sv[k], k=1,...,n-1;
394: net[k], k=0, ...,n-1, are in ascending order */
395: PetscMalloc1(nta,&sedges);
396: for (nsv = 0; nsv < nta; nsv++) {
397: /* for a single svtx, put shared vertices in ascending order of gidx */
398: PetscTableGetCount(svtas[nsv],&n);
399: PetscCalloc1(2*n,&sv);
400: sedges[nsv].sv = sv;
401: sedges[nsv].n = n;
402: sedges[nsv].gidx = -1; /* initialization */
404: PetscTableGetHeadPosition(svtas[nsv],&ppos);
405: for (k=0; k<n; k++) { /* gidx is sorted in ascending order */
406: PetscTableGetNext(svtas[nsv],&ppos,&gidx,&i);
407: gidx--; i--;
409: j = ta2sv[nsv][i]; /* maps i to index of sv_wk */
410: sv[2*k] = sv_wk[2*j];
411: sv[2*k+1] = sv_wk[2*j + 1];
412: }
413: }
415: for (j=0; j<nta; j++) {
416: PetscTableDestroy(&svtas[j]);
417: PetscFree(ta2sv[j]);
418: }
419: PetscFree4(svtas,tdata,sv_wk,ta2sv);
421: *Nsvtx = nta;
422: *svtx = sedges;
423: return(0);
424: }
426: static PetscErrorCode DMNetworkLayoutSetUp_Coupling(DM dm)
427: {
429: DM_Network *network = (DM_Network*)dm->data;
430: PetscInt i,j,ctr,np,*edges,*subnetvtx,*subnetedge,vStart;
431: PetscInt k,*vidxlTog,Nsv=0,Nsubnet=network->Nsubnet;
432: PetscInt *sedgelist=network->sedgelist;
433: const PetscInt *cone;
434: MPI_Comm comm;
435: PetscMPIInt size,rank,*recvcounts=NULL,*displs=NULL;
436: PetscInt net,idx,gidx,nmerged,e,v,vfrom,vto,*vrange,*eowners,gidx_from,net_from,sv_idx;
437: SVtxType svtype = SVNONE;
438: SVtx *svtx=NULL;
439: PetscSection sectiong;
442: /* This implementation requires user input each subnet by a single processor, thus subnet[net].nvtx=subnet[net].Nvtx */
443: for (net=0; net<Nsubnet; net++) {
444: if (network->subnet[net].nvtx && network->subnet[net].nvtx != network->subnet[net].Nvtx) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"subnetwork %D local num of vertices %D != %D global num",net,network->subnet[net].nvtx,network->subnet[net].Nvtx);
445: }
447: PetscObjectGetComm((PetscObject)dm,&comm);
448: MPI_Comm_rank(comm,&rank);
449: MPI_Comm_size(comm,&size);
451: /* (1) Create svtx[] from sedgelist */
452: /* -------------------------------- */
453: /* Nsv: global number of SVtx; svtx: shared vertices, see SVtx in dmnetworkimpl.h */
454: SVtxCreate(dm,network->Nsvtx,sedgelist,&Nsv,&svtx);
456: /* (2) Setup svtx; Shared vto vertices are merged to their vfrom vertex with same global vetex index (gidx) */
457: /* -------------------------------------------------------------------------------------------------------- */
458: /* (2.1) compute vrage[rank]: global index of 1st local vertex in proc[rank] */
459: PetscMalloc3(size+1,&vrange,size,&displs,size,&recvcounts);
460: for (i=0; i<size; i++) {displs[i] = i; recvcounts[i] = 1;}
462: vrange[0] = 0;
463: MPI_Allgatherv(&network->nVertices,1,MPIU_INT,vrange+1,recvcounts,displs,MPIU_INT,comm);
464: for (i=2; i<size+1; i++) {
465: vrange[i] += vrange[i-1];
466: }
468: /* (2.2) Create vidxlTog: maps UN-MERGED local vertex index i to global index gidx (plex, excluding ghost vertices) */
469: PetscMalloc1(network->nVertices,&vidxlTog);
470: i = 0; gidx = 0;
471: nmerged = 0; /* local num of merged vertices */
472: network->nsvtx = 0;
473: for (net=0; net<Nsubnet; net++) {
474: for (idx=0; idx<network->subnet[net].Nvtx; idx++) {
475: gidx_from = gidx;
476: sv_idx = -1;
478: SVtxSetUp(Nsv,svtx,net,idx,&gidx_from,&svtype,&sv_idx);
479: if (svtype == SVTO) {
480: if (network->subnet[net].nvtx) {/* this proc owns sv_to */
481: net_from = svtx[sv_idx].sv[0]; /* subnet num of its shared vertex */
482: if (network->subnet[net_from].nvtx == 0) {
483: /* this proc does not own v_from, thus a new local coupling vertex */
484: network->nsvtx++;
485: }
486: vidxlTog[i++] = gidx_from;
487: nmerged++; /* a coupling vertex -- merged */
488: }
489: } else {
490: if (svtype == SVFROM) {
491: if (network->subnet[net].nvtx) {
492: /* this proc owns this v_from, a new local coupling vertex */
493: network->nsvtx++;
494: }
495: }
496: if (network->subnet[net].nvtx) vidxlTog[i++] = gidx;
497: gidx++;
498: }
499: }
500: }
501: #if defined(PETSC_USE_DEBUG)
502: if (i != network->nVertices) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"%D != %D nVertices",i,network->nVertices);
503: #endif
505: /* (2.3) Setup svtable for querry shared vertices */
506: for (v=0; v<Nsv; v++) {
507: gidx = svtx[v].gidx;
508: PetscTableAdd(network->svtable,gidx+1,v+1,INSERT_VALUES);
509: }
511: /* (2.4) Shared vertices in the subnetworks are merged, update global NVertices: np = sum(local nmerged) */
512: MPI_Allreduce(&nmerged,&np,1,MPIU_INT,MPI_SUM,comm);
513: network->NVertices -= np;
515: PetscCalloc1(2*network->nEdges,&edges);
517: ctr = 0;
518: for (net=0; net<Nsubnet; net++) {
519: for (j = 0; j < network->subnet[net].nedge; j++) {
520: /* vfrom: */
521: i = network->subnet[net].edgelist[2*j] + (network->subnet[net].vStart - vrange[rank]);
522: edges[2*ctr] = vidxlTog[i];
524: /* vto */
525: i = network->subnet[net].edgelist[2*j+1] + (network->subnet[net].vStart - vrange[rank]);
526: edges[2*ctr+1] = vidxlTog[i];
527: ctr++;
528: }
529: }
530: PetscFree3(vrange,displs,recvcounts);
531: PetscFree(vidxlTog);
533: /* (3) Create network->plex */
534: DMCreate(comm,&network->plex);
535: DMSetType(network->plex,DMPLEX);
536: DMSetDimension(network->plex,1);
537: if (size == 1) {
538: DMPlexBuildFromCellList(network->plex,network->nEdges,network->nVertices-nmerged,2,edges);
539: } else {
540: DMPlexBuildFromCellListParallel(network->plex,network->nEdges,network->nVertices-nmerged,PETSC_DECIDE,2,edges,NULL);
541: }
543: DMPlexGetChart(network->plex,&network->pStart,&network->pEnd);
544: DMPlexGetHeightStratum(network->plex,0,&network->eStart,&network->eEnd);
545: DMPlexGetHeightStratum(network->plex,1,&network->vStart,&network->vEnd);
546: vStart = network->vStart;
548: PetscSectionCreate(comm,&network->DataSection);
549: PetscSectionCreate(comm,&network->DofSection);
550: PetscSectionSetChart(network->DataSection,network->pStart,network->pEnd);
551: PetscSectionSetChart(network->DofSection,network->pStart,network->pEnd);
553: np = network->pEnd - network->pStart;
554: PetscCalloc2(np,&network->header,np,&network->cvalue);
555: for (i=0; i<np; i++) {
556: network->header[i].maxcomps = 1;
557: SetUpNetworkHeaderComponentValue(dm,&network->header[i],&network->cvalue[i]);
558: }
560: /* (4) Create edge and vertex arrays for the subnetworks */
561: PetscCalloc2(network->nEdges,&subnetedge,network->nVertices+network->nsvtx,&subnetvtx); /* Maps local edge/vertex to local subnetwork's edge/vertex */
562: network->subnetedge = subnetedge;
563: network->subnetvtx = subnetvtx;
564: for (j=0; j < Nsubnet; j++) {
565: network->subnet[j].edges = subnetedge;
566: subnetedge += network->subnet[j].nedge;
568: network->subnet[j].vertices = subnetvtx;
569: subnetvtx += network->subnet[j].nvtx;
570: }
571: network->svertices = subnetvtx;
573: /* Get edge ownership */
574: PetscMalloc1(size+1,&eowners);
575: np = network->eEnd - network->eStart; /* num of local edges */
576: MPI_Allgather(&np,1,MPIU_INT,eowners+1,1,MPIU_INT,comm);
577: eowners[0] = 0;
578: for (i=2; i<=size; i++) eowners[i] += eowners[i-1];
580: /* Setup edge and vertex arrays for subnetworks */
581: e = 0;
582: for (i=0; i < Nsubnet; i++) {
583: ctr = 0;
584: for (j = 0; j < network->subnet[i].nedge; j++) {
585: /* edge e */
586: network->header[e].index = e + eowners[rank]; /* Global edge index */
587: network->header[e].subnetid = i; /* Subnetwork id */
588: network->subnet[i].edges[j] = e;
589: network->header[e].ndata = 0;
590: network->header[e].offset[0] = 0;
591: network->header[e].offsetvarrel[0] = 0;
592: PetscSectionAddDof(network->DataSection,e,network->header[e].hsize);
594: /* connected vertices */
595: DMPlexGetCone(network->plex,e,&cone);
597: /* vertex cone[0] */
598: v = cone[0];
599: network->header[v].index = edges[2*e]; /* Global vertex index */
600: network->header[v].subnetid = i; /* Subnetwork id */
601: vfrom = network->subnet[i].edgelist[2*ctr];
602: network->subnet[i].vertices[vfrom] = v; /* user's subnet[].dix = petsc's v */
604: /* vertex cone[1] */
605: v = cone[1];
606: network->header[v].index = edges[2*e+1]; /* Global vertex index */
607: network->header[v].subnetid = i;
608: vto = network->subnet[i].edgelist[2*ctr+1];
609: network->subnet[i].vertices[vto]= v;
611: e++; ctr++;
612: }
613: }
614: PetscFree(eowners);
615: PetscFree(edges);
617: /* Set vertex array for the subnetworks */
618: k = 0;
619: for (v=vStart; v<network->vEnd; v++) { /* local vertices, including ghosts */
620: network->header[v].ndata = 0;
621: network->header[v].offset[0] = 0;
622: network->header[v].offsetvarrel[0] = 0;
623: PetscSectionAddDof(network->DataSection,v,network->header[v].hsize);
625: /* shared vertex */
626: PetscTableFind(network->svtable,network->header[v].index+1,&i);
627: if (i) network->svertices[k++] = v;
628: }
629: #if defined(PETSC_USE_DEBUG)
630: if (k != network->nsvtx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"k %D != %D nsvtx",k,network->nsvtx);
631: #endif
633: network->svtx = svtx;
634: network->Nsvtx = Nsv;
635: PetscFree(sedgelist);
637: /* Create a global section to be used by DMNetworkIsGhostVertex() which is a non-collective routine */
638: /* see snes_tutorials_network-ex1_4 */
639: DMGetGlobalSection(network->plex,§iong);
640: return(0);
641: }
643: /*@
644: DMNetworkLayoutSetUp - Sets up the bare layout (graph) for the network
646: Collective on dm
648: Input Parameters:
649: . dm - the dmnetwork object
651: Notes:
652: This routine should be called after the network sizes and edgelists have been provided. It creates
653: the bare layout of the network and sets up the network to begin insertion of components.
655: All the components should be registered before calling this routine.
657: Level: beginner
659: .seealso: DMNetworkSetNumSubNetworks(), DMNetworkAddSubnetwork()
660: @*/
661: PetscErrorCode DMNetworkLayoutSetUp(DM dm)
662: {
664: DM_Network *network = (DM_Network*)dm->data;
665: PetscInt i,j,ctr,Nsubnet=network->Nsubnet,*eowners,np,*edges,*subnetvtx,*subnetedge,e,v,vfrom,vto;
666: const PetscInt *cone;
667: MPI_Comm comm;
668: PetscMPIInt size,rank;
671: if (network->nsubnet != network->Nsubnet) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Must call DMNetworkAddSubnetwork() %D times",network->Nsubnet);
673: /* Create svtable for querry shared vertices */
674: PetscTableCreate(network->Nsvtx,network->NVertices+1,&network->svtable);
676: if (network->Nsvtx) {
677: DMNetworkLayoutSetUp_Coupling(dm);
678: return(0);
679: }
681: PetscObjectGetComm((PetscObject)dm,&comm);
682: MPI_Comm_rank(comm,&rank);
683: MPI_Comm_size(comm,&size);
685: /* Create LOCAL edgelist in global vertex ordering for the network by concatenating local input edgelists of the subnetworks */
686: PetscCalloc1(2*network->nEdges,&edges);
687: ctr = 0;
688: for (i=0; i < Nsubnet; i++) {
689: for (j = 0; j < network->subnet[i].nedge; j++) {
690: edges[2*ctr] = network->subnet[i].vStart + network->subnet[i].edgelist[2*j];
691: edges[2*ctr+1] = network->subnet[i].vStart + network->subnet[i].edgelist[2*j+1];
692: ctr++;
693: }
694: }
696: /* Create network->plex; One dimensional network, numCorners=2 */
697: DMCreate(comm,&network->plex);
698: DMSetType(network->plex,DMPLEX);
699: DMSetDimension(network->plex,1);
700: if (size == 1) {
701: DMPlexBuildFromCellList(network->plex,network->nEdges,network->nVertices,2,edges);
702: } else {
703: DMPlexBuildFromCellListParallel(network->plex,network->nEdges,network->nVertices,PETSC_DECIDE,2,edges,NULL);
704: }
706: DMPlexGetChart(network->plex,&network->pStart,&network->pEnd);
707: DMPlexGetHeightStratum(network->plex,0,&network->eStart,&network->eEnd);
708: DMPlexGetHeightStratum(network->plex,1,&network->vStart,&network->vEnd);
710: PetscSectionCreate(comm,&network->DataSection);
711: PetscSectionCreate(comm,&network->DofSection);
712: PetscSectionSetChart(network->DataSection,network->pStart,network->pEnd);
713: PetscSectionSetChart(network->DofSection,network->pStart,network->pEnd);
715: np = network->pEnd - network->pStart;
716: PetscCalloc2(np,&network->header,np,&network->cvalue);
717: for (i=0; i < np; i++) {
718: network->header[i].maxcomps = 1;
719: SetUpNetworkHeaderComponentValue(dm,&network->header[i],&network->cvalue[i]);
720: }
722: /* Create edge and vertex arrays for the subnetworks
723: This implementation assumes that DMNetwork reads
724: (1) a single subnetwork in parallel; or
725: (2) n subnetworks using n processors, one subnetwork/processor.
726: */
727: PetscCalloc2(network->nEdges,&subnetedge,network->nVertices,&subnetvtx); /* Maps local edge/vertex to local subnetwork's edge/vertex */
729: network->subnetedge = subnetedge;
730: network->subnetvtx = subnetvtx;
731: for (j=0; j < network->Nsubnet; j++) {
732: network->subnet[j].edges = subnetedge;
733: subnetedge += network->subnet[j].nedge;
735: network->subnet[j].vertices = subnetvtx;
736: subnetvtx += network->subnet[j].nvtx;
737: }
739: /* Get edge ownership */
740: PetscMalloc1(size+1,&eowners);
741: np = network->eEnd - network->eStart;
742: MPI_Allgather(&np,1,MPIU_INT,eowners+1,1,MPIU_INT,comm);
743: eowners[0] = 0;
744: for (i=2; i<=size; i++) eowners[i] += eowners[i-1];
746: /* Setup edge and vertex arrays for subnetworks */
747: e = 0;
748: for (i=0; i < Nsubnet; i++) {
749: ctr = 0;
750: for (j = 0; j < network->subnet[i].nedge; j++) {
751: /* edge e */
752: network->header[e].index = e + eowners[rank]; /* Global edge index */
753: network->header[e].subnetid = i;
754: network->subnet[i].edges[j] = e;
756: network->header[e].ndata = 0;
757: network->header[e].offset[0] = 0;
758: network->header[e].offsetvarrel[0] = 0;
759: PetscSectionAddDof(network->DataSection,e,network->header[e].hsize);
761: /* connected vertices */
762: DMPlexGetCone(network->plex,e,&cone);
764: /* vertex cone[0] */
765: v = cone[0];
766: network->header[v].index = edges[2*e]; /* Global vertex index */
767: network->header[v].subnetid = i; /* Subnetwork id */
768: if (Nsubnet == 1) {
769: network->subnet[i].vertices[v - network->vStart] = v; /* user's subnet[].idx = petsc's v */
770: } else {
771: vfrom = network->subnet[i].edgelist[2*ctr]; /* =subnet[i].idx, Global index! */
772: network->subnet[i].vertices[vfrom] = v; /* user's subnet[].dix = petsc's v */
773: }
775: /* vertex cone[1] */
776: v = cone[1];
777: network->header[v].index = edges[2*e+1]; /* Global vertex index */
778: network->header[v].subnetid = i; /* Subnetwork id */
779: if (Nsubnet == 1) {
780: network->subnet[i].vertices[v - network->vStart] = v; /* user's subnet[].idx = petsc's v */
781: } else {
782: vto = network->subnet[i].edgelist[2*ctr+1]; /* =subnet[i].idx, Global index! */
783: network->subnet[i].vertices[vto] = v; /* user's subnet[].dix = petsc's v */
784: }
786: e++; ctr++;
787: }
788: }
789: PetscFree(eowners);
790: PetscFree(edges); /* local edge list with global idx used by DMPlexBuildFromCellList() */
792: for (v = network->vStart; v < network->vEnd; v++) {
793: network->header[v].ndata = 0;
794: network->header[v].offset[0] = 0;
795: network->header[v].offsetvarrel[0] = 0;
796: PetscSectionAddDof(network->DataSection,v,network->header[e].hsize);
797: }
798: return(0);
799: }
801: /*@C
802: DMNetworkGetSubnetwork - Returns the information about a requested subnetwork
804: Not collective
806: Input Parameters:
807: + dm - the DM object
808: - netnum - the global index of the subnetwork
810: Output Parameters:
811: + nv - number of vertices (local)
812: . ne - number of edges (local)
813: . vtx - local vertices of the subnetwork
814: - edge - local edges of the subnetwork
816: Notes:
817: Cannot call this routine before DMNetworkLayoutSetup()
819: Level: intermediate
821: .seealso: DMNetworkCreate(), DMNetworkAddSubnetwork(), DMNetworkLayoutSetUp()
822: @*/
823: PetscErrorCode DMNetworkGetSubnetwork(DM dm,PetscInt netnum,PetscInt *nv,PetscInt *ne,const PetscInt **vtx,const PetscInt **edge)
824: {
825: DM_Network *network = (DM_Network*)dm->data;
828: if (netnum >= network->Nsubnet) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subnet index %D exceeds the num of subnets %D",netnum,network->Nsubnet);
829: if (nv) *nv = network->subnet[netnum].nvtx;
830: if (ne) *ne = network->subnet[netnum].nedge;
831: if (vtx) *vtx = network->subnet[netnum].vertices;
832: if (edge) *edge = network->subnet[netnum].edges;
833: return(0);
834: }
836: /*@
837: DMNetworkAddSharedVertices - Add shared vertices that connect two given subnetworks
839: Collective on dm
841: Input Parameters:
842: + dm - the dm object
843: . anetnum - first subnetwork global numbering returned by DMNetworkAddSubnetwork()
844: . bnetnum - second subnetwork global numbering returned by DMNetworkAddSubnetwork()
845: . nsvtx - number of vertices that are shared by the two subnetworks
846: . asvtx - vertex index in the first subnetwork
847: - bsvtx - vertex index in the second subnetwork
849: Level: beginner
851: .seealso: DMNetworkCreate(), DMNetworkAddSubnetwork(), DMNetworkGetSharedVertices()
852: @*/
853: PetscErrorCode DMNetworkAddSharedVertices(DM dm,PetscInt anetnum,PetscInt bnetnum,PetscInt nsvtx,PetscInt asvtx[],PetscInt bsvtx[])
854: {
856: DM_Network *network = (DM_Network*)dm->data;
857: PetscInt i,nsubnet = network->Nsubnet,*sedgelist,Nsvtx=network->Nsvtx;
860: if (anetnum == bnetnum) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Subnetworks must have different netnum");
861: if (anetnum < 0 || bnetnum < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"netnum cannot be negative");
862: if (!Nsvtx) {
863: /* allocate network->sedgelist to hold at most 2*nsubnet pairs of shared vertices */
864: PetscMalloc1(2*4*nsubnet,&network->sedgelist);
865: }
867: sedgelist = network->sedgelist;
868: for (i=0; i<nsvtx; i++) {
869: sedgelist[4*Nsvtx] = anetnum; sedgelist[4*Nsvtx+1] = asvtx[i];
870: sedgelist[4*Nsvtx+2] = bnetnum; sedgelist[4*Nsvtx+3] = bsvtx[i];
871: Nsvtx++;
872: }
873: if (Nsvtx > 2*nsubnet) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"allocate more space for coupling edgelist");
874: network->Nsvtx = Nsvtx;
875: return(0);
876: }
878: /*@C
879: DMNetworkGetSharedVertices - Returns the info for the shared vertices
881: Not collective
883: Input Parameter:
884: . dm - the DM object
886: Output Parameters:
887: + nsv - number of local shared vertices
888: - svtx - local shared vertices
890: Notes:
891: Cannot call this routine before DMNetworkLayoutSetup()
893: Level: intermediate
895: .seealso: DMNetworkGetSubnetwork(), DMNetworkLayoutSetUp(), DMNetworkAddSharedVertices()
896: @*/
897: PetscErrorCode DMNetworkGetSharedVertices(DM dm,PetscInt *nsv,const PetscInt **svtx)
898: {
899: DM_Network *net = (DM_Network*)dm->data;
902: if (net->Nsvtx) {
903: *nsv = net->nsvtx;
904: *svtx = net->svertices;
905: } else {
906: *nsv = 0;
907: *svtx = NULL;
908: }
909: return(0);
910: }
912: /*@C
913: DMNetworkRegisterComponent - Registers the network component
915: Logically collective on dm
917: Input Parameters:
918: + dm - the network object
919: . name - the component name
920: - size - the storage size in bytes for this component data
922: Output Parameters:
923: . key - an integer key that defines the component
925: Notes
926: This routine should be called by all processors before calling DMNetworkLayoutSetup().
928: Level: beginner
930: .seealso: DMNetworkCreate(), DMNetworkLayoutSetUp()
931: @*/
932: PetscErrorCode DMNetworkRegisterComponent(DM dm,const char *name,size_t size,PetscInt *key)
933: {
934: PetscErrorCode ierr;
935: DM_Network *network = (DM_Network*) dm->data;
936: DMNetworkComponent *component=NULL,*newcomponent=NULL;
937: PetscBool flg=PETSC_FALSE;
938: PetscInt i;
941: if (!network->component) {
942: PetscCalloc1(network->max_comps_registered,&network->component);
943: }
945: for (i=0; i < network->ncomponent; i++) {
946: PetscStrcmp(network->component[i].name,name,&flg);
947: if (flg) {
948: *key = i;
949: return(0);
950: }
951: }
953: if (network->ncomponent == network->max_comps_registered) {
954: /* Reached max allowed so resize component */
955: network->max_comps_registered += 2;
956: PetscCalloc1(network->max_comps_registered,&newcomponent);
957: /* Copy over the previous component info */
958: for (i=0; i < network->ncomponent; i++) {
959: PetscStrcpy(newcomponent[i].name,network->component[i].name);
960: newcomponent[i].size = network->component[i].size;
961: }
962: /* Free old one */
963: PetscFree(network->component);
964: /* Update pointer */
965: network->component = newcomponent;
966: }
968: component = &network->component[network->ncomponent];
970: PetscStrcpy(component->name,name);
971: component->size = size/sizeof(DMNetworkComponentGenericDataType);
972: *key = network->ncomponent;
973: network->ncomponent++;
974: return(0);
975: }
977: /*@
978: DMNetworkGetVertexRange - Get the bounds [start, end) for the local vertices
980: Not Collective
982: Input Parameter:
983: . dm - the DMNetwork object
985: Output Parameters:
986: + vStart - the first vertex point
987: - vEnd - one beyond the last vertex point
989: Level: beginner
991: .seealso: DMNetworkGetEdgeRange()
992: @*/
993: PetscErrorCode DMNetworkGetVertexRange(DM dm,PetscInt *vStart,PetscInt *vEnd)
994: {
995: DM_Network *network = (DM_Network*)dm->data;
998: if (vStart) *vStart = network->vStart;
999: if (vEnd) *vEnd = network->vEnd;
1000: return(0);
1001: }
1003: /*@
1004: DMNetworkGetEdgeRange - Get the bounds [start, end) for the local edges
1006: Not Collective
1008: Input Parameter:
1009: . dm - the DMNetwork object
1011: Output Parameters:
1012: + eStart - The first edge point
1013: - eEnd - One beyond the last edge point
1015: Level: beginner
1017: .seealso: DMNetworkGetVertexRange()
1018: @*/
1019: PetscErrorCode DMNetworkGetEdgeRange(DM dm,PetscInt *eStart,PetscInt *eEnd)
1020: {
1021: DM_Network *network = (DM_Network*)dm->data;
1024: if (eStart) *eStart = network->eStart;
1025: if (eEnd) *eEnd = network->eEnd;
1026: return(0);
1027: }
1029: static PetscErrorCode DMNetworkGetIndex(DM dm,PetscInt p,PetscInt *index)
1030: {
1031: PetscErrorCode ierr;
1032: DM_Network *network = (DM_Network*)dm->data;
1033: PetscInt offsetp;
1034: DMNetworkComponentHeader header;
1037: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
1038: PetscSectionGetOffset(network->DataSection,p,&offsetp);
1039: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetp);
1040: *index = header->index;
1041: return(0);
1042: }
1044: /*@
1045: DMNetworkGetGlobalEdgeIndex - Get the global numbering for the edge on the network
1047: Not Collective
1049: Input Parameters:
1050: + dm - DMNetwork object
1051: - p - edge point
1053: Output Parameters:
1054: . index - the global numbering for the edge
1056: Level: intermediate
1058: .seealso: DMNetworkGetGlobalVertexIndex()
1059: @*/
1060: PetscErrorCode DMNetworkGetGlobalEdgeIndex(DM dm,PetscInt p,PetscInt *index)
1061: {
1065: DMNetworkGetIndex(dm,p,index);
1066: return(0);
1067: }
1069: /*@
1070: DMNetworkGetGlobalVertexIndex - Get the global numbering for the vertex on the network
1072: Not Collective
1074: Input Parameters:
1075: + dm - DMNetwork object
1076: - p - vertex point
1078: Output Parameters:
1079: . index - the global numbering for the vertex
1081: Level: intermediate
1083: .seealso: DMNetworkGetGlobalEdgeIndex(), DMNetworkGetLocalVertexIndex()
1084: @*/
1085: PetscErrorCode DMNetworkGetGlobalVertexIndex(DM dm,PetscInt p,PetscInt *index)
1086: {
1090: DMNetworkGetIndex(dm,p,index);
1091: return(0);
1092: }
1094: /*@
1095: DMNetworkGetNumComponents - Get the number of components at a vertex/edge
1097: Not Collective
1099: Input Parameters:
1100: + dm - the DMNetwork object
1101: - p - vertex/edge point
1103: Output Parameters:
1104: . numcomponents - Number of components at the vertex/edge
1106: Level: beginner
1108: .seealso: DMNetworkRegisterComponent(), DMNetworkAddComponent()
1109: @*/
1110: PetscErrorCode DMNetworkGetNumComponents(DM dm,PetscInt p,PetscInt *numcomponents)
1111: {
1113: PetscInt offset;
1114: DM_Network *network = (DM_Network*)dm->data;
1117: PetscSectionGetOffset(network->DataSection,p,&offset);
1118: *numcomponents = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata;
1119: return(0);
1120: }
1122: /*@
1123: DMNetworkGetLocalVecOffset - Get the offset for accessing the variables associated with a component at the given vertex/edge from the local vector
1125: Not Collective
1127: Input Parameters:
1128: + dm - the DMNetwork object
1129: . p - the edge or vertex point
1130: - compnum - component number; use ALL_COMPONENTS if no specific component is requested
1132: Output Parameters:
1133: . offset - the local offset
1135: Notes:
1136: These offsets can be passed to MatSetValuesLocal() for matrices obtained with DMCreateMatrix().
1138: For vectors obtained with DMCreateLocalVector() the offsets can be used with VecSetValues().
1140: For vectors obtained with DMCreateLocalVector() and the array obtained with VecGetArray(vec,&array) you can access or set
1141: the vector values with array[offset].
1143: For vectors obtained with DMCreateGlobalVector() the offsets can be used with VecSetValuesLocal().
1145: Level: intermediate
1147: .seealso: DMGetLocalVector(), DMNetworkGetComponent(), DMNetworkGetGlobalVecOffset(), DMCreateGlobalVector(), VecGetArray(), VecSetValuesLocal(), MatSetValuesLocal()
1148: @*/
1149: PetscErrorCode DMNetworkGetLocalVecOffset(DM dm,PetscInt p,PetscInt compnum,PetscInt *offset)
1150: {
1152: DM_Network *network = (DM_Network*)dm->data;
1153: PetscInt offsetp,offsetd;
1154: DMNetworkComponentHeader header;
1157: PetscSectionGetOffset(network->plex->localSection,p,&offsetp);
1158: if (compnum == ALL_COMPONENTS) {
1159: *offset = offsetp;
1160: return(0);
1161: }
1163: PetscSectionGetOffset(network->DataSection,p,&offsetd);
1164: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetd);
1165: *offset = offsetp + header->offsetvarrel[compnum];
1166: return(0);
1167: }
1169: /*@
1170: DMNetworkGetGlobalVecOffset - Get the global offset for accessing the variables associated with a component for the given vertex/edge from the global vector
1172: Not Collective
1174: Input Parameters:
1175: + dm - the DMNetwork object
1176: . p - the edge or vertex point
1177: - compnum - component number; use ALL_COMPONENTS if no specific component is requested
1179: Output Parameters:
1180: . offsetg - the global offset
1182: Notes:
1183: These offsets can be passed to MatSetValues() for matrices obtained with DMCreateMatrix().
1185: For vectors obtained with DMCreateGlobalVector() the offsets can be used with VecSetValues().
1187: For vectors obtained with DMCreateGlobalVector() and the array obtained with VecGetArray(vec,&array) you can access or set
1188: the vector values with array[offset - rstart] where restart is obtained with VecGetOwnershipRange(v,&rstart,NULL);
1190: Level: intermediate
1192: .seealso: DMNetworkGetLocalVecOffset(), DMGetGlobalVector(), DMNetworkGetComponent(), DMCreateGlobalVector(), VecGetArray(), VecSetValues(), MatSetValues()
1193: @*/
1194: PetscErrorCode DMNetworkGetGlobalVecOffset(DM dm,PetscInt p,PetscInt compnum,PetscInt *offsetg)
1195: {
1197: DM_Network *network = (DM_Network*)dm->data;
1198: PetscInt offsetp,offsetd;
1199: DMNetworkComponentHeader header;
1202: PetscSectionGetOffset(network->plex->globalSection,p,&offsetp);
1203: if (offsetp < 0) offsetp = -(offsetp + 1); /* Convert to actual global offset for ghost vertex */
1205: if (compnum == ALL_COMPONENTS) {
1206: *offsetg = offsetp;
1207: return(0);
1208: }
1209: PetscSectionGetOffset(network->DataSection,p,&offsetd);
1210: header = (DMNetworkComponentHeader)(network->componentdataarray+offsetd);
1211: *offsetg = offsetp + header->offsetvarrel[compnum];
1212: return(0);
1213: }
1215: /*@
1216: DMNetworkGetEdgeOffset - Get the offset for accessing the variables associated with the given edge from the local subvector
1218: Not Collective
1220: Input Parameters:
1221: + dm - the DMNetwork object
1222: - p - the edge point
1224: Output Parameters:
1225: . offset - the offset
1227: Level: intermediate
1229: .seealso: DMNetworkGetLocalVecOffset(), DMGetLocalVector()
1230: @*/
1231: PetscErrorCode DMNetworkGetEdgeOffset(DM dm,PetscInt p,PetscInt *offset)
1232: {
1234: DM_Network *network = (DM_Network*)dm->data;
1237: PetscSectionGetOffset(network->edge.DofSection,p,offset);
1238: return(0);
1239: }
1241: /*@
1242: DMNetworkGetVertexOffset - Get the offset for accessing the variables associated with the given vertex from the local subvector
1244: Not Collective
1246: Input Parameters:
1247: + dm - the DMNetwork object
1248: - p - the vertex point
1250: Output Parameters:
1251: . offset - the offset
1253: Level: intermediate
1255: .seealso: DMNetworkGetEdgeOffset(), DMGetLocalVector()
1256: @*/
1257: PetscErrorCode DMNetworkGetVertexOffset(DM dm,PetscInt p,PetscInt *offset)
1258: {
1260: DM_Network *network = (DM_Network*)dm->data;
1263: p -= network->vStart;
1264: PetscSectionGetOffset(network->vertex.DofSection,p,offset);
1265: return(0);
1266: }
1268: /*@
1269: DMNetworkAddComponent - Adds a network component and number of variables at the given point (vertex/edge)
1271: Not Collective
1273: Input Parameters:
1274: + dm - the DMNetwork
1275: . p - the vertex/edge point
1276: . componentkey - component key returned while registering the component; ignored if compvalue=NULL
1277: . compvalue - pointer to the data structure for the component, or NULL if not required.
1278: - nvar - number of variables for the component at the vertex/edge point
1280: Level: beginner
1282: .seealso: DMNetworkGetComponent()
1283: @*/
1284: PetscErrorCode DMNetworkAddComponent(DM dm,PetscInt p,PetscInt componentkey,void* compvalue,PetscInt nvar)
1285: {
1286: PetscErrorCode ierr;
1287: DM_Network *network = (DM_Network*)dm->data;
1288: DMNetworkComponent *component = &network->component[componentkey];
1289: DMNetworkComponentHeader header;
1290: DMNetworkComponentValue cvalue;
1291: PetscBool sharedv=PETSC_FALSE;
1292: PetscInt compnum;
1293: PetscInt *compsize,*compkey,*compoffset,*compnvar,*compoffsetvarrel;
1294: void* *compdata;
1297: PetscSectionAddDof(network->DofSection,p,nvar);
1298: if (!compvalue) return(0);
1300: DMNetworkIsSharedVertex(dm,p,&sharedv);
1301: if (sharedv) {
1302: PetscBool ghost;
1303: DMNetworkIsGhostVertex(dm,p,&ghost);
1304: if (ghost) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Adding a component at a leaf(ghost) shared vertex is not supported");
1305: }
1307: header = &network->header[p];
1308: cvalue = &network->cvalue[p];
1309: if (header->ndata == header->maxcomps) {
1310: PetscInt additional_size;
1312: /* Reached limit so resize header component arrays */
1313: header->maxcomps += 2;
1315: /* Allocate arrays for component information and value */
1316: PetscCalloc5(header->maxcomps,&compsize,header->maxcomps,&compkey,header->maxcomps,&compoffset,header->maxcomps,&compnvar,header->maxcomps,&compoffsetvarrel);
1317: PetscMalloc1(header->maxcomps,&compdata);
1319: /* Recalculate header size */
1320: header->hsize = sizeof(struct _p_DMNetworkComponentHeader) + 5*header->maxcomps*sizeof(PetscInt);
1322: header->hsize /= sizeof(DMNetworkComponentGenericDataType);
1324: /* Copy over component info */
1325: PetscMemcpy(compsize,header->size,header->ndata*sizeof(PetscInt));
1326: PetscMemcpy(compkey,header->key,header->ndata*sizeof(PetscInt));
1327: PetscMemcpy(compoffset,header->offset,header->ndata*sizeof(PetscInt));
1328: PetscMemcpy(compnvar,header->nvar,header->ndata*sizeof(PetscInt));
1329: PetscMemcpy(compoffsetvarrel,header->offsetvarrel,header->ndata*sizeof(PetscInt));
1331: /* Copy over component data pointers */
1332: PetscMemcpy(compdata,cvalue->data,header->ndata*sizeof(void*));
1334: /* Free old arrays */
1335: PetscFree5(header->size,header->key,header->offset,header->nvar,header->offsetvarrel);
1336: PetscFree(cvalue->data);
1338: /* Update pointers */
1339: header->size = compsize;
1340: header->key = compkey;
1341: header->offset = compoffset;
1342: header->nvar = compnvar;
1343: header->offsetvarrel = compoffsetvarrel;
1345: cvalue->data = compdata;
1347: /* Update DataSection Dofs */
1348: /* The dofs for datasection point p equals sizeof the header (i.e. header->hsize) + sizes of the components added at point p. With the resizing of the header, we need to update the dofs for point p. Hence, we add the extra size added for the header */
1349: additional_size = (5*(header->maxcomps - header->ndata)*sizeof(PetscInt))/sizeof(DMNetworkComponentGenericDataType);
1350: PetscSectionAddDof(network->DataSection,p,additional_size);
1351: }
1352: header = &network->header[p];
1353: cvalue = &network->cvalue[p];
1355: compnum = header->ndata;
1357: header->size[compnum] = component->size;
1358: PetscSectionAddDof(network->DataSection,p,component->size);
1359: header->key[compnum] = componentkey;
1360: if (compnum != 0) header->offset[compnum] = header->offset[compnum-1] + header->size[compnum-1];
1361: cvalue->data[compnum] = (void*)compvalue;
1363: /* variables */
1364: header->nvar[compnum] += nvar;
1365: if (compnum != 0) header->offsetvarrel[compnum] = header->offsetvarrel[compnum-1] + header->nvar[compnum-1];
1367: header->ndata++;
1368: return(0);
1369: }
1371: /*@
1372: DMNetworkGetComponent - Gets the component key, the component data, and the number of variables at a given network point
1374: Not Collective
1376: Input Parameters:
1377: + dm - the DMNetwork object
1378: . p - vertex/edge point
1379: - compnum - component number; use ALL_COMPONENTS if sum up all the components
1381: Output Parameters:
1382: + compkey - the key obtained when registering the component (use NULL if not required)
1383: . component - the component data (use NULL if not required)
1384: - nvar - number of variables (use NULL if not required)
1386: Level: beginner
1388: .seealso: DMNetworkAddComponent(), DMNetworkGetNumComponents()
1389: @*/
1390: PetscErrorCode DMNetworkGetComponent(DM dm,PetscInt p,PetscInt compnum,PetscInt *compkey,void **component,PetscInt *nvar)
1391: {
1393: DM_Network *network = (DM_Network*)dm->data;
1394: PetscInt offset = 0;
1395: DMNetworkComponentHeader header;
1398: if (compnum == ALL_COMPONENTS) {
1399: PetscSectionGetDof(network->DofSection,p,nvar);
1400: return(0);
1401: }
1403: PetscSectionGetOffset(network->DataSection,p,&offset);
1404: header = (DMNetworkComponentHeader)(network->componentdataarray+offset);
1406: if (compnum >= 0) {
1407: if (compkey) *compkey = header->key[compnum];
1408: if (component) {
1409: offset += header->hsize+header->offset[compnum];
1410: *component = network->componentdataarray+offset;
1411: }
1412: }
1414: if (nvar) *nvar = header->nvar[compnum];
1416: return(0);
1417: }
1419: /*
1420: Sets up the array that holds the data for all components and its associated section.
1421: It copies the data for all components in a contiguous array called componentdataarray. The component data is stored pointwise with an additional header (metadata) stored for each point. The header has metadata information such as number of components at each point, number of variables for each component, offsets for the components data, etc.
1422: */
1423: PetscErrorCode DMNetworkComponentSetUp(DM dm)
1424: {
1425: PetscErrorCode ierr;
1426: DM_Network *network = (DM_Network*)dm->data;
1427: PetscInt arr_size,p,offset,offsetp,ncomp,i,*headerarr;
1428: DMNetworkComponentHeader header;
1429: DMNetworkComponentValue cvalue;
1430: DMNetworkComponentHeader headerinfo;
1431: DMNetworkComponentGenericDataType *componentdataarray;
1434: PetscSectionSetUp(network->DataSection);
1435: PetscSectionGetStorageSize(network->DataSection,&arr_size);
1436: /* arr_size+1 fixes pipeline test of opensolaris-misc for src/dm/tests/ex10.c -- Do not know why */
1437: PetscCalloc1(arr_size+1,&network->componentdataarray);
1438: componentdataarray = network->componentdataarray;
1439: for (p = network->pStart; p < network->pEnd; p++) {
1440: PetscSectionGetOffset(network->DataSection,p,&offsetp);
1441: /* Copy header */
1442: header = &network->header[p];
1443: headerinfo = (DMNetworkComponentHeader)(componentdataarray+offsetp);
1444: PetscMemcpy(headerinfo,header,sizeof(struct _p_DMNetworkComponentHeader));
1445: headerarr = (PetscInt*)(headerinfo+1);
1446: PetscMemcpy(headerarr,header->size,header->maxcomps*sizeof(PetscInt));
1447: headerarr += header->maxcomps;
1448: PetscMemcpy(headerarr,header->key,header->maxcomps*sizeof(PetscInt));
1449: headerarr += header->maxcomps;
1450: PetscMemcpy(headerarr,header->offset,header->maxcomps*sizeof(PetscInt));
1451: headerarr += header->maxcomps;
1452: PetscMemcpy(headerarr,header->nvar,header->maxcomps*sizeof(PetscInt));
1453: headerarr += header->maxcomps;
1454: PetscMemcpy(headerarr,header->offsetvarrel,header->maxcomps*sizeof(PetscInt));
1456: /* Copy data */
1457: cvalue = &network->cvalue[p];
1458: ncomp = header->ndata;
1460: for (i = 0; i < ncomp; i++) {
1461: offset = offsetp + header->hsize + header->offset[i];
1462: PetscMemcpy(componentdataarray+offset,cvalue->data[i],header->size[i]*sizeof(DMNetworkComponentGenericDataType));
1463: }
1464: }
1465: return(0);
1466: }
1468: /* Sets up the section for dofs. This routine is called during DMSetUp() */
1469: static PetscErrorCode DMNetworkVariablesSetUp(DM dm)
1470: {
1472: DM_Network *network = (DM_Network*)dm->data;
1473: MPI_Comm comm;
1474: PetscMPIInt size;
1477: PetscObjectGetComm((PetscObject)dm,&comm);
1478: MPI_Comm_size(comm,&size);
1480: if (size > 1) { /* Sync nvar at shared vertices for all processes */
1481: PetscSF sf = network->plex->sf;
1482: PetscInt *local_nvar, *remote_nvar,nroots,nleaves,p=-1,i,nsv;
1483: const PetscInt *svtx;
1484: PetscBool ghost;
1485: /*
1486: PetscMPIInt rank;
1487: const PetscInt *ilocal;
1488: const PetscSFNode *iremote;
1489: MPI_Comm_rank(comm,&rank);
1490: PetscSFGetGraph(sf,&nroots,&nleaves,&ilocal,&iremote);
1491: */
1492: PetscSFGetGraph(sf,&nroots,&nleaves,NULL,NULL);
1493: PetscCalloc2(nroots,&local_nvar,nroots,&remote_nvar);
1495: /* Leaves copy user's nvar to local_nvar */
1496: DMNetworkGetSharedVertices(dm,&nsv,&svtx);
1497: for (i=0; i<nsv; i++) {
1498: p = svtx[i];
1499: DMNetworkIsGhostVertex(dm,p,&ghost);
1500: if (!ghost) continue;
1501: PetscSectionGetDof(network->DofSection,p,&local_nvar[p]);
1502: /* printf("[%d] Before SFReduce: leaf local_nvar[%d] = %d\n",rank,p,local_nvar[p]); */
1503: }
1505: /* Leaves add local_nvar to root remote_nvar */
1506: PetscSFReduceBegin(sf, MPIU_INT, local_nvar, remote_nvar, MPI_SUM);
1507: PetscSFReduceEnd(sf, MPIU_INT, local_nvar, remote_nvar, MPI_SUM);
1508: /* printf("[%d] remote_nvar[%d] = %d\n",rank,p,remote_nvar[p]); */
1510: /* Update roots' local_nvar */
1511: for (i=0; i<nsv; i++) {
1512: p = svtx[i];
1513: DMNetworkIsGhostVertex(dm,p,&ghost);
1514: if (ghost) continue;
1516: PetscSectionAddDof(network->DofSection,p,remote_nvar[p]);
1517: PetscSectionGetDof(network->DofSection,p,&local_nvar[p]);
1518: /* printf("[%d] After SFReduce: root local_nvar[%d] = %d\n",rank,p,local_nvar[p]); */
1519: }
1521: /* Roots Bcast nvar to leaves */
1522: PetscSFBcastBegin(sf, MPIU_INT, local_nvar, remote_nvar,MPI_REPLACE);
1523: PetscSFBcastEnd(sf, MPIU_INT, local_nvar, remote_nvar,MPI_REPLACE);
1525: /* Leaves reset receved/remote nvar to dm */
1526: for (i=0; i<nsv; i++) {
1527: DMNetworkIsGhostVertex(dm,p,&ghost);
1528: if (!ghost) continue;
1529: p = svtx[i];
1530: /* printf("[%d] leaf reset nvar %d at p= %d \n",rank,remote_nvar[p],p); */
1531: /* DMNetworkSetNumVariables(dm,p,remote_nvar[p]) */
1532: PetscSectionSetDof(network->DofSection,p,remote_nvar[p]);
1533: }
1535: PetscFree2(local_nvar,remote_nvar);
1536: }
1538: PetscSectionSetUp(network->DofSection);
1539: return(0);
1540: }
1542: /* Get a subsection from a range of points */
1543: static PetscErrorCode DMNetworkGetSubSection_private(PetscSection main,PetscInt pstart,PetscInt pend,PetscSection *subsection)
1544: {
1546: PetscInt i, nvar;
1549: PetscSectionCreate(PetscObjectComm((PetscObject)main), subsection);
1550: PetscSectionSetChart(*subsection, 0, pend - pstart);
1551: for (i = pstart; i < pend; i++) {
1552: PetscSectionGetDof(main,i,&nvar);
1553: PetscSectionSetDof(*subsection, i - pstart, nvar);
1554: }
1556: PetscSectionSetUp(*subsection);
1557: return(0);
1558: }
1560: /* Create a submap of points with a GlobalToLocal structure */
1561: static PetscErrorCode DMNetworkSetSubMap_private(PetscInt pstart, PetscInt pend, ISLocalToGlobalMapping *map)
1562: {
1564: PetscInt i, *subpoints;
1567: /* Create index sets to map from "points" to "subpoints" */
1568: PetscMalloc1(pend - pstart, &subpoints);
1569: for (i = pstart; i < pend; i++) {
1570: subpoints[i - pstart] = i;
1571: }
1572: ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,pend-pstart,subpoints,PETSC_COPY_VALUES,map);
1573: PetscFree(subpoints);
1574: return(0);
1575: }
1577: /*@
1578: DMNetworkAssembleGraphStructures - Assembles vertex and edge data structures. Must be called after DMNetworkDistribute
1580: Collective on dm
1582: Input Parameters:
1583: . dm - the DMNetworkObject
1585: Note: the routine will create alternative orderings for the vertices and edges. Assume global network points are:
1587: points = [0 1 2 3 4 5 6]
1589: 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]).
1591: With this new ordering a local PetscSection, global PetscSection and PetscSF will be created specific to the subset.
1593: Level: intermediate
1595: @*/
1596: PetscErrorCode DMNetworkAssembleGraphStructures(DM dm)
1597: {
1599: MPI_Comm comm;
1600: PetscMPIInt size;
1601: DM_Network *network = (DM_Network*)dm->data;
1604: PetscObjectGetComm((PetscObject)dm,&comm);
1605: MPI_Comm_size(comm, &size);
1607: /* Create maps for vertices and edges */
1608: DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);
1609: DMNetworkSetSubMap_private(network->eStart,network->eEnd,&network->edge.mapping);
1611: /* Create local sub-sections */
1612: DMNetworkGetSubSection_private(network->DofSection,network->vStart,network->vEnd,&network->vertex.DofSection);
1613: DMNetworkGetSubSection_private(network->DofSection,network->eStart,network->eEnd,&network->edge.DofSection);
1615: if (size > 1) {
1616: PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);
1618: PetscSectionCreateGlobalSection(network->vertex.DofSection, network->vertex.sf, PETSC_FALSE, PETSC_FALSE, &network->vertex.GlobalDofSection);
1619: PetscSFGetSubSF(network->plex->sf, network->edge.mapping, &network->edge.sf);
1620: PetscSectionCreateGlobalSection(network->edge.DofSection, network->edge.sf, PETSC_FALSE, PETSC_FALSE, &network->edge.GlobalDofSection);
1621: } else {
1622: /* create structures for vertex */
1623: PetscSectionClone(network->vertex.DofSection,&network->vertex.GlobalDofSection);
1624: /* create structures for edge */
1625: PetscSectionClone(network->edge.DofSection,&network->edge.GlobalDofSection);
1626: }
1628: /* Add viewers */
1629: PetscObjectSetName((PetscObject)network->edge.GlobalDofSection,"Global edge dof section");
1630: PetscObjectSetName((PetscObject)network->vertex.GlobalDofSection,"Global vertex dof section");
1631: PetscSectionViewFromOptions(network->edge.GlobalDofSection, NULL, "-edge_global_section_view");
1632: PetscSectionViewFromOptions(network->vertex.GlobalDofSection, NULL, "-vertex_global_section_view");
1633: return(0);
1634: }
1636: /*
1637: Add all subnetid for the input vertex v in this process to the btable
1638: vertex_subnetid = supportingedge_subnetid
1639: */
1640: PETSC_STATIC_INLINE PetscErrorCode SetSubnetIdLookupBT(DM dm,PetscInt v,PetscInt Nsubnet,PetscBT btable)
1641: {
1643: PetscInt e,nedges,offset;
1644: const PetscInt *edges;
1645: DM_Network *newDMnetwork = (DM_Network*)dm->data;
1646: DMNetworkComponentHeader header;
1649: DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
1650: PetscBTMemzero(Nsubnet,btable);
1651: for (e=0; e<nedges; e++) {
1652: PetscSectionGetOffset(newDMnetwork->DataSection,edges[e],&offset);
1653: header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1654: PetscBTSet(btable,header->subnetid);
1655: }
1656: return(0);
1657: }
1659: /*@
1660: DMNetworkDistribute - Distributes the network and moves associated component data
1662: Collective
1664: Input Parameters:
1665: + DM - the DMNetwork object
1666: - overlap - the overlap of partitions, 0 is the default
1668: Options Database Key:
1669: + -dmnetwork_view - Calls DMView() at the conclusion of DMSetUp()
1670: - -dmnetwork_view_distributed - Calls DMView() at the conclusion of DMNetworkDistribute()
1672: Notes:
1673: Distributes the network with <overlap>-overlapping partitioning of the edges.
1675: Level: intermediate
1677: .seealso: DMNetworkCreate()
1678: @*/
1679: PetscErrorCode DMNetworkDistribute(DM *dm,PetscInt overlap)
1680: {
1681: MPI_Comm comm;
1683: PetscMPIInt size;
1684: DM_Network *oldDMnetwork = (DM_Network*)((*dm)->data);
1685: DM_Network *newDMnetwork;
1686: PetscSF pointsf=NULL;
1687: DM newDM;
1688: PetscInt j,e,v,offset,*subnetvtx,*subnetedge,Nsubnet,gidx,svtx_idx,nv;
1689: PetscInt to_net,from_net,*svto;
1690: PetscBT btable;
1691: PetscPartitioner part;
1692: DMNetworkComponentHeader header;
1695: PetscObjectGetComm((PetscObject)*dm,&comm);
1696: MPI_Comm_size(comm, &size);
1697: if (size == 1) return(0);
1699: if (overlap) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"overlap %D != 0 is not supported yet",overlap);
1701: /* This routine moves the component data to the appropriate processors. It makes use of the DataSection and the componentdataarray to move the component data to appropriate processors and returns a new DataSection and new componentdataarray. */
1702: DMNetworkCreate(PetscObjectComm((PetscObject)*dm),&newDM);
1703: newDMnetwork = (DM_Network*)newDM->data;
1704: newDMnetwork->max_comps_registered = oldDMnetwork->max_comps_registered;
1705: PetscMalloc1(newDMnetwork->max_comps_registered,&newDMnetwork->component);
1707: /* Enable runtime options for petscpartitioner */
1708: DMPlexGetPartitioner(oldDMnetwork->plex,&part);
1709: PetscPartitionerSetFromOptions(part);
1711: /* Distribute plex dm */
1712: DMPlexDistribute(oldDMnetwork->plex,overlap,&pointsf,&newDMnetwork->plex);
1714: /* Distribute dof section */
1715: PetscSectionCreate(comm,&newDMnetwork->DofSection);
1716: PetscSFDistributeSection(pointsf,oldDMnetwork->DofSection,NULL,newDMnetwork->DofSection);
1718: /* Distribute data and associated section */
1719: PetscSectionCreate(comm,&newDMnetwork->DataSection);
1720: DMPlexDistributeData(newDMnetwork->plex,pointsf,oldDMnetwork->DataSection,MPIU_INT,(void*)oldDMnetwork->componentdataarray,newDMnetwork->DataSection,(void**)&newDMnetwork->componentdataarray);
1722: PetscSectionGetChart(newDMnetwork->DataSection,&newDMnetwork->pStart,&newDMnetwork->pEnd);
1723: DMPlexGetHeightStratum(newDMnetwork->plex,0, &newDMnetwork->eStart,&newDMnetwork->eEnd);
1724: DMPlexGetHeightStratum(newDMnetwork->plex,1,&newDMnetwork->vStart,&newDMnetwork->vEnd);
1725: newDMnetwork->nEdges = newDMnetwork->eEnd - newDMnetwork->eStart;
1726: newDMnetwork->nVertices = newDMnetwork->vEnd - newDMnetwork->vStart;
1727: newDMnetwork->NVertices = oldDMnetwork->NVertices;
1728: newDMnetwork->NEdges = oldDMnetwork->NEdges;
1729: newDMnetwork->svtable = oldDMnetwork->svtable; /* global table! */
1730: oldDMnetwork->svtable = NULL;
1732: /* Set Dof section as the section for dm */
1733: DMSetLocalSection(newDMnetwork->plex,newDMnetwork->DofSection);
1734: DMGetGlobalSection(newDMnetwork->plex,&newDMnetwork->GlobalDofSection);
1736: /* Setup subnetwork info in the newDM */
1737: newDMnetwork->Nsubnet = oldDMnetwork->Nsubnet;
1738: newDMnetwork->Nsvtx = oldDMnetwork->Nsvtx;
1739: oldDMnetwork->Nsvtx = 0;
1740: newDMnetwork->svtx = oldDMnetwork->svtx; /* global vertices! */
1741: oldDMnetwork->svtx = NULL;
1742: PetscCalloc1(newDMnetwork->Nsubnet,&newDMnetwork->subnet);
1744: /* Copy over the global number of vertices and edges in each subnetwork.
1745: Note: these are calculated in DMNetworkLayoutSetUp()
1746: */
1747: Nsubnet = newDMnetwork->Nsubnet;
1748: for (j = 0; j < Nsubnet; j++) {
1749: newDMnetwork->subnet[j].Nvtx = oldDMnetwork->subnet[j].Nvtx;
1750: newDMnetwork->subnet[j].Nedge = oldDMnetwork->subnet[j].Nedge;
1751: }
1753: /* Count local nedges for subnetworks */
1754: for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++) {
1755: PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);
1756: header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1757: /* Update pointers */
1758: header->size = (PetscInt*)(header + 1);
1759: header->key = header->size + header->maxcomps;
1760: header->offset = header->key + header->maxcomps;
1761: header->nvar = header->offset + header->maxcomps;
1762: header->offsetvarrel = header->nvar + header->maxcomps;
1764: newDMnetwork->subnet[header->subnetid].nedge++;
1765: }
1767: /* Count local nvtx for subnetworks */
1768: PetscBTCreate(Nsubnet,&btable);
1769: for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++) {
1770: PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);
1771: header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1772: /* Update pointers */
1773: header->size = (PetscInt*)(header + 1);
1774: header->key = header->size + header->maxcomps;
1775: header->offset = header->key + header->maxcomps;
1776: header->nvar = header->offset + header->maxcomps;
1777: header->offsetvarrel = header->nvar + header->maxcomps;
1779: /* shared vertices: use gidx=header->index to check if v is a shared vertex */
1780: gidx = header->index;
1781: PetscTableFind(newDMnetwork->svtable,gidx+1,&svtx_idx);
1782: svtx_idx--;
1784: if (svtx_idx < 0) { /* not a shared vertex */
1785: newDMnetwork->subnet[header->subnetid].nvtx++;
1786: } else { /* a shared vertex belongs to more than one subnetworks, it is being counted by multiple subnets */
1787: SetSubnetIdLookupBT(newDM,v,Nsubnet,btable);
1789: from_net = newDMnetwork->svtx[svtx_idx].sv[0];
1790: if (PetscBTLookup(btable,from_net)) newDMnetwork->subnet[from_net].nvtx++; /* sv is on from_net */
1792: for (j=1; j<newDMnetwork->svtx[svtx_idx].n; j++) {
1793: svto = newDMnetwork->svtx[svtx_idx].sv + 2*j;
1794: to_net = svto[0];
1795: if (PetscBTLookup(btable,to_net)) newDMnetwork->subnet[to_net].nvtx++; /* sv is on to_net */
1796: }
1797: }
1798: }
1800: /* Get total local nvtx for subnetworks */
1801: nv = 0;
1802: for (j=0; j<Nsubnet; j++) nv += newDMnetwork->subnet[j].nvtx;
1803: nv += newDMnetwork->Nsvtx;
1805: /* Now create the vertices and edge arrays for the subnetworks */
1806: PetscCalloc2(newDMnetwork->nEdges,&subnetedge,nv,&subnetvtx); /* Maps local vertex to local subnetwork's vertex */
1807: newDMnetwork->subnetedge = subnetedge;
1808: newDMnetwork->subnetvtx = subnetvtx;
1809: for (j=0; j < newDMnetwork->Nsubnet; j++) {
1810: newDMnetwork->subnet[j].edges = subnetedge;
1811: subnetedge += newDMnetwork->subnet[j].nedge;
1813: newDMnetwork->subnet[j].vertices = subnetvtx;
1814: subnetvtx += newDMnetwork->subnet[j].nvtx;
1816: /* Temporarily setting nvtx and nedge to 0 so we can use them as counters in the below for loop. These get updated when the vertices and edges are added. */
1817: newDMnetwork->subnet[j].nvtx = newDMnetwork->subnet[j].nedge = 0;
1818: }
1819: newDMnetwork->svertices = subnetvtx;
1821: /* Set the edges and vertices in each subnetwork */
1822: for (e = newDMnetwork->eStart; e < newDMnetwork->eEnd; e++) {
1823: PetscSectionGetOffset(newDMnetwork->DataSection,e,&offset);
1824: header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1825: newDMnetwork->subnet[header->subnetid].edges[newDMnetwork->subnet[header->subnetid].nedge++] = e;
1826: }
1828: nv = 0;
1829: for (v = newDMnetwork->vStart; v < newDMnetwork->vEnd; v++) {
1830: PetscSectionGetOffset(newDMnetwork->DataSection,v,&offset);
1831: header = (DMNetworkComponentHeader)(newDMnetwork->componentdataarray+offset);
1833: /* coupling vertices: use gidx = header->index to check if v is a coupling vertex */
1834: PetscTableFind(newDMnetwork->svtable,header->index+1,&svtx_idx);
1835: svtx_idx--;
1836: if (svtx_idx < 0) {
1837: newDMnetwork->subnet[header->subnetid].vertices[newDMnetwork->subnet[header->subnetid].nvtx++] = v;
1838: } else { /* a shared vertex */
1839: newDMnetwork->svertices[nv++] = v;
1841: /* add all subnetid for this shared vertex in this process to btable */
1842: SetSubnetIdLookupBT(newDM,v,Nsubnet,btable);
1844: from_net = newDMnetwork->svtx[svtx_idx].sv[0];
1845: if (PetscBTLookup(btable,from_net))
1846: newDMnetwork->subnet[from_net].vertices[newDMnetwork->subnet[from_net].nvtx++] = v;
1848: for (j=1; j<newDMnetwork->svtx[svtx_idx].n; j++) {
1849: svto = newDMnetwork->svtx[svtx_idx].sv + 2*j;
1850: to_net = svto[0];
1851: if (PetscBTLookup(btable,to_net))
1852: newDMnetwork->subnet[to_net].vertices[newDMnetwork->subnet[to_net].nvtx++] = v;
1853: }
1854: }
1855: }
1856: newDMnetwork->nsvtx = nv; /* num of local shared vertices */
1858: newDM->setupcalled = (*dm)->setupcalled;
1859: newDMnetwork->distributecalled = PETSC_TRUE;
1861: /* Free spaces */
1862: PetscSFDestroy(&pointsf);
1863: DMDestroy(dm);
1864: PetscBTDestroy(&btable);
1866: /* View distributed dmnetwork */
1867: DMViewFromOptions(newDM,NULL,"-dmnetwork_view_distributed");
1869: *dm = newDM;
1870: return(0);
1871: }
1873: /*@C
1874: PetscSFGetSubSF - Returns an SF for a specific subset of points. Leaves are re-numbered to reflect the new ordering
1876: Collective
1878: Input Parameters:
1879: + mainSF - the original SF structure
1880: - map - a ISLocalToGlobal mapping that contains the subset of points
1882: Output Parameter:
1883: . subSF - a subset of the mainSF for the desired subset.
1885: Level: intermediate
1886: @*/
1887: PetscErrorCode PetscSFGetSubSF(PetscSF mainsf,ISLocalToGlobalMapping map,PetscSF *subSF)
1888: {
1889: PetscErrorCode ierr;
1890: PetscInt nroots, nleaves, *ilocal_sub;
1891: PetscInt i, *ilocal_map, nroots_sub, nleaves_sub = 0;
1892: PetscInt *local_points, *remote_points;
1893: PetscSFNode *iremote_sub;
1894: const PetscInt *ilocal;
1895: const PetscSFNode *iremote;
1898: PetscSFGetGraph(mainsf,&nroots,&nleaves,&ilocal,&iremote);
1900: /* Look for leaves that pertain to the subset of points. Get the local ordering */
1901: PetscMalloc1(nleaves,&ilocal_map);
1902: ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nleaves,ilocal,NULL,ilocal_map);
1903: for (i = 0; i < nleaves; i++) {
1904: if (ilocal_map[i] != -1) nleaves_sub += 1;
1905: }
1906: /* Re-number ilocal with subset numbering. Need information from roots */
1907: PetscMalloc2(nroots,&local_points,nroots,&remote_points);
1908: for (i = 0; i < nroots; i++) local_points[i] = i;
1909: ISGlobalToLocalMappingApply(map,IS_GTOLM_MASK,nroots,local_points,NULL,local_points);
1910: PetscSFBcastBegin(mainsf, MPIU_INT, local_points, remote_points,MPI_REPLACE);
1911: PetscSFBcastEnd(mainsf, MPIU_INT, local_points, remote_points,MPI_REPLACE);
1912: /* Fill up graph using local (that is, local to the subset) numbering. */
1913: PetscMalloc1(nleaves_sub,&ilocal_sub);
1914: PetscMalloc1(nleaves_sub,&iremote_sub);
1915: nleaves_sub = 0;
1916: for (i = 0; i < nleaves; i++) {
1917: if (ilocal_map[i] != -1) {
1918: ilocal_sub[nleaves_sub] = ilocal_map[i];
1919: iremote_sub[nleaves_sub].rank = iremote[i].rank;
1920: iremote_sub[nleaves_sub].index = remote_points[ilocal[i]];
1921: nleaves_sub += 1;
1922: }
1923: }
1924: PetscFree2(local_points,remote_points);
1925: ISLocalToGlobalMappingGetSize(map,&nroots_sub);
1927: /* Create new subSF */
1928: PetscSFCreate(PETSC_COMM_WORLD,subSF);
1929: PetscSFSetFromOptions(*subSF);
1930: PetscSFSetGraph(*subSF,nroots_sub,nleaves_sub,ilocal_sub,PETSC_OWN_POINTER,iremote_sub,PETSC_COPY_VALUES);
1931: PetscFree(ilocal_map);
1932: PetscFree(iremote_sub);
1933: return(0);
1934: }
1936: /*@C
1937: DMNetworkGetSupportingEdges - Return the supporting edges for this vertex point
1939: Not Collective
1941: Input Parameters:
1942: + dm - the DMNetwork object
1943: - p - the vertex point
1945: Output Parameters:
1946: + nedges - number of edges connected to this vertex point
1947: - edges - list of edge points
1949: Level: beginner
1951: Fortran Notes:
1952: Since it returns an array, this routine is only available in Fortran 90, and you must
1953: include petsc.h90 in your code.
1955: .seealso: DMNetworkCreate(), DMNetworkGetConnectedVertices()
1956: @*/
1957: PetscErrorCode DMNetworkGetSupportingEdges(DM dm,PetscInt vertex,PetscInt *nedges,const PetscInt *edges[])
1958: {
1960: DM_Network *network = (DM_Network*)dm->data;
1963: DMPlexGetSupportSize(network->plex,vertex,nedges);
1964: if (edges) {DMPlexGetSupport(network->plex,vertex,edges);}
1965: return(0);
1966: }
1968: /*@C
1969: DMNetworkGetConnectedVertices - Return the connected vertices for this edge point
1971: Not Collective
1973: Input Parameters:
1974: + dm - the DMNetwork object
1975: - p - the edge point
1977: Output Parameters:
1978: . vertices - vertices connected to this edge
1980: Level: beginner
1982: Fortran Notes:
1983: Since it returns an array, this routine is only available in Fortran 90, and you must
1984: include petsc.h90 in your code.
1986: .seealso: DMNetworkCreate(), DMNetworkGetSupportingEdges()
1987: @*/
1988: PetscErrorCode DMNetworkGetConnectedVertices(DM dm,PetscInt edge,const PetscInt *vertices[])
1989: {
1991: DM_Network *network = (DM_Network*)dm->data;
1994: DMPlexGetCone(network->plex,edge,vertices);
1995: return(0);
1996: }
1998: /*@
1999: DMNetworkIsSharedVertex - Returns TRUE if the vertex is shared by subnetworks
2001: Not Collective
2003: Input Parameters:
2004: + dm - the DMNetwork object
2005: - p - the vertex point
2007: Output Parameter:
2008: . flag - TRUE if the vertex is shared by subnetworks
2010: Level: beginner
2012: .seealso: DMNetworkAddSharedVertices(), DMNetworkIsGhostVertex()
2013: @*/
2014: PetscErrorCode DMNetworkIsSharedVertex(DM dm,PetscInt p,PetscBool *flag)
2015: {
2017: PetscInt i;
2020: *flag = PETSC_FALSE;
2022: if (dm->setupcalled) { /* DMNetworkGetGlobalVertexIndex() requires DMSetUp() be called */
2023: DM_Network *network = (DM_Network*)dm->data;
2024: PetscInt gidx;
2025: DMNetworkGetGlobalVertexIndex(dm,p,&gidx);
2026: PetscTableFind(network->svtable,gidx+1,&i);
2027: if (i) *flag = PETSC_TRUE;
2028: } else { /* would be removed? */
2029: PetscInt nv;
2030: const PetscInt *vtx;
2031: DMNetworkGetSharedVertices(dm,&nv,&vtx);
2032: for (i=0; i<nv; i++) {
2033: if (p == vtx[i]) {
2034: *flag = PETSC_TRUE;
2035: break;
2036: }
2037: }
2038: }
2039: return(0);
2040: }
2042: /*@
2043: DMNetworkIsGhostVertex - Returns TRUE if the vertex is a ghost vertex
2045: Not Collective
2047: Input Parameters:
2048: + dm - the DMNetwork object
2049: - p - the vertex point
2051: Output Parameter:
2052: . isghost - TRUE if the vertex is a ghost point
2054: Level: beginner
2056: .seealso: DMNetworkGetConnectedVertices(), DMNetworkGetVertexRange(), DMNetworkIsSharedVertex()
2057: @*/
2058: PetscErrorCode DMNetworkIsGhostVertex(DM dm,PetscInt p,PetscBool *isghost)
2059: {
2061: DM_Network *network = (DM_Network*)dm->data;
2062: PetscInt offsetg;
2063: PetscSection sectiong;
2066: *isghost = PETSC_FALSE;
2067: DMGetGlobalSection(network->plex,§iong);
2068: PetscSectionGetOffset(sectiong,p,&offsetg);
2069: if (offsetg < 0) *isghost = PETSC_TRUE;
2070: return(0);
2071: }
2073: PetscErrorCode DMSetUp_Network(DM dm)
2074: {
2076: DM_Network *network=(DM_Network*)dm->data;
2079: DMNetworkComponentSetUp(dm);
2080: DMNetworkVariablesSetUp(dm);
2082: DMSetLocalSection(network->plex,network->DofSection);
2083: DMGetGlobalSection(network->plex,&network->GlobalDofSection);
2085: dm->setupcalled = PETSC_TRUE;
2087: /* View dmnetwork */
2088: DMViewFromOptions(dm,NULL,"-dmnetwork_view");
2089: return(0);
2090: }
2092: /*@
2093: DMNetworkHasJacobian - Sets global flag for using user's sub Jacobian matrices
2094: -- replaced by DMNetworkSetOption(network,userjacobian,PETSC_TURE)?
2096: Collective
2098: Input Parameters:
2099: + dm - the DMNetwork object
2100: . eflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for edges
2101: - vflg - turn the option on (PETSC_TRUE) or off (PETSC_FALSE) if user provides Jacobian for vertices
2103: Level: intermediate
2105: @*/
2106: PetscErrorCode DMNetworkHasJacobian(DM dm,PetscBool eflg,PetscBool vflg)
2107: {
2108: DM_Network *network=(DM_Network*)dm->data;
2110: PetscInt nVertices = network->nVertices;
2113: network->userEdgeJacobian = eflg;
2114: network->userVertexJacobian = vflg;
2116: if (eflg && !network->Je) {
2117: PetscCalloc1(3*network->nEdges,&network->Je);
2118: }
2120: if (vflg && !network->Jv && nVertices) {
2121: PetscInt i,*vptr,nedges,vStart=network->vStart;
2122: PetscInt nedges_total;
2123: const PetscInt *edges;
2125: /* count nvertex_total */
2126: nedges_total = 0;
2127: PetscMalloc1(nVertices+1,&vptr);
2129: vptr[0] = 0;
2130: for (i=0; i<nVertices; i++) {
2131: DMNetworkGetSupportingEdges(dm,i+vStart,&nedges,&edges);
2132: nedges_total += nedges;
2133: vptr[i+1] = vptr[i] + 2*nedges + 1;
2134: }
2136: PetscCalloc1(2*nedges_total+nVertices,&network->Jv);
2137: network->Jvptr = vptr;
2138: }
2139: return(0);
2140: }
2142: /*@
2143: DMNetworkEdgeSetMatrix - Sets user-provided Jacobian matrices for this edge to the network
2145: Not Collective
2147: Input Parameters:
2148: + dm - the DMNetwork object
2149: . p - the edge point
2150: - J - array (size = 3) of Jacobian submatrices for this edge point:
2151: J[0]: this edge
2152: J[1] and J[2]: connected vertices, obtained by calling DMNetworkGetConnectedVertices()
2154: Level: advanced
2156: .seealso: DMNetworkVertexSetMatrix()
2157: @*/
2158: PetscErrorCode DMNetworkEdgeSetMatrix(DM dm,PetscInt p,Mat J[])
2159: {
2160: DM_Network *network=(DM_Network*)dm->data;
2163: if (!network->Je) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkEdgeSetMatrix");
2165: if (J) {
2166: network->Je[3*p] = J[0];
2167: network->Je[3*p+1] = J[1];
2168: network->Je[3*p+2] = J[2];
2169: }
2170: return(0);
2171: }
2173: /*@
2174: DMNetworkVertexSetMatrix - Sets user-provided Jacobian matrix for this vertex to the network
2176: Not Collective
2178: Input Parameters:
2179: + dm - The DMNetwork object
2180: . p - the vertex point
2181: - J - array of Jacobian (size = 2*(num of supporting edges) + 1) submatrices for this vertex point:
2182: J[0]: this vertex
2183: J[1+2*i]: i-th supporting edge
2184: J[1+2*i+1]: i-th connected vertex
2186: Level: advanced
2188: .seealso: DMNetworkEdgeSetMatrix()
2189: @*/
2190: PetscErrorCode DMNetworkVertexSetMatrix(DM dm,PetscInt p,Mat J[])
2191: {
2193: DM_Network *network=(DM_Network*)dm->data;
2194: PetscInt i,*vptr,nedges,vStart=network->vStart;
2195: const PetscInt *edges;
2198: if (!network->Jv) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"Must call DMNetworkHasJacobian() collectively before calling DMNetworkVertexSetMatrix");
2200: if (J) {
2201: vptr = network->Jvptr;
2202: network->Jv[vptr[p-vStart]] = J[0]; /* Set Jacobian for this vertex */
2204: /* Set Jacobian for each supporting edge and connected vertex */
2205: DMNetworkGetSupportingEdges(dm,p,&nedges,&edges);
2206: for (i=1; i<=2*nedges; i++) network->Jv[vptr[p-vStart]+i] = J[i];
2207: }
2208: return(0);
2209: }
2211: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
2212: {
2214: PetscInt j;
2215: PetscScalar val=(PetscScalar)ncols;
2218: if (!ghost) {
2219: for (j=0; j<nrows; j++) {
2220: VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
2221: }
2222: } else {
2223: for (j=0; j<nrows; j++) {
2224: VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
2225: }
2226: }
2227: return(0);
2228: }
2230: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
2231: {
2233: PetscInt j,ncols_u;
2234: PetscScalar val;
2237: if (!ghost) {
2238: for (j=0; j<nrows; j++) {
2239: MatGetRow(Ju,j,&ncols_u,NULL,NULL);
2240: val = (PetscScalar)ncols_u;
2241: VecSetValues(vdnz,1,&rows[j],&val,ADD_VALUES);
2242: MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
2243: }
2244: } else {
2245: for (j=0; j<nrows; j++) {
2246: MatGetRow(Ju,j,&ncols_u,NULL,NULL);
2247: val = (PetscScalar)ncols_u;
2248: VecSetValues(vonz,1,&rows[j],&val,ADD_VALUES);
2249: MatRestoreRow(Ju,j,&ncols_u,NULL,NULL);
2250: }
2251: }
2252: return(0);
2253: }
2255: PETSC_STATIC_INLINE PetscErrorCode MatSetPreallocationblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscBool ghost,Vec vdnz,Vec vonz)
2256: {
2260: if (Ju) {
2261: MatSetPreallocationUserblock_private(Ju,nrows,rows,ncols,ghost,vdnz,vonz);
2262: } else {
2263: MatSetPreallocationDenseblock_private(nrows,rows,ncols,ghost,vdnz,vonz);
2264: }
2265: return(0);
2266: }
2268: PETSC_STATIC_INLINE PetscErrorCode MatSetDenseblock_private(PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
2269: {
2271: PetscInt j,*cols;
2272: PetscScalar *zeros;
2275: PetscCalloc2(ncols,&cols,nrows*ncols,&zeros);
2276: for (j=0; j<ncols; j++) cols[j] = j+ cstart;
2277: MatSetValues(*J,nrows,rows,ncols,cols,zeros,INSERT_VALUES);
2278: PetscFree2(cols,zeros);
2279: return(0);
2280: }
2282: PETSC_STATIC_INLINE PetscErrorCode MatSetUserblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
2283: {
2285: PetscInt j,M,N,row,col,ncols_u;
2286: const PetscInt *cols;
2287: PetscScalar zero=0.0;
2290: MatGetSize(Ju,&M,&N);
2291: if (nrows != M || ncols != N) SETERRQ4(PetscObjectComm((PetscObject)Ju),PETSC_ERR_USER,"%D by %D must equal %D by %D",nrows,ncols,M,N);
2293: for (row=0; row<nrows; row++) {
2294: MatGetRow(Ju,row,&ncols_u,&cols,NULL);
2295: for (j=0; j<ncols_u; j++) {
2296: col = cols[j] + cstart;
2297: MatSetValues(*J,1,&rows[row],1,&col,&zero,INSERT_VALUES);
2298: }
2299: MatRestoreRow(Ju,row,&ncols_u,&cols,NULL);
2300: }
2301: return(0);
2302: }
2304: PETSC_STATIC_INLINE PetscErrorCode MatSetblock_private(Mat Ju,PetscInt nrows,PetscInt *rows,PetscInt ncols,PetscInt cstart,Mat *J)
2305: {
2309: if (Ju) {
2310: MatSetUserblock_private(Ju,nrows,rows,ncols,cstart,J);
2311: } else {
2312: MatSetDenseblock_private(nrows,rows,ncols,cstart,J);
2313: }
2314: return(0);
2315: }
2317: /* 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.
2318: */
2319: PetscErrorCode CreateSubGlobalToLocalMapping_private(PetscSection globalsec, PetscSection localsec, ISLocalToGlobalMapping *ltog)
2320: {
2322: PetscInt i,size,dof;
2323: PetscInt *glob2loc;
2326: PetscSectionGetStorageSize(localsec,&size);
2327: PetscMalloc1(size,&glob2loc);
2329: for (i = 0; i < size; i++) {
2330: PetscSectionGetOffset(globalsec,i,&dof);
2331: dof = (dof >= 0) ? dof : -(dof + 1);
2332: glob2loc[i] = dof;
2333: }
2335: ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,size,glob2loc,PETSC_OWN_POINTER,ltog);
2336: #if 0
2337: PetscIntView(size,glob2loc,PETSC_VIEWER_STDOUT_WORLD);
2338: #endif
2339: return(0);
2340: }
2342: #include <petsc/private/matimpl.h>
2344: PetscErrorCode DMCreateMatrix_Network_Nest(DM dm,Mat *J)
2345: {
2347: DM_Network *network = (DM_Network*)dm->data;
2348: PetscInt eDof,vDof;
2349: Mat j11,j12,j21,j22,bA[2][2];
2350: MPI_Comm comm;
2351: ISLocalToGlobalMapping eISMap,vISMap;
2354: PetscObjectGetComm((PetscObject)dm,&comm);
2356: PetscSectionGetConstrainedStorageSize(network->edge.GlobalDofSection,&eDof);
2357: PetscSectionGetConstrainedStorageSize(network->vertex.GlobalDofSection,&vDof);
2359: MatCreate(comm, &j11);
2360: MatSetSizes(j11, eDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
2361: MatSetType(j11, MATMPIAIJ);
2363: MatCreate(comm, &j12);
2364: MatSetSizes(j12, eDof, vDof, PETSC_DETERMINE ,PETSC_DETERMINE);
2365: MatSetType(j12, MATMPIAIJ);
2367: MatCreate(comm, &j21);
2368: MatSetSizes(j21, vDof, eDof, PETSC_DETERMINE, PETSC_DETERMINE);
2369: MatSetType(j21, MATMPIAIJ);
2371: MatCreate(comm, &j22);
2372: MatSetSizes(j22, vDof, vDof, PETSC_DETERMINE, PETSC_DETERMINE);
2373: MatSetType(j22, MATMPIAIJ);
2375: bA[0][0] = j11;
2376: bA[0][1] = j12;
2377: bA[1][0] = j21;
2378: bA[1][1] = j22;
2380: CreateSubGlobalToLocalMapping_private(network->edge.GlobalDofSection,network->edge.DofSection,&eISMap);
2381: CreateSubGlobalToLocalMapping_private(network->vertex.GlobalDofSection,network->vertex.DofSection,&vISMap);
2383: MatSetLocalToGlobalMapping(j11,eISMap,eISMap);
2384: MatSetLocalToGlobalMapping(j12,eISMap,vISMap);
2385: MatSetLocalToGlobalMapping(j21,vISMap,eISMap);
2386: MatSetLocalToGlobalMapping(j22,vISMap,vISMap);
2388: MatSetUp(j11);
2389: MatSetUp(j12);
2390: MatSetUp(j21);
2391: MatSetUp(j22);
2393: MatCreateNest(comm,2,NULL,2,NULL,&bA[0][0],J);
2394: MatSetUp(*J);
2395: MatNestSetVecType(*J,VECNEST);
2396: MatDestroy(&j11);
2397: MatDestroy(&j12);
2398: MatDestroy(&j21);
2399: MatDestroy(&j22);
2401: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
2402: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
2403: MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2405: /* Free structures */
2406: ISLocalToGlobalMappingDestroy(&eISMap);
2407: ISLocalToGlobalMappingDestroy(&vISMap);
2408: return(0);
2409: }
2411: PetscErrorCode DMCreateMatrix_Network(DM dm,Mat *J)
2412: {
2414: DM_Network *network = (DM_Network*)dm->data;
2415: PetscInt eStart,eEnd,vStart,vEnd,rstart,nrows,*rows,localSize;
2416: PetscInt cstart,ncols,j,e,v;
2417: PetscBool ghost,ghost_vc,ghost2,isNest;
2418: Mat Juser;
2419: PetscSection sectionGlobal;
2420: PetscInt nedges,*vptr=NULL,vc,*rows_v; /* suppress maybe-uninitialized warning */
2421: const PetscInt *edges,*cone;
2422: MPI_Comm comm;
2423: MatType mtype;
2424: Vec vd_nz,vo_nz;
2425: PetscInt *dnnz,*onnz;
2426: PetscScalar *vdnz,*vonz;
2429: mtype = dm->mattype;
2430: PetscStrcmp(mtype,MATNEST,&isNest);
2431: if (isNest) {
2432: DMCreateMatrix_Network_Nest(dm,J);
2433: MatSetDM(*J,dm);
2434: return(0);
2435: }
2437: if (!network->userEdgeJacobian && !network->userVertexJacobian) {
2438: /* user does not provide Jacobian blocks */
2439: DMCreateMatrix_Plex(network->plex,J);
2440: MatSetDM(*J,dm);
2441: return(0);
2442: }
2444: MatCreate(PetscObjectComm((PetscObject)dm),J);
2445: DMGetGlobalSection(network->plex,§ionGlobal);
2446: PetscSectionGetConstrainedStorageSize(sectionGlobal,&localSize);
2447: MatSetSizes(*J,localSize,localSize,PETSC_DETERMINE,PETSC_DETERMINE);
2449: MatSetType(*J,MATAIJ);
2450: MatSetFromOptions(*J);
2452: /* (1) Set matrix preallocation */
2453: /*------------------------------*/
2454: PetscObjectGetComm((PetscObject)dm,&comm);
2455: VecCreate(comm,&vd_nz);
2456: VecSetSizes(vd_nz,localSize,PETSC_DECIDE);
2457: VecSetFromOptions(vd_nz);
2458: VecSet(vd_nz,0.0);
2459: VecDuplicate(vd_nz,&vo_nz);
2461: /* Set preallocation for edges */
2462: /*-----------------------------*/
2463: DMNetworkGetEdgeRange(dm,&eStart,&eEnd);
2465: PetscMalloc1(localSize,&rows);
2466: for (e=eStart; e<eEnd; e++) {
2467: /* Get row indices */
2468: DMNetworkGetGlobalVecOffset(dm,e,ALL_COMPONENTS,&rstart);
2469: PetscSectionGetDof(network->DofSection,e,&nrows);
2470: if (nrows) {
2471: for (j=0; j<nrows; j++) rows[j] = j + rstart;
2473: /* Set preallocation for connected vertices */
2474: DMNetworkGetConnectedVertices(dm,e,&cone);
2475: for (v=0; v<2; v++) {
2476: PetscSectionGetDof(network->DofSection,cone[v],&ncols);
2478: if (network->Je) {
2479: Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
2480: } else Juser = NULL;
2481: DMNetworkIsGhostVertex(dm,cone[v],&ghost);
2482: MatSetPreallocationblock_private(Juser,nrows,rows,ncols,ghost,vd_nz,vo_nz);
2483: }
2485: /* Set preallocation for edge self */
2486: cstart = rstart;
2487: if (network->Je) {
2488: Juser = network->Je[3*e]; /* Jacobian(e,e) */
2489: } else Juser = NULL;
2490: MatSetPreallocationblock_private(Juser,nrows,rows,nrows,PETSC_FALSE,vd_nz,vo_nz);
2491: }
2492: }
2494: /* Set preallocation for vertices */
2495: /*--------------------------------*/
2496: DMNetworkGetVertexRange(dm,&vStart,&vEnd);
2497: if (vEnd - vStart) vptr = network->Jvptr;
2499: for (v=vStart; v<vEnd; v++) {
2500: /* Get row indices */
2501: DMNetworkGetGlobalVecOffset(dm,v,ALL_COMPONENTS,&rstart);
2502: PetscSectionGetDof(network->DofSection,v,&nrows);
2503: if (!nrows) continue;
2505: DMNetworkIsGhostVertex(dm,v,&ghost);
2506: if (ghost) {
2507: PetscMalloc1(nrows,&rows_v);
2508: } else {
2509: rows_v = rows;
2510: }
2512: for (j=0; j<nrows; j++) rows_v[j] = j + rstart;
2514: /* Get supporting edges and connected vertices */
2515: DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
2517: for (e=0; e<nedges; e++) {
2518: /* Supporting edges */
2519: DMNetworkGetGlobalVecOffset(dm,edges[e],ALL_COMPONENTS,&cstart);
2520: PetscSectionGetDof(network->DofSection,edges[e],&ncols);
2522: if (network->Jv) {
2523: Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */
2524: } else Juser = NULL;
2525: MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost,vd_nz,vo_nz);
2527: /* Connected vertices */
2528: DMNetworkGetConnectedVertices(dm,edges[e],&cone);
2529: vc = (v == cone[0]) ? cone[1]:cone[0];
2530: DMNetworkIsGhostVertex(dm,vc,&ghost_vc);
2532: PetscSectionGetDof(network->DofSection,vc,&ncols);
2534: if (network->Jv) {
2535: Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
2536: } else Juser = NULL;
2537: if (ghost_vc||ghost) {
2538: ghost2 = PETSC_TRUE;
2539: } else {
2540: ghost2 = PETSC_FALSE;
2541: }
2542: MatSetPreallocationblock_private(Juser,nrows,rows_v,ncols,ghost2,vd_nz,vo_nz);
2543: }
2545: /* Set preallocation for vertex self */
2546: DMNetworkIsGhostVertex(dm,v,&ghost);
2547: if (!ghost) {
2548: DMNetworkGetGlobalVecOffset(dm,v,ALL_COMPONENTS,&cstart);
2549: if (network->Jv) {
2550: Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
2551: } else Juser = NULL;
2552: MatSetPreallocationblock_private(Juser,nrows,rows_v,nrows,PETSC_FALSE,vd_nz,vo_nz);
2553: }
2554: if (ghost) {
2555: PetscFree(rows_v);
2556: }
2557: }
2559: VecAssemblyBegin(vd_nz);
2560: VecAssemblyBegin(vo_nz);
2562: PetscMalloc2(localSize,&dnnz,localSize,&onnz);
2564: VecAssemblyEnd(vd_nz);
2565: VecAssemblyEnd(vo_nz);
2567: VecGetArray(vd_nz,&vdnz);
2568: VecGetArray(vo_nz,&vonz);
2569: for (j=0; j<localSize; j++) {
2570: dnnz[j] = (PetscInt)PetscRealPart(vdnz[j]);
2571: onnz[j] = (PetscInt)PetscRealPart(vonz[j]);
2572: }
2573: VecRestoreArray(vd_nz,&vdnz);
2574: VecRestoreArray(vo_nz,&vonz);
2575: VecDestroy(&vd_nz);
2576: VecDestroy(&vo_nz);
2578: MatSeqAIJSetPreallocation(*J,0,dnnz);
2579: MatMPIAIJSetPreallocation(*J,0,dnnz,0,onnz);
2580: MatSetOption(*J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2582: PetscFree2(dnnz,onnz);
2584: /* (2) Set matrix entries for edges */
2585: /*----------------------------------*/
2586: for (e=eStart; e<eEnd; e++) {
2587: /* Get row indices */
2588: DMNetworkGetGlobalVecOffset(dm,e,ALL_COMPONENTS,&rstart);
2589: PetscSectionGetDof(network->DofSection,e,&nrows);
2590: if (nrows) {
2591: for (j=0; j<nrows; j++) rows[j] = j + rstart;
2593: /* Set matrix entries for connected vertices */
2594: DMNetworkGetConnectedVertices(dm,e,&cone);
2595: for (v=0; v<2; v++) {
2596: DMNetworkGetGlobalVecOffset(dm,cone[v],ALL_COMPONENTS,&cstart);
2597: PetscSectionGetDof(network->DofSection,cone[v],&ncols);
2599: if (network->Je) {
2600: Juser = network->Je[3*e+1+v]; /* Jacobian(e,v) */
2601: } else Juser = NULL;
2602: MatSetblock_private(Juser,nrows,rows,ncols,cstart,J);
2603: }
2605: /* Set matrix entries for edge self */
2606: cstart = rstart;
2607: if (network->Je) {
2608: Juser = network->Je[3*e]; /* Jacobian(e,e) */
2609: } else Juser = NULL;
2610: MatSetblock_private(Juser,nrows,rows,nrows,cstart,J);
2611: }
2612: }
2614: /* Set matrix entries for vertices */
2615: /*---------------------------------*/
2616: for (v=vStart; v<vEnd; v++) {
2617: /* Get row indices */
2618: DMNetworkGetGlobalVecOffset(dm,v,ALL_COMPONENTS,&rstart);
2619: PetscSectionGetDof(network->DofSection,v,&nrows);
2620: if (!nrows) continue;
2622: DMNetworkIsGhostVertex(dm,v,&ghost);
2623: if (ghost) {
2624: PetscMalloc1(nrows,&rows_v);
2625: } else {
2626: rows_v = rows;
2627: }
2628: for (j=0; j<nrows; j++) rows_v[j] = j + rstart;
2630: /* Get supporting edges and connected vertices */
2631: DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);
2633: for (e=0; e<nedges; e++) {
2634: /* Supporting edges */
2635: DMNetworkGetGlobalVecOffset(dm,edges[e],ALL_COMPONENTS,&cstart);
2636: PetscSectionGetDof(network->DofSection,edges[e],&ncols);
2638: if (network->Jv) {
2639: Juser = network->Jv[vptr[v-vStart]+2*e+1]; /* Jacobian(v,e) */
2640: } else Juser = NULL;
2641: MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);
2643: /* Connected vertices */
2644: DMNetworkGetConnectedVertices(dm,edges[e],&cone);
2645: vc = (v == cone[0]) ? cone[1]:cone[0];
2647: DMNetworkGetGlobalVecOffset(dm,vc,ALL_COMPONENTS,&cstart);
2648: PetscSectionGetDof(network->DofSection,vc,&ncols);
2650: if (network->Jv) {
2651: Juser = network->Jv[vptr[v-vStart]+2*e+2]; /* Jacobian(v,vc) */
2652: } else Juser = NULL;
2653: MatSetblock_private(Juser,nrows,rows_v,ncols,cstart,J);
2654: }
2656: /* Set matrix entries for vertex self */
2657: if (!ghost) {
2658: DMNetworkGetGlobalVecOffset(dm,v,ALL_COMPONENTS,&cstart);
2659: if (network->Jv) {
2660: Juser = network->Jv[vptr[v-vStart]]; /* Jacobian(v,v) */
2661: } else Juser = NULL;
2662: MatSetblock_private(Juser,nrows,rows_v,nrows,cstart,J);
2663: }
2664: if (ghost) {
2665: PetscFree(rows_v);
2666: }
2667: }
2668: PetscFree(rows);
2670: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
2671: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
2673: MatSetDM(*J,dm);
2674: return(0);
2675: }
2677: PetscErrorCode DMDestroy_Network(DM dm)
2678: {
2680: DM_Network *network = (DM_Network*)dm->data;
2681: PetscInt j,np;
2684: if (--network->refct > 0) return(0);
2685: if (network->Je) {
2686: PetscFree(network->Je);
2687: }
2688: if (network->Jv) {
2689: PetscFree(network->Jvptr);
2690: PetscFree(network->Jv);
2691: }
2693: ISLocalToGlobalMappingDestroy(&network->vertex.mapping);
2694: PetscSectionDestroy(&network->vertex.DofSection);
2695: PetscSectionDestroy(&network->vertex.GlobalDofSection);
2696: if (network->vltog) {
2697: PetscFree(network->vltog);
2698: }
2699: if (network->vertex.sf) {
2700: PetscSFDestroy(&network->vertex.sf);
2701: }
2702: /* edge */
2703: ISLocalToGlobalMappingDestroy(&network->edge.mapping);
2704: PetscSectionDestroy(&network->edge.DofSection);
2705: PetscSectionDestroy(&network->edge.GlobalDofSection);
2706: if (network->edge.sf) {
2707: PetscSFDestroy(&network->edge.sf);
2708: }
2709: DMDestroy(&network->plex);
2710: PetscSectionDestroy(&network->DataSection);
2711: PetscSectionDestroy(&network->DofSection);
2713: for (j=0; j<network->Nsvtx; j++) {
2714: PetscFree(network->svtx[j].sv);
2715: }
2716: if (network->svtx) {PetscFree(network->svtx);}
2717: PetscFree2(network->subnetedge,network->subnetvtx);
2719: PetscTableDestroy(&network->svtable);
2720: PetscFree(network->subnet);
2721: PetscFree(network->component);
2722: PetscFree(network->componentdataarray);
2724: if (network->header) {
2725: np = network->pEnd - network->pStart;
2726: for (j=0; j < np; j++) {
2727: PetscFree5(network->header[j].size,network->header[j].key,network->header[j].offset,network->header[j].nvar,network->header[j].offsetvarrel);
2728: PetscFree(network->cvalue[j].data);
2729: }
2730: PetscFree2(network->header,network->cvalue);
2731: }
2732: PetscFree(network);
2733: return(0);
2734: }
2736: PetscErrorCode DMView_Network(DM dm,PetscViewer viewer)
2737: {
2739: PetscBool iascii;
2740: PetscMPIInt rank;
2743: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() first");
2744: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
2747: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2748: if (iascii) {
2749: const PetscInt *cone,*vtx,*edges;
2750: PetscInt vfrom,vto,i,j,nv,ne,ncv,p,nsubnet;
2751: DM_Network *network = (DM_Network*)dm->data;
2753: nsubnet = network->Nsubnet; /* num of subnetworks */
2754: if (rank == 0) {
2755: PetscPrintf(PETSC_COMM_SELF," NSubnets: %D; NEdges: %D; NVertices: %D; NSharedVertices: %D.\n",nsubnet,network->NEdges,network->NVertices,network->Nsvtx);
2756: }
2758: DMNetworkGetSharedVertices(dm,&ncv,&vtx);
2759: PetscViewerASCIIPushSynchronized(viewer);
2760: PetscViewerASCIISynchronizedPrintf(viewer, " [%d] nEdges: %D; nVertices: %D; nSharedVertices: %D\n",rank,network->nEdges,network->nVertices,ncv);
2762: for (i=0; i<nsubnet; i++) {
2763: DMNetworkGetSubnetwork(dm,i,&nv,&ne,&vtx,&edges);
2764: if (ne) {
2765: PetscViewerASCIISynchronizedPrintf(viewer, " Subnet %D: nEdges %D, nVertices(include shared vertices) %D\n",i,ne,nv);
2766: for (j=0; j<ne; j++) {
2767: p = edges[j];
2768: DMNetworkGetConnectedVertices(dm,p,&cone);
2769: DMNetworkGetGlobalVertexIndex(dm,cone[0],&vfrom);
2770: DMNetworkGetGlobalVertexIndex(dm,cone[1],&vto);
2771: DMNetworkGetGlobalEdgeIndex(dm,edges[j],&p);
2772: PetscViewerASCIISynchronizedPrintf(viewer, " edge %D: %D ----> %D\n",p,vfrom,vto);
2773: }
2774: }
2775: }
2777: /* Shared vertices */
2778: DMNetworkGetSharedVertices(dm,&ncv,&vtx);
2779: if (ncv) {
2780: SVtx *svtx = network->svtx;
2781: PetscInt gidx,svtx_idx,nvto,vfrom_net,vfrom_idx,*svto;
2782: PetscBool ghost;
2783: PetscViewerASCIISynchronizedPrintf(viewer, " SharedVertices:\n");
2784: for (i=0; i<ncv; i++) {
2785: DMNetworkIsGhostVertex(dm,vtx[i],&ghost);
2786: if (ghost) continue;
2788: DMNetworkGetGlobalVertexIndex(dm,vtx[i],&gidx);
2789: PetscTableFind(network->svtable,gidx+1,&svtx_idx);
2790: svtx_idx--;
2791: nvto = svtx[svtx_idx].n;
2793: vfrom_net = svtx[svtx_idx].sv[0];
2794: vfrom_idx = svtx[svtx_idx].sv[1];
2795: PetscViewerASCIISynchronizedPrintf(viewer, " svtx %D: global index %D, subnet[%D].%D ---->\n",i,gidx,vfrom_net,vfrom_idx);
2796: for (j=1; j<nvto; j++) {
2797: svto = svtx[svtx_idx].sv + 2*j;
2798: PetscViewerASCIISynchronizedPrintf(viewer, " ----> subnet[%D].%D\n",svto[0],svto[1]);
2799: }
2800: }
2801: }
2802: PetscViewerFlush(viewer);
2803: PetscViewerASCIIPopSynchronized(viewer);
2804: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMNetwork writing", ((PetscObject)viewer)->type_name);
2805: return(0);
2806: }
2808: PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l)
2809: {
2811: DM_Network *network = (DM_Network*)dm->data;
2814: DMGlobalToLocalBegin(network->plex,g,mode,l);
2815: return(0);
2816: }
2818: PetscErrorCode DMGlobalToLocalEnd_Network(DM dm, Vec g, InsertMode mode, Vec l)
2819: {
2821: DM_Network *network = (DM_Network*)dm->data;
2824: DMGlobalToLocalEnd(network->plex,g,mode,l);
2825: return(0);
2826: }
2828: PetscErrorCode DMLocalToGlobalBegin_Network(DM dm, Vec l, InsertMode mode, Vec g)
2829: {
2831: DM_Network *network = (DM_Network*)dm->data;
2834: DMLocalToGlobalBegin(network->plex,l,mode,g);
2835: return(0);
2836: }
2838: PetscErrorCode DMLocalToGlobalEnd_Network(DM dm, Vec l, InsertMode mode, Vec g)
2839: {
2841: DM_Network *network = (DM_Network*)dm->data;
2844: DMLocalToGlobalEnd(network->plex,l,mode,g);
2845: return(0);
2846: }
2848: /*@
2849: DMNetworkGetVertexLocalToGlobalOrdering - Get vertex global index
2851: Not collective
2853: Input Parameters:
2854: + dm - the dm object
2855: - vloc - local vertex ordering, start from 0
2857: Output Parameters:
2858: . vg - global vertex ordering, start from 0
2860: Level: advanced
2862: .seealso: DMNetworkSetVertexLocalToGlobalOrdering()
2863: @*/
2864: PetscErrorCode DMNetworkGetVertexLocalToGlobalOrdering(DM dm,PetscInt vloc,PetscInt *vg)
2865: {
2866: DM_Network *network = (DM_Network*)dm->data;
2867: PetscInt *vltog = network->vltog;
2870: if (!vltog) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkSetVertexLocalToGlobalOrdering() first");
2871: *vg = vltog[vloc];
2872: return(0);
2873: }
2875: /*@
2876: DMNetworkSetVertexLocalToGlobalOrdering - Create and setup vertex local to global map
2878: Collective
2880: Input Parameters:
2881: . dm - the dm object
2883: Level: advanced
2885: .seealso: DMNetworkGetGlobalVertexIndex()
2886: @*/
2887: PetscErrorCode DMNetworkSetVertexLocalToGlobalOrdering(DM dm)
2888: {
2889: PetscErrorCode ierr;
2890: DM_Network *network=(DM_Network*)dm->data;
2891: MPI_Comm comm;
2892: PetscMPIInt rank,size,*displs=NULL,*recvcounts=NULL,remoterank;
2893: PetscBool ghost;
2894: PetscInt *vltog,nroots,nleaves,i,*vrange,k,N,lidx;
2895: const PetscSFNode *iremote;
2896: PetscSF vsf;
2897: Vec Vleaves,Vleaves_seq;
2898: VecScatter ctx;
2899: PetscScalar *varr,val;
2900: const PetscScalar *varr_read;
2903: PetscObjectGetComm((PetscObject)dm,&comm);
2904: MPI_Comm_size(comm,&size);
2905: MPI_Comm_rank(comm,&rank);
2907: if (size == 1) {
2908: nroots = network->vEnd - network->vStart;
2909: PetscMalloc1(nroots, &vltog);
2910: for (i=0; i<nroots; i++) vltog[i] = i;
2911: network->vltog = vltog;
2912: return(0);
2913: }
2915: if (!network->distributecalled) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE,"Must call DMNetworkDistribute() first");
2916: if (network->vltog) {
2917: PetscFree(network->vltog);
2918: }
2920: DMNetworkSetSubMap_private(network->vStart,network->vEnd,&network->vertex.mapping);
2921: PetscSFGetSubSF(network->plex->sf, network->vertex.mapping, &network->vertex.sf);
2922: vsf = network->vertex.sf;
2924: PetscMalloc3(size+1,&vrange,size,&displs,size,&recvcounts);
2925: PetscSFGetGraph(vsf,&nroots,&nleaves,NULL,&iremote);
2927: for (i=0; i<size; i++) { displs[i] = i; recvcounts[i] = 1;}
2929: i = nroots - nleaves; /* local number of vertices, excluding ghosts */
2930: vrange[0] = 0;
2931: MPI_Allgatherv(&i,1,MPIU_INT,vrange+1,recvcounts,displs,MPIU_INT,comm);
2932: for (i=2; i<size+1; i++) {vrange[i] += vrange[i-1];}
2934: PetscMalloc1(nroots, &vltog);
2935: network->vltog = vltog;
2937: /* Set vltog for non-ghost vertices */
2938: k = 0;
2939: for (i=0; i<nroots; i++) {
2940: DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);
2941: if (ghost) continue;
2942: vltog[i] = vrange[rank] + k++;
2943: }
2944: PetscFree3(vrange,displs,recvcounts);
2946: /* Set vltog for ghost vertices */
2947: /* (a) create parallel Vleaves and sequential Vleaves_seq to convert local iremote[*].index to global index */
2948: VecCreate(comm,&Vleaves);
2949: VecSetSizes(Vleaves,2*nleaves,PETSC_DETERMINE);
2950: VecSetFromOptions(Vleaves);
2951: VecGetArray(Vleaves,&varr);
2952: for (i=0; i<nleaves; i++) {
2953: varr[2*i] = (PetscScalar)(iremote[i].rank); /* rank of remote process */
2954: varr[2*i+1] = (PetscScalar)(iremote[i].index); /* local index in remote process */
2955: }
2956: VecRestoreArray(Vleaves,&varr);
2958: /* (b) scatter local info to remote processes via VecScatter() */
2959: VecScatterCreateToAll(Vleaves,&ctx,&Vleaves_seq);
2960: VecScatterBegin(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);
2961: VecScatterEnd(ctx,Vleaves,Vleaves_seq,INSERT_VALUES,SCATTER_FORWARD);
2963: /* (c) convert local indices to global indices in parallel vector Vleaves */
2964: VecGetSize(Vleaves_seq,&N);
2965: VecGetArrayRead(Vleaves_seq,&varr_read);
2966: for (i=0; i<N; i+=2) {
2967: remoterank = (PetscMPIInt)PetscRealPart(varr_read[i]);
2968: if (remoterank == rank) {
2969: k = i+1; /* row number */
2970: lidx = (PetscInt)PetscRealPart(varr_read[i+1]);
2971: val = (PetscScalar)vltog[lidx]; /* global index for non-ghost vertex computed above */
2972: VecSetValues(Vleaves,1,&k,&val,INSERT_VALUES);
2973: }
2974: }
2975: VecRestoreArrayRead(Vleaves_seq,&varr_read);
2976: VecAssemblyBegin(Vleaves);
2977: VecAssemblyEnd(Vleaves);
2979: /* (d) Set vltog for ghost vertices by copying local values of Vleaves */
2980: VecGetArrayRead(Vleaves,&varr_read);
2981: k = 0;
2982: for (i=0; i<nroots; i++) {
2983: DMNetworkIsGhostVertex(dm,i+network->vStart,&ghost);
2984: if (!ghost) continue;
2985: vltog[i] = (PetscInt)PetscRealPart(varr_read[2*k+1]); k++;
2986: }
2987: VecRestoreArrayRead(Vleaves,&varr_read);
2989: VecDestroy(&Vleaves);
2990: VecDestroy(&Vleaves_seq);
2991: VecScatterDestroy(&ctx);
2992: return(0);
2993: }
2995: PETSC_STATIC_INLINE PetscErrorCode DMISAddSize_private(DM_Network *network,PetscInt p,PetscInt numkeys,PetscInt keys[],PetscInt blocksize[],PetscInt nselectedvar[],PetscInt *nidx)
2996: {
2997: PetscErrorCode ierr;
2998: PetscInt i,j,ncomps,nvar,key,offset=0;
2999: DMNetworkComponentHeader header;
3002: PetscSectionGetOffset(network->DataSection,p,&offset);
3003: ncomps = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata;
3004: header = (DMNetworkComponentHeader)(network->componentdataarray+offset);
3006: for (i=0; i<ncomps; i++) {
3007: key = header->key[i];
3008: nvar = header->nvar[i];
3009: for (j=0; j<numkeys; j++) {
3010: if (key == keys[j]) {
3011: if (!blocksize || blocksize[j] == -1) {
3012: *nidx += nvar;
3013: } else {
3014: *nidx += nselectedvar[j]*nvar/blocksize[j];
3015: }
3016: }
3017: }
3018: }
3019: return(0);
3020: }
3022: PETSC_STATIC_INLINE PetscErrorCode DMISComputeIdx_private(DM dm,PetscInt p,PetscInt numkeys,PetscInt keys[],PetscInt blocksize[],PetscInt nselectedvar[],PetscInt *selectedvar[],PetscInt *ii,PetscInt *idx)
3023: {
3024: PetscErrorCode ierr;
3025: PetscInt i,j,ncomps,nvar,key,offsetg,k,k1,offset=0;
3026: DM_Network *network = (DM_Network*)dm->data;
3027: DMNetworkComponentHeader header;
3030: PetscSectionGetOffset(network->DataSection,p,&offset);
3031: ncomps = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata;
3032: header = (DMNetworkComponentHeader)(network->componentdataarray+offset);
3034: for (i=0; i<ncomps; i++) {
3035: key = header->key[i];
3036: nvar = header->nvar[i];
3037: for (j=0; j<numkeys; j++) {
3038: if (key != keys[j]) continue;
3040: DMNetworkGetGlobalVecOffset(dm,p,i,&offsetg);
3041: if (!blocksize || blocksize[j] == -1) {
3042: for (k=0; k<nvar; k++) idx[(*ii)++] = offsetg + k;
3043: } else {
3044: for (k=0; k<nvar; k+=blocksize[j]) {
3045: for (k1=0; k1<nselectedvar[j]; k1++) idx[(*ii)++] = offsetg + k + selectedvar[j][k1];
3046: }
3047: }
3048: }
3049: }
3050: return(0);
3051: }
3053: /*@
3054: DMNetworkCreateIS - Create an index set object from the global vector of the network
3056: Collective
3058: Input Parameters:
3059: + dm - DMNetwork object
3060: . numkeys - number of keys
3061: . keys - array of keys that define the components of the variables you wish to extract
3062: . blocksize - block size of the variables associated to the component
3063: . nselectedvar - number of variables in each block to select
3064: - selectedvar - the offset into the block of each variable in each block to select
3066: Output Parameters:
3067: . is - the index set
3069: Level: Advanced
3071: Notes:
3072: Use blocksize[i] of -1 to indicate select all the variables of the i-th component, for which nselectvar[i] and selectedvar[i] are ignored. Use NULL, NULL, NULL to indicate for all selected components one wishes to obtain all the values of that component. For example, DMNetworkCreateIS(dm,1,&key,NULL,NULL,NULL,&is) will return an is that extracts all the variables for the 0-th component.
3074: .seealso: DMNetworkCreate(), ISCreateGeneral(), DMNetworkCreateLocalIS()
3075: @*/
3076: PetscErrorCode DMNetworkCreateIS(DM dm,PetscInt numkeys,PetscInt keys[],PetscInt blocksize[],PetscInt nselectedvar[],PetscInt *selectedvar[],IS *is)
3077: {
3079: MPI_Comm comm;
3080: DM_Network *network = (DM_Network*)dm->data;
3081: PetscInt i,p,estart,eend,vstart,vend,nidx,*idx;
3082: PetscBool ghost;
3085: PetscObjectGetComm((PetscObject)dm,&comm);
3087: /* Check input parameters */
3088: for (i=0; i<numkeys; i++) {
3089: if (!blocksize || blocksize[i] == -1) continue;
3090: if (nselectedvar[i] > blocksize[i]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"number of selectedvariables %D cannot be larger than blocksize %D",nselectedvar[i],blocksize[i]);
3091: }
3093: DMNetworkGetEdgeRange(dm,&estart,&eend);
3094: DMNetworkGetVertexRange(dm,&vstart,&vend);
3096: /* Get local number of idx */
3097: nidx = 0;
3098: for (p=estart; p<eend; p++) {
3099: DMISAddSize_private(network,p,numkeys,keys,blocksize,nselectedvar,&nidx);
3100: }
3101: for (p=vstart; p<vend; p++) {
3102: DMNetworkIsGhostVertex(dm,p,&ghost);
3103: if (ghost) continue;
3104: DMISAddSize_private(network,p,numkeys,keys,blocksize,nselectedvar,&nidx);
3105: }
3107: /* Compute idx */
3108: PetscMalloc1(nidx,&idx);
3109: i = 0;
3110: for (p=estart; p<eend; p++) {
3111: DMISComputeIdx_private(dm,p,numkeys,keys,blocksize,nselectedvar,selectedvar,&i,idx);
3112: }
3113: for (p=vstart; p<vend; p++) {
3114: DMNetworkIsGhostVertex(dm,p,&ghost);
3115: if (ghost) continue;
3116: DMISComputeIdx_private(dm,p,numkeys,keys,blocksize,nselectedvar,selectedvar,&i,idx);
3117: }
3119: /* Create is */
3120: ISCreateGeneral(comm,nidx,idx,PETSC_COPY_VALUES,is);
3121: PetscFree(idx);
3122: return(0);
3123: }
3125: PETSC_STATIC_INLINE PetscErrorCode DMISComputeLocalIdx_private(DM dm,PetscInt p,PetscInt numkeys,PetscInt keys[],PetscInt blocksize[],PetscInt nselectedvar[],PetscInt *selectedvar[],PetscInt *ii,PetscInt *idx)
3126: {
3127: PetscErrorCode ierr;
3128: PetscInt i,j,ncomps,nvar,key,offsetl,k,k1,offset=0;
3129: DM_Network *network = (DM_Network*)dm->data;
3130: DMNetworkComponentHeader header;
3133: PetscSectionGetOffset(network->DataSection,p,&offset);
3134: ncomps = ((DMNetworkComponentHeader)(network->componentdataarray+offset))->ndata;
3135: header = (DMNetworkComponentHeader)(network->componentdataarray+offset);
3137: for (i=0; i<ncomps; i++) {
3138: key = header->key[i];
3139: nvar = header->nvar[i];
3140: for (j=0; j<numkeys; j++) {
3141: if (key != keys[j]) continue;
3143: DMNetworkGetLocalVecOffset(dm,p,i,&offsetl);
3144: if (!blocksize || blocksize[j] == -1) {
3145: for (k=0; k<nvar; k++) idx[(*ii)++] = offsetl + k;
3146: } else {
3147: for (k=0; k<nvar; k+=blocksize[j]) {
3148: for (k1=0; k1<nselectedvar[j]; k1++) idx[(*ii)++] = offsetl + k + selectedvar[j][k1];
3149: }
3150: }
3151: }
3152: }
3153: return(0);
3154: }
3156: /*@
3157: DMNetworkCreateLocalIS - Create an index set object from the local vector of the network
3159: Not collective
3161: Input Parameters:
3162: + dm - DMNetwork object
3163: . numkeys - number of keys
3164: . keys - array of keys that define the components of the variables you wish to extract
3165: . blocksize - block size of the variables associated to the component
3166: . nselectedvar - number of variables in each block to select
3167: - selectedvar - the offset into the block of each variable in each block to select
3169: Output Parameters:
3170: . is - the index set
3172: Level: Advanced
3174: Notes:
3175: Use blocksize[i] of -1 to indicate select all the variables of the i-th component, for which nselectvar[i] and selectedvar[i] are ignored. Use NULL, NULL, NULL to indicate for all selected components one wishes to obtain all the values of that component. For example, DMNetworkCreateLocalIS(dm,1,&key,NULL,NULL,NULL,&is) will return an is that extracts all the variables for the 0-th component.
3177: .seealso: DMNetworkCreate(), DMNetworkCreateIS, ISCreateGeneral()
3178: @*/
3179: PetscErrorCode DMNetworkCreateLocalIS(DM dm,PetscInt numkeys,PetscInt keys[],PetscInt blocksize[],PetscInt nselectedvar[],PetscInt *selectedvar[],IS *is)
3180: {
3182: DM_Network *network = (DM_Network*)dm->data;
3183: PetscInt i,p,pstart,pend,nidx,*idx;
3186: /* Check input parameters */
3187: for (i=0; i<numkeys; i++) {
3188: if (!blocksize || blocksize[i] == -1) continue;
3189: if (nselectedvar[i] > blocksize[i]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"number of selectedvariables %D cannot be larger than blocksize %D",nselectedvar[i],blocksize[i]);
3190: }
3192: pstart = network->pStart;
3193: pend = network->pEnd;
3195: /* Get local number of idx */
3196: nidx = 0;
3197: for (p=pstart; p<pend; p++) {
3198: DMISAddSize_private(network,p,numkeys,keys,blocksize,nselectedvar,&nidx);
3199: }
3201: /* Compute local idx */
3202: PetscMalloc1(nidx,&idx);
3203: i = 0;
3204: for (p=pstart; p<pend; p++) {
3205: DMISComputeLocalIdx_private(dm,p,numkeys,keys,blocksize,nselectedvar,selectedvar,&i,idx);
3206: }
3208: /* Create is */
3209: ISCreateGeneral(PETSC_COMM_SELF,nidx,idx,PETSC_COPY_VALUES,is);
3210: PetscFree(idx);
3211: return(0);
3212: }