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: }