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