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, k;
84: const PetscScalar *xv;
85: MPI_Comm comm;
86: PetscMPIInt size, rank, tag = ((PetscObject)viewer)->tag, len;
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(MPI_IN_PLACE, &len_loc, 1, MPIU_INT, MPI_MAX, comm));
111: PetscCall(PetscMPIIntCast(len_loc, &len));
112: PetscCall(PetscCalloc1(len, &values));
114: if (rank == 0) PetscCall(PetscViewerASCIIPrintf(viewer, "Process [%d]\n", rank));
116: /* iterate over edges */
117: k = 2;
118: for (e = eStart; e < eEnd; e++) {
119: PetscCall(DMNetworkGetComponent(networkdm, e, ALL_COMPONENTS, NULL, NULL, &nvar));
120: if (!nvar) continue;
122: PetscCall(DMNetworkGetLocalVecOffset(networkdm, e, ALL_COMPONENTS, &offset));
123: PetscCall(DMNetworkGetGlobalEdgeIndex(networkdm, e, &id));
125: if (rank == 0) { /* print its own entries */
126: PetscCall(PetscViewerASCIIPrintf(viewer, " Edge %" PetscInt_FMT ":\n", id));
127: PetscCall(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: PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghostvtex));
139: if (ghostvtex) continue;
140: PetscCall(DMNetworkGetComponent(networkdm, v, ALL_COMPONENTS, NULL, NULL, &nvar));
141: if (!nvar) continue;
143: PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, ALL_COMPONENTS, &offset));
144: PetscCall(DMNetworkGetGlobalVertexIndex(networkdm, v, &id));
146: if (rank == 0) {
147: PetscCall(PetscViewerASCIIPrintf(viewer, " Vertex %" PetscInt_FMT ":\n", id));
148: PetscCall(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 (PetscMPIInt j = 1; j < size; j++) {
160: PetscCall(PetscViewerASCIIPrintf(viewer, "Process [%d]\n", j));
162: PetscCallMPI(MPI_Recv(values, 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: PetscCall(PetscViewerASCIIPrintf(viewer, " Edge %" PetscInt_FMT ":\n", id));
173: PetscCall(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: PetscCall(PetscViewerASCIIPrintf(viewer, " Vertex %" PetscInt_FMT ":\n", id));
182: PetscCall(VecArrayPrint_private(viewer, nvar, values + k));
183: k += nvar;
184: }
185: }
186: } else {
187: /* sends values to proc[0] */
188: PetscCallMPI(MPIU_Send((void *)values, k, MPIU_SCALAR, 0, tag, comm));
189: }
191: PetscCall(PetscFree(values));
192: PetscCall(VecRestoreArrayRead(localX, &xv));
193: PetscCall(DMRestoreLocalVector(networkdm, &localX));
194: PetscFunctionReturn(PETSC_SUCCESS);
195: }
197: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
199: static PetscErrorCode VecView_Network(Vec v, PetscViewer viewer)
200: {
201: DM dm;
202: PetscBool isseq;
203: PetscBool iascii;
205: PetscFunctionBegin;
206: PetscCall(VecGetDM(v, &dm));
207: PetscCheck(dm, PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
208: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
209: PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq));
211: /* Use VecView_Network if the viewer is ASCII; use VecView_Seq/MPI for other viewer formats */
212: if (iascii) {
213: if (isseq) PetscCall(VecView_Network_Seq(dm, v, viewer));
214: else PetscCall(VecView_Network_MPI(dm, v, viewer));
215: } else {
216: if (isseq) PetscCall(VecView_Seq(v, viewer));
217: else PetscCall(VecView_MPI(v, viewer));
218: }
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: static PetscErrorCode DMCreateGlobalVector_Network(DM dm, Vec *vec)
223: {
224: DM_Network *network = (DM_Network *)dm->data;
226: PetscFunctionBegin;
227: PetscCall(DMCreateGlobalVector(network->plex, vec));
228: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Network));
229: PetscCall(VecSetDM(*vec, dm));
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: static PetscErrorCode DMCreateLocalVector_Network(DM dm, Vec *vec)
234: {
235: DM_Network *network = (DM_Network *)dm->data;
237: PetscFunctionBegin;
238: PetscCall(DMCreateLocalVector(network->plex, vec));
239: PetscCall(VecSetDM(*vec, dm));
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: PetscErrorCode DMNetworkInitializeToDefault_NonShared(DM dm)
244: {
245: DM_Network *network = (DM_Network *)dm->data;
247: PetscFunctionBegin;
248: network->Je = NULL;
249: network->Jv = NULL;
250: network->Jvptr = NULL;
251: network->userEdgeJacobian = PETSC_FALSE;
252: network->userVertexJacobian = PETSC_FALSE;
254: network->vertex.DofSection = NULL;
255: network->vertex.GlobalDofSection = NULL;
256: network->vertex.mapping = NULL;
257: network->vertex.sf = NULL;
259: network->edge.DofSection = NULL;
260: network->edge.GlobalDofSection = NULL;
261: network->edge.mapping = NULL;
262: network->edge.sf = NULL;
264: network->DataSection = NULL;
265: network->DofSection = NULL;
266: network->GlobalDofSection = NULL;
267: network->componentsetup = PETSC_FALSE;
269: network->plex = NULL;
271: network->component = NULL;
272: network->ncomponent = 0;
274: network->header = NULL;
275: network->cvalue = NULL;
276: network->componentdataarray = NULL;
278: network->max_comps_registered = DMNETWORK_MAX_COMP_REGISTERED_DEFAULT; /* return to default */
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
281: /* Default values for the parameters in DMNetwork */
282: PetscErrorCode DMNetworkInitializeToDefault(DM dm)
283: {
284: DM_Network *network = (DM_Network *)dm->data;
285: DMNetworkCloneShared cloneshared = network->cloneshared;
287: PetscFunctionBegin;
288: PetscCall(DMNetworkInitializeToDefault_NonShared(dm));
289: /* Default values for shared data */
290: cloneshared->refct = 1;
291: cloneshared->NVertices = 0;
292: cloneshared->NEdges = 0;
293: cloneshared->nVertices = 0;
294: cloneshared->nEdges = 0;
295: cloneshared->nsubnet = 0;
296: cloneshared->pStart = -1;
297: cloneshared->pEnd = -1;
298: cloneshared->vStart = -1;
299: cloneshared->vEnd = -1;
300: cloneshared->eStart = -1;
301: cloneshared->eEnd = -1;
302: cloneshared->vltog = NULL;
303: cloneshared->distributecalled = PETSC_FALSE;
305: cloneshared->subnet = NULL;
306: cloneshared->subnetvtx = NULL;
307: cloneshared->subnetedge = NULL;
308: cloneshared->svtx = NULL;
309: cloneshared->nsvtx = 0;
310: cloneshared->Nsvtx = 0;
311: cloneshared->svertices = NULL;
312: cloneshared->sedgelist = NULL;
313: cloneshared->svtable = NULL;
314: PetscFunctionReturn(PETSC_SUCCESS);
315: }
317: static PetscErrorCode DMInitialize_Network(DM dm)
318: {
319: PetscFunctionBegin;
320: PetscCall(DMSetDimension(dm, 1));
321: dm->ops->view = DMView_Network;
322: dm->ops->setfromoptions = DMSetFromOptions_Network;
323: dm->ops->clone = DMClone_Network;
324: dm->ops->setup = DMSetUp_Network;
325: dm->ops->createglobalvector = DMCreateGlobalVector_Network;
326: dm->ops->createlocalvector = DMCreateLocalVector_Network;
327: dm->ops->getlocaltoglobalmapping = NULL;
328: dm->ops->createfieldis = NULL;
329: dm->ops->createcoordinatedm = DMCreateCoordinateDM_Network;
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: }