1 #define PETSCDM_DLL 2 #include <petsc/private/dmnetworkimpl.h> /*I "petscdmnetwork.h" I*/ 3 #include <petscdmda.h> 4 #include <petsc/private/vecimpl.h> 5 6 PetscErrorCode DMSetFromOptions_Network(PetscOptionItems *PetscOptionsObject,DM dm) 7 { 8 PetscErrorCode ierr; 9 10 PetscFunctionBegin; 11 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 12 ierr = PetscOptionsHead(PetscOptionsObject,"DMNetwork Options");CHKERRQ(ierr); 13 ierr = PetscOptionsTail();CHKERRQ(ierr); 14 PetscFunctionReturn(0); 15 } 16 17 /* External function declarations here */ 18 extern PetscErrorCode DMCreateMatrix_Network(DM, Mat*); 19 extern PetscErrorCode DMDestroy_Network(DM); 20 extern PetscErrorCode DMView_Network(DM, PetscViewer); 21 extern PetscErrorCode DMGlobalToLocalBegin_Network(DM, Vec, InsertMode, Vec); 22 extern PetscErrorCode DMGlobalToLocalEnd_Network(DM, Vec, InsertMode, Vec); 23 extern PetscErrorCode DMLocalToGlobalBegin_Network(DM, Vec, InsertMode, Vec); 24 extern PetscErrorCode DMLocalToGlobalEnd_Network(DM, Vec, InsertMode, Vec); 25 extern PetscErrorCode DMSetUp_Network(DM); 26 extern PetscErrorCode DMClone_Network(DM, DM*); 27 28 static PetscErrorCode VecArrayPrint_private(PetscViewer viewer,PetscInt n,const PetscScalar *xv) 29 { 30 PetscErrorCode ierr; 31 PetscInt i; 32 33 PetscFunctionBegin; 34 for (i=0; i<n; i++) { 35 #if defined(PETSC_USE_COMPLEX) 36 if (PetscImaginaryPart(xv[i]) > 0.0) { 37 ierr = PetscViewerASCIIPrintf(viewer," %g + %g i\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));CHKERRQ(ierr); 38 } else if (PetscImaginaryPart(xv[i]) < 0.0) { 39 ierr = PetscViewerASCIIPrintf(viewer," %g - %g i\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));CHKERRQ(ierr); 40 } else { 41 ierr = PetscViewerASCIIPrintf(viewer," %g\n",(double)PetscRealPart(xv[i]));CHKERRQ(ierr); 42 } 43 #else 44 ierr = PetscViewerASCIIPrintf(viewer," %g\n",(double)xv[i]);CHKERRQ(ierr); 45 #endif 46 } 47 PetscFunctionReturn(0); 48 } 49 50 static PetscErrorCode VecView_Network_Seq(DM networkdm,Vec X,PetscViewer viewer) 51 { 52 PetscErrorCode ierr; 53 PetscInt e,v,Start,End,offset,nvar,id; 54 const PetscScalar *xv; 55 56 PetscFunctionBegin; 57 ierr = VecGetArrayRead(X,&xv);CHKERRQ(ierr); 58 59 /* iterate over edges */ 60 ierr = DMNetworkGetEdgeRange(networkdm,&Start,&End);CHKERRQ(ierr); 61 for (e=Start; e<End; e++) { 62 ierr = DMNetworkGetNumVariables(networkdm,e,&nvar);CHKERRQ(ierr); 63 if (!nvar) continue; 64 65 ierr = DMNetworkGetVariableOffset(networkdm,e,&offset);CHKERRQ(ierr); 66 ierr = DMNetworkGetGlobalEdgeIndex(networkdm,e,&id);CHKERRQ(ierr); 67 68 ierr = PetscViewerASCIIPrintf(viewer," Edge %D:\n",id);CHKERRQ(ierr); 69 ierr = VecArrayPrint_private(viewer,nvar,xv+offset);CHKERRQ(ierr); 70 } 71 72 /* iterate over vertices */ 73 ierr = DMNetworkGetVertexRange(networkdm,&Start,&End);CHKERRQ(ierr); 74 for (v=Start; v<End; v++) { 75 ierr = DMNetworkGetNumVariables(networkdm,v,&nvar);CHKERRQ(ierr); 76 if (!nvar) continue; 77 78 ierr = DMNetworkGetVariableOffset(networkdm,v,&offset);CHKERRQ(ierr); 79 ierr = DMNetworkGetGlobalVertexIndex(networkdm,v,&id);CHKERRQ(ierr); 80 81 ierr = PetscViewerASCIIPrintf(viewer," Vertex %D:\n",id);CHKERRQ(ierr); 82 ierr = VecArrayPrint_private(viewer,nvar,xv+offset);CHKERRQ(ierr); 83 } 84 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 85 ierr = VecRestoreArrayRead(X,&xv);CHKERRQ(ierr); 86 PetscFunctionReturn(0); 87 } 88 89 static PetscErrorCode VecView_Network_MPI(DM networkdm,Vec X,PetscViewer viewer) 90 { 91 PetscErrorCode ierr; 92 PetscInt i,e,v,eStart,eEnd,vStart,vEnd,offset,nvar,len_loc,len,k; 93 const PetscScalar *xv; 94 MPI_Comm comm; 95 PetscMPIInt size,rank,tag = ((PetscObject)viewer)->tag; 96 Vec localX; 97 PetscBool ghostvtex; 98 PetscScalar *values; 99 PetscInt j,ne,nv,id; 100 MPI_Status status; 101 102 PetscFunctionBegin; 103 ierr = PetscObjectGetComm((PetscObject)networkdm,&comm);CHKERRQ(ierr); 104 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 105 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 106 107 ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr); 108 ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 109 ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 110 ierr = VecGetArrayRead(localX,&xv);CHKERRQ(ierr); 111 112 ierr = VecGetLocalSize(localX,&len_loc);CHKERRQ(ierr); 113 114 ierr = DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);CHKERRQ(ierr); 115 ierr = DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);CHKERRQ(ierr); 116 len_loc += 2*(1 + eEnd-eStart + vEnd-vStart); 117 118 /* values = [nedges, nvertices; id, nvar, xedge; ...; id, nvars, xvertex;...], to be sent to proc[0] */ 119 ierr = MPI_Allreduce(&len_loc,&len,1,MPIU_INT,MPI_MAX,comm);CHKERRQ(ierr); 120 ierr = PetscCalloc1(len,&values);CHKERRQ(ierr); 121 122 if (!rank) { 123 ierr = PetscViewerASCIIPrintf(viewer,"Process [%d]\n",rank);CHKERRQ(ierr); 124 } 125 126 /* iterate over edges */ 127 k = 2; 128 for (e=eStart; e<eEnd; e++) { 129 ierr = DMNetworkGetNumVariables(networkdm,e,&nvar);CHKERRQ(ierr); 130 if (!nvar) continue; 131 132 ierr = DMNetworkGetVariableOffset(networkdm,e,&offset);CHKERRQ(ierr); 133 ierr = DMNetworkGetGlobalEdgeIndex(networkdm,e,&id);CHKERRQ(ierr); 134 135 if (!rank) { /* print its own entries */ 136 ierr = PetscViewerASCIIPrintf(viewer," Edge %D:\n",id);CHKERRQ(ierr); 137 ierr = VecArrayPrint_private(viewer,nvar,xv+offset);CHKERRQ(ierr); 138 } else { 139 values[0] += 1; /* number of edges */ 140 values[k++] = id; 141 values[k++] = nvar; 142 for (i=offset; i< offset+nvar; i++) values[k++] = xv[i]; 143 } 144 } 145 146 /* iterate over vertices */ 147 for (v=vStart; v<vEnd; v++) { 148 ierr = DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);CHKERRQ(ierr); 149 if (ghostvtex) continue; 150 ierr = DMNetworkGetNumVariables(networkdm,v,&nvar);CHKERRQ(ierr); 151 if (!nvar) continue; 152 153 ierr = DMNetworkGetVariableOffset(networkdm,v,&offset);CHKERRQ(ierr); 154 ierr = DMNetworkGetGlobalVertexIndex(networkdm,v,&id);CHKERRQ(ierr); 155 156 if (!rank) { 157 ierr = PetscViewerASCIIPrintf(viewer," Vertex %D:\n",id);CHKERRQ(ierr); 158 ierr = VecArrayPrint_private(viewer,nvar,xv+offset);CHKERRQ(ierr); 159 } else { 160 values[1] += 1; /* number of vertices */ 161 values[k++] = id; 162 values[k++] = nvar; 163 for (i=offset; i< offset+nvar; i++) values[k++] = xv[i]; 164 } 165 } 166 167 if (!rank) { 168 /* proc[0] receives and prints messages */ 169 for (j=1; j<size; j++) { 170 ierr = PetscViewerASCIIPrintf(viewer,"Process [%d]\n",j);CHKERRQ(ierr); 171 172 ierr = MPI_Recv(values,(PetscMPIInt)len,MPIU_SCALAR,j,tag,comm,&status);CHKERRQ(ierr); 173 174 ne = (PetscInt)PetscAbsScalar(values[0]); 175 nv = (PetscInt)PetscAbsScalar(values[1]); 176 177 /* print received edges */ 178 k = 2; 179 for (i=0; i<ne; i++) { 180 id = (PetscInt)PetscAbsScalar(values[k++]); 181 nvar = (PetscInt)PetscAbsScalar(values[k++]); 182 ierr = PetscViewerASCIIPrintf(viewer," Edge %D:\n",id);CHKERRQ(ierr); 183 ierr = VecArrayPrint_private(viewer,nvar,values+k);CHKERRQ(ierr); 184 k += nvar; 185 } 186 187 /* print received vertices */ 188 for (i=0; i<nv; i++) { 189 id = (PetscInt)PetscAbsScalar(values[k++]); 190 nvar = (PetscInt)PetscAbsScalar(values[k++]); 191 ierr = PetscViewerASCIIPrintf(viewer," Vertex %D:\n",id);CHKERRQ(ierr); 192 ierr = VecArrayPrint_private(viewer,nvar,values+k);CHKERRQ(ierr); 193 k += nvar; 194 } 195 } 196 } else { 197 /* sends values to proc[0] */ 198 ierr = MPI_Send((void*)values,k,MPIU_SCALAR,0,tag,comm);CHKERRQ(ierr); 199 } 200 201 ierr = PetscFree(values);CHKERRQ(ierr); 202 ierr = VecRestoreArrayRead(localX,&xv);CHKERRQ(ierr); 203 ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr); 204 PetscFunctionReturn(0); 205 } 206 207 PETSC_EXTERN PetscErrorCode VecView_MPI(Vec,PetscViewer); 208 209 PetscErrorCode VecView_Network(Vec v,PetscViewer viewer) 210 { 211 DM dm; 212 PetscErrorCode ierr; 213 PetscBool isseq; 214 PetscBool iascii; 215 216 PetscFunctionBegin; 217 ierr = VecGetDM(v,&dm);CHKERRQ(ierr); 218 if (!dm) SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_WRONG,"Vector not generated from a DM"); 219 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 220 ierr = PetscObjectTypeCompare((PetscObject)v,VECSEQ,&isseq);CHKERRQ(ierr); 221 222 /* Use VecView_Network if the viewer is ASCII; use VecView_Seq/MPI for other viewer formats */ 223 if (iascii) { 224 if (isseq) { 225 ierr = VecView_Network_Seq(dm,v,viewer);CHKERRQ(ierr); 226 } else { 227 ierr = VecView_Network_MPI(dm,v,viewer);CHKERRQ(ierr); 228 } 229 } else { 230 if (isseq) { 231 ierr = VecView_Seq(v,viewer);CHKERRQ(ierr); 232 } else { 233 ierr = VecView_MPI(v,viewer);CHKERRQ(ierr); 234 } 235 } 236 PetscFunctionReturn(0); 237 } 238 239 static PetscErrorCode DMCreateGlobalVector_Network(DM dm,Vec *vec) 240 { 241 PetscErrorCode ierr; 242 DM_Network *network = (DM_Network*) dm->data; 243 244 PetscFunctionBegin; 245 ierr = DMCreateGlobalVector(network->plex,vec);CHKERRQ(ierr); 246 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Network);CHKERRQ(ierr); 247 ierr = VecSetDM(*vec,dm);CHKERRQ(ierr); 248 PetscFunctionReturn(0); 249 } 250 251 static PetscErrorCode DMCreateLocalVector_Network(DM dm,Vec *vec) 252 { 253 PetscErrorCode ierr; 254 DM_Network *network = (DM_Network*) dm->data; 255 256 PetscFunctionBegin; 257 ierr = DMCreateLocalVector(network->plex,vec);CHKERRQ(ierr); 258 ierr = VecSetDM(*vec,dm);CHKERRQ(ierr); 259 PetscFunctionReturn(0); 260 } 261 262 PetscErrorCode DMInitialize_Network(DM dm) 263 { 264 265 PetscFunctionBegin; 266 267 dm->ops->view = NULL; 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 = 0; 277 dm->ops->creatematrix = DMCreateMatrix_Network; 278 dm->ops->createinterpolation = 0; 279 dm->ops->getaggregates = 0; 280 dm->ops->getinjection = 0; 281 dm->ops->refine = 0; 282 dm->ops->coarsen = 0; 283 dm->ops->refinehierarchy = 0; 284 dm->ops->coarsenhierarchy = 0; 285 dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Network; 286 dm->ops->globaltolocalend = DMGlobalToLocalEnd_Network; 287 dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Network; 288 dm->ops->localtoglobalend = DMLocalToGlobalEnd_Network; 289 dm->ops->destroy = DMDestroy_Network; 290 dm->ops->createsubdm = NULL; 291 dm->ops->locatepoints = NULL; 292 PetscFunctionReturn(0); 293 } 294 295 PetscErrorCode DMClone_Network(DM dm, DM *newdm) 296 { 297 DM_Network *network = (DM_Network *) dm->data; 298 PetscErrorCode ierr; 299 300 PetscFunctionBegin; 301 network->refct++; 302 (*newdm)->data = network; 303 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMNETWORK);CHKERRQ(ierr); 304 ierr = DMInitialize_Network(*newdm);CHKERRQ(ierr); 305 PetscFunctionReturn(0); 306 } 307 308 /*MC 309 DMNETWORK = "network" - A DM object that encapsulates an unstructured network. The implementation is based on the DM object 310 DMPlex that manages unstructured grids. Distributed networks use a non-overlapping partitioning of 311 the edges. In the local representation, Vecs contain all unknowns in the interior and shared boundary. 312 This is specified by a PetscSection object. Ownership in the global representation is determined by 313 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 314 315 Level: intermediate 316 317 .seealso: DMType, DMNetworkCreate(), DMCreate(), DMSetType() 318 M*/ 319 320 PETSC_EXTERN PetscErrorCode DMCreate_Network(DM dm) 321 { 322 DM_Network *network; 323 PetscErrorCode ierr; 324 325 PetscFunctionBegin; 326 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 327 ierr = PetscNewLog(dm,&network);CHKERRQ(ierr); 328 dm->data = network; 329 330 network->refct = 1; 331 network->NVertices = 0; 332 network->NEdges = 0; 333 network->nVertices = 0; 334 network->nEdges = 0; 335 network->nsubnet = 0; 336 337 338 ierr = DMInitialize_Network(dm);CHKERRQ(ierr); 339 PetscFunctionReturn(0); 340 } 341 342 /*@ 343 DMNetworkCreate - Creates a DMNetwork object, which encapsulates an unstructured network. 344 345 Collective on MPI_Comm 346 347 Input Parameter: 348 . comm - The communicator for the DMNetwork object 349 350 Output Parameter: 351 . network - The DMNetwork object 352 353 Level: beginner 354 355 .keywords: DMNetwork, create 356 @*/ 357 PetscErrorCode DMNetworkCreate(MPI_Comm comm, DM *network) 358 { 359 PetscErrorCode ierr; 360 361 PetscFunctionBegin; 362 PetscValidPointer(network,2); 363 ierr = DMCreate(comm, network);CHKERRQ(ierr); 364 ierr = DMSetType(*network, DMNETWORK);CHKERRQ(ierr); 365 PetscFunctionReturn(0); 366 } 367