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