Actual source code: networkcreate.c

      1: #define PETSCDM_DLL
      2: #include <petsc/private/dmnetworkimpl.h>
      3: #include <petsc/private/vecimpl.h>

      5: PetscErrorCode  DMSetFromOptions_Network(PetscOptionItems *PetscOptionsObject,DM dm)
      6: {

     11:   PetscOptionsHead(PetscOptionsObject,"DMNetwork Options");
     12:   PetscOptionsTail();
     13:   return(0);
     14: }

     16: /* External function declarations here */
     17: extern PetscErrorCode DMCreateMatrix_Network(DM, Mat*);
     18: extern PetscErrorCode DMDestroy_Network(DM);
     19: extern PetscErrorCode DMView_Network(DM, PetscViewer);
     20: extern PetscErrorCode DMGlobalToLocalBegin_Network(DM, Vec, InsertMode, Vec);
     21: extern PetscErrorCode DMGlobalToLocalEnd_Network(DM, Vec, InsertMode, Vec);
     22: extern PetscErrorCode DMLocalToGlobalBegin_Network(DM, Vec, InsertMode, Vec);
     23: extern PetscErrorCode DMLocalToGlobalEnd_Network(DM, Vec, InsertMode, Vec);
     24: extern PetscErrorCode DMSetUp_Network(DM);
     25: extern PetscErrorCode DMClone_Network(DM, DM*);

     27: static PetscErrorCode VecArrayPrint_private(PetscViewer viewer,PetscInt n,const PetscScalar *xv)
     28: {
     30:   PetscInt       i;

     33:   for (i=0; i<n; i++) {
     34: #if defined(PETSC_USE_COMPLEX)
     35:     if (PetscImaginaryPart(xv[i]) > 0.0) {
     36:       PetscViewerASCIIPrintf(viewer,"    %g + %g i\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
     37:     } else if (PetscImaginaryPart(xv[i]) < 0.0) {
     38:       PetscViewerASCIIPrintf(viewer,"    %g - %g i\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
     39:     } else {
     40:       PetscViewerASCIIPrintf(viewer,"    %g\n",(double)PetscRealPart(xv[i]));
     41:     }
     42: #else
     43:     PetscViewerASCIIPrintf(viewer,"    %g\n",(double)xv[i]);
     44: #endif
     45:   }
     46:   return(0);
     47: }

     49: static PetscErrorCode VecView_Network_Seq(DM networkdm,Vec X,PetscViewer viewer)
     50: {
     51:   PetscErrorCode    ierr;
     52:   PetscInt          e,v,Start,End,offset,nvar,id;
     53:   const PetscScalar *xv;

     56:   VecGetArrayRead(X,&xv);

     58:   /* iterate over edges */
     59:   DMNetworkGetEdgeRange(networkdm,&Start,&End);
     60:   for (e=Start; e<End; e++) {
     61:     DMNetworkGetComponent(networkdm,e,ALL_COMPONENTS,NULL,NULL,&nvar);
     62:     if (!nvar) continue;

     64:     DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&offset);
     65:     DMNetworkGetGlobalEdgeIndex(networkdm,e,&id);

     67:     PetscViewerASCIIPrintf(viewer,"  Edge %D:\n",id);
     68:     VecArrayPrint_private(viewer,nvar,xv+offset);
     69:   }

     71:   /* iterate over vertices */
     72:   DMNetworkGetVertexRange(networkdm,&Start,&End);
     73:   for (v=Start; v<End; v++) {
     74:     DMNetworkGetComponent(networkdm,v,ALL_COMPONENTS,NULL,NULL,&nvar);
     75:     if (!nvar) continue;

     77:     DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&offset);
     78:     DMNetworkGetGlobalVertexIndex(networkdm,v,&id);

     80:     PetscViewerASCIIPrintf(viewer,"  Vertex %D:\n",id);
     81:     VecArrayPrint_private(viewer,nvar,xv+offset);
     82:   }
     83:   PetscViewerFlush(viewer);
     84:   VecRestoreArrayRead(X,&xv);
     85:   return(0);
     86: }

     88: static PetscErrorCode VecView_Network_MPI(DM networkdm,Vec X,PetscViewer viewer)
     89: {
     90:   PetscErrorCode    ierr;
     91:   PetscInt          i,e,v,eStart,eEnd,vStart,vEnd,offset,nvar,len_loc,len,k;
     92:   const PetscScalar *xv;
     93:   MPI_Comm          comm;
     94:   PetscMPIInt       size,rank,tag = ((PetscObject)viewer)->tag;
     95:   Vec               localX;
     96:   PetscBool         ghostvtex;
     97:   PetscScalar       *values;
     98:   PetscInt          j,ne,nv,id;
     99:   MPI_Status        status;

    102:   PetscObjectGetComm((PetscObject)networkdm,&comm);
    103:   MPI_Comm_size(comm,&size);
    104:   MPI_Comm_rank(comm,&rank);

    106:   DMGetLocalVector(networkdm,&localX);
    107:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
    108:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
    109:   VecGetArrayRead(localX,&xv);

    111:   VecGetLocalSize(localX,&len_loc);

    113:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
    114:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
    115:   len_loc += 2*(1 + eEnd-eStart + vEnd-vStart);

    117:   /* values = [nedges, nvertices; id, nvar, xedge; ...; id, nvars, xvertex;...], to be sent to proc[0] */
    118:   MPI_Allreduce(&len_loc,&len,1,MPIU_INT,MPI_MAX,comm);
    119:   PetscCalloc1(len,&values);

    121:   if (!rank) {
    122:     PetscViewerASCIIPrintf(viewer,"Process [%d]\n",rank);
    123:   }

    125:   /* iterate over edges */
    126:   k = 2;
    127:   for (e=eStart; e<eEnd; e++) {
    128:     DMNetworkGetComponent(networkdm,e,ALL_COMPONENTS,NULL,NULL,&nvar);
    129:     if (!nvar) continue;

    131:     DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&offset);
    132:     DMNetworkGetGlobalEdgeIndex(networkdm,e,&id);

    134:     if (!rank) { /* print its own entries */
    135:       PetscViewerASCIIPrintf(viewer,"  Edge %D:\n",id);
    136:       VecArrayPrint_private(viewer,nvar,xv+offset);
    137:     } else {
    138:       values[0]  += 1; /* number of edges */
    139:       values[k++] = id;
    140:       values[k++] = nvar;
    141:       for (i=offset; i< offset+nvar; i++) values[k++] = xv[i];
    142:     }
    143:   }

    145:   /* iterate over vertices */
    146:   for (v=vStart; v<vEnd; v++) {
    147:     DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
    148:     if (ghostvtex) continue;
    149:     DMNetworkGetComponent(networkdm,v,ALL_COMPONENTS,NULL,NULL,&nvar);
    150:     if (!nvar) continue;

    152:     DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&offset);
    153:     DMNetworkGetGlobalVertexIndex(networkdm,v,&id);

    155:     if (!rank) {
    156:       PetscViewerASCIIPrintf(viewer,"  Vertex %D:\n",id);
    157:       VecArrayPrint_private(viewer,nvar,xv+offset);
    158:     } else {
    159:       values[1]  += 1; /* number of vertices */
    160:       values[k++] = id;
    161:       values[k++] = nvar;
    162:       for (i=offset; i< offset+nvar; i++) values[k++] = xv[i];
    163:     }
    164:   }

    166:   if (!rank) {
    167:     /* proc[0] receives and prints messages */
    168:     for (j=1; j<size; j++) {
    169:       PetscViewerASCIIPrintf(viewer,"Process [%d]\n",j);

    171:       MPI_Recv(values,(PetscMPIInt)len,MPIU_SCALAR,j,tag,comm,&status);

    173:       ne = (PetscInt)PetscAbsScalar(values[0]);
    174:       nv = (PetscInt)PetscAbsScalar(values[1]);

    176:       /* print received edges */
    177:       k = 2;
    178:       for (i=0; i<ne; i++) {
    179:         id   = (PetscInt)PetscAbsScalar(values[k++]);
    180:         nvar = (PetscInt)PetscAbsScalar(values[k++]);
    181:         PetscViewerASCIIPrintf(viewer,"  Edge %D:\n",id);
    182:         VecArrayPrint_private(viewer,nvar,values+k);
    183:         k   += nvar;
    184:       }

    186:       /* print received vertices */
    187:       for (i=0; i<nv; i++) {
    188:         id   = (PetscInt)PetscAbsScalar(values[k++]);
    189:         nvar = (PetscInt)PetscAbsScalar(values[k++]);
    190:         PetscViewerASCIIPrintf(viewer,"  Vertex %D:\n",id);
    191:         VecArrayPrint_private(viewer,nvar,values+k);
    192:         k   += nvar;
    193:       }
    194:     }
    195:   } else {
    196:     /* sends values to proc[0] */
    197:     MPI_Send((void*)values,k,MPIU_SCALAR,0,tag,comm);
    198:   }

    200:   PetscFree(values);
    201:   VecRestoreArrayRead(localX,&xv);
    202:   DMRestoreLocalVector(networkdm,&localX);
    203:   return(0);
    204: }

    206: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec,PetscViewer);

    208: PetscErrorCode VecView_Network(Vec v,PetscViewer viewer)
    209: {
    210:   DM             dm;
    212:   PetscBool      isseq;
    213:   PetscBool      iascii;

    216:   VecGetDM(v,&dm);
    217:   if (!dm) SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_WRONG,"Vector not generated from a DM");
    218:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
    219:   PetscObjectTypeCompare((PetscObject)v,VECSEQ,&isseq);

    221:   /* Use VecView_Network if the viewer is ASCII; use VecView_Seq/MPI for other viewer formats */
    222:   if (iascii) {
    223:     if (isseq) {
    224:       VecView_Network_Seq(dm,v,viewer);
    225:     } else {
    226:       VecView_Network_MPI(dm,v,viewer);
    227:     }
    228:   } else {
    229:     if (isseq) {
    230:       VecView_Seq(v,viewer);
    231:     } else {
    232:       VecView_MPI(v,viewer);
    233:     }
    234:   }
    235:   return(0);
    236: }

    238: static PetscErrorCode DMCreateGlobalVector_Network(DM dm,Vec *vec)
    239: {
    241:   DM_Network     *network = (DM_Network*) dm->data;

    244:   DMCreateGlobalVector(network->plex,vec);
    245:   VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Network);
    246:   VecSetDM(*vec,dm);
    247:   return(0);
    248: }

    250: static PetscErrorCode DMCreateLocalVector_Network(DM dm,Vec *vec)
    251: {
    253:   DM_Network     *network = (DM_Network*) dm->data;

    256:   DMCreateLocalVector(network->plex,vec);
    257:   VecSetDM(*vec,dm);
    258:   return(0);
    259: }

    261: PetscErrorCode DMInitialize_Network(DM dm)
    262: {

    266:   DMSetDimension(dm,1);
    267:   dm->ops->view                            = DMView_Network;
    268:   dm->ops->setfromoptions                  = DMSetFromOptions_Network;
    269:   dm->ops->clone                           = DMClone_Network;
    270:   dm->ops->setup                           = DMSetUp_Network;
    271:   dm->ops->createglobalvector              = DMCreateGlobalVector_Network;
    272:   dm->ops->createlocalvector               = DMCreateLocalVector_Network;
    273:   dm->ops->getlocaltoglobalmapping         = NULL;
    274:   dm->ops->createfieldis                   = NULL;
    275:   dm->ops->createcoordinatedm              = NULL;
    276:   dm->ops->getcoloring                     = NULL;
    277:   dm->ops->creatematrix                    = DMCreateMatrix_Network;
    278:   dm->ops->createinterpolation             = NULL;
    279:   dm->ops->createinjection                 = NULL;
    280:   dm->ops->refine                          = NULL;
    281:   dm->ops->coarsen                         = NULL;
    282:   dm->ops->refinehierarchy                 = NULL;
    283:   dm->ops->coarsenhierarchy                = NULL;
    284:   dm->ops->globaltolocalbegin              = DMGlobalToLocalBegin_Network;
    285:   dm->ops->globaltolocalend                = DMGlobalToLocalEnd_Network;
    286:   dm->ops->localtoglobalbegin              = DMLocalToGlobalBegin_Network;
    287:   dm->ops->localtoglobalend                = DMLocalToGlobalEnd_Network;
    288:   dm->ops->destroy                         = DMDestroy_Network;
    289:   dm->ops->createsubdm                     = NULL;
    290:   dm->ops->locatepoints                    = NULL;
    291:   return(0);
    292: }

    294: PetscErrorCode DMClone_Network(DM dm, DM *newdm)
    295: {
    296:   DM_Network     *network = (DM_Network *) dm->data;

    300:   network->refct++;
    301:   (*newdm)->data = network;
    302:   PetscObjectChangeTypeName((PetscObject) *newdm, DMNETWORK);
    303:   DMInitialize_Network(*newdm);
    304:   return(0);
    305: }

    307: /*MC
    308:   DMNETWORK = "network" - A DM object that encapsulates an unstructured network. The implementation is based on the DM object
    309:                           DMPlex that manages unstructured grids. Distributed networks use a non-overlapping partitioning of
    310:                           the edges. In the local representation, Vecs contain all unknowns in the interior and shared boundary.
    311:                           This is specified by a PetscSection object. Ownership in the global representation is determined by
    312:                           ownership of the underlying DMPlex points. This is specified by another PetscSection object.

    314:   Level: intermediate

    316: .seealso: DMType, DMNetworkCreate(), DMCreate(), DMSetType()
    317: M*/

    319: PETSC_EXTERN PetscErrorCode DMCreate_Network(DM dm)
    320: {
    321:   DM_Network     *network;

    326:   PetscNewLog(dm,&network);
    327:   dm->data = network;

    329:   network->refct     = 1;
    330:   network->NVertices = 0;
    331:   network->NEdges    = 0;
    332:   network->nVertices = 0;
    333:   network->nEdges    = 0;
    334:   network->nsubnet   = 0;

    337:   DMInitialize_Network(dm);
    338:   return(0);
    339: }

    341: /*@
    342:   DMNetworkCreate - Creates a DMNetwork object, which encapsulates an unstructured network.

    344:   Collective

    346:   Input Parameter:
    347: . comm - The communicator for the DMNetwork object

    349:   Output Parameter:
    350: . network  - The DMNetwork object

    352:   Level: beginner

    354: @*/
    355: PetscErrorCode DMNetworkCreate(MPI_Comm comm, DM *network)
    356: {

    361:   DMCreate(comm, network);
    362:   DMSetType(*network, DMNETWORK);
    363:   return(0);
    364: }