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(DM dm, PetscOptionItems *PetscOptionsObject)
6: {
7: PetscOptionsHeadBegin(PetscOptionsObject, "DMNetwork Options");
8: PetscOptionsHeadEnd();
9: return 0;
10: }
12: /* External function declarations here */
13: extern PetscErrorCode DMCreateMatrix_Network(DM, Mat *);
14: extern PetscErrorCode DMDestroy_Network(DM);
15: extern PetscErrorCode DMView_Network(DM, PetscViewer);
16: extern PetscErrorCode DMGlobalToLocalBegin_Network(DM, Vec, InsertMode, Vec);
17: extern PetscErrorCode DMGlobalToLocalEnd_Network(DM, Vec, InsertMode, Vec);
18: extern PetscErrorCode DMLocalToGlobalBegin_Network(DM, Vec, InsertMode, Vec);
19: extern PetscErrorCode DMLocalToGlobalEnd_Network(DM, Vec, InsertMode, Vec);
20: extern PetscErrorCode DMSetUp_Network(DM);
21: extern PetscErrorCode DMClone_Network(DM, DM *);
23: static PetscErrorCode VecArrayPrint_private(PetscViewer viewer, PetscInt n, const PetscScalar *xv)
24: {
25: PetscInt i;
27: for (i = 0; i < n; i++) {
28: #if defined(PETSC_USE_COMPLEX)
29: if (PetscImaginaryPart(xv[i]) > 0.0) {
30: PetscViewerASCIIPrintf(viewer, " %g + %g i\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i]));
31: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
32: PetscViewerASCIIPrintf(viewer, " %g - %g i\n", (double)PetscRealPart(xv[i]), -(double)PetscImaginaryPart(xv[i]));
33: } else PetscViewerASCIIPrintf(viewer, " %g\n", (double)PetscRealPart(xv[i]));
34: #else
35: PetscViewerASCIIPrintf(viewer, " %g\n", (double)xv[i]);
36: #endif
37: }
38: return 0;
39: }
41: static PetscErrorCode VecView_Network_Seq(DM networkdm, Vec X, PetscViewer viewer)
42: {
43: PetscInt e, v, Start, End, offset, nvar, id;
44: const PetscScalar *xv;
46: VecGetArrayRead(X, &xv);
48: /* iterate over edges */
49: DMNetworkGetEdgeRange(networkdm, &Start, &End);
50: for (e = Start; e < End; e++) {
51: DMNetworkGetComponent(networkdm, e, ALL_COMPONENTS, NULL, NULL, &nvar);
52: if (!nvar) continue;
54: DMNetworkGetLocalVecOffset(networkdm, e, ALL_COMPONENTS, &offset);
55: DMNetworkGetGlobalEdgeIndex(networkdm, e, &id);
57: PetscViewerASCIIPrintf(viewer, " Edge %" PetscInt_FMT ":\n", id);
58: VecArrayPrint_private(viewer, nvar, xv + offset);
59: }
61: /* iterate over vertices */
62: DMNetworkGetVertexRange(networkdm, &Start, &End);
63: for (v = Start; v < End; v++) {
64: DMNetworkGetComponent(networkdm, v, ALL_COMPONENTS, NULL, NULL, &nvar);
65: if (!nvar) continue;
67: DMNetworkGetLocalVecOffset(networkdm, v, ALL_COMPONENTS, &offset);
68: DMNetworkGetGlobalVertexIndex(networkdm, v, &id);
70: PetscViewerASCIIPrintf(viewer, " Vertex %" PetscInt_FMT ":\n", id);
71: VecArrayPrint_private(viewer, nvar, xv + offset);
72: }
73: PetscViewerFlush(viewer);
74: VecRestoreArrayRead(X, &xv);
75: return 0;
76: }
78: static PetscErrorCode VecView_Network_MPI(DM networkdm, Vec X, PetscViewer viewer)
79: {
80: PetscInt i, e, v, eStart, eEnd, vStart, vEnd, offset, nvar, len_loc, len, k;
81: const PetscScalar *xv;
82: MPI_Comm comm;
83: PetscMPIInt size, rank, tag = ((PetscObject)viewer)->tag;
84: Vec localX;
85: PetscBool ghostvtex;
86: PetscScalar *values;
87: PetscInt j, ne, nv, id;
88: MPI_Status status;
90: PetscObjectGetComm((PetscObject)networkdm, &comm);
91: MPI_Comm_size(comm, &size);
92: MPI_Comm_rank(comm, &rank);
94: DMGetLocalVector(networkdm, &localX);
95: DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX);
96: DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX);
97: VecGetArrayRead(localX, &xv);
99: VecGetLocalSize(localX, &len_loc);
101: DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd);
102: DMNetworkGetVertexRange(networkdm, &vStart, &vEnd);
103: len_loc += 2 * (1 + eEnd - eStart + vEnd - vStart);
105: /* values = [nedges, nvertices; id, nvar, xedge; ...; id, nvars, xvertex;...], to be sent to proc[0] */
106: MPI_Allreduce(&len_loc, &len, 1, MPIU_INT, MPI_MAX, comm);
107: PetscCalloc1(len, &values);
109: if (rank == 0) PetscViewerASCIIPrintf(viewer, "Process [%d]\n", rank);
111: /* iterate over edges */
112: k = 2;
113: for (e = eStart; e < eEnd; e++) {
114: DMNetworkGetComponent(networkdm, e, ALL_COMPONENTS, NULL, NULL, &nvar);
115: if (!nvar) continue;
117: DMNetworkGetLocalVecOffset(networkdm, e, ALL_COMPONENTS, &offset);
118: DMNetworkGetGlobalEdgeIndex(networkdm, e, &id);
120: if (rank == 0) { /* print its own entries */
121: PetscViewerASCIIPrintf(viewer, " Edge %" PetscInt_FMT ":\n", id);
122: VecArrayPrint_private(viewer, nvar, xv + offset);
123: } else {
124: values[0] += 1; /* number of edges */
125: values[k++] = id;
126: values[k++] = nvar;
127: for (i = offset; i < offset + nvar; i++) values[k++] = xv[i];
128: }
129: }
131: /* iterate over vertices */
132: for (v = vStart; v < vEnd; v++) {
133: DMNetworkIsGhostVertex(networkdm, v, &ghostvtex);
134: if (ghostvtex) continue;
135: DMNetworkGetComponent(networkdm, v, ALL_COMPONENTS, NULL, NULL, &nvar);
136: if (!nvar) continue;
138: DMNetworkGetLocalVecOffset(networkdm, v, ALL_COMPONENTS, &offset);
139: DMNetworkGetGlobalVertexIndex(networkdm, v, &id);
141: if (rank == 0) {
142: PetscViewerASCIIPrintf(viewer, " Vertex %" PetscInt_FMT ":\n", id);
143: VecArrayPrint_private(viewer, nvar, xv + offset);
144: } else {
145: values[1] += 1; /* number of vertices */
146: values[k++] = id;
147: values[k++] = nvar;
148: for (i = offset; i < offset + nvar; i++) values[k++] = xv[i];
149: }
150: }
152: if (rank == 0) {
153: /* proc[0] receives and prints messages */
154: for (j = 1; j < size; j++) {
155: PetscViewerASCIIPrintf(viewer, "Process [%" PetscInt_FMT "]\n", j);
157: MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, comm, &status);
159: ne = (PetscInt)PetscAbsScalar(values[0]);
160: nv = (PetscInt)PetscAbsScalar(values[1]);
162: /* print received edges */
163: k = 2;
164: for (i = 0; i < ne; i++) {
165: id = (PetscInt)PetscAbsScalar(values[k++]);
166: nvar = (PetscInt)PetscAbsScalar(values[k++]);
167: PetscViewerASCIIPrintf(viewer, " Edge %" PetscInt_FMT ":\n", id);
168: VecArrayPrint_private(viewer, nvar, values + k);
169: k += nvar;
170: }
172: /* print received vertices */
173: for (i = 0; i < nv; i++) {
174: id = (PetscInt)PetscAbsScalar(values[k++]);
175: nvar = (PetscInt)PetscAbsScalar(values[k++]);
176: PetscViewerASCIIPrintf(viewer, " Vertex %" PetscInt_FMT ":\n", id);
177: VecArrayPrint_private(viewer, nvar, values + k);
178: k += nvar;
179: }
180: }
181: } else {
182: /* sends values to proc[0] */
183: MPI_Send((void *)values, k, MPIU_SCALAR, 0, tag, comm);
184: }
186: PetscFree(values);
187: VecRestoreArrayRead(localX, &xv);
188: DMRestoreLocalVector(networkdm, &localX);
189: return 0;
190: }
192: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
194: PetscErrorCode VecView_Network(Vec v, PetscViewer viewer)
195: {
196: DM dm;
197: PetscBool isseq;
198: PetscBool iascii;
200: VecGetDM(v, &dm);
202: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
203: PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq);
205: /* Use VecView_Network if the viewer is ASCII; use VecView_Seq/MPI for other viewer formats */
206: if (iascii) {
207: if (isseq) VecView_Network_Seq(dm, v, viewer);
208: else VecView_Network_MPI(dm, v, viewer);
209: } else {
210: if (isseq) VecView_Seq(v, viewer);
211: else VecView_MPI(v, viewer);
212: }
213: return 0;
214: }
216: static PetscErrorCode DMCreateGlobalVector_Network(DM dm, Vec *vec)
217: {
218: DM_Network *network = (DM_Network *)dm->data;
220: DMCreateGlobalVector(network->plex, vec);
221: VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Network);
222: VecSetDM(*vec, dm);
223: return 0;
224: }
226: static PetscErrorCode DMCreateLocalVector_Network(DM dm, Vec *vec)
227: {
228: DM_Network *network = (DM_Network *)dm->data;
230: DMCreateLocalVector(network->plex, vec);
231: VecSetDM(*vec, dm);
232: return 0;
233: }
235: PetscErrorCode DMNetworkInitializeToDefault_NonShared(DM dm)
236: {
237: DM_Network *network = (DM_Network *)dm->data;
239: network->Je = NULL;
240: network->Jv = NULL;
241: network->Jvptr = NULL;
242: network->userEdgeJacobian = PETSC_FALSE;
243: network->userVertexJacobian = PETSC_FALSE;
245: network->vertex.DofSection = NULL;
246: network->vertex.GlobalDofSection = NULL;
247: network->vertex.mapping = NULL;
248: network->vertex.sf = NULL;
250: network->edge.DofSection = NULL;
251: network->edge.GlobalDofSection = NULL;
252: network->edge.mapping = NULL;
253: network->edge.sf = NULL;
255: network->DataSection = NULL;
256: network->DofSection = NULL;
257: network->GlobalDofSection = NULL;
258: network->componentsetup = PETSC_FALSE;
260: network->plex = NULL;
262: network->component = NULL;
263: network->ncomponent = 0;
265: network->header = NULL;
266: network->cvalue = NULL;
267: network->componentdataarray = NULL;
269: network->max_comps_registered = DMNETWORK_MAX_COMP_REGISTERED_DEFAULT; /* return to default */
270: return 0;
271: }
272: /* Default values for the parameters in DMNetwork */
273: PetscErrorCode DMNetworkInitializeToDefault(DM dm)
274: {
275: DM_Network *network = (DM_Network *)dm->data;
276: DMNetworkCloneShared cloneshared = network->cloneshared;
278: DMNetworkInitializeToDefault_NonShared(dm);
279: /* Default values for shared data */
280: cloneshared->refct = 1;
281: cloneshared->NVertices = 0;
282: cloneshared->NEdges = 0;
283: cloneshared->nVertices = 0;
284: cloneshared->nEdges = 0;
285: cloneshared->nsubnet = 0;
286: cloneshared->pStart = -1;
287: cloneshared->pEnd = -1;
288: cloneshared->vStart = -1;
289: cloneshared->vEnd = -1;
290: cloneshared->eStart = -1;
291: cloneshared->eEnd = -1;
292: cloneshared->vltog = NULL;
293: cloneshared->distributecalled = PETSC_FALSE;
295: cloneshared->subnet = NULL;
296: cloneshared->subnetvtx = NULL;
297: cloneshared->subnetedge = NULL;
298: cloneshared->svtx = NULL;
299: cloneshared->nsvtx = 0;
300: cloneshared->Nsvtx = 0;
301: cloneshared->svertices = NULL;
302: cloneshared->sedgelist = NULL;
303: cloneshared->svtable = NULL;
304: return 0;
305: }
307: PetscErrorCode DMInitialize_Network(DM dm)
308: {
309: DMSetDimension(dm, 1);
310: dm->ops->view = DMView_Network;
311: dm->ops->setfromoptions = DMSetFromOptions_Network;
312: dm->ops->clone = DMClone_Network;
313: dm->ops->setup = DMSetUp_Network;
314: dm->ops->createglobalvector = DMCreateGlobalVector_Network;
315: dm->ops->createlocalvector = DMCreateLocalVector_Network;
316: dm->ops->getlocaltoglobalmapping = NULL;
317: dm->ops->createfieldis = NULL;
318: dm->ops->createcoordinatedm = NULL;
319: dm->ops->getcoloring = NULL;
320: dm->ops->creatematrix = DMCreateMatrix_Network;
321: dm->ops->createinterpolation = NULL;
322: dm->ops->createinjection = NULL;
323: dm->ops->refine = NULL;
324: dm->ops->coarsen = NULL;
325: dm->ops->refinehierarchy = NULL;
326: dm->ops->coarsenhierarchy = NULL;
327: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Network;
328: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Network;
329: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Network;
330: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Network;
331: dm->ops->destroy = DMDestroy_Network;
332: dm->ops->createsubdm = NULL;
333: dm->ops->locatepoints = NULL;
334: return 0;
335: }
336: /*
337: copies over the subnetid and index portions of the DMNetworkComponentHeader from orignal dm to the newdm
338: */
339: static PetscErrorCode DMNetworkCopyHeaderTopological(DM dm, DM newdm)
340: {
341: DM_Network *network = (DM_Network *)dm->data, *newnetwork = (DM_Network *)newdm->data;
342: PetscInt p, i, np, index, subnetid;
344: np = network->cloneshared->pEnd - network->cloneshared->pStart;
345: PetscCalloc2(np, &newnetwork->header, np, &newnetwork->cvalue);
346: for (i = 0; i < np; i++) {
347: p = i + network->cloneshared->pStart;
348: DMNetworkGetSubnetID(dm, p, &subnetid);
349: DMNetworkGetIndex(dm, p, &index);
350: newnetwork->header[i].index = index;
351: newnetwork->header[i].subnetid = subnetid;
352: newnetwork->header[i].size = NULL;
353: newnetwork->header[i].key = NULL;
354: newnetwork->header[i].offset = NULL;
355: newnetwork->header[i].nvar = NULL;
356: newnetwork->header[i].offsetvarrel = NULL;
357: newnetwork->header[i].ndata = 0;
358: newnetwork->header[i].maxcomps = DMNETWORK_MAX_COMP_AT_POINT_DEFAULT;
359: newnetwork->header[i].hsize = sizeof(struct _p_DMNetworkComponentHeader) / sizeof(sizeof(DMNetworkComponentGenericDataType));
360: }
361: return 0;
362: }
364: PetscErrorCode DMClone_Network(DM dm, DM *newdm)
365: {
366: DM_Network *network = (DM_Network *)dm->data, *newnetwork = NULL;
368: network->cloneshared->refct++;
369: PetscNew(&newnetwork);
370: (*newdm)->data = newnetwork;
371: DMNetworkInitializeToDefault_NonShared(*newdm);
372: newnetwork->cloneshared = network->cloneshared; /* Share all data that can be cloneshared */
373: DMClone(network->plex, &newnetwork->plex);
374: DMNetworkCopyHeaderTopological(dm, *newdm);
375: DMNetworkInitializeNonTopological(*newdm); /* initialize all non-topological data to the state after DMNetworkLayoutSetUp as been called */
376: PetscObjectChangeTypeName((PetscObject)*newdm, DMNETWORK);
377: DMInitialize_Network(*newdm);
378: return 0;
379: }
381: /*MC
382: DMNETWORK = "network" - A DM object that encapsulates an unstructured network. The implementation is based on the DM object
383: DMPlex that manages unstructured grids. Distributed networks use a non-overlapping partitioning of
384: the edges. In the local representation, Vecs contain all unknowns in the interior and shared boundary.
385: This is specified by a PetscSection object. Ownership in the global representation is determined by
386: ownership of the underlying DMPlex points. This is specified by another PetscSection object.
388: Level: intermediate
390: .seealso: `DMType`, `DMNetworkCreate()`, `DMCreate()`, `DMSetType()`
391: M*/
393: PETSC_EXTERN PetscErrorCode DMCreate_Network(DM dm)
394: {
395: DM_Network *network;
398: PetscNew(&network);
399: PetscNew(&network->cloneshared);
400: dm->data = network;
402: DMNetworkInitializeToDefault(dm);
403: DMInitialize_Network(dm);
404: return 0;
405: }
407: /*@
408: DMNetworkCreate - Creates a DMNetwork object, which encapsulates an unstructured network.
410: Collective
412: Input Parameter:
413: . comm - The communicator for the DMNetwork object
415: Output Parameter:
416: . network - The DMNetwork object
418: Level: beginner
420: @*/
421: PetscErrorCode DMNetworkCreate(MPI_Comm comm, DM *network)
422: {
424: DMCreate(comm, network);
425: DMSetType(*network, DMNETWORK);
426: return 0;
427: }