Actual source code: networkcreate.c
1: #include <petsc/private/dmnetworkimpl.h>
2: #include <petsc/private/vecimpl.h>
4: static PetscErrorCode DMSetFromOptions_Network(DM dm, PetscOptionItems PetscOptionsObject)
5: {
6: PetscFunctionBegin;
7: PetscOptionsHeadBegin(PetscOptionsObject, "DMNetwork Options");
8: PetscOptionsHeadEnd();
9: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
28: for (i = 0; i < n; i++) {
29: if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(xv[i]) > 0.0) {
30: PetscCall(PetscViewerASCIIPrintf(viewer, " %g + %g i\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i])));
31: } else if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(xv[i]) < 0.0) {
32: PetscCall(PetscViewerASCIIPrintf(viewer, " %g - %g i\n", (double)PetscRealPart(xv[i]), -(double)PetscImaginaryPart(xv[i])));
33: } else PetscCall(PetscViewerASCIIPrintf(viewer, " %g\n", (double)PetscRealPart(xv[i])));
34: }
35: PetscFunctionReturn(PETSC_SUCCESS);
36: }
38: static PetscErrorCode VecView_Network_Seq(DM networkdm, Vec X, PetscViewer viewer)
39: {
40: PetscInt e, v, Start, End, offset, nvar, id;
41: const PetscScalar *xv;
43: PetscFunctionBegin;
44: PetscCall(VecGetArrayRead(X, &xv));
46: /* iterate over edges */
47: PetscCall(DMNetworkGetEdgeRange(networkdm, &Start, &End));
48: for (e = Start; e < End; e++) {
49: PetscCall(DMNetworkGetComponent(networkdm, e, ALL_COMPONENTS, NULL, NULL, &nvar));
50: if (!nvar) continue;
52: PetscCall(DMNetworkGetLocalVecOffset(networkdm, e, ALL_COMPONENTS, &offset));
53: PetscCall(DMNetworkGetGlobalEdgeIndex(networkdm, e, &id));
55: PetscCall(PetscViewerASCIIPrintf(viewer, " Edge %" PetscInt_FMT ":\n", id));
56: PetscCall(VecArrayPrint_private(viewer, nvar, xv + offset));
57: }
59: /* iterate over vertices */
60: PetscCall(DMNetworkGetVertexRange(networkdm, &Start, &End));
61: for (v = Start; v < End; v++) {
62: PetscCall(DMNetworkGetComponent(networkdm, v, ALL_COMPONENTS, NULL, NULL, &nvar));
63: if (!nvar) continue;
65: PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, ALL_COMPONENTS, &offset));
66: PetscCall(DMNetworkGetGlobalVertexIndex(networkdm, v, &id));
68: PetscCall(PetscViewerASCIIPrintf(viewer, " Vertex %" PetscInt_FMT ":\n", id));
69: PetscCall(VecArrayPrint_private(viewer, nvar, xv + offset));
70: }
71: PetscCall(PetscViewerFlush(viewer));
72: PetscCall(VecRestoreArrayRead(X, &xv));
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: static PetscErrorCode VecView_Network_MPI(DM networkdm, Vec X, PetscViewer viewer)
77: {
78: PetscInt i, e, v, eStart, eEnd, vStart, vEnd, offset, nvar, len_loc, k;
79: const PetscScalar *xv;
80: MPI_Comm comm;
81: PetscMPIInt size, rank, tag = ((PetscObject)viewer)->tag, len;
82: Vec localX;
83: PetscBool ghostvtex;
84: PetscScalar *values;
85: PetscInt ne, nv, id;
86: MPI_Status status;
88: PetscFunctionBegin;
89: PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm));
90: PetscCallMPI(MPI_Comm_size(comm, &size));
91: PetscCallMPI(MPI_Comm_rank(comm, &rank));
93: PetscCall(DMGetLocalVector(networkdm, &localX));
94: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
95: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
96: PetscCall(VecGetArrayRead(localX, &xv));
98: PetscCall(VecGetLocalSize(localX, &len_loc));
100: PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
101: PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
102: len_loc += 2 * (1 + eEnd - eStart + vEnd - vStart);
104: /* values = [nedges, nvertices; id, nvar, xedge; ...; id, nvars, xvertex;...], to be sent to proc[0] */
105: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &len_loc, 1, MPIU_INT, MPI_MAX, comm));
106: PetscCall(PetscMPIIntCast(len_loc, &len));
107: PetscCall(PetscCalloc1(len, &values));
109: if (rank == 0) PetscCall(PetscViewerASCIIPrintf(viewer, "Process [%d]\n", rank));
111: /* iterate over edges */
112: k = 2;
113: for (e = eStart; e < eEnd; e++) {
114: PetscCall(DMNetworkGetComponent(networkdm, e, ALL_COMPONENTS, NULL, NULL, &nvar));
115: if (!nvar) continue;
117: PetscCall(DMNetworkGetLocalVecOffset(networkdm, e, ALL_COMPONENTS, &offset));
118: PetscCall(DMNetworkGetGlobalEdgeIndex(networkdm, e, &id));
120: if (rank == 0) { /* print its own entries */
121: PetscCall(PetscViewerASCIIPrintf(viewer, " Edge %" PetscInt_FMT ":\n", id));
122: PetscCall(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: PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghostvtex));
134: if (ghostvtex) continue;
135: PetscCall(DMNetworkGetComponent(networkdm, v, ALL_COMPONENTS, NULL, NULL, &nvar));
136: if (!nvar) continue;
138: PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, ALL_COMPONENTS, &offset));
139: PetscCall(DMNetworkGetGlobalVertexIndex(networkdm, v, &id));
141: if (rank == 0) {
142: PetscCall(PetscViewerASCIIPrintf(viewer, " Vertex %" PetscInt_FMT ":\n", id));
143: PetscCall(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 (PetscMPIInt j = 1; j < size; j++) {
155: PetscCall(PetscViewerASCIIPrintf(viewer, "Process [%d]\n", j));
157: PetscCallMPI(MPI_Recv(values, 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: PetscCall(PetscViewerASCIIPrintf(viewer, " Edge %" PetscInt_FMT ":\n", id));
168: PetscCall(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: PetscCall(PetscViewerASCIIPrintf(viewer, " Vertex %" PetscInt_FMT ":\n", id));
177: PetscCall(VecArrayPrint_private(viewer, nvar, values + k));
178: k += nvar;
179: }
180: }
181: } else {
182: /* sends values to proc[0] */
183: PetscCallMPI(MPIU_Send((void *)values, k, MPIU_SCALAR, 0, tag, comm));
184: }
186: PetscCall(PetscFree(values));
187: PetscCall(VecRestoreArrayRead(localX, &xv));
188: PetscCall(DMRestoreLocalVector(networkdm, &localX));
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: PETSC_SINGLE_LIBRARY_INTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
194: static PetscErrorCode VecView_Network(Vec v, PetscViewer viewer)
195: {
196: DM dm;
197: PetscBool isseq;
198: PetscBool isascii;
200: PetscFunctionBegin;
201: PetscCall(VecGetDM(v, &dm));
202: PetscCheck(dm, PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
203: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
204: PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq));
206: /* Use VecView_Network if the viewer is ASCII; use VecView_Seq/MPI for other viewer formats */
207: if (isascii) {
208: if (isseq) PetscCall(VecView_Network_Seq(dm, v, viewer));
209: else PetscCall(VecView_Network_MPI(dm, v, viewer));
210: } else {
211: if (isseq) PetscCall(VecView_Seq(v, viewer));
212: else PetscCall(VecView_MPI(v, viewer));
213: }
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: static PetscErrorCode DMCreateGlobalVector_Network(DM dm, Vec *vec)
218: {
219: DM_Network *network = (DM_Network *)dm->data;
221: PetscFunctionBegin;
222: PetscCall(DMCreateGlobalVector(network->plex, vec));
223: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (PetscErrorCodeFn *)VecView_Network));
224: PetscCall(VecSetDM(*vec, dm));
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: static PetscErrorCode DMCreateLocalVector_Network(DM dm, Vec *vec)
229: {
230: DM_Network *network = (DM_Network *)dm->data;
232: PetscFunctionBegin;
233: PetscCall(DMCreateLocalVector(network->plex, vec));
234: PetscCall(VecSetDM(*vec, dm));
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: PetscErrorCode DMNetworkInitializeToDefault_NonShared(DM dm)
239: {
240: DM_Network *network = (DM_Network *)dm->data;
242: PetscFunctionBegin;
243: network->Je = NULL;
244: network->Jv = NULL;
245: network->Jvptr = NULL;
246: network->userEdgeJacobian = PETSC_FALSE;
247: network->userVertexJacobian = PETSC_FALSE;
249: network->vertex.DofSection = NULL;
250: network->vertex.GlobalDofSection = NULL;
251: network->vertex.mapping = NULL;
252: network->vertex.sf = NULL;
254: network->edge.DofSection = NULL;
255: network->edge.GlobalDofSection = NULL;
256: network->edge.mapping = NULL;
257: network->edge.sf = NULL;
259: network->DataSection = NULL;
260: network->DofSection = NULL;
261: network->GlobalDofSection = NULL;
262: network->componentsetup = PETSC_FALSE;
264: network->plex = NULL;
266: network->component = NULL;
267: network->ncomponent = 0;
269: network->header = NULL;
270: network->cvalue = NULL;
271: network->componentdataarray = NULL;
273: network->max_comps_registered = DMNETWORK_MAX_COMP_REGISTERED_DEFAULT; /* return to default */
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
276: /* Default values for the parameters in DMNetwork */
277: PetscErrorCode DMNetworkInitializeToDefault(DM dm)
278: {
279: DM_Network *network = (DM_Network *)dm->data;
280: DMNetworkCloneShared cloneshared = network->cloneshared;
282: PetscFunctionBegin;
283: PetscCall(DMNetworkInitializeToDefault_NonShared(dm));
284: /* Default values for shared data */
285: cloneshared->refct = 1;
286: cloneshared->NVertices = 0;
287: cloneshared->NEdges = 0;
288: cloneshared->nVertices = 0;
289: cloneshared->nEdges = 0;
290: cloneshared->nsubnet = 0;
291: cloneshared->pStart = -1;
292: cloneshared->pEnd = -1;
293: cloneshared->vStart = -1;
294: cloneshared->vEnd = -1;
295: cloneshared->eStart = -1;
296: cloneshared->eEnd = -1;
297: cloneshared->vltog = NULL;
298: cloneshared->distributecalled = PETSC_FALSE;
300: cloneshared->subnet = NULL;
301: cloneshared->subnetvtx = NULL;
302: cloneshared->subnetedge = NULL;
303: cloneshared->svtx = NULL;
304: cloneshared->nsvtx = 0;
305: cloneshared->Nsvtx = 0;
306: cloneshared->svertices = NULL;
307: cloneshared->sedgelist = NULL;
308: cloneshared->svtable = NULL;
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: static PetscErrorCode DMInitialize_Network(DM dm)
313: {
314: PetscFunctionBegin;
315: PetscCall(DMSetDimension(dm, 1));
316: dm->ops->view = DMView_Network;
317: dm->ops->setfromoptions = DMSetFromOptions_Network;
318: dm->ops->clone = DMClone_Network;
319: dm->ops->setup = DMSetUp_Network;
320: dm->ops->createglobalvector = DMCreateGlobalVector_Network;
321: dm->ops->createlocalvector = DMCreateLocalVector_Network;
322: dm->ops->getlocaltoglobalmapping = NULL;
323: dm->ops->createfieldis = NULL;
324: dm->ops->createcoordinatedm = DMCreateCoordinateDM_Network;
325: dm->ops->createcellcoordinatedm = NULL;
326: dm->ops->getcoloring = NULL;
327: dm->ops->creatematrix = DMCreateMatrix_Network;
328: dm->ops->createinterpolation = NULL;
329: dm->ops->createinjection = NULL;
330: dm->ops->refine = NULL;
331: dm->ops->coarsen = NULL;
332: dm->ops->refinehierarchy = NULL;
333: dm->ops->coarsenhierarchy = NULL;
334: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Network;
335: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Network;
336: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Network;
337: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Network;
338: dm->ops->destroy = DMDestroy_Network;
339: dm->ops->createsubdm = NULL;
340: dm->ops->locatepoints = NULL;
341: PetscFunctionReturn(PETSC_SUCCESS);
342: }
343: /*
344: copies over the subnetid and index portions of the DMNetworkComponentHeader from original dm to the newdm
345: */
346: static PetscErrorCode DMNetworkCopyHeaderTopological(DM dm, DM newdm)
347: {
348: DM_Network *network = (DM_Network *)dm->data, *newnetwork = (DM_Network *)newdm->data;
349: PetscInt p, i, np, index, subnetid;
351: PetscFunctionBegin;
352: np = network->cloneshared->pEnd - network->cloneshared->pStart;
353: PetscCall(PetscCalloc2(np, &newnetwork->header, np, &newnetwork->cvalue));
354: for (i = 0; i < np; i++) {
355: p = i + network->cloneshared->pStart;
356: PetscCall(DMNetworkGetSubnetID(dm, p, &subnetid));
357: PetscCall(DMNetworkGetIndex(dm, p, &index));
358: newnetwork->header[i].index = index;
359: newnetwork->header[i].subnetid = subnetid;
360: newnetwork->header[i].size = NULL;
361: newnetwork->header[i].key = NULL;
362: newnetwork->header[i].offset = NULL;
363: newnetwork->header[i].nvar = NULL;
364: newnetwork->header[i].offsetvarrel = NULL;
365: newnetwork->header[i].ndata = 0;
366: newnetwork->header[i].maxcomps = DMNETWORK_MAX_COMP_AT_POINT_DEFAULT;
367: newnetwork->header[i].hsize = sizeof(struct _p_DMNetworkComponentHeader) / sizeof(sizeof(DMNetworkComponentGenericDataType));
368: }
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: PetscErrorCode DMClone_Network(DM dm, DM *newdm)
373: {
374: DM_Network *network = (DM_Network *)dm->data, *newnetwork = NULL;
376: PetscFunctionBegin;
377: network->cloneshared->refct++;
378: PetscCall(PetscNew(&newnetwork));
379: (*newdm)->data = newnetwork;
380: PetscCall(DMNetworkInitializeToDefault_NonShared(*newdm));
381: newnetwork->cloneshared = network->cloneshared; /* Share all data that can be cloneshared */
383: PetscCheck(network->plex, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_NULL, "Must call DMNetworkLayoutSetUp() first");
384: PetscCall(DMClone(network->plex, &newnetwork->plex));
385: PetscCall(DMNetworkCopyHeaderTopological(dm, *newdm));
386: PetscCall(DMNetworkInitializeNonTopological(*newdm)); /* initialize all non-topological data to the state after DMNetworkLayoutSetUp as been called */
387: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMNETWORK));
388: PetscCall(DMInitialize_Network(*newdm));
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: /* Developer Note: Be aware that the plex inside of the network does not have a coordinate plex.
393: */
394: PetscErrorCode DMCreateCoordinateDM_Network(DM dm, DM *cdm)
395: {
396: DM_Network *newnetwork = NULL;
397: PetscInt Nf;
398: const char *prefix;
400: PetscFunctionBegin;
401: PetscCall(DMClone(dm, cdm));
402: newnetwork = (DM_Network *)(*cdm)->data;
403: PetscCall(DMGetNumFields(newnetwork->plex, &Nf));
404: PetscCall(DMSetNumFields(*cdm, Nf)); /* consistency with the coordinate plex */
405: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
406: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*cdm, prefix));
407: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)*cdm, "cdm_"));
408: PetscFunctionReturn(PETSC_SUCCESS);
409: }
411: /*MC
412: DMNETWORK = "network" - A DM object that encapsulates an unstructured network. The implementation is based on the DM object
413: DMPlex that manages unstructured grids. Distributed networks use a non-overlapping partitioning of
414: the edges. In the local representation, Vecs contain all unknowns in the interior and shared boundary.
415: This is specified by a PetscSection object. Ownership in the global representation is determined by
416: ownership of the underlying DMPlex points. This is specified by another PetscSection object.
418: Level: intermediate
420: .seealso: `DMType`, `DMNetworkCreate()`, `DMCreate()`, `DMSetType()`
421: M*/
423: PETSC_EXTERN PetscErrorCode DMCreate_Network(DM dm)
424: {
425: DM_Network *network;
427: PetscFunctionBegin;
429: PetscCall(PetscNew(&network));
430: PetscCall(PetscNew(&network->cloneshared));
431: dm->data = network;
433: PetscCall(DMNetworkInitializeToDefault(dm));
434: PetscCall(DMInitialize_Network(dm));
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: /*@
439: DMNetworkCreate - Creates a DMNetwork object, which encapsulates an unstructured network.
441: Collective
443: Input Parameter:
444: . comm - The communicator for the DMNetwork object
446: Output Parameter:
447: . network - The DMNetwork object
449: Level: beginner
451: .seealso: `DMCreate()`
452: @*/
453: PetscErrorCode DMNetworkCreate(MPI_Comm comm, DM *network)
454: {
455: PetscFunctionBegin;
456: PetscAssertPointer(network, 2);
457: PetscCall(DMCreate(comm, network));
458: PetscCall(DMSetType(*network, DMNETWORK));
459: PetscFunctionReturn(PETSC_SUCCESS);
460: }