Actual source code: networkcreate.c
1: #define PETSCDM_DLL
2: #include <petsc/private/dmnetworkimpl.h>
3: #include <petsc/private/vecimpl.h>
5: PetscErrorCode DMSetFromOptions_Network(PetscOptionItems *PetscOptionsObject,DM dm)
6: {
11: PetscOptionsHead(PetscOptionsObject,"DMNetwork Options");
12: PetscOptionsTail();
13: return(0);
14: }
16: /* External function declarations here */
17: extern PetscErrorCode DMCreateMatrix_Network(DM, Mat*);
18: extern PetscErrorCode DMDestroy_Network(DM);
19: extern PetscErrorCode DMView_Network(DM, PetscViewer);
20: extern PetscErrorCode DMGlobalToLocalBegin_Network(DM, Vec, InsertMode, Vec);
21: extern PetscErrorCode DMGlobalToLocalEnd_Network(DM, Vec, InsertMode, Vec);
22: extern PetscErrorCode DMLocalToGlobalBegin_Network(DM, Vec, InsertMode, Vec);
23: extern PetscErrorCode DMLocalToGlobalEnd_Network(DM, Vec, InsertMode, Vec);
24: extern PetscErrorCode DMSetUp_Network(DM);
25: extern PetscErrorCode DMClone_Network(DM, DM*);
27: static PetscErrorCode VecArrayPrint_private(PetscViewer viewer,PetscInt n,const PetscScalar *xv)
28: {
30: PetscInt i;
33: for (i=0; i<n; i++) {
34: #if defined(PETSC_USE_COMPLEX)
35: if (PetscImaginaryPart(xv[i]) > 0.0) {
36: PetscViewerASCIIPrintf(viewer," %g + %g i\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
37: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
38: PetscViewerASCIIPrintf(viewer," %g - %g i\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
39: } else {
40: PetscViewerASCIIPrintf(viewer," %g\n",(double)PetscRealPart(xv[i]));
41: }
42: #else
43: PetscViewerASCIIPrintf(viewer," %g\n",(double)xv[i]);
44: #endif
45: }
46: return(0);
47: }
49: static PetscErrorCode VecView_Network_Seq(DM networkdm,Vec X,PetscViewer viewer)
50: {
51: PetscErrorCode ierr;
52: PetscInt e,v,Start,End,offset,nvar,id;
53: const PetscScalar *xv;
56: VecGetArrayRead(X,&xv);
58: /* iterate over edges */
59: DMNetworkGetEdgeRange(networkdm,&Start,&End);
60: for (e=Start; e<End; e++) {
61: DMNetworkGetComponent(networkdm,e,ALL_COMPONENTS,NULL,NULL,&nvar);
62: if (!nvar) continue;
64: DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&offset);
65: DMNetworkGetGlobalEdgeIndex(networkdm,e,&id);
67: PetscViewerASCIIPrintf(viewer," Edge %D:\n",id);
68: VecArrayPrint_private(viewer,nvar,xv+offset);
69: }
71: /* iterate over vertices */
72: DMNetworkGetVertexRange(networkdm,&Start,&End);
73: for (v=Start; v<End; v++) {
74: DMNetworkGetComponent(networkdm,v,ALL_COMPONENTS,NULL,NULL,&nvar);
75: if (!nvar) continue;
77: DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&offset);
78: DMNetworkGetGlobalVertexIndex(networkdm,v,&id);
80: PetscViewerASCIIPrintf(viewer," Vertex %D:\n",id);
81: VecArrayPrint_private(viewer,nvar,xv+offset);
82: }
83: PetscViewerFlush(viewer);
84: VecRestoreArrayRead(X,&xv);
85: return(0);
86: }
88: static PetscErrorCode VecView_Network_MPI(DM networkdm,Vec X,PetscViewer viewer)
89: {
90: PetscErrorCode ierr;
91: PetscInt i,e,v,eStart,eEnd,vStart,vEnd,offset,nvar,len_loc,len,k;
92: const PetscScalar *xv;
93: MPI_Comm comm;
94: PetscMPIInt size,rank,tag = ((PetscObject)viewer)->tag;
95: Vec localX;
96: PetscBool ghostvtex;
97: PetscScalar *values;
98: PetscInt j,ne,nv,id;
99: MPI_Status status;
102: PetscObjectGetComm((PetscObject)networkdm,&comm);
103: MPI_Comm_size(comm,&size);
104: MPI_Comm_rank(comm,&rank);
106: DMGetLocalVector(networkdm,&localX);
107: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
108: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
109: VecGetArrayRead(localX,&xv);
111: VecGetLocalSize(localX,&len_loc);
113: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
114: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
115: len_loc += 2*(1 + eEnd-eStart + vEnd-vStart);
117: /* values = [nedges, nvertices; id, nvar, xedge; ...; id, nvars, xvertex;...], to be sent to proc[0] */
118: MPI_Allreduce(&len_loc,&len,1,MPIU_INT,MPI_MAX,comm);
119: PetscCalloc1(len,&values);
121: if (!rank) {
122: PetscViewerASCIIPrintf(viewer,"Process [%d]\n",rank);
123: }
125: /* iterate over edges */
126: k = 2;
127: for (e=eStart; e<eEnd; e++) {
128: DMNetworkGetComponent(networkdm,e,ALL_COMPONENTS,NULL,NULL,&nvar);
129: if (!nvar) continue;
131: DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&offset);
132: DMNetworkGetGlobalEdgeIndex(networkdm,e,&id);
134: if (!rank) { /* print its own entries */
135: PetscViewerASCIIPrintf(viewer," Edge %D:\n",id);
136: VecArrayPrint_private(viewer,nvar,xv+offset);
137: } else {
138: values[0] += 1; /* number of edges */
139: values[k++] = id;
140: values[k++] = nvar;
141: for (i=offset; i< offset+nvar; i++) values[k++] = xv[i];
142: }
143: }
145: /* iterate over vertices */
146: for (v=vStart; v<vEnd; v++) {
147: DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
148: if (ghostvtex) continue;
149: DMNetworkGetComponent(networkdm,v,ALL_COMPONENTS,NULL,NULL,&nvar);
150: if (!nvar) continue;
152: DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&offset);
153: DMNetworkGetGlobalVertexIndex(networkdm,v,&id);
155: if (!rank) {
156: PetscViewerASCIIPrintf(viewer," Vertex %D:\n",id);
157: VecArrayPrint_private(viewer,nvar,xv+offset);
158: } else {
159: values[1] += 1; /* number of vertices */
160: values[k++] = id;
161: values[k++] = nvar;
162: for (i=offset; i< offset+nvar; i++) values[k++] = xv[i];
163: }
164: }
166: if (!rank) {
167: /* proc[0] receives and prints messages */
168: for (j=1; j<size; j++) {
169: PetscViewerASCIIPrintf(viewer,"Process [%d]\n",j);
171: MPI_Recv(values,(PetscMPIInt)len,MPIU_SCALAR,j,tag,comm,&status);
173: ne = (PetscInt)PetscAbsScalar(values[0]);
174: nv = (PetscInt)PetscAbsScalar(values[1]);
176: /* print received edges */
177: k = 2;
178: for (i=0; i<ne; i++) {
179: id = (PetscInt)PetscAbsScalar(values[k++]);
180: nvar = (PetscInt)PetscAbsScalar(values[k++]);
181: PetscViewerASCIIPrintf(viewer," Edge %D:\n",id);
182: VecArrayPrint_private(viewer,nvar,values+k);
183: k += nvar;
184: }
186: /* print received vertices */
187: for (i=0; i<nv; i++) {
188: id = (PetscInt)PetscAbsScalar(values[k++]);
189: nvar = (PetscInt)PetscAbsScalar(values[k++]);
190: PetscViewerASCIIPrintf(viewer," Vertex %D:\n",id);
191: VecArrayPrint_private(viewer,nvar,values+k);
192: k += nvar;
193: }
194: }
195: } else {
196: /* sends values to proc[0] */
197: MPI_Send((void*)values,k,MPIU_SCALAR,0,tag,comm);
198: }
200: PetscFree(values);
201: VecRestoreArrayRead(localX,&xv);
202: DMRestoreLocalVector(networkdm,&localX);
203: return(0);
204: }
206: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec,PetscViewer);
208: PetscErrorCode VecView_Network(Vec v,PetscViewer viewer)
209: {
210: DM dm;
212: PetscBool isseq;
213: PetscBool iascii;
216: VecGetDM(v,&dm);
217: if (!dm) SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_WRONG,"Vector not generated from a DM");
218: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
219: PetscObjectTypeCompare((PetscObject)v,VECSEQ,&isseq);
221: /* Use VecView_Network if the viewer is ASCII; use VecView_Seq/MPI for other viewer formats */
222: if (iascii) {
223: if (isseq) {
224: VecView_Network_Seq(dm,v,viewer);
225: } else {
226: VecView_Network_MPI(dm,v,viewer);
227: }
228: } else {
229: if (isseq) {
230: VecView_Seq(v,viewer);
231: } else {
232: VecView_MPI(v,viewer);
233: }
234: }
235: return(0);
236: }
238: static PetscErrorCode DMCreateGlobalVector_Network(DM dm,Vec *vec)
239: {
241: DM_Network *network = (DM_Network*) dm->data;
244: DMCreateGlobalVector(network->plex,vec);
245: VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Network);
246: VecSetDM(*vec,dm);
247: return(0);
248: }
250: static PetscErrorCode DMCreateLocalVector_Network(DM dm,Vec *vec)
251: {
253: DM_Network *network = (DM_Network*) dm->data;
256: DMCreateLocalVector(network->plex,vec);
257: VecSetDM(*vec,dm);
258: return(0);
259: }
261: PetscErrorCode DMInitialize_Network(DM dm)
262: {
266: DMSetDimension(dm,1);
267: dm->ops->view = DMView_Network;
268: dm->ops->setfromoptions = DMSetFromOptions_Network;
269: dm->ops->clone = DMClone_Network;
270: dm->ops->setup = DMSetUp_Network;
271: dm->ops->createglobalvector = DMCreateGlobalVector_Network;
272: dm->ops->createlocalvector = DMCreateLocalVector_Network;
273: dm->ops->getlocaltoglobalmapping = NULL;
274: dm->ops->createfieldis = NULL;
275: dm->ops->createcoordinatedm = NULL;
276: dm->ops->getcoloring = NULL;
277: dm->ops->creatematrix = DMCreateMatrix_Network;
278: dm->ops->createinterpolation = NULL;
279: dm->ops->createinjection = NULL;
280: dm->ops->refine = NULL;
281: dm->ops->coarsen = NULL;
282: dm->ops->refinehierarchy = NULL;
283: dm->ops->coarsenhierarchy = NULL;
284: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Network;
285: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Network;
286: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Network;
287: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Network;
288: dm->ops->destroy = DMDestroy_Network;
289: dm->ops->createsubdm = NULL;
290: dm->ops->locatepoints = NULL;
291: return(0);
292: }
294: PetscErrorCode DMClone_Network(DM dm, DM *newdm)
295: {
296: DM_Network *network = (DM_Network *) dm->data;
300: network->refct++;
301: (*newdm)->data = network;
302: PetscObjectChangeTypeName((PetscObject) *newdm, DMNETWORK);
303: DMInitialize_Network(*newdm);
304: return(0);
305: }
307: /*MC
308: DMNETWORK = "network" - A DM object that encapsulates an unstructured network. The implementation is based on the DM object
309: DMPlex that manages unstructured grids. Distributed networks use a non-overlapping partitioning of
310: the edges. In the local representation, Vecs contain all unknowns in the interior and shared boundary.
311: This is specified by a PetscSection object. Ownership in the global representation is determined by
312: ownership of the underlying DMPlex points. This is specified by another PetscSection object.
314: Level: intermediate
316: .seealso: DMType, DMNetworkCreate(), DMCreate(), DMSetType()
317: M*/
319: PETSC_EXTERN PetscErrorCode DMCreate_Network(DM dm)
320: {
321: DM_Network *network;
326: PetscNewLog(dm,&network);
327: dm->data = network;
329: network->refct = 1;
330: network->NVertices = 0;
331: network->NEdges = 0;
332: network->nVertices = 0;
333: network->nEdges = 0;
334: network->nsubnet = 0;
337: DMInitialize_Network(dm);
338: return(0);
339: }
341: /*@
342: DMNetworkCreate - Creates a DMNetwork object, which encapsulates an unstructured network.
344: Collective
346: Input Parameter:
347: . comm - The communicator for the DMNetwork object
349: Output Parameter:
350: . network - The DMNetwork object
352: Level: beginner
354: @*/
355: PetscErrorCode DMNetworkCreate(MPI_Comm comm, DM *network)
356: {
361: DMCreate(comm, network);
362: DMSetType(*network, DMNETWORK);
363: return(0);
364: }