Actual source code: networkcreate.c
1: #define PETSCDM_DLL
2: #include <petsc/private/dmnetworkimpl.h>
3: #include <petsc/private/vecimpl.h>
5: static PetscErrorCode DMSetFromOptions_Network(DM dm, PetscOptionItems *PetscOptionsObject)
6: {
7: PetscFunctionBegin;
8: PetscOptionsHeadBegin(PetscOptionsObject, "DMNetwork Options");
9: PetscOptionsHeadEnd();
10: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
29: for (i = 0; i < n; i++) {
30: #if defined(PETSC_USE_COMPLEX)
31: if (PetscImaginaryPart(xv[i]) > 0.0) {
32: PetscCall(PetscViewerASCIIPrintf(viewer, " %g + %g i\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i])));
33: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
34: PetscCall(PetscViewerASCIIPrintf(viewer, " %g - %g i\n", (double)PetscRealPart(xv[i]), -(double)PetscImaginaryPart(xv[i])));
35: } else PetscCall(PetscViewerASCIIPrintf(viewer, " %g\n", (double)PetscRealPart(xv[i])));
36: #else
37: PetscCall(PetscViewerASCIIPrintf(viewer, " %g\n", (double)xv[i]));
38: #endif
39: }
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: static PetscErrorCode VecView_Network_Seq(DM networkdm, Vec X, PetscViewer viewer)
44: {
45: PetscInt e, v, Start, End, offset, nvar, id;
46: const PetscScalar *xv;
48: PetscFunctionBegin;
49: PetscCall(VecGetArrayRead(X, &xv));
51: /* iterate over edges */
52: PetscCall(DMNetworkGetEdgeRange(networkdm, &Start, &End));
53: for (e = Start; e < End; e++) {
54: PetscCall(DMNetworkGetComponent(networkdm, e, ALL_COMPONENTS, NULL, NULL, &nvar));
55: if (!nvar) continue;
57: PetscCall(DMNetworkGetLocalVecOffset(networkdm, e, ALL_COMPONENTS, &offset));
58: PetscCall(DMNetworkGetGlobalEdgeIndex(networkdm, e, &id));
60: PetscCall(PetscViewerASCIIPrintf(viewer, " Edge %" PetscInt_FMT ":\n", id));
61: PetscCall(VecArrayPrint_private(viewer, nvar, xv + offset));
62: }
64: /* iterate over vertices */
65: PetscCall(DMNetworkGetVertexRange(networkdm, &Start, &End));
66: for (v = Start; v < End; v++) {
67: PetscCall(DMNetworkGetComponent(networkdm, v, ALL_COMPONENTS, NULL, NULL, &nvar));
68: if (!nvar) continue;
70: PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, ALL_COMPONENTS, &offset));
71: PetscCall(DMNetworkGetGlobalVertexIndex(networkdm, v, &id));
73: PetscCall(PetscViewerASCIIPrintf(viewer, " Vertex %" PetscInt_FMT ":\n", id));
74: PetscCall(VecArrayPrint_private(viewer, nvar, xv + offset));
75: }
76: PetscCall(PetscViewerFlush(viewer));
77: PetscCall(VecRestoreArrayRead(X, &xv));
78: PetscFunctionReturn(PETSC_SUCCESS);
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 ne, nv, id;
91: MPI_Status status;
93: PetscFunctionBegin;
94: PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm));
95: PetscCallMPI(MPI_Comm_size(comm, &size));
96: PetscCallMPI(MPI_Comm_rank(comm, &rank));
98: PetscCall(DMGetLocalVector(networkdm, &localX));
99: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
100: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
101: PetscCall(VecGetArrayRead(localX, &xv));
103: PetscCall(VecGetLocalSize(localX, &len_loc));
105: PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
106: PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
107: len_loc += 2 * (1 + eEnd - eStart + vEnd - vStart);
109: /* values = [nedges, nvertices; id, nvar, xedge; ...; id, nvars, xvertex;...], to be sent to proc[0] */
110: PetscCallMPI(MPIU_Allreduce(&len_loc, &len, 1, MPIU_INT, MPI_MAX, comm));
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, (PetscMPIInt)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_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
198: static PetscErrorCode VecView_Network(Vec v, PetscViewer viewer)
199: {
200: DM dm;
201: PetscBool isseq;
202: PetscBool iascii;
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, &iascii));
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 (iascii) {
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, (void (*)(void))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->getcoloring = NULL;
330: dm->ops->creatematrix = DMCreateMatrix_Network;
331: dm->ops->createinterpolation = NULL;
332: dm->ops->createinjection = NULL;
333: dm->ops->refine = NULL;
334: dm->ops->coarsen = NULL;
335: dm->ops->refinehierarchy = NULL;
336: dm->ops->coarsenhierarchy = NULL;
337: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Network;
338: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Network;
339: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Network;
340: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Network;
341: dm->ops->destroy = DMDestroy_Network;
342: dm->ops->createsubdm = NULL;
343: dm->ops->locatepoints = NULL;
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
346: /*
347: copies over the subnetid and index portions of the DMNetworkComponentHeader from original dm to the newdm
348: */
349: static PetscErrorCode DMNetworkCopyHeaderTopological(DM dm, DM newdm)
350: {
351: DM_Network *network = (DM_Network *)dm->data, *newnetwork = (DM_Network *)newdm->data;
352: PetscInt p, i, np, index, subnetid;
354: PetscFunctionBegin;
355: np = network->cloneshared->pEnd - network->cloneshared->pStart;
356: PetscCall(PetscCalloc2(np, &newnetwork->header, np, &newnetwork->cvalue));
357: for (i = 0; i < np; i++) {
358: p = i + network->cloneshared->pStart;
359: PetscCall(DMNetworkGetSubnetID(dm, p, &subnetid));
360: PetscCall(DMNetworkGetIndex(dm, p, &index));
361: newnetwork->header[i].index = index;
362: newnetwork->header[i].subnetid = subnetid;
363: newnetwork->header[i].size = NULL;
364: newnetwork->header[i].key = NULL;
365: newnetwork->header[i].offset = NULL;
366: newnetwork->header[i].nvar = NULL;
367: newnetwork->header[i].offsetvarrel = NULL;
368: newnetwork->header[i].ndata = 0;
369: newnetwork->header[i].maxcomps = DMNETWORK_MAX_COMP_AT_POINT_DEFAULT;
370: newnetwork->header[i].hsize = sizeof(struct _p_DMNetworkComponentHeader) / sizeof(sizeof(DMNetworkComponentGenericDataType));
371: }
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: PetscErrorCode DMClone_Network(DM dm, DM *newdm)
376: {
377: DM_Network *network = (DM_Network *)dm->data, *newnetwork = NULL;
379: PetscFunctionBegin;
380: network->cloneshared->refct++;
381: PetscCall(PetscNew(&newnetwork));
382: (*newdm)->data = newnetwork;
383: PetscCall(DMNetworkInitializeToDefault_NonShared(*newdm));
384: newnetwork->cloneshared = network->cloneshared; /* Share all data that can be cloneshared */
386: PetscCheck(network->plex, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_NULL, "Must call DMNetworkLayoutSetUp() first");
387: PetscCall(DMClone(network->plex, &newnetwork->plex));
388: PetscCall(DMNetworkCopyHeaderTopological(dm, *newdm));
389: PetscCall(DMNetworkInitializeNonTopological(*newdm)); /* initialize all non-topological data to the state after DMNetworkLayoutSetUp as been called */
390: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMNETWORK));
391: PetscCall(DMInitialize_Network(*newdm));
392: PetscFunctionReturn(PETSC_SUCCESS);
393: }
395: /* Developer Note: Be aware that the plex inside of the network does not have a coordinate plex.
396: */
397: PetscErrorCode DMCreateCoordinateDM_Network(DM dm, DM *cdm)
398: {
399: DM_Network *newnetwork = NULL;
400: PetscInt Nf;
401: const char *prefix;
403: PetscFunctionBegin;
404: PetscCall(DMClone(dm, cdm));
405: newnetwork = (DM_Network *)(*cdm)->data;
406: PetscCall(DMGetNumFields(newnetwork->plex, &Nf));
407: PetscCall(DMSetNumFields(*cdm, Nf)); /* consistency with the coordinate plex */
408: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
409: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*cdm, prefix));
410: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)*cdm, "cdm_"));
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: /*MC
415: DMNETWORK = "network" - A DM object that encapsulates an unstructured network. The implementation is based on the DM object
416: DMPlex that manages unstructured grids. Distributed networks use a non-overlapping partitioning of
417: the edges. In the local representation, Vecs contain all unknowns in the interior and shared boundary.
418: This is specified by a PetscSection object. Ownership in the global representation is determined by
419: ownership of the underlying DMPlex points. This is specified by another PetscSection object.
421: Level: intermediate
423: .seealso: `DMType`, `DMNetworkCreate()`, `DMCreate()`, `DMSetType()`
424: M*/
426: PETSC_EXTERN PetscErrorCode DMCreate_Network(DM dm)
427: {
428: DM_Network *network;
430: PetscFunctionBegin;
432: PetscCall(PetscNew(&network));
433: PetscCall(PetscNew(&network->cloneshared));
434: dm->data = network;
436: PetscCall(DMNetworkInitializeToDefault(dm));
437: PetscCall(DMInitialize_Network(dm));
438: PetscFunctionReturn(PETSC_SUCCESS);
439: }
441: /*@
442: DMNetworkCreate - Creates a DMNetwork object, which encapsulates an unstructured network.
444: Collective
446: Input Parameter:
447: . comm - The communicator for the DMNetwork object
449: Output Parameter:
450: . network - The DMNetwork object
452: Level: beginner
454: .seealso: `DMCreate()`
455: @*/
456: PetscErrorCode DMNetworkCreate(MPI_Comm comm, DM *network)
457: {
458: PetscFunctionBegin;
459: PetscAssertPointer(network, 2);
460: PetscCall(DMCreate(comm, network));
461: PetscCall(DMSetType(*network, DMNETWORK));
462: PetscFunctionReturn(PETSC_SUCCESS);
463: }