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 defined(PETSC_USE_COMPLEX)
30: if (PetscImaginaryPart(xv[i]) > 0.0) {
31: PetscCall(PetscViewerASCIIPrintf(viewer, " %g + %g i\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i])));
32: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
33: PetscCall(PetscViewerASCIIPrintf(viewer, " %g - %g i\n", (double)PetscRealPart(xv[i]), -(double)PetscImaginaryPart(xv[i])));
34: } else PetscCall(PetscViewerASCIIPrintf(viewer, " %g\n", (double)PetscRealPart(xv[i])));
35: #else
36: PetscCall(PetscViewerASCIIPrintf(viewer, " %g\n", (double)xv[i]));
37: #endif
38: }
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: static PetscErrorCode VecView_Network_Seq(DM networkdm, Vec X, PetscViewer viewer)
43: {
44: PetscInt e, v, Start, End, offset, nvar, id;
45: const PetscScalar *xv;
47: PetscFunctionBegin;
48: PetscCall(VecGetArrayRead(X, &xv));
50: /* iterate over edges */
51: PetscCall(DMNetworkGetEdgeRange(networkdm, &Start, &End));
52: for (e = Start; e < End; e++) {
53: PetscCall(DMNetworkGetComponent(networkdm, e, ALL_COMPONENTS, NULL, NULL, &nvar));
54: if (!nvar) continue;
56: PetscCall(DMNetworkGetLocalVecOffset(networkdm, e, ALL_COMPONENTS, &offset));
57: PetscCall(DMNetworkGetGlobalEdgeIndex(networkdm, e, &id));
59: PetscCall(PetscViewerASCIIPrintf(viewer, " Edge %" PetscInt_FMT ":\n", id));
60: PetscCall(VecArrayPrint_private(viewer, nvar, xv + offset));
61: }
63: /* iterate over vertices */
64: PetscCall(DMNetworkGetVertexRange(networkdm, &Start, &End));
65: for (v = Start; v < End; v++) {
66: PetscCall(DMNetworkGetComponent(networkdm, v, ALL_COMPONENTS, NULL, NULL, &nvar));
67: if (!nvar) continue;
69: PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, ALL_COMPONENTS, &offset));
70: PetscCall(DMNetworkGetGlobalVertexIndex(networkdm, v, &id));
72: PetscCall(PetscViewerASCIIPrintf(viewer, " Vertex %" PetscInt_FMT ":\n", id));
73: PetscCall(VecArrayPrint_private(viewer, nvar, xv + offset));
74: }
75: PetscCall(PetscViewerFlush(viewer));
76: PetscCall(VecRestoreArrayRead(X, &xv));
77: PetscFunctionReturn(PETSC_SUCCESS);
78: }
80: static PetscErrorCode VecView_Network_MPI(DM networkdm, Vec X, PetscViewer viewer)
81: {
82: PetscInt i, e, v, eStart, eEnd, vStart, vEnd, offset, nvar, len_loc, k;
83: const PetscScalar *xv;
84: MPI_Comm comm;
85: PetscMPIInt size, rank, tag = ((PetscObject)viewer)->tag, len;
86: Vec localX;
87: PetscBool ghostvtex;
88: PetscScalar *values;
89: PetscInt ne, nv, id;
90: MPI_Status status;
92: PetscFunctionBegin;
93: PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm));
94: PetscCallMPI(MPI_Comm_size(comm, &size));
95: PetscCallMPI(MPI_Comm_rank(comm, &rank));
97: PetscCall(DMGetLocalVector(networkdm, &localX));
98: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
99: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
100: PetscCall(VecGetArrayRead(localX, &xv));
102: PetscCall(VecGetLocalSize(localX, &len_loc));
104: PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
105: PetscCall(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: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &len_loc, 1, MPIU_INT, MPI_MAX, comm));
110: PetscCall(PetscMPIIntCast(len_loc, &len));
111: PetscCall(PetscCalloc1(len, &values));
113: if (rank == 0) PetscCall(PetscViewerASCIIPrintf(viewer, "Process [%d]\n", rank));
115: /* iterate over edges */
116: k = 2;
117: for (e = eStart; e < eEnd; e++) {
118: PetscCall(DMNetworkGetComponent(networkdm, e, ALL_COMPONENTS, NULL, NULL, &nvar));
119: if (!nvar) continue;
121: PetscCall(DMNetworkGetLocalVecOffset(networkdm, e, ALL_COMPONENTS, &offset));
122: PetscCall(DMNetworkGetGlobalEdgeIndex(networkdm, e, &id));
124: if (rank == 0) { /* print its own entries */
125: PetscCall(PetscViewerASCIIPrintf(viewer, " Edge %" PetscInt_FMT ":\n", id));
126: PetscCall(VecArrayPrint_private(viewer, nvar, xv + offset));
127: } else {
128: values[0] += 1; /* number of edges */
129: values[k++] = id;
130: values[k++] = nvar;
131: for (i = offset; i < offset + nvar; i++) values[k++] = xv[i];
132: }
133: }
135: /* iterate over vertices */
136: for (v = vStart; v < vEnd; v++) {
137: PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghostvtex));
138: if (ghostvtex) continue;
139: PetscCall(DMNetworkGetComponent(networkdm, v, ALL_COMPONENTS, NULL, NULL, &nvar));
140: if (!nvar) continue;
142: PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, ALL_COMPONENTS, &offset));
143: PetscCall(DMNetworkGetGlobalVertexIndex(networkdm, v, &id));
145: if (rank == 0) {
146: PetscCall(PetscViewerASCIIPrintf(viewer, " Vertex %" PetscInt_FMT ":\n", id));
147: PetscCall(VecArrayPrint_private(viewer, nvar, xv + offset));
148: } else {
149: values[1] += 1; /* number of vertices */
150: values[k++] = id;
151: values[k++] = nvar;
152: for (i = offset; i < offset + nvar; i++) values[k++] = xv[i];
153: }
154: }
156: if (rank == 0) {
157: /* proc[0] receives and prints messages */
158: for (PetscMPIInt j = 1; j < size; j++) {
159: PetscCall(PetscViewerASCIIPrintf(viewer, "Process [%d]\n", j));
161: PetscCallMPI(MPI_Recv(values, len, MPIU_SCALAR, j, tag, comm, &status));
163: ne = (PetscInt)PetscAbsScalar(values[0]);
164: nv = (PetscInt)PetscAbsScalar(values[1]);
166: /* print received edges */
167: k = 2;
168: for (i = 0; i < ne; i++) {
169: id = (PetscInt)PetscAbsScalar(values[k++]);
170: nvar = (PetscInt)PetscAbsScalar(values[k++]);
171: PetscCall(PetscViewerASCIIPrintf(viewer, " Edge %" PetscInt_FMT ":\n", id));
172: PetscCall(VecArrayPrint_private(viewer, nvar, values + k));
173: k += nvar;
174: }
176: /* print received vertices */
177: for (i = 0; i < nv; i++) {
178: id = (PetscInt)PetscAbsScalar(values[k++]);
179: nvar = (PetscInt)PetscAbsScalar(values[k++]);
180: PetscCall(PetscViewerASCIIPrintf(viewer, " Vertex %" PetscInt_FMT ":\n", id));
181: PetscCall(VecArrayPrint_private(viewer, nvar, values + k));
182: k += nvar;
183: }
184: }
185: } else {
186: /* sends values to proc[0] */
187: PetscCallMPI(MPIU_Send((void *)values, k, MPIU_SCALAR, 0, tag, comm));
188: }
190: PetscCall(PetscFree(values));
191: PetscCall(VecRestoreArrayRead(localX, &xv));
192: PetscCall(DMRestoreLocalVector(networkdm, &localX));
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
196: PETSC_SINGLE_LIBRARY_INTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
198: static PetscErrorCode VecView_Network(Vec v, PetscViewer viewer)
199: {
200: DM dm;
201: PetscBool isseq;
202: PetscBool isascii;
204: PetscFunctionBegin;
205: PetscCall(VecGetDM(v, &dm));
206: PetscCheck(dm, PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
207: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
208: PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq));
210: /* Use VecView_Network if the viewer is ASCII; use VecView_Seq/MPI for other viewer formats */
211: if (isascii) {
212: if (isseq) PetscCall(VecView_Network_Seq(dm, v, viewer));
213: else PetscCall(VecView_Network_MPI(dm, v, viewer));
214: } else {
215: if (isseq) PetscCall(VecView_Seq(v, viewer));
216: else PetscCall(VecView_MPI(v, viewer));
217: }
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: static PetscErrorCode DMCreateGlobalVector_Network(DM dm, Vec *vec)
222: {
223: DM_Network *network = (DM_Network *)dm->data;
225: PetscFunctionBegin;
226: PetscCall(DMCreateGlobalVector(network->plex, vec));
227: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (PetscErrorCodeFn *)VecView_Network));
228: PetscCall(VecSetDM(*vec, dm));
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: static PetscErrorCode DMCreateLocalVector_Network(DM dm, Vec *vec)
233: {
234: DM_Network *network = (DM_Network *)dm->data;
236: PetscFunctionBegin;
237: PetscCall(DMCreateLocalVector(network->plex, vec));
238: PetscCall(VecSetDM(*vec, dm));
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: PetscErrorCode DMNetworkInitializeToDefault_NonShared(DM dm)
243: {
244: DM_Network *network = (DM_Network *)dm->data;
246: PetscFunctionBegin;
247: network->Je = NULL;
248: network->Jv = NULL;
249: network->Jvptr = NULL;
250: network->userEdgeJacobian = PETSC_FALSE;
251: network->userVertexJacobian = PETSC_FALSE;
253: network->vertex.DofSection = NULL;
254: network->vertex.GlobalDofSection = NULL;
255: network->vertex.mapping = NULL;
256: network->vertex.sf = NULL;
258: network->edge.DofSection = NULL;
259: network->edge.GlobalDofSection = NULL;
260: network->edge.mapping = NULL;
261: network->edge.sf = NULL;
263: network->DataSection = NULL;
264: network->DofSection = NULL;
265: network->GlobalDofSection = NULL;
266: network->componentsetup = PETSC_FALSE;
268: network->plex = NULL;
270: network->component = NULL;
271: network->ncomponent = 0;
273: network->header = NULL;
274: network->cvalue = NULL;
275: network->componentdataarray = NULL;
277: network->max_comps_registered = DMNETWORK_MAX_COMP_REGISTERED_DEFAULT; /* return to default */
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
280: /* Default values for the parameters in DMNetwork */
281: PetscErrorCode DMNetworkInitializeToDefault(DM dm)
282: {
283: DM_Network *network = (DM_Network *)dm->data;
284: DMNetworkCloneShared cloneshared = network->cloneshared;
286: PetscFunctionBegin;
287: PetscCall(DMNetworkInitializeToDefault_NonShared(dm));
288: /* Default values for shared data */
289: cloneshared->refct = 1;
290: cloneshared->NVertices = 0;
291: cloneshared->NEdges = 0;
292: cloneshared->nVertices = 0;
293: cloneshared->nEdges = 0;
294: cloneshared->nsubnet = 0;
295: cloneshared->pStart = -1;
296: cloneshared->pEnd = -1;
297: cloneshared->vStart = -1;
298: cloneshared->vEnd = -1;
299: cloneshared->eStart = -1;
300: cloneshared->eEnd = -1;
301: cloneshared->vltog = NULL;
302: cloneshared->distributecalled = PETSC_FALSE;
304: cloneshared->subnet = NULL;
305: cloneshared->subnetvtx = NULL;
306: cloneshared->subnetedge = NULL;
307: cloneshared->svtx = NULL;
308: cloneshared->nsvtx = 0;
309: cloneshared->Nsvtx = 0;
310: cloneshared->svertices = NULL;
311: cloneshared->sedgelist = NULL;
312: cloneshared->svtable = NULL;
313: PetscFunctionReturn(PETSC_SUCCESS);
314: }
316: static PetscErrorCode DMInitialize_Network(DM dm)
317: {
318: PetscFunctionBegin;
319: PetscCall(DMSetDimension(dm, 1));
320: dm->ops->view = DMView_Network;
321: dm->ops->setfromoptions = DMSetFromOptions_Network;
322: dm->ops->clone = DMClone_Network;
323: dm->ops->setup = DMSetUp_Network;
324: dm->ops->createglobalvector = DMCreateGlobalVector_Network;
325: dm->ops->createlocalvector = DMCreateLocalVector_Network;
326: dm->ops->getlocaltoglobalmapping = NULL;
327: dm->ops->createfieldis = NULL;
328: dm->ops->createcoordinatedm = DMCreateCoordinateDM_Network;
329: dm->ops->createcellcoordinatedm = NULL;
330: dm->ops->getcoloring = NULL;
331: dm->ops->creatematrix = DMCreateMatrix_Network;
332: dm->ops->createinterpolation = NULL;
333: dm->ops->createinjection = NULL;
334: dm->ops->refine = NULL;
335: dm->ops->coarsen = NULL;
336: dm->ops->refinehierarchy = NULL;
337: dm->ops->coarsenhierarchy = NULL;
338: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Network;
339: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Network;
340: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Network;
341: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Network;
342: dm->ops->destroy = DMDestroy_Network;
343: dm->ops->createsubdm = NULL;
344: dm->ops->locatepoints = NULL;
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
347: /*
348: copies over the subnetid and index portions of the DMNetworkComponentHeader from original dm to the newdm
349: */
350: static PetscErrorCode DMNetworkCopyHeaderTopological(DM dm, DM newdm)
351: {
352: DM_Network *network = (DM_Network *)dm->data, *newnetwork = (DM_Network *)newdm->data;
353: PetscInt p, i, np, index, subnetid;
355: PetscFunctionBegin;
356: np = network->cloneshared->pEnd - network->cloneshared->pStart;
357: PetscCall(PetscCalloc2(np, &newnetwork->header, np, &newnetwork->cvalue));
358: for (i = 0; i < np; i++) {
359: p = i + network->cloneshared->pStart;
360: PetscCall(DMNetworkGetSubnetID(dm, p, &subnetid));
361: PetscCall(DMNetworkGetIndex(dm, p, &index));
362: newnetwork->header[i].index = index;
363: newnetwork->header[i].subnetid = subnetid;
364: newnetwork->header[i].size = NULL;
365: newnetwork->header[i].key = NULL;
366: newnetwork->header[i].offset = NULL;
367: newnetwork->header[i].nvar = NULL;
368: newnetwork->header[i].offsetvarrel = NULL;
369: newnetwork->header[i].ndata = 0;
370: newnetwork->header[i].maxcomps = DMNETWORK_MAX_COMP_AT_POINT_DEFAULT;
371: newnetwork->header[i].hsize = sizeof(struct _p_DMNetworkComponentHeader) / sizeof(sizeof(DMNetworkComponentGenericDataType));
372: }
373: PetscFunctionReturn(PETSC_SUCCESS);
374: }
376: PetscErrorCode DMClone_Network(DM dm, DM *newdm)
377: {
378: DM_Network *network = (DM_Network *)dm->data, *newnetwork = NULL;
380: PetscFunctionBegin;
381: network->cloneshared->refct++;
382: PetscCall(PetscNew(&newnetwork));
383: (*newdm)->data = newnetwork;
384: PetscCall(DMNetworkInitializeToDefault_NonShared(*newdm));
385: newnetwork->cloneshared = network->cloneshared; /* Share all data that can be cloneshared */
387: PetscCheck(network->plex, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_NULL, "Must call DMNetworkLayoutSetUp() first");
388: PetscCall(DMClone(network->plex, &newnetwork->plex));
389: PetscCall(DMNetworkCopyHeaderTopological(dm, *newdm));
390: PetscCall(DMNetworkInitializeNonTopological(*newdm)); /* initialize all non-topological data to the state after DMNetworkLayoutSetUp as been called */
391: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMNETWORK));
392: PetscCall(DMInitialize_Network(*newdm));
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: /* Developer Note: Be aware that the plex inside of the network does not have a coordinate plex.
397: */
398: PetscErrorCode DMCreateCoordinateDM_Network(DM dm, DM *cdm)
399: {
400: DM_Network *newnetwork = NULL;
401: PetscInt Nf;
402: const char *prefix;
404: PetscFunctionBegin;
405: PetscCall(DMClone(dm, cdm));
406: newnetwork = (DM_Network *)(*cdm)->data;
407: PetscCall(DMGetNumFields(newnetwork->plex, &Nf));
408: PetscCall(DMSetNumFields(*cdm, Nf)); /* consistency with the coordinate plex */
409: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
410: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*cdm, prefix));
411: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)*cdm, "cdm_"));
412: PetscFunctionReturn(PETSC_SUCCESS);
413: }
415: /*MC
416: DMNETWORK = "network" - A DM object that encapsulates an unstructured network. The implementation is based on the DM object
417: DMPlex that manages unstructured grids. Distributed networks use a non-overlapping partitioning of
418: the edges. In the local representation, Vecs contain all unknowns in the interior and shared boundary.
419: This is specified by a PetscSection object. Ownership in the global representation is determined by
420: ownership of the underlying DMPlex points. This is specified by another PetscSection object.
422: Level: intermediate
424: .seealso: `DMType`, `DMNetworkCreate()`, `DMCreate()`, `DMSetType()`
425: M*/
427: PETSC_EXTERN PetscErrorCode DMCreate_Network(DM dm)
428: {
429: DM_Network *network;
431: PetscFunctionBegin;
433: PetscCall(PetscNew(&network));
434: PetscCall(PetscNew(&network->cloneshared));
435: dm->data = network;
437: PetscCall(DMNetworkInitializeToDefault(dm));
438: PetscCall(DMInitialize_Network(dm));
439: PetscFunctionReturn(PETSC_SUCCESS);
440: }
442: /*@
443: DMNetworkCreate - Creates a DMNetwork object, which encapsulates an unstructured network.
445: Collective
447: Input Parameter:
448: . comm - The communicator for the DMNetwork object
450: Output Parameter:
451: . network - The DMNetwork object
453: Level: beginner
455: .seealso: `DMCreate()`
456: @*/
457: PetscErrorCode DMNetworkCreate(MPI_Comm comm, DM *network)
458: {
459: PetscFunctionBegin;
460: PetscAssertPointer(network, 2);
461: PetscCall(DMCreate(comm, network));
462: PetscCall(DMSetType(*network, DMNETWORK));
463: PetscFunctionReturn(PETSC_SUCCESS);
464: }