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: {
8: PetscOptionsHead(PetscOptionsObject,"DMNetwork Options");
9: PetscOptionsTail();
10: return 0;
11: }
13: /* External function declarations here */
14: extern PetscErrorCode DMCreateMatrix_Network(DM, Mat*);
15: extern PetscErrorCode DMDestroy_Network(DM);
16: extern PetscErrorCode DMView_Network(DM, PetscViewer);
17: extern PetscErrorCode DMGlobalToLocalBegin_Network(DM, Vec, InsertMode, Vec);
18: extern PetscErrorCode DMGlobalToLocalEnd_Network(DM, Vec, InsertMode, Vec);
19: extern PetscErrorCode DMLocalToGlobalBegin_Network(DM, Vec, InsertMode, Vec);
20: extern PetscErrorCode DMLocalToGlobalEnd_Network(DM, Vec, InsertMode, Vec);
21: extern PetscErrorCode DMSetUp_Network(DM);
22: extern PetscErrorCode DMClone_Network(DM, DM*);
24: static PetscErrorCode VecArrayPrint_private(PetscViewer viewer,PetscInt n,const PetscScalar *xv)
25: {
26: PetscInt i;
28: for (i=0; i<n; i++) {
29: #if defined(PETSC_USE_COMPLEX)
30: if (PetscImaginaryPart(xv[i]) > 0.0) {
31: PetscViewerASCIIPrintf(viewer," %g + %g i\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
32: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
33: PetscViewerASCIIPrintf(viewer," %g - %g i\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
34: } else {
35: PetscViewerASCIIPrintf(viewer," %g\n",(double)PetscRealPart(xv[i]));
36: }
37: #else
38: PetscViewerASCIIPrintf(viewer," %g\n",(double)xv[i]);
39: #endif
40: }
41: return 0;
42: }
44: static PetscErrorCode VecView_Network_Seq(DM networkdm,Vec X,PetscViewer viewer)
45: {
46: PetscInt e,v,Start,End,offset,nvar,id;
47: const PetscScalar *xv;
49: VecGetArrayRead(X,&xv);
51: /* iterate over edges */
52: DMNetworkGetEdgeRange(networkdm,&Start,&End);
53: for (e=Start; e<End; e++) {
54: DMNetworkGetComponent(networkdm,e,ALL_COMPONENTS,NULL,NULL,&nvar);
55: if (!nvar) continue;
57: DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&offset);
58: DMNetworkGetGlobalEdgeIndex(networkdm,e,&id);
60: PetscViewerASCIIPrintf(viewer," Edge %" PetscInt_FMT ":\n",id);
61: VecArrayPrint_private(viewer,nvar,xv+offset);
62: }
64: /* iterate over vertices */
65: DMNetworkGetVertexRange(networkdm,&Start,&End);
66: for (v=Start; v<End; v++) {
67: DMNetworkGetComponent(networkdm,v,ALL_COMPONENTS,NULL,NULL,&nvar);
68: if (!nvar) continue;
70: DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&offset);
71: DMNetworkGetGlobalVertexIndex(networkdm,v,&id);
73: PetscViewerASCIIPrintf(viewer," Vertex %" PetscInt_FMT ":\n",id);
74: VecArrayPrint_private(viewer,nvar,xv+offset);
75: }
76: PetscViewerFlush(viewer);
77: VecRestoreArrayRead(X,&xv);
78: return 0;
79: }
81: static PetscErrorCode VecView_Network_MPI(DM networkdm,Vec X,PetscViewer viewer)
82: {
83: PetscInt i,e,v,eStart,eEnd,vStart,vEnd,offset,nvar,len_loc,len,k;
84: const PetscScalar *xv;
85: MPI_Comm comm;
86: PetscMPIInt size,rank,tag = ((PetscObject)viewer)->tag;
87: Vec localX;
88: PetscBool ghostvtex;
89: PetscScalar *values;
90: PetscInt j,ne,nv,id;
91: MPI_Status status;
93: PetscObjectGetComm((PetscObject)networkdm,&comm);
94: MPI_Comm_size(comm,&size);
95: MPI_Comm_rank(comm,&rank);
97: DMGetLocalVector(networkdm,&localX);
98: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
99: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
100: VecGetArrayRead(localX,&xv);
102: VecGetLocalSize(localX,&len_loc);
104: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
105: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
106: len_loc += 2*(1 + eEnd-eStart + vEnd-vStart);
108: /* values = [nedges, nvertices; id, nvar, xedge; ...; id, nvars, xvertex;...], to be sent to proc[0] */
109: MPI_Allreduce(&len_loc,&len,1,MPIU_INT,MPI_MAX,comm);
110: PetscCalloc1(len,&values);
112: if (rank == 0) {
113: PetscViewerASCIIPrintf(viewer,"Process [%d]\n",rank);
114: }
116: /* iterate over edges */
117: k = 2;
118: for (e=eStart; e<eEnd; e++) {
119: DMNetworkGetComponent(networkdm,e,ALL_COMPONENTS,NULL,NULL,&nvar);
120: if (!nvar) continue;
122: DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&offset);
123: DMNetworkGetGlobalEdgeIndex(networkdm,e,&id);
125: if (rank == 0) { /* print its own entries */
126: PetscViewerASCIIPrintf(viewer," Edge %" PetscInt_FMT ":\n",id);
127: VecArrayPrint_private(viewer,nvar,xv+offset);
128: } else {
129: values[0] += 1; /* number of edges */
130: values[k++] = id;
131: values[k++] = nvar;
132: for (i=offset; i< offset+nvar; i++) values[k++] = xv[i];
133: }
134: }
136: /* iterate over vertices */
137: for (v=vStart; v<vEnd; v++) {
138: DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
139: if (ghostvtex) continue;
140: DMNetworkGetComponent(networkdm,v,ALL_COMPONENTS,NULL,NULL,&nvar);
141: if (!nvar) continue;
143: DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&offset);
144: DMNetworkGetGlobalVertexIndex(networkdm,v,&id);
146: if (rank == 0) {
147: PetscViewerASCIIPrintf(viewer," Vertex %" PetscInt_FMT ":\n",id);
148: VecArrayPrint_private(viewer,nvar,xv+offset);
149: } else {
150: values[1] += 1; /* number of vertices */
151: values[k++] = id;
152: values[k++] = nvar;
153: for (i=offset; i< offset+nvar; i++) values[k++] = xv[i];
154: }
155: }
157: if (rank == 0) {
158: /* proc[0] receives and prints messages */
159: for (j=1; j<size; j++) {
160: PetscViewerASCIIPrintf(viewer,"Process [%d]\n",j);
162: MPI_Recv(values,(PetscMPIInt)len,MPIU_SCALAR,j,tag,comm,&status);
164: ne = (PetscInt)PetscAbsScalar(values[0]);
165: nv = (PetscInt)PetscAbsScalar(values[1]);
167: /* print received edges */
168: k = 2;
169: for (i=0; i<ne; i++) {
170: id = (PetscInt)PetscAbsScalar(values[k++]);
171: nvar = (PetscInt)PetscAbsScalar(values[k++]);
172: PetscViewerASCIIPrintf(viewer," Edge %" PetscInt_FMT ":\n",id);
173: VecArrayPrint_private(viewer,nvar,values+k);
174: k += nvar;
175: }
177: /* print received vertices */
178: for (i=0; i<nv; i++) {
179: id = (PetscInt)PetscAbsScalar(values[k++]);
180: nvar = (PetscInt)PetscAbsScalar(values[k++]);
181: PetscViewerASCIIPrintf(viewer," Vertex %" PetscInt_FMT ":\n",id);
182: VecArrayPrint_private(viewer,nvar,values+k);
183: k += nvar;
184: }
185: }
186: } else {
187: /* sends values to proc[0] */
188: MPI_Send((void*)values,k,MPIU_SCALAR,0,tag,comm);
189: }
191: PetscFree(values);
192: VecRestoreArrayRead(localX,&xv);
193: DMRestoreLocalVector(networkdm,&localX);
194: return 0;
195: }
197: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec,PetscViewer);
199: PetscErrorCode VecView_Network(Vec v,PetscViewer viewer)
200: {
201: DM dm;
202: PetscBool isseq;
203: PetscBool iascii;
205: VecGetDM(v,&dm);
207: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
208: PetscObjectTypeCompare((PetscObject)v,VECSEQ,&isseq);
210: /* Use VecView_Network if the viewer is ASCII; use VecView_Seq/MPI for other viewer formats */
211: if (iascii) {
212: if (isseq) {
213: VecView_Network_Seq(dm,v,viewer);
214: } else {
215: VecView_Network_MPI(dm,v,viewer);
216: }
217: } else {
218: if (isseq) {
219: VecView_Seq(v,viewer);
220: } else {
221: VecView_MPI(v,viewer);
222: }
223: }
224: return 0;
225: }
227: static PetscErrorCode DMCreateGlobalVector_Network(DM dm,Vec *vec)
228: {
229: DM_Network *network = (DM_Network*) dm->data;
231: DMCreateGlobalVector(network->plex,vec);
232: VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Network);
233: VecSetDM(*vec,dm);
234: return 0;
235: }
237: static PetscErrorCode DMCreateLocalVector_Network(DM dm,Vec *vec)
238: {
239: DM_Network *network = (DM_Network*) dm->data;
241: DMCreateLocalVector(network->plex,vec);
242: VecSetDM(*vec,dm);
243: return 0;
244: }
246: PetscErrorCode DMInitialize_Network(DM dm)
247: {
248: DMSetDimension(dm,1);
249: dm->ops->view = DMView_Network;
250: dm->ops->setfromoptions = DMSetFromOptions_Network;
251: dm->ops->clone = DMClone_Network;
252: dm->ops->setup = DMSetUp_Network;
253: dm->ops->createglobalvector = DMCreateGlobalVector_Network;
254: dm->ops->createlocalvector = DMCreateLocalVector_Network;
255: dm->ops->getlocaltoglobalmapping = NULL;
256: dm->ops->createfieldis = NULL;
257: dm->ops->createcoordinatedm = NULL;
258: dm->ops->getcoloring = NULL;
259: dm->ops->creatematrix = DMCreateMatrix_Network;
260: dm->ops->createinterpolation = NULL;
261: dm->ops->createinjection = NULL;
262: dm->ops->refine = NULL;
263: dm->ops->coarsen = NULL;
264: dm->ops->refinehierarchy = NULL;
265: dm->ops->coarsenhierarchy = NULL;
266: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Network;
267: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Network;
268: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Network;
269: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Network;
270: dm->ops->destroy = DMDestroy_Network;
271: dm->ops->createsubdm = NULL;
272: dm->ops->locatepoints = NULL;
273: return 0;
274: }
276: PetscErrorCode DMClone_Network(DM dm, DM *newdm)
277: {
278: DM_Network *network = (DM_Network *) dm->data;
280: network->refct++;
281: (*newdm)->data = network;
282: PetscObjectChangeTypeName((PetscObject) *newdm, DMNETWORK);
283: DMInitialize_Network(*newdm);
284: return 0;
285: }
287: /*MC
288: DMNETWORK = "network" - A DM object that encapsulates an unstructured network. The implementation is based on the DM object
289: DMPlex that manages unstructured grids. Distributed networks use a non-overlapping partitioning of
290: the edges. In the local representation, Vecs contain all unknowns in the interior and shared boundary.
291: This is specified by a PetscSection object. Ownership in the global representation is determined by
292: ownership of the underlying DMPlex points. This is specified by another PetscSection object.
294: Level: intermediate
296: .seealso: DMType, DMNetworkCreate(), DMCreate(), DMSetType()
297: M*/
299: PETSC_EXTERN PetscErrorCode DMCreate_Network(DM dm)
300: {
301: DM_Network *network;
304: PetscNewLog(dm,&network);
305: dm->data = network;
307: network->refct = 1;
308: network->NVertices = 0;
309: network->NEdges = 0;
310: network->nVertices = 0;
311: network->nEdges = 0;
312: network->nsubnet = 0;
314: network->max_comps_registered = 20;
315: network->component = NULL;
316: network->header = NULL;
317: network->cvalue = NULL;
319: DMInitialize_Network(dm);
320: return 0;
321: }
323: /*@
324: DMNetworkCreate - Creates a DMNetwork object, which encapsulates an unstructured network.
326: Collective
328: Input Parameter:
329: . comm - The communicator for the DMNetwork object
331: Output Parameter:
332: . network - The DMNetwork object
334: Level: beginner
336: @*/
337: PetscErrorCode DMNetworkCreate(MPI_Comm comm, DM *network)
338: {
340: DMCreate(comm, network);
341: DMSetType(*network, DMNETWORK);
342: return 0;
343: }